Support wetting-and-drying with RK4 and subglacial thin film

components/mpas-ocean/src/Registry.xml

@@ -1339,6 +1339,10 @@
 					description="If true, use take into account the ssh gradient across previously flux-limited edges in the wetting and drying algorithm"
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_land_ice_overburden_height" type="real" default_value="0.0" units="m"
+					description="The amount of pressure overburden allowed in grounded regions. If ice must float at level of critical thickness, set equal to -1 * critical thickness"
+					possible_values="Any real"
+		/>
 	</nml_record>
 	<nml_record name="ocean_constants">
 		<nml_option name="config_density0" type="real" default_value="1026.0" units="kg m^-3"

components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F

@@ -76,6 +76,7 @@ module ocn_diagnostics
    !
    !--------------------------------------------------------------------
 
+   logical :: landIcePressureInitialized
    integer :: ke_cell_flag, ke_vertex_flag
    real (kind=RKIND) ::  fCoef
    real (kind=RKIND) ::  landIceTopDragCoeff
@@ -806,6 +807,33 @@ subroutine ocn_diagnostic_solve_wctEdge(btrVelocity, deltaSSH, &
          !$acc exit data copyout(wctEdge) delete(baroclinicThickness, deltaSSH, btrVelocity)
 #endif
 
+      case (thickEdgeFluxConstant)
+
+         ! Use centered (mean) thickness as flux value
+
+         ! Compute mean layer thickness at edge
+#ifdef MPAS_OPENACC
+         !$acc enter data copyin(baroclinicThickness, deltaSSH, wctEdge)
+         !$acc parallel loop &
+         !$acc    present(baroclinicThickness, deltaSSH, wctEdge, cellsOnEdge) &
+         !$acc    private(iEdge, cell1, cell2)
+#else
+         !$omp parallel
+         !$omp do schedule(runtime)
+#endif
+         do iEdge = 1, nEdgesAll
+            if ( maxLevelEdgeTop(iEdge) == 0 ) cycle
+            wctEdge(iEdge) = baroclinicThickness(iEdge)
+         end do
+#ifndef MPAS_OPENACC
+         !$omp end do
+         !$omp end parallel
+#endif
+
+#ifdef MPAS_OPENACC
+         !$acc exit data copyout(wctEdge) delete(baroclinicThickness, deltaSSH)
+#endif
+
       case default
          ! Should have been caught on init
          call mpas_log_write('Thickness flux option unknown', &
@@ -2594,6 +2622,8 @@ subroutine ocn_diagnostic_solve_surface_pressure(forcingPool, &
 
       integer :: iCell ! cell loop index
 
+      real (kind=RKIND) :: h_crit, ssh_crit, landIcePressureMax
+      real (kind=RKIND), dimension(:), allocatable :: landIcePressureApplied
       real (kind=RKIND), dimension(:), pointer :: &
          frazilSurfacePressure, &! frazil ice pressure from forcing
          landIcePressure         ! land ice pressure
@@ -2609,6 +2639,13 @@ subroutine ocn_diagnostic_solve_surface_pressure(forcingPool, &
       call mpas_pool_get_array(forcingPool, 'landIcePressure', &
                                              landIcePressure)
 
+      allocate(landIcePressureApplied(nCellsAll))
+      ! First initialize landIcePressure unchanged, then modify in grounded regions
+      if (landIcePressureOn .and. .not. landIcePressureInitialized) then
+         landIcePressureApplied(:) = landIcePressure(:)
+         landIcePressureInitialized = .true.
+      endif
+
       ! Compute pressure for generalized coordinates.
       ! Pressure at top surface may be due to atmospheric pressure,
       ! sea ice, frazil ice or the weight of an ice shelf
@@ -2652,6 +2689,31 @@ subroutine ocn_diagnostic_solve_surface_pressure(forcingPool, &
       ! Add land ice pressure if needed
 
       if (landIcePressureOn) then
+#ifdef MPAS_OPENACC
+         !$acc parallel loop &
+         !$acc    present(surfacePressure)
+         !! need this eventually, but there is currently
+         !! an issue, probably with the re-retrieval of pointer
+         !!acc    present(surfacePressure, landIcePressure)
+#else
+         !$omp parallel
+         !$omp do schedule(runtime) private(h_crit, ssh_crit, landIcePressureMax)
+#endif
+         do iCell = 1, nCellsAll
+            ! the maximum pressure is the one that results in ssh at the critical
+            ! thickness (gravity is positive), but we also allow some amount of overburden
+            ! because rho may deviate significantly from rho_sw
+            h_crit = config_drying_min_cell_height * (maxLevelCell(iCell) - minLevelCell(iCell) + 1)
+            ssh_crit = bottomDepth(iCell) - h_crit + config_land_ice_overburden_height
+            landIcePressureMax = gravity * rho_sw * ssh_crit
+            landIcePressureApplied(iCell) = &
+               min(landIcePressureMax, landIcePressure(iCell))
+         end do
+#ifndef M_OPENACC
+         !$omp end do
+         !$omp end parallel
+#endif
+
 #ifdef MPAS_OPENACC
          !$acc parallel loop &
          !$acc    present(surfacePressure)
@@ -2663,7 +2725,7 @@ subroutine ocn_diagnostic_solve_surface_pressure(forcingPool, &
 #endif
          do iCell = 1, nCellsAll
             surfacePressure(iCell) = surfacePressure(iCell) &
-                                   + landIcePressure(iCell)
+                                   + landIcePressureApplied(iCell)
          end do
 #ifndef MPAS_OPENACC
          !$omp end do
@@ -4714,6 +4776,8 @@ subroutine ocn_diagnostics_init(domain, err)!{{{
                  & 'Set config_disable_vel_hadv to true for a linear test case. ', &
                  MPAS_LOG_CRIT)
          end if
+         call mpas_log_write('Thickness flux option set to ' //&
+              & trim(config_thickness_flux_type))
       elseif (trim(config_thickness_flux_type) == 'upwind') then
          thickEdgeFluxChoice = thickEdgeFluxUpwind
          call mpas_log_write('Thickness flux option set to ' //&
@@ -4750,6 +4814,12 @@ subroutine ocn_diagnostics_init(domain, err)!{{{
         call mpas_log_write("config_sgr_flux_vertical_location not one of 'bottom', 'uniform', 'top'.", MPAS_LOG_CRIT)
       end if
 
+      if (config_do_restart) then
+         landIcePressureInitialized = .true.
+      else
+         landIcePressureInitialized = .false.
+      endif
+
     end subroutine ocn_diagnostics_init!}}}
 
 !***********************************************************************

components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F

@@ -113,6 +113,7 @@ module ocn_diagnostics_variables
    real (kind=RKIND), dimension(:), pointer :: surfaceBuoyancyForcing
    real (kind=RKIND), dimension(:), pointer :: surfaceFrictionVelocity
    real (kind=RKIND), dimension(:), pointer :: surfacePressure
+   real (kind=RKIND), dimension(:), pointer :: landIcePressureApplied
 
    real (kind=RKIND), dimension(:,:), pointer :: &
       transportVelocityX, transportVelocityY, transportVelocityZ, transportVelocityZonal, &
@@ -714,6 +715,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
       call mpas_pool_get_array(diagnosticsPool, 'dThreshMLD', dThreshMLD)
       call mpas_pool_get_array(diagnosticsPool, 'surfacePressure', &
                 surfacePressure)
+      call mpas_pool_get_array(diagnosticsPool, 'landIcePressureApplied', &
+                landIcePressureApplied)
 
       ! Semi-implicit Array Pointer retrievals
       if (trim(config_time_integrator) == 'split_implicit') then
@@ -966,6 +969,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
       !$acc                   thermExpCoeff,                           &
       !$acc                   sfcFlxAttCoeff,                          &
       !$acc                   surfacePressure,                         &
+      !$acc                   landIcePressureApplied,                  &
       !$acc                   potentialDensity,                        &
       !$acc                   displacedDensity,                        &
       !$acc                   boundaryLayerDepth,                      &
@@ -1228,6 +1232,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
       !$acc                   thermExpCoeff,                           &
       !$acc                   sfcFlxAttCoeff,                          &
       !$acc                   surfacePressure,                         &
+      !$acc                   landIcePressureApplied,                  &
       !$acc                   potentialDensity,                        &
       !$acc                   displacedDensity,                        &
       !$acc                   boundaryLayerDepth,                      &
@@ -1433,6 +1438,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
               thermExpCoeff,                           &
               sfcFlxAttCoeff,                          &
               surfacePressure,                         &
+              landIcePressureApplied,                  &
               potentialDensity,                        &
               displacedDensity,                        &
               boundaryLayerDepth,                      &