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