diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index b16681156b..7a15d58bb1 100644 --- a/src/core/MOM.F90 +++ b/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)") diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 9dd3791211..13af5a936a 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/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 ! 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 diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 0d5342d9be..f95192f5f8 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/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) 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 )