Add wettingVelocityFactor ramp as in O'Dea

fixup ramp
E3SM-Ocean-Discussion · Jul 13, 2022 · fd86b92 · fd86b92
1 parent 412b4d5
commit fd86b92
Show file tree

Hide file tree

Showing 2 changed files with 81 additions and 38 deletions.
diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
@@ -1089,6 +1089,10 @@
 					description="If true, just zero out velocity that is contributing to drying for cell that is drying.  This option can be used to estimate acceptable minimum thicknesses for a run."
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_zero_drying_velocity_ramp" type="logical" default_value=".false." units="unitless"
+					description="If true, ramp velocities and tendencies to zero rather than applying a simple on/off switch."
+					possible_values=".true. or .false."
+		/>
 		<nml_option name="config_verify_not_dry" type="logical" default_value=".false." units="unitless"
 					description="If true, verify that cells are at least config_min_cell_height thick."
 					possible_values=".true. or .false."

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F
@@ -235,6 +235,7 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying
       type (mpas_pool_type), pointer :: tendPool
       type (mpas_pool_type), pointer :: statePool
       type (mpas_pool_type), pointer :: provisStatePool
+      real (kind=RKIND), dimension(:), pointer :: ssh
       real (kind=RKIND), dimension(:, :), pointer :: layerThicknessCur
       real (kind=RKIND), dimension(:, :), pointer :: layerThicknessProvis
       real (kind=RKIND), dimension(:, :), pointer :: normalVelocity
@@ -247,6 +248,7 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying
      call mpas_pool_get_subpool(block % structs, 'state', statePool)
      call mpas_pool_get_subpool(block % structs, 'provis_state', provisStatePool)
 
+     call mpas_pool_get_array(statePool, 'ssh', ssh, 1)
      call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 1)
      ! use thickness at n because constraint is h_n + dt*T_h > h_min
      call mpas_pool_get_array(statePool, 'layerThickness', layerThicknessCur, 1)
@@ -261,10 +263,12 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying
      !$omp end do
      !$omp end parallel
 
-     ! ensure cells stay wet by selectively damping cells with a damping tendency to make sure tendency doesn't dry cells
+     ! ensure cells stay wet by selectively damping cells with a damping tendency to make 
+     ! sure tendency doesn't dry cells
 
-     call ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
-                                             normalTransportVelocity, rkSubstepWeight, wettingVelocityFactor, err)
+     call ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, &
+                                             layerThicknessProvis, normalTransportVelocity, &
+                                             ssh, rkSubstepWeight, wettingVelocityFactor, err)
 
      ! prevent drying from happening with selective wettingVelocityFactor
      if (config_zero_drying_velocity) then
@@ -274,8 +278,10 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying
          do k = minLevelEdgeTop(iEdge), maxLevelEdgeBot(iEdge)
 
            if (abs(wettingVelocityFactor(k, iEdge)) > 0.0_RKIND) then
-             normalTransportVelocity(k, iEdge) = 0.0_RKIND
-             normalVelocity(k, iEdge) = 0.0_RKIND
+             normalTransportVelocity(k, iEdge) = (1.0_RKIND - &
+               wettingVelocityFactor(k, iEdge)) * normalTransportVelocity(k, iEdge)
+             normalVelocity(k, iEdge) = (1.0_RKIND - &
+               wettingVelocityFactor(k, iEdge)) * normalVelocity(k, iEdge)
            end if
 
          end do
@@ -301,7 +307,7 @@ end subroutine ocn_prevent_drying_rk4 !}}}
 !-----------------------------------------------------------------------
    subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
 
-       normalVelocity, dt, wettingVelocityFactor, err)!{{{
+       normalVelocity, ssh, dt, wettingVelocityFactor, err)!{{{
 
       !-----------------------------------------------------------------
       !
@@ -321,6 +327,8 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness
       real (kind=RKIND), dimension(:,:), intent(in) :: &
          normalVelocity     !< Input: transport
 
+      real (kind=RKIND), dimension(:), intent(in) :: ssh
+
       real (kind=RKIND), intent(in) :: &
          dt     !< Input: time step
 
@@ -347,55 +355,86 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness
       !
       !-----------------------------------------------------------------
 
-      integer :: iEdge, iCell, k, i
+      integer :: cell1, cell2, iEdge, iCell, k, i
 
       real (kind=RKIND) :: divOutFlux
       real (kind=RKIND) :: layerThickness
+      real (kind=RKIND) :: hCrit, hEdgeTotal
+
+      character (len=100) :: log_string
 
       err = 0
 
-      if (.not. config_zero_drying_velocity ) return
+      if (.not. config_zero_drying_velocity) return
 
-      ! need predicted transport velocity to limit drying flux
-      !$omp parallel
-      !$omp do schedule(runtime) private(i, iEdge, k, divOutFlux, layerThickness)
-      do iCell = 1, nCellsAll
-        ! can switch with maxLevelEdgeBot(iEdge)
-        do k = minLevelCell(iCell), maxLevelCell(iCell)
-          divOutFlux = 0.0_RKIND
-          layerThickness = min(layerThicknessProvis(k, iCell), layerThicknessCur(k, iCell))
-          do i = 1, nEdgesOnCell(iCell)
-            iEdge = edgesOnCell(i, iCell)
-            if (k <= maxLevelEdgeBot(iEdge) .or. k >= minLevelEdgeTop(iEdge)) then
-              ! only consider divergence flux leaving the cell
-              if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) < 0.0_RKIND ) then
-                divOutFlux = divOutFlux &
-                             + normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) * &
-                               layerThickEdgeFlux(k, iEdge) * dvEdge(iEdge) * &
-                               invAreaCell(iCell)
-              end if
-            end if
-          end do
+      if (config_zero_drying_velocity_ramp) then
+
+        ! Following O'Dea et al. (2021), if total upwinded wct is less than
+        ! 2*critical thickness, apply damping throughout the water column
+
+        !$omp parallel
+        !$omp do schedule(runtime) private(cell1, cell2, iEdge, k, hEdgeTotal)
+        do iEdge = 1, nEdgesAll
+
+          cell1 = cellsOnEdge(1,iEdge)
+          cell2 = cellsOnEdge(2,iEdge)
+          hEdgeTotal = 0.0_RKIND
+          hCrit = 0.0_RKIND
+          do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
+             hEdgeTotal = hEdgeTotal + layerThickEdgeFlux(k, iEdge)
+             hCrit = hCrit + config_drying_min_cell_height
+          enddo
 
-          ! if layer thickness is too small, limit divergence flux outwards with opposite velocity
-          if ((layerThickness + dt * divOutFlux ) <= &
-              (config_drying_min_cell_height + config_drying_safety_height)) then
+          if (hEdgeTotal <= hCrit) then
+             wettingVelocityFactor(:, iEdge) = 1.0_RKIND
+          elseif (hEdgeTotal > hCrit .and. hEdgeTotal <= 2.0_RKIND*hCrit) then
+             wettingVelocityFactor(:, iEdge) = 1.0_RKIND - &
+               TANH(50.0_RKIND * (hEdgeTotal - hCrit)/hCrit)
+          endif
+        enddo
+        !$omp end do
+        !$omp end parallel
 
+      else
+
+        ! need predicted transport velocity to limit drying flux
+        !$omp parallel
+        !$omp do schedule(runtime) private(i, iEdge, k, divOutFlux, layerThickness)
+        do iCell = 1, nCellsAll
+          do k = minLevelCell(iCell), maxLevelCell(iCell)
+            divOutFlux = 0.0_RKIND
+            layerThickness = min(layerThicknessProvis(k, iCell), layerThicknessCur(k, iCell))
             do i = 1, nEdgesOnCell(iCell)
               iEdge = edgesOnCell(i, iCell)
               if (k <= maxLevelEdgeBot(iEdge) .or. k >= minLevelEdgeTop(iEdge)) then
-                if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then
-                  ! just go with simple boolean approach for zero wetting velocity for debugging purposes
-                  wettingVelocityFactor(k, iEdge) = 1.0_RKIND
+                ! only consider divergence flux leaving the cell
+                if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) < 0.0_RKIND ) then
+                  divOutFlux = divOutFlux + &
+                               normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) * &
+                               layerThickEdgeFlux(k, iEdge) * dvEdge(iEdge) * &
+                               invAreaCell(iCell)
                 end if
               end if
             end do
-          end if
 
+            ! if layer thickness is too small, limit divergence flux outwards with opposite velocity
+            if ((layerThickness + dt * divOutFlux) <= &
+                (config_drying_min_cell_height + config_drying_safety_height)) then
+              do i = 1, nEdgesOnCell(iCell)
+                iEdge = edgesOnCell(i, iCell)
+                if (k <= maxLevelEdgeBot(iEdge) .or. k >= minLevelEdgeTop(iEdge)) then
+                  if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then
+                    wettingVelocityFactor(k, iEdge) = 1.0_RKIND
+                  end if
+                end if
+              end do
+            end if
+
+          end do
         end do
-      end do
-      !$omp end do
-      !$omp end parallel
+        !$omp end do
+        !$omp end parallel
+      end if
 
    end subroutine ocn_wetting_drying_wettingVelocity !}}}