Change wettingVelocity to wettingVelocityFactor

E3SM-Ocean-Discussion · Apr 11, 2022 · 0a4d363 · 0a4d363
1 parent 4e95577
commit 0a4d363
Show file tree

Hide file tree

Showing 5 changed files with 21 additions and 33 deletions.
diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
@@ -2537,8 +2537,8 @@
 			 description="horizontal velocity used to transport mass and tracers"
 			 packages="forwardMode;analysisMode"
 		/>
-		<var name="wettingVelocity" type="real" dimensions="nVertLevels nEdges Time" units="m s^{-1}"
-			 description="Velocity to prevent drying of cell."
+		<var name="wettingVelocityFactor" type="real" dimensions="nVertLevels nEdges Time" units="unitless"
+			 description="Scaling factor for normalVelocity and its tendency to prevent drying of cell between 0 (no scaling) and 1 (normalVelocity and tendency set to 0)."
 			 packages="forwardMode"
 		/>
 		<var name="vertAleTransportTop" type="real" dimensions="nVertLevelsP1 nCells Time" units="m s^{-1}"

diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F
@@ -1205,7 +1205,7 @@ subroutine ocn_time_integrator_rk4_diagnostic_update(block, dt, rkSubstepWeight,
          do k = 1, maxLevelEdgeTop(iEdge)
             normalVelocityProvis(k, iEdge) = normalVelocityCur(k, iEdge) + rkSubstepWeight &
                                            * normalVelocityTend(k, iEdge)
-            normalVelocityProvis(k, iEdge) = normalVelocityProvis(k, iEdge) * (1.0_RKIND - wettingVelocity(k, iEdge))
+            normalVelocityProvis(k, iEdge) = normalVelocityProvis(k, iEdge) * (1.0_RKIND - wettingVelocityFactor(k, iEdge))
          end do
       end do
       !$omp end do
@@ -1441,7 +1441,7 @@ subroutine ocn_time_integrator_rk4_accumulate_update(block, rkWeight, err)!{{{
       !$omp do schedule(runtime)
       do iEdge = 1, nEdges
          normalVelocityNew(:, iEdge) = normalVelocityNew(:, iEdge) + rkWeight * normalVelocityTend(:, iEdge)
-         normalVelocityNew(:, iEdge) = normalVelocityNew(:, iEdge) * (1.0_RKIND - wettingVelocity(:, iEdge))
+         normalVelocityNew(:, iEdge) = normalVelocityNew(:, iEdge) * (1.0_RKIND - wettingVelocityFactor(:, iEdge))
       end do
       !$omp end do
       !$omp end parallel

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
@@ -131,7 +131,7 @@ module ocn_diagnostics_variables
    real (kind=RKIND), dimension(:,:), pointer :: viscosity
    real (kind=RKIND), dimension(:,:), pointer :: vertViscTopOfEdge
    real (kind=RKIND), dimension(:,:), pointer :: vertViscTopOfCell
-   real (kind=RKIND), dimension(:,:), pointer :: wettingVelocity
+   real (kind=RKIND), dimension(:,:), pointer :: wettingVelocityFactor
    type (field2DReal), pointer :: wettingVelocityField
 
    real (kind=RKIND), dimension(:,:), pointer :: vertDiffTopOfCell
@@ -239,9 +239,9 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
       end if
 
       if ( trim(config_ocean_run_mode) == 'forward' ) then
-         call mpas_pool_get_array(diagnosticsPool, 'wettingVelocity', &
-                   wettingVelocity)
-         call mpas_pool_get_field(diagnosticsPool, 'wettingVelocity', &
+         call mpas_pool_get_array(diagnosticsPool, 'wettingVelocityFactor', &
+                   wettingVelocityFactor)
+         call mpas_pool_get_field(diagnosticsPool, 'wettingVelocityFactor', &
                    wettingVelocityField)
       end if
 
@@ -734,8 +734,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
          !$acc                   )
       end if
       if ( trim(config_ocean_run_mode) == 'forward' ) then
-         !$acc enter data create( &
-         !$acc                   wettingVelocity              &
+         !$acc enter data copyin( &
+         !$acc                   wettingVelocityFactor            &
          !$acc                   )
       end if
       if ( config_use_GM ) then
@@ -973,7 +973,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
       end if
       if ( trim(config_ocean_run_mode) == 'forward' ) then
          !$acc exit data delete(                              &
-         !$acc                   wettingVelocity              &
+         !$acc                   wettingVelocityFactor        &
          !$acc                   )
       end if
       if ( config_use_GM ) then
@@ -1170,7 +1170,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
                  normalizedRelativeVorticityCell)
       end if
       if ( trim(config_ocean_run_mode) == 'forward' ) then
-         nullify(wettingVelocity)
+         nullify(wettingVelocityFactor)
       end if
       if ( config_use_GM ) then
          nullify(RediKappa,         &

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F
@@ -346,7 +346,7 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, &
       !$acc               density, normRelVortEdge, montgomeryPotential, pressure, &
       !$acc               thermExpCoeff, salineContractCoeff, tangentialVelocity, &
       !$acc               layerThickEdgeFlux, kineticEnergyCell, sfcFlxAttCoeff, potentialDensity, &
-      !$acc               vertAleTransportTop, vertViscTopOfEdge, wettingVelocity)
+      !$acc               vertAleTransportTop, vertViscTopOfEdge, wettingVelocityFactor)
       !$acc enter data &
       !$acc    copyin(tendVel, sfcStress, sfcStressMag, &
       !$acc    ssh, normalVelocity, &
@@ -437,15 +437,14 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, &
 
 #ifdef MPAS_OPENACC
       !$acc parallel loop collapse(2) &
-      !$acc    present(tendVel, wettingVelocity)
+      !$acc    present(tendVel, wettingVelocityFactor)
 #else
       !$omp parallel
       !$omp do schedule(runtime) private(k)
 #endif
       do iEdge = 1, nEdgesAll
       do k=1,nVertLevels
-         tendVel(k,iEdge) = tendVel(k,iEdge)* &
-                            (1.0_RKIND - wettingVelocity(k,iEdge))
+         tendVel(k,iEdge) = tendVel(k,iEdge) * (1.0_RKIND - wettingVelocityFactor(k,iEdge))
       end do
       end do
 #ifndef MPAS_OPENACC

diff --git a/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F b/components/mpas-ocean/src/shared/mpas_ocn_wetting_drying.F
@@ -253,7 +253,7 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying
      !$omp parallel
      !$omp do schedule(runtime)
      do iEdge = 1, nEdgesAll
-       wettingVelocity(:, iEdge) = 0.0_RKIND
+       wettingVelocityFactor(:, iEdge) = 0.0_RKIND
      end do
      !$omp end do
      !$omp end parallel
@@ -263,21 +263,13 @@ subroutine ocn_prevent_drying_rk4(block, dt, rkSubstepWeight, config_zero_drying
      call ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
                                              normalTransportVelocity, rkSubstepWeight, wettingVelocityFactor, err)
 
-     ! prevent drying from happening with selective wettingVelocity
+     ! prevent drying from happening with selective wettingVelocityFactor
      !$omp parallel
      !$omp do schedule(runtime) private(k)
      do iEdge = 1, nEdgesAll
        do k = minLevelEdgeTop(iEdge), maxLevelEdgeBot(iEdge)
-         if (abs(normalTransportVelocity(k,iEdge) + wettingVelocity(k,iEdge)) < eps)  then
-           ! prevent spurious flux for close to zero values
-           normalTransportVelocity(k, iEdge) = 0.0_RKIND
-           normalVelocity(k, iEdge) = 0.0_RKIND
-         else if (abs(normalTransportVelocity(k,iEdge) + wettingVelocity(k,iEdge)) <= abs(normalTransportVelocity(k,iEdge))) then
-           normalTransportVelocity(k, iEdge) = normalTransportVelocity(k, iEdge) + wettingVelocity(k, iEdge)
-           normalVelocity(k, iEdge) = normalVelocity(k, iEdge) + wettingVelocity(k, iEdge)
-         end if
 
-         if (abs(wettingVelocity(k, iEdge)) > 0.0_RKIND .and. config_zero_drying_velocity) then
+         if (abs(wettingVelocityFactor(k, iEdge)) > 0.0_RKIND .and. config_zero_drying_velocity) then
            normalTransportVelocity(k, iEdge) = 0.0_RKIND
            normalVelocity(k, iEdge) = 0.0_RKIND
          end if
@@ -304,7 +296,7 @@ end subroutine ocn_prevent_drying_rk4 !}}}
 !-----------------------------------------------------------------------
    subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThicknessCur, layerThicknessProvis, &
 
-       normalVelocity, dt, wettingVelocity, err)!{{{
+       normalVelocity, dt, wettingVelocityFactor, err)!{{{
 
       !-----------------------------------------------------------------
       !
@@ -334,7 +326,7 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness
       !-----------------------------------------------------------------
 
       real (kind=RKIND), dimension(:,:), intent(inout) :: &
-         wettingVelocity          !< Input/Output: velocity wettingVelocityency
+         wettingVelocityFactor          !< Input/Output: velocity wettingVelocityFactor
 
       !-----------------------------------------------------------------
       !
@@ -386,12 +378,9 @@ subroutine ocn_wetting_drying_wettingVelocity(layerThickEdgeFlux, layerThickness
               iEdge = edgesOnCell(i, iCell)
               if (k <= maxLevelEdgeBot(iEdge) .or. k >= minLevelEdgeTop(iEdge)) then
                 if ( normalVelocity(k, iEdge) * edgeSignOnCell(i, iCell) <= 0.0_RKIND ) then
-                  ! each outgoing velocity is penalized (but not the incoming, wetting velocities)
-                  ! square the fractional term to make values near zero go to zero much quicker (to prevent threshold from being hit)
-                  wettingVelocity(k, iEdge) =  - (min(max(0.0_RKIND, 1.0_RKIND - (divOutFlux*divOutFlux)), 1.0_RKIND)) * normalVelocity(k, iEdge)
                   ! just go with simple boolean approach for zero wetting velocity for debugging purposes
                   if (config_zero_drying_velocity) then
-                    wettingVelocity(k, iEdge) = 1.0_RKIND
+                    wettingVelocityFactor(k, iEdge) = 1.0_RKIND
                   end if
                 end if
               end if