Hydrostatic initialization in ice cavities (NOAA-EMC#41)

* Hydrostatic initialization in ice cavities

    - Iteratively solve for the initial ice shelf displacement
      in cavities by calculating the pressure at the current
      displacement depth using the unperturbed profile.
    - This change should obsolete TRIM_IC_FOR_PSURF and DEPRESS_INITIAL_SURFACE for
      ice shelf applications and should work for arbitrary equations of state.
    - Existing implementations (e.g. ISOMIP) should turn off the above options
      in order to exercise this feature.
    - This code change should not impact non ice-shelf configurations
      or those with either of the above two options.

* Addresses a number of issues identified in code review.

* add appropriate intent to ice_shelf_query

* fix unit scaling comments

Co-authored-by: Marshall Ward <marshall.ward@noaa.gov>

Author: MJHarrison-GFDL

src/core/MOM.F90

@@ -283,7 +283,8 @@ module MOM
   type(time_type) :: dtbt_reset_interval !< A time_time representation of dtbt_reset_period.
   type(time_type) :: dtbt_reset_time !< The next time DTBT should be calculated.
   real, dimension(:,:), pointer :: frac_shelf_h => NULL() !< fraction of total area occupied
-                                     !! by ice shelf [nondim]
+  !! by ice shelf [nondim]
+  real, dimension(:,:), pointer :: mass_shelf => NULL() !< Mass of ice shelf [R Z ~> kg m-2]
   real, dimension(:,:,:), pointer :: &
     h_pre_dyn => NULL(), &      !< The thickness before the transports [H ~> m or kg m-2].
     T_pre_dyn => NULL(), &      !< Temperature before the transports [degC].
@@ -1748,9 +1749,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   real, allocatable         :: v_in(:,:,:) ! Initial meridional velocities [L T-1 ~> m s-1]
   real, allocatable         :: h_in(:,:,:) ! Initial layer thicknesses [H ~> m or kg m-2]
   real, allocatable, target :: frac_shelf_in(:,:) ! Initial fraction of the total cell area occupied
-                                           ! by an ice shelf [nondim]
+                                                  ! by an ice shelf [nondim]
+  real, allocatable, target :: mass_shelf_in(:,:) ! Initial mass of ice shelf contained within a grid cell
+                                                  ! [R Z ~> kg m-2]
   real, allocatable, target :: T_in(:,:,:) ! Initial temperatures [degC]
   real, allocatable, target :: S_in(:,:,:) ! Initial salinities [ppt]
+
   type(ocean_OBC_type), pointer :: OBC_in => NULL()
   type(sponge_CS), pointer :: sponge_in_CSp => NULL()
   type(ALE_sponge_CS), pointer :: ALE_sponge_in_CSp => NULL()
@@ -2523,14 +2527,17 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
       ! for legacy reasons. The actual ice shelf diag CS is internal to the ice shelf
       call initialize_ice_shelf(param_file, G_in, Time, ice_shelf_CSp, diag_ptr)
       allocate(frac_shelf_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed), source=0.0)
+      allocate(mass_shelf_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed), source=0.0)
       allocate(CS%frac_shelf_h(isd:ied, jsd:jed), source=0.0)
-      call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h)
+      allocate(CS%mass_shelf(isd:ied, jsd:jed), source=0.0)
+      call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h, CS%mass_shelf)
       ! MOM_initialize_state is using the  unrotated metric
       call rotate_array(CS%frac_shelf_h, -turns, frac_shelf_in)
+      call rotate_array(CS%mass_shelf, -turns, mass_shelf_in)
       call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, &
           param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, &
           sponge_in_CSp, ALE_sponge_in_CSp, oda_incupd_in_CSp, OBC_in, Time_in, &
-          frac_shelf_h=frac_shelf_in)
+          frac_shelf_h=frac_shelf_in, mass_shelf = mass_shelf_in)
     else
       call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, &
           param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, &
@@ -2574,16 +2581,17 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
       deallocate(S_in)
     endif
     if (use_ice_shelf) &
-      deallocate(frac_shelf_in)
+      deallocate(frac_shelf_in,mass_shelf_in)
   else
     if (use_ice_shelf) then
       call initialize_ice_shelf(param_file, G, Time, ice_shelf_CSp, diag_ptr)
       allocate(CS%frac_shelf_h(isd:ied, jsd:jed), source=0.0)
-      call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h)
+      allocate(CS%mass_shelf(isd:ied, jsd:jed), source=0.0)
+      call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h, CS%mass_shelf)
       call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, &
           param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, &
           CS%sponge_CSp, CS%ALE_sponge_CSp,CS%oda_incupd_CSp, CS%OBC, Time_in, &
-          frac_shelf_h=CS%frac_shelf_h)
+          frac_shelf_h=CS%frac_shelf_h, mass_shelf=CS%mass_shelf)
     else
       call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, &
           param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, &
@@ -2598,8 +2606,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
     endif
   endif
 
-  if (use_ice_shelf .and. CS%debug) &
+  if (use_ice_shelf .and. CS%debug) then
     call hchksum(CS%frac_shelf_h, "MOM:frac_shelf_h", G%HI, haloshift=0)
+    call hchksum(CS%mass_shelf, "MOM:mass_shelf", G%HI, haloshift=0,scale=US%RZ_to_kg_m2)
+  endif
 
   call cpu_clock_end(id_clock_MOM_init)
   call callTree_waypoint("returned from MOM_initialize_state() (initialize_MOM)")

src/ice_shelf/MOM_ice_shelf.F90

@@ -2032,11 +2032,12 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time)
 end subroutine update_shelf_mass
 
 !> Save the ice shelf restart file
-subroutine ice_shelf_query(CS, G, frac_shelf_h)
+subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf)
   type(ice_shelf_CS),         pointer    :: CS !< ice shelf control structure
   type(ocean_grid_type), intent(in)      :: G  !< A pointer to an ocean grid control structure.
-  real, optional, dimension(SZI_(G),SZJ_(G))  :: frac_shelf_h !<
-                                      !< Ice shelf area fraction [nodim].
+  real, optional, dimension(SZI_(G),SZJ_(G)), intent(out)  :: frac_shelf_h !< Ice shelf area fraction [nodim].
+  real, optional, dimension(SZI_(G),SZJ_(G)), intent(out)  :: mass_shelf !<Ice shelf mass [R Z -> kg m-2]
+
 
   integer :: i, j
 
@@ -2047,6 +2048,13 @@ subroutine ice_shelf_query(CS, G, frac_shelf_h)
     enddo ; enddo
   endif
 
+  if (present(mass_shelf)) then
+    do j=G%jsd,G%jed ; do i=G%isd,G%ied
+      mass_shelf(i,j) = 0.0
+      if (G%areaT(i,j)>0.) mass_shelf(i,j) = CS%ISS%mass_shelf(i,j)
+    enddo ; enddo
+  endif
+
 end subroutine ice_shelf_query
 
 !> Save the ice shelf restart file

src/initialization/MOM_state_initialization.F90

@@ -117,7 +117,7 @@ module MOM_state_initialization
 !! conditions or by reading them from a restart (or saves) file.
 subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
                                 restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, &
-                                ALE_sponge_CSp, oda_incupd_CSp, OBC, Time_in, frac_shelf_h)
+                                ALE_sponge_CSp, oda_incupd_CSp, OBC, Time_in, frac_shelf_h, mass_shelf)
   type(ocean_grid_type),      intent(inout) :: G    !< The ocean's grid structure.
   type(verticalGrid_type),    intent(in)    :: GV   !< The ocean's vertical grid structure.
   type(unit_scale_type),      intent(in)    :: US   !< A dimensional unit scaling type
@@ -147,6 +147,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
   real, dimension(SZI_(G),SZJ_(G)), &
                      optional, intent(in)   :: frac_shelf_h    !< The fraction of the grid cell covered
                                                                !! by a floating ice shelf [nondim].
+  real, dimension(SZI_(G),SZJ_(G)), &
+                     optional, intent(in)   :: mass_shelf      !< The mass per unit area of the overlying
+                                                               !! ice shelf [ R Z ~> kg m-2 ]
   ! Local variables
   real :: depth_tot(SZI_(G),SZJ_(G))  ! The nominal total depth of the ocean [Z ~> m]
   character(len=200) :: filename   ! The name of an input file.
@@ -158,6 +161,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
   real :: vel_rescale ! A rescaling factor for velocities from the representation in
                       ! a restart file to the internal representation in this run.
   real :: dt          ! The baroclinic dynamics timestep for this run [T ~> s].
+
   logical :: from_Z_file, useALE
   logical :: new_sim
   integer :: write_geom
@@ -404,6 +408,23 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
   if (use_temperature .and. use_OBC) &
     call fill_temp_salt_segments(G, GV, OBC, tv)
 
+  ! Calculate the initial surface displacement under ice shelf
+
+  call get_param(PF, mdl, "DEPRESS_INITIAL_SURFACE", depress_sfc, &
+       "If true,  depress the initial surface to avoid huge "//&
+       "tsunamis when a large surface pressure is applied.", &
+       default=.false., do_not_log=just_read)
+  call get_param(PF, mdl, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, &
+       "If true, cuts way the top of the column for initial conditions "//&
+       "at the depth where the hydrostatic pressure matches the imposed "//&
+       "surface pressure which is read from file.", default=.false., &
+       do_not_log=just_read)
+
+  if (new_sim) then
+    if (use_ice_shelf .and. present(mass_shelf) .and. .not. (trim_ic_for_p_surf .or. depress_sfc)) &
+         call calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h)
+  endif
+
   ! The thicknesses in halo points might be needed to initialize the velocities.
   if (new_sim) call pass_var(h, G%Domain)
 
@@ -458,15 +479,6 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
     call convert_thickness(h, G, GV, US, tv)
 
   ! Remove the mass that would be displaced by an ice shelf or inverse barometer.
-  call get_param(PF, mdl, "DEPRESS_INITIAL_SURFACE", depress_sfc, &
-               "If true,  depress the initial surface to avoid huge "//&
-               "tsunamis when a large surface pressure is applied.", &
-               default=.false., do_not_log=just_read)
-  call get_param(PF, mdl, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, &
-               "If true, cuts way the top of the column for initial conditions "//&
-               "at the depth where the hydrostatic pressure matches the imposed "//&
-               "surface pressure which is read from file.", default=.false., &
-               do_not_log=just_read)
   if (depress_sfc .and. trim_ic_for_p_surf) call MOM_error(FATAL, "MOM_initialize_state: "//&
            "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True")
   if (new_sim .and. debug .and. (depress_sfc .or. trim_ic_for_p_surf)) &
@@ -1035,7 +1047,7 @@ subroutine convert_thickness(h, G, GV, US, tv)
 end subroutine convert_thickness
 
 !> Depress the sea-surface based on an initial condition file
-subroutine depress_surface(h, G, GV, US, param_file, tv, just_read)
+subroutine depress_surface(h, G, GV, US, param_file, tv, just_read, z_top_shelf)
   type(ocean_grid_type),   intent(in)    :: G    !< The ocean's grid structure
   type(verticalGrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure
   type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
@@ -1045,6 +1057,8 @@ subroutine depress_surface(h, G, GV, US, param_file, tv, just_read)
   type(thermo_var_ptrs),   intent(in)    :: tv   !< A structure pointing to various thermodynamic variables
   logical,                 intent(in)    :: just_read !< If true, this call will only read
                                                       !! parameters without changing h.
+  real, dimension(SZI_(G),SZJ_(G)), &
+                 optional, intent(in)    :: z_top_shelf    !< Top interface position under ice shelf [Z ~> m]
   ! Local variables
   real, dimension(SZI_(G),SZJ_(G)) :: &
     eta_sfc  ! The free surface height that the model should use [Z ~> m].
@@ -1057,30 +1071,40 @@ subroutine depress_surface(h, G, GV, US, param_file, tv, just_read)
   character(len=200) :: inputdir, eta_srf_file ! Strings for file/path
   character(len=200) :: filename, eta_srf_var  ! Strings for file/path
   integer :: i, j, k, is, ie, js, je, nz
+  logical :: use_z_shelf
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
-  ! Read the surface height (or pressure) from a file.
+  use_z_shelf = present(z_top_shelf)
 
-  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
-  inputdir = slasher(inputdir)
-  call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_FILE", eta_srf_file,&
-                 "The initial condition file for the surface height.", &
-                 fail_if_missing=.not.just_read, do_not_log=just_read)
-  call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_VAR", eta_srf_var, &
-                 "The initial condition variable for the surface height.",&
-                 default="SSH", do_not_log=just_read)
-  filename = trim(inputdir)//trim(eta_srf_file)
-  if (.not.just_read) &
-    call log_param(param_file,  mdl, "INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
 
-  call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_SCALE", scale_factor, &
-                 "A scaling factor to convert SURFACE_HEIGHT_IC_VAR into units of m", &
-                 units="variable", default=1.0, scale=US%m_to_Z, do_not_log=just_read)
+  if (.not. use_z_shelf) then
+  ! Read the surface height (or pressure) from a file.
 
-  if (just_read) return ! All run-time parameters have been read, so return.
+    call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
+    inputdir = slasher(inputdir)
+    call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_FILE", eta_srf_file,&
+                   "The initial condition file for the surface height.", &
+                   fail_if_missing=.not.just_read, do_not_log=just_read)
+    call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_VAR", eta_srf_var, &
+                   "The initial condition variable for the surface height.",&
+                   default="SSH", do_not_log=just_read)
+    filename = trim(inputdir)//trim(eta_srf_file)
+    if (.not.just_read) &
+      call log_param(param_file,  mdl, "INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
+
+    call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_SCALE", scale_factor, &
+                   "A scaling factor to convert SURFACE_HEIGHT_IC_VAR into units of m", &
+                   units="variable", default=1.0, scale=US%m_to_Z, do_not_log=just_read)
+
+    if (just_read) return ! All run-time parameters have been read, so return.
 
-  call MOM_read_data(filename, eta_srf_var, eta_sfc, G%Domain, scale=scale_factor)
+    call MOM_read_data(filename, eta_srf_var, eta_sfc, G%Domain, scale=scale_factor)
+  else
+    do j=js,je ; do i=is,ie
+      eta_sfc(i,j) = z_top_shelf(i,j)
+    enddo; enddo
+  endif
 
   ! Convert thicknesses to interface heights.
   call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref)
@@ -1201,6 +1225,88 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read)
 
 end subroutine trim_for_ice
 
+!> Calculate the hydrostatic equilibrium position of the surface under an ice shelf
+subroutine calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h)
+  type(param_file_type),   intent(in)    :: PF !< Parameter file structure
+  type(ocean_grid_type),   intent(in)    :: G  !< Ocean grid structure
+  type(verticalGrid_type), intent(in)    :: GV !< Vertical grid structure
+  type(unit_scale_type),   intent(in)    :: US !< A dimensional unit scaling type
+  real, dimension(SZI_(G),SZJ_(G)), &
+                           intent(in)    :: mass_shelf  !< Ice shelf mass [R Z ~> kg m-2]
+  type(thermo_var_ptrs),   intent(inout) :: tv !< Thermodynamics structure
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                           intent(inout) :: h  !< Layer thickness [H ~> m or kg m-2]
+
+  real :: z_top_shelf(SZI_(G),SZJ_(G))  ! The depth of the top interface under ice shelves [Z ~> m]
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
+                                   eta  ! The free surface height that the model should use [Z ~> m].
+  ! temporary arrays
+  real, dimension(SZK_(GV)) :: rho_col   ! potential density in the column for use in ice
+  real, dimension(SZK_(GV)) :: rho_h     ! potential density multiplied by thickness [R Z ~> kg m-2 ]
+  real, dimension(SZK_(GV)) :: h_tmp     ! temporary storage for thicknesses [H ~> m]
+  real, dimension(SZK_(GV)) :: p_ref     ! pressure for density [R Z ~> kg m-2]
+  real, dimension(SZK_(GV)+1) :: ei_tmp, ei_orig ! temporary storage for interface positions [Z ~> m]
+  real :: z_top, z_col, mass_disp, residual, tol
+  integer :: is, ie, js, je, k, nz, i, j, max_iter, iter
+
+  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+
+
+  tol = 0.001 ! The initialization tolerance for ice shelf initialization (m)
+  call get_param(PF, mdl, "ICE_SHELF_INITIALIZATION_Z_TOLERANCE", tol, &
+                "A initialization tolerance for the calculation of the static "// &
+                "ice shelf displacement (m) using initial temperature and salinity profile.",&
+                 default=tol, units="m", scale=US%m_to_Z)
+  max_iter = 1e3
+  call MOM_mesg("Started calculating initial interface position under ice shelf ")
+  ! Convert thicknesses to interface heights.
+  call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref)
+  do j = js, je ; do i = is, ie
+    iter = 1
+    z_top_shelf(i,j) = 0.0
+    p_ref(:) = tv%p_ref
+    if (G%mask2dT(i,j) .gt. 0. .and. mass_shelf(i,j) .gt. 0.) then
+      call calculate_density(tv%T(i,j,:), tv%S(i,j,:), P_Ref, rho_col, tv%eqn_of_state)
+      z_top = min(max(-1.0*mass_shelf(i,j)/rho_col(1),-G%bathyT(i,j)),0.)
+      h_tmp = 0.0
+      z_col = 0.0
+      ei_tmp(1:nz+1)=eta(i,j,1:nz+1)
+      ei_orig(1:nz+1)=eta(i,j,1:nz+1)
+      do k=1,nz+1
+        if (ei_tmp(k)<z_top) ei_tmp(k)=z_top
+      enddo
+      mass_disp = 0.0
+      do k=1,nz
+        h_tmp(k) = max(ei_tmp(k)-ei_tmp(k+1),GV%Angstrom_H)
+        rho_h(k) = h_tmp(k) * rho_col(k)
+        mass_disp = mass_disp + rho_h(k)
+      enddo
+      residual = mass_shelf(i,j) - mass_disp
+      do while (abs(residual) .gt. tol .and. z_top .gt. -G%bathyT(i,j) .and. iter .lt. max_iter)
+        z_top=min(max(z_top-(residual*0.5e-3),-G%bathyT(i,j)),0.0)
+        h_tmp = 0.0
+        z_col = 0.0
+        ei_tmp(1:nz+1) = ei_orig(1:nz+1)
+        do k=1,nz+1
+          if (ei_tmp(k)<z_top) ei_tmp(k)=z_top
+        enddo
+        mass_disp = 0.0
+        do k=1,nz
+          h_tmp(k) = max(ei_tmp(k)-ei_tmp(k+1),GV%Angstrom_H)
+          rho_h(k) = h_tmp(k) * rho_col(k)
+          mass_disp = mass_disp + rho_h(k)
+        enddo
+        residual = mass_shelf(i,j) - mass_disp
+        iter = iter+1
+      end do
+      if (iter .ge. max_iter) call MOM_mesg("Warning: calc_sfc_displacement too many iterations.")
+      z_top_shelf(i,j) = z_top
+    endif
+  enddo; enddo
+  call MOM_mesg("Calling depress_surface ")
+  call depress_surface(h, G, GV, US, PF, tv, just_read=.false.,z_top_shelf=z_top_shelf)
+  call MOM_mesg("Finishing calling depress_surface ")
+end subroutine calc_sfc_displacement
 
 !> Adjust the layer thicknesses by removing the top of the water column above the
 !! depth where the hydrostatic pressure matches p_surf
@@ -2597,6 +2703,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just
                           old_remap=remap_old_alg, answers_2018=answers_2018 )
     call ALE_remap_scalar(remapCS, G, GV, nkd, h1, tmpS1dIn, h, tv%S, all_cells=remap_full_column, &
                           old_remap=remap_old_alg, answers_2018=answers_2018 )
+
     deallocate( h1 )
     deallocate( tmpT1dIn )
     deallocate( tmpS1dIn )