MOM_ice_shelf.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Implements the thermodynamic aspects of ocean / ice-shelf interactions,
6!! along with a crude placeholder for a later implementation of full
7!! ice shelf dynamics, all using the MOM framework and coding style.
8module mom_ice_shelf
9
10use mom_array_transform, only : rotate_array
11use mom_constants, only : hlf
12use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
13use mom_cpu_clock, only : clock_component, clock_routine
14use mom_coms, only : num_pes, reproducing_sum
15use mom_data_override, only : data_override
16use mom_diag_mediator, only : mom_diag_ctrl=>diag_ctrl
17use mom_is_diag_mediator, only : post_data=>post_is_data, post_scalar_data=>post_is_data_0d
18use mom_is_diag_mediator, only : register_diag_field=>register_mom_is_diag_field, safe_alloc_ptr
19use mom_is_diag_mediator, only : register_scalar_field=>register_mom_is_scalar_field
26use mom_domains, only : mom_domains_init, pass_var, pass_vector, clone_mom_domain
27use mom_domains, only : to_all, cgrid_ne, bgrid_ne, corner
29use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
32use mom_file_parser, only : read_param, get_param, log_param, log_version, param_file_type
38use mom_fixed_initialization, only : mom_initialize_rotation
40use mom_io, only : field_exists, file_exists, mom_read_data, write_version_number
41use mom_io, only : slasher, fieldtype, vardesc, var_desc
42use mom_io, only : close_file, single_file, multiple
43use mom_restart, only : register_restart_field, save_restart
44use mom_restart, only : restart_init, restore_state, mom_restart_cs, register_restart_pair
45use mom_time_manager, only : time_type, time_to_real, real_to_time, operator(>), operator(-)
47use mom_transcribe_grid, only : rotate_dyngrid
51use mom_forcing_type, only : forcing, allocate_forcing_type, deallocate_forcing_type, mom_forcing_chksum
55use mom_eos, only : calculate_density, calculate_density_derivs, calculate_tfreeze, eos_domain
56use mom_eos, only : eos_type, eos_init
63!MJH use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary
68use mom_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum
69use mom_interpolate, only : init_external_field, time_interp_external, time_interp_external_init
70use mom_interpolate, only : external_field
71
72implicit none ; private
73
74#include <MOM_memory.h>
75#ifdef SYMMETRIC_MEMORY_
76# define GRID_SYM_ .true.
77#else
78# define GRID_SYM_ .false.
79#endif
80
87public update_ice_smb
88
89! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
90! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
91! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
92! vary with the Boussinesq approximation, the Boussinesq variant is given first.
93
94!> Control structure that contains ice shelf parameters and diagnostics handles
95type, public :: ice_shelf_cs ; private
96 ! Parameters
97 type(mom_restart_cs), pointer :: restart_csp => null() !< A pointer to the restart control
98 !! structure for the ice shelves
99 type(ocean_grid_type), pointer :: grid_in => null() !< un-rotated input grid metric
100 logical :: rotate_index = .false. !< True if index map is rotated
101 integer :: turns !< The number of quarter turns for rotation testing.
102 type(ocean_grid_type), pointer :: grid => null() !< Grid for the ice-shelf model
103 type(unit_scale_type), pointer :: &
104 us => null() !< A structure containing various unit conversion factors
105 type(ocean_grid_type), pointer :: ocn_grid => null() !< A pointer to the ocean model grid
106 !! The rest is private
107 real :: flux_factor = 1.0 !< A factor that can be used to turn off ice shelf
108 !! melting (flux_factor = 0) [nondim].
109 character(len=128) :: restart_output_dir = ' ' !< The directory in which to write restart files
110 type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
111 !! the ice-shelf state
112 type(ice_shelf_dyn_cs), pointer :: dcs => null() !< The control structure for the ice-shelf dynamics.
113
114 real, pointer, dimension(:,:) :: &
115 utide => null() !< An unresolved tidal velocity [L T-1 ~> m s-1]
116
117 real :: ustar_bg !< A minimum value for ustar under ice shelves [Z T-1 ~> m s-1].
118 real :: ustar_max !< A maximum value for ustar under ice shelves, or a negative value to
119 !! have no limit [Z T-1 ~> m s-1].
120 real :: cdrag !< drag coefficient under ice shelves [nondim].
121 real :: g_earth !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
122 real :: cp !< The heat capacity of sea water [Q C-1 ~> J kg-1 degC-1].
123 real :: rho_ocn !< A reference ocean density [R ~> kg m-3].
124 real :: cp_ice !< The heat capacity of fresh ice [Q C-1 ~> J kg-1 degC-1].
125 real :: gamma_t !< The (fixed) turbulent exchange velocity in the
126 !< 2-equation formulation [Z T-1 ~> m s-1].
127 real :: salin_ice !< The salinity of shelf ice [S ~> ppt].
128 real :: temp_ice !< The core temperature of shelf ice [C ~> degC].
129 real :: kv_ice !< The viscosity of ice [L4 Z-2 T-1 ~> m2 s-1].
130 real :: density_ice !< A typical density of ice [R ~> kg m-3].
131 real :: kv_molec !< The molecular kinematic viscosity of sea water [Z2 T-1 ~> m2 s-1].
132 real :: kd_molec_salt!< The molecular diffusivity of salt [Z2 T-1 ~> m2 s-1].
133 real :: kd_molec_temp!< The molecular diffusivity of heat [Z2 T-1 ~> m2 s-1].
134 real :: lat_fusion !< The latent heat of fusion [Q ~> J kg-1].
135 real :: gamma_t_3eq !< Nondimensional heat-transfer coefficient, used in the 3Eq. formulation [nondim]
136 real :: gamma_s_3eq !< Nondimensional salt-transfer coefficient, used in the 3Eq. formulation [nondim]
137 !< This number should be specified by the user.
138 real :: col_mass_melt_threshold !< An ocean column mass below the iceshelf below which melting
139 !! does not occur [R Z ~> kg m-2]
140 logical :: mass_from_file !< Read the ice shelf mass from a file every dt
141 logical :: ustar_shelf_from_vel !< If true, use the surface velocities, and not the previous
142 !! values of the stresses to set ustar.
143
144 !!!! PHYSICAL AND NUMERICAL PARAMETERS FOR ICE DYNAMICS !!!!!!
145
146 real :: time_step !< this is the shortest timestep that the ice shelf sees [T ~> s], and
147 !! is equal to the forcing timestep (it is passed in when the shelf
148 !! is initialized - so need to reorganize MOM driver.
149 !! it will be the prognostic timestep ... maybe.
150
151 logical :: solo_ice_sheet !< whether the ice model is running without being
152 !! coupled to the ocean
153 logical :: gl_regularize !< whether to regularize the floatation condition
154 !! at the grounding line a la Goldberg Holland Schoof 2009
155 logical :: gl_couple !< whether to let the floatation condition be
156 !!determined by ocean column thickness means update_OD_ffrac
157 !! will be called (note: GL_regularize and GL_couple
158 !! should be exclusive)
159 logical :: calve_to_mask !< If true, calve any ice that passes outside of a masked area
160 logical :: calve_ice_shelf_bergs=.false. !< If true, flux through a static ice front is converted
161 !! to point bergs
162 real :: min_thickness_simple_calve !< min. ice shelf thickness criteria for calving [Z ~> m].
163 real :: t0 !< temperature at ocean surface in the restoring region [C ~> degC]
164 real :: s0 !< Salinity at ocean surface in the restoring region [S ~> ppt].
165 real :: input_flux !< The vertically integrated inward ice thickness flux per
166 !! unit face length at an upstream boundary [Z L T-1 ~> m2 s-1]
167 real :: input_thickness !< Ice thickness at an upstream open boundary [Z ~> m].
168
169 type(time_type) :: time !< The component's time.
170 type(eos_type) :: eqn_of_state !< Type that indicates the equation of state to use.
171 logical :: active_shelf_dynamics !< True if the ice shelf mass changes as a result
172 !! the dynamic ice-shelf model.
173 logical :: shelf_mass_is_dynamic !< True if ice shelf mass changes over time. If true, ice
174 !! shelf dynamics will be initialized
175 logical :: data_override_shelf_fluxes !< True if the ice shelf surface mass fluxes can be
176 !! written using the data_override feature (only for MOSAIC grids)
177 logical :: override_shelf_movement !< If true, user code specifies the shelf movement
178 !! instead of using the dynamic ice-shelf mode.
179 logical :: isthermo !< True if the ice shelf can exchange heat and
180 !! mass with the underlying ocean.
181 logical :: threeeq !< If true, the 3 equation consistency equations are
182 !! used to calculate the flux at the ocean-ice
183 !! interface.
184 logical :: insulator !< If true, ice shelf is a perfect insulator
185 logical :: const_gamma !< If true, gamma_T is specified by the user.
186 logical :: constant_sea_level !< if true, apply an evaporative, heat and salt
187 !! fluxes. It will avoid large increase in sea level.
188 logical :: constant_sea_level_misomip !< If true, constant_sea_level fluxes are applied only over
189 !! the surface sponge cells from the ISOMIP/MISOMIP configuration
190 logical :: smb_diag !< If true, calculate diagnostics related to surface mass balance
191 logical :: bmb_diag !< If true, calculate diagnostics related to basal mass balance
192 real :: min_ocean_mass_float !< The minimum ocean mass per unit area before the ice
193 !! shelf is considered to float when constant_sea_level
194 !! is used [R Z ~> kg m-2]
195 real :: cutoff_depth !< Depth above which melt is set to zero (>= 0) [Z ~> m].
196 logical :: find_salt_root !< If true, if true find Sbdry using a quadratic eq.
197 real :: tfr_0_0 !< The freezing point at 0 pressure and 0 salinity [C ~> degC]
198 real :: dtfr_ds !< Partial derivative of freezing temperature with
199 !! salinity [C S-1 ~> degC ppt-1]
200 real :: dtfr_dp !< Partial derivative of freezing temperature with
201 !! pressure [C T2 R-1 L-2 ~> degC Pa-1]
202 real :: zeta_n !< The stability constant xi_N = 0.052 from Holland & Jenkins '99
203 !! divided by the von Karman constant VK [nondim]. Was 1/8.
204 real :: vk !< Von Karman's constant [nondim]
205 real :: rc !< critical flux Richardson number [nondim]
206 logical :: ustar_from_vel_bugfix !< If true, fixes ustar from ocean velocity bug
207 logical :: buoy_flux_itt_bugfix !< If true, fixes buoyancy iteration bug
208 logical :: salt_flux_itt_bugfix !< If true, fixes salt iteration bug
209 real :: buoy_flux_tol !< Fractional buoyancy iteration tolerance for convergence [nondim]
210
211 !>@{ Diagnostic handles
212 integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, &
213 id_tfreeze = -1, id_tfl_shelf = -1, &
214 id_thermal_driving = -1, id_haline_driving = -1, &
215 id_u_ml = -1, id_v_ml = -1, id_sbdry = -1, &
216 id_h_shelf = -1, id_dhdt_shelf = -1, id_h_mask = -1, id_frazil = -1, &
217 id_surf_elev = -1, id_bathym = -1, &
218 id_area_shelf_h = -1, &
219 id_ustar_shelf = -1, id_shelf_mass = -1, id_mass_flux = -1, &
220 id_shelf_sfc_mass_flux = -1, &
221 id_vaf = -1, id_g_adott = -1, id_f_adott = -1, id_adott = -1, &
222 id_bdott_melt = -1, id_bdott_accum = -1, id_bdott = -1, &
223 id_dvafdt = -1, id_g_adot = -1, id_f_adot = -1, id_adot = -1, &
224 id_bdot_melt = -1, id_bdot_accum = -1, id_bdot = -1, &
225 id_t_area = -1, id_g_area = -1, id_f_area = -1, &
226 id_ant_vaf = -1, id_ant_g_adott = -1, id_ant_f_adott = -1, id_ant_adott = -1, &
227 id_ant_bdott_melt = -1, id_ant_bdott_accum = -1, id_ant_bdott = -1, &
228 id_ant_dvafdt = -1, id_ant_g_adot = -1, id_ant_f_adot = -1, id_ant_adot = -1, &
229 id_ant_bdot_melt = -1, id_ant_bdot_accum = -1, id_ant_bdot = -1, &
230 id_ant_t_area = -1, id_ant_g_area = -1, id_ant_f_area = -1, &
231 id_gr_vaf = -1, id_gr_g_adott = -1, id_gr_f_adott = -1, id_gr_adott = -1, &
232 id_gr_bdott_melt = -1, id_gr_bdott_accum = -1, id_gr_bdott = -1, &
233 id_gr_dvafdt = -1, id_gr_g_adot = -1, id_gr_f_adot = -1, id_gr_adot = -1, &
234 id_gr_bdot_melt = -1, id_gr_bdot_accum = -1, id_gr_bdot = -1, &
235 id_gr_t_area = -1, id_gr_g_area = -1, id_gr_f_area = -1
236 !>@}
237
238 type(external_field) :: smb_file
239 !< Handle for reading the time interpolated smb from a file
240 logical :: time_varying_smb
241 !< logical flag set true if reading time-varying smb
242 type(external_field) :: mass_handle
243 !< Handle for reading the time interpolated ice shelf mass from a file
244 type(external_field) :: area_handle
245 !< Handle for reading the time interpolated ice shelf area from a file
246
247 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to control diagnostic output.
248 type(user_ice_shelf_cs), pointer :: user_cs => null() !< A pointer to the control structure for
249 !! user-supplied modifications to the ice shelf code.
250
251 logical :: debug !< If true, write verbose checksums for debugging purposes
252 !! and use reproducible sums
253end type ice_shelf_cs
254
255!>@{ CPU time clock IDs
256integer :: id_clock_shelf=-1 !< CPU Clock for the ice shelf code
257integer :: id_clock_pass=-1 !< CPU Clock for ice shelf group pass calls
258!>@}
259
260contains
261
262!> Calculates fluxes between the ocean and ice-shelf using the three-equations
263!! formulation (optional to use just two equations).
264!! See \ref section_ICE_SHELF_equations
265subroutine shelf_calc_flux(sfc_state_in, fluxes_in, Time, time_step_in, CS)
266 type(surface), target, intent(inout) :: sfc_state_in !< A structure containing fields that
267 !! describe the surface state of the ocean. The
268 !! intent is only inout to allow for halo updates.
269 type(forcing), target, intent(inout) :: fluxes_in !< structure containing pointers to any
270 !! possible thermodynamic or mass-flux forcing fields.
271 type(time_type), intent(in) :: time !< Start time of the fluxes.
272 real, intent(in) :: time_step_in !< Length of time over which these fluxes
273 !! will be applied [T ~> s].
274 type(ice_shelf_cs), pointer :: cs !< A pointer to the control structure returned
275 !! by a previous call to initialize_ice_shelf.
276
277 ! Local variables
278 type(ocean_grid_type), pointer :: g => null() !< The grid structure used by the ice shelf.
279 type(unit_scale_type), pointer :: us => null() !< Pointer to a structure containing
280 !! various unit conversion factors
281 type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
282 !! the ice-shelf state
283
284 type(surface), pointer :: sfc_state => null()
285 type(forcing), pointer :: fluxes => null()
286
287 real, dimension(SZI_(CS%grid)) :: &
288 rhoml, & !< Ocean mixed layer density [R ~> kg m-3].
289 dr0_dt, & !< Partial derivative of the mixed layer density
290 !< with temperature [R C-1 ~> kg m-3 degC-1].
291 dr0_ds, & !< Partial derivative of the mixed layer density
292 !< with salinity [R S-1 ~> kg m-3 ppt-1].
293 p_int !< The pressure at the ice-ocean interface [R L2 T-2 ~> Pa].
294
295 real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: &
296 exch_vel_t, & !< Sub-shelf thermal exchange velocity [Z T-1 ~> m s-1]
297 exch_vel_s, & !< Sub-shelf salt exchange velocity [Z T-1 ~> m s-1]
298 dh_bdott, & !< Basal melt/accumulation over a time step, used for diagnostics [Z ~> m]
299 dh_adott !< Surface melt/accumulation over a time step, used for diagnostics [Z ~> m]
300 real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
301 mass_flux !< Total mass flux of freshwater across the ice-ocean interface. [R Z L2 T-1 ~> kg s-1]
302 real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
303 haline_driving !< (SSS - S_boundary) ice-ocean
304 !! interface, positive for melting and negative for freezing [S ~> ppt].
305 !! This is computed as part of the ISOMIP diagnostics.
306 real :: time_step !< Length of time over which these fluxes will be applied [T ~> s].
307 real :: itime_step !< Inverse of the length of time over which these fluxes will be applied [T-1 ~> s-1]
308 real :: vk !< Von Karman's constant [nondim]
309 real :: zeta_n !< This is the stability constant xi_N = 0.052 from Holland & Jenkins '99
310 !! divided by the von Karman constant VK. Was 1/8. [nondim]
311 real :: rf_crit !< critical flux Richardson number [nondim]
312 real :: i_2zeta_n !< Half the inverse of Zeta_N [nondim].
313 real :: i_lf !< The inverse of the latent heat of fusion [Q-1 ~> kg J-1].
314 real :: i_dt_lhf ! The inverse of the timestep times the latent heat of fusion [Q-1 T-1 ~> kg J-1 s-1].
315 real :: i_vk !< The inverse of the Von Karman constant [nondim].
316 real :: pr, sc !< The Prandtl number and Schmidt number [nondim].
317
318 ! 3 equations formulation variables
319 real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
320 sbdry !< Salinities in the ocean at the interface with the ice shelf [S ~> ppt].
321 real :: sbdry_it ! The boundary salinity at an iteration [S ~> ppt]
322 real :: s_a ! A variable used to find salt roots [S-1 ~> ppt-1]
323 real :: s_b ! A variable used to find salt roots [nondim]
324 real :: s_c ! A variable used to find salt roots [S ~> ppt]
325 real :: ds_it !< The interface salinity change during an iteration [S ~> ppt].
326 real :: hbl_neut !< The neutral boundary layer thickness [Z ~> m].
327 real :: hbl_neut_h_molec !< The ratio of the neutral boundary layer thickness
328 !! to the molecular boundary layer thickness [nondim].
329 real :: wt_flux !< The downward vertical flux of heat just inside the ocean [C Z T-1 ~> degC m s-1].
330 real :: wb_flux !< The downward vertical flux of buoyancy just inside the ocean [Z2 T-3 ~> m2 s-3].
331 real :: db_ds !< The derivative of buoyancy with salinity [Z T-2 S-1 ~> m s-2 ppt-1].
332 real :: db_dt !< The derivative of buoyancy with temperature [Z T-2 C-1 ~> m s-2 degC-1].
333 real :: i_n_star ! The inverse of the ratio of working boundary layer thickness
334 ! to the neutral thickness [nondim]
335 real :: n_star_term ! A term in the expression for nstar [T3 Z-2 ~> s3 m-2]
336 real :: absf ! The absolute value of the Coriolis parameter [T-1 ~> s-1]
337 real :: dins_dwb !< The partial derivative of I_n_star with wB_flux, in [T3 Z-2 ~> s3 m-2]
338 real :: dt_ustar ! The difference between the freezing point and the ocean boundary layer
339 ! temperature times the friction velocity [C Z T-1 ~> degC m s-1]
340 real :: ds_ustar ! The difference between the salinity at the ice-ocean interface and the ocean
341 ! boundary layer salinity times the friction velocity [S Z T-1 ~> ppt m s-1]
342 real :: ustar_h ! The friction velocity in the water below the ice shelf [Z T-1 ~> m s-1]
343 real :: gam_turb ! A relative turbluent diffusivity [nondim]
344 real :: gam_mol_t, gam_mol_s ! Relative coefficients of molecular diffusivities [nondim]
345 real :: rhocp ! A typical ocean density times the heat capacity of water [Q R C-1 ~> J m-3 degC-1]
346 real :: ln_neut ! The log of the ratio of the neutral boundary layer thickness to the molecular
347 ! boundary layer thickness if it is greater than 1 or 0 otherwise [nondim]
348 real :: mass_exch ! A mass exchange rate [R Z T-1 ~> kg m-2 s-1]
349 real :: sb_min, sb_max ! Minimum and maximum boundary salinities [S ~> ppt]
350 real :: ds_min, ds_max ! Minimum and maximum salinity changes [S ~> ppt]
351 ! Variables used in iterating for wB_flux.
352 real :: wb_flux_next ! The next interation's guess for wB_flux [Z2 T-3 ~> m2 s-3]
353 real :: wb_flux_new ! An updated value of wB_flux when Gam_turb is based on wB_flux [Z2 T-3 ~> m2 s-3]
354 real :: wb_flux_max ! The upper bound on wB_flux [Z2 T-3 ~> m2 s-3]
355 real :: wb_flux_min ! The lower bound on wB_flux [Z2 T-3 ~> m2 s-3]
356 real :: ddwb_dwb ! The slope of the change in wB_flux between iterations with wB_flux [nondim]
357 real :: dwb_max ! The change in wB_flux when it is wB_flux_max [Z2 T-3 ~> m2 s-3]
358 real :: dwb_min ! The change in wB_flux when it is wB_flux_min [Z2 T-3 ~> m2 s-3]
359 real :: i_gam_t, i_gam_s ! Terms that vary inversely with Gam_mol_T or Gam_mol_S and Gam_turb [nondim]
360 real :: dg_dwb ! The derivative of Gam_turb with wB [T3 Z-2 ~> s3 m-2]
361 real :: taux2, tauy2 ! The squared surface stresses [R2 L2 Z2 T-4 ~> Pa2].
362 real :: u2_av, v2_av ! The ice-area weighted average squared ocean velocities [L2 T-2 ~> m2 s-2]
363 real :: asu1, asu2 ! Ocean areas covered by ice shelves at neighboring u-points [L2 ~> m2]
364 real :: asv1, asv2 ! Ocean areas covered by ice shelves at neighboring v-points [L2 ~> m2]
365 real :: i_au, i_av ! The Adcroft reciprocals of the ice shelf areas at adjacent points [L-2 ~> m-2]
366 real :: irho0 ! The inverse of the mean density times a unit conversion factor [R-1 L Z-1 ~> m3 kg-1]
367 logical :: sb_min_set, sb_max_set
368 logical :: root_found
369 logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
370 logical :: coupled_gl ! If true, the grounding line position is determined based on
371 ! coupled ice-ocean dynamics.
372 logical :: add_frazil ! If true, allow frazil formation to modify ice-shelf water flux
373 real, parameter :: c2_3 = 2.0/3.0 ! Two thirds [nondim]
374 character(len=320) :: mesg ! The text of an error message
375 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
376 integer :: i, j, is, ie, js, je, ied, jed, it1, it3
377 real :: vaf0, vaf0_a, vaf0_g ! The previous volumes above floatation [Z L2 ~> m3]
378 ! for all ice sheets, Antarctica only, or Greenland only
379
380 if (.not. associated(cs)) call mom_error(fatal, "shelf_calc_flux: "// &
381 "initialize_ice_shelf must be called before shelf_calc_flux.")
382 call cpu_clock_begin(id_clock_shelf)
383
384 g => cs%grid ; us => cs%US
385 iss => cs%ISS
386 time_step = time_step_in
387 itime_step = 1./time_step
388
389 dh_adott(:,:) = 0.0 ; dh_bdott(:,:) = 0.0
390
391 if (cs%active_shelf_dynamics) then
392 !calculate previous volumes above floatation
393 if (cs%id_dvafdt > 0) call volume_above_floatation(cs%dCS, g, iss, vaf0) !all ice sheet
394 if (cs%id_Ant_dvafdt > 0) call volume_above_floatation(cs%dCS, g, iss, vaf0_a, hemisphere=0) !Antarctica only
395 if (cs%id_Gr_dvafdt > 0) call volume_above_floatation(cs%dCS, g, iss, vaf0_g, hemisphere=1) !Greenland only
396 endif
397
398 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; ied = g%ied ; jed = g%jed
399 if (cs%data_override_shelf_fluxes .and. cs%active_shelf_dynamics) then
400 call data_override(g%Domain, 'shelf_sfc_mass_flux', fluxes_in%shelf_sfc_mass_flux(is:ie,js:je), cs%Time, &
401 scale=us%kg_m2s_to_RZ_T)
402 call pass_var(fluxes_in%shelf_sfc_mass_flux, g%domain, complete=.true.)
403 endif
404
405 if (cs%rotate_index) then
406 allocate(sfc_state)
407 call rotate_surface_state(sfc_state_in, sfc_state, cs%Grid, cs%turns)
408 allocate(fluxes)
409 call allocate_forcing_type(fluxes_in, g, fluxes, turns=cs%turns)
410 call rotate_forcing(fluxes_in, fluxes, cs%turns)
411 else
412 sfc_state => sfc_state_in
413 fluxes => fluxes_in
414 endif
415 ! useful parameters
416 zeta_n = cs%Zeta_N
417 vk = cs%Vk
418 rf_crit = cs%Rc
419 i_2zeta_n = 0.5 / cs%Zeta_N
420 i_lf = 1.0 / cs%Lat_fusion
421 i_dt_lhf = 1.0 / (time_step * cs%Lat_fusion)
422 sc = cs%kv_molec/cs%kd_molec_salt
423 pr = cs%kv_molec/cs%kd_molec_temp
424 i_vk = 1.0/vk
425 rhocp = cs%Rho_ocn * cs%Cp
426
427 !first calculate molecular component
428 gam_mol_t = 12.5 * (pr**c2_3) - 6.0
429 gam_mol_s = 12.5 * (sc**c2_3) - 6.0
430
431 ! GMM, zero some fields of the ice shelf structure (ice_shelf_CS)
432 ! these fields are already set to zero during initialization
433 ! However, they seem to be changed somewhere and, for diagnostic
434 ! reasons, it is better to set them to zero again.
435 exch_vel_t(:,:) = 0.0 ; exch_vel_s(:,:) = 0.0
436 iss%tflux_shelf(:,:) = 0.0 ; iss%water_flux(:,:) = 0.0
437 iss%salt_flux(:,:) = 0.0 ; iss%tflux_ocn(:,:) = 0.0 ; iss%tfreeze(:,:) = 0.0
438 ! define Sbdry to avoid Run-Time Check Failure, when melt is not computed.
439 haline_driving(:,:) = 0.0
440 sbdry(:,:) = sfc_state%sss(:,:)
441
442 !update time
443 cs%Time = time
444
445 if (cs%override_shelf_movement) then
446 cs%time_step = time_step
447 ! update shelf mass
448 if (cs%mass_from_file) call update_shelf_mass(g, us, cs, iss, time)
449 endif
450
451 if (cs%debug) then
452 call hchksum(fluxes_in%frac_shelf_h, "frac_shelf_h before apply melting", cs%Grid_in%HI, haloshift=0)
453 call hchksum(sfc_state_in%sst, "sst before apply melting", cs%Grid_in%HI, haloshift=0, unscale=us%C_to_degC)
454 call hchksum(sfc_state_in%sss, "sss before apply melting", cs%Grid_in%HI, haloshift=0, unscale=us%S_to_ppt)
455 call uvchksum("[uv]_ml before apply melting", sfc_state_in%u, sfc_state_in%v, &
456 cs%Grid_in%HI, haloshift=0, unscale=us%L_T_to_m_s)
457 call hchksum(sfc_state_in%ocean_mass, "ocean_mass before apply melting", cs%Grid_in%HI, haloshift=0, &
458 unscale=us%RZ_to_kg_m2)
459 endif
460
461 ! Calculate the friction velocity under ice shelves, using taux_shelf and tauy_shelf if possible.
462 if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then
463 call pass_vector(sfc_state%taux_shelf, sfc_state%tauy_shelf, g%domain, to_all, cgrid_ne)
464 endif
465 irho0 = us%Z_to_L / cs%Rho_ocn
466 do j=js,je ; do i=is,ie ; if (fluxes%frac_shelf_h(i,j) > 0.0) then
467 taux2 = 0.0 ; tauy2 = 0.0 ; u2_av = 0.0 ; v2_av = 0.0
468 asu1 = (iss%area_shelf_h(i-1,j) + iss%area_shelf_h(i,j))
469 asu2 = (iss%area_shelf_h(i,j) + iss%area_shelf_h(i+1,j))
470 asv1 = (iss%area_shelf_h(i,j-1) + iss%area_shelf_h(i,j))
471 asv2 = (iss%area_shelf_h(i,j) + iss%area_shelf_h(i,j+1))
472 i_au = 0.0 ; if (asu1 + asu2 > 0.0) i_au = 1.0 / (asu1 + asu2)
473 i_av = 0.0 ; if (asv1 + asv2 > 0.0) i_av = 1.0 / (asv1 + asv2)
474 if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then
475 taux2 = (((asu1 * (sfc_state%taux_shelf(i-1,j)**2)) + (asu2 * (sfc_state%taux_shelf(i,j)**2)) ) * i_au)
476 tauy2 = (((asv1 * (sfc_state%tauy_shelf(i,j-1)**2)) + (asv2 * (sfc_state%tauy_shelf(i,j)**2)) ) * i_av)
477 endif
478 u2_av = (((asu1 * (sfc_state%u(i-1,j)**2)) + (asu2 * sfc_state%u(i,j)**2)) * i_au)
479 if (cs%ustar_from_vel_bugfix) then
480 v2_av = (((asv1 * (sfc_state%v(i,j-1)**2)) + (asv2 * sfc_state%v(i,j)**2)) * i_av)
481 else
482 v2_av = (((asv1 * (sfc_state%v(i,j-1)**2)) + (asu2 * sfc_state%v(i,j)**2)) * i_av)
483 endif
484
485 if ((taux2 + tauy2 > 0.0) .and. .not.cs%ustar_shelf_from_vel) then
486 if (cs%ustar_max >= 0.0) then
487 fluxes%ustar_shelf(i,j) = min(cs%ustar_max, max(cs%ustar_bg, us%L_to_Z * &
488 sqrt(irho0 * sqrt(taux2 + tauy2) + cs%cdrag*cs%utide(i,j)**2)))
489 else
490 fluxes%ustar_shelf(i,j) = max(cs%ustar_bg, us%L_to_Z * &
491 sqrt(irho0 * sqrt(taux2 + tauy2) + cs%cdrag*cs%utide(i,j)**2))
492 endif
493 else ! Take care of the cases when taux_shelf is not set or not allocated.
494 fluxes%ustar_shelf(i,j) = max(cs%ustar_bg, us%L_TO_Z * &
495 sqrt(cs%cdrag*((u2_av + v2_av) + cs%utide(i,j)**2)))
496 endif
497 else ! There is no shelf here.
498 fluxes%ustar_shelf(i,j) = 0.0
499 endif ; enddo ; enddo
500
501 eosdom(:) = eos_domain(g%HI)
502 do j=js,je
503 ! Find the pressure at the ice-ocean interface, averaged only over the
504 ! part of the cell covered by ice shelf.
505 do i=is,ie ; p_int(i) = cs%g_Earth * iss%mass_shelf(i,j) ; enddo
506
507 ! Calculate insitu densities and expansion coefficients
508 call calculate_density(sfc_state%sst(:,j), sfc_state%sss(:,j), p_int, rhoml(:), &
509 cs%eqn_of_state, eosdom)
510 call calculate_density_derivs(sfc_state%sst(:,j), sfc_state%sss(:,j), p_int, &
511 dr0_dt, dr0_ds, cs%eqn_of_state, eosdom)
512
513 do i=is,ie
514 if ((sfc_state%ocean_mass(i,j) > cs%col_mass_melt_threshold) .and. &
515 (iss%area_shelf_h(i,j) > 0.0) .and. cs%isthermo &
516 .and. iss%melt_mask(i,j)>0.0) then
517
518 if (cs%threeeq) then
519 ! Iteratively determine a self-consistent set of fluxes, with the ocean
520 ! salinity just below the ice-shelf as the variable that is being
521 ! iterated for.
522
523 ustar_h = fluxes%ustar_shelf(i,j)
524
525 ! Estimate the neutral ocean boundary layer thickness as the minimum of the
526 ! reported ocean mixed layer thickness and the neutral Ekman depth.
527 absf = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
528 (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
529 if (absf*sfc_state%Hml(i,j) <= vk*ustar_h) then ; hbl_neut = sfc_state%Hml(i,j)
530 else ; hbl_neut = (vk*ustar_h) / absf ; endif
531 hbl_neut_h_molec = zeta_n * ((hbl_neut * ustar_h) / (5.0 * cs%kv_molec))
532 ln_neut = 0.0 ; if (hbl_neut_h_molec > 1.0) ln_neut = log(hbl_neut_h_molec)
533 n_star_term = (zeta_n * hbl_neut * vk) / (rf_crit * ustar_h**3)
534
535 ! Determine the mixed layer buoyancy flux, wB_flux.
536 db_ds = (us%L_to_Z**2*cs%g_Earth / rhoml(i)) * dr0_ds(i)
537 db_dt = (us%L_to_Z**2*cs%g_Earth / rhoml(i)) * dr0_dt(i)
538
539 if (cs%find_salt_root) then
540 ! Solve for the skin salinity using the linearized liquidus parameters and
541 ! balancing the turbulent fresh water flux in the near-boundary layer with
542 ! the net fresh water or salt added by melting:
543 ! (Cp/Lat_fusion)*Gamma_T_3Eq*(TFr_skin-T_ocn) = Gamma_S_3Eq*(S_skin-S_ocn)/S_skin
544
545 ! S_a is always < 0.0 with a realistic expression for the freezing point.
546 s_a = cs%dTFr_dS * cs%Gamma_T_3EQ * cs%Cp
547 s_b = cs%Gamma_T_3EQ*cs%Cp*(cs%TFr_0_0 + cs%dTFr_dp*p_int(i) - sfc_state%sst(i,j)) - &
548 cs%Lat_fusion * cs%Gamma_S_3EQ ! S_b Can take either sign, but is usually negative.
549 s_c = cs%Lat_fusion * cs%Gamma_S_3EQ * sfc_state%sss(i,j) ! Always >= 0
550
551 if (s_c == 0.0) then ! The solution for fresh water.
552 sbdry(i,j) = 0.0
553 elseif (s_a < 0.0) then ! This is the usual ocean case
554 if (s_b < 0.0) then ! This is almost always the case
555 sbdry(i,j) = 2.0*s_c / (-s_b + sqrt(s_b*s_b - 4.*s_a*s_c))
556 else
557 sbdry(i,j) = (s_b + sqrt(s_b*s_b - 4.*s_a*s_c)) / (-2.*s_a)
558 endif
559 elseif ((s_a == 0.0) .and. (s_b < 0.0)) then ! It should be the case that S_b < 0.
560 sbdry(i,j) = -s_c / s_b
561 else
562 call mom_error(fatal, "Impossible conditions found in 3-equation skin salinity calculation.")
563 endif
564
565 ! Safety check
566 if (sbdry(i,j) < 0.) then
567 write(mesg,*) 'sfc_state%sss(i,j) = ',us%S_to_ppt*sfc_state%sss(i,j), &
568 'S_a, S_b, S_c', us%ppt_to_S*s_a, s_b, us%S_to_ppt*s_c
569 call mom_error(warning, mesg, .true.)
570 call mom_error(fatal, "shelf_calc_flux: Negative salinity (Sbdry).")
571 endif
572 else
573 ! Guess sss as the iteration starting point for the boundary salinity.
574 sbdry(i,j) = sfc_state%sss(i,j) ; sb_max_set = .false.
575 sb_min_set = .false.
576 endif !find_salt_root
577
578 do it1 = 1,20
579 ! Determine the potential temperature at the ice-ocean interface.
580 ! The following two lines are equivalent:
581 ! call calculate_TFreeze(Sbdry(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state, scale_from_EOS=.true.)
582 call calculate_tfreeze(sbdry(i:i,j), p_int(i:i), iss%tfreeze(i:i,j), cs%eqn_of_state)
583
584 dt_ustar = (iss%tfreeze(i,j) - sfc_state%sst(i,j)) * ustar_h
585 ds_ustar = (sbdry(i,j) - sfc_state%sss(i,j)) * ustar_h
586
587 if (cs%const_gamma) then
588 ! If using a constant gamma_T, there are no effects of the buoyancy flux on the turbulence.
589 i_gam_t = cs%Gamma_T_3EQ
590 i_gam_s = cs%Gamma_S_3EQ
591 wt_flux = dt_ustar * cs%Gamma_T_3EQ
592 wb_flux = db_ds * (ds_ustar * cs%Gamma_S_3EQ) + db_dt * wt_flux
593 elseif (.not.cs%buoy_flux_itt_bugfix) then
594 ! Gamma_T and gamma_S are a function of the buoyancy flux, and there should have been
595 ! iteration to find the root where wB_flux is consistent with the values of gamma with
596 ! that flux, but it was omitted.
597 gam_turb = i_vk * (ln_neut + (i_2zeta_n - 1.0))
598 i_gam_t = 1.0 / (gam_mol_t + gam_turb)
599 i_gam_s = 1.0 / (gam_mol_s + gam_turb)
600 wb_flux = db_ds * (ds_ustar * i_gam_s) + db_dt * (dt_ustar * i_gam_t)
601
602 if (wb_flux < 0.0) then ! The stabilising buoyancy flux reduces the turbulent fluxes.
603 i_n_star = sqrt(1.0 - n_star_term * wb_flux)
604 if (hbl_neut_h_molec > i_n_star**2) then
605 gam_turb = i_vk * ((ln_neut - 2.0*log(i_n_star)) + (i_2zeta_n*i_n_star - 1.0))
606 else ! The layer dominated by molecular viscosity is smaller than the boundary layer.
607 gam_turb = i_vk * (i_2zeta_n*i_n_star - 1.0)
608 endif
609 i_gam_t = 1.0 / (gam_mol_t + gam_turb)
610 i_gam_s = 1.0 / (gam_mol_s + gam_turb)
611 endif
612 wt_flux = dt_ustar * i_gam_t
613 else ! gamma_T and gamma_S are a function of the buoyancy flux with proper iteration.
614 ! Find the root where wB_flux is consistent with the values of gamma with that flux.
615
616 ! First, determine the buoyancy flux assuming no effects of stability
617 ! on the turbulence. Following H & J '99, this limit also applies
618 ! when the buoyancy flux is destabilizing.
619 gam_turb = i_vk * (ln_neut + (i_2zeta_n - 1.0))
620 i_gam_t = 1.0 / (gam_mol_t + gam_turb)
621 i_gam_s = 1.0 / (gam_mol_s + gam_turb)
622 wb_flux = (db_ds * ds_ustar) * i_gam_s + (db_dt * dt_ustar) * i_gam_t
623
624 if (wb_flux < 0.0) then
625 ! The buoyancy flux is stabilizing and will reduce the turbulent
626 ! fluxes, and iteration is required.
627
628 ! n_star <= 1.0 is the ratio of working boundary layer thickness
629 ! to the neutral thickness. I_n_star is its inverse.
630 i_n_star = sqrt(1.0 - n_star_term * wb_flux)
631 if (hbl_neut_h_molec > i_n_star**2) then
632 gam_turb = i_vk * ((ln_neut - 2.0*log(i_n_star)) + (i_2zeta_n*i_n_star - 1.0))
633 else ! The layer dominated by molecular viscosity is smaller than the boundary layer.
634 gam_turb = i_vk * (i_2zeta_n*i_n_star - 1.0)
635 endif
636 i_gam_t = 1.0 / (gam_mol_t + gam_turb)
637 i_gam_s = 1.0 / (gam_mol_s + gam_turb)
638
639 wb_flux_new = (db_ds * ds_ustar) * i_gam_s + (db_dt * dt_ustar) * i_gam_t
640 root_found = (abs(wb_flux_new - wb_flux) < cs%buoy_flux_tol*(abs(wb_flux_new) + abs(wb_flux)))
641 ! Do not update the flux if its maagnitude would be increased by the otherwise
642 ! stabilizing buoyancy fluxes. This can happen when the buoyancy flux
643 ! is stabilizing when one of the heat or salt fluxes are destabilizing due
644 ! to their different molecular properties.
645 if (wb_flux_new <= wb_flux) root_found = .true.
646
647 if (.not.root_found) then
648 wb_flux_max = 0.0 ; dwb_max = wb_flux
649 wb_flux_min = wb_flux ; dwb_min = wb_flux_new - wb_flux
650
651 if ((wb_flux_min*n_star_term < (1.0 - hbl_neut_h_molec)) .and. &
652 ((1.0 - hbl_neut_h_molec) < wb_flux_max*n_star_term)) then
653 ! The derivative of Gam_turb with wB_flux has a discontinuous change within the
654 ! bracketed range of values. Take this discontinous slope value for a first
655 ! guess, because Newton's method and the false position method may not converge
656 ! quickly when this discontinuity is between a guess and the solution.
657 wb_flux = (1.0 - hbl_neut_h_molec) / n_star_term
658 i_n_star = sqrt(hbl_neut_h_molec)
659 gam_turb = i_vk * (i_2zeta_n*i_n_star - 1.0)
660 i_gam_t = 1.0 / (gam_mol_t + gam_turb)
661 i_gam_s = 1.0 / (gam_mol_s + gam_turb)
662 wb_flux_new = (db_ds * ds_ustar) * i_gam_s + (db_dt * dt_ustar) * i_gam_t
663
664 if (abs(wb_flux_new - wb_flux) <= cs%buoy_flux_tol*(abs(wb_flux_new) + abs(wb_flux))) then
665 ! The root has been found to within the tolerance at the kink. This should be very rare.
666 root_found = .true.
667 elseif (wb_flux_new > wb_flux) then
668 ! The solution is in the limit where abs(wB_flux) is small and
669 ! Gam_turb = I_VK * ((ln_neut - 2.0*log(I_n_star)) + (I_2Zeta_N*I_n_star - 1.0))
670 wb_flux_min = wb_flux ; dwb_min = wb_flux_new - wb_flux
671 else
672 ! The solution is in the limt where abs(wB_flux) is large and
673 ! Gam_turb = I_VK * (I_2Zeta_N*I_n_star - 1.0)
674 wb_flux_max = wb_flux ; dwb_max = wb_flux_new - wb_flux
675 endif
676 endif
677 endif
678
679 if (.not.root_found) then
680 ! Use the false position for the next guess.
681 wb_flux = wb_flux_min + (wb_flux_max-wb_flux_min) * (dwb_min / (dwb_min - dwb_max))
682
683 do it3 = 1,30
684 ! Iterate using Newton's method with bounds or the false position method to find the root.
685
686 i_n_star = sqrt(1.0 - n_star_term * wb_flux)
687 dins_dwb = -0.5 * n_star_term / i_n_star
688 if (hbl_neut_h_molec > i_n_star**2) then
689 gam_turb = i_vk * ((ln_neut - 2.0*log(i_n_star)) + (i_2zeta_n*i_n_star - 1.0))
690 dg_dwb = i_vk * (( -2.0 / i_n_star + i_2zeta_n) * dins_dwb)
691 else
692 ! The layer dominated by molecular viscosity is smaller than the boundary layer.
693 gam_turb = i_vk * (i_2zeta_n*i_n_star - 1.0)
694 dg_dwb = i_vk * (i_2zeta_n * dins_dwb)
695 endif
696 i_gam_t = 1.0 / (gam_mol_t + gam_turb)
697 i_gam_s = 1.0 / (gam_mol_s + gam_turb)
698 wb_flux_new = (db_ds * ds_ustar) * i_gam_s + (db_dt * dt_ustar) * i_gam_t
699
700 ! Test for convergence to within tolerance at the point where wB_flux_new = wB_flux.
701 if (abs(wb_flux_new - wb_flux) <= cs%buoy_flux_tol*(abs(wb_flux_new) + abs(wb_flux))) &
702 root_found = .true.
703 if (root_found) exit
704
705 ddwb_dwb = -dg_dwb * ((db_ds * ds_ustar) * i_gam_s**2 + &
706 (db_dt * dt_ustar) * i_gam_t**2) - 1.0
707 if ((ddwb_dwb >= 0.0) .or. &
708 ( wb_flux - wb_flux_new >= abs(ddwb_dwb)*(wb_flux_max - wb_flux)) .or. &
709 ( wb_flux - wb_flux_new <= abs(ddwb_dwb)*(wb_flux_min - wb_flux)) ) then
710 ! Use the False position method to determine the guess for the next iteration when
711 ! Newton's method would go out of bounds
712 wb_flux_next = wb_flux_min + (wb_flux_max-wb_flux_min) * (dwb_min / (dwb_min - dwb_max))
713 else
714 ! Use Newton's method for the next guess.
715 wb_flux_next = wb_flux - (wb_flux_new - wb_flux) / ddwb_dwb
716 endif
717
718 ! Reset one of the bounds inward.
719 if (wb_flux_new - wb_flux > 0) then
720 wb_flux_min = wb_flux ; dwb_min = wb_flux_new - wb_flux
721 else
722 wb_flux_max = wb_flux ; dwb_max = wb_flux_new - wb_flux
723 endif
724
725 ! Update wB_flux
726 wb_flux = wb_flux_next
727 enddo ! it3
728 endif
729
730 endif ! End of test for first guess of wB_flux < 0.
731 wt_flux = dt_ustar * i_gam_t
732 endif ! End of test for CS%const_gamma
733
734 iss%tflux_ocn(i,j) = rhocp * wt_flux
735 exch_vel_t(i,j) = ustar_h * i_gam_t
736 exch_vel_s(i,j) = ustar_h * i_gam_s
737
738 ! Calculate the heat flux inside the ice shelf.
739 ! Vertical adv/diff as in H+J 1999, equations (26) & approx from (31).
740 ! Q_ice = density_ice * CS%Cp_ice * K_ice * dT/dz (at interface)
741 ! vertical adv/diff as in H+J 1999, equations (31) & (26)...
742 ! dT/dz ~= min( (lprec/(density_ice*K_ice))*(CS%Temp_Ice-T_freeze) , 0.0 )
743 ! If this approximation is not made, iterations are required... See H+J Fig 3.
744
745 if (iss%tflux_ocn(i,j) >= 0.0) then
746 ! Freezing occurs due to downward ocean heat flux, so zero iout ce heat flux.
747 iss%water_flux(i,j) = -i_lf * iss%tflux_ocn(i,j)
748 iss%tflux_shelf(i,j) = 0.0
749 else
750 if (cs%insulator) then
751 !no conduction/perfect insulator
752 iss%tflux_shelf(i,j) = 0.0
753 iss%water_flux(i,j) = i_lf * (iss%tflux_shelf(i,j) - iss%tflux_ocn(i,j))
754
755 else
756 ! With melting, from H&J 1999, eqs (31) & (26)...
757 ! Q_ice ~= Cp_ice * (CS%Temp_Ice-T_freeze) * lprec
758 ! RhoLF*lprec = Q_ice - ISS%tflux_ocn(i,j)
759 ! lprec = -(ISS%tflux_ocn(i,j)) / (CS%Lat_fusion + Cp_ice * (T_freeze-CS%Temp_Ice))
760 iss%water_flux(i,j) = -iss%tflux_ocn(i,j) / &
761 (cs%Lat_fusion + cs%Cp_ice * (iss%tfreeze(i,j) - cs%Temp_Ice))
762
763 iss%tflux_shelf(i,j) = iss%tflux_ocn(i,j) + cs%Lat_fusion*iss%water_flux(i,j)
764 endif
765
766 endif
767 !other options: dTi/dz linear through shelf, with draft in [Z ~> m], KTI in [Z2 T-1 ~> m2 s-1]
768 ! dTi_dz = (CS%Temp_Ice - ISS%tfreeze(i,j)) / draft(i,j)
769 ! ISS%tflux_shelf(i,j) = Rho_Ice * CS%Cp_ice * KTI * dTi_dz
770
771
772 if (cs%find_salt_root) then
773 exit ! no need to do interaction, so exit loop
774 else
775
776 mass_exch = exch_vel_s(i,j) * cs%Rho_ocn
777 sbdry_it = (sfc_state%sss(i,j) * mass_exch + cs%Salin_ice * iss%water_flux(i,j)) / &
778 (mass_exch + iss%water_flux(i,j))
779 ds_it = sbdry_it - sbdry(i,j)
780 if (abs(ds_it) < 1.0e-4*(0.5*(sfc_state%sss(i,j) + sbdry(i,j) + 1.0e-10*us%ppt_to_S))) exit
781
782 if (ds_it < 0.0) then ! Sbdry is now the upper bound.
783 if (sb_max_set) then
784 if (sbdry(i,j) > sb_max) &
785 call mom_error(fatal,"shelf_calc_flux: Irregular iteration for Sbdry (max).")
786 endif
787 sb_max = sbdry(i,j) ; ds_max = ds_it ; sb_max_set = .true.
788 else ! Sbdry is now the lower bound.
789 if (sb_min_set) then
790 if (sbdry(i,j) < sb_min) &
791 call mom_error(fatal, "shelf_calc_flux: Irregular iteration for Sbdry (min).")
792 endif
793 sb_min = sbdry(i,j) ; ds_min = ds_it ; sb_min_set = .true.
794 endif ! dS_it < 0.0
795
796 if (sb_min_set .and. sb_max_set) then
797 ! Use the false position method for the next iteration.
798 sbdry(i,j) = sb_min + (sb_max-sb_min) * (ds_min / (ds_min - ds_max))
799 else
800 sbdry(i,j) = sbdry_it
801 endif ! Sb_min_set
802
803 if (.not.cs%salt_flux_itt_bugfix) sbdry(i,j) = sbdry_it
804
805 endif ! CS%find_salt_root
806
807 enddo !it1
808 ! Check for non-convergence and/or non-boundedness?
809
810 else
811 ! In the 2-equation form, the mixed layer turbulent exchange velocity
812 ! is specified and large enough that the ocean salinity at the interface
813 ! is about the same as the boundary layer salinity.
814 ! The following two lines are equivalent:
815 ! call calculate_TFreeze(Sbdry(i,j), p_int(i), ISS%tfreeze(i,j), CS%eqn_of_state, scale_from_EOS=.true.)
816 call calculate_tfreeze(sfc_state%SSS(i:i,j), p_int(i:i), iss%tfreeze(i:i,j), cs%eqn_of_state)
817
818 exch_vel_t(i,j) = cs%gamma_t
819 iss%tflux_ocn(i,j) = rhocp * exch_vel_t(i,j) * (iss%tfreeze(i,j) - sfc_state%sst(i,j))
820 iss%tflux_shelf(i,j) = 0.0
821 iss%water_flux(i,j) = -i_lf * iss%tflux_ocn(i,j)
822 sbdry(i,j) = 0.0
823 endif
824 elseif (iss%area_shelf_h(i,j) > 0.0) then ! This is an ice-sheet, not a floating shelf.
825 iss%tflux_ocn(i,j) = 0.0
826 else ! There is no ice shelf or sheet here.
827 iss%tflux_ocn(i,j) = 0.0
828 endif
829
830! haline_driving(i,j) = sfc_state%sss(i,j) - Sbdry(i,j)
831
832 enddo ! i-loop
833 enddo ! j-loop
834
835 if (allocated(sfc_state%frazil)) then
836 add_frazil = .true.
837 else
838 add_frazil = .false.
839 endif
840
841 do j=js,je ; do i=is,ie
842 ! ISS%water_flux = net liquid water into the ocean [R Z T-1 ~> kg m-2 s-1]
843 if (cs%flux_factor/=1.0) then
844 iss%water_flux(i,j) = iss%water_flux(i,j) * cs%flux_factor
845 iss%tflux_ocn(i,j) = iss%tflux_ocn(i,j) * cs%flux_factor
846 if (cs%threeeq .and. iss%tflux_ocn(i,j) < 0.0 .and. (.not. cs%insulator)) &
847 iss%tflux_shelf(i,j)=iss%tflux_ocn(i,j) + cs%Lat_fusion * iss%water_flux(i,j)
848 endif
849
850 if ((sfc_state%ocean_mass(i,j) > cs%col_mass_melt_threshold) .and. &
851 (iss%area_shelf_h(i,j) > 0.0) .and. (cs%isthermo)) then
852
853 ! Set melt to zero above a cutoff pressure (CS%Rho_ocn*CS%cutoff_depth*CS%g_Earth).
854 ! This is needed for the ISOMIP test case.
855 if (iss%mass_shelf(i,j) < cs%Rho_ocn*cs%cutoff_depth) then
856 iss%water_flux(i,j) = 0.0
857 endif
858 ! Compute haline driving, which is one of the diags. used in ISOMIP
859 if (exch_vel_s(i,j)>0.) haline_driving(i,j) = (iss%water_flux(i,j) * sbdry(i,j)) / (cs%Rho_ocn * exch_vel_s(i,j))
860
861 !!!!!!!!!!!!!!!!!!!!!!!!!!!!Safety checks !!!!!!!!!!!!!!!!!!!!!!!!!
862 !1)Check if haline_driving computed above is consistent with
863 ! haline_driving = sfc_state%sss - Sbdry
864 !if (ISS%water_flux(i,j) /= 0.0) then
865 ! if (haline_driving(i,j) /= (sfc_state%sss(i,j) - Sbdry(i,j))) then
866 ! write(mesg,*) 'at i,j=',i,j,' haline_driving, sss-Sbdry',US%S_to_ppt*haline_driving(i,j), &
867 ! US%S_to_ppt*(sfc_state%sss(i,j) - Sbdry(i,j))
868 ! call MOM_error(FATAL, &
869 ! "shelf_calc_flux: Inconsistency in melt and haline_driving"//trim(mesg))
870 ! endif
871 !endif
872
873 ! 2) check if |melt| > 0 when ustar_shelf = 0.
874 ! this should never happen
875 if ((abs(iss%water_flux(i,j))>0.0) .and. (fluxes%ustar_shelf(i,j) == 0.0)) then
876 write(mesg,*) "|melt| = ",iss%water_flux(i,j)," > 0 and ustar_shelf = 0. at i,j", i, j
877 call mom_error(fatal, "shelf_calc_flux: "//trim(mesg))
878 endif
879 !!!!!!!!!!!!!!!!!!!!!!!!!!!!End of safety checks !!!!!!!!!!!!!!!!!!!
880 elseif (iss%area_shelf_h(i,j) > 0.0) then
881 ! This is grounded ice, that could be modified to melt if a geothermal heat flux were used.
882 haline_driving(i,j) = 0.0
883 iss%water_flux(i,j) = 0.0
884 endif ! area_shelf_h
885
886 ! mass flux [R Z L2 T-1 ~> kg s-1], part of ISOMIP diags.
887 mass_flux(i,j) = iss%water_flux(i,j) * iss%area_shelf_h(i,j)
888
889 !Add frazil formation
890 if (add_frazil .and. (iss%hmask(i,j) == 1 .or. iss%hmask(i,j) == 2)) &
891 iss%water_flux(i,j) = iss%water_flux(i,j) - iss%frazil(i,j) * i_dt_lhf
892 fluxes%iceshelf_melt(i,j) = iss%water_flux(i,j)
893 enddo ; enddo ! i- and j-loops
894
895 if (cs%active_shelf_dynamics .or. cs%override_shelf_movement) then
896 call cpu_clock_begin(id_clock_pass)
897 call pass_var(iss%area_shelf_h, g%domain, complete=.false.)
898 call pass_var(iss%mass_shelf, g%domain)
899 call cpu_clock_end(id_clock_pass)
900 endif
901
902 ! Melting has been computed, now is time to update thickness and mass
903 if ( cs%override_shelf_movement .and. (.not.cs%mass_from_file)) then
904 if (cs%bmb_diag) dh_bdott(is:ie,js:je) = iss%h_shelf(is:ie,js:je)
905 call change_thickness_using_melt(iss, g, us, time_step, fluxes, cs%density_ice, cs%debug)
906 if (cs%bmb_diag) dh_bdott(is:ie,js:je) = iss%h_shelf(is:ie,js:je) - dh_bdott(is:ie,js:je)
907
908 if (cs%debug) then
909 call hchksum(iss%h_shelf, "h_shelf after change thickness using melt", g%HI, haloshift=0, unscale=us%Z_to_m)
910 call hchksum(iss%mass_shelf, "mass_shelf after change thickness using melt", g%HI, haloshift=0, &
911 unscale=us%RZ_to_kg_m2)
912 endif
913 endif
914
915 ! Melting has been computed, now is time to update thickness and mass with dynamic ice shelf
916 if (cs%active_shelf_dynamics) then
917
918 iss%dhdt_shelf(:,:) = iss%h_shelf(:,:)
919
920 if (cs%bmb_diag) dh_bdott(is:ie,js:je) = iss%h_shelf(is:ie,js:je)
921 call change_thickness_using_melt(iss, g, us, time_step, fluxes, cs%density_ice, cs%debug)
922 if (cs%bmb_diag) dh_bdott(is:ie,js:je) = iss%h_shelf(is:ie,js:je) - dh_bdott(is:ie,js:je)
923
924 if (cs%debug) then
925 call hchksum(iss%h_shelf, "h_shelf after change thickness using melt", g%HI, haloshift=0, unscale=us%Z_to_m)
926 call hchksum(iss%mass_shelf, "mass_shelf after change thickness using melt", g%HI, haloshift=0, &
927 unscale=us%RZ_to_kg_m2)
928 endif
929
930 if (cs%smb_diag) dh_adott(is:ie,js:je) = iss%h_shelf(is:ie,js:je)
931 call change_thickness_using_precip(cs, iss, g, us, fluxes, time_step, time)
932 if (cs%smb_diag) dh_adott(is:ie,js:je) = iss%h_shelf(is:ie,js:je) - dh_adott(is:ie,js:je)
933
934 if (cs%debug) then
935 call hchksum(iss%h_shelf, "h_shelf after change thickness using surf acc", g%HI, haloshift=0, unscale=us%Z_to_m)
936 call hchksum(iss%mass_shelf, "mass_shelf after change thickness using surf acc", g%HI, haloshift=0, &
937 unscale=us%RZ_to_kg_m2)
938 endif
939
940 update_ice_vel = .false.
941 coupled_gl = (cs%GL_couple .and. .not. cs%solo_ice_sheet)
942
943 ! advect the ice shelf, and advance the front. Calving will be in here somewhere as well..
944 ! when we decide on how to do it
945 call update_ice_shelf(cs%dCS, iss, g, us, time_step, time, cs%calve_ice_shelf_bergs, &
946 sfc_state%ocean_mass, coupled_gl)
947
948 do j=js,je ; do i=is,ie
949 iss%dhdt_shelf(i,j) = (iss%h_shelf(i,j) - iss%dhdt_shelf(i,j))*itime_step
950 enddo ; enddo
951
952 call is_dynamics_post_data(time_step, time, cs%dCS, iss, g)
953 endif
954
955 if (cs%shelf_mass_is_dynamic) &
956 call write_ice_shelf_energy(cs%dCS, g, us, iss%mass_shelf, iss%area_shelf_h, time, &
957 time_step=real_to_time(time_step, unscale=us%T_to_s) )
958
959 if (cs%debug) call mom_forcing_chksum("Before add shelf flux", fluxes, g, cs%US, haloshift=0)
960
961 ! pass on the updated ice sheet geometry (for pressure on ocean) and thermodynamic data
962 call add_shelf_flux(g, us, cs, sfc_state, fluxes, time_step)
963
964 call enable_averages(time_step, time, cs%diag)
965 if (cs%id_shelf_mass > 0) call post_data(cs%id_shelf_mass, iss%mass_shelf, cs%diag)
966 if (cs%id_area_shelf_h > 0) call post_data(cs%id_area_shelf_h, iss%area_shelf_h, cs%diag)
967 if (cs%id_ustar_shelf > 0) call post_data(cs%id_ustar_shelf, fluxes%ustar_shelf, cs%diag)
968 if (cs%id_shelf_sfc_mass_flux > 0) call post_data(cs%id_shelf_sfc_mass_flux, fluxes%shelf_sfc_mass_flux, cs%diag)
969
970 if (cs%id_melt > 0) call post_data(cs%id_melt, fluxes%iceshelf_melt, cs%diag)
971 if (cs%id_thermal_driving > 0) call post_data(cs%id_thermal_driving, (sfc_state%sst-iss%tfreeze), cs%diag)
972 if (cs%id_Sbdry > 0) call post_data(cs%id_Sbdry, sbdry, cs%diag)
973 if (cs%id_haline_driving > 0) call post_data(cs%id_haline_driving, haline_driving, cs%diag)
974 if (cs%id_mass_flux > 0) call post_data(cs%id_mass_flux, mass_flux, cs%diag)
975 if (cs%id_u_ml > 0) call post_data(cs%id_u_ml, sfc_state%u, cs%diag)
976 if (cs%id_v_ml > 0) call post_data(cs%id_v_ml, sfc_state%v, cs%diag)
977 if (cs%id_tfreeze > 0) call post_data(cs%id_tfreeze, iss%tfreeze, cs%diag)
978 if (cs%id_tfl_shelf > 0) call post_data(cs%id_tfl_shelf, iss%tflux_shelf, cs%diag)
979 if (cs%id_exch_vel_t > 0) call post_data(cs%id_exch_vel_t, exch_vel_t, cs%diag)
980 if (cs%id_exch_vel_s > 0) call post_data(cs%id_exch_vel_s, exch_vel_s, cs%diag)
981 if (cs%id_h_shelf > 0) call post_data(cs%id_h_shelf, iss%h_shelf, cs%diag)
982 if (cs%id_dhdt_shelf > 0) call post_data(cs%id_dhdt_shelf, iss%dhdt_shelf, cs%diag)
983 if (cs%id_h_mask > 0) call post_data(cs%id_h_mask,iss%hmask,cs%diag)
984 if (cs%id_frazil > 0) call post_data(cs%id_frazil,iss%frazil,cs%diag)
985 if (cs%active_shelf_dynamics) &
986 call process_and_post_scalar_data(cs, vaf0, vaf0_a, vaf0_g, itime_step, dh_adott, dh_bdott)
987 call disable_averaging(cs%diag)
988
989 !reset used frazil
990 if (add_frazil) iss%frazil(:,:) = 0.0
991
992 call cpu_clock_end(id_clock_shelf)
993
994 if (cs%debug) call mom_forcing_chksum("End of shelf calc flux", fluxes, g, cs%US, haloshift=0)
995
996 if (cs%rotate_index) then
997! call rotate_surface_state(sfc_state, sfc_state_in, CS%Grid_in, -CS%turns)
998 call rotate_forcing(fluxes, fluxes_in, -cs%turns)
999 call deallocate_surface_state(sfc_state)
1000 deallocate(sfc_state)
1001 call deallocate_forcing_type(fluxes)
1002 deallocate(fluxes)
1003 endif
1004
1005end subroutine shelf_calc_flux
1006
1007!> Copies frazil from the ocean surface state to the ice sheet state. Removes frazil that will
1008!! be used by the ice sheet from the ocean surface state
1009subroutine adjust_ice_sheet_frazil(sfc_state_in, fluxes_in, CS)
1010 type(surface), target, intent(inout) :: sfc_state_in !< A structure containing fields that
1011 !! describe the surface state of the ocean. The
1012 !! intent is only inout to allow for halo updates.
1013 type(forcing), target, intent(in) :: fluxes_in !< structure containing pointers to any
1014 !! possible thermodynamic or mass-flux forcing fields.
1015 type(ice_shelf_cs), pointer :: cs !< A pointer to the control structure returned
1016 !! by a previous call to initialize_ice_shelf.
1017 ! Local variables
1018 type(ocean_grid_type), pointer :: g => null() !< The grid structure used by the ice shelf.
1019 type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
1020 !! the ice-shelf state
1021 type(surface), pointer :: sfc_state => null()
1022 type(forcing), pointer :: fluxes => null()
1023 integer :: i,j,is,ie,js,je
1024
1025 g => cs%grid ; iss => cs%ISS
1026
1027 if (cs%rotate_index) then
1028 allocate(sfc_state)
1029 call rotate_surface_state(sfc_state_in, sfc_state, g, cs%turns)
1030 allocate(fluxes)
1031 call allocate_forcing_type(fluxes_in, g, fluxes, turns=cs%turns)
1032 call rotate_forcing(fluxes_in, fluxes, cs%turns)
1033 else
1034 sfc_state => sfc_state_in
1035 fluxes => fluxes_in
1036 endif
1037
1038 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1039
1040 do j=js,je ; do i=is,ie
1041 !Copy frazil to the ice sheet module where ice sheet is present.
1042 !No scaling to account for partial ice-sheet cells is necessary here, as
1043 !this is taken care of when applied to the ice sheet.
1044 if (fluxes%frac_shelf_h(i,j)>0.0) iss%frazil(i,j) = sfc_state%frazil(i,j)
1045 !Remove the frazil that is used by the ice sheet from sfc_state%frazil
1046 !The sfc_state%frazil is sent to the sea-ice module
1047 sfc_state%frazil(i,j) = sfc_state%frazil(i,j) * (1.0-fluxes%frac_shelf_h(i,j))
1048 enddo ; enddo
1049
1050 if (cs%rotate_index) then
1051 call rotate_surface_state(sfc_state, sfc_state_in, g, -cs%turns)
1052 ! call rotate_forcing(fluxes, fluxes_in, -CS%turns)
1053 call deallocate_surface_state(sfc_state)
1054 deallocate(sfc_state)
1055 call deallocate_forcing_type(fluxes)
1056 deallocate(fluxes)
1057 endif
1058end subroutine adjust_ice_sheet_frazil
1059
1060function integrate_over_ice_sheet_area(G, ISS, var, unscale, hemisphere) result(var_out)
1061 type(ocean_grid_type), intent(in) :: g !< The grid structure used by the ice shelf.
1062 type(ice_shelf_state), intent(in) :: iss !< A structure with elements that describe the ice-shelf state
1063 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: var !< Ice variable to integrate in arbitrary units [A ~> a]
1064 real, intent(in) :: unscale !< Dimensional scaling for variable to integrate [a A-1 ~> 1]
1065 integer, optional, intent(in) :: hemisphere !< 0 for Antarctica only, 1 for Greenland only. Otherwise, all ice sheets
1066 real :: var_out !< Variable integrated over the area of the ice sheet in arbitrary scaled units [A L2 ~> a m2]
1067
1068 ! Local variables
1069 integer :: is_id ! local copy of hemisphere
1070 real, dimension(SZI_(G),SZJ_(G)) :: var_cell !< Variable integrated over the ice-sheet area of each cell
1071 !! in arbitrary units [A L2 ~> a m2]
1072 integer, dimension(SZI_(G),SZJ_(G)) :: mask ! a mask for active cells depending on hemisphere indicated
1073 integer :: i, j
1074
1075 if (present(hemisphere)) then
1076 is_id = hemisphere
1077 else
1078 is_id = -1
1079 endif
1080
1081 mask(:,:) = 0
1082 if (is_id==0) then !Antarctica (S. Hemisphere) only
1083 do j = g%jsc,g%jec ; do i = g%isc,g%iec
1084 if (iss%hmask(i,j)>0 .and. g%geoLatT(i,j)<=0.0) mask(i,j)=1
1085 enddo ; enddo
1086 elseif (is_id==1) then !Greenland (N. Hemisphere) only
1087 do j = g%jsc,g%jec ; do i = g%isc,g%iec
1088 if (iss%hmask(i,j)>0 .and. g%geoLatT(i,j)>0.0) mask(i,j)=1
1089 enddo ; enddo
1090 else !All ice sheets
1091 mask(g%isc:g%iec,g%jsc:g%jec) = iss%hmask(g%isc:g%iec,g%jsc:g%jec)
1092 endif
1093
1094 var_cell(:,:) = 0.0
1095 do j = g%jsc,g%jec ; do i = g%isc,g%iec
1096 if (mask(i,j)>0) var_cell(i,j) = var(i,j) * iss%area_shelf_h(i,j)
1097 enddo ; enddo
1098
1099 var_out = reproducing_sum(var_cell, unscale=unscale*g%US%L_to_m**2)
1101
1102!> Converts the ice-shelf-to-ocean calving and calving_hflx variables from the ice-shelf state (ISS) type
1103!! to the ocean public type
1104subroutine ice_sheet_calving_to_ocean_sfc(CS,US,calving,calving_hflx)
1105 type(ice_shelf_cs), pointer :: cs !< A pointer to the ice shelf control structure
1106 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1107 real, dimension(:,:), intent(inout) :: calving !< The mass flux per unit area of the ice shelf
1108 !! to convert to bergs [R Z T-1 ~> kg m-2 s-1].
1109 real, dimension(:,:), intent(inout) :: calving_hflx !< Calving heat flux [Q R Z T-1 ~> W m-2].
1110 ! Local variables
1111 type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
1112 !! the ice-shelf state
1113 type(ocean_grid_type), pointer :: g => null() !< A pointer to the ocean grid metric.
1114 integer :: is, ie, js, je
1115
1116 g=>cs%Grid
1117 iss => cs%ISS
1118 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1119
1120 calving = us%RZ_T_to_kg_m2s * iss%calving(is:ie,js:je)
1121 calving_hflx = us%QRZ_T_to_W_m2 * iss%calving_hflx(is:ie,js:je)
1122
1123 !CS%calve_ice_shelf_bergs=.true.
1124
1125end subroutine ice_sheet_calving_to_ocean_sfc
1126
1127!> Changes the thickness (mass) of the ice shelf based on sub-ice-shelf melting
1128subroutine change_thickness_using_melt(ISS, G, US, time_step, fluxes, density_ice, debug)
1129 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1130 type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
1131 !! the ice-shelf state
1132 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1133 real, intent(in) :: time_step !< The time step for this update [T ~> s].
1134 type(forcing), intent(inout) :: fluxes !< structure containing pointers to any possible
1135 !! thermodynamic or mass-flux forcing fields.
1136 real, intent(in) :: density_ice !< The density of ice-shelf ice [R ~> kg m-3].
1137 logical, optional, intent(in) :: debug !< If present and true, write chksums
1138
1139 ! locals
1140 real :: I_rho_ice ! Ice specific volume [R-1 ~> m3 kg-1]
1141 integer :: i, j
1142
1143 i_rho_ice = 1.0 / density_ice
1144
1145
1146 do j=g%jsc,g%jec ; do i=g%isc,g%iec
1147 if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
1148 ! first, zero out fluxes applied during previous time step
1149 if (associated(fluxes%lprec)) fluxes%lprec(i,j) = 0.0
1150 if (associated(fluxes%sens)) fluxes%sens(i,j) = 0.0
1151 if (associated(fluxes%frac_shelf_h)) fluxes%frac_shelf_h(i,j) = 0.0
1152 if (associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = 0.0
1153
1154 if (iss%water_flux(i,j) * time_step / density_ice < iss%h_shelf(i,j)) then
1155 iss%h_shelf(i,j) = iss%h_shelf(i,j) - iss%water_flux(i,j) * time_step / density_ice
1156 else
1157 ! the ice is about to melt away, so set thickness, area, and mask to zero
1158 ! NOTE: this is not mass conservative should maybe scale salt & heat flux for this cell
1159 iss%h_shelf(i,j) = 0.0
1160 iss%hmask(i,j) = 0.0
1161 iss%area_shelf_h(i,j) = 0.0
1162 endif
1163 iss%mass_shelf(i,j) = iss%h_shelf(i,j) * density_ice
1164 endif
1165 enddo ; enddo
1166
1167 call pass_var(iss%area_shelf_h, g%domain, complete=.false.)
1168 call pass_var(iss%h_shelf, g%domain, complete=.false.)
1169 call pass_var(iss%hmask, g%domain, complete=.false.)
1170 call pass_var(iss%mass_shelf, g%domain)
1171
1172end subroutine change_thickness_using_melt
1173
1174!> This subroutine adds the mechanical forcing fields and perhaps shelf areas, based on
1175!! the ice state in ice_shelf_CS.
1176subroutine add_shelf_forces(Ocn_grid, US, CS, forces_in, do_shelf_area, external_call)
1177 type(ocean_grid_type), intent(in) :: ocn_grid !< The ocean's grid structure.
1178 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1179 type(ice_shelf_cs), pointer :: cs !< This module's control structure.
1180 type(mech_forcing), target, intent(inout) :: forces_in !< A structure with the
1181 !! driving mechanical forces
1182 logical, optional, intent(in) :: do_shelf_area !< If true find the shelf-covered areas.
1183 logical, optional, intent(in) :: external_call !< If true the incoming forcing type
1184 !! is using the unrotated input grid and may need
1185 !! to be rotated.
1186 type(ocean_grid_type), pointer :: g => null() !< A pointer to the ocean grid metric.
1187 type(mech_forcing), pointer :: forces !< A structure with the driving mechanical forces
1188 real :: kv_rho_ice ! The viscosity of ice divided by its density [L4 T-1 R-1 Z-2 ~> m5 kg-1 s-1].
1189 real :: press_ice ! The pressure of the ice shelf per unit area of ocean (not ice) [R L2 T-2 ~> Pa].
1190 logical :: find_area ! If true find the shelf areas at u & v points.
1191 logical :: rotate = .false.
1192 type(ice_shelf_state), pointer :: iss => null() ! A structure with elements that describe
1193 ! the ice-shelf state
1194
1195 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
1196
1197 rotate = .false. ; if (present(external_call)) rotate = external_call
1198
1199 if (cs%rotate_index .and. rotate) then
1200 if ((ocn_grid%isc /= cs%Grid_in%isc) .or. (ocn_grid%iec /= cs%Grid_in%iec) .or. &
1201 (ocn_grid%jsc /= cs%Grid_in%jsc) .or. (ocn_grid%jec /= cs%Grid_in%jec)) &
1202 call mom_error(fatal,"add_shelf_forces: Incompatible Ocean and external Ice shelf grids.")
1203 allocate(forces)
1204 call allocate_mech_forcing(forces_in, cs%Grid, forces)
1205 call rotate_mech_forcing(forces_in, cs%turns, forces)
1206 else
1207 if ((ocn_grid%isc /= cs%Grid%isc) .or. (ocn_grid%iec /= cs%Grid%iec) .or. &
1208 (ocn_grid%jsc /= cs%Grid%jsc) .or. (ocn_grid%jec /= cs%Grid%jec)) &
1209 call mom_error(fatal,"add_shelf_forces: Incompatible Ocean and internal Ice shelf grids.")
1210
1211 forces=>forces_in
1212 endif
1213
1214 g=>cs%Grid
1215
1216 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1217 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
1218
1219 iss => cs%ISS
1220
1221 find_area = .true. ; if (present(do_shelf_area)) find_area = do_shelf_area
1222
1223 if (find_area) then
1224 ! The frac_shelf is set over the widest possible area. Could it be smaller?
1225 do j=jsd,jed ; do i=isd,ied-1
1226 forces%frac_shelf_u(i,j) = 0.0
1227 if ((g%areaT(i,j) + g%areaT(i+1,j) > 0.0)) & ! .and. (G%areaCu(I,j) > 0.0)) &
1228 forces%frac_shelf_u(i,j) = (iss%area_shelf_h(i,j) + iss%area_shelf_h(i+1,j)) / &
1229 (g%areaT(i,j) + g%areaT(i+1,j))
1230 enddo ; enddo
1231 do j=jsd,jed-1 ; do i=isd,ied
1232 forces%frac_shelf_v(i,j) = 0.0
1233 if ((g%areaT(i,j) + g%areaT(i,j+1) > 0.0)) & ! .and. (G%areaCv(i,J) > 0.0)) &
1234 forces%frac_shelf_v(i,j) = (iss%area_shelf_h(i,j) + iss%area_shelf_h(i,j+1)) / &
1235 (g%areaT(i,j) + g%areaT(i,j+1))
1236 enddo ; enddo
1237 call pass_vector(forces%frac_shelf_u, forces%frac_shelf_v, g%domain, to_all, cgrid_ne)
1238 endif
1239
1240 do j=js,je ; do i=is,ie
1241 press_ice = (iss%area_shelf_h(i,j) * g%IareaT(i,j)) * (cs%g_Earth * iss%mass_shelf(i,j))
1242 if (associated(forces%p_surf)) then
1243 if (.not.forces%accumulate_p_surf) forces%p_surf(i,j) = 0.0
1244 forces%p_surf(i,j) = forces%p_surf(i,j) + press_ice
1245 endif
1246 if (associated(forces%p_surf_full)) then
1247 if (.not.forces%accumulate_p_surf) forces%p_surf_full(i,j) = 0.0
1248 forces%p_surf_full(i,j) = forces%p_surf_full(i,j) + press_ice
1249 endif
1250 enddo ; enddo
1251
1252 ! For various reasons, forces%rigidity_ice_[uv] is always updated here. Note
1253 ! that it may have been zeroed out where IOB is translated to forces and
1254 ! contributions from icebergs and the sea-ice pack added subsequently.
1255 !### THE RIGIDITY SHOULD ALSO INCORPORATE AREAL-COVERAGE INFORMATION.
1256 kv_rho_ice = cs%kv_ice / cs%density_ice
1257 do j=js,je ; do i=is-1,ie
1258 if (.not.forces%accumulate_rigidity) forces%rigidity_ice_u(i,j) = 0.0
1259 forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + &
1260 kv_rho_ice * min(iss%mass_shelf(i,j), iss%mass_shelf(i+1,j))
1261 enddo ; enddo
1262 do j=js-1,je ; do i=is,ie
1263 if (.not.forces%accumulate_rigidity) forces%rigidity_ice_v(i,j) = 0.0
1264 forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + &
1265 kv_rho_ice * min(iss%mass_shelf(i,j), iss%mass_shelf(i,j+1))
1266 enddo ; enddo
1267
1268 if (cs%debug) then
1269 call uvchksum("rigidity_ice_[uv]", forces%rigidity_ice_u, &
1270 forces%rigidity_ice_v, cs%Grid%HI, symmetric=.true., &
1271 unscale=us%L_to_m**3*us%L_to_Z*us%s_to_T, scalar_pair=.true.)
1272 call uvchksum("frac_shelf_[uv]", forces%frac_shelf_u, &
1273 forces%frac_shelf_v, cs%Grid%HI, symmetric=.true., &
1274 scalar_pair=.true.)
1275 endif
1276
1277 if (cs%rotate_index .and. rotate) then
1278 call rotate_mech_forcing(forces, -cs%turns, forces_in)
1279 call deallocate_mech_forcing(forces)
1280 endif
1281
1282end subroutine add_shelf_forces
1283
1284!> This subroutine adds the ice shelf pressure to the fluxes type.
1285subroutine add_shelf_pressure(Ocn_grid, US, CS, fluxes)
1286 type(ocean_grid_type), intent(in) :: Ocn_grid !< The ocean's grid structure.
1287 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1288 type(ice_shelf_cs), intent(in) :: CS !< This module's control structure.
1289 type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be updated.
1290
1291 type(ocean_grid_type), pointer :: G => null() ! A pointer to ocean's grid structure.
1292 real :: press_ice !< The pressure of the ice shelf per unit area of ocean (not ice) [R L2 T-2 ~> Pa].
1293 integer :: i, j, is, ie, js, je
1294
1295 g=>cs%Grid
1296 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1297
1298 if ((cs%grid%isc /= g%isc) .or. (cs%grid%iec /= g%iec) .or. &
1299 (cs%grid%jsc /= g%jsc) .or. (cs%grid%jec /= g%jec)) &
1300 call mom_error(fatal,"add_shelf_pressure: Incompatible ocean and ice shelf grids.")
1301
1302 do j=js,je ; do i=is,ie
1303 press_ice = (cs%ISS%area_shelf_h(i,j) * g%IareaT(i,j)) * (cs%g_Earth * cs%ISS%mass_shelf(i,j))
1304 if (associated(fluxes%p_surf)) then
1305 if (.not.fluxes%accumulate_p_surf) fluxes%p_surf(i,j) = 0.0
1306 fluxes%p_surf(i,j) = fluxes%p_surf(i,j) + press_ice
1307 endif
1308 if (associated(fluxes%p_surf_full)) then
1309 if (.not.fluxes%accumulate_p_surf) fluxes%p_surf_full(i,j) = 0.0
1310 fluxes%p_surf_full(i,j) = fluxes%p_surf_full(i,j) + press_ice
1311 endif
1312 enddo ; enddo
1313
1314end subroutine add_shelf_pressure
1315
1316!> Updates surface fluxes that are influenced by sub-ice-shelf melting
1317subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes, time_step)
1318 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1319 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1320 type(ice_shelf_cs), pointer :: CS !< This module's control structure.
1321 type(surface), intent(inout) :: sfc_state !< Surface ocean state
1322 type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be used/updated.
1323 real, intent(in) :: time_step !< Time step over which fluxes are applied [T ~> s]
1324 ! local variables
1325 real :: frac_shelf !< The fractional area covered by the ice shelf [nondim].
1326 real :: frac_open !< The fractional area of the ocean that is not covered by the ice shelf [nondim].
1327 real :: delta_mass_shelf !< Change in ice shelf mass over one time step [R Z L2 T-1 ~> kg s-1]
1328 real :: balancing_flux !< The fresh water flux that balances the integrated melt flux [R Z T-1 ~> kg m-2 s-1]
1329 real :: balancing_area !< total area where the balancing flux is applied [L2 ~> m2]
1330 type(time_type) :: dTime !< The time step as a time_type
1331 type(time_type) :: Time0 !< The previous time (Time-dt)
1332 real, dimension(SZDI_(G),SZDJ_(G)) :: bal_frac !< Fraction of the cell where the mass flux
1333 !! balancing the net melt flux occurs, 0 to 1 [nondim]
1334 real, dimension(SZDI_(G),SZDJ_(G)) :: last_mass_shelf !< Ice shelf mass
1335 !! at at previous time (Time-dt) [R Z ~> kg m-2]
1336 real, dimension(SZDI_(G),SZDJ_(G)) :: delta_float_mass !< The change in the floating mass between
1337 !! the two timesteps at (Time) and (Time-dt) [R Z ~> kg m-2].
1338 real, dimension(SZDI_(G),SZDJ_(G)) :: last_h_shelf !< Ice shelf thickness [Z ~> m]
1339 !! at at previous time (Time-dt)
1340 real, dimension(SZDI_(G),SZDJ_(G)) :: last_hmask !< Ice shelf mask [nondim]
1341 !! at at previous time (Time-dt)
1342 real, dimension(SZDI_(G),SZDJ_(G)) :: last_area_shelf_h !< Ice shelf area [L2 ~> m2]
1343 !! at at previous time (Time-dt)
1344 real, dimension(SZDI_(G),SZDJ_(G)) :: delta_draft !< change in ice shelf draft thickness [L ~> m]
1345 !! since previous time (Time-dt)
1346 type(ice_shelf_state), pointer :: ISS => null() !< A structure with elements that describe
1347 !! the ice-shelf state
1348
1349 character(len=160) :: mesg ! The text of an error message
1350 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
1351 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1352 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
1353
1354 if ((cs%grid%isc /= g%isc) .or. (cs%grid%iec /= g%iec) .or. &
1355 (cs%grid%jsc /= g%jsc) .or. (cs%grid%jec /= g%jec)) &
1356 call mom_error(fatal,"add_shelf_flux: Incompatible ocean and ice shelf grids.")
1357
1358 iss => cs%ISS
1359
1360
1361 call add_shelf_pressure(g, us, cs, fluxes)
1362
1363 ! Determine ustar and the square magnitude of the velocity in the
1364 ! bottom boundary layer. Together these give the TKE source and
1365 ! vertical decay scale.
1366
1367 if (cs%debug) then
1368 if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then
1369 call uvchksum("tau[xy]_shelf", sfc_state%taux_shelf, sfc_state%tauy_shelf, &
1370 g%HI, haloshift=0, unscale=us%RZ_T_to_kg_m2s*us%L_T_to_m_s)
1371 endif
1372 endif
1373
1374 if (cs%active_shelf_dynamics .or. cs%override_shelf_movement) then
1375 do j=jsd,jed ; do i=isd,ied
1376 if (g%areaT(i,j) > 0.0) &
1377 fluxes%frac_shelf_h(i,j) = min(1.0, iss%area_shelf_h(i,j) * g%IareaT(i,j))
1378 enddo ; enddo
1379 endif
1380
1381 if (cs%debug) then
1382 call mom_forcing_chksum("Before adding shelf fluxes", fluxes, g, cs%US, haloshift=0)
1383 endif
1384
1385 do j=js,je ; do i=is,ie ; if (iss%area_shelf_h(i,j) > 0.0) then
1386 ! Replace fluxes intercepted by the ice shelf with fluxes from the ice shelf
1387 frac_shelf = min(1.0, iss%area_shelf_h(i,j) * g%IareaT(i,j))
1388 frac_open = max(0.0, 1.0 - frac_shelf)
1389
1390 if (associated(fluxes%sw)) fluxes%sw(i,j) = frac_open * fluxes%sw(i,j)
1391 if (associated(fluxes%sw_vis_dir)) fluxes%sw_vis_dir(i,j) = frac_open * fluxes%sw_vis_dir(i,j)
1392 if (associated(fluxes%sw_vis_dif)) fluxes%sw_vis_dif(i,j) = frac_open * fluxes%sw_vis_dif(i,j)
1393 if (associated(fluxes%sw_nir_dir)) fluxes%sw_nir_dir(i,j) = frac_open * fluxes%sw_nir_dir(i,j)
1394 if (associated(fluxes%sw_nir_dif)) fluxes%sw_nir_dif(i,j) = frac_open * fluxes%sw_nir_dif(i,j)
1395 if (associated(fluxes%lw)) fluxes%lw(i,j) = frac_open * fluxes%lw(i,j)
1396 if (associated(fluxes%latent)) fluxes%latent(i,j) = frac_open * fluxes%latent(i,j)
1397 if (associated(fluxes%evap)) fluxes%evap(i,j) = frac_open * fluxes%evap(i,j)
1398 if (associated(fluxes%lprec)) then
1399 if (iss%water_flux(i,j) > 0.0) then
1400 fluxes%lprec(i,j) = frac_shelf*iss%water_flux(i,j) + frac_open * fluxes%lprec(i,j)
1401 else
1402 fluxes%lprec(i,j) = frac_open * fluxes%lprec(i,j)
1403 fluxes%evap(i,j) = fluxes%evap(i,j) + frac_shelf*iss%water_flux(i,j)
1404 endif
1405 endif
1406
1407 if (associated(fluxes%sens)) &
1408 fluxes%sens(i,j) = frac_shelf*iss%tflux_ocn(i,j) + frac_open * fluxes%sens(i,j)
1409 ! The salt flux should be mostly from sea ice, so perhaps none should be intercepted and this should be changed.
1410 if (associated(fluxes%salt_flux)) &
1411 fluxes%salt_flux(i,j) = frac_shelf * iss%salt_flux(i,j)*cs%flux_factor + frac_open * fluxes%salt_flux(i,j)
1412 endif ; enddo ; enddo
1413
1414 if (cs%debug) then
1415 call hchksum(iss%water_flux, "water_flux add shelf fluxes", g%HI, haloshift=0, unscale=us%RZ_T_to_kg_m2s)
1416 call hchksum(iss%tflux_ocn, "tflux_ocn add shelf fluxes", g%HI, haloshift=0, unscale=us%QRZ_T_to_W_m2)
1417 call mom_forcing_chksum("After adding shelf fluxes", fluxes, g, cs%US, haloshift=0)
1418 endif
1419
1420 ! Keep sea level constant by removing mass via a balancing flux that might be applied
1421 ! in the open ocean or the sponge region (via virtual precip, vprec). Apply additional
1422 ! salt/heat fluxes so that the resultant surface buoyancy forcing is ~ 0.
1423 ! This is needed for some of the ISOMIP+ experiments.
1424
1425 if (cs%constant_sea_level) then
1426 if (.not. associated(fluxes%salt_flux)) allocate(fluxes%salt_flux(ie,je))
1427 if (.not. associated(fluxes%vprec)) allocate(fluxes%vprec(ie,je))
1428 fluxes%salt_flux(:,:) = 0.0 ; fluxes%vprec(:,:) = 0.0
1429
1430 ! take into account changes in mass (or thickness) when imposing ice shelf mass
1431 if (cs%override_shelf_movement .and. cs%mass_from_file) then
1432 dtime = real_to_time(cs%time_step, unscale=us%T_to_s)
1433
1434 ! Compute changes in mass after at least one full time step
1435 if (cs%Time > dtime) then
1436 time0 = cs%Time - dtime
1437 do j=js,je ; do i=is,ie
1438 last_hmask(i,j) = iss%hmask(i,j) ; last_area_shelf_h(i,j) = iss%area_shelf_h(i,j)
1439 enddo ; enddo
1440 call time_interp_external(cs%mass_handle, time0, last_mass_shelf, scale=us%kg_m3_to_R*us%m_to_Z)
1441 do j=js,je ; do i=is,ie
1442 last_h_shelf(i,j) = last_mass_shelf(i,j) / cs%density_ice
1443 enddo ; enddo
1444
1445 ! apply calving
1446 if (cs%min_thickness_simple_calve > 0.0) then
1447 call ice_shelf_min_thickness_calve(g, last_h_shelf, last_area_shelf_h, last_hmask, &
1448 cs%min_thickness_simple_calve, halo=0)
1449 ! convert to mass again
1450 do j=js,je ; do i=is,ie
1451 last_mass_shelf(i,j) = last_h_shelf(i,j) * cs%density_ice
1452 enddo ; enddo
1453 endif
1454
1455 ! get total ice shelf mass at (Time-dt) and (Time), in kg
1456 do j=js,je ; do i=is,ie
1457 ! Just consider the change in the mass of the floating shelf.
1458 if ((sfc_state%ocean_mass(i,j) > cs%min_ocean_mass_float) .and. &
1459 (iss%area_shelf_h(i,j) > 0.0)) then
1460 delta_float_mass(i,j) = iss%mass_shelf(i,j) - last_mass_shelf(i,j)
1461 else
1462 delta_float_mass(i,j) = 0.0
1463 endif
1464 enddo ; enddo
1465 delta_mass_shelf = global_area_integral(delta_float_mass, g, tmp_scale=us%RZ_to_kg_m2, &
1466 area=iss%area_shelf_h) / cs%time_step
1467 else! first time step
1468 delta_mass_shelf = 0.0
1469 endif
1470 else
1471 if (cs%active_shelf_dynamics) then ! change in ice_shelf draft
1472 do j=js,je ; do i=is,ie
1473 last_h_shelf(i,j) = iss%h_shelf(i,j) - time_step * iss%dhdt_shelf(i,j)
1474 enddo ; enddo
1475 call change_in_draft(cs%dCS, g, last_h_shelf, iss%h_shelf, delta_draft)
1476
1477 !this currently assumes area_shelf_h is constant over the time step
1478 delta_mass_shelf = global_area_integral(delta_draft, g, tmp_scale=us%RZ_to_kg_m2, &
1479 area=iss%area_shelf_h) &
1480 * cs%Rho_ocn / cs%time_step
1481 else ! ice shelf mass does not change
1482 delta_mass_shelf = 0.0
1483 endif
1484 endif
1485
1486 ! average total melt flux over sponge area (ISOMIP/MISOMIP only) or open ocean (general case)
1487 do j=js,je ; do i=is,ie
1488 if (cs%constant_sea_level_misomip) then !for ismip/misomip only
1489 if (g%geoLonT(i,j) >= 790.0) then
1490 bal_frac(i,j) = max(1.0 - iss%area_shelf_h(i,j) * g%IareaT(i,j), 0.0)
1491 else
1492 bal_frac(i,j) = 0.0
1493 endif
1494 elseif ((g%mask2dT(i,j) > 0.0) .and. (iss%area_shelf_h(i,j) * g%IareaT(i,j) < 1.0)) then !general case
1495 bal_frac(i,j) = max(1.0 - iss%area_shelf_h(i,j) * g%IareaT(i,j), 0.0)
1496 else
1497 bal_frac(i,j) = 0.0
1498 endif
1499 enddo ; enddo
1500
1501 balancing_area = global_area_integral(bal_frac, g, area=g%areaT, tmp_scale=1.0)
1502 if (balancing_area > 0.0) then
1503 balancing_flux = ( global_area_integral(iss%water_flux, g, tmp_scale=us%RZ_T_to_kg_m2s, &
1504 area=iss%area_shelf_h) + &
1505 delta_mass_shelf ) / balancing_area
1506 else
1507 balancing_flux = 0.0
1508 endif
1509
1510 ! apply fluxes
1511 do j=js,je ; do i=is,ie
1512 if (bal_frac(i,j) > 0.0) then
1513 ! evap is negative, and vprec has units of [R Z T-1 ~> kg m-2 s-1]
1514 fluxes%vprec(i,j) = -balancing_flux
1515 fluxes%sens(i,j) = fluxes%vprec(i,j) * cs%Cp * cs%T0 ! [Q R Z T-1 ~> W m-2]
1516 fluxes%salt_flux(i,j) = fluxes%vprec(i,j) * cs%S0*1.0e-3*us%S_to_ppt ! [1e-3 S R Z T-1 ~> kgSalt m-2 s-1]
1517 endif
1518 enddo ; enddo
1519
1520 if (cs%debug) then
1521 write(mesg,*) 'Balancing flux (kg/(m^2 s)), dt = ', balancing_flux*us%RZ_T_to_kg_m2s, us%T_to_s*cs%time_step
1522 call mom_mesg(mesg)
1523 call mom_forcing_chksum("After constant sea level", fluxes, g, cs%US, haloshift=0)
1524 endif
1525
1526 endif ! constant_sea_level
1527
1528end subroutine add_shelf_flux
1529
1530
1531!> Initializes shelf model data, parameters and diagnostics
1532subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, Time_init, directory, forces_in, &
1533 fluxes_in, sfc_state_in, solo_ice_sheet_in, calve_ice_shelf_bergs)
1534 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1535 type(ocean_grid_type), pointer :: ocn_grid !< The calling ocean model's horizontal grid structure
1536 type(time_type), intent(inout) :: time !< The clock that that will indicate the model time
1537 type(ice_shelf_cs), pointer :: cs !< A pointer to the ice shelf control structure
1538 type(mom_diag_ctrl), pointer :: diag !< This is a pointer to the MOM diag CS
1539 !! which will be discarded
1540 type(time_type), intent(in) :: time_init !< The time at initialization.
1541 character(len=*), intent(in) :: directory !< The directory where the energy file goes.
1542
1543 type(mech_forcing), optional, target, intent(inout) :: forces_in !< A structure with the driving mechanical forces
1544 type(forcing), optional, target, intent(inout) :: fluxes_in !< A structure containing pointers to any
1545 !! possible thermodynamic or mass-flux forcing fields.
1546 type(surface), target, optional, intent(inout) :: sfc_state_in !< A structure containing fields that
1547 !! describe the surface state of the ocean. The
1548 !! intent is only inout to allow for halo updates.
1549 logical, optional, intent(in) :: solo_ice_sheet_in !< If present, this indicates whether
1550 !! a solo ice-sheet driver.
1551 logical, optional :: calve_ice_shelf_bergs !< If true, will add point iceberg calving variables to the ice
1552 !! shelf restart
1553
1554 type(ocean_grid_type), pointer :: g => null(), og => null() ! Pointers to grids for convenience.
1555 type(unit_scale_type), pointer :: us => null() ! Pointer to a structure containing
1556 ! various unit conversion factors
1557 type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
1558 !! the ice-shelf state
1559 type(directories) :: dirs
1560 type(dyn_horgrid_type), pointer :: dg => null()
1561 type(dyn_horgrid_type), pointer :: dg_in => null()
1562 real :: meltrate_conversion ! The conversion factor to use for in the melt rate diagnostic
1563 ! [T kg R-1 Z-1 m-2 s-1 ~> nondim]
1564 real :: dz_ocean_min_float ! The minimum ocean thickness above which the ice shelf is considered
1565 ! to be floating when CONST_SEA_LEVEL = True [Z ~> m].
1566 real :: cdrag ! The drag coefficient at the ice-ocean interface [nondim]
1567 real :: drag_bg_vel ! A background velocity used in the quadratic drag [Z T-1 ~> m s-1]
1568 logical :: new_sim, save_ic
1569 !This include declares and sets the variable "version".
1570# include "version_variable.h"
1571 character(len=200) :: ic_file, inputdir ! Input file names or paths
1572 character(len=40) :: mdl = "MOM_ice_shelf" ! This module's name.
1573 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, isdq, iedq, jsdq, jedq
1574 integer :: wd_halos(2)
1575 logical :: showcalltree
1576 logical :: read_tideamp, debug
1577 logical :: global_indexing
1578 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
1579 ! recreate the bugs, or if false bugs are only used if actively selected.
1580 character(len=240) :: tideamp_file ! Input file names
1581 character(len=80) :: tideamp_var ! Input file variable names
1582 real :: utide ! A tidal velocity [L T-1 ~> m s-1]
1583 real :: col_thick_melt_thresh ! An ocean column thickness below which iceshelf melting
1584 ! does not occur [Z ~> m]
1585 real, allocatable, dimension(:,:) :: tmp2d ! Temporary array for ice shelf input data [L T-1 ~> m s-1]
1586 real, allocatable, dimension(:,:) :: maskt ! Temporary array for the tracer points masks [nondim]
1587
1588 type(surface), pointer :: sfc_state => null()
1589 type(vardesc) :: u_desc, v_desc
1590
1591 if (associated(cs)) then
1592 call mom_error(fatal, "MOM_ice_shelf.F90, initialize_ice_shelf: "// &
1593 "called with an associated control structure.")
1594 return
1595 endif
1596 allocate(cs)
1597
1598 ! Go through all of the infrastructure initialization calls, since this is
1599 ! being treated as an independent component that just happens to use the
1600 ! MOM's grid and infrastructure.
1601 call get_mom_input(dirs=dirs)
1602
1603 call mom_is_diag_mediator_infrastructure_init()
1604
1605 ! Determining the internal unit scaling factors for this run.
1606 call unit_scaling_init(param_file, cs%US)
1607
1608 call get_param(param_file, mdl, "ROTATE_INDEX", cs%rotate_index, &
1609 "Enable rotation of the horizontal indices.", default=.false., &
1610 debuggingparam=.true.)
1611
1612 call get_param(param_file, "MOM", "GLOBAL_INDEXING", global_indexing, &
1613 "If true, use a global lateral indexing convention, so "//&
1614 "that corresponding points on different processors have "//&
1615 "the same index. This does not work with static memory.", &
1616 default=.false., layoutparam=.true.)
1617
1618 ! Set up the ice-shelf domain and grid
1619 wd_halos(:)=0
1620 allocate(cs%Grid_in)
1621 call mom_domains_init(cs%Grid_in%domain, param_file, min_halo=wd_halos, symmetric=grid_sym_,&
1622 domain_name='MOM_Ice_Shelf_in', us=cs%US)
1623 !allocate(CS%Grid_in%HI)
1624 !call hor_index_init(CS%Grid%Domain, CS%Grid%HI, param_file, &
1625 ! local_indexing=.not.global_indexing)
1626 call mom_grid_init(cs%Grid_in, param_file, cs%US)
1627
1628 if (cs%rotate_index) then
1629 ! ! TODO: Index rotation currently only works when index rotation does not
1630 ! ! change the MPI rank of each domain. Resolving this will require a
1631 ! ! modification to FMS PE assignment.
1632 ! ! For now, we only permit single-core runs.
1633
1634 if (num_pes() /= 1) call mom_error(fatal, "Index rotation is only supported on one PE.")
1635
1636 call get_param(param_file, mdl, "INDEX_TURNS", cs%turns, &
1637 "Number of counterclockwise quarter-turn index rotations.", &
1638 default=1, debuggingparam=.true.)
1639 ! NOTE: If indices are rotated, then CS%Grid and CS%Grid_in must both be initialized.
1640 ! If not rotated, then CS%Grid_in and CS%Ggrid are the same grid.
1641 call create_dyn_horgrid(dg_in, cs%Grid_in%HI)
1642 call clone_mom_domain(cs%Grid_in%Domain, dg_in%Domain)
1643 call set_grid_metrics(dg_in, param_file, cs%US)
1644 ! Set up the bottom depth, dG_in%bathyT, either analytically or from file
1645 call mom_initialize_topography(dg_in%bathyT, cs%Grid_in%max_depth, dg_in, param_file, cs%US)
1646
1647 ! The use of maskT here sets all ice shelf points to be unmasked.
1648 allocate(maskt(dg_in%isd:dg_in%ied,dg_in%jsd:dg_in%jed), source=1.0)
1649 call initialize_masks(dg_in, param_file, cs%US, maskt=maskt)
1650 deallocate(maskt)
1651
1652 call copy_dyngrid_to_mom_grid(dg_in, cs%Grid_in, cs%US)
1653
1654 ! Now set up the rotated ice-shelf grid.
1655 allocate(cs%Grid)
1656 call clone_mom_domain(cs%Grid_in%Domain, cs%Grid%Domain, turns=cs%turns)
1657 call rotate_hor_index(cs%Grid_in%HI, cs%turns, cs%Grid%HI)
1658 call mom_grid_init(cs%Grid, param_file, cs%US, cs%Grid%HI)
1659 call create_dyn_horgrid(dg, cs%Grid%HI)
1660 call rotate_dyngrid(dg_in, dg, cs%US, cs%turns)
1661 call copy_dyngrid_to_mom_grid(dg, cs%Grid, cs%US)
1662
1663 call destroy_dyn_horgrid(dg_in)
1664 call destroy_dyn_horgrid(dg)
1665 else
1666 cs%Grid => cs%Grid_in
1667 dg => null()
1668 call create_dyn_horgrid(dg, cs%Grid%HI)
1669 call clone_mom_domain(cs%Grid%Domain, dg%Domain)
1670 call set_grid_metrics(dg, param_file, cs%US)
1671 ! Set up the bottom depth, dG%bathyT, either analytically or from file
1672 call mom_initialize_topography(dg%bathyT, cs%Grid%max_depth, dg, param_file, cs%US)
1673
1674 ! The use of maskT here sets all ice shelf points to be unmasked.
1675 allocate(maskt(dg%isd:dg%ied,dg%jsd:dg%jed), source=1.0)
1676 call initialize_masks(dg, param_file, cs%US, maskt=maskt)
1677 deallocate(maskt)
1678
1679 call copy_dyngrid_to_mom_grid(dg, cs%Grid, cs%US)
1680 call destroy_dyn_horgrid(dg)
1681 endif
1682 g => cs%Grid
1683
1684 allocate(cs%diag)
1685 call mom_is_diag_mediator_init(g, cs%US, param_file, cs%diag, component='MOM_IceShelf')
1686 ! This call sets up the diagnostic axes. These are needed,
1687 ! e.g. to generate the target grids below.
1688 call set_is_axes_info(g, cs%diag)
1689
1690
1691 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1692 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
1693 isdq = g%IsdB ; iedq = g%IedB ; jsdq = g%JsdB ; jedq = g%JedB
1694
1695 ! The ocean grid possibly uses different symmetry.
1696 if (associated(ocn_grid)) then ; cs%ocn_grid => ocn_grid
1697 else ; cs%ocn_grid => cs%grid ; endif
1698
1699 ! Convenience pointers
1700 og => cs%ocn_grid
1701 us => cs%US
1702
1703 ! Are we being called from the solo ice-sheet driver? When called by the ocean
1704 ! model solo_ice_sheet_in is not preset.
1705 cs%solo_ice_sheet = .false.
1706 if (present(solo_ice_sheet_in)) cs%solo_ice_sheet = solo_ice_sheet_in
1707
1708 !if (present(Time_in)) Time = Time_in
1709
1710
1711 cs%override_shelf_movement = .false. ; cs%active_shelf_dynamics = .false.
1712
1713 call log_version(param_file, mdl, version, "")
1714 call get_param(param_file, mdl, "DEBUG", debug, default=.false.)
1715 call get_param(param_file, mdl, "DEBUG_IS", cs%debug, &
1716 "If true, write verbose debugging messages for the ice shelf.", &
1717 default=debug)
1718 call get_param(param_file, mdl, "DYNAMIC_SHELF_MASS", cs%shelf_mass_is_dynamic, &
1719 "If true, the ice sheet mass can evolve with time.", &
1720 default=.false.)
1721 if (cs%shelf_mass_is_dynamic) then
1722 call get_param(param_file, mdl, "OVERRIDE_SHELF_MOVEMENT", cs%override_shelf_movement, &
1723 "If true, user provided code specifies the ice-shelf "//&
1724 "movement instead of the dynamic ice model.", default=.false.)
1725 cs%active_shelf_dynamics = .not.cs%override_shelf_movement
1726 call get_param(param_file, mdl, "DATA_OVERRIDE_SHELF_FLUXES", &
1727 cs%data_override_shelf_fluxes, &
1728 "If true, the data override feature is used to write "//&
1729 "the surface mass flux deposition. This option is only "//&
1730 "available for MOSAIC grid types.", default=.false.)
1731 call get_param(param_file, mdl, "GROUNDING_LINE_INTERPOLATE", cs%GL_regularize, &
1732 "If true, regularize the floatation condition at the "//&
1733 "grounding line as in Goldberg Holland Schoof 2009.", default=.false.)
1734 call get_param(param_file, mdl, "GROUNDING_LINE_COUPLE", cs%GL_couple, &
1735 "If true, let the floatation condition be determined by "//&
1736 "ocean column thickness. This means that update_OD_ffrac "//&
1737 "will be called. GL_REGULARIZE and GL_COUPLE are exclusive.", &
1738 default=.false., do_not_log=cs%GL_regularize)
1739 if (cs%GL_regularize) cs%GL_couple = .false.
1740 if (cs%solo_ice_sheet) cs%GL_couple = .false.
1741 endif
1742
1743 call get_param(param_file, mdl, "SHELF_THERMO", cs%isthermo, &
1744 "If true, use a thermodynamically interactive ice shelf.", &
1745 default=.false.)
1746 call get_param(param_file, mdl, "LATENT_HEAT_FUSION", cs%Lat_fusion, &
1747 "The latent heat of fusion.", units="J/kg", default=hlf, scale=us%J_kg_to_Q)
1748 call get_param(param_file, mdl, "SHELF_THREE_EQN", cs%threeeq, &
1749 "If true, use the three equation expression of "//&
1750 "consistency to calculate the fluxes at the ice-ocean "//&
1751 "interface.", default=.true.)
1752 call get_param(param_file, mdl, "SHELF_INSULATOR", cs%insulator, &
1753 "If true, the ice shelf is a perfect insulator "//&
1754 "(no conduction).", default=.false.)
1755 call get_param(param_file, mdl, "MELTING_CUTOFF_DEPTH", cs%cutoff_depth, &
1756 "Depth above which the melt is set to zero (it must be >= 0) "//&
1757 "Default value won't affect the solution.", units="m", default=0.0, scale=us%m_to_Z)
1758 if (cs%cutoff_depth < 0.) &
1759 call mom_error(warning,"Initialize_ice_shelf: MELTING_CUTOFF_DEPTH must be >= 0.")
1760
1761 call get_param(param_file, mdl, "CONST_SEA_LEVEL", cs%constant_sea_level, &
1762 "If true, apply evaporative, heat and salt fluxes in "//&
1763 "the sponge region. This will avoid a large increase "//&
1764 "in sea level. This option is needed for some of the "//&
1765 "ISOMIP+ experiments (Ocean3 and Ocean4). "//&
1766 "IMPORTANT: it is not currently possible to do "//&
1767 "prefect restarts using this flag.", default=.false.)
1768 call get_param(param_file, mdl, "CONST_SEA_LEVEL_MISOMIP", cs%constant_sea_level_misomip, &
1769 "If true, constant_sea_level fluxes are applied only over "//&
1770 "the surface sponge cells from the ISOMIP/MISOMIP configuration", default=.false.)
1771 call get_param(param_file, mdl, "MIN_OCEAN_FLOAT_THICK", dz_ocean_min_float, &
1772 "The minimum ocean thickness above which the ice shelf is considered to be "//&
1773 "floating when CONST_SEA_LEVEL = True.", &
1774 default=0.1, units="m", scale=us%m_to_Z, do_not_log=.not.cs%constant_sea_level)
1775
1776 call get_param(param_file, mdl, "ISOMIP_S_SUR_SPONGE", cs%S0, &
1777 "Surface salinity in the restoring region.", &
1778 default=33.8, units='ppt', scale=us%ppt_to_S, do_not_log=.true.)
1779
1780 call get_param(param_file, mdl, "ISOMIP_T_SUR_SPONGE", cs%T0, &
1781 "Surface temperature in the restoring region.", &
1782 default=-1.9, units='degC', scale=us%degC_to_C, do_not_log=.true.)
1783
1784 call get_param(param_file, mdl, "SHELF_3EQ_GAMMA", cs%const_gamma, &
1785 "If true, user specifies a constant nondimensional heat-transfer coefficient "//&
1786 "(GAMMA_T_3EQ), from which the default salt-transfer coefficient is set "//&
1787 "as GAMMA_T_3EQ/35. This is used with SHELF_THREE_EQN.", default=.false.)
1788 call get_param(param_file, mdl, "SHELF_S_ROOT", cs%find_salt_root, &
1789 "If SHELF_S_ROOT = True, salinity at the ice/ocean interface (Sbdry) "//&
1790 "is computed from a quadratic equation. Otherwise, the previous "//&
1791 "interactive method to estimate Sbdry is used.", &
1792 default=.false., do_not_log=.not.cs%threeeq)
1793 if (.not.cs%threeeq) then
1794 call get_param(param_file, mdl, "SHELF_2EQ_GAMMA_T", cs%gamma_t, &
1795 "If SHELF_THREE_EQN is false, this the fixed turbulent "//&
1796 "exchange velocity at the ice-ocean interface.", &
1797 units="m s-1", scale=us%m_to_Z*us%T_to_s, fail_if_missing=.true.)
1798 endif
1799 if (cs%const_gamma .or. cs%find_salt_root) then
1800 call get_param(param_file, mdl, "SHELF_3EQ_GAMMA_T", cs%Gamma_T_3EQ, &
1801 "Nondimensional heat-transfer coefficient.", &
1802 units="nondim", default=2.2e-2)
1803 call get_param(param_file, mdl, "SHELF_3EQ_GAMMA_S", cs%Gamma_S_3EQ, &
1804 "Nondimensional salt-transfer coefficient.", &
1805 default=cs%Gamma_T_3EQ/35.0, units="nondim")
1806 endif
1807
1808 call get_param(param_file, mdl, "ICE_SHELF_MASS_FROM_FILE", &
1809 cs%mass_from_file, "Read the mass of the "//&
1810 "ice shelf (every time step) from a file.", default=.false.)
1811
1812 if (cs%find_salt_root) then ! read liquidus coeffs.
1813 call get_param(param_file, mdl, "TFREEZE_S0_P0", cs%TFr_0_0, &
1814 "this is the freezing potential temperature at "//&
1815 "S=0, P=0.", units="degC", default=0.0, scale=us%degC_to_C, do_not_log=.true.)
1816 call get_param(param_file, mdl, "DTFREEZE_DS", cs%dTFr_dS, &
1817 "this is the derivative of the freezing potential temperature with salinity.", &
1818 units="degC psu-1", default=-0.054, scale=us%degC_to_C*us%S_to_ppt, do_not_log=.true.)
1819 call get_param(param_file, mdl, "DTFREEZE_DP", cs%dTFr_dp, &
1820 "this is the derivative of the freezing potential temperature with pressure.", &
1821 units="degC Pa-1", default=0.0, scale=us%degC_to_C*us%RL2_T2_to_Pa, do_not_log=.true.)
1822 endif
1823
1824 call get_param(param_file, mdl, "G_EARTH", cs%g_Earth, &
1825 "The gravitational acceleration of the Earth.", &
1826 units="m s-2", default=9.80, scale=us%m_s_to_L_T**2*us%Z_to_m)
1827 call get_param(param_file, mdl, "C_P", cs%Cp, &
1828 "The heat capacity of sea water, approximated as a constant. "//&
1829 "The default value is from the TEOS-10 definition of conservative temperature.", &
1830 units="J kg-1 K-1", default=3991.86795711963, scale=us%J_kg_to_Q*us%C_to_degC)
1831 call get_param(param_file, mdl, "RHO_0", cs%Rho_ocn, &
1832 "The mean ocean density used with BOUSSINESQ true to "//&
1833 "calculate accelerations and the mass for conservation "//&
1834 "properties, or with BOUSSINESQ false to convert some "//&
1835 "parameters from vertical units of m to kg m-2.", &
1836 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)
1837 call get_param(param_file, mdl, "C_P_ICE", cs%Cp_ice, &
1838 "The heat capacity of ice.", units="J kg-1 K-1", scale=us%J_kg_to_Q*us%C_to_degC, &
1839 default=2.10e3)
1840 if (cs%constant_sea_level) cs%min_ocean_mass_float = dz_ocean_min_float*cs%Rho_ocn
1841
1842 call get_param(param_file, mdl, "ICE_SHELF_FLUX_FACTOR", cs%flux_factor, &
1843 "Non-dimensional factor applied to shelf thermodynamic fluxes.", &
1844 units="none", default=1.0)
1845
1846 call get_param(param_file, mdl, "KV_ICE", cs%kv_ice, &
1847 "The viscosity of the ice.", &
1848 units="m2 s-1", default=1.0e10, scale=us%Z_to_L**2*us%m_to_L**2*us%T_to_s)
1849 call get_param(param_file, mdl, "KV_MOLECULAR", cs%kv_molec, &
1850 "The molecular kinematic viscosity of sea water at the freezing temperature.", &
1851 units="m2 s-1", default=1.95e-6, scale=us%m2_s_to_Z2_T)
1852 call get_param(param_file, mdl, "ICE_SHELF_SALINITY", cs%Salin_ice, &
1853 "The salinity of the ice inside the ice shelf.", &
1854 units="psu", default=0.0, scale=us%ppt_to_S)
1855 call get_param(param_file, mdl, "ICE_SHELF_TEMPERATURE", cs%Temp_ice, &
1856 "The temperature at the center of the ice shelf.", &
1857 units="degC", default=-15.0, scale=us%degC_to_C)
1858 call get_param(param_file, mdl, "KD_SALT_MOLECULAR", cs%kd_molec_salt, &
1859 "The molecular diffusivity of salt in sea water at the "//&
1860 "freezing point.", units="m2 s-1", default=8.02e-10, scale=us%m2_s_to_Z2_T)
1861 call get_param(param_file, mdl, "KD_TEMP_MOLECULAR", cs%kd_molec_temp, &
1862 "The molecular diffusivity of heat in sea water at the "//&
1863 "freezing point.", units="m2 s-1", default=1.41e-7, scale=us%m2_s_to_Z2_T)
1864 call get_param(param_file, mdl, "DT_FORCING", cs%time_step, &
1865 "The time step for changing forcing, coupling with other "//&
1866 "components, or potentially writing certain diagnostics. "//&
1867 "The default value is given by DT.", units="s", default=0.0, scale=us%s_to_T)
1868
1869 call get_param(param_file, mdl, "COL_THICK_MELT_THRESHOLD", col_thick_melt_thresh, &
1870 "The minimum ocean column thickness where melting is allowed.", &
1871 units="m", scale=us%m_to_Z, default=0.0)
1872 cs%col_mass_melt_threshold = cs%Rho_ocn * col_thick_melt_thresh
1873
1874 call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, &
1875 "If true, read a file (given by TIDEAMP_FILE) containing "//&
1876 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
1877 call get_param(param_file, mdl, "ICE_SHELF_LINEAR_SHELF_FRAC", cs%Zeta_N, &
1878 "Ratio of HJ99 stability constant xi_N (ratio of maximum "//&
1879 "mixing length to planetary boundary layer depth in "//&
1880 "neutrally stable conditions) to the von Karman constant", &
1881 units="nondim", default=0.13)
1882 call get_param(param_file, mdl, "ICE_SHELF_VK_CNST", cs%Vk, &
1883 "Von Karman constant.", &
1884 units="nondim", default=0.40)
1885 call get_param(param_file, mdl, "ICE_SHELF_RC", cs%Rc, &
1886 "Critical flux Richardson number for ice melt ", &
1887 units="nondim", default=0.20)
1888 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
1889 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
1890 call get_param(param_file, mdl, "ICE_SHELF_USTAR_FROM_VEL_BUGFIX", cs%ustar_from_vel_bugfix, &
1891 "Bug fix for ice-area weighting of squared ocean velocities "//&
1892 "used to calculate friction velocity under ice shelves", default=.not.enable_bugs)
1893 call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_BUGFIX", cs%buoy_flux_itt_bugfix, &
1894 "Bug fix of buoyancy iteration", default=.true., old_name="ICE_SHELF_BUOYANCY_FLUX_ITT_BUG")
1895 call get_param(param_file, mdl, "ICE_SHELF_SALT_FLUX_ITT_BUGFIX", cs%salt_flux_itt_bugfix, &
1896 "Bug fix of salt iteration", default=.true., old_name="ICE_SHELF_SALT_FLUX_ITT_BUG")
1897 call get_param(param_file, mdl, "ICE_SHELF_BUOYANCY_FLUX_ITT_THRESHOLD", cs%buoy_flux_tol, &
1898 "Convergence criterion of Newton's method for ice shelf "//&
1899 "buoyancy iteration.", units="nondim", default=1.0e-4)
1900
1901 if (PRESENT(sfc_state_in)) then
1902 ! assuming frazil is enabled in ocean. This could break some configurations?
1903 call allocate_surface_state(sfc_state_in, cs%Grid_in, use_temperature=.true., &
1904 do_integrals=.true., omit_frazil=.false., use_iceshelves=.true.)
1905 if (cs%rotate_index) then
1906 allocate(sfc_state)
1907 call rotate_surface_state(sfc_state_in, sfc_state, cs%Grid, cs%turns)
1908 else
1909 sfc_state => sfc_state_in
1910 endif
1911 endif
1912
1913
1914 call safe_alloc_ptr(cs%utide,isd,ied,jsd,jed) ; cs%utide(:,:) = 0.0
1915
1916 if (read_tideamp) then
1917 call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
1918 "The path to the file containing the spatially varying tidal amplitudes.", &
1919 default="tideamp.nc")
1920 call get_param(param_file, mdl, "TIDEAMP_VARNAME", tideamp_var, &
1921 "The name of the tidal amplitude variable in the input file.", &
1922 default="tideamp")
1923 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1924 inputdir = slasher(inputdir)
1925 tideamp_file = trim(inputdir) // trim(tideamp_file)
1926 if (cs%rotate_index) then
1927 allocate(tmp2d(cs%Grid_in%isd:cs%Grid_in%ied,cs%Grid_in%jsd:cs%Grid_in%jed), source=0.0)
1928 call mom_read_data(tideamp_file, tideamp_var, tmp2d, cs%Grid_in%domain, timelevel=1, scale=us%m_s_to_L_T)
1929 call rotate_array(tmp2d, cs%turns, cs%utide)
1930 deallocate(tmp2d)
1931 else
1932 call mom_read_data(tideamp_file, tideamp_var, cs%utide, cs%Grid%domain, timelevel=1, scale=us%m_s_to_L_T)
1933 endif
1934 else
1935 call get_param(param_file, mdl, "UTIDE", utide, &
1936 "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
1937 units="m s-1", default=0.0 , scale=us%m_s_to_L_T)
1938 cs%utide(:,:) = utide
1939 endif
1940
1941 call eos_init(param_file, cs%eqn_of_state, us)
1942
1943 !! new parameters that need to be in MOM_input
1944
1945 if (cs%active_shelf_dynamics) then
1946
1947 call get_param(param_file, mdl, "DENSITY_ICE", cs%density_ice, &
1948 "A typical density of ice.", units="kg m-3", default=917.0, scale=us%kg_m3_to_R)
1949
1950 call get_param(param_file, mdl, "INPUT_FLUX_ICE_SHELF", cs%input_flux, &
1951 "volume flux at upstream boundary", units="m2 s-1", default=0., scale=us%m_to_Z*us%m_s_to_L_T)
1952 call get_param(param_file, mdl, "INPUT_THICK_ICE_SHELF", cs%input_thickness, &
1953 "flux thickness at upstream boundary", units="m", default=1000., scale=us%m_to_Z)
1954 else
1955 ! This is here because of inconsistent defaults. I don't know why. RWH
1956 call get_param(param_file, mdl, "DENSITY_ICE", cs%density_ice, &
1957 "A typical density of ice.", units="kg m-3", default=900.0, scale=us%kg_m3_to_R)
1958 endif
1959 call get_param(param_file, mdl, "MIN_THICKNESS_SIMPLE_CALVE", &
1960 cs%min_thickness_simple_calve, &
1961 "Min thickness rule for the very simple calving law",&
1962 units="m", default=0.0, scale=us%m_to_Z)
1963
1964 call get_param(param_file, mdl, "USTAR_SHELF_BG", cs%ustar_bg, &
1965 "The minimum value of ustar under ice shelves.", &
1966 units="m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
1967 call get_param(param_file, mdl, "CDRAG_SHELF", cdrag, &
1968 "CDRAG is the drag coefficient relating the magnitude of "//&
1969 "the velocity field to the surface stress.", units="nondim", &
1970 default=0.003)
1971 cs%cdrag = cdrag
1972 if (cs%ustar_bg <= 0.0) then
1973 call get_param(param_file, mdl, "DRAG_BG_VEL_SHELF", drag_bg_vel, &
1974 "DRAG_BG_VEL is either the assumed bottom velocity (with "//&
1975 "LINEAR_DRAG) or an unresolved velocity that is "//&
1976 "combined with the resolved velocity to estimate the "//&
1977 "velocity magnitude.", units="m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
1978 if (cs%cdrag*drag_bg_vel > 0.0) cs%ustar_bg = sqrt(cs%cdrag)*drag_bg_vel
1979 endif
1980 call get_param(param_file, mdl, "USTAR_SHELF_FROM_VEL", cs%ustar_shelf_from_vel, &
1981 "If true, use the surface velocities to set the friction velocity under ice "//&
1982 "shelves instead of using the previous values of the stresses.", &
1983 default=.true.)
1984 call get_param(param_file, mdl, "USTAR_SHELF_MAX", cs%ustar_max, &
1985 "The maximum value of ustar under ice shelves, or a negative value for no limit.", &
1986 units="m s-1", default=-1.0, scale=us%m_to_Z*us%T_to_s, &
1987 do_not_log=cs%ustar_shelf_from_vel)
1988
1989 if (present(calve_ice_shelf_bergs)) cs%calve_ice_shelf_bergs=calve_ice_shelf_bergs
1990
1991 ! Allocate and initialize state variables to default values
1992 call ice_shelf_state_init(cs%ISS, cs%grid)
1993 iss => cs%ISS
1994
1995 new_sim = .false.
1996 if ((dirs%input_filename(1:1) == 'n') .and. &
1997 (len_trim(dirs%input_filename) == 1)) new_sim = .true.
1998
1999 iss%area_shelf_h(:,:)=0.0
2000 iss%h_shelf(:,:)=0.0
2001 iss%hmask(:,:)=0.0
2002 iss%mass_shelf(:,:)=0.0
2003
2004 if (cs%override_shelf_movement .and. cs%mass_from_file) then
2005
2006 ! initialize the ids for reading shelf mass from a netCDF
2007 call initialize_shelf_mass(g, param_file, cs, iss)
2008
2009 if (new_sim) then
2010 ! new simulation, initialize ice thickness as in the static case
2011 call initialize_ice_thickness(iss%h_shelf, iss%area_shelf_h, iss%hmask, iss%melt_mask, cs%Grid, cs%Grid_in, &
2012 us, param_file, cs%rotate_index, cs%turns)
2013
2014 ! next make sure mass is consistent with thickness
2015 do j=g%jsd,g%jed ; do i=g%isd,g%ied
2016 if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2) .or. (iss%hmask(i,j)==3)) then
2017 iss%mass_shelf(i,j) = iss%h_shelf(i,j)*cs%density_ice
2018 endif
2019 enddo ; enddo
2020
2021 if (cs%min_thickness_simple_calve > 0.0) &
2022 call ice_shelf_min_thickness_calve(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, &
2023 cs%min_thickness_simple_calve)
2024 endif
2025 endif
2026
2027 if (cs%active_shelf_dynamics) then
2028 ! the only reason to initialize boundary conds is if the shelf is dynamic - MJH
2029
2030 ! call initialize_ice_shelf_boundary ( CS%u_face_mask_bdry, CS%v_face_mask_bdry, &
2031 ! CS%u_flux_bdry_val, CS%v_flux_bdry_val, &
2032 ! CS%u_bdry_val, CS%v_bdry_val, CS%h_bdry_val, &
2033 ! ISS%hmask, G, param_file)
2034
2035 endif
2036
2037 ! Set up the restarts.
2038
2039 call restart_init(param_file, cs%restart_CSp, "Shelf.res")
2040 call register_restart_field(iss%mass_shelf, "shelf_mass", .true., cs%restart_CSp, &
2041 "Ice shelf mass", "kg m-2", conversion=us%RZ_to_kg_m2)
2042 call register_restart_field(iss%area_shelf_h, "shelf_area", .true., cs%restart_CSp, &
2043 "Ice shelf area in cell", "m2", conversion=us%L_to_m**2)
2044 call register_restart_field(iss%h_shelf, "h_shelf", .true., cs%restart_CSp, &
2045 "ice sheet/shelf thickness", "m", conversion=us%Z_to_m)
2046 call register_restart_field(iss%melt_mask, "melt_mask", .false., cs%restart_CSp, &
2047 "Mask that is >0 where ice-shelf melting is allowed", "none")
2048 if (cs%calve_ice_shelf_bergs) then
2049 call register_restart_field(iss%calving, "shelf_calving", .true., cs%restart_CSp, &
2050 "Calving flux from ice shelf into icebergs", "kg m-2", conversion=us%RZ_to_kg_m2)
2051 call register_restart_field(iss%calving_hflx, "shelf_calving_hflx", .true., cs%restart_CSp, &
2052 "Calving heat flux from ice shelf into icebergs", "W m-2", conversion=us%QRZ_T_to_W_m2)
2053 endif
2054
2055 if (PRESENT(sfc_state_in)) then
2056 if (allocated(sfc_state%taux_shelf) .and. allocated(sfc_state%tauy_shelf)) then
2057 u_desc = var_desc("taux_shelf", "Pa", "the zonal stress on the ocean under ice shelves", &
2058 hor_grid='Cu',z_grid='1')
2059 v_desc = var_desc("tauy_shelf", "Pa", "the meridional stress on the ocean under ice shelves", &
2060 hor_grid='Cv',z_grid='1')
2061 call register_restart_pair(sfc_state%taux_shelf, sfc_state%tauy_shelf, u_desc, v_desc, &
2062 .false., cs%restart_CSp, conversion=us%RLZ_T2_to_Pa)
2063 endif
2064 endif
2065
2066 if (cs%active_shelf_dynamics) then
2067 call register_restart_field(iss%hmask, "h_mask", .true., cs%restart_CSp, &
2068 "ice sheet/shelf thickness mask" ,"none")
2069 endif
2070
2071 if (cs%active_shelf_dynamics) then
2072 ! Allocate CS%dCS and specify additional restarts for ice shelf dynamics
2073 call register_ice_shelf_dyn_restarts(cs%Grid_in, us, param_file, cs%dCS, cs%restart_CSp)
2074 endif
2075
2076 !GMM - I think we do not need to save ustar_shelf and iceshelf_melt in the restart file
2077 !if (.not. CS%solo_ice_sheet) then
2078 ! call register_restart_field(fluxes%ustar_shelf, "ustar_shelf", .false., CS%restart_CSp, &
2079 ! "Friction velocity under ice shelves", "m s-1", conversion=US%Z_to_m*US%s_to_T)
2080 !endif
2081
2082 cs%restart_output_dir = dirs%restart_output_dir
2083
2084 if (present(fluxes_in)) then
2085 call initialize_ice_shelf_fluxes(cs, ocn_grid, us, fluxes_in)
2086 call register_restart_field(fluxes_in%shelf_sfc_mass_flux, "sfc_mass_flux", .true., cs%restart_CSp, &
2087 "ice shelf surface mass flux deposition from atmosphere", &
2088 'kg m-2 s-1', conversion=us%RZ_T_to_kg_m2s)
2089 endif
2090
2091 if (new_sim .and. (.not. (cs%override_shelf_movement .and. cs%mass_from_file))) then
2092 ! This model is initialized internally or from a file.
2093 call initialize_ice_thickness(iss%h_shelf, iss%area_shelf_h, iss%hmask, iss%melt_mask, cs%Grid, cs%Grid_in, &
2094 us, param_file, cs%rotate_index, cs%turns)
2095 ! next make sure mass is consistent with thickness
2096 do j=g%jsd,g%jed ; do i=g%isd,g%ied
2097 if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2) .or. (iss%hmask(i,j) == 3)) then
2098 iss%mass_shelf(i,j) = iss%h_shelf(i,j)*cs%density_ice
2099 endif
2100 enddo ; enddo
2101 if (cs%debug) then
2102 call hchksum(iss%mass_shelf, "IS init: mass_shelf", g%HI, haloshift=0, unscale=us%RZ_to_kg_m2)
2103 call hchksum(iss%area_shelf_h, "IS init: area_shelf", g%HI, haloshift=0, unscale=us%L_to_m*us%L_to_m)
2104 call hchksum(iss%hmask, "IS init: hmask", g%HI, haloshift=0)
2105 endif
2106
2107 ! else ! Previous block for new_sim=.T., this block restores the state.
2108 elseif (.not.new_sim) then
2109 ! This line calls a subroutine that reads the initial conditions from a restart file.
2110 call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: Restoring ice shelf from file.")
2111 call restore_state(dirs%input_filename, dirs%restart_input_dir, time, g, cs%restart_CSp)
2112
2113 endif ! .not. new_sim
2114
2115! do j=G%jsc,G%jec ; do i=G%isc,G%iec
2116! ISS%area_shelf_h(i,j) = ISS%area_shelf_h(i,j)*G%mask2dT(i,j)
2117! enddo ; enddo
2118
2119 id_clock_shelf = cpu_clock_id('Ice shelf', grain=clock_component)
2120 id_clock_pass = cpu_clock_id(' Ice shelf halo updates', grain=clock_routine)
2121
2122 call cpu_clock_begin(id_clock_pass)
2123 call pass_var(iss%area_shelf_h, g%domain, complete=.false.)
2124 call pass_var(iss%h_shelf, g%domain, complete=.false.)
2125 call pass_var(iss%mass_shelf, g%domain, complete=.false.)
2126 call pass_var(iss%hmask, g%domain, complete=.false.)
2127 call pass_var(g%bathyT, g%domain)
2128 call cpu_clock_end(id_clock_pass)
2129
2130 do j=jsd,jed ; do i=isd,ied
2131 if (iss%area_shelf_h(i,j) > g%areaT(i,j)) then
2132 call mom_error(warning,"Initialize_ice_shelf: area_shelf_h exceeds G%areaT.")
2133 iss%area_shelf_h(i,j) = g%areaT(i,j)
2134 endif
2135 enddo ; enddo
2136
2137 if (cs%debug) then
2138 call hchksum(iss%area_shelf_h, "IS init: area_shelf_h", g%HI, haloshift=0, unscale=us%L_to_m*us%L_to_m)
2139 endif
2140
2141 cs%Time = time
2142
2143 if (cs%active_shelf_dynamics .and. .not.cs%isthermo) then
2144 iss%water_flux(:,:) = 0.0
2145 endif
2146
2147 if (cs%shelf_mass_is_dynamic) &
2148 call initialize_ice_shelf_dyn(param_file, time, iss, cs%dCS, g, us, cs%diag, new_sim, cs%Cp_ice, &
2149 time_init, directory, solo_ice_sheet_in)
2150
2151 call fix_restart_unit_scaling(us, unscaled=.true.)
2152
2153 call get_param(param_file, mdl, "SAVE_INITIAL_CONDS", save_ic, &
2154 "If true, save the ice shelf initial conditions.", default=.false.)
2155 if (save_ic) call get_param(param_file, mdl, "SHELF_IC_OUTPUT_FILE", ic_file,&
2156 "The name-root of the output file for the ice shelf initial conditions.", &
2157 default="MOM_Shelf_IC")
2158
2159 if (save_ic .and. .not.((dirs%input_filename(1:1) == 'r') .and. &
2160 (len_trim(dirs%input_filename) == 1))) then
2161 showcalltree = calltree_showquery()
2162 if (showcalltree) call calltree_waypoint("About to call save_restart (MOM_ice_shelf)")
2163 call save_restart(dirs%output_directory, cs%Time, cs%Grid_in, cs%restart_CSp, &
2164 filename=ic_file, write_ic=.true.)
2165 if (showcalltree) call calltree_waypoint("Done with call to save_restart (MOM_ice_shelf)")
2166 endif
2167
2168 cs%id_area_shelf_h = register_diag_field('ice_shelf_model', 'area_shelf_h', cs%diag%axesT1, cs%Time, &
2169 'Ice Shelf Area in cell', 'meter2', conversion=us%L_to_m**2)
2170 cs%id_shelf_mass = register_diag_field('ice_shelf_model', 'shelf_mass', cs%diag%axesT1, cs%Time, &
2171 'mass of shelf', 'kg/m^2', conversion=us%RZ_to_kg_m2)
2172 cs%id_h_shelf = register_diag_field('ice_shelf_model', 'h_shelf', cs%diag%axesT1, cs%Time, &
2173 'ice shelf thickness', 'm', conversion=us%Z_to_m)
2174 cs%id_dhdt_shelf = register_diag_field('ice_shelf_model', 'dhdt_shelf', cs%diag%axesT1, cs%Time, &
2175 'change in ice shelf thickness over time', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
2176 cs%id_mass_flux = register_diag_field('ice_shelf_model', 'mass_flux', cs%diag%axesT1,&
2177 cs%Time, 'Total mass flux of freshwater across the ice-ocean interface.', &
2178 'kg/s', conversion=us%RZ_T_to_kg_m2s*us%L_to_m**2)
2179
2180 if (cs%const_gamma) then ! use ISOMIP+ eq. with rho_fw = 1000. kg m-3
2181 meltrate_conversion = 86400.0*365.0*us%Z_to_m*us%s_to_T / (1000.0*us%kg_m3_to_R)
2182 else ! use original eq.
2183 meltrate_conversion = 86400.0*365.0*us%Z_to_m*us%s_to_T / cs%density_ice
2184 endif
2185 cs%id_melt = register_diag_field('ice_shelf_model', 'melt', cs%diag%axesT1, cs%Time, &
2186 'Ice Shelf Melt Rate', 'm yr-1', conversion=meltrate_conversion)
2187 cs%id_thermal_driving = register_diag_field('ice_shelf_model', 'thermal_driving', cs%diag%axesT1, cs%Time, &
2188 'pot. temp. in the boundary layer minus freezing pot. temp. at the ice-ocean interface.', &
2189 'Celsius', conversion=us%C_to_degC)
2190 cs%id_haline_driving = register_diag_field('ice_shelf_model', 'haline_driving', cs%diag%axesT1, cs%Time, &
2191 'salinity in the boundary layer minus salinity at the ice-ocean interface.', &
2192 'psu', conversion=us%S_to_ppt)
2193 cs%id_Sbdry = register_diag_field('ice_shelf_model', 'sbdry', cs%diag%axesT1, cs%Time, &
2194 'salinity at the ice-ocean interface.', 'psu', conversion=us%S_to_ppt)
2195 cs%id_u_ml = register_diag_field('ice_shelf_model', 'u_ml', cs%diag%axesCu1, cs%Time, &
2196 'Eastward vel. in the boundary layer (used to compute ustar)', 'm s-1', conversion=us%L_T_to_m_s)
2197 cs%id_v_ml = register_diag_field('ice_shelf_model', 'v_ml', cs%diag%axesCv1, cs%Time, &
2198 'Northward vel. in the boundary layer (used to compute ustar)', 'm s-1', conversion=us%L_T_to_m_s)
2199 cs%id_exch_vel_s = register_diag_field('ice_shelf_model', 'exch_vel_s', cs%diag%axesT1, cs%Time, &
2200 'Sub-shelf salinity exchange velocity', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
2201 cs%id_exch_vel_t = register_diag_field('ice_shelf_model', 'exch_vel_t', cs%diag%axesT1, cs%Time, &
2202 'Sub-shelf thermal exchange velocity', 'm s-1' , conversion=us%Z_to_m*us%s_to_T)
2203 cs%id_tfreeze = register_diag_field('ice_shelf_model', 'tfreeze', cs%diag%axesT1, cs%Time, &
2204 'In Situ Freezing point at ice shelf interface', 'degC', conversion=us%C_to_degC)
2205 cs%id_tfl_shelf = register_diag_field('ice_shelf_model', 'tflux_shelf', cs%diag%axesT1, cs%Time, &
2206 'Heat conduction into ice shelf', 'W m-2', conversion=-us%QRZ_T_to_W_m2)
2207 cs%id_ustar_shelf = register_diag_field('ice_shelf_model', 'ustar_shelf', cs%diag%axesT1, cs%Time, &
2208 'Fric vel under shelf', 'm/s', conversion=us%Z_to_m*us%s_to_T)
2209 cs%id_frazil = register_diag_field('ice_shelf_model', 'frazil', cs%diag%axesT1, cs%Time, &
2210 'Frazil heat rejected by the ocean', 'J m-2', conversion=us%Q_to_J_kg*us%RZ_to_kg_m2)
2211 if (cs%active_shelf_dynamics) then
2212 cs%id_h_mask = register_diag_field('ice_shelf_model', 'h_mask', cs%diag%axesT1, cs%Time, &
2213 'ice shelf thickness mask', 'none', conversion=1.0)
2214 endif
2215
2216 cs%id_shelf_sfc_mass_flux = register_diag_field('ice_shelf_model', 'sfc_mass_flux', cs%diag%axesT1, cs%Time, &
2217 'ice shelf surface mass flux deposition from atmosphere', &
2218 'kg m-2 s-1', conversion=us%RZ_T_to_kg_m2s)
2219
2220 ! Scalars (area integrated over all ice sheets)
2221 cs%id_vaf = register_scalar_field('ice_shelf_model', 'int_vaf', cs%diag%axesT1, cs%Time, &
2222 'Area integrated ice sheet volume above floatation', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2223 cs%id_adott = register_scalar_field('ice_shelf_model', 'int_a', cs%diag%axesT1, cs%Time, &
2224 'Area integrated change in ice-sheet thickness ' //&
2225 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2226 cs%id_g_adott = register_scalar_field('ice_shelf_model', 'int_a_ground', cs%diag%axesT1, cs%Time, &
2227 'Area integrated change in grounded ice-sheet thickness ' //&
2228 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2229 cs%id_f_adott = register_scalar_field('ice_shelf_model', 'int_a_float', cs%diag%axesT1, cs%Time, &
2230 'Area integrated change in floating ice-shelf thickness ' //&
2231 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2232 cs%id_bdott = register_scalar_field('ice_shelf_model', 'int_b', cs%diag%axesT1, cs%Time, &
2233 'Area integrated change in floating ice-shelf thickness '//&
2234 'due to basal accum+melt during a DT_THERM time step', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2235 cs%id_bdott_melt = register_scalar_field('ice_shelf_model', 'int_b_melt', cs%diag%axesT1, cs%Time, &
2236 'Area integrated basal melt over ice shelves during a DT_THERM time step', &
2237 units='m3', conversion=us%Z_to_m*us%L_to_m**2)
2238 cs%id_bdott_accum = register_scalar_field('ice_shelf_model', 'int_b_accum', cs%diag%axesT1, cs%Time, &
2239 'Area integrated basal accumulation over ice shelves during a DT_THERM a time step', &
2240 units='m3', conversion=us%Z_to_m*us%L_to_m**2)
2241 cs%id_t_area = register_scalar_field('ice_shelf_model', 'tot_area', cs%diag%axesT1, cs%Time, &
2242 'Total ice-sheet area', 'm2', conversion=us%L_to_m**2)
2243 cs%id_f_area = register_scalar_field('ice_shelf_model', 'tot_area_float', cs%diag%axesT1, cs%Time, &
2244 'Total area of floating ice shelves', 'm2', conversion=us%L_to_m**2)
2245 cs%id_g_area = register_scalar_field('ice_shelf_model', 'tot_area_ground', cs%diag%axesT1, cs%Time, &
2246 'Total area of grounded ice sheets', 'm2', conversion=us%L_to_m**2)
2247 !scalars (area integrated rates over all ice sheets)
2248 cs%id_dvafdt = register_scalar_field('ice_shelf_model', 'int_vafdot', cs%diag%axesT1, cs%Time, &
2249 'Area integrated rate of change in ice-sheet volume above floatation', &
2250 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2251 cs%id_adot = register_scalar_field('ice_shelf_model', 'int_adot', cs%diag%axesT1, cs%Time, &
2252 'Area integrated rate of change in ice-sheet thickness due to surface accum+melt', &
2253 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2254 cs%id_g_adot = register_scalar_field('ice_shelf_model', 'int_adot_ground', cs%diag%axesT1, cs%Time, &
2255 'Area integrated rate of change in grounded ice-sheet thickness due to surface accum+melt', &
2256 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2257 cs%id_f_adot = register_scalar_field('ice_shelf_model', 'int_adot_float', cs%diag%axesT1, cs%Time, &
2258 'Area integrated rate of change in floating ice-shelf thickness due to surface accum+melt', &
2259 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2260 cs%id_bdot = register_scalar_field('ice_shelf_model', 'int_bdot', cs%diag%axesT1, cs%Time, &
2261 'Area integrated rate of change in ice-shelf thickness due to basal accum+melt', &
2262 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2263 cs%id_bdot_melt = register_scalar_field('ice_shelf_model', 'int_bdot_melt', cs%diag%axesT1, cs%Time, &
2264 'Area integrated basal melt rate over ice shelves', &
2265 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2266 cs%id_bdot_accum = register_scalar_field('ice_shelf_model', 'int_bdot_accum', cs%diag%axesT1, cs%Time, &
2267 'Area integrated basal accumulation rate over ice shelves', &
2268 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2269
2270 !scalars (area integrated over the Antarctic ice sheet)
2271 cs%id_Ant_vaf = register_scalar_field('ice_shelf_model', 'int_vaf_A', cs%diag%axesT1, cs%Time, &
2272 'Area integrated Antarctic ice sheet volume above floatation', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2273 cs%id_Ant_adott = register_scalar_field('ice_shelf_model', 'int_a_A', cs%diag%axesT1, cs%Time, &
2274 'Area integrated (Antarctic ice sheet) change in ice-sheet thickness ' //&
2275 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2276 cs%id_Ant_g_adott = register_scalar_field('ice_shelf_model', 'int_a_ground_A', cs%diag%axesT1, cs%Time, &
2277 'Area integrated change in Antarctic grounded ice-sheet thickness ' //&
2278 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2279 cs%id_Ant_f_adott = register_scalar_field('ice_shelf_model', 'int_a_float_A', cs%diag%axesT1, cs%Time, &
2280 'Area integrated change in Antarctic floating ice-shelf thickness ' //&
2281 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2282 cs%id_Ant_bdott = register_scalar_field('ice_shelf_model', 'int_b_A', cs%diag%axesT1, cs%Time, &
2283 'Area integrated change in Antarctic floating ice-shelf thickness '//&
2284 'due to basal accum+melt during a DT_THERM time step', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2285 cs%id_Ant_bdott_melt = register_scalar_field('ice_shelf_model', 'int_b_melt_A', cs%diag%axesT1, cs%Time, &
2286 'Area integrated basal melt over Antarctic ice shelves during a DT_THERM time step', &
2287 units='m3', conversion=us%Z_to_m*us%L_to_m**2)
2288 cs%id_Ant_bdott_accum = register_scalar_field('ice_shelf_model', 'int_b_accum_A', cs%diag%axesT1, cs%Time, &
2289 'Area integrated basal accumulation over Antarctic ice shelves during a DT_THERM a time step', &
2290 units='m3', conversion=us%Z_to_m*us%L_to_m**2)
2291 cs%id_Ant_t_area = register_scalar_field('ice_shelf_model', 'tot_area_A', cs%diag%axesT1, cs%Time, &
2292 'Total area of Antarctic ice sheet', 'm2', conversion=us%L_to_m**2)
2293 cs%id_Ant_f_area = register_scalar_field('ice_shelf_model', 'tot_area_float_A', cs%diag%axesT1, cs%Time, &
2294 'Total area of Antarctic floating ice shelves', 'm2', conversion=us%L_to_m**2)
2295 cs%id_Ant_g_area = register_scalar_field('ice_shelf_model', 'tot_area_ground_A', cs%diag%axesT1, cs%Time, &
2296 'Total area of Antarctic grounded ice sheet', 'm2', conversion=us%L_to_m**2)
2297 !scalars (area integrated rates over the Antarctic ice sheet)
2298 cs%id_Ant_dvafdt = register_scalar_field('ice_shelf_model', 'int_vafdot_A', cs%diag%axesT1, cs%Time, &
2299 'Area integrated rate of change in Antarctic ice-sheet volume above floatation', &
2300 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2301 cs%id_Ant_adot = register_scalar_field('ice_shelf_model', 'int_adot_A', cs%diag%axesT1, cs%Time, &
2302 'Area integrated rate of change in Antarctic ice-sheet thickness due to surface accum+melt', &
2303 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2304 cs%id_Ant_g_adot = register_scalar_field('ice_shelf_model', 'int_adot_ground_A', cs%diag%axesT1, cs%Time, &
2305 'Area integrated rate of change in Antarctic grounded ice-sheet thickness due to surface accum+melt', &
2306 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2307 cs%id_Ant_f_adot = register_scalar_field('ice_shelf_model', 'int_adot_float_A', cs%diag%axesT1, cs%Time, &
2308 'Area integrated rate of change in Antarctic floating ice-shelf thickness due to surface accum+melt', &
2309 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2310 cs%id_Ant_bdot = register_scalar_field('ice_shelf_model', 'int_bdot_A', cs%diag%axesT1, cs%Time, &
2311 'Area integrated rate of change in Antarctic ice-shelf thickness due to basal accum+melt', &
2312 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2313 cs%id_Ant_bdot_melt = register_scalar_field('ice_shelf_model', 'int_bdot_melt_A', cs%diag%axesT1, cs%Time, &
2314 'Area integrated basal melt rate over Antarctic ice shelves', &
2315 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2316 cs%id_Ant_bdot_accum = register_scalar_field('ice_shelf_model', 'int_bdot_accum_A', cs%diag%axesT1, cs%Time, &
2317 'Area integrated basal accumulation rate over Antarctic ice shelves', &
2318 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2319
2320 !scalars (area integrated over the Greenland ice sheet)
2321 cs%id_Gr_vaf = register_scalar_field('ice_shelf_model', 'int_vaf_G', cs%diag%axesT1, cs%Time, &
2322 'Area integrated Greenland ice sheet volume above floatation', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2323 cs%id_Gr_adott = register_scalar_field('ice_shelf_model', 'int_a_G', cs%diag%axesT1, cs%Time, &
2324 'Area integrated (Greenland ice sheet) change in ice-sheet thickness ' //&
2325 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2326 cs%id_Gr_g_adott = register_scalar_field('ice_shelf_model', 'int_a_ground_G', cs%diag%axesT1, cs%Time, &
2327 'Area integrated change in Greenland grounded ice-sheet thickness ' //&
2328 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2329 cs%id_Gr_f_adott = register_scalar_field('ice_shelf_model', 'int_a_float_G', cs%diag%axesT1, cs%Time, &
2330 'Area integrated change in Greenland floating ice-shelf thickness ' //&
2331 'due to surface accum+melt during a DT_THERM time step', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2332 cs%id_Gr_bdott = register_scalar_field('ice_shelf_model', 'int_b_G', cs%diag%axesT1, cs%Time, &
2333 'Area integrated change in Greenland floating ice-shelf thickness '//&
2334 'due to basal accum+melt during a DT_THERM time step', 'm3', conversion=us%Z_to_m*us%L_to_m**2)
2335 cs%id_Gr_bdott_melt = register_scalar_field('ice_shelf_model', 'int_b_melt_G', cs%diag%axesT1, cs%Time, &
2336 'Area integrated basal melt over Greenland ice shelves during a DT_THERM time step', &
2337 units='m3', conversion=us%Z_to_m*us%L_to_m**2)
2338 cs%id_Gr_bdott_accum = register_scalar_field('ice_shelf_model', 'int_b_accum_G', cs%diag%axesT1, cs%Time, &
2339 'Area integrated basal accumulation over Greenland ice shelves during a DT_THERM a time step', &
2340 units='m3', conversion=us%Z_to_m*us%L_to_m**2)
2341 cs%id_Gr_t_area = register_scalar_field('ice_shelf_model', 'tot_area_G', cs%diag%axesT1, cs%Time, &
2342 'Total area of Greenland ice sheet', 'm2', conversion=us%L_to_m**2)
2343 cs%id_Gr_f_area = register_scalar_field('ice_shelf_model', 'tot_area_float_G', cs%diag%axesT1, cs%Time, &
2344 'Total area of Greenland floating ice shelves', 'm2', conversion=us%L_to_m**2)
2345 cs%id_Gr_g_area = register_scalar_field('ice_shelf_model', 'tot_area_ground_G', cs%diag%axesT1, cs%Time, &
2346 'Total area of Greenland grounded ice sheet', 'm2', conversion=us%L_to_m**2)
2347 !scalars (area integrated rates over the Greenland ice sheet)
2348 cs%id_Gr_dvafdt = register_scalar_field('ice_shelf_model', 'int_vafdot_G', cs%diag%axesT1, cs%Time, &
2349 'Area integrated rate of change in Greenland ice-sheet volume above floatation', &
2350 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2351 cs%id_Gr_adot = register_scalar_field('ice_shelf_model', 'int_adot_G', cs%diag%axesT1, cs%Time, &
2352 'Area integrated rate of change in Greenland ice-sheet thickness due to surface accum+melt', &
2353 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2354 cs%id_Gr_g_adot = register_scalar_field('ice_shelf_model', 'int_adot_ground_G', cs%diag%axesT1, cs%Time, &
2355 'Area integrated rate of change in Greenland grounded ice-sheet thickness due to surface accum+melt', &
2356 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2357 cs%id_Gr_f_adot = register_scalar_field('ice_shelf_model', 'int_adot_float_G', cs%diag%axesT1, cs%Time, &
2358 'Area integrated rate of change in Greenland floating ice-shelf thickness due to surface accum+melt', &
2359 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2360 cs%id_Gr_bdot = register_scalar_field('ice_shelf_model', 'int_bdot_G', cs%diag%axesT1, cs%Time, &
2361 'Area integrated rate of change in Greenland ice-shelf thickness due to basal accum+melt', &
2362 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2363 cs%id_Gr_bdot_melt = register_scalar_field('ice_shelf_model', 'int_bdot_melt_G', cs%diag%axesT1, cs%Time, &
2364 'Area integrated basal melt rate over Greenland ice shelves', &
2365 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2366 cs%id_Gr_bdot_accum = register_scalar_field('ice_shelf_model', 'int_bdot_accum_G', cs%diag%axesT1, cs%Time, &
2367 'Area integrated basal accumulation rate over Greenland ice shelves', &
2368 units='m3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2369
2370 !Flags to calculate diagnostics related to surface/basal mass balance
2371 if (cs%id_adott>0 .or. cs%id_g_adott>0 .or. cs%id_f_adott>0 .or. &
2372 cs%id_adot >0 .or. cs%id_g_adot >0 .or. cs%id_f_adot >0 .or. &
2373 cs%id_Ant_adott>0 .or. cs%id_Ant_g_adott>0 .or. cs%id_Ant_f_adott>0 .or. &
2374 cs%id_Ant_adot >0 .or. cs%id_Ant_g_adot >0 .or. cs%id_Ant_f_adot >0 .or. &
2375 cs%id_Gr_adott>0 .or. cs%id_Gr_g_adott>0 .or. cs%id_Gr_f_adott>0 .or. &
2376 cs%id_Gr_adot >0 .or. cs%id_Gr_g_adot >0 .or. cs%id_Gr_f_adot >0) then
2377 cs%smb_diag=.true.
2378 else
2379 cs%smb_diag=.false.
2380 endif
2381
2382 if (cs%id_bdott>0 .or. cs%id_bdott_melt>0 .or. cs%id_bdott_accum>0 .or. &
2383 cs%id_bdot >0 .or. cs%id_bdot_melt >0 .or. cs%id_bdot_accum >0 .or. &
2384 cs%id_Ant_bdott>0 .or. cs%id_Ant_bdott_melt>0 .or. cs%id_Ant_bdott_accum>0 .or. &
2385 cs%id_Ant_bdot >0 .or. cs%id_Ant_bdot_melt >0 .or. cs%id_Ant_bdot_accum >0 .or. &
2386 cs%id_Gr_bdott>0 .or. cs%id_Gr_bdott_melt>0 .or. cs%id_Gr_bdott_accum>0 .or. &
2387 cs%id_Gr_bdot >0 .or. cs%id_Gr_bdot_melt >0 .or. cs%id_Gr_bdot_accum >0) then
2388 cs%bmb_diag=.true.
2389 else
2390 cs%bmb_diag=.false.
2391 endif
2392
2394
2395 if (present(forces_in)) call initialize_ice_shelf_forces(cs, ocn_grid, us, forces_in)
2396
2397end subroutine initialize_ice_shelf
2398
2399subroutine initialize_ice_shelf_fluxes(CS, ocn_grid, US, fluxes_in)
2400 type(ice_shelf_cs), pointer :: cs !< A pointer to the ice shelf control structure
2401 type(ocean_grid_type), pointer :: ocn_grid !< The calling ocean model's horizontal grid structure
2402 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
2403 type(forcing), target, intent(inout) :: fluxes_in !< A structure containing pointers to any
2404 !! possible thermodynamic or mass-flux forcing fields.
2405
2406 ! Local variables
2407 type(ocean_grid_type), pointer :: g => null() ! Pointers to grids for convenience.
2408 type(forcing), pointer :: fluxes => null()
2409 integer :: i, j, isd, ied, jsd, jed
2410
2411 g => cs%Grid
2412 isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
2413
2414 ! Allocate the arrays for passing ice-shelf data through the forcing type.
2415 if (.not. cs%solo_ice_sheet) then
2416 call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes.")
2417 ! GMM: the following assures that water/heat fluxes are just allocated
2418 ! when SHELF_THERMO = True. These fluxes are necessary if one wants to
2419 ! use either ENERGETICS_SFC_PBL (ALE mode) or BULKMIXEDLAYER (layer mode).
2420 call allocate_forcing_type(cs%Grid_in, fluxes_in, ustar=.true., shelf=.true., &
2421 press=.true., water=cs%isthermo, heat=cs%isthermo, shelf_sfc_accumulation=cs%active_shelf_dynamics, &
2422 tau_mag=.true.)
2423 else
2424 call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes in solo mode.")
2425 call allocate_forcing_type(cs%Grid_in, fluxes_in, ustar=.true., shelf=.true., &
2426 press=.true., shelf_sfc_accumulation=cs%active_shelf_dynamics, tau_mag=.true.)
2427 endif
2428 if (cs%rotate_index) then
2429 allocate(fluxes)
2430 call allocate_forcing_type(fluxes_in, cs%Grid, fluxes, turns=cs%turns)
2431 call rotate_forcing(fluxes_in, fluxes, cs%turns)
2432 else
2433 fluxes => fluxes_in
2434 endif
2435
2436 do j=jsd,jed ; do i=isd,ied
2437 if (g%areaT(i,j)>0.) fluxes%frac_shelf_h(i,j) = cs%ISS%area_shelf_h(i,j) / g%areaT(i,j)
2438 enddo ; enddo
2439 if (cs%debug) call hchksum(fluxes%frac_shelf_h, "IS init: frac_shelf_h", g%HI, haloshift=0)
2440 call add_shelf_pressure(ocn_grid, us, cs, fluxes)
2441
2442 if (cs%rotate_index) then
2443 call rotate_forcing(fluxes, fluxes_in, -cs%turns)
2444 call deallocate_forcing_type(fluxes)
2445 deallocate(fluxes)
2446 endif
2447
2448end subroutine initialize_ice_shelf_fluxes
2449
2450!> Allocate and initialize the ice-shelf forcing elements of a mechanical forcing type.
2451!! This forcing type is on the unrotated grid that is used outside of the MOM6 ice shelf code.
2452subroutine initialize_ice_shelf_forces(CS, ocn_grid, US, forces_in)
2453 type(ice_shelf_cs), pointer :: cs !< A pointer to the ice shelf control structure
2454 type(ocean_grid_type), pointer :: ocn_grid !< The calling ocean model's horizontal grid structure
2455 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
2456 type(mech_forcing), target, intent(inout) :: forces_in !< A structure with the driving mechanical forces
2457
2458 ! Local variables
2459 type(mech_forcing), pointer :: forces => null()
2460
2461 call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating forces.")
2462
2463 if ((ocn_grid%isc /= cs%Grid_in%isc) .or. (ocn_grid%iec /= cs%Grid_in%iec) .or. &
2464 (ocn_grid%jsc /= cs%Grid_in%jsc) .or. (ocn_grid%jec /= cs%Grid_in%jec)) &
2465 call mom_error(fatal,"initialize_ice_shelf_forces: Incompatible ocean and external ice shelf grids.")
2466
2467 call allocate_mech_forcing(cs%Grid_in, forces_in, ustar=.true., shelf=.true., press=.true., tau_mag=.true.)
2468 if (cs%rotate_index) then
2469 allocate(forces)
2470 call allocate_mech_forcing(forces_in, cs%Grid, forces)
2471 call rotate_mech_forcing(forces_in, cs%turns, forces)
2472 else
2473 forces=>forces_in
2474 endif
2475
2476 call add_shelf_forces(cs%grid, us, cs, forces, do_shelf_area=.not.cs%solo_ice_sheet, &
2477 external_call=.false.)
2478
2479 if (cs%rotate_index) then
2480 call rotate_mech_forcing(forces, -cs%turns, forces_in)
2481 call deallocate_mech_forcing(forces)
2482 endif
2483
2484end subroutine initialize_ice_shelf_forces
2485
2486!> Initializes shelf mass based on three options (file, zero and user)
2487subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim)
2488
2489 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2490 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
2491 type(ice_shelf_cs), pointer :: CS !< A pointer to the ice shelf control structure
2492 type(ice_shelf_state), intent(inout) :: ISS !< The ice shelf state type that is being updated
2493 logical, optional, intent(in) :: new_sim !< If present and false, this run is being restarted
2494
2495 integer :: i, j, is, ie, js, je
2496 logical :: read_shelf_area, new_sim_2
2497 character(len=240) :: config, inputdir, shelf_file, filename
2498 character(len=120) :: shelf_mass_var ! The name of shelf mass in the file.
2499 character(len=120) :: shelf_area_var ! The name of shelf area in the file.
2500 character(len=40) :: mdl = "MOM_ice_shelf"
2501 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2502
2503 new_sim_2 = .true. ; if (present(new_sim)) new_sim_2 = new_sim
2504
2505 call get_param(param_file, mdl, "ICE_SHELF_CONFIG", config, &
2506 "A string that specifies how the ice shelf is "//&
2507 "initialized. Valid options include:\n"//&
2508 " \tfile\t Read from a file.\n"//&
2509 " \tzero\t Set shelf mass to 0 everywhere.\n"//&
2510 " \tUSER\t Call USER_initialize_shelf_mass.\n", &
2511 fail_if_missing=.true.)
2512
2513 select case ( trim(config) )
2514 case ("file")
2515
2516 call time_interp_external_init()
2517
2518 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
2519 inputdir = slasher(inputdir)
2520
2521 call get_param(param_file, mdl, "SHELF_FILE", shelf_file, &
2522 "If DYNAMIC_SHELF_MASS = True, OVERRIDE_SHELF_MOVEMENT = True "//&
2523 "and ICE_SHELF_MASS_FROM_FILE = True, this is the file from "//&
2524 "which to read the shelf mass and area.", &
2525 default="shelf_mass.nc")
2526 call get_param(param_file, mdl, "SHELF_MASS_VAR", shelf_mass_var, &
2527 "The variable in SHELF_FILE with the shelf mass.", &
2528 default="shelf_mass")
2529 call get_param(param_file, mdl, "READ_SHELF_AREA", read_shelf_area, &
2530 "If true, also read the area covered by ice-shelf from SHELF_FILE.", &
2531 default=.false.)
2532
2533 filename = trim(slasher(inputdir))//trim(shelf_file)
2534 call log_param(param_file, mdl, "INPUTDIR/SHELF_FILE", filename)
2535
2536 cs%mass_handle = init_external_field(filename, shelf_mass_var, &
2537 mom_domain=cs%Grid_in%Domain, verbose=cs%debug)
2538
2539 if (read_shelf_area) then
2540 call get_param(param_file, mdl, "SHELF_AREA_VAR", shelf_area_var, &
2541 "The variable in SHELF_FILE with the shelf area.", &
2542 default="shelf_area")
2543
2544 cs%area_handle = init_external_field(filename, shelf_area_var, &
2545 mom_domain=cs%Grid_in%Domain)
2546 endif
2547
2548 if (.not.file_exists(filename, cs%Grid_in%Domain)) call mom_error(fatal, &
2549 " initialize_shelf_mass: Unable to open "//trim(filename))
2550
2551 case ("zero")
2552 do j=js,je ; do i=is,ie
2553 iss%mass_shelf(i,j) = 0.0
2554 iss%area_shelf_h(i,j) = 0.0
2555 enddo ; enddo
2556
2557 case ("USER")
2558 call user_initialize_shelf_mass(iss%mass_shelf, iss%area_shelf_h, &
2559 iss%h_shelf, iss%hmask, g, cs%US, cs%user_CS, param_file, new_sim_2)
2560
2561 case default ; call mom_error(fatal,"initialize_ice_shelf: "// &
2562 "Unrecognized ice shelf setup "//trim(config))
2563 end select
2564
2565end subroutine initialize_shelf_mass
2566
2567!> This subroutine applies net accumulation/ablation at the top surface to the dynamic ice shelf.
2568!! acc_rate[m-s]=surf_mass_flux/density_ice is ablation/accumulation rate
2569!! positive for accumulation negative for ablation
2570subroutine change_thickness_using_precip(CS, ISS, G, US, fluxes, time_step, Time)
2571 type(ice_shelf_cs), intent(in) :: CS !< A pointer to the ice shelf control structure
2572 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
2573 type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
2574 !! the ice-shelf state
2575 type(forcing), intent(in) :: fluxes !< A structure of surface fluxes that
2576 !! includes surface mass flux
2577 type(time_type), intent(in) :: Time !< The current model time
2578 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2579 real, intent(in) :: time_step !< The time step for this update [T ~> s].
2580
2581 ! locals
2582 integer :: i, j
2583 real :: I_rho_ice ! The specific volume of ice [R-1 ~> m3 kg-1]
2584
2585 i_rho_ice = 1.0 / cs%density_ice
2586
2587 !update time
2588! CS%Time = Time
2589
2590! CS%time_step = time_step
2591 ! update surface mass flux rate
2592! if (CS%surf_mass_flux_from_file) call update_surf_mass_flux(G, US, CS, ISS, Time)
2593
2594 do j=g%jsc,g%jec ; do i=g%isc,g%iec
2595 if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
2596
2597 if (-fluxes%shelf_sfc_mass_flux(i,j) * time_step * i_rho_ice < iss%h_shelf(i,j)) then
2598 iss%h_shelf(i,j) = iss%h_shelf(i,j) + fluxes%shelf_sfc_mass_flux(i,j) * time_step * i_rho_ice
2599 else
2600 ! the ice is about to ablate, so set thickness, area, and mask to zero
2601 ! NOTE: this is not mass conservative should maybe scale salt & heat flux for this cell
2602 iss%h_shelf(i,j) = 0.0
2603 iss%hmask(i,j) = 0.0
2604 iss%area_shelf_h(i,j) = 0.0
2605 endif
2606 iss%mass_shelf(i,j) = iss%h_shelf(i,j) * cs%density_ice
2607 endif
2608 enddo ; enddo
2609
2610 call pass_var(iss%area_shelf_h, g%domain, complete=.false.)
2611 call pass_var(iss%h_shelf, g%domain, complete=.false.)
2612 call pass_var(iss%hmask, g%domain, complete=.false.)
2613 call pass_var(iss%mass_shelf, g%domain, complete=.true.)
2614
2615end subroutine change_thickness_using_precip
2616
2617
2618!> Updates the ice shelf mass using data from a file.
2619subroutine update_shelf_mass(G, US, CS, ISS, Time)
2620 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
2621 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2622 type(ice_shelf_cs), intent(in) :: CS !< A pointer to the ice shelf control structure
2623 type(ice_shelf_state), intent(inout) :: ISS !< The ice shelf state type that is being updated
2624 type(time_type), intent(in) :: Time !< The current model time
2625
2626 ! local variables
2627 integer :: i, j, is, ie, js, je
2628 real, allocatable, dimension(:,:) :: tmp2d ! Temporary array for storing ice shelf input data [R Z ~> kg m-2]
2629
2630 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2631
2632
2633 if (cs%rotate_index) then
2634 allocate(tmp2d(cs%Grid_in%isc:cs%Grid_in%iec,cs%Grid_in%jsc:cs%Grid_in%jec), source=0.0)
2635 else
2636 allocate(tmp2d(is:ie,js:je), source=0.0)
2637 endif
2638
2639 call time_interp_external(cs%mass_handle, time, tmp2d, scale=us%kg_m3_to_R*us%m_to_Z)
2640 call rotate_array(tmp2d, cs%turns, iss%mass_shelf)
2641 deallocate(tmp2d)
2642
2643 do j=js,je ; do i=is,ie
2644 iss%area_shelf_h(i,j) = 0.0
2645 iss%hmask(i,j) = 0.
2646 if (iss%mass_shelf(i,j) > 0.0) then
2647 iss%area_shelf_h(i,j) = g%areaT(i,j)
2648 iss%h_shelf(i,j) = iss%mass_shelf(i,j) / cs%density_ice
2649 iss%hmask(i,j) = 1.
2650 endif
2651 enddo ; enddo
2652
2653 !call USER_update_shelf_mass(ISS%mass_shelf, ISS%area_shelf_h, ISS%h_shelf, &
2654 ! ISS%hmask, CS%grid, CS%user_CS, Time, .true.)
2655
2656 if (cs%min_thickness_simple_calve > 0.0) then
2657 call ice_shelf_min_thickness_calve(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, &
2658 cs%min_thickness_simple_calve, halo=0)
2659 endif
2660
2661 call pass_var(iss%area_shelf_h, g%domain, complete=.false.)
2662 call pass_var(iss%h_shelf, g%domain, complete=.false.)
2663 call pass_var(iss%hmask, g%domain, complete=.false.)
2664 call pass_var(iss%mass_shelf, g%domain, complete=.true.)
2665
2666end subroutine update_shelf_mass
2667
2668!> Update the ice-shelf surface mass balance (SMB) field
2669subroutine update_ice_smb(CS, G, SMB, Time)
2670 type(ice_shelf_cs), pointer :: cs !< A pointer to the control structure returned
2671 !! by a previous call to initialize_ice_shelf.
2672 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
2673 real, dimension(SZDI_(G),SZDJ_(G)), &
2674 intent(inout) :: smb !< Ice surface mass balance parameter, often in [R Z T-1 ~> kg m-2 s-1]
2675 type(time_type), intent(in) :: time !< Current model time
2676
2677
2678 if (cs%time_varying_smb) then
2679 call time_interp_external(cs%smb_file, time, smb)
2680 if (cs%debug) call hchksum(smb, "update_ice_SMB", cs%Grid_in%HI, haloshift=0)
2681 endif
2682
2683
2684end subroutine update_ice_smb
2685
2686
2687!> Save the ice shelf restart file
2688subroutine ice_shelf_query(CS, G, frac_shelf_h, mass_shelf, data_override_shelf_fluxes)
2689 type(ice_shelf_cs), pointer :: cs !< ice shelf control structure
2690 type(ocean_grid_type), intent(in) :: g !< A pointer to an ocean grid control structure.
2691 real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: frac_shelf_h !< Ice shelf area fraction [nondim].
2692 real, optional, dimension(SZI_(G),SZJ_(G)), intent(out) :: mass_shelf !< Ice shelf mass [R Z ~> kg m-2]
2693 logical, optional :: data_override_shelf_fluxes !< If true, shelf fluxes can be written using
2694 !! the data_override capability (only for MOSAIC grids)
2695
2696 integer :: i, j
2697
2698 if (present(frac_shelf_h)) then
2699 do j=g%jsd,g%jed ; do i=g%isd,g%ied
2700 frac_shelf_h(i,j) = 0.0
2701 if (g%areaT(i,j)>0.) frac_shelf_h(i,j) = cs%ISS%area_shelf_h(i,j) / g%areaT(i,j)
2702 enddo ; enddo
2703 endif
2704
2705 if (present(mass_shelf)) then
2706 do j=g%jsd,g%jed ; do i=g%isd,g%ied
2707 mass_shelf(i,j) = 0.0
2708 if (g%areaT(i,j)>0.) mass_shelf(i,j) = cs%ISS%mass_shelf(i,j)
2709 enddo ; enddo
2710 endif
2711
2712 if (present(data_override_shelf_fluxes)) then
2713 data_override_shelf_fluxes=.false.
2714 if (cs%active_shelf_dynamics) data_override_shelf_fluxes = cs%data_override_shelf_fluxes
2715 endif
2716
2717end subroutine ice_shelf_query
2718
2719!> Save the ice shelf restart file
2720subroutine ice_shelf_save_restart(CS, Time, directory, time_stamped, filename_suffix)
2721 type(ice_shelf_cs), pointer :: cs !< ice shelf control structure
2722 type(time_type), intent(in) :: time !< model time at this call
2723 character(len=*), optional, intent(in) :: directory !< An optional directory into which to write
2724 !! these restart files.
2725 logical, optional, intent(in) :: time_stamped !< f true, the restart file names include
2726 !! a unique time stamp. The default is false.
2727 character(len=*), optional, intent(in) :: filename_suffix !< An optional suffix (e.g., a
2728 !! time-stamp) to append to the restart file names.
2729 ! local variables
2730 type(ocean_grid_type), pointer :: g => null()
2731 character(len=200) :: restart_dir
2732
2733 g => cs%grid
2734
2735 if (present(directory)) then ; restart_dir = directory
2736 else ; restart_dir = cs%restart_output_dir ; endif
2737
2738 call save_restart(restart_dir, time, cs%grid_in, cs%restart_CSp, time_stamped)
2739
2740end subroutine ice_shelf_save_restart
2741
2742!> Deallocates all memory associated with this module
2743subroutine ice_shelf_end(CS)
2744 type(ice_shelf_cs), pointer :: cs !< A pointer to the ice shelf control structure
2745
2746 if (.not.associated(cs)) return
2747
2748 call ice_shelf_state_end(cs%ISS)
2749
2750 if (cs%active_shelf_dynamics) call ice_shelf_dyn_end(cs%dCS)
2751
2752 call mom_is_diag_mediator_end(cs%diag)
2753 deallocate(cs)
2754
2755end subroutine ice_shelf_end
2756
2757!> This routine is for stepping a stand-alone ice shelf model without an ocean.
2758subroutine solo_step_ice_shelf(CS, time_interval, nsteps, Time, min_time_step_in, fluxes_in)
2759 type(ice_shelf_cs), pointer :: cs !< A pointer to the ice shelf control structure
2760 type(time_type), intent(in) :: time_interval !< The time interval for this update [s].
2761 integer, intent(inout) :: nsteps !< The running number of ice shelf steps.
2762 type(time_type), intent(inout) :: time !< The current model time
2763 real, optional, intent(in) :: min_time_step_in !< The minimum permitted time step [T ~> s].
2764 type(forcing), optional, target, intent(inout) :: fluxes_in !< A structure containing pointers to any
2765 !! possible thermodynamic or mass-flux forcing fields.
2766 type(ocean_grid_type), pointer :: g => null() ! A pointer to the ocean's grid structure
2767 type(unit_scale_type), pointer :: us => null() ! Pointer to a structure containing
2768 ! various unit conversion factors
2769 type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
2770 !! the ice-shelf state
2771 real :: remaining_time ! The remaining time in this call [T ~> s]
2772 real :: time_step ! The internal time step during this call [T ~> s]
2773 real :: full_time_step ! The external time step (sum of internal time steps) during this call [T ~> s]
2774 real :: ifull_time_step ! The inverse of the external time step [T-1 ~> s-1]
2775 real :: min_time_step ! The minimal required timestep that would indicate a fatal problem [T ~> s]
2776 character(len=240) :: mesg
2777 logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
2778 logical :: coupled_gl ! If true the grounding line position is determined based on
2779 ! coupled ice-ocean dynamics.
2780 integer :: is, ie, js, je, i, j
2781 real :: vaf0, vaf0_a, vaf0_g !The previous volumes above floatation
2782 !for all ice sheets, Antarctica only, or Greenland only [Z L2 ~> m3]
2783 real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: &
2784 dh_adott_sum, & ! Surface melt/accumulation over a full time step, used for diagnostics [Z ~> m]
2785 dh_adott ! Surface melt/accumulation over a partial time step, used for diagnostics [Z ~> m]
2786
2787 g => cs%grid
2788 us => cs%US
2789 iss => cs%ISS
2790 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2791
2792 remaining_time = time_to_real(time_interval, scale=us%s_to_T)
2793 full_time_step = remaining_time
2794 ifull_time_step = 1./full_time_step
2795
2796 if (present (min_time_step_in)) then
2797 min_time_step = min_time_step_in
2798 else
2799 min_time_step = 1000.0*us%s_to_T ! At 1 km resolution this would imply ice is moving at ~1 meter per second
2800 endif
2801
2802 write (mesg,*) "TIME in ice shelf call, yrs: ", time_to_real(time)/(365. * 86400.)
2803 call mom_mesg("solo_step_ice_shelf: "//mesg, 5)
2804
2805 iss%dhdt_shelf(:,:) = iss%h_shelf(:,:)
2806
2807 dh_adott(:,:)=0.0
2808
2809 if (cs%smb_diag) dh_adott_sum(:,:) = 0.0
2810
2811 !calculate previous volumes above floatation
2812 if (cs%id_dvafdt > 0) call volume_above_floatation(cs%dCS, g, iss, vaf0) !all ice sheet
2813 if (cs%id_Ant_dvafdt > 0) call volume_above_floatation(cs%dCS, g, iss, vaf0_a, hemisphere=0) !Antarctica only
2814 if (cs%id_Gr_dvafdt > 0) call volume_above_floatation(cs%dCS, g, iss, vaf0_g, hemisphere=1) !Greenland only
2815
2816 do while (remaining_time > 0.0)
2817 nsteps = nsteps+1
2818
2819 ! If time_interval is not too long, this is unnecessary.
2820 time_step = min(ice_time_step_cfl(cs%dCS, iss, g), remaining_time)
2821
2822 write (mesg,*) "Ice model timestep = ", us%T_to_s*time_step, " seconds"
2823 if ((time_step < min_time_step) .and. (time_step < remaining_time)) then
2824 call mom_error(fatal, "MOM_ice_shelf:solo_step_ice_shelf: abnormally small timestep "//mesg)
2825 else
2826 call mom_mesg("solo_step_ice_shelf: "//mesg, 5)
2827 endif
2828
2829 if (cs%smb_diag) dh_adott(is:ie,js:je) = iss%h_shelf(is:ie,js:je)
2830 call change_thickness_using_precip(cs, iss, g, us, fluxes_in, time_step, time)
2831 if (cs%smb_diag) dh_adott_sum(is:ie,js:je) = dh_adott_sum(is:ie,js:je) + &
2832 (iss%h_shelf(is:ie,js:je) - dh_adott(is:ie,js:je))
2833 remaining_time = remaining_time - time_step
2834
2835 ! If the last mini-timestep is a day or less, we cannot expect velocities to change by much.
2836 ! Do not update the velocities if the last step is very short.
2837 update_ice_vel = ((time_step > min_time_step) .or. (remaining_time > 0.0))
2838 coupled_gl = .false.
2839
2840 call update_ice_shelf(cs%dCS, iss, g, us, time_step, time, cs%calve_ice_shelf_bergs, &
2841 must_update_vel=update_ice_vel)
2842
2843 enddo
2844
2845 call write_ice_shelf_energy(cs%dCS, g, us, iss%mass_shelf, iss%area_shelf_h, time, &
2846 time_step=time_interval)
2847 do j=js,je ; do i=is,ie
2848 iss%dhdt_shelf(i,j) = (iss%h_shelf(i,j) - iss%dhdt_shelf(i,j)) * ifull_time_step
2849 enddo ; enddo
2850
2851 call enable_averages(full_time_step, time, cs%diag)
2852 if (cs%id_area_shelf_h > 0) call post_data(cs%id_area_shelf_h ,iss%area_shelf_h,cs%diag)
2853 if (cs%id_h_shelf > 0) call post_data(cs%id_h_shelf ,iss%h_shelf ,cs%diag)
2854 if (cs%id_dhdt_shelf > 0) call post_data(cs%id_dhdt_shelf ,iss%dhdt_shelf ,cs%diag)
2855 if (cs%id_h_mask > 0) call post_data(cs%id_h_mask ,iss%hmask ,cs%diag)
2856 if (cs%id_shelf_sfc_mass_flux > 0) call post_data(cs%id_shelf_sfc_mass_flux, fluxes_in%shelf_sfc_mass_flux, cs%diag)
2857 call process_and_post_scalar_data(cs, vaf0, vaf0_a, vaf0_g, ifull_time_step, dh_adott, dh_adott*0.0)
2858 call disable_averaging(cs%diag)
2859
2860 call is_dynamics_post_data(full_time_step, time, cs%dCS, iss, g)
2861end subroutine solo_step_ice_shelf
2862
2863!> Post_data calls for ice-sheet scalars
2864subroutine process_and_post_scalar_data(CS, vaf0, vaf0_A, vaf0_G, Itime_step, dh_adott, dh_bdott)
2865 type(ice_shelf_cs), pointer :: CS !< A pointer to the ice shelf control structure
2866 real :: vaf0 !< The previous volumes above floatation for all ice sheets [Z L2 ~> m3]
2867 real :: vaf0_A !< The previous volumes above floatation for the Antarctic ice sheet [Z L2 ~> m3]
2868 real :: vaf0_G !< The previous volumes above floatation for the Greenland ice sheet [Z L2 ~> m3]
2869 real :: Itime_step !< Inverse of the time step [T-1 ~> s-1]
2870 real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: dh_adott !< Surface (plus basal if solo shelf mode)
2871 !! melt/accumulation over a time step [Z ~> m]
2872 real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: dh_bdott !< Surface (plus basal if solo shelf mode)
2873 !! melt/accumulation over a time step [Z ~> m]
2874
2875 ! Local variables
2876 real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: tmp ! Temporary field used when calculating diagnostics [various]
2877 real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: ones ! Temporary field used when calculating diagnostics [various]
2878 real :: vaf ! The current ice-sheet volume above floatation [Z L2 ~> m3]
2879 real :: val ! Temporary value when calculating scalar diagnostics [various]
2880 type(ocean_grid_type), pointer :: G => null() ! A pointer to the ocean's grid structure
2881 type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing various unit conversion factors
2882 type(ice_shelf_state), pointer :: ISS => null() ! A structure with elements that describe the ice-shelf state
2883 integer :: is, ie, js, je, i, j
2884
2885 g => cs%grid
2886 us => cs%US
2887 iss => cs%ISS
2888 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2889
2890 !---ALL ICE SHEET---!
2891 if (cs%id_vaf > 0 .or. cs%id_dvafdt > 0) & !calculate current volume above floatation (vaf)
2892 call volume_above_floatation(cs%dCS, g, iss, vaf)
2893 if (cs%id_vaf > 0) call post_scalar_data(cs%id_vaf ,vaf ,cs%diag) !current vaf
2894 if (cs%id_dvafdt > 0) call post_scalar_data(cs%id_dvafdt,(vaf-vaf0)*itime_step,cs%diag) !d(vaf)/dt
2895 if (cs%id_adott > 0 .or. cs%id_adot > 0) then !surface accumulation - surface melt
2896 val = integrate_over_ice_sheet_area(g, iss, dh_adott, unscale=us%Z_to_m)
2897 if (cs%id_adott > 0) call post_scalar_data(cs%id_adott,val ,cs%diag)
2898 if (cs%id_adot > 0) call post_scalar_data(cs%id_adot ,val*itime_step,cs%diag)
2899 endif
2900 if (cs%id_g_adott > 0 .or. cs%id_g_adot > 0) then !grounded only: surface accumulation - surface melt
2901 call masked_var_grounded(g,cs%dCS,dh_adott,tmp)
2902 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=us%Z_to_m)
2903 if (cs%id_g_adott > 0) call post_scalar_data(cs%id_g_adott,val ,cs%diag)
2904 if (cs%id_g_adot > 0) call post_scalar_data(cs%id_g_adot ,val*itime_step,cs%diag)
2905 endif
2906 if (cs%id_f_adott > 0 .or. cs%id_f_adot > 0) then !floating only: surface accumulation - surface melt
2907 call masked_var_grounded(g,cs%dCS,dh_adott,tmp)
2908 do j=js,je ; do i=is,ie
2909 tmp(i,j) = dh_adott(i,j) - tmp(i,j)
2910 enddo ; enddo
2911 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=us%Z_to_m)
2912 if (cs%id_f_adott > 0) call post_scalar_data(cs%id_f_adott,val ,cs%diag)
2913 if (cs%id_f_adot > 0) call post_scalar_data(cs%id_f_adot ,val*itime_step,cs%diag)
2914 endif
2915 if (cs%id_bdott > 0 .or. cs%id_bdot > 0) then !bottom accumulation - bottom melt
2916 val = integrate_over_ice_sheet_area(g, iss, dh_bdott, unscale=us%Z_to_m)
2917 if (cs%id_bdott > 0) call post_scalar_data(cs%id_bdott,val ,cs%diag)
2918 if (cs%id_bdot > 0) call post_scalar_data(cs%id_bdot ,val*itime_step,cs%diag)
2919 endif
2920 if (cs%id_bdott_melt > 0 .or. cs%id_bdot_melt > 0) then !bottom melt
2921 tmp(:,:)=0.0
2922 do j=js,je ; do i=is,ie
2923 if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j)
2924 enddo ; enddo
2925 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=us%Z_to_m)
2926 if (cs%id_bdott_melt > 0) call post_scalar_data(cs%id_bdott_melt,val ,cs%diag)
2927 if (cs%id_bdot_melt > 0) call post_scalar_data(cs%id_bdot_melt ,val*itime_step,cs%diag)
2928 endif
2929 if (cs%id_bdott_accum > 0 .or. cs%id_bdot_accum > 0) then !bottom accumulation
2930 tmp(:,:)=0.0
2931 do j=js,je ; do i=is,ie
2932 if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j)
2933 enddo ; enddo
2934 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=us%Z_to_m)
2935 if (cs%id_bdott_accum > 0) call post_scalar_data(cs%id_bdott_accum,val ,cs%diag)
2936 if (cs%id_bdot_accum > 0) call post_scalar_data(cs%id_bdot_accum ,val*itime_step,cs%diag)
2937 endif
2938 if (cs%id_t_area > 0) then !ice sheet area
2939 tmp(:,:) = 1.0 ; val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=1.0)
2940 call post_scalar_data(cs%id_t_area,val,cs%diag)
2941 endif
2942 if (cs%id_g_area > 0 .or. cs%id_f_area > 0) then
2943 ones(:,:) = 1.0 ; call masked_var_grounded(g, cs%dCS, ones, tmp)
2944 if (cs%id_g_area > 0) then !grounded only ice sheet area
2945 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=1.0)
2946 call post_scalar_data(cs%id_g_area,val,cs%diag)
2947 endif
2948 if (cs%id_f_area > 0) then !floating only ice sheet area (ice shelf area)
2949 val = integrate_over_ice_sheet_area(g, iss, 1.0-tmp, unscale=1.0)
2950 call post_scalar_data(cs%id_f_area,val,cs%diag)
2951 endif
2952 endif
2953
2954 !---ANTARCTICA ONLY---!
2955 if (cs%id_Ant_vaf > 0 .or. cs%id_Ant_dvafdt > 0) & !calculate current volume above floatation (vaf)
2956 call volume_above_floatation(cs%dCS, g, iss, vaf, hemisphere=0)
2957 if (cs%id_Ant_vaf > 0) call post_scalar_data(cs%id_Ant_vaf ,vaf ,cs%diag) !current vaf
2958 if (cs%id_Ant_dvafdt > 0) call post_scalar_data(cs%id_Ant_dvafdt,(vaf-vaf0_a)*itime_step,cs%diag) !d(vaf)/dt
2959 if (cs%id_Ant_adott > 0 .or. cs%id_Ant_adot > 0) then !surface accumulation - surface melt
2960 val = integrate_over_ice_sheet_area(g, iss, dh_adott, unscale=us%Z_to_m, hemisphere=0)
2961 if (cs%id_Ant_adott > 0) call post_scalar_data(cs%id_Ant_adott,val ,cs%diag)
2962 if (cs%id_Ant_adot > 0) call post_scalar_data(cs%id_Ant_adot ,val*itime_step,cs%diag)
2963 endif
2964 if (cs%id_Ant_g_adott > 0 .or. cs%id_Ant_g_adot > 0) then !grounded only: surface accumulation - surface melt
2965 call masked_var_grounded(g,cs%dCS,dh_adott,tmp)
2966 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=us%Z_to_m, hemisphere=0)
2967 if (cs%id_Ant_g_adott > 0) call post_scalar_data(cs%id_Ant_g_adott,val ,cs%diag)
2968 if (cs%id_Ant_g_adot > 0) call post_scalar_data(cs%id_Ant_g_adot ,val*itime_step,cs%diag)
2969 endif
2970 if (cs%id_Ant_f_adott > 0 .or. cs%id_Ant_f_adot > 0) then !floating only: surface accumulation - surface melt
2971 call masked_var_grounded(g,cs%dCS,dh_adott,tmp)
2972 do j=js,je ; do i=is,ie
2973 tmp(i,j) = dh_adott(i,j) - tmp(i,j)
2974 enddo ; enddo
2975 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=us%Z_to_m, hemisphere=0)
2976 if (cs%id_Ant_f_adott > 0) call post_scalar_data(cs%id_Ant_f_adott,val ,cs%diag)
2977 if (cs%id_Ant_f_adot > 0) call post_scalar_data(cs%id_Ant_f_adot ,val*itime_step,cs%diag)
2978 endif
2979 if (cs%id_Ant_bdott > 0 .or. cs%id_Ant_bdot > 0) then !bottom accumulation - bottom melt
2980 val = integrate_over_ice_sheet_area(g, iss, dh_bdott, unscale=us%Z_to_m, hemisphere=0)
2981 if (cs%id_Ant_bdott > 0) call post_scalar_data(cs%id_Ant_bdott,val ,cs%diag)
2982 if (cs%id_Ant_bdot > 0) call post_scalar_data(cs%id_Ant_bdot ,val*itime_step,cs%diag)
2983 endif
2984 if (cs%id_Ant_bdott_melt > 0 .or. cs%id_Ant_bdot_melt > 0) then !bottom melt
2985 tmp(:,:)=0.0
2986 do j=js,je ; do i=is,ie
2987 if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j)
2988 enddo ; enddo
2989 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=us%Z_to_m, hemisphere=0)
2990 if (cs%id_Ant_bdott_melt > 0) call post_scalar_data(cs%id_Ant_bdott_melt,val ,cs%diag)
2991 if (cs%id_Ant_bdot_melt > 0) call post_scalar_data(cs%id_Ant_bdot_melt ,val*itime_step,cs%diag)
2992 endif
2993 if (cs%id_Ant_bdott_accum > 0 .or. cs%id_Ant_bdot_accum > 0) then !bottom accumulation
2994 tmp(:,:)=0.0
2995 do j=js,je ; do i=is,ie
2996 if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j)
2997 enddo ; enddo
2998 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=us%Z_to_m, hemisphere=0)
2999 if (cs%id_Ant_bdott_accum > 0) call post_scalar_data(cs%id_Ant_bdott_accum,val ,cs%diag)
3000 if (cs%id_Ant_bdot_accum > 0) call post_scalar_data(cs%id_Ant_bdot_accum ,val*itime_step,cs%diag)
3001 endif
3002 if (cs%id_Ant_t_area > 0) then !ice sheet area
3003 tmp(:,:) = 1.0 ; val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=1.0, hemisphere=0)
3004 call post_scalar_data(cs%id_Ant_t_area,val,cs%diag)
3005 endif
3006 if (cs%id_Ant_g_area > 0 .or. cs%id_Ant_f_area > 0) then
3007 ones(:,:) = 1.0 ; call masked_var_grounded(g, cs%dCS, ones, tmp)
3008 if (cs%id_Ant_g_area > 0) then !grounded only ice sheet area
3009 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=1.0, hemisphere=0)
3010 call post_scalar_data(cs%id_Ant_g_area,val,cs%diag)
3011 endif
3012 if (cs%id_Ant_f_area > 0) then !floating only ice sheet area (ice shelf area)
3013 val = integrate_over_ice_sheet_area(g, iss, 1.0-tmp, unscale=1.0, hemisphere=0)
3014 call post_scalar_data(cs%id_Ant_f_area,val,cs%diag)
3015 endif
3016 endif
3017
3018 !---GREENLAND ONLY---!
3019 if (cs%id_Gr_vaf > 0 .or. cs%id_Gr_dvafdt > 0) & !calculate current volume above floatation (vaf)
3020 call volume_above_floatation(cs%dCS, g, iss, vaf, hemisphere=1)
3021 if (cs%id_Gr_vaf > 0) call post_scalar_data(cs%id_Gr_vaf ,vaf ,cs%diag) !current vaf
3022 if (cs%id_Gr_dvafdt > 0) call post_scalar_data(cs%id_Gr_dvafdt,(vaf-vaf0_a)*itime_step,cs%diag) !d(vaf)/dt
3023 if (cs%id_Gr_adott > 0 .or. cs%id_Gr_adot > 0) then !surface accumulation - surface melt
3024 val = integrate_over_ice_sheet_area(g, iss, dh_adott, unscale=us%Z_to_m, hemisphere=1)
3025 if (cs%id_Gr_adott > 0) call post_scalar_data(cs%id_Gr_adott,val ,cs%diag)
3026 if (cs%id_Gr_adot > 0) call post_scalar_data(cs%id_Gr_adot ,val*itime_step,cs%diag)
3027 endif
3028 if (cs%id_Gr_g_adott > 0 .or. cs%id_Gr_g_adot > 0) then !grounded only: surface accumulation - surface melt
3029 call masked_var_grounded(g,cs%dCS,dh_adott,tmp)
3030 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=us%Z_to_m, hemisphere=1)
3031 if (cs%id_Gr_g_adott > 0) call post_scalar_data(cs%id_Gr_g_adott,val ,cs%diag)
3032 if (cs%id_Gr_g_adot > 0) call post_scalar_data(cs%id_Gr_g_adot ,val*itime_step,cs%diag)
3033 endif
3034 if (cs%id_Gr_f_adott > 0 .or. cs%id_Gr_f_adot > 0) then !floating only: surface accumulation - surface melt
3035 call masked_var_grounded(g,cs%dCS,dh_adott,tmp)
3036 do j=js,je ; do i=is,ie
3037 tmp(i,j) = dh_adott(i,j) - tmp(i,j)
3038 enddo ; enddo
3039 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=us%Z_to_m, hemisphere=1)
3040 if (cs%id_Gr_f_adott > 0) call post_scalar_data(cs%id_Gr_f_adott,val ,cs%diag)
3041 if (cs%id_Gr_f_adot > 0) call post_scalar_data(cs%id_Gr_f_adot ,val*itime_step,cs%diag)
3042 endif
3043 if (cs%id_Gr_bdott > 0 .or. cs%id_Gr_bdot > 0) then !bottom accumulation - bottom melt
3044 val = integrate_over_ice_sheet_area(g, iss, dh_bdott, unscale=us%Z_to_m, hemisphere=1)
3045 if (cs%id_Gr_bdott > 0) call post_scalar_data(cs%id_Gr_bdott,val ,cs%diag)
3046 if (cs%id_Gr_bdot > 0) call post_scalar_data(cs%id_Gr_bdot ,val*itime_step,cs%diag)
3047 endif
3048 if (cs%id_Gr_bdott_melt > 0 .or. cs%id_Gr_bdot_melt > 0) then !bottom melt
3049 tmp(:,:)=0.0
3050 do j=js,je ; do i=is,ie
3051 if (dh_bdott(i,j) < 0) tmp(i,j) = -dh_bdott(i,j)
3052 enddo ; enddo
3053 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=us%Z_to_m, hemisphere=1)
3054 if (cs%id_Gr_bdott_melt > 0) call post_scalar_data(cs%id_Gr_bdott_melt,val ,cs%diag)
3055 if (cs%id_Gr_bdot_melt > 0) call post_scalar_data(cs%id_Gr_bdot_melt ,val*itime_step,cs%diag)
3056 endif
3057 if (cs%id_Gr_bdott_accum > 0 .or. cs%id_Gr_bdot_accum > 0) then !bottom accumulation
3058 tmp(:,:)=0.0
3059 do j=js,je ; do i=is,ie
3060 if (dh_bdott(i,j) > 0) tmp(i,j) = dh_bdott(i,j)
3061 enddo ; enddo
3062 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=us%Z_to_m, hemisphere=1)
3063 if (cs%id_Gr_bdott_accum > 0) call post_scalar_data(cs%id_Gr_bdott_accum,val ,cs%diag)
3064 if (cs%id_Gr_bdot_accum > 0) call post_scalar_data(cs%id_Gr_bdot_accum ,val*itime_step,cs%diag)
3065 endif
3066 if (cs%id_Gr_t_area > 0) then !ice sheet area
3067 tmp(:,:) = 1.0 ; val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=1.0, hemisphere=1)
3068 call post_scalar_data(cs%id_Gr_t_area,val,cs%diag)
3069 endif
3070 if (cs%id_Gr_g_area > 0 .or. cs%id_Gr_f_area > 0) then
3071 ones(:,:) = 1.0 ; call masked_var_grounded(g, cs%dCS, ones, tmp)
3072 if (cs%id_Gr_g_area > 0) then !grounded only ice sheet area
3073 val = integrate_over_ice_sheet_area(g, iss, tmp, unscale=1.0, hemisphere=1)
3074 call post_scalar_data(cs%id_Gr_g_area,val,cs%diag)
3075 endif
3076 if (cs%id_Gr_f_area > 0) then !floating only ice sheet area (ice shelf area)
3077 val = integrate_over_ice_sheet_area(g, iss, 1.0-tmp, unscale=1.0, hemisphere=1)
3078 call post_scalar_data(cs%id_Gr_f_area,val,cs%diag)
3079 endif
3080 endif
3081end subroutine process_and_post_scalar_data
3082
3083!> Initialize ice surface mass balance field that is held constant over time
3084subroutine initialize_ice_smb(CS, SMB, G, US, PF)
3085 type(ice_shelf_cs), pointer :: cs !< A pointer to the control structure returned
3086 !! by a previous call to initialize_ice_shelf.
3087 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
3088 real, dimension(SZDI_(G),SZDJ_(G)), &
3089 intent(inout) :: smb !< Ice surface mass balance parameter, often in [R Z T-1 ~> kg m-2 s-1]
3090 type(unit_scale_type), intent(in) :: us !< A structure containing unit conversion factors
3091 type(param_file_type), intent(in) :: pf !< A structure to parse for run-time parameters
3092
3093 real :: smb_val ! Constant ice surface mass balance parameter, often in [R Z T-1 ~> kg m-2 s-1]
3094 character(len=40) :: mdl = "initialize_ice_SMB" ! This subroutine's name.
3095 character(len=200) :: config
3096 character(len=200) :: varname
3097 character(len=200) :: inputdir, filename, smb_file
3098 logical :: smb_file_has_time
3099
3100 call get_param(pf, mdl, "ICE_SMB_CONFIG", config, &
3101 "This specifies how the initial ice surface mass balance parameter is specified. "//&
3102 "Valid values are: CONSTANT and FILE.", &
3103 default="CONSTANT")
3104
3105 if (trim(config)=="CONSTANT") then
3106 call get_param(pf, mdl, "SMB", smb_val, &
3107 "Surface mass balance.", units="kg m-2 s-1", default=0.0, scale=us%kg_m2s_to_RZ_T)
3108
3109 smb(:,:) = smb_val
3110
3111 elseif (trim(config)=="FILE") then
3112 call mom_mesg(" MOM_ice_shelf.F90, initialize_ice_shelf: reading SMB parameter")
3113 call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
3114 inputdir = slasher(inputdir)
3115
3116 call get_param(pf, mdl, "ICE_SMB_FILE", smb_file, &
3117 "The file from which the ice surface mass balance is read.", &
3118 default="ice_SMB.nc")
3119 filename = trim(inputdir)//trim(smb_file)
3120 call log_param(pf, mdl, "INPUTDIR/ICE_SMB_FILE", filename)
3121 call get_param(pf, mdl, "ICE_SMB_VARNAME", varname, &
3122 "The variable to use as surface mass balance.", &
3123 default="SMB")
3124 call get_param(pf, mdl, "ICE_SMB_TIME_VARYING", smb_file_has_time, &
3125 "The file from which the ice surface mass balance is read has a time axis.", &
3126 default=.false.)
3127 call log_param(pf, mdl, "ICE_SMB_TIME_VARYING", smb_file_has_time)
3128 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
3129 " initialize_ice_SMV_from_file: Unable to open "//trim(filename))
3130 if (smb_file_has_time) then
3131 cs%smb_file = init_external_field(filename, varname, mom_domain=g%Domain,correct_leap_year_inconsistency=.true.)
3132 cs%time_varying_smb = .true.
3133 else
3134 call mom_read_data(filename,trim(varname), smb, g%Domain, scale=us%kg_m2s_to_RZ_T)
3135 cs%time_varying_smb = .false.
3136 endif
3137 endif
3138end subroutine initialize_ice_smb
3139
3140!> \namespace mom_ice_shelf
3141!!
3142!! \section section_ICE_SHELF
3143!!
3144!! This module implements the thermodynamic aspects of ocean/ice-shelf
3145!! inter-actions using the MOM framework and coding style.
3146!!
3147!! Derived from code by Chris Little, early 2010.
3148!!
3149!! The ice-sheet dynamics subroutines do the following:
3150!! initialize_shelf_mass - Initializes the ice shelf mass distribution.
3151!! - Initializes h_shelf, h_mask, area_shelf_h
3152!! - CURRENTLY: initializes mass_shelf as well, but this is unnecessary, as mass_shelf is initialized based on
3153!! h_shelf and density_ice immediately afterwards. Possibly subroutine should be renamed
3154!! update_shelf_mass - updates ice shelf mass via netCDF file
3155!! USER_update_shelf_mass (TODO).
3156!! solo_step_ice_shelf - called only in ice-only mode.
3157!! shelf_calc_flux - after melt rate & fluxes are calculated, ice dynamics are done. Currently mass_shelf is
3158!! updated immediately after ice_shelf_advect in fully dynamic mode.
3159!!
3160!! NOTES: be aware that hmask(:,:) has a number of functions; it is used for front advancement,
3161!! for subroutines in the velocity solve, and for thickness boundary conditions (this last one may be removed).
3162!! in other words, interfering with its updates will have implications you might not expect.
3163!!
3164!! Overall issues: Many variables need better documentation and units and the
3165!! subgrid on which they are discretized.
3166!!
3167!! \subsection section_ICE_SHELF_equations ICE_SHELF equations
3168!!
3169!! The three fundamental equations are:
3170!! Heat flux
3171!! \f[ \qquad \rho_w C_{pw} \gamma_T (T_w - T_b) = \rho_i \dot{m} L_f \f]
3172!! Salt flux
3173!! \f[ \qquad \rho_w \gamma_s (S_w - S_b) = \rho_i \dot{m} S_b \f]
3174!! Freezing temperature
3175!! \f[ \qquad T_b = a S_b + b + c P \f]
3176!!
3177!! where ....
3178!!
3179!! \subsection section_ICE_SHELF_references References
3180!!
3181!! Asay-Davis, Xylar S., Stephen L. Cornford, Benjamin K. Galton-Fenzi, Rupert M. Gladstone, G. Hilmar Gudmundsson,
3182!! David M. Holland, Paul R. Holland, and Daniel F. Martin. Experimental design for three interrelated marine ice sheet
3183!! and ocean model intercomparison projects: MISMIP v. 3 (MISMIP+), ISOMIP v. 2 (ISOMIP+) and MISOMIP v. 1 (MISOMIP1).
3184!! Geoscientific Model Development 9, no. 7 (2016): 2471.
3185!!
3186!! Goldberg, D. N., et al. Investigation of land ice-ocean interaction with a fully coupled ice-ocean model: 1.
3187!! Model description and behavior. Journal of Geophysical Research: Earth Surface 117.F2 (2012).
3188!!
3189!! Goldberg, D. N., et al. Investigation of land ice-ocean interaction with a fully coupled ice-ocean model: 2.
3190!! Sensitivity to external forcings. Journal of Geophysical Research: Earth Surface 117.F2 (2012).
3191!!
3192!! Holland, David M., and Adrian Jenkins. Modeling thermodynamic ice-ocean interactions at the base of an ice shelf.
3193!! Journal of Physical Oceanography 29.8 (1999): 1787-1800.
3194
3195end module mom_ice_shelf