Skip to content

Static ice shelf initialize #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
MJHarrison-GFDL wants to merge 2 commits into from
Closed

Static ice shelf initialize #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d168ede
Hydrostatic initialization in ice cavities
MJHarrison-GFDL Dec 13, 2021
47e13d7
max number of iterations for ice shelf initialization
MJHarrison-GFDL Dec 13, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 17 additions & 8 deletions src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -277,7 +277,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 @@ -1703,7 +1704,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &

! Initial state on the input index map
real, allocatable, dimension(:,:,:) :: u_in, v_in, h_in
real, allocatable, dimension(:,:), target :: frac_shelf_in
real, allocatable, dimension(:,:), target :: frac_shelf_in, mass_shelf_in
real, allocatable, dimension(:,:,:), target :: T_in, S_in
type(ocean_OBC_type), pointer :: OBC_in => NULL()
type(sponge_CS), pointer :: sponge_in_CSp => NULL()
Expand Down Expand Up @@ -2457,14 +2458,18 @@ 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)
CS%mass_shelf = CS%mass_shelf*US%kg_m3_to_R*US%m_to_Z
! 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 @@ -2507,25 +2512,29 @@ 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)
CS%mass_shelf = CS%mass_shelf*US%kg_m3_to_R*US%m_to_Z
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, &
CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp, CS%OBC, Time_in)
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)
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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2034,11 +2034,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)) :: frac_shelf_h !< Ice shelf area fraction [nodim].
real, optional, dimension(SZI_(G),SZJ_(G)) :: mass_shelf !<Ice shelf mass [RZ -> kg m-2]


integer :: i, j

Expand All @@ -2049,6 +2050,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: 133 additions & 30 deletions src/initialization/MOM_state_initialization.F90
View file
Open in desktop
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 @@ -148,6 +148,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 [ RZ ~> 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 @@ -159,6 +162,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 @@ -257,7 +261,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
"use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true")

call MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, &
just_read_params=just_read, frac_shelf_h=frac_shelf_h)
just_read_params=just_read, frac_shelf_h=frac_shelf_h)
else
! Initialize thickness, h.
call get_param(PF, mdl, "THICKNESS_CONFIG", config, &
Expand Down Expand Up @@ -405,6 +409,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

if (new_sim) then
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 (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 @@ -459,15 +480,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 @@ -1047,7 +1059,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_params)
subroutine depress_surface(h, G, GV, US, param_file, tv, just_read_params, 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 @@ -1056,7 +1068,9 @@ subroutine depress_surface(h, G, GV, US, param_file, tv, just_read_params)
type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
logical, optional, intent(in) :: just_read_params !< If present and true, this call will
!! only read parameters without changing h.
!! 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 @@ -1070,31 +1084,41 @@ subroutine depress_surface(h, G, GV, US, param_file, tv, just_read_params)
character(len=200) :: filename, eta_srf_var ! Strings for file/path
logical :: just_read ! If true, just read parameters but set nothing.
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

just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

! 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.

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.
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 @@ -1218,6 +1242,84 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read_params)

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 [RZ ~> 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 [RZ ~> kg m-2 ]
real, dimension(SZK_(GV)) :: h_tmp ! temporary storage for thicknesses [H ~> m]
real, dimension(SZK_(GV)) :: p_ref ! pressure for density [RZ ~> 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 = US%m_to_Z*0.001 ! 1mm tolerance
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) > 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)>tol .and. z_top>-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, 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 @@ -2640,6 +2742,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