!KGEN-generated Fortran source file !Generated at : 2020-04-03 14:07:28 !KGEN version : 0.9.0 ! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) ! and the University Corporation for Atmospheric Research (UCAR). ! Unless noted otherwise source code is licensed under the BSD license. ! Additional copyright and license information can be found in the LICENSE file ! distributed with this code, or at http://mpas-dev.github.com/license.html ! ! module ocn_gm USE mpas_constants USE ocn_constants USE kgen_utils_mod, ONLY: kgen_dp, kgen_array_sumcheck USE tprof_mod, ONLY: tstart, tstop, tnull, tprnt IMPLICIT NONE PRIVATE SAVE !-------------------------------------------------------------------- ! Public parameters !-------------------------------------------------------------------- !-------------------------------------------------------------------- ! Public member functions !-------------------------------------------------------------------- ! ! ! ! PUBLIC ocn_gm_compute_bolus_velocity !-------------------------------------------------------------------- ! Private module variables !-------------------------------------------------------------------- ! ! PRIVATE tridiagonal_solve ! Config options real (kind=RKIND), pointer :: config_gravWaveSpeed_trunc real (kind=RKIND), pointer :: config_max_relative_slope logical, pointer :: config_disable_redi_k33 logical, pointer :: config_use_Redi_surface_layer_tapering logical, pointer :: config_use_Redi_bottom_layer_tapering real (kind=RKIND), pointer :: config_Redi_surface_layer_tapering_extent real (kind=RKIND), pointer :: config_Redi_bottom_layer_tapering_depth logical, pointer :: config_gm_lat_variable_c2 logical, pointer :: config_gm_kappa_lat_depth_variable real (kind=RKIND), pointer :: config_gm_min_stratification_ratio real (kind=RKIND), pointer :: config_gm_min_phase_speed real (kind=RKIND), parameter :: epsGM = 1.0e-12_RKIND !*********************************************************************** #ifdef _MPI include "mpif.h" #endif PUBLIC kr_externs_in_ocn_gm PUBLIC kr_externs_out_ocn_gm contains !*********************************************************************** ! routine ocn_gm_compute_Bolus_velocity !> \brief Computes GM Bolus velocity !> \author Qingshan Chen, Mark Petersen, Todd Ringler !> \date January 2013 !> \details !> This routine is the main driver for the Gent-McWilliams (GM) parameterization. !> It computes horizontal and vertical density gradients, the slope !> of isopycnal surfaces, and solves a boundary value problem in each column !> for the stream function, which is used to compute the Bolus velocity. !----------------------------------------------------------------------- ! ! ! SUBROUTINE ocn_gm_compute_bolus_velocity(kgen_unit, kgen_measure, kgen_isverified, kgen_filepath) !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- USE kgen_utils_mod, ONLY: kgen_dp, kgen_array_sumcheck USE kgen_utils_mod, ONLY: kgen_perturb_real USE kgen_utils_mod, ONLY: check_t, kgen_init_check, kgen_init_verify, kgen_tolerance, kgen_minvalue, kgen_verboselevel, & &CHECK_IDENTICAL, CHECK_IN_TOL, CHECK_OUT_TOL ! ! !----------------------------------------------------------------- ! input/output variables !----------------------------------------------------------------- ! ! !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- ! ! REAL(KIND=rkind), dimension(:,:), pointer :: density, displaceddensity, zmid, normalgmbolusvelocity, layerthicknessedge, & &graddensityedge, graddensitytopofedge, graddensityconstztopofedge, gradzmidedge, gradzmidtopofedge, relativeslopetopofedge, & &relativeslopetopofcell, k33, gmstreamfunctopofedge, bruntvaisalafreqtop, gmstreamfunctopofcell, ddensitydztopofedge, & &ddensitydztopofcell, relativeslopetapering, relativeslopetaperingcell, areacellsum, kappagm3d REAL(KIND=rkind), dimension(:), pointer :: boundarylayerdepth, gmboluskappa, cgmphasespeed REAL(KIND=rkind), dimension(:), pointer :: areacell, dcedge, dvedge, tridiaga, tridiagb, tridiagc, righthandside INTEGER, dimension(:), pointer :: maxleveledgetop, maxlevelcell, nedgesoncell INTEGER, dimension(:,:), pointer :: cellsonedge, edgesoncell INTEGER :: i, k, iedge, cell1, cell2, icell, n, iter REAL(KIND=rkind) :: h1, h2, areaedge, bruntvaisalafreqtopedge, rtmp, stmp REAL(KIND=rkind) :: sumn2, countn2, maxn, ltsum ! Dimensions INTEGER :: ncells, nedges INTEGER, pointer :: nvertlevels INTEGER, dimension(:), pointer :: ncellsarray, nedgesarray INTEGER, INTENT(IN) :: kgen_unit REAL(KIND=kgen_dp), INTENT(OUT) :: kgen_measure LOGICAL, INTENT(OUT) :: kgen_isverified CHARACTER(LEN=*), INTENT(IN) :: kgen_filepath LOGICAL :: kgen_istrue REAL(KIND=8) :: kgen_array_sum INTEGER :: kgen_intvar, kgen_ierr INTEGER :: kgen_mpirank, kgen_openmptid, kgen_kernelinvoke LOGICAL :: kgen_evalstage, kgen_warmupstage, kgen_mainstage COMMON / state / kgen_mpirank, kgen_openmptid, kgen_kernelinvoke, kgen_evalstage, kgen_warmupstage, kgen_mainstage INTEGER, PARAMETER :: KGEN_MAXITER = NUM_REPEAT TYPE(check_t) :: check_status INTEGER*8 :: kgen_start_clock, kgen_stop_clock, kgen_rate_clock REAL(KIND=kgen_dp) :: gkgen_measure REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_graddensityedge REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_gradzmidedge REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_normalgmbolusvelocity REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_graddensitytopofedge REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_ddensitydztopofedge REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_gradzmidtopofedge REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_relativeslopetopofedge REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_relativeslopetapering REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_ddensitydztopofcell REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_k33 REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_relativeslopetopofcell REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_relativeslopetaperingcell REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_graddensityconstztopofedge REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_areacellsum REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_kappagm3d REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_gmstreamfunctopofedge REAL(KIND=rkind), pointer, dimension(:,:) :: kgenref_gmstreamfunctopofcell REAL(KIND=rkind), pointer, dimension(:) :: kgenref_cgmphasespeed REAL(KIND=rkind), pointer, dimension(:) :: kgenref_tridiagb REAL(KIND=rkind), pointer, dimension(:) :: kgenref_tridiagc REAL(KIND=rkind), pointer, dimension(:) :: kgenref_righthandside REAL(KIND=rkind), pointer, dimension(:) :: kgenref_tridiaga INTEGER, pointer, dimension(:) :: kgenref_maxleveledgetop INTEGER :: kgenref_iedge INTEGER :: kgenref_k INTEGER :: kgenref_icell INTEGER :: kgenref_cell1 INTEGER :: kgenref_cell2 INTEGER :: kgenref_i INTEGER :: kgenref_iter INTEGER :: kgenref_n REAL(KIND=rkind) :: kgenref_rtmp REAL(KIND=rkind) :: kgenref_h1 REAL(KIND=rkind) :: kgenref_h2 REAL(KIND=rkind) :: kgenref_areaedge REAL(KIND=rkind) :: kgenref_stmp REAL(KIND=rkind) :: kgenref_bruntvaisalafreqtopedge REAL(KIND=rkind) :: kgenref_sumn2 REAL(KIND=rkind) :: kgenref_ltsum REAL(KIND=rkind) :: kgenref_countn2 REAL(KIND=rkind) :: kgenref_maxn INTEGER :: kgenref_ncells INTEGER :: kgenref_nedges !parent block preprocessing #ifdef _MPI call mpi_comm_rank(mpi_comm_world, kgen_mpirank, kgen_ierr) #else kgen_mpirank = 0 #endif !local input variables CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(displaceddensity, kgen_unit, "displaceddensity", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(density, kgen_unit, "density", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(zmid, kgen_unit, "zmid", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(ddensitydztopofcell, kgen_unit, "ddensitydztopofcell", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(layerthicknessedge, kgen_unit, "layerthicknessedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(graddensityedge, kgen_unit, "graddensityedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(gradzmidedge, kgen_unit, "gradzmidedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(graddensitytopofedge, kgen_unit, "graddensitytopofedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(ddensitydztopofedge, kgen_unit, "ddensitydztopofedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(gradzmidtopofedge, kgen_unit, "gradzmidtopofedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(relativeslopetopofedge, kgen_unit, "relativeslopetopofedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(relativeslopetopofcell, kgen_unit, "relativeslopetopofcell", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(areacellsum, kgen_unit, "areacellsum", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(relativeslopetaperingcell, kgen_unit, "relativeslopetaperingcell", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(bruntvaisalafreqtop, kgen_unit, "bruntvaisalafreqtop", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kappagm3d, kgen_unit, "kappagm3d", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(graddensityconstztopofedge, kgen_unit, "graddensityconstztopofedge", & &.FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(gmstreamfunctopofedge, kgen_unit, "gmstreamfunctopofedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(gmstreamfunctopofcell, kgen_unit, "gmstreamfunctopofcell", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(normalgmbolusvelocity, kgen_unit, "normalgmbolusvelocity", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(relativeslopetapering, kgen_unit, "relativeslopetapering", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(k33, kgen_unit, "k33", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(boundarylayerdepth, kgen_unit, "boundarylayerdepth", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(gmboluskappa, kgen_unit, "gmboluskappa", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(cgmphasespeed, kgen_unit, "cgmphasespeed", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(dcedge, kgen_unit, "dcedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(dvedge, kgen_unit, "dvedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(tridiaga, kgen_unit, "tridiaga", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(tridiagb, kgen_unit, "tridiagb", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(tridiagc, kgen_unit, "tridiagc", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(righthandside, kgen_unit, "righthandside", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(areacell, kgen_unit, "areacell", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp2(maxlevelcell, kgen_unit, "maxlevelcell", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp2(maxleveledgetop, kgen_unit, "maxleveledgetop", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp2(nedgesoncell, kgen_unit, "nedgesoncell", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp3(cellsonedge, kgen_unit, "cellsonedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp3(edgesoncell, kgen_unit, "edgesoncell", .FALSE.) READ (UNIT = kgen_unit) k READ (UNIT = kgen_unit) iedge READ (UNIT = kgen_unit) icell READ (UNIT = kgen_unit) cell1 READ (UNIT = kgen_unit) cell2 READ (UNIT = kgen_unit) i READ (UNIT = kgen_unit) n READ (UNIT = kgen_unit) iter READ (UNIT = kgen_unit) rtmp READ (UNIT = kgen_unit) h2 READ (UNIT = kgen_unit) h1 READ (UNIT = kgen_unit) areaedge READ (UNIT = kgen_unit) stmp READ (UNIT = kgen_unit) bruntvaisalafreqtopedge READ (UNIT = kgen_unit) sumn2 READ (UNIT = kgen_unit) ltsum READ (UNIT = kgen_unit) countn2 READ (UNIT = kgen_unit) maxn READ (UNIT = kgen_unit) nedges READ (UNIT = kgen_unit) ncells CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp5(nvertlevels, kgen_unit, "nvertlevels", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp2(ncellsarray, kgen_unit, "ncellsarray", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp2(nedgesarray, kgen_unit, "nedgesarray", .FALSE.) !extern output variables CALL kr_externs_out_ocn_gm(kgen_unit) !local output variables CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_graddensityedge, kgen_unit, "kgenref_graddensityedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_gradzmidedge, kgen_unit, "kgenref_gradzmidedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_normalgmbolusvelocity, kgen_unit, "kgenref_normalgmbolusvelocity", & &.FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_graddensitytopofedge, kgen_unit, "kgenref_graddensitytopofedge", & &.FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_ddensitydztopofedge, kgen_unit, "kgenref_ddensitydztopofedge", & &.FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_gradzmidtopofedge, kgen_unit, "kgenref_gradzmidtopofedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_relativeslopetopofedge, kgen_unit, "kgenref_relativeslopetopofedge", & &.FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_relativeslopetapering, kgen_unit, "kgenref_relativeslopetapering", & &.FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_ddensitydztopofcell, kgen_unit, "kgenref_ddensitydztopofcell", & &.FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_k33, kgen_unit, "kgenref_k33", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_relativeslopetopofcell, kgen_unit, "kgenref_relativeslopetopofcell", & &.FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_relativeslopetaperingcell, kgen_unit, & &"kgenref_relativeslopetaperingcell", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_graddensityconstztopofedge, kgen_unit, & &"kgenref_graddensityconstztopofedge", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_areacellsum, kgen_unit, "kgenref_areacellsum", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_kappagm3d, kgen_unit, "kgenref_kappagm3d", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_gmstreamfunctopofedge, kgen_unit, "kgenref_gmstreamfunctopofedge", & &.FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp0(kgenref_gmstreamfunctopofcell, kgen_unit, "kgenref_gmstreamfunctopofcell", & &.FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(kgenref_cgmphasespeed, kgen_unit, "kgenref_cgmphasespeed", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(kgenref_tridiagb, kgen_unit, "kgenref_tridiagb", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(kgenref_tridiagc, kgen_unit, "kgenref_tridiagc", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(kgenref_righthandside, kgen_unit, "kgenref_righthandside", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp1(kgenref_tridiaga, kgen_unit, "kgenref_tridiaga", .FALSE.) CALL kr_kgen_ocn_gm_compute_bolus_velocity_subp2(kgenref_maxleveledgetop, kgen_unit, "kgenref_maxleveledgetop", .FALSE.) READ (UNIT = kgen_unit) kgenref_iedge READ (UNIT = kgen_unit) kgenref_k READ (UNIT = kgen_unit) kgenref_icell READ (UNIT = kgen_unit) kgenref_cell1 READ (UNIT = kgen_unit) kgenref_cell2 READ (UNIT = kgen_unit) kgenref_i READ (UNIT = kgen_unit) kgenref_iter READ (UNIT = kgen_unit) kgenref_n READ (UNIT = kgen_unit) kgenref_rtmp READ (UNIT = kgen_unit) kgenref_h1 READ (UNIT = kgen_unit) kgenref_h2 READ (UNIT = kgen_unit) kgenref_areaedge READ (UNIT = kgen_unit) kgenref_stmp READ (UNIT = kgen_unit) kgenref_bruntvaisalafreqtopedge READ (UNIT = kgen_unit) kgenref_sumn2 READ (UNIT = kgen_unit) kgenref_ltsum READ (UNIT = kgen_unit) kgenref_countn2 READ (UNIT = kgen_unit) kgenref_maxn READ (UNIT = kgen_unit) kgenref_ncells READ (UNIT = kgen_unit) kgenref_nedges !$kgen begin_callsite ocn_gm_compute_Bolus_velocity IF (kgen_evalstage) THEN END IF IF (kgen_warmupstage) THEN END IF IF (kgen_mainstage) THEN END IF !Uncomment following call statement to turn on perturbation experiment. !Adjust perturbation value and/or kind parameter if required. !CALL kgen_perturb_real( your_variable, 1.0D-15_8 ) !call to kgen kernel !$acc data copyin(nvertlevels) nCells = nCellsArray( size(nCellsArray) ) nEdges = nEdgesArray( size(nEdgesArray) ) ! Assign a huge value to the scratch variables which may manifest itself when ! there is a bug. !$omp do schedule(runtime) private(k) !print *, "NEDGES", nEdges !print *, "NVERTLEVELS", nvertlevels !$acc parallel loop gang do iEdge = 1, nEdges !$acc loop vector do k = 1, nVertLevels gradDensityEdge(k, iEdge) = huge(0D0) gradZMidEdge(k, iEdge) = huge(0D0) normalGMBolusVelocity(k, iEdge) = 0.0_RKIND end do end do !$acc end parallel !$omp end do !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iEdge = 1, nEdges !$acc loop vector do k = 1, nVertLevels + 1 gradDensityTopOfEdge(k, iEdge) = huge(0D0) dDensityDzTopOfEdge(k, iEdge) = huge(0D0) gradZMidTopOfEdge(k, iEdge) = huge(0D0) relativeSlopeTopOfEdge(k, iEdge) = 0.0_RKIND relativeSlopeTapering(k, iEdge) = 0.0_RKIND end do end do !$acc end parallel !$omp end do !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iCell = 1, nCells + 1 !$acc loop vector do k = 1, nVertLevels dDensityDzTopOfCell(k, iCell) = huge(0D0) k33(k, iCell) = 0.0_RKIND relativeSlopeTopOfCell(k, iCell) = 0.0_RKIND relativeSlopeTaperingCell(k, iCell) = 0.0_RKIND end do end do !$acc end parallel !$omp end do !-------------------------------------------------------------------- ! ! Compute vertical derivative of density at top of cell, interpolate to top of edge ! This is required for Redi and Bolus parts. ! !-------------------------------------------------------------------- nCells = nCellsArray( 3 ) ! Compute vertical derivative of density (dDensityDzTopOfCell) at cell center and layer interface ! Note that displacedDensity is used from the upper cell, so that the EOS reference level for ! pressure is the same for both displacedDensity(k-1,iCell) and density(k,iCell). !$omp do schedule(runtime) private(k, rtmp) ! verification failure !!$acc parallel loop gang do iCell = 1, nCells !!$acc loop vector do k = 2, maxLevelCell(iCell) rtmp = (displacedDensity(k-1,iCell) - density(k,iCell)) / (zMid(k-1,iCell) - zMid(k,iCell)) dDensityDzTopOfCell(k,iCell) = min(rtmp, -epsGM) end do ! Approximation of dDensityDzTopOfCell on the top and bottom interfaces through the idea of having ! ghost cells above the top and below the bottom layers of the same depths and density. ! Essentially, this enforces the boundary condition (d density)/dz = 0 at the top and bottom. dDensityDzTopOfCell(1,iCell) = 0.0_RKIND dDensityDzTopOfCell(maxLevelCell(iCell)+1,iCell) = 0.0_RKIND end do !!$acc end parallel !$omp end do nEdges = nEdgesArray( 3 ) ! Interpolate dDensityDzTopOfCell to edge and layer interface !$omp do schedule(runtime) private(k, cell1, cell2) !$acc parallel loop gang do iEdge = 1, nEdges !$acc loop vector do k = 1, maxLevelEdgeTop(iEdge)+1 cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) dDensityDzTopOfEdge(k,iEdge) = 0.5_RKIND * (dDensityDzTopOfCell(k,cell1) + dDensityDzTopOfCell(k,cell2)) end do end do !$acc end parallel !$omp end do !-------------------------------------------------------------------- ! ! Compute horizontal gradient and mid-layer of edge, interpolate to top of edge ! This is required for Redi and Bolus parts. ! !-------------------------------------------------------------------- nEdges = nEdgesArray( 3 ) ! Compute density gradient (gradDensityEdge) and gradient of zMid (gradZMidEdge) ! along the constant coordinate surface. ! The computed variables lives at edge and mid-layer depth !$omp do schedule(runtime) private(cell1, cell2, k) ! verification failure !!$acc parallel loop gang do iEdge = 1, nEdges cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) !!$acc loop vector do k=1,maxLevelEdgeTop(iEdge) gradDensityEdge(k,iEdge) = (density(k,cell2) - density(k,cell1)) / dcEdge(iEdge) gradZMidEdge(k,iEdge) = (zMid(k,cell2) - zMid(k,cell1)) / dcEdge(iEdge) end do end do !!$acc end parallel !$omp end do nEdges = nEdgesArray( 3 ) ! Interpolate gradDensityEdge and gradZMidEdge to layer interface !$omp do schedule(runtime) private(k, h1, h2) ! verification failure !!$acc parallel loop gang do iEdge = 1, nEdges ! The interpolation can only be carried out on non-boundary edges if (maxLevelEdgeTop(iEdge) .GE. 1) then !!$acc loop vector do k = 2, maxLevelEdgeTop(iEdge) h1 = layerThicknessEdge(k-1,iEdge) h2 = layerThicknessEdge(k,iEdge) ! Using second-order interpolation below gradDensityTopOfEdge(k,iEdge) = (h2 * gradDensityEdge(k-1,iEdge) + h1 * gradDensityEdge(k,iEdge)) / (h1 + h2) gradZMidTopOfEdge(k,iEdge) = (h2 * gradZMidEdge(k-1,iEdge) + h1 * gradZMidEdge(k,iEdge)) / (h1 + h2) end do ! Approximation of values on the top and bottom interfaces through the idea of having ghost cells above ! the top and below the bottom layers of the same depths and density. gradDensityTopOfEdge(1,iEdge) = gradDensityEdge(1,iEdge) gradDensityTopOfEdge(maxLevelEdgeTop(iEdge)+1,iEdge) = gradDensityEdge(maxLevelEdgeTop(iEdge),iEdge) gradZMidTopOfEdge(1,iEdge) = gradZMidEdge(1,iEdge) gradZMidTopOfEdge(maxLevelEdgeTop(iEdge)+1,iEdge) = gradZMidEdge(maxLevelEdgeTop(iEdge),iEdge) end if end do !!$acc end parallel !$omp end do !-------------------------------------------------------------------- ! ! Compute horizontal gradient required for Bolus part (along constant z) ! !-------------------------------------------------------------------- nEdges = nEdgesArray( 3 ) !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iEdge = 1, nEdges if (maxLevelEdgeTop(iEdge) .GE. 1) then !$acc loop vector do k = 1, maxLevelEdgeTop(iEdge)+1 gradDensityConstZTopOfEdge(k,iEdge) = gradDensityTopOfEdge(k,iEdge) - dDensityDzTopOfEdge(k,iEdge) & * gradZMidTopOfEdge(k,iEdge) end do end if end do !$acc end parallel !$omp end do !-------------------------------------------------------------------- ! ! Compute relative slope and k33 for Redi part of GM. ! These variables are used in del2 velocity tendency routines. ! !-------------------------------------------------------------------- nEdges = nEdgesArray( 3 ) ! Compute relativeSlopeTopOfEdge at edge and layer interface ! set relativeSlopeTopOfEdge to zero for horizontal land/water edges. !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iEdge = 1, nEdges relativeSlopeTopOfEdge(:, iEdge) = 0.0_RKIND ! Beside a full land cell (e.g. missing cell) maxLevelEdgeTop=0, so relativeSlopeTopOfEdge at that edge will remain zero. !$acc loop vector do k = 2, maxLevelEdgeTop(iEdge) relativeSlopeTopOfEdge(k,iEdge) = - gradDensityTopOfEdge(k,iEdge) / min(dDensityDzTopOfEdge(k,iEdge),-epsGM) end do ! Since dDensityDzTopOfEdge is guaranteed to be zero on the top surface, relativeSlopeTopOfEdge on the top ! surface is identified with its value on the second interface. relativeSlopeTopOfEdge(1,iEdge) = relativeSlopeTopOfEdge(2,iEdge) ! dDensityDzTopOfEdge may or may not equal zero on the bottom surface, depending on whether ! maxLevelEdgeTop(iEdge) = maxLevelEdgeBottom(iEdge). But here we ! take a simplistic approach and identify relativeSlopeTopOfEdge on the bottom surface with its value on ! the interface just above. relativeSlopeTopOfEdge( maxLevelEdgeTop(iEdge)+1, iEdge ) = relativeSlopeTopOfEdge( max(1,maxLevelEdgeTop(iEdge)), iEdge ) end do !$acc end parallel !$omp end do nEdges = nEdgesArray( 3 ) ! slope can be unbounded in regions of neutral stability, reset to the large, but bounded, value ! values is hardwrite to 1.0, this is equivalent to a slope of 45 degrees !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iEdge = 1, nEdges !$acc loop vector do k = 1, nVertLevels relativeSlopeTopOfEdge(k, iEdge) = max( min( relativeSlopeTopOfEdge(k, iEdge), 1.0_RKIND), -1.0_RKIND) end do end do !$acc end parallel !$omp end do ! average relative slope to cell centers ! do this by computing (relative slope)^2, then taking sqrt nCells = nCellsArray( 2 ) !$omp do schedule(runtime) private(i, iEdge, areaEdge, rtmp, k) !$acc parallel loop gang do iCell = 1, nCells areaCellSum(:, iCell) = 1.0e-34_RKIND do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) !contribution of cell area from this edge * 2.0 areaEdge = 0.5_RKIND * dcEdge(iEdge) * dvEdge(iEdge) !$acc loop vector do k = 1, maxLevelEdgeTop(iEdge) rtmp = areaEdge * relativeSlopeTopOfEdge(k, iEdge)**2 relativeSlopeTopOfCell(k, iCell) = relativeSlopeTopOfCell(k, iCell) + rtmp areaCellSum(k, iCell) = areaCellSum(k, iCell) + areaEdge end do end do end do !$acc end parallel !$omp end do nCells = nCellsArray( 2 ) !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iCell=1,nCells !$acc loop vector do k = 1, maxLevelCell(iCell) relativeSlopeTopOfCell(k,iCell) = sqrt( relativeSlopeTopOfCell(k,iCell)/areaCellSum(k,iCell) ) end do end do !$acc end parallel !$omp end do ! Compute tapering function ! Compute k33 at cell center and layer interface nCells = nCellsArray( size(nCellsArray) ) !$omp do schedule(runtime) ! verification failure !!$acc parallel loop gang do iCell = 1, nCells k33(:, iCell) = 0.0_RKIND end do !!$acc end parallel !$omp end do ! use relativeSlopeTaperingCell as a temporary space for smoothing of relativeSlopeTopOfCell relativeSlopeTaperingCell = relativeSlopeTopOfCell !!$acc parallel loop gang do iter = 1, 5 nCells = nCellsArray( 2 ) !$omp do schedule(runtime) do iCell=1,nCells relativeSlopeTaperingCell(1, iCell) = 0.0_RKIND relativeSlopeTaperingCell(maxLevelCell(iCell):nVertLevels, iCell) = 0.0_RKIND !!$acc loop vector do k = 2, maxLevelCell(iCell)-1 rtmp = relativeSlopeTopOfCell(k-1,iCell) + relativeSlopeTopOfCell(k+1,iCell) stmp = 2.0_RKIND*relativeSlopeTopOfCell(k,iCell) relativeSlopeTaperingCell(k,iCell) = (rtmp+stmp)/4.0_RKIND end do relativeSlopeTopOfCell(:, iCell) = relativeSlopeTaperingCell(:, iCell) end do !$omp end do end do ! iter !!$acc end parallel nCells = nCellsArray ( 2 ) ! first, compute tapering across full domain based on a maximum allowable slope !$omp do schedule(runtime) private(k) !!$acc parallel loop gang do iCell=1,nCells ! compilation error !!$acc loop vector do k = 1, maxLevelCell(iCell) relativeSlopeTaperingCell(k,iCell) = min(1.0_RKIND, config_max_relative_slope / (relativeSlopeTopOfCell(k,iCell)+epsGM)) end do end do !!$acc end parallel !$omp end do ! now further taper in the boundary layer ! vertical (k33) tapering starts at 2*OBL, increases linearly to OBL and is held uniform across OBL ! rtmp = 1 @ zMid = -2.0*OBL, rtmp = 0 @ zMid = -OBL if(config_use_Redi_surface_layer_tapering) then nCells = nCellsArray ( 2 ) !$omp do schedule(runtime) private(k, rtmp) !!$acc parallel loop gang do iCell=1,nCells ! compilation error !!$acc loop vector do k = 1, maxLevelCell(iCell) rtmp = -zMid(k,iCell)/max(config_Redi_surface_layer_tapering_extent,boundaryLayerDepth(iCell)+epsGM) rtmp = max(0.0_RKIND,rtmp) rtmp = min(1.0_RKIND,rtmp) relativeSlopeTaperingCell(k,iCell) = rtmp*relativeSlopeTaperingCell(k,iCell) end do end do !!$acc end parallel !$omp end do endif ! config_use_Redi_surface_layer_tapering ! now further taper in the boundary layer ! vertical (k33) tapering starts at 2*OBL, increases linearly to OBL and is held uniform across OBL ! rtmp = 1 @ zMid = zMid(maxLevelCell) + config_Redi_bottom_layer_tapering_depth, rtmp = 0 @ zMid = zMid(maxLevelCell) if(config_use_Redi_bottom_layer_tapering) then nCells = nCellsArray ( 2 ) !$omp do schedule(runtime) private(k, rtmp) do iCell=1,nCells do k = 1, maxLevelCell(iCell) rtmp = (zMid(k,iCell)-zMid(maxLevelCell(iCell),iCell))/(config_Redi_bottom_layer_tapering_depth+epsGM) rtmp = max(0.0_RKIND,rtmp) rtmp = min(1.0_RKIND,rtmp) relativeSlopeTaperingCell(k,iCell) = rtmp*relativeSlopeTaperingCell(k,iCell) end do end do !$omp end do endif ! config_use_Redi_bottom_layer_tapering nCells = nCellsArray( 2 ) !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iCell=1,nCells k33(:, iCell) = 0.0_RKIND !$acc loop vector do k = 2, maxLevelCell(iCell) k33(k,iCell) = ( relativeSlopeTaperingCell(k,iCell) * relativeSlopeTopOfCell(k,iCell) )**2 end do end do !$acc end parallel !$omp end do nEdges = nEdgesArray( 3 ) ! average tapering function to layer edges !$omp do schedule(runtime) private(cell1, cell2, k) !$acc parallel loop gang do iEdge = 1, nEdges cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) !$acc loop vector do k = 1, maxLevelEdgeTop(iEdge) relativeSlopeTapering(k,iEdge) = 0.5_RKIND * (relativeSlopeTaperingCell(k,cell1) + relativeSlopeTaperingCell(k,cell2)) enddo enddo !$acc end parallel !$omp end do ! allow disabling of K33 for testing if(config_disable_redi_k33) then nCells = nCellsArray( size(nCellsArray) ) !$omp do schedule(runtime) !$acc parallel loop gang do iCell = 1, nCells k33(:, iCell) = 0.0_RKIND end do !$acc end parallel !$omp end do end if !-------------------------------------------------------------------- ! ! Compute stream function and Bolus velocity for Bolus part of GM ! !-------------------------------------------------------------------- if (config_gm_lat_variable_c2) then !$omp do schedule(runtime) private(cell1, cell2, sumN2, ltSum, countN2, BruntVaisalaFreqTopEdge) ! compilation error !!$acc parallel loop gang do iEdge = 1, nEdges cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) sumN2 = 0.0 ltSum = 0.0 countN2 = 0 !!$acc loop vector do k=2,maxLevelEdgeTop(iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) sumN2 = sumN2 + BruntVaisalaFreqTopEdge*layerThicknessEdge(k,iEdge) ltSum = ltSum + layerThicknessEdge(k,iEdge) countN2 = countN2 +1 enddo if(countN2 > 0) cGMphaseSpeed(iEdge) = max(config_gm_min_phase_speed ,sqrt(sumN2/ltSum)*ltSum / 3.141592_RKIND) enddo !!$acc end parallel !$omp end do else !$omp do schedule(runtime) ! compilation error !!$acc parallel loop gang do iEdge = 1, nEdges cGMphaseSpeed(iEdge) = config_gravWaveSpeed_trunc enddo !!$acc end parallel !$omp end do endif !$omp do schedule(runtime) ! runtime failure !!$acc parallel loop gang do iEdge=1,nEdges kappaGM3D(:,iEdge) = gmBolusKappa(iEdge) enddo !!$acc end parallel !$omp end do if (config_gm_kappa_lat_depth_variable) then !$omp do schedule(runtime) private(cell1, cell2, k, BruntVaisalaFreqTopEdge, maxN) ! compilation error !!$acc parallel loop gang do iEdge = 1,nEdges cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) maxN = -1.0_RKIND do k=2,maxLevelEdgeTop(iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) maxN = max(maxN,BruntVaisalaFreqTopEdge) enddo do k=2,maxLevelEdgeTop(iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) kappaGM3D(k,iEdge) = gmBolusKappa(iEdge)*max(config_gm_min_stratification_ratio, & BruntVaisalaFreqTopEdge / (maxN + 1.0E-10_RKIND)) enddo enddo !!$acc end parallel !$omp end do endif nEdges = nEdgesArray( 3 ) !$omp do schedule(runtime) ! compilation error !!$acc parallel loop gang do iEdge = 1, nEdges cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) gmStreamFuncTopOfEdge(:, iEdge) = 0.0_RKIND ! Construct the tridiagonal matrix if (maxLevelEdgeTop(iEdge) .GE. 3) then ! First row k = 2 BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1,iEdge) & * layerThicknessEdge(k,iEdge)) - BruntVaisalaFreqTopEdge tridiagC(k-1) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k, iEdge) & / (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge)) rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge) ! Second to next to the last rows do k = 3, maxLevelEdgeTop(iEdge)-1 BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) tridiagA(k-2) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k-1, iEdge) & / (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge)) tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1, iEdge) & * layerThicknessEdge(k, iEdge) ) - BruntVaisalaFreqTopEdge tridiagC(k-1) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k, iEdge) & / (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge)) rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge) end do ! Last row k = maxLevelEdgeTop(iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) tridiagA(k-2) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k-1,iEdge) & / (layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge)) tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1, iEdge) & * layerThicknessEdge(k, iEdge)) - BruntVaisalaFreqTopEdge rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge) ! Total number of rows N = maxLevelEdgeTop(iEdge) - 1 ! Call the tridiagonal solver call tridiagonal_solve(tridiagA, tridiagB, tridiagC, rightHandSide, & gmStreamFuncTopOfEdge(2:maxLevelEdgeTop(iEdge), iEdge), N) end if end do !!$acc end parallel !$omp end do nEdges = nEdgesArray( 3 ) ! Compute normalGMBolusVelocity from the stream function !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iEdge = 1, nEdges !$acc loop vector do k = 1, maxLevelEdgeTop(iEdge) normalGMBolusVelocity(k,iEdge) = (gmStreamFuncTopOfEdge(k,iEdge) - gmStreamFuncTopOfEdge(k+1,iEdge)) & / layerThicknessEdge(k,iEdge) end do end do !$acc end parallel !$omp end do nCells = nCellsArray( 1 ) ! Interpolate gmStreamFuncTopOfEdge to cell centers for visualization !$omp do schedule(runtime) private(i, iEdge, areaEdge, k, rtmp) ! verification failure !!$acc parallel loop gang do iCell = 1, nCells gmStreamFuncTopOfCell(:, iCell) = 0.0_RKIND do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) areaEdge = 0.25_RKIND * dcEdge(iEdge) * dvEdge(iEdge) !!$acc loop vector do k = 1, maxLevelEdgeTop(iEdge) rtmp = 0.5_RKIND * ( gmStreamFuncTopOfEdge(k, iEdge) + gmStreamFuncTopOfEdge(k+1, iEdge) ) * areaEdge gmStreamFuncTopOfCell(k, iCell) = gmStreamFuncTopOfCell(k, iCell) + rtmp end do end do end do !!$acc end parallel !$omp end do !$omp do schedule(runtime) ! execution failture !!$acc parallel loop gang do iCell = 1, nCells gmStreamFuncTopOfCell(:, iCell) = gmStreamFuncTopOfCell(:,iCell) / areaCell(iCell) end do !!$acc end parallel IF (kgen_mainstage) THEN !verify init CALL kgen_init_verify(tolerance=MAX_TOL, minvalue=1.D-14, verboseLevel=VERBOSITY) CALL kgen_init_check(check_status, rank=kgen_mpirank) !extern verify variables !local verify variables CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("graddensityedge", check_status, graddensityedge, & &kgenref_graddensityedge) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("relativeslopetopofedge", check_status, relativeslopetopofedge, & &kgenref_relativeslopetopofedge) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("relativeslopetaperingcell", check_status, relativeslopetaperingcell, & &kgenref_relativeslopetaperingcell) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("k33", check_status, k33, kgenref_k33) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("graddensitytopofedge", check_status, graddensitytopofedge, & &kgenref_graddensitytopofedge) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("normalgmbolusvelocity", check_status, normalgmbolusvelocity, & &kgenref_normalgmbolusvelocity) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("graddensityconstztopofedge", check_status, & &graddensityconstztopofedge, kgenref_graddensityconstztopofedge) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("gmstreamfunctopofedge", check_status, gmstreamfunctopofedge, & &kgenref_gmstreamfunctopofedge) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("kappagm3d", check_status, kappagm3d, kgenref_kappagm3d) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("gradzmidtopofedge", check_status, gradzmidtopofedge, & &kgenref_gradzmidtopofedge) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("ddensitydztopofedge", check_status, ddensitydztopofedge, & &kgenref_ddensitydztopofedge) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("ddensitydztopofcell", check_status, ddensitydztopofcell, & &kgenref_ddensitydztopofcell) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("gmstreamfunctopofcell", check_status, gmstreamfunctopofcell, & &kgenref_gmstreamfunctopofcell) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("relativeslopetapering", check_status, relativeslopetapering, & &kgenref_relativeslopetapering) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("relativeslopetopofcell", check_status, relativeslopetopofcell, & &kgenref_relativeslopetopofcell) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("areacellsum", check_status, areacellsum, kgenref_areacellsum) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp0("gradzmidedge", check_status, gradzmidedge, kgenref_gradzmidedge) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp1("cgmphasespeed", check_status, cgmphasespeed, kgenref_cgmphasespeed) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp1("tridiagb", check_status, tridiagb, kgenref_tridiagb) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp1("tridiagc", check_status, tridiagc, kgenref_tridiagc) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp1("tridiaga", check_status, tridiaga, kgenref_tridiaga) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp1("righthandside", check_status, righthandside, kgenref_righthandside) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp2("maxleveledgetop", check_status, maxleveledgetop, & &kgenref_maxleveledgetop) CALL kv_ocn_gm_compute_bolus_velocity_integer__("i", check_status, i, kgenref_i) CALL kv_ocn_gm_compute_bolus_velocity_integer__("icell", check_status, icell, kgenref_icell) CALL kv_ocn_gm_compute_bolus_velocity_integer__("k", check_status, k, kgenref_k) CALL kv_ocn_gm_compute_bolus_velocity_integer__("cell2", check_status, cell2, kgenref_cell2) CALL kv_ocn_gm_compute_bolus_velocity_integer__("iter", check_status, iter, kgenref_iter) CALL kv_ocn_gm_compute_bolus_velocity_integer__("n", check_status, n, kgenref_n) CALL kv_ocn_gm_compute_bolus_velocity_integer__("cell1", check_status, cell1, kgenref_cell1) CALL kv_ocn_gm_compute_bolus_velocity_integer__("iedge", check_status, iedge, kgenref_iedge) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp3("h2", check_status, h2, kgenref_h2) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp3("stmp", check_status, stmp, kgenref_stmp) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp3("h1", check_status, h1, kgenref_h1) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp3("bruntvaisalafreqtopedge", check_status, bruntvaisalafreqtopedge, & &kgenref_bruntvaisalafreqtopedge) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp3("rtmp", check_status, rtmp, kgenref_rtmp) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp3("areaedge", check_status, areaedge, kgenref_areaedge) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp3("sumn2", check_status, sumn2, kgenref_sumn2) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp3("maxn", check_status, maxn, kgenref_maxn) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp3("ltsum", check_status, ltsum, kgenref_ltsum) CALL kv_kgen_ocn_gm_compute_bolus_velocity_subp3("countn2", check_status, countn2, kgenref_countn2) CALL kv_ocn_gm_compute_bolus_velocity_integer__("nedges", check_status, nedges, kgenref_nedges) CALL kv_ocn_gm_compute_bolus_velocity_integer__("ncells", check_status, ncells, kgenref_ncells) IF (check_status%rank == 0) THEN WRITE (*, *) "" END IF IF (kgen_verboseLevel > 0) THEN IF (check_status%rank == 0) THEN WRITE (*, *) "Number of output variables: ", check_status%numTotal WRITE (*, *) "Number of identical variables: ", check_status%numIdentical WRITE (*, *) "Number of non-identical variables within tolerance: ", check_status%numInTol WRITE (*, *) "Number of non-identical variables out of tolerance: ", check_status%numOutTol WRITE (*, *) "Tolerance: ", kgen_tolerance END IF END IF IF (check_status%rank == 0) THEN WRITE (*, *) "" END IF IF (check_status%numOutTol > 0) THEN IF (check_status%rank == 0) THEN WRITE (*, *) "Verification FAILED with" // TRIM(ADJUSTL(kgen_filepath)) END IF check_status%Passed = .FALSE. kgen_isverified = .FALSE. ELSE IF (check_status%rank == 0) THEN WRITE (*, *) "Verification PASSED with " // TRIM(ADJUSTL(kgen_filepath)) END IF check_status%Passed = .TRUE. kgen_isverified = .TRUE. END IF IF (check_status%rank == 0) THEN WRITE (*, *) "" END IF #ifdef _MPI call mpi_barrier(mpi_comm_world, kgen_ierr) #endif CALL SYSTEM_CLOCK(kgen_start_clock, kgen_rate_clock) DO kgen_intvar = 1, KGEN_MAXITER nCells = nCellsArray( size(nCellsArray) ) nEdges = nEdgesArray( size(nEdgesArray) ) ! Assign a huge value to the scratch variables which may manifest itself when ! there is a bug. !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iEdge = 1, nEdges !$acc loop vector do k = 1, nVertLevels gradDensityEdge(k, iEdge) = huge(0D0) gradZMidEdge(k, iEdge) = huge(0D0) normalGMBolusVelocity(k, iEdge) = 0.0_RKIND end do end do !$acc end parallel !$omp end do !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iEdge = 1, nEdges !$acc loop vector do k = 1, nVertLevels + 1 gradDensityTopOfEdge(k, iEdge) = huge(0D0) dDensityDzTopOfEdge(k, iEdge) = huge(0D0) gradZMidTopOfEdge(k, iEdge) = huge(0D0) relativeSlopeTopOfEdge(k, iEdge) = 0.0_RKIND relativeSlopeTapering(k, iEdge) = 0.0_RKIND end do end do !$acc end parallel !$omp end do !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iCell = 1, nCells + 1 !$acc loop vector do k = 1, nVertLevels dDensityDzTopOfCell(k, iCell) = huge(0D0) k33(k, iCell) = 0.0_RKIND relativeSlopeTopOfCell(k, iCell) = 0.0_RKIND relativeSlopeTaperingCell(k, iCell) = 0.0_RKIND end do end do !$acc end parallel !$omp end do !-------------------------------------------------------------------- ! ! Compute vertical derivative of density at top of cell, interpolate to top of edge ! This is required for Redi and Bolus parts. ! !-------------------------------------------------------------------- nCells = nCellsArray( 3 ) ! Compute vertical derivative of density (dDensityDzTopOfCell) at cell center and layer interface ! Note that displacedDensity is used from the upper cell, so that the EOS reference level for ! pressure is the same for both displacedDensity(k-1,iCell) and density(k,iCell). !$omp do schedule(runtime) private(k, rtmp) do iCell = 1, nCells do k = 2, maxLevelCell(iCell) rtmp = (displacedDensity(k-1,iCell) - density(k,iCell)) / (zMid(k-1,iCell) - zMid(k,iCell)) dDensityDzTopOfCell(k,iCell) = min(rtmp, -epsGM) end do ! Approximation of dDensityDzTopOfCell on the top and bottom interfaces through the idea of having ! ghost cells above the top and below the bottom layers of the same depths and density. ! Essentially, this enforces the boundary condition (d density)/dz = 0 at the top and bottom. dDensityDzTopOfCell(1,iCell) = 0.0_RKIND dDensityDzTopOfCell(maxLevelCell(iCell)+1,iCell) = 0.0_RKIND end do !$omp end do nEdges = nEdgesArray( 3 ) ! Interpolate dDensityDzTopOfCell to edge and layer interface !$omp do schedule(runtime) private(k, cell1, cell2) !$acc parallel loop gang do iEdge = 1, nEdges !$acc loop vector do k = 1, maxLevelEdgeTop(iEdge)+1 cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) dDensityDzTopOfEdge(k,iEdge) = 0.5_RKIND * (dDensityDzTopOfCell(k,cell1) + dDensityDzTopOfCell(k,cell2)) end do end do !$acc end parallel !$omp end do !-------------------------------------------------------------------- ! ! Compute horizontal gradient and mid-layer of edge, interpolate to top of edge ! This is required for Redi and Bolus parts. ! !-------------------------------------------------------------------- nEdges = nEdgesArray( 3 ) ! Compute density gradient (gradDensityEdge) and gradient of zMid (gradZMidEdge) ! along the constant coordinate surface. ! The computed variables lives at edge and mid-layer depth !$omp do schedule(runtime) private(cell1, cell2, k) do iEdge = 1, nEdges cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) do k=1,maxLevelEdgeTop(iEdge) gradDensityEdge(k,iEdge) = (density(k,cell2) - density(k,cell1)) / dcEdge(iEdge) gradZMidEdge(k,iEdge) = (zMid(k,cell2) - zMid(k,cell1)) / dcEdge(iEdge) end do end do !$omp end do nEdges = nEdgesArray( 3 ) ! Interpolate gradDensityEdge and gradZMidEdge to layer interface !$omp do schedule(runtime) private(k, h1, h2) do iEdge = 1, nEdges ! The interpolation can only be carried out on non-boundary edges if (maxLevelEdgeTop(iEdge) .GE. 1) then do k = 2, maxLevelEdgeTop(iEdge) h1 = layerThicknessEdge(k-1,iEdge) h2 = layerThicknessEdge(k,iEdge) ! Using second-order interpolation below gradDensityTopOfEdge(k,iEdge) = (h2 * gradDensityEdge(k-1,iEdge) + h1 * gradDensityEdge(k,iEdge)) / (h1 + h2) gradZMidTopOfEdge(k,iEdge) = (h2 * gradZMidEdge(k-1,iEdge) + h1 * gradZMidEdge(k,iEdge)) / (h1 + h2) end do ! Approximation of values on the top and bottom interfaces through the idea of having ghost cells above ! the top and below the bottom layers of the same depths and density. gradDensityTopOfEdge(1,iEdge) = gradDensityEdge(1,iEdge) gradDensityTopOfEdge(maxLevelEdgeTop(iEdge)+1,iEdge) = gradDensityEdge(maxLevelEdgeTop(iEdge),iEdge) gradZMidTopOfEdge(1,iEdge) = gradZMidEdge(1,iEdge) gradZMidTopOfEdge(maxLevelEdgeTop(iEdge)+1,iEdge) = gradZMidEdge(maxLevelEdgeTop(iEdge),iEdge) end if end do !$omp end do !-------------------------------------------------------------------- ! ! Compute horizontal gradient required for Bolus part (along constant z) ! !-------------------------------------------------------------------- nEdges = nEdgesArray( 3 ) !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iEdge = 1, nEdges if (maxLevelEdgeTop(iEdge) .GE. 1) then !$acc loop vector do k = 1, maxLevelEdgeTop(iEdge)+1 gradDensityConstZTopOfEdge(k,iEdge) = gradDensityTopOfEdge(k,iEdge) - dDensityDzTopOfEdge(k,iEdge) & * gradZMidTopOfEdge(k,iEdge) end do end if end do !$acc end parallel !$omp end do !-------------------------------------------------------------------- ! ! Compute relative slope and k33 for Redi part of GM. ! These variables are used in del2 velocity tendency routines. ! !-------------------------------------------------------------------- nEdges = nEdgesArray( 3 ) ! Compute relativeSlopeTopOfEdge at edge and layer interface ! set relativeSlopeTopOfEdge to zero for horizontal land/water edges. !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iEdge = 1, nEdges relativeSlopeTopOfEdge(:, iEdge) = 0.0_RKIND ! Beside a full land cell (e.g. missing cell) maxLevelEdgeTop=0, so relativeSlopeTopOfEdge at that edge will remain zero. !$acc loop vector do k = 2, maxLevelEdgeTop(iEdge) relativeSlopeTopOfEdge(k,iEdge) = - gradDensityTopOfEdge(k,iEdge) / min(dDensityDzTopOfEdge(k,iEdge),-epsGM) end do ! Since dDensityDzTopOfEdge is guaranteed to be zero on the top surface, relativeSlopeTopOfEdge on the top ! surface is identified with its value on the second interface. relativeSlopeTopOfEdge(1,iEdge) = relativeSlopeTopOfEdge(2,iEdge) ! dDensityDzTopOfEdge may or may not equal zero on the bottom surface, depending on whether ! maxLevelEdgeTop(iEdge) = maxLevelEdgeBottom(iEdge). But here we ! take a simplistic approach and identify relativeSlopeTopOfEdge on the bottom surface with its value on ! the interface just above. relativeSlopeTopOfEdge( maxLevelEdgeTop(iEdge)+1, iEdge ) = relativeSlopeTopOfEdge( max(1,maxLevelEdgeTop(iEdge)), iEdge ) end do !$acc end parallel !$omp end do nEdges = nEdgesArray( 3 ) ! slope can be unbounded in regions of neutral stability, reset to the large, but bounded, value ! values is hardwrite to 1.0, this is equivalent to a slope of 45 degrees !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iEdge = 1, nEdges !$acc loop vector do k = 1, nVertLevels relativeSlopeTopOfEdge(k, iEdge) = max( min( relativeSlopeTopOfEdge(k, iEdge), 1.0_RKIND), -1.0_RKIND) end do end do !$acc end parallel !$omp end do ! average relative slope to cell centers ! do this by computing (relative slope)^2, then taking sqrt nCells = nCellsArray( 2 ) !$omp do schedule(runtime) private(i, iEdge, areaEdge, rtmp, k) !$acc parallel loop gang do iCell = 1, nCells areaCellSum(:, iCell) = 1.0e-34_RKIND do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) !contribution of cell area from this edge * 2.0 areaEdge = 0.5_RKIND * dcEdge(iEdge) * dvEdge(iEdge) !$acc loop vector do k = 1, maxLevelEdgeTop(iEdge) rtmp = areaEdge * relativeSlopeTopOfEdge(k, iEdge)**2 relativeSlopeTopOfCell(k, iCell) = relativeSlopeTopOfCell(k, iCell) + rtmp areaCellSum(k, iCell) = areaCellSum(k, iCell) + areaEdge end do end do end do !$acc end parallel !$omp end do nCells = nCellsArray( 2 ) !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iCell=1,nCells !$acc loop vector do k = 1, maxLevelCell(iCell) relativeSlopeTopOfCell(k,iCell) = sqrt( relativeSlopeTopOfCell(k,iCell)/areaCellSum(k,iCell) ) end do end do !$acc end parallel !$omp end do ! Compute tapering function ! Compute k33 at cell center and layer interface nCells = nCellsArray( size(nCellsArray) ) !$omp do schedule(runtime) do iCell = 1, nCells k33(:, iCell) = 0.0_RKIND end do !$omp end do ! use relativeSlopeTaperingCell as a temporary space for smoothing of relativeSlopeTopOfCell relativeSlopeTaperingCell = relativeSlopeTopOfCell do iter = 1, 5 nCells = nCellsArray( 2 ) !$omp do schedule(runtime) do iCell=1,nCells relativeSlopeTaperingCell(1, iCell) = 0.0_RKIND relativeSlopeTaperingCell(maxLevelCell(iCell):nVertLevels, iCell) = 0.0_RKIND do k = 2, maxLevelCell(iCell)-1 rtmp = relativeSlopeTopOfCell(k-1,iCell) + relativeSlopeTopOfCell(k+1,iCell) stmp = 2.0_RKIND*relativeSlopeTopOfCell(k,iCell) relativeSlopeTaperingCell(k,iCell) = (rtmp+stmp)/4.0_RKIND end do relativeSlopeTopOfCell(:, iCell) = relativeSlopeTaperingCell(:, iCell) end do !$omp end do end do ! iter nCells = nCellsArray ( 2 ) ! first, compute tapering across full domain based on a maximum allowable slope !$omp do schedule(runtime) private(k) do iCell=1,nCells do k = 1, maxLevelCell(iCell) relativeSlopeTaperingCell(k,iCell) = min(1.0_RKIND, config_max_relative_slope / (relativeSlopeTopOfCell(k,iCell)+epsGM)) end do end do !$omp end do ! now further taper in the boundary layer ! vertical (k33) tapering starts at 2*OBL, increases linearly to OBL and is held uniform across OBL ! rtmp = 1 @ zMid = -2.0*OBL, rtmp = 0 @ zMid = -OBL if(config_use_Redi_surface_layer_tapering) then nCells = nCellsArray ( 2 ) !$omp do schedule(runtime) private(k, rtmp) do iCell=1,nCells do k = 1, maxLevelCell(iCell) rtmp = -zMid(k,iCell)/max(config_Redi_surface_layer_tapering_extent,boundaryLayerDepth(iCell)+epsGM) rtmp = max(0.0_RKIND,rtmp) rtmp = min(1.0_RKIND,rtmp) relativeSlopeTaperingCell(k,iCell) = rtmp*relativeSlopeTaperingCell(k,iCell) end do end do !$omp end do endif ! config_use_Redi_surface_layer_tapering ! now further taper in the boundary layer ! vertical (k33) tapering starts at 2*OBL, increases linearly to OBL and is held uniform across OBL ! rtmp = 1 @ zMid = zMid(maxLevelCell) + config_Redi_bottom_layer_tapering_depth, rtmp = 0 @ zMid = zMid(maxLevelCell) if(config_use_Redi_bottom_layer_tapering) then nCells = nCellsArray ( 2 ) !$omp do schedule(runtime) private(k, rtmp) do iCell=1,nCells do k = 1, maxLevelCell(iCell) rtmp = (zMid(k,iCell)-zMid(maxLevelCell(iCell),iCell))/(config_Redi_bottom_layer_tapering_depth+epsGM) rtmp = max(0.0_RKIND,rtmp) rtmp = min(1.0_RKIND,rtmp) relativeSlopeTaperingCell(k,iCell) = rtmp*relativeSlopeTaperingCell(k,iCell) end do end do !$omp end do endif ! config_use_Redi_bottom_layer_tapering nCells = nCellsArray( 2 ) !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iCell=1,nCells k33(:, iCell) = 0.0_RKIND !$acc loop vector do k = 2, maxLevelCell(iCell) k33(k,iCell) = ( relativeSlopeTaperingCell(k,iCell) * relativeSlopeTopOfCell(k,iCell) )**2 end do end do !$acc end parallel !$omp end do nEdges = nEdgesArray( 3 ) ! average tapering function to layer edges !$omp do schedule(runtime) private(cell1, cell2, k) !$acc parallel loop gang do iEdge = 1, nEdges cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) !$acc loop vector do k = 1, maxLevelEdgeTop(iEdge) relativeSlopeTapering(k,iEdge) = 0.5_RKIND * (relativeSlopeTaperingCell(k,cell1) + relativeSlopeTaperingCell(k,cell2)) enddo enddo !$acc end parallel !$omp end do ! allow disabling of K33 for testing if(config_disable_redi_k33) then nCells = nCellsArray( size(nCellsArray) ) !$omp do schedule(runtime) !$acc parallel loop gang do iCell = 1, nCells k33(:, iCell) = 0.0_RKIND end do !$acc end parallel !$omp end do end if !-------------------------------------------------------------------- ! ! Compute stream function and Bolus velocity for Bolus part of GM ! !-------------------------------------------------------------------- if (config_gm_lat_variable_c2) then !$omp do schedule(runtime) private(cell1, cell2, sumN2, ltSum, countN2, BruntVaisalaFreqTopEdge) do iEdge = 1, nEdges cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) sumN2 = 0.0 ltSum = 0.0 countN2 = 0 do k=2,maxLevelEdgeTop(iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) sumN2 = sumN2 + BruntVaisalaFreqTopEdge*layerThicknessEdge(k,iEdge) ltSum = ltSum + layerThicknessEdge(k,iEdge) countN2 = countN2 +1 enddo if(countN2 > 0) cGMphaseSpeed(iEdge) = max(config_gm_min_phase_speed ,sqrt(sumN2/ltSum)*ltSum / 3.141592_RKIND) enddo !$omp end do else !$omp do schedule(runtime) do iEdge = 1, nEdges cGMphaseSpeed(iEdge) = config_gravWaveSpeed_trunc enddo !$omp end do endif !$omp do schedule(runtime) do iEdge=1,nEdges kappaGM3D(:,iEdge) = gmBolusKappa(iEdge) enddo !$omp end do if (config_gm_kappa_lat_depth_variable) then !$omp do schedule(runtime) private(cell1, cell2, k, BruntVaisalaFreqTopEdge, maxN) do iEdge = 1,nEdges cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) maxN = -1.0_RKIND do k=2,maxLevelEdgeTop(iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) maxN = max(maxN,BruntVaisalaFreqTopEdge) enddo do k=2,maxLevelEdgeTop(iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) kappaGM3D(k,iEdge) = gmBolusKappa(iEdge)*max(config_gm_min_stratification_ratio, & BruntVaisalaFreqTopEdge / (maxN + 1.0E-10_RKIND)) enddo enddo !$omp end do endif nEdges = nEdgesArray( 3 ) !$omp do schedule(runtime) do iEdge = 1, nEdges cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) gmStreamFuncTopOfEdge(:, iEdge) = 0.0_RKIND ! Construct the tridiagonal matrix if (maxLevelEdgeTop(iEdge) .GE. 3) then ! First row k = 2 BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1,iEdge) & * layerThicknessEdge(k,iEdge)) - BruntVaisalaFreqTopEdge tridiagC(k-1) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k, iEdge) & / (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge)) rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge) ! Second to next to the last rows do k = 3, maxLevelEdgeTop(iEdge)-1 BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) tridiagA(k-2) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k-1, iEdge) & / (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge)) tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1, iEdge) & * layerThicknessEdge(k, iEdge) ) - BruntVaisalaFreqTopEdge tridiagC(k-1) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k, iEdge) & / (layerThicknessEdge(k-1, iEdge) + layerThicknessEdge(k, iEdge)) rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge) end do ! Last row k = maxLevelEdgeTop(iEdge) BruntVaisalaFreqTopEdge = 0.5_RKIND * (BruntVaisalaFreqTop(k,cell1) + BruntVaisalaFreqTop(k,cell2)) BruntVaisalaFreqTopEdge = max(BruntVaisalaFreqTopEdge, 0.0_RKIND) tridiagA(k-2) = 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / layerThicknessEdge(k-1,iEdge) & / (layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge)) tridiagB(k-1) = - 2.0_RKIND * cGMphaseSpeed(iEdge)**2 / (layerThicknessEdge(k-1, iEdge) & * layerThicknessEdge(k, iEdge)) - BruntVaisalaFreqTopEdge rightHandSide(k-1) = kappaGM3D(k-1,iEdge) * gravity / rho_sw * gradDensityConstZTopOfEdge(k,iEdge) ! Total number of rows N = maxLevelEdgeTop(iEdge) - 1 ! Call the tridiagonal solver call tridiagonal_solve(tridiagA, tridiagB, tridiagC, rightHandSide, & gmStreamFuncTopOfEdge(2:maxLevelEdgeTop(iEdge), iEdge), N) end if end do !$omp end do nEdges = nEdgesArray( 3 ) ! Compute normalGMBolusVelocity from the stream function !$omp do schedule(runtime) private(k) !$acc parallel loop gang do iEdge = 1, nEdges !$acc loop vector do k = 1, maxLevelEdgeTop(iEdge) normalGMBolusVelocity(k,iEdge) = (gmStreamFuncTopOfEdge(k,iEdge) - gmStreamFuncTopOfEdge(k+1,iEdge)) & / layerThicknessEdge(k,iEdge) end do end do !$acc end parallel !$omp end do nCells = nCellsArray( 1 ) ! Interpolate gmStreamFuncTopOfEdge to cell centers for visualization !$omp do schedule(runtime) private(i, iEdge, areaEdge, k, rtmp) do iCell = 1, nCells gmStreamFuncTopOfCell(:, iCell) = 0.0_RKIND do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) areaEdge = 0.25_RKIND * dcEdge(iEdge) * dvEdge(iEdge) do k = 1, maxLevelEdgeTop(iEdge) rtmp = 0.5_RKIND * ( gmStreamFuncTopOfEdge(k, iEdge) + gmStreamFuncTopOfEdge(k+1, iEdge) ) * areaEdge gmStreamFuncTopOfCell(k, iCell) = gmStreamFuncTopOfCell(k, iCell) + rtmp end do end do end do !$omp end do !$omp do schedule(runtime) do iCell = 1, nCells gmStreamFuncTopOfCell(:, iCell) = gmStreamFuncTopOfCell(:,iCell) / areaCell(iCell) end do END DO CALL SYSTEM_CLOCK(kgen_stop_clock, kgen_rate_clock) kgen_measure = 1.0D6*(kgen_stop_clock - kgen_start_clock)/DBLE(kgen_rate_clock*KGEN_MAXITER) #ifdef _MPI CALL mpi_allreduce(kgen_measure, gkgen_measure, 1, mpi_real8, mpi_max, mpi_comm_world, kgen_ierr) kgen_measure = gkgen_measure #endif IF (check_status%rank==0) THEN WRITE (*, *) "ocn_gm_compute_Bolus_velocity : Time per call (usec): ", kgen_measure END IF END IF IF (kgen_warmupstage) THEN END IF IF (kgen_evalstage) THEN END IF !$omp end do !$acc end data !$kgen end_callsite ! Deallocate scratch variables CONTAINS !read state subroutine for kr_kgen_ocn_gm_compute_bolus_velocity_subp0 SUBROUTINE kr_kgen_ocn_gm_compute_bolus_velocity_subp0(var, kgen_unit, printname, printvar) REAL(KIND=rkind), INTENT(INOUT), POINTER, DIMENSION(:,:) :: var INTEGER, INTENT(IN) :: kgen_unit CHARACTER(LEN=*), INTENT(IN) :: printname LOGICAL, INTENT(IN), OPTIONAL :: printvar LOGICAL :: kgen_istrue REAL(KIND=8) :: kgen_array_sum INTEGER :: idx1, idx2 INTEGER, DIMENSION(2,2) :: kgen_bound READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( var )) THEN NULLIFY (var) END IF READ (UNIT = kgen_unit) kgen_array_sum READ (UNIT = kgen_unit) kgen_bound(1, 1) READ (UNIT = kgen_unit) kgen_bound(2, 1) READ (UNIT = kgen_unit) kgen_bound(1, 2) READ (UNIT = kgen_unit) kgen_bound(2, 2) ALLOCATE (var(kgen_bound(1,1):kgen_bound(2,1), kgen_bound(1,2):kgen_bound(2,2))) READ (UNIT = kgen_unit) var CALL kgen_array_sumcheck(printname, kgen_array_sum, DBLE(SUM(var, mask=(var .eq. var))), .TRUE.) IF (PRESENT( printvar ) .AND. printvar) THEN WRITE (*, *) "KGEN DEBUG: DBLE(SUM(" // printname // ")) = ", DBLE(SUM(var, mask=(var .eq. var))) END IF END IF END SUBROUTINE kr_kgen_ocn_gm_compute_bolus_velocity_subp0 !read state subroutine for kr_kgen_ocn_gm_compute_bolus_velocity_subp1 SUBROUTINE kr_kgen_ocn_gm_compute_bolus_velocity_subp1(var, kgen_unit, printname, printvar) REAL(KIND=rkind), INTENT(INOUT), POINTER, DIMENSION(:) :: var INTEGER, INTENT(IN) :: kgen_unit CHARACTER(LEN=*), INTENT(IN) :: printname LOGICAL, INTENT(IN), OPTIONAL :: printvar LOGICAL :: kgen_istrue REAL(KIND=8) :: kgen_array_sum INTEGER :: idx1 INTEGER, DIMENSION(2,1) :: kgen_bound READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( var )) THEN NULLIFY (var) END IF READ (UNIT = kgen_unit) kgen_array_sum READ (UNIT = kgen_unit) kgen_bound(1, 1) READ (UNIT = kgen_unit) kgen_bound(2, 1) ALLOCATE (var(kgen_bound(1,1):kgen_bound(2,1))) READ (UNIT = kgen_unit) var CALL kgen_array_sumcheck(printname, kgen_array_sum, DBLE(SUM(var, mask=(var .eq. var))), .TRUE.) IF (PRESENT( printvar ) .AND. printvar) THEN WRITE (*, *) "KGEN DEBUG: DBLE(SUM(" // printname // ")) = ", DBLE(SUM(var, mask=(var .eq. var))) END IF END IF END SUBROUTINE kr_kgen_ocn_gm_compute_bolus_velocity_subp1 !read state subroutine for kr_kgen_ocn_gm_compute_bolus_velocity_subp2 SUBROUTINE kr_kgen_ocn_gm_compute_bolus_velocity_subp2(var, kgen_unit, printname, printvar) INTEGER, INTENT(INOUT), POINTER, DIMENSION(:) :: var INTEGER, INTENT(IN) :: kgen_unit CHARACTER(LEN=*), INTENT(IN) :: printname LOGICAL, INTENT(IN), OPTIONAL :: printvar LOGICAL :: kgen_istrue REAL(KIND=8) :: kgen_array_sum INTEGER :: idx1 INTEGER, DIMENSION(2,1) :: kgen_bound READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( var )) THEN NULLIFY (var) END IF READ (UNIT = kgen_unit) kgen_array_sum READ (UNIT = kgen_unit) kgen_bound(1, 1) READ (UNIT = kgen_unit) kgen_bound(2, 1) ALLOCATE (var(kgen_bound(1,1):kgen_bound(2,1))) READ (UNIT = kgen_unit) var CALL kgen_array_sumcheck(printname, kgen_array_sum, DBLE(SUM(var, mask=(var .eq. var))), .TRUE.) IF (PRESENT( printvar ) .AND. printvar) THEN WRITE (*, *) "KGEN DEBUG: DBLE(SUM(" // printname // ")) = ", DBLE(SUM(var, mask=(var .eq. var))) END IF END IF END SUBROUTINE kr_kgen_ocn_gm_compute_bolus_velocity_subp2 !read state subroutine for kr_kgen_ocn_gm_compute_bolus_velocity_subp3 SUBROUTINE kr_kgen_ocn_gm_compute_bolus_velocity_subp3(var, kgen_unit, printname, printvar) INTEGER, INTENT(INOUT), POINTER, DIMENSION(:,:) :: var INTEGER, INTENT(IN) :: kgen_unit CHARACTER(LEN=*), INTENT(IN) :: printname LOGICAL, INTENT(IN), OPTIONAL :: printvar LOGICAL :: kgen_istrue REAL(KIND=8) :: kgen_array_sum INTEGER :: idx1, idx2 INTEGER, DIMENSION(2,2) :: kgen_bound READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( var )) THEN NULLIFY (var) END IF READ (UNIT = kgen_unit) kgen_array_sum READ (UNIT = kgen_unit) kgen_bound(1, 1) READ (UNIT = kgen_unit) kgen_bound(2, 1) READ (UNIT = kgen_unit) kgen_bound(1, 2) READ (UNIT = kgen_unit) kgen_bound(2, 2) ALLOCATE (var(kgen_bound(1,1):kgen_bound(2,1), kgen_bound(1,2):kgen_bound(2,2))) READ (UNIT = kgen_unit) var CALL kgen_array_sumcheck(printname, kgen_array_sum, DBLE(SUM(var, mask=(var .eq. var))), .TRUE.) IF (PRESENT( printvar ) .AND. printvar) THEN WRITE (*, *) "KGEN DEBUG: DBLE(SUM(" // printname // ")) = ", DBLE(SUM(var, mask=(var .eq. var))) END IF END IF END SUBROUTINE kr_kgen_ocn_gm_compute_bolus_velocity_subp3 !read state subroutine for kr_kgen_ocn_gm_compute_bolus_velocity_subp5 SUBROUTINE kr_kgen_ocn_gm_compute_bolus_velocity_subp5(var, kgen_unit, printname, printvar) INTEGER, INTENT(INOUT), POINTER :: var INTEGER, INTENT(IN) :: kgen_unit CHARACTER(LEN=*), INTENT(IN) :: printname LOGICAL, INTENT(IN), OPTIONAL :: printvar LOGICAL :: kgen_istrue REAL(KIND=8) :: kgen_array_sum READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( var )) THEN NULLIFY (var) END IF ALLOCATE (var) READ (UNIT = kgen_unit) var IF (PRESENT( printvar ) .AND. printvar) THEN WRITE (*, *) "KGEN DEBUG: " // printname // " = ", var END IF END IF END SUBROUTINE kr_kgen_ocn_gm_compute_bolus_velocity_subp5 !verify state subroutine for kv_kgen_ocn_gm_compute_bolus_velocity_subp0 RECURSIVE SUBROUTINE kv_kgen_ocn_gm_compute_bolus_velocity_subp0(varname, check_status, var, kgenref_var) CHARACTER(LEN=*), INTENT(IN) :: varname TYPE(check_t), INTENT(INOUT) :: check_status REAL(KIND=rkind), pointer, INTENT(IN), DIMENSION(:,:) :: var, kgenref_var INTEGER :: check_result LOGICAL :: is_print = .FALSE. INTEGER :: idx1, idx2 INTEGER :: n real(KIND=rkind) :: nrmsdiff, rmsdiff real(KIND=rkind), ALLOCATABLE :: buf1(:,:), buf2(:,:) IF (ASSOCIATED(var)) THEN check_status%numTotal = check_status%numTotal + 1 IF (ALL(var == kgenref_var)) THEN check_status%numIdentical = check_status%numIdentical + 1 IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is IDENTICAL." END IF END IF check_result = CHECK_IDENTICAL ELSE ALLOCATE (buf1(SIZE(var,dim=1),SIZE(var,dim=2))) ALLOCATE (buf2(SIZE(var,dim=1),SIZE(var,dim=2))) n = SIZE(var) WHERE ( ABS(kgenref_var) > kgen_minvalue ) buf1 = ((var-kgenref_var)/kgenref_var)**2 buf2 = (var-kgenref_var)**2 ELSEWHERE buf1 = (var-kgenref_var)**2 buf2 = buf1 END WHERE nrmsdiff = SQRT(SUM(buf1)/DBLE(n)) rmsdiff = SQRT(SUM(buf2)/DBLE(n)) IF (rmsdiff > kgen_tolerance) THEN check_status%numOutTol = check_status%numOutTol + 1 IF (kgen_verboseLevel > 0) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is NOT IDENTICAL(out of tolerance)." END IF END IF check_result = CHECK_OUT_TOL ELSE check_status%numInTol = check_status%numInTol + 1 IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is NOT IDENTICAL(within tolerance)." END IF END IF check_result = CHECK_IN_TOL END IF END IF IF (check_result == CHECK_IDENTICAL) THEN IF (kgen_verboseLevel > 2) THEN IF (check_status%rank == 0) THEN WRITE (*, *) count( var /= kgenref_var), " of ", size( var ), " elements are different." WRITE (*, *) "Average - kernel ", sum(var)/real(size(var)) WRITE (*, *) "Average - reference ", sum(kgenref_var)/real(size(kgenref_var)) WRITE (*, *) "RMS of difference is ", 0 WRITE (*, *) "Normalized RMS of difference is ", 0 WRITE (*, *) "" END IF END IF ELSE IF (check_result == CHECK_OUT_TOL) THEN IF (kgen_verboseLevel > 0) THEN IF (check_status%rank == 0) THEN WRITE (*, *) count( var /= kgenref_var), " of ", size( var ), " elements are different." WRITE (*, *) "Average - kernel ", sum(var)/real(size(var)) WRITE (*, *) "Average - reference ", sum(kgenref_var)/real(size(kgenref_var)) WRITE (*, *) "RMS of difference is ", rmsdiff WRITE (*, *) "Normalized RMS of difference is ", nrmsdiff WRITE (*, *) "" END IF END IF ELSE IF (check_result == CHECK_IN_TOL) THEN IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) count( var /= kgenref_var), " of ", size( var ), " elements are different." WRITE (*, *) "Average - kernel ", sum(var)/real(size(var)) WRITE (*, *) "Average - reference ", sum(kgenref_var)/real(size(kgenref_var)) WRITE (*, *) "RMS of difference is ", rmsdiff WRITE (*, *) "Normalized RMS of difference is ", nrmsdiff WRITE (*, *) "" END IF END IF END IF END IF END SUBROUTINE kv_kgen_ocn_gm_compute_bolus_velocity_subp0 !verify state subroutine for kv_kgen_ocn_gm_compute_bolus_velocity_subp1 RECURSIVE SUBROUTINE kv_kgen_ocn_gm_compute_bolus_velocity_subp1(varname, check_status, var, kgenref_var) CHARACTER(LEN=*), INTENT(IN) :: varname TYPE(check_t), INTENT(INOUT) :: check_status REAL(KIND=rkind), pointer, INTENT(IN), DIMENSION(:) :: var, kgenref_var INTEGER :: check_result LOGICAL :: is_print = .FALSE. INTEGER :: idx1 INTEGER :: n real(KIND=rkind) :: nrmsdiff, rmsdiff real(KIND=rkind), ALLOCATABLE :: buf1(:), buf2(:) IF (ASSOCIATED(var)) THEN check_status%numTotal = check_status%numTotal + 1 IF (ALL(var == kgenref_var)) THEN check_status%numIdentical = check_status%numIdentical + 1 IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is IDENTICAL." END IF END IF check_result = CHECK_IDENTICAL ELSE ALLOCATE (buf1(SIZE(var,dim=1))) ALLOCATE (buf2(SIZE(var,dim=1))) n = SIZE(var) WHERE ( ABS(kgenref_var) > kgen_minvalue ) buf1 = ((var-kgenref_var)/kgenref_var)**2 buf2 = (var-kgenref_var)**2 ELSEWHERE buf1 = (var-kgenref_var)**2 buf2 = buf1 END WHERE nrmsdiff = SQRT(SUM(buf1)/DBLE(n)) rmsdiff = SQRT(SUM(buf2)/DBLE(n)) IF (rmsdiff > kgen_tolerance) THEN check_status%numOutTol = check_status%numOutTol + 1 IF (kgen_verboseLevel > 0) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is NOT IDENTICAL(out of tolerance)." END IF END IF check_result = CHECK_OUT_TOL ELSE check_status%numInTol = check_status%numInTol + 1 IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is NOT IDENTICAL(within tolerance)." END IF END IF check_result = CHECK_IN_TOL END IF END IF IF (check_result == CHECK_IDENTICAL) THEN IF (kgen_verboseLevel > 2) THEN IF (check_status%rank == 0) THEN WRITE (*, *) count( var /= kgenref_var), " of ", size( var ), " elements are different." WRITE (*, *) "Average - kernel ", sum(var)/real(size(var)) WRITE (*, *) "Average - reference ", sum(kgenref_var)/real(size(kgenref_var)) WRITE (*, *) "RMS of difference is ", 0 WRITE (*, *) "Normalized RMS of difference is ", 0 WRITE (*, *) "" END IF END IF ELSE IF (check_result == CHECK_OUT_TOL) THEN IF (kgen_verboseLevel > 0) THEN IF (check_status%rank == 0) THEN WRITE (*, *) count( var /= kgenref_var), " of ", size( var ), " elements are different." WRITE (*, *) "Average - kernel ", sum(var)/real(size(var)) WRITE (*, *) "Average - reference ", sum(kgenref_var)/real(size(kgenref_var)) WRITE (*, *) "RMS of difference is ", rmsdiff WRITE (*, *) "Normalized RMS of difference is ", nrmsdiff WRITE (*, *) "" END IF END IF ELSE IF (check_result == CHECK_IN_TOL) THEN IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) count( var /= kgenref_var), " of ", size( var ), " elements are different." WRITE (*, *) "Average - kernel ", sum(var)/real(size(var)) WRITE (*, *) "Average - reference ", sum(kgenref_var)/real(size(kgenref_var)) WRITE (*, *) "RMS of difference is ", rmsdiff WRITE (*, *) "Normalized RMS of difference is ", nrmsdiff WRITE (*, *) "" END IF END IF END IF END IF END SUBROUTINE kv_kgen_ocn_gm_compute_bolus_velocity_subp1 !verify state subroutine for kv_kgen_ocn_gm_compute_bolus_velocity_subp2 RECURSIVE SUBROUTINE kv_kgen_ocn_gm_compute_bolus_velocity_subp2(varname, check_status, var, kgenref_var) CHARACTER(LEN=*), INTENT(IN) :: varname TYPE(check_t), INTENT(INOUT) :: check_status INTEGER, pointer, INTENT(IN), DIMENSION(:) :: var, kgenref_var INTEGER :: check_result LOGICAL :: is_print = .FALSE. INTEGER :: idx1 INTEGER :: n integer :: nrmsdiff, rmsdiff integer, ALLOCATABLE :: buf1(:), buf2(:) IF (ASSOCIATED(var)) THEN check_status%numTotal = check_status%numTotal + 1 IF (ALL(var == kgenref_var)) THEN check_status%numIdentical = check_status%numIdentical + 1 IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is IDENTICAL." END IF END IF check_result = CHECK_IDENTICAL ELSE ALLOCATE (buf1(SIZE(var,dim=1))) ALLOCATE (buf2(SIZE(var,dim=1))) n = SIZE(var) WHERE ( ABS(kgenref_var) > kgen_minvalue ) buf1 = ((var-kgenref_var)/kgenref_var)**2 buf2 = (var-kgenref_var)**2 ELSEWHERE buf1 = (var-kgenref_var)**2 buf2 = buf1 END WHERE nrmsdiff = SQRT(SUM(buf1)/DBLE(n)) rmsdiff = SQRT(SUM(buf2)/DBLE(n)) IF (rmsdiff > kgen_tolerance) THEN check_status%numOutTol = check_status%numOutTol + 1 IF (kgen_verboseLevel > 0) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is NOT IDENTICAL(out of tolerance)." END IF END IF check_result = CHECK_OUT_TOL ELSE check_status%numInTol = check_status%numInTol + 1 IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is NOT IDENTICAL(within tolerance)." END IF END IF check_result = CHECK_IN_TOL END IF END IF IF (check_result == CHECK_IDENTICAL) THEN IF (kgen_verboseLevel > 2) THEN IF (check_status%rank == 0) THEN WRITE (*, *) count( var /= kgenref_var), " of ", size( var ), " elements are different." WRITE (*, *) "Average - kernel ", sum(var)/real(size(var)) WRITE (*, *) "Average - reference ", sum(kgenref_var)/real(size(kgenref_var)) WRITE (*, *) "RMS of difference is ", 0 WRITE (*, *) "Normalized RMS of difference is ", 0 WRITE (*, *) "" END IF END IF ELSE IF (check_result == CHECK_OUT_TOL) THEN IF (kgen_verboseLevel > 0) THEN IF (check_status%rank == 0) THEN WRITE (*, *) count( var /= kgenref_var), " of ", size( var ), " elements are different." WRITE (*, *) "Average - kernel ", sum(var)/real(size(var)) WRITE (*, *) "Average - reference ", sum(kgenref_var)/real(size(kgenref_var)) WRITE (*, *) "RMS of difference is ", rmsdiff WRITE (*, *) "Normalized RMS of difference is ", nrmsdiff WRITE (*, *) "" END IF END IF ELSE IF (check_result == CHECK_IN_TOL) THEN IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) count( var /= kgenref_var), " of ", size( var ), " elements are different." WRITE (*, *) "Average - kernel ", sum(var)/real(size(var)) WRITE (*, *) "Average - reference ", sum(kgenref_var)/real(size(kgenref_var)) WRITE (*, *) "RMS of difference is ", rmsdiff WRITE (*, *) "Normalized RMS of difference is ", nrmsdiff WRITE (*, *) "" END IF END IF END IF END IF END SUBROUTINE kv_kgen_ocn_gm_compute_bolus_velocity_subp2 !verify state subroutine for kv_ocn_gm_compute_bolus_velocity_integer__ RECURSIVE SUBROUTINE kv_ocn_gm_compute_bolus_velocity_integer__(varname, check_status, var, kgenref_var) CHARACTER(LEN=*), INTENT(IN) :: varname TYPE(check_t), INTENT(INOUT) :: check_status INTEGER, INTENT(IN) :: var, kgenref_var INTEGER :: check_result LOGICAL :: is_print = .FALSE. integer :: diff check_status%numTotal = check_status%numTotal + 1 IF (var == kgenref_var) THEN check_status%numIdentical = check_status%numIdentical + 1 IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is IDENTICAL." END IF END IF check_result = CHECK_IDENTICAL ELSE diff = ABS(var - kgenref_var) IF (diff <= kgen_tolerance) THEN check_status%numInTol = check_status%numInTol + 1 IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is NOT IDENTICAL(within tolerance)." END IF END IF check_result = CHECK_IN_TOL ELSE check_status%numOutTol = check_status%numOutTol + 1 IF (kgen_verboseLevel > 0) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is NOT IDENTICAL(out of tolerance)." END IF END IF check_result = CHECK_OUT_TOL END IF END IF IF (check_result == CHECK_IDENTICAL) THEN IF (kgen_verboseLevel > 2) THEN IF (check_status%rank == 0) THEN WRITE (*, *) "Difference is ", 0 WRITE (*, *) "" END IF END IF ELSE IF (check_result == CHECK_OUT_TOL) THEN IF (kgen_verboseLevel > 0) THEN IF (check_status%rank == 0) THEN WRITE (*, *) "Difference is ", diff WRITE (*, *) "" END IF END IF ELSE IF (check_result == CHECK_IN_TOL) THEN IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) "Difference is ", diff WRITE (*, *) "" END IF END IF END IF END SUBROUTINE kv_ocn_gm_compute_bolus_velocity_integer__ !verify state subroutine for kv_kgen_ocn_gm_compute_bolus_velocity_subp3 RECURSIVE SUBROUTINE kv_kgen_ocn_gm_compute_bolus_velocity_subp3(varname, check_status, var, kgenref_var) CHARACTER(LEN=*), INTENT(IN) :: varname TYPE(check_t), INTENT(INOUT) :: check_status REAL(KIND=rkind), INTENT(IN) :: var, kgenref_var INTEGER :: check_result LOGICAL :: is_print = .FALSE. real(KIND=rkind) :: diff check_status%numTotal = check_status%numTotal + 1 IF (var == kgenref_var) THEN check_status%numIdentical = check_status%numIdentical + 1 IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is IDENTICAL." END IF END IF check_result = CHECK_IDENTICAL ELSE diff = ABS(var - kgenref_var) IF (diff <= kgen_tolerance) THEN check_status%numInTol = check_status%numInTol + 1 IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is NOT IDENTICAL(within tolerance)." END IF END IF check_result = CHECK_IN_TOL ELSE check_status%numOutTol = check_status%numOutTol + 1 IF (kgen_verboseLevel > 0) THEN IF (check_status%rank == 0) THEN WRITE (*, *) trim(adjustl(varname)), " is NOT IDENTICAL(out of tolerance)." END IF END IF check_result = CHECK_OUT_TOL END IF END IF IF (check_result == CHECK_IDENTICAL) THEN IF (kgen_verboseLevel > 2) THEN IF (check_status%rank == 0) THEN WRITE (*, *) "Difference is ", 0 WRITE (*, *) "" END IF END IF ELSE IF (check_result == CHECK_OUT_TOL) THEN IF (kgen_verboseLevel > 0) THEN IF (check_status%rank == 0) THEN WRITE (*, *) "Difference is ", diff WRITE (*, *) "" END IF END IF ELSE IF (check_result == CHECK_IN_TOL) THEN IF (kgen_verboseLevel > 1) THEN IF (check_status%rank == 0) THEN WRITE (*, *) "Difference is ", diff WRITE (*, *) "" END IF END IF END IF END SUBROUTINE kv_kgen_ocn_gm_compute_bolus_velocity_subp3 END SUBROUTINE ocn_gm_compute_bolus_velocity !*********************************************************************** ! routine tridiagonal_solve !> \brief Solve the matrix equation Ax=r for x, where A is tridiagonal. !> \author Mark Petersen !> \date September 2011 !> \details !> Solve the matrix equation Ax=r for x, where A is tridiagonal. !> A is an nxn matrix, with: !> a sub-diagonal, filled from 1:n-1 (a(1) appears on row 2) !> b diagonal, filled from 1:n !> c sup-diagonal, filled from 1:n-1 (c(1) apears on row 1) !----------------------------------------------------------------------- ! mrp note: This subroutine also appears in vmix and should really be put in the framework. ! ! ! subroutine tridiagonal_solve(a,b,c,r,x,n) !{{{ !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- ! ! integer,intent(in) :: n real (KIND=RKIND), dimension(n), intent(in) :: a,b,c,r !----------------------------------------------------------------- ! output variables !----------------------------------------------------------------- ! ! real (KIND=RKIND), dimension(n), intent(out) :: x !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- ! ! real (KIND=RKIND), dimension(n) :: bTemp,rTemp real (KIND=RKIND) :: m integer i ! Use work variables for b and r bTemp(1) = b(1) rTemp(1) = r(1) ! First pass: set the coefficients do i = 2,n m = a(i-1)/bTemp(i-1) bTemp(i) = b(i) - m*c(i-1) rTemp(i) = r(i) - m*rTemp(i-1) end do x(n) = rTemp(n)/bTemp(n) ! Second pass: back-substition do i = n-1, 1, -1 x(i) = (rTemp(i) - c(i)*x(i+1))/bTemp(i) end do end subroutine tridiagonal_solve !}}} !*********************************************************************** ! routine ocn_gm_init !> \brief Initializes ocean momentum horizontal pressure gradient !> \author Mark Petersen !> \date September 2011 !> \details !> This routine initializes parameters required for the computation of the !> horizontal pressure gradient. !----------------------------------------------------------------------- ! ! ! !*********************************************************************** !read state subroutine for kr_externs_in_ocn_gm SUBROUTINE kr_externs_in_ocn_gm(kgen_unit) INTEGER, INTENT(IN) :: kgen_unit LOGICAL :: kgen_istrue REAL(KIND=8) :: kgen_array_sum READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( config_gravwavespeed_trunc )) THEN NULLIFY (config_gravwavespeed_trunc) END IF ALLOCATE (config_gravwavespeed_trunc) READ (UNIT = kgen_unit) config_gravwavespeed_trunc END IF READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( config_max_relative_slope )) THEN NULLIFY (config_max_relative_slope) END IF ALLOCATE (config_max_relative_slope) READ (UNIT = kgen_unit) config_max_relative_slope END IF READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( config_disable_redi_k33 )) THEN NULLIFY (config_disable_redi_k33) END IF ALLOCATE (config_disable_redi_k33) READ (UNIT = kgen_unit) config_disable_redi_k33 END IF READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( config_use_redi_surface_layer_tapering )) THEN NULLIFY (config_use_redi_surface_layer_tapering) END IF ALLOCATE (config_use_redi_surface_layer_tapering) READ (UNIT = kgen_unit) config_use_redi_surface_layer_tapering END IF READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( config_use_redi_bottom_layer_tapering )) THEN NULLIFY (config_use_redi_bottom_layer_tapering) END IF ALLOCATE (config_use_redi_bottom_layer_tapering) READ (UNIT = kgen_unit) config_use_redi_bottom_layer_tapering END IF READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( config_redi_surface_layer_tapering_extent )) THEN NULLIFY (config_redi_surface_layer_tapering_extent) END IF ALLOCATE (config_redi_surface_layer_tapering_extent) READ (UNIT = kgen_unit) config_redi_surface_layer_tapering_extent END IF READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( config_redi_bottom_layer_tapering_depth )) THEN NULLIFY (config_redi_bottom_layer_tapering_depth) END IF ALLOCATE (config_redi_bottom_layer_tapering_depth) READ (UNIT = kgen_unit) config_redi_bottom_layer_tapering_depth END IF READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( config_gm_lat_variable_c2 )) THEN NULLIFY (config_gm_lat_variable_c2) END IF ALLOCATE (config_gm_lat_variable_c2) READ (UNIT = kgen_unit) config_gm_lat_variable_c2 END IF READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( config_gm_kappa_lat_depth_variable )) THEN NULLIFY (config_gm_kappa_lat_depth_variable) END IF ALLOCATE (config_gm_kappa_lat_depth_variable) READ (UNIT = kgen_unit) config_gm_kappa_lat_depth_variable END IF READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( config_gm_min_stratification_ratio )) THEN NULLIFY (config_gm_min_stratification_ratio) END IF ALLOCATE (config_gm_min_stratification_ratio) READ (UNIT = kgen_unit) config_gm_min_stratification_ratio END IF READ (UNIT = kgen_unit) kgen_istrue IF (kgen_istrue) THEN IF (ASSOCIATED( config_gm_min_phase_speed )) THEN NULLIFY (config_gm_min_phase_speed) END IF ALLOCATE (config_gm_min_phase_speed) READ (UNIT = kgen_unit) config_gm_min_phase_speed END IF END SUBROUTINE kr_externs_in_ocn_gm !read state subroutine for kr_externs_out_ocn_gm SUBROUTINE kr_externs_out_ocn_gm(kgen_unit) INTEGER, INTENT(IN) :: kgen_unit LOGICAL :: kgen_istrue REAL(KIND=8) :: kgen_array_sum END SUBROUTINE kr_externs_out_ocn_gm end module ocn_gm !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||