From 77105d25f012522ad891d1f93b2f96eb3ab5669e Mon Sep 17 00:00:00 2001 From: Phil Jones Date: Wed, 16 Dec 2020 15:26:40 -0500 Subject: [PATCH 1/5] GPU implementation of ocean baroclinic velocity tendencies - phase 1 - added OpenACC loop and data directives - removed unnecessary pool and pointer accesses - did a lot of general cleanup - retained current loop structure so CPU is bit for bit - will optimize the loops for CPU and GPU in a separate branch --- .../mpas_ocn_time_integration_rk4.F | 6 +- .../mpas_ocn_time_integration_split.F | 3 +- .../shared/mpas_ocn_surface_bulk_forcing.F | 141 ++-- .../shared/mpas_ocn_surface_land_ice_fluxes.F | 189 +++-- src/core_ocean/shared/mpas_ocn_tendency.F | 385 ++++++--- src/core_ocean/shared/mpas_ocn_vel_forcing.F | 91 +- ...pas_ocn_vel_forcing_explicit_bottom_drag.F | 115 +-- .../mpas_ocn_vel_forcing_surface_stress.F | 157 ++-- .../shared/mpas_ocn_vel_hadv_coriolis.F | 247 ++++-- src/core_ocean/shared/mpas_ocn_vel_hmix.F | 100 +-- .../shared/mpas_ocn_vel_hmix_del2.F | 160 ++-- .../shared/mpas_ocn_vel_hmix_del4.F | 288 ++++--- .../shared/mpas_ocn_vel_hmix_leith.F | 180 ++-- .../shared/mpas_ocn_vel_pressure_grad.F | 798 ++++++++---------- .../shared/mpas_ocn_vel_tidal_potential.F | 538 ++++++------ src/core_ocean/shared/mpas_ocn_vel_vadv.F | 168 ++-- 16 files changed, 1956 insertions(+), 1610 deletions(-) diff --git a/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F b/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F index 070244ae77..e24ffef1c5 100644 --- a/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/src/core_ocean/mode_forward/mpas_ocn_time_integration_rk4.F @@ -858,7 +858,8 @@ subroutine ocn_time_integrator_rk4_compute_vel_tends(block, dt, & vertAleTransportTop, err) endif - call ocn_tend_vel(tendPool, provisStatePool, forcingPool, diagnosticsPool, meshPool, scratchPool, 1, dt) + call ocn_tend_vel(tendPool, provisStatePool, forcingPool, & + diagnosticsPool, 1, dt) end subroutine ocn_time_integrator_rk4_compute_vel_tends!}}} @@ -1052,7 +1053,8 @@ subroutine ocn_time_integrator_rk4_compute_tends(block, dt, rkWeight, err)!{{{ vertAleTransportTop, err) endif - call ocn_tend_vel(tendPool, provisStatePool, forcingPool, diagnosticsPool, meshPool, scratchPool, 1, dt) + call ocn_tend_vel(tendPool, provisStatePool, forcingPool, & + diagnosticsPool, 1, dt) if (associated(highFreqThicknessProvis)) then call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, scratchPool, & diff --git a/src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F b/src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F index 01e27c2dfe..7b9fe57bb6 100644 --- a/src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F +++ b/src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F @@ -534,7 +534,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ sshCur, dt, vertAleTransportTop, err) endif - call ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshPool, scratchPool, stage1_tend_time, dt) + call ocn_tend_vel(tendPool, statePool, forcingPool, & + diagnosticsPool, stage1_tend_time, dt) block => block % next end do diff --git a/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F b/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F index bc3ba013a3..cbf268a19e 100644 --- a/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F +++ b/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F @@ -27,6 +27,7 @@ module ocn_surface_bulk_forcing use mpas_timekeeping use ocn_constants use ocn_config + use ocn_mesh use ocn_equation_of_state implicit none @@ -56,7 +57,9 @@ module ocn_surface_bulk_forcing ! !-------------------------------------------------------------------- - logical :: bulkWindStressOn, bulkThicknessFluxOn + logical :: & + bulkWindStressOff, &! on/off switch for bulk wind stress + bulkThicknessFluxOff ! on/off switch for bulk thickness flux !*********************************************************************** @@ -133,95 +136,121 @@ end subroutine ocn_surface_bulk_forcing_tracers!}}} !> \author Doug Jacobsen !> \date 04/25/12 !> \details -!> This routine computes the velocity forcing arrays used later in MPAS. +!> This routine computes the surface stress for velocity forcing +!> based on bulk flux formulae. ! !----------------------------------------------------------------------- - subroutine ocn_surface_bulk_forcing_vel(meshPool, forcingPool, surfaceStress, surfaceStressMagnitude, err)!{{{ + subroutine ocn_surface_bulk_forcing_vel( & + windStressZonal, windStressMeridional, & + sfcStress, sfcStressMag, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- - type (mpas_pool_type), intent(in) :: meshPool !< Input: mesh information - type (mpas_pool_type), intent(in) :: forcingPool !< Input: Forcing information + + real (kind=RKIND), dimension(:), intent(in) :: & + windStressZonal, &!< [in] zonal wind stress cell center + windStressMeridional !< [in] meridional wind stress cell center !----------------------------------------------------------------- - ! ! input/output variables - ! !----------------------------------------------------------------- - real (kind=RKIND), dimension(:), intent(inout) :: surfaceStress, & !< Input/Output: Array for surface stress - surfaceStressMagnitude !< Input/Output: Array for magnitude of surface stress + + real (kind=RKIND), dimension(:), intent(inout) :: & + sfcStress, & !< [inout] Array for surface stress + sfcStressMag !< [inout] Array for magnitude of surface stress !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: Error flag + integer, intent(out) :: err !< [out] Error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer :: iEdge, cell1, cell2, iCell, nCells, nEdges - integer, dimension(:), pointer :: nCellsArray, nEdgesArray + integer :: & + iEdge, iCell, &! loop indices for edge, cell loops + nCells, nEdges, &! number of cells, edges participating + cell1, cell2 ! neighbor cell indices across edge - integer, dimension(:,:), pointer :: cellsOnEdge + real (kind=RKIND) :: & + meridWSEdge, &! meridional wind stress avg to edge + zonalWSEdge ! zonal wind stress averaged to edge - real (kind=RKIND) :: meridionalAverage, zonalAverage - real (kind=RKIND), dimension(:), pointer :: angleEdge - real (kind=RKIND), dimension(:), pointer :: windStressZonal, windStressMeridional + ! End preamble + !----------------------------------------------------------------- + ! Begin code - err = 0 + !*** Set error code and exit early if bulk wind stress not on + !*** Otherwise start timer - if ( .not. bulkWindStressOn ) return + err = 0 + if (bulkWindStressOff) return call mpas_timer_start("bulk_ws", .false.) - call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) - call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) + !*** Forcing needed over halo regions - call mpas_pool_get_array(meshPool, 'angleEdge', angleEdge) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + nEdges = nEdgesHalo(3) + nCells = nCellsHalo(2) - call mpas_pool_get_array(forcingPool, 'windStressZonal', windStressZonal) - call mpas_pool_get_array(forcingPool, 'windStressMeridional', windStressMeridional) + !*** Convert coupled climate model wind stress to internal + !*** ocean model wind stress by averaging cell-centered wind + !*** stress to edges and then rotating vectors - nEdges = nEdgesArray( 4 ) - nCells = nCellsArray( 3 ) - - ! Convert coupled climate model wind stress to MPAS-O wind stress + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, angleEdge, sfcStress, & + !$acc windStressZonal, windStressMeridional) & + !$acc private(cell1, cell2, zonalWSEdge, meridWSEdge) + #else !$omp parallel - !$omp do schedule(runtime) private(cell1, cell2, zonalAverage, meridionalAverage) + !$omp do schedule(runtime) & + !$omp private(cell1, cell2, zonalWSEdge, meridWSEdge) + #endif do iEdge = 1, nEdges cell1 = cellsOnEdge(1, iEdge) cell2 = cellsOnEdge(2, iEdge) - zonalAverage = 0.5_RKIND * (windStressZonal(cell1) + windStressZonal(cell2)) - meridionalAverage = 0.5_RKIND * (windStressMeridional(cell1) + windStressMeridional(cell2)) + zonalWSEdge = 0.5_RKIND * (windStressZonal(cell1) + & + windStressZonal(cell2)) + meridWSEdge = 0.5_RKIND * (windStressMeridional(cell1) + & + windStressMeridional(cell2)) - surfaceStress(iEdge) = surfaceStress(iEdge) + cos(angleEdge(iEdge)) * zonalAverage + sin(angleEdge(iEdge)) & - * meridionalAverage + sfcStress(iEdge) = sfcStress(iEdge) + & + cos(angleEdge(iEdge))*zonalWSEdge + & + sin(angleEdge(iEdge))*meridWSEdge end do + #ifndef MPAS_OPENACC !$omp end do + #endif - ! Build surface fluxes at cell centers + ! Calculate surface flux magnitude at cell centers + + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(sfcStressMag, windStressZonal, & + !$acc windStressMeridional) + #else !$omp do schedule(runtime) + #endif do iCell = 1, nCells - surfaceStressMagnitude(iCell) = surfaceStressMagnitude(iCell) + sqrt( windStressZonal(iCell)**2 & - + windStressMeridional(iCell)**2 ) + sfcStressMag(iCell) = sfcStressMag(iCell) + & + sqrt( windStressZonal(iCell)**2 & + + windStressMeridional(iCell)**2 ) end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif call mpas_timer_stop("bulk_ws") + !----------------------------------------------------------------- + end subroutine ocn_surface_bulk_forcing_vel!}}} !*********************************************************************** @@ -280,7 +309,7 @@ subroutine ocn_surface_bulk_forcing_thick(meshPool, forcingPool, surfaceThicknes err = 0 - if ( .not. bulkThicknessFluxOn ) return + if (bulkThicknessFluxOff) return call mpas_timer_start("bulk_thick", .false.) @@ -325,12 +354,30 @@ end subroutine ocn_surface_bulk_forcing_thick!}}} subroutine ocn_surface_bulk_forcing_init(err)!{{{ - integer, intent(out) :: err !< Output: error flag + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< [out] error flag + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + !*** Initialize error flag and module defaults err = 0 + bulkWindStressOff = .true. + bulkThicknessFluxOff = .true. + + !*** Set on/off flags based in user input configuration - bulkWindStressOn = config_use_bulk_wind_stress - bulkThicknessFluxOn = config_use_bulk_thickness_flux + if (config_use_bulk_wind_stress) & + bulkWindStressOff = .false. + if (config_use_bulk_thickness_flux) & + bulkThicknessFluxOff = .false. + + !----------------------------------------------------------------- end subroutine ocn_surface_bulk_forcing_init!}}} @@ -470,7 +517,7 @@ subroutine ocn_surface_bulk_forcing_active_tracers(meshPool, forcingPool, tracer ! Assume surface fluxes of water have zero salinity. So the RHS forcing is zero for salinity. ! Only include this heat forcing when bulk thickness is turned on ! indices on tracerGroup are (iTracer, iLevel, iCell) - if (bulkThicknessFluxOn) then + if (.not. bulkThicknessFluxOff) then !$omp parallel !$omp do schedule(runtime) do iCell = 1, nCells diff --git a/src/core_ocean/shared/mpas_ocn_surface_land_ice_fluxes.F b/src/core_ocean/shared/mpas_ocn_surface_land_ice_fluxes.F index a74a5e34dd..069231789f 100644 --- a/src/core_ocean/shared/mpas_ocn_surface_land_ice_fluxes.F +++ b/src/core_ocean/shared/mpas_ocn_surface_land_ice_fluxes.F @@ -27,6 +27,7 @@ module ocn_surface_land_ice_fluxes use ocn_constants use ocn_config + use ocn_mesh use ocn_equation_of_state implicit none @@ -57,9 +58,18 @@ module ocn_surface_land_ice_fluxes ! !-------------------------------------------------------------------- - logical :: landIceFluxesOn, standaloneOn, isomipOn, jenkinsOn, hollandJenkinsOn + logical :: & + landIceFluxesOff, &! on/off switch for land ice fluxes + standaloneOff, &! flag to determine coupled or standalone + isomipOn, &! use ISOMIP flux formulation + jenkinsOn, &! use Jenkins flux formulation + hollandJenkinsOn, &! use Holland-Jenkins flux formulation + useHollandJenkinsAdvDiff ! use Holland-Jenkins advection-diffusion - real (kind=RKIND) :: cp_land_ice, rho_land_ice + real (kind=RKIND) :: & + cpLandIce, &! specific heat for land ice + rhoLandIce, &! density of land ice + ISOMIPgammaT ! parameter for ISOMIP formulation !*********************************************************************** @@ -112,7 +122,7 @@ subroutine ocn_surface_land_ice_fluxes_tracers(meshPool, groupName, forcingPool, err = 0 - if ( .not. landIceFluxesOn ) return + if (landIceFluxesOff) return call mpas_timer_start("land_ice_" // trim(groupName)) @@ -132,73 +142,84 @@ end subroutine ocn_surface_land_ice_fluxes_tracers!}}} !> \author Xylar Asay-Davis !> \date 9 September 2015 !> \details -!> This routine computes the top-drag tendency for momentum -!> based on current state. +!> This routine computes the top-drag tendency under land ice for +!> momentum based on current state. ! !----------------------------------------------------------------------- - subroutine ocn_surface_land_ice_fluxes_vel(meshPool, diagnosticsPool, surfaceStress, surfaceStressMagnitude, err)!{{{ + subroutine ocn_surface_land_ice_fluxes_vel(topDrag, topDragMag, & + sfcStress, sfcStressMag, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- - type (mpas_pool_type), intent(in) :: meshPool !< Input: mesh information - type (mpas_pool_type), intent(in) :: diagnosticsPool !< Input: Diagnostics information + + real (kind=RKIND), dimension(:), intent(in) :: & + topDrag, &!< [in] drag at top of ocean on edge + topDragMag !< [in] magnitude of drag at top of ocean at center !----------------------------------------------------------------- - ! ! input/output variables - ! !----------------------------------------------------------------- - real (kind=RKIND), dimension(:), intent(inout) :: surfaceStress, & !< Input/Output: Array for total surface stress - surfaceStressMagnitude !< Input/Output: Array for magnitude of surface stress + + real (kind=RKIND), dimension(:), intent(inout) :: & + sfcStress, &!< [inout] accumulated surface stress + sfcStressMag !< [inout] accumulated surface stress magnitude !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: Error flag + integer, intent(out) :: err !< [out] Error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer :: iEdge, iCell - integer, pointer :: nCells, nEdges - real (kind=RKIND), dimension(:), pointer :: topDrag, topDragMagnitude + integer :: & + iEdge, iCell ! loop indices for edge, cell loops - err = 0 + ! End preamble + !----------------------------------------------------------------- + ! Begin code - if ( .not. landIceFluxesOn ) return + !*** set error flag and return early if land ice fluxes not on + !*** otherwise, start timer - call mpas_timer_start("top_drag", .false.) + err = 0 + if (landIceFluxesOff) return - call mpas_pool_get_dimension(meshPool, 'nCells', nCells) - call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) + call mpas_timer_start("top_drag", .false.) - call mpas_pool_get_array(diagnosticsPool, 'topDrag', topDrag) - call mpas_pool_get_array(diagnosticsPool, 'topDragMagnitude', topDragMagnitude) + !*** simply add input values to accumulated stresses + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(sfcStress, topDrag) + #else !$omp parallel !$omp do schedule(runtime) - do iEdge = 1, nEdges - surfaceStress(iEdge) = surfaceStress(iEdge) + topDrag(iEdge) + #endif + do iEdge = 1, nEdgesAll + sfcStress(iEdge) = sfcStress(iEdge) + topDrag(iEdge) end do + #ifndef MPAS_OPENACC !$omp end do + #endif - ! Build surface stress magnitude at cell centers + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(sfcStressMag, topDragMag) + #else !$omp do schedule(runtime) - do iCell = 1, nCells - surfaceStressMagnitude(iCell) = surfaceStressMagnitude(iCell) + topDragMagnitude(iCell) + #endif + do iCell = 1, nCellsAll + sfcStressMag(iCell) = sfcStressMag(iCell) + topDragMag(iCell) end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif call mpas_timer_stop("top_drag") @@ -257,7 +278,7 @@ subroutine ocn_surface_land_ice_fluxes_thick(meshPool, forcingPool, surfaceThick err = 0 - if ( .not. landIceFluxesOn ) return + if (landIceFluxesOff) return call mpas_timer_start("land_ice_thick") @@ -329,7 +350,7 @@ subroutine ocn_surface_land_ice_fluxes_active_tracers(meshPool, forcingPool, tra err = 0 - if ( .not. landIceFluxesOn ) return + if (landIceFluxesOff) return call mpas_pool_get_dimension(meshPool, 'nCells', nCells) @@ -436,7 +457,7 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, diagnosticsPool, & err = 0 - if ( .not. standaloneOn ) return + if (standaloneOff) return call mpas_timer_start("land_ice_build_arrays") @@ -470,7 +491,7 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, diagnosticsPool, & call mpas_pool_get_array(statePool, 'accumulatedLandIceHeat', accumulatedLandIceHeatNew, 2) call mpas_pool_get_array(statePool, 'accumulatedLandIceHeat', accumulatedLandIceHeatOld, 1) - if(config_land_ice_flux_useHollandJenkinsAdvDiff) then + if (useHollandJenkinsAdvDiff) then call mpas_pool_get_array(forcingPool, 'landIceSurfaceTemperature', landIceSurfaceTemperature) allocate(freezeInterfaceSalinity(nCells), & @@ -482,7 +503,8 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, diagnosticsPool, & nCells = nCellsArray( size(nCellsArray) ) - if(isomipOn) then + if (isomipOn) then !*** ISOMIP formulation + !$omp parallel !$omp do schedule(runtime) private(freshwaterFlux, heatFlux) do iCell = 1, nCells @@ -498,7 +520,7 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, diagnosticsPool, & ! or (7) from Jenkins et al. (2001) if gamma constant ! and no heat flux into ice ! freshwater flux = density * melt rate is in kg/m^2/s - freshwaterFlux = -rho_sw * config_land_ice_flux_ISOMIP_gammaT * (cp_sw/latent_heat_fusion_mks) & + freshwaterFlux = -rho_sw * ISOMIPgammaT * (cp_sw/latent_heat_fusion_mks) & * (landIceInterfaceTracers(indexIT,iCell)-landIceBoundaryLayerTracers(indexBLT,iCell)) landIceFreshwaterFlux(iCell) = landIceFraction(iCell)*freshwaterFlux @@ -506,7 +528,7 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, diagnosticsPool, & ! Using (13) from Jenkins et al. (2001) ! heat flux is in W/s heatFlux = cp_sw*(freshwaterFlux*landIceInterfaceTracers(indexIT,iCell) & - + rho_sw*config_land_ice_flux_ISOMIP_gammaT & + + rho_sw*ISOMIPgammaT & * (landIceInterfaceTracers(indexIT,iCell)-landIceBoundaryLayerTracers(indexBLT,iCell))) landIceHeatFlux(iCell) = landIceFraction(iCell)*heatFlux @@ -515,10 +537,10 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, diagnosticsPool, & end do !$omp end do !$omp end parallel - end if - - if(jenkinsOn .or. hollandJenkinsOn) then - if(config_land_ice_flux_useHollandJenkinsAdvDiff) then + endif ! isomipOn + + if (jenkinsOn .or. hollandJenkinsOn) then + if(useHollandJenkinsAdvDiff) then ! melting solution call compute_HJ99_melt_fluxes( & landIceMask, & @@ -602,9 +624,9 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, diagnosticsPool, & heatFluxToLandIce(iCell) = landIceFraction(iCell)*heatFluxToLandIce(iCell) end do - end if + endif ! jenkinsOn or hollandJenkinsOn - if(config_land_ice_flux_useHollandJenkinsAdvDiff) then + if(useHollandJenkinsAdvDiff) then deallocate(freezeInterfaceSalinity, & freezeInterfaceTemperature, & freezeFreshwaterFlux, & @@ -639,35 +661,70 @@ end subroutine ocn_surface_land_ice_fluxes_build_arrays!}}} subroutine ocn_surface_land_ice_fluxes_init(err)!{{{ - integer, intent(out) :: err !< Output: error flag + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< [out] error flag - err = 0 - isomipOn = .false. - jenkinsOn = .false. - hollandJenkinsOn = .false. + ! End preamble + !----------------------------------------------------------------- + ! Begin code - landIceFluxesOn = (trim(config_land_ice_flux_mode) == 'standalone') & - .or. (trim(config_land_ice_flux_mode) == 'coupled') - if(.not. landIceFluxesOn) return + !*** Set error code and default values for module variables - standaloneOn = trim(config_land_ice_flux_mode) == 'standalone' + err = 0 - if ( trim(config_land_ice_flux_formulation) == 'ISOMIP' ) then + landIceFluxesOff = .true. + standaloneOff = .true. + isomipOn = .false. + jenkinsOn = .false. + hollandJenkinsOn = .false. + cpLandIce = 0.0_RKIND + rhoLandIce = 0.0_RKIND + ISOMIPgammaT = 0.0_RKIND + + !*** Determine whether land ice fluxes are on + + select case (trim(config_land_ice_flux_mode)) + case ('standalone','Standalone','STANDALONE') + landIceFluxesOff = .false. + standaloneOff = .false. + case ('coupled','Coupled','COUPLED') + landIceFluxesOff = .false. + standaloneOff = .true. + case default + landIceFluxesOff = .true. + standaloneOff = .true. + return + end select + + !*** Determine which flux formulation to use + + select case (trim(config_land_ice_flux_formulation)) + + case ('isomip','Isomip','ISOMIP') isomipOn = .true. - else if ( trim(config_land_ice_flux_formulation) == 'Jenkins' ) then + case ('jenkins', 'Jenkins','JENKINS') jenkinsOn = .true. - else if ( trim(config_land_ice_flux_formulation) == 'HollandJenkins' ) then + case ('hollandjenkins','HollandJenkins','HOLLANDJENKINS') hollandJenkinsOn = .true. - else + case default call mpas_log_write( & "config_land_ice_flux_formulation not one of 'ISOMIP', 'Jenkins', " & // "or 'HollandJenkins'.", & MPAS_LOG_CRIT) err = 1 - end if + end select + + useHollandJenkinsAdvDiff = & + config_land_ice_flux_useHollandJenkinsAdvDiff + + !*** Set physical values - cp_land_ice = config_land_ice_flux_cp_ice - rho_land_ice = config_land_ice_flux_rho_ice + cpLandIce = config_land_ice_flux_cp_ice + rhoLandIce = config_land_ice_flux_rho_ice + ISOMIPgammaT = config_land_ice_flux_ISOMIP_gammaT !-------------------------------------------------------------------- @@ -807,7 +864,7 @@ subroutine compute_melt_fluxes( & ! 2. the diffusion (if any) of heat into the ice, based on temperature difference between ! the reference point in the ice (either the surface or the middle of the bottom layer) ! and the interface - outIceHeatFlux(iCell) = -cp_land_ice*outFreshwaterFlux(iCell)*outInterfaceTemperature(iCell) + outIceHeatFlux(iCell) = -cpLandIce*outFreshwaterFlux(iCell)*outInterfaceTemperature(iCell) end do !$omp end do !$omp end parallel @@ -911,7 +968,7 @@ subroutine compute_HJ99_melt_fluxes( & integer :: iCell err = 0 - cpRatio = cp_land_ice/cp_sw + cpRatio = cpLandIce/cp_sw !$omp parallel !$omp do schedule(runtime) & !$omp private(T0, dTf_dS, transferVelocityRatio, Tlatent, eta, TlatentStar, a, b, c, err) @@ -953,7 +1010,7 @@ subroutine compute_HJ99_melt_fluxes( & ! Since we're considering only melting and ignoring diffusion, ! the ice loses heat simply by the loss of ice mass at the prescribed ! (surface?) ice temperature - outIceHeatFlux(iCell) = -cp_land_ice*outFreshwaterFlux(iCell)*iceTemperature(iCell) + outIceHeatFlux(iCell) = -cpLandIce*outFreshwaterFlux(iCell)*iceTemperature(iCell) end do !$omp end do !$omp end parallel diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index 7b6977ec54..a9c94f770e 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -28,6 +28,7 @@ module ocn_tendency use ocn_diagnostics use ocn_constants use ocn_config + use ocn_mesh use ocn_surface_bulk_forcing use ocn_surface_land_ice_fluxes @@ -109,7 +110,6 @@ module ocn_tendency !----------------------------------------------------------------------- subroutine ocn_tend_thick(tendPool, forcingPool, diagnosticsPool, meshPool)!{{{ - implicit none type (mpas_pool_type), intent(inout) :: tendPool !< Input/Output: Tendency structure type (mpas_pool_type), intent(inout) :: forcingPool !< Input: Forcing information @@ -208,175 +208,304 @@ end subroutine ocn_tend_thick!}}} ! !----------------------------------------------------------------------- - subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshPool, scratchPool, timeLevelIn, dt)!{{{ - implicit none + subroutine ocn_tend_vel(tendPool, statePool, forcingPool, & + diagnosticsPool, timeLevelIn, dt)!{{{ - type (mpas_pool_type), intent(inout) :: tendPool !< Input/Output: Tendency structure - type (mpas_pool_type), intent(in) :: statePool !< Input: State information - type (mpas_pool_type), intent(inout) :: forcingPool !< Input: Forcing information - type (mpas_pool_type), intent(in) :: diagnosticsPool !< Input: Diagnostic information - type (mpas_pool_type), intent(in) :: meshPool !< Input: Mesh information - real (kind=RKIND), intent(in) :: dt - type (mpas_pool_type), intent(inout) :: scratchPool !< Input: Scratch structure - integer, intent(in), optional :: timeLevelIn !< Input: Time level for state fields + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + real (kind=RKIND), intent(in) :: & + dt !< [in] time step + + integer, intent(in), optional :: & + timeLevelIn !< [in] Time level for state fields + + ! Input structures with lots of fields + type (mpas_pool_type), intent(in) :: statePool ! replace + type (mpas_pool_type), intent(in) :: diagnosticsPool ! replace + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + ! currently updated as pool variables - change to args? + + type (mpas_pool_type), intent(inout) :: & + tendPool, &!< [out] Tendency structure w/ vel tend + forcingPool !< [out] Forcing structure w/ sfc stresses + + real (kind=RKIND), dimension(:,:), pointer :: & + tendVel !< [out] normal velocity tendency at edges + + real (kind=RKIND), dimension(:), pointer :: & + sfcStress, &!< [out] surface stress at edges + sfcStressMag !< [out] surface stress magnitude (on cell) + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: & + err, &! local error flag + iEdge, iCell, k, &! loop indices for edge, cell and vertical + indxTemp, &! tracer index for temperature + indxSalt, &! tracer index for salinity + timeLevel ! time level to use for state variables + + ! various pointers for array retrievals + integer, pointer :: indexTemperature, indexSalinity type (mpas_pool_type), pointer :: tracersPool - real (kind=RKIND), dimension(:), pointer :: surfaceStress, surfaceStressMagnitude, surfaceFluxAttenuationCoefficient, ssh - real (kind=RKIND), dimension(:), pointer :: tidalPotentialEta + real (kind=RKIND), dimension(:), pointer :: & + windStressZonal, &! zonal wind stress from forcing + windStressMerid, &! meridional wind stress from forcing + topDrag, &! drag at top surface + topDragMag, &! magnitude of drag at top + sfcFlxAttCoeff, &! surface flux attenuation coefficient + ssh ! sea surface height + real (kind=RKIND), dimension(:,:), pointer :: & - layerThicknessEdge, normalVelocity, tangentialVelocity, density, potentialDensity, zMid, pressure, & - layerThickness, & - tend_normalVelocity, circulation, relativeVorticity, viscosity, kineticEnergyCell, & - normalizedRelativeVorticityEdge, normalizedPlanetaryVorticityEdge, & - montgomeryPotential, vertAleTransportTop, divergence, vertViscTopOfEdge, & - inSituThermalExpansionCoeff, inSituSalineContractionCoeff - real (kind=RKIND), dimension(:,:), pointer :: tidalPotentialZMid - real (kind=RKIND), dimension(:,:), pointer :: wettingVelocity - real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers + layerThickEdge, &! layer thicknesss at edge + normalVelocity, &! normal velocity + tangentialVelocity, &! tangential velocity + density, & + potentialDensity, & + zMid, & + pressure, & + layerThickness, & + circulation, & + relativeVorticity, & + kineticEnergyCell, & + normRelVortEdge, &! normalized relative vorticity on edge + normPlanetVortEdge, &! normalized planetary vorticity on edge + montgomeryPotential, & + vertAleTransportTop, & + divergence, & + vertViscTopOfEdge, & + thermExpCoeff, &! in situ thermal expansion coeff + salineContractCoeff, &! in situ saline contraction coeff + wettingVelocity - integer :: timeLevel + real (kind=RKIND), dimension(:,:,:), pointer :: & + activeTracers - integer :: err, iEdge, iCell, k - integer, pointer :: indexTemperature, indexSalinity, nEdges, nCells - integer, dimension(:), pointer :: maxLevelCell + ! End preamble + !----------------------------------------------------------------- + ! Begin code + !*** Initialize error code and determine time level + + err = 0 if (present(timeLevelIn)) then timeLevel = timeLevelIn else timeLevel = 1 end if - call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) - call mpas_pool_get_dimension(meshPool, 'nCells', nCells) - + !*** Retrieve pool variables call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) - call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, timeLevel) - call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevel) - call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevel) - call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracers, timeLevel) - call mpas_pool_get_dimension(tracersPool, 'index_temperature', indexTemperature) - call mpas_pool_get_dimension(tracersPool, 'index_salinity', indexSalinity) - - call mpas_pool_get_array(diagnosticsPool, 'kineticEnergyCell', kineticEnergyCell) - call mpas_pool_get_array(diagnosticsPool, 'layerThicknessEdge', layerThicknessEdge) - call mpas_pool_get_array(diagnosticsPool, 'vertAleTransportTop', vertAleTransportTop) + call mpas_pool_get_array(statePool, 'normalVelocity', & + normalVelocity, timeLevel) + call mpas_pool_get_array(statePool, 'layerThickness', & + layerThickness, timeLevel) + call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevel) + call mpas_pool_get_array(tracersPool, 'activeTracers', & + activeTracers, timeLevel) + + call mpas_pool_get_dimension(tracersPool, 'index_temperature', & + indexTemperature) + call mpas_pool_get_dimension(tracersPool, 'index_salinity', & + indexSalinity) + indxTemp = indexTemperature + indxSalt = indexSalinity + + call mpas_pool_get_array(diagnosticsPool, 'kineticEnergyCell', & + kineticEnergyCell) + call mpas_pool_get_array(diagnosticsPool, 'layerThicknessEdge', & + layerThickEdge) + call mpas_pool_get_array(diagnosticsPool, 'vertAleTransportTop', & + vertAleTransportTop) call mpas_pool_get_array(diagnosticsPool, 'zMid', zMid) - call mpas_pool_get_array(diagnosticsPool, 'relativeVorticity', relativeVorticity) - call mpas_pool_get_array(diagnosticsPool, 'normalizedRelativeVorticityEdge', normalizedRelativeVorticityEdge) - call mpas_pool_get_array(diagnosticsPool, 'normalizedPlanetaryVorticityEdge', normalizedPlanetaryVorticityEdge) - call mpas_pool_get_array(diagnosticsPool, 'divergence', divergence) - call mpas_pool_get_array(diagnosticsPool, 'viscosity', viscosity) - call mpas_pool_get_array(diagnosticsPool, 'montgomeryPotential', montgomeryPotential) - call mpas_pool_get_array(diagnosticsPool, 'pressure', pressure) - call mpas_pool_get_array(diagnosticsPool, 'vertViscTopOfEdge', vertViscTopOfEdge) - call mpas_pool_get_array(diagnosticsPool, 'density', density) - call mpas_pool_get_array(diagnosticsPool, 'potentialDensity', potentialDensity) - call mpas_pool_get_array(diagnosticsPool, 'tangentialVelocity', tangentialVelocity) - call mpas_pool_get_array(diagnosticsPool, 'surfaceFluxAttenuationCoefficient', surfaceFluxAttenuationCoefficient) - - call mpas_pool_get_array(tendPool, 'normalVelocity', tend_normalVelocity) - - call mpas_pool_get_array(forcingPool, 'surfaceStress', surfaceStress) - call mpas_pool_get_array(forcingPool, 'surfaceStressMagnitude', surfaceStressMagnitude) - - call mpas_pool_get_array(diagnosticsPool, 'wettingVelocity', wettingVelocity) - ! - ! velocity tendency: start accumulating tendency terms - ! - !$omp Parallel - !$omp do schedule(runtime) - do iEdge = 1, nEdges - tend_normalVelocity(:, iEdge) = 0.0_RKIND - surfaceStress(iEdge) = 0.0_RKIND + call mpas_pool_get_array(diagnosticsPool, 'relativeVorticity', & + relativeVorticity) + call mpas_pool_get_array(diagnosticsPool, 'normalizedRelativeVorticityEdge', & + normRelVortEdge) + call mpas_pool_get_array(diagnosticsPool, 'normalizedPlanetaryVorticityEdge', & + normPlanetVortEdge) + call mpas_pool_get_array(diagnosticsPool, 'divergence', & + divergence) + call mpas_pool_get_array(diagnosticsPool, 'montgomeryPotential', & + montgomeryPotential) + call mpas_pool_get_array(diagnosticsPool, 'pressure', & + pressure) + call mpas_pool_get_array(diagnosticsPool, 'vertViscTopOfEdge', & + vertViscTopOfEdge) + call mpas_pool_get_array(diagnosticsPool, 'density', & + density) + call mpas_pool_get_array(diagnosticsPool, 'potentialDensity', & + potentialDensity) + call mpas_pool_get_array(diagnosticsPool, 'tangentialVelocity', & + tangentialVelocity) + call mpas_pool_get_array(diagnosticsPool, 'surfaceFluxAttenuationCoefficient', & + sfcFlxAttCoeff) + call mpas_pool_get_array(diagnosticsPool, 'inSituThermalExpansionCoeff', & + thermExpCoeff) + call mpas_pool_get_array(diagnosticsPool, 'inSituSalineContractionCoeff', & + salineContractCoeff) + call mpas_pool_get_array(diagnosticsPool, 'topDrag', & + topDrag) + call mpas_pool_get_array(diagnosticsPool, 'topDragMagnitude', & + topDragMag) + call mpas_pool_get_array(diagnosticsPool, 'wettingVelocity', & + wettingVelocity) + call mpas_pool_get_array(forcingPool, 'windStressZonal', & + windStressZonal) + call mpas_pool_get_array(forcingPool, 'windStressMeridional', & + windStressMerid) + + call mpas_pool_get_array(tendPool, 'normalVelocity', tendVel) + + call mpas_pool_get_array(forcingPool, 'surfaceStress', & + sfcStress) + call mpas_pool_get_array(forcingPool, 'surfaceStressMagnitude', & + sfcStressMag) + + !*** Transfer data to device + !$acc enter data & + !$acc copyin(tendVel, sfcStress, sfcStressMag, & + !$acc windStressZonal, windStressMerid, topDrag, topDragMag, & + !$acc sfcFlxAttCoeff, ssh, layerThickEdge, normalVelocity, & + !$acc tangentialVelocity, density, potentialDensity, zMid, & + !$acc pressure, layerThickness, circulation, activeTracers, & + !$acc relativeVorticity, kineticEnergyCell, normRelVortEdge, & + !$acc normPlanetVortEdge, montgomeryPotential, & + !$acc vertAleTransportTop, divergence, vertViscTopOfEdge, & + !$acc thermExpCoeff, salineContractCoeff, wettingVelocity) + + !*** Set output variables to zero before accumulating results + + #ifdef MPAS_OPENACC + !$acc parallel loop present(sfcStress, tendVel) private(k) + #else + !$omp parallel + !$omp do schedule(runtime) private(k) + #endif + do iEdge = 1, nEdgesAll + sfcStress(iEdge) = 0.0_RKIND + do k=1,nVertLevels + tendVel(k,iEdge) = 0.0_RKIND + end do end do + #ifndef MPAS_OPENACC !$omp end do + #endif + #ifdef MPAS_OPENACC + !$acc parallel loop present(sfcStressMag) + #else + !$omp end do !$omp do schedule(runtime) - do iCell = 1, nCells - surfaceStressMagnitude(iCell) = 0.0_RKIND + #endif + do iCell = 1, nCellsAll + sfcStressMag(iCell) = 0.0_RKIND end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif - if(config_disable_vel_all_tend) return + !*** Return early if all vel tendencies disabled + !*** Otherwise, start timer + if (config_disable_vel_all_tend) return call mpas_timer_start("ocn_tend_vel") - ! Build bulk forcing surface stress - call ocn_surface_bulk_forcing_vel(meshPool, forcingPool, surfaceStress, surfaceStressMagnitude, err) + ! Compute bulk forcing surface stress + call ocn_surface_bulk_forcing_vel( & + windStressZonal, windStressMerid, & + sfcStress, sfcStressMag, err) ! Add top drag to surface stress - call ocn_surface_land_ice_fluxes_vel(meshPool, diagnosticsPool, surfaceStress, surfaceStressMagnitude, err) - - ! - ! velocity tendency: nonlinear Coriolis term and grad of kinetic energy - ! - call ocn_vel_hadv_coriolis_tend(meshPool, normalizedRelativeVorticityEdge, normalizedPlanetaryVorticityEdge, layerThicknessEdge, & - normalVelocity, kineticEnergyCell, tend_normalVelocity, err) - - ! - ! velocity tendency: vertical advection term -w du/dz - ! - call ocn_vel_vadv_tend(meshPool, normalVelocity, layerThicknessEdge, vertAleTransportTop, tend_normalVelocity, err) - - ! - ! velocity tendency: pressure gradient - ! - if (config_pressure_gradient_type.eq.'Jacobian_from_TS') then - ! only pass EOS derivatives if needed. - call mpas_pool_get_array(diagnosticsPool, 'inSituThermalExpansionCoeff',inSituThermalExpansionCoeff) - call mpas_pool_get_array(diagnosticsPool, 'inSituSalineContractionCoeff', inSituSalineContractionCoeff) - call ocn_vel_pressure_grad_tend(meshPool, ssh, pressure, montgomeryPotential, zMid, density, potentialDensity, & - indexTemperature, indexSalinity, activeTracers, tend_normalVelocity, err, & - inSituThermalExpansionCoeff,inSituSalineContractionCoeff) - else - call ocn_vel_pressure_grad_tend(meshPool, ssh, pressure, montgomeryPotential, zMid, density, potentialDensity, & - indexTemperature, indexSalinity, activeTracers, tend_normalVelocity, err, & - inSituThermalExpansionCoeff,inSituSalineContractionCoeff) - endif - - ! - ! velocity tendency: tidal potential (if needed) - ! - call ocn_compute_tidal_potential_forcing(meshPool, forcingPool, diagnosticsPool, err) + call ocn_surface_land_ice_fluxes_vel(topDrag, topDragMag, & + sfcStress, sfcStressMag, err) + + ! Add nonlinear Coriolis term and horizontal advection of + ! momentum, formulated as grad of kinetic energy + call ocn_vel_hadv_coriolis_tend(normRelVortEdge, & + normPlanetVortEdge, & + layerThickEdge, normalVelocity, & + kineticEnergyCell, tendVel, err) + + ! Add vertical advection term -w du/dz + call ocn_vel_vadv_tend(normalVelocity, layerThickEdge, & + vertAleTransportTop, tendVel, err) + + ! Add pressure gradient + call ocn_vel_pressure_grad_tend(ssh, pressure, & + montgomeryPotential, zMid, & + density, potentialDensity, & + indxTemp, indxSalt, activeTracers, & + thermExpCoeff, salineContractCoeff, & + tendVel, err) + + ! Add tidal potential (if needed) + call ocn_compute_tidal_potential_forcing(diagnosticsPool, err) if (config_time_integrator == 'RK4') then - ! for split explicit, tidal forcing is added in barotropic subsycles - call ocn_vel_tidal_potential_tend(meshPool, forcingPool, ssh, tend_normalVelocity, err) + ! for split explicit, tidal forcing is added in barotropic subcycles + call ocn_vel_tidal_potential_tend(ssh, tendVel, err) endif - ! - ! velocity tendency: del2 dissipation, \nu_2 \nabla^2 u - ! computed as \nu( \nabla divergence + k \times \nabla relativeVorticity ) - ! strictly only valid for config_mom_del2 == constant - ! - call ocn_vel_hmix_tend(meshPool, scratchPool, divergence, relativeVorticity, normalVelocity, tangentialVelocity, viscosity, & - tend_normalVelocity, err) + ! Add horizontal mixing + call ocn_vel_hmix_tend(divergence, relativeVorticity, tendVel,err) - ! - ! velocity tendency: forcing and bottom drag - ! - call ocn_vel_forcing_tend(meshPool, normalVelocity, surfaceFluxAttenuationCoefficient, & - surfaceStress, kineticEnergyCell, layerThicknessEdge, & - tend_normalVelocity, err) + ! Add forcing and bottom drag + call ocn_vel_forcing_tend(normalVelocity, sfcFlxAttCoeff, & + sfcStress, kineticEnergyCell, & + layerThickEdge, tendVel, err) - ! - ! velocity tendency: zero if drying - ! + ! vertical mixing treated implicitly in a later routine + ! adjust total velocity tendency based on wetting/drying + + #ifdef MPAS_OPENACC + !$acc parallel loop collapse(2) & + !$acc present(tendVel, wettingVelocity) + #else !$omp parallel - !$omp do schedule(runtime) - do iEdge = 1, nEdges - tend_normalVelocity(:, iEdge) = tend_normalVelocity(:, iEdge) * (1.0_RKIND - wettingVelocity(:, iEdge)) + !$omp do schedule(runtime) private(k) + #endif + do iEdge = 1, nEdgesAll + do k=1,nVertLevels + tendVel(k,iEdge) = tendVel(k,iEdge)* & + (1.0_RKIND - wettingVelocity(k,iEdge)) + end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif + + !*** Transfer any needed data back to host and delete inputs + !$acc exit data & + !$acc copyout(tendVel, sfcStress, sfcStressMag) & + !$acc delete( & + !$acc windStressZonal, windStressMerid, topDrag, topDragMag, & + !$acc sfcFlxAttCoeff, ssh, layerThickEdge, normalVelocity, & + !$acc tangentialVelocity, density, potentialDensity, zMid, & + !$acc pressure, layerThickness, circulation, activeTracers, & + !$acc relativeVorticity, kineticEnergyCell, normRelVortEdge, & + !$acc normPlanetVortEdge, montgomeryPotential, & + !$acc vertAleTransportTop, divergence, vertViscTopOfEdge, & + !$acc thermExpCoeff, salineContractCoeff, wettingVelocity) + + ! stop the timer and exit - ! - ! velocity tendency: vertical mixing d/dz( nu_v du/dz)) - ! call mpas_timer_stop("ocn_tend_vel") + !-------------------------------------------------------------------- + end subroutine ocn_tend_vel!}}} !*********************************************************************** diff --git a/src/core_ocean/shared/mpas_ocn_vel_forcing.F b/src/core_ocean/shared/mpas_ocn_vel_forcing.F index 58d0d55eb5..ef82187bf6 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_forcing.F +++ b/src/core_ocean/shared/mpas_ocn_vel_forcing.F @@ -20,11 +20,10 @@ module ocn_vel_forcing + use mpas_kind_types use mpas_derived_types - + use mpas_log use ocn_constants - use ocn_forcing - use ocn_vel_forcing_surface_stress use ocn_vel_forcing_explicit_bottom_drag @@ -75,71 +74,56 @@ module ocn_vel_forcing ! !----------------------------------------------------------------------- - subroutine ocn_vel_forcing_tend(meshPool, normalVelocity, surfaceFluxAttenuationCoefficient, & - surfaceStress, kineticEnergyCell, layerThicknessEdge, tend, err)!{{{ + subroutine ocn_vel_forcing_tend(normalVelocity, sfcFlxAttCoeff, & + surfaceStress, kineticEnergyCell, & + layerThicknessEdge, tend, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(in) :: & - normalVelocity, & !< Input: Normal velocity at edges - kineticEnergyCell !< Input: kinetic energy at cell + normalVelocity, &!< [in] Normal velocity at edges + kineticEnergyCell, &!< [in] kinetic energy at cell + layerThicknessEdge !< [in] thickness at edge real (kind=RKIND), dimension(:), intent(in) :: & - surfaceFluxAttenuationCoefficient, & !< Input: attenuation coefficient for surface fluxes at cell centers - surfaceStress !< Input: surface stress at edges - - real (kind=RKIND), dimension(:,:), intent(in) :: & - layerThicknessEdge !< Input: thickness at edge - - type (mpas_pool_type), intent(in) :: & - meshPool !< Input: mesh information + sfcFlxAttCoeff, &!< [in] attenuation coefficient for sfc fluxes + surfaceStress !< [in] surface stress at edges !----------------------------------------------------------------- - ! ! input/output variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(inout) :: & - tend !< Input/Output: velocity tendency + tend !< [inout] accumulated velocity tendency !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer :: err1 + integer :: err1 ! local error flag + ! End preamble !----------------------------------------------------------------- - ! - ! call relevant routines for computing tendencies - ! note that the user can choose multiple options and the - ! tendencies will be added together - ! - !----------------------------------------------------------------- + ! Begin code + + !*** initialize error flag, then call individual forcing routines err = 0 - call ocn_vel_forcing_surface_stress_tend(meshPool, surfaceFluxAttenuationCoefficient, & - surfaceStress, layerThicknessEdge, tend, err1) + call ocn_vel_forcing_surface_stress_tend(sfcFlxAttCoeff, & + surfaceStress, layerThicknessEdge, tend, err1) err = ior(err, err1) - call ocn_vel_forcing_explicit_bottom_drag_tend(meshPool, normalVelocity, & - kineticEnergyCell, layerThicknessEdge, tend, err1) - + call ocn_vel_forcing_explicit_bottom_drag_tend(normalVelocity, & + kineticEnergyCell, layerThicknessEdge, tend, err1) err = ior(err, err1) !-------------------------------------------------------------------- @@ -163,22 +147,39 @@ end subroutine ocn_vel_forcing_tend!}}} subroutine ocn_vel_forcing_init(err)!{{{ - !-------------------------------------------------------------------- + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! - ! call individual init routines for each parameterization - ! + ! local variables !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer :: err1, err2 ! local error flags - integer :: err1, err2 + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + err = 0 call ocn_vel_forcing_surface_stress_init(err1) - call ocn_vel_forcing_explicit_bottom_drag_init(err2) + if (err1 /= 0) then + call mpas_log_write( & + 'vel_forcing_init encountered error in sfc stress init', & + MPAS_LOG_ERR) + err = err1 + endif - err = ior(err1, err2) + call ocn_vel_forcing_explicit_bottom_drag_init(err2) + if (err2 /= 0) then + call mpas_log_write( & + 'vel_forcing_init encountered error in bot drag init', & + MPAS_LOG_ERR) + err = err2 + endif !-------------------------------------------------------------------- diff --git a/src/core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F b/src/core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F index 69dc8c18af..aed856924f 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F +++ b/src/core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F @@ -20,12 +20,11 @@ module ocn_vel_forcing_explicit_bottom_drag - use mpas_derived_types - use mpas_pool_routines use mpas_timer use ocn_constants use ocn_config + use ocn_mesh use ocn_forcing implicit none @@ -53,8 +52,11 @@ module ocn_vel_forcing_explicit_bottom_drag ! !-------------------------------------------------------------------- - logical :: explicitBottomDragOn - real (kind=RKIND) :: explicitBottomDragCoef + logical :: & + explicitBottomDragOff ! on/off switch for explicit bottom drag + + real (kind=RKIND) :: & + dragCoeff ! drag coefficient !*********************************************************************** @@ -73,83 +75,81 @@ module ocn_vel_forcing_explicit_bottom_drag ! !----------------------------------------------------------------------- - subroutine ocn_vel_forcing_explicit_bottom_drag_tend(meshPool, normalVelocity, & !{{{ - kineticEnergyCell, layerThicknessEdge, tend, err) + subroutine ocn_vel_forcing_explicit_bottom_drag_tend(normVelocity, & + KECell, layerThickEdge, tend, err) !{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(in) :: & - normalVelocity, &!< Input: velocity - kineticEnergyCell, &!< Input: kinetic energy at cell - layerThicknessEdge !< Input: thickness at edge - - type (mpas_pool_type), intent(in) :: & - meshPool !< Input: mesh information + normVelocity, &!< [in] normal velocity + KECell, &!< [in] kinetic energy at cell + layerThickEdge !< [in] layer thickness at edge !----------------------------------------------------------------- - ! ! input/output variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(inout) :: & - tend !< Input/Output: velocity tendency + tend !< [inout] accumulated velocity tendency !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer :: iEdge, k, cell1, cell2, nEdges - integer, dimension(:), pointer :: nEdgesArray - integer, dimension(:), pointer :: maxLevelEdgeTop - integer, dimension(:,:), pointer :: cellsOnEdge + integer :: & + iEdge, &! loop index for edge loop + k, &! vertical index of lowest active layer at edge + cell1, cell2 ! neighbor cell addresses across edge - err = 0 + ! End preamble + !----------------------------------------------------------------- + ! Begin code - if ( .not. explicitBottomDragOn ) return + !*** Initialize error code and return if not turned on + !*** Otherwise start timer + err = 0 + if (explicitBottomDragOff) return call mpas_timer_start('vel explicit bottom drag') - call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) - call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - - nEdges = nEdgesArray( 1 ) - + ! Explicit bottom drag term: + ! du/dt = ... - c |u| u / h + ! appied to bottom layer only. + ! This term comes from the bottom boundary condition in the vertical + ! momentum mixing, and is explicit if both |u| and u are chosen to be at + ! time level n. + + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, maxLevelEdgeTop, KECell, & + !$acc tend, normVelocity, layerThickEdge) & + !$acc private(k, cell1, cell2) + #else !$omp parallel !$omp do schedule(runtime) private(k, cell1, cell2) - do iEdge = 1, nEdges + #endif + do iEdge = 1, nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) k = maxLevelEdgeTop(iEdge) - ! Explicit bottom drag term: - ! du/dt = ... - c |u| u / h - ! appied to bottom layer only. - ! This term comes from the bottom boundary condition in the vertical - ! momentum mixing, and is explicit if both |u| and u are chosen to be at - ! time level n. - - tend(k,iEdge) = tend(k,iEdge) - explicitBottomDragCoef * & - sqrt(kineticEnergyCell(k,cell1) + kineticEnergyCell(k,cell2)) * normalVelocity(k,iEdge) / layerThicknessEdge(k,iEdge) + tend(k,iEdge) = tend(k,iEdge) - dragCoeff* & + sqrt(KECell(k,cell1) + KECell(k,cell2))* & + normVelocity(k,iEdge)/layerThickEdge(k,iEdge) enddo + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif call mpas_timer_stop('vel explicit bottom drag') @@ -172,26 +172,33 @@ end subroutine ocn_vel_forcing_explicit_bottom_drag_tend!}}} subroutine ocn_vel_forcing_explicit_bottom_drag_init(err)!{{{ - !-------------------------------------------------------------------- - !----------------------------------------------------------------- - ! - ! call individual init routines for each parameterization - ! + ! output variables !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + !*** Initialize return error code and set module defaults err = 0 - explicitBottomDragCoef = 0.0_RKIND + explicitBottomDragOff = .true. + dragCoeff = 0.0_RKIND + + !*** Reset values based on input model configuration if (config_use_explicit_bottom_drag) then - explicitBottomDragOn = .true. - explicitBottomDragCoef = config_explicit_bottom_drag_coeff + explicitBottomDragOff = .false. + dragCoeff = config_explicit_bottom_drag_coeff endif - if (config_disable_vel_explicit_bottom_drag) explicitBottomDragOn = .false. + if (config_disable_vel_explicit_bottom_drag) then + explicitBottomDragOff = .true. + dragCoeff = 0.0_RKIND + endif !-------------------------------------------------------------------- diff --git a/src/core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F b/src/core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F index ca79b604c4..4d30e92241 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F +++ b/src/core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F @@ -20,13 +20,12 @@ module ocn_vel_forcing_surface_stress - use mpas_derived_types - use mpas_pool_routines use mpas_timer use ocn_constants use ocn_config use ocn_forcing + use ocn_mesh implicit none private @@ -53,7 +52,12 @@ module ocn_vel_forcing_surface_stress ! !-------------------------------------------------------------------- - logical :: surfaceStressOn + logical :: & + surfaceStressOff !< on/off switch for surface stress + + real (kind=RKIND), parameter :: & + maxAttDepth = -100.0_RKIND ! max attenuation depth to prevent + ! underflow in exponential !*********************************************************************** @@ -72,112 +76,118 @@ module ocn_vel_forcing_surface_stress ! !----------------------------------------------------------------------- - subroutine ocn_vel_forcing_surface_stress_tend(meshPool, surfaceFluxAttenuationCoefficient, surfaceStress, & !{{{ - layerThicknessEdge, tend, err) + subroutine ocn_vel_forcing_surface_stress_tend(sfcFlxAttCoeff, & + sfcStress, layerThickEdge, tend, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:), intent(in) :: & - surfaceStress, & !< Input: Wind stress at surface - surfaceFluxAttenuationCoefficient !< Input: attenuation coefficient for surface fluxes + sfcStress, & !< [in] Wind stress at surface + sfcFlxAttCoeff !< [in] attenuation coefficient for sfc fluxes real (kind=RKIND), dimension(:,:), intent(in) :: & - layerThicknessEdge !< Input: thickness at edge - - type (mpas_pool_type), intent(in) :: & - meshPool !< Input: mesh information + layerThickEdge !< [in] thickness at edge !----------------------------------------------------------------- - ! ! input/output variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(inout) :: & - tend !< Input/Output: velocity tendency + tend !< [inout] accumulated velocity tendency !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer :: iEdge, k, cell1, cell2, nEdges - integer, dimension(:), pointer :: nEdgesArray - integer, dimension(:), pointer :: maxLevelEdgeTop - integer, dimension(:,:), pointer :: edgeMask, cellsOnEdge - - real (kind=RKIND) :: transmissionCoeffTop, transmissionCoeffBot, zTop, zBot, remainingStress, & - attenuationCoeff - - !----------------------------------------------------------------- - ! - ! call relevant routines for computing tendencies - ! note that the user can choose multiple options and the - ! tendencies will be added together - ! + integer :: & + iEdge, k, &! loop indices for edge, vertical loops + kmax, &! index of deepest active edge + cell1, cell2 ! neighbor cell addresses across edge + + real (kind=RKIND) :: & + transCoeffTop, &! transmission coefficent at top of edge + transCoeffBot, &! transmission coefficent at bottom of edge + zTop, zBot, &! depths at top and bottom of edge + remainingStress, &! stress remaining at bottom of edge + absorb, &! absorption fraction within layer + attDepth, &! depth scaled by attenuation length + attCoeff ! attenuation coefficient on edge + + ! End preamble !----------------------------------------------------------------- + ! Begin code - err = 0 - - if ( .not. surfaceStressOn ) return + !*** Initialize error code and return if not on + !*** If turned on, start timer. + err = 0 + if (surfaceStressOff) return call mpas_timer_start('vel surface stress') - call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) - - call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'edgeMask', edgeMask) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - - nEdges = nEdgesArray ( 1 ) - + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, maxLevelEdgeTop, edgeMask, & + !$acc sfcFlxAttCoeff, sfcStress, layerThickEdge, & + !$acc tend) & + !$acc private(k, kmax, cell1, cell2, zBot, zTop, & + !$acc attCoeff, attDepth, remainingStress, & + !$acc transCoeffTop, transCoeffBot, absorb) + #else !$omp parallel !$omp do schedule(runtime) & - !$omp private(zTop, cell1, cell2, attenuationCoeff, transmissionCoeffTop, & - !$omp remainingStress, k, zBot, transmissionCoeffBot) - do iEdge = 1, nEdges - zTop = 0.0_RKIND + !$omp private(k, kmax, cell1, cell2, zBot, zTop, & + !$omp attCoeff, attDepth, remainingStress, & + !$omp transCoeffTop, transCoeffBot, absorb) + #endif + do iEdge = 1, nEdgesOwned + zTop = 0.0_RKIND cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) - attenuationCoeff = 0.5_RKIND * (surfaceFluxAttenuationCoefficient(cell1) & - + surfaceFluxAttenuationCoefficient(cell2)) - transmissionCoeffTop = ocn_forcing_transmission(zTop, attenuationCoeff) + attCoeff = 0.5_RKIND*(sfcFlxAttCoeff(cell1) & + + sfcFlxAttCoeff(cell2)) + attDepth = 0.0_RKIND + transCoeffTop = 1.0_RKIND remainingStress = 1.0_RKIND - do k = 1, maxLevelEdgeTop(iEdge) - zBot = zTop - layerThicknessEdge(k, iEdge) + kmax = maxLevelEdgeTop(iEdge) + do k = 1, kmax + zBot = zTop - layerThickEdge(k,iEdge) + attDepth = max(zBot/attCoeff, maxAttDepth) - transmissionCoeffBot = ocn_forcing_transmission(zBot, attenuationCoeff) + transCoeffBot = exp(attDepth) - remainingStress = remainingStress - (transmissionCoeffTop - transmissionCoeffBot) + absorb = transCoeffTop - transCoeffBot + remainingStress = remainingStress - absorb - tend(k,iEdge) = tend(k,iEdge) + edgeMask(k, iEdge) * surfaceStress(iEdge) & - * (transmissionCoeffTop - transmissionCoeffBot) / rho_sw / layerThicknessEdge(k,iEdge) + tend(k,iEdge) = tend(k,iEdge) + & + edgeMask(k,iEdge)*sfcStress(iEdge) * & + absorb/rho_sw/layerThickEdge(k,iEdge) zTop = zBot - transmissionCoeffTop = transmissionCoeffBot + transCoeffTop = transCoeffBot enddo - if ( maxLevelEdgeTop(iEdge) > 0 .and. remainingStress > 0.0_RKIND) then - tend(maxLevelEdgeTop(iEdge), iEdge) = tend(maxLevelEdgeTop(iEdge), iEdge) & - + edgeMask(maxLevelEdgeTop(iEdge), iEdge) * surfaceStress(iEdge) * remainingStress & - / rho_sw / layerThicknessEdge(maxLevelEdgeTop(iEdge), iEdge) + !*** if there is any remaining stress at the bottom, add it + !*** to the bottom layer + + if ( kmax > 0 .and. remainingStress > 0.0_RKIND) then + tend(kmax,iEdge) = tend(kmax,iEdge) & + + edgeMask(kmax,iEdge)*sfcStress(iEdge)* & + remainingStress/rho_sw/ & + layerThickEdge(kmax,iEdge) end if enddo + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif call mpas_timer_stop('vel surface stress') @@ -200,21 +210,24 @@ end subroutine ocn_vel_forcing_surface_stress_tend!}}} subroutine ocn_vel_forcing_surface_stress_init(err)!{{{ - !-------------------------------------------------------------------- - !----------------------------------------------------------------- - ! - ! call individual init routines for each parameterization - ! + ! output variables !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag - surfaceStressOn = .true. + ! End preamble + !----------------------------------------------------------------- + ! Begin code - if(config_disable_vel_surface_stress) surfaceStressOn = .false. + !*** Initialize return error code and set module defaults err = 0 + surfaceStressOff = .false. + + !*** Reset based on input model configuration + + if (config_disable_vel_surface_stress) surfaceStressOff = .true. !-------------------------------------------------------------------- diff --git a/src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F b/src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F index 049ee2ba6b..f0e3029abb 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F @@ -22,9 +22,8 @@ module ocn_vel_hadv_coriolis use mpas_timer - use mpas_derived_types - use mpas_pool_routines use ocn_constants + use ocn_mesh use ocn_config implicit none @@ -52,8 +51,13 @@ module ocn_vel_hadv_coriolis ! !-------------------------------------------------------------------- - logical :: hadvAndCoriolisOn - integer :: RK4On + logical :: & + hadvCoriolisDisabled ! on/off switch for hadv/Coriolis tend + ! use disabled since default is on + + logical :: & + usePlanetVorticity ! multiplicative mask for including + ! planetary vorticity term !*********************************************************************** @@ -66,100 +70,132 @@ module ocn_vel_hadv_coriolis !> \brief Computes tendency term for horizontal advection and coriolis force !> \author Mark Petersen, Doug Jacobsen, Todd Ringler !> \date September 2011 -!> \details - -!> This routine computes the horizontal advection and coriolis and -!> advection tendencies for momentum based on current state. - +!> \details This routine computes the horizontal advection and +!> coriolis advection tendencies for momentum based on +!> current state. ! !----------------------------------------------------------------------- - subroutine ocn_vel_hadv_coriolis_tend(meshPool, normalizedRelativeVorticityEdge, normalizedPlanetaryVorticityEdge, & - layerThicknessEdge, normalVelocity, kineticEnergyCell, tend, err)!{{{ + subroutine ocn_vel_hadv_coriolis_tend(normRelVortEdge, & + normPlanetVortEdge, & + layerThickEdge, normalVelocity, & + kineticEnergyCell, tend, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(in) :: & - normalizedRelativeVorticityEdge, &!< Input: relative vorticity over thickness, on an edge - normalizedPlanetaryVorticityEdge,&!< Input: planetary vorticity over thickness, on an edge - layerThicknessEdge,&!< Input: Thickness on edge - normalVelocity,& !< Input: Horizontal velocity - kineticEnergyCell !< Input: Kinetic Energy - - type (mpas_pool_type), intent(in) :: & - meshPool !< Input: mesh information + normRelVortEdge, &!< [in] relative vorticity/thickness at edge + normPlanetVortEdge, &!< [in] planetary vorticity/thickness at edge + layerThickEdge, &!< [in] layer thickness on edge + normalVelocity, &!< [in] Horizontal velocity + kineticEnergyCell !< [in] Kinetic Energy !----------------------------------------------------------------- - ! ! input/output variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(inout) :: & - tend !< Input/Output: velocity tendency + tend !< [inout] accumulated velocity tendency !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnEdge - integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnEdge, edgeMask - real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge - real (kind=RKIND), dimension(:), pointer :: dcEdge + integer :: & + iEdge, j, k, &! loop indices for edges, vertical + cell1, cell2, &! neighbor cell indices + eoe ! edge on edge index - integer :: j, k - integer :: cell1, cell2, iEdge, eoe, nEdges - integer, pointer :: nVertLevels - integer, dimension(:), pointer :: nEdgesArray - real (kind=RKIND) :: workVorticity, invLength, edgeWeight, r_tmp - real (kind=RKIND), dimension(:), allocatable :: qArr + real (kind=RKIND) :: & + avgVorticity, &! vorticity averaged across edge + invLength, &! temp variable for 1/dcedge + edgeWeight ! weight for summing over edges - err = 0 - - if ( .not. hadvAndCoriolisOn ) return + real (kind=RKIND), dimension(:,:), allocatable :: & + qArr, &! temporary vorticity array + tmpVorticity ! vorticity on edge - call mpas_timer_start("coriolis") + ! End preamble + !----------------------------------------------------------------- + ! Begin code - call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'nEdgesOnEdge', nEdgesOnEdge) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - call mpas_pool_get_array(meshPool, 'edgesOnEdge', edgesOnEdge) - call mpas_pool_get_array(meshPool, 'weightsOnEdge', weightsOnEdge) - call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge) - call mpas_pool_get_array(meshPool, 'edgeMask', edgeMask) + !*** Set error code and exit early if disabled + !*** Otherwise, start timer - call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) - call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + err = 0 + if (hadvCoriolisDisabled) return + call mpas_timer_start("coriolis") - nEdges = nEdgesArray( 1 ) + !*** allocate and transfer temporary arrays - allocate( qArr(nVertLevels) ) + allocate(tmpVorticity(nVertLevels,nEdgesAll), & + qArr(nVertLevels,nEdgesAll)) + !$acc enter data create(tmpVorticity, qArr) + #ifndef MPAS_OPENACC + !$omp parallel + #endif + if (usePlanetVorticity) then + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(tmpVorticity, normRelVortEdge, & + !$acc normPlanetVortEdge) & + !$acc private(k) + #else + !$omp do schedule(runtime) & + !$omp private(k) + #endif + do iEdge = 1, nEdgesAll + do k=1,nVertLevels + tmpVorticity(k,iEdge) = normRelVortEdge(k,iEdge) + & + normPlanetVortEdge(k,iEdge) + end do + end do + #ifndef MPAS_OPENACC + !$omp end do + #endif + else + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(tmpVorticity, normRelVortEdge) & + !$acc private(k) + #else + !$omp do schedule(runtime) & + !$omp private(k) + #endif + do iEdge = 1, nEdgesAll + do k=1,nVertLevels + tmpVorticity(k,iEdge) = normRelVortEdge(k,iEdge) + end do + end do + #ifndef MPAS_OPENACC + !$omp end do + #endif + endif + + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(nEdgesOnEdge, edgesOnEdge, weightsOnEdge, & + !$acc maxLevelEdgeTop, qArr, tmpVorticity, & + !$acc normalVelocity, layerThickEdge) & + !$acc private(k, j, eoe, edgeWeight, avgVorticity) + #else !$omp parallel !$omp do schedule(runtime) & - !$omp private(cell1, cell2, invLength, k, qArr, j, eoe, edgeWeight, workVorticity) - do iEdge = 1, nEdges - cell1 = cellsOnEdge(1,iEdge) - cell2 = cellsOnEdge(2,iEdge) + !$omp private(k, j, eoe, edgeWeight, avgVorticity) + #endif + do iEdge = 1, nEdgesOwned - invLength = 1.0_RKIND / dcEdge(iEdge) - - do k = 1, maxLevelEdgeTop(iEdge) - qArr(k) = 0.0_RKIND + do k = 1, nVertLevels + qArr(k,iEdge) = 0.0_RKIND end do do j = 1, nEdgesOnEdge(iEdge) @@ -167,24 +203,48 @@ subroutine ocn_vel_hadv_coriolis_tend(meshPool, normalizedRelativeVorticityEdge, edgeWeight = weightsOnEdge(j, iEdge) do k = 1, maxLevelEdgeTop(iEdge) - workVorticity = 0.5_RKIND & - * ( normalizedRelativeVorticityEdge(k, iEdge) + RK4On * normalizedPlanetaryVorticityEdge(k, iEdge) & - + normalizedRelativeVorticityEdge(k, eoe) + RK4On * normalizedPlanetaryVorticityEdge(k, eoe)) - qArr(k) = qArr(k) + edgeWeight * normalVelocity(k, eoe) * workVorticity * layerThicknessEdge(k, eoe) + avgVorticity = 0.5_RKIND* & + (tmpVorticity(k,iEdge) + & + tmpVorticity(k,eoe)) + qArr(k,iEdge) = qArr(k,iEdge) + & + edgeWeight*normalVelocity(k,eoe)* & + avgVorticity*layerThickEdge(k,eoe) end do end do + end do + #ifndef MPAS_OPENACC + !$omp end do + #endif + + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, maxLevelEdgeTop, edgeMask, & + !$acc dcEdge, tend, qArr, kineticEnergyCell) & + !$acc private(cell1, cell2, invLength, k) + #else + !$omp do schedule(runtime) & + !$omp private(cell1, cell2, invLength, k) + #endif + do iEdge = 1, nEdgesOwned + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + invLength = 1.0_RKIND / dcEdge(iEdge) do k = 1, maxLevelEdgeTop(iEdge) - tend(k, iEdge) = tend(k, iEdge) + edgeMask(k, iEdge) * ( qArr(k) - ( kineticEnergyCell(k, cell2) & - - kineticEnergyCell(k, cell1) ) * invLength ) + tend(k,iEdge) = tend(k,iEdge) + & + edgeMask(k,iEdge)* (qArr(k,iEdge) - & + (kineticEnergyCell(k,cell2) & + - kineticEnergyCell(k,cell1))*invLength) end do - end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif - deallocate( qArr ) + !$acc exit data delete(tmpVorticity, qArr) + deallocate(qArr,tmpVorticity) call mpas_timer_stop("coriolis") @@ -197,7 +257,7 @@ end subroutine ocn_vel_hadv_coriolis_tend!}}} ! routine ocn_vel_hadv_coriolis_init ! !> \brief Initializes ocean momentum horizontal advection and -!> coriolis tendencies +!> coriolis tendencies !> \author Mark Petersen, Doug Jacobsen, Todd Ringler !> \date September 2011 ! @@ -208,30 +268,47 @@ subroutine ocn_vel_hadv_coriolis_init(err)!{{{ !-------------------------------------------------------------------- !----------------------------------------------------------------- - ! ! Output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + !*** Initialize returned error code and default flags err = 0 - hadvAndCoriolisOn = .true. + hadvCoriolisDisabled = .false. + usePlanetVorticity = .false. - if ( config_disable_vel_coriolis ) hadvAndCoriolisOn = .false. + !*** Reset module variables based in input configuration - if ( trim( config_time_integrator ) == 'RK4') then + if ( config_disable_vel_coriolis ) hadvCoriolisDisabled = .true. + + select case (trim(config_time_integrator)) + case ('RK4','rk4') ! For RK4, coriolis tendency term includes f: (eta+f)/h. - RK4On = 1 - elseif ( trim( config_time_integrator ) == 'split_explicit' & - .or. trim( config_time_integrator ) == 'unsplit_explicit' & - .or. trim( config_time_integrator ) == 'semi_implicit') then - ! For split explicit, Coriolis tendency uses eta/h because the Coriolis term - ! is added separately to the momentum tendencies. - RK4On = 0 - end if + usePlanetVorticity = .true. + + case ('split_explicit') + ! For split explicit, Coriolis tendency uses eta/h because + ! the Coriolis term is added separately to the momentum + ! tendencies. + usePlanetVorticity = .false. + + case ('unsplit_explicit') + ! For split explicit, Coriolis tendency uses eta/h because + ! the Coriolis term is added separately to the momentum + ! tendencies. + usePlanetVorticity = .false. + + case ('semi_implicit') + usePlanetaryVorticity = .false. + end select !-------------------------------------------------------------------- diff --git a/src/core_ocean/shared/mpas_ocn_vel_hmix.F b/src/core_ocean/shared/mpas_ocn_vel_hmix.F index 039e330128..0333fe6680 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hmix.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hmix.F @@ -24,8 +24,7 @@ module ocn_vel_hmix use mpas_timer use mpas_derived_types - use mpas_pool_routines - use mpas_threading + use mpas_log use ocn_vel_hmix_del2 use ocn_vel_hmix_leith use ocn_vel_hmix_del4 @@ -57,7 +56,8 @@ module ocn_vel_hmix ! !-------------------------------------------------------------------- - logical :: hmixOn + logical :: & + hmixOff ! on/off switch for horizontal velocity mixing !*********************************************************************** @@ -80,82 +80,58 @@ module ocn_vel_hmix ! !----------------------------------------------------------------------- - subroutine ocn_vel_hmix_tend(meshPool, scratchPool, divergence, relativeVorticity, normalVelocity, tangentialVelocity, & - viscosity, tend, err)!{{{ + subroutine ocn_vel_hmix_tend(div, relVort, tend, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- - type (mpas_pool_type), intent(in) :: & - meshPool !< Input: mesh information - - type (mpas_pool_type), intent(inout) :: scratchPool !< Input: scratch variables - - real (kind=RKIND), dimension(:,:), intent(in) :: & - divergence !< Input: velocity divergence - - real (kind=RKIND), dimension(:,:), intent(in) :: & - relativeVorticity !< Input: relative vorticity - - real (kind=RKIND), dimension(:,:), intent(in) :: & - normalVelocity !< Input: velocity normal to an edge - real (kind=RKIND), dimension(:,:), intent(in) :: & - tangentialVelocity !< Input: velocity, tangent to an edge + div, &!< [in] velocity divergence + relVort !< [in] relative vorticity !----------------------------------------------------------------- - ! ! input/output variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(inout) :: & - viscosity !< Input: viscosity - - real (kind=RKIND), dimension(:,:), intent(inout) :: & - tend !< Input/Output: velocity tendency + tend !< [inout] accumulated velocity tendency !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer :: err1 + integer :: & + err1 ! local error flag + ! End preamble !----------------------------------------------------------------- - ! - ! call relevant routines for computing tendencies - ! note that the user can choose multiple options and the - ! tendencies will be added together - ! - !----------------------------------------------------------------- + ! Begin code - if(.not.hmixOn) return + !*** Initialize return error code and viscosity + !*** exit early if not turned on, otherwise start relevant timer + err = 0 + if (hmixOff) return call mpas_timer_start("vel hmix") - viscosity = 0.0_RKIND - err = 0 + ! call relevant routines for computing tendencies + ! note that the user can choose multiple options and the + ! tendencies will be added together - call ocn_vel_hmix_del2_tend(meshPool, divergence, relativeVorticity, viscosity, tend, err1) + call ocn_vel_hmix_del2_tend(div, relVort, tend, err1) err = ior(err1, err) - call ocn_vel_hmix_leith_tend(meshPool, divergence, relativeVorticity, viscosity, tend, err1) + call ocn_vel_hmix_leith_tend(div, relVort, tend, err1) err = ior(err1, err) - call ocn_vel_hmix_del4_tend(meshPool, divergence, relativeVorticity, tend, err1) + call ocn_vel_hmix_del4_tend(div, relVort, tend, err1) err = ior(err1, err) call mpas_timer_stop("vel hmix") @@ -181,27 +157,43 @@ end subroutine ocn_vel_hmix_tend!}}} subroutine ocn_vel_hmix_init(err)!{{{ - !-------------------------------------------------------------------- + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! - ! call individual init routines for each parameterization - ! + ! local variables !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer :: err1, err2, err3 ! local error flags - integer :: err1, err2, err3 + ! End preamble + !----------------------------------------------------------------- + ! Begin code - hmixOn = .true. + !*** set error code and defaults + err = 0 + hmixOff = .false. + + !*** call individual init routines call ocn_vel_hmix_del2_init(err1) + if (err1 /= 0) call mpas_log_write( & + 'Error encountered initializing vel_hmix_del2', MPAS_LOG_ERR) + call ocn_vel_hmix_leith_init(err2) + if (err2 /= 0) call mpas_log_write( & + 'Error encountered initializing vel_hmix_leith', MPAS_LOG_ERR) + call ocn_vel_hmix_del4_init(err3) + if (err3 /= 0) call mpas_log_write( & + 'Error encountered initializing vel_hmix_del4', MPAS_LOG_ERR) err = ior(ior(err1, err2),err3) - if(config_disable_vel_hmix) hmixOn = .false. + if (config_disable_vel_hmix) hmixOff = .true. !-------------------------------------------------------------------- diff --git a/src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F b/src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F index b3148813b7..a5c5adbc16 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F @@ -21,13 +21,9 @@ module ocn_vel_hmix_del2 use mpas_timer - use mpas_derived_types - use mpas_pool_routines - use mpas_threading - use mpas_vector_operations - use mpas_matrix_operations use ocn_constants use ocn_config + use ocn_mesh implicit none private @@ -54,8 +50,12 @@ module ocn_vel_hmix_del2 ! !-------------------------------------------------------------------- - logical :: hmixDel2On !< integer flag to determine whether del2 chosen + logical :: & + hmixDel2Off !< on/off flag to determine whether del2 chosen + real (kind=RKIND) :: & + viscDel2 !< viscosity coefficient for del2 horz mixing + !*********************************************************************** contains @@ -64,70 +64,56 @@ module ocn_vel_hmix_del2 ! ! routine ocn_vel_hmix_del2_tend ! -!> \brief Computes tendency term for Laplacian horizontal momentum mixing +!> \brief Computes tendency for Laplacian horizontal momentum mixing !> \author Mark Petersen, Doug Jacobsen, Todd Ringler !> \date 22 August 2011 !> \details !> This routine computes the horizontal mixing tendency for momentum !> based on a Laplacian form for the mixing, \f$\nu_2 \nabla^2 u\f$ -!> This tendency takes the -!> form \f$\nu( \nabla divergence + k \times \nabla relativeVorticity )\f$, -!> where \f$\nu\f$ is a viscosity and \f$k\f$ is the vertical unit vector. -!> This form is strictly only valid for constant \f$\nu\f$ . +!> This tendency takes the form +!> \f$\nu( \nabla divergence + k \times \nabla relativeVorticity )\f$, +!> where \f$\nu\f$ is a viscosity and \f$k\f$ is the vertical unit +!> vector. This form is strictly only valid for constant \f$\nu\f$ . ! !----------------------------------------------------------------------- - subroutine ocn_vel_hmix_del2_tend(meshPool, divergence, relativeVorticity, viscosity, tend, err)!{{{ + subroutine ocn_vel_hmix_del2_tend(div, relVort, tend, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(in) :: & - divergence !< Input: velocity divergence - - real (kind=RKIND), dimension(:,:), intent(in) :: & - relativeVorticity !< Input: relative vorticity - - type (mpas_pool_type), intent(in) :: & - meshPool !< Input: mesh information + div, &!< [in] velocity divergence + relVort !< [in] relative vorticity - !------ ----------------------------------------------------------- - ! + !----------------------------------------------------------------- ! input /output variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(inout) :: & - viscosity !< Input: viscosity - - real (kind=RKIND), dimension(:,:), intent(inout) :: & - tend !< Input/Output: velocity tendency + tend !< [inout] accumulated velocity tendency !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer :: iEdge, cell1, cell2, vertex1, vertex2, k, nEdges - integer, dimension(:), pointer :: nEdgesArray - integer, dimension(:), pointer :: maxLevelEdgeTop - integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask + integer :: & + iEdge, k, &! loop indices for edge, vertical loops + cell1, cell2, &! neighbor cell addresses across edge + vertex1, vertex2 ! neighbor vertex addresses along edge - real (kind=RKIND) :: u_diffusion, invLength1, invLength2, visc2 - real (kind=RKIND), dimension(:), pointer :: meshScalingDel2, & - dcEdge, dvEdge + real (kind=RKIND) :: & + uDiff, &! velocity diffusion + dcEdgeInv, &! 1/dcEdge + dvEdgeInv, &! 1/dvEdge + visc2 ! scaled viscosity coefficient !----------------------------------------------------------------- ! @@ -137,52 +123,55 @@ subroutine ocn_vel_hmix_del2_tend(meshPool, divergence, relativeVorticity, visco err = 0 - if(.not.hmixDel2On) return + if (hmixDel2Off) return call mpas_timer_start("vel del2") - call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) - - call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge) - call mpas_pool_get_array(meshPool, 'meshScalingDel2', meshScalingDel2) - call mpas_pool_get_array(meshPool, 'edgeMask', edgeMask) - call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge) - call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) - - nEdges = nEdgesArray( 1 ) - + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, maxLevelEdgeTop, verticesOnEdge, & + !$acc edgeMask, dcEdge, dvEdge, meshScalingDel2, & + !$acc div, relVort, tend) & + !$acc private(k, cell1, cell2, uDiff, vertex1, vertex2, & + !$acc dcEdgeInv, dvEdgeInv, visc2) + #else !$omp parallel - !$omp do schedule(runtime) private(cell1, cell2, vertex1, vertex2, invLength1, invLength2, k, u_diffusion, visc2) - do iEdge = 1, nEdges + !$omp do schedule(runtime) & + !$omp private(k, cell1, cell2, uDiff, vertex1, vertex2, & + !$omp dcEdgeInv, dvEdgeInv, visc2) + #endif + do iEdge = 1, nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) vertex1 = verticesOnEdge(1,iEdge) vertex2 = verticesOnEdge(2,iEdge) - invLength1 = 1.0_RKIND / dcEdge(iEdge) - invLength2 = 1.0_RKIND / dvEdge(iEdge) - - do k = 1, maxLevelEdgeTop(iEdge) + dcEdgeInv = 1.0_RKIND / dcEdge(iEdge) + dvEdgeInv = 1.0_RKIND / dvEdge(iEdge) - ! Here -( relativeVorticity(k,vertex2) - relativeVorticity(k,vertex1) ) / dvEdge(iEdge) - ! is - \nabla relativeVorticity pointing from vertex 2 to vertex 1, or equivalently - ! + k \times \nabla relativeVorticity pointing from cell1 to cell2. + visc2 = viscDel2*meshScalingDel2(iEdge) - u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) * invLength1 & - -( relativeVorticity(k,vertex2) - relativeVorticity(k,vertex1) ) * invLength2 + do k = 1, maxLevelEdgeTop(iEdge) - visc2 = config_mom_del2 * meshScalingDel2(iEdge) + ! Here -( relativeVorticity(k,vertex2) - + ! relativeVorticity(k,vertex1) ) / dvEdge(iEdge) + ! is - \nabla relativeVorticity pointing from vertex 2 + ! to vertex 1, or equivalently + ! + k \times \nabla relativeVorticity pointing from cell1 + ! to cell2. - tend(k,iEdge) = tend(k,iEdge) + edgeMask(k, iEdge) * visc2 * u_diffusion + uDiff = (div(k,cell2) - div(k,cell1))*dcEdgeInv & + -(relVort(k,vertex2)-relVort(k,vertex1))*dvEdgeInv - viscosity(k,iEdge) = viscosity(k,iEdge) + visc2 + tend(k,iEdge) = tend(k,iEdge) + & + edgeMask(k,iEdge)*visc2*uDiff end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif call mpas_timer_stop("vel del2") @@ -205,27 +194,36 @@ end subroutine ocn_vel_hmix_del2_tend!}}} subroutine ocn_vel_hmix_del2_init(err)!{{{ + !----------------------------------------------------------------- + ! Output variables + !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag - - !-------------------------------------------------------------------- - ! - ! set some local module variables based on input config choices - ! - !-------------------------------------------------------------------- + integer, intent(out) :: err !< [out] error flag - err = 0 + ! End preamble + !----------------------------------------------------------------- + ! Begin code - hmixDel2On = .false. + !*** Initialize return code and default variables + err = 0 + hmixDel2Off = .true. + viscDel2 = 0.0_RKIND - if ( config_mom_del2 > 0.0_RKIND ) then - hmixDel2On = .true. - endif + !*** reset values based on input configuration - if ( .not. config_use_mom_del2 ) hmixDel2On = .false. + if (config_use_mom_del2) then + hmixDel2Off = .false. + viscDel2 = config_mom_del2 + if (viscDel2 <= 0.0_RKIND) then + call mpas_log_write( & + 'vel_hmix_del2_init: viscosity must be > 0', & + MPAS_LOG_ERR) + err = -1 + endif + endif - !-------------------------------------------------------------------- + !----------------------------------------------------------------- end subroutine ocn_vel_hmix_del2_init!}}} diff --git a/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F b/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F index 1e3dc8a12b..f9838d98a7 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F @@ -21,13 +21,10 @@ module ocn_vel_hmix_del4 use mpas_timer - use mpas_derived_types - use mpas_pool_routines - use mpas_threading - use mpas_vector_operations - use mpas_matrix_operations + use mpas_log use ocn_constants use ocn_config + use ocn_mesh implicit none private @@ -54,7 +51,12 @@ module ocn_vel_hmix_del4 ! !-------------------------------------------------------------------- - logical :: hmixDel4On !< local flag to determine whether del4 chosen + logical :: & + hmixDel4Off !< on/off flag to determine whether del4 chosen + + real (kind=RKIND) :: & + viscDel4, &!< biharmonic viscosity coefficient + divFactor !< factor for divergence term !*********************************************************************** @@ -64,7 +66,7 @@ module ocn_vel_hmix_del4 ! ! routine ocn_vel_hmix_del4_tend ! -!> \brief Computes tendency term for biharmonic horizontal momentum mixing +!> \brief Computes tendency for biharmonic horizontal momentum mixing !> \author Mark Petersen, Doug Jacobsen, Todd Ringler !> \date September 2011 !> \details @@ -78,197 +80,214 @@ module ocn_vel_hmix_del4 ! !----------------------------------------------------------------------- - subroutine ocn_vel_hmix_del4_tend(meshPool, divergence, relativeVorticity, tend, err)!{{{ + subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(in) :: & - divergence !< Input: velocity divergence - - real (kind=RKIND), dimension(:,:), intent(in) :: & - relativeVorticity !< Input: relative vorticity - - type (mpas_pool_type), intent(in) :: & - meshPool !< Input: mesh information + div, &!< [in] velocity divergence + relVort !< [in] relative vorticity !----------------------------------------------------------------- - ! ! input/output variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(inout) :: & - tend !< Input/Output: velocity tendency + tend !< [inout] accumulated velocity tendency !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer :: iEdge, cell1, cell2, vertex1, vertex2, k, i - integer :: iCell, iVertex, nEdges, nCells, nVertices - integer, pointer :: nVertLevels, vertexDegree - integer, dimension(:), pointer :: nEdgesArray, nCellsArray, nVerticesArray + integer :: & + iEdge, iCell, iVertex, k, i, &! loop counters + cell1, cell2, &! neighbor cell addresses across edge + vertex1, vertex2, &! neighbor vertex addresses along edge + nEdges, nCells, nVertices ! num edges,cells,vertices incl halos - integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelVertexTop, & - maxLevelCell, nEdgesOnCell - integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask, edgesOnVertex, edgesOnCell, edgeSignOnVertex, & - edgeSignOnCell + real (kind=RKIND) :: & + uDiff, &! diffusion operator temporary + areaCellInv, &! 1/area of cell + areaTriInv, &! 1/area of triangle + dcEdgeInv, &! 1/dcEdge + dvEdgeInv, &! 1/dvEdge + visc4 ! scaled biharmonic viscosity coeff + ! Scratch Arrays + real (kind=RKIND), dimension(:,:), allocatable :: & + del2div, &! + del2relVort, &! + del2u - real (kind=RKIND) :: u_diffusion, invAreaCell1, invAreaCell2, invAreaTri1, & - invAreaTri2, invDcEdge, invDvEdge, r_tmp - real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaTriangle, & - meshScalingDel4, areaCell + ! End preamble + !----------------------------------------------------------------- + ! Begin code - ! Scratch Arrays - real (kind=RKIND), dimension(:,:), allocatable :: delsq_divergence - real (kind=RKIND), dimension(:,:), allocatable :: delsq_relativeVorticity - real (kind=RKIND), dimension(:,:), allocatable :: delsq_u + !*** initialize error code and return if not selected + !*** otherwise start timer err = 0 + if (hmixDel4Off) return + call mpas_timer_start("vel del4") - if(.not.hmixDel4On) return + !*** allocate temporaries - call mpas_timer_start("vel del4") + allocate(del2u(nVertLevels, nEdgesAll + 1), & + del2div(nVertLevels, nCellsAll + 1), & + del2relVort(nVertLevels, nVerticesAll + 1)) + !$acc enter data create(del2u, del2div, del2relVort) + + nEdges = nEdgesHalo(2) - call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) - call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) - call mpas_pool_get_dimension(meshPool, 'nVerticesArray', nVerticesArray) - call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) - call mpas_pool_get_dimension(meshPool, 'vertexDegree', vertexDegree) - - call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'maxLevelVertexTop', maxLevelVertexTop) - call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge) - call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge) - call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) - call mpas_pool_get_array(meshPool, 'areaTriangle', areaTriangle) - call mpas_pool_get_array(meshPool, 'areaCell', areaCell) - call mpas_pool_get_array(meshPool, 'meshScalingDel4', meshScalingDel4) - call mpas_pool_get_array(meshPool, 'edgeMask', edgeMask) - call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) - call mpas_pool_get_array(meshPool, 'edgesOnVertex', edgesOnVertex) - call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) - call mpas_pool_get_array(meshPool, 'edgeSignOnVertex', edgeSignOnVertex) - call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell) - - nEdges = nEdgesArray(size(nEdgesArray)) - nCells = nCellsArray(size(nCellsArray)) - nVertices = nVerticesArray(size(nVerticesArray)) - allocate(delsq_u(nVertLevels, nEdges + 1), & - delsq_divergence(nVertLevels, nCells + 1), & - delsq_relativeVorticity(nVertLevels, nVertices + 1)) - - nEdges = nEdgesArray( 3 ) - - !Compute delsq_u + !Compute Laplacian of velocity + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, verticesOnEdge, maxLevelEdgeTop, & + !$acc dvEdge, dcEdge, del2u, div, relVort) & + !$acc private(k, cell1, cell2, vertex1, vertex2, & + !$acc dcEdgeInv, dvEdgeInv) + #else !$omp parallel !$omp do schedule(runtime) & - !$omp private(cell1, cell2, vertex1, vertex2, invDcEdge, invDvEdge, k) + !$omp private(k, cell1, cell2, vertex1, vertex2, & + !$omp dcEdgeInv, dvEdgeInv) + #endif do iEdge = 1, nEdges - delsq_u(:, iEdge) = 0.0_RKIND + del2u(:, iEdge) = 0.0_RKIND cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) vertex1 = verticesOnEdge(1,iEdge) vertex2 = verticesOnEdge(2,iEdge) - invDcEdge = 1.0_RKIND / dcEdge(iEdge) - invDvEdge = 1.0_RKIND / max(dvEdge(iEdge), 0.25_RKIND*dcEdge(iEdge)) + dcEdgeInv = 1.0_RKIND / dcEdge(iEdge) + dvEdgeInv = 1.0_RKIND / max(dvEdge(iEdge), 0.25_RKIND*dcEdge(iEdge)) do k=1,maxLevelEdgeTop(iEdge) ! Compute \nabla^2 u = \nabla divergence + k \times \nabla relativeVorticity - delsq_u(k, iEdge) = ( divergence(k,cell2) - divergence(k,cell1) ) * invDcEdge & - -( relativeVorticity(k,vertex2) - relativeVorticity(k,vertex1)) * invDvEdge + del2u(k,iEdge) = (div(k,cell2) - div(k,cell1))*dcEdgeInv & + -(relVort(k,vertex2) - relVort(k,vertex1))*& + dvEdgeInv end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel - - nVertices = nVerticesArray( 2 ) - - ! Compute delsq_relativeVorticity + #endif + + nVertices = nVerticesHalo(1) + + ! Compute del2 of relative vorticity + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(edgesOnVertex, maxLevelVertexTop, dcEdge, & + !$acc areaTriangle, edgeSignOnVertex, del2u, & + !$acc del2relVort) & + !$acc private(i, k, iEdge, areaTriInv) + #else !$omp parallel - !$omp do schedule(runtime) private(invAreaTri1, i, iEdge, k) + !$omp do schedule(runtime) private(i, k, iEdge, areaTriInv) + #endif do iVertex = 1, nVertices - delsq_relativeVorticity(:, iVertex) = 0.0_RKIND - invAreaTri1 = 1.0_RKIND / areaTriangle(iVertex) + del2relVort(:, iVertex) = 0.0_RKIND + areaTriInv = 1.0_RKIND/areaTriangle(iVertex) do i = 1, vertexDegree iEdge = edgesOnVertex(i, iVertex) do k = 1, maxLevelVertexTop(iVertex) - delsq_relativeVorticity(k, iVertex) = delsq_relativeVorticity(k, iVertex) + edgeSignOnVertex(i, iVertex) & - * dcEdge(iEdge) * delsq_u(k, iEdge) * invAreaTri1 + del2relVort(k,iVertex) = del2relVort(k,iVertex) & + + edgeSignOnVertex(i,iVertex) & + *dcEdge(iEdge)*del2u(k,iEdge)*areaTriInv end do end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif - nCells = nCellsArray( 2 ) + nCells = nCellsHalo(1) - ! Compute delsq_divergence + ! Compute del2 of divergence + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(nEdgesOnCell, edgesOnCell, maxLevelCell, dvEdge,& + !$acc edgeSignOnCell, areaCell, del2u, del2div) & + !$acc private(i, k, iEdge, areaCellInv) + #else !$omp parallel - !$omp do schedule(runtime) private(invAreaCell1, i, iEdge, k) + !$omp do schedule(runtime) private(i, k, iEdge, areaCellInv) + #endif do iCell = 1, nCells - delsq_divergence(:, iCell) = 0.0_RKIND - invAreaCell1 = 1.0_RKIND / areaCell(iCell) + del2div(:, iCell) = 0.0_RKIND + areaCellInv = 1.0_RKIND / areaCell(iCell) do i = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(i, iCell) do k = 1, maxLevelCell(iCell) - delsq_divergence(k, iCell) = delsq_divergence(k, iCell) - edgeSignOnCell(i, iCell) * dvEdge(iEdge) & - * delsq_u(k, iEdge) * invAreaCell1 + del2div(k,iCell) = del2div(k,iCell) - & + edgeSignOnCell(i,iCell)*dvEdge(iEdge) & + *del2u(k,iEdge)*areaCellInv end do end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel - - nEdges = nEdgesArray( 1 ) - - ! Compute - \kappa \nabla^4 u - ! as \nabla div(\nabla^2 u) + k \times \nabla ( k \cross curl(\nabla^2 u) ) + #endif + + ! Compute final tendency \kappa \nabla^4 u + ! as \nabla div(\nabla^2 u) + k \times + ! \nabla ( k \cross curl(\nabla^2 u) ) + + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, verticesOnEdge, maxLevelEdgeTop, & + !$acc dcEdge, dvEdge, meshScalingDel4, edgeMask, & + !$acc del2div, del2relVort, tend) & + !$acc private(k, cell1, cell2, vertex1, vertex2, & + !$acc dcEdgeInv, dvEdgeInv, visc4, uDiff) + #else !$omp parallel !$omp do schedule(runtime) & - !$omp private(cell1, cell2, vertex1, vertex2, invDcEdge, invDvEdge, r_tmp, k, u_diffusion) - do iEdge = 1, nEdges + !$omp private(k, cell1, cell2, vertex1, vertex2, & + !$omp dcEdgeInv, dvEdgeInv, visc4, uDiff) + #endif + do iEdge = 1, nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) vertex1 = verticesOnEdge(1,iEdge) vertex2 = verticesOnEdge(2,iEdge) - invDcEdge = 1.0_RKIND / dcEdge(iEdge) - invDvEdge = 1.0_RKIND / dvEdge(iEdge) - r_tmp = config_mom_del4 * meshScalingDel4(iEdge) + dcEdgeInv = 1.0_RKIND / dcEdge(iEdge) + dvEdgeInv = 1.0_RKIND / dvEdge(iEdge) + visc4 = viscDel4*meshScalingDel4(iEdge) - do k=1,maxLevelEdgeTop(iEdge) - u_diffusion = config_mom_del4_div_factor*(delsq_divergence(k,cell2) - delsq_divergence(k,cell1)) * invDcEdge & - - (delsq_relativeVorticity(k,vertex2) - delsq_relativeVorticity(k,vertex1) ) * invDvEdge + do k=1, maxLevelEdgeTop(iEdge) + uDiff = divFactor*(del2div(k,cell2) - del2div(k,cell1))* & + dcEdgeInv & + - (del2relVort(k,vertex2) - del2relVort(k,vertex1))* & + dvEdgeInv - tend(k,iEdge) = tend(k,iEdge) - edgeMask(k, iEdge) * u_diffusion * r_tmp + tend(k,iEdge) = tend(k,iEdge) - & + edgeMask(k,iEdge)*uDiff*visc4 end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif - deallocate(delsq_u, & - delsq_divergence, & - delsq_relativeVorticity) + !$acc exit data delete(del2u, del2div, del2relVort) + deallocate(del2u, & + del2div, & + del2relVort) call mpas_timer_stop("vel del4") @@ -291,23 +310,36 @@ end subroutine ocn_vel_hmix_del4_tend!}}} subroutine ocn_vel_hmix_del4_init(err)!{{{ - integer, intent(out) :: err !< Output: error flag - - !-------------------------------------------------------------------- - ! - ! set some local module variables based on input config choices - ! - !-------------------------------------------------------------------- + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- - err = 0 + integer, intent(out) :: err !< [out] error flag - hmixDel4On = .false. + ! End preamble + !----------------------------------------------------------------- + ! Begin code - if ( config_mom_del4 > 0.0_RKIND ) then - hmixDel4On = .true. - endif + !*** initialize return error code and set module variable defaults - if(.not.config_use_mom_del4) hmixDel4On = .false. + err = 0 + hmixDel4Off = .true. + viscDel4 = 0.0_RKIND + divFactor = 1.0_RKIND + + !*** reset module variables based on input configuration + + if (config_use_mom_del4) then + hmixDel4Off = .false. + viscDel4 = config_mom_del4 + divFactor = config_mom_del4_div_factor + if ( config_mom_del4 <= 0.0_RKIND ) then + call mpas_log_write( & + 'vel_hmix_del4_init: config_mom_del4 must be > 0 ', & + MPAS_LOG_ERR) + err = -1 + endif + endif !-------------------------------------------------------------------- diff --git a/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F b/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F index 04e614446b..0d37a839cf 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F @@ -21,11 +21,10 @@ module ocn_vel_hmix_leith use mpas_timer - use mpas_derived_types - use mpas_pool_routines use mpas_constants use ocn_constants use ocn_config + use ocn_mesh implicit none private @@ -52,7 +51,14 @@ module ocn_vel_hmix_leith ! !-------------------------------------------------------------------- - logical :: hmixLeithOn !< integer flag to determine whether leith chosen + logical :: & + hmixLeithOff !< on/off switch to determine whether leith chosen + + real (kind=RKIND) :: & + leithParam, &!< Leith parameter + dxLeith, &!< Leith length scale + viscMaxLeith, &!< maximum viscosity + sqrt3fact !< sqrt(3) !*********************************************************************** @@ -62,137 +68,128 @@ module ocn_vel_hmix_leith ! ! routine ocn_vel_hmix_leith_tend ! -!> \brief Computes tendency term for horizontal momentum mixing with Leith parameterization +!> \brief Computes horizontal momentum mixing with Leith formulation !> \author Mark Petersen, Todd Ringler !> \date 22 October 2012 !> \details !> This routine computes the horizontal mixing tendency for momentum !> based on the Leith closure. The Leith closure is the !> enstrophy-cascade analogy to the Smagorinsky (1963) energy-cascade -!> closure, i.e. Leith (1996) assumes an inertial range of enstrophy flux -!> moving toward the mesh scale. The assumption of an enstrophy cascade -!> and dimensional analysis produces right-hand-side dissipation, -!> $\bf{D}$, of velocity of the form +!> closure, i.e. Leith (1996) assumes an inertial range of enstrophy +!> flux moving toward the mesh scale. The assumption of an enstrophy +!> cascade and dimensional analysis produces right-hand-side +!> dissipation, $\bf{D}$, of velocity of the form !> $ {\bf D} = \nabla \cdot \left( \nu_\ast \nabla {\bf u} \right) !> = \nabla \cdot \left( \gamma \left| \nabla \omega \right| !> \left( \Delta x \right)^3 \nabla \bf{u} \right) -!> where $\omega$ is the relative vorticity and $\gamma$ is a non-dimensional, -!> $O(1)$ parameter. We set $\gamma=1$. - +!> where $\omega$ is the relative vorticity and $\gamma$ is a +!> non-dimensional, $O(1)$ parameter. We set $\gamma=1$. ! !----------------------------------------------------------------------- - subroutine ocn_vel_hmix_leith_tend(meshPool, divergence, relativeVorticity, viscosity, tend, err)!{{{ + subroutine ocn_vel_hmix_leith_tend(div, relVort, tend, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(in) :: & - divergence !< Input: velocity divergence - - real (kind=RKIND), dimension(:,:), intent(in) :: & - relativeVorticity !< Input: relative vorticity - - type (mpas_pool_type), intent(in) :: & - meshPool !< Input: mesh information + div, &!< [in] velocity divergence + relVort !< [in] relative vorticity !----------------------------------------------------------------- - ! ! input/output variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(inout) :: & - tend !< Input/Output: velocity tendency - - real (kind=RKIND), dimension(:,:), intent(inout) :: & - viscosity !< Input: viscosity + tend !< [inout] accumulated velocity tendency !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer :: iEdge, cell1, cell2, vertex1, vertex2, k, nEdges - integer, dimension(:), pointer :: nEdgesArray - integer, dimension(:), pointer :: maxLevelEdgeTop - integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask + integer :: & + iEdge, k, &! loop indices for edge, vertical loops + cell1, cell2, &! neighbor cell addresses across edge + vertex1, vertex2 ! neighbor vertex addresses along edge - real (kind=RKIND) :: u_diffusion, invLength_dc, invLength_dv, visc2 - real (kind=RKIND), dimension(:), pointer :: meshScalingDel2, & - dcEdge, dvEdge + real (kind=RKIND) :: &! + uDiff, &! velocity diffusion operator + dcEdgeInv, &! 1/dcEdge + dvEdgeInv, &! 1/dvEdge + visc2tmp, &! common factor for visc2 + visc2 ! scaled viscosity coeff + ! End preamble !----------------------------------------------------------------- - ! - ! exit if this mixing is not selected - ! - !----------------------------------------------------------------- - - err = 0 + ! Begin code - if(.not.hmixLeithOn) return + !*** initialize return code and exit if Leith not chosen + !*** start timer if chosen + err = 0 + if (hmixLeithOff) return call mpas_timer_start("vel leith") - call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) - - call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge) - call mpas_pool_get_array(meshPool, 'meshScalingDel2', meshScalingDel2) - call mpas_pool_get_array(meshPool, 'edgeMask', edgeMask) - call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge) - call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) - - nEdges = nEdgesArray( 1 ) - + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, verticesOnEdge, maxLevelEdgeTop, & + !$acc dcEdge, dvEdge, meshScalingDel2, div, relVort, & + !$acc tend, edgeMask) & + !$acc private(k, cell1, cell2, vertex1, vertex2, dcEdgeInv, & + !$acc dvEdgeInv, uDiff, visc2, visc2tmp) + #else !$omp parallel - !$omp do schedule(runtime) private(cell1, cell2, vertex1, vertex2, invLength_dc, invLength_dv, k, u_diffusion, visc2) - do iEdge = 1, nEdges + !$omp do schedule(runtime) & + !$omp private(k, cell1, cell2, vertex1, vertex2, dcEdgeInv, & + !$omp dvEdgeInv, uDiff, visc2, visc2tmp) + #endif + do iEdge = 1, nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) vertex1 = verticesOnEdge(1,iEdge) vertex2 = verticesOnEdge(2,iEdge) - invLength_dc = 1.0_RKIND / dcEdge(iEdge) - invLength_dv = 1.0_RKIND / dvEdge(iEdge) + dcEdgeInv = 1.0_RKIND / dcEdge(iEdge) + dvEdgeInv = 1.0_RKIND / dvEdge(iEdge) + + visc2tmp = (leithParam*dxLeith*meshScalingDel2(iEdge)/pii)**3 do k = 1, maxLevelEdgeTop(iEdge) - ! Here -( relativeVorticity(k,vertex2) - relativeVorticity(k,vertex1) ) / dvEdge(iEdge) - ! is - \nabla relativeVorticity pointing from vertex 2 to vertex 1, or equivalently - ! + k \times \nabla relativeVorticity pointing from cell1 to cell2. + ! Here -( relativeVorticity(k,vertex2) - + ! relativeVorticity(k,vertex1) ) / dvEdge(iEdge) + ! is - \nabla relativeVorticity pointing from vertex 2 to + ! vertex 1, or equivalently + ! + k \times \nabla relativeVorticity pointing from cell1 + ! to cell2. - u_diffusion = ( divergence(k,cell2) - divergence(k,cell1) ) * invLength_dc & - -( relativeVorticity(k,vertex2) - relativeVorticity(k,vertex1) ) * invLength_dv + uDiff = (div(k,cell2) - div(k,cell1))*dcEdgeInv & + -(relVort(k,vertex2) - relVort(k,vertex1))*dvEdgeInv ! Here the first line is (\delta x)^3 ! the second line is |\nabla \omega| ! and u_diffusion is \nabla^2 u (see formula for $\bf{D}$ above). - visc2 = ( config_leith_parameter * config_leith_dx * meshScalingDel2(iEdge) / pii)**3 & - * abs( relativeVorticity(k,vertex2) - relativeVorticity(k,vertex1) ) * invLength_dc * sqrt(3.0_RKIND) - visc2 = min(visc2, config_leith_visc2_max) - - tend(k,iEdge) = tend(k,iEdge) + edgeMask(k, iEdge) * visc2 * u_diffusion + visc2 = visc2tmp & + *abs(relVort(k,vertex2) - relVort(k,vertex1))* & + dcEdgeInv*sqrt3fact + visc2 = min(visc2, viscMaxLeith) - viscosity(k,iEdge) = viscosity(k,iEdge) + visc2 + tend(k,iEdge) = tend(k,iEdge) + & + edgeMask(k,iEdge)*visc2*uDiff end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif call mpas_timer_stop("vel leith") @@ -204,7 +201,7 @@ end subroutine ocn_vel_hmix_leith_tend!}}} ! ! routine ocn_vel_hmix_leith_init ! -!> \brief Initializes ocean momentum horizontal mixing with Leith parameterization +!> \brief Initializes ocean momentum horiz mixing in Leith formulation !> \author Mark Petersen !> \date 22 October 2012 !> \details @@ -215,24 +212,35 @@ end subroutine ocn_vel_hmix_leith_tend!}}} subroutine ocn_vel_hmix_leith_init(err)!{{{ + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag - !-------------------------------------------------------------------- - ! - ! set some local module variables based on input config choices - ! - !-------------------------------------------------------------------- + ! End preamble + !----------------------------------------------------------------- + ! Begin code - err = 0 + !*** initialize return flag and set default values - hmixLeithOn = .false. + err = 0 + hmixLeithOff = .true. + leithParam = 0.0_RKIND + dxLeith = 0.0_RKIND + viscMaxLeith = 0.0_RKIND + sqrt3fact = sqrt(3.0_RKIND) - if (config_use_leith_del2) then - hmixLeithOn = .true. - endif + !*** reset values based on input configuration - !-------------------------------------------------------------------- + if (config_use_leith_del2) then + hmixLeithOff = .false. + leithParam = config_leith_parameter + dxLeith = config_leith_dx + viscMaxLeith = config_leith_visc2_max + endif + + !----------------------------------------------------------------- end subroutine ocn_vel_hmix_leith_init!}}} diff --git a/src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F b/src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F index e02c474516..095da74024 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F +++ b/src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F @@ -22,13 +22,12 @@ module ocn_vel_pressure_grad use mpas_timer - use mpas_derived_types - use mpas_pool_routines use mpas_constants use mpas_log use ocn_constants use ocn_config + use ocn_mesh implicit none private @@ -55,11 +54,26 @@ module ocn_vel_pressure_grad ! !-------------------------------------------------------------------- -! character (len=StrKIND), pointer :: config_pressure_gradient_type -! real (kind=RKIND), pointer :: config_common_level_weight - logical :: pgradOn - real (kind=RKIND) :: density0Inv, gdensity0Inv, inv12 - + logical :: & + pGradOff ! flag for turning pressure gradient on/off + + integer :: & + pGradType ! id for pressure radient method selected + + integer, parameter :: &! ids for supported methods + pGradTypeNone = 0, &! none selected + pGradTypeSSHgrad = 1, &! ssh gradient + pGradTypePZmid = 2, &! pressure and zMid + pGradTypeMontPot = 3, &! Montgomery potential + pGradTypeMontPotDens = 4, &! Montgomery potential and density + pGradTypeJacobDens = 5, &! Jacobian from density + pGradTypeJacobTS = 6 ! Jacobian from T,S + + real (kind=RKIND) :: &! precomputed constants for efficiency + density0Inv, &! 1/density0 + gdensity0Inv, &! g/density0 + inv12, &! 1/12 + pGradLvlWgt ! weighting for levels in Jacobian formulations !*********************************************************************** @@ -78,282 +92,456 @@ module ocn_vel_pressure_grad ! !----------------------------------------------------------------------- - subroutine ocn_vel_pressure_grad_tend(meshPool, ssh, pressure, montgomeryPotential, zMid, density, potentialDensity, & - indexT, indexS, tracers, tend, err, inSituThermalExpansionCoeff,inSituSalineContractionCoeff)!{{{ + subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & + montgomeryPotential, zMid, & + density, potentialDensity, & + indxT, indxS, tracers, & + thermExpCoeff, salineContractCoeff, & + tend, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- - integer, intent(in) :: indexT, indexS + integer, intent(in) :: & + indxT, &!< [in] tracer array index for temperature + indxS !< [in] tracer array index for salt real (kind=RKIND), dimension(:), intent(in) :: & - ssh !< Input: Sea surface height + ssh !< [in] sea surface height real (kind=RKIND), dimension(:,:), intent(in) :: & - pressure, & !< Input: Pressure field - montgomeryPotential, & !< Input: Mongomery potential - zMid, & !< Input: z-coordinate at mid-depth of layer - density, & !< Input: density - potentialDensity !< Input: potentialDensity + pressure, &!< [in] Pressure field + montgomeryPotential, &!< [in] Mongomery potential + zMid, &!< [in] z-coordinate at layer mid-depth + density, &!< [in] density + potentialDensity !< [in] potentialDensity - real (kind=RKIND), dimension(:,:), intent(in), optional :: & - inSituThermalExpansionCoeff, & - inSituSalineContractionCoeff - - real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers + real (kind=RKIND), dimension(:,:), intent(in) :: & + thermExpCoeff, &!< [in] in situ thermal expansion coeff + salineContractCoeff !< [in] in situ saline contraction coeff - type (mpas_pool_type), intent(in) :: & - meshPool !< Input: mesh information + real (kind=RKIND), dimension(:,:,:), intent(in) :: & + tracers !< [in] array of active tracers !----------------------------------------------------------------- - ! ! input/output variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(inout) :: & - tend !< Input/Output: velocity tendency + tend !< [inout] accumulated velocity tendency !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer :: iEdge, k, cell1, cell2, iCell, kMax, nEdges - integer, pointer :: nVertLevels - integer, dimension(:), pointer :: nEdgesArray - integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell - integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask + integer :: & + iEdge, iCell, k, &! loop indices for edge, cell, vertical loops + cell1, cell2, &! neighbor cell indices across edge + kMax ! deepest active layer - real (kind=RKIND), dimension(:), pointer :: dcEdge - real (kind=RKIND), dimension(:), allocatable :: JacobianDxDs,JacobianTz,JacobianSz,T1,T2,S1,S2 - real (kind=RKIND), dimension(:,:), allocatable :: FXTop, work1, work2 - real (kind=RKIND) :: invdcEdge, pGrad, sumAJTop, AJTop, FC, FCPrev, alpha, beta + real (kind=RKIND) :: & + invdcEdge, &! temporary for 1/dcEdge + pGrad, &! vertically integrated pgrad + alpha, beta, &! T,S expansion factors for Jacobian + Area, &! Area for Jacobian calculation + zStar, zC, zGamma,&! z coordinate combination for Jacobian + rhoL, rhoR, &! density neighbors for Jacobian + TL, TR, &! temperature neighbors for Jacobian + SL, SR ! salinity neighbors for Jacobian - err = 0 + real (kind=RKIND), dimension(:,:), allocatable :: & + JacobianDxDs, &! Jacobian in density * Dx * DS + JacobianTz, &! Jacobian associated with temperature + JacobianSz ! Jacobian associated with salinity - if (.not. pgradOn) return + ! End preamble + !----------------------------------------------------------------- + ! Begin code - call mpas_timer_start("pressure grad") + !*** set error code and exit if turned off + !*** start timer if turned on - call mpas_pool_get_dimension(meshPool, 'nVertLevels',nVertLevels) - call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) - call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge) - call mpas_pool_get_array(meshPool, 'edgeMask', edgeMask) + err = 0 + if (pgradOff) return + call mpas_timer_start("pressure grad") - nEdges = nEdgesArray( 1 ) + !*** Compute pGrad based on method selected + select case (pGradType) - if (config_pressure_gradient_type.eq.'ssh_gradient') then + case (pGradTypeSSHgrad) - ! pressure for sea surface height - ! - g grad ssh + ! pressure for sea surface height = - g grad ssh + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, & + !$acc tend, edgeMask, ssh) & + !$acc private(cell1, cell2, invdcEdge, k, kMax) + #else !$omp parallel - !$omp do schedule(runtime) private(cell1, cell2, invdcEdge, k) - do iEdge=1,nEdges + !$omp do schedule(runtime) & + !$omp private(cell1, cell2, invdcEdge, k, kMax) + #endif + do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) invdcEdge = 1.0_RKIND / dcEdge(iEdge) + kMax = maxLevelEdgeTop(iEdge) - do k=1,maxLevelEdgeTop(iEdge) - tend(k,iEdge) = tend(k,iEdge) - gravity * edgeMask(k,iEdge) * invdcEdge * ( ssh(cell2) - ssh(cell1) ) + do k=1,kMax + tend(k,iEdge) = tend(k,iEdge) - & + gravity*edgeMask(k,iEdge)*invdcEdge* & + (ssh(cell2) - ssh(cell1)) end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif - elseif (config_pressure_gradient_type.eq.'pressure_and_zmid') then + case (pGradTypePZmid) ! pressure for generalized coordinates ! -1/density_0 (grad p_k + density g grad z_k^{mid}) + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, zMid, & + !$acc tend, edgeMask, pressure, density) & + !$acc private(cell1, cell2, invdcEdge, k, kMax) + #else !$omp parallel - !$omp do schedule(runtime) private(cell1, cell2, invdcEdge, k) - do iEdge=1,nEdges + !$omp do schedule(runtime) & + !$omp private(cell1, cell2, invdcEdge, k, kMax) + #endif + do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) invdcEdge = 1.0_RKIND / dcEdge(iEdge) + kMax = maxLevelEdgeTop(iEdge) - do k=1,maxLevelEdgeTop(iEdge) - tend(k,iEdge) = tend(k,iEdge) + edgeMask(k,iEdge) * invdcEdge * ( & - - density0Inv * ( pressure(k,cell2) - pressure(k,cell1) ) & - - gdensity0Inv * 0.5_RKIND*(density(k,cell1)+density(k,cell2)) * ( zMid(k,cell2) - zMid(k,cell1) ) ) + do k=1,kMax + tend(k,iEdge) = tend(k,iEdge) + & + edgeMask(k,iEdge)*invdcEdge*( & + - density0Inv*(pressure(k,cell2) - & + pressure(k,cell1)) & + - gdensity0Inv*0.5_RKIND* & + (density(k,cell1)+density(k,cell2))* & + (zMid(k,cell2)-zMid(k,cell1))) end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif - elseif (config_pressure_gradient_type.eq.'MontgomeryPotential') then + case (pGradTypeMontPot) ! For pure isopycnal coordinates, this is just grad(M), ! the gradient of Montgomery Potential + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, & + !$acc tend, edgeMask, montgomeryPotential) & + !$acc private(cell1, cell2, invdcEdge, k, kMax) + #else !$omp parallel - !$omp do schedule(runtime) private(cell1, cell2, invdcEdge, k) - do iEdge=1,nEdges + !$omp do schedule(runtime) & + !$omp private(cell1, cell2, invdcEdge, k, kMax) + #endif + do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) invdcEdge = 1.0_RKIND / dcEdge(iEdge) + kMax = maxLevelEdgeTop(iEdge) - do k=1,maxLevelEdgeTop(iEdge) - tend(k,iEdge) = tend(k,iEdge) + edgeMask(k,iEdge) * invdcEdge * ( & - - ( montgomeryPotential(k,cell2) - montgomeryPotential(k,cell1) ) ) + do k=1,kMax + tend(k,iEdge) = tend(k,iEdge) + & + edgeMask(k,iEdge)*invdcEdge* ( & + - (montgomeryPotential(k,cell2) - & + montgomeryPotential(k,cell1))) end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif - elseif (config_pressure_gradient_type.eq.'MontgomeryPotential_and_density') then + case (pGradTypeMontPotDens) - ! This formulation has not been extensively tested and is not supported at this time. + ! This formulation has not been extensively tested and is not + ! supported at this time. ! This is -grad(M)+p grad(1/rho) ! Where rho is the potential density. ! See Bleck (2002) equation 1, and last equation in Appendix A. + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, dcEdge, maxLevelEdgeTop, & + !$acc edgeMask, tend, montgomeryPotential, & + !$acc pressure, potentialDensity) & + !$acc private(cell1, cell2, invdcEdge, k, kMax) + #else !$omp parallel - !$omp do schedule(runtime) private(cell1, cell2, invdcEdge, k) - do iEdge=1,nEdges + !$omp do schedule(runtime) & + !$omp private(cell1, cell2, invdcEdge, k, kMax) + #endif + do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) invdcEdge = 1.0_RKIND / dcEdge(iEdge) + kMax = maxLevelEdgeTop(iEdge) - do k=1,maxLevelEdgeTop(iEdge) - tend(k,iEdge) = tend(k,iEdge) + edgeMask(k,iEdge) * invdcEdge * ( & - - ( montgomeryPotential(k,cell2) - montgomeryPotential(k,cell1) ) & - + 0.5_RKIND*(pressure(k,cell1)+pressure(k,cell2)) * ( 1.0_RKIND/potentialDensity(k,cell2) & - - 1.0_RKIND/potentialDensity(k,cell1) ) ) + do k=1,kMax + tend(k,iEdge) = tend(k,iEdge) + & + edgeMask(k,iEdge)*invdcEdge*( & + - (montgomeryPotential(k,cell2) - & + montgomeryPotential(k,cell1)) & + + 0.5_RKIND*(pressure(k,cell1) + & + pressure(k,cell2))* & + ( 1.0_RKIND/potentialDensity(k,cell2) & + - 1.0_RKIND/potentialDensity(k,cell1))) end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif - elseif (config_pressure_gradient_type.eq.'Jacobian_from_density') then + case (pGradTypeJacobDens) - allocate(JacobianDxDs(nVertLevels)) + allocate(JacobianDxDs(nVertLevels,nEdgesAll)) + !$acc enter data create (JacobianDxDs) + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, & + !$acc edgeMask, pressure, density, zMid, tend, & + !$acc JacobianDxDs) & + !$acc private(cell1, cell2, invdcEdge, k, kMax, pGrad, & + !$acc Area, zStar, zC, zGamma, rhoL, rhoR) + #else !$omp parallel - !$omp do schedule(runtime) private(cell1, cell2, invdcEdge, k, pGrad, JacobianDxDs) - do iEdge=1,nEdges + !$omp do schedule(runtime) & + !$omp private(cell1, cell2, invdcEdge, k, kMax, pGrad, & + !$omp Area, zStar, zC, zGamma, rhoL, rhoR) + #endif + do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) invdcEdge = 1.0_RKIND / dcEdge(iEdge) + kMax = maxLevelEdgeTop(iEdge) - call pGrad_Jacobian_common_level(density(:,cell1),density(:,cell2),zMid(:,cell1),zMid(:,cell2), & - maxLevelEdgeTop(iEdge), config_common_level_weight, JacobianDxDs) + ! Compute the density-Jacobian in common_level form. + ! See Shchepetkin and McWilliams (2003) Ocean Modeling, + ! sections 2-4 + + JacobianDxDs(1,iEdge) = 0.0_RKIND + + do k=2,kMax + + ! eqn 2.7 in Shchepetkin and McWilliams (2003) + ! Note delta x was removed. It must be an error in the + ! paper, ! as it makes the units incorrect. + Area = 0.5_RKIND*(zMid(k-1,cell1) - zMid(k,cell1) + & + zMid(k-1,cell2) - zMid(k,cell2)) + + ! eqn 2.8 + zStar = ( zMid(k-1,cell2)*zMid(k-1,cell1) - & + zMid(k ,cell2)*zMid(k ,cell1) )/ & + ( zMid(k-1,cell2)-zMid(k ,cell2) + & + zMid(k-1,cell1)-zMid(k,cell1)) + + ! eqn 3.2 + zC = 0.25_RKIND*(zMid(k,cell1) + zMid(k-1,cell1) + & + zMid(k,cell2) + zMid(k-1,cell2)) + + ! eqn 4.1 + zGamma = (1.0_RKIND - pGradLvlWgt)*zStar + pGradLvlWgt*zC + + rhoL = (density(k ,cell1)*(zMid(k-1,cell1)-zGamma) + & + density(k-1,cell1)*(zGamma-zMid(k ,cell1)))/ & + (zMid(k-1,cell1) - zMid(k,cell1)) + rhoR = (density(k ,cell2)*(zMid(k-1,cell2)-zGamma) + & + density(k-1,cell2)*(zGamma-zMid(k ,cell2)))/ & + (zMid(k-1,cell2) - zMid(k,cell2)) + + ! eqn 2.6 in Shchepetkin and McWilliams (2003) + JacobianDxDs(k,iEdge) = Area * (rhoL - rhoR) + end do ! In layer 1, use pressure for generalized coordinates ! pGrad = -1/density_0 (grad p_k + density g grad z_k^{mid}) k = 1 - pGrad = edgeMask(k,iEdge) * invdcEdge * ( & - - density0Inv * ( pressure(k,cell2) - pressure(k,cell1) ) & - - gdensity0Inv * 0.5_RKIND*(density(k,cell1)+density(k,cell2)) * ( zMid(k,cell2) - zMid(k,cell1) ) ) + pGrad = edgeMask(k,iEdge)*invdcEdge*( & + - density0Inv*(pressure(k,cell2) - & + pressure(k,cell1)) & + - gdensity0Inv*0.5_RKIND* & + (density(k,cell1)+density(k,cell2))* & + (zMid(k,cell2)- zMid(k,cell1) ) ) tend(k,iEdge) = tend(k,iEdge) + pGrad - do k=2,maxLevelEdgeTop(iEdge) + do k=2,kMax ! note JacobianDxDs includes negative sign, so ! pGrad is - g/rho_0 dP/dx - pGrad = pGrad + gdensity0Inv * JacobianDxDs(k) * invdcEdge + pGrad = pGrad + gdensity0Inv*JacobianDxDs(k,iEdge)*invdcEdge tend(k,iEdge) = tend(k,iEdge) + pGrad end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif + !$acc exit data delete(JacobianDxDs) deallocate(JacobianDxDs) - elseif (config_pressure_gradient_type.eq.'Jacobian_from_TS') then - - allocate(JacobianDxDs(nVertLevels),JacobianTz(nVertLevels),JacobianSz(nVertLevels), T1(nVertLevels)) - allocate(T2(nVertLevels), S1(nVertLevels), S2(nVertLevels)) - + case (pGradTypeJacobTS) + + allocate(JacobianDxDs(nVertLevels,nEdgesAll), & + JacobianTz(nVertLevels,nEdgesAll), & + JacobianSz(nVertLevels,nEdgesAll)) + !$acc enter data create(JacobianDxDs, & + !$acc JacobianTz, JacobianSz) + + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, & + !$acc zMid, tracers, edgeMask, pressure, & + !$acc density, tend, thermExpCoeff, & + !$acc salineContractCoeff, & + !$acc JacobianTz, JacobianSz, JacobianDxDs) & + !$acc private(cell1, cell2, invdcEdge, k,kMax, alpha, beta,& + !$acc TL, TR, SL, SR, Area, zStar, zC, zGamma, pGrad) + #else !$omp parallel !$omp do schedule(runtime) & - !$omp private(cell1, cell2, invdcEdge, kMax, T1, T2, S1, S2, k, alpha, beta, & - !$omp pGrad, JacobianDxDs, JacobianTz, JacobianSz) - do iEdge=1,nEdges + !$omp private(cell1, cell2, invdcEdge, k, kMax, alpha, beta, & + !$omp TL, TR, SL, SR, Area, zStar, zC, zGamma, pGrad) + #endif + do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) invdcEdge = 1.0_RKIND / dcEdge(iEdge) kMax = maxLevelEdgeTop(iEdge) - ! copy T and S to local column arrays - T1(1:kMax) = tracers(indexT,1:kMax,cell1) - T2(1:kMax) = tracers(indexT,1:kMax,cell2) - S1(1:kMax) = tracers(indexS,1:kMax,cell1) - S2(1:kMax) = tracers(indexS,1:kMax,cell2) + ! compute J(T,z) and J(S,z) + ! in Shchepetkin and McWilliams (2003) (7.16) - ! compute J(T,z) and J(S,z) in Shchepetkin and McWilliams (2003) (7.16) - call pGrad_Jacobian_common_level(T1, T2 ,zMid(:,cell1),zMid(:,cell2),kMax,config_common_level_weight, JacobianTz) - call pGrad_Jacobian_common_level(S1, S2 ,zMid(:,cell1),zMid(:,cell2),kMax,config_common_level_weight, JacobianSz) + JacobianTz(1,iEdge) = 0.0_RKIND + JacobianSz(1,iEdge) = 0.0_RKIND + + do k=2,kMax + + ! eqn 2.7 in Shchepetkin and McWilliams (2003) + ! Note delta x was removed. It must be an error in the + ! paper, ! as it makes the units incorrect. + Area = 0.5_RKIND*(zMid(k-1,cell1) - zMid(k,cell1) + & + zMid(k-1,cell2) - zMid(k,cell2)) + + ! eqn 2.8 + zStar = ( zMid(k-1,cell2)*zMid(k-1,cell1) - & + zMid(k ,cell2)*zMid(k ,cell1) )/ & + ( zMid(k-1,cell2)-zMid(k ,cell2) + & + zMid(k-1,cell1)-zMid(k,cell1)) + + ! eqn 3.2 + zC = 0.25_RKIND*(zMid(k,cell1) + zMid(k-1,cell1) + & + zMid(k,cell2) + zMid(k-1,cell2)) + + ! eqn 4.1 + zGamma = (1.0_RKIND - pGradLvlWgt)*zStar + pGradLvlWgt*zC + + + TL = (tracers(indxT,k ,cell1)*(zMid(k-1,cell1)-zGamma) + & + tracers(indxT,k-1,cell1)*(zGamma-zMid(k ,cell1)))/ & + (zMid(k-1,cell1) - zMid(k,cell1)) + TR = (tracers(indxT,k ,cell2)*(zMid(k-1,cell2)-zGamma) + & + tracers(indxT,k-1,cell2)*(zGamma-zMid(k ,cell2)))/ & + (zMid(k-1,cell2) - zMid(k,cell2)) + + SL = (tracers(indxS,k ,cell1)*(zMid(k-1,cell1)-zGamma) + & + tracers(indxS,k-1,cell1)*(zGamma-zMid(k ,cell1)))/ & + (zMid(k-1,cell1) - zMid(k,cell1)) + SR = (tracers(indxS,k ,cell2)*(zMid(k-1,cell2)-zGamma) + & + tracers(indxS,k-1,cell2)*(zGamma-zMid(k ,cell2)))/ & + (zMid(k-1,cell2) - zMid(k,cell2)) + + + ! eqn 2.6 in Shchepetkin and McWilliams (2003) + JacobianTz(k,iEdge) = Area*(TL - TR) + JacobianSz(k,iEdge) = Area*(SL - SR) + end do ! In layer 1, use pressure for generalized coordinates ! pGrad = -1/density_0 (grad p_k + density g grad z_k^{mid}) + k = 1 - pGrad = edgeMask(k,iEdge) * invdcEdge * ( & - - density0Inv * ( pressure(k,cell2) - pressure(k,cell1) ) & - - gdensity0Inv * 0.5_RKIND*(density(k,cell1)+density(k,cell2)) * ( zMid(k,cell2) - zMid(k,cell1) ) ) + pGrad = edgeMask(k,iEdge)*invdcEdge*( & + - density0Inv*(pressure(k,cell2) - & + pressure(k,cell1)) & + - gdensity0Inv*0.5_RKIND* & + (density(k,cell1)+density(k,cell2))* & + (zMid(k,cell2)- zMid(k,cell1) ) ) tend(k,iEdge) = tend(k,iEdge) + pGrad do k=2,kMax ! Average alpha and beta over four data points of the Jacobian cell. - ! Note that inSituThermalExpansionCoeff and inSituSalineContractionCoeff include a 1/density factor, + ! Note that thermExpCoeff and salineContractCoeff include a 1/density factor, ! so must multiply by density here. - alpha = 0.25_RKIND*( density(k,cell1)*inSituThermalExpansionCoeff (k,cell1) & - + density(k-1,cell1)*inSituThermalExpansionCoeff (k-1,cell1) & - + density(k,cell2)*inSituThermalExpansionCoeff (k,cell2) & - + density(k-1,cell2)*inSituThermalExpansionCoeff (k-1,cell2) ) - beta = 0.25_RKIND*( density(k,cell1)*inSituSalineContractionCoeff(k,cell1) & - + density(k-1,cell1)*inSituSalineContractionCoeff(k-1,cell1) & - + density(k,cell2)*inSituSalineContractionCoeff(k,cell2) & - + density(k-1,cell2)*inSituSalineContractionCoeff(k-1,cell2) ) + alpha = 0.25_RKIND*( & + density(k ,cell1)*thermExpCoeff (k ,cell1) & + + density(k-1,cell1)*thermExpCoeff (k-1,cell1) & + + density(k ,cell2)*thermExpCoeff (k ,cell2) & + + density(k-1,cell2)*thermExpCoeff (k-1,cell2) ) + beta = 0.25_RKIND*( & + density(k ,cell1)*salineContractCoeff(k ,cell1) & + + density(k-1,cell1)*salineContractCoeff(k-1,cell1) & + + density(k ,cell2)*salineContractCoeff(k ,cell2) & + + density(k-1,cell2)*salineContractCoeff(k-1,cell2) ) ! Shchepetkin and McWilliams (2003) (7.16) - JacobianDxDs(k) = -alpha*JacobianTz(k) + beta*JacobianSz(k) + JacobianDxDs(k,iEdge) = -alpha*JacobianTz(k,iEdge) + & + beta*JacobianSz(k,iEdge) ! note JacobianDxDs includes negative sign, so ! pGrad is - g/rho_0 dP/dx - pGrad = pGrad + gdensity0Inv * JacobianDxDs(k) * invdcEdge + pGrad = pGrad + gdensity0Inv * JacobianDxDs(k,iEdge) * invdcEdge tend(k,iEdge) = tend(k,iEdge) + pGrad end do end do + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel + #endif - deallocate(JacobianDxDs,JacobianTz,JacobianSz, T1, T2, S1, S2) + !$acc exit data delete (JacobianDxDs, JacobianTz, JacobianSz) + deallocate(JacobianDxDs, JacobianTz, JacobianSz) - else + case default - call mpas_log_write(' Pressure type is: '//trim(config_pressure_gradient_type)) - call mpas_log_write(' Incorrect choice of config_pressure_gradient_type.',MPAS_LOG_CRIT) - err = 1 + !*** error already trapped during init - endif + end select ! pGradType call mpas_timer_stop("pressure grad") @@ -363,331 +551,84 @@ end subroutine ocn_vel_pressure_grad_tend!}}} !*********************************************************************** ! -! routine pGrad_Jacobian_common_level -! -!> \brief Computes density-Jacobian -!> \author Mark Petersen -!> \date February 2014 -!> \details -!> This routine computes the density-Jacobian in common_level form. -!> See Shchepetkin and McWilliams (2003) Ocean Modeling, sections 2-4 -! -!----------------------------------------------------------------------- - - subroutine pGrad_Jacobian_common_level(rho1,rho2,z1,z2,kMax,gamma,JacobianDxDs) - - !----------------------------------------------------------------- - ! - ! input variables - ! - !----------------------------------------------------------------- - - real (kind=RKIND), dimension(:), intent(in) :: & - rho1, & ! density of column 1 - rho2, & ! density of column 2 - z1, & ! z-coordinate at middle of cell, column 1 - z2 ! z-coordinate at middle of cell, column 2 - - real (kind=RKIND), intent(in) :: & - gamma ! weight between zStar (original Jacobian) and z_C (weighted Jacobian) - - integer, intent(in) :: & - kMax - - !----------------------------------------------------------------- - ! - ! input/output variables - ! - !----------------------------------------------------------------- - - !----------------------------------------------------------------- - ! - ! output variables - ! - !----------------------------------------------------------------- - - real (kind=RKIND), dimension(:), intent(out) :: & - JacobianDxDs ! - Delta x Delta s J(rho,z) - - !----------------------------------------------------------------- - ! - ! local variables - ! - !----------------------------------------------------------------- - - integer :: k - real (kind=RKIND) :: Area, zStar, rhoL, rhoR, zC, zGamma - - JacobianDxDs = 0.0_RKIND - - do k=2,kMax - - ! eqn 2.7 in Shchepetkin and McWilliams (2003) - ! Note delta x was removed. It must be an error in the paper, - ! as it makes the units incorrect. - Area = 0.5_RKIND*(z1(k-1) - z1(k) + z2(k-1) - z2(k) ) - - ! eqn 2.8 - zStar = ( z2(k-1)*z1(k-1) - z2(k)*z1(k) )/(z2(k-1)-z2(k) + z1(k-1)-z1(k)) - - ! eqn 3.2 - zC = 0.25_RKIND*( z1(k) + z1(k-1) + z2(k) + z2(k-1) ) - - ! eqn 4.1 - zGamma = (1.0_RKIND - gamma)*zStar + gamma*zC - - rhoL = (rho1(k)*(z1(k-1)-zGamma) + rho1(k-1)*(zGamma-z1(k)))/(z1(k-1) - z1(k)) - rhoR = (rho2(k)*(z2(k-1)-zGamma) + rho2(k-1)*(zGamma-z2(k)))/(z2(k-1) - z2(k)) - - ! eqn 2.6 in Shchepetkin and McWilliams (2003) - JacobianDxDs(k) = Area * (rhoL - rhoR) - end do - - end subroutine pGrad_Jacobian_common_level - -!*********************************************************************** -! -! routine pGrad_Jacobian_POM_SCRUM -! -!> \brief Computes density-Jacobian -!> \author Mark Petersen -!> \date February 2014 -!> \details -!> This routine computes the density-Jacobian in POM/SCRUM form. -!> See Shchepetkin and McWilliams (2003) Ocean Modeling, section 2. -! -!----------------------------------------------------------------------- - - subroutine pGrad_Jacobian_POM_SCRUM(rho1,rho2,z1,z2,kMax,JacobianDxDs) - - !----------------------------------------------------------------- - ! - ! input variables - ! - !----------------------------------------------------------------- - - real (kind=RKIND), dimension(:), intent(in) :: & - rho1, & ! density of column 1 - rho2, & ! density of column 2 - z1, & ! z-coordinate at middle of cell, column 1 - z2 ! z-coordinate at middle of cell, column 2 - - integer, intent(in) :: & - kMax ! maximum level - - !----------------------------------------------------------------- - ! - ! input/output variables - ! - !----------------------------------------------------------------- - - !----------------------------------------------------------------- - ! - ! output variables - ! - !----------------------------------------------------------------- - - real (kind=RKIND), dimension(:), intent(out) :: & - JacobianDxDs ! - Delta x Delta s J(rho,z) - - !----------------------------------------------------------------- - ! - ! local variables - ! - !----------------------------------------------------------------- - - integer :: k - - JacobianDxDs = 0.0_RKIND - - do k=2,kMax - - ! eqn 2.3 in Shchepetkin and McWilliams (2003) - JacobianDxDs(k) = 0.25_RKIND*(& - (rho1(k) + rho1(k-1) - rho2(k) - rho2(k-1) )*(z1(k-1) - z1(k) + z2(k-1) - z2(k) ) & - - (rho1(k-1) - rho1(k) + rho2(k-1) - rho2(k) )*(z1(k) + z1(k-1) - z2(k) - z2(k-1) ) ) - end do - - end subroutine pGrad_Jacobian_POM_SCRUM - -!*********************************************************************** -! -! routine pGrad_Jacobian_diagonal -! -!> \brief Computes density-Jacobian -!> \author Mark Petersen -!> \date February 2014 -!> \details -!> This routine computes the density-Jacobian in diagonal form. -!> See Shchepetkin and McWilliams (2003) Ocean Modeling, section 2. -! -!----------------------------------------------------------------------- - - subroutine pGrad_Jacobian_diagonal(rho1,rho2,z1,z2,kMax,JacobianDxDs) - - !----------------------------------------------------------------- - ! - ! input variables - ! - !----------------------------------------------------------------- - - real (kind=RKIND), dimension(:), intent(in) :: & - rho1, & ! density of column 1 - rho2, & ! density of column 2 - z1, & ! z-coordinate at middle of cell, column 1 - z2 ! z-coordinate at middle of cell, column 2 - - integer, intent(in) :: & - kMax - - !----------------------------------------------------------------- - ! - ! input/output variables - ! - !----------------------------------------------------------------- - - !----------------------------------------------------------------- - ! - ! output variables - ! - !----------------------------------------------------------------- - - real (kind=RKIND), dimension(:), intent(out) :: & - JacobianDxDs ! - Delta x Delta s J(rho,z) - - !----------------------------------------------------------------- - ! - ! local variables - ! - !----------------------------------------------------------------- - - integer :: k - - JacobianDxDs = 0.0_RKIND - - do k=2,kMax - - ! eqn 2.5 in Shchepetkin and McWilliams (2003) - JacobianDxDs(k) = 0.5_RKIND*( & - (rho1(k-1) - rho2(k))*(z2(k-1) - z1(k) ) & - + (rho1(k) - rho2(k-1))*(z1(k-1) - z2(k)) ) - end do - - end subroutine pGrad_Jacobian_diagonal - -!*********************************************************************** -! -! routine pGrad_Jacobian_pseudo_flux +! routine ocn_vel_pressure_grad_init ! -!> \brief Computes density-Jacobian +!> \brief Initializes ocean momentum horizontal pressure gradient !> \author Mark Petersen -!> \date February 2014 +!> \date September 2011 !> \details -!> This routine computes the density-Jacobian in pseudo_flux form. -!> See Shchepetkin and McWilliams (2003) Ocean Modeling, section 2. +!> This routine initializes parameters required for the computation of +!> the horizontal pressure gradient. ! !----------------------------------------------------------------------- - subroutine pGrad_Jacobian_pseudo_flux(rho1,rho2,z1,z2,kMax,JacobianDxDs) - - !----------------------------------------------------------------- - ! - ! input variables - ! - !----------------------------------------------------------------- - - real (kind=RKIND), dimension(:), intent(in) :: & - rho1, & ! density of column 1 - rho2, & ! density of column 2 - z1, & ! z-coordinate at middle of cell, column 1 - z2 ! z-coordinate at middle of cell, column 2 - - integer, intent(in) :: & - kMax - - !----------------------------------------------------------------- - ! - ! input/output variables - ! - !----------------------------------------------------------------- + subroutine ocn_vel_pressure_grad_init(err)!{{{ !----------------------------------------------------------------- - ! - ! output variables - ! + ! Output Variables !----------------------------------------------------------------- - real (kind=RKIND), dimension(:), intent(out) :: & - JacobianDxDs ! - Delta x Delta s J(rho,z) + integer, intent(out) :: err !< [out] error flag + ! End preamble !----------------------------------------------------------------- - ! - ! local variables - ! - !----------------------------------------------------------------- - - integer :: k - real (kind=RKIND) :: FLeft, FTop, FRight, FBottom + ! Begin code - JacobianDxDs = 0.0_RKIND + !*** Initialize error code and default values - do k=2,kMax - - FLeft = 0.5_RKIND*( rho1(k) + rho1(k-1) ) * (z1(k-1) - z1(k)) - FTop = 0.5_RKIND*( rho1(k-1) + rho2(k-1) ) * (z2(k-1) - z1(k-1)) - FRight = 0.5_RKIND*( rho2(k) + rho2(k-1) ) * (z2(k-1) - z2(k)) - FBottom = 0.5_RKIND*( rho1(k) + rho2(k) ) * (z2(k) - z1(k)) - - ! eqn 2.11 in Shchepetkin and McWilliams (2003) - JacobianDxDs(k) = FLeft + FTop - FRight - FBottom - end do - - end subroutine pGrad_Jacobian_pseudo_flux - -!*********************************************************************** -! -! routine ocn_vel_pressure_grad_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. -! -!----------------------------------------------------------------------- + err = 0 + pGradOff = .true. + pGradType = pGradTypeNone + pGradLvlWgt = 0.0_RKIND - subroutine ocn_vel_pressure_grad_init(err)!{{{ + density0Inv = 1.0_RKIND / rho_sw + gdensity0Inv = gravity / rho_sw + inv12 = 1.0_RKIND / 12.0_RKIND - !-------------------------------------------------------------------- + !*** Now reset variables based on input configuration + if (config_disable_vel_pgrad) then + pGradOff = .true. + return + else + pGradOff = .false. + pGradLvlWgt = config_common_level_weight + endif - !----------------------------------------------------------------- - ! - ! Output Variables - ! - !----------------------------------------------------------------- + select case (trim(config_pressure_gradient_type)) - integer, intent(out) :: err !< Output: error flag + case ('ssh_gradient') + pGradType = pGradTypeSSHgrad + call mpas_log_write(' Pressure type is: ssh_gradient') + case ('pressure_and_zmid') + pGradType = pGradTypePZmid + call mpas_log_write(' Pressure type is: pressure_and_zmid') - !----------------------------------------------------------------- - ! - ! call individual init routines for each parameterization - ! - !----------------------------------------------------------------- + case ('MontgomeryPotential') + pGradType = pGradTypeMontPot + call mpas_log_write(' Pressure type is: MontgomeryPotential') - err = 0 + case ('MontgomeryPotential_and_density') + pGradType = pGradTypeMontPotDens + call mpas_log_write( & + ' Pressure type is: MontgomeryPotential_and_density') - pgradOn = .true. + case ('Jacobian_from_density') + pGradType = pGradTypeJacobDens + call mpas_log_write(' Pressure type is: Jacobian_from_density') - density0Inv = 1.0_RKIND / rho_sw - gdensity0Inv = gravity / rho_sw - inv12 = 1.0_RKIND / 12.0_RKIND + case ('Jacobian_from_TS') + pGradType = pGradTypeJacobTS + call mpas_log_write(' Pressure type is: Jacobian_from_TS') - if (config_disable_vel_pgrad) pgradOn = .false. + case default + call mpas_log_write( & + ' Incorrect choice of config_pressure_gradient_type.', & + MPAS_LOG_CRIT) + err = 1 - call mpas_log_write(' Pressure type is: '//trim(config_pressure_gradient_type)) + end select !-------------------------------------------------------------------- @@ -698,3 +639,4 @@ end subroutine ocn_vel_pressure_grad_init!}}} end module ocn_vel_pressure_grad !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F b/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F index b6c5bb2115..67387101a4 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F +++ b/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F @@ -28,6 +28,8 @@ module ocn_vel_tidal_potential use mpas_timekeeping use mpas_timer use ocn_constants + use ocn_config + use ocn_mesh implicit none private @@ -56,7 +58,34 @@ module ocn_vel_tidal_potential ! !-------------------------------------------------------------------- - logical :: tidalPotentialOn + logical :: & + tidalPotentialOff ! on/off switch for tidal potential + + !*** physical constants + real (kind=RKIND) :: & + betaSelfAttrLoad, &! self-attraction and loading beta + tidalPotRamp ! scale for ramping tidal potential + + !*** eventually these will be private module arrays and + !*** not in the forcing pool so retain pointers here as placeholders + real (kind=RKIND), dimension(:), pointer :: & + tidalPotEta, &! + tidalConstituentAmplitude, &! + tidalConstituentFrequency, &! + tidalConstituentLoveNumbers, &! + tidalConstituentNodalAmplitude, &! + tidalConstituentAstronomical, &! + tidalConstituentNodalPhase ! + + real (kind=RKIND), dimension(:,:), pointer :: & + latitudeFunction + + integer, dimension(:), pointer :: & + tidalConstituentType + + integer :: & + nTidalConstituents ! number of tidal constituents + type :: char_array character(:), allocatable :: constituent end type @@ -80,84 +109,77 @@ module ocn_vel_tidal_potential ! !----------------------------------------------------------------------- - subroutine ocn_vel_tidal_potential_tend(meshPool, forcingPool, ssh, tend, err)!{{{ + subroutine ocn_vel_tidal_potential_tend(ssh, tend, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:), intent(in) :: & - ssh !< Input: Sea surface height - - type (mpas_pool_type), intent(in) :: meshPool !< Input: mesh information - type (mpas_pool_type), intent(in) :: forcingPool !< Input: forcinginformation + ssh !< [in] Sea surface height !----------------------------------------------------------------- - ! ! input/output variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(inout) :: & - tend !< Input/Output: velocity tendency + tend !< [inout] accumulated velocity tendency !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer :: iEdge, k, cell1, cell2, nEdges - integer, dimension(:), pointer :: nEdgesArray - integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell - integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask - - real (kind=RKIND), dimension(:), pointer :: dcEdge - real (kind=RKIND), dimension(:), pointer :: tidalPotentialEta - real (kind=RKIND), pointer :: config_self_attraction_and_loading_beta - real (kind=RKIND) :: invdcEdge - real (kind=RKIND) :: potentialGrad - logical, pointer :: config_use_tidal_potential_forcing - - err = 0 - - if (.not. tidalPotentialOn) return + integer :: & + iEdge, k, &! loop indices for edge, vertical loops + cell1, cell2, &! neighbor cell indices across edge + kMax ! deepest active layer on edge - call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) - call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge) - call mpas_pool_get_array(meshPool, 'edgeMask', edgeMask) - call mpas_pool_get_array(forcingPool, 'tidalPotentialEta', tidalPotentialEta) - call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing', config_use_tidal_potential_forcing) - call mpas_pool_get_config(ocnConfigs, 'config_self_attraction_and_loading_beta', config_self_attraction_and_loading_beta) + real (kind=RKIND) :: & + invdcEdge, &! 1/dcEdge + potentialGrad ! tmp potential gradient term - nEdges = nEdgesArray( 1 ) + ! End preamble + !----------------------------------------------------------------- + ! Begin code - !$omp do schedule(runtime) private(cell1, cell2, invdcEdge, potentialGrad, k) - do iEdge=1,nEdges + !*** Initialize error code and return early if not turned on + err = 0 + if (tidalPotentialOff) return + + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, edgeMask, & + !$acc tidalPotEta, ssh, tend) & + !$acc private(cell1, cell2, invdcEdge, potentialGrad, k, kMax) + #else + !$omp parallel do schedule(runtime) & + !$omp private(cell1, cell2, invdcEdge, potentialGrad, k, kMax) + #endif + do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) invdcEdge = 1.0_RKIND / dcEdge(iEdge) + kMax = maxLevelEdgeTop(iEdge) - potentialGrad = - gravity * invdcEdge * ( tidalPotentialEta(cell2) - tidalPotentialEta(cell1) & - + config_self_attraction_and_loading_beta*(ssh(cell2) - ssh(cell1))) + potentialGrad = - gravity*invdcEdge* & + (tidalPotEta(cell2) - tidalPotEta(cell1) + & + betaSelfAttrLoad* & + (ssh(cell2) - ssh(cell1))) - do k=1,maxLevelEdgeTop(iEdge) - tend(k,iEdge) = tend(k,iEdge) - edgeMask(k,iEdge) * potentialGrad + do k=1,kMax + tend(k,iEdge) = tend(k,iEdge) - & + edgeMask(k,iEdge)*potentialGrad end do end do + #ifndef MPAS_OPENACC !$omp end do + #endif end subroutine ocn_vel_tidal_potential_tend!}}} @@ -174,96 +196,74 @@ end subroutine ocn_vel_tidal_potential_tend!}}} ! !----------------------------------------------------------------------- - subroutine ocn_compute_tidal_potential_forcing(meshPool, forcingPool, diagnosticsPool, err)!{{{ + subroutine ocn_compute_tidal_potential_forcing(diagnosticsPool, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- - type (mpas_pool_type), intent(in) :: meshPool !< Input: mesh information - type (mpas_pool_type), intent(in) :: diagnosticsPool - !----------------------------------------------------------------- - ! - ! input/output variables - ! - !----------------------------------------------------------------- - type (mpas_pool_type), intent(inout) :: forcingPool + type (mpas_pool_type), intent(in) :: diagnosticsPool !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- integer, intent(out) :: err !< Output: Error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer, pointer :: nCells - real (kind=RKIND), pointer :: tidalPotentialRamp - real (kind=RKIND), dimension(:), pointer :: lonCell - integer, dimension(:), pointer :: maxLevelCell - - integer, pointer :: nTidalConstituents - real (kind=RKIND), dimension(:), pointer :: amplitude, frequency, loveNumbers - real (kind=RKIND), dimension(:), pointer :: nodalFactorAmplitude, nodalFactorPhase, astronomicalArgument - real (kind=RKIND), dimension(:,:), pointer :: latitudeFunction - integer, dimension(:), pointer :: constituentType - real (kind=RKIND), dimension(:), pointer :: eta - real (kind=RKIND), dimension(:,:), pointer :: zMid, tidalPotentialZMid real (kind=RKIND), pointer :: daysSinceStartOfSim - integer :: iCell - integer :: jCon, conType - real (kind=RKIND) :: lon,tArg,ramp,t - real (kind=RKIND) :: nCycles,period + integer :: & + iCell, jCon, &! loop indices for cell, constituent loops + conType ! id for constituent type + + real (kind=RKIND) :: & + lon,& + tArg,& + ramp,& + t, & + nCycles,& + period + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + !*** Set error flag and return early if not turned on err = 0 + if (tidalPotentialOff) return + + !*** Get pool variables (eventually private module variables?) + call mpas_pool_get_array(diagnosticsPool, "daysSinceStartOfSim", & + daysSinceStartOfSim) - if ( .not. tidalPotentialOn ) return - - call mpas_pool_get_config(ocnConfigs, 'config_tidal_potential_ramp', tidalPotentialRamp) - call mpas_pool_get_dimension(meshPool, 'nCells', nCells) - call mpas_pool_get_array(meshPool, 'lonCell', lonCell) - call mpas_pool_get_array(forcingPool, 'nTidalPotentialConstituents', nTidalConstituents) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentAmplitude', amplitude) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentFrequency', frequency) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentLoveNumbers', loveNumbers) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentNodalAmplitude', nodalFactorAmplitude) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentAstronomical', astronomicalArgument) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentNodalPhase', nodalFactorPhase) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentType', constituentType) - call mpas_pool_get_array(forcingPool, 'tidalPotentialLatitudeFunction', latitudeFunction) - call mpas_pool_get_array(forcingPool, 'tidalPotentialEta', eta) - call mpas_pool_get_array(diagnosticsPool, "daysSinceStartOfSim", daysSinceStartOfSim) - - ramp = tanh((2.0_RKIND*daysSinceStartOfSim)/tidalPotentialRamp) + ramp = tanh((2.0_RKIND*daysSinceStartOfSim)/tidalPotRamp) t = daysSinceStartOfSim*86400.0_RKIND - do iCell = 1, nCells - eta(iCell) = 0.0_RKIND + do iCell = 1, nCellsAll + tidalPotEta(iCell) = 0.0_RKIND end do + !*** Compute eta by summing all constituent contributions do jCon = 1, nTidalConstituents - period = 2.0_RKIND*pii/frequency(jCon) - nCycles = real(int(t/period),RKIND) - targ = frequency(jCon)*(t - nCycles*period) + nodalFactorPhase(jCon) + astronomicalArgument(jCon) - conType = constituentType(jCon) - do iCell = 1, nCells - lon = lonCell(iCell) - eta(iCell) = eta(iCell) + ramp & - * amplitude(jCon) & - * nodalFactorAmplitude(jCon) & - * loveNumbers(jCon) & - * latitudeFunction(iCell,conType+1) & - * cos(tArg + real(conType,RKIND)*lon) - end do + period = 2.0_RKIND*pii/tidalConstituentFrequency(jCon) + nCycles = real(int(t/period),RKIND) + targ = tidalConstituentFrequency(jCon)*(t - nCycles*period) + & + tidalConstituentNodalPhase(jCon) + & + tidalConstituentAstronomical(jCon) + conType = tidalConstituentType(jCon) + do iCell = 1, nCellsAll + lon = lonCell(iCell) + tidalPotEta(iCell) = tidalPotEta(iCell) + & + ramp * tidalConstituentAmplitude(jCon) & + * tidalConstituentNodalAmplitude(jCon) & + * tidalConstituentLoveNumbers(jCon) & + * latitudeFunction(iCell,conType+1) & + * cos(tArg + real(conType,RKIND)*lon) + end do end do @@ -285,139 +285,149 @@ end subroutine ocn_compute_tidal_potential_forcing!}}} subroutine ocn_vel_tidal_potential_init(domain,err)!{{{ + !----------------------------------------------------------------- + ! input/output variables + !----------------------------------------------------------------- + type (domain_type), intent(inout) :: domain - integer, intent(out) :: err !< Output: error flag - - logical, pointer :: config_use_tidal_potential_forcing - logical, pointer :: config_use_tidal_potential_forcing_M2 - logical, pointer :: config_use_tidal_potential_forcing_S2 - logical, pointer :: config_use_tidal_potential_forcing_N2 - logical, pointer :: config_use_tidal_potential_forcing_K2 - logical, pointer :: config_use_tidal_potential_forcing_K1 - logical, pointer :: config_use_tidal_potential_forcing_O1 - logical, pointer :: config_use_tidal_potential_forcing_Q1 - logical, pointer :: config_use_tidal_potential_forcing_P1 - character (len=StrKIND), pointer :: config_tidal_potential_reference_time + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< [out] error flag + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + type (MPAS_Time_Type) :: & + refTime ! internal reference time + + integer :: & + iCell, iCon ! loop indices for cell, constituent loops + + real (kind=RKIND) :: & + lat ! latitude temporary + + ! pointers for pool variables type (block_type), pointer :: block_ptr type (mpas_pool_type), pointer :: forcingPool type (mpas_pool_type), pointer :: meshPool type (mpas_pool_type), pointer :: diagnosticsPool - integer, pointer :: nTidalConstituents - integer, pointer :: nCells - real (kind=RKIND), dimension(:), pointer :: tidalConstituentAmplitude - real (kind=RKIND), dimension(:), pointer :: tidalConstituentFreqency - real (kind=RKIND), dimension(:), pointer :: tidalConstituentLoveNumbers - real (kind=RKIND), dimension(:), pointer :: tidalConstituentNodalAmplitude - real (kind=RKIND), dimension(:), pointer :: tidalConstituentAstronomical - real (kind=RKIND), dimension(:), pointer :: tidalConstituentNodalPhase - real (kind=RKIND), dimension(:), pointer :: latCell - real (kind=RKIND), dimension(:), pointer :: eta - real (kind=RKIND), dimension(:,:), pointer :: latitudeFunction - integer, dimension(:), pointer :: tidalConstituentType - type (MPAS_Time_Type) :: refTime - integer :: iCell,iCon - real (kind=RKIND) :: lat + integer, pointer :: nTidalPotentialConstituents + ! End preamble + !----------------------------------------------------------------- + ! Begin code + !*** Set error flags and default values err = 0 + tidalPotentialOff = .true. + betaSelfAttrLoad = 0.0_RKIND + tidalPotRamp = 1.0e20_RKIND + + !*** If tidal potential turned on, set all relevant variables - call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing', config_use_tidal_potential_forcing) - call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_M2', config_use_tidal_potential_forcing_M2) - call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_S2', config_use_tidal_potential_forcing_S2) - call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_N2', config_use_tidal_potential_forcing_N2) - call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_K2', config_use_tidal_potential_forcing_K2) - call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_K1', config_use_tidal_potential_forcing_K1) - call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_O1', config_use_tidal_potential_forcing_O1) - call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_Q1', config_use_tidal_potential_forcing_Q1) - call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_P1', config_use_tidal_potential_forcing_P1) - call mpas_pool_get_config(ocnConfigs, 'config_tidal_potential_reference_time', config_tidal_potential_reference_time) - - block_ptr => domain % blocklist - do while(associated(block_ptr)) - call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool) - call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool) - call mpas_pool_get_dimension(meshPool, 'nCells', nCells) - call mpas_pool_get_array(forcingPool, 'tidalPotentialEta', eta) - - do iCell = 1,nCells - eta(iCell) = 0.0_RKIND - end do - - block_ptr => block_ptr % next - end do - - tidalPotentialOn = .false. if (config_use_tidal_potential_forcing) then - tidalPotentialOn = .true. - else - return - end if + tidalPotentialOff = .false. + + betaSelfAttrLoad = config_self_attraction_and_loading_beta + tidalPotRamp = config_tidal_potential_ramp + + !*** retrieve variables from pool + !*** eventually replace with allocated module private vars + + block_ptr => domain % blocklist + call mpas_pool_get_subpool(block_ptr%structs, 'forcing', & + forcingPool) + + call mpas_pool_get_array(forcingPool, 'tidalPotentialEta', & + tidalPotEta) + + call mpas_pool_get_array(forcingPool, & + 'nTidalPotentialConstituents', & + nTidalPotentialConstituents) + call mpas_pool_get_array(forcingPool, & + 'tidalPotentialConstituentAmplitude', & + tidalConstituentAmplitude) + call mpas_pool_get_array(forcingPool, & + 'tidalPotentialConstituentFrequency', & + tidalConstituentFrequency) + call mpas_pool_get_array(forcingPool, & + 'tidalPotentialConstituentLoveNumbers', & + tidalConstituentLoveNumbers) + call mpas_pool_get_array(forcingPool, & + 'tidalPotentialConstituentNodalAmplitude', & + tidalConstituentNodalAmplitude) + call mpas_pool_get_array(forcingPool, & + 'tidalPotentialConstituentAstronomical', & + tidalConstituentAstronomical) + call mpas_pool_get_array(forcingPool, & + 'tidalPotentialConstituentNodalPhase', & + tidalConstituentNodalPhase) + call mpas_pool_get_array(forcingPool, & + 'tidalPotentialConstituentType', & + tidalConstituentType) + call mpas_pool_get_array(forcingPool, & + 'tidalPotentialLatitudeFunction', & + latitudeFunction) + + call mpas_set_time(refTime, & + dateTimeString=config_tidal_potential_reference_time) + + do iCell = 1,nCellsAll + tidalPotEta(iCell) = 0.0_RKIND + end do - block_ptr => domain % blocklist - do while(associated(block_ptr)) - - call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool) - call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool) - call mpas_pool_get_subpool(block_ptr % structs, 'diagnostics', diagnosticsPool) - call mpas_pool_get_array(forcingPool, 'nTidalPotentialConstituents', nTidalConstituents) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentAmplitude', tidalConstituentAmplitude) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentFrequency', tidalConstituentFreqency) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentLoveNumbers', tidalConstituentLoveNumbers) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentNodalAmplitude', tidalConstituentNodalAmplitude) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentAstronomical', tidalConstituentAstronomical) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentNodalPhase', tidalConstituentNodalPhase) - call mpas_pool_get_array(forcingPool, 'tidalPotentialConstituentType', tidalConstituentType) - call mpas_pool_get_array(forcingPool, 'tidalPotentialLatitudeFunction', latitudeFunction) - call mpas_pool_get_array(meshPool, 'latCell', latCell) - call mpas_pool_get_dimension(meshPool, 'nCells', nCells) - - call mpas_set_time(refTime, dateTimeString=config_tidal_potential_reference_time) - - nTidalConstituents = 0 - if (config_use_tidal_potential_forcing_M2) then - nTidalConstituents = nTidalConstituents + 1 - constituentList(nTidalConstituents)%constituent = 'M2' - end if - - if (config_use_tidal_potential_forcing_S2) then - nTidalConstituents = nTidalConstituents + 1 - constituentList(nTidalConstituents)%constituent = 'S2' - end if + nTidalConstituents = 0 + if (config_use_tidal_potential_forcing_M2) then + nTidalConstituents = nTidalConstituents + 1 + constituentList(nTidalConstituents)%constituent = 'M2' + end if + + if (config_use_tidal_potential_forcing_S2) then + nTidalConstituents = nTidalConstituents + 1 + constituentList(nTidalConstituents)%constituent = 'S2' + end if - if (config_use_tidal_potential_forcing_N2) then - nTidalConstituents = nTidalConstituents + 1 - constituentList(nTidalConstituents)%constituent = 'N2' - end if - - if (config_use_tidal_potential_forcing_K2) then - nTidalConstituents = nTidalConstituents + 1 - constituentList(nTidalConstituents)%constituent = 'K2' - end if - - if (config_use_tidal_potential_forcing_K1) then - nTidalConstituents = nTidalConstituents + 1 - constituentList(nTidalConstituents)%constituent = 'K1' - end if - - if (config_use_tidal_potential_forcing_O1) then - nTidalConstituents = nTidalConstituents + 1 - constituentList(nTidalConstituents)%constituent = 'O1' - end if - - if (config_use_tidal_potential_forcing_Q1) then - nTidalConstituents = nTidalConstituents + 1 - constituentList(nTidalConstituents)%constituent = 'Q1' - end if - - if (config_use_tidal_potential_forcing_P1) then - nTidalConstituents = nTidalConstituents + 1 - constituentList(nTidalConstituents)%constituent = 'P1' - end if - - call tidal_constituent_factors(constituentList,nTidalConstituents,refTime, & - tidalConstituentFreqency, & + if (config_use_tidal_potential_forcing_N2) then + nTidalConstituents = nTidalConstituents + 1 + constituentList(nTidalConstituents)%constituent = 'N2' + end if + + if (config_use_tidal_potential_forcing_K2) then + nTidalConstituents = nTidalConstituents + 1 + constituentList(nTidalConstituents)%constituent = 'K2' + end if + + if (config_use_tidal_potential_forcing_K1) then + nTidalConstituents = nTidalConstituents + 1 + constituentList(nTidalConstituents)%constituent = 'K1' + end if + + if (config_use_tidal_potential_forcing_O1) then + nTidalConstituents = nTidalConstituents + 1 + constituentList(nTidalConstituents)%constituent = 'O1' + end if + + if (config_use_tidal_potential_forcing_Q1) then + nTidalConstituents = nTidalConstituents + 1 + constituentList(nTidalConstituents)%constituent = 'Q1' + end if + + if (config_use_tidal_potential_forcing_P1) then + nTidalConstituents = nTidalConstituents + 1 + constituentList(nTidalConstituents)%constituent = 'P1' + end if + + ! set point value in forcing pool just in case + nTidalPotentialConstituents = nTidalConstituents + + call tidal_constituent_factors(constituentList, & + nTidalConstituents, refTime, & + tidalConstituentFrequency, & tidalConstituentAmplitude, & tidalConstituentLoveNumbers, & tidalConstituentNodalAmplitude, & @@ -426,28 +436,36 @@ subroutine ocn_vel_tidal_potential_init(domain,err)!{{{ tidalConstituentType, & err) - do iCell = 1,nCells - lat = latCell(iCell) - latitudeFunction(iCell,1) = 3.0_RKIND*sin(lat)**2 - 1.0_RKIND - latitudeFunction(iCell,2) = sin(2.0_RKIND*lat) - latitudeFunction(iCell,3) = cos(lat)**2 - end do + do iCell = 1,nCellsAll + lat = latCell(iCell) + latitudeFunction(iCell,1) = 3.0_RKIND*sin(lat)**2 -1.0_RKIND + latitudeFunction(iCell,2) = sin(2.0_RKIND*lat) + latitudeFunction(iCell,3) = cos(lat)**2 + end do - do iCon = 1,nTidalConstituents - call mpas_log_write('Constituent '//constituentList(iCon)%constituent) - call mpas_log_write(' Frequency $r', realArgs=(/ tidalConstituentFreqency(iCon) /)) - call mpas_log_write(' Amplitude $r', realArgs=(/ tidalConstituentAmplitude(iCon) /)) - call mpas_log_write(' LoveNumbers $r', realArgs=(/ tidalConstituentLoveNumbers(iCon) /)) - call mpas_log_write(' NodalAmplitude $r', realArgs=(/ tidalConstituentNodalAmplitude(iCon) /)) - call mpas_log_write(' Astronomical argument $r', realArgs=(/ tidalConstituentAstronomical(iCon) /)) - call mpas_log_write(' NodalPhase $r', realArgs=(/ tidalConstituentNodalPhase(iCon) /)) - call mpas_log_write(' Type $i', intArgs=(/ tidalConstituentType(iCon) /)) - call mpas_log_write(' ') - end do - - block_ptr => block_ptr % next - end do + do iCon = 1,nTidalConstituents + call mpas_log_write( & + 'Constituent '//constituentList(iCon)%constituent) + call mpas_log_write(' Frequency $r', & + realArgs=(/ tidalConstituentFrequency(iCon) /)) + call mpas_log_write(' Amplitude $r', & + realArgs=(/ tidalConstituentAmplitude(iCon) /)) + call mpas_log_write(' LoveNumbers $r', & + realArgs=(/ tidalConstituentLoveNumbers(iCon) /)) + call mpas_log_write(' NodalAmplitude $r', & + realArgs=(/ tidalConstituentNodalAmplitude(iCon) /)) + call mpas_log_write(' Astronomical argument $r', & + realArgs=(/ tidalConstituentAstronomical(iCon) /)) + call mpas_log_write(' NodalPhase $r', & + realArgs=(/ tidalConstituentNodalPhase(iCon) /)) + call mpas_log_write(' Type $i', & + intArgs=(/ tidalConstituentType(iCon) /)) + call mpas_log_write(' ') + end do + end if ! tidal potential on + + !----------------------------------------------------------------- end subroutine ocn_vel_tidal_potential_init!}}} @@ -460,7 +478,7 @@ end subroutine ocn_vel_tidal_potential_init!}}} !> \author Steven Brus !> \date September 2019 !> \details -!> This routine initializes the ampiltude, frequency, love numbers, +!> This routine initializes the amplitude, frequency, love numbers, !> astronomical argument and nodal factors for each tidal constituent !> Nodal factor equations are from: !> "Manual of Harmonic Analysis and Prediction of Tides" @@ -470,7 +488,7 @@ end subroutine ocn_vel_tidal_potential_init!}}} !----------------------------------------------------------------------- subroutine tidal_constituent_factors(constituentList,nTidalConstituents,refTime, & - tidalConstituentFreqency, & + tidalConstituentFrequency, & tidalConstituentAmplitude, & tidalConstituentLoveNumbers, & tidalConstituentNodalAmplitude, & @@ -483,7 +501,7 @@ subroutine tidal_constituent_factors(constituentList,nTidalConstituents,refTime, integer, intent(in) :: nTidalConstituents type (MPAS_Time_Type), intent(in) :: refTime real (kind=RKIND), dimension(:), intent(out) :: tidalConstituentAmplitude - real (kind=RKIND), dimension(:), intent(out) :: tidalConstituentFreqency + real (kind=RKIND), dimension(:), intent(out) :: tidalConstituentFrequency real (kind=RKIND), dimension(:), intent(out) :: tidalConstituentLoveNumbers real (kind=RKIND), dimension(:), intent(out) :: tidalConstituentNodalAmplitude real (kind=RKIND), dimension(:), intent(out) :: tidalConstituentAstronomical @@ -552,7 +570,7 @@ subroutine tidal_constituent_factors(constituentList,nTidalConstituents,refTime, do j = 1,nTidalConstituents if (constituentList(j)%constituent == 'M2') then tidalConstituentAmplitude(j) = 0.242334_RKIND - tidalConstituentFreqency(j) = 1.405189e-4_RKIND + tidalConstituentFrequency(j) = 1.405189e-4_RKIND tidalConstituentLoveNumbers(j) = 0.693_RKIND tidalConstituentAstronomical(j) = 0.0_RKIND tidalConstituentNodalAmplitude(j) = cos(0.5_RKIND*I*deg2rad)**4/0.91544_RKIND @@ -561,7 +579,7 @@ subroutine tidal_constituent_factors(constituentList,nTidalConstituents,refTime, else if (constituentList(j)%constituent == 'S2') then tidalConstituentAmplitude(j) = 0.112743_RKIND - tidalConstituentFreqency(j) = 1.454441e-4_RKIND + tidalConstituentFrequency(j) = 1.454441e-4_RKIND tidalConstituentLoveNumbers(j) = 0.693_RKIND tidalConstituentAstronomical(j) = 0.0_RKIND tidalConstituentNodalAmplitude(j) = 1.0_RKIND @@ -570,7 +588,7 @@ subroutine tidal_constituent_factors(constituentList,nTidalConstituents,refTime, else if (constituentList(j)%constituent == 'N2') then tidalConstituentAmplitude(j) = 0.046397_RKIND - tidalConstituentFreqency(j) = 1.378797e-4_RKIND + tidalConstituentFrequency(j) = 1.378797e-4_RKIND tidalConstituentLoveNumbers(j) = 0.693_RKIND tidalConstituentAstronomical(j) = 0.0_RKIND tidalConstituentNodalAmplitude(j) = cos(0.5_RKIND*I*deg2rad)**4/0.91544_RKIND @@ -579,7 +597,7 @@ subroutine tidal_constituent_factors(constituentList,nTidalConstituents,refTime, else if (constituentList(j)%constituent == 'K2') then tidalConstituentAmplitude(j) = 0.030684_RKIND - tidalConstituentFreqency(j) = 1.458423e-4_RKIND + tidalConstituentFrequency(j) = 1.458423e-4_RKIND tidalConstituentLoveNumbers(j) = 0.693_RKIND tidalConstituentAstronomical(j) = 0.0_RKIND tidalConstituentNodalAmplitude(j) = 0.001_RKIND+sqrt(19.0444_RKIND*sin(I*deg2rad)**4 + & @@ -590,7 +608,7 @@ subroutine tidal_constituent_factors(constituentList,nTidalConstituents,refTime, else if (constituentList(j)%constituent == 'K1') then tidalConstituentAmplitude(j) = 0.141565_RKIND - tidalConstituentFreqency(j) = 0.7292117e-4_RKIND + tidalConstituentFrequency(j) = 0.7292117e-4_RKIND tidalConstituentLoveNumbers(j) = 0.736_RKIND tidalConstituentAstronomical(j) = 0.0_RKIND tidalConstituentNodalAmplitude(j) = sqrt(0.8965_RKIND*sin(2.0_RKIND*I*deg2rad)**2 + & @@ -601,7 +619,7 @@ subroutine tidal_constituent_factors(constituentList,nTidalConstituents,refTime, else if (constituentList(j)%constituent == 'O1') then tidalConstituentAmplitude(j) = 0.100661_RKIND - tidalConstituentFreqency(j) = 0.6759774e-4_RKIND + tidalConstituentFrequency(j) = 0.6759774e-4_RKIND tidalConstituentLoveNumbers(j) = 0.695_RKIND tidalConstituentAstronomical(j) = 0.0_RKIND tidalConstituentNodalAmplitude(j) = sin(I*deg2rad)*cos(0.5_RKIND*I*deg2rad)**2/0.37988_RKIND @@ -610,7 +628,7 @@ subroutine tidal_constituent_factors(constituentList,nTidalConstituents,refTime, else if (constituentList(j)%constituent == 'Q1') then tidalConstituentAmplitude(j) = 0.019273_RKIND - tidalConstituentFreqency(j) = 0.6495854e-4_RKIND + tidalConstituentFrequency(j) = 0.6495854e-4_RKIND tidalConstituentLoveNumbers(j) = 0.695_RKIND tidalConstituentAstronomical(j) = 0.0_RKIND tidalConstituentNodalAmplitude(j) = sin(I*deg2rad)*cos(0.5_RKIND*I*deg2rad)**2/0.37988_RKIND @@ -619,7 +637,7 @@ subroutine tidal_constituent_factors(constituentList,nTidalConstituents,refTime, else if (constituentList(j)%constituent == 'P1') then tidalConstituentAmplitude(j) = 0.046848_RKIND - tidalConstituentFreqency(j) = 0.7252295e-4_RKIND + tidalConstituentFrequency(j) = 0.7252295e-4_RKIND tidalConstituentLoveNumbers(j) = 0.706_RKIND tidalConstituentAstronomical(j) = 0.0_RKIND tidalConstituentNodalAmplitude(j) = 1.0_RKIND diff --git a/src/core_ocean/shared/mpas_ocn_vel_vadv.F b/src/core_ocean/shared/mpas_ocn_vel_vadv.F index d590b1b754..8af4caba81 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_vadv.F +++ b/src/core_ocean/shared/mpas_ocn_vel_vadv.F @@ -22,10 +22,9 @@ module ocn_vel_vadv use mpas_timer - use mpas_derived_types - use mpas_pool_routines use ocn_constants use ocn_config + use ocn_mesh implicit none private @@ -52,7 +51,8 @@ module ocn_vel_vadv ! !-------------------------------------------------------------------- - logical :: velVadvOn + logical :: & + velVadvOff ! on/off switch for vertical advection !*********************************************************************** @@ -72,102 +72,119 @@ module ocn_vel_vadv ! !----------------------------------------------------------------------- - subroutine ocn_vel_vadv_tend(meshPool, normalVelocity, layerThicknessEdge, vertAleTransportTop, tend, err)!{{{ + subroutine ocn_vel_vadv_tend(normalVelocity, layerThickEdge, & + vertAleTransportTop, tend, err)!{{{ !----------------------------------------------------------------- - ! ! input variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(in) :: & - normalVelocity !< Input: Horizontal velocity - real (kind=RKIND), dimension(:,:), intent(in) :: & - layerThicknessEdge,&!< Input: thickness at edge - vertAleTransportTop !< Input: Vertical velocity on top layer - - type (mpas_pool_type), intent(in) :: & - meshPool !< Input: mesh information + normalVelocity, &!< [in] Horizontal velocity + layerThickEdge, &!< [in] Layer thickness at edge + vertAleTransportTop !< [in] Vertical velocity on top layer !----------------------------------------------------------------- - ! ! input/output variables - ! !----------------------------------------------------------------- real (kind=RKIND), dimension(:,:), intent(inout) :: & - tend !< Input/Output: velocity tendency + tend !< [inout] accumulated velocity tendency !----------------------------------------------------------------- - ! ! output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag !----------------------------------------------------------------- - ! ! local variables - ! !----------------------------------------------------------------- - integer :: iEdge, cell1, cell2, k, nEdges - integer, pointer :: nVertLevels - integer, dimension(:), pointer :: nEdgesArray - integer, dimension(:), pointer :: maxLevelEdgeTop - integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask + integer :: & + iEdge, k, &! loop indices for edge, vertical loops + kmax, &! deepest active layer on edge + cell1, cell2 ! neighbor cell indices across edge - real (kind=RKIND) :: vertAleTransportTopEdge - real (kind=RKIND), dimension(:), allocatable :: w_dudzTopEdge + real (kind=RKIND) :: & + wAvg ! ALE transport velocity across top edge - if (.not. velVadvOn) return + real (kind=RKIND), dimension(:,:), allocatable :: & + w_dudzTopEdge ! w*du/dz at top of edge - call mpas_timer_start("vel vadv") + ! End preamble + !----------------------------------------------------------------- + ! Begin code - err = 0 + !*** Set return error code and return early if not turned on + !*** Start relevant timer - call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) - call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) - call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) - call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) - call mpas_pool_get_array(meshPool, 'edgeMask', edgeMask) + err = 0 + if (velVadvOff) return + call mpas_timer_start("vel vadv") - allocate(w_dudzTopEdge(nVertLevels+1)) - w_dudzTopEdge = 0.0_RKIND - nEdges = nEdgesArray( 1 ) + allocate(w_dudzTopEdge(nVertLevels+1,nEdgesAll)) + !$acc enter data create(w_dudzTopEdge) -#ifndef CPRPGI - !$omp parallel + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(cellsOnEdge, maxLevelEdgeTop, w_dudzTopEdge, & + !$acc vertAleTransportTop, normalVelocity, & + !$acc layerThickEdge) & + !$acc private(cell1, cell2, k, kmax, wAvg) + #else !$omp do schedule(runtime) & - !$omp private(iEdge, cell1, cell2, k, vertAleTransportTopEdge) & - !$omp firstprivate(w_dudzTopEdge) -#endif - do iEdge = 1, nEdges - cell1 = cellsOnEdge(1,iEdge) - cell2 = cellsOnEdge(2,iEdge) - - do k = 2, maxLevelEdgeTop(iEdge) - ! Average w from cell center to edge - vertAleTransportTopEdge = 0.5_RKIND*(vertAleTransportTop(k,cell1) + vertAleTransportTop(k,cell2)) - - ! compute dudz at vertical interface with first order derivative. - w_dudzTopEdge(k) = vertAleTransportTopEdge * (normalVelocity(k-1,iEdge)-normalVelocity(k,iEdge)) & - / (0.5_RKIND*(layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge))) - end do - w_dudzTopEdge(maxLevelEdgeTop(iEdge)+1) = 0.0_RKIND - ! Average w*du/dz from vertical interface to vertical middle of cell - do k = 1, maxLevelEdgeTop(iEdge) - - tend(k,iEdge) = tend(k,iEdge) - edgeMask(k, iEdge) * 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1)) - enddo + !$omp private(cell1, cell2, k, kmax, wAvg) + #endif + do iEdge = 1, nEdgesOwned + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + kmax = maxLevelEdgeTop(iEdge) + + w_dudzTopEdge(1,iEdge) = 0.0_RKIND ! flux is zero at top + do k = 2,kmax + + ! Average w from cell center to edge + wAvg = 0.5_RKIND*(vertAleTransportTop(k,cell1) + & + vertAleTransportTop(k,cell2)) + + ! compute dudz at vertical interface with first order derivative. + w_dudzTopEdge(k,iEdge) = wAvg* & + (normalVelocity(k-1,iEdge) - & + normalVelocity(k, iEdge))/ & + (0.5_RKIND*(layerThickEdge(k-1,iEdge) + & + layerThickEdge(k ,iEdge))) + end do + ! transport is zero at bottom + w_dudzTopEdge(kmax+1,iEdge) = 0.0_RKIND + end do + #ifndef MPAS_OPENACC + !$omp end do + #endif + + ! Average w*du/dz from vertical interface to vertical middle of cell + #ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(maxLevelEdgeTop, tend, edgeMask, w_dudzTopEdge)& + !$acc private(k) + #else + !$omp do schedule(runtime) & + !$omp private(k) + #endif + do iEdge = 1, nEdgesOwned + do k = 1, maxLevelEdgeTop(iEdge) + tend(k,iEdge) = tend(k,iEdge) - edgeMask(k,iEdge)* & + 0.5_RKIND*(w_dudzTopEdge(k ,iEdge) + & + w_dudzTopEdge(k+1,iEdge)) + enddo enddo -#ifndef CPRPGI + #ifndef MPAS_OPENACC !$omp end do !$omp end parallel -#endif + #endif + !$acc exit data delete(w_dudzTopEdge) deallocate(w_dudzTopEdge) call mpas_timer_stop("vel vadv") @@ -191,25 +208,30 @@ end subroutine ocn_vel_vadv_tend!}}} subroutine ocn_vel_vadv_init(err)!{{{ - !-------------------------------------------------------------------- - !----------------------------------------------------------------- - ! ! Output variables - ! !----------------------------------------------------------------- - integer, intent(out) :: err !< Output: error flag + integer, intent(out) :: err !< [out] error flag + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + !*** initialize error code and default values err = 0 + velVadvOff = .true. - velVadvOn = .false. + !*** Set values based on input configuration - if (config_vert_coord_movement .ne.'impermeable_interfaces') then - velVadvOn = .true. + if (config_vert_coord_movement == 'impermeable_interfaces') then + velVadvOff = .true. + else + velVadvOff = .false. end if - if ( config_disable_vel_vadv ) velVadvOn = .false. + if (config_disable_vel_vadv ) velVadvOff = .true. !-------------------------------------------------------------------- From 70bcb4956b5035dc358b007e23e38d984597dd90 Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Thu, 17 Dec 2020 16:43:27 -0700 Subject: [PATCH 2/5] Re-add the passing of forcingPool to sfc forcing --- .../shared/mpas_ocn_surface_bulk_forcing.F | 62 ++++++++++++------- .../shared/mpas_ocn_surface_land_ice_fluxes.F | 41 +++++++----- src/core_ocean/shared/mpas_ocn_tendency.F | 51 ++++++--------- 3 files changed, 83 insertions(+), 71 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F b/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F index cbf268a19e..bfacf1d8c1 100644 --- a/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F +++ b/src/core_ocean/shared/mpas_ocn_surface_bulk_forcing.F @@ -142,16 +142,15 @@ end subroutine ocn_surface_bulk_forcing_tracers!}}} !----------------------------------------------------------------------- subroutine ocn_surface_bulk_forcing_vel( & - windStressZonal, windStressMeridional, & - sfcStress, sfcStressMag, err)!{{{ + forcingPool, sfcStress, & + sfcStressMag, err)!{{{ !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- - real (kind=RKIND), dimension(:), intent(in) :: & - windStressZonal, &!< [in] zonal wind stress cell center - windStressMeridional !< [in] meridional wind stress cell center + type (mpas_pool_type), intent(in) :: & + forcingPool !< [in] Forcing structure w/ sfc stresses !----------------------------------------------------------------- ! input/output variables @@ -180,6 +179,10 @@ subroutine ocn_surface_bulk_forcing_vel( & meridWSEdge, &! meridional wind stress avg to edge zonalWSEdge ! zonal wind stress averaged to edge + real (kind=RKIND), dimension(:), pointer :: & + windStressZonal, &! zonal wind stress from forcing + windStressMerid ! meridional wind stress from forcing + ! End preamble !----------------------------------------------------------------- ! Begin code @@ -192,60 +195,71 @@ subroutine ocn_surface_bulk_forcing_vel( & call mpas_timer_start("bulk_ws", .false.) + call mpas_pool_get_array(forcingPool, 'windStressZonal', windStressZonal) + call mpas_pool_get_array(forcingPool, 'windStressMeridional', windStressMerid) + + !*** Transfer data to device + !$acc enter data & + !$acc copyin(windStressZonal, windStressMerid) + !*** Forcing needed over halo regions nEdges = nEdgesHalo(3) nCells = nCellsHalo(2) - !*** Convert coupled climate model wind stress to internal - !*** ocean model wind stress by averaging cell-centered wind + !*** Convert coupled climate model wind stress to internal + !*** ocean model wind stress by averaging cell-centered wind !*** stress to edges and then rotating vectors - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, angleEdge, sfcStress, & - !$acc windStressZonal, windStressMeridional) & - !$acc private(cell1, cell2, zonalWSEdge, meridWSEdge) - #else + !$acc windStressZonal, windStressMerid) & + !$acc private(cell1, cell2, zonalWSEdge, meridWSEdge) +#else !$omp parallel !$omp do schedule(runtime) & - !$omp private(cell1, cell2, zonalWSEdge, meridWSEdge) - #endif + !$omp private(cell1, cell2, zonalWSEdge, meridWSEdge) +#endif do iEdge = 1, nEdges cell1 = cellsOnEdge(1, iEdge) cell2 = cellsOnEdge(2, iEdge) zonalWSEdge = 0.5_RKIND * (windStressZonal(cell1) + & windStressZonal(cell2)) - meridWSEdge = 0.5_RKIND * (windStressMeridional(cell1) + & - windStressMeridional(cell2)) + meridWSEdge = 0.5_RKIND * (windStressMerid(cell1) + & + windStressMerid(cell2)) sfcStress(iEdge) = sfcStress(iEdge) + & cos(angleEdge(iEdge))*zonalWSEdge + & sin(angleEdge(iEdge))*meridWSEdge end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do - #endif +#endif ! Calculate surface flux magnitude at cell centers - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(sfcStressMag, windStressZonal, & - !$acc windStressMeridional) - #else + !$acc windStressMerid) +#else !$omp do schedule(runtime) - #endif +#endif do iCell = 1, nCells sfcStressMag(iCell) = sfcStressMag(iCell) + & sqrt( windStressZonal(iCell)**2 & - + windStressMeridional(iCell)**2 ) + + windStressMerid(iCell)**2 ) end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif + + !*** Transfer any needed data back to host and delete inputs + !$acc exit data & + !$acc delete(windStressZonal, windStressMerid) call mpas_timer_stop("bulk_ws") diff --git a/src/core_ocean/shared/mpas_ocn_surface_land_ice_fluxes.F b/src/core_ocean/shared/mpas_ocn_surface_land_ice_fluxes.F index 069231789f..7c094c83da 100644 --- a/src/core_ocean/shared/mpas_ocn_surface_land_ice_fluxes.F +++ b/src/core_ocean/shared/mpas_ocn_surface_land_ice_fluxes.F @@ -147,16 +147,14 @@ end subroutine ocn_surface_land_ice_fluxes_tracers!}}} ! !----------------------------------------------------------------------- - subroutine ocn_surface_land_ice_fluxes_vel(topDrag, topDragMag, & + subroutine ocn_surface_land_ice_fluxes_vel(diagnosticsPool, & sfcStress, sfcStressMag, err)!{{{ !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- - real (kind=RKIND), dimension(:), intent(in) :: & - topDrag, &!< [in] drag at top of ocean on edge - topDragMag !< [in] magnitude of drag at top of ocean at center + type(mpas_pool_type), intent(in) :: diagnosticsPool !----------------------------------------------------------------- ! input/output variables @@ -179,6 +177,10 @@ subroutine ocn_surface_land_ice_fluxes_vel(topDrag, topDragMag, & integer :: & iEdge, iCell ! loop indices for edge, cell loops + real (kind=RKIND), dimension(:), pointer :: & + topDrag, &!< [in] drag at top of ocean on edge + topDragMag !< [in] magnitude of drag at top of ocean at center + ! End preamble !----------------------------------------------------------------- ! Begin code @@ -191,35 +193,46 @@ subroutine ocn_surface_land_ice_fluxes_vel(topDrag, topDragMag, & call mpas_timer_start("top_drag", .false.) + call mpas_pool_get_array(diagnosticsPool, 'topDrag', topDrag) + call mpas_pool_get_array(diagnosticsPool, 'topDragMagnitude', topDragMag) + + !*** Transfer data to device + !$acc enter data & + !$acc copyin(topDrag, topDragMag) + !*** simply add input values to accumulated stresses - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(sfcStress, topDrag) - #else +#else !$omp parallel !$omp do schedule(runtime) - #endif +#endif do iEdge = 1, nEdgesAll sfcStress(iEdge) = sfcStress(iEdge) + topDrag(iEdge) end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do - #endif +#endif - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(sfcStressMag, topDragMag) - #else +#else !$omp do schedule(runtime) - #endif +#endif do iCell = 1, nCellsAll sfcStressMag(iCell) = sfcStressMag(iCell) + topDragMag(iCell) end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif + + !*** Transfer any needed data back to host and delete inputs + !$acc exit data & + !$acc delete(topDrag, topDragMag) call mpas_timer_stop("top_drag") diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index a9c94f770e..b8062d116d 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -259,10 +259,6 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, & type (mpas_pool_type), pointer :: tracersPool real (kind=RKIND), dimension(:), pointer :: & - windStressZonal, &! zonal wind stress from forcing - windStressMerid, &! meridional wind stress from forcing - topDrag, &! drag at top surface - topDragMag, &! magnitude of drag at top sfcFlxAttCoeff, &! surface flux attenuation coefficient ssh ! sea surface height @@ -355,16 +351,8 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, & thermExpCoeff) call mpas_pool_get_array(diagnosticsPool, 'inSituSalineContractionCoeff', & salineContractCoeff) - call mpas_pool_get_array(diagnosticsPool, 'topDrag', & - topDrag) - call mpas_pool_get_array(diagnosticsPool, 'topDragMagnitude', & - topDragMag) call mpas_pool_get_array(diagnosticsPool, 'wettingVelocity', & wettingVelocity) - call mpas_pool_get_array(forcingPool, 'windStressZonal', & - windStressZonal) - call mpas_pool_get_array(forcingPool, 'windStressMeridional', & - windStressMerid) call mpas_pool_get_array(tendPool, 'normalVelocity', tendVel) @@ -376,7 +364,6 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, & !*** Transfer data to device !$acc enter data & !$acc copyin(tendVel, sfcStress, sfcStressMag, & - !$acc windStressZonal, windStressMerid, topDrag, topDragMag, & !$acc sfcFlxAttCoeff, ssh, layerThickEdge, normalVelocity, & !$acc tangentialVelocity, density, potentialDensity, zMid, & !$acc pressure, layerThickness, circulation, activeTracers, & @@ -387,35 +374,34 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, & !*** Set output variables to zero before accumulating results - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop present(sfcStress, tendVel) private(k) - #else +#else !$omp parallel !$omp do schedule(runtime) private(k) - #endif +#endif do iEdge = 1, nEdgesAll sfcStress(iEdge) = 0.0_RKIND do k=1,nVertLevels tendVel(k,iEdge) = 0.0_RKIND end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do - #endif +#endif - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop present(sfcStressMag) - #else - !$omp end do +#else !$omp do schedule(runtime) - #endif +#endif do iCell = 1, nCellsAll sfcStressMag(iCell) = 0.0_RKIND end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif !*** Return early if all vel tendencies disabled !*** Otherwise, start timer @@ -425,11 +411,11 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, & ! Compute bulk forcing surface stress call ocn_surface_bulk_forcing_vel( & - windStressZonal, windStressMerid, & - sfcStress, sfcStressMag, err) + forcingPool, sfcStress, & + sfcStressMag, err) ! Add top drag to surface stress - call ocn_surface_land_ice_fluxes_vel(topDrag, topDragMag, & + call ocn_surface_land_ice_fluxes_vel(diagnosticsPool, & sfcStress, sfcStressMag, err) ! Add nonlinear Coriolis term and horizontal advection of @@ -469,29 +455,28 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, & ! vertical mixing treated implicitly in a later routine ! adjust total velocity tendency based on wetting/drying - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop collapse(2) & !$acc present(tendVel, wettingVelocity) - #else +#else !$omp parallel !$omp do schedule(runtime) private(k) - #endif +#endif do iEdge = 1, nEdgesAll do k=1,nVertLevels tendVel(k,iEdge) = tendVel(k,iEdge)* & (1.0_RKIND - wettingVelocity(k,iEdge)) end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif !*** Transfer any needed data back to host and delete inputs !$acc exit data & !$acc copyout(tendVel, sfcStress, sfcStressMag) & !$acc delete( & - !$acc windStressZonal, windStressMerid, topDrag, topDragMag, & !$acc sfcFlxAttCoeff, ssh, layerThickEdge, normalVelocity, & !$acc tangentialVelocity, density, potentialDensity, zMid, & !$acc pressure, layerThickness, circulation, activeTracers, & From a824fc3978370064d1e7e5e237a890445497b995 Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Thu, 17 Dec 2020 16:57:17 -0700 Subject: [PATCH 3/5] Move preprocessor statments to first column --- ...pas_ocn_vel_forcing_explicit_bottom_drag.F | 8 ++-- .../mpas_ocn_vel_forcing_surface_stress.F | 8 ++-- .../shared/mpas_ocn_vel_hadv_coriolis.F | 47 +++++++++--------- .../shared/mpas_ocn_vel_hmix_del2.F | 8 ++-- .../shared/mpas_ocn_vel_hmix_del4.F | 32 ++++++------- .../shared/mpas_ocn_vel_hmix_leith.F | 8 ++-- .../shared/mpas_ocn_vel_pressure_grad.F | 48 +++++++++---------- .../shared/mpas_ocn_vel_tidal_potential.F | 12 ++--- src/core_ocean/shared/mpas_ocn_vel_vadv.F | 21 ++++---- 9 files changed, 98 insertions(+), 94 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F b/src/core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F index aed856924f..ce45e34cc0 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F +++ b/src/core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F @@ -127,7 +127,7 @@ subroutine ocn_vel_forcing_explicit_bottom_drag_tend(normVelocity, & ! momentum mixing, and is explicit if both |u| and u are chosen to be at ! time level n. - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, KECell, & !$acc tend, normVelocity, layerThickEdge) & @@ -135,7 +135,7 @@ subroutine ocn_vel_forcing_explicit_bottom_drag_tend(normVelocity, & #else !$omp parallel !$omp do schedule(runtime) private(k, cell1, cell2) - #endif +#endif do iEdge = 1, nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -146,10 +146,10 @@ subroutine ocn_vel_forcing_explicit_bottom_drag_tend(normVelocity, & normVelocity(k,iEdge)/layerThickEdge(k,iEdge) enddo - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif call mpas_timer_stop('vel explicit bottom drag') diff --git a/src/core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F b/src/core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F index 4d30e92241..c1c336a865 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F +++ b/src/core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F @@ -132,7 +132,7 @@ subroutine ocn_vel_forcing_surface_stress_tend(sfcFlxAttCoeff, & if (surfaceStressOff) return call mpas_timer_start('vel surface stress') - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, edgeMask, & !$acc sfcFlxAttCoeff, sfcStress, layerThickEdge, & @@ -146,7 +146,7 @@ subroutine ocn_vel_forcing_surface_stress_tend(sfcFlxAttCoeff, & !$omp private(k, kmax, cell1, cell2, zBot, zTop, & !$omp attCoeff, attDepth, remainingStress, & !$omp transCoeffTop, transCoeffBot, absorb) - #endif +#endif do iEdge = 1, nEdgesOwned zTop = 0.0_RKIND cell1 = cellsOnEdge(1,iEdge) @@ -184,10 +184,10 @@ subroutine ocn_vel_forcing_surface_stress_tend(sfcFlxAttCoeff, & layerThickEdge(kmax,iEdge) end if enddo - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif call mpas_timer_stop('vel surface stress') diff --git a/src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F b/src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F index f0e3029abb..bb06b2a720 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F @@ -140,58 +140,61 @@ subroutine ocn_vel_hadv_coriolis_tend(normRelVortEdge, & qArr(nVertLevels,nEdgesAll)) !$acc enter data create(tmpVorticity, qArr) - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp parallel - #endif +#endif if (usePlanetVorticity) then - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(tmpVorticity, normRelVortEdge, & !$acc normPlanetVortEdge) & !$acc private(k) - #else +#else !$omp do schedule(runtime) & !$omp private(k) - #endif +#endif do iEdge = 1, nEdgesAll do k=1,nVertLevels tmpVorticity(k,iEdge) = normRelVortEdge(k,iEdge) + & normPlanetVortEdge(k,iEdge) end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do - #endif +#endif else - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(tmpVorticity, normRelVortEdge) & !$acc private(k) - #else +#else !$omp do schedule(runtime) & !$omp private(k) - #endif +#endif do iEdge = 1, nEdgesAll do k=1,nVertLevels tmpVorticity(k,iEdge) = normRelVortEdge(k,iEdge) end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do - #endif +#endif endif +#ifndef MPAS_OPENACC + !$omp end parallel +#endif - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(nEdgesOnEdge, edgesOnEdge, weightsOnEdge, & !$acc maxLevelEdgeTop, qArr, tmpVorticity, & !$acc normalVelocity, layerThickEdge) & !$acc private(k, j, eoe, edgeWeight, avgVorticity) - #else +#else !$omp parallel !$omp do schedule(runtime) & !$omp private(k, j, eoe, edgeWeight, avgVorticity) - #endif +#endif do iEdge = 1, nEdgesOwned do k = 1, nVertLevels @@ -213,19 +216,19 @@ subroutine ocn_vel_hadv_coriolis_tend(normRelVortEdge, & end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do - #endif +#endif - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, edgeMask, & !$acc dcEdge, tend, qArr, kineticEnergyCell) & !$acc private(cell1, cell2, invLength, k) - #else +#else !$omp do schedule(runtime) & !$omp private(cell1, cell2, invLength, k) - #endif +#endif do iEdge = 1, nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -238,10 +241,10 @@ subroutine ocn_vel_hadv_coriolis_tend(normRelVortEdge, & - kineticEnergyCell(k,cell1))*invLength) end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif !$acc exit data delete(tmpVorticity, qArr) deallocate(qArr,tmpVorticity) diff --git a/src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F b/src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F index a5c5adbc16..6bf21ace5d 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F @@ -127,7 +127,7 @@ subroutine ocn_vel_hmix_del2_tend(div, relVort, tend, err)!{{{ call mpas_timer_start("vel del2") - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, verticesOnEdge, & !$acc edgeMask, dcEdge, dvEdge, meshScalingDel2, & @@ -139,7 +139,7 @@ subroutine ocn_vel_hmix_del2_tend(div, relVort, tend, err)!{{{ !$omp do schedule(runtime) & !$omp private(k, cell1, cell2, uDiff, vertex1, vertex2, & !$omp dcEdgeInv, dvEdgeInv, visc2) - #endif +#endif do iEdge = 1, nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -168,10 +168,10 @@ subroutine ocn_vel_hmix_del2_tend(div, relVort, tend, err)!{{{ end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif call mpas_timer_stop("vel del2") diff --git a/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F b/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F index f9838d98a7..18ae73eff3 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F @@ -148,7 +148,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ nEdges = nEdgesHalo(2) !Compute Laplacian of velocity - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, verticesOnEdge, maxLevelEdgeTop, & !$acc dvEdge, dcEdge, del2u, div, relVort) & @@ -159,7 +159,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ !$omp do schedule(runtime) & !$omp private(k, cell1, cell2, vertex1, vertex2, & !$omp dcEdgeInv, dvEdgeInv) - #endif +#endif do iEdge = 1, nEdges del2u(:, iEdge) = 0.0_RKIND cell1 = cellsOnEdge(1,iEdge) @@ -178,15 +178,15 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ dvEdgeInv end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif nVertices = nVerticesHalo(1) ! Compute del2 of relative vorticity - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(edgesOnVertex, maxLevelVertexTop, dcEdge, & !$acc areaTriangle, edgeSignOnVertex, del2u, & @@ -195,7 +195,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ #else !$omp parallel !$omp do schedule(runtime) private(i, k, iEdge, areaTriInv) - #endif +#endif do iVertex = 1, nVertices del2relVort(:, iVertex) = 0.0_RKIND areaTriInv = 1.0_RKIND/areaTriangle(iVertex) @@ -208,15 +208,15 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ end do end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif nCells = nCellsHalo(1) ! Compute del2 of divergence - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(nEdgesOnCell, edgesOnCell, maxLevelCell, dvEdge,& !$acc edgeSignOnCell, areaCell, del2u, del2div) & @@ -224,7 +224,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ #else !$omp parallel !$omp do schedule(runtime) private(i, k, iEdge, areaCellInv) - #endif +#endif do iCell = 1, nCells del2div(:, iCell) = 0.0_RKIND areaCellInv = 1.0_RKIND / areaCell(iCell) @@ -237,16 +237,16 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ end do end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif ! Compute final tendency \kappa \nabla^4 u ! as \nabla div(\nabla^2 u) + k \times ! \nabla ( k \cross curl(\nabla^2 u) ) - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, verticesOnEdge, maxLevelEdgeTop, & !$acc dcEdge, dvEdge, meshScalingDel4, edgeMask, & @@ -258,7 +258,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ !$omp do schedule(runtime) & !$omp private(k, cell1, cell2, vertex1, vertex2, & !$omp dcEdgeInv, dvEdgeInv, visc4, uDiff) - #endif +#endif do iEdge = 1, nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -279,10 +279,10 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ edgeMask(k,iEdge)*uDiff*visc4 end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif !$acc exit data delete(del2u, del2div, del2relVort) deallocate(del2u, & diff --git a/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F b/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F index 0d37a839cf..65d7c2653f 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F @@ -137,7 +137,7 @@ subroutine ocn_vel_hmix_leith_tend(div, relVort, tend, err)!{{{ if (hmixLeithOff) return call mpas_timer_start("vel leith") - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, verticesOnEdge, maxLevelEdgeTop, & !$acc dcEdge, dvEdge, meshScalingDel2, div, relVort, & @@ -149,7 +149,7 @@ subroutine ocn_vel_hmix_leith_tend(div, relVort, tend, err)!{{{ !$omp do schedule(runtime) & !$omp private(k, cell1, cell2, vertex1, vertex2, dcEdgeInv, & !$omp dvEdgeInv, uDiff, visc2, visc2tmp) - #endif +#endif do iEdge = 1, nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -186,10 +186,10 @@ subroutine ocn_vel_hmix_leith_tend(div, relVort, tend, err)!{{{ end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif call mpas_timer_stop("vel leith") diff --git a/src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F b/src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F index 095da74024..d41083bc39 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F +++ b/src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F @@ -179,7 +179,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & ! pressure for sea surface height = - g grad ssh - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, & !$acc tend, edgeMask, ssh) & @@ -188,7 +188,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, k, kMax) - #endif +#endif do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -201,17 +201,17 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & (ssh(cell2) - ssh(cell1)) end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif case (pGradTypePZmid) ! pressure for generalized coordinates ! -1/density_0 (grad p_k + density g grad z_k^{mid}) - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, zMid, & !$acc tend, edgeMask, pressure, density) & @@ -220,7 +220,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, k, kMax) - #endif +#endif do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -237,17 +237,17 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & (zMid(k,cell2)-zMid(k,cell1))) end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif case (pGradTypeMontPot) ! For pure isopycnal coordinates, this is just grad(M), ! the gradient of Montgomery Potential - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, & !$acc tend, edgeMask, montgomeryPotential) & @@ -256,7 +256,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, k, kMax) - #endif +#endif do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -270,10 +270,10 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & montgomeryPotential(k,cell1))) end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif case (pGradTypeMontPotDens) @@ -284,7 +284,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & ! Where rho is the potential density. ! See Bleck (2002) equation 1, and last equation in Appendix A. - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, dcEdge, maxLevelEdgeTop, & !$acc edgeMask, tend, montgomeryPotential, & @@ -294,7 +294,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, k, kMax) - #endif +#endif do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -312,17 +312,17 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & - 1.0_RKIND/potentialDensity(k,cell1))) end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif case (pGradTypeJacobDens) allocate(JacobianDxDs(nVertLevels,nEdgesAll)) !$acc enter data create (JacobianDxDs) - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, & !$acc edgeMask, pressure, density, zMid, tend, & @@ -334,7 +334,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$omp do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, k, kMax, pGrad, & !$omp Area, zStar, zC, zGamma, rhoL, rhoR) - #endif +#endif do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -402,10 +402,10 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif !$acc exit data delete(JacobianDxDs) deallocate(JacobianDxDs) @@ -418,7 +418,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$acc enter data create(JacobianDxDs, & !$acc JacobianTz, JacobianSz) - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, & !$acc zMid, tracers, edgeMask, pressure, & @@ -432,7 +432,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$omp do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, k, kMax, alpha, beta, & !$omp TL, TR, SL, SR, Area, zStar, zC, zGamma, pGrad) - #endif +#endif do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -529,10 +529,10 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & end do end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif !$acc exit data delete (JacobianDxDs, JacobianTz, JacobianSz) deallocate(JacobianDxDs, JacobianTz, JacobianSz) diff --git a/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F b/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F index 67387101a4..bbac972c58 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F +++ b/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F @@ -152,15 +152,15 @@ subroutine ocn_vel_tidal_potential_tend(ssh, tend, err)!{{{ err = 0 if (tidalPotentialOff) return - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, edgeMask, & !$acc tidalPotEta, ssh, tend) & !$acc private(cell1, cell2, invdcEdge, potentialGrad, k, kMax) - #else +#else !$omp parallel do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, potentialGrad, k, kMax) - #endif +#endif do iEdge=1,nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -177,9 +177,9 @@ subroutine ocn_vel_tidal_potential_tend(ssh, tend, err)!{{{ edgeMask(k,iEdge)*potentialGrad end do end do - #ifndef MPAS_OPENACC - !$omp end do - #endif +#ifndef MPAS_OPENACC + !$omp end parallel do +#endif end subroutine ocn_vel_tidal_potential_tend!}}} diff --git a/src/core_ocean/shared/mpas_ocn_vel_vadv.F b/src/core_ocean/shared/mpas_ocn_vel_vadv.F index 8af4caba81..0c0170e03b 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_vadv.F +++ b/src/core_ocean/shared/mpas_ocn_vel_vadv.F @@ -127,16 +127,17 @@ subroutine ocn_vel_vadv_tend(normalVelocity, layerThickEdge, & allocate(w_dudzTopEdge(nVertLevels+1,nEdgesAll)) !$acc enter data create(w_dudzTopEdge) - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(cellsOnEdge, maxLevelEdgeTop, w_dudzTopEdge, & !$acc vertAleTransportTop, normalVelocity, & !$acc layerThickEdge) & !$acc private(cell1, cell2, k, kmax, wAvg) - #else +#else + !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1, cell2, k, kmax, wAvg) - #endif +#endif do iEdge = 1, nEdgesOwned cell1 = cellsOnEdge(1,iEdge) cell2 = cellsOnEdge(2,iEdge) @@ -159,19 +160,19 @@ subroutine ocn_vel_vadv_tend(normalVelocity, layerThickEdge, & ! transport is zero at bottom w_dudzTopEdge(kmax+1,iEdge) = 0.0_RKIND end do - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do - #endif +#endif ! Average w*du/dz from vertical interface to vertical middle of cell - #ifdef MPAS_OPENACC +#ifdef MPAS_OPENACC !$acc parallel loop & !$acc present(maxLevelEdgeTop, tend, edgeMask, w_dudzTopEdge)& !$acc private(k) - #else +#else !$omp do schedule(runtime) & !$omp private(k) - #endif +#endif do iEdge = 1, nEdgesOwned do k = 1, maxLevelEdgeTop(iEdge) tend(k,iEdge) = tend(k,iEdge) - edgeMask(k,iEdge)* & @@ -179,10 +180,10 @@ subroutine ocn_vel_vadv_tend(normalVelocity, layerThickEdge, & w_dudzTopEdge(k+1,iEdge)) enddo enddo - #ifndef MPAS_OPENACC +#ifndef MPAS_OPENACC !$omp end do !$omp end parallel - #endif +#endif !$acc exit data delete(w_dudzTopEdge) deallocate(w_dudzTopEdge) From 90115d002eb49bbf077a4613802f8d3bfaab6a00 Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Thu, 17 Dec 2020 16:59:08 -0700 Subject: [PATCH 4/5] Move preprocessor else statments to first column --- .../mpas_ocn_vel_forcing_explicit_bottom_drag.F | 2 +- .../shared/mpas_ocn_vel_forcing_surface_stress.F | 2 +- src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F | 2 +- src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F | 8 ++++---- src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F | 2 +- src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F | 12 ++++++------ 6 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F b/src/core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F index ce45e34cc0..3e8cf7b71a 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F +++ b/src/core_ocean/shared/mpas_ocn_vel_forcing_explicit_bottom_drag.F @@ -132,7 +132,7 @@ subroutine ocn_vel_forcing_explicit_bottom_drag_tend(normVelocity, & !$acc present(cellsOnEdge, maxLevelEdgeTop, KECell, & !$acc tend, normVelocity, layerThickEdge) & !$acc private(k, cell1, cell2) - #else +#else !$omp parallel !$omp do schedule(runtime) private(k, cell1, cell2) #endif diff --git a/src/core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F b/src/core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F index c1c336a865..eb1eb09e6c 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F +++ b/src/core_ocean/shared/mpas_ocn_vel_forcing_surface_stress.F @@ -140,7 +140,7 @@ subroutine ocn_vel_forcing_surface_stress_tend(sfcFlxAttCoeff, & !$acc private(k, kmax, cell1, cell2, zBot, zTop, & !$acc attCoeff, attDepth, remainingStress, & !$acc transCoeffTop, transCoeffBot, absorb) - #else +#else !$omp parallel !$omp do schedule(runtime) & !$omp private(k, kmax, cell1, cell2, zBot, zTop, & diff --git a/src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F b/src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F index 6bf21ace5d..da72525eee 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hmix_del2.F @@ -134,7 +134,7 @@ subroutine ocn_vel_hmix_del2_tend(div, relVort, tend, err)!{{{ !$acc div, relVort, tend) & !$acc private(k, cell1, cell2, uDiff, vertex1, vertex2, & !$acc dcEdgeInv, dvEdgeInv, visc2) - #else +#else !$omp parallel !$omp do schedule(runtime) & !$omp private(k, cell1, cell2, uDiff, vertex1, vertex2, & diff --git a/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F b/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F index 18ae73eff3..94ed413f50 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hmix_del4.F @@ -154,7 +154,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ !$acc dvEdge, dcEdge, del2u, div, relVort) & !$acc private(k, cell1, cell2, vertex1, vertex2, & !$acc dcEdgeInv, dvEdgeInv) - #else +#else !$omp parallel !$omp do schedule(runtime) & !$omp private(k, cell1, cell2, vertex1, vertex2, & @@ -192,7 +192,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ !$acc areaTriangle, edgeSignOnVertex, del2u, & !$acc del2relVort) & !$acc private(i, k, iEdge, areaTriInv) - #else +#else !$omp parallel !$omp do schedule(runtime) private(i, k, iEdge, areaTriInv) #endif @@ -221,7 +221,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ !$acc present(nEdgesOnCell, edgesOnCell, maxLevelCell, dvEdge,& !$acc edgeSignOnCell, areaCell, del2u, del2div) & !$acc private(i, k, iEdge, areaCellInv) - #else +#else !$omp parallel !$omp do schedule(runtime) private(i, k, iEdge, areaCellInv) #endif @@ -253,7 +253,7 @@ subroutine ocn_vel_hmix_del4_tend(div, relVort, tend, err)!{{{ !$acc del2div, del2relVort, tend) & !$acc private(k, cell1, cell2, vertex1, vertex2, & !$acc dcEdgeInv, dvEdgeInv, visc4, uDiff) - #else +#else !$omp parallel !$omp do schedule(runtime) & !$omp private(k, cell1, cell2, vertex1, vertex2, & diff --git a/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F b/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F index 65d7c2653f..467ab15f19 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hmix_leith.F @@ -144,7 +144,7 @@ subroutine ocn_vel_hmix_leith_tend(div, relVort, tend, err)!{{{ !$acc tend, edgeMask) & !$acc private(k, cell1, cell2, vertex1, vertex2, dcEdgeInv, & !$acc dvEdgeInv, uDiff, visc2, visc2tmp) - #else +#else !$omp parallel !$omp do schedule(runtime) & !$omp private(k, cell1, cell2, vertex1, vertex2, dcEdgeInv, & diff --git a/src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F b/src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F index d41083bc39..5f3287a25c 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F +++ b/src/core_ocean/shared/mpas_ocn_vel_pressure_grad.F @@ -184,7 +184,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, & !$acc tend, edgeMask, ssh) & !$acc private(cell1, cell2, invdcEdge, k, kMax) - #else +#else !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, k, kMax) @@ -216,7 +216,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, zMid, & !$acc tend, edgeMask, pressure, density) & !$acc private(cell1, cell2, invdcEdge, k, kMax) - #else +#else !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, k, kMax) @@ -252,7 +252,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$acc present(cellsOnEdge, maxLevelEdgeTop, dcEdge, & !$acc tend, edgeMask, montgomeryPotential) & !$acc private(cell1, cell2, invdcEdge, k, kMax) - #else +#else !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, k, kMax) @@ -290,7 +290,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$acc edgeMask, tend, montgomeryPotential, & !$acc pressure, potentialDensity) & !$acc private(cell1, cell2, invdcEdge, k, kMax) - #else +#else !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, k, kMax) @@ -329,7 +329,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$acc JacobianDxDs) & !$acc private(cell1, cell2, invdcEdge, k, kMax, pGrad, & !$acc Area, zStar, zC, zGamma, rhoL, rhoR) - #else +#else !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, k, kMax, pGrad, & @@ -427,7 +427,7 @@ subroutine ocn_vel_pressure_grad_tend(ssh, pressure, & !$acc JacobianTz, JacobianSz, JacobianDxDs) & !$acc private(cell1, cell2, invdcEdge, k,kMax, alpha, beta,& !$acc TL, TR, SL, SR, Area, zStar, zC, zGamma, pGrad) - #else +#else !$omp parallel !$omp do schedule(runtime) & !$omp private(cell1, cell2, invdcEdge, k, kMax, alpha, beta, & From a30fd34006d149d41321e0c0927ff8dd25e18768 Mon Sep 17 00:00:00 2001 From: Matt Turner Date: Thu, 14 Jan 2021 11:09:15 -0700 Subject: [PATCH 5/5] Updates for recent SI merge --- src/core_ocean/mode_forward/mpas_ocn_time_integration_si.F | 2 +- src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/core_ocean/mode_forward/mpas_ocn_time_integration_si.F b/src/core_ocean/mode_forward/mpas_ocn_time_integration_si.F index c5f0270a3e..b08f926e9e 100644 --- a/src/core_ocean/mode_forward/mpas_ocn_time_integration_si.F +++ b/src/core_ocean/mode_forward/mpas_ocn_time_integration_si.F @@ -510,7 +510,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ call mpas_pool_get_array(diagnosticsPool, 'layerThicknessEdge', layerThicknessEdge) - call ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshPool, scratchPool, stage1_tend_time, dt) + call ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, stage1_tend_time, dt) block => block % next end do diff --git a/src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F b/src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F index bb06b2a720..6fc5d147a4 100644 --- a/src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F +++ b/src/core_ocean/shared/mpas_ocn_vel_hadv_coriolis.F @@ -303,13 +303,16 @@ subroutine ocn_vel_hadv_coriolis_init(err)!{{{ usePlanetVorticity = .false. case ('unsplit_explicit') - ! For split explicit, Coriolis tendency uses eta/h because + ! For unsplit explicit, Coriolis tendency uses eta/h because ! the Coriolis term is added separately to the momentum ! tendencies. usePlanetVorticity = .false. case ('semi_implicit') - usePlanetaryVorticity = .false. + ! For semi-implicit, Coriolis tendency uses eta/h because + ! the Coriolis term is added separately to the momentum + ! tendencies. + usePlanetVorticity = .false. end select