diff --git a/cime_config/tests.py b/cime_config/tests.py index 3ebac841239f..744947b26867 100644 --- a/cime_config/tests.py +++ b/cime_config/tests.py @@ -298,7 +298,9 @@ "e3sm_ocnice_stealth_features" : { "tests" : ( "SMS_D_Ld1.T62_oQU240wLI.GMPAS-IAF-PISMF.mpaso-impl_top_drag", + "SMS_D_Ld1.T62_oQU240wLI.GMPAS-IAF-PISMF.mpaso-vertical_remap", "SMS_D_Ld1.T62_oQU240.GMPAS-IAF.mpaso-harmonic_mean_drag", + "SMS_D_Ld1.T62_oQU240.GMPAS-IAF.mpaso-vertical_remap", "SMS_D_Ld1.T62_oQU240.GMPAS-IAF.mpaso-upwind_advection", "SMS_D_Ld1.T62_oQU240.GMPAS-IAF.mpaso-freshwater_tracers", "ERS_Ld5_D.T62_oQU240.GMPAS-IAF.mpaso-conservation_check", diff --git a/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/vertical_remap/README b/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/vertical_remap/README new file mode 100644 index 000000000000..74a87a24ec97 --- /dev/null +++ b/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/vertical_remap/README @@ -0,0 +1 @@ +Test vertical Lagrangian-remapping capability with a remapping interval of every third timestep. diff --git a/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/vertical_remap/user_nl_mpaso b/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/vertical_remap/user_nl_mpaso new file mode 100644 index 000000000000..e9b0e07ee6d7 --- /dev/null +++ b/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/vertical_remap/user_nl_mpaso @@ -0,0 +1,2 @@ +config_vert_advection_method = 'remap' +config_vert_remap_interval = 3 diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index dd218f789e4f..69a7cfe9d6b1 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -3095,7 +3095,7 @@ packages="forwardMode;analysisMode" /> + domain % blocklist do while (associated(block)) - call ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, rk_substep_weights(rk_step), domain % dminfo, err ) - call ocn_time_integrator_rk4_compute_thick_tends( block, dt, rk_substep_weights(rk_step), err ) + call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) + allocate(div_hu(nVertLevels, nCells)) + call ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, rk_substep_weights(rk_step), domain % dminfo, div_hu, err ) + call ocn_time_integrator_rk4_compute_thick_tends( block, dt, div_hu, rk_substep_weights(rk_step), err ) + deallocate(div_hu) block => block % next end do @@ -525,7 +531,11 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ block => domain % blocklist do while (associated(block)) - call ocn_time_integrator_rk4_compute_tracer_tends( block, dt, rk_substep_weights(rk_step), err ) + call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) + allocate(div_hu(nVertLevels, nCells)) + div_hu = 0.0_RKIND ! This is only needed for vertical remap and we do not support it for RK4 + call ocn_time_integrator_rk4_compute_tracer_tends( block, dt, div_hu, rk_substep_weights(rk_step), err ) + deallocate(div_hu) block => block % next end do @@ -869,15 +879,20 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ end subroutine ocn_time_integrator_rk4!}}} subroutine ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, & - rkSubstepWeight, dminfo, err)!{{{ + rkSubstepWeight, dminfo, div_hu, err)!{{{ type (domain_type), intent(inout) :: domain !< Input/Output: domain information type (block_type), intent(in) :: block real (kind=RKIND), intent(in) :: dt real (kind=RKIND), intent(in) :: rkSubstepWeight + real (kind=RKIND), intent(out), dimension(:,:) :: div_hu integer, intent(out) :: err type (dm_info), intent(in) :: dminfo + real (kind=RKIND), dimension(:), allocatable :: & + projectedSSH ! SSH used to determine layerInterfaceVelocity + real (kind=RKIND), dimension(:,:), allocatable:: layerInterfaceVelocity + type (mpas_pool_type), pointer :: meshPool, verticalMeshPool type (mpas_pool_type), pointer :: statePool, forcingPool type (mpas_pool_type), pointer :: scratchPool, tendPool, provisStatePool @@ -911,28 +926,46 @@ subroutine ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, & call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(verticalMeshPool, & - layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) + ! The diagnostic variable layerThicknessTarget is updated + allocate(layerInterfaceVelocity(nVertLevels, nCellsAll), & + projectedSSH(nCellsAll)) +#ifdef MPAS_OPENACC + !$acc enter data create(projectedSSH, div_hu) +#endif + call ocn_diagnostic_solve_divergence_hu(sshCur, layerThickEdgeFlux, & + normalVelocityProvis, dt, div_hu, projectedSSH) + if ( TRIM(config_vert_advection_method) == 'remap' ) then + vertAleTransportTop = 0.0_RKIND else - call ocn_vert_transport_velocity_top(verticalMeshPool, & - layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err) + if (associated(highFreqThicknessProvis)) then + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH, layerThicknessCur, rkSubstepWeight, & + layerInterfaceVelocity, err, highFreqThicknessProvis) + else + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH, layerThicknessCur, & + rkSubstepWeight, layerInterfaceVelocity, err) + endif + call ocn_vert_transport_velocity_top(div_hu, layerInterfaceVelocity, & + vertAleTransportTop, err) endif call ocn_tend_vel(domain, tendPool, provisStatePool, forcingPool, 1, dminfo, dt) end subroutine ocn_time_integrator_rk4_compute_vel_tends!}}} - subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeight, err)!{{{ + subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, div_hu, rkSubstepWeight, err)!{{{ type (block_type), intent(in) :: block real (kind=RKIND), intent(in) :: dt + real (kind=RKIND), intent(in), dimension(:,:) :: div_hu real (kind=RKIND), intent(in) :: rkSubstepWeight integer, intent(out) :: err + real (kind=RKIND), dimension(:), allocatable :: & + projectedSSH_tr ! SSH used to determine layerInterfaceVelocity + real (kind=RKIND), dimension(:,:), allocatable:: & + div_hu_tr, layerInterfaceVelocity + type (mpas_pool_type), pointer :: meshPool, verticalMeshPool type (mpas_pool_type), pointer :: statePool, forcingPool type (mpas_pool_type), pointer :: scratchPool, tendPool, provisStatePool @@ -942,7 +975,7 @@ subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeigh real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur, normalVelocityCur real (kind=RKIND), dimension(:, :), pointer :: normalVelocityProvis, highFreqThicknessProvis - integer :: iEdge + integer :: iEdge, k, iCell integer, pointer :: nEdges, nVertLevels logical, pointer :: config_filter_btr_mode @@ -984,26 +1017,48 @@ subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeigh call ocn_MLE_add_to_transport_vel(normalTransportVelocity, nEdges) ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) + ! The diagnostic variable layerThicknessTarget is updated + allocate(layerInterfaceVelocity(nVertLevels, nCellsAll), & + div_hu_tr(nVertLevels, nCellsAll), & + projectedSSH_tr(nCellsAll)) + call ocn_diagnostic_solve_divergence_hu(sshCur, layerThickEdgeFlux, & + normalTransportVelocity, dt, div_hu_tr, projectedSSH_tr) + if ( TRIM(config_vert_advection_method) == 'remap' ) then + ! This is the layerInterfaceVelocity implicit in w=0 during Stage 1 + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k) + do iCell = 1, nCellsAll + do k = minLevelCell(iCell), maxLevelCell(iCell) + layerInterfaceVelocity(k,iCell) = div_hu(k,iCell) + enddo + enddo + !$omp end do + !$omp end parallel else - call ocn_vert_transport_velocity_top(verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err) + if (associated(highFreqThicknessProvis)) then + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH_tr, layerThicknessCur, dt, & + layerInterfaceVelocity, err, highFreqThicknessProvis) + else + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH_tr, layerThicknessCur, & + dt, layerInterfaceVelocity, err) + endif endif + call ocn_vert_transport_velocity_top(div_hu_tr, layerInterfaceVelocity, & + vertAleTransportTop, err) + deallocate(layerInterfaceVelocity, div_hu_tr, projectedSSH_tr) call ocn_tend_thick(tendPool, forcingPool) end subroutine ocn_time_integrator_rk4_compute_thick_tends!}}} - subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, rkSubstepWeight, err)!{{{ + subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, div_hu, rkSubstepWeight, err)!{{{ type (block_type), intent(in) :: block real (kind=RKIND), intent(in) :: dt + real (kind=RKIND), intent(in), dimension(:,:) :: div_hu real (kind=RKIND), intent(in) :: rkSubstepWeight integer, intent(out) :: err @@ -1012,13 +1067,15 @@ subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, rkSubstepWeig type (mpas_pool_type), pointer :: scratchPool, tendPool, provisStatePool type (mpas_pool_type), pointer :: swForcingPool, tracersPool + real (kind=RKIND), dimension(:), allocatable :: projectedSSH_tr + real (kind=RKIND), dimension(:, :), allocatable :: div_hu_tr, layerInterfaceVelocity real (kind=RKIND), dimension(:), pointer :: sshCur real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur, normalVelocityCur real (kind=RKIND), dimension(:, :), pointer :: normalVelocityProvis, highFreqThicknessProvis logical, pointer :: config_filter_btr_mode - integer :: iEdge + integer :: k, iCell, iEdge integer, pointer :: nEdges, nVertLevels err = 0 @@ -1060,17 +1117,37 @@ subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, rkSubstepWeig call ocn_MLE_add_to_transport_vel(normalTransportVelocity, nEdges) ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) + ! The diagnostic variable layerThicknessTarget is updated + allocate(div_hu_tr(nVertLevels, nCellsAll), & + projectedSSH_tr(nCellsAll)) + call ocn_diagnostic_solve_divergence_hu(sshCur, layerThickEdgeFlux, & + normalTransportVelocity, dt, div_hu_tr, projectedSSH_tr) + if ( TRIM(config_vert_advection_method) == 'remap' ) then + ! This is the layerInterfaceVelocity implicit in w=0 during Stage 1 + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k) + do iCell = 1, nCellsAll + do k = minLevelCell(iCell), maxLevelCell(iCell) + layerInterfaceVelocity(k,iCell) = div_hu(k,iCell) + enddo + enddo + !$omp end do + !$omp end parallel else - call ocn_vert_transport_velocity_top(verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkSubstepWeight, & - vertAleTransportTop, err) + if (associated(highFreqThicknessProvis)) then + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH_tr, layerThicknessCur, dt, & + layerInterfaceVelocity, err, highFreqThicknessProvis) + else + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH_tr, layerThicknessCur, & + dt, layerInterfaceVelocity, err) + endif endif + call ocn_vert_transport_velocity_top(div_hu_tr, layerInterfaceVelocity, & + vertAleTransportTop, err) + deallocate(layerInterfaceVelocity, div_hu_tr, projectedSSH_tr) if (config_filter_btr_mode) then call ocn_filter_btr_mode_tend_vel(tendPool, provisStatePool, meshPool, 1) @@ -1098,6 +1175,10 @@ subroutine ocn_time_integrator_rk4_compute_tends(domain, block, dt, rkWeight, dm real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur, normalVelocityCur real (kind=RKIND), dimension(:, :), pointer :: normalVelocityProvis, highFreqThicknessProvis + real (kind=RKIND), dimension(:), allocatable :: projectedSSH, projectedSSH_tr + real (kind=RKIND), dimension(:, :), allocatable :: div_hu, div_hu_tr, layerInterfaceVelocity + + integer :: k, iCell logical, pointer :: config_filter_btr_mode err = 0 @@ -1123,31 +1204,63 @@ subroutine ocn_time_integrator_rk4_compute_tends(domain, block, dt, rkWeight, dm call mpas_pool_get_array(provisStatePool, 'highFreqThickness', highFreqThicknessProvis, 1) ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity. - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(verticalMeshPool, & - layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) + ! The diagnostic variable layerThicknessTarget is updated + allocate(layerInterfaceVelocity(nVertLevels, nCellsAll), & + div_hu(nVertLevels, nCellsAll), & + projectedSSH(nCellsAll)) +#ifdef MPAS_OPENACC + !$acc enter data create(projectedSSH, div_hu) +#endif + call ocn_diagnostic_solve_divergence_hu(sshCur, layerThickEdgeFlux, & + normalVelocityProvis, dt, div_hu, projectedSSH) + if ( TRIM(config_vert_advection_method) == 'remap' ) then + vertAleTransportTop = 0.0_RKIND else - call ocn_vert_transport_velocity_top(verticalMeshPool, & - layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, & - sshCur, rkWeight, & - vertAleTransportTop, err) + if (associated(highFreqThicknessProvis)) then + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH, layerThicknessCur, dt, & + layerInterfaceVelocity, err, highFreqThicknessProvis) + else + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH, layerThicknessCur, & + dt, layerInterfaceVelocity, err) + endif + call ocn_vert_transport_velocity_top(div_hu, layerInterfaceVelocity, & + vertAleTransportTop, err) endif call ocn_tend_vel(domain, tendPool, provisStatePool, forcingPool, 1, dminfo, dt) - if (associated(highFreqThicknessProvis)) then - call ocn_vert_transport_velocity_top(verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkWeight, & - vertAleTransportTop, err, highFreqThicknessProvis) + allocate(div_hu_tr(nVertLevels, nCellsAll), & + projectedSSH_tr(nCellsAll)) + call ocn_diagnostic_solve_divergence_hu(sshCur, layerThickEdgeFlux, & + normalTransportVelocity, dt, div_hu_tr, projectedSSH_tr) + if ( TRIM(config_vert_advection_method) == 'remap' ) then + ! This is the layerInterfaceVelocity implicit in w=0 during Stage 1 + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k) + do iCell = 1, nCellsAll + do k = minLevelCell(iCell), maxLevelCell(iCell) + layerInterfaceVelocity(k,iCell) = div_hu(k,iCell) + enddo + enddo + !$omp end do + !$omp end parallel else - call ocn_vert_transport_velocity_top(verticalMeshPool, & - layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, & - sshCur, rkWeight, & - vertAleTransportTop, err) + if (associated(highFreqThicknessProvis)) then + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH_tr, layerThicknessCur, dt, & + layerInterfaceVelocity, err, highFreqThicknessProvis) + else + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH_tr, layerThicknessCur, & + dt, layerInterfaceVelocity, err) + endif endif + call ocn_vert_transport_velocity_top(div_hu_tr, layerInterfaceVelocity, & + vertAleTransportTop, err) + deallocate(layerInterfaceVelocity, div_hu_tr, projectedSSH_tr) call ocn_tend_thick(tendPool, forcingPool) diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F index 77bf1331d3d1..e8830ceb6418 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F @@ -298,10 +298,11 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ fluxAx real (kind=RKIND), dimension(:,:), allocatable:: & - uTemp + div_hu, div_hu_tr, uTemp, layerInterfaceVelocity real (kind=RKIND), dimension(:), allocatable :: & btrvel_temp, & + projectedSSH, projectedSSH_tr, & ! SSH used to determine layerInterfaceVelocity bottomDepthEdge ! State Array Pointers @@ -384,7 +385,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ ! Remap variables real (kind=RKIND), dimension(:,:), pointer :: & - layerThicknessLagNew + layerThicknessLagCur, layerThicknessLagNew integer :: & siLargeIter ! Index of the large btr system iteration @@ -477,6 +478,9 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ call mpas_pool_get_array(statePool, 'layerThickness', & layerThicknessNew, 2) + call mpas_pool_get_array(statePool, 'layerThicknessLag', layerThicknessLagCur, 1) + call mpas_pool_get_array(statePool, 'layerThicknessLag', layerThicknessLagNew, 2) + call mpas_pool_get_array(statePool, 'highFreqThickness', & highFreqThicknessCur, 1) call mpas_pool_get_array(statePool, 'highFreqThickness', & @@ -540,6 +544,29 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ #endif + if ( configVertAdvMethod == vertAdvRemap ) then + if ( associated(layerThicknessLagCur) ) then +#ifdef MPAS_OPENACC + !$acc loop collapse(2) +#else + !$omp do schedule(runtime) private(k) +#endif + do iCell= 1,nCellsAll + do k = 1,nVertLevels + ! Store the current, Lagrangian location for computation of the + ! Eulerian vertical velocity induced by remapping + layerThicknessLagCur(k,iCell) = layerThicknessCur(k,iCell) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do +#endif + endif + ! Perform the remapping + ! The updated layerThickness goes in layerThicknessCur which is assigned to layerThicknessNew below + call ocn_remap_vert_state(block, err, 1) + end if + ! Initialize * variables that are used to compute baroclinic ! tendencies below. @@ -828,16 +855,28 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ ! compute vertAleTransportTop. Use u (rather than & ! normalTransportVelocity) for momentum advection. ! Use the most recent time level available. - - if (associated(highFreqThicknessNew)) then - call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalVelocityCur, sshCur, dt, & - vertAleTransportTop, err, highFreqThicknessNew) + ! The diagnostic variable layerThicknessTarget is updated + allocate(layerInterfaceVelocity(nVertLevels, nCellsAll), & + div_hu(nVertLevels, nCellsAll), & + projectedSSH(nCellsAll)) +#ifdef MPAS_OPENACC + !$acc enter data create(projectedSSH, div_hu) +#endif + call ocn_diagnostic_solve_divergence_hu(sshCur, layerThickEdgeFlux, & + normalVelocityCur, dt, div_hu, projectedSSH) + if ( configVertAdvMethod == vertAdvRemap ) then + vertAleTransportTop = 0.0_RKIND else - call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalVelocityCur, sshCur, dt, & + if (associated(highFreqThicknessNew)) then + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH, layerThicknessCur, & + dt, layerInterfaceVelocity, err, highFreqThicknessNew) + else + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH, layerThicknessCur, & + dt, layerInterfaceVelocity, err) + endif + call ocn_vert_transport_velocity_top(div_hu, layerInterfaceVelocity, & vertAleTransportTop, err) endif @@ -2127,27 +2166,51 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$acc enter data copyin(layerThicknessCur, sshCur) !$acc update device(layerThickEdgeFlux, normalTransportVelocity) #endif - if (associated(highFreqThicknessNew)) then + allocate(div_hu_tr(nVertLevels, nCellsAll), & + projectedSSH_tr(nCellsAll)) #ifdef MPAS_OPENACC - !$acc enter data copyin(highFreqThicknessNew) + !$acc enter data create(projectedSSH_tr, div_hu_tr) #endif - call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalTransportVelocity, sshCur, & - dt, vertAleTransportTop, err, highFreqThicknessNew) + call ocn_diagnostic_solve_divergence_hu(sshCur, layerThickEdgeFlux, & + normalTransportVelocity, dt, div_hu_tr, projectedSSH_tr) + if ( configVertAdvMethod == vertAdvRemap ) then + ! This is the layerInterfaceVelocity implicit in w=0 during Stage 1 + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k) + do iCell = 1, nCellsAll + do k = minLevelCell(iCell), maxLevelCell(iCell) + layerInterfaceVelocity(k,iCell) = div_hu(k,iCell) + enddo + enddo + !$omp end do + !$omp end parallel + else + if (associated(highFreqThicknessNew)) then #ifdef MPAS_OPENACC - !$acc exit data delete(highFreqThicknessNew) + !$acc enter data copyin(highFreqThicknessNew) #endif - else - call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalTransportVelocity, sshCur, & - dt, vertAleTransportTop, err) + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH_tr, layerThicknessCur, & + dt, layerInterfaceVelocity, err, highFreqThicknessNew) +#ifdef MPAS_OPENACC + !$acc exit data delete(highFreqThicknessNew) +#endif + else + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH_tr, layerThicknessCur, & + dt, layerInterfaceVelocity, err) + endif endif + call ocn_vert_transport_velocity_top(div_hu_tr, layerInterfaceVelocity, & + vertAleTransportTop, err) #ifdef MPAS_OPENACC !$acc exit data delete(layerThicknessCur, sshCur) + !$acc exit data delete(projectedSSH_tr, div_hu) !$acc update host(vertAleTransportTop, normalTransportVelocity) #endif + deallocate(layerInterfaceVelocity, div_hu, projectedSSH) + deallocate(div_hu_tr, projectedSSH_tr) call mpas_timer_stop('thick vert trans vel top') call ocn_tend_thick(tendPool, forcingPool) @@ -2420,6 +2483,24 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ !$omp end parallel #endif + if ( associated(layerThicknessLagNew) ) then +#ifdef MPAS_OPENACC + !$acc loop collapse(2) +#else + !$omp do schedule(runtime) private(k) +#endif + do iCell= 1,nCellsAll + do k = 1,nVertLevels + ! Store the current, Lagrangian location for computation of the + ! Eulerian vertical velocity induced by remapping + layerThicknessLagNew(k,iCell) = layerThicknessLagCur(k,iCell) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do +#endif + endif + if (config_compute_active_tracer_budgets) then #ifdef MPAS_OPENACC @@ -2946,23 +3027,6 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ call mpas_timer_stop("si implicit vert mix") - if ( configVertAdvMethod == vertAdvRemap ) then - - call mpas_timer_start("vertical remap") - - call mpas_pool_get_array(statePool, 'layerThicknessLag', layerThicknessLagNew, 2) - - ! Store the current, Lagrangian location for computation of the - ! Eulerian vertical velocity induced by remapping - layerThicknessLagNew = layerThicknessNew - - ! Perform the remapping - call ocn_remap_vert_state(block, err) - - call mpas_timer_stop("vertical remap") - - end if - call mpas_timer_start('si fini') #ifdef MPAS_OPENACC diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index 5c1069f4e3b4..80400ff80f6c 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -194,7 +194,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ sshGrad ! temp for holding ssh gradient real (kind=RKIND), dimension(:,:), allocatable:: & - uTemp + uTemp, layerInterfaceVelocity real (kind=RKIND), dimension(:), allocatable :: & btrvel_temp, & @@ -223,6 +223,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ normalVelocityNew, &! full velocity at new time layerThicknessCur, &! layer thick at current time layerThicknessNew, &! layer thick at new time + layerThicknessLagCur, &! layer thick at current time + layerThicknessLagNew, &! layer thick at new time highFreqThicknessCur, &! high frq thick at current time highFreqThicknessNew, &! high frq thick at new time lowFreqDivergenceCur, &! low frq div at current time @@ -282,10 +284,6 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ real (kind=RKIND), pointer :: forcingTimeIncrement - ! Remap variables - real (kind=RKIND), dimension(:,:), pointer :: & - layerThicknessLagNew - ! These are public variables used from the diagnostics module ! indexSurfaceVelocityZonal, indexSurfaceVelocityMeridional ! indexSSHGradientZonal, indexSSHGradientMeridional @@ -366,6 +364,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call mpas_pool_get_array(statePool, 'layerThickness', & layerThicknessNew, 2) + call mpas_pool_get_array(statePool, 'layerThicknessLag', & + layerThicknessLagCur, 1) + call mpas_pool_get_array(statePool, 'layerThicknessLag', & + layerThicknessLagNew, 2) + call mpas_pool_get_array(statePool, 'highFreqThickness', & highFreqThicknessCur, 1) call mpas_pool_get_array(statePool, 'highFreqThickness', & @@ -488,8 +491,21 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ !$acc end parallel #else !$omp end do - !$omp end parallel #endif + if ( configVertAdvMethod == vertAdvRemap ) then +#ifndef MPAS_OPENACC + !$omp do schedule(runtime) private(k) +#endif + do iCell = 1, nCellsAll + do k = 1,nVertLevels + layerThicknessLagNew(k,iCell) = layerThicknessLagCur(k,iCell) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do +#endif + endif + !$omp end parallel call mpas_pool_begin_iteration(tracersPool) do while ( mpas_pool_get_next_member(tracersPool, groupItr)) @@ -705,16 +721,27 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ ! compute vertAleTransportTop. Use u (rather than & ! normalTransportVelocity) for momentum advection. ! Use the most recent time level available. - - if (associated(highFreqThicknessNew)) then - call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalVelocityCur, sshCur, dt, & - vertAleTransportTop, err, highFreqThicknessNew) + ! The diagnostic variable layerThicknessTarget is updated + nCells = nCellsAll + allocate(layerInterfaceVelocity(nVertLevels, nCells)) + ! This is the Eulerian vertical velocity + ! We update the diagnostic variable vertVelocityTop + call mpas_log_write('ocn_diagnostic_solve_vertVel') + call ocn_diagnostic_solve_vertVel(layerThickEdgeFlux, & + normalVelocityCur, vertVelocityTop) + if ( configVertAdvMethod == vertAdvRemap ) then + vertAleTransportTop(:,:) = 0.0_RKIND else - call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalVelocityCur, sshCur, dt, & + if (associated(highFreqThicknessNew)) then + call ocn_layer_interface_velocity(verticalMeshPool, & + sshCur, layerThicknessCur, vertVelocityTop, & + dt, layerInterfaceVelocity, err, highFreqThicknessNew) + else + call ocn_layer_interface_velocity(verticalMeshPool, & + sshCur, layerThicknessCur, vertVelocityTop, & + dt, layerInterfaceVelocity, err) + endif + call ocn_vert_transport_velocity_top(vertVelocityTop, layerInterfaceVelocity, & vertAleTransportTop, err) endif @@ -1766,27 +1793,46 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ !$acc enter data copyin(layerThicknessCur, sshCur) !$acc update device(layerThickEdgeFlux, normalTransportVelocity) #endif - if (associated(highFreqThicknessNew)) then + call ocn_diagnostic_solve_vertVel(layerThickEdgeFlux, & + normalTransportVelocity, vertTransportVelocityTop) + if ( configVertAdvMethod == vertAdvRemap ) then + call mpas_log_write('set layerInterfaceVelocity = vertVelocityTop') + ! This is the layerInterfaceVelocity implicit in w=0 during Stage 1 + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k) + do iCell = 1, nCellsAll + do k = minLevelCell(iCell), maxLevelCell(iCell) + layerInterfaceVelocity(k,iCell) = vertVelocityTop(k,iCell) + enddo + enddo + !$omp end do + !$omp end parallel + else + if (associated(highFreqThicknessNew)) then #ifdef MPAS_OPENACC - !$acc enter data copyin(highFreqThicknessNew) + !$acc enter data copyin(highFreqThicknessNew) #endif - call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalTransportVelocity, sshCur, & - dt, vertAleTransportTop, err, highFreqThicknessNew) + call ocn_layer_interface_velocity(verticalMeshPool, & + sshCur, layerThicknessCur, vertTransportVelocityTop, & + dt, layerInterfaceVelocity, err, highFreqThicknessNew) #ifdef MPAS_OPENACC - !$acc exit data delete(highFreqThicknessNew) + !$acc exit data delete(highFreqThicknessNew) #endif - else - call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalTransportVelocity, sshCur, & - dt, vertAleTransportTop, err) + else + call ocn_layer_interface_velocity(verticalMeshPool, & + sshCur, layerThicknessCur, vertTransportVelocityTop, & + dt, layerInterfaceVelocity, err) + endif endif + call ocn_vert_transport_velocity_top(vertTransportVelocityTop, layerInterfaceVelocity, & + vertAleTransportTop, err) #ifdef MPAS_OPENACC !$acc exit data delete(layerThicknessCur, sshCur) !$acc update host(vertAleTransportTop, normalTransportVelocity) #endif + deallocate(layerInterfaceVelocity) + call mpas_timer_stop('thick vert trans vel top') call ocn_tend_thick(tendPool, forcingPool) @@ -2577,20 +2623,10 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call mpas_timer_stop("se implicit vert mix") if ( configVertAdvMethod == vertAdvRemap ) then - - call mpas_timer_start("vertical remap") - - call mpas_pool_get_array(statePool, 'layerThicknessLag', layerThicknessLagNew, 2) - - ! Store the current, Lagrangian location for computation of the - ! Eulerian vertical velocity induced by remapping - layerThicknessLagNew = layerThicknessNew - - ! Perform the remapping - call ocn_remap_vert_state(block, err) - - call mpas_timer_stop("vertical remap") - + ! Perform the remapping + ! The time level specified here is the one that will get overwritten during remapping + call ocn_remap_vert_state(block, err, 2) + ! layerThickEdgeFlux will be updated in diagnostic solve end if call mpas_timer_start('se fini') diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F index 0693d00a72fd..5506db7a92cb 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F @@ -203,10 +203,11 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ sshGrad ! temp for holding ssh gradient real (kind=RKIND), dimension(:,:), allocatable:: & - uTemp + div_hu, div_hu_tr, uTemp, layerInterfaceVelocity real (kind=RKIND), dimension(:), allocatable :: & btrvel_temp, & + projectedSSH, projectedSSH_tr, & ! SSH used to determine layerInterfaceVelocity sshTemp, & wctEdge, &! flux water column thickness at edge bottomDepthEdge @@ -286,7 +287,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ ! Remap variables real (kind=RKIND), dimension(:,:), pointer :: & - layerThicknessLagNew + layerThicknessLagCur, layerThicknessLagNew ! These are public variables used from the diagnostics module ! indexSurfaceVelocityZonal, indexSurfaceVelocityMeridional @@ -384,6 +385,9 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ call mpas_pool_get_array(statePool, 'lowFreqDivergence', & lowFreqDivergenceNew, 2) + call mpas_pool_get_array(statePool, 'layerThicknessLag', layerThicknessLagCur, 1) + call mpas_pool_get_array(statePool, 'layerThicknessLag', layerThicknessLagNew, 2) + call mpas_pool_get_dimension_scalar(tracersPool, 'index_salinity', & indexSalinity) @@ -430,6 +434,29 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ endif #endif + if ( configVertAdvMethod == vertAdvRemap ) then + if ( associated(layerThicknessLagCur) ) then +#ifdef MPAS_OPENACC + !$acc loop collapse(2) +#else + !$omp do schedule(runtime) private(k) +#endif + do iCell= 1,nCellsAll + do k = 1,nVertLevels + ! Store the current, Lagrangian location for computation of the + ! Eulerian vertical velocity induced by remapping + layerThicknessLagCur(k,iCell) = layerThicknessCur(k,iCell) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do +#endif + endif + ! Perform the remapping + ! The updated layerThickness goes in layerThicknessCur which is assigned to layerThicknessNew below + call ocn_remap_vert_state(block, err, 1) + end if + ! Initialize * variables that are used to compute baroclinic ! tendencies below. @@ -629,7 +656,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ normalVelocityCur, stage1_tend_time) #ifdef MPAS_OPENACC - !$acc enter data copyin(layerThicknessCur, normalVelocityCur, sshCur) + !$acc enter data copyin(layerThicknessCur, normalVelocityCur, sshCur, layerThicknessTarget) !$acc update device(layerThickEdgeFlux) if (config_use_freq_filtered_thickness .or. associated(highFreqThicknessNew) ) then !$acc enter data create(highFreqThicknessNew) @@ -708,16 +735,27 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ ! compute vertAleTransportTop. Use u (rather than & ! normalTransportVelocity) for momentum advection. ! Use the most recent time level available. - - if (associated(highFreqThicknessNew)) then - call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalVelocityCur, sshCur, dt, & - vertAleTransportTop, err, highFreqThicknessNew) + allocate(layerInterfaceVelocity(nVertLevels, nCellsAll), & + div_hu(nVertLevels, nCellsAll), & + projectedSSH(nCellsAll)) +#ifdef MPAS_OPENACC + !$acc enter data create(projectedSSH, div_hu) +#endif + call ocn_diagnostic_solve_divergence_hu(sshCur, layerThickEdgeFlux, & + normalVelocityCur, dt, div_hu, projectedSSH) + if ( configVertAdvMethod == vertAdvRemap ) then + vertAleTransportTop = 0.0_RKIND else - call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalVelocityCur, sshCur, dt, & + if (associated(highFreqThicknessNew)) then + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH, layerThicknessCur, & + dt, layerInterfaceVelocity, err, highFreqThicknessNew) + else + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH, layerThicknessCur, & + dt, layerInterfaceVelocity, err) + endif + call ocn_vert_transport_velocity_top(div_hu, layerInterfaceVelocity, & vertAleTransportTop, err) endif @@ -1749,27 +1787,51 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ !$acc enter data copyin(layerThicknessCur, sshCur) !$acc update device(layerThickEdgeFlux, normalTransportVelocity) #endif - if (associated(highFreqThicknessNew)) then + allocate(div_hu_tr(nVertLevels, nCellsAll), & + projectedSSH_tr(nCellsAll)) #ifdef MPAS_OPENACC - !$acc enter data copyin(highFreqThicknessNew) + !$acc enter data create(projectedSSH_tr, div_hu_tr) #endif - call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalTransportVelocity, sshCur, & - dt, vertAleTransportTop, err, highFreqThicknessNew) + call ocn_diagnostic_solve_divergence_hu(sshCur, layerThickEdgeFlux, & + normalTransportVelocity, dt, div_hu_tr, projectedSSH_tr) + if ( configVertAdvMethod == vertAdvRemap ) then + ! This is the layerInterfaceVelocity implicit in w=0 during Stage 1 + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k) + do iCell = 1, nCellsAll + do k = minLevelCell(iCell), maxLevelCell(iCell) + layerInterfaceVelocity(k,iCell) = div_hu(k,iCell) + enddo + enddo + !$omp end do + !$omp end parallel + else + if (associated(highFreqThicknessNew)) then #ifdef MPAS_OPENACC - !$acc exit data delete(highFreqThicknessNew) + !$acc enter data copyin(highFreqThicknessNew) #endif - else - call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalTransportVelocity, sshCur, & - dt, vertAleTransportTop, err) + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH_tr, layerThicknessCur, & + dt, layerInterfaceVelocity, err, highFreqThicknessNew) +#ifdef MPAS_OPENACC + !$acc exit data delete(highFreqThicknessNew) +#endif + else + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH_tr, layerThicknessCur, & + dt, layerInterfaceVelocity, err) + endif endif + call ocn_vert_transport_velocity_top(div_hu_tr, layerInterfaceVelocity, & + vertAleTransportTop, err) #ifdef MPAS_OPENACC !$acc exit data delete(layerThicknessCur, sshCur) + !$acc exit data delete(projectedSSH_tr, div_hu) !$acc update host(vertAleTransportTop, normalTransportVelocity) #endif + deallocate(layerInterfaceVelocity, div_hu, projectedSSH) + deallocate(div_hu_tr, projectedSSH_tr) call mpas_timer_stop('thick vert trans vel top') call ocn_tend_thick(tendPool, forcingPool) @@ -2177,10 +2239,48 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ !$acc update device(normalTransportVelocity) #endif +#ifdef MPAS_OPENACC + !$acc enter data copyin(layerThicknessCur, sshCur) + !$acc update device(layerThickEdgeFlux, normalTransportVelocity) +#endif + allocate(div_hu_tr(nVertLevels, nCellsAll), & + projectedSSH_tr(nCellsAll)) +#ifdef MPAS_OPENACC + !$acc enter data create(projectedSSH_tr, div_hu_tr) +#endif + call ocn_diagnostic_solve_divergence_hu(sshCur, layerThickEdgeFlux, & + normalTransportVelocity, dt, div_hu_tr, projectedSSH_tr) + if ( configVertAdvMethod == vertAdvRemap ) then + ! This is the layerInterfaceVelocity implicit in w=0 during Stage 1 + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k) + do iCell = 1, nCellsAll + do k = minLevelCell(iCell), maxLevelCell(iCell) + layerInterfaceVelocity(k,iCell) = div_hu(k,iCell) + enddo + enddo + !$omp end do + !$omp end parallel + else + if (associated(highFreqThicknessNew)) then +#ifdef MPAS_OPENACC + !$acc enter data copyin(highFreqThicknessNew) +#endif + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH_tr, layerThicknessCur, & + dt, layerInterfaceVelocity, err, highFreqThicknessNew) +#ifdef MPAS_OPENACC + !$acc exit data delete(highFreqThicknessNew) +#endif + else + call ocn_layer_interface_velocity(verticalMeshPool, & + projectedSSH_tr, layerThicknessCur, & + dt, layerInterfaceVelocity, err) + endif + endif call ocn_vert_transport_velocity_top( & - verticalMeshPool, layerThicknessCur, & - layerThickEdgeFlux, normalTransportVelocity,sshCur, dt, & - vertAleTransportTop, err) + div_hu_tr, layerInterfaceVelocity, vertAleTransportTop, err) call ocn_tend_thick(tendPool, forcingPool) @@ -2211,6 +2311,24 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ !$omp end do !$omp end parallel #endif + if ( associated(layerThicknessLagNew) ) then +#ifdef MPAS_OPENACC + !$acc loop collapse(2) +#else + !$omp do schedule(runtime) private(k) +#endif + do iCell= 1,nCellsAll + do k = 1,nVertLevels + ! Store the current, Lagrangian location for computation of the + ! Eulerian vertical velocity induced by remapping + layerThicknessLagNew(k,iCell) = layerThicknessLagCur(k,iCell) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do +#endif + endif + endif ! nowTimeStamp @@ -2749,23 +2867,6 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ call mpas_timer_stop("se implicit vert mix") - if ( configVertAdvMethod == vertAdvRemap ) then - - call mpas_timer_start("vertical remap") - - call mpas_pool_get_array(statePool, 'layerThicknessLag', layerThicknessLagNew, 2) - - ! Store the current, Lagrangian location for computation of the - ! Eulerian vertical velocity induced by remapping - layerThicknessLagNew = layerThicknessNew - - ! Perform the remapping - call ocn_remap_vert_state(block, err) - - call mpas_timer_stop("vertical remap") - - end if - call mpas_timer_start('se fini') #ifdef MPAS_OPENACC diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile index 92019f57efe8..a0ae75d06bc3 100644 --- a/components/mpas-ocean/src/shared/Makefile +++ b/components/mpas-ocean/src/shared/Makefile @@ -231,7 +231,7 @@ mpas_ocn_vel_tidal_potential.o: mpas_ocn_diagnostics_variables.o mpas_ocn_vertical_advection.o: mpas_ocn_config.o -mpas_ocn_vertical_regrid.o: mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_mesh.o +mpas_ocn_vertical_regrid.o: mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_mesh.o mpas_ocn_thick_ale.o mpas_ocn_vertical_remap.o: mpas_ocn_config.o mpas_ocn_constants.o mpas_ocn_diagnostics_variables.o mpas_ocn_mesh.o mpas_ocn_vertical_regrid.o diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index a86c24386c8f..22adafa5179c 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -59,8 +59,11 @@ module ocn_diagnostics public :: ocn_diagnostic_solve, & ocn_diagnostic_solve_layerThicknessEdge, & + ocn_diagnostic_solve_divergence_hu, & ocn_diagnostic_solve_wctEdge, & + ocn_diagnostic_solve_vertVel, & ocn_relativeVorticity_circulation, & + ocn_layer_interface_velocity, & ocn_vert_transport_velocity_top, & ocn_fuperp, & ocn_filter_btr_mode_tend_vel, & @@ -281,49 +284,6 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMe call ocn_diagnostic_solve_MLEvel(layerThickEdgeFlux, normalMLEVelocity, vertMLEBolusVelocityTop) end if - if ( configVertAdvMethod == vertAdvRemap ) then - - ! Overwrite vertVelocityTop with exact Eulerian solution from remapping - call mpas_pool_get_array(statePool, 'layerThicknessLag', & - layerThicknessLag, timeLevel) - call ocn_diagnostic_solve_vertVel_remap(layerThickness, & - layerThicknessLag, dt, vertVelocityTop) - - ! Overwrite vertical tracer budget terms with new vertVelocityTop - if ( config_compute_active_tracer_budgets ) then - numTracers = size(activeTracers, dim=1) - do iTracer = 1, numTracers - !$omp parallel - !$omp do schedule(runtime) private(k,kmin,kmax) - do iCell = 1, nCells - kmin = minLevelCell(iCell) - kmax = maxLevelCell(iCell) - - activeTracerVerticalAdvectionTopFlux(iTracer, :, iCell) = 0.0_RKIND - do k = kmin+1, kmax - activeTracerVerticalAdvectionTopFlux(iTracer, k, iCell) = & - min(0.0_RKIND,vertVelocityTop(k, iCell))* & - activeTracers(iTracer, k-1, iCell) & - + max(0.0_RKIND,vertVelocityTop(k,iCell))* & - activeTracers(iTracer, k, iCell) - enddo - - do k = kmin, kmax-1 - activeTracerVerticalAdvectionTendency(iTracer, k, iCell) = ( & - activeTracerVerticalAdvectionTopFlux(iTracer, k+1, iCell) & - - activeTracerVerticalAdvectionTopFlux(iTracer, k, iCell) & - ) / layerThickness(k, iCell) - enddo - activeTracerVerticalAdvectionTendency(iTracer, kmax, iCell) = & - - activeTracerVerticalAdvectionTopFlux(iTracer, kmax, iCell) & - / layerThickness(kmax, iCell) - enddo - !$omp end do - !$omp end parallel - enddo - endif - endif - ! ! Compute kinetic energy ! @@ -2265,6 +2225,317 @@ end subroutine ocn_diagnostic_solve_MLEvel!}}} ! !----------------------------------------------------------------------- + subroutine ocn_diagnostic_solve_vorticityCell( & + relativeVorticity, & + relativeVorticityCell)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + relativeVorticity !< [in] relative vorticity in cells + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + relativeVorticityCell !< [out] relative vorticity in cell + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: & + i, j, k, &! loop indices + iCell, iVertex, &! indices for cells, vertices + nCells, &! number of cells + kmin, kmax ! min, max active vert level index + + real(kind=RKIND) :: & + invAreaCell1, &! 1/cell area + r_tmp ! tmp for common factors + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + ! Compute relative vorticity in cells + +#ifdef MPAS_OPENACC + !$acc parallel loop gang vector & + !$acc present(relativeVorticityCell, relativeVorticity, & + !$acc invAreaCell, kiteAreasOnVertex, & + !$acc nEdgesOnCell, kiteIndexOnCell, verticesOnCell, & + !$acc minLevelCell, maxLevelCell) & + !$acc private(i, iVertex, j, k, kmin, kmax, invAreaCell1) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(i, iVertex, j, k, kmin, kmax, invAreaCell1) +#endif + do iCell = 1, nCellsOwned + relativeVorticityCell(:,iCell) = 0.0_RKIND + invAreaCell1 = invAreaCell(iCell) + kmin = minLevelCell(iCell) + kmax = maxLevelCell(iCell) + + do i = 1, nEdgesOnCell(iCell) + j = kiteIndexOnCell(i, iCell) + iVertex = verticesOnCell(i, iCell) + do k = kmin,kmax + relativeVorticityCell(k,iCell) = & + relativeVorticityCell(k,iCell) + & + kiteAreasOnVertex(j,iVertex)* & + relativeVorticity(k,iVertex)*invAreaCell1 + end do + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + end subroutine ocn_diagnostic_solve_vorticityCell!}}} + + subroutine ocn_diagnostic_solve_KE( & + normalVelocity, kineticEnergyCell)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + normalVelocity !< [in] velocity normal to edge + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + kineticEnergyCell !< [out] kinetic energy in cell + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: & + i, j, k, &! loop indices + iCell, iVertex, iEdge, &! indices for cells, vertices, edges + nCells, nEdges, &! number of cells, edges + kmin, kmax, &! min, max active vert level index + eoe, &! index for edge on edge + edgeSignOnCell_temp ! local tmp for sign across edge + + real(kind=RKIND) :: & + invAreaCell1, &! 1/cell area + dcEdge_temp, &! local val of dcEdge + dvEdge_temp, &! local val of dvEdge + r_tmp, &! tmp for common factors + weightsOnEdge_temp ! local tmp for weights on edge + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + ! + ! Compute kinetic energy + ! +#ifdef MPAS_OPENACC + !$acc parallel loop gang vector & + !$acc present(kineticEnergyCell, & + !$acc normalVelocity, & + !$acc dcEdge, dvEdge, & + !$acc invAreaCell, nEdgesOnCell, edgesOnCell, & + !$acc edgeSignOnCell, minLevelCell, maxLevelCell) & + !$acc private(i, k, kmin, kmax, & + !$acc iEdge, invAreaCell1, r_tmp, dcEdge_temp, & + !$acc dvEdge_temp, edgeSignOnCell_temp) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(i, k, kmin, kmax, & + !$omp iEdge, invAreaCell1, r_tmp, dcEdge_temp, & + !$omp dvEdge_temp, edgeSignOnCell_temp) +#endif + do iCell = 1, nCells + kineticEnergyCell(:,iCell) = 0.0_RKIND + invAreaCell1 = invAreaCell(iCell) + kmin = minLevelCell(iCell) + kmax = maxLevelCell(iCell) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + edgeSignOnCell_temp = edgeSignOnCell(i, iCell) + dcEdge_temp = dcEdge(iEdge) + dvEdge_temp = dvEdge(iEdge) + do k = kmin,kmax + r_tmp = dvEdge_temp*normalVelocity(k,iEdge)*invAreaCell1 + kineticEnergyCell(k,iCell) = kineticEnergyCell(k,iCell) & + + 0.25*r_tmp*dcEdge_temp* & + normalVelocity(k,iEdge) + end do + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + end subroutine ocn_diagnostic_solve_KE!}}} + + subroutine ocn_diagnostic_solve_tangentialVel( & + normalVelocity, tangentialVelocity)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + normalVelocity !< [in] velocity normal to edge + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tangentialVelocity !< [out] tangential velocity at edge + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: & + i, k, &! loop indices + iEdge, &! indices for edges + nEdges, &! number of edges + kmin, kmax, &! min, max active vert level index + eoe ! index for edge on edge + + real(kind=RKIND) :: & + weightsOnEdge_temp ! local tmp for weights on edge + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + ! Compute tangential velocity + + nEdges = nEdgesHalo(2) + +#ifdef MPAS_OPENACC + !$acc parallel loop gang vector & + !$acc present(tangentialVelocity, normalVelocity, & + !$acc weightsOnEdge, nEdgesOnEdge, edgesOnEdge, & + !$acc minLevelEdgeBot, maxLevelEdgeTop) & + !$acc private(i, k, kmin, kmax, eoe, weightsOnEdge_temp) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(i, k, kmin, kmax, eoe, weightsOnEdge_temp) +#endif + do iEdge = 1, nEdges + tangentialVelocity(:,iEdge) = 0.0_RKIND + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + do i = 1, nEdgesOnEdge(iEdge) + eoe = edgesOnEdge(i,iEdge) + weightsOnEdge_temp = weightsOnEdge(i, iEdge) + do k = kmin,kmax + tangentialVelocity(k,iEdge) = & + tangentialVelocity(k,iEdge) + weightsOnEdge_temp* & + normalVelocity(k, eoe) + end do + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + end subroutine ocn_diagnostic_solve_tangentialVel!}}} + + subroutine ocn_diagnostic_solve_divergence_hu(layerThicknessEdgeFlux, & + normalVelocity, div_hu)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThicknessEdgeFlux, &!< [in] layerThickness at edge + normalVelocity !< [in] transport + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + real (kind=RKIND), dimension(:,:), intent(out) :: & + div_hu ! divergence of (thickness*velocity) + + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + + real (kind=RKIND) :: & + r_tmp, &! tmp for common factors + flux, &! + invAreaCell1, &! 1/cell area + div_hu_btr ! divergence of h*u for barotropic u + + integer :: & + iEdge, iCell, k, i, &! edge, cell, vert, nbr indices + nCells, &! number of cells + kmin, kmax ! min/max layer indices for active layers + + nCells = nCellsHalo(2) + + ! + ! thickness-weighted divergence and barotropic divergence + ! + ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3. +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(oldSSH, normalVelocity, div_hu, & + !$acc layerThicknessEdgeFlux, dvEdge, invAreaCell, & + !$acc nEdgesOnCell, edgesOnCell, edgeSignOnCell, & + !$acc minLevelEdgeBot, maxLevelEdgeTop) & + !$acc private(i, iEdge, k, kmin, kmax, div_hu_btr, & + !$acc invAreaCell1, flux) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(i, iEdge, k, kmin, kmax, div_hu_btr, & + !$omp invAreaCell1, flux) +#endif + do iCell = 1, nCells + div_hu(:,iCell) = 0.0_RKIND + div_hu_btr = 0.0_RKIND + invAreaCell1 = invAreaCell(iCell) + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + + do k = kmin, kmax + r_tmp = dvEdge(iEdge)*normalVelocity(k,iEdge)*invAreaCell1 + + divergence(k,iCell) = divergence(k,iCell) & + - edgeSignOnCell(i,iCell) * r_tmp + + flux = layerThicknessEdgeFlux(k,iEdge) * & + r_tmp * edgeSignOnCell(i,iCell) + + div_hu(k,iCell) = div_hu(k,iCell) - flux + end do + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + end subroutine ocn_diagnostic_solve_divergence_hu!}}} + subroutine ocn_diagnostic_solve_vortVel( & relativeVorticity, layerThicknessEdgeFlux, & normalVelocity, normalTransportVelocity, & @@ -2473,6 +2744,97 @@ subroutine ocn_diagnostic_solve_vortVel( & end subroutine ocn_diagnostic_solve_vortVel!}}} + subroutine ocn_diagnostic_solve_vertVel( & + layerThicknessEdgeFlux, & + normalVelocity, & + vertVelTop)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThicknessEdgeFlux, &!< [in] thickness flux at edge + normalVelocity !< [in] velocity normal to edge + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + vertVelTop !< [out] vertical velocity top of cell + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: & + i, j, k, &! loop indices + iCell, iVertex, iEdge, &! indices for cells, vertices, edges + nCells, nEdges, &! number of cells, edges + kmin, kmax, &! min, max active vert level index + eoe, &! index for edge on edge + edgeSignOnCell_temp ! local tmp for sign across edge + + real(kind=RKIND) :: & + invAreaCell1, &! 1/cell area + dcEdge_temp, &! local val of dcEdge + dvEdge_temp, &! local val of dvEdge + r_tmp, &! tmp for common factors + weightsOnEdge_temp ! local tmp for weights on edge + + real (kind=RKIND), dimension(:,:), allocatable :: & + div_hu ! local divergence of thick-weighted velocity + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + + ! Need 0,1,2 halo cells + nCells = nCellsHalo(2) + + ! + ! Compute divergence, kinetic energy, and vertical velocity + ! + allocate(div_hu(nVertLevels, nCells)) + !$acc enter data create(div_hu) + + call ocn_diagnostic_solve_divergence_hu(layerThicknessEdgeFlux, & + normalVelocity, div_hu)!{{{ + +#ifdef MPAS_OPENACC + !$acc parallel loop gang vector & + !$acc present(vertVelTop, & + !$acc normalVelocity, & + !$acc minLevelCell, maxLevelCell) & + !$acc private(k, kmin, kmax) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, kmin, kmax) +#endif + do iCell = 1, nCells + kmin = minLevelCell(iCell) + kmax = maxLevelCell(iCell) + ! Vertical velocity at bottom is zero, initialized above. + vertVelTop(1:kmin-1,iCell) = 0.0_RKIND + vertVelTop(kmax+1 ,iCell) = 0.0_RKIND + do k = kmax, 1, -1 + vertVelTop(k,iCell) = & + vertVelTop(k+1,iCell) - div_hu(k, iCell) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + !$acc exit data delete(div_hu) + deallocate(div_hu) + + end subroutine ocn_diagnostic_solve_vertVel!}}} + !*********************************************************************** ! ! routine ocn_diagnostic_solve_vertVel_remap @@ -3049,32 +3411,104 @@ end subroutine ocn_diagnostic_solve_ssh!}}} ! !----------------------------------------------------------------------- - subroutine ocn_vert_transport_velocity_top(verticalMeshPool, & - oldLayerThickness, layerThicknessEdgeFlux, & - normalVelocity, oldSSH, dt, vertAleTransportTop, & - err, newHighFreqThickness)!{{{ - - !----------------------------------------------------------------- - ! input variables - !----------------------------------------------------------------- + subroutine ocn_layer_interface_velocity(verticalMeshPool, & + oldSSH, oldLayerThickness, vertVelTop, dt, & + layerInterfaceVelocity, err, newHighFreqThickness)!{{{ type (mpas_pool_type), intent(in) :: & - verticalMeshPool !< [in] vertical mesh information + verticalMeshPool !< [in] vertical mesh information real (kind=RKIND), dimension(:), intent(in) :: & - oldSSH !< [in] sea surface height at old time + oldSSH !< [in] SSH at old time real (kind=RKIND), dimension(:,:), intent(in) :: & - oldLayerThickness, &!< [in] layer thickness at old time - layerThicknessEdgeFlux, &!< [in] layerThickness at edge - normalVelocity !< [in] transport + oldLayerThickness, &!< [in] layer thickness at old time + vertVelTop !< [in] Eulerian vertical velocity + + real (kind=RKIND), intent(in) :: & + dt !< Input: time step ! alters ALE thickness when freq-filtering real (kind=RKIND), dimension(:,:), intent(in), optional :: & newHighFreqThickness !< [in] high frequency thickness - real (kind=RKIND), intent(in) :: & - dt !< Input: time step + real (kind=RKIND), dimension(:,:), intent(out) :: & + layerInterfaceVelocity !< [out] layer interface motion over the time step + + integer, intent(out) :: err !< [out] error flag + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: & + nCells, &! number of cells + iCell, k ! cell, vert indices + + real (kind=RKIND), dimension(:), allocatable :: projectedSSH + !real (kind=RKIND), dimension(:,:) :: & + ! layerThicknessTarget ! This is already a diagnostic variable + + err = 0 + + ! Only need to compute over 0 and 1 halos + nCells = nCellsAll + allocate(projectedSSH(nCells)) + ! + ! Compute desired thickness at new time + ! + if (configVertAdvMethod == vertAdvRemap) then + ! The layer does not move here + layerInterfaceVelocity = 0.0_RKIND + else +#ifndef MPAS_OPENACC + !$omp parallel + !$omp do schedule(runtime) private(k) +#endif + do iCell = 1,nCells + projectedSSH(iCell) = oldSSH(iCell) + & + vertVelTop(minLevelCell(iCell), iCell) * dt + enddo +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + if (present(newHighFreqThickness)) then + call ocn_ALE_thickness(verticalMeshPool, projectedSSH, & + layerThicknessTarget, err, newHighFreqThickness) + else + call ocn_ALE_thickness(verticalMeshPool, projectedSSH, & + layerThicknessTarget, err) + endif +#ifndef MPAS_OPENACC + !$omp parallel + !$omp do schedule(runtime) private(k) +#endif + do iCell = 1,nCells + do k = minLevelCell(iCell), maxLevelCell(iCell) + layerInterfaceVelocity(k,iCell) = (layerThicknessTarget(k,iCell) - & + oldLayerThickness(k,iCell))/dt + enddo + enddo +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + endif + + end subroutine ocn_layer_interface_velocity!}}} + + subroutine ocn_vert_transport_velocity_top(& + vertEulerVelocityTop, layerInterfaceVelocity, & + vertAleTransportTop, err)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertEulerVelocityTop, & + layerInterfaceVelocity !< [in] layer interface motion over the time step !----------------------------------------------------------------- ! output variables @@ -3090,21 +3524,8 @@ subroutine ocn_vert_transport_velocity_top(verticalMeshPool, & !----------------------------------------------------------------- integer :: & - iEdge, iCell, k, i, &! edge, cell, vert, nbr indices - nCells, &! number of cells - kmin, kmax ! min/max layer indices for active layers - - real (kind=RKIND) :: & - flux, &! - invAreaCell1, &! 1/cell area - div_hu_btr ! divergence of h*u for barotropic u - - real (kind=RKIND), dimension(:), allocatable :: & - projectedSSH ! projectedSSH: projected SSH at a new time - - real (kind=RKIND), dimension(:,:), allocatable :: & - div_hu, &! divergence of (thickness*velocity) - ALE_Thickness ! ALE thickness at new time + nCells, &! number of cells + iCell, k ! cell, vert indices ! End preamble !----------------------------------------------------------------- @@ -3112,73 +3533,15 @@ subroutine ocn_vert_transport_velocity_top(verticalMeshPool, & err = 0 - if (config_vert_coord_movement == 'impermeable_interfaces' .or. & - configVertAdvMethod == vertAdvRemap) then + if (config_vert_coord_movement == 'impermeable_interfaces') then vertAleTransportTop = 0.0_RKIND !$acc update device(vertAleTransportTop) return end if - allocate(div_hu(nVertLevels, nCellsAll), & - projectedSSH(nCellsAll), & - ALE_Thickness(nVertLevels, nCellsAll)) - !$acc enter data create(projectedSSH, ALE_Thickness, div_hu) ! Only need to compute over 0 and 1 halos - nCells = nCellsHalo( 1 ) - - ! - ! thickness-weighted divergence and barotropic divergence - ! - ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3. -#ifdef MPAS_OPENACC - !$acc parallel loop & - !$acc present(projectedSSH, oldSSH, normalVelocity, div_hu, & - !$acc layerThicknessEdgeFlux, dvEdge, invAreaCell, & - !$acc nEdgesOnCell, edgesOnCell, edgeSignOnCell, & - !$acc minLevelEdgeBot, maxLevelEdgeTop) & - !$acc private(i, iEdge, k, kmin, kmax, div_hu_btr, & - !$acc invAreaCell1, flux) -#else - !$omp parallel - !$omp do schedule(runtime) & - !$omp private(i, iEdge, k, kmin, kmax, div_hu_btr, & - !$omp invAreaCell1, flux) -#endif - do iCell = 1, nCells - div_hu(:,iCell) = 0.0_RKIND - div_hu_btr = 0.0_RKIND - invAreaCell1 = invAreaCell(iCell) - do i = 1, nEdgesOnCell(iCell) - iEdge = edgesOnCell(i, iCell) - kmin = minLevelEdgeBot(iEdge) - kmax = maxLevelEdgeTop(iEdge) - - do k = kmin, kmax - flux = layerThicknessEdgeFlux(k,iEdge)* & - normalVelocity(k,iEdge)*dvEdge(iEdge)* & - edgeSignOnCell(i,iCell) * invAreaCell1 - div_hu(k,iCell) = div_hu(k,iCell) - flux - div_hu_btr = div_hu_btr - flux - end do - end do - projectedSSH(iCell) = oldSSH(iCell) - dt*div_hu_btr - end do -#ifndef MPAS_OPENACC - !$omp end do - !$omp end parallel -#endif - - ! - ! Compute desired thickness at new time - ! - if (present(newHighFreqThickness)) then - call ocn_ALE_thickness(verticalMeshPool, projectedSSH, & - ALE_thickness, err, newHighFreqThickness) - else - call ocn_ALE_thickness(verticalMeshPool, projectedSSH, & - ALE_thickness, err) - endif + nCells = nCellsAll ! ! Vertical transport through layer interfaces @@ -3186,12 +3549,12 @@ subroutine ocn_vert_transport_velocity_top(verticalMeshPool, & ! Vertical transport through layer interface at top and ! bottom is zero. Here we are using solving the continuity ! equation for vertAleTransportTop ($w^t$), and using - ! ALE_Thickness for thickness at the new time. + ! layerThicknessTarget for thickness at the new time. #ifdef MPAS_OPENACC !$acc parallel loop gang vector & - !$acc present(vertAleTransportTop, ALE_Thickness, & - !$acc oldLayerThickness, div_hu, & + !$acc present(vertAleTransportTop, & + !$acc oldLayerThickness, vertEulerVelocityTop, & !$acc minLevelCell, maxLevelCell) #else !$omp parallel @@ -3202,9 +3565,8 @@ subroutine ocn_vert_transport_velocity_top(verticalMeshPool, & vertAleTransportTop(maxLevelCell(iCell)+1,iCell) = 0.0_RKIND do k = maxLevelCell(iCell), minLevelCell(iCell)+1, -1 vertAleTransportTop(k,iCell) = & - vertAleTransportTop(k+1,iCell) - div_hu(k,iCell) & - - (ALE_Thickness(k,iCell) - & - oldLayerThickness(k,iCell))/dt + vertAleTransportTop(k+1,iCell) + vertEulerVelocityTop(k,iCell) & + - layerInterfaceVelocity(k,iCell) end do end do #ifndef MPAS_OPENACC @@ -3212,11 +3574,6 @@ subroutine ocn_vert_transport_velocity_top(verticalMeshPool, & !$omp end parallel #endif - !$acc exit data delete(projectedSSH, ALE_Thickness, div_hu) - deallocate(div_hu, & - projectedSSH, & - ALE_Thickness) - !-------------------------------------------------------------------- end subroutine ocn_vert_transport_velocity_top!}}} diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F index abdc7c1ff9b4..9004d1299c67 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F @@ -56,6 +56,7 @@ module ocn_diagnostics_variables real (kind=RKIND), dimension(:,:), pointer :: salineContractCoeff real (kind=RKIND), dimension(:,:), pointer :: BruntVaisalaFreqTop real (kind=RKIND), dimension(:,:), pointer :: tangentialVelocity + real (kind=RKIND), dimension(:,:), pointer :: layerThicknessTarget real (kind=RKIND), dimension(:,:), pointer :: layerThickEdgeDrag real (kind=RKIND), dimension(:,:), pointer :: layerThickEdgeFlux real (kind=RKIND), dimension(:,:), pointer :: layerThickEdgeMean @@ -291,6 +292,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e layerThicknessCellWetDry) call mpas_pool_get_array(diagnosticsPool, 'sshCellWetDry', & sshCellWetDry) + call mpas_pool_get_array(diagnosticsPool, 'layerThicknessTarget', & + layerThicknessTarget) end if if (config_use_spatially_variable_upwind) then call mpas_pool_get_array(diagnosticsPool, 'upwindFactor', & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_thick_ale.F b/components/mpas-ocean/src/shared/mpas_ocn_thick_ale.F index a4e861f6ed67..d12d71ae4d72 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_thick_ale.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_thick_ale.F @@ -157,7 +157,7 @@ subroutine ocn_ALE_thickness(verticalMeshPool, SSH, ALE_thickness, & restingThickness) !xacc enter data copyin(restingThickness) - nCells = nCellsHalo(1) + nCells = nCellsAll ! ! ALE thickness alteration due to SSH (z-star) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_regrid.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_regrid.F index feeb585ab6f8..df4f4b90750b 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vertical_regrid.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_regrid.F @@ -33,6 +33,7 @@ module ocn_vertical_regrid use ocn_constants use ocn_config use ocn_mesh + use ocn_thick_ale implicit none private @@ -68,40 +69,54 @@ module ocn_vertical_regrid ! !----------------------------------------------------------------------- - subroutine ocn_vert_regrid(restingThickness, layerThicknessLag, & - layerThicknessTarget, err) + subroutine ocn_vert_regrid(verticalMeshPool, layerThicknessLag, & + layerThicknessTarget, err, newHighFreqThickness) + + type (mpas_pool_type), intent(in) :: & + verticalMeshPool ! vertical mesh information real (kind=RKIND), dimension(:, :), intent(in) :: & - layerThicknessLag, & ! layerThickness after the lagrangian step - restingThickness - real (kind=RKIND), dimension(:, :), intent(out) :: & + layerThicknessLag ! layerThickness after the lagrangian step + + ! alters ALE thickness when freq-filtering + real (kind=RKIND), dimension(:,:), intent(in), optional :: & + newHighFreqThickness !< [in] high frequency thickness + + real (kind=RKIND), dimension(:,:), intent(inout) :: & layerThicknessTarget ! adjusted target layerThickness for remapping + ! generally set equal to layerThicknessTarget diagnostics variable + integer, intent(out) :: err !< Output: Error flag integer :: iCell, k, kmin, kmax - real (kind=RKIND) :: totalThickness + real (kind=RKIND), dimension(:), allocatable :: projectedSSH err = 0 - ! Calculate z-star layer locations - layerThicknessTarget = 0.0_RKIND + allocate(projectedSSH(nCellsAll)) !$omp parallel - !$omp do schedule(runtime) private(k,totalThickness) + !$omp do schedule(runtime) private(k) do iCell = 1, nCellsAll - totalThickness = 0.0_RKIND - do k = maxLevelCell(iCell), minLevelCell(iCell), -1 - totalThickness = totalThickness + layerThicknessLag(k,iCell) - end do + projectedSSH(iCell) = -bottomDepth(iCell) do k = maxLevelCell(iCell), minLevelCell(iCell), -1 - layerThicknessTarget(k, iCell) = restingThickness(k, iCell) * & - totalThickness / bottomDepth(iCell) + projectedSSH(iCell) = projectedSSH(iCell) + layerThicknessLag(k,iCell) end do end do !$omp end do !$omp end parallel + ! + ! Compute desired thickness at new time + ! + if (present(newHighFreqThickness)) then + call ocn_ALE_thickness(verticalMeshPool, projectedSSH, & + layerThicknessTarget, err, newHighFreqThickness) + else + call ocn_ALE_thickness(verticalMeshPool, projectedSSH, & + layerThicknessTarget, err) + endif end subroutine ocn_vert_regrid diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F index eaf40063e860..51f5c73fedfc 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F @@ -75,11 +75,13 @@ module ocn_vertical_remap ! !----------------------------------------------------------------------- - subroutine ocn_remap_vert_state(block, err) + subroutine ocn_remap_vert_state(block, err, timeLevel) type (block_type), intent(in) :: block + integer, intent(in) :: timeLevel integer, intent(out) :: err + real (kind=RKIND) :: predSSH integer :: nCells, nEdges integer :: nLayers, nLevels, nVars, nDoFs, nTracers integer :: cell1, cell2, iCell, jCell, iEdge, iTrac @@ -95,12 +97,11 @@ subroutine ocn_remap_vert_state(block, err) real (kind=RKIND), dimension(:, :, :), allocatable :: dataNow, dataNew real (kind=RKIND), dimension(:, :), allocatable :: & - layerThicknessNew, & ! layerThickness (new target grid) layerThickEdgeNew, & ! layerThickness at edges (new target grid) heightCellNow, heightEdgeNow, & ! depth at layer interfaces (lagrangian grid) heightCellNew, heightEdgeNew ! depth at layer interfaces (new target grid) - real (kind=RKIND), dimension(:, :), pointer :: layerThickness, restingThickness + real (kind=RKIND), dimension(:, :), pointer :: layerThicknessLag, restingThickness, layerThicknessNew real (kind=RKIND), dimension(:, :), pointer :: normalVelocity real (kind=RKIND), dimension(:, :), pointer :: highFreqThickness,lowFreqDivergence real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroup @@ -111,6 +112,7 @@ subroutine ocn_remap_vert_state(block, err) if (itimestepLastRemap < config_vert_remap_interval) return itimestepLastRemap = 0 + call mpas_timer_start("vertical remap") ! Remapping currently only supports one block for the whole domain rather ! than sub-blocks call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool) @@ -120,11 +122,12 @@ subroutine ocn_remap_vert_state(block, err) call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) - call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 2) - call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 2) + call mpas_pool_get_array(statePool, 'layerThicknessLag', layerThicknessLag, timeLevel) + call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessNew, timeLevel) + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, timeLevel) if (config_use_freq_filtered_thickness) then - call mpas_pool_get_array(statePool, 'highFreqThickness', highFreqThickness, 2) - call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergence, 2) + call mpas_pool_get_array(statePool, 'highFreqThickness', highFreqThickness, timeLevel) + call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergence, timeLevel) end if nCells = nCellsAll @@ -136,49 +139,79 @@ subroutine ocn_remap_vert_state(block, err) allocate(heightCellNew(nVertLevels + 1, nCells)) allocate(heightEdgeNow(nVertLevels + 1, nEdges)) allocate(heightEdgeNew(nVertLevels + 1, nEdges)) - allocate(layerThicknessNew(nVertLevels, nCells + 1)) allocate(layerThickEdgeNew(nVertLevels, nEdges)) + if ( associated(layerThicknessLag) ) then +#ifdef MPAS_OPENACC + !$acc loop collapse(2) +#else + !$omp do schedule(runtime) private(k) +#endif + do iCell= 1,nCells + do k = 1,nVertLevels + ! Store the current, Lagrangian location for computation of the + ! Eulerian vertical velocity induced by remapping + layerThicknessLag(k,iCell) = layerThicknessNew(k,iCell) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do +#endif + endif ! Compute new layer thicknesses based on vertical coordinate choice ! Currently, zstar coordinate is hard-coded - call ocn_vert_regrid(restingThickness, layerThickness, & - layerThicknessNew, err) - - ! Compute the layer interface locations between bottomDepth and -ssh + ! The diagnostic variable layerThicknessTarget is updated + if (associated(highFreqThickness)) then + call ocn_vert_regrid(verticalMeshPool, layerThicknessLag, & + layerThicknessTarget, err, highFreqThickness) + else + call ocn_vert_regrid(verticalMeshPool, layerThicknessLag, & + layerThicknessTarget, err) + endif + ! Compute the depth of the top of each layer between bottomDepth and -ssh ! (rather than between -bottomDepth and ssh due to PPR limitations) !$omp parallel - !$omp do schedule(runtime) private(k) + !$omp do schedule(runtime) private(k, predSSH) do iCell = 1, nCells - heightCellNow(:, iCell) = bottomDepth(iCell) - heightCellNew(:, iCell) = bottomDepth(iCell) + if (maxLevelCell(iCell) < 1) cycle + heightCellNow(maxLevelCell(iCell) + 1:, iCell) = bottomDepth(iCell) + heightCellNew(maxLevelCell(iCell) + 1:, iCell) = bottomDepth(iCell) ! reconstruct "now" heights do k = maxLevelCell(iCell), minLevelCell(iCell), -1 heightCellNow(k, iCell) = & - heightCellNow(k + 1, iCell) - layerThickness(k, iCell) + heightCellNow(k + 1, iCell) - layerThicknessLag(k, iCell) end do ! reconstruct "new" heights - heightCellNew(minLevelCell(iCell), iCell) = heightCellNow(minLevelCell(iCell), iCell) - do k = maxLevelCell(iCell), minLevelCell(iCell)+1, -1 + do k = maxLevelCell(iCell), minLevelCell(iCell) + 1, -1 heightCellNew(k, iCell) = & - heightCellNew(k + 1, iCell) - layerThicknessNew(k, iCell) + heightCellNew(k + 1, iCell) - layerThicknessTarget(k, iCell) end do + + ! We force the reconstruction to agree on the total baroclinic depth so we don't have propagagtion of floating point errors + heightCellNew(minLevelCell(iCell), iCell) = heightCellNow(minLevelCell(iCell), iCell) + + predSSH = heightCellNew(minLevelCell(iCell)+1, iCell) - layerThicknessTarget(minLevelCell(iCell), iCell) end do !$omp end do !$omp end parallel + ! Assign new layerThicknesses + + layerThicknessNew = layerThicknessTarget + #ifdef MPAS_OPENACC !$acc parallel loop & - !$acc present(layerThickness, layerThickEdgeMean, & + !$acc present(layerThicknessNew, layerThickEdgeMean, & !$acc minLevelEdgeBot, maxLevelEdgeTop, cellsOnEdge) & !$acc private(k, kmin, kmax, cell1, cell2) #else !$omp parallel !$omp do schedule(runtime) private(k, kmin, kmax, cell1, cell2) #endif - do iEdge = 1, nEdgesAll + do iEdge = 1, nEdges kmin = minLevelEdgeBot(iEdge) kmax = maxLevelEdgeTop(iEdge) cell1 = cellsOnEdge(1,iEdge) @@ -212,10 +245,10 @@ subroutine ocn_remap_vert_state(block, err) ! We only need to define the kBot+1 depth for one active cell ! as the absolute kBot+1 depth has no effect on remapping - if (cell1 <= nCellsAll .and. maxLevelCell(cell1) > 0) then + if (cell1 <= nCells .and. maxLevelCell(cell1) > 0) then heightEdgeNow(:, iEdge) = heightCellNow(kBot+1, cell1) heightEdgeNew(:, iEdge) = heightCellNew(kBot+1, cell1) - elseif (cell2 <= nCellsAll .and. maxLevelCell(cell2) > 0) then + elseif (cell2 <= nCells .and. maxLevelCell(cell2) > 0) then heightEdgeNow(:, iEdge) = heightCellNow(kBot+1, cell2) heightEdgeNew(:, iEdge) = heightCellNew(kBot+1, cell2) else @@ -233,9 +266,6 @@ subroutine ocn_remap_vert_state(block, err) !$omp end do !$omp end parallel - ! Assign new layerThicknesses - layerThickness = layerThicknessNew - !---------------------------------------------------------------------------- ! ACTUAL REMAPPING FROM HERE @@ -330,7 +360,7 @@ subroutine ocn_remap_vert_state(block, err) call mpas_pool_begin_iteration(tracersPool) do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) if ( groupItr % memberType == MPAS_POOL_FIELD ) then - call mpas_pool_get_array(tracersPool, groupItr % memberName, tracersGroup, 2) + call mpas_pool_get_array(tracersPool, groupItr % memberName, tracersGroup, timeLevel) if ( associated(tracersGroup) ) then nTracers = size(tracersGroup, dim=1) @@ -389,6 +419,7 @@ subroutine ocn_remap_vert_state(block, err) deallocate(heightEdgeNow) deallocate(heightCellNew) deallocate(heightEdgeNew) + call mpas_timer_stop("vertical remap") end subroutine ocn_remap_vert_state