MOM6
MOM_ice_shelf.F90
Go to the documentation of this file.
1 !> Implements the thermodynamic aspects of ocean / ice-shelf interactions,
2 !! along with a crude placeholder for a later implementation of full
3 !! ice shelf dynamics, all using the MOM framework and coding style.
4 module mom_ice_shelf
5 
6 ! This file is part of MOM6. See LICENSE.md for the license.
7 
8 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
9 use mom_cpu_clock, only : clock_component, clock_routine
10 use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
11 use mom_diag_mediator, only : diag_mediator_init, set_diag_mediator_grid
12 use mom_diag_mediator, only : diag_ctrl, time_type, enable_averaging, disable_averaging
13 use mom_domains, only : mom_domains_init, clone_mom_domain
14 use mom_domains, only : pass_var, pass_vector, to_all, cgrid_ne, bgrid_ne, corner
15 use mom_dyn_horgrid, only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid
16 use mom_dyn_horgrid, only : rescale_dyn_horgrid_bathymetry
17 use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
18 use mom_file_parser, only : read_param, get_param, log_param, log_version, param_file_type
19 use mom_grid, only : mom_grid_init, ocean_grid_type
20 use mom_grid_initialize, only : set_grid_metrics
21 use mom_fixed_initialization, only : mom_initialize_topography
22 use mom_fixed_initialization, only : mom_initialize_rotation
23 use user_initialization, only : user_initialize_topography
24 use mom_io, only : field_exists, file_exists, mom_read_data, write_version_number
25 use mom_io, only : slasher, fieldtype
26 use mom_io, only : write_field, close_file, single_file, multiple
27 use mom_restart, only : register_restart_field, query_initialized, save_restart
28 use mom_restart, only : restart_init, restore_state, mom_restart_cs
29 use mom_time_manager, only : time_type, time_type_to_real, time_type_to_real, real_to_time
30 use mom_transcribe_grid, only : copy_dyngrid_to_mom_grid, copy_mom_grid_to_dyngrid
31 use mom_unit_scaling, only : unit_scale_type, unit_scaling_init, fix_restart_unit_scaling
32 use mom_variables, only : surface
33 use mom_forcing_type, only : forcing, allocate_forcing_type, mom_forcing_chksum
34 use mom_forcing_type, only : mech_forcing, allocate_mech_forcing, mom_mech_forcing_chksum
35 use mom_forcing_type, only : copy_common_forcing_fields
36 use mom_get_input, only : directories, get_mom_input
37 use mom_eos, only : calculate_density, calculate_density_derivs, calculate_tfreeze
38 use mom_eos, only : eos_type, eos_init
39 use mom_ice_shelf_dynamics, only : ice_shelf_dyn_cs, update_ice_shelf
40 use mom_ice_shelf_dynamics, only : register_ice_shelf_dyn_restarts, initialize_ice_shelf_dyn
41 use mom_ice_shelf_dynamics, only : ice_shelf_min_thickness_calve
42 use mom_ice_shelf_dynamics, only : ice_time_step_cfl, ice_shelf_dyn_end
43 use mom_ice_shelf_initialize, only : initialize_ice_thickness
44 !MJH use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary
45 use mom_ice_shelf_state, only : ice_shelf_state, ice_shelf_state_end, ice_shelf_state_init
46 use user_shelf_init, only : user_initialize_shelf_mass, user_update_shelf_mass
47 use user_shelf_init, only : user_ice_shelf_cs
48 use mom_coms, only : reproducing_sum, sum_across_pes
49 use mom_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum
50 use time_interp_external_mod, only : init_external_field, time_interp_external
51 use time_interp_external_mod, only : time_interp_external_init
52 use time_manager_mod, only : print_time
53 implicit none ; private
54 
55 #include <MOM_memory.h>
56 #ifdef SYMMETRIC_LAND_ICE
57 # define GRID_SYM_ .true.
58 #else
59 # define GRID_SYM_ .false.
60 #endif
61 
62 public shelf_calc_flux, add_shelf_flux, initialize_ice_shelf, ice_shelf_end
63 public ice_shelf_save_restart, solo_time_step, add_shelf_forces
64 
65 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
66 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
67 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
68 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
69 
70 !> Control structure that contains ice shelf parameters and diagnostics handles
71 type, public :: ice_shelf_cs ; private
72  ! Parameters
73  type(mom_restart_cs), pointer :: restart_csp => null() !< A pointer to the restart control
74  !! structure for the ice shelves
75  type(ocean_grid_type) :: grid !< Grid for the ice-shelf model
76  type(unit_scale_type), pointer :: &
77  us => null() !< A structure containing various unit conversion factors
78  !type(dyn_horgrid_type), pointer :: dG !< Dynamic grid for the ice-shelf model
79  type(ocean_grid_type), pointer :: ocn_grid => null() !< A pointer to the ocean model grid
80  !! The rest is private
81  real :: flux_factor = 1.0 !< A factor that can be used to turn off ice shelf
82  !! melting (flux_factor = 0) [nondim].
83  character(len=128) :: restart_output_dir = ' ' !< The directory in which to write restart files
84  type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
85  !! the ice-shelf state
86  type(ice_shelf_dyn_cs), pointer :: dcs => null() !< The control structure for the ice-shelf dynamics.
87 
88  real, pointer, dimension(:,:) :: &
89  utide => null() !< tidal velocity [m s-1]
90 
91  real :: ustar_bg !< A minimum value for ustar under ice shelves [Z T-1 ~> m s-1].
92  real :: cdrag !< drag coefficient under ice shelves [nondim].
93  real :: g_earth !< The gravitational acceleration [m s-2]
94  real :: cp !< The heat capacity of sea water [J kg-1 degC-1].
95  real :: rho0 !< A reference ocean density [kg m-3].
96  real :: cp_ice !< The heat capacity of fresh ice [J kg-1 degC-1].
97  real :: gamma_t !< The (fixed) turbulent exchange velocity in the
98  !< 2-equation formulation [m s-1].
99  real :: salin_ice !< The salinity of shelf ice [ppt].
100  real :: temp_ice !< The core temperature of shelf ice [degC].
101  real :: kv_ice !< The viscosity of ice [m2 s-1].
102  real :: density_ice !< A typical density of ice [kg m-3].
103  real :: rho_ice !< Nominal ice density [kg m-2 Z-1 ~> kg m-3].
104  real :: kv_molec !< The molecular kinematic viscosity of sea water [m2 s-1].
105  real :: kd_molec_salt!< The molecular diffusivity of salt [m2 s-1].
106  real :: kd_molec_temp!< The molecular diffusivity of heat [m2 s-1].
107  real :: lat_fusion !< The latent heat of fusion [J kg-1].
108  real :: gamma_t_3eq !< Nondimensional heat-transfer coefficient, used in the 3Eq. formulation
109  !< This number should be specified by the user.
110  real :: col_thick_melt_threshold !< if the mixed layer is below this threshold, melt rate
111  logical :: mass_from_file !< Read the ice shelf mass from a file every dt
112 
113  !!!! PHYSICAL AND NUMERICAL PARAMETERS FOR ICE DYNAMICS !!!!!!
114 
115  real :: time_step !< this is the shortest timestep that the ice shelf sees, and
116  !! is equal to the forcing timestep (it is passed in when the shelf
117  !! is initialized - so need to reorganize MOM driver.
118  !! it will be the prognistic timestep ... maybe.
119 
120  logical :: solo_ice_sheet !< whether the ice model is running without being
121  !! coupled to the ocean
122  logical :: gl_regularize !< whether to regularize the floatation condition
123  !! at the grounding line a la Goldberg Holland Schoof 2009
124  logical :: gl_couple !< whether to let the floatation condition be
125  !!determined by ocean column thickness means update_OD_ffrac
126  !! will be called (note: GL_regularize and GL_couple
127  !! should be exclusive)
128  real :: density_ocean_avg !< this does not affect ocean circulation OR thermodynamics
129  !! it is to estimate the gravitational driving force at the
130  !! shelf front (until we think of a better way to do it,
131  !! but any difference will be negligible)
132  logical :: calve_to_mask !< If true, calve any ice that passes outside of a masked area
133  real :: min_thickness_simple_calve !< min. ice shelf thickness criteria for calving [Z ~> m].
134  real :: t0 !< temperature at ocean surface in the restoring region [degC]
135  real :: s0 !< Salinity at ocean surface in the restoring region [ppt].
136  real :: input_flux !< Ice volume flux at an upstream open boundary [m3 s-1].
137  real :: input_thickness !< Ice thickness at an upstream open boundary [m].
138 
139  type(time_type) :: time !< The component's time.
140  type(eos_type), pointer :: eqn_of_state => null() !< Type that indicates the
141  !! equation of state to use.
142  logical :: active_shelf_dynamics !< True if the ice shelf mass changes as a result
143  !! the dynamic ice-shelf model.
144  logical :: override_shelf_movement !< If true, user code specifies the shelf movement
145  !! instead of using the dynamic ice-shelf mode.
146  logical :: isthermo !< True if the ice shelf can exchange heat and
147  !! mass with the underlying ocean.
148  logical :: threeeq !< If true, the 3 equation consistency equations are
149  !! used to calculate the flux at the ocean-ice
150  !! interface.
151  logical :: insulator !< If true, ice shelf is a perfect insulator
152  logical :: const_gamma !< If true, gamma_T is specified by the user.
153  logical :: find_salt_root !< If true, if true find Sbdry using a quadratic eq.
154  logical :: constant_sea_level !< if true, apply an evaporative, heat and salt
155  !! fluxes. It will avoid large increase in sea level.
156  real :: cutoff_depth !< depth above which melt is set to zero (>= 0).
157  real :: lambda1 !< liquidus coeff., Needed if find_salt_root = true
158  real :: lambda2 !< liquidus coeff., Needed if find_salt_root = true
159  real :: lambda3 !< liquidus coeff., Needed if find_salt_root = true
160  !>@{ Diagnostic handles
161  integer :: id_melt = -1, id_exch_vel_s = -1, id_exch_vel_t = -1, &
162  id_tfreeze = -1, id_tfl_shelf = -1, &
163  id_thermal_driving = -1, id_haline_driving = -1, &
164  id_u_ml = -1, id_v_ml = -1, id_sbdry = -1, &
165  id_h_shelf = -1, id_h_mask = -1, &
166  id_surf_elev = -1, id_bathym = -1, &
167  id_area_shelf_h = -1, &
168  id_ustar_shelf = -1, id_shelf_mass = -1, id_mass_flux = -1
169  !>@}
170 
171  integer :: id_read_mass !< An integer handle used in time interpolation of
172  !! the ice shelf mass read from a file
173  integer :: id_read_area !< An integer handle used in time interpolation of
174  !! the ice shelf mass read from a file
175 
176  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to control diagnostic output.
177  type(user_ice_shelf_cs), pointer :: user_cs => null() !< A pointer to the control structure for
178  !! user-supplied modifications to the ice shelf code.
179 
180  logical :: debug !< If true, write verbose checksums for debugging purposes
181  !! and use reproducible sums
182 end type ice_shelf_cs
183 
184 integer :: id_clock_shelf !< CPU Clock for the ice shelf code
185 integer :: id_clock_pass !< CPU Clock for group pass calls
186 
187 contains
188 
189 !> Calculates fluxes between the ocean and ice-shelf using the three-equations
190 !! formulation (optional to use just two equations).
191 !! See \ref section_ICE_SHELF_equations
192 subroutine shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
193  type(surface), intent(inout) :: state !< A structure containing fields that
194  !! describe the surface state of the ocean. The
195  !! intent is only inout to allow for halo updates.
196  type(forcing), intent(inout) :: fluxes !< structure containing pointers to any possible
197  !! thermodynamic or mass-flux forcing fields.
198  type(time_type), intent(in) :: time !< Start time of the fluxes.
199  real, intent(in) :: time_step !< Length of time over which
200  !! these fluxes will be applied [s].
201  type(ice_shelf_cs), pointer :: cs !< A pointer to the control structure
202  !! returned by a previous call to
203  !! initialize_ice_shelf.
204  type(mech_forcing), optional, intent(inout) :: forces !< A structure with the driving mechanical forces
205 
206  type(ocean_grid_type), pointer :: g => null() ! The grid structure used by the ice shelf.
207  type(unit_scale_type), pointer :: us => null() ! Pointer to a structure containing
208  ! various unit conversion factors
209  type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
210  !! the ice-shelf state
211 
212  real, dimension(SZI_(CS%grid)) :: &
213  rhoml, & !< Ocean mixed layer density [kg m-3].
214  dr0_dt, & !< Partial derivative of the mixed layer density
215  !< with temperature [kg m-3 degC-1].
216  dr0_ds, & !< Partial derivative of the mixed layer density
217  !< with salinity [kg m-3 ppt-1].
218  p_int !< The pressure at the ice-ocean interface [Pa].
219 
220  real, dimension(SZI_(CS%grid),SZJ_(CS%grid)) :: &
221  exch_vel_t, & !< Sub-shelf thermal exchange velocity [m s-1]
222  exch_vel_s !< Sub-shelf salt exchange velocity [m s-1]
223 
224  real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
225  mass_flux !< total mass flux of freshwater across
226  real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
227  haline_driving !< (SSS - S_boundary) ice-ocean
228  !! interface, positive for melting and negative for freezing.
229  !! This is computed as part of the ISOMIP diagnostics.
230  real, parameter :: vk = 0.40 !< Von Karman's constant - dimensionless
231  real :: zeta_n = 0.052 !> The fraction of the boundary layer over which the
232  !! viscosity is linearly increasing. (Was 1/8. Why?)
233  real, parameter :: rc = 0.20 ! critical flux Richardson number.
234  real :: i_zeta_n !< The inverse of ZETA_N.
235  real :: lf, i_lf !< Latent Heat of fusion [J kg-1] and its inverse.
236  real :: i_vk !< The inverse of VK.
237  real :: pr, sc !< The Prandtl number and Schmidt number [nondim].
238 
239  ! 3 equations formulation variables
240  real, dimension(SZDI_(CS%grid),SZDJ_(CS%grid)) :: &
241  sbdry !< Salinities in the ocean at the interface with the ice shelf [ppt].
242  real :: sbdry_it
243  real :: sbdry1, sbdry2, s_a, s_b, s_c ! use to find salt roots
244  real :: ds_it !< The interface salinity change during an iteration [ppt].
245  real :: hbl_neut !< The neutral boundary layer thickness [m].
246  real :: hbl_neut_h_molec !< The ratio of the neutral boundary layer thickness
247  !! to the molecular boundary layer thickness [nondim].
248  !### THESE ARE CURRENTLY POSITIVE UPWARD.
249  real :: wt_flux !< The vertical flux of heat just inside the ocean [degC m s-1].
250  real :: wb_flux !< The vertical flux of heat just inside the ocean [m2 s-3].
251  real :: db_ds !< The derivative of buoyancy with salinity [m s-2 ppt-1].
252  real :: db_dt !< The derivative of buoyancy with temperature [m s-2 degC-1].
253  real :: i_n_star, n_star_term, absf
254  real :: dins_dwb !< The partial derivative of I_n_star with wB_flux, in ???.
255  real :: dt_ustar, ds_ustar
256  real :: ustar_h
257  real :: gam_turb
258  real :: gam_mol_t, gam_mol_s
259  real :: rhocp
260  real :: i_rholf
261  real :: ln_neut
262  real :: mass_exch
263  real :: sb_min, sb_max
264  real :: ds_min, ds_max
265  ! Variables used in iterating for wB_flux.
266  real :: wb_flux_new, dwb, ddwb_dwb_in
267  real :: i_gam_t, i_gam_s, dg_dwb, idens
268  real :: u_at_h, v_at_h, isqrt2
269  logical :: sb_min_set, sb_max_set
270  logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
271  logical :: coupled_gl ! If true, the grouding line position is determined based on
272  ! coupled ice-ocean dynamics.
273 
274  real, parameter :: c2_3 = 2.0/3.0
275  character(len=160) :: mesg ! The text of an error message
276  integer :: i, j, is, ie, js, je, ied, jed, it1, it3
277  real, parameter :: rho_fw = 1000.0 ! fresh water density
278 
279  if (.not. associated(cs)) call mom_error(fatal, "shelf_calc_flux: "// &
280  "initialize_ice_shelf must be called before shelf_calc_flux.")
281  call cpu_clock_begin(id_clock_shelf)
282 
283  g => cs%grid ; us => cs%US
284  iss => cs%ISS
285 
286  ! useful parameters
287  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; ied = g%ied ; jed = g%jed
288  i_zeta_n = 1.0 / zeta_n
289  lf = cs%Lat_fusion
290  i_rholf = 1.0/(cs%Rho0*lf)
291  i_lf = 1.0 / lf
292  sc = cs%kv_molec/cs%kd_molec_salt
293  pr = cs%kv_molec/cs%kd_molec_temp
294  i_vk = 1.0/vk
295  rhocp = cs%Rho0 * cs%Cp
296  isqrt2 = 1.0/sqrt(2.0)
297 
298  !first calculate molecular component
299  gam_mol_t = 12.5 * (pr**c2_3) - 6
300  gam_mol_s = 12.5 * (sc**c2_3) - 6
301 
302  idens = 1.0/cs%density_ocean_avg
303 
304  ! GMM, zero some fields of the ice shelf structure (ice_shelf_CS)
305  ! these fields are already set to zero during initialization
306  ! However, they seem to be changed somewhere and, for diagnostic
307  ! reasons, it is better to set them to zero again.
308  exch_vel_t(:,:) = 0.0 ; exch_vel_s(:,:) = 0.0
309  iss%tflux_shelf(:,:) = 0.0 ; iss%water_flux(:,:) = 0.0
310  iss%salt_flux(:,:) = 0.0; iss%tflux_ocn(:,:) = 0.0
311  iss%tfreeze(:,:) = 0.0
312  ! define Sbdry to avoid Run-Time Check Failure, when melt is not computed.
313  haline_driving(:,:) = 0.0
314  sbdry(:,:) = state%sss(:,:)
315 
316  !update time
317  cs%Time = time
318 
319  if (cs%override_shelf_movement) then
320  cs%time_step = time_step
321  ! update shelf mass
322  if (cs%mass_from_file) call update_shelf_mass(g, cs, iss, time)
323  endif
324 
325  if (cs%debug) then
326  call hchksum(fluxes%frac_shelf_h, "frac_shelf_h before apply melting", g%HI, haloshift=0)
327  call hchksum(state%sst, "sst before apply melting", g%HI, haloshift=0)
328  call hchksum(state%sss, "sss before apply melting", g%HI, haloshift=0)
329  call hchksum(state%u, "u_ml before apply melting", g%HI, haloshift=0)
330  call hchksum(state%v, "v_ml before apply melting", g%HI, haloshift=0)
331  call hchksum(state%ocean_mass, "ocean_mass before apply melting", g%HI, haloshift=0)
332  endif
333 
334  do j=js,je
335  ! Find the pressure at the ice-ocean interface, averaged only over the
336  ! part of the cell covered by ice shelf.
337  do i=is,ie ; p_int(i) = cs%g_Earth * iss%mass_shelf(i,j) ; enddo
338 
339  ! Calculate insitu densities and expansion coefficients
340  call calculate_density(state%sst(:,j), state%sss(:,j), p_int, &
341  rhoml(:), is, ie-is+1, cs%eqn_of_state)
342  call calculate_density_derivs(state%sst(:,j), state%sss(:,j), p_int, &
343  dr0_dt, dr0_ds, is, ie-is+1, cs%eqn_of_state)
344 
345  do i=is,ie
346  ! set ustar_shelf to zero. This is necessary if shelf_mass_is_dynamic
347  ! but it won't make a difference otherwise.
348  fluxes%ustar_shelf(i,j)= 0.0
349 
350  ! DNG - to allow this everywhere Hml>0.0 allows for melting under grounded cells
351  ! propose instead to allow where Hml > [some threshold]
352 
353  if ((idens*state%ocean_mass(i,j) > cs%col_thick_melt_threshold) .and. &
354  (iss%area_shelf_h(i,j) > 0.0) .and. &
355  (cs%isthermo) .and. (state%Hml(i,j) > 0.0) ) then
356 
357  if (cs%threeeq) then
358  ! Iteratively determine a self-consistent set of fluxes, with the ocean
359  ! salinity just below the ice-shelf as the variable that is being
360  ! iterated for.
361  ! ### SHOULD I SET USTAR_SHELF YET?
362 
363  u_at_h = state%u(i,j)
364  v_at_h = state%v(i,j)
365 
366  !### I think that CS%utide**1 should be CS%utide**2
367  ! Also I think that if taux_shelf and tauy_shelf have been calculated by the
368  ! ocean stress calculation, they should be used here or later to set ustar_shelf. - RWH
369  fluxes%ustar_shelf(i,j) = max(cs%ustar_bg, us%m_to_Z*us%T_to_s * &
370  sqrt(cs%cdrag*((u_at_h**2 + v_at_h**2) + cs%utide(i,j)**1)))
371 
372  ustar_h = us%Z_to_m*us%s_to_T*fluxes%ustar_shelf(i,j)
373 
374  ! I think that the following can be deleted without causing any problems.
375  ! if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then
376  ! ! These arrays are supposed to be stress components at C-grid points, which is
377  ! ! inconsistent with what is coded up here.
378  ! state%taux_shelf(i,j) = ustar_h*ustar_h*CS%Rho0*Isqrt2
379  ! state%tauy_shelf(i,j) = state%taux_shelf(i,j)
380  ! endif
381 
382  ! Estimate the neutral ocean boundary layer thickness as the minimum of the
383  ! reported ocean mixed layer thickness and the neutral Ekman depth.
384  absf = 0.25*us%s_to_T*((abs(us%s_to_T*g%CoriolisBu(i,j)) + abs(us%s_to_T*g%CoriolisBu(i-1,j-1))) + &
385  (abs(us%s_to_T*g%CoriolisBu(i,j-1)) + abs(us%s_to_T*g%CoriolisBu(i-1,j))))
386  if (absf*state%Hml(i,j) <= vk*ustar_h) then ; hbl_neut = state%Hml(i,j)
387  else ; hbl_neut = (vk*ustar_h) / absf ; endif
388  hbl_neut_h_molec = zeta_n * ((hbl_neut * ustar_h) / (5.0 * cs%Kv_molec))
389 
390  ! Determine the mixed layer buoyancy flux, wB_flux.
391  db_ds = (cs%g_Earth / rhoml(i)) * dr0_ds(i)
392  db_dt = (cs%g_Earth / rhoml(i)) * dr0_dt(i)
393  ln_neut = 0.0 ; if (hbl_neut_h_molec > 1.0) ln_neut = log(hbl_neut_h_molec)
394 
395  if (cs%find_salt_root) then
396  ! read liquidus parameters
397 
398  s_a = cs%lambda1 * cs%Gamma_T_3EQ * cs%Cp
399 ! S_b = -CS%Gamma_T_3EQ*(CS%lambda2-CS%lambda3*p_int(i)-state%sst(i,j)) &
400 ! -LF*CS%Gamma_T_3EQ/35.0
401 
402  s_b = cs%Gamma_T_3EQ*cs%Cp*(cs%lambda2+cs%lambda3*p_int(i)- &
403  state%sst(i,j))-lf*cs%Gamma_T_3EQ/35.0
404  s_c = lf*(cs%Gamma_T_3EQ/35.0)*state%sss(i,j)
405 
406  !### Depending on the sign of S_b, one of these will be inaccurate!
407  sbdry1 = (-s_b + sqrt(s_b*s_b-4*s_a*s_c))/(2*s_a)
408  sbdry2 = (-s_b - sqrt(s_b*s_b-4*s_a*s_c))/(2*s_a)
409  sbdry(i,j) = max(sbdry1, sbdry2)
410  ! Safety check
411  if (sbdry(i,j) < 0.) then
412  write(mesg,*) 'state%sss(i,j) = ',state%sss(i,j), 'S_a, S_b, S_c', s_a, s_b, s_c
413  call mom_error(warning, mesg, .true.)
414  write(mesg,*) 'I,J,Sbdry1,Sbdry2',i,j,sbdry1,sbdry2
415  call mom_error(warning, mesg, .true.)
416  call mom_error(fatal, "shelf_calc_flux: Negative salinity (Sbdry).")
417  endif
418  else
419  ! Guess sss as the iteration starting point for the boundary salinity.
420  sbdry(i,j) = state%sss(i,j) ; sb_max_set = .false.
421  sb_min_set = .false.
422  endif !find_salt_root
423 
424  do it1 = 1,20
425  ! Determine the potential temperature at the ice-ocean interface.
426  call calculate_tfreeze(sbdry(i,j), p_int(i), iss%tfreeze(i,j), cs%eqn_of_state)
427 
428  dt_ustar = (state%sst(i,j) - iss%tfreeze(i,j)) * ustar_h
429  ds_ustar = (state%sss(i,j) - sbdry(i,j)) * ustar_h
430 
431  ! First, determine the buoyancy flux assuming no effects of stability
432  ! on the turbulence. Following H & J '99, this limit also applies
433  ! when the buoyancy flux is destabilizing.
434 
435  if (cs%const_gamma) then ! if using a constant gamma_T
436  ! note the different form, here I_Gam_T is NOT 1/Gam_T!
437  i_gam_t = cs%Gamma_T_3EQ
438  i_gam_s = cs%Gamma_T_3EQ/35.
439  else
440  gam_turb = i_vk * (ln_neut + (0.5 * i_zeta_n - 1.0))
441  i_gam_t = 1.0 / (gam_mol_t + gam_turb)
442  i_gam_s = 1.0 / (gam_mol_s + gam_turb)
443  endif
444 
445  wt_flux = dt_ustar * i_gam_t
446  wb_flux = db_ds * (ds_ustar * i_gam_s) + db_dt * wt_flux
447 
448  if (wb_flux > 0.0) then
449  ! The buoyancy flux is stabilizing and will reduce the tubulent
450  ! fluxes, and iteration is required.
451  n_star_term = (zeta_n/rc) * (hbl_neut * vk) / ustar_h**3
452  do it3 = 1,30
453  ! n_star <= 1.0 is the ratio of working boundary layer thickness
454  ! to the neutral thickness.
455  ! hBL = n_star*hBL_neut ; hSub = 1/8*n_star*hBL
456 
457  i_n_star = sqrt(1.0 + n_star_term * wb_flux)
458  dins_dwb = 0.5 * n_star_term / i_n_star
459  if (hbl_neut_h_molec > i_n_star**2) then
460  gam_turb = i_vk * ((ln_neut - 2.0*log(i_n_star)) + &
461  (0.5*i_zeta_n*i_n_star - 1.0))
462  dg_dwb = i_vk * ( -2.0 / i_n_star + (0.5 * i_zeta_n)) * dins_dwb
463  else
464  ! The layer dominated by molecular viscosity is smaller than
465  ! the assumed boundary layer. This should be rare!
466  gam_turb = i_vk * (0.5 * i_zeta_n*i_n_star - 1.0)
467  dg_dwb = i_vk * (0.5 * i_zeta_n) * dins_dwb
468  endif
469 
470  if (cs%const_gamma) then ! if using a constant gamma_T
471  ! note the different form, here I_Gam_T is NOT 1/Gam_T!
472  i_gam_t = cs%Gamma_T_3EQ
473  i_gam_s = cs%Gamma_T_3EQ/35.
474  else
475  i_gam_t = 1.0 / (gam_mol_t + gam_turb)
476  i_gam_s = 1.0 / (gam_mol_s + gam_turb)
477  endif
478 
479  wt_flux = dt_ustar * i_gam_t
480  wb_flux_new = db_ds * (ds_ustar * i_gam_s) + db_dt * wt_flux
481 
482  ! Find the root where dwB = 0.0
483  dwb = wb_flux_new - wb_flux
484  if (abs(wb_flux_new - wb_flux) < &
485  1e-4*(abs(wb_flux_new) + abs(wb_flux))) exit
486 
487  ddwb_dwb_in = -dg_dwb * (db_ds * (ds_ustar * i_gam_s**2) + &
488  db_dt * (dt_ustar * i_gam_t**2)) - 1.0
489  ! This is Newton's method without any bounds.
490  ! ### SHOULD BOUNDS BE NEEDED?
491  wb_flux_new = wb_flux - dwb / ddwb_dwb_in
492  enddo !it3
493  endif
494 
495  iss%tflux_ocn(i,j) = rhocp * wt_flux
496  exch_vel_t(i,j) = ustar_h * i_gam_t
497  exch_vel_s(i,j) = ustar_h * i_gam_s
498 
499  !Calculate the heat flux inside the ice shelf.
500 
501  !vertical adv/diff as in H+J 1999, eqns (26) & approx from (31).
502  ! Q_ice = rho_ice * CS%CP_Ice * K_ice * dT/dz (at interface)
503  !vertical adv/diff as in H+J 199, eqs (31) & (26)...
504  ! dT/dz ~= min( (lprec/(rho_ice*K_ice))*(CS%Temp_Ice-T_freeze) , 0.0 )
505  !If this approximation is not made, iterations are required... See H+J Fig 3.
506 
507  if (iss%tflux_ocn(i,j) <= 0.0) then ! Freezing occurs, so zero ice heat flux.
508  iss%water_flux(i,j) = i_lf * iss%tflux_ocn(i,j)
509  iss%tflux_shelf(i,j) = 0.0
510  else
511  if (cs%insulator) then
512  !no conduction/perfect insulator
513  iss%tflux_shelf(i,j) = 0.0
514  iss%water_flux(i,j) = i_lf * (- iss%tflux_shelf(i,j) + iss%tflux_ocn(i,j))
515 
516  else
517  ! With melting, from H&J 1999, eqs (31) & (26)...
518  ! Q_ice ~= cp_ice * (CS%Temp_Ice-T_freeze) * lprec
519  ! RhoLF*lprec = Q_ice + ISS%tflux_ocn(i,j)
520  ! lprec = (ISS%tflux_ocn(i,j)) / (LF + cp_ice * (T_freeze-CS%Temp_Ice))
521  iss%water_flux(i,j) = iss%tflux_ocn(i,j) / &
522  (lf + cs%CP_Ice * (iss%tfreeze(i,j) - cs%Temp_Ice))
523 
524  iss%tflux_shelf(i,j) = iss%tflux_ocn(i,j) - lf*iss%water_flux(i,j)
525  endif
526 
527  endif
528  !other options: dTi/dz linear through shelf
529  ! dTi_dz = (CS%Temp_Ice - ISS%tfreeze(i,j))/G%draft(i,j)
530  ! ISS%tflux_shelf(i,j) = - Rho_Ice * CS%CP_Ice * KTI * dTi_dz
531 
532 
533  if (cs%find_salt_root) then
534  exit ! no need to do interaction, so exit loop
535  else
536 
537  mass_exch = exch_vel_s(i,j) * cs%Rho0
538  sbdry_it = (state%sss(i,j) * mass_exch + cs%Salin_ice * &
539  iss%water_flux(i,j)) / (mass_exch + iss%water_flux(i,j))
540  ds_it = sbdry_it - sbdry(i,j)
541  if (abs(ds_it) < 1e-4*(0.5*(state%sss(i,j) + sbdry(i,j) + 1.e-10))) exit
542 
543 
544  if (ds_it < 0.0) then ! Sbdry is now the upper bound.
545  if (sb_max_set .and. (sbdry(i,j) > sb_max)) &
546  call mom_error(fatal,"shelf_calc_flux: Irregular iteration for Sbdry (max).")
547  sb_max = sbdry(i,j) ; ds_max = ds_it ; sb_max_set = .true.
548  else ! Sbdry is now the lower bound.
549  if (sb_min_set .and. (sbdry(i,j) < sb_min)) &
550  call mom_error(fatal, &
551  "shelf_calc_flux: Irregular iteration for Sbdry (min).")
552  sb_min = sbdry(i,j) ; ds_min = ds_it ; sb_min_set = .true.
553  endif ! dS_it < 0.0
554 
555  if (sb_min_set .and. sb_max_set) then
556  ! Use the false position method for the next iteration.
557  sbdry(i,j) = sb_min + (sb_max-sb_min) * &
558  (ds_min / (ds_min - ds_max))
559  else
560  sbdry(i,j) = sbdry_it
561  endif ! Sb_min_set
562 
563  sbdry(i,j) = sbdry_it
564  endif ! CS%find_salt_root
565 
566  enddo !it1
567  ! Check for non-convergence and/or non-boundedness?
568 
569  else
570  ! In the 2-equation form, the mixed layer turbulent exchange velocity
571  ! is specified and large enough that the ocean salinity at the interface
572  ! is about the same as the boundary layer salinity.
573 
574  call calculate_tfreeze(state%sss(i,j), p_int(i), iss%tfreeze(i,j), cs%eqn_of_state)
575 
576  exch_vel_t(i,j) = cs%gamma_t
577  iss%tflux_ocn(i,j) = rhocp * exch_vel_t(i,j) * (state%sst(i,j) - iss%tfreeze(i,j))
578  iss%tflux_shelf(i,j) = 0.0
579  iss%water_flux(i,j) = i_lf * iss%tflux_ocn(i,j)
580  sbdry(i,j) = 0.0
581  endif
582  else !not shelf
583  iss%tflux_ocn(i,j) = 0.0
584  endif
585 
586 ! haline_driving(:,:) = state%sss(i,j) - Sbdry(i,j)
587 
588  enddo ! i-loop
589  enddo ! j-loop
590 
591  ! ISS%water_flux = net liquid water into the ocean ( kg/(m^2 s) )
592  ! We want melt in m/year
593  if (cs%const_gamma) then ! use ISOMIP+ eq. with rho_fw
594  fluxes%iceshelf_melt = iss%water_flux * (86400.0*365.0/rho_fw) * cs%flux_factor
595  else ! use original eq.
596  fluxes%iceshelf_melt = iss%water_flux * (86400.0*365.0/cs%density_ice) * cs%flux_factor
597  endif
598 
599  do j=js,je ; do i=is,ie
600  if ((idens*state%ocean_mass(i,j) > cs%col_thick_melt_threshold) .and. &
601  (iss%area_shelf_h(i,j) > 0.0) .and. &
602  (cs%isthermo) .and. (state%Hml(i,j) > 0.0) ) then
603 
604  ! Set melt to zero above a cutoff pressure
605  ! (CS%Rho0*CS%cutoff_depth*CS%g_Earth) this is needed for the isomip
606  ! test case.
607  if ((cs%g_Earth * iss%mass_shelf(i,j)) < cs%Rho0*cs%cutoff_depth* &
608  cs%g_Earth) then
609  iss%water_flux(i,j) = 0.0
610  fluxes%iceshelf_melt(i,j) = 0.0
611  endif
612  ! Compute haline driving, which is one of the diags. used in ISOMIP
613  haline_driving(i,j) = (iss%water_flux(i,j) * sbdry(i,j)) / &
614  (cs%Rho0 * exch_vel_s(i,j))
615 
616  !!!!!!!!!!!!!!!!!!!!!!!!!!!!Safety checks !!!!!!!!!!!!!!!!!!!!!!!!!
617  !1)Check if haline_driving computed above is consistent with
618  ! haline_driving = state%sss - Sbdry
619  !if (fluxes%iceshelf_melt(i,j) /= 0.0) then
620  ! if (haline_driving(i,j) /= (state%sss(i,j) - Sbdry(i,j))) then
621  ! write(mesg,*) 'at i,j=',i,j,' haline_driving, sss-Sbdry',haline_driving(i,j), &
622  ! (state%sss(i,j) - Sbdry(i,j))
623  ! call MOM_error(FATAL, &
624  ! "shelf_calc_flux: Inconsistency in melt and haline_driving"//trim(mesg))
625  ! endif
626  !endif
627 
628  ! 2) check if |melt| > 0 when ustar_shelf = 0.
629  ! this should never happen
630  if ((abs(fluxes%iceshelf_melt(i,j))>0.0) .and. (fluxes%ustar_shelf(i,j) == 0.0)) then
631  write(mesg,*) "|melt| = ",fluxes%iceshelf_melt(i,j)," > 0 and ustar_shelf = 0. at i,j", i, j
632  call mom_error(fatal, "shelf_calc_flux: "//trim(mesg))
633  endif
634  endif ! area_shelf_h
635  !!!!!!!!!!!!!!!!!!!!!!!!!!!!End of safety checks !!!!!!!!!!!!!!!!!!!
636  enddo ; enddo ! i- and j-loops
637 
638  ! mass flux [kg s-1], part of ISOMIP diags.
639  mass_flux(:,:) = 0.0
640  mass_flux(:,:) = iss%water_flux(:,:) * iss%area_shelf_h(:,:)
641 
642  if (cs%active_shelf_dynamics .or. cs%override_shelf_movement) then
643  call cpu_clock_begin(id_clock_pass)
644  call pass_var(iss%area_shelf_h, g%domain, complete=.false.)
645  call pass_var(iss%mass_shelf, g%domain)
646  call cpu_clock_end(id_clock_pass)
647  endif
648 
649  ! Melting has been computed, now is time to update thickness and mass
650  if ( cs%override_shelf_movement .and. (.not.cs%mass_from_file)) then
651  call change_thickness_using_melt(iss, g, time_step, fluxes, cs%rho_ice, cs%debug)
652 
653  if (cs%debug) then
654  call hchksum(iss%h_shelf, "h_shelf after change thickness using melt", g%HI, haloshift=0, scale=us%Z_to_m)
655  call hchksum(iss%mass_shelf, "mass_shelf after change thickness using melt", g%HI, haloshift=0)
656  endif
657  endif
658 
659  if (cs%debug) call mom_forcing_chksum("Before add shelf flux", fluxes, g, cs%US, haloshift=0)
660 
661  call add_shelf_flux(g, us, cs, state, fluxes)
662 
663  ! now the thermodynamic data is passed on... time to update the ice dynamic quantities
664 
665  if (cs%active_shelf_dynamics) then
666  update_ice_vel = .false.
667  coupled_gl = (cs%GL_couple .and. .not.cs%solo_ice_sheet)
668 
669  ! advect the ice shelf, and advance the front. Calving will be in here somewhere as well..
670  ! when we decide on how to do it
671  call update_ice_shelf(cs%dCS, iss, g, us, time_step, time, state%ocean_mass, coupled_gl)
672 
673  endif
674 
675  call enable_averaging(time_step,time,cs%diag)
676  if (cs%id_shelf_mass > 0) call post_data(cs%id_shelf_mass, iss%mass_shelf, cs%diag)
677  if (cs%id_area_shelf_h > 0) call post_data(cs%id_area_shelf_h, iss%area_shelf_h, cs%diag)
678  if (cs%id_ustar_shelf > 0) call post_data(cs%id_ustar_shelf, fluxes%ustar_shelf, cs%diag)
679  if (cs%id_melt > 0) call post_data(cs%id_melt, fluxes%iceshelf_melt, cs%diag)
680  if (cs%id_thermal_driving > 0) call post_data(cs%id_thermal_driving, (state%sst-iss%tfreeze), cs%diag)
681  if (cs%id_Sbdry > 0) call post_data(cs%id_Sbdry, sbdry, cs%diag)
682  if (cs%id_haline_driving > 0) call post_data(cs%id_haline_driving, haline_driving, cs%diag)
683  if (cs%id_mass_flux > 0) call post_data(cs%id_mass_flux, mass_flux, cs%diag)
684  if (cs%id_u_ml > 0) call post_data(cs%id_u_ml, state%u, cs%diag)
685  if (cs%id_v_ml > 0) call post_data(cs%id_v_ml, state%v, cs%diag)
686  if (cs%id_tfreeze > 0) call post_data(cs%id_tfreeze, iss%tfreeze, cs%diag)
687  if (cs%id_tfl_shelf > 0) call post_data(cs%id_tfl_shelf, iss%tflux_shelf, cs%diag)
688  if (cs%id_exch_vel_t > 0) call post_data(cs%id_exch_vel_t, exch_vel_t, cs%diag)
689  if (cs%id_exch_vel_s > 0) call post_data(cs%id_exch_vel_s, exch_vel_s, cs%diag)
690  if (cs%id_h_shelf > 0) call post_data(cs%id_h_shelf,iss%h_shelf,cs%diag)
691  if (cs%id_h_mask > 0) call post_data(cs%id_h_mask,iss%hmask,cs%diag)
692  call disable_averaging(cs%diag)
693 
694  if (present(forces)) then
695  call add_shelf_forces(g, us, cs, forces, do_shelf_area=(cs%active_shelf_dynamics .or. &
696  cs%override_shelf_movement))
697  endif
698 
699  call cpu_clock_end(id_clock_shelf)
700 
701  if (cs%debug) call mom_forcing_chksum("End of shelf calc flux", fluxes, g, cs%US, haloshift=0)
702 
703 end subroutine shelf_calc_flux
704 
705 !> Changes the thickness (mass) of the ice shelf based on sub-ice-shelf melting
706 subroutine change_thickness_using_melt(ISS, G, time_step, fluxes, rho_ice, debug)
707  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
708  type(ice_shelf_state), intent(inout) :: ISS !< A structure with elements that describe
709  !! the ice-shelf state
710  real, intent(in) :: time_step !< The time step for this update [s].
711  type(forcing), intent(inout) :: fluxes !< structure containing pointers to any possible
712  !! thermodynamic or mass-flux forcing fields.
713  real, intent(in) :: rho_ice !< The density of ice-shelf ice [kg m-2 Z-1 ~> kg m-3].
714  logical, optional, intent(in) :: debug !< If present and true, write chksums
715 
716  ! locals
717  real :: I_rho_ice
718  integer :: i, j
719 
720  i_rho_ice = 1.0 / rho_ice
721 
722  do j=g%jsc,g%jec ; do i=g%isc,g%iec
723  if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
724  ! first, zero out fluxes applied during previous time step
725  if (associated(fluxes%lprec)) fluxes%lprec(i,j) = 0.0
726  if (associated(fluxes%sens)) fluxes%sens(i,j) = 0.0
727  if (associated(fluxes%frac_shelf_h)) fluxes%frac_shelf_h(i,j) = 0.0
728  if (associated(fluxes%salt_flux)) fluxes%salt_flux(i,j) = 0.0
729 
730  if (iss%water_flux(i,j) / rho_ice * time_step < iss%h_shelf(i,j)) then
731  iss%h_shelf(i,j) = iss%h_shelf(i,j) - iss%water_flux(i,j) / rho_ice * time_step
732  else
733  ! the ice is about to melt away, so set thickness, area, and mask to zero
734  ! NOTE: this is not mass conservative should maybe scale salt & heat flux for this cell
735  iss%h_shelf(i,j) = 0.0
736  iss%hmask(i,j) = 0.0
737  iss%area_shelf_h(i,j) = 0.0
738  endif
739  endif
740  enddo ; enddo
741 
742  call pass_var(iss%area_shelf_h, g%domain)
743  call pass_var(iss%h_shelf, g%domain)
744  call pass_var(iss%hmask, g%domain)
745 
746  !### combine this with the loops above.
747  do j=g%jsd,g%jed ; do i=g%isd,g%ied
748  if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
749  iss%mass_shelf(i,j) = iss%h_shelf(i,j)*rho_ice
750  endif
751  enddo ; enddo
752 
753  call pass_var(iss%mass_shelf, g%domain)
754 
755 end subroutine change_thickness_using_melt
756 
757 !> This subroutine adds the mechanical forcing fields and perhaps shelf areas, based on
758 !! the ice state in ice_shelf_CS.
759 subroutine add_shelf_forces(G, US, CS, forces, do_shelf_area)
760  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
761  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
762  type(ice_shelf_cs), pointer :: cs !< This module's control structure.
763  type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
764  logical, optional, intent(in) :: do_shelf_area !< If true find the shelf-covered areas.
765 
766  real :: kv_rho_ice ! The viscosity of ice divided by its density [m5 kg-1 s-1].
767  real :: press_ice ! The pressure of the ice shelf per unit area of ocean (not ice) [Pa].
768  logical :: find_area ! If true find the shelf areas at u & v points.
769  type(ice_shelf_state), pointer :: iss => null() ! A structure with elements that describe
770  ! the ice-shelf state
771 
772  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
773  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
774  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
775 
776  if ((cs%grid%isc /= g%isc) .or. (cs%grid%iec /= g%iec) .or. &
777  (cs%grid%jsc /= g%jsc) .or. (cs%grid%jec /= g%jec)) &
778  call mom_error(fatal,"add_shelf_forces: Incompatible ocean and ice shelf grids.")
779 
780  iss => cs%ISS
781 
782  find_area = .true. ; if (present(do_shelf_area)) find_area = do_shelf_area
783 
784  if (find_area) then
785  ! The frac_shelf is set over the widest possible area. Could it be smaller?
786  do j=jsd,jed ; do i=isd,ied-1
787  forces%frac_shelf_u(i,j) = 0.0
788  if ((g%areaT(i,j) + g%areaT(i+1,j) > 0.0)) & ! .and. (G%areaCu(I,j) > 0.0)) &
789  forces%frac_shelf_u(i,j) = ((iss%area_shelf_h(i,j) + iss%area_shelf_h(i+1,j)) / &
790  (us%L_to_m**2*g%areaT(i,j) + us%L_to_m**2*g%areaT(i+1,j)))
791  enddo ; enddo
792  do j=jsd,jed-1 ; do i=isd,ied
793  forces%frac_shelf_v(i,j) = 0.0
794  if ((g%areaT(i,j) + g%areaT(i,j+1) > 0.0)) & ! .and. (G%areaCv(i,J) > 0.0)) &
795  forces%frac_shelf_v(i,j) = ((iss%area_shelf_h(i,j) + iss%area_shelf_h(i,j+1)) / &
796  (us%L_to_m**2*g%areaT(i,j) + us%L_to_m**2*g%areaT(i,j+1)))
797  enddo ; enddo
798  call pass_vector(forces%frac_shelf_u, forces%frac_shelf_v, g%domain, to_all, cgrid_ne)
799  endif
800 
801  !### Consider working over a smaller array range.
802  do j=jsd,jed ; do i=isd,ied
803  press_ice = (iss%area_shelf_h(i,j) * us%m_to_L**2*g%IareaT(i,j)) * (cs%g_Earth * iss%mass_shelf(i,j))
804  if (associated(forces%p_surf)) then
805  if (.not.forces%accumulate_p_surf) forces%p_surf(i,j) = 0.0
806  forces%p_surf(i,j) = forces%p_surf(i,j) + press_ice
807  endif
808  if (associated(forces%p_surf_full)) then
809  if (.not.forces%accumulate_p_surf) forces%p_surf_full(i,j) = 0.0
810  forces%p_surf_full(i,j) = forces%p_surf_full(i,j) + press_ice
811  endif
812  enddo ; enddo
813 
814  ! For various reasons, forces%rigidity_ice_[uv] is always updated here. Note
815  ! that it may have been zeroed out where IOB is translated to forces and
816  ! contributions from icebergs and the sea-ice pack added subsequently.
817  !### THE RIGIDITY SHOULD ALSO INCORPORATE AREAL-COVERAGE INFORMATION.
818  kv_rho_ice = cs%kv_ice / cs%density_ice
819  do j=js,je ; do i=is-1,ie
820  if (.not.forces%accumulate_rigidity) forces%rigidity_ice_u(i,j) = 0.0
821  forces%rigidity_ice_u(i,j) = forces%rigidity_ice_u(i,j) + &
822  kv_rho_ice * min(iss%mass_shelf(i,j), iss%mass_shelf(i+1,j))
823  enddo ; enddo
824  do j=js-1,je ; do i=is,ie
825  if (.not.forces%accumulate_rigidity) forces%rigidity_ice_v(i,j) = 0.0
826  forces%rigidity_ice_v(i,j) = forces%rigidity_ice_v(i,j) + &
827  kv_rho_ice * min(iss%mass_shelf(i,j), iss%mass_shelf(i,j+1))
828  enddo ; enddo
829 
830  if (cs%debug) then
831  call uvchksum("rigidity_ice_[uv]", forces%rigidity_ice_u, forces%rigidity_ice_v, &
832  g%HI, symmetric=.true.)
833  call uvchksum("frac_shelf_[uv]", forces%frac_shelf_u, forces%frac_shelf_v, &
834  g%HI, symmetric=.true.)
835  endif
836 
837 end subroutine add_shelf_forces
838 
839 !> This subroutine adds the ice shelf pressure to the fluxes type.
840 subroutine add_shelf_pressure(G, US, CS, fluxes)
841  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
842  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
843  type(ice_shelf_cs), intent(in) :: CS !< This module's control structure.
844  type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be updated.
845 
846  real :: press_ice !< The pressure of the ice shelf per unit area of ocean (not ice) [Pa].
847  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
848  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
849 
850  if ((cs%grid%isc /= g%isc) .or. (cs%grid%iec /= g%iec) .or. &
851  (cs%grid%jsc /= g%jsc) .or. (cs%grid%jec /= g%jec)) &
852  call mom_error(fatal,"add_shelf_pressure: Incompatible ocean and ice shelf grids.")
853 
854  do j=js,je ; do i=is,ie
855  press_ice = (cs%ISS%area_shelf_h(i,j) * us%m_to_L**2*g%IareaT(i,j)) * (cs%g_Earth * cs%ISS%mass_shelf(i,j))
856  if (associated(fluxes%p_surf)) then
857  if (.not.fluxes%accumulate_p_surf) fluxes%p_surf(i,j) = 0.0
858  fluxes%p_surf(i,j) = fluxes%p_surf(i,j) + press_ice
859  endif
860  if (associated(fluxes%p_surf_full)) then
861  if (.not.fluxes%accumulate_p_surf) fluxes%p_surf_full(i,j) = 0.0
862  fluxes%p_surf_full(i,j) = fluxes%p_surf_full(i,j) + press_ice
863  endif
864  enddo ; enddo
865 
866 end subroutine add_shelf_pressure
867 
868 !> Updates surface fluxes that are influenced by sub-ice-shelf melting
869 subroutine add_shelf_flux(G, US, CS, state, fluxes)
870  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
871  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
872  type(ice_shelf_cs), pointer :: cs !< This module's control structure.
873  type(surface), intent(inout) :: state!< Surface ocean state
874  type(forcing), intent(inout) :: fluxes !< A structure of surface fluxes that may be used/updated.
875 
876  ! local variables
877  real :: irho0 !< The inverse of the mean density [m3 kg-1].
878  real :: frac_area !< The fractional area covered by the ice shelf [nondim].
879  real :: shelf_mass0 !< Total ice shelf mass at previous time (Time-dt).
880  real :: shelf_mass1 !< Total ice shelf mass at current time (Time).
881  real :: delta_mass_shelf!< Change in ice shelf mass over one time step [kg s-1]
882  real :: taux2, tauy2 !< The squared surface stresses [Pa].
883  real :: press_ice !< The pressure of the ice shelf per unit area of ocean (not ice) [Pa].
884  real :: asu1, asu2 !< Ocean areas covered by ice shelves at neighboring u-
885  real :: asv1, asv2 !< and v-points [m2].
886  real :: fraz !< refreezing rate [kg m-2 s-1]
887  real :: mean_melt_flux !< spatial mean melt flux [kg s-1] or [kg m-2 s-1] at various points in the code.
888  real :: sponge_area !< total area of sponge region [m2]
889  real :: t0 !< The previous time (Time-dt) [s].
890  type(time_type) :: time0!< The previous time (Time-dt)
891  real, dimension(SZDI_(G),SZDJ_(G)) :: last_mass_shelf !< Ice shelf mass
892  !! at at previous time (Time-dt) [kg m-2]
893  real, dimension(SZDI_(G),SZDJ_(G)) :: last_h_shelf !< Ice shelf thickness [Z ~> m]
894  !! at at previous time (Time-dt)
895  real, dimension(SZDI_(G),SZDJ_(G)) :: last_hmask !< Ice shelf mask
896  !! at at previous time (Time-dt)
897  real, dimension(SZDI_(G),SZDJ_(G)) :: last_area_shelf_h !< Ice shelf area [m2]
898  !! at at previous time (Time-dt)
899  type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
900  !! the ice-shelf state
901 
902  real :: kv_rho_ice ! The viscosity of ice divided by its density [m5 kg-1 s-1]
903  real, parameter :: rho_fw = 1000.0 ! fresh water density
904  character(len=160) :: mesg ! The text of an error message
905  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
906  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
907  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
908 
909  if ((cs%grid%isc /= g%isc) .or. (cs%grid%iec /= g%iec) .or. &
910  (cs%grid%jsc /= g%jsc) .or. (cs%grid%jec /= g%jec)) &
911  call mom_error(fatal,"add_shelf_flux: Incompatible ocean and ice shelf grids.")
912 
913  iss => cs%ISS
914 
915  call add_shelf_pressure(g, us, cs, fluxes)
916 
917  ! Determine ustar and the square magnitude of the velocity in the
918  ! bottom boundary layer. Together these give the TKE source and
919  ! vertical decay scale.
920 
921  if (cs%debug) then
922  if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then
923  call uvchksum("tau[xy]_shelf", state%taux_shelf, state%tauy_shelf, &
924  g%HI, haloshift=0)
925  endif
926  endif
927 
928  if (allocated(state%taux_shelf) .and. allocated(state%tauy_shelf)) then
929  call pass_vector(state%taux_shelf, state%tauy_shelf, g%domain, to_all, cgrid_ne)
930  ! GMM: melting is computed using ustar_shelf (and not ustar), which has already
931  ! been passed, I so believe we do not need to update fluxes%ustar.
932 ! Irho0 = 1.0 / CS%Rho0
933 ! do j=js,je ; do i=is,ie ; if (fluxes%frac_shelf_h(i,j) > 0.0) then
934  ! ### THIS SHOULD BE AN AREA WEIGHTED AVERAGE OF THE ustar_shelf POINTS.
935  ! taux2 = 0.0 ; tauy2 = 0.0
936  ! asu1 = (ISS%area_shelf_h(i-1,j) + ISS%area_shelf_h(i,j))
937  ! asu2 = (ISS%area_shelf_h(i,j) + ISS%area_shelf_h(i+1,j))
938  ! asv1 = (ISS%area_shelf_h(i,j-1) + ISS%area_shelf_h(i,j))
939  ! asv2 = (ISS%area_shelf_h(i,j) + ISS%area_shelf_h(i,j+1))
940  ! if ((asu1 + asu2 > 0.0) .and. associated(state%taux_shelf)) &
941  ! taux2 = (asu1 * state%taux_shelf(I-1,j)**2 + &
942  ! asu2 * state%taux_shelf(I,j)**2 ) / (asu1 + asu2)
943  ! if ((asv1 + asv2 > 0.0) .and. associated(state%tauy_shelf)) &
944  ! tauy2 = (asv1 * state%tauy_shelf(i,J-1)**2 + &
945  ! asv2 * state%tauy_shelf(i,J)**2 ) / (asv1 + asv2)
946 
947  ! fluxes%ustar(i,j) = MAX(CS%ustar_bg, US%m_to_Z*US%T_to_s*sqrt(Irho0 * sqrt(taux2 + tauy2)))
948 ! endif ; enddo ; enddo
949  endif
950 
951  if (cs%active_shelf_dynamics .or. cs%override_shelf_movement) then
952  do j=jsd,jed ; do i=isd,ied
953  if (g%areaT(i,j) > 0.0) &
954  fluxes%frac_shelf_h(i,j) = iss%area_shelf_h(i,j) * us%m_to_L**2*g%IareaT(i,j)
955  enddo ; enddo
956  endif
957 
958  do j=js,je ; do i=is,ie ; if (iss%area_shelf_h(i,j) > 0.0) then
959  frac_area = fluxes%frac_shelf_h(i,j) !### Should this be 1-frac_shelf_h?
960  if (associated(fluxes%sw)) fluxes%sw(i,j) = 0.0
961  if (associated(fluxes%sw_vis_dir)) fluxes%sw_vis_dir(i,j) = 0.0
962  if (associated(fluxes%sw_vis_dif)) fluxes%sw_vis_dif(i,j) = 0.0
963  if (associated(fluxes%sw_nir_dir)) fluxes%sw_nir_dir(i,j) = 0.0
964  if (associated(fluxes%sw_nir_dif)) fluxes%sw_nir_dif(i,j) = 0.0
965  if (associated(fluxes%lw)) fluxes%lw(i,j) = 0.0
966  if (associated(fluxes%latent)) fluxes%latent(i,j) = 0.0
967  if (associated(fluxes%evap)) fluxes%evap(i,j) = 0.0
968  if (associated(fluxes%lprec)) then
969  if (iss%water_flux(i,j) > 0.0) then
970  fluxes%lprec(i,j) = us%kg_m3_to_R*us%m_to_Z*us%T_to_s*frac_area*iss%water_flux(i,j)*cs%flux_factor
971  else
972  fluxes%lprec(i,j) = 0.0
973  fluxes%evap(i,j) = us%kg_m3_to_R*us%m_to_Z*us%T_to_s*frac_area*iss%water_flux(i,j)*cs%flux_factor
974  endif
975  endif
976 
977  if (associated(fluxes%sens)) &
978  fluxes%sens(i,j) = -frac_area*iss%tflux_ocn(i,j)*cs%flux_factor
979  if (associated(fluxes%salt_flux)) &
980  fluxes%salt_flux(i,j) = frac_area * iss%salt_flux(i,j)*cs%flux_factor
981  endif ; enddo ; enddo
982 
983  ! keep sea level constant by removing mass in the sponge
984  ! region (via virtual precip, vprec). Apply additional
985  ! salt/heat fluxes so that the resultant surface buoyancy
986  ! forcing is ~ 0.
987  ! This is needed for some of the ISOMIP+ experiments.
988 
989  if (cs%constant_sea_level) then
990  !### This code has lots of problems with hard coded constants and the use of
991  !### of non-reproducing sums. It needs to be refactored. -RWH
992 
993  if (.not. associated(fluxes%salt_flux)) allocate(fluxes%salt_flux(ie,je))
994  if (.not. associated(fluxes%vprec)) allocate(fluxes%vprec(ie,je))
995  fluxes%salt_flux(:,:) = 0.0 ; fluxes%vprec(:,:) = 0.0
996 
997  mean_melt_flux = 0.0; sponge_area = 0.0
998  do j=js,je ; do i=is,ie
999  frac_area = fluxes%frac_shelf_h(i,j)
1000  if (frac_area > 0.0) &
1001  mean_melt_flux = mean_melt_flux + (iss%water_flux(i,j)) * iss%area_shelf_h(i,j)
1002 
1003  !### These hard-coded limits need to be corrected. They are inappropriate here.
1004  if (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0) then
1005  sponge_area = sponge_area + us%L_to_m**2*g%areaT(i,j)
1006  endif
1007  enddo ; enddo
1008 
1009  ! take into account changes in mass (or thickness) when imposing ice shelf mass
1010  if (cs%override_shelf_movement .and. cs%mass_from_file) then
1011  t0 = time_type_to_real(cs%Time) - cs%time_step
1012 
1013  ! just compute changes in mass after first time step
1014  if (t0>0.0) then
1015  time0 = real_to_time(t0)
1016  last_hmask(:,:) = iss%hmask(:,:) ; last_area_shelf_h(:,:) = iss%area_shelf_h(:,:)
1017  call time_interp_external(cs%id_read_mass, time0, last_mass_shelf)
1018  last_h_shelf(:,:) = last_mass_shelf(:,:) / cs%rho_ice
1019 
1020  ! apply calving
1021  if (cs%min_thickness_simple_calve > 0.0) then
1022  call ice_shelf_min_thickness_calve(g, last_h_shelf, last_area_shelf_h, last_hmask, &
1023  cs%min_thickness_simple_calve)
1024  ! convert to mass again
1025  last_mass_shelf(:,:) = last_h_shelf(:,:) * cs%rho_ice
1026  endif
1027 
1028  shelf_mass0 = 0.0; shelf_mass1 = 0.0
1029  ! get total ice shelf mass at (Time-dt) and (Time), in kg
1030  do j=js,je ; do i=is,ie
1031  ! just floating shelf (0.1 is a threshold for min ocean thickness)
1032  if (((1.0/cs%density_ocean_avg)*state%ocean_mass(i,j) > 0.1) .and. &
1033  (iss%area_shelf_h(i,j) > 0.0)) then
1034  shelf_mass0 = shelf_mass0 + (last_mass_shelf(i,j) * iss%area_shelf_h(i,j))
1035  shelf_mass1 = shelf_mass1 + (iss%mass_shelf(i,j) * iss%area_shelf_h(i,j))
1036  endif
1037  enddo ; enddo
1038  call sum_across_pes(shelf_mass0); call sum_across_pes(shelf_mass1)
1039  delta_mass_shelf = (shelf_mass1 - shelf_mass0)/cs%time_step
1040 ! delta_mass_shelf = (shelf_mass1 - shelf_mass0)* &
1041 ! (rho_fw/CS%density_ice)/CS%time_step
1042 ! write(mesg,*)'delta_mass_shelf = ',delta_mass_shelf
1043 ! call MOM_mesg(mesg,5)
1044  else! first time step
1045  delta_mass_shelf = 0.0
1046  endif
1047  else ! ice shelf mass does not change
1048  delta_mass_shelf = 0.0
1049  endif
1050 
1051  call sum_across_pes(mean_melt_flux)
1052  call sum_across_pes(sponge_area)
1053 
1054  ! average total melt flux over sponge area
1055  mean_melt_flux = (mean_melt_flux+delta_mass_shelf) / sponge_area !kg/(m^2 s)
1056 
1057  ! apply fluxes
1058  do j=js,je ; do i=is,ie
1059  ! Note the following is hard coded for ISOMIP
1060  if (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0) then
1061  fluxes%vprec(i,j) = -mean_melt_flux * cs%density_ice/1000. ! evap is negative
1062  fluxes%sens(i,j) = fluxes%vprec(i,j) * cs%Cp * cs%T0 ! W /m^2
1063  ! Rescale fluxes%vprec to the proper units.
1064  fluxes%vprec(i,j) = us%kg_m3_to_R*us%m_to_Z*us%T_to_s * fluxes%vprec(i,j)
1065  fluxes%salt_flux(i,j) = fluxes%vprec(i,j) * cs%S0*1.0e-3 ! kg (salt)/(m^2 s)
1066  endif
1067  enddo ; enddo
1068 
1069  if (cs%debug) then
1070  write(mesg,*) 'Mean melt flux (kg/(m^2 s)), dt = ', mean_melt_flux, cs%time_step
1071  call mom_mesg(mesg)
1072  call mom_forcing_chksum("After constant sea level", fluxes, g, cs%US, haloshift=0)
1073  endif
1074 
1075  endif !constant_sea_level
1076 
1077 end subroutine add_shelf_flux
1078 
1079 
1080 !> Initializes shelf model data, parameters and diagnostics
1081 subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fluxes, Time_in, solo_ice_sheet_in)
1082  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1083  type(ocean_grid_type), pointer :: ocn_grid !< The calling ocean model's horizontal grid structure
1084  type(time_type), intent(inout) :: time !< The clock that that will indicate the model time
1085  type(ice_shelf_cs), pointer :: cs !< A pointer to the ice shelf control structure
1086  type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate the diagnostic output.
1087  type(forcing), optional, intent(inout) :: fluxes !< A structure containing pointers to any possible
1088  !! thermodynamic or mass-flux forcing fields.
1089  type(mech_forcing), optional, intent(inout) :: forces !< A structure with the driving mechanical forces
1090  type(time_type), optional, intent(in) :: time_in !< The time at initialization.
1091  logical, optional, intent(in) :: solo_ice_sheet_in !< If present, this indicates whether
1092  !! a solo ice-sheet driver.
1093 
1094  type(ocean_grid_type), pointer :: g => null(), og => null() ! Pointers to grids for convenience.
1095  type(unit_scale_type), pointer :: us => null() ! Pointer to a structure containing
1096  ! various unit conversion factors
1097  type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
1098  !! the ice-shelf state
1099  type(directories) :: dirs
1100  type(dyn_horgrid_type), pointer :: dg => null()
1101  real :: z_rescale ! A rescaling factor for heights from the representation in
1102  ! a reastart fole to the internal representation in this run.
1103  real :: cdrag, drag_bg_vel
1104  logical :: new_sim, save_ic, var_force
1105  !This include declares and sets the variable "version".
1106 #include "version_variable.h"
1107  character(len=200) :: config
1108  character(len=200) :: ic_file,filename,inputdir
1109  character(len=40) :: mdl = "MOM_ice_shelf" ! This module's name.
1110  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed, isdq, iedq, jsdq, jedq
1111  integer :: wd_halos(2)
1112  logical :: read_tideamp, shelf_mass_is_dynamic, debug
1113  character(len=240) :: tideamp_file
1114  real :: utide
1115  if (associated(cs)) then
1116  call mom_error(fatal, "MOM_ice_shelf.F90, initialize_ice_shelf: "// &
1117  "called with an associated control structure.")
1118  return
1119  endif
1120  allocate(cs)
1121 
1122  ! Go through all of the infrastructure initialization calls, since this is
1123  ! being treated as an independent component that just happens to use the
1124  ! MOM's grid and infrastructure.
1125  call get_mom_input(dirs=dirs)
1126 
1127  ! Determining the internal unit scaling factors for this run.
1128  call unit_scaling_init(param_file, cs%US)
1129 
1130  ! Set up the ice-shelf domain and grid
1131  wd_halos(:)=0
1132  call mom_domains_init(cs%grid%domain, param_file, min_halo=wd_halos, symmetric=grid_sym_)
1133  ! call diag_mediator_init(CS%grid,param_file,CS%diag)
1134  ! this needs to be fixed - will probably break when not using coupled driver 0
1135  call mom_grid_init(cs%grid, param_file, cs%US)
1136 
1137  call create_dyn_horgrid(dg, cs%grid%HI)
1138  call clone_mom_domain(cs%grid%Domain, dg%Domain)
1139 
1140  call set_grid_metrics(dg, param_file, cs%US)
1141  ! call set_diag_mediator_grid(CS%grid, CS%diag)
1142 
1143  ! The ocean grid possibly uses different symmetry.
1144  if (associated(ocn_grid)) then ; cs%ocn_grid => ocn_grid
1145  else ; cs%ocn_grid => cs%grid ; endif
1146 
1147  ! Convenience pointers
1148  g => cs%grid
1149  og => cs%ocn_grid
1150  us => cs%US
1151 
1152  if (is_root_pe()) then
1153  write(0,*) 'OG: ', og%isd, og%isc, og%iec, og%ied, og%jsd, og%jsc, og%jsd, og%jed
1154  write(0,*) 'IG: ', g%isd, g%isc, g%iec, g%ied, g%jsd, g%jsc, g%jsd, g%jed
1155  endif
1156 
1157  cs%Time = time ! ### This might not be in the right place?
1158  cs%diag => diag
1159 
1160  ! Are we being called from the solo ice-sheet driver? When called by the ocean
1161  ! model solo_ice_sheet_in is not preset.
1162  cs%solo_ice_sheet = .false.
1163  if (present(solo_ice_sheet_in)) cs%solo_ice_sheet = solo_ice_sheet_in
1164 
1165  if (present(time_in)) time = time_in
1166 
1167  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1168  isd = g%isd ; jsd = g%jsd ; ied = g%ied ; jed = g%jed
1169  isdq = g%IsdB ; iedq = g%IedB ; jsdq = g%JsdB ; jedq = g%JedB
1170 
1171  cs%Lat_fusion = 3.34e5
1172  cs%override_shelf_movement = .false. ; cs%active_shelf_dynamics = .false.
1173 
1174  call log_version(param_file, mdl, version, "")
1175  call get_param(param_file, mdl, "DEBUG", debug, default=.false.)
1176  call get_param(param_file, mdl, "DEBUG_IS", cs%debug, &
1177  "If true, write verbose debugging messages for the ice shelf.", &
1178  default=debug)
1179  call get_param(param_file, mdl, "DYNAMIC_SHELF_MASS", shelf_mass_is_dynamic, &
1180  "If true, the ice sheet mass can evolve with time.", &
1181  default=.false.)
1182  if (shelf_mass_is_dynamic) then
1183  call get_param(param_file, mdl, "OVERRIDE_SHELF_MOVEMENT", cs%override_shelf_movement, &
1184  "If true, user provided code specifies the ice-shelf "//&
1185  "movement instead of the dynamic ice model.", default=.false.)
1186  cs%active_shelf_dynamics = .not.cs%override_shelf_movement
1187  call get_param(param_file, mdl, "GROUNDING_LINE_INTERPOLATE", cs%GL_regularize, &
1188  "If true, regularize the floatation condition at the "//&
1189  "grounding line as in Goldberg Holland Schoof 2009.", default=.false.)
1190  call get_param(param_file, mdl, "GROUNDING_LINE_COUPLE", cs%GL_couple, &
1191  "If true, let the floatation condition be determined by "//&
1192  "ocean column thickness. This means that update_OD_ffrac "//&
1193  "will be called. GL_REGULARIZE and GL_COUPLE are exclusive.", &
1194  default=.false., do_not_log=cs%GL_regularize)
1195  if (cs%GL_regularize) cs%GL_couple = .false.
1196  endif
1197 
1198  call get_param(param_file, mdl, "SHELF_THERMO", cs%isthermo, &
1199  "If true, use a thermodynamically interactive ice shelf.", &
1200  default=.false.)
1201  call get_param(param_file, mdl, "SHELF_THREE_EQN", cs%threeeq, &
1202  "If true, use the three equation expression of "//&
1203  "consistency to calculate the fluxes at the ice-ocean "//&
1204  "interface.", default=.true.)
1205  call get_param(param_file, mdl, "SHELF_INSULATOR", cs%insulator, &
1206  "If true, the ice shelf is a perfect insulatior "//&
1207  "(no conduction).", default=.false.)
1208  call get_param(param_file, mdl, "MELTING_CUTOFF_DEPTH", cs%cutoff_depth, &
1209  "Depth above which the melt is set to zero (it must be >= 0) "//&
1210  "Default value won't affect the solution.", default=0.0)
1211  if (cs%cutoff_depth < 0.) &
1212  call mom_error(warning,"Initialize_ice_shelf: MELTING_CUTOFF_DEPTH must be >= 0.")
1213 
1214  call get_param(param_file, mdl, "CONST_SEA_LEVEL", cs%constant_sea_level, &
1215  "If true, apply evaporative, heat and salt fluxes in "//&
1216  "the sponge region. This will avoid a large increase "//&
1217  "in sea level. This option is needed for some of the "//&
1218  "ISOMIP+ experiments (Ocean3 and Ocean4). "//&
1219  "IMPORTANT: it is not currently possible to do "//&
1220  "prefect restarts using this flag.", default=.false.)
1221 
1222  call get_param(param_file, mdl, "ISOMIP_S_SUR_SPONGE", &
1223  cs%S0, "Surface salinity in the resoring region.", &
1224  default=33.8, do_not_log=.true.)
1225 
1226  call get_param(param_file, mdl, "ISOMIP_T_SUR_SPONGE", &
1227  cs%T0, "Surface temperature in the resoring region.", &
1228  default=-1.9, do_not_log=.true.)
1229 
1230  call get_param(param_file, mdl, "SHELF_3EQ_GAMMA", cs%const_gamma, &
1231  "If true, user specifies a constant nondimensional heat-transfer coefficient "//&
1232  "(GAMMA_T_3EQ), from which the salt-transfer coefficient is then computed "//&
1233  " as GAMMA_T_3EQ/35. This is used with SHELF_THREE_EQN.", default=.false.)
1234  if (cs%const_gamma) call get_param(param_file, mdl, "SHELF_3EQ_GAMMA_T", cs%Gamma_T_3EQ, &
1235  "Nondimensional heat-transfer coefficient.",default=2.2e-2, &
1236  units="nondim.", fail_if_missing=.true.)
1237 
1238  call get_param(param_file, mdl, "ICE_SHELF_MASS_FROM_FILE", &
1239  cs%mass_from_file, "Read the mass of the "//&
1240  "ice shelf (every time step) from a file.", default=.false.)
1241 
1242  if (cs%threeeq) &
1243  call get_param(param_file, mdl, "SHELF_S_ROOT", cs%find_salt_root, &
1244  "If SHELF_S_ROOT = True, salinity at the ice/ocean interface (Sbdry) "//&
1245  "is computed from a quadratic equation. Otherwise, the previous "//&
1246  "interactive method to estimate Sbdry is used.", default=.false.)
1247  if (cs%find_salt_root) then ! read liquidus coeffs.
1248  call get_param(param_file, mdl, "TFREEZE_S0_P0",cs%lambda1, &
1249  "this is the freezing potential temperature at "//&
1250  "S=0, P=0.", units="degC", default=0.0, do_not_log=.true.)
1251  call get_param(param_file, mdl, "DTFREEZE_DS",cs%lambda1, &
1252  "this is the derivative of the freezing potential "//&
1253  "temperature with salinity.", &
1254  units="degC psu-1", default=-0.054, do_not_log=.true.)
1255  call get_param(param_file, mdl, "DTFREEZE_DP",cs%lambda3, &
1256  "this is the derivative of the freezing potential "//&
1257  "temperature with pressure.", &
1258  units="degC Pa-1", default=0.0, do_not_log=.true.)
1259 
1260  endif
1261 
1262  if (.not.cs%threeeq) &
1263  call get_param(param_file, mdl, "SHELF_2EQ_GAMMA_T", cs%gamma_t, &
1264  "If SHELF_THREE_EQN is false, this the fixed turbulent "//&
1265  "exchange velocity at the ice-ocean interface.", &
1266  units="m s-1", fail_if_missing=.true.)
1267 
1268  call get_param(param_file, mdl, "G_EARTH", cs%g_Earth, &
1269  "The gravitational acceleration of the Earth.", &
1270  units="m s-2", default = 9.80)
1271  call get_param(param_file, mdl, "C_P", cs%Cp, &
1272  "The heat capacity of sea water.", units="J kg-1 K-1", &
1273  fail_if_missing=.true.)
1274  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
1275  "The mean ocean density used with BOUSSINESQ true to "//&
1276  "calculate accelerations and the mass for conservation "//&
1277  "properties, or with BOUSSINSEQ false to convert some "//&
1278  "parameters from vertical units of m to kg m-2.", &
1279  units="kg m-3", default=1035.0) !### MAKE THIS A SEPARATE PARAMETER.
1280  call get_param(param_file, mdl, "C_P_ICE", cs%Cp_ice, &
1281  "The heat capacity of ice.", units="J kg-1 K-1", &
1282  default=2.10e3)
1283 
1284  call get_param(param_file, mdl, "ICE_SHELF_FLUX_FACTOR", cs%flux_factor, &
1285  "Non-dimensional factor applied to shelf thermodynamic "//&
1286  "fluxes.", units="none", default=1.0)
1287 
1288  call get_param(param_file, mdl, "KV_ICE", cs%kv_ice, &
1289  "The viscosity of the ice.", units="m2 s-1", default=1.0e10)
1290  call get_param(param_file, mdl, "KV_MOLECULAR", cs%kv_molec, &
1291  "The molecular kinimatic viscosity of sea water at the "//&
1292  "freezing temperature.", units="m2 s-1", default=1.95e-6)
1293  call get_param(param_file, mdl, "ICE_SHELF_SALINITY", cs%Salin_ice, &
1294  "The salinity of the ice inside the ice shelf.", units="psu", &
1295  default=0.0)
1296  call get_param(param_file, mdl, "ICE_SHELF_TEMPERATURE", cs%Temp_ice, &
1297  "The temperature at the center of the ice shelf.", &
1298  units = "degC", default=-15.0)
1299  call get_param(param_file, mdl, "KD_SALT_MOLECULAR", cs%kd_molec_salt, &
1300  "The molecular diffusivity of salt in sea water at the "//&
1301  "freezing point.", units="m2 s-1", default=8.02e-10)
1302  call get_param(param_file, mdl, "KD_TEMP_MOLECULAR", cs%kd_molec_temp, &
1303  "The molecular diffusivity of heat in sea water at the "//&
1304  "freezing point.", units="m2 s-1", default=1.41e-7)
1305  call get_param(param_file, mdl, "RHO_0", cs%density_ocean_avg, &
1306  "avg ocean density used in floatation cond", &
1307  units="kg m-3", default=1035.)
1308  call get_param(param_file, mdl, "DT_FORCING", cs%time_step, &
1309  "The time step for changing forcing, coupling with other "//&
1310  "components, or potentially writing certain diagnostics. "//&
1311  "The default value is given by DT.", units="s", default=0.0)
1312 
1313  call get_param(param_file, mdl, "COL_THICK_MELT_THRESHOLD", cs%col_thick_melt_threshold, &
1314  "The minimum ML thickness where melting is allowed.", units="m", &
1315  default=0.0)
1316 
1317  call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, &
1318  "If true, read a file (given by TIDEAMP_FILE) containing "//&
1319  "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
1320 
1321  call safe_alloc_ptr(cs%utide,isd,ied,jsd,jed) ; cs%utide(:,:) = 0.0
1322 
1323  if (read_tideamp) then
1324  call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
1325  "The path to the file containing the spatially varying "//&
1326  "tidal amplitudes.", &
1327  default="tideamp.nc")
1328  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1329  inputdir = slasher(inputdir)
1330  tideamp_file = trim(inputdir) // trim(tideamp_file)
1331  call mom_read_data(tideamp_file,'tideamp',cs%utide,g%domain,timelevel=1)
1332  else
1333  call get_param(param_file, mdl, "UTIDE", utide, &
1334  "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
1335  units="m s-1", default=0.0)
1336  cs%utide(:,:) = utide
1337  endif
1338 
1339  call eos_init(param_file, cs%eqn_of_state)
1340 
1341  !! new parameters that need to be in MOM_input
1342 
1343  if (cs%active_shelf_dynamics) then
1344 
1345  call get_param(param_file, mdl, "DENSITY_ICE", cs%density_ice, &
1346  "A typical density of ice.", units="kg m-3", default=917.0)
1347 
1348  call get_param(param_file, mdl, "INPUT_FLUX_ICE_SHELF", cs%input_flux, &
1349  "volume flux at upstream boundary", units="m2 s-1", default=0.)
1350  call get_param(param_file, mdl, "INPUT_THICK_ICE_SHELF", cs%input_thickness, &
1351  "flux thickness at upstream boundary", units="m", default=1000.)
1352  else
1353  ! This is here because of inconsistent defaults. I don't know why. RWH
1354  call get_param(param_file, mdl, "DENSITY_ICE", cs%density_ice, &
1355  "A typical density of ice.", units="kg m-3", default=900.0)
1356  endif
1357  cs%rho_ice = cs%density_ice*us%Z_to_m
1358  call get_param(param_file, mdl, "MIN_THICKNESS_SIMPLE_CALVE", &
1359  cs%min_thickness_simple_calve, &
1360  "Min thickness rule for the very simple calving law",&
1361  units="m", default=0.0, scale=us%m_to_Z)
1362 
1363  call get_param(param_file, mdl, "USTAR_SHELF_BG", cs%ustar_bg, &
1364  "The minimum value of ustar under ice sheves.", &
1365  units="m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
1366  call get_param(param_file, mdl, "CDRAG_SHELF", cdrag, &
1367  "CDRAG is the drag coefficient relating the magnitude of "//&
1368  "the velocity field to the surface stress.", units="nondim", &
1369  default=0.003)
1370  cs%cdrag = cdrag
1371  if (cs%ustar_bg <= 0.0) then
1372  call get_param(param_file, mdl, "DRAG_BG_VEL_SHELF", drag_bg_vel, &
1373  "DRAG_BG_VEL is either the assumed bottom velocity (with "//&
1374  "LINEAR_DRAG) or an unresolved velocity that is "//&
1375  "combined with the resolved velocity to estimate the "//&
1376  "velocity magnitude.", units="m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
1377  if (cs%cdrag*drag_bg_vel > 0.0) cs%ustar_bg = sqrt(cs%cdrag)*drag_bg_vel
1378  endif
1379 
1380  ! Allocate and initialize state variables to default values
1381  call ice_shelf_state_init(cs%ISS, cs%grid)
1382  iss => cs%ISS
1383 
1384  ! Allocate the arrays for passing ice-shelf data through the forcing type.
1385  if (.not. cs%solo_ice_sheet) then
1386  call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes.")
1387  ! GMM: the following assures that water/heat fluxes are just allocated
1388  ! when SHELF_THERMO = True. These fluxes are necessary if one wants to
1389  ! use either ENERGETICS_SFC_PBL (ALE mode) or BULKMIXEDLAYER (layer mode).
1390  if (present(fluxes)) &
1391  call allocate_forcing_type(cs%ocn_grid, fluxes, ustar=.true., shelf=.true., &
1392  press=.true., water=cs%isthermo, heat=cs%isthermo)
1393  if (present(forces)) &
1394  call allocate_mech_forcing(cs%ocn_grid, forces, ustar=.true., shelf=.true., press=.true.)
1395  else
1396  call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: allocating fluxes in solo mode.")
1397  if (present(fluxes)) &
1398  call allocate_forcing_type(g, fluxes, ustar=.true., shelf=.true., press=.true.)
1399  if (present(forces)) &
1400  call allocate_mech_forcing(g, forces, ustar=.true., shelf=.true., press=.true.)
1401  endif
1402 
1403  ! Set up the bottom depth, G%D either analytically or from file
1404  call mom_initialize_topography(dg%bathyT, g%max_depth, dg, param_file)
1405  call rescale_dyn_horgrid_bathymetry(dg, us%Z_to_m)
1406 
1407  ! Set up the Coriolis parameter, G%f, usually analytically.
1408  call mom_initialize_rotation(dg%CoriolisBu, dg, param_file, us)
1409  ! This copies grid elements, including bathyT and CoriolisBu from dG to CS%grid.
1410  call copy_dyngrid_to_mom_grid(dg, cs%grid, us)
1411 
1412  call destroy_dyn_horgrid(dg)
1413 
1414  ! Set up the restarts.
1415  call restart_init(param_file, cs%restart_CSp, "Shelf.res")
1416  call register_restart_field(iss%mass_shelf, "shelf_mass", .true., cs%restart_CSp, &
1417  "Ice shelf mass", "kg m-2")
1418  call register_restart_field(iss%area_shelf_h, "shelf_area", .true., cs%restart_CSp, &
1419  "Ice shelf area in cell", "m2")
1420  call register_restart_field(iss%h_shelf, "h_shelf", .true., cs%restart_CSp, &
1421  "ice sheet/shelf thickness", "m")
1422  call register_restart_field(us%m_to_Z_restart, "m_to_Z", .false., cs%restart_CSp, &
1423  "Height unit conversion factor", "Z meter-1")
1424  if (cs%active_shelf_dynamics) then
1425  call register_restart_field(iss%hmask, "h_mask", .true., cs%restart_CSp, &
1426  "ice sheet/shelf thickness mask" ,"none")
1427  endif
1428 
1429  ! if (CS%active_shelf_dynamics) then !### Consider adding an ice shelf dynamics switch.
1430  ! Allocate CS%dCS and specify additional restarts for ice shelf dynamics
1431  call register_ice_shelf_dyn_restarts(g, param_file, cs%dCS, cs%restart_CSp)
1432  ! endif
1433 
1434  !GMM - I think we do not need to save ustar_shelf and iceshelf_melt in the restart file
1435  !if (.not. CS%solo_ice_sheet) then
1436  ! call register_restart_field(fluxes%ustar_shelf, "ustar_shelf", .false., CS%restart_CSp, &
1437  ! "Friction velocity under ice shelves", "m s-1")
1438  ! call register_restart_field(fluxes%iceshelf_melt, "iceshelf_melt", .false., CS%restart_CSp, &
1439  ! "Ice Shelf Melt Rate", "m year-1")
1440  !endif
1441 
1442  cs%restart_output_dir = dirs%restart_output_dir
1443 
1444  new_sim = .false.
1445  if ((dirs%input_filename(1:1) == 'n') .and. &
1446  (len_trim(dirs%input_filename) == 1)) new_sim = .true.
1447 
1448  if (cs%override_shelf_movement .and. cs%mass_from_file) then
1449 
1450  ! initialize the ids for reading shelf mass from a netCDF
1451  call initialize_shelf_mass(g, param_file, cs, iss)
1452 
1453  if (new_sim) then
1454  ! new simulation, initialize ice thickness as in the static case
1455  call initialize_ice_thickness(iss%h_shelf, iss%area_shelf_h, iss%hmask, g, us, param_file)
1456 
1457  ! next make sure mass is consistent with thickness
1458  do j=g%jsd,g%jed ; do i=g%isd,g%ied
1459  if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
1460  iss%mass_shelf(i,j) = iss%h_shelf(i,j)*cs%rho_ice
1461  endif
1462  enddo ; enddo
1463 
1464  if (cs%min_thickness_simple_calve > 0.0) &
1465  call ice_shelf_min_thickness_calve(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, &
1466  cs%min_thickness_simple_calve)
1467  endif
1468  endif
1469 
1470  if (cs%active_shelf_dynamics) then
1471  ! the only reason to initialize boundary conds is if the shelf is dynamic - MJH
1472 
1473  ! call initialize_ice_shelf_boundary ( CS%u_face_mask_bdry, CS%v_face_mask_bdry, &
1474  ! CS%u_flux_bdry_val, CS%v_flux_bdry_val, &
1475  ! CS%u_bdry_val, CS%v_bdry_val, CS%h_bdry_val, &
1476  ! ISS%hmask, G, param_file)
1477 
1478  endif
1479 
1480  if (new_sim .and. (.not. (cs%override_shelf_movement .and. cs%mass_from_file))) then
1481 
1482  ! This model is initialized internally or from a file.
1483  call initialize_ice_thickness(iss%h_shelf, iss%area_shelf_h, iss%hmask, g, us, param_file)
1484 
1485  ! next make sure mass is consistent with thickness
1486  do j=g%jsd,g%jed ; do i=g%isd,g%ied
1487  if ((iss%hmask(i,j) == 1) .or. (iss%hmask(i,j) == 2)) then
1488  iss%mass_shelf(i,j) = iss%h_shelf(i,j)*cs%rho_ice
1489  endif
1490  enddo ; enddo
1491 
1492  ! else ! Previous block for new_sim=.T., this block restores the state.
1493  elseif (.not.new_sim) then
1494  ! This line calls a subroutine that reads the initial conditions from a restart file.
1495  call mom_mesg("MOM_ice_shelf.F90, initialize_ice_shelf: Restoring ice shelf from file.")
1496  call restore_state(dirs%input_filename, dirs%restart_input_dir, time, &
1497  g, cs%restart_CSp)
1498 
1499  if ((us%m_to_Z_restart /= 0.0) .and. (us%m_to_Z_restart /= us%m_to_Z)) then
1500  z_rescale = us%m_to_Z / us%m_to_Z_restart
1501  do j=g%jsc,g%jec ; do i=g%isc,g%iec
1502  iss%h_shelf(i,j) = z_rescale * iss%h_shelf(i,j)
1503  enddo ; enddo
1504  endif
1505 
1506  endif ! .not. new_sim
1507 
1508  cs%Time = time
1509 
1510  call cpu_clock_begin(id_clock_pass)
1511  call pass_var(iss%area_shelf_h, g%domain)
1512  call pass_var(iss%h_shelf, g%domain)
1513  call pass_var(iss%mass_shelf, g%domain)
1514  call pass_var(iss%hmask, g%domain)
1515  call pass_var(g%bathyT, g%domain)
1516  call cpu_clock_end(id_clock_pass)
1517 
1518  do j=jsd,jed ; do i=isd,ied
1519  if (iss%area_shelf_h(i,j) > us%L_to_m**2*g%areaT(i,j)) then
1520  call mom_error(warning,"Initialize_ice_shelf: area_shelf_h exceeds G%areaT.")
1521  iss%area_shelf_h(i,j) = us%L_to_m**2*g%areaT(i,j)
1522  endif
1523  enddo ; enddo
1524  if (present(fluxes)) then ; do j=jsd,jed ; do i=isd,ied
1525  if (g%areaT(i,j) > 0.0) fluxes%frac_shelf_h(i,j) = iss%area_shelf_h(i,j) / (us%L_to_m**2*g%areaT(i,j))
1526  enddo ; enddo ; endif
1527 
1528  if (cs%debug) then
1529  call hchksum(fluxes%frac_shelf_h, "IS init: frac_shelf_h", g%HI, haloshift=0)
1530  endif
1531 
1532  if (present(forces)) &
1533  call add_shelf_forces(g, us, cs, forces, do_shelf_area=.not.cs%solo_ice_sheet)
1534 
1535  if (present(fluxes)) call add_shelf_pressure(g, us, cs, fluxes)
1536 
1537  if (cs%active_shelf_dynamics .and. .not.cs%isthermo) then
1538  iss%water_flux(:,:) = 0.0
1539  endif
1540 
1541  if (shelf_mass_is_dynamic) &
1542  call initialize_ice_shelf_dyn(param_file, time, iss, cs%dCS, g, us, diag, new_sim, solo_ice_sheet_in)
1543 
1544  call fix_restart_unit_scaling(us)
1545 
1546  call get_param(param_file, mdl, "SAVE_INITIAL_CONDS", save_ic, &
1547  "If true, save the ice shelf initial conditions.", &
1548  default=.false.)
1549  if (save_ic) call get_param(param_file, mdl, "SHELF_IC_OUTPUT_FILE", ic_file,&
1550  "The name-root of the output file for the ice shelf "//&
1551  "initial conditions.", default="MOM_Shelf_IC")
1552 
1553  if (save_ic .and. .not.((dirs%input_filename(1:1) == 'r') .and. &
1554  (len_trim(dirs%input_filename) == 1))) then
1555  call save_restart(dirs%output_directory, cs%Time, g, &
1556  cs%restart_CSp, filename=ic_file)
1557  endif
1558 
1559 
1560  cs%id_area_shelf_h = register_diag_field('ocean_model', 'area_shelf_h', cs%diag%axesT1, cs%Time, &
1561  'Ice Shelf Area in cell', 'meter-2')
1562  cs%id_shelf_mass = register_diag_field('ocean_model', 'shelf_mass', cs%diag%axesT1, cs%Time, &
1563  'mass of shelf', 'kg/m^2')
1564  cs%id_h_shelf = register_diag_field('ocean_model', 'h_shelf', cs%diag%axesT1, cs%Time, &
1565  'ice shelf thickness', 'm', conversion=us%Z_to_m)
1566  cs%id_mass_flux = register_diag_field('ocean_model', 'mass_flux', cs%diag%axesT1,&
1567  cs%Time,'Total mass flux of freshwater across the ice-ocean interface.', 'kg/s')
1568  cs%id_melt = register_diag_field('ocean_model', 'melt', cs%diag%axesT1, cs%Time, &
1569  'Ice Shelf Melt Rate', 'm yr-1')
1570  cs%id_thermal_driving = register_diag_field('ocean_model', 'thermal_driving', cs%diag%axesT1, cs%Time, &
1571  'pot. temp. in the boundary layer minus freezing pot. temp. at the ice-ocean interface.', 'Celsius')
1572  cs%id_haline_driving = register_diag_field('ocean_model', 'haline_driving', cs%diag%axesT1, cs%Time, &
1573  'salinity in the boundary layer minus salinity at the ice-ocean interface.', 'psu')
1574  cs%id_Sbdry = register_diag_field('ocean_model', 'sbdry', cs%diag%axesT1, cs%Time, &
1575  'salinity at the ice-ocean interface.', 'psu')
1576  cs%id_u_ml = register_diag_field('ocean_model', 'u_ml', cs%diag%axesCu1, cs%Time, &
1577  'Eastward vel. in the boundary layer (used to compute ustar)', 'm s-1')
1578  cs%id_v_ml = register_diag_field('ocean_model', 'v_ml', cs%diag%axesCv1, cs%Time, &
1579  'Northward vel. in the boundary layer (used to compute ustar)', 'm s-1')
1580  cs%id_exch_vel_s = register_diag_field('ocean_model', 'exch_vel_s', cs%diag%axesT1, cs%Time, &
1581  'Sub-shelf salinity exchange velocity', 'm s-1')
1582  cs%id_exch_vel_t = register_diag_field('ocean_model', 'exch_vel_t', cs%diag%axesT1, cs%Time, &
1583  'Sub-shelf thermal exchange velocity', 'm s-1')
1584  cs%id_tfreeze = register_diag_field('ocean_model', 'tfreeze', cs%diag%axesT1, cs%Time, &
1585  'In Situ Freezing point at ice shelf interface', 'degC')
1586  cs%id_tfl_shelf = register_diag_field('ocean_model', 'tflux_shelf', cs%diag%axesT1, cs%Time, &
1587  'Heat conduction into ice shelf', 'W m-2')
1588  cs%id_ustar_shelf = register_diag_field('ocean_model', 'ustar_shelf', cs%diag%axesT1, cs%Time, &
1589  'Fric vel under shelf', 'm/s', conversion=us%Z_to_m*us%s_to_T)
1590  if (cs%active_shelf_dynamics) then
1591  cs%id_h_mask = register_diag_field('ocean_model', 'h_mask', cs%diag%axesT1, cs%Time, &
1592  'ice shelf thickness mask', 'none')
1593  endif
1594 
1595  id_clock_shelf = cpu_clock_id('Ice shelf', grain=clock_component)
1596  id_clock_pass = cpu_clock_id(' Ice shelf halo updates', grain=clock_routine)
1597 
1598 end subroutine initialize_ice_shelf
1599 
1600 !> Initializes shelf mass based on three options (file, zero and user)
1601 subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim)
1602 
1603  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1604  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1605  type(ice_shelf_cs), pointer :: CS !< A pointer to the ice shelf control structure
1606  type(ice_shelf_state), intent(inout) :: ISS !< The ice shelf state type that is being updated
1607  logical, optional, intent(in) :: new_sim !< If present and false, this run is being restarted
1608 
1609  integer :: i, j, is, ie, js, je
1610  logical :: read_shelf_area, new_sim_2
1611  character(len=240) :: config, inputdir, shelf_file, filename
1612  character(len=120) :: shelf_mass_var ! The name of shelf mass in the file.
1613  character(len=120) :: shelf_area_var ! The name of shelf area in the file.
1614  character(len=40) :: mdl = "MOM_ice_shelf"
1615  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1616 
1617  new_sim_2 = .true. ; if (present(new_sim)) new_sim_2 = new_sim
1618 
1619  call get_param(param_file, mdl, "ICE_SHELF_CONFIG", config, &
1620  "A string that specifies how the ice shelf is "//&
1621  "initialized. Valid options include:\n"//&
1622  " \tfile\t Read from a file.\n"//&
1623  " \tzero\t Set shelf mass to 0 everywhere.\n"//&
1624  " \tUSER\t Call USER_initialize_shelf_mass.\n", &
1625  fail_if_missing=.true.)
1626 
1627  select case ( trim(config) )
1628  case ("file")
1629 
1630  call time_interp_external_init()
1631 
1632  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1633  inputdir = slasher(inputdir)
1634 
1635  call get_param(param_file, mdl, "SHELF_FILE", shelf_file, &
1636  "If DYNAMIC_SHELF_MASS = True, OVERRIDE_SHELF_MOVEMENT = True "//&
1637  "and ICE_SHELF_MASS_FROM_FILE = True, this is the file from "//&
1638  "which to read the shelf mass and area.", &
1639  default="shelf_mass.nc")
1640  call get_param(param_file, mdl, "SHELF_MASS_VAR", shelf_mass_var, &
1641  "The variable in SHELF_FILE with the shelf mass.", &
1642  default="shelf_mass")
1643  call get_param(param_file, mdl, "READ_SHELF_AREA", read_shelf_area, &
1644  "If true, also read the area covered by ice-shelf from SHELF_FILE.", &
1645  default=.false.)
1646 
1647  filename = trim(slasher(inputdir))//trim(shelf_file)
1648  call log_param(param_file, mdl, "INPUTDIR/SHELF_FILE", filename)
1649 
1650  cs%id_read_mass = init_external_field(filename, shelf_mass_var, &
1651  domain=g%Domain%mpp_domain, verbose=cs%debug)
1652 
1653  if (read_shelf_area) then
1654  call get_param(param_file, mdl, "SHELF_AREA_VAR", shelf_area_var, &
1655  "The variable in SHELF_FILE with the shelf area.", &
1656  default="shelf_area")
1657 
1658  cs%id_read_area = init_external_field(filename,shelf_area_var, &
1659  domain=g%Domain%mpp_domain)
1660  endif
1661 
1662  if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
1663  " initialize_shelf_mass: Unable to open "//trim(filename))
1664 
1665  case ("zero")
1666  do j=js,je ; do i=is,ie
1667  iss%mass_shelf(i,j) = 0.0
1668  iss%area_shelf_h(i,j) = 0.0
1669  enddo ; enddo
1670 
1671  case ("USER")
1672  call user_initialize_shelf_mass(iss%mass_shelf, iss%area_shelf_h, &
1673  iss%h_shelf, iss%hmask, g, cs%US, cs%user_CS, param_file, new_sim_2)
1674 
1675  case default ; call mom_error(fatal,"initialize_ice_shelf: "// &
1676  "Unrecognized ice shelf setup "//trim(config))
1677  end select
1678 
1679 end subroutine initialize_shelf_mass
1680 
1681 !> Updates the ice shelf mass using data from a file.
1682 subroutine update_shelf_mass(G, CS, ISS, Time)
1683  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1684  type(ice_shelf_cs), intent(in) :: CS !< A pointer to the ice shelf control structure
1685  type(ice_shelf_state), intent(inout) :: ISS !< The ice shelf state type that is being updated
1686  type(time_type), intent(in) :: Time !< The current model time
1687 
1688  ! local variables
1689  integer :: i, j, is, ie, js, je
1690  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1691 
1692  call time_interp_external(cs%id_read_mass, time, iss%mass_shelf)
1693 
1694  do j=js,je ; do i=is,ie
1695  iss%area_shelf_h(i,j) = 0.0
1696  iss%hmask(i,j) = 0.
1697  if (iss%mass_shelf(i,j) > 0.0) then
1698  iss%area_shelf_h(i,j) = g%US%L_to_m**2*g%areaT(i,j)
1699  iss%h_shelf(i,j) = iss%mass_shelf(i,j) / cs%rho_ice
1700  iss%hmask(i,j) = 1.
1701  endif
1702  enddo ; enddo
1703 
1704  !call USER_update_shelf_mass(ISS%mass_shelf, ISS%area_shelf_h, ISS%h_shelf, &
1705  ! ISS%hmask, CS%grid, CS%user_CS, Time, .true.)
1706 
1707  if (cs%min_thickness_simple_calve > 0.0) then
1708  call ice_shelf_min_thickness_calve(g, iss%h_shelf, iss%area_shelf_h, iss%hmask, &
1709  cs%min_thickness_simple_calve)
1710  endif
1711 
1712  call pass_var(iss%area_shelf_h, g%domain)
1713  call pass_var(iss%h_shelf, g%domain)
1714  call pass_var(iss%hmask, g%domain)
1715  call pass_var(iss%mass_shelf, g%domain)
1716 
1717 end subroutine update_shelf_mass
1718 
1719 !> Save the ice shelf restart file
1720 subroutine ice_shelf_save_restart(CS, Time, directory, time_stamped, filename_suffix)
1721  type(ice_shelf_cs), pointer :: cs !< ice shelf control structure
1722  type(time_type), intent(in) :: time !< model time at this call
1723  character(len=*), optional, intent(in) :: directory !< An optional directory into which to write
1724  !! these restart files.
1725  logical, optional, intent(in) :: time_stamped !< f true, the restart file names include
1726  !! a unique time stamp. The default is false.
1727  character(len=*), optional, intent(in) :: filename_suffix !< An optional suffix (e.g., a
1728  !! time-stamp) to append to the restart file names.
1729  ! local variables
1730  type(ocean_grid_type), pointer :: g => null()
1731  character(len=200) :: restart_dir
1732 
1733  g => cs%grid
1734 
1735  if (present(directory)) then ; restart_dir = directory
1736  else ; restart_dir = cs%restart_output_dir ; endif
1737 
1738  call save_restart(restart_dir, time, cs%grid, cs%restart_CSp, time_stamped)
1739 
1740 end subroutine ice_shelf_save_restart
1741 
1742 !> Deallocates all memory associated with this module
1743 subroutine ice_shelf_end(CS)
1744  type(ice_shelf_cs), pointer :: cs !< A pointer to the ice shelf control structure
1745 
1746  if (.not.associated(cs)) return
1747 
1748  call ice_shelf_state_end(cs%ISS)
1749 
1750  if (cs%active_shelf_dynamics) call ice_shelf_dyn_end(cs%dCS)
1751 
1752  deallocate(cs)
1753 
1754 end subroutine ice_shelf_end
1755 
1756 !> This routine is for stepping a stand-alone ice shelf model without an ocean.
1757 subroutine solo_time_step(CS, time_step, nsteps, Time, min_time_step_in)
1758  type(ice_shelf_cs), pointer :: cs !< A pointer to the ice shelf control structure
1759  real, intent(in) :: time_step !< The time interval for this update [s].
1760  integer, intent(inout) :: nsteps !< The running number of ice shelf steps.
1761  type(time_type), intent(inout) :: time !< The current model time
1762  real, optional, intent(in) :: min_time_step_in !< The minimum permitted time step [s].
1763 
1764  type(ocean_grid_type), pointer :: g => null()
1765  type(unit_scale_type), pointer :: us => null() ! Pointer to a structure containing
1766  ! various unit conversion factors
1767  type(ice_shelf_state), pointer :: iss => null() !< A structure with elements that describe
1768  !! the ice-shelf state
1769  integer :: is, iec, js, jec, i, j
1770  real :: time_step_remain
1771  real :: time_step_int, min_time_step
1772  character(len=240) :: mesg
1773  logical :: update_ice_vel ! If true, it is time to update the ice shelf velocities.
1774  logical :: coupled_gl ! If true the grouding line position is determined based on
1775  ! coupled ice-ocean dynamics.
1776 
1777  g => cs%grid
1778  us => cs%US
1779  iss => cs%ISS
1780  is = g%isc ; iec = g%iec ; js = g%jsc ; jec = g%jec
1781 
1782  time_step_remain = time_step
1783  if (present (min_time_step_in)) then
1784  min_time_step = min_time_step_in
1785  else
1786  min_time_step = 1000.0 ! This is in seconds - at 1 km resolution it would imply ice is moving at ~1 meter per second
1787  endif
1788 
1789  write (mesg,*) "TIME in ice shelf call, yrs: ", time_type_to_real(time)/(365. * 86400.)
1790  call mom_mesg("solo_time_step: "//mesg)
1791 
1792  do while (time_step_remain > 0.0)
1793  nsteps = nsteps+1
1794 
1795  ! If time_step is not too long, this is unnecessary.
1796  time_step_int = min(ice_time_step_cfl(cs%dCS, iss, g), time_step)
1797 
1798  write (mesg,*) "Ice model timestep = ", time_step_int, " seconds"
1799  if (time_step_int < min_time_step) then
1800  call mom_error(fatal, "MOM_ice_shelf:solo_time_step: abnormally small timestep "//mesg)
1801  else
1802  call mom_mesg("solo_time_step: "//mesg)
1803  endif
1804 
1805  if (time_step_int >= time_step_remain) then
1806  time_step_int = time_step_remain
1807  time_step_remain = 0.0
1808  else
1809  time_step_remain = time_step_remain - time_step_int
1810  endif
1811 
1812  ! If the last mini-timestep is a day or less, we cannot expect velocities to change by much.
1813  ! Do not update the velocities if the last step is very short.
1814  update_ice_vel = ((time_step_int > min_time_step) .or. (time_step_int >= time_step))
1815  coupled_gl = .false.
1816 
1817  call update_ice_shelf(cs%dCS, iss, g, us, time_step_int, time, must_update_vel=update_ice_vel)
1818 
1819  call enable_averaging(time_step,time,cs%diag)
1820  if (cs%id_area_shelf_h > 0) call post_data(cs%id_area_shelf_h, iss%area_shelf_h, cs%diag)
1821  if (cs%id_h_shelf > 0) call post_data(cs%id_h_shelf, iss%h_shelf, cs%diag)
1822  if (cs%id_h_mask > 0) call post_data(cs%id_h_mask, iss%hmask, cs%diag)
1823  call disable_averaging(cs%diag)
1824 
1825  enddo
1826 
1827 end subroutine solo_time_step
1828 
1829 !> \namespace mom_ice_shelf
1830 !!
1831 !! \section section_ICE_SHELF
1832 !!
1833 !! This module implements the thermodynamic aspects of ocean/ice-shelf
1834 !! inter-actions using the MOM framework and coding style.
1835 !!
1836 !! Derived from code by Chris Little, early 2010.
1837 !!
1838 !! The ice-sheet dynamics subroutines do the following:
1839 !! initialize_shelf_mass - Initializes the ice shelf mass distribution.
1840 !! - Initializes h_shelf, h_mask, area_shelf_h
1841 !! - CURRENTLY: initializes mass_shelf as well, but this is unnecessary, as mass_shelf is initialized based on
1842 !! h_shelf and density_ice immediately afterwards. Possibly subroutine should be renamed
1843 !! update_shelf_mass - updates ice shelf mass via netCDF file
1844 !! USER_update_shelf_mass (TODO).
1845 !! solo_time_step - called only in ice-only mode.
1846 !! shelf_calc_flux - after melt rate & fluxes are calculated, ice dynamics are done. currently mass_shelf is
1847 !! updated immediately after ice_shelf_advect in fully dynamic mode.
1848 !!
1849 !! NOTES: be aware that hmask(:,:) has a number of functions; it is used for front advancement,
1850 !! for subroutines in the velocity solve, and for thickness boundary conditions (this last one may be removed).
1851 !! in other words, interfering with its updates will have implications you might not expect.
1852 !!
1853 !! Overall issues: Many variables need better documentation and units and the
1854 !! subgrid on which they are discretized.
1855 !!
1856 !! \subsection section_ICE_SHELF_equations ICE_SHELF equations
1857 !!
1858 !! The three fundamental equations are:
1859 !! Heat flux
1860 !! \f[ \qquad \rho_w C_{pw} \gamma_T (T_w - T_b) = \rho_i \dot{m} L_f \f]
1861 !! Salt flux
1862 !! \f[ \qquad \rho_w \gamma_s (S_w - S_b) = \rho_i \dot{m} S_b \f]
1863 !! Freezing temperature
1864 !! \f[ \qquad T_b = a S_b + b + c P \f]
1865 !!
1866 !! where ....
1867 !!
1868 !! \subsection section_ICE_SHELF_references References
1869 !!
1870 !! Asay-Davis, Xylar S., Stephen L. Cornford, Benjamin K. Galton-Fenzi, Rupert M. Gladstone, G. Hilmar Gudmundsson,
1871 !! David M. Holland, Paul R. Holland, and Daniel F. Martin. Experimental design for three interrelated marine ice sheet
1872 !! and ocean model intercomparison projects: MISMIP v. 3 (MISMIP+), ISOMIP v. 2 (ISOMIP+) and MISOMIP v. 1 (MISOMIP1).
1873 !! Geoscientific Model Development 9, no. 7 (2016): 2471.
1874 !!
1875 !! Goldberg, D. N., et al. Investigation of land ice-ocean interaction with a fully coupled ice-ocean model: 1.
1876 !! Model description and behavior. Journal of Geophysical Research: Earth Surface 117.F2 (2012).
1877 !!
1878 !! Goldberg, D. N., et al. Investigation of land ice-ocean interaction with a fully coupled ice-ocean model: 2.
1879 !! Sensitivity to external forcings. Journal of Geophysical Research: Earth Surface 117.F2 (2012).
1880 !!
1881 !! Holland, David M., and Adrian Jenkins. Modeling thermodynamic ice-ocean interactions at the base of an ice shelf.
1882 !! Journal of Physical Oceanography 29.8 (1999): 1787-1800.
1883 
1884 end module mom_ice_shelf
mom_time_manager
Wraps the FMS time manager functions.
Definition: MOM_time_manager.F90:2
mom_restart::restore_state
subroutine, public restore_state(filename, directory, day, G, CS)
restore_state reads the model state from previously generated files. All restart variables are read f...
Definition: MOM_restart.F90:984
mom_forcing_type::allocate_forcing_type
subroutine, public allocate_forcing_type(G, fluxes, water, heat, ustar, press, shelf, iceberg, salt)
Conditionally allocate fields within the forcing type.
Definition: MOM_forcing_type.F90:2811
mom_ice_shelf::ice_shelf_end
subroutine, public ice_shelf_end(CS)
Deallocates all memory associated with this module.
Definition: MOM_ice_shelf.F90:1744
mom_ice_shelf_dynamics
Implements a crude placeholder for a later implementation of full ice shelf dynamics.
Definition: MOM_ice_shelf_dynamics.F90:3
mom_forcing_type::mech_forcing
Structure that contains pointers to the mechanical forcing at the surface used to drive the liquid oc...
Definition: MOM_forcing_type.F90:187
mom_grid_initialize::set_grid_metrics
subroutine, public set_grid_metrics(G, param_file, US)
set_grid_metrics is used to set the primary values in the model's horizontal grid....
Definition: MOM_grid_initialize.F90:63
mom_variables::surface
Pointers to various fields which may be used describe the surface state of MOM, and which will be ret...
Definition: MOM_variables.F90:38
mom_forcing_type::mom_mech_forcing_chksum
subroutine, public mom_mech_forcing_chksum(mesg, forces, G, US, haloshift)
Write out chksums for the driving mechanical forces.
Definition: MOM_forcing_type.F90:1104
mom_ice_shelf_dynamics::update_ice_shelf
subroutine, public update_ice_shelf(CS, ISS, G, US, time_step, Time, ocean_mass, coupled_grounding, must_update_vel)
This subroutine updates the ice shelf velocities, mass, stresses and properties due to the ice shelf ...
Definition: MOM_ice_shelf_dynamics.F90:609
mom_ice_shelf::id_clock_pass
integer id_clock_pass
CPU Clock for group pass calls.
Definition: MOM_ice_shelf.F90:185
mom_ice_shelf
Implements the thermodynamic aspects of ocean / ice-shelf interactions, along with a crude placeholde...
Definition: MOM_ice_shelf.F90:4
mom_dyn_horgrid::destroy_dyn_horgrid
subroutine, public destroy_dyn_horgrid(G)
Release memory used by the dyn_horgrid_type and related structures.
Definition: MOM_dyn_horgrid.F90:377
mom_ice_shelf_state::ice_shelf_state_init
subroutine, public ice_shelf_state_init(ISS, G)
Deallocates all memory associated with this module.
Definition: MOM_ice_shelf_state.F90:58
mom_ice_shelf::change_thickness_using_melt
subroutine change_thickness_using_melt(ISS, G, time_step, fluxes, rho_ice, debug)
Changes the thickness (mass) of the ice shelf based on sub-ice-shelf melting.
Definition: MOM_ice_shelf.F90:707
mom_checksums::vchksum
Checksums an array (2d or 3d) staggered at C-grid v points.
Definition: MOM_checksums.F90:38
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_ice_shelf_dynamics::register_ice_shelf_dyn_restarts
subroutine, public register_ice_shelf_dyn_restarts(G, param_file, CS, restart_CS)
This subroutine is used to register any fields related to the ice shelf dynamics that should be writt...
Definition: MOM_ice_shelf_dynamics.F90:203
mom_checksums::qchksum
This is an older interface that has been renamed Bchksum.
Definition: MOM_checksums.F90:58
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_ice_shelf_initialize
Initialize ice shelf variables.
Definition: MOM_ice_shelf_initialize.F90:2
mom_get_input::directories
Container for paths and parameter file names.
Definition: MOM_get_input.F90:20
mom_forcing_type::allocate_mech_forcing
subroutine, public allocate_mech_forcing(G, forces, stress, ustar, shelf, press, iceberg)
Conditionally allocate fields within the mechanical forcing type.
Definition: MOM_forcing_type.F90:2879
mom_dyn_horgrid
Contains a shareable dynamic type for describing horizontal grids and metric data and utilty routines...
Definition: MOM_dyn_horgrid.F90:3
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_ice_shelf::add_shelf_pressure
subroutine add_shelf_pressure(G, US, CS, fluxes)
This subroutine adds the ice shelf pressure to the fluxes type.
Definition: MOM_ice_shelf.F90:841
mom_ice_shelf_dynamics::ice_shelf_dyn_end
subroutine, public ice_shelf_dyn_end(CS)
Deallocates all memory associated with the ice shelf dynamics module.
Definition: MOM_ice_shelf_dynamics.F90:3454
mom_diag_mediator::diag_mediator_init
subroutine, public diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir)
diag_mediator_init initializes the MOM diag_mediator and opens the available diagnostics file,...
Definition: MOM_diag_mediator.F90:2972
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_checksums::uvchksum
Checksums a pair velocity arrays (2d or 3d) staggered at C-grid locations.
Definition: MOM_checksums.F90:28
mom_domains::pass_var
Do a halo update on an array.
Definition: MOM_domains.F90:49
user_shelf_init::user_update_shelf_mass
subroutine, public user_update_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, G, CS, Time, new_sim)
This subroutine updates the ice shelf mass, as specified by user-provided code.
Definition: user_shelf_init.F90:126
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_checksums::chksum
This is an older interface for 1-, 2-, or 3-D checksums.
Definition: MOM_checksums.F90:63
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
mom_fixed_initialization::mom_initialize_topography
subroutine, public mom_initialize_topography(D, max_depth, G, PF, US)
MOM_initialize_topography makes the appropriate call to set up the bathymetry. At this point the topo...
Definition: MOM_fixed_initialization.F90:174
mom_restart::mom_restart_cs
A restart registry and the control structure for restarts.
Definition: MOM_restart.F90:72
mom_transcribe_grid::copy_mom_grid_to_dyngrid
subroutine, public copy_mom_grid_to_dyngrid(oG, dG, US)
Copies information from an ocean_grid_type into a dynamic (shared) horizontal grid type.
Definition: MOM_transcribe_grid.F90:167
mom_get_input
Reads the only Fortran name list needed to boot-strap the model.
Definition: MOM_get_input.F90:6
mom_ice_shelf::solo_time_step
subroutine, public solo_time_step(CS, time_step, nsteps, Time, min_time_step_in)
This routine is for stepping a stand-alone ice shelf model without an ocean.
Definition: MOM_ice_shelf.F90:1758
mom_ice_shelf_dynamics::ice_time_step_cfl
real function, public ice_time_step_cfl(CS, ISS, G)
This function returns the global maximum timestep that can be taken based on the current ice velociti...
Definition: MOM_ice_shelf_dynamics.F90:576
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
user_shelf_init::user_initialize_shelf_mass
subroutine, public user_initialize_shelf_mass(mass_shelf, area_shelf_h, h_shelf, hmask, G, US, CS, param_file, new_sim)
This subroutine sets up the initial mass and area covered by the ice shelf, based on user-provided co...
Definition: user_shelf_init.F90:44
mom_eos::calculate_tfreeze
Calculates the freezing point of sea water from T, S and P.
Definition: MOM_EOS.F90:81
mom_coms
Interfaces to non-domain-oriented communication subroutines, including the MOM6 reproducing sums faci...
Definition: MOM_coms.F90:3
mom_ice_shelf_dynamics::ice_shelf_dyn_cs
The control structure for the ice shelf dynamics.
Definition: MOM_ice_shelf_dynamics.F90:41
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_domains::pass_vector
Do a halo update on a pair of arrays representing the two components of a vector.
Definition: MOM_domains.F90:54
mom_ice_shelf::ice_shelf_save_restart
subroutine, public ice_shelf_save_restart(CS, Time, directory, time_stamped, filename_suffix)
Save the ice shelf restart file.
Definition: MOM_ice_shelf.F90:1721
mom_transcribe_grid
Module with routines for copying information from a shared dynamic horizontal grid to an ocean-specif...
Definition: MOM_transcribe_grid.F90:3
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:91
mom_domains::clone_mom_domain
Copy one MOM_domain_type into another.
Definition: MOM_domains.F90:94
mom_forcing_type
This module implements boundary forcing for MOM6.
Definition: MOM_forcing_type.F90:2
mom_domains::to_all
integer, parameter, public to_all
A flag for passing in all directions.
Definition: MOM_domains.F90:132
mom_ice_shelf_state
Implements the thermodynamic aspects of ocean / ice-shelf interactions, along with a crude placeholde...
Definition: MOM_ice_shelf_state.F90:4
mom_grid::mom_grid_init
subroutine, public mom_grid_init(G, param_file, US, HI, global_indexing, bathymetry_at_vel)
MOM_grid_init initializes the ocean grid array sizes and grid memory.
Definition: MOM_grid.F90:184
mom_ice_shelf::ice_shelf_cs
Control structure that contains ice shelf parameters and diagnostics handles.
Definition: MOM_ice_shelf.F90:71
mom_ice_shelf::add_shelf_flux
subroutine, public add_shelf_flux(G, US, CS, state, fluxes)
Updates surface fluxes that are influenced by sub-ice-shelf melting.
Definition: MOM_ice_shelf.F90:870
mom_checksums
Routines to calculate checksums of various array and vector types.
Definition: MOM_checksums.F90:2
mom_restart
The MOM6 facility for reading and writing restart files, and querying what has been read.
Definition: MOM_restart.F90:2
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_dyn_horgrid::create_dyn_horgrid
subroutine, public create_dyn_horgrid(G, HI, bathymetry_at_vel)
Allocate memory used by the dyn_horgrid_type and related structures.
Definition: MOM_dyn_horgrid.F90:176
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_ice_shelf_initialize::initialize_ice_thickness
subroutine, public initialize_ice_thickness(h_shelf, area_shelf_h, hmask, G, US, PF)
Initialize ice shelf thickness.
Definition: MOM_ice_shelf_initialize.F90:29
mom_forcing_type::mom_forcing_chksum
subroutine, public mom_forcing_chksum(mesg, fluxes, G, US, haloshift)
Write out chksums for thermodynamic fluxes.
Definition: MOM_forcing_type.F90:1012
mom_io::mom_read_data
Read a data field from a file.
Definition: MOM_io.F90:74
mom_time_manager::real_to_time
type(time_type) function, public real_to_time(x, err_msg)
This is an alternate implementation of the FMS function real_to_time_type that is accurate over a lar...
Definition: MOM_time_manager.F90:47
mom_ice_shelf_state::ice_shelf_state
Structure that describes the ice shelf state.
Definition: MOM_ice_shelf_state.F90:24
mom_fixed_initialization
Initializes fixed aspects of the model, such as horizontal grid metrics, topography and Coriolis.
Definition: MOM_fixed_initialization.F90:3
mom_ice_shelf::update_shelf_mass
subroutine update_shelf_mass(G, CS, ISS, Time)
Updates the ice shelf mass using data from a file.
Definition: MOM_ice_shelf.F90:1683
mom_cpu_clock
Wraps the MPP cpu clock functions.
Definition: MOM_cpu_clock.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
user_initialization::user_initialize_topography
subroutine, public user_initialize_topography(D, G, param_file, max_depth, US)
Initialize topography.
Definition: user_initialization.F90:65
mom_get_input::get_mom_input
subroutine, public get_mom_input(param_file, dirs, check_params, default_input_filename, ensemble_num)
Get the names of the I/O directories and initialization file. Also calls the subroutine that opens ru...
Definition: MOM_get_input.F90:34
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:70
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_diag_mediator::disable_averaging
subroutine, public disable_averaging(diag_cs)
Call this subroutine to avoid averaging any offered fields.
Definition: MOM_diag_mediator.F90:1840
mom_restart::register_restart_field
Register fields for restarts.
Definition: MOM_restart.F90:107
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_restart::restart_init
subroutine, public restart_init(param_file, CS, restart_root)
Initialize this module and set up a restart control structure.
Definition: MOM_restart.F90:1421
mom_domains::mom_domains_init
subroutine, public mom_domains_init(MOM_dom, param_file, symmetric, static_memory, NIHALO, NJHALO, NIGLOBAL, NJGLOBAL, NIPROC, NJPROC, min_halo, domain_name, include_name, param_suffix)
MOM_domains_init initalizes a MOM_domain_type variable, based on the information read in from a param...
Definition: MOM_domains.F90:1155
mom_grid_initialize
Initializes horizontal grid.
Definition: MOM_grid_initialize.F90:2
mom_forcing_type::forcing
Structure that contains pointers to the boundary forcing used to drive the liquid ocean simulated by ...
Definition: MOM_forcing_type.F90:50
mom_ice_shelf::id_clock_shelf
integer id_clock_shelf
CPU Clock for the ice shelf code.
Definition: MOM_ice_shelf.F90:184
mom_unit_scaling::fix_restart_unit_scaling
subroutine, public fix_restart_unit_scaling(US)
Set the unit scaling factors for output to restart files to the unit scaling factors for this run.
Definition: MOM_unit_scaling.F90:123
mom_io::file_exists
Indicate whether a file exists, perhaps with domain decomposition.
Definition: MOM_io.F90:68
mom_transcribe_grid::copy_dyngrid_to_mom_grid
subroutine, public copy_dyngrid_to_mom_grid(dG, oG, US)
Copies information from a dynamic (shared) horizontal grid type into an ocean_grid_type.
Definition: MOM_transcribe_grid.F90:23
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_restart::save_restart
subroutine, public save_restart(directory, time, G, CS, time_stamped, filename, GV)
save_restart saves all registered variables to restart files.
Definition: MOM_restart.F90:781
user_shelf_init
This module specifies the initial values and evolving properties of the MOM6 ice shelf,...
Definition: user_shelf_init.F90:3
mom_checksums::uchksum
Checksums an array (2d or 3d) staggered at C-grid u points.
Definition: MOM_checksums.F90:33
mom_forcing_type::copy_common_forcing_fields
subroutine, public copy_common_forcing_fields(forces, fluxes, G, skip_pres)
This subroutine copies the computational domains of common forcing fields from a mech_forcing type to...
Definition: MOM_forcing_type.F90:2046
mom_ice_shelf::initialize_ice_shelf
subroutine, public initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fluxes, Time_in, solo_ice_sheet_in)
Initializes shelf model data, parameters and diagnostics.
Definition: MOM_ice_shelf.F90:1082
mom_dyn_horgrid::rescale_dyn_horgrid_bathymetry
subroutine, public rescale_dyn_horgrid_bathymetry(G, m_in_new_units)
rescale_dyn_horgrid_bathymetry permits a change in the internal units for the bathymetry on the grid,...
Definition: MOM_dyn_horgrid.F90:285
mom_ice_shelf::add_shelf_forces
subroutine, public add_shelf_forces(G, US, CS, forces, do_shelf_area)
This subroutine adds the mechanical forcing fields and perhaps shelf areas, based on the ice state in...
Definition: MOM_ice_shelf.F90:760
mom_coms::reproducing_sum
Find an accurate and order-invariant sum of distributed 2d or 3d fields.
Definition: MOM_coms.F90:54
mom_checksums::hchksum
Checksums an array (2d or 3d) staggered at tracer points.
Definition: MOM_checksums.F90:48
mom_eos::eos_init
subroutine, public eos_init(param_file, EOS)
Initializes EOS_type by allocating and reading parameters.
Definition: MOM_EOS.F90:806
mom_diag_mediator::set_diag_mediator_grid
subroutine, public set_diag_mediator_grid(G, diag_cs)
Set up the array extents for doing diagnostics.
Definition: MOM_diag_mediator.F90:1187
mom_file_parser::log_param
An overloaded interface to log the values of various types of parameters.
Definition: MOM_file_parser.F90:96
mom_diag_mediator::enable_averaging
subroutine, public enable_averaging(time_int_in, time_end_in, diag_cs)
This subroutine enables the accumulation of time averages over the specified time interval.
Definition: MOM_diag_mediator.F90:1805
mom_ice_shelf_state::ice_shelf_state_end
subroutine, public ice_shelf_state_end(ISS)
Deallocates all memory associated with this module.
Definition: MOM_ice_shelf_state.F90:87
mom_restart::query_initialized
Indicate whether a field has been read from a restart file.
Definition: MOM_restart.F90:116
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_ice_shelf::shelf_calc_flux
subroutine, public shelf_calc_flux(state, fluxes, Time, time_step, CS, forces)
Calculates fluxes between the ocean and ice-shelf using the three-equations formulation (optional to ...
Definition: MOM_ice_shelf.F90:193
mom_ice_shelf::initialize_shelf_mass
subroutine initialize_shelf_mass(G, param_file, CS, ISS, new_sim)
Initializes shelf mass based on three options (file, zero and user)
Definition: MOM_ice_shelf.F90:1602
mom_dyn_horgrid::dyn_horgrid_type
Describes the horizontal ocean grid with only dynamic memory arrays.
Definition: MOM_dyn_horgrid.F90:23
user_shelf_init::user_ice_shelf_cs
The control structure for the user_ice_shelf module.
Definition: user_shelf_init.F90:29
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_ice_shelf_dynamics::ice_shelf_min_thickness_calve
subroutine, public ice_shelf_min_thickness_calve(G, h_shelf, area_shelf_h, hmask, thickness_calve)
Apply a very simple calving law using a minimum thickness rule.
Definition: MOM_ice_shelf_dynamics.F90:2039
mom_unit_scaling::unit_scaling_init
subroutine, public unit_scaling_init(param_file, US)
Allocates and initializes the ocean model unit scaling type.
Definition: MOM_unit_scaling.F90:44
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_ice_shelf_dynamics::initialize_ice_shelf_dyn
subroutine, public initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_sim, solo_ice_sheet_in)
Initializes shelf model data, parameters and diagnostics.
Definition: MOM_ice_shelf_dynamics.F90:263
user_initialization
A template of a user to code up customized initial conditions.
Definition: user_initialization.F90:2
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60
mom_file_parser::read_param
An overloaded interface to read various types of parameters.
Definition: MOM_file_parser.F90:90