Skip to content

Commit

Permalink
Hydrostatic initialization in ice cavities (NOAA-EMC#41)
Browse files Browse the repository at this point in the history 
* 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>
  • Loading branch information
@MJHarrison-GFDL @marshallward
MJHarrison-GFDL and marshallward authored Jan 6, 2022
1 parent 2d32631 commit df46be4
Show file tree
Hide file tree
Showing 3 changed files with 164 additions and 39 deletions.
26 changes: 18 additions & 8 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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].
Expand Down Expand Up @@ -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()
Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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, &
Expand All @@ -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)")
Expand Down
14 changes: 11 additions & 3 deletions src/ice_shelf/MOM_ice_shelf.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand Down
163 changes: 135 additions & 28 deletions src/initialization/MOM_state_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)) &
Expand Down Expand Up @@ -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
Expand All @@ -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].
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 )
Expand Down

0 comments on commit df46be4

Please sign in to comment.