Fix bugs related to vertical lagrangian remapping

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",

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.

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

components/mpas-ocean/src/Registry.xml

@@ -3095,7 +3095,7 @@
 			 packages="forwardMode;analysisMode"
 		/>
 		<var name="vertVelocityTop" type="real" dimensions="nVertLevelsP1 nCells Time" units="m s^-1"
-			 description="vertical velocity defined at center (horizontally) and top (vertically) of cell"
+			 description="vertical velocity defined at center (horizontally) and top (vertically) of cell.  This is not the vertical ALE transport, but is Eulerian (fixed-frame) in the vertical, and computed from the continuity equation from the horizontal velocity."
 			 packages="forwardMode;analysisMode"
 		/>
 		<var name="vertTransportVelocityTop" type="real" dimensions="nVertLevelsP1 nCells Time" units="m s^-1"
@@ -3111,6 +3111,9 @@
 			 description="vertical tracer-transport velocity defined at center (horizontally) and top (vertically) of cell.  This is not the vertical ALE transport, but is Eulerian (fixed-frame) in the vertical, and computed from the continuity equation from the horizontal MLE (submesoscale) Bolus velocity."
 			 packages="submeso"
 		/>
+		<var name="layerThicknessTarget" type="real" dimensions="nVertLevels nCells Time" units="m"
+			 description="The layerThickness from ALE routine, prior to remapping."
+		/>
 		<var name="tangentialVelocity" type="real" dimensions="nVertLevels nEdges Time" units="m s^-1"
 			 description="horizontal velocity, tangential to an edge"
 			 missing_value="FILLVAL" missing_value_mask="edgeMask"

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