MOM.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> The central module of the MOM6 ocean model
6module mom
7
8! Infrastructure modules
9use mom_array_transform, only : rotate_array, rotate_vector
10use mom_debugging, only : mom_debugging_init, hchksum, uvchksum, totaltands
11use mom_debugging, only : check_redundant, query_debugging_checks
12use mom_checksum_packages, only : mom_thermo_chksum, mom_state_chksum
14use mom_coms, only : num_pes
15use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
16use mom_cpu_clock, only : clock_component, clock_subcomponent
17use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine
21use mom_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
28use mom_domains, only : mom_domains_init, mom_domain_type
29use mom_domains, only : sum_across_pes, pass_var, pass_vector
30use mom_domains, only : clone_mom_domain, deallocate_mom_domain
31use mom_domains, only : to_north, to_east, to_south, to_west
32use mom_domains, only : to_all, omit_corners, cgrid_ne, scalar_pair
33use mom_domains, only : create_group_pass, do_group_pass, group_pass_type
34use mom_domains, only : start_group_pass, complete_group_pass, omit_corners
35use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
38use mom_file_parser, only : read_param, get_param, log_version, param_file_type
39use mom_forcing_type, only : forcing, mech_forcing, find_ustar
43use mom_io, only : slasher, file_exists, mom_read_data
45use mom_restart, only : register_restart_field, register_restart_pair, save_restart
46use mom_restart, only : query_initialized, set_initialized, restart_registry_lock
49use mom_time_manager, only : time_type, real_to_time, operator(+)
50use mom_time_manager, only : operator(-), operator(>), operator(*), operator(/)
51use mom_time_manager, only : operator(>=), operator(==), increment_date
52use mom_unit_tests, only : unit_tests
53
54! MOM core modules
62use mom_ale_sponge, only : rotate_ale_sponge, update_ale_sponge_field
92use mom_eos, only : eos_init, calculate_density, calculate_tfreeze, eos_domain
94use mom_forcing_type, only : allocate_forcing_type, allocate_mech_forcing
104use mom_interface_heights, only : find_eta, calc_derived_thermo, thickness_to_dz
113use mom_meke_types, only : meke_type
140use mom_ale_sponge, only : init_ale_sponge_diags, ale_sponge_cs
147use mom_tracer_registry, only : tracer_registry_type, register_tracer, tracer_registry_init
149use mom_tracer_registry, only : post_tracer_transport_diagnostics, mom_tracer_chksum
165
166! Database client used for machine-learning interface
168
169! ODA modules
173
174! Offline modules
184implicit none ; private
185
186#include <MOM_memory.h>
187
188! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
189! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
190! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
191! vary with the Boussinesq approximation, the Boussinesq variant is given first.
192
193!> A structure with diagnostic IDs of the state variables
194type mom_diag_ids
195 !>@{ 3-d state field diagnostic IDs
196 integer :: id_u = -1, id_v = -1, id_h = -1
197 !>@}
198 !> 2-d state field diagnostic ID
199 integer :: id_ssh_inst = -1
200end type mom_diag_ids
201
202!> Control structure for the MOM module, including the variables that describe
203!! the state of the ocean.
204type, public :: mom_control_struct ; private
205 real allocable_, dimension(NIMEM_,NJMEM_,NKMEM_) :: &
206 h, & !< layer thickness [H ~> m or kg m-2]
207 t, & !< potential temperature [C ~> degC]
208 s !< salinity [S ~> ppt]
209 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
210 u, & !< zonal velocity component [L T-1 ~> m s-1]
211 uh, & !< uh = u * h * dy at u grid points [H L2 T-1 ~> m3 s-1 or kg s-1]
212 uhtr !< accumulated zonal thickness fluxes to advect tracers [H L2 ~> m3 or kg]
213 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
214 v, & !< meridional velocity [L T-1 ~> m s-1]
215 vh, & !< vh = v * h * dx at v grid points [H L2 T-1 ~> m3 s-1 or kg s-1]
216 vhtr !< accumulated meridional thickness fluxes to advect tracers [H L2 ~> m3 or kg]
217 real allocable_, dimension(NIMEM_,NJMEM_) :: ssh_rint
218 !< A running time integral of the sea surface height [T Z ~> s m].
219 real allocable_, dimension(NIMEM_,NJMEM_) :: ave_ssh_ibc
220 !< time-averaged (over a forcing time step) sea surface height
221 !! with a correction for the inverse barometer [Z ~> m]
222 real allocable_, dimension(NIMEM_,NJMEM_) :: eta_av_bc
223 !< free surface height or column mass time averaged over the last
224 !! baroclinic dynamics time step [H ~> m or kg m-2]
225 real, dimension(:,:), pointer :: hml => null()
226 !< active mixed layer depth, or 0 if there is no boundary layer scheme [Z ~> m]
227 real :: time_in_cycle !< The running time of the current time-stepping cycle
228 !! in calls that step the dynamics, and also the length of
229 !! the time integral of ssh_rint [T ~> s].
230 real :: time_in_thermo_cycle !< The running time of the current time-stepping
231 !! cycle in calls that step the thermodynamics [T ~> s].
232
233 type(ocean_grid_type) :: g_in !< Input grid metric
234 type(ocean_grid_type), pointer :: g => null() !< Model grid metric
235 logical :: rotate_index = .false. !< True if index map is rotated
236 logical :: homogenize_forcings = .false. !< True if all inputs are homogenized
237 logical :: update_ustar = .false. !< True to update ustar from homogenized tau
238 logical :: vertex_shear = .false. !< True if vertex shear is on
239
240 type(verticalgrid_type), pointer :: &
241 gv => null() !< structure containing vertical grid info
242 type(unit_scale_type), pointer :: &
243 us => null() !< structure containing various unit conversion factors
244 type(thermo_var_ptrs) :: tv !< structure containing pointers to available thermodynamic fields
245 real :: t_dyn_rel_adv !< The time of the dynamics relative to tracer advection and lateral mixing
246 !! [T ~> s], or equivalently the elapsed time since advectively updating the
247 !! tracers. t_dyn_rel_adv is invariably positive and may span multiple coupling timesteps.
248 integer :: n_dyn_steps_in_adv !< The number of dynamics time steps that contributed to uhtr
249 !! and vhtr since the last time tracer advection occured.
250 real :: t_dyn_rel_thermo !< The time of the dynamics relative to diabatic processes and remapping
251 !! [T ~> s]. t_dyn_rel_thermo can be negative or positive depending on whether
252 !! the diabatic processes are applied before or after the dynamics and may span
253 !! multiple coupling timesteps.
254 real :: t_dyn_rel_diag !< The time of the diagnostics relative to diabatic processes and remapping
255 !! [T ~> s]. t_dyn_rel_diag is always positive, since the diagnostics must lag.
256 logical :: preadv_h_stored = .false. !< If true, the thicknesses from before the advective cycle
257 !! have been stored for use in diagnostics.
258
259 type(diag_ctrl) :: diag !< structure to regulate diagnostic output timing
260 type(vertvisc_type) :: visc !< structure containing vertical viscosities,
261 !! bottom drag viscosities, and related fields
262 type(meke_type) :: meke !< Fields related to the Mesoscale Eddy Kinetic Energy
263 logical :: adiabatic !< If true, there are no diapycnal mass fluxes, and no calls
264 !! to routines to calculate or apply diapycnal fluxes.
265 logical :: diabatic_first !< If true, apply diabatic and thermodynamic processes before time
266 !! stepping the dynamics.
267 logical :: use_ale_algorithm !< If true, use the ALE algorithm rather than layered
268 !! isopycnal/stacked shallow water mode. This logical is set by calling the
269 !! function useRegridding() from the MOM_regridding module.
270 logical :: remap_aux_vars !< If true, apply ALE remapping to all of the auxiliary 3-D
271 !! variables that are needed to reproduce across restarts,
272 !! similarly to what is done with the primary state variables.
273 logical :: remap_uv_using_old_alg !< If true, use the old "remapping via a delta z" method for
274 !! velocities. If false, remap between two grids described by thicknesses.
275
276 type(mom_stoch_eos_cs) :: stoch_eos_cs !< structure containing random pattern for stoch EOS
277 logical :: alternate_first_direction !< If true, alternate whether the x- or y-direction
278 !! updates occur first in directionally split parts of the calculation.
279 real :: first_dir_restart = -1.0 !< A real copy of G%first_direction for use in restart files [nondim]
280 logical :: offline_tracer_mode = .false.
281 !< If true, step_offline() is called instead of step_MOM().
282 !! This is intended for running MOM6 in offline tracer mode
283 logical :: meke_in_dynamics !< If .true. (default), MEKE is called in the dynamics routine otherwise
284 !! it is called during the tracer dynamics
285
286 type(time_type), pointer :: time !< pointer to the ocean clock
287 real :: dt !< (baroclinic) dynamics time step [T ~> s]
288 real :: dt_therm !< diabatic time step [T ~> s]
289 real :: dt_tr_adv !< tracer advection time step [T ~> s]
290 logical :: thermo_spans_coupling !< If true, thermodynamic and tracer time
291 !! steps can span multiple coupled time steps.
292 logical :: tradv_spans_coupling !< If true, thermodynamic and tracer time
293 integer :: nstep_tot = 0 !< The total number of dynamic timesteps taken
294 !! so far in this run segment
295 logical :: count_calls = .false. !< If true, count the calls to step_MOM, rather than the
296 !! number of dynamics steps in nstep_tot
297 logical :: debug !< If true, write verbose checksums for debugging purposes.
298 logical :: debug_obcs !< If true, write verbose OBC values for debugging purposes.
299 integer :: ntrunc !< number u,v truncations since last call to write_energy
300
301 integer :: cont_stencil !< The stencil for thickness from the continuity solver.
302 integer :: dyn_h_stencil !< The stencil for thickness for the dynamics based on
303 !! the continuity solver and Coriolis schemes.
304 ! These elements are used to control the dynamics updates.
305 logical :: do_dynamics !< If false, does not call step_MOM_dyn_*. This is an
306 !! undocumented run-time flag that is fragile.
307 logical :: split !< If true, use the split time stepping scheme.
308 logical :: use_alt_split !< If true, use a version of the split explicit time stepping
309 !! scheme that exchanges velocities with step_MOM that have the
310 !! average barotropic phase over a baroclinic timestep rather
311 !! than the instantaneous barotropic phase.
312 logical :: use_rk2 !< If true, use RK2 instead of RK3 in unsplit mode
313 !! (i.e., no split between barotropic and baroclinic).
314 logical :: interface_filter !< If true, apply an interface height filter immediately
315 !! after any calls to thickness_diffuse.
316 logical :: thickness_diffuse !< If true, diffuse interface height w/ a diffusivity KHTH.
317 logical :: thickness_diffuse_first !< If true, diffuse thickness before dynamics.
318 logical :: interface_filter_dt_bug !< If true, uses the wrong time interval in
319 !! calls to interface_filter and thickness_diffuse.
320 logical :: mixedlayer_restrat !< If true, use submesoscale mixed layer restratifying scheme.
321 logical :: usemeke !< If true, call the MEKE parameterization.
322 logical :: use_stochastic_eos !< If true, use the stochastic EOS parameterizations.
323 logical :: usewaves !< If true, update Stokes drift
324 real :: dtbt_reset_period !< The time interval between dynamic recalculation of the
325 !! barotropic time step [T ~> s]. If this is negative dtbt is never
326 !! calculated, and if it is 0, dtbt is calculated every step.
327 type(time_type) :: dtbt_reset_interval !< A time_time representation of dtbt_reset_period.
328 type(time_type) :: dtbt_reset_time !< The next time DTBT should be calculated.
329 real :: dt_obc_seg_period !< The time interval between OBC segment updates for OBGC
330 !! tracers [T ~> s], or a negative value if the segment
331 !! data are time-invarant, or zero to update the OBGC
332 !! segment data with every call to update_OBC_segment_data.
333 type(time_type) :: dt_obc_seg_interval !< A time_time representation of dt_obc_seg_period.
334 type(time_type) :: dt_obc_seg_time !< The next time OBC segment update is applied to OBGC tracers.
335
336 real, dimension(:,:), pointer :: frac_shelf_h => null() !< fraction of total area occupied
337 !! by ice shelf [nondim]
338 real, dimension(:,:), pointer :: mass_shelf => null() !< Mass of ice shelf [R Z ~> kg m-2]
339 type(accel_diag_ptrs) :: adp !< structure containing pointers to accelerations,
340 !! for derived diagnostics (e.g., energy budgets)
341 type(cont_diag_ptrs) :: cdp !< structure containing pointers to continuity equation
342 !! terms, for derived diagnostics (e.g., energy budgets)
343 real, dimension(:,:,:), pointer :: &
344 u_prev => null(), & !< previous value of u stored for diagnostics [L T-1 ~> m s-1]
345 v_prev => null() !< previous value of v stored for diagnostics [L T-1 ~> m s-1]
346
347 logical :: interp_p_surf !< If true, linearly interpolate surface pressure
348 !! over the coupling time step, using specified value
349 !! at the end of the coupling step. False by default.
350 logical :: p_surf_prev_set !< If true, p_surf_prev has been properly set from
351 !! a previous time-step or the ocean restart file.
352 !! This is only valid when interp_p_surf is true.
353 real, dimension(:,:), pointer :: &
354 p_surf_prev => null(), & !< surface pressure [R L2 T-2 ~> Pa] at end previous call to step_MOM
355 p_surf_begin => null(), & !< surface pressure [R L2 T-2 ~> Pa] at start of step_MOM_dyn_...
356 p_surf_end => null() !< surface pressure [R L2 T-2 ~> Pa] at end of step_MOM_dyn_...
357
358 ! Variables needed to reach between start and finish phases of initialization
359 logical :: write_ic !< If true, then the initial conditions will be written to file
360 character(len=120) :: ic_file !< A file into which the initial conditions are
361 !! written in a new run if SAVE_INITIAL_CONDS is true.
362
363 logical :: calc_rho_for_sea_lev !< If true, calculate rho to convert pressure to sea level
364
365 ! These elements are used to control the calculation and error checking of the surface state
366 real :: hmix !< Diagnostic mixed layer thickness over which to
367 !! average surface tracer properties when a bulk
368 !! mixed layer is not used [H ~> m or kg m-2], or a negative value
369 !! if a bulk mixed layer is being used.
370 real :: hfrz !< If HFrz > 0, the nominal depth over which melt potential is computed
371 !! [H ~> m or kg m-2]. The actual depth over which melt potential is
372 !! computed is min(HFrz, OBLD), where OBLD is the boundary layer depth.
373 !! If HFrz <= 0 (default), melt potential will not be computed.
374 real :: hmix_uv !< Depth scale over which to average surface flow to
375 !! feedback to the coupler/driver [H ~> m or kg m-2] when
376 !! bulk mixed layer is not used, or a negative value
377 !! if a bulk mixed layer is being used.
378 logical :: check_bad_sfc_vals !< If true, scan surface state for ridiculous values.
379 real :: bad_val_ssh_max !< Maximum SSH before triggering bad value message [Z ~> m]
380 real :: bad_val_sst_max !< Maximum SST before triggering bad value message [C ~> degC]
381 real :: bad_val_sst_min !< Minimum SST before triggering bad value message [C ~> degC]
382 real :: bad_val_sss_max !< Maximum SSS before triggering bad value message [S ~> ppt]
383 real :: bad_val_col_thick !< Minimum column thickness before triggering bad value message [Z ~> m]
384 integer :: answer_date !< The vintage of the expressions for the surface properties. Values
385 !! below 20190101 recover the answers from the end of 2018, while
386 !! higher values use more appropriate expressions that differ at
387 !! roundoff for non-Boussinesq cases.
388 logical :: use_particles !< Turns on the particles package
389 logical :: use_uh_particles !< particles are advected by uh/h
390 logical :: uh_particles_bug !< If true, uses an inconsistent timestep for particle advection
391 logical :: use_dbclient !< Turns on the database client used for ML inference/analysis
392 character(len=10) :: particle_type !< Particle types include: surface(default), profiling and sail drone.
393
394 type(mom_diag_ids) :: ids !< Handles used for diagnostics.
395 type(transport_diag_ids) :: transport_ids !< Handles used for transport diagnostics.
396 type(surface_diag_ids) :: sfc_ids !< Handles used for surface diagnostics.
397 type(diag_grid_storage) :: diag_pre_sync !< The grid (thicknesses) before remapping
398 type(diag_grid_storage) :: diag_pre_dyn !< The grid (thicknesses) before dynamics
399
400 ! The remainder of this type provides pointers to child module control structures.
401
402 type(mom_dyn_unsplit_cs), pointer :: dyn_unsplit_csp => null()
403 !< Pointer to the control structure used for the unsplit dynamics
404 type(mom_dyn_unsplit_rk2_cs), pointer :: dyn_unsplit_rk2_csp => null()
405 !< Pointer to the control structure used for the unsplit RK2 dynamics
406 type(mom_dyn_split_rk2_cs), pointer :: dyn_split_rk2_csp => null()
407 !< Pointer to the control structure used for the mode-split RK2 dynamics
408 type(mom_dyn_split_rk2b_cs), pointer :: dyn_split_rk2b_csp => null()
409 !< Pointer to the control structure used for an alternate version of the mode-split RK2 dynamics
410 type(harmonic_analysis_cs), pointer :: ha_csp => null()
411 !< Pointer to the control structure for harmonic analysis
412 type(thickness_diffuse_cs) :: thickness_diffuse_csp
413 !< Pointer to the control structure used for the isopycnal height diffusive transport.
414 !! This is also common referred to as Gent-McWilliams diffusion
415 type(interface_filter_cs) :: interface_filter_csp
416 !< Control structure used for the interface height smoothing operator.
417 type(mixedlayer_restrat_cs) :: mixedlayer_restrat_csp
418 !< Pointer to the control structure used for the mixed layer restratification
419 type(set_visc_cs) :: set_visc_csp
420 !< Pointer to the control structure used to set viscosities
421 type(diabatic_cs), pointer :: diabatic_csp => null()
422 !< Pointer to the control structure for the diabatic driver
423 type(meke_cs) :: meke_csp
424 !< Pointer to the control structure for the MEKE updates
425 type(varmix_cs) :: varmix
426 !< Control structure for the variable mixing module
427 type(tracer_registry_type), pointer :: tracer_reg => null()
428 !< Pointer to the MOM tracer registry
429 type(tracer_advect_cs), pointer :: tracer_adv_csp => null()
430 !< Pointer to the MOM tracer advection control structure
431 type(tracer_hor_diff_cs), pointer :: tracer_diff_csp => null()
432 !< Pointer to the MOM along-isopycnal tracer diffusion control structure
433 type(tracer_flow_control_cs), pointer :: tracer_flow_csp => null()
434 !< Pointer to the control structure that orchestrates the calling of tracer packages
435 ! Although update_OBC_CS is not used directly outside of initialization, other modules
436 ! set pointers to this type, so it should be kept for the duration of the run.
437 type(update_obc_cs), pointer :: update_obc_csp => null()
438 !< Pointer to the control structure for updating open boundary condition properties
439 type(ocean_obc_type), pointer :: obc => null()
440 !< Pointer to the MOM open boundary condition type
441 type(sponge_cs), pointer :: sponge_csp => null()
442 !< Pointer to the layered-mode sponge control structure
443 type(ale_sponge_cs), pointer :: ale_sponge_csp => null()
444 !< Pointer to the ALE-mode sponge control structure
445 type(oda_incupd_cs), pointer :: oda_incupd_csp => null()
446 !< Pointer to the oda incremental update control structure
447 type(int_tide_cs), pointer :: int_tide_csp => null()
448 !< Pointer to the internal tides control structure
449 type(ale_cs), pointer :: ale_csp => null()
450 !< Pointer to the Arbitrary Lagrangian Eulerian (ALE) vertical coordinate control structure
451
452 ! Pointers to control structures used for diagnostics
453 type(sum_output_cs), pointer :: sum_output_csp => null()
454 !< Pointer to the globally summed output control structure
455 type(diagnostics_cs) :: diagnostics_csp
456 !< Pointer to the MOM diagnostics control structure
457 type(offline_transport_cs), pointer :: offline_csp => null()
458 !< Pointer to the offline tracer transport control structure
459 type(porous_barrier_cs) :: por_bar_cs
460 !< Control structure for porous barrier
461
462 logical :: ensemble_ocean !< if true, this run is part of a
463 !! larger ensemble for the purpose of data assimilation
464 !! or statistical analysis.
465 type(oda_cs), pointer :: odacs => null() !< a pointer to the control structure for handling
466 !! ensemble model state vectors and data assimilation
467 !! increments and priors
468 type(dbcomms_cs_type) :: dbcomms_cs !< Control structure for database client used for online ML/AI
469 logical :: use_porbar !< If true, use porous barrier to constrain the widths and face areas
470 !! at the edges of the grid cells.
471 type(porous_barrier_type) :: pbv !< porous barrier fractional cell metrics
472 type(particles), pointer :: particles => null() !<Lagrangian particles
473 type(stochastic_cs), pointer :: stoch_cs => null() !< a pointer to the stochastics control structure
474 type(mom_restart_cs), pointer :: restart_cs => null()
475 !< Pointer to MOM's restart control structure
476end type mom_control_struct
477
482public allocate_surface_state, deallocate_surface_state
483public save_mom_restart
484
485!>@{ CPU time clock IDs
486integer :: id_clock_ocean
487integer :: id_clock_dynamics
488integer :: id_clock_thermo
489integer :: id_clock_mom_end
490integer :: id_clock_remap
491integer :: id_clock_tracer
492integer :: id_clock_diabatic
493integer :: id_clock_adiabatic
494integer :: id_clock_continuity ! also in dynamics s/r
495integer :: id_clock_thick_diff
496integer :: id_clock_int_filter
497integer :: id_clock_bbl_visc
498integer :: id_clock_ml_restrat
499integer :: id_clock_diagnostics
500integer :: id_clock_z_diag
501integer :: id_clock_init
502integer :: id_clock_mom_init
503integer :: id_clock_pass ! also in dynamics d/r
504integer :: id_clock_pass_init ! also in dynamics d/r
505integer :: id_clock_ale
506integer :: id_clock_other
507integer :: id_clock_offline_tracer
508integer :: id_clock_save_restart
509integer :: id_clock_unit_tests
510integer :: id_clock_stoch
511integer :: id_clock_vart
512!>@}
513
514contains
515
516!> This subroutine orchestrates the time stepping of MOM. The adiabatic
517!! dynamics are stepped by calls to one of the step_MOM_dyn_...routines.
518!! The action of lateral processes on tracers occur in calls to
519!! advect_tracer and tracer_hordiff. Vertical mixing and possibly remapping
520!! occur inside of diabatic.
521subroutine step_mom(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS, &
522 Waves, do_dynamics, do_thermodynamics, start_cycle, &
523 end_cycle, cycle_length, reset_therm)
524 type(mech_forcing), target, intent(inout) :: forces_in !< A structure with the driving mechanical forces
525 type(forcing), target, intent(inout) :: fluxes_in !< A structure with pointers to themodynamic,
526 !! tracer and mass exchange forcing fields
527 type(surface), target, intent(inout) :: sfc_state !< surface ocean state
528 type(time_type), intent(in) :: time_start !< starting time of a segment, as a time type
529 real, intent(in) :: time_int_in !< time interval covered by this run segment [T ~> s].
530 type(mom_control_struct), intent(inout), target :: cs !< control structure from initialize_MOM
531 type(wave_parameters_cs), &
532 optional, pointer :: waves !< An optional pointer to a wave property CS
533 logical, optional, intent(in) :: do_dynamics !< Present and false, do not do updates due
534 !! to the dynamics.
535 logical, optional, intent(in) :: do_thermodynamics !< Present and false, do not do updates due
536 !! to the thermodynamics or remapping.
537 logical, optional, intent(in) :: start_cycle !< This indicates whether this call is to be
538 !! treated as the first call to step_MOM in a
539 !! time-stepping cycle; missing is like true.
540 logical, optional, intent(in) :: end_cycle !< This indicates whether this call is to be
541 !! treated as the last call to step_MOM in a
542 !! time-stepping cycle; missing is like true.
543 real, optional, intent(in) :: cycle_length !< The amount of time in a coupled time
544 !! stepping cycle [T ~> s].
545 logical, optional, intent(in) :: reset_therm !< This indicates whether the running sums of
546 !! thermodynamic quantities should be reset.
547 !! If missing, this is like start_cycle.
548
549 ! local variables
550 type(ocean_grid_type), pointer :: g => null() ! pointer to a structure containing
551 ! metrics and related information
552 type(ocean_grid_type), pointer :: g_in => null() ! Input grid metric
553 type(verticalgrid_type), pointer :: gv => null() ! Pointer to the vertical grid structure
554 type(unit_scale_type), pointer :: us => null() ! Pointer to a structure containing
555 ! various unit conversion factors
556 integer :: ntstep ! number of time steps between diabatic forcing updates
557 integer :: ntastep ! number of time steps between tracer advection updates
558 integer :: n_max ! number of steps to take in this call
559 integer :: halo_sz, dynamics_stencil
560
561 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, n
562 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
563
564 real :: time_interval ! time interval covered by this run segment [T ~> s].
565 real :: dt ! baroclinic time step [T ~> s]
566 real :: dtdia ! time step for diabatic processes [T ~> s]
567 real :: dt_tr_adv ! time step for tracer advection [T ~> s]
568 real :: dt_therm ! a limited and quantized version of CS%dt_therm [T ~> s]
569 real :: dt_tradv_here ! a further limited value of dt_tr_adv [T ~> s]
570
571 real :: wt_end, wt_beg ! Fractional weights of the future pressure at the end
572 ! and beginning of the current time step [nondim]
573 real :: bbl_time_int ! The amount of time over which the calculated BBL
574 ! properties will apply, for use in diagnostics, or 0
575 ! if it is not to be calculated anew [T ~> s].
576 real :: rel_time = 0.0 ! relative time since start of this call [T ~> s].
577
578 logical :: do_advection ! If true, do tracer advection.
579 logical :: do_diabatic ! If true, do diabatic update.
580 logical :: thermo_does_span_coupling ! If true,thermodynamic (diabatic) forcing spans
581 ! multiple coupling timesteps.
582 logical :: tradv_does_span_coupling ! If true, tracer advection spans
583 ! multiple coupling timesteps.
584 logical :: do_dyn ! If true, dynamics are updated with this call.
585 logical :: do_thermo ! If true, thermodynamics and remapping may be applied with this call.
586 logical :: debug_redundant ! If true, check redundant values on PE boundaries when debugging.
587 logical :: nonblocking_p_surf_update ! A flag to indicate whether surface properties
588 ! can use nonblocking halo updates
589 logical :: cycle_start ! If true, do calculations that are only done at the start of
590 ! a stepping cycle (whatever that may mean).
591 logical :: cycle_end ! If true, do calculations and diagnostics that are only done at
592 ! the end of a stepping cycle (whatever that may mean).
593 logical :: therm_reset ! If true, reset running sums of thermodynamic quantities.
594 real :: cycle_time ! The length of the coupled time-stepping cycle [T ~> s].
595 real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
596 u_star ! The wind friction velocity, calculated using the Boussinesq reference density or
597 ! the time-evolving surface density in non-Boussinesq mode [Z T-1 ~> m s-1]
598 real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
599 ssh ! sea surface height, which may be based on eta_av [Z ~> m]
600 real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%GV)) :: &
601 dz ! Vertical distance across layers [Z ~> m]
602
603 real, dimension(:,:,:), pointer :: &
604 u => null(), & ! u : zonal velocity component [L T-1 ~> m s-1]
605 v => null(), & ! v : meridional velocity component [L T-1 ~> m s-1]
606 h => null() ! h : layer thickness [H ~> m or kg m-2]
607 real, dimension(:,:), pointer :: &
608 p_surf => null() ! A pointer to the ocean surface pressure [R L2 T-2 ~> Pa].
609 real :: i_wt_ssh ! The inverse of the time weights [T-1 ~> s-1]
610
611 type(time_type) :: time_local, end_time_thermo
612 type(time_type) :: time_end_diag ! End time of a diagnostic segment, as a time type
613
614 type(group_pass_type) :: pass_tau_ustar_psurf
615 logical :: showcalltree
616
617 ! External forcing fields on the model index map
618 type(mech_forcing), pointer :: forces ! Mechanical forcing
619 type(forcing), pointer :: fluxes ! Boundary fluxes
620 type(surface), pointer :: sfc_state_diag ! Surface boundary fields
621 integer :: turns ! Number of quarter turns from input to model indexing
622
623 g => cs%G ; g_in => cs%G_in ; gv => cs%GV ; us => cs%US
624 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
625 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
626 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
627 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
628 u => cs%u ; v => cs%v ; h => cs%h
629
630 time_interval = time_int_in
631 do_dyn = .true. ; if (present(do_dynamics)) do_dyn = do_dynamics
632 do_thermo = .true. ; if (present(do_thermodynamics)) do_thermo = do_thermodynamics
633 if (.not.(do_dyn .or. do_thermo)) call mom_error(fatal,"Step_MOM: "//&
634 "Both do_dynamics and do_thermodynamics are false, which makes no sense.")
635 cycle_start = .true. ; if (present(start_cycle)) cycle_start = start_cycle
636 cycle_end = .true. ; if (present(end_cycle)) cycle_end = end_cycle
637 cycle_time = time_interval ; if (present(cycle_length)) cycle_time = cycle_length
638 therm_reset = cycle_start ; if (present(reset_therm)) therm_reset = reset_therm
639
640 call cpu_clock_begin(id_clock_ocean)
641 call cpu_clock_begin(id_clock_other)
642
643 if (cs%debug) then
644 call query_debugging_checks(do_redundant=debug_redundant)
645 call mom_state_chksum("Beginning of step_MOM ", u, v, h, cs%uh, cs%vh, g, gv, us)
646 endif
647
648 showcalltree = calltree_showquery()
649 if (showcalltree) call calltree_enter("step_MOM(), MOM.F90")
650
651 ! Rotate the forces from G_in to G
652 if (cs%rotate_index) then
653 turns = g%HI%turns
654 allocate(forces)
655 call allocate_mech_forcing(forces_in, g, forces)
656 call rotate_mech_forcing(forces_in, turns, forces)
657
658 allocate(fluxes)
659 call allocate_forcing_type(fluxes_in, g, fluxes, turns=turns)
660 call rotate_forcing(fluxes_in, fluxes, turns)
661 else
662 forces => forces_in
663 fluxes => fluxes_in
664 endif
665
666 ! Homogenize the forces
667 if (cs%homogenize_forcings) then
668 ! Homogenize all forcing and fluxes fields.
669 call homogenize_mech_forcing(forces, g, us, gv%Rho0, cs%update_ustar)
670 ! Note the following computes the mean ustar as the mean of ustar rather than
671 ! ustar of the mean of tau.
672 call homogenize_forcing(fluxes, g, gv, us)
673 if (cs%update_ustar) then
674 ! These calls corrects the ustar values
675 call copy_common_forcing_fields(forces, fluxes, g)
676 call set_derived_forcing_fields(forces, fluxes, g, us, gv%Rho0)
677 endif
678 endif
679
680 ! This will be replaced later with the pressures from forces or fluxes if they are available.
681 if (associated(cs%tv%p_surf)) cs%tv%p_surf(:,:) = 0.0
682
683 ! First determine the time step that is consistent with this call and an
684 ! integer fraction of time_interval.
685 if (do_dyn) then
686 n_max = 1
687 if (time_interval > cs%dt) n_max = ceiling(time_interval/cs%dt - 0.001)
688
689 dt = time_interval / real(n_max)
690 thermo_does_span_coupling = (cs%thermo_spans_coupling .and. &
691 (cs%dt_therm > 1.5*cycle_time))
692 tradv_does_span_coupling = (cs%tradv_spans_coupling .and. &
693 (cs%dt_tr_adv > 1.5*cycle_time))
694 if (thermo_does_span_coupling) then
695 ! Set dt_therm to be an integer multiple of the coupling time step.
696 dt_therm = cycle_time * floor(cs%dt_therm / cycle_time + 0.001)
697 ntstep = floor(dt_therm/dt + 0.001)
698 elseif (.not.do_thermo) then
699 dt_therm = cs%dt_therm
700 if (present(cycle_length)) dt_therm = min(cs%dt_therm, cycle_length)
701 ntstep = 1 ! ntstep is initialized to avoid an error in a secondary logical test,
702 ! but the nonzero value of ntstep does not matter when do_thermo is false.
703 else
704 ntstep = max(1, min(n_max, floor(cs%dt_therm/dt + 0.001)))
705 dt_therm = dt*ntstep
706 endif
707 if (tradv_does_span_coupling) then
708 ! Set dt_tr_adv to be an integer multiple of the coupling time step.
709 dt_tr_adv = cycle_time * floor(cs%dt_tr_adv / cycle_time + 0.001)
710 ntastep = floor(dt_tr_adv/dt + 0.001)
711 elseif (.not.do_thermo) then
712 dt_tr_adv = cs%dt_tr_adv
713 if (present(cycle_length)) dt_tr_adv = min(cs%dt_tr_adv, cycle_length)
714 ! ntastep is not used.
715 else
716 ntastep = max(1, min(n_max, floor(cs%dt_tr_adv/dt + 0.001)))
717 dt_tr_adv = dt*ntastep
718 endif
719
720 !---------- Initiate group halo pass of the forcing fields
721 call cpu_clock_begin(id_clock_pass)
722 ! Halo updates for surface pressure need to be completed before calling calc_resoln_function
723 ! among other routines if the surface pressure is used in the equation of state.
724 nonblocking_p_surf_update = g%nonblocking_updates .and. &
725 .not.(associated(cs%tv%p_surf) .and. associated(forces%p_surf) .and. &
726 allocated(cs%tv%SpV_avg) .and. associated(cs%tv%T))
727 if (.not.associated(forces%taux) .or. .not.associated(forces%tauy)) &
728 call mom_error(fatal,'step_MOM:forces%taux,tauy not associated')
729 call create_group_pass(pass_tau_ustar_psurf, forces%taux, forces%tauy, g%Domain)
730 if (associated(forces%ustar)) &
731 call create_group_pass(pass_tau_ustar_psurf, forces%ustar, g%Domain)
732 if (associated(forces%tau_mag)) &
733 call create_group_pass(pass_tau_ustar_psurf, forces%tau_mag, g%Domain)
734 if (associated(forces%p_surf)) &
735 call create_group_pass(pass_tau_ustar_psurf, forces%p_surf, g%Domain)
736 if (nonblocking_p_surf_update) then
737 call start_group_pass(pass_tau_ustar_psurf, g%Domain)
738 else
739 call do_group_pass(pass_tau_ustar_psurf, g%Domain)
740 endif
741 call cpu_clock_end(id_clock_pass)
742
743 if (associated(forces%p_surf)) p_surf => forces%p_surf
744 if (.not.associated(forces%p_surf)) cs%interp_p_surf = .false.
745 if (associated(cs%tv%p_surf) .and. associated(forces%p_surf)) then
746 do j=jsd,jed ; do i=isd,ied ; cs%tv%p_surf(i,j) = forces%p_surf(i,j) ; enddo ; enddo
747
748 if (allocated(cs%tv%SpV_avg) .and. associated(cs%tv%T)) then
749 ! The internal ocean state depends on the surface pressues, so update SpV_avg.
750 dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
751 call calc_derived_thermo(cs%tv, h, g, gv, us, halo=dynamics_stencil, debug=cs%debug)
752 endif
753 endif
754
755 else
756 ! This step only updates the thermodynamics so setting timesteps is simpler.
757 n_max = 1
758 if ((time_interval > cs%dt_therm) .and. (cs%dt_therm > 0.0)) &
759 n_max = ceiling(time_interval/cs%dt_therm - 0.001)
760
761 dt = time_interval / real(n_max)
762 dt_therm = dt ; ntstep = 1
763
764 if (cs%UseWaves .and. associated(fluxes%ustar)) &
765 call pass_var(fluxes%ustar, g%Domain, clock=id_clock_pass, halo=1)
766 if (cs%UseWaves .and. associated(fluxes%tau_mag)) &
767 call pass_var(fluxes%tau_mag, g%Domain, clock=id_clock_pass, halo=1)
768
769 if (associated(fluxes%p_surf)) p_surf => fluxes%p_surf
770 if (associated(cs%tv%p_surf) .and. associated(fluxes%p_surf)) then
771 do j=js,je ; do i=is,ie ; cs%tv%p_surf(i,j) = fluxes%p_surf(i,j) ; enddo ; enddo
772 if (allocated(cs%tv%SpV_avg)) then
773 call pass_var(cs%tv%p_surf, g%Domain, clock=id_clock_pass)
774 ! The internal ocean state depends on the surface pressues, so update SpV_avg.
775 call extract_diabatic_member(cs%diabatic_CSp, diabatic_halo=halo_sz)
776 halo_sz = max(halo_sz, 1)
777 call calc_derived_thermo(cs%tv, h, g, gv, us, halo=halo_sz, debug=cs%debug)
778 endif
779 endif
780 endif
781
782 if (therm_reset) then
783 cs%time_in_thermo_cycle = 0.0
784 if (associated(cs%tv%frazil)) then
785 cs%tv%frazil(:,:) = 0.0
786 cs%tv%frazil_was_reset = .true.
787 endif
788 if (associated(cs%tv%salt_deficit)) cs%tv%salt_deficit(:,:) = 0.0
789 if (associated(cs%tv%TempxPmE)) cs%tv%TempxPmE(:,:) = 0.0
790 if (associated(cs%tv%internal_heat)) cs%tv%internal_heat(:,:) = 0.0
791 endif
792
793 if (cycle_start) then
794 cs%time_in_cycle = 0.0
795 do j=js,je ; do i=is,ie ; cs%ssh_rint(i,j) = 0.0 ; enddo ; enddo
796
797 if (cs%VarMix%use_variable_mixing) then
798 time_end_diag = time_start + real_to_time(cycle_time, unscale=us%T_to_s)
799 call enable_averages(cycle_time, time_end_diag, cs%diag)
800 call calc_resoln_function(h, cs%tv, g, gv, us, cs%VarMix, cs%MEKE, cs%OBC, dt)
801 call calc_depth_function(g, cs%VarMix)
802 call disable_averaging(cs%diag)
803 endif
804 endif
805 ! advance the random pattern if stochastic physics is active
806 if (cs%stoch_CS%do_sppt .OR. cs%stoch_CS%pert_epbl .OR. cs%stoch_CS%do_skeb) &
807 call update_stochastics(cs%stoch_CS)
808
809 if (do_dyn) then
810 if (nonblocking_p_surf_update) &
811 call complete_group_pass(pass_tau_ustar_psurf, g%Domain, clock=id_clock_pass)
812
813 if (cs%interp_p_surf) then
814 if (.not.associated(cs%p_surf_end)) allocate(cs%p_surf_end(isd:ied,jsd:jed))
815 if (.not.associated(cs%p_surf_begin)) allocate(cs%p_surf_begin(isd:ied,jsd:jed))
816 if (.not.cs%p_surf_prev_set) then
817 do j=jsd,jed ; do i=isd,ied
818 cs%p_surf_prev(i,j) = forces%p_surf(i,j)
819 enddo ; enddo
820 cs%p_surf_prev_set = .true.
821 endif
822 else
823 cs%p_surf_end => forces%p_surf
824 endif
825 if (cs%UseWaves) then
826 ! Update wave information, which is presently kept static over each call to step_mom
827 time_end_diag = time_start + real_to_time(time_interval, unscale=us%T_to_s)
828 call enable_averages(time_interval, time_end_diag, cs%diag)
829 call find_ustar(forces, cs%tv, u_star, g, gv, us, halo=1)
830 call thickness_to_dz(h, cs%tv, dz, g, gv, us, halo_size=1)
831 call update_stokes_drift(g, gv, us, waves, dz, u_star, time_interval, do_dyn)
832 call disable_averaging(cs%diag)
833 endif
834 else ! not do_dyn.
835 if (cs%UseWaves) then ! Diagnostics are not enabled in this call.
836 call find_ustar(fluxes, cs%tv, u_star, g, gv, us, halo=1)
837 call thickness_to_dz(h, cs%tv, dz, g, gv, us, halo_size=1)
838 call update_stokes_drift(g, gv, us, waves, dz, u_star, time_interval, do_dyn)
839 endif
840 endif
841
842 if (cs%debug) then
843 if (cycle_start) &
844 call mom_state_chksum("Before steps ", u, v, h, cs%uh, cs%vh, g, gv, us)
845 if (cycle_start .and. debug_redundant) &
846 call check_redundant("Before steps ", u, v, g, unscale=us%L_T_to_m_s)
847 if (do_dyn) call mom_mech_forcing_chksum("Before steps", forces, g, us, haloshift=0)
848 if (do_dyn .and. debug_redundant) &
849 call check_redundant("Before steps ", forces%taux, forces%tauy, g, &
850 unscale=us%RZ_T_to_kg_m2s*us%L_T_to_m_s)
851 endif
852 call cpu_clock_end(id_clock_other)
853
854 rel_time = 0.0
855 do n=1,n_max
856 rel_time = rel_time + dt ! The relative time at the end of the step.
857 ! Set the universally visible time to the middle of the time step.
858 cs%Time = time_start + real_to_time(rel_time - 0.5*dt, unscale=us%T_to_s)
859 ! Set the local time to the end of the time step.
860 time_local = time_start + real_to_time(rel_time, unscale=us%T_to_s)
861
862 if (showcalltree) call calltree_enter("DT cycles (step_MOM) n=",n)
863
864 ! Update the vertically extensive diagnostic grids so that they are
865 ! referenced to the beginning timestep
866 call diag_update_remap_grids(cs%diag, update_intensive = .false., update_extensive = .true. )
867
868 !===========================================================================
869 ! This is the first place where the diabatic processes and remapping could occur.
870 if (cs%diabatic_first .and. (cs%t_dyn_rel_adv==0.0) .and. do_thermo) then ! do thermodynamics.
871
872 if (.not.do_dyn) then
873 dtdia = dt
874 elseif (thermo_does_span_coupling) then
875 dtdia = dt_therm
876 if ((fluxes%dt_buoy_accum > 0.0) .and. (dtdia > time_interval) .and. &
877 (abs(fluxes%dt_buoy_accum - dtdia) > 1e-6*dtdia)) then
878 call mom_error(fatal, "step_MOM: Mismatch between long thermodynamic "//&
879 "timestep and time over which buoyancy fluxes have been accumulated.")
880 endif
881 call mom_error(fatal, "MOM is not yet set up to have restarts that work "//&
882 "with THERMO_SPANS_COUPLING and DIABATIC_FIRST.")
883 else
884 dtdia = dt*min(ntstep,n_max-(n-1))
885 endif
886
887 end_time_thermo = time_local
888 if (dtdia > dt) then
889 ! If necessary, temporarily reset CS%Time to the center of the period covered
890 ! by the call to step_MOM_thermo, noting that they begin at the same time.
891 cs%Time = cs%Time + real_to_time(0.5*(dtdia-dt), unscale=us%T_to_s)
892 ! The end-time of the diagnostic interval needs to be set ahead if there
893 ! are multiple dynamic time steps worth of thermodynamics applied here.
894 end_time_thermo = time_local + real_to_time(dtdia-dt, unscale=us%T_to_s)
895 endif
896
897 ! Apply diabatic forcing, do mixing, and regrid.
898 call step_mom_thermo(cs, g, gv, us, u, v, h, cs%tv, fluxes, dtdia, &
899 end_time_thermo, .true., waves=waves)
900 if ( cs%use_ALE_algorithm ) &
901 call ale_regridding_and_remapping(cs, g, gv, us, u, v, h, cs%tv, dtdia, time_local)
902 call post_diabatic_halo_updates(cs, g, gv, us, u, v, h, cs%tv)
903 cs%time_in_thermo_cycle = cs%time_in_thermo_cycle + dtdia
904
905 ! The diabatic processes are now ahead of the dynamics by dtdia.
906 cs%t_dyn_rel_thermo = -dtdia
907 if (showcalltree) call calltree_waypoint("finished diabatic_first (step_MOM)")
908
909 if (dtdia > dt) & ! Reset CS%Time to its previous value.
910 cs%Time = time_start + real_to_time(rel_time - 0.5*dt, unscale=us%T_to_s)
911 endif ! end of block "(CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0))"
912
913 if (do_dyn) then
914 ! Store pre-dynamics thicknesses for proper diagnostic remapping for transports or
915 ! advective tendencies. If there are more than one dynamics steps per advective
916 ! step (i.e DT_THERM > DT), this needs to be stored at the first dynamics call.
917 if (.not.cs%preadv_h_stored .and. (cs%t_dyn_rel_adv == 0.)) then
918 call diag_copy_diag_to_storage(cs%diag_pre_dyn, h, cs%diag)
919 cs%preadv_h_stored = .true.
920 endif
921
922 ! The pre-dynamics velocities might be stored for debugging truncations.
923 if (associated(cs%u_prev) .and. associated(cs%v_prev)) then
924 do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb
925 cs%u_prev(i,j,k) = u(i,j,k)
926 enddo ; enddo ; enddo
927 do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
928 cs%v_prev(i,j,k) = v(i,j,k)
929 enddo ; enddo ; enddo
930 endif
931
932 if (cs%interface_filter_dt_bug) then
933 dt_tradv_here = dt_therm
934 if (do_thermo .and. do_dyn .and. .not.thermo_does_span_coupling) &
935 dt_tradv_here = dt*min(ntstep, n_max-n+1)
936 else
937 dt_tradv_here = dt_tr_adv
938 if (do_thermo .and. do_dyn .and. .not.tradv_does_span_coupling) &
939 dt_tradv_here = dt*min(ntstep, n_max-n+1)
940 endif
941
942 ! Indicate whether the bottom boundary layer properties need to be
943 ! recalculated, and if so for how long an interval they are valid.
944 bbl_time_int = 0.0
945 if (do_thermo) then
946 if ((cs%t_dyn_rel_adv == 0.0) .or. (n==1)) &
947 bbl_time_int = max(dt, min(dt_therm - cs%t_dyn_rel_adv, dt*(1+n_max-n)) )
948 else
949 if ((cs%t_dyn_rel_adv == 0.0) .or. ((n==1) .and. cycle_start)) &
950 bbl_time_int = min(dt_therm, cycle_time)
951 endif
952
953 if (cs%interp_p_surf) then
954 wt_end = real(n) / real(n_max)
955 wt_beg = real(n-1) / real(n_max)
956 do j=jsd,jed ; do i=isd,ied
957 cs%p_surf_end(i,j) = wt_end * forces%p_surf(i,j) + &
958 (1.0-wt_end) * cs%p_surf_prev(i,j)
959 cs%p_surf_begin(i,j) = wt_beg * forces%p_surf(i,j) + &
960 (1.0-wt_beg) * cs%p_surf_prev(i,j)
961 enddo ; enddo
962 endif
963
964 call step_mom_dynamics(forces, cs%p_surf_begin, cs%p_surf_end, dt, &
965 dt_tradv_here, bbl_time_int, cs, &
966 time_local, waves=waves)
967
968 !===========================================================================
969 ! This is the start of the tracer advection part of the algorithm.
970 if (tradv_does_span_coupling .or. .not.do_thermo) then
971 do_advection = ((cs%t_dyn_rel_adv + 0.5*dt > dt_tr_adv) .or. &
972 (cs%t_dyn_rel_thermo + 0.5*dt > dt_therm))
973 else
974 do_advection = ((mod(n,ntastep) == 0) .or. (n==n_max))
975 endif
976
977 if (do_advection) then ! Do advective transport and lateral tracer mixing.
978 call step_mom_tracer_dyn(cs, g, gv, us, h, time_local)
979 if (cs%diabatic_first .and. abs(cs%t_dyn_rel_thermo) > 1e-6*dt) call mom_error(fatal, &
980 "step_MOM: Mismatch between the dynamics and diabatic times "//&
981 "with DIABATIC_FIRST.")
982 endif
983 endif ! end of (do_dyn)
984
985 !===========================================================================
986 ! This is the second place where the diabatic processes and remapping could occur.
987 if (thermo_does_span_coupling .or. .not.do_dyn) then
988 do_diabatic = (do_thermo .and. (cs%t_dyn_rel_thermo + 0.5*dt > dt_therm))
989 else
990 do_diabatic = (do_thermo .and. ((mod(n,ntstep) == 0) .or. (n==n_max)))
991 endif
992 if ((cs%t_dyn_rel_adv==0.0) .and. (.not.cs%diabatic_first) .and. do_diabatic) then
993
994 dtdia = cs%t_dyn_rel_thermo
995 ! If the MOM6 dynamic and thermodynamic time stepping is being orchestrated
996 ! by the coupler, the value of diabatic_first does not matter.
997 if ((cs%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) dtdia = dt
998
999 if (cs%thermo_spans_coupling .and. (cs%dt_therm > 1.5*cycle_time) .and. &
1000 (abs(dt_therm - dtdia) > 1e-6*dt_therm)) then
1001 call mom_error(fatal, "step_MOM: Mismatch between dt_therm and dtdia "//&
1002 "before call to diabatic.")
1003 endif
1004
1005 ! If necessary, temporarily reset CS%Time to the center of the period covered
1006 ! by the call to step_MOM_thermo, noting that they end at the same time.
1007 if (dtdia > dt) &
1008 cs%Time = cs%Time - real_to_time(0.5*(dtdia-dt), unscale=us%T_to_s)
1009
1010 ! Apply diabatic forcing, do mixing, and regrid.
1011 call step_mom_thermo(cs, g, gv, us, u, v, h, cs%tv, fluxes, dtdia, &
1012 time_local, .false., waves=waves)
1013 if ( cs%use_ALE_algorithm ) &
1014 call ale_regridding_and_remapping(cs, g, gv, us, u, v, h, cs%tv, dtdia, time_local)
1015 call post_diabatic_halo_updates(cs, g, gv, us, u, v, h, cs%tv)
1016 cs%time_in_thermo_cycle = cs%time_in_thermo_cycle + dtdia
1017
1018 if ((cs%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) then
1019 ! The diabatic processes are now ahead of the dynamics by dtdia.
1020 cs%t_dyn_rel_thermo = -dtdia
1021 else ! The diabatic processes and the dynamics are synchronized.
1022 cs%t_dyn_rel_thermo = 0.0
1023 endif
1024
1025 ! Reset CS%Time to its previous value.
1026 if (dtdia > dt) &
1027 cs%Time = time_start + real_to_time(rel_time - 0.5*dt, unscale=us%T_to_s)
1028 endif
1029
1030 if (do_dyn) then
1031 call cpu_clock_begin(id_clock_dynamics)
1032 ! Determining the time-average sea surface height is part of the algorithm.
1033 ! This may be eta_av if Boussinesq, or need to be diagnosed if not.
1034 cs%time_in_cycle = cs%time_in_cycle + dt
1035 call find_eta(h, cs%tv, g, gv, us, ssh, cs%eta_av_bc, dzref=g%Z_ref)
1036 do j=js,je ; do i=is,ie
1037 cs%ssh_rint(i,j) = cs%ssh_rint(i,j) + dt*ssh(i,j)
1038 enddo ; enddo
1039 if (cs%IDs%id_ssh_inst > 0) then
1040 call enable_averages(dt, time_local, cs%diag)
1041 call post_data(cs%IDs%id_ssh_inst, ssh, cs%diag)
1042 call disable_averaging(cs%diag)
1043 endif
1044 call cpu_clock_end(id_clock_dynamics)
1045 endif
1046
1047 !===========================================================================
1048 ! Calculate diagnostics at the end of the time step if the state is self-consistent.
1049 if (mom_state_is_synchronized(cs)) then
1050 !### Perhaps this should be if (CS%t_dyn_rel_thermo == 0.0)
1051 call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
1052 ! Diagnostics that require the complete state to be up-to-date can be calculated.
1053
1054 call enable_averages(cs%t_dyn_rel_diag, time_local, cs%diag)
1055 call calculate_diagnostic_fields(u, v, h, cs%uh, cs%vh, cs%tv, cs%ADp, &
1056 cs%CDp, p_surf, cs%t_dyn_rel_diag, cs%diag_pre_sync,&
1057 g, gv, us, cs%diagnostics_CSp)
1058 call post_tracer_diagnostics_at_sync(cs%Tracer_reg, h, cs%diag_pre_sync, cs%diag, g, gv, cs%t_dyn_rel_diag)
1059 call diag_copy_diag_to_storage(cs%diag_pre_sync, h, cs%diag)
1060 if (showcalltree) call calltree_waypoint("finished calculate_diagnostic_fields (step_MOM)")
1061 call disable_averaging(cs%diag)
1062 cs%t_dyn_rel_diag = 0.0
1063
1064 call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
1065 endif
1066
1067 if (do_dyn .and. .not.cs%count_calls) cs%nstep_tot = cs%nstep_tot + 1
1068 if (showcalltree) call calltree_leave("DT cycles (step_MOM)")
1069
1070 enddo ! complete the n loop
1071
1072 if (cs%count_calls .and. cycle_start) cs%nstep_tot = cs%nstep_tot + 1
1073
1074 call cpu_clock_begin(id_clock_other)
1075
1076 if (cs%time_in_cycle > 0.0) then
1077 i_wt_ssh = 1.0/cs%time_in_cycle
1078 do j=js,je ; do i=is,ie
1079 ssh(i,j) = cs%ssh_rint(i,j)*i_wt_ssh
1080 cs%ave_ssh_ibc(i,j) = ssh(i,j)
1081 enddo ; enddo
1082 if (associated(cs%HA_CSp)) call ha_accum('ssh', ssh, time_local, g, cs%HA_CSp)
1083 if (do_dyn) then
1084 call adjust_ssh_for_p_atm(cs%tv, g, gv, us, cs%ave_ssh_ibc, forces%p_surf_SSH, &
1085 cs%calc_rho_for_sea_lev)
1086 elseif (do_thermo) then
1087 call adjust_ssh_for_p_atm(cs%tv, g, gv, us, cs%ave_ssh_ibc, fluxes%p_surf_SSH, &
1088 cs%calc_rho_for_sea_lev)
1089 endif
1090 endif
1091
1092 if (do_dyn .and. cs%interp_p_surf) then ; do j=jsd,jed ; do i=isd,ied
1093 cs%p_surf_prev(i,j) = forces%p_surf(i,j)
1094 enddo ; enddo ; endif
1095
1096 if (cs%ensemble_ocean) then
1097 ! store ensemble vector in odaCS
1098 call set_prior_tracer(cs%Time, g, gv, cs%h, cs%tv, cs%odaCS)
1099 ! call DA interface
1100 call oda(cs%Time,cs%odaCS)
1101 ! update the time for the next analysis step if needed
1102 call set_analysis_time(cs%Time,cs%odaCS)
1103 endif
1104
1105 if (showcalltree) call calltree_waypoint("calling extract_surface_state (step_MOM)")
1106 ! NOTE: sfc_state uses input indexing, since it is also used by drivers.
1107 call extract_surface_state(cs, sfc_state)
1108
1109 ! Do diagnostics that only occur at the end of a complete forcing step.
1110 if (cycle_end) then
1111 if (showcalltree) call calltree_waypoint("Do cycle end diagnostics (step_MOM)")
1112 if (cs%rotate_index) then
1113 allocate(sfc_state_diag)
1114 call rotate_surface_state(sfc_state, sfc_state_diag, g, turns)
1115 else
1116 sfc_state_diag => sfc_state
1117 endif
1118
1119 call cpu_clock_begin(id_clock_diagnostics)
1120 if (cs%time_in_cycle > 0.0) then
1121 call enable_averages(cs%time_in_cycle, time_local, cs%diag)
1122 call post_surface_dyn_diags(cs%sfc_IDs, g, cs%diag, sfc_state_diag, ssh)
1123 endif
1124 if (cs%time_in_thermo_cycle > 0.0) then
1125 call enable_averages(cs%time_in_thermo_cycle, time_local, cs%diag)
1126 call post_surface_thermo_diags(cs%sfc_IDs, g, gv, us, cs%diag, cs%time_in_thermo_cycle, &
1127 sfc_state_diag, cs%tv, ssh, cs%ave_ssh_ibc)
1128 endif
1129 call disable_averaging(cs%diag)
1130 call cpu_clock_end(id_clock_diagnostics)
1131 if (cs%rotate_index) then
1132 call deallocate_surface_state(sfc_state_diag)
1133 endif
1134 if (showcalltree) call calltree_waypoint("Done with end cycle diagnostics (step_MOM)")
1135 endif
1136
1137 ! Accumulate the surface fluxes for assessing conservation
1138 if (do_thermo .and. fluxes%fluxes_used) &
1139 call accumulate_net_input(fluxes, sfc_state, cs%tv, fluxes%dt_buoy_accum, &
1140 g, us, cs%sum_output_CSp)
1141
1142 if (mom_state_is_synchronized(cs)) &
1143 call write_energy(cs%u, cs%v, cs%h, cs%tv, time_local, cs%nstep_tot, &
1144 g, gv, us, cs%sum_output_CSp, cs%tracer_flow_CSp, &
1145 dt_forcing=real_to_time(time_interval, unscale=us%T_to_s) )
1146
1147 call cpu_clock_end(id_clock_other)
1148
1149 ! De-rotate fluxes and copy back to the input, since they can be changed.
1150 if (cs%rotate_index) then
1151 call rotate_forcing(fluxes, fluxes_in, -turns)
1152 call rotate_mech_forcing(forces, -turns, forces_in)
1153 call deallocate_mech_forcing(forces)
1154 deallocate(forces)
1155 call deallocate_forcing_type(fluxes)
1156 deallocate(fluxes)
1157 endif
1158
1159 if (showcalltree) call calltree_leave("step_MOM()")
1160 call cpu_clock_end(id_clock_ocean)
1161
1162end subroutine step_mom
1163
1164!> Time step the ocean dynamics, including the momentum and continuity equations
1165subroutine step_mom_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_tr_adv, &
1166 bbl_time_int, CS, Time_local, Waves)
1167 type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1168 real, dimension(:,:), pointer :: p_surf_begin !< A pointer (perhaps NULL) to the surface
1169 !! pressure at the beginning of this dynamic
1170 !! step, intent in [R L2 T-2 ~> Pa].
1171 real, dimension(:,:), pointer :: p_surf_end !< A pointer (perhaps NULL) to the surface
1172 !! pressure at the end of this dynamic step,
1173 !! intent in [R L2 T-2 ~> Pa].
1174 real, intent(in) :: dt !< time interval covered by this call [T ~> s].
1175 real, intent(in) :: dt_tr_adv !< time interval covered by any updates that may
1176 !! span multiple dynamics steps [T ~> s].
1177 real, intent(in) :: bbl_time_int !< time interval over which updates to the
1178 !! bottom boundary layer properties will apply [T ~> s],
1179 !! or zero not to update the properties.
1180 type(mom_control_struct), intent(inout), target :: CS !< control structure from initialize_MOM
1181 type(time_type), intent(in) :: Time_local !< End time of a segment, as a time type
1182 type(wave_parameters_cs), &
1183 optional, pointer :: Waves !< Container for wave related parameters; the
1184 !! fields in Waves are intent in here.
1185
1186 ! local variables
1187 type(ocean_grid_type), pointer :: G => null() ! pointer to a structure containing
1188 ! metrics and related information
1189 type(verticalgrid_type), pointer :: GV => null() ! Pointer to the vertical grid structure
1190 type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing
1191 ! various unit conversion factors
1192 type(mom_diag_ids), pointer :: IDs => null() ! A structure with the diagnostic IDs.
1193 real, dimension(:,:,:), pointer :: &
1194 u => null(), & ! u : zonal velocity component [L T-1 ~> m s-1]
1195 v => null(), & ! v : meridional velocity component [L T-1 ~> m s-1]
1196 h => null() ! h : layer thickness [H ~> m or kg m-2]
1197
1198 type(time_type) :: Time_end_diag ! End time of a diagnostic segment, as a time type
1199 logical :: calc_dtbt ! Indicates whether the dynamically adjusted
1200 ! barotropic time step needs to be updated.
1201 logical :: showCallTree
1202
1203 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1204 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
1205
1206 g => cs%G ; gv => cs%GV ; us => cs%US ; ids => cs%IDs
1207 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1208 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1209 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1210 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1211 u => cs%u ; v => cs%v ; h => cs%h
1212 showcalltree = calltree_showquery()
1213
1214 call cpu_clock_begin(id_clock_dynamics)
1215 call cpu_clock_begin(id_clock_stoch)
1216 if (cs%use_stochastic_EOS) call mom_stoch_eos_run(g, u, v, dt, time_local, cs%stoch_eos_CS)
1217 call cpu_clock_end(id_clock_stoch)
1218 call cpu_clock_begin(id_clock_vart)
1219 if (cs%use_stochastic_EOS) then
1220 call mom_calc_vart(g, gv, us, h, cs%tv, cs%stoch_eos_CS, dt)
1221 if (associated(cs%tv%varT)) call pass_var(cs%tv%varT, g%Domain, clock=id_clock_pass, halo=1)
1222 endif
1223 call cpu_clock_end(id_clock_vart)
1224
1225 if ((cs%t_dyn_rel_adv == 0.0) .and. cs%thickness_diffuse_first .and. &
1226 (cs%thickness_diffuse .or. cs%interface_filter)) then
1227
1228 time_end_diag = time_local + real_to_time(dt_tr_adv - dt, unscale=us%T_to_s)
1229 call enable_averages(dt_tr_adv, time_end_diag, cs%diag)
1230 if (cs%thickness_diffuse) then
1231 call cpu_clock_begin(id_clock_thick_diff)
1232 if (cs%VarMix%use_variable_mixing) &
1233 call calc_slope_functions(h, cs%tv, dt, g, gv, us, cs%VarMix, obc=cs%OBC)
1234 call thickness_diffuse(h, cs%uhtr, cs%vhtr, cs%tv, dt_tr_adv, g, gv, us, &
1235 cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp, &
1236 cs%stoch_CS, u, v)
1237 call cpu_clock_end(id_clock_thick_diff)
1238 call pass_var(h, g%Domain, clock=id_clock_pass, halo=cs%dyn_h_stencil)
1239 if (showcalltree) call calltree_waypoint("finished thickness_diffuse_first (step_MOM)")
1240 endif
1241
1242 if (cs%interface_filter) then
1243 if (allocated(cs%tv%SpV_avg)) call pass_var(cs%tv%SpV_avg, g%Domain, clock=id_clock_pass)
1244 cs%tv%valid_SpV_halo = min(g%Domain%nihalo, g%Domain%njhalo)
1245 call cpu_clock_begin(id_clock_int_filter)
1246 call interface_filter(h, cs%uhtr, cs%vhtr, cs%tv, dt_tr_adv, g, gv, us, &
1247 cs%CDp, cs%interface_filter_CSp)
1248 call cpu_clock_end(id_clock_int_filter)
1249 call pass_var(h, g%Domain, clock=id_clock_pass, halo=cs%dyn_h_stencil)
1250 if (showcalltree) call calltree_waypoint("finished interface_filter_first (step_MOM)")
1251 endif
1252
1253 call disable_averaging(cs%diag)
1254 ! Whenever thickness changes let the diag manager know, target grids
1255 ! for vertical remapping may need to be regenerated.
1256 call diag_update_remap_grids(cs%diag)
1257 endif
1258
1259 ! Update porous barrier fractional cell metrics
1260 if (cs%use_porbar) then
1261 call enable_averages(dt, time_local, cs%diag)
1262 call porous_widths_layer(h, cs%tv, g, gv, us, cs%pbv, cs%por_bar_CS)
1263 call disable_averaging(cs%diag)
1264 call pass_vector(cs%pbv%por_face_areaU, cs%pbv%por_face_areaV, &
1265 g%Domain, direction=to_all+scalar_pair, clock=id_clock_pass, halo=cs%cont_stencil)
1266 endif
1267
1268 ! The bottom boundary layer properties need to be recalculated.
1269 if (bbl_time_int > 0.0) then
1270 time_end_diag = time_local + real_to_time(bbl_time_int - dt, unscale=us%T_to_s)
1271 call enable_averages(bbl_time_int, time_end_diag, cs%diag)
1272 ! Calculate the BBL properties and store them inside visc (u,h).
1273 call cpu_clock_begin(id_clock_bbl_visc)
1274 call set_viscous_bbl(cs%u, cs%v, cs%h, cs%tv, cs%visc, g, gv, us, cs%set_visc_CSp, cs%pbv)
1275 call cpu_clock_end(id_clock_bbl_visc)
1276 if (showcalltree) call calltree_waypoint("done with set_viscous_BBL (step_MOM)")
1277 call disable_averaging(cs%diag)
1278 endif
1279
1280 !OBC segment data update for some fields can be less frequent than others
1281 if (associated(cs%OBC)) then
1282 cs%OBC%update_OBC_seg_data = .false.
1283 if (cs%dt_obc_seg_period == 0.0) cs%OBC%update_OBC_seg_data = .true.
1284 if (cs%dt_obc_seg_period > 0.0) then
1285 if (time_local >= cs%dt_obc_seg_time) then
1286 cs%OBC%update_OBC_seg_data = .true.
1287 cs%dt_obc_seg_time = cs%dt_obc_seg_time + cs%dt_obc_seg_interval
1288 endif
1289 endif
1290 endif
1291 ! if (CS%debug_OBCs .and. associated(CS%OBC)) call chksum_OBC_segments(CS%OBC, G, GV, US, 3)
1292
1293 if (cs%do_dynamics .and. cs%split) then !--------------------------- start SPLIT
1294 ! This section uses a split time stepping scheme for the dynamic equations,
1295 ! basically the stacked shallow water equations with viscosity.
1296
1297 calc_dtbt = .false.
1298 if (cs%dtbt_reset_period == 0.0) calc_dtbt = .true.
1299 if (cs%dtbt_reset_period > 0.0) then
1300 if (time_local >= cs%dtbt_reset_time) then !### Change >= to > here.
1301 calc_dtbt = .true.
1302 cs%dtbt_reset_time = cs%dtbt_reset_time + cs%dtbt_reset_interval
1303 endif
1304 endif
1305
1306 if (cs%use_alt_split) then
1307 call step_mom_dyn_split_rk2b(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
1308 p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
1309 cs%eta_av_bc, g, gv, us, cs%dyn_split_RK2b_CSp, calc_dtbt, cs%VarMix, &
1310 cs%MEKE, cs%thickness_diffuse_CSp, cs%pbv, waves=waves)
1311 else
1312 call step_mom_dyn_split_rk2(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
1313 p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
1314 cs%eta_av_bc, g, gv, us, cs%dyn_split_RK2_CSp, calc_dtbt, cs%VarMix, &
1315 cs%MEKE, cs%thickness_diffuse_CSp, cs%pbv, cs%stoch_CS, waves=waves)
1316 endif
1317 if (showcalltree) call calltree_waypoint("finished step_MOM_dyn_split (step_MOM)")
1318
1319 elseif (cs%do_dynamics) then ! ------------------------------------ not SPLIT
1320 ! This section uses an unsplit stepping scheme for the dynamic
1321 ! equations; basically the stacked shallow water equations with viscosity.
1322 ! Because the time step is limited by CFL restrictions on the external
1323 ! gravity waves, the unsplit is usually much less efficient that the split
1324 ! approaches. But because of its simplicity, the unsplit method is very
1325 ! useful for debugging purposes.
1326
1327 if (cs%use_RK2) then
1328 call step_mom_dyn_unsplit_rk2(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
1329 p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
1330 cs%eta_av_bc, g, gv, us, cs%dyn_unsplit_RK2_CSp, cs%VarMix, cs%MEKE, cs%pbv, &
1331 cs%stoch_CS)
1332 else
1333 call step_mom_dyn_unsplit(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
1334 p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
1335 cs%eta_av_bc, g, gv, us, cs%dyn_unsplit_CSp, cs%VarMix, cs%MEKE, cs%pbv, &
1336 cs%stoch_CS, waves=waves)
1337 endif
1338 if (showcalltree) call calltree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)")
1339
1340 endif ! -------------------------------------------------- end SPLIT
1341
1342 if (cs%use_particles .and. cs%do_dynamics .and. (.not. cs%use_uh_particles)) then
1343 if (cs%thickness_diffuse_first) call mom_error(warning,"particles_run: "//&
1344 "Thickness_diffuse_first is true and use_uh_particles is false. "//&
1345 "This is usually a bad combination.")
1346 !Run particles using unweighted velocity
1347 call particles_run(cs%particles, time_local, cs%u, cs%v, cs%h, &
1348 cs%tv, dt, cs%use_uh_particles)
1349 call particles_to_z_space(cs%particles, h)
1350 endif
1351
1352 ! Update the model's current to reflect wind-wave growth
1353 if (waves%Stokes_DDT .and. (.not.waves%Passive_Stokes_DDT)) then
1354 do j=jsq,jeq ; do i=is,ie
1355 v(i,j,:) = v(i,j,:) + waves%ddt_us_y(i,j,:)*dt
1356 enddo ; enddo
1357 do j=js,je ; do i=isq,ieq
1358 u(i,j,:) = u(i,j,:) + waves%ddt_us_x(i,j,:)*dt
1359 enddo ; enddo
1360 call pass_vector(u, v, g%Domain)
1361 endif
1362 ! Added an additional output to track Stokes drift time tendency.
1363 ! It is mostly for debugging, and perhaps doesn't need to hang
1364 ! around permanently.
1365 if (waves%Stokes_DDT .and. (waves%id_3dstokes_y_from_ddt>0)) then
1366 do j=jsq,jeq ; do i=is,ie
1367 waves%us_y_from_ddt(i,j,:) = waves%us_y_from_ddt(i,j,:) + waves%ddt_us_y(i,j,:)*dt
1368 enddo ; enddo
1369 endif
1370 if (waves%Stokes_DDT .and. (waves%id_3dstokes_x_from_ddt>0)) then
1371 do j=js,je ; do i=isq,ieq
1372 waves%us_x_from_ddt(i,j,:) = waves%us_x_from_ddt(i,j,:) + waves%ddt_us_x(i,j,:)*dt
1373 enddo ; enddo
1374 endif
1375
1376
1377 if ((cs%thickness_diffuse .or. cs%interface_filter) .and. &
1378 .not.cs%thickness_diffuse_first) then
1379
1380 if (cs%debug) call hchksum(h,"Pre-thickness_diffuse h", g%HI, haloshift=0, unscale=gv%H_to_MKS)
1381
1382 if (cs%thickness_diffuse) then
1383 call cpu_clock_begin(id_clock_thick_diff)
1384 if (cs%VarMix%use_variable_mixing) &
1385 call calc_slope_functions(h, cs%tv, dt, g, gv, us, cs%VarMix, obc=cs%OBC)
1386 call thickness_diffuse(h, cs%uhtr, cs%vhtr, cs%tv, dt, g, gv, us, &
1387 cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp, &
1388 cs%stoch_CS, u, v)
1389 call cpu_clock_end(id_clock_thick_diff)
1390 call pass_var(h, g%Domain, clock=id_clock_pass, halo=cs%dyn_h_stencil)
1391 if (cs%debug) call hchksum(h,"Post-thickness_diffuse h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
1392 if (showcalltree) call calltree_waypoint("finished thickness_diffuse (step_MOM)")
1393 endif
1394
1395 if (cs%interface_filter) then
1396 if (allocated(cs%tv%SpV_avg)) call pass_var(cs%tv%SpV_avg, g%Domain, clock=id_clock_pass)
1397 cs%tv%valid_SpV_halo = min(g%Domain%nihalo, g%Domain%njhalo)
1398 call cpu_clock_begin(id_clock_int_filter)
1399 if (cs%interface_filter_dt_bug) then
1400 call interface_filter(h, cs%uhtr, cs%vhtr, cs%tv, dt_tr_adv, g, gv, us, &
1401 cs%CDp, cs%interface_filter_CSp)
1402 else
1403 call interface_filter(h, cs%uhtr, cs%vhtr, cs%tv, dt, g, gv, us, &
1404 cs%CDp, cs%interface_filter_CSp)
1405 endif
1406 call cpu_clock_end(id_clock_int_filter)
1407 call pass_var(h, g%Domain, clock=id_clock_pass, halo=cs%dyn_h_stencil)
1408 if (showcalltree) call calltree_waypoint("finished interface_filter (step_MOM)")
1409 endif
1410 endif
1411
1412 ! apply the submesoscale mixed layer restratification parameterization
1413 if (cs%mixedlayer_restrat) then
1414 if (cs%debug) then
1415 call hchksum(h,"Pre-mixedlayer_restrat h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
1416 call uvchksum("Pre-mixedlayer_restrat uhtr", &
1417 cs%uhtr, cs%vhtr, g%HI, haloshift=0, unscale=gv%H_to_MKS*us%L_to_m**2)
1418 endif
1419 call cpu_clock_begin(id_clock_ml_restrat)
1420 call mixedlayer_restrat(h, cs%uhtr, cs%vhtr, cs%tv, forces, dt, cs%visc%MLD, cs%visc%h_ML, &
1421 cs%visc%sfc_buoy_flx, cs%VarMix, g, gv, us, cs%mixedlayer_restrat_CSp)
1422 call cpu_clock_end(id_clock_ml_restrat)
1423 call pass_var(h, g%Domain, clock=id_clock_pass, halo=cs%dyn_h_stencil)
1424 if (cs%debug) then
1425 call hchksum(h,"Post-mixedlayer_restrat h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
1426 call uvchksum("Post-mixedlayer_restrat [uv]htr", &
1427 cs%uhtr, cs%vhtr, g%HI, haloshift=0, unscale=gv%H_to_MKS*us%L_to_m**2)
1428 endif
1429 endif
1430
1431 ! Whenever thickness changes let the diag manager know, target grids
1432 ! for vertical remapping may need to be regenerated.
1433 call diag_update_remap_grids(cs%diag)
1434
1435 if (cs%useMEKE .and. cs%MEKE_in_dynamics) then
1436 call step_forward_meke(cs%MEKE, h, cs%VarMix%SN_u, cs%VarMix%SN_v, &
1437 cs%visc, dt, g, gv, us, cs%MEKE_CSp, cs%uhtr, cs%vhtr, &
1438 cs%u, cs%v, cs%tv, time_local)
1439 endif
1440 call disable_averaging(cs%diag)
1441
1442 ! Advance the dynamics time by dt.
1443 cs%t_dyn_rel_adv = cs%t_dyn_rel_adv + dt
1444
1445 if (cs%use_particles .and. cs%do_dynamics .and. cs%use_uh_particles .and. &
1446 cs%uh_particles_bug) then
1447 ! Run particles using thickness-weighted velocity
1448 call particles_run(cs%particles, time_local, cs%uhtr, cs%vhtr, cs%h, &
1449 cs%tv, cs%t_dyn_rel_adv, cs%use_uh_particles)
1450 endif
1451
1452 cs%n_dyn_steps_in_adv = cs%n_dyn_steps_in_adv + 1
1453 if (cs%alternate_first_direction) then
1454 call set_first_direction(g, modulo(g%first_direction+1,2))
1455 cs%first_dir_restart = real(g%first_direction)
1456 elseif (cs%use_particles .and. cs%do_dynamics .and. (.not.cs%use_uh_particles)) then
1457 call particles_to_k_space(cs%particles, h)
1458 endif
1459 cs%t_dyn_rel_thermo = cs%t_dyn_rel_thermo + dt
1460 if (abs(cs%t_dyn_rel_thermo) < 1e-6*dt) cs%t_dyn_rel_thermo = 0.0
1461 cs%t_dyn_rel_diag = cs%t_dyn_rel_diag + dt
1462
1463 call cpu_clock_end(id_clock_dynamics)
1464
1465 call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
1466 call enable_averages(dt, time_local, cs%diag)
1467 ! These diagnostics are available after every time dynamics step.
1468 if (ids%id_u > 0) call post_data(ids%id_u, u, cs%diag)
1469 if (ids%id_v > 0) call post_data(ids%id_v, v, cs%diag)
1470 if (ids%id_h > 0) call post_data(ids%id_h, h, cs%diag)
1471 if (cs%use_stochastic_EOS) call post_stoch_eos_diags(cs%stoch_eos_CS, cs%tv, cs%diag)
1472 call disable_averaging(cs%diag)
1473 call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
1474
1475end subroutine step_mom_dynamics
1476
1477!> step_MOM_tracer_dyn does tracer advection and lateral diffusion, bringing the
1478!! tracers up to date with the changes in state due to the dynamics. Surface
1479!! sources and sinks and remapping are handled via step_MOM_thermo.
1480subroutine step_mom_tracer_dyn(CS, G, GV, US, h, Time_local)
1481 type(mom_control_struct), intent(inout) :: CS !< control structure
1482 type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1483 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
1484 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1485 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1486 intent(in) :: h !< layer thicknesses after the transports [H ~> m or kg m-2]
1487 type(time_type), intent(in) :: Time_local !< The model time at the end
1488 !! of the time step.
1489 type(group_pass_type) :: pass_T_S
1490 integer :: halo_sz ! The size of a halo where data must be valid.
1491 logical :: x_first ! If true, advect tracers first in the x-direction, then y.
1492 logical :: showCallTree
1493 showcalltree = calltree_showquery()
1494
1495 if (cs%debug) then
1496 call cpu_clock_begin(id_clock_other)
1497 call hchksum(h,"Pre-advection h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
1498 call uvchksum("Pre-advection uhtr", cs%uhtr, cs%vhtr, g%HI, &
1499 haloshift=0, unscale=gv%H_to_MKS*us%L_to_m**2)
1500 if (associated(cs%tv%T)) call hchksum(cs%tv%T, "Pre-advection T", g%HI, haloshift=1, unscale=us%C_to_degC)
1501 if (associated(cs%tv%S)) call hchksum(cs%tv%S, "Pre-advection S", g%HI, haloshift=1, unscale=us%S_to_ppt)
1502 if (associated(cs%tv%frazil)) call hchksum(cs%tv%frazil, "Pre-advection frazil", g%HI, haloshift=0, &
1503 unscale=us%Q_to_J_kg*us%RZ_to_kg_m2)
1504 if (associated(cs%tv%salt_deficit)) call hchksum(cs%tv%salt_deficit, &
1505 "Pre-advection salt deficit", g%HI, haloshift=0, unscale=us%S_to_ppt*us%RZ_to_kg_m2)
1506 ! call MOM_thermo_chksum("Pre-advection ", CS%tv, G, US)
1507 call cpu_clock_end(id_clock_other)
1508 endif
1509
1510 call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer)
1511 call enable_averages(cs%t_dyn_rel_adv, time_local, cs%diag)
1512
1513 if (cs%use_particles .and. cs%use_uh_particles .and. (.not. cs%uh_particles_bug)) then
1514 ! Run particles using thickness-weighted velocity
1515 call particles_run(cs%particles, time_local, cs%uhtr, cs%vhtr, cs%h, &
1516 cs%tv, cs%t_dyn_rel_adv, cs%use_uh_particles)
1517 endif
1518
1519
1520 if (cs%alternate_first_direction) then
1521 ! This calculation of the value of G%first_direction from the start of the accumulation of
1522 ! mass transports for use by the tracers is the equivalent to adding 2*n_dyn_steps before
1523 ! subtracting n_dyn_steps so that the mod will be taken of a non-negative number.
1524 x_first = (modulo(g%first_direction+cs%n_dyn_steps_in_adv,2) == 0)
1525 else
1526 x_first = (modulo(g%first_direction,2) == 0)
1527 endif
1528 if (cs%debug) call mom_tracer_chksum("Pre-advect ", cs%tracer_Reg, g)
1529 call advect_tracer(h, cs%uhtr, cs%vhtr, cs%OBC, cs%t_dyn_rel_adv, g, gv, us, &
1530 cs%tracer_adv_CSp, cs%tracer_Reg, x_first_in=x_first)
1531 if (cs%debug) call mom_tracer_chksum("Post-advect ", cs%tracer_Reg, g)
1532 call tracer_hordiff(h, cs%t_dyn_rel_adv, cs%MEKE, cs%VarMix, cs%visc, g, gv, us, &
1533 cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1534 if (cs%debug) call mom_tracer_chksum("Post-diffuse ", cs%tracer_Reg, g)
1535 if (showcalltree) call calltree_waypoint("finished tracer advection/diffusion (step_MOM)")
1536 if (associated(cs%OBC)) then
1537 call pass_vector(cs%uhtr, cs%vhtr, g%Domain)
1538 call update_segment_tracer_reservoirs(g, gv, cs%uhtr, cs%vhtr, h, cs%OBC, &
1539 cs%tracer_Reg)
1540 endif
1541 call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo)
1542
1543 call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
1544 call post_transport_diagnostics(g, gv, us, cs%uhtr, cs%vhtr, h, cs%transport_IDs, &
1545 cs%diag_pre_dyn, cs%diag, cs%t_dyn_rel_adv, cs%tracer_reg)
1546 ! Rebuild the remap grids now that we've posted the fields which rely on thicknesses
1547 ! from before the dynamics calls
1548 call diag_update_remap_grids(cs%diag)
1549
1550 call disable_averaging(cs%diag)
1551 call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
1552
1553 ! Reset the accumulated transports to 0 and record that the dynamics
1554 ! and advective times now agree.
1555 call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer)
1556 cs%uhtr(:,:,:) = 0.0
1557 cs%vhtr(:,:,:) = 0.0
1558 cs%n_dyn_steps_in_adv = 0
1559 cs%t_dyn_rel_adv = 0.0
1560 call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo)
1561
1562 if (cs%useMEKE .and. (.not. cs%MEKE_in_dynamics)) then
1563 call step_forward_meke(cs%MEKE, h, cs%VarMix%SN_u, cs%VarMix%SN_v, &
1564 cs%visc, cs%t_dyn_rel_adv, g, gv, us, cs%MEKE_CSp, cs%uhtr, cs%vhtr, &
1565 cs%u, cs%v, cs%tv, time_local)
1566 endif
1567
1568 if (associated(cs%tv%T)) then
1569 call extract_diabatic_member(cs%diabatic_CSp, diabatic_halo=halo_sz)
1570 ! The bottom boundary layer calculation may need halo values of SpV_avg, including the corners.
1571 if (allocated(cs%tv%SpV_avg)) halo_sz = max(halo_sz, 1)
1572 if (halo_sz > 0) then
1573 call create_group_pass(pass_t_s, cs%tv%T, g%Domain, to_all, halo=halo_sz)
1574 call create_group_pass(pass_t_s, cs%tv%S, g%Domain, to_all, halo=halo_sz)
1575 call do_group_pass(pass_t_s, g%Domain, clock=id_clock_pass)
1576 elseif (cs%diabatic_first) then
1577 ! Temperature and salinity need halo updates because they will be used
1578 ! in the dynamics before they are changed again.
1579 call create_group_pass(pass_t_s, cs%tv%T, g%Domain, to_all+omit_corners, halo=1)
1580 call create_group_pass(pass_t_s, cs%tv%S, g%Domain, to_all+omit_corners, halo=1)
1581 call do_group_pass(pass_t_s, g%Domain, clock=id_clock_pass)
1582 halo_sz = 1
1583 endif
1584
1585 ! Update derived thermodynamic quantities.
1586 if (allocated(cs%tv%SpV_avg)) then
1587 call calc_derived_thermo(cs%tv, h, g, gv, us, halo=halo_sz, debug=cs%debug)
1588 endif
1589 endif
1590
1591 cs%preadv_h_stored = .false.
1592
1593end subroutine step_mom_tracer_dyn
1594
1595!> MOM_step_thermo orchestrates the thermodynamic time stepping and vertical
1596!! remapping, via calls to diabatic (or adiabatic).
1597subroutine step_mom_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
1598 Time_end_thermo, update_BBL, Waves)
1599 type(mom_control_struct), intent(inout) :: CS !< Master MOM control structure
1600 type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1601 type(verticalgrid_type), intent(inout) :: GV !< ocean vertical grid structure
1602 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1603 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1604 intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1605 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1606 intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
1607 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1608 intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
1609 type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
1610 type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
1611 real, intent(in) :: dtdia !< The time interval over which to advance [T ~> s]
1612 type(time_type), intent(in) :: Time_end_thermo !< End of averaging interval for thermo diags
1613 logical, intent(in) :: update_BBL !< If true, calculate the bottom boundary layer properties.
1614 type(wave_parameters_cs), &
1615 optional, pointer :: Waves !< Container for wave related parameters
1616 !! the fields in Waves are intent in here.
1617
1618 logical :: debug_redundant ! If true, check redundant values on PE boundaries when debugging.
1619 logical :: showCallTree
1620 type(group_pass_type) :: pass_T_S
1621 integer :: dynamics_stencil ! The computational stencil for the calculations
1622 ! in the dynamic core.
1623 integer :: halo_sz ! The size of a halo where data must be valid.
1624
1625 showcalltree = calltree_showquery()
1626 if (showcalltree) call calltree_enter("step_MOM_thermo(), MOM.F90")
1627 if (cs%debug) call query_debugging_checks(do_redundant=debug_redundant)
1628
1629 call enable_averages(dtdia, time_end_thermo, cs%diag)
1630
1631 if (associated(cs%odaCS)) then
1632 if (cs%debug) then
1633 call mom_thermo_chksum("Pre-oda ", tv, g, us, haloshift=0)
1634 endif
1635 call apply_oda_tracer_increments(dtdia, time_end_thermo, g, gv, tv, h, cs%odaCS)
1636 if (cs%debug) then
1637 call mom_thermo_chksum("Post-oda ", tv, g, us, haloshift=0)
1638 endif
1639 endif
1640
1641 if (associated(fluxes%p_surf) .or. associated(fluxes%p_surf_full)) then
1642 call extract_diabatic_member(cs%diabatic_CSp, diabatic_halo=halo_sz)
1643 if (halo_sz > 0) then
1644 if (associated(fluxes%p_surf_full)) &
1645 call pass_var(fluxes%p_surf_full, g%Domain, &
1646 clock=id_clock_pass, halo=halo_sz, complete=.not.associated(fluxes%p_surf))
1647 call pass_var(fluxes%p_surf, g%Domain, clock=id_clock_pass, halo=halo_sz, complete=.true.)
1648 endif
1649 endif
1650
1651 if (update_bbl) then
1652 ! Calculate the BBL properties and store them inside visc (u,h).
1653 ! This is here so that CS%visc is updated before diabatic() when
1654 ! DIABATIC_FIRST=True. Otherwise diabatic() is called after the dynamics
1655 ! and set_viscous_BBL is called as a part of the dynamic stepping.
1656 call cpu_clock_begin(id_clock_bbl_visc)
1657 !update porous barrier fractional cell metrics
1658 if (cs%use_porbar) then
1659 call porous_widths_interface(h, cs%tv, g, gv, us, cs%pbv, cs%por_bar_CS)
1660 call pass_vector(cs%pbv%por_layer_widthU, cs%pbv%por_layer_widthV, &
1661 g%Domain, direction=to_all+scalar_pair, clock=id_clock_pass, halo=cs%cont_stencil)
1662 endif
1663 call set_viscous_bbl(u, v, h, tv, cs%visc, g, gv, us, cs%set_visc_CSp, cs%pbv)
1664 call cpu_clock_end(id_clock_bbl_visc)
1665 if (showcalltree) call calltree_waypoint("done with set_viscous_BBL (step_MOM_thermo)")
1666 endif
1667
1668 call cpu_clock_begin(id_clock_thermo)
1669 if (.not.cs%adiabatic) then
1670 if (cs%debug) then
1671 call uvchksum("Pre-diabatic [uv]", u, v, g%HI, haloshift=2, unscale=us%L_T_to_m_s)
1672 call hchksum(h,"Pre-diabatic h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
1673 call uvchksum("Pre-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1674 haloshift=0, unscale=gv%H_to_MKS*us%L_to_m**2)
1675 ! call MOM_state_chksum("Pre-diabatic ", u, v, h, CS%uhtr, CS%vhtr, G, GV, vel_scale=1.0)
1676 call mom_thermo_chksum("Pre-diabatic ", tv, g, us, haloshift=0)
1677 if (debug_redundant) &
1678 call check_redundant("Pre-diabatic ", u, v, g, unscale=us%L_T_to_m_s)
1679 call mom_forcing_chksum("Pre-diabatic", fluxes, g, us, haloshift=0)
1680 endif
1681
1682 call cpu_clock_begin(id_clock_diabatic)
1683
1684 call diabatic(u, v, h, tv, cs%Hml, fluxes, cs%visc, cs%ADp, cs%CDp, dtdia, &
1685 time_end_thermo, g, gv, us, cs%diabatic_CSp, cs%stoch_CS, cs%OBC, waves)
1686 fluxes%fluxes_used = .true.
1687
1688 if (cs%stoch_CS%do_skeb) then
1689 call apply_skeb(cs%G,cs%GV,cs%stoch_CS,cs%u,cs%v,cs%h,cs%tv,dtdia,time_end_thermo)
1690 endif
1691
1692 if (showcalltree) call calltree_waypoint("finished diabatic (step_MOM_thermo)")
1693
1694 if (cs%debug) then
1695 call uvchksum("Post-diabatic u", u, v, g%HI, haloshift=2, unscale=us%L_T_to_m_s)
1696 call hchksum(h, "Post-diabatic h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
1697 call uvchksum("Post-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1698 haloshift=0, unscale=gv%H_to_MKS*us%L_to_m**2)
1699 ! call MOM_state_chksum("Post-diabatic ", u, v, &
1700 ! h, CS%uhtr, CS%vhtr, G, GV, haloshift=1)
1701 if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", g%HI, haloshift=1, unscale=us%C_to_degC)
1702 if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", g%HI, haloshift=1, unscale=us%S_to_ppt)
1703 if (associated(tv%frazil)) call hchksum(tv%frazil, "Post-diabatic frazil", g%HI, haloshift=0, &
1704 unscale=us%Q_to_J_kg*us%RZ_to_kg_m2)
1705 if (associated(tv%salt_deficit)) call hchksum(tv%salt_deficit, &
1706 "Post-diabatic salt deficit", g%HI, haloshift=0, unscale=us%RZ_to_kg_m2)
1707 ! call MOM_thermo_chksum("Post-diabatic ", tv, G, US)
1708 if (debug_redundant) &
1709 call check_redundant("Post-diabatic ", u, v, g, unscale=us%L_T_to_m_s)
1710 endif
1711 call disable_averaging(cs%diag)
1712
1713 call cpu_clock_end(id_clock_diabatic)
1714 else ! complement of "if (.not.CS%adiabatic)"
1715
1716 call cpu_clock_begin(id_clock_adiabatic)
1717 call adiabatic(h, tv, fluxes, dtdia, g, gv, us, cs%diabatic_CSp)
1718 fluxes%fluxes_used = .true.
1719 call cpu_clock_end(id_clock_adiabatic)
1720
1721 if (associated(tv%T)) then
1722 dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
1723 call create_group_pass(pass_t_s, tv%T, g%Domain, to_all+omit_corners, halo=dynamics_stencil)
1724 call create_group_pass(pass_t_s, tv%S, g%Domain, to_all+omit_corners, halo=dynamics_stencil)
1725 call do_group_pass(pass_t_s, g%Domain, clock=id_clock_pass)
1726 if (cs%debug) then
1727 if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", g%HI, haloshift=1, unscale=us%C_to_degC)
1728 if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", g%HI, haloshift=1, unscale=us%S_to_ppt)
1729 endif
1730
1731 ! Update derived thermodynamic quantities.
1732 if (allocated(tv%SpV_avg)) then
1733 call calc_derived_thermo(tv, h, g, gv, us, halo=dynamics_stencil, debug=cs%debug)
1734 endif
1735 endif
1736
1737 endif ! endif for the block "if (.not.CS%adiabatic)"
1738 call cpu_clock_end(id_clock_thermo)
1739
1740 call disable_averaging(cs%diag)
1741
1742! This works in general:
1743! if (associated(tv%T)) &
1744! call totalTandS(G%HI, h, G%areaT, tv%T, tv%S, "End of step_MOM", US, GV%H_to_mks)
1745! This works only if there is no rescaling being used:
1746! if (associated(tv%T)) &
1747! call totalTandS(G%HI, h, G%areaT, tv%T, tv%S, "End of step_MOM")
1748
1749 if (showcalltree) call calltree_leave("step_MOM_thermo(), MOM.F90")
1750
1751end subroutine step_mom_thermo
1752
1753!> ALE_regridding_and_remapping does regridding (the generation of a new grid) and remapping
1754!! (from the old grid to the new grid). This is done after the themrodynamic step.
1755subroutine ale_regridding_and_remapping(CS, G, GV, US, u, v, h, tv, dtdia, Time_end_thermo)
1756 type(mom_control_struct), intent(inout) :: CS !< Master MOM control structure
1757 type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1758 type(verticalgrid_type), intent(inout) :: GV !< ocean vertical grid structure
1759 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1760 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1761 intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1762 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1763 intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
1764 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1765 intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
1766 type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
1767 real, intent(in) :: dtdia !< The time interval over which to advance [T ~> s]
1768 type(time_type), intent(in) :: Time_end_thermo !< End of averaging interval for thermo diags
1769
1770 real :: h_new(SZI_(G),SZJ_(G),SZK_(GV)) ! Layer thicknesses after regridding [H ~> m or kg m-2]
1771 real :: dzRegrid(SZI_(G),SZJ_(G),SZK_(GV)+1) ! The change in grid interface positions due to regridding,
1772 ! in the same units as thicknesses [H ~> m or kg m-2]
1773 real :: h_old_u(SZIB_(G),SZJ_(G),SZK_(GV)) ! Source grid thickness at zonal
1774 ! velocity points [H ~> m or kg m-2]
1775 real :: h_old_v(SZI_(G),SZJB_(G),SZK_(GV)) ! Source grid thickness at meridional
1776 ! velocity points [H ~> m or kg m-2]
1777 real :: h_new_u(SZIB_(G),SZJ_(G),SZK_(GV)) ! Destination grid thickness at zonal
1778 ! velocity points [H ~> m or kg m-2]
1779 real :: h_new_v(SZI_(G),SZJB_(G),SZK_(GV)) ! Destination grid thickness at meridional
1780 ! velocity points [H ~> m or kg m-2]
1781 logical :: PCM_cell(SZI_(G),SZJ_(G),SZK_(GV)) ! If true, PCM remapping should be used in a cell.
1782 logical :: use_ice_shelf ! Needed for selecting the right ALE interface.
1783 logical :: debug_redundant ! If true, check redundant values on PE boundaries when debugging.
1784 logical :: showCallTree
1785 type(group_pass_type) :: pass_T_S_h
1786 integer :: i, j, k, is, ie, js, je, nz
1787
1788 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1789 use_ice_shelf = .false.
1790 if (associated(cs%frac_shelf_h)) use_ice_shelf = .true.
1791 showcalltree = calltree_showquery()
1792 if (showcalltree) call calltree_enter("ALE_regridding_and_remapping(), MOM.F90")
1793 if (cs%debug) call query_debugging_checks(do_redundant=debug_redundant)
1794
1795 call cpu_clock_begin(id_clock_remap)
1796
1797 ! Regridding/remapping is done here, at end of thermodynamics time step
1798 ! (that may comprise several dynamical time steps)
1799 ! The routine 'ALE_regrid' can be found in 'MOM_ALE.F90'.
1800 call enable_averages(dtdia, time_end_thermo, cs%diag)
1801
1802 call cpu_clock_begin(id_clock_pass)
1803 if (associated(tv%T)) &
1804 call create_group_pass(pass_t_s_h, tv%T, g%Domain, to_all+omit_corners, halo=1)
1805 if (associated(tv%S)) &
1806 call create_group_pass(pass_t_s_h, tv%S, g%Domain, to_all+omit_corners, halo=1)
1807 call create_group_pass(pass_t_s_h, h, g%Domain, to_all+omit_corners, halo=1)
1808 call do_group_pass(pass_t_s_h, g%Domain)
1809 call cpu_clock_end(id_clock_pass)
1810
1811 call preale_tracer_diagnostics(cs%tracer_Reg, g, gv)
1812
1813 if (cs%use_particles) then
1814 call particles_to_z_space(cs%particles, h)
1815 endif
1816
1817 if (cs%debug) then
1818 call mom_state_chksum("Pre-ALE ", u, v, h, cs%uh, cs%vh, g, gv, us, omit_corners=.true.)
1819 call hchksum(tv%T,"Pre-ALE T", g%HI, haloshift=1, omit_corners=.true., unscale=us%C_to_degC)
1820 call hchksum(tv%S,"Pre-ALE S", g%HI, haloshift=1, omit_corners=.true., unscale=us%S_to_ppt)
1821 if (debug_redundant) &
1822 call check_redundant("Pre-ALE ", u, v, g, unscale=us%L_T_to_m_s)
1823 endif
1824 call cpu_clock_begin(id_clock_ale)
1825
1826 call pre_ale_diagnostics(g, gv, us, h, u, v, tv, cs%ALE_CSp)
1827 call ale_update_regrid_weights(dtdia, cs%ALE_CSp)
1828 ! Do any necessary adjustments ot the state prior to remapping.
1829 call pre_ale_adjustments(g, gv, us, h, tv, cs%tracer_Reg, cs%ALE_CSp, u, v)
1830 ! Adjust the target grids for diagnostics, in case there have been thickness adjustments.
1831 call diag_update_remap_grids(cs%diag)
1832
1833 if (use_ice_shelf) then
1834 call ale_regrid(g, gv, us, h, h_new, dzregrid, tv, cs%ALE_CSp, cs%frac_shelf_h, pcm_cell)
1835 else
1836 call ale_regrid(g, gv, us, h, h_new, dzregrid, tv, cs%ALE_CSp, pcm_cell=pcm_cell)
1837 endif
1838
1839 if (showcalltree) call calltree_waypoint("new grid generated")
1840 ! Remap all variables from the old grid h onto the new grid h_new
1841 call ale_remap_tracers(cs%ALE_CSp, g, gv, h, h_new, cs%tracer_Reg, showcalltree, dtdia, pcm_cell)
1842
1843 ! Determine the old and new grid thicknesses at velocity points.
1844 call ale_remap_set_h_vel(cs%ALE_CSp, g, gv, h, h_old_u, h_old_v, cs%OBC, debug=showcalltree)
1845 if (cs%remap_uv_using_old_alg) then
1846 call ale_remap_set_h_vel_via_dz(cs%ALE_CSp, g, gv, h_new, h_new_u, h_new_v, cs%OBC, h, dzregrid, showcalltree)
1847 else
1848 call ale_remap_set_h_vel(cs%ALE_CSp, g, gv, h_new, h_new_u, h_new_v, cs%OBC, debug=showcalltree)
1849 endif
1850
1851 ! Remap the velocity components.
1852 call ale_remap_velocities(cs%ALE_CSp, g, gv, h_old_u, h_old_v, h_new_u, h_new_v, u, v, showcalltree, &
1853 dtdia, allow_preserve_variance=.true.)
1854
1855 if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid.
1856
1857 if (cs%remap_aux_vars) then
1858 if (cs%split .and. cs%use_alt_split) then
1859 call remap_dyn_split_rk2b_aux_vars(g, gv, cs%dyn_split_RK2b_CSp, h_old_u, h_old_v, &
1860 h_new_u, h_new_v, cs%ALE_CSp)
1861 elseif (cs%split) then
1862 call remap_dyn_split_rk2_aux_vars(g, gv, cs%dyn_split_RK2_CSp, h_old_u, h_old_v, h_new_u, h_new_v, cs%ALE_CSp)
1863 endif
1864
1865 if (associated(cs%OBC) .or. associated(cs%visc%Kv_shear_Bu)) then
1866 call pass_var(h, g%Domain, complete=.false.)
1867 call pass_var(h_new, g%Domain, complete=.true.)
1868 endif
1869
1870 if (associated(cs%OBC)) &
1871 call remap_obc_fields(g, gv, h, h_new, cs%OBC, pcm_cell=pcm_cell)
1872
1873 call remap_vertvisc_aux_vars(g, gv, cs%visc, h, h_new, cs%ALE_CSp, cs%OBC)
1874 if (associated(cs%visc%Kv_shear)) &
1875 call pass_var(cs%visc%Kv_shear, g%Domain, to_all+omit_corners, clock=id_clock_pass, halo=1)
1876 endif
1877
1878 ! Replace the old grid with new one. All remapping must be done by this point in the code.
1879 !$OMP parallel do default(shared)
1880 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
1881 h(i,j,k) = h_new(i,j,k)
1882 enddo ; enddo ; enddo
1883
1884 if (showcalltree) call calltree_waypoint("finished ALE_regrid (ALE_regridding_and_remapping)")
1885 call cpu_clock_end(id_clock_ale)
1886
1887 ! Update derived thermodynamic quantities.
1888 if (allocated(cs%tv%SpV_avg)) then
1889 call calc_derived_thermo(cs%tv, cs%h, g, gv, us, halo=1, debug=cs%debug)
1890 endif
1891
1892 ! Whenever thickness changes let the diag manager know, target grids
1893 ! for vertical remapping may need to be regenerated. In non-Boussinesq mode,
1894 ! calc_derived_thermo needs to be called before diag_update_remap_grids.
1895 ! This needs to happen after the H update and before the next post_data.
1896 call diag_update_remap_grids(cs%diag)
1897
1898 call postale_tracer_diagnostics(cs%tracer_Reg, g, gv, cs%diag, dtdia)
1899
1900 if (cs%debug .and. cs%use_ALE_algorithm) then
1901 call mom_state_chksum("Post-ALE ", u, v, h, cs%uh, cs%vh, g, gv, us)
1902 call hchksum(tv%T, "Post-ALE T", g%HI, haloshift=1, unscale=us%C_to_degC)
1903 call hchksum(tv%S, "Post-ALE S", g%HI, haloshift=1, unscale=us%S_to_ppt)
1904 if (debug_redundant) &
1905 call check_redundant("Post-ALE ", u, v, g, unscale=us%L_T_to_m_s)
1906 endif
1907 if (cs%debug) then
1908 call uvchksum("Post-ALE, Post-diabatic u", u, v, g%HI, haloshift=2, unscale=us%L_T_to_m_s)
1909 call hchksum(h, "Post-ALE, Post-diabatic h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
1910 call uvchksum("Post-ALE, Post-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1911 haloshift=0, unscale=gv%H_to_MKS*us%L_to_m**2)
1912 ! call MOM_state_chksum("Post-diabatic ", u, v, &
1913 ! h, CS%uhtr, CS%vhtr, G, GV, haloshift=1)
1914 if (associated(tv%T)) call hchksum(tv%T, "Post-ALE, Post-diabatic T", g%HI, haloshift=1, unscale=us%C_to_degC)
1915 if (associated(tv%S)) call hchksum(tv%S, "Post-ALE, Post-diabatic S", g%HI, haloshift=1, unscale=us%S_to_ppt)
1916 if (associated(tv%frazil)) call hchksum(tv%frazil, "Post-ALE, Post-diabatic frazil", g%HI, haloshift=0, &
1917 unscale=us%Q_to_J_kg*us%RZ_to_kg_m2)
1918 if (associated(tv%salt_deficit)) call hchksum(tv%salt_deficit, &
1919 "Post-ALE, Post-diabatic salt deficit", g%HI, haloshift=0, unscale=us%RZ_to_kg_m2)
1920 ! call MOM_thermo_chksum("Post-diabatic ", tv, G, US)
1921 if (debug_redundant) &
1922 call check_redundant("Post-ALE, Post-diabatic ", u, v, g, unscale=us%L_T_to_m_s)
1923 endif
1924 call disable_averaging(cs%diag)
1925
1926 call cpu_clock_end(id_clock_remap)
1927
1928 if (showcalltree) call calltree_leave("ALE_regridding_and_remapping(), MOM.F90")
1929
1930end subroutine ale_regridding_and_remapping
1931
1932!> post_diabatic_halo_updates does halo updates and calculates derived thermodynamic quantities
1933!! (e.g. specific volume). This must be done after the diabatic step regardless of is ALE
1934!! cooridinates are used or not.
1935subroutine post_diabatic_halo_updates(CS, G, GV, US, u, v, h, tv)
1936 type(mom_control_struct), intent(inout) :: CS !< Master MOM control structure
1937 type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1938 type(verticalgrid_type), intent(inout) :: GV !< ocean vertical grid structure
1939 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1940 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1941 intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1942 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1943 intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
1944 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1945 intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
1946 type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
1947
1948 logical :: debug_redundant ! If true, check redundant values on PE boundaries when debugging.
1949 logical :: showCallTree
1950 type(group_pass_type) :: pass_uv_T_S_h
1951 integer :: dynamics_stencil ! The computational stencil for the calculations
1952 ! in the dynamic core.
1953
1954 showcalltree = calltree_showquery()
1955 if (showcalltree) call calltree_enter("post_diabatic_halo_updates, MOM.F90")
1956 if (cs%debug) call query_debugging_checks(do_redundant=debug_redundant)
1957
1958 if (cs%use_particles) then
1959 call particles_to_k_space(cs%particles, h)
1960 endif
1961
1962 dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
1963 call create_group_pass(pass_uv_t_s_h, u, v, g%Domain, halo=dynamics_stencil)
1964 if (associated(tv%T)) &
1965 call create_group_pass(pass_uv_t_s_h, tv%T, g%Domain, halo=dynamics_stencil)
1966 if (associated(tv%S)) &
1967 call create_group_pass(pass_uv_t_s_h, tv%S, g%Domain, halo=dynamics_stencil)
1968 call create_group_pass(pass_uv_t_s_h, h, g%Domain, halo=dynamics_stencil)
1969 call do_group_pass(pass_uv_t_s_h, g%Domain, clock=id_clock_pass)
1970
1971 if (associated(tv%frazil) .and. (.not.tv%frazil_was_reset) .and. cs%vertex_shear) &
1972 call pass_var(tv%frazil, g%Domain, halo=1)
1973
1974 ! Update derived thermodynamic quantities.
1975 if (allocated(tv%SpV_avg)) then
1976 call calc_derived_thermo(tv, h, g, gv, us, halo=dynamics_stencil, debug=cs%debug)
1977 endif
1978 if (showcalltree) call calltree_leave("post_diabatic_halo_updates, MOM.F90")
1979end subroutine post_diabatic_halo_updates
1980
1981!> step_offline is the main driver for running tracers offline in MOM6. This has been primarily
1982!! developed with ALE configurations in mind. Some work has been done in isopycnal configuration, but
1983!! the work is very preliminary. Some more detail about this capability along with some of the subroutines
1984!! called here can be found in tracers/MOM_offline_control.F90
1985subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS)
1986 type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1987 type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
1988 type(surface), intent(inout) :: sfc_state !< surface ocean state
1989 type(time_type), intent(in) :: time_start !< starting time of a segment, as a time type
1990 real, intent(in) :: time_interval !< time interval [T ~> s]
1991 type(mom_control_struct), intent(inout) :: cs !< control structure from initialize_MOM
1992
1993 ! Local pointers
1994 type(ocean_grid_type), pointer :: g => null() ! Pointer to a structure containing
1995 ! metrics and related information
1996 type(verticalgrid_type), pointer :: gv => null() ! Pointer to structure containing information
1997 ! about the vertical grid
1998 type(unit_scale_type), pointer :: us => null() ! Pointer to a structure containing
1999 ! various unit conversion factors
2000
2001 logical :: first_iter !< True if this is the first time step_offline has been called in a given interval
2002 logical :: last_iter !< True if this is the last time step_tracer is to be called in an offline interval
2003 logical :: do_vertical !< If enough time has elapsed, do the diabatic tracer sources/sinks
2004 logical :: adv_converged !< True if all the horizontal fluxes have been used
2005
2006 real, allocatable, dimension(:,:,:) :: h_new ! Layer thicknesses after regridding [H ~> m or kg m-2]
2007 real, allocatable, dimension(:,:,:) :: dzregrid ! The change in grid interface positions due to regridding,
2008 ! in the same units as thicknesses [H ~> m or kg m-2]
2009 real :: dt_offline ! The offline timestep for advection [T ~> s]
2010 real :: dt_offline_vertical ! The offline timestep for vertical fluxes and remapping [T ~> s]
2011 logical :: skip_diffusion
2012
2013 type(time_type), pointer :: accumulated_time => null()
2014 type(time_type), pointer :: vertical_time => null()
2015 integer :: dynamics_stencil ! The computational stencil for the calculations
2016 ! in the dynamic core.
2017 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
2018
2019 ! 3D pointers
2020 real, dimension(:,:,:), pointer :: &
2021 uhtr => null(), & ! Accumulated zonal thickness fluxes to advect tracers [H L2 ~> m3 or kg]
2022 vhtr => null(), & ! Accumulated meridional thickness fluxes to advect tracers [H L2 ~> m3 or kg]
2023 eatr => null(), & ! Layer entrainment rates across the interface above [H ~> m or kg m-2]
2024 ebtr => null(), & ! Layer entrainment rates across the interface below [H ~> m or kg m-2]
2025 h_end => null() ! Layer thicknesses at the end of a step [H ~> m or kg m-2]
2026
2027 type(time_type) :: time_end ! End time of a segment, as a time type
2028
2029 ! Grid-related pointer assignments
2030 g => cs%G ; gv => cs%GV ; us => cs%US
2031
2032 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
2033 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2034
2035 call cpu_clock_begin(id_clock_offline_tracer)
2036 call extract_offline_main(cs%offline_CSp, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, &
2037 vertical_time, dt_offline, dt_offline_vertical, skip_diffusion)
2038 time_end = increment_date(time_start, seconds=floor(us%T_to_s*time_interval+0.001))
2039
2040 call enable_averages(time_interval, time_end, cs%diag)
2041
2042 ! Check to see if this is the first iteration of the offline interval
2043 first_iter = (accumulated_time == real_to_time(0.0))
2044
2045 ! Check to see if vertical tracer functions should be done
2046 do_vertical = (first_iter .or. (accumulated_time >= vertical_time))
2047 if (do_vertical) vertical_time = accumulated_time + real_to_time(dt_offline_vertical, unscale=us%T_to_s)
2048
2049 ! Increment the amount of time elapsed since last read and check if it's time to roll around
2050 accumulated_time = accumulated_time + real_to_time(time_interval, unscale=us%T_to_s)
2051
2052 last_iter = (accumulated_time >= real_to_time(dt_offline, unscale=us%T_to_s))
2053
2054 if (cs%use_ALE_algorithm) then
2055 ! If this is the first iteration in the offline timestep, then we need to read in fields and
2056 ! perform the main advection.
2057 if (first_iter) then
2058 call mom_mesg("Reading in new offline fields")
2059 ! Read in new transport and other fields
2060 ! call update_transport_from_files(G, GV, CS%offline_CSp, h_end, eatr, ebtr, uhtr, vhtr, &
2061 ! CS%tv%T, CS%tv%S, fluxes, CS%use_ALE_algorithm)
2062 ! call update_transport_from_arrays(CS%offline_CSp)
2063 call update_offline_fields(cs%offline_CSp, g, gv, us, cs%h, fluxes, cs%use_ALE_algorithm)
2064
2065 ! Apply any fluxes into the ocean
2066 call offline_fw_fluxes_into_ocean(g, gv, cs%offline_CSp, fluxes, cs%h)
2067
2068 if (.not.cs%diabatic_first) then
2069 call offline_advection_ale(fluxes, time_start, time_interval, g, gv, us, cs%offline_CSp, &
2070 id_clock_ale, cs%h, uhtr, vhtr, converged=adv_converged)
2071
2072 ! Redistribute any remaining transport
2073 call offline_redistribute_residual(cs%offline_CSp, g, gv, us, cs%h, uhtr, vhtr, adv_converged)
2074
2075 ! Perform offline diffusion if requested
2076 if (.not. skip_diffusion) then
2077 if (cs%VarMix%use_variable_mixing) then
2078 call pass_var(cs%h, g%Domain)
2079 call calc_resoln_function(cs%h, cs%tv, g, gv, us, cs%VarMix, cs%MEKE, cs%OBC, dt_offline)
2080 call calc_depth_function(g, cs%VarMix)
2081 call calc_slope_functions(cs%h, cs%tv, dt_offline, g, gv, us, cs%VarMix, obc=cs%OBC)
2082 endif
2083 call tracer_hordiff(cs%h, dt_offline, cs%MEKE, cs%VarMix, cs%visc, g, gv, us, &
2084 cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
2085 endif
2086 endif
2087 endif
2088 ! The functions related to column physics of tracers is performed separately in ALE mode
2089 if (do_vertical) then
2090 call offline_diabatic_ale(fluxes, time_start, time_end, g, gv, us, cs%offline_CSp, &
2091 cs%h, cs%tv, eatr, ebtr)
2092 endif
2093
2094 ! Last thing that needs to be done is the final ALE remapping
2095 if (last_iter) then
2096 if (cs%diabatic_first) then
2097 call offline_advection_ale(fluxes, time_start, time_interval, g, gv, us, cs%offline_CSp, &
2098 id_clock_ale, cs%h, uhtr, vhtr, converged=adv_converged)
2099
2100 ! Redistribute any remaining transport and perform the remaining advection
2101 call offline_redistribute_residual(cs%offline_CSp, g, gv, us, cs%h, uhtr, vhtr, adv_converged)
2102 ! Perform offline diffusion if requested
2103 if (.not. skip_diffusion) then
2104 if (cs%VarMix%use_variable_mixing) then
2105 call pass_var(cs%h, g%Domain)
2106 call calc_resoln_function(cs%h, cs%tv, g, gv, us, cs%VarMix, cs%MEKE, cs%OBC, dt_offline)
2107 call calc_depth_function(g, cs%VarMix)
2108 call calc_slope_functions(cs%h, cs%tv, dt_offline, g, gv, us, cs%VarMix, obc=cs%OBC)
2109 endif
2110 call tracer_hordiff(cs%h, dt_offline, cs%MEKE, cs%VarMix, cs%visc, g, gv, us, &
2111 cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
2112 endif
2113 endif
2114
2115 call mom_mesg("Last iteration of offline interval")
2116
2117 ! Apply freshwater fluxes out of the ocean
2118 call offline_fw_fluxes_out_ocean(g, gv, cs%offline_CSp, fluxes, cs%h)
2119 ! These diagnostic can be used to identify which grid points did not converge within
2120 ! the specified number of advection sub iterations
2121 call post_offline_convergence_diags(g, gv, cs%offline_CSp, cs%h, h_end, uhtr, vhtr)
2122
2123 ! Call ALE one last time to make sure that tracers are remapped onto the layer thicknesses
2124 ! stored from the forward run
2125 call cpu_clock_begin(id_clock_ale)
2126
2127 ! Do any necessary adjustments ot the state prior to remapping.
2128 call pre_ale_adjustments(g, gv, us, h_end, cs%tv, cs%tracer_Reg, cs%ALE_CSp)
2129
2130 allocate(h_new(isd:ied, jsd:jed, nz), source=0.0)
2131 allocate(dzregrid(isd:ied, jsd:jed, nz+1), source=0.0)
2132
2133 ! Generate the new grid based on the tracer grid at the end of the interval.
2134 call ale_regrid(g, gv, us, h_end, h_new, dzregrid, cs%tv, cs%ALE_CSp)
2135
2136 ! Remap the tracers from the previous tracer grid onto the new grid. The thicknesses that
2137 ! are used are intended to ensure that in the case where transports don't quite conserve,
2138 ! the offline layer thicknesses do not drift too far away from the online model.
2139 call ale_remap_tracers(cs%ALE_CSp, g, gv, cs%h, h_new, cs%tracer_Reg, debug=cs%debug)
2140 if (allocated(cs%tv%SpV_avg)) cs%tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid.
2141
2142 ! Update the tracer grid.
2143 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
2144 cs%h(i,j,k) = h_new(i,j,k)
2145 enddo ; enddo ; enddo
2146
2147 deallocate(h_new, dzregrid)
2148
2149 call cpu_clock_end(id_clock_ale)
2150 call pass_var(cs%h, g%Domain)
2151 endif
2152 else ! NON-ALE MODE...NOT WELL TESTED
2153 call mom_error(warning, &
2154 "Offline tracer mode in non-ALE configuration has not been thoroughly tested")
2155 ! Note that for the layer mode case, the calls to tracer sources and sinks is embedded in
2156 ! main_offline_advection_layer. Warning: this may not be appropriate for tracers that
2157 ! exchange with the atmosphere
2158 if (abs(time_interval - dt_offline) > 1.0e-6*us%s_to_T) then
2159 call mom_error(fatal, &
2160 "For offline tracer mode in a non-ALE configuration, dt_offline must equal time_interval")
2161 endif
2162 call update_offline_fields(cs%offline_CSp, g, gv, us, cs%h, fluxes, cs%use_ALE_algorithm)
2163 call offline_advection_layer(fluxes, time_start, time_interval, g, gv, us, cs%offline_CSp, &
2164 cs%h, eatr, ebtr, uhtr, vhtr)
2165 ! Perform offline diffusion if requested
2166 if (.not. skip_diffusion) then
2167 call tracer_hordiff(h_end, dt_offline, cs%MEKE, cs%VarMix, cs%visc, g, gv, us, &
2168 cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
2169 endif
2170
2171 cs%h = h_end
2172
2173 call pass_var(cs%tv%T, g%Domain)
2174 call pass_var(cs%tv%S, g%Domain)
2175 call pass_var(cs%h, g%Domain)
2176
2177 endif
2178
2179 call adjust_ssh_for_p_atm(cs%tv, g, gv, us, cs%ave_ssh_ibc, forces%p_surf_SSH, &
2180 cs%calc_rho_for_sea_lev)
2181 call extract_surface_state(cs, sfc_state)
2182
2183 call disable_averaging(cs%diag)
2184 call pass_var(cs%tv%T, g%Domain)
2185 call pass_var(cs%tv%S, g%Domain)
2186 call pass_var(cs%h, g%Domain)
2187
2188 fluxes%fluxes_used = .true.
2189
2190 ! Update derived thermodynamic quantities.
2191 if (allocated(cs%tv%SpV_avg)) then
2192 dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
2193 call calc_derived_thermo(cs%tv, cs%h, g, gv, us, halo=dynamics_stencil)
2194 endif
2195
2196 if (last_iter) then
2197 accumulated_time = real_to_time(0.0)
2198 endif
2199
2200 call cpu_clock_end(id_clock_offline_tracer)
2201
2202end subroutine step_offline
2203
2204!> Initialize MOM, including memory allocation, setting up parameters and diagnostics,
2205!! initializing the ocean state variables, and initializing subsidiary modules
2206subroutine initialize_mom(Time, Time_init, param_file, dirs, CS, &
2207 Time_in, offline_tracer_mode, input_restart_file, diag_ptr, &
2208 count_calls, tracer_flow_CSp, ice_shelf_CSp, waves_CSp, ensemble_num, &
2209 calve_ice_shelf_bergs)
2210 type(time_type), target, intent(inout) :: time !< model time, set in this routine
2211 type(time_type), intent(in) :: time_init !< The start time for the coupled model's calendar
2212 type(param_file_type), intent(out) :: param_file !< structure indicating parameter file to parse
2213 type(directories), intent(out) :: dirs !< structure with directory paths
2214 type(mom_control_struct), intent(inout), target :: cs !< pointer set in this routine to MOM control structure
2215 type(time_type), optional, intent(in) :: time_in !< time passed to MOM_initialize_state when
2216 !! model is not being started from a restart file
2217 logical, optional, intent(out) :: offline_tracer_mode !< True is returned if tracers are being run offline
2218 character(len=*),optional, intent(in) :: input_restart_file !< If present, name of restart file to read
2219 type(diag_ctrl), optional, pointer :: diag_ptr !< A pointer set in this routine to the diagnostic
2220 !! regulatory structure
2221 type(tracer_flow_control_cs), &
2222 optional, pointer :: tracer_flow_csp !< A pointer set in this routine to
2223 !! the tracer flow control structure.
2224 logical, optional, intent(in) :: count_calls !< If true, nstep_tot counts the number of
2225 !! calls to step_MOM instead of the number of
2226 !! dynamics timesteps.
2227 type(ice_shelf_cs), optional, pointer :: ice_shelf_csp !< A pointer to an ice shelf control structure
2228 type(wave_parameters_cs), &
2229 optional, pointer :: waves_csp !< An optional pointer to a wave property CS
2230 integer, optional :: ensemble_num !< Ensemble index provided by the cap (instead of FMS
2231 !! ensemble manager)
2232 logical, optional :: calve_ice_shelf_bergs !< If true, will add point iceberg calving variables to the ice
2233 !! shelf restart
2234 ! local variables
2235 type(ocean_grid_type), pointer :: g => null() ! A pointer to the metric grid use for the run
2236 type(ocean_grid_type), pointer :: g_in => null() ! Pointer to the input grid
2237 type(hor_index_type), pointer :: hi => null() ! A hor_index_type for array extents
2238 type(hor_index_type), target :: hi_in ! HI on the input grid
2239 type(hor_index_type) :: hi_in_unmasked ! HI on the unmasked input grid
2240 type(verticalgrid_type), pointer :: gv => null()
2241 type(dyn_horgrid_type), pointer :: dg => null(), test_dg => null()
2242 type(dyn_horgrid_type), pointer :: dg_in => null()
2243 type(dyn_horgrid_type), pointer :: dg_unmasked_in => null()
2244 type(diag_ctrl), pointer :: diag => null()
2245 type(unit_scale_type), pointer :: us => null()
2246 type(mom_restart_cs), pointer :: restart_csp => null()
2247 character(len=4), parameter :: vers_num = 'v2.0'
2248 integer :: turns ! Number of grid quarter-turns
2249 logical :: point_calving
2250
2251 ! Initial state on the input index map
2252 real, allocatable :: u_in(:,:,:) ! Initial zonal velocities [L T-1 ~> m s-1]
2253 real, allocatable :: v_in(:,:,:) ! Initial meridional velocities [L T-1 ~> m s-1]
2254 real, allocatable :: h_in(:,:,:) ! Initial layer thicknesses [H ~> m or kg m-2]
2255 real, allocatable, target :: frac_shelf_in(:,:) ! Initial fraction of the total cell area occupied
2256 ! by an ice shelf [nondim]
2257 real, allocatable, target :: mass_shelf_in(:,:) ! Initial mass of ice shelf contained within a grid cell
2258 ! [R Z ~> kg m-2]
2259 real, allocatable, target :: t_in(:,:,:) ! Initial temperatures [C ~> degC]
2260 real, allocatable, target :: s_in(:,:,:) ! Initial salinities [S ~> ppt]
2261
2262 type(ocean_obc_type), pointer :: obc_in => null()
2263 type(sponge_cs), pointer :: sponge_in_csp => null()
2264 type(ale_sponge_cs), pointer :: ale_sponge_in_csp => null()
2265 type(oda_incupd_cs),pointer :: oda_incupd_in_csp => null()
2266 ! This include declares and sets the variable "version".
2267# include "version_variable.h"
2268
2269 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
2270 integer :: isdb, iedb, jsdb, jedb
2271 real :: dtbt ! If negative, this specifies the barotropic timestep as a fraction
2272 ! of the maximum stable value [nondim].
2273
2274 real, allocatable, dimension(:,:) :: eta ! free surface height or column mass [H ~> m or kg m-2]
2275 real, allocatable, dimension(:,:,:) :: h_new ! Layer thicknesses after regridding [H ~> m or kg m-2]
2276 real, allocatable, dimension(:,:,:) :: dzregrid ! The change in grid interface positions due to regridding,
2277 ! in the same units as thicknesses [H ~> m or kg m-2]
2278 real, allocatable, dimension(:,:,:) :: h_old_u ! Source grid thickness at zonal velocity points [H ~> m or kg m-2]
2279 real, allocatable, dimension(:,:,:) :: h_old_v ! Source grid thickness at meridional velocity
2280 ! points [H ~> m or kg m-2]
2281 real, allocatable, dimension(:,:,:) :: h_new_u ! Destination grid thickness at zonal
2282 ! velocity points [H ~> m or kg m-2]
2283 real, allocatable, dimension(:,:,:) :: h_new_v ! Destination grid thickness at meridional
2284 ! velocity points [H ~> m or kg m-2]
2285 logical, allocatable, dimension(:,:,:) :: pcm_cell ! If true, PCM remapping should be used in a cell.
2286 type(group_pass_type) :: tmp_pass_uv_t_s_h, pass_uv_t_s_h
2287
2288 real :: hmix_z, hmix_uv_z ! Temporary variables with averaging depths [Z ~> m]
2289 real :: hfrz_z ! Temporary variable with the melt potential depth [Z ~> m]
2290 real :: default_val ! The default value for DTBT_RESET_PERIOD [s]
2291 logical :: write_geom_files ! If true, write out the grid geometry files.
2292 logical :: new_sim ! If true, this has been determined to be a new simulation
2293 logical :: use_geothermal ! If true, apply geothermal heating.
2294 logical :: use_eos ! If true, density calculated from T & S using an equation of state.
2295 logical :: symmetric ! If true, use symmetric memory allocation.
2296 logical :: save_ic ! If true, save the initial conditions.
2297 logical :: do_unit_tests ! If true, call unit tests.
2298 logical :: fpmix ! Needed to decide if BLD should be passed to RK2.
2299 logical :: test_grid_copy = .false.
2300
2301 logical :: bulkmixedlayer ! If true, a refined bulk mixed layer scheme is used
2302 ! with nkml sublayers and nkbl buffer layer.
2303 logical :: use_temperature ! If true, temperature and salinity used as state variables.
2304 logical :: use_p_surf_in_eos ! If true, always include the surface pressure contributions
2305 ! in equation of state calculations.
2306 logical :: use_frazil ! If true, liquid seawater freezes if temp below freezing,
2307 ! with accumulated heat deficit returned to surface ocean.
2308 logical :: bound_salinity ! If true, salt is added to keep salinity above
2309 ! a minimum value, and the deficit is reported.
2310 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
2311 logical :: use_cont_abss ! If true, the prognostics T & S are conservative temperature
2312 ! and absolute salinity. Care should be taken to convert them
2313 ! to potential temperature and practical salinity before
2314 ! exchanging them with the coupler and/or reporting T&S diagnostics.
2315 logical :: advect_ts ! If false, then no horizontal advection of temperature
2316 ! and salnity is performed
2317 logical :: use_ice_shelf ! Needed for ALE
2318 logical :: global_indexing ! If true use global horizontal index values instead
2319 ! of having the data domain on each processor start at 1.
2320 logical :: bathy_at_vel ! If true, also define bathymetric fields at the
2321 ! the velocity points.
2322 logical :: calc_dtbt ! Indicates whether the dynamically adjusted barotropic
2323 ! time step needs to be updated before it is used.
2324 logical :: debug_truncations ! If true, turn on diagnostics useful for debugging truncations.
2325 integer :: first_direction ! An integer that indicates which direction is to be
2326 ! updated first in directionally split parts of the
2327 ! calculation.
2328 logical :: enable_bugs ! If true, the defaults for certain recently added bug-fix flags are
2329 ! set to recreate the bugs so that the code can be moved forward
2330 ! without changing answers for existing configurations. When this is
2331 ! false, bugs are only used if they are actively selected.
2332 logical :: non_bous ! If true, this run is fully non-Boussinesq
2333 logical :: boussinesq ! If true, this run is fully Boussinesq
2334 logical :: semi_boussinesq ! If true, this run is partially non-Boussinesq
2335 logical :: use_kpp ! If true, diabatic is using KPP vertical mixing
2336 logical :: mle_use_pbl_mld ! If true, use stored boundary layer depths for submesoscale restratification.
2337 logical :: obc_reservoir_init_bug
2338 logical :: obc_bgc_time_ref_bug ! If true, use the start of the current run (not the overall
2339 ! start time) as the reference for OBC BGC tracer update schedule.
2340 integer :: nkml, nkbl, verbosity, write_geom, number_of_obc_segments
2341 integer :: dynamics_stencil ! The computational stencil for the calculations
2342 ! in the dynamic core.
2343 real :: salin_underflow ! A tiny value of salinity below which the it is set to 0 [S ~> ppt]
2344 real :: temp_underflow ! A tiny magnitude of temperatures below which they are set to 0 [C ~> degC]
2345 real :: conv2watt ! A conversion factor from temperature fluxes to heat
2346 ! fluxes [J m-2 H-1 C-1 ~> J m-3 degC-1 or J kg-1 degC-1]
2347 real :: conv2salt ! A conversion factor for salt fluxes [m H-1 ~> 1] or [kg m-2 H-1 ~> 1]
2348 character(len=48) :: s_flux_units
2349
2350 type(vardesc) :: vd_t, vd_s ! Structures describing temperature and salinity variables.
2351 type(time_type) :: start_time
2352 type(ocean_internal_state) :: mom_internal_state
2353 type(mom_domain_type), pointer :: mom_dom_unmasked => null() ! Unmasked MOM domain instance
2354 ! (To be used for writing out ocean geometry)
2355 character(len=240) :: geom_file ! Name of the ocean geometry file
2356
2357 cs%Time => time
2358
2359 id_clock_ocean = cpu_clock_id('Ocean', grain=clock_component)
2360 id_clock_init = cpu_clock_id('Ocean Initialization', grain=clock_subcomponent)
2361 call cpu_clock_begin(id_clock_ocean) ; call cpu_clock_begin(id_clock_init)
2362
2363 start_time = time ; if (present(time_in)) start_time = time_in
2364
2365 ! Read paths and filenames from namelist and store in "dirs".
2366 ! Also open the parsed input parameter file(s) and setup param_file.
2367 call get_mom_input(param_file, dirs, default_input_filename=input_restart_file, ensemble_num=ensemble_num)
2368
2369 verbosity = 2 ; call read_param(param_file, "VERBOSITY", verbosity)
2370 call mom_set_verbosity(verbosity, .true.)
2371 call calltree_enter("initialize_MOM(), MOM.F90")
2372
2373 call find_obsolete_params(param_file)
2374
2375 ! Determining the internal unit scaling factors for this run.
2376 call unit_scaling_init(param_file, cs%US)
2377 us => cs%US
2378
2379 ! Read relevant parameters and write them to the model log.
2380 call log_version(param_file, "MOM", version, "", log_to_all=.true., layout=.true., debugging=.true.)
2381 call get_param(param_file, "MOM", "VERBOSITY", verbosity, &
2382 "Integer controlling level of messaging\n" // &
2383 "\t0 = Only FATAL messages\n" // &
2384 "\t2 = Only FATAL, WARNING, NOTE [default]\n" // &
2385 "\t6 = Above plus call tree messages\n" //&
2386 "\t9 = All)", default=2, debuggingparam=.true.)
2387 call get_param(param_file, "MOM", "DO_UNIT_TESTS", do_unit_tests, &
2388 "If True, exercises unit tests at model start up.", &
2389 default=.false., debuggingparam=.true.)
2390 if (do_unit_tests) then
2391 id_clock_unit_tests = cpu_clock_id('(Ocean unit tests)', grain=clock_module)
2392 call cpu_clock_begin(id_clock_unit_tests)
2393 call unit_tests(verbosity)
2394 call cpu_clock_end(id_clock_unit_tests)
2395 endif
2396
2397 call get_param(param_file, "MOM", "SPLIT", cs%split, &
2398 "Use the split time stepping if true.", default=.true.)
2399 call get_param(param_file, "MOM", "SPLIT_RK2B", cs%use_alt_split, &
2400 "If true, use a version of the split explicit time stepping scheme that "//&
2401 "exchanges velocities with step_MOM that have the average barotropic phase over "//&
2402 "a baroclinic timestep rather than the instantaneous barotropic phase.", &
2403 default=.false., do_not_log=.not.cs%split)
2404 if (cs%split) then
2405 cs%use_RK2 = .false.
2406 else
2407 call get_param(param_file, "MOM", "USE_RK2", cs%use_RK2, &
2408 "If true, use RK2 instead of RK3 in the unsplit time stepping.", &
2409 default=.false.)
2410 endif
2411
2412 ! FPMIX is needed to decide if boundary layer depth should be passed to RK2
2413 call get_param(param_file, '', "FPMIX", fpmix, &
2414 "If true, add non-local momentum flux increments and diffuse down the Eulerian gradient.", &
2415 default=.false., do_not_log=.true.)
2416
2417 if (fpmix .and. .not. cs%split) then
2418 call mom_error(fatal, "initialize_MOM: "//&
2419 "FPMIX=True only works when SPLIT=True.")
2420 endif
2421
2422 call get_param(param_file, "MOM", "BOUSSINESQ", boussinesq, &
2423 "If true, make the Boussinesq approximation.", default=.true., do_not_log=.true.)
2424 call get_param(param_file, "MOM", "SEMI_BOUSSINESQ", semi_boussinesq, &
2425 "If true, do non-Boussinesq pressure force calculations and use mass-based "//&
2426 "thicknesses, but use RHO_0 to convert layer thicknesses into certain "//&
2427 "height changes. This only applies if BOUSSINESQ is false.", &
2428 default=.true., do_not_log=.true.)
2429 non_bous = .not.(boussinesq .or. semi_boussinesq)
2430 call get_param(param_file, "MOM", "CALC_RHO_FOR_SEA_LEVEL", cs%calc_rho_for_sea_lev, &
2431 "If true, the in-situ density is used to calculate the "//&
2432 "effective sea level that is returned to the coupler. If false, "//&
2433 "the Boussinesq parameter RHO_0 is used.", default=non_bous)
2434 call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, &
2435 "If true, Temperature and salinity are used as state "//&
2436 "variables.", default=.true.)
2437 call get_param(param_file, "MOM", "USE_EOS", use_eos, &
2438 "If true, density is calculated from temperature and "//&
2439 "salinity with an equation of state. If USE_EOS is "//&
2440 "true, ENABLE_THERMODYNAMICS must be true as well.", &
2441 default=use_temperature)
2442 call get_param(param_file, "MOM", "DIABATIC_FIRST", cs%diabatic_first, &
2443 "If true, apply diabatic and thermodynamic processes, "//&
2444 "including buoyancy forcing and mass gain or loss, "//&
2445 "before stepping the dynamics forward.", default=.false.)
2446 call get_param(param_file, "MOM", "USE_CONTEMP_ABSSAL", use_cont_abss, &
2447 "If true, the prognostics T&S are the conservative temperature "//&
2448 "and absolute salinity. Care should be taken to convert them "//&
2449 "to potential temperature and practical salinity before "//&
2450 "exchanging them with the coupler and/or reporting T&S diagnostics.", &
2451 default=.false.)
2452 cs%tv%T_is_conT = use_cont_abss ; cs%tv%S_is_absS = use_cont_abss
2453 call get_param(param_file, "MOM", "ADIABATIC", cs%adiabatic, &
2454 "There are no diapycnal mass fluxes if ADIABATIC is true. "//&
2455 "This assumes that KD = 0.0 and that there is no buoyancy forcing, "//&
2456 "but makes the model faster by eliminating subroutine calls.", default=.false.)
2457 call get_param(param_file, "MOM", "DO_DYNAMICS", cs%do_dynamics, &
2458 "If False, skips the dynamics calls that update u & v, as well as "//&
2459 "the gravity wave adjustment to h. This may be a fragile feature, "//&
2460 "but can be useful during development", default=.true.)
2461 call get_param(param_file, "MOM", "ADVECT_TS", advect_ts, &
2462 "If True, advect temperature and salinity horizontally "//&
2463 "If False, T/S are registered for advection. "//&
2464 "This is intended only to be used in offline tracer mode "//&
2465 "and is by default false in that case.", &
2466 do_not_log=.true., default=.true.)
2467 if (present(offline_tracer_mode)) then ! Only read this parameter in enabled modes
2468 call get_param(param_file, "MOM", "OFFLINE_TRACER_MODE", cs%offline_tracer_mode, &
2469 "If true, barotropic and baroclinic dynamics, thermodynamics "//&
2470 "are all bypassed with all the fields necessary to integrate "//&
2471 "the tracer advection and diffusion equation are read in from "//&
2472 "files stored from a previous integration of the prognostic model. "//&
2473 "NOTE: This option only used in the ocean_solo_driver.", default=.false.)
2474 if (cs%offline_tracer_mode) then
2475 call get_param(param_file, "MOM", "ADVECT_TS", advect_ts, &
2476 "If True, advect temperature and salinity horizontally "//&
2477 "If False, T/S are registered for advection. "//&
2478 "This is intended only to be used in offline tracer mode, "//&
2479 "and is by default false in that case", &
2480 default=.false. )
2481 endif
2482 endif
2483 call get_param(param_file, "MOM", "USE_REGRIDDING", cs%use_ALE_algorithm, &
2484 "If True, use the ALE algorithm (regridding/remapping). "//&
2485 "If False, use the layered isopycnal algorithm.", default=.false. )
2486 call get_param(param_file, "MOM", "REMAP_UV_USING_OLD_ALG", cs%remap_uv_using_old_alg, &
2487 "If true, uses the old remapping-via-a-delta-z method for "//&
2488 "remapping u and v. If false, uses the new method that remaps "//&
2489 "between grids described by an old and new thickness.", &
2490 default=.false., do_not_log=.not.cs%use_ALE_algorithm)
2491 call get_param(param_file, "MOM", "REMAP_AUXILIARY_VARS", cs%remap_aux_vars, &
2492 "If true, apply ALE remapping to all of the auxiliary 3-dimensional "//&
2493 "variables that are needed to reproduce across restarts, similarly to "//&
2494 "what is already being done with the primary state variables. "//&
2495 "The default should be changed to true.", default=.false., &
2496 do_not_log=.not.cs%use_ALE_algorithm)
2497 call get_param(param_file, "MOM", "BULKMIXEDLAYER", bulkmixedlayer, &
2498 "If true, use a Kraus-Turner-like bulk mixed layer "//&
2499 "with transitional buffer layers. Layers 1 through "//&
2500 "NKML+NKBL have variable densities. There must be at "//&
2501 "least NKML+NKBL+1 layers if BULKMIXEDLAYER is true. "//&
2502 "BULKMIXEDLAYER can not be used with USE_REGRIDDING. "//&
2503 "The default is influenced by ENABLE_THERMODYNAMICS.", &
2504 default=use_temperature .and. .not.cs%use_ALE_algorithm)
2505 call get_param(param_file, "MOM", "USE_POROUS_BARRIER", cs%use_porbar, &
2506 "If true, use porous barrier to constrain the widths "//&
2507 "and face areas at the edges of the grid cells. ", &
2508 default=.false.)
2509 call get_param(param_file, "MOM", "BATHYMETRY_AT_VEL", bathy_at_vel, &
2510 "If true, there are separate values for the basin depths "//&
2511 "at velocity points. Otherwise the effects of topography "//&
2512 "are entirely determined from thickness points.", &
2513 default=.false.)
2514 call get_param(param_file, "MOM", "USE_WAVES", cs%UseWaves, default=.false., &
2515 do_not_log=.true.)
2516
2517 call get_param(param_file, "MOM", "DEBUG", cs%debug, &
2518 "If true, write out verbose debugging data.", &
2519 default=.false., debuggingparam=.true.)
2520 call get_param(param_file, "MOM", "DEBUG_TRUNCATIONS", debug_truncations, &
2521 "If true, calculate all diagnostics that are useful for "//&
2522 "debugging truncations.", default=.false., debuggingparam=.true.)
2523 call get_param(param_file, "MOM", "OBC_NUMBER_OF_SEGMENTS", number_of_obc_segments, &
2524 default=0, do_not_log=.true.)
2525 call get_param(param_file, "MOM", "DEBUG_OBCS", cs%debug_OBCs, &
2526 "If true, write out verbose debugging data about OBCs.", &
2527 default=.false., debuggingparam=.true., do_not_log=(number_of_obc_segments<=0))
2528 call get_param(param_file, "MOM", "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
2529 "If true, the defaults for certain recently added bug-fix flags are set to "//&
2530 "recreate the bugs so that the code can be moved forward without changing "//&
2531 "answers for existing configurations. The defaults for groups of bug-fix "//&
2532 "flags are periodically changed to correct the bugs, at which point this "//&
2533 "parameter will no longer be used to set their default. Setting this to false "//&
2534 "means that bugs are only used if they are actively selected, but it also "//&
2535 "means that answers may change when code is updated due to newly found bugs.", &
2536 default=.true.)
2537
2538 call get_param(param_file, "MOM", "DT", cs%dt, &
2539 "The (baroclinic) dynamics time step. The time-step that "//&
2540 "is actually used will be an integer fraction of the "//&
2541 "forcing time-step (DT_FORCING in ocean-only mode or the "//&
2542 "coupling timestep in coupled mode.)", units="s", scale=us%s_to_T, &
2543 fail_if_missing=.true.)
2544 call get_param(param_file, "MOM", "DT_THERM", cs%dt_therm, &
2545 "The thermodynamic time step. Ideally DT_THERM should be an "//&
2546 "integer multiple of DT and of DT_TRACER_ADVECT "//&
2547 "and less than the forcing or coupling time-step. However, if "//&
2548 "THERMO_SPANS_COUPLING is true, DT_THERM can be an integer multiple "//&
2549 "of the coupling timestep. By default DT_THERM is set to DT.", &
2550 units="s", scale=us%s_to_T, default=us%T_to_s*cs%dt)
2551 call get_param(param_file, "MOM", "THERMO_SPANS_COUPLING", cs%thermo_spans_coupling, &
2552 "If true, the MOM will take thermodynamic "//&
2553 "timesteps that can be longer than the coupling timestep. "//&
2554 "The actual thermodynamic timestep that is used in this "//&
2555 "case is the largest integer multiple of the coupling "//&
2556 "timestep that is less than or equal to DT_THERM.", default=.false.)
2557 call get_param(param_file, "MOM", "DT_TRACER_ADVECT", cs%dt_tr_adv, &
2558 "The tracer advection time step. Ideally DT_TRACER_ADVECT should be an "//&
2559 "integer multiple of DT, less than DT_THERM, and less than the forcing "//&
2560 "or coupling time-step. However, if TRADV_SPANS_COUPLING is true, "//&
2561 "DT_TRACER_ADVECT can be longer than the coupling timestep. By "//&
2562 "default DT_TRACER_ADVECT is set to DT_THERM.", &
2563 units="s", scale=us%s_to_T, default=us%T_to_s*cs%dt_therm)
2564 call get_param(param_file, "MOM", "TRADV_SPANS_COUPLING", cs%tradv_spans_coupling, &
2565 "If true, the MOM will take tracer advection "//&
2566 "timesteps that can be longer than the coupling timestep. "//&
2567 "The actual tracer advection timestep that is used in this "//&
2568 "case is the largest integer multiple of the coupling "//&
2569 "timestep that is less than or equal to DT_TRACER_ADVECT.", &
2570 default=cs%thermo_spans_coupling)
2571 if ( cs%diabatic_first .and. (cs%dt_tr_adv /= cs%dt_therm) ) then
2572 call mom_error(fatal,"MOM: If using DIABATIC_FIRST, DT_TRACER_ADVECT must equal DT_THERM.")
2573 endif
2574 call get_param(param_file, "MOM", "THICKNESSDIFFUSE", cs%thickness_diffuse, &
2575 "If true, isopycnal surfaces are diffused with a Laplacian "//&
2576 "coefficient of KHTH.", default=.false.)
2577 call get_param(param_file, "MOM", "APPLY_INTERFACE_FILTER", cs%interface_filter, &
2578 "If true, model interface heights are subjected to a grid-scale "//&
2579 "dependent spatial smoothing, often with biharmonic filter.", default=.false.)
2580 call get_param(param_file, "MOM", "THICKNESSDIFFUSE_FIRST", cs%thickness_diffuse_first, &
2581 "If true, do thickness diffusion or interface height smoothing before dynamics. "//&
2582 "This is only used if THICKNESSDIFFUSE or APPLY_INTERFACE_FILTER is true.", &
2583 default=.false., do_not_log=.not.(cs%thickness_diffuse.or.cs%interface_filter))
2584 cs%interface_filter_dt_bug = .false.
2585 if ((.not.cs%thickness_diffuse_first .and. cs%interface_filter) .or. &
2586 (cs%thickness_diffuse_first .and. (cs%thickness_diffuse .or. cs%interface_filter) &
2587 .and. (cs%dt_tr_adv /= cs%dt_therm))) then
2588 call get_param(param_file, "MOM", "INTERFACE_FILTER_DT_BUG", cs%interface_filter_dt_bug, &
2589 "If true, uses the wrong time interval in calls to interface_filter "//&
2590 "and thickness_diffuse. Has no effect when THICKNESSDIFFUSE_FIRST is "//&
2591 "true and DT_TRACER_ADVECT = DT_THERMO or when THICKNESSDIFFUSE_FIRST "//&
2592 "is false and APPLY_INTERFACE_FILTER is false. ", default=.false.)
2593 endif
2594
2595 if (bulkmixedlayer) then
2596 cs%Hmix = -1.0 ; cs%Hmix_UV = -1.0
2597 else
2598 call get_param(param_file, "MOM", "HMIX_SFC_PROP", hmix_z, &
2599 "If BULKMIXEDLAYER is false, HMIX_SFC_PROP is the depth "//&
2600 "over which to average to find surface properties like "//&
2601 "SST and SSS or density (but not surface velocities).", &
2602 units="m", default=1.0, scale=us%m_to_Z)
2603 call get_param(param_file, "MOM", "HMIX_UV_SFC_PROP", hmix_uv_z, &
2604 "If BULKMIXEDLAYER is false, HMIX_UV_SFC_PROP is the depth "//&
2605 "over which to average to find surface flow properties, "//&
2606 "SSU, SSV. A non-positive value indicates no averaging.", &
2607 units="m", default=0.0, scale=us%m_to_Z)
2608 endif
2609 call get_param(param_file, "MOM", "HFREEZE", hfrz_z, &
2610 "If HFREEZE > 0, melt potential will be computed. The actual depth "//&
2611 "over which melt potential is computed will be min(HFREEZE, OBLD), "//&
2612 "where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), "//&
2613 "melt potential will not be computed.", &
2614 units="m", default=-1.0, scale=us%m_to_Z)
2615 call get_param(param_file, "MOM", "INTERPOLATE_P_SURF", cs%interp_p_surf, &
2616 "If true, linearly interpolate the surface pressure "//&
2617 "over the coupling time step, using the specified value "//&
2618 "at the end of the step.", default=.false.)
2619
2620 if (cs%split) then
2621 call get_param(param_file, "MOM", "DTBT", dtbt, units="s or nondim", default=-0.98)
2622 default_val = us%T_to_s*cs%dt_therm ; if (dtbt > 0.0) default_val = -1.0
2623 cs%dtbt_reset_period = -1.0
2624 call get_param(param_file, "MOM", "DTBT_RESET_PERIOD", cs%dtbt_reset_period, &
2625 "The period between recalculations of DTBT (if DTBT <= 0). If "//&
2626 "DTBT_RESET_PERIOD is negative, DTBT is set based only on information "//&
2627 "available at initialization. If 0, DTBT will be set every dynamics time "//&
2628 "step. Values between 0 and DT are treated as 0. The default is set by "//&
2629 "DT_THERM. This is only used if SPLIT is true.", &
2630 units="s", default=default_val, scale=us%s_to_T, do_not_read=(dtbt > 0.0))
2631 endif
2632
2633 ! This is here in case these values are used inappropriately.
2634 use_frazil = .false. ; bound_salinity = .false. ; use_p_surf_in_eos = .false.
2635 cs%tv%P_Ref = 2.0e7*us%Pa_to_RL2_T2
2636 if (use_temperature) then
2637 call get_param(param_file, "MOM", "FRAZIL", use_frazil, &
2638 "If true, water freezes if it gets too cold, and the "//&
2639 "accumulated heat deficit is returned in the "//&
2640 "surface state. FRAZIL is only used if "//&
2641 "ENABLE_THERMODYNAMICS is true.", default=.false.)
2642 call get_param(param_file, "MOM", "DO_GEOTHERMAL", use_geothermal, &
2643 "If true, apply geothermal heating.", default=.false.)
2644 call get_param(param_file, "MOM", "BOUND_SALINITY", bound_salinity, &
2645 "If true, limit salinity to being positive. (The sea-ice "//&
2646 "model may ask for more salt than is available and "//&
2647 "drive the salinity negative otherwise.)", default=.false.)
2648 call get_param(param_file, "MOM", "MIN_SALINITY", cs%tv%min_salinity, &
2649 "The minimum value of salinity when BOUND_SALINITY=True.", &
2650 units="PPT", default=0.0, scale=us%ppt_to_S, do_not_log=.not.bound_salinity)
2651 call get_param(param_file, "MOM", "SALINITY_UNDERFLOW", salin_underflow, &
2652 "A tiny value of salinity below which the it is set to 0. For reference, "//&
2653 "one molecule of salt per square meter of ocean is of order 1e-29 ppt.", &
2654 units="PPT", default=0.0, scale=us%ppt_to_S)
2655 call get_param(param_file, "MOM", "TEMPERATURE_UNDERFLOW", temp_underflow, &
2656 "A tiny magnitude of temperatures below which they are set to 0.", &
2657 units="degC", default=0.0, scale=us%degC_to_C)
2658 call get_param(param_file, "MOM", "C_P", cs%tv%C_p, &
2659 "The heat capacity of sea water, approximated as a constant. "//&
2660 "This is only used if ENABLE_THERMODYNAMICS is true. The default "//&
2661 "value is from the TEOS-10 definition of conservative temperature.", &
2662 units="J kg-1 K-1", default=3991.86795711963, scale=us%J_kg_to_Q*us%C_to_degC)
2663 call get_param(param_file, "MOM", "USE_PSURF_IN_EOS", use_p_surf_in_eos, &
2664 "If true, always include the surface pressure contributions "//&
2665 "in equation of state calculations.", default=.true.)
2666 endif
2667 if (use_eos) call get_param(param_file, "MOM", "P_REF", cs%tv%P_Ref, &
2668 "The pressure that is used for calculating the coordinate "//&
2669 "density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
2670 "This is only used if USE_EOS and ENABLE_THERMODYNAMICS are true.", &
2671 units="Pa", default=2.0e7, scale=us%Pa_to_RL2_T2)
2672
2673 if (bulkmixedlayer) then
2674 call get_param(param_file, "MOM", "NKML", nkml, &
2675 "The number of sublayers within the mixed layer if "//&
2676 "BULKMIXEDLAYER is true.", units="nondim", default=2)
2677 call get_param(param_file, "MOM", "NKBL", nkbl, &
2678 "The number of layers that are used as variable density buffer "//&
2679 "layers if BULKMIXEDLAYER is true.", units="nondim", default=2)
2680 endif
2681
2682 call get_param(param_file, "MOM", "GLOBAL_INDEXING", global_indexing, &
2683 "If true, use a global lateral indexing convention, so "//&
2684 "that corresponding points on different processors have "//&
2685 "the same index. This does not work with static memory.", &
2686 default=.false., layoutparam=.true.)
2687#ifdef STATIC_MEMORY_
2688 if (global_indexing) call mom_error(fatal, "initialize_MOM: "//&
2689 "GLOBAL_INDEXING can not be true with STATIC_MEMORY.")
2690#endif
2691 call get_param(param_file, "MOM", "FIRST_DIRECTION", first_direction, &
2692 "An integer that indicates which direction goes first "//&
2693 "in parts of the code that use directionally split "//&
2694 "updates, with even numbers (or 0) used for x- first "//&
2695 "and odd numbers used for y-first.", default=0)
2696 call get_param(param_file, "MOM", "ALTERNATE_FIRST_DIRECTION", cs%alternate_first_direction, &
2697 "If true, after every dynamic timestep alternate whether the x- or y- "//&
2698 "direction updates occur first in directionally split parts of the calculation. "//&
2699 "If this is true, FIRST_DIRECTION applies at the start of a new run or if "//&
2700 "the next first direction can not be found in the restart file.", default=.false.)
2701 call get_param(param_file, "MOM", "CHECK_BAD_SURFACE_VALS", cs%check_bad_sfc_vals, &
2702 "If true, check the surface state for ridiculous values.", &
2703 default=.false.)
2704 if (cs%check_bad_sfc_vals) then
2705 call get_param(param_file, "MOM", "BAD_VAL_SSH_MAX", cs%bad_val_ssh_max, &
2706 "The value of SSH above which a bad value message is "//&
2707 "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
2708 units="m", default=20.0, scale=us%m_to_Z)
2709 call get_param(param_file, "MOM", "BAD_VAL_SSS_MAX", cs%bad_val_sss_max, &
2710 "The value of SSS above which a bad value message is "//&
2711 "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
2712 units="PPT", default=45.0, scale=us%ppt_to_S)
2713 call get_param(param_file, "MOM", "BAD_VAL_SST_MAX", cs%bad_val_sst_max, &
2714 "The value of SST above which a bad value message is "//&
2715 "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
2716 units="deg C", default=45.0, scale=us%degC_to_C)
2717 call get_param(param_file, "MOM", "BAD_VAL_SST_MIN", cs%bad_val_sst_min, &
2718 "The value of SST below which a bad value message is "//&
2719 "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
2720 units="deg C", default=-2.1, scale=us%degC_to_C)
2721 call get_param(param_file, "MOM", "BAD_VAL_COLUMN_THICKNESS", cs%bad_val_col_thick, &
2722 "The value of column thickness below which a bad value message is "//&
2723 "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
2724 units="m", default=0.0, scale=us%m_to_Z)
2725 endif
2726 call get_param(param_file, "MOM", "DEFAULT_ANSWER_DATE", default_answer_date, &
2727 "This sets the default value for the various _ANSWER_DATE parameters.", &
2728 default=99991231)
2729 call get_param(param_file, "MOM", "SURFACE_ANSWER_DATE", cs%answer_date, &
2730 "The vintage of the expressions for the surface properties. Values below "//&
2731 "20190101 recover the answers from the end of 2018, while higher values "//&
2732 "use updated and more robust forms of the same expressions.", &
2733 default=default_answer_date, do_not_log=non_bous)
2734 if (non_bous) cs%answer_date = 99991231
2735
2736 call get_param(param_file, "MOM", "SAVE_INITIAL_CONDS", save_ic, &
2737 "If true, write the initial conditions to a file given "//&
2738 "by IC_OUTPUT_FILE.", default=.false.)
2739 call get_param(param_file, "MOM", "IC_OUTPUT_FILE", cs%IC_file, &
2740 "The file into which to write the initial conditions.", &
2741 default="MOM_IC")
2742 call get_param(param_file, "MOM", "WRITE_GEOM", write_geom, &
2743 "If =0, never write the geometry and vertical grid files. "//&
2744 "If =1, write the geometry and vertical grid files only for "//&
2745 "a new simulation. If =2, always write the geometry and "//&
2746 "vertical grid files. Other values are invalid.", default=1)
2747 if (write_geom<0 .or. write_geom>2) call mom_error(fatal,"MOM: "//&
2748 "WRITE_GEOM must be equal to 0, 1 or 2.")
2749 call get_param(param_file, "MOM", "GEOM_FILE", geom_file, &
2750 "The file into which to write the ocean geometry.", &
2751 default="ocean_geometry")
2752 call get_param(param_file, "MOM", "USE_DBCLIENT", cs%use_dbclient, &
2753 "If true, initialize a client to a remote database that can "//&
2754 "be used for online analysis and machine-learning inference.",&
2755 default=.false.)
2756
2757 ! Check for inconsistent parameter settings.
2758 if (cs%use_ALE_algorithm .and. bulkmixedlayer) call mom_error(fatal, &
2759 "MOM: BULKMIXEDLAYER can not currently be used with the ALE algorithm.")
2760 if (cs%use_ALE_algorithm .and. .not.use_temperature) call mom_error(fatal, &
2761 "MOM: At this time, USE_EOS should be True when using the ALE algorithm.")
2762 if (cs%adiabatic .and. use_temperature) call mom_error(warning, &
2763 "MOM: ADIABATIC and ENABLE_THERMODYNAMICS both defined is usually unwise.")
2764 if (use_eos .and. .not.use_temperature) call mom_error(fatal, &
2765 "MOM: ENABLE_THERMODYNAMICS must be defined to use USE_EOS.")
2766 if (cs%adiabatic .and. bulkmixedlayer) call mom_error(fatal, &
2767 "MOM: ADIABATIC and BULKMIXEDLAYER can not both be defined.")
2768 if (bulkmixedlayer .and. .not.use_eos) call mom_error(fatal, &
2769 "initialize_MOM: A bulk mixed layer can only be used with T & S as "//&
2770 "state variables. Add USE_EOS = True to MOM_input.")
2771
2772 use_ice_shelf = .false.
2773 if (present(ice_shelf_csp)) then
2774 call get_param(param_file, "MOM", "ICE_SHELF", use_ice_shelf, &
2775 "If true, enables the ice shelf model.", default=.false.)
2776 endif
2777
2778 call get_param(param_file, "MOM", "USE_PARTICLES", cs%use_particles, &
2779 "If true, use the particles package.", default=.false.)
2780 call get_param(param_file, "MOM", "USE_UH_PARTICLES", cs%use_uh_particles, &
2781 "If true, use the uh velocity in the particles package.", &
2782 default=.false., do_not_log=.not.cs%use_particles)
2783 call get_param(param_file, "MOM", "UH_PARTICLES_BUG", cs%uh_particles_bug, &
2784 "If true, use a bug in which the particles are advected inconsistently"//&
2785 "with the dynamics timestep instead of the tracer timestep.", &
2786 default=enable_bugs, do_not_log=.not.cs%use_uh_particles)
2787 cs%ensemble_ocean=.false.
2788 call get_param(param_file, "MOM", "ENSEMBLE_OCEAN", cs%ensemble_ocean, &
2789 "If False, The model is being run in serial mode as a single realization. "//&
2790 "If True, The current model realization is part of a larger ensemble "//&
2791 "and at the end of step MOM, we will perform a gather of the ensemble "//&
2792 "members for statistical evaluation and/or data assimilation.", default=.false.)
2793
2794 call calltree_waypoint("MOM parameters read (initialize_MOM)")
2795
2796 call get_param(param_file, "MOM", "HOMOGENIZE_FORCINGS", cs%homogenize_forcings, &
2797 "If True, homogenize the forces and fluxes.", default=.false.)
2798 call get_param(param_file, "MOM", "UPDATE_USTAR",cs%update_ustar, &
2799 "If True, update ustar from homogenized tau when using the "//&
2800 "HOMOGENIZE_FORCINGS option. Note that this will not work "//&
2801 "with a non-zero gustiness factor.", default=.false., &
2802 do_not_log=.not.cs%homogenize_forcings)
2803
2804 ! Grid rotation test
2805 call get_param(param_file, "MOM", "ROTATE_INDEX", cs%rotate_index, &
2806 "Enable rotation of the horizontal indices.", default=.false., &
2807 debuggingparam=.true.)
2808 if (cs%rotate_index) then
2809 ! TODO: Index rotation currently only works when index rotation does not
2810 ! change the MPI rank of each domain. Resolving this will require a
2811 ! modification to FMS PE assignment.
2812 ! For now, we only permit single-core runs.
2813
2814 if (num_pes() /= 1) &
2815 call mom_error(fatal, "Index rotation is only supported on one PE.")
2816
2817 ! Alternate_first_direction is not permitted with index rotation.
2818 ! This feature can be added later in the future if needed.
2819 if (cs%alternate_first_direction) &
2820 call mom_error(fatal, "Alternating_first_direction is not compatible with index rotation.")
2821
2822 call get_param(param_file, "MOM", "INDEX_TURNS", turns, &
2823 "Number of counterclockwise quarter-turn index rotations.", &
2824 default=1, debuggingparam=.true.)
2825 else
2826 turns = 0
2827 endif
2828
2829 ! Set up the model domain and grids.
2830#ifdef SYMMETRIC_MEMORY_
2831 symmetric = .true.
2832#else
2833 symmetric = .false.
2834#endif
2835 g_in => cs%G_in
2836#ifdef STATIC_MEMORY_
2837 call mom_domains_init(g_in%domain, param_file, symmetric=symmetric, &
2838 static_memory=.true., nihalo=nihalo_, njhalo=njhalo_, &
2839 niglobal=niglobal_, njglobal=njglobal_, niproc=niproc_, &
2840 njproc=njproc_, us=us, mom_dom_unmasked=mom_dom_unmasked)
2841#else
2842 call mom_domains_init(g_in%domain, param_file, symmetric=symmetric, &
2843 domain_name="MOM_in", us=us, mom_dom_unmasked=mom_dom_unmasked)
2844#endif
2845
2846 ! Copy input grid (G_in) domain to active grid G
2847 ! Swap axes for quarter and 3-quarter turns
2848 if (cs%rotate_index) then
2849 allocate(cs%G)
2850 call clone_mom_domain(g_in%Domain, cs%G%Domain, turns=turns, domain_name="MOM_rot")
2851 else
2852 cs%G => g_in
2853 endif
2854
2855 ! TODO: It is unlikely that test_grid_copy and rotate_index would work at the
2856 ! same time. It may be possible to enable both but for now we prevent it.
2857 if (test_grid_copy .and. cs%rotate_index) &
2858 call mom_error(fatal, "Grid cannot be copied during index rotation.")
2859
2860 if (test_grid_copy) then ; allocate(g)
2861 else ; g => cs%G ; endif
2862
2863 call calltree_waypoint("domains initialized (initialize_MOM)")
2864
2865 call mom_debugging_init(param_file)
2866 call diag_mediator_infrastructure_init()
2867 call mom_io_init(param_file)
2868
2869 ! Create HI and dG on the input index map.
2870 call hor_index_init(g_in%Domain, hi_in, param_file, &
2871 local_indexing=.not.global_indexing)
2872 call create_dyn_horgrid(dg_in, hi_in, bathymetry_at_vel=bathy_at_vel)
2873 call clone_mom_domain(g_in%Domain, dg_in%Domain)
2874 ! Also allocate the input ocean_grid_type type at this point based on the same information.
2875 call mom_grid_init(g_in, param_file, us, hi_in, bathymetry_at_vel=bathy_at_vel)
2876
2877 ! Allocate initialize time-invariant MOM variables.
2878 call mom_initialize_fixed(dg_in, us, obc_in, param_file)
2879
2880 call get_param(param_file, "MOM", "DT_OBC_SEG_UPDATE_OBGC", cs%dt_obc_seg_period, &
2881 "The time between OBC segment data updates for OBGC tracers. This must be an "//&
2882 "integer multiple of DT and DT_THERM. The default is set to DT.", units="s", &
2883 default=us%T_to_s*cs%dt, scale=us%s_to_T, do_not_log=.not.associated(obc_in))
2884 call get_param(param_file, "MOM", "OBC_BGC_TIME_REF_BUG", obc_bgc_time_ref_bug, &
2885 "If true, recover a bug that the BGC OBC segment update schedule is "//&
2886 "referenced to the start of the current run rather than the overall start "//&
2887 "time, which can lead to restart reproducibility failures.", &
2888 default=.true., do_not_log=.not.associated(obc_in))
2889
2890 ! Copy the grid metrics and bathymetry to the ocean_grid_type
2891 call copy_dyngrid_to_mom_grid(dg_in, g_in, us)
2892
2893 call calltree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)")
2894
2895 call verticalgridinit( param_file, cs%GV, us )
2896 gv => cs%GV
2897
2898 ! Now that the vertical grid has been initialized, rescale parameters that depend on factors
2899 ! that are set with the vertical grid to their desired units. This added rescaling step would
2900 ! be unnecessary if the vertical grid were initialized earlier in this routine.
2901 if (.not.bulkmixedlayer) then
2902 cs%Hmix = (us%Z_to_m * gv%m_to_H) * hmix_z
2903 cs%Hmix_UV = (us%Z_to_m * gv%m_to_H) * hmix_uv_z
2904 endif
2905 cs%HFrz = (us%Z_to_m * gv%m_to_H) * hfrz_z
2906
2907 ! Shift from using the temporary dynamic grid type to using the final (potentially static)
2908 ! and properly rotated ocean-specific grid type and horizontal index type.
2909 if (cs%rotate_index) then
2910 allocate(hi)
2911 call rotate_hor_index(hi_in, turns, hi)
2912 ! NOTE: If indices are rotated, then G and G_in must both be initialized separately, and
2913 ! the dynamic grid must be created to handle the grid rotation. G%domain has already been
2914 ! initialized above.
2915 call mom_grid_init(g, param_file, us, hi, bathymetry_at_vel=bathy_at_vel)
2916 call create_dyn_horgrid(dg, hi, bathymetry_at_vel=bathy_at_vel)
2917 call clone_mom_domain(g%Domain, dg%Domain)
2918 call rotate_dyn_horgrid(dg_in, dg, us, turns)
2919 call copy_dyngrid_to_mom_grid(dg, g, us)
2920
2921 if (associated(obc_in)) then
2922 allocate(cs%OBC)
2923 call rotate_obc_config(obc_in, dg_in, cs%OBC, dg, turns)
2924 endif
2925
2926 call destroy_dyn_horgrid(dg)
2927 else
2928 ! If not rotated, then G_in and G are the same grid.
2929 hi => hi_in
2930 g => g_in
2931 cs%OBC => obc_in
2932 endif
2933 ! dG_in is retained for now so that it can be used with write_ocean_geometry_file() below.
2934
2935 if (is_root_pe()) call check_mom6_scaling_factors(cs%GV, us)
2936
2937 call calltree_waypoint("grids initialized (initialize_MOM)")
2938
2939 call mom_timing_init(cs)
2940
2941 call tracer_registry_init(param_file, cs%tracer_Reg)
2942
2943 ! Allocate and initialize space for the primary time-varying MOM variables.
2944 is = hi%isc ; ie = hi%iec ; js = hi%jsc ; je = hi%jec ; nz = gv%ke
2945 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed
2946 isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
2947 alloc_(cs%u(isdb:iedb,jsd:jed,nz)) ; cs%u(:,:,:) = 0.0
2948 alloc_(cs%v(isd:ied,jsdb:jedb,nz)) ; cs%v(:,:,:) = 0.0
2949 alloc_(cs%h(isd:ied,jsd:jed,nz)) ; cs%h(:,:,:) = gv%Angstrom_H
2950 alloc_(cs%uh(isdb:iedb,jsd:jed,nz)) ; cs%uh(:,:,:) = 0.0
2951 alloc_(cs%vh(isd:ied,jsdb:jedb,nz)) ; cs%vh(:,:,:) = 0.0
2952 if (use_temperature) then
2953 alloc_(cs%T(isd:ied,jsd:jed,nz)) ; cs%T(:,:,:) = 0.0
2954 alloc_(cs%S(isd:ied,jsd:jed,nz)) ; cs%S(:,:,:) = 0.0
2955 cs%tv%T => cs%T ; cs%tv%S => cs%S
2956 if (cs%tv%T_is_conT) then
2957 vd_t = var_desc(name="contemp", units="Celsius", longname="Conservative Temperature", &
2958 cmor_field_name="bigthetao", cmor_longname="Sea Water Conservative Temperature", &
2959 conversion=us%C_to_degC)
2960 else
2961 vd_t = var_desc(name="temp", units="degC", longname="Potential Temperature", &
2962 cmor_field_name="thetao", cmor_longname="Sea Water Potential Temperature", &
2963 conversion=us%C_to_degC)
2964 endif
2965 if (cs%tv%S_is_absS) then
2966 vd_s = var_desc(name="abssalt", units="g kg-1", longname="Absolute Salinity", &
2967 cmor_field_name="absso", cmor_longname="Sea Water Absolute Salinity", &
2968 conversion=us%S_to_ppt)
2969 else
2970 vd_s = var_desc(name="salt", units="psu", longname="Salinity", &
2971 cmor_field_name="so", cmor_longname="Sea Water Salinity", &
2972 conversion=us%S_to_ppt)
2973 endif
2974
2975 if (advect_ts) then
2976 s_flux_units = get_tr_flux_units(gv, "psu") ! Could change to "kg m-2 s-1"?
2977 conv2watt = gv%H_to_kg_m2 * us%Q_to_J_kg*cs%tv%C_p
2978 if (gv%Boussinesq) then
2979 conv2salt = us%S_to_ppt*gv%H_to_m ! Could change to US%S_to_ppt*GV%H_to_kg_m2 * 0.001?
2980 else
2981 conv2salt = us%S_to_ppt*gv%H_to_kg_m2
2982 endif
2983 call register_tracer(cs%tv%T, cs%tracer_Reg, param_file, hi, gv, &
2984 tr_desc=vd_t, registry_diags=.true., conc_scale=us%C_to_degC, &
2985 flux_nameroot='T', flux_units='W', flux_longname='Heat', &
2986 net_surfflux_name='KPP_QminusSW', nlt_budget_name='KPP_NLT_temp_budget', &
2987 net_surfflux_longname='Net temperature flux ignoring short-wave, as used by [CVMix] KPP', &
2988 flux_scale=conv2watt, convergence_units='W m-2', &
2989 convergence_scale=conv2watt, cmor_tendprefix="opottemp", &
2990 diag_form=2, underflow_conc=temp_underflow, tr_out=cs%tv%tr_T)
2991 call register_tracer(cs%tv%S, cs%tracer_Reg, param_file, hi, gv, &
2992 tr_desc=vd_s, registry_diags=.true., conc_scale=us%S_to_ppt, &
2993 flux_nameroot='S', flux_units=s_flux_units, flux_longname='Salt', &
2994 net_surfflux_name='KPP_netSalt', nlt_budget_name='KPP_NLT_saln_budget', &
2995 flux_scale=conv2salt, convergence_units='kg m-2 s-1', &
2996 convergence_scale=0.001*us%S_to_ppt*gv%H_to_kg_m2, cmor_tendprefix="osalt", &
2997 diag_form=2, underflow_conc=salin_underflow, tr_out=cs%tv%tr_S)
2998 endif
2999 endif
3000
3001 if (use_p_surf_in_eos) allocate(cs%tv%p_surf(isd:ied,jsd:jed), source=0.0)
3002 if (use_frazil) then
3003 allocate(cs%tv%frazil(isd:ied,jsd:jed), source=0.0)
3004 cs%tv%frazil_was_reset = .true.
3005 endif
3006 if (bound_salinity) allocate(cs%tv%salt_deficit(isd:ied,jsd:jed), source=0.0)
3007
3008 allocate(cs%Hml(isd:ied,jsd:jed), source=0.0)
3009
3010 if (bulkmixedlayer) then
3011 gv%nkml = nkml ; gv%nk_rho_varies = nkml + nkbl
3012 else
3013 gv%nkml = 0 ; gv%nk_rho_varies = 0
3014 endif
3015 if (cs%use_ALE_algorithm) then
3016 call get_param(param_file, "MOM", "NK_RHO_VARIES", gv%nk_rho_varies, default=0) ! Will default to nz later... -AJA
3017 endif
3018
3019 alloc_(cs%uhtr(isdb:iedb,jsd:jed,nz)) ; cs%uhtr(:,:,:) = 0.0
3020 alloc_(cs%vhtr(isd:ied,jsdb:jedb,nz)) ; cs%vhtr(:,:,:) = 0.0
3021 cs%t_dyn_rel_adv = 0.0 ; cs%t_dyn_rel_thermo = 0.0 ; cs%t_dyn_rel_diag = 0.0
3022 cs%n_dyn_steps_in_adv = 0
3023
3024 if (debug_truncations) then
3025 allocate(cs%u_prev(isdb:iedb,jsd:jed,nz), source=0.0)
3026 allocate(cs%v_prev(isd:ied,jsdb:jedb,nz), source=0.0)
3027 mom_internal_state%u_prev => cs%u_prev
3028 mom_internal_state%v_prev => cs%v_prev
3029 call safe_alloc_ptr(cs%ADp%du_dt_visc,isdb,iedb,jsd,jed,nz)
3030 call safe_alloc_ptr(cs%ADp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
3031 if (.not.cs%adiabatic) then
3032 call safe_alloc_ptr(cs%ADp%du_dt_dia,isdb,iedb,jsd,jed,nz)
3033 call safe_alloc_ptr(cs%ADp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
3034 endif
3035 endif
3036
3037 mom_internal_state%u => cs%u ; mom_internal_state%v => cs%v
3038 mom_internal_state%h => cs%h
3039 mom_internal_state%uh => cs%uh ; mom_internal_state%vh => cs%vh
3040 if (use_temperature) then
3041 mom_internal_state%T => cs%T ; mom_internal_state%S => cs%S
3042 endif
3043
3044 cs%CDp%uh => cs%uh ; cs%CDp%vh => cs%vh
3045
3046 if (cs%interp_p_surf) allocate(cs%p_surf_prev(isd:ied,jsd:jed), source=0.0)
3047
3048 alloc_(cs%ssh_rint(isd:ied,jsd:jed)) ; cs%ssh_rint(:,:) = 0.0
3049 alloc_(cs%ave_ssh_ibc(isd:ied,jsd:jed)) ; cs%ave_ssh_ibc(:,:) = 0.0
3050 alloc_(cs%eta_av_bc(isd:ied,jsd:jed)) ; cs%eta_av_bc(:,:) = 0.0 ! -G%Z_ref
3051 cs%time_in_cycle = 0.0 ; cs%time_in_thermo_cycle = 0.0
3052
3053 !allocate porous topography variables
3054 allocate(cs%pbv%por_face_areaU(isdb:iedb,jsd:jed,nz), source=1.0)
3055 allocate(cs%pbv%por_face_areaV(isd:ied,jsdb:jedb,nz), source=1.0)
3056 allocate(cs%pbv%por_layer_widthU(isdb:iedb,jsd:jed,nz+1), source=1.0)
3057 allocate(cs%pbv%por_layer_widthV(isd:ied,jsdb:jedb,nz+1), source=1.0)
3058
3059 ! Use the Wright equation of state by default, unless otherwise specified
3060 ! Note: this line and the following block ought to be in a separate
3061 ! initialization routine for tv.
3062 if (use_eos) then
3063 allocate(cs%tv%eqn_of_state)
3064 call eos_init(param_file, cs%tv%eqn_of_state, us, use_cont_abss)
3065 endif
3066 if (use_temperature) then
3067 allocate(cs%tv%TempxPmE(isd:ied,jsd:jed), source=0.0)
3068 if (use_geothermal) then
3069 allocate(cs%tv%internal_heat(isd:ied,jsd:jed), source=0.0)
3070 endif
3071 endif
3072 call calltree_waypoint("state variables allocated (initialize_MOM)")
3073
3074 ! Set the fields that are needed for bitwise identical restarting
3075 ! the time stepping scheme.
3076 call restart_init(param_file, cs%restart_CS)
3077 restart_csp => cs%restart_CS
3078
3079 call set_restart_fields(gv, us, param_file, cs, restart_csp)
3080 if (cs%split .and. cs%use_alt_split) then
3081 call register_restarts_dyn_split_rk2b(hi, gv, us, param_file, &
3082 cs%dyn_split_RK2b_CSp, restart_csp, cs%uh, cs%vh)
3083 elseif (cs%split) then
3084 call register_restarts_dyn_split_rk2(hi, gv, us, param_file, &
3085 cs%dyn_split_RK2_CSp, restart_csp, cs%uh, cs%vh)
3086 elseif (cs%use_RK2) then
3087 call register_restarts_dyn_unsplit_rk2(hi, gv, param_file, &
3088 cs%dyn_unsplit_RK2_CSp)
3089 else
3090 call register_restarts_dyn_unsplit(hi, gv, param_file, &
3091 cs%dyn_unsplit_CSp)
3092 endif
3093
3094 ! This subroutine calls user-specified tracer registration routines.
3095 ! Additional calls can be added to MOM_tracer_flow_control.F90.
3096 call call_tracer_register(g, gv, us, param_file, cs%tracer_flow_CSp, &
3097 cs%tracer_Reg, restart_csp)
3098
3099 call meke_alloc_register_restart(hi, us, param_file, cs%MEKE, restart_csp)
3100 call set_visc_register_restarts(hi, g, gv, us, param_file, cs%visc, restart_csp, use_ice_shelf)
3101 call mixedlayer_restrat_register_restarts(hi, gv, us, param_file, &
3102 cs%mixedlayer_restrat_CSp, restart_csp)
3103
3104 if (associated(cs%OBC)) then
3105 ! This call initializes the relevant vertical remapping structures.
3106 call open_boundary_setup_vert(gv, us, cs%OBC)
3107
3108 ! Set up remaining information about open boundary conditions that is needed for OBCs.
3109 ! Package specific changes to OBCs occur here.
3110 call call_obc_register(g, gv, us, param_file, cs%update_OBC_CSp, cs%OBC, cs%tracer_Reg)
3111
3112 ! This is the equivalent to 2 calls to register_segment_tracer (per segment), which
3113 ! could occur with the call to update_OBC_data or after the main initialization.
3114 if (use_temperature) &
3115 call register_temp_salt_segments(gv, us, cs%OBC, cs%tracer_Reg, param_file)
3116 ! This is the equivalent call to register_temp_salt_segments for external tracers with OBC
3117 call call_tracer_register_obc_segments(gv, param_file, cs%tracer_flow_CSp, cs%tracer_Reg, cs%OBC)
3118
3119 ! Set up the thickness reservoirs if using them.
3120 if (cs%OBC%use_h_res) &
3121 call segment_thickness_reservoir_init(gv, us, cs%OBC, param_file)
3122
3123 ! This needs the number of tracers and to have called any code that sets whether
3124 ! reservoirs are used.
3125 call open_boundary_register_restarts(hi, gv, us, cs%OBC, cs%tracer_Reg, &
3126 param_file, restart_csp, use_temperature)
3127
3128 ! This call allocates the arrays on the segments for open boundary data, but it must occur
3129 ! after any calls to call_tracer_register_obc_segments.
3130 call initialize_segment_data(gv, us, cs%OBC, param_file, turns, use_temperature)
3131
3132 if (cs%debug_OBCs) call write_obc_info(cs%OBC, g, gv, us)
3133 endif
3134
3135 if (present(waves_csp)) then
3136 call waves_register_restarts(waves_csp, hi, gv, us, param_file, restart_csp)
3137 endif
3138
3139 if (use_temperature) then
3140 call stoch_eos_register_restarts(hi, param_file, cs%stoch_eos_CS, restart_csp)
3141 endif
3142
3143 if (.not. cs%adiabatic) then
3144 call register_diabatic_restarts(g, gv, us, param_file, cs%int_tide_CSp, restart_csp, cs%diabatic_CSp)
3145 endif
3146
3147 call calltree_waypoint("restart registration complete (initialize_MOM)")
3148 call restart_registry_lock(restart_csp)
3149
3150 ! Write out all of the grid data used by this run.
3151 new_sim = determine_is_new_run(dirs%input_filename, dirs%restart_input_dir, g_in, restart_csp)
3152 write_geom_files = ((write_geom==2) .or. ((write_geom==1) .and. new_sim))
3153 if (write_geom_files) then
3154 if (associated(mom_dom_unmasked)) then
3155 call hor_index_init(mom_dom_unmasked, hi_in_unmasked, param_file, &
3156 local_indexing=.not.global_indexing)
3157 call create_dyn_horgrid(dg_unmasked_in, hi_in_unmasked, bathymetry_at_vel=bathy_at_vel)
3158 call clone_mom_domain(mom_dom_unmasked, dg_unmasked_in%Domain)
3159 call mom_initialize_fixed(dg_unmasked_in, us, obc_in, param_file)
3160 call write_ocean_geometry_file(dg_unmasked_in, param_file, dirs%output_directory, us=us, geom_file=geom_file)
3161 call deallocate_mom_domain(mom_dom_unmasked)
3162 call destroy_dyn_horgrid(dg_unmasked_in)
3163 else
3164 call write_ocean_geometry_file(dg_in, param_file, dirs%output_directory, us=us, geom_file=geom_file)
3165 endif
3166 endif
3167 call destroy_dyn_horgrid(dg_in)
3168
3169 ! Initialize dynamically evolving fields, perhaps from restart files.
3170 call cpu_clock_begin(id_clock_mom_init)
3171 call mom_initialize_coord(gv, us, param_file, cs%tv, g%max_depth)
3172 call calltree_waypoint("returned from MOM_initialize_coord() (initialize_MOM)")
3173
3174 if (cs%use_ALE_algorithm) then
3175 call ale_init(param_file, g, gv, us, g%max_depth, cs%ALE_CSp)
3176 call calltree_waypoint("returned from ALE_init() (initialize_MOM)")
3177 endif
3178
3179 ! Set a few remaining fields that are specific to the ocean grid type.
3180 if (cs%rotate_index) then
3181 call set_first_direction(g, modulo(first_direction + turns, 2))
3182 else
3183 call set_first_direction(g, modulo(first_direction, 2))
3184 endif
3185 ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
3186 if (cs%debug .or. g%symmetric) then
3187 call clone_mom_domain(g%Domain, g%Domain_aux, symmetric=.false.)
3188 else ; g%Domain_aux => g%Domain ; endif
3189 ! Copy common variables from the vertical grid to the horizontal grid.
3190 ! Consider removing this later?
3191 g%ke = gv%ke
3192
3193 if (use_ice_shelf) then
3194 point_calving = .false. ; if (present(calve_ice_shelf_bergs)) point_calving = calve_ice_shelf_bergs
3195 endif
3196
3197 if (cs%rotate_index) then
3198 g_in%ke = gv%ke
3199
3200 ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
3201 if (cs%debug .or. g_in%symmetric) then
3202 call clone_mom_domain(g_in%Domain, g_in%Domain_aux, symmetric=.false.)
3203 else ; g_in%Domain_aux => g_in%Domain ; endif
3204
3205 allocate(u_in(g_in%IsdB:g_in%IedB, g_in%jsd:g_in%jed, nz), source=0.0)
3206 allocate(v_in(g_in%isd:g_in%ied, g_in%JsdB:g_in%JedB, nz), source=0.0)
3207 allocate(h_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed, nz), source=gv%Angstrom_H)
3208
3209 if (use_temperature) then
3210 allocate(t_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed, nz), source=0.0)
3211 allocate(s_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed, nz), source=0.0)
3212
3213 cs%tv%T => t_in
3214 cs%tv%S => s_in
3215
3216 if (associated(cs%OBC)) then
3217 ! Log this parameter in MOM_initialize_state
3218 call get_param(param_file, "MOM", "OBC_RESERVOIR_INIT_BUG", obc_reservoir_init_bug, &
3219 "If true, set the OBC tracer reservoirs at the startup of a new run from the "//&
3220 "interior tracer concentrations regardless of properties that may be explicitly "//&
3221 "specified for the reservoir concentrations.", default=enable_bugs, do_not_log=.true.)
3222 if (obc_reservoir_init_bug .and. (allocated(cs%OBC%tres_x) .or. allocated(cs%OBC%tres_y))) &
3223 call mom_error(fatal, "OBC_RESERVOIR_INIT_BUG can not be set to true with grid rotation.")
3224 endif
3225 endif
3226
3227 if (use_ice_shelf) then
3228 ! These arrays are not initialized in most solo cases, but are needed
3229 ! when using an ice shelf. Passing the ice shelf diagnostics CS from MOM
3230 ! for legacy reasons. The actual ice shelf diag CS is internal to the ice shelf
3231 call initialize_ice_shelf(param_file, g, time, ice_shelf_csp, diag_ptr, &
3232 time_init, dirs%output_directory, calve_ice_shelf_bergs=point_calving)
3233 allocate(frac_shelf_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed), source=0.0)
3234 allocate(mass_shelf_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed), source=0.0)
3235 allocate(cs%frac_shelf_h(isd:ied, jsd:jed), source=0.0)
3236 allocate(cs%mass_shelf(isd:ied, jsd:jed), source=0.0)
3237 call ice_shelf_query(ice_shelf_csp, g, cs%frac_shelf_h, cs%mass_shelf)
3238 ! MOM_initialize_state is using the unrotated metric
3239 call rotate_array(cs%frac_shelf_h, -turns, frac_shelf_in)
3240 call rotate_array(cs%mass_shelf, -turns, mass_shelf_in)
3241 call mom_initialize_state(u_in, v_in, h_in, cs%tv, time, g_in, gv, us, &
3242 param_file, dirs, restart_csp, cs%ALE_CSp, cs%tracer_Reg, &
3243 sponge_in_csp, ale_sponge_in_csp, oda_incupd_in_csp, obc_in, time_in, &
3244 frac_shelf_h=frac_shelf_in, mass_shelf=mass_shelf_in)
3245 else
3246 call mom_initialize_state(u_in, v_in, h_in, cs%tv, time, g_in, gv, us, &
3247 param_file, dirs, restart_csp, cs%ALE_CSp, cs%tracer_Reg, &
3248 sponge_in_csp, ale_sponge_in_csp, oda_incupd_in_csp, obc_in, time_in)
3249 endif
3250
3251 if (use_temperature) then
3252 cs%tv%T => cs%T
3253 cs%tv%S => cs%S
3254 endif
3255
3256 ! Reset the first direction if it was found in a restart file
3257 if (cs%first_dir_restart > -1.0) then
3258 call set_first_direction(g, modulo(nint(cs%first_dir_restart) + turns, 2))
3259 else
3260 cs%first_dir_restart = real(modulo(first_direction, 2))
3261 endif
3262
3263 call rotate_initial_state(u_in, v_in, h_in, t_in, s_in, use_temperature, &
3264 turns, cs%u, cs%v, cs%h, cs%T, cs%S)
3265
3266 if (associated(sponge_in_csp)) then
3267 ! TODO: Implementation and testing of non-ALE sponge rotation
3268 call mom_error(fatal, "Index rotation of non-ALE sponge is not yet implemented.")
3269 endif
3270
3271 if (associated(ale_sponge_in_csp)) then
3272 call rotate_ale_sponge(ale_sponge_in_csp, g_in, cs%ALE_sponge_CSp, g, gv, us, turns, param_file)
3273 call update_ale_sponge_field(cs%ALE_sponge_CSp, t_in, g, gv, cs%T)
3274 call update_ale_sponge_field(cs%ALE_sponge_CSp, s_in, g, gv, cs%S)
3275 endif
3276
3277 ! Deallocate the unrotated arrays and types that are no longer needed.
3278 deallocate(u_in)
3279 deallocate(v_in)
3280 deallocate(h_in)
3281 if (use_temperature) then
3282 deallocate(t_in)
3283 deallocate(s_in)
3284 endif
3285 if (use_ice_shelf) deallocate(frac_shelf_in, mass_shelf_in)
3286 if (associated(obc_in)) call open_boundary_end(obc_in)
3287
3288 else ! The model is being run without grid rotation. This is true of all production runs.
3289 if (use_ice_shelf) then
3290 call initialize_ice_shelf(param_file, g, time, ice_shelf_csp, diag_ptr, time_init, &
3291 dirs%output_directory, calve_ice_shelf_bergs=point_calving)
3292 allocate(cs%frac_shelf_h(isd:ied, jsd:jed), source=0.0)
3293 allocate(cs%mass_shelf(isd:ied, jsd:jed), source=0.0)
3294 call ice_shelf_query(ice_shelf_csp,g,cs%frac_shelf_h, cs%mass_shelf)
3295 call mom_initialize_state(cs%u, cs%v, cs%h, cs%tv, time, g, gv, us, &
3296 param_file, dirs, restart_csp, cs%ALE_CSp, cs%tracer_Reg, &
3297 cs%sponge_CSp, cs%ALE_sponge_CSp, cs%oda_incupd_CSp, cs%OBC, time_in, &
3298 frac_shelf_h=cs%frac_shelf_h, mass_shelf=cs%mass_shelf, obc_for_bug=cs%OBC)
3299 else
3300 call mom_initialize_state(cs%u, cs%v, cs%h, cs%tv, time, g, gv, us, &
3301 param_file, dirs, restart_csp, cs%ALE_CSp, cs%tracer_Reg, &
3302 cs%sponge_CSp, cs%ALE_sponge_CSp, cs%oda_incupd_CSp, cs%OBC, time_in, obc_for_bug=cs%OBC)
3303 endif
3304
3305 ! Reset the first direction if it was found in a restart file.
3306 if (cs%first_dir_restart > -1.0) then
3307 call set_first_direction(g, nint(cs%first_dir_restart))
3308 else
3309 cs%first_dir_restart = real(modulo(first_direction, 2))
3310 endif
3311 endif
3312
3313 ! Allocate any derived densities or other equation of state derived fields.
3314 if (.not.(gv%Boussinesq .or. gv%semi_Boussinesq)) then
3315 allocate(cs%tv%SpV_avg(isd:ied,jsd:jed,nz), source=0.0)
3316 cs%tv%valid_SpV_halo = -1 ! This array does not yet have any valid data.
3317 endif
3318
3319 if (associated(cs%OBC)) then
3320 call mom_initialize_obcs(cs%h, cs%tv, cs%OBC, time, g, gv, us, param_file, restart_csp, cs%tracer_Reg)
3321
3322 if (use_temperature) then
3323 call pass_var(cs%tv%T, g%Domain, complete=.false.)
3324 call pass_var(cs%tv%S, g%Domain, complete=.true.)
3325 endif
3326 call calc_derived_thermo(cs%tv, cs%h, g, gv, us)
3327
3328 ! Call this during initialization to fill boundary arrays from fixed values
3329 call read_obc_segment_data(g, gv, us, cs%OBC, cs%tv, cs%h, time)
3330 call update_obc_segment_data(g, gv, us, cs%OBC, cs%h, time)
3331 call initialize_obc_segment_reservoirs(gv, cs%OBC)
3332 endif
3333
3334 if (use_ice_shelf .and. cs%debug) then
3335 call hchksum(cs%frac_shelf_h, "MOM:frac_shelf_h", g%HI, haloshift=0)
3336 call hchksum(cs%mass_shelf, "MOM:mass_shelf", g%HI, haloshift=0, unscale=us%RZ_to_kg_m2)
3337 endif
3338
3339 call cpu_clock_end(id_clock_mom_init)
3340 call calltree_waypoint("returned from MOM_initialize_state() (initialize_MOM)")
3341
3342 ! From this point, there may be pointers being set, so the final grid type
3343 ! that will persist throughout the run has to be used.
3344
3345 if (test_grid_copy) then
3346 ! Copy the data from the temporary grid to the dyn_hor_grid to CS%G.
3347 call create_dyn_horgrid(test_dg, g%HI)
3348 call clone_mom_domain(g%Domain, test_dg%Domain)
3349
3350 call clone_mom_domain(g%Domain, cs%G%Domain)
3351 call mom_grid_init(cs%G, param_file, us)
3352
3353 call copy_mom_grid_to_dyngrid(g, test_dg, us)
3354 call copy_dyngrid_to_mom_grid(test_dg, cs%G, us)
3355
3356 call destroy_dyn_horgrid(test_dg)
3357 call mom_grid_end(g) ; deallocate(g)
3358
3359 g => cs%G
3360 if (cs%debug .or. cs%G%symmetric) then
3361 call clone_mom_domain(cs%G%Domain, cs%G%Domain_aux, symmetric=.false.)
3362 else ; cs%G%Domain_aux => cs%G%Domain ; endif
3363 g%ke = gv%ke
3364 endif
3365
3366 ! At this point, all user-modified initialization code has been called. The
3367 ! remainder of this subroutine is controlled by the parameters that have
3368 ! have already been set.
3369
3370 if (ale_remap_init_conds(cs%ALE_CSp) .and. .not. query_initialized(cs%h,"h",restart_csp)) then
3371 ! This block is controlled by the ALE parameter REMAP_AFTER_INITIALIZATION.
3372 ! \todo This block exists for legacy reasons and we should phase it out of all examples. !###
3373 if (cs%debug) then
3374 call uvchksum("Pre ALE adjust init cond [uv]", cs%u, cs%v, g%HI, haloshift=1, unscale=us%L_T_to_m_s)
3375 call hchksum(cs%h,"Pre ALE adjust init cond h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
3376 endif
3377 call calltree_waypoint("Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)")
3378 call adjustgridforintegrity(cs%ALE_CSp, g, gv, cs%h )
3379 if (allocated(cs%tv%SpV_avg)) call calc_derived_thermo(cs%tv, cs%h, g, gv, us, halo=1)
3380 call pre_ale_adjustments(g, gv, us, cs%h, cs%tv, cs%tracer_Reg, cs%ALE_CSp, cs%u, cs%v)
3381
3382 call calltree_waypoint("Calling ALE_regrid() to remap initial conditions (initialize_MOM)")
3383 allocate(h_new(isd:ied, jsd:jed, nz), source=0.0)
3384 allocate(dzregrid(isd:ied, jsd:jed, nz+1), source=0.0)
3385 allocate(pcm_cell(isd:ied, jsd:jed, nz), source=.false.)
3386 allocate(h_old_u(isdb:iedb, jsd:jed, nz), source=0.0)
3387 allocate(h_new_u(isdb:iedb, jsd:jed, nz), source=0.0)
3388 allocate(h_old_v(isd:ied, jsdb:jedb, nz), source=0.0)
3389 allocate(h_new_v(isd:ied, jsdb:jedb, nz), source=0.0)
3390 if (use_ice_shelf) then
3391 call ale_regrid(g, gv, us, cs%h, h_new, dzregrid, cs%tv, cs%ALE_CSp, cs%frac_shelf_h, pcm_cell)
3392 else
3393 call ale_regrid(g, gv, us, cs%h, h_new, dzregrid, cs%tv, cs%ALE_CSp, pcm_cell=pcm_cell)
3394 endif
3395
3396 if (calltree_showquery()) call calltree_waypoint("new grid generated")
3397 ! Remap all variables from the old grid h onto the new grid h_new
3398 call ale_remap_tracers(cs%ALE_CSp, g, gv, cs%h, h_new, cs%tracer_Reg, cs%debug, pcm_cell=pcm_cell)
3399
3400 ! Determine the old and new grid thicknesses at velocity points.
3401 call ale_remap_set_h_vel(cs%ALE_CSp, g, gv, cs%h, h_old_u, h_old_v, cs%OBC, debug=cs%debug)
3402 if (cs%remap_uv_using_old_alg) then
3403 call ale_remap_set_h_vel_via_dz(cs%ALE_CSp, g, gv, h_new, h_new_u, h_new_v, cs%OBC, cs%h, dzregrid, cs%debug)
3404 else
3405 call ale_remap_set_h_vel(cs%ALE_CSp, g, gv, h_new, h_new_u, h_new_v, cs%OBC, debug=cs%debug)
3406 endif
3407
3408 ! Remap the velocity components.
3409 call ale_remap_velocities(cs%ALE_CSp, g, gv, h_old_u, h_old_v, h_new_u, h_new_v, cs%u, cs%v, cs%debug)
3410
3411 if (allocated(cs%tv%SpV_avg)) cs%tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid.
3412
3413 ! Replace the old grid with new one. All remapping must be done at this point.
3414 !$OMP parallel do default(shared)
3415 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
3416 cs%h(i,j,k) = h_new(i,j,k)
3417 enddo ; enddo ; enddo
3418
3419 deallocate(h_new, dzregrid, pcm_cell, h_old_u, h_new_u, h_old_v, h_new_v)
3420
3421 call cpu_clock_begin(id_clock_pass_init)
3422 call create_group_pass(tmp_pass_uv_t_s_h, cs%u, cs%v, g%Domain)
3423 if (use_temperature) then
3424 call create_group_pass(tmp_pass_uv_t_s_h, cs%tv%T, g%Domain)
3425 call create_group_pass(tmp_pass_uv_t_s_h, cs%tv%S, g%Domain)
3426 endif
3427 call create_group_pass(tmp_pass_uv_t_s_h, cs%h, g%Domain)
3428 call do_group_pass(tmp_pass_uv_t_s_h, g%Domain)
3429 call cpu_clock_end(id_clock_pass_init)
3430
3431 if (cs%debug) then
3432 call uvchksum("Post ALE adjust init cond [uv]", cs%u, cs%v, g%HI, haloshift=1, unscale=us%L_T_to_m_s)
3433 call hchksum(cs%h, "Post ALE adjust init cond h", g%HI, haloshift=2, unscale=gv%H_to_MKS)
3434 if (use_temperature) then
3435 call hchksum(cs%tv%T, "Post ALE adjust init cond T", g%HI, haloshift=2, unscale=us%C_to_degC)
3436 call hchksum(cs%tv%S, "Post ALE adjust init cond S", g%HI, haloshift=2, unscale=us%S_to_ppt)
3437 endif
3438 endif
3439 endif
3440 if ( cs%use_ALE_algorithm ) then
3441 call ale_set_extrap_boundaries (param_file, cs%ALE_CSp)
3442 call calltree_waypoint("returned from ALE_init() (initialize_MOM)")
3443 call ale_updateverticalgridtype( cs%ALE_CSp, gv )
3444 endif
3445 ! The basic state variables have now been fully initialized, so update their halos and
3446 ! calculate any derived thermodynmics quantities.
3447
3448 !--- set up group pass for u,v,T,S and h. pass_uv_T_S_h also is used in step_MOM
3449 call cpu_clock_begin(id_clock_pass_init)
3450 dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
3451 call create_group_pass(pass_uv_t_s_h, cs%u, cs%v, g%Domain, halo=dynamics_stencil)
3452 if (use_temperature) then
3453 call create_group_pass(pass_uv_t_s_h, cs%tv%T, g%Domain, halo=dynamics_stencil)
3454 call create_group_pass(pass_uv_t_s_h, cs%tv%S, g%Domain, halo=dynamics_stencil)
3455 endif
3456 call create_group_pass(pass_uv_t_s_h, cs%h, g%Domain, halo=dynamics_stencil)
3457
3458 call do_group_pass(pass_uv_t_s_h, g%Domain)
3459 if (associated(cs%tv%p_surf)) call pass_var(cs%tv%p_surf, g%Domain, halo=dynamics_stencil)
3460 call cpu_clock_end(id_clock_pass_init)
3461
3462 ! Update derived thermodynamic quantities.
3463 if (allocated(cs%tv%SpV_avg)) then
3464 call calc_derived_thermo(cs%tv, cs%h, g, gv, us, halo=dynamics_stencil, debug=cs%debug)
3465 endif
3466
3467
3468 diag => cs%diag
3469 ! Initialize the diag mediator.
3470 call diag_mediator_init(g, gv, us, gv%ke, param_file, diag, doc_file_dir=dirs%output_directory)
3471 if (associated(cs%OBC)) then
3472 call diag_mediator_set_obc_info(g, cs%OBC%segnum_u, cs%OBC%segnum_v, diag)
3473 endif
3474 if (present(diag_ptr)) diag_ptr => cs%diag
3475
3476 ! Initialize the diagnostics masks for native arrays.
3477 ! This step has to be done after call to MOM_initialize_state
3478 ! and before MOM_diagnostics_init
3479 call diag_masks_set(g, gv%ke, diag)
3480
3481 ! Set up pointers within diag mediator control structure,
3482 ! this needs to occur _after_ CS%h etc. have been allocated.
3483 call diag_set_state_ptrs(cs%h, cs%tv, diag)
3484
3485 ! This call sets up the diagnostic axes. These are needed,
3486 ! e.g. to generate the target grids below.
3487 call set_axes_info(g, gv, us, param_file, diag)
3488
3489 ! Whenever thickness/T/S changes let the diag manager know, target grids
3490 ! for vertical remapping may need to be regenerated. In non-Boussinesq mode,
3491 ! calc_derived_thermo needs to be called before diag_update_remap_grids.
3492 call diag_update_remap_grids(diag)
3493
3494 ! Setup the diagnostic grid storage types
3495 call diag_grid_storage_init(cs%diag_pre_sync, g, gv, diag)
3496 call diag_grid_storage_init(cs%diag_pre_dyn, g, gv, diag)
3497
3498 ! Calculate masks for diagnostics arrays in non-native coordinates
3499 ! This step has to be done after set_axes_info() because the axes needed
3500 ! to be configured, and after diag_update_remap_grids() because the grids
3501 ! must be defined.
3502 call set_masks_for_axes(g, diag)
3503
3504 ! Register the volume cell measure (must be one of first diagnostics)
3505 call register_cell_measure(g, cs%diag, time)
3506
3507 call cpu_clock_begin(id_clock_mom_init)
3508 ! Diagnose static fields AND associate areas/volumes with axes
3509 call write_static_fields(g, gv, us, cs%tv, cs%diag)
3510 call calltree_waypoint("static fields written (initialize_MOM)")
3511
3512 if (cs%use_ALE_algorithm) then
3513 call ale_writecoordinatefile( cs%ALE_CSp, gv, dirs%output_directory )
3514 call calltree_waypoint("ALE initialized (initialize_MOM)")
3515 elseif (write_geom_files) then
3516 call write_vertgrid_file(gv, us, param_file, dirs%output_directory)
3517 endif
3518 call cpu_clock_end(id_clock_mom_init)
3519
3520 if (cs%use_dbclient) call database_comms_init(param_file, cs%dbcomms_CS)
3521 cs%useMEKE = meke_init(time, g, gv, us, param_file, diag, cs%dbcomms_CS, cs%MEKE_CSp, cs%MEKE, &
3522 restart_csp, cs%MEKE_in_dynamics)
3523
3524 call varmix_init(time, g, gv, us, param_file, diag, cs%VarMix)
3525 call set_visc_init(time, g, gv, us, param_file, diag, cs%visc, cs%set_visc_CSp, restart_csp, cs%OBC)
3526 call thickness_diffuse_init(time, g, gv, us, param_file, diag, cs%CDp, cs%thickness_diffuse_CSp)
3527 if (cs%interface_filter) &
3528 call interface_filter_init(time, g, gv, us, param_file, diag, cs%CDp, cs%interface_filter_CSp)
3529
3530 new_sim = is_new_run(restart_csp)
3531 if (use_temperature) then
3532 cs%use_stochastic_EOS = mom_stoch_eos_init(time, g, gv, us, param_file, diag, cs%stoch_eos_CS, restart_csp)
3533 else
3534 cs%use_stochastic_EOS = .false.
3535 endif
3536
3537 if (cs%use_porbar) &
3538 call porous_barriers_init(time, gv, us, param_file, diag, cs%por_bar_CS)
3539
3540 if (cs%split) then
3541 allocate(eta(szi_(g),szj_(g)), source=0.0)
3542 if (cs%use_alt_split) then
3543 call initialize_dyn_split_rk2b(cs%u, cs%v, cs%h, cs%tv, cs%uh, cs%vh, eta, time, &
3544 g, gv, us, param_file, diag, cs%dyn_split_RK2b_CSp, cs%HA_CSp, restart_csp, &
3545 cs%dt, cs%ADp, cs%CDp, mom_internal_state, cs%VarMix, cs%MEKE, &
3546 cs%thickness_diffuse_CSp, cs%OBC, cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, &
3547 cs%visc, dirs, cs%ntrunc, cs%pbv, calc_dtbt=calc_dtbt, &
3548 cont_stencil=cs%cont_stencil, dyn_h_stencil=cs%dyn_h_stencil)
3549 else
3550 call initialize_dyn_split_rk2(cs%u, cs%v, cs%h, cs%tv, cs%uh, cs%vh, eta, time, &
3551 g, gv, us, param_file, diag, cs%dyn_split_RK2_CSp, cs%HA_CSp, restart_csp, &
3552 cs%dt, cs%ADp, cs%CDp, mom_internal_state, cs%VarMix, cs%MEKE, &
3553 cs%thickness_diffuse_CSp, cs%OBC, cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, &
3554 cs%visc, dirs, cs%ntrunc, cs%pbv, calc_dtbt=calc_dtbt, &
3555 cont_stencil=cs%cont_stencil, dyn_h_stencil=cs%dyn_h_stencil)
3556 endif
3557 ! A reset period no longer than dt is equivalent to recalculating every step.
3558 if (cs%dtbt_reset_period > 0.0 .and. cs%dtbt_reset_period <= cs%dt) &
3559 cs%dtbt_reset_period = 0.0
3560 if (cs%dtbt_reset_period > 0.0) then
3561 cs%dtbt_reset_interval = real_to_time(cs%dtbt_reset_period, unscale=us%T_to_s)
3562 if (calc_dtbt) then
3563 ! No restart DTBT (not found or new run) or DTBT_RESTART_BUG=True: set to the most recent
3564 ! multiple of the interval before or equal to current time, so the set_dtbt is called on
3565 ! the first step.
3566 cs%dtbt_reset_time = time_init + cs%dtbt_reset_interval * &
3567 ((time - time_init) / cs%dtbt_reset_interval)
3568 else
3569 ! Restart DTBT available: defer to the next multiple after current time so the first step
3570 ! uses the restart value unless a reset is naturally due.
3571 cs%dtbt_reset_time = time_init + cs%dtbt_reset_interval * &
3572 ((time - time_init) / cs%dtbt_reset_interval + 1)
3573 endif
3574 endif
3575 elseif (cs%use_RK2) then
3576 call initialize_dyn_unsplit_rk2(cs%u, cs%v, cs%h, cs%tv, time, g, gv, &
3577 us, param_file, diag, cs%dyn_unsplit_RK2_CSp, &
3578 cs%ADp, cs%CDp, mom_internal_state, cs%OBC, &
3579 cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, cs%visc, dirs, &
3580 cs%ntrunc, cont_stencil=cs%cont_stencil, dyn_h_stencil=cs%dyn_h_stencil)
3581 else
3582 call initialize_dyn_unsplit(cs%u, cs%v, cs%h, cs%tv, time, g, gv, &
3583 us, param_file, diag, cs%dyn_unsplit_CSp, &
3584 cs%ADp, cs%CDp, mom_internal_state, cs%OBC, &
3585 cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, cs%visc, dirs, &
3586 cs%ntrunc, cont_stencil=cs%cont_stencil, dyn_h_stencil=cs%dyn_h_stencil)
3587 endif
3588 cs%dyn_h_stencil = max(2, cs%dyn_h_stencil)
3589
3590 ! Set the next time to update OBC segment BGC tracer data
3591 if (associated(cs%OBC) .and. (cs%dt_obc_seg_period > 0.0)) then
3592 if (cs%dt_obc_seg_period > cs%dt_tr_adv) then
3593 if (new_sim) then
3594 call mom_error(warning, "DT_OBC_SEG_UPDATE_OBGC > DT_TRACER_ADVECT: this run will "//&
3595 "proceed normally, but any restart from it will fail with a FATAL error because "//&
3596 "BGC OBC external data is not saved in restart files. "//&
3597 "Set DT_OBC_SEG_UPDATE_OBGC = DT_TRACER_ADVECT or 0 until this is fixed.")
3598 else
3599 call mom_error(fatal, "DT_OBC_SEG_UPDATE_OBGC > DT_TRACER_ADVECT is not supported "//&
3600 "for restart runs: BGC OBC external data is not saved in restart files, so the "//&
3601 "first tracer advection step after restart may use incorrect boundary data. "//&
3602 "Set DT_OBC_SEG_UPDATE_OBGC = DT_TRACER_ADVECT or 0 until this is fixed.")
3603 endif
3604 endif
3605 cs%dt_obc_seg_interval = real_to_time(cs%dt_obc_seg_period, unscale=us%T_to_s)
3606 if (obc_bgc_time_ref_bug) then
3607 cs%dt_obc_seg_time = time + cs%dt_obc_seg_interval
3608 else
3609 ! Set to the next update point after current time, so that %t read from initialization is
3610 ! not overwritten. Note that even though this line by itself is correct, there are still
3611 ! issue with restart runs, as external data %t is not saved in restart files.
3612 cs%dt_obc_seg_time = time_init + cs%dt_obc_seg_interval * &
3613 ((time - time_init) / cs%dt_obc_seg_interval + 1)
3614 endif
3615 endif
3616
3617 call calltree_waypoint("dynamics initialized (initialize_MOM)")
3618
3619 cs%mixedlayer_restrat = mixedlayer_restrat_init(time, g, gv, us, param_file, diag, &
3620 cs%mixedlayer_restrat_CSp, restart_csp)
3621
3622 if (gv%Boussinesq .and. associated(cs%visc%h_ML)) then
3623 ! This is here to allow for a transition of restart files between model versions.
3624 call get_param(param_file, "MOM", "MLE_USE_PBL_MLD", mle_use_pbl_mld, &
3625 default=.false., do_not_log=.true.)
3626 if (mle_use_pbl_mld .and. .not.query_initialized(cs%visc%h_ML, "h_ML", restart_csp) .and. &
3627 associated(cs%visc%MLD)) then
3628 do j=js,je ; do i=is,ie ; cs%visc%h_ML(i,j) = gv%Z_to_H * cs%visc%MLD(i,j) ; enddo ; enddo
3629 endif
3630 endif
3631
3632 if (cs%mixedlayer_restrat) then
3633 if (.not.(bulkmixedlayer .or. cs%use_ALE_algorithm)) &
3634 call mom_error(fatal, "MOM: MIXEDLAYER_RESTRAT true requires a boundary layer scheme.")
3635 ! When DIABATIC_FIRST=False and using CS%visc%ML in mixedlayer_restrat we need to update after a restart
3636 if (.not. cs%diabatic_first .and. associated(cs%visc%MLD)) &
3637 call pass_var(cs%visc%MLD, g%domain, halo=1)
3638 if (.not. cs%diabatic_first .and. associated(cs%visc%h_ML)) &
3639 call pass_var(cs%visc%h_ML, g%domain, halo=1)
3640 endif
3641
3642 call mom_diagnostics_init(mom_internal_state, cs%ADp, cs%CDp, time, g, gv, us, &
3643 param_file, diag, cs%diagnostics_CSp, cs%tv)
3644 call diag_copy_diag_to_storage(cs%diag_pre_sync, cs%h, cs%diag)
3645
3646
3647 if (cs%adiabatic) then
3648 call adiabatic_driver_init(time, g, param_file, diag, cs%diabatic_CSp, &
3649 cs%tracer_flow_CSp)
3650 else
3651 call diabatic_driver_init(time, g, gv, us, param_file, cs%use_ALE_algorithm, diag, &
3652 cs%ADp, cs%CDp, cs%diabatic_CSp, cs%tracer_flow_CSp, &
3653 cs%sponge_CSp, cs%ALE_sponge_CSp, cs%oda_incupd_CSp, cs%int_tide_CSp)
3654 endif
3655
3656 cs%vertex_shear = kappa_shear_at_vertex(param_file)
3657
3658 ! GMM, the following is needed to get BLDs into the dynamics module
3659 if (cs%split .and. fpmix) then
3660 call init_dyn_split_rk2_diabatic(cs%diabatic_CSp, cs%dyn_split_RK2_CSp)
3661 endif
3662
3663 if (associated(cs%sponge_CSp)) &
3664 call init_sponge_diags(time, g, gv, us, diag, cs%sponge_CSp)
3665
3666 if (associated(cs%oda_incupd_CSp)) &
3667 call init_oda_incupd_diags(time, g, gv, diag, cs%oda_incupd_CSp, us)
3668
3669 call tracer_advect_init(time, g, us, param_file, diag, cs%tracer_adv_CSp)
3670 call tracer_hor_diff_init(time, g, gv, us, param_file, diag, cs%tv%eqn_of_state, cs%diabatic_CSp, &
3671 cs%tracer_diff_CSp)
3672
3673 call lock_tracer_registry(cs%tracer_Reg)
3674 call calltree_waypoint("tracer registry now locked (initialize_MOM)")
3675
3676 ! now register some diagnostics since the tracer registry is now locked
3677 call register_surface_diags(time, g, us, cs%sfc_IDs, cs%diag, cs%tv)
3678 call register_diags(time, g, gv, us, cs%IDs, cs%diag)
3679 call register_transport_diags(time, g, gv, us, cs%transport_IDs, cs%diag)
3680 call extract_diabatic_member(cs%diabatic_CSp, use_kpp=use_kpp)
3681 call register_tracer_diagnostics(cs%tracer_Reg, cs%h, time, diag, g, gv, us, &
3682 cs%use_ALE_algorithm, use_kpp)
3683 if (cs%use_ALE_algorithm) then
3684 call ale_register_diags(time, g, gv, us, diag, cs%ALE_CSp)
3685 endif
3686
3687 ! Do any necessary halo updates on any auxiliary variables that have been initialized.
3688 call cpu_clock_begin(id_clock_pass_init)
3689 if (associated(cs%visc%Kv_shear)) &
3690 call pass_var(cs%visc%Kv_shear, g%Domain, to_all+omit_corners, halo=1)
3691
3692 if (associated(cs%visc%Kv_slow)) &
3693 call pass_var(cs%visc%Kv_slow, g%Domain, to_all+omit_corners, halo=1)
3694 call cpu_clock_end(id_clock_pass_init)
3695
3696 ! This subroutine initializes any tracer packages.
3697 call tracer_flow_control_init(.not.new_sim, time, g, gv, us, cs%h, param_file, &
3698 cs%diag, cs%OBC, cs%tracer_flow_CSp, cs%sponge_CSp, &
3699 cs%ALE_sponge_CSp, cs%tv)
3700 if (present(tracer_flow_csp)) tracer_flow_csp => cs%tracer_flow_CSp
3701
3702 if (associated(cs%ALE_sponge_CSp)) &
3703 call init_ale_sponge_diags(time, g, diag, cs%ALE_sponge_CSp, us)
3704
3705 ! If running in offline tracer mode, initialize the necessary control structure and
3706 ! parameters
3707 if (present(offline_tracer_mode)) offline_tracer_mode=cs%offline_tracer_mode
3708
3709 if (cs%offline_tracer_mode) then
3710 ! Setup some initial parameterizations and also assign some of the subtypes
3711 call offline_transport_init(param_file, cs%offline_CSp, cs%diabatic_CSp, g, gv, us)
3712 call insert_offline_main( cs=cs%offline_CSp, ale_csp=cs%ALE_CSp, diabatic_csp=cs%diabatic_CSp, &
3713 diag=cs%diag, obc=cs%OBC, tracer_adv_csp=cs%tracer_adv_CSp, &
3714 tracer_flow_csp=cs%tracer_flow_CSp, tracer_reg=cs%tracer_Reg, &
3715 tv=cs%tv, x_before_y=(modulo(first_direction,2)==0), debug=cs%debug )
3716 call register_diags_offline_transport(time, cs%diag, cs%offline_CSp, gv, us)
3717 endif
3718
3719 if (associated(cs%OBC)) then
3720 ! At this point any information related to the tracer reservoirs has either been read from
3721 ! the restart file or has been specified in the segments. Initialize the tracer reservoir
3722 ! values from the segments if they have not been set via the restart file.
3723 call setup_obc_tracer_reservoirs(g, gv, cs%OBC, restart_csp)
3724 call setup_obc_thickness_reservoirs(g, gv, cs%OBC, restart_csp)
3725 call open_boundary_halo_update(g, cs%OBC)
3726 endif
3727
3728 call register_obsolete_diagnostics(param_file, cs%diag)
3729
3730 if (use_frazil) then
3731 if (.not.query_initialized(cs%tv%frazil, "frazil", restart_csp)) then
3732 cs%tv%frazil(:,:) = 0.0
3733 call set_initialized(cs%tv%frazil, "frazil", restart_csp)
3734 endif
3735 endif
3736
3737 if (cs%interp_p_surf) then
3738 cs%p_surf_prev_set = query_initialized(cs%p_surf_prev, "p_surf_prev", restart_csp)
3739
3740 if (cs%p_surf_prev_set) then
3741 call pass_var(cs%p_surf_prev, g%domain)
3742 endif
3743 endif
3744
3745 if (.not.query_initialized(cs%ave_ssh_ibc, "ave_ssh", restart_csp)) then
3746 if (cs%split) then
3747 call find_eta(cs%h, cs%tv, g, gv, us, cs%ave_ssh_ibc, eta, dzref=g%Z_ref)
3748 else
3749 call find_eta(cs%h, cs%tv, g, gv, us, cs%ave_ssh_ibc, dzref=g%Z_ref)
3750 endif
3751 call set_initialized(cs%ave_ssh_ibc, "ave_ssh", restart_csp)
3752 endif
3753 if (cs%split) deallocate(eta)
3754
3755 cs%nstep_tot = 0
3756 if (present(count_calls)) cs%count_calls = count_calls
3757 call mom_sum_output_init(g_in, gv, us, param_file, dirs%output_directory, &
3758 cs%ntrunc, time_init, cs%sum_output_CSp)
3759
3760 ! Flag whether to save initial conditions in finish_MOM_initialization() or not.
3761 cs%write_IC = save_ic .and. &
3762 .not.((dirs%input_filename(1:1) == 'r') .and. &
3763 (len_trim(dirs%input_filename) == 1))
3764
3765 if (cs%ensemble_ocean) then
3766 call init_oda(time, g, gv, us, cs%diag, cs%odaCS)
3767 endif
3768
3769 ! initialize stochastic physics
3770 call stochastics_init(cs%dt_therm, cs%G, cs%GV, cs%stoch_CS, param_file, diag, time)
3771
3772 call calltree_leave("initialize_MOM()")
3773 call cpu_clock_end(id_clock_init) ; call cpu_clock_end(id_clock_ocean)
3774
3775end subroutine initialize_mom
3776
3777!> Finishes initializing MOM and writes out the initial conditions.
3778subroutine finish_mom_initialization(Time, dirs, CS)
3779 type(time_type), intent(in) :: time !< model time, used in this routine
3780 type(directories), intent(in) :: dirs !< structure with directory paths
3781 type(mom_control_struct), intent(inout) :: cs !< MOM control structure
3782
3783 type(ocean_grid_type), pointer :: g => null() ! pointer to a structure containing
3784 ! metrics and related information
3785 type(verticalgrid_type), pointer :: gv => null() ! Pointer to the vertical grid structure
3786 type(unit_scale_type), pointer :: us => null() ! Pointer to a structure containing
3787 ! various unit conversion factors
3788 type(mom_restart_cs), pointer :: restart_csp_tmp => null()
3789 real, allocatable :: z_interface(:,:,:) ! Interface heights [Z ~> m]
3790
3791 call cpu_clock_begin(id_clock_init)
3792 call calltree_enter("finish_MOM_initialization()")
3793
3794 ! Pointers for convenience
3795 g => cs%G ; gv => cs%GV ; us => cs%US
3796
3797 if (cs%use_particles) then
3798 call particles_init(cs%particles, g, cs%Time, cs%dt_therm, cs%u, cs%v, cs%h)
3799 endif
3800
3801 ! Write initial conditions
3802 if (cs%write_IC) then
3803 allocate(restart_csp_tmp)
3804 restart_csp_tmp = cs%restart_CS
3805 call restart_registry_lock(restart_csp_tmp, unlocked=.true.)
3806 allocate(z_interface(szi_(g),szj_(g),szk_(gv)+1))
3807 call find_eta(cs%h, cs%tv, g, gv, us, z_interface, dzref=g%Z_ref)
3808 call register_restart_field(z_interface, "eta", .true., restart_csp_tmp, &
3809 "Interface heights", "meter", z_grid='i', conversion=us%Z_to_m)
3810 ! NOTE: write_ic=.true. routes routine to fms2 IO write_initial_conditions interface
3811 call save_restart(dirs%output_directory, time, cs%G_in, &
3812 restart_csp_tmp, filename=cs%IC_file, gv=gv, write_ic=.true.)
3813 deallocate(z_interface)
3814 deallocate(restart_csp_tmp)
3815 endif
3816
3817 call write_energy(cs%u, cs%v, cs%h, cs%tv, time, 0, g, gv, us, &
3818 cs%sum_output_CSp, cs%tracer_flow_CSp)
3819
3820 call calltree_leave("finish_MOM_initialization()")
3821 call cpu_clock_end(id_clock_init)
3822
3823end subroutine finish_mom_initialization
3824
3825!> Register certain diagnostics
3826subroutine register_diags(Time, G, GV, US, IDs, diag)
3827 type(time_type), intent(in) :: Time !< current model time
3828 type(ocean_grid_type), intent(in) :: G !< ocean grid structure
3829 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
3830 type(unit_scale_type), intent(inout) :: US !< A dimensional unit scaling type
3831 type(mom_diag_ids), intent(inout) :: IDs !< A structure with the diagnostic IDs.
3832 type(diag_ctrl), intent(inout) :: diag !< regulates diagnostic output
3833
3834 character(len=48) :: thickness_units
3835
3836 thickness_units = get_thickness_units(gv)
3837
3838 ! Diagnostics of the rapidly varying dynamic state
3839 ids%id_u = register_diag_field('ocean_model', 'u_dyn', diag%axesCuL, time, &
3840 'Zonal velocity after the dynamics update', 'm s-1', conversion=us%L_T_to_m_s)
3841 ids%id_v = register_diag_field('ocean_model', 'v_dyn', diag%axesCvL, time, &
3842 'Meridional velocity after the dynamics update', 'm s-1', conversion=us%L_T_to_m_s)
3843 ids%id_h = register_diag_field('ocean_model', 'h_dyn', diag%axesTL, time, &
3844 'Layer Thickness after the dynamics update', thickness_units, conversion=gv%H_to_MKS, &
3845 v_extensive=.true.)
3846 ids%id_ssh_inst = register_diag_field('ocean_model', 'SSH_inst', diag%axesT1, &
3847 time, 'Instantaneous Sea Surface Height', 'm', conversion=us%Z_to_m)
3848
3849end subroutine register_diags
3850
3851!> Set up CPU clock IDs for timing various subroutines.
3852subroutine mom_timing_init(CS)
3853 type(mom_control_struct), intent(in) :: CS !< control structure set up by initialize_MOM.
3854
3855 id_clock_dynamics = cpu_clock_id('Ocean dynamics', grain=clock_subcomponent)
3856 id_clock_thermo = cpu_clock_id('Ocean thermodynamics and tracers', grain=clock_subcomponent)
3857 id_clock_remap = cpu_clock_id('Ocean grid generation and remapping', grain=clock_subcomponent)
3858 id_clock_other = cpu_clock_id('Ocean Other', grain=clock_subcomponent)
3859 id_clock_mom_end = cpu_clock_id('Ocean MOM_end', grain=clock_subcomponent)
3860 id_clock_tracer = cpu_clock_id('(Ocean tracer advection)', grain=clock_module_driver)
3861 if (.not.cs%adiabatic) then
3862 id_clock_diabatic = cpu_clock_id('(Ocean diabatic driver)', grain=clock_module_driver)
3863 else
3864 id_clock_adiabatic = cpu_clock_id('(Ocean adiabatic driver)', grain=clock_module_driver)
3865 endif
3866
3867 id_clock_continuity = cpu_clock_id('(Ocean continuity equation *)', grain=clock_module)
3868 id_clock_bbl_visc = cpu_clock_id('(Ocean set BBL viscosity)', grain=clock_module)
3869 id_clock_pass = cpu_clock_id('(Ocean message passing *)', grain=clock_module)
3870 id_clock_mom_init = cpu_clock_id('(Ocean MOM_initialize_state)', grain=clock_module)
3871 id_clock_pass_init = cpu_clock_id('(Ocean init message passing *)', grain=clock_routine)
3872 if (cs%thickness_diffuse) &
3873 id_clock_thick_diff = cpu_clock_id('(Ocean thickness diffusion *)', grain=clock_module)
3874 if (cs%interface_filter) &
3875 id_clock_int_filter = cpu_clock_id('(Ocean interface height filter *)', grain=clock_module)
3876 !if (CS%mixedlayer_restrat) &
3877 id_clock_ml_restrat = cpu_clock_id('(Ocean mixed layer restrat)', grain=clock_module)
3878 id_clock_diagnostics = cpu_clock_id('(Ocean collective diagnostics)', grain=clock_module)
3879 id_clock_z_diag = cpu_clock_id('(Ocean Z-space diagnostics)', grain=clock_module)
3880 id_clock_ale = cpu_clock_id('(Ocean ALE)', grain=clock_module)
3881 if (cs%offline_tracer_mode) then
3882 id_clock_offline_tracer = cpu_clock_id('Ocean offline tracers', grain=clock_subcomponent)
3883 endif
3884 id_clock_stoch = cpu_clock_id('(Stochastic EOS)', grain=clock_module)
3885 id_clock_vart = cpu_clock_id('(SGS Temperature Variance)', grain=clock_module)
3886
3887 id_clock_save_restart = cpu_clock_id('(Ocean MOM save_restart)', grain=clock_module)
3888
3889end subroutine mom_timing_init
3890
3891!> Set the fields that are needed for bitwise identical restarting
3892!! the time stepping scheme. In addition to those specified here
3893!! directly, there may be fields related to the forcing or to the
3894!! barotropic solver that are needed; these are specified in sub-
3895!! routines that are called from this one.
3896!!
3897!! This routine should be altered if there are any changes to the
3898!! time stepping scheme. The CHECK_RESTART facility may be used to
3899!! confirm that all needed restart fields have been included.
3900subroutine set_restart_fields(GV, US, param_file, CS, restart_CSp)
3901 type(verticalgrid_type), intent(inout) :: GV !< ocean vertical grid structure
3902 type(unit_scale_type), intent(inout) :: US !< A dimensional unit scaling type
3903 type(param_file_type), intent(in) :: param_file !< opened file for parsing to get parameters
3904 type(mom_control_struct), intent(in) :: CS !< control structure set up by initialize_MOM
3905 type(mom_restart_cs), pointer :: restart_CSp !< pointer to the restart control
3906 !! structure that will be used for MOM.
3907 ! Local variables
3908 logical :: use_ice_shelf ! Needed to determine whether to add CS%Hml to restarts
3909 character(len=48) :: thickness_units, flux_units
3910 type(vardesc) :: u_desc, v_desc
3911
3912 thickness_units = get_thickness_units(gv)
3913 flux_units = get_flux_units(gv)
3914
3915 if (associated(cs%tv%T)) &
3916 call register_restart_field(cs%tv%T, "Temp", .true., restart_csp, &
3917 "Potential Temperature", "degC", conversion=us%C_to_degC)
3918 if (associated(cs%tv%S)) &
3919 call register_restart_field(cs%tv%S, "Salt", .true., restart_csp, &
3920 "Salinity", "PPT", conversion=us%S_to_ppt)
3921
3922 call register_restart_field(cs%h, "h", .true., restart_csp, &
3923 "Layer Thickness", thickness_units, conversion=gv%H_to_MKS)
3924
3925 u_desc = var_desc("u", "m s-1", "Zonal velocity", hor_grid='Cu')
3926 v_desc = var_desc("v", "m s-1", "Meridional velocity", hor_grid='Cv')
3927 call register_restart_pair(cs%u, cs%v, u_desc, v_desc, .true., restart_csp, conversion=us%L_T_to_m_s)
3928
3929 if (associated(cs%tv%frazil)) &
3930 call register_restart_field(cs%tv%frazil, "frazil", .false., restart_csp, &
3931 "Frazil heat flux into ocean", &
3932 "J m-2", conversion=us%Q_to_J_kg*us%RZ_to_kg_m2)
3933
3934 if (cs%interp_p_surf) then
3935 call register_restart_field(cs%p_surf_prev, "p_surf_prev", .false., restart_csp, &
3936 "Previous ocean surface pressure", "Pa", conversion=us%RL2_T2_to_Pa)
3937 endif
3938
3939 if (associated(cs%tv%p_surf)) &
3940 call register_restart_field(cs%tv%p_surf, "p_surf_EOS", .false., restart_csp, &
3941 "Ocean surface pressure used in EoS", "Pa", conversion=us%RL2_T2_to_Pa)
3942
3943 call register_restart_field(cs%ave_ssh_ibc, "ave_ssh", .false., restart_csp, &
3944 "Time average sea surface height", "meter", conversion=us%Z_to_m)
3945
3946 ! hML is needed when using the ice shelf module
3947 call get_param(param_file, '', "ICE_SHELF", use_ice_shelf, default=.false., &
3948 do_not_log=.true.)
3949 if (use_ice_shelf .and. associated(cs%Hml)) then
3950 call register_restart_field(cs%Hml, "hML", .false., restart_csp, &
3951 "Mixed layer thickness", "m", conversion=us%Z_to_m)
3952 endif
3953
3954 ! Register scalar unit conversion factors.
3955 call register_restart_field(cs%first_dir_restart, "First_direction", .false., restart_csp, &
3956 "Indicator of the first direction in split calculations.", "nondim")
3957
3958end subroutine set_restart_fields
3959
3960!> Apply a correction to the sea surface height to compensate
3961!! for the atmospheric pressure (the inverse barometer).
3962subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS)
3963 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
3964 type(ocean_grid_type), intent(in) :: G !< ocean grid structure
3965 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
3966 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3967 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: ssh !< time mean surface height [Z ~> m]
3968 real, dimension(:,:), pointer :: p_atm !< Ocean surface pressure [R L2 T-2 ~> Pa]
3969 logical, intent(in) :: use_EOS !< If true, calculate the density for
3970 !! the SSH correction using the equation of state.
3971
3972 real :: Rho_conv(SZI_(G)) ! The density used to convert surface pressure to
3973 ! a corrected effective SSH [R ~> kg m-3].
3974 real :: IgR0 ! The SSH conversion factor from R L2 T-2 to Z [Z T2 R-1 L-2 ~> m Pa-1].
3975 logical :: calc_rho
3976 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
3977 integer :: i, j, is, ie, js, je
3978
3979 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3980 eosdom(:) = eos_domain(g%HI)
3981 if (associated(p_atm)) then
3982 calc_rho = use_eos .and. associated(tv%eqn_of_state)
3983 ! Correct the output sea surface height for the contribution from the ice pressure.
3984 do j=js,je
3985 if (calc_rho) then
3986 call calculate_density(tv%T(:,j,1), tv%S(:,j,1), 0.5*p_atm(:,j), rho_conv, &
3987 tv%eqn_of_state, eosdom)
3988 do i=is,ie
3989 igr0 = 1.0 / (rho_conv(i) * gv%g_Earth)
3990 ssh(i,j) = ssh(i,j) + p_atm(i,j) * igr0
3991 enddo
3992 else
3993 igr0 = 1.0 / (gv%Rho0 * gv%g_Earth)
3994 do i=is,ie
3995 ssh(i,j) = ssh(i,j) + p_atm(i,j) * igr0
3996 enddo
3997 endif
3998 enddo
3999 endif
4000
4001end subroutine adjust_ssh_for_p_atm
4002
4003!> Set the surface (return) properties of the ocean model by
4004!! setting the appropriate fields in sfc_state. Unused fields
4005!! are set to NULL or are unallocated.
4006subroutine extract_surface_state(CS, sfc_state_in)
4007 type(mom_control_struct), intent(inout), target :: cs !< Master MOM control structure
4008 type(surface), target, intent(inout) :: sfc_state_in !< transparent ocean surface state
4009 !! structure shared with the calling routine
4010 !! data in this structure is intent out.
4011
4012 ! Local variables
4013 real :: hu, hv ! Thicknesses interpolated to velocity points [H ~> m or kg m-2]
4014 type(ocean_grid_type), pointer :: g => null() !< pointer to a structure containing
4015 !! metrics and related information
4016 type(ocean_grid_type), pointer :: g_in => null() !< Input grid metric
4017 type(verticalgrid_type), pointer :: gv => null() !< structure containing vertical grid info
4018 type(unit_scale_type), pointer :: us => null() !< structure containing various unit conversion factors
4019 type(surface), pointer :: sfc_state => null() ! surface state on the model grid
4020 real, dimension(:,:,:), pointer :: h => null() !< h : layer thickness [H ~> m or kg m-2]
4021 real :: depth(szi_(cs%g)) !< Distance from the surface in depth units [Z ~> m] or [H ~> m or kg m-2]
4022 real :: depth_ml !< Depth over which to average to determine mixed
4023 !! layer properties [Z ~> m] or [H ~> m or kg m-2]
4024 real :: dh !< Thickness of a layer within the mixed layer [Z ~> m] or [H ~> m or kg m-2]
4025 real :: mass !< Mass per unit area of a layer [R Z ~> kg m-2]
4026 real :: i_depth !< The inverse of depth [Z-1 ~> m-1] or [H-1 ~> m-1 or m2 kg-1]
4027 real :: missing_depth !< The portion of depth_ml that can not be found in a column [H ~> m or kg m-2]
4028 real :: h_rescale !< A conversion factor from thickness units to the units used in the
4029 !! calculation of properties of the uppermost ocean [nondim] or [Z H-1 ~> 1 or m3 kg-1]
4030 ! After the ANSWERS_2018 flag has been obsoleted, H_rescale will be 1.
4031 real :: t_freeze(szi_(cs%g)) !< freezing temperature [C ~> degC]
4032 real :: pres(szi_(cs%g)) !< Pressure to use for the freezing temperature calculation [R L2 T-2 ~> Pa]
4033 real :: delt(szi_(cs%g)) !< Depth integral of T-T_freeze [H C ~> m degC or degC kg m-2]
4034 logical :: use_temperature !< If true, temperature and salinity are used as state variables.
4035 integer :: i, j, k, is, ie, js, je, nz, numberoferrors, ig, jg
4036 integer :: isd, ied, jsd, jed
4037 integer :: iscb, iecb, jscb, jecb, isdb, iedb, jsdb, jedb
4038 logical :: localerror
4039 logical :: use_iceshelves
4040 character(240) :: msg
4041 integer :: turns ! Number of quarter turns
4042
4043 call calltree_enter("extract_surface_state(), MOM.F90")
4044 g => cs%G ; g_in => cs%G_in ; gv => cs%GV ; us => cs%US
4045 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
4046 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
4047 iscb = g%iscB ; iecb = g%iecB ; jscb = g%jscB ; jecb = g%jecB
4048 isdb = g%isdB ; iedb = g%iedB ; jsdb = g%jsdB ; jedb = g%jedB
4049 h => cs%h
4050
4051 use_temperature = associated(cs%tv%T)
4052
4053 use_iceshelves=.false.
4054 if (associated(cs%frac_shelf_h)) use_iceshelves = .true.
4055
4056 turns = 0
4057 if (cs%rotate_index) &
4058 turns = g%HI%turns
4059
4060 if (.not.sfc_state_in%arrays_allocated) then
4061 ! Consider using a run-time flag to determine whether to do the vertical
4062 ! integrals, since the 3-d sums are not negligible in cost.
4063 call allocate_surface_state(sfc_state_in, g_in, use_temperature, &
4064 do_integrals=.true., omit_frazil=.not.associated(cs%tv%frazil),&
4065 use_iceshelves=use_iceshelves)
4066 endif
4067
4068 if (cs%rotate_index) then
4069 allocate(sfc_state)
4070 call allocate_surface_state(sfc_state, g, use_temperature, &
4071 do_integrals=.true., omit_frazil=.not.associated(cs%tv%frazil),&
4072 use_iceshelves=use_iceshelves, sfc_state_in=sfc_state_in, turns=turns)
4073 else
4074 sfc_state => sfc_state_in
4075 endif
4076
4077 sfc_state%T_is_conT = cs%tv%T_is_conT
4078 sfc_state%S_is_absS = cs%tv%S_is_absS
4079
4080 do j=js,je ; do i=is,ie
4081 sfc_state%sea_lev(i,j) = cs%ave_ssh_ibc(i,j)
4082 enddo ; enddo
4083
4084 if (allocated(sfc_state%frazil) .and. associated(cs%tv%frazil)) then ; do j=js,je ; do i=is,ie
4085 sfc_state%frazil(i,j) = cs%tv%frazil(i,j)
4086 enddo ; enddo ; endif
4087
4088 ! copy Hml into sfc_state, so that caps can access it
4089 do j=js,je ; do i=is,ie
4090 sfc_state%Hml(i,j) = cs%Hml(i,j)
4091 enddo ; enddo
4092
4093 if (cs%Hmix < 0.0) then ! A bulk mixed layer is in use, so layer 1 has the properties
4094 if (use_temperature) then ; do j=js,je ; do i=is,ie
4095 sfc_state%SST(i,j) = cs%tv%T(i,j,1)
4096 sfc_state%SSS(i,j) = cs%tv%S(i,j,1)
4097 enddo ; enddo ; endif
4098 do j=js,je ; do i=is-1,ie
4099 sfc_state%u(i,j) = cs%u(i,j,1)
4100 enddo ; enddo
4101 do j=js-1,je ; do i=is,ie
4102 sfc_state%v(i,j) = cs%v(i,j,1)
4103 enddo ; enddo
4104
4105 else ! (CS%Hmix >= 0.0)
4106 h_rescale = 1.0
4107 depth_ml = cs%Hmix
4108 if (cs%answer_date < 20190101) then
4109 h_rescale = gv%H_to_Z
4110 depth_ml = gv%H_to_Z*cs%Hmix
4111 endif
4112 ! Determine the mean tracer properties of the uppermost depth_ml fluid.
4113
4114 !$OMP parallel do default(shared) private(depth,dh)
4115 do j=js,je
4116 do i=is,ie
4117 depth(i) = 0.0
4118 if (use_temperature) then
4119 sfc_state%SST(i,j) = 0.0 ; sfc_state%SSS(i,j) = 0.0
4120 else
4121 sfc_state%sfc_density(i,j) = 0.0
4122 endif
4123 enddo
4124
4125 do k=1,nz ; do i=is,ie
4126 if (depth(i) + h(i,j,k)*h_rescale < depth_ml) then
4127 dh = h(i,j,k)*h_rescale
4128 elseif (depth(i) < depth_ml) then
4129 dh = depth_ml - depth(i)
4130 else
4131 dh = 0.0
4132 endif
4133 if (use_temperature) then
4134 sfc_state%SST(i,j) = sfc_state%SST(i,j) + dh * cs%tv%T(i,j,k)
4135 sfc_state%SSS(i,j) = sfc_state%SSS(i,j) + dh * cs%tv%S(i,j,k)
4136 else
4137 sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) + dh * gv%Rlay(k)
4138 endif
4139 depth(i) = depth(i) + dh
4140 enddo ; enddo
4141 ! Calculate the average properties of the mixed layer depth.
4142 do i=is,ie
4143 if (cs%answer_date < 20190101) then
4144 if (depth(i) < gv%H_subroundoff*h_rescale) &
4145 depth(i) = gv%H_subroundoff*h_rescale
4146 if (use_temperature) then
4147 sfc_state%SST(i,j) = sfc_state%SST(i,j) / depth(i)
4148 sfc_state%SSS(i,j) = sfc_state%SSS(i,j) / depth(i)
4149 else
4150 sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) / depth(i)
4151 endif
4152 else
4153 if (depth(i) < gv%H_subroundoff*h_rescale) then
4154 i_depth = 1.0 / (gv%H_subroundoff*h_rescale)
4155 missing_depth = gv%H_subroundoff*h_rescale - depth(i)
4156 if (use_temperature) then
4157 sfc_state%SST(i,j) = (sfc_state%SST(i,j) + missing_depth*cs%tv%T(i,j,1)) * i_depth
4158 sfc_state%SSS(i,j) = (sfc_state%SSS(i,j) + missing_depth*cs%tv%S(i,j,1)) * i_depth
4159 else
4160 sfc_state%sfc_density(i,j) = (sfc_state%sfc_density(i,j) + &
4161 missing_depth*gv%Rlay(1)) * i_depth
4162 endif
4163 else
4164 i_depth = 1.0 / depth(i)
4165 if (use_temperature) then
4166 sfc_state%SST(i,j) = sfc_state%SST(i,j) * i_depth
4167 sfc_state%SSS(i,j) = sfc_state%SSS(i,j) * i_depth
4168 else
4169 sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) * i_depth
4170 endif
4171 endif
4172 endif
4173 enddo
4174 enddo ! end of j loop
4175
4176! Determine the mean velocities in the uppermost depth_ml fluid.
4177 ! NOTE: Velocity loops start on `[ij]s-1` in order to update halo values
4178 ! required by the speed diagnostic on the non-symmetric grid.
4179 ! This assumes that u and v halos have already been updated.
4180 if (cs%Hmix_UV>0.) then
4181 depth_ml = cs%Hmix_UV
4182 if (cs%answer_date < 20190101) depth_ml = gv%H_to_Z*cs%Hmix_UV
4183 !$OMP parallel do default(shared) private(depth,dh,hv)
4184 do j=js-1,je
4185 do i=is,ie
4186 depth(i) = 0.0
4187 sfc_state%v(i,j) = 0.0
4188 enddo
4189 do k=1,nz ; do i=is,ie
4190 hv = 0.5 * (h(i,j,k) + h(i,j+1,k)) * h_rescale
4191 if (depth(i) + hv < depth_ml) then
4192 dh = hv
4193 elseif (depth(i) < depth_ml) then
4194 dh = depth_ml - depth(i)
4195 else
4196 dh = 0.0
4197 endif
4198 sfc_state%v(i,j) = sfc_state%v(i,j) + dh * cs%v(i,j,k)
4199 depth(i) = depth(i) + dh
4200 enddo ; enddo
4201 ! Calculate the average properties of the mixed layer depth.
4202 do i=is,ie
4203 sfc_state%v(i,j) = sfc_state%v(i,j) / max(depth(i), gv%H_subroundoff*h_rescale)
4204 enddo
4205 enddo ! end of j loop
4206
4207 !$OMP parallel do default(shared) private(depth,dh,hu)
4208 do j=js,je
4209 do i=is-1,ie
4210 depth(i) = 0.0
4211 sfc_state%u(i,j) = 0.0
4212 enddo
4213 do k=1,nz ; do i=is-1,ie
4214 hu = 0.5 * (h(i,j,k) + h(i+1,j,k)) * h_rescale
4215 if (depth(i) + hu < depth_ml) then
4216 dh = hu
4217 elseif (depth(i) < depth_ml) then
4218 dh = depth_ml - depth(i)
4219 else
4220 dh = 0.0
4221 endif
4222 sfc_state%u(i,j) = sfc_state%u(i,j) + dh * cs%u(i,j,k)
4223 depth(i) = depth(i) + dh
4224 enddo ; enddo
4225 ! Calculate the average properties of the mixed layer depth.
4226 do i=is-1,ie
4227 sfc_state%u(i,j) = sfc_state%u(i,j) / max(depth(i), gv%H_subroundoff*h_rescale)
4228 enddo
4229 enddo ! end of j loop
4230 else ! Hmix_UV<=0.
4231 do j=js,je ; do i=is-1,ie
4232 sfc_state%u(i,j) = cs%u(i,j,1)
4233 enddo ; enddo
4234 do j=js-1,je ; do i=is,ie
4235 sfc_state%v(i,j) = cs%v(i,j,1)
4236 enddo ; enddo
4237 endif
4238 endif ! (CS%Hmix >= 0.0)
4239
4240
4241 if (allocated(sfc_state%melt_potential)) then
4242 !$OMP parallel do default(shared) private(depth_ml, dh, T_freeze, depth, pres, delT)
4243 do j=js,je
4244 do i=is,ie
4245 depth(i) = 0.0
4246 delt(i) = 0.0
4247 pres(i) = 0.0
4248 ! Here it is assumed that p=0 is OK, since HFrz ~ 10 to 20m, but under ice-shelves this
4249 ! can be a very bad assumption. ###To fix this, uncomment the following...
4250 ! pres(i) = p_surface(i) + 0.5*(GV%g_Earth*GV%H_to_RZ)*h(i,j,1)
4251 enddo
4252
4253 do k=1,nz
4254 call calculate_tfreeze(cs%tv%S(is:ie,j,k), pres(is:ie), t_freeze(is:ie), cs%tv%eqn_of_state)
4255 do i=is,ie
4256 depth_ml = min(cs%HFrz, cs%visc%h_ML(i,j))
4257 if (depth(i) + h(i,j,k) < depth_ml) then
4258 dh = h(i,j,k)
4259 elseif (depth(i) < depth_ml) then
4260 dh = depth_ml - depth(i)
4261 else
4262 dh = 0.0
4263 endif
4264
4265 depth(i) = depth(i) + dh
4266 delt(i) = delt(i) + dh * (cs%tv%T(i,j,k) - t_freeze(i))
4267 enddo
4268 ! If there is a pressure-dependent freezing point calculation uncomment the following.
4269 ! if (k<nz) then ; do i=is,ie
4270 ! pres(i) = pres(i) + 0.5*(GV%g_Earth*GV%H_to_RZ) * (h(i,j,k) + h(i,j,k+1))
4271 ! enddo ; endif
4272 enddo
4273
4274 do i=is,ie
4275 ! set melt_potential to zero to avoid passing previous values
4276 sfc_state%melt_potential(i,j) = 0.0
4277
4278 if (g%mask2dT(i,j)>0.) then
4279 ! instantaneous melt_potential [Q R Z ~> J m-2]
4280 sfc_state%melt_potential(i,j) = cs%tv%C_p * gv%H_to_RZ * delt(i)
4281 endif
4282 enddo
4283 enddo ! end of j loop
4284 endif ! melt_potential
4285
4286 if (allocated(sfc_state%taux_shelf) .and. allocated(cs%visc%taux_shelf)) then
4287 !$OMP parallel do default(shared)
4288 do j=js,je ; do i=is-1,ie
4289 sfc_state%taux_shelf(i,j) = cs%visc%taux_shelf(i,j)
4290 enddo ; enddo
4291 endif
4292 if (allocated(sfc_state%tauy_shelf) .and. allocated(cs%visc%tauy_shelf)) then
4293 !$OMP parallel do default(shared)
4294 do j=js-1,je ; do i=is,ie
4295 sfc_state%tauy_shelf(i,j) = cs%visc%tauy_shelf(i,j)
4296 enddo ; enddo
4297 endif
4298
4299 if (allocated(sfc_state%ocean_mass) .and. allocated(sfc_state%ocean_heat) .and. &
4300 allocated(sfc_state%ocean_salt)) then
4301 !$OMP parallel do default(shared)
4302 do j=js,je ; do i=is,ie
4303 sfc_state%ocean_mass(i,j) = 0.0
4304 sfc_state%ocean_heat(i,j) = 0.0 ; sfc_state%ocean_salt(i,j) = 0.0
4305 enddo ; enddo
4306 !$OMP parallel do default(shared) private(mass)
4307 do j=js,je ; do k=1,nz ; do i=is,ie
4308 mass = gv%H_to_RZ*h(i,j,k)
4309 sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + mass
4310 sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass * cs%tv%T(i,j,k)
4311 sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + mass * (1.0e-3*cs%tv%S(i,j,k))
4312 enddo ; enddo ; enddo
4313 else
4314 if (allocated(sfc_state%ocean_mass)) then
4315 !$OMP parallel do default(shared)
4316 do j=js,je ; do i=is,ie ; sfc_state%ocean_mass(i,j) = 0.0 ; enddo ; enddo
4317 !$OMP parallel do default(shared)
4318 do j=js,je ; do k=1,nz ; do i=is,ie
4319 sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + gv%H_to_RZ*h(i,j,k)
4320 enddo ; enddo ; enddo
4321 endif
4322 if (allocated(sfc_state%ocean_heat)) then
4323 !$OMP parallel do default(shared)
4324 do j=js,je ; do i=is,ie ; sfc_state%ocean_heat(i,j) = 0.0 ; enddo ; enddo
4325 !$OMP parallel do default(shared) private(mass)
4326 do j=js,je ; do k=1,nz ; do i=is,ie
4327 mass = gv%H_to_RZ*h(i,j,k)
4328 sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass * cs%tv%T(i,j,k)
4329 enddo ; enddo ; enddo
4330 endif
4331 if (allocated(sfc_state%ocean_salt)) then
4332 !$OMP parallel do default(shared)
4333 do j=js,je ; do i=is,ie ; sfc_state%ocean_salt(i,j) = 0.0 ; enddo ; enddo
4334 !$OMP parallel do default(shared) private(mass)
4335 do j=js,je ; do k=1,nz ; do i=is,ie
4336 mass = gv%H_to_RZ*h(i,j,k)
4337 sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + mass * (1.0e-3*cs%tv%S(i,j,k))
4338 enddo ; enddo ; enddo
4339 endif
4340 endif
4341
4342 if (associated(cs%tracer_flow_CSp)) then
4343 call call_tracer_surface_state(sfc_state, h, g, gv, us, cs%tracer_flow_CSp)
4344 endif
4345
4346 if (cs%check_bad_sfc_vals) then
4347 numberoferrors=0 ! count number of errors
4348 do j=js,je ; do i=is,ie
4349 if (g%mask2dT(i,j)>0.) then
4350 localerror = sfc_state%sea_lev(i,j) < -g%bathyT(i,j) - g%Z_ref &
4351 .or. sfc_state%sea_lev(i,j) >= cs%bad_val_ssh_max + (g%meanSL(i,j) - g%Z_ref) &
4352 .or. sfc_state%sea_lev(i,j) <= -cs%bad_val_ssh_max + (g%meanSL(i,j) - g%Z_ref) &
4353 .or. sfc_state%sea_lev(i,j) + g%bathyT(i,j) + g%Z_ref < cs%bad_val_col_thick
4354 if (use_temperature) localerror = localerror &
4355 .or. sfc_state%SSS(i,j)<0. &
4356 .or. sfc_state%SSS(i,j)>=cs%bad_val_sss_max &
4357 .or. sfc_state%SST(i,j)< cs%bad_val_sst_min &
4358 .or. sfc_state%SST(i,j)>=cs%bad_val_sst_max
4359 if (localerror) then
4360 numberoferrors=numberoferrors+1
4361 if (numberoferrors<9) then ! Only report details for the first few errors
4362 ig = i + g%HI%idg_offset ! Global i-index
4363 jg = j + g%HI%jdg_offset ! Global j-index
4364 if (use_temperature) then
4365 write(msg(1:240),'(2(a,I0,1x),4(a,f8.3,1x),8(a,es11.4,1x))') &
4366 'Extreme surface sfc_state detected: i=',ig,'j=',jg, &
4367 'lon=',g%geoLonT(i,j), 'lat=',g%geoLatT(i,j), &
4368 'x=',g%gridLonT(ig), 'y=',g%gridLatT(jg), &
4369 'D=',us%Z_to_m*(g%bathyT(i,j)+g%Z_ref), 'SSH=',us%Z_to_m*sfc_state%sea_lev(i,j), &
4370 'SST=',us%C_to_degC*sfc_state%SST(i,j), 'SSS=',us%S_to_ppt*sfc_state%SSS(i,j), &
4371 'U-=',us%L_T_to_m_s*sfc_state%u(i-1,j), 'U+=',us%L_T_to_m_s*sfc_state%u(i,j), &
4372 'V-=',us%L_T_to_m_s*sfc_state%v(i,j-1), 'V+=',us%L_T_to_m_s*sfc_state%v(i,j)
4373 else
4374 write(msg(1:240),'(2(a,I0,1x),4(a,f8.3,1x),6(a,es11.4))') &
4375 'Extreme surface sfc_state detected: i=',ig,'j=',jg, &
4376 'lon=',g%geoLonT(i,j), 'lat=',g%geoLatT(i,j), &
4377 'x=',g%gridLonT(ig), 'y=',g%gridLatT(jg), &
4378 'D=',us%Z_to_m*(g%bathyT(i,j)+g%Z_ref), 'SSH=',us%Z_to_m*sfc_state%sea_lev(i,j), &
4379 'U-=',us%L_T_to_m_s*sfc_state%u(i-1,j), 'U+=',us%L_T_to_m_s*sfc_state%u(i,j), &
4380 'V-=',us%L_T_to_m_s*sfc_state%v(i,j-1), 'V+=',us%L_T_to_m_s*sfc_state%v(i,j)
4381 endif
4382 call mom_error(warning, trim(msg), all_print=.true.)
4383 elseif (numberoferrors==9) then ! Indicate once that there are more errors
4384 call mom_error(warning, 'There were more unreported extreme events!', all_print=.true.)
4385 endif ! numberOfErrors
4386 endif ! localError
4387 endif ! mask2dT
4388 enddo ; enddo
4389 call sum_across_pes(numberoferrors)
4390 if (numberoferrors>0) then
4391 write(msg(1:240),'(a,i0,a)') 'There were a total of ',numberoferrors, &
4392 ' locations detected with extreme surface values!'
4393 call mom_error(fatal, trim(msg))
4394 endif
4395 endif
4396
4397 if (cs%debug) call mom_surface_chksum("Post extract_sfc", sfc_state, g, us, haloshift=0, symmetric=.true.)
4398
4399 ! Rotate sfc_state back onto the input grid, sfc_state_in
4400 if (cs%rotate_index) then
4401 call rotate_surface_state(sfc_state, sfc_state_in, g_in, -turns)
4402 call deallocate_surface_state(sfc_state)
4403 endif
4404
4405 call calltree_leave("extract_surface_sfc_state()")
4406end subroutine extract_surface_state
4407
4408!> Rotate initialization fields from input to rotated arrays.
4409subroutine rotate_initial_state(u_in, v_in, h_in, T_in, S_in, &
4410 use_temperature, turns, u, v, h, T, S)
4411 real, dimension(:,:,:), intent(in) :: u_in !< Zonal velocity on the initial grid [L T-1 ~> m s-1]
4412 real, dimension(:,:,:), intent(in) :: v_in !< Meridional velocity on the initial grid [L T-1 ~> m s-1]
4413 real, dimension(:,:,:), intent(in) :: h_in !< Layer thickness on the initial grid [H ~> m or kg m-2]
4414 real, dimension(:,:,:), intent(in) :: T_in !< Temperature on the initial grid [C ~> degC]
4415 real, dimension(:,:,:), intent(in) :: S_in !< Salinity on the initial grid [S ~> ppt]
4416 logical, intent(in) :: use_temperature !< If true, temperature and salinity are active
4417 integer, intent(in) :: turns !< The number quarter-turns to apply
4418 real, dimension(:,:,:), intent(out) :: u !< Zonal velocity on the rotated grid [L T-1 ~> m s-1]
4419 real, dimension(:,:,:), intent(out) :: v !< Meridional velocity on the rotated grid [L T-1 ~> m s-1]
4420 real, dimension(:,:,:), intent(out) :: h !< Layer thickness on the rotated grid [H ~> m or kg m-2]
4421 real, dimension(:,:,:), intent(out) :: T !< Temperature on the rotated grid [C ~> degC]
4422 real, dimension(:,:,:), intent(out) :: S !< Salinity on the rotated grid [S ~> ppt]
4423
4424 call rotate_vector(u_in, v_in, turns, u, v)
4425 call rotate_array(h_in, turns, h)
4426 if (use_temperature) then
4427 call rotate_array(t_in, turns, t)
4428 call rotate_array(s_in, turns, s)
4429 endif
4430end subroutine rotate_initial_state
4431
4432!> Return true if all phases of step_MOM are at the same point in time.
4433function mom_state_is_synchronized(CS, adv_dyn) result(in_synch)
4434 type(mom_control_struct), intent(inout) :: cs !< MOM control structure
4435 logical, optional, intent(in) :: adv_dyn !< If present and true, only check
4436 !! whether the advection is up-to-date with
4437 !! the dynamics.
4438 logical :: in_synch !< True if all phases of the update are synchronized.
4439
4440 logical :: adv_only
4441
4442 adv_only = .false. ; if (present(adv_dyn)) adv_only = adv_dyn
4443
4444 if (adv_only) then
4445 in_synch = (cs%t_dyn_rel_adv == 0.0)
4446 else
4447 in_synch = ((cs%t_dyn_rel_adv == 0.0) .and. (cs%t_dyn_rel_thermo == 0.0))
4448 endif
4449
4450end function mom_state_is_synchronized
4451
4452!> This subroutine offers access to values or pointers to other types from within
4453!! the MOM_control_struct, allowing the MOM_control_struct to be opaque.
4454subroutine get_mom_state_elements(CS, G, GV, US, C_p, C_p_scaled, use_temp)
4455 type(mom_control_struct), intent(inout), target :: cs !< MOM control structure
4456 type(ocean_grid_type), optional, pointer :: g !< structure containing metrics and grid info
4457 type(verticalgrid_type), optional, pointer :: gv !< structure containing vertical grid info
4458 type(unit_scale_type), optional, pointer :: us !< A dimensional unit scaling type
4459 real, optional, intent(out) :: c_p !< The heat capacity [J kg degC-1]
4460 real, optional, intent(out) :: c_p_scaled !< The heat capacity in scaled
4461 !! units [Q C-1 ~> J kg-1 degC-1]
4462 logical, optional, intent(out) :: use_temp !< True if temperature is a state variable
4463
4464 if (present(g)) g => cs%G_in
4465 if (present(gv)) gv => cs%GV
4466 if (present(us)) us => cs%US
4467 if (present(c_p)) c_p = cs%US%Q_to_J_kg*us%degC_to_C * cs%tv%C_p
4468 if (present(c_p_scaled)) c_p_scaled = cs%tv%C_p
4469 if (present(use_temp)) use_temp = associated(cs%tv%T)
4470end subroutine get_mom_state_elements
4471
4472!> Find the global integrals of various quantities.
4473subroutine get_ocean_stocks(CS, mass, heat, salt, on_PE_only)
4474 type(mom_control_struct), intent(inout) :: cs !< MOM control structure
4475 real, optional, intent(out) :: heat !< The globally integrated integrated ocean heat [J].
4476 real, optional, intent(out) :: salt !< The globally integrated integrated ocean salt [kg].
4477 real, optional, intent(out) :: mass !< The globally integrated integrated ocean mass [kg].
4478 logical, optional, intent(in) :: on_pe_only !< If present and true, only sum on the local PE.
4479
4480 if (present(mass)) &
4481 mass = global_mass_integral(cs%h, cs%G, cs%GV, on_pe_only=on_pe_only)
4482 if (present(heat)) &
4483 heat = cs%US%Q_to_J_kg*cs%US%RZL2_to_kg * cs%tv%C_p * &
4484 global_mass_integral(cs%h, cs%G, cs%GV, cs%tv%T, on_pe_only=on_pe_only, tmp_scale=cs%US%C_to_degC)
4485 if (present(salt)) &
4486 salt = 1.0e-3 * global_mass_integral(cs%h, cs%G, cs%GV, cs%tv%S, on_pe_only=on_pe_only, unscale=cs%US%S_to_ppt)
4487
4488end subroutine get_ocean_stocks
4489
4490
4491!> Save restart/pickup files required to initialize the MOM6 internal state.
4492subroutine save_mom_restart(CS, directory, time, G, time_stamped, filename, &
4493 GV, num_rest_files, write_IC)
4494 type(mom_control_struct), intent(inout) :: cs
4495 !< MOM control structure
4496 character(len=*), intent(in) :: directory
4497 !< The directory where the restart files are to be written
4498 type(time_type), intent(in) :: time
4499 !< The current model time
4500 type(ocean_grid_type), intent(inout) :: g
4501 !< The ocean's grid structure
4502 logical, optional, intent(in) :: time_stamped
4503 !< If present and true, add time-stamp to the restart file names
4504 character(len=*), optional, intent(in) :: filename
4505 !< A filename that overrides the name in CS%restartfile
4506 type(verticalgrid_type), optional, intent(in) :: gv
4507 !< The ocean's vertical grid structure
4508 integer, optional, intent(out) :: num_rest_files
4509 !< number of restart files written
4510 logical, optional, intent(in) :: write_ic
4511 !< If present and true, initial conditions are being written
4512
4513 logical :: showcalltree
4514 showcalltree = calltree_showquery()
4515
4516 call cpu_clock_begin(id_clock_ocean) ; call cpu_clock_begin(id_clock_save_restart)
4517 if (showcalltree) call calltree_waypoint("About to call save_restart (step_MOM)")
4518 call save_restart(directory, time, g, cs%restart_CS, &
4519 time_stamped=time_stamped, filename=filename, gv=gv, &
4520 num_rest_files=num_rest_files, write_ic=write_ic)
4521 if (showcalltree) call calltree_waypoint("Done with call to save_restart (step_MOM)")
4522
4523 if (cs%use_particles) call particles_save_restart(cs%particles, cs%h, directory, time, time_stamped)
4524 call cpu_clock_end(id_clock_save_restart) ; call cpu_clock_end(id_clock_ocean)
4525end subroutine save_mom_restart
4526
4527
4528!> End of ocean model, including memory deallocation
4529subroutine mom_end(CS)
4530 type(mom_control_struct), intent(inout) :: cs !< MOM control structure
4531
4532 call cpu_clock_begin(id_clock_ocean) ; call cpu_clock_begin(id_clock_mom_end)
4533
4534 call mom_sum_output_end(cs%sum_output_CSp)
4535
4536 if (cs%use_ALE_algorithm) call ale_end(cs%ALE_CSp)
4537
4538 !deallocate porous topography variables
4539 deallocate(cs%pbv%por_face_areaU) ; deallocate(cs%pbv%por_face_areaV)
4540 deallocate(cs%pbv%por_layer_widthU) ; deallocate(cs%pbv%por_layer_widthV)
4541
4542 ! NOTE: Allocated in PressureForce_FV_Bouss
4543 if (associated(cs%tv%varT)) deallocate(cs%tv%varT)
4544
4545 call tracer_advect_end(cs%tracer_adv_CSp)
4546 call tracer_hor_diff_end(cs%tracer_diff_CSp)
4547 call tracer_registry_end(cs%tracer_Reg)
4548 call tracer_flow_control_end(cs%tracer_flow_CSp)
4549
4550 if (.not. cs%adiabatic) then
4551 call diabatic_driver_end(cs%diabatic_CSp)
4552 deallocate(cs%diabatic_CSp)
4553 endif
4554
4555 call mom_diagnostics_end(cs%diagnostics_CSp, cs%ADp, cs%CDp)
4556
4557 if (cs%offline_tracer_mode) call offline_transport_end(cs%offline_CSp)
4558
4559 if (cs%split .and. cs%use_alt_split) then
4560 call end_dyn_split_rk2b(cs%dyn_split_RK2b_CSp)
4561 elseif (cs%split) then
4562 call end_dyn_split_rk2(cs%dyn_split_RK2_CSp)
4563 elseif (cs%use_RK2) then
4564 call end_dyn_unsplit_rk2(cs%dyn_unsplit_RK2_CSp)
4565 else
4566 call end_dyn_unsplit(cs%dyn_unsplit_CSp)
4567 endif
4568
4569 if (cs%use_particles) then
4570 call particles_end(cs%particles, cs%h)
4571 deallocate(cs%particles)
4572 endif
4573
4574 call thickness_diffuse_end(cs%thickness_diffuse_CSp, cs%CDp)
4575 if (cs%interface_filter) call interface_filter_end(cs%interface_filter_CSp, cs%CDp)
4576 call varmix_end(cs%VarMix)
4577 call set_visc_end(cs%visc, cs%set_visc_CSp)
4578 call meke_end(cs%MEKE)
4579
4580 if (associated(cs%tv%internal_heat)) deallocate(cs%tv%internal_heat)
4581 if (associated(cs%tv%TempxPmE)) deallocate(cs%tv%TempxPmE)
4582
4583 dealloc_(cs%ave_ssh_ibc) ; dealloc_(cs%ssh_rint) ; dealloc_(cs%eta_av_bc)
4584
4585 ! TODO: debug_truncations deallocation
4586
4587 dealloc_(cs%uhtr) ; dealloc_(cs%vhtr)
4588
4589 if (associated(cs%Hml)) deallocate(cs%Hml)
4590 if (associated(cs%tv%salt_deficit)) deallocate(cs%tv%salt_deficit)
4591 if (associated(cs%tv%frazil)) deallocate(cs%tv%frazil)
4592 if (allocated(cs%tv%SpV_avg)) deallocate(cs%tv%SpV_avg)
4593
4594 if (associated(cs%tv%T)) then
4595 dealloc_(cs%T) ; cs%tv%T => null() ; dealloc_(cs%S) ; cs%tv%S => null()
4596 endif
4597
4598 dealloc_(cs%u) ; dealloc_(cs%v) ; dealloc_(cs%h)
4599 dealloc_(cs%uh) ; dealloc_(cs%vh)
4600
4601 if (associated(cs%update_OBC_CSp)) call obc_register_end(cs%update_OBC_CSp)
4602 if (associated(cs%OBC)) call open_boundary_end(cs%OBC)
4603
4604 call verticalgridend(cs%GV)
4605 call mom_grid_end(cs%G)
4606
4607 if (cs%debug .or. cs%G%symmetric) &
4608 call deallocate_mom_domain(cs%G%Domain_aux)
4609
4610 if (cs%rotate_index) &
4611 call deallocate_mom_domain(cs%G%Domain)
4612
4613 ! The MPP domains may be needed by an external coupler, so use `cursory`.
4614 ! TODO: This may create a domain memory leak, and needs investigation.
4615 call deallocate_mom_domain(cs%G_in%domain, cursory=.true.)
4616
4617 call unit_scaling_end(cs%US)
4618
4619 call cpu_clock_end(id_clock_mom_end) ; call cpu_clock_end(id_clock_ocean)
4620
4621end subroutine mom_end
4622
4623!> \namespace mom
4624!!
4625!! Modular Ocean Model (MOM) Version 6.0 (MOM6)
4626!!
4627!! \authors Alistair Adcroft, Robert Hallberg, and Stephen Griffies
4628!!
4629!! Additional contributions from:
4630!! * Whit Anderson
4631!! * Brian Arbic
4632!! * Will Cooke
4633!! * Anand Gnanadesikan
4634!! * Matthew Harrison
4635!! * Mehmet Ilicak
4636!! * Laura Jackson
4637!! * Jasmine John
4638!! * John Krasting
4639!! * Zhi Liang
4640!! * Bonnie Samuels
4641!! * Harper Simmons
4642!! * Laurent White
4643!! * Niki Zadeh
4644!!
4645!! MOM ice-shelf code was developed by
4646!! * Daniel Goldberg
4647!! * Robert Hallberg
4648!! * Chris Little
4649!! * Olga Sergienko
4650!!
4651!! \section section_overview Overview of MOM
4652!!
4653!! This program (MOM) simulates the ocean by numerically solving
4654!! the hydrostatic primitive equations in generalized Lagrangian
4655!! vertical coordinates, typically tracking stretched pressure (p*)
4656!! surfaces or following isopycnals in the ocean's interior, and
4657!! general orthogonal horizontal coordinates. Unlike earlier versions
4658!! of MOM, in MOM6 these equations are horizontally discretized on an
4659!! Arakawa C-grid. (It remains to be seen whether a B-grid dynamic
4660!! core will be revived in MOM6 at a later date; for now applications
4661!! requiring a B-grid discretization should use MOM5.1.) MOM6 offers
4662!! a range of options for the physical parameterizations, from those
4663!! most appropriate to highly idealized models for geophysical fluid
4664!! dynamics studies to a rich suite of processes appropriate for
4665!! realistic ocean simulations. The thermodynamic options typically
4666!! use conservative temperature and preformed salinity as conservative
4667!! state variables and a full nonlinear equation of state, but there
4668!! are also idealized adiabatic configurations of the model that use
4669!! fixed density layers. Version 6.0 of MOM continues in the long
4670!! tradition of a commitment to climate-quality ocean simulations
4671!! embodied in previous versions of MOM, even as it draws extensively
4672!! on the lessons learned in the development of the Generalized Ocean
4673!! Layered Dynamics (GOLD) ocean model, which was also primarily
4674!! developed at NOAA/GFDL. MOM has also benefited tremendously from
4675!! the FMS infrastructure, which it utilizes and shares with other
4676!! component models developed at NOAA/GFDL.
4677!!
4678!! When run is isopycnal-coordinate mode, the uppermost few layers
4679!! are often used to describe a bulk mixed layer, including the
4680!! effects of penetrating shortwave radiation. Either a split-
4681!! explicit time stepping scheme or a non-split scheme may be used
4682!! for the dynamics, while the time stepping may be split (and use
4683!! different numbers of steps to cover the same interval) for the
4684!! forcing, the thermodynamics, and for the dynamics. Most of the
4685!! numerics are second order accurate in space. MOM can run with an
4686!! absurdly thin minimum layer thickness. A variety of non-isopycnal
4687!! vertical coordinate options are under development, but all exploit
4688!! the advantages of a Lagrangian vertical coordinate, as discussed
4689!! in detail by Adcroft and Hallberg (Ocean Modelling, 2006).
4690!!
4691!! Details of the numerics and physical parameterizations are
4692!! provided in the appropriate source files. All of the available
4693!! options are selected at run-time by parsing the input files,
4694!! usually MOM_input and MOM_override, and the options choices are
4695!! then documented for each run in MOM_param_docs.
4696!!
4697!! MOM6 integrates the equations forward in time in three distinct
4698!! phases. In one phase, the dynamic equations for the velocities
4699!! and layer thicknesses are advanced, capturing the propagation of
4700!! external and internal inertia-gravity waves, Rossby waves, and
4701!! other strictly adiabatic processes, including lateral stresses,
4702!! vertical viscosity and momentum forcing, and interface height
4703!! diffusion (commonly called Gent-McWilliams diffusion in depth-
4704!! coordinate models). In the second phase, all tracers are advected
4705!! and diffused along the layers. The third phase applies diabatic
4706!! processes, vertical mixing of water properties, and perhaps
4707!! vertical remapping to cause the layers to track the desired
4708!! vertical coordinate.
4709!!
4710!! The present file (MOM.F90) orchestrates the main time stepping
4711!! loops. One time integration option for the dynamics uses a split
4712!! explicit time stepping scheme to rapidly step the barotropic
4713!! pressure and velocity fields. The barotropic velocities are
4714!! averaged over the baroclinic time step before they are used to
4715!! advect thickness and determine the baroclinic accelerations. As
4716!! described in Hallberg and Adcroft (2009), a barotropic correction
4717!! is applied to the time-mean layer velocities to ensure that the
4718!! sum of the layer transports agrees with the time-mean barotropic
4719!! transport, thereby ensuring that the estimates of the free surface
4720!! from the sum of the layer thicknesses agrees with the final free
4721!! surface height as calculated by the barotropic solver. The
4722!! barotropic and baroclinic velocities are kept consistent by
4723!! recalculating the barotropic velocities from the baroclinic
4724!! transports each time step. This scheme is described in Hallberg,
4725!! 1997, J. Comp. Phys. 135, 54-65 and in Hallberg and Adcroft, 2009,
4726!! Ocean Modelling, 29, 15-26.
4727!!
4728!! The other time integration options use non-split time stepping
4729!! schemes based on the 3-step third order Runge-Kutta scheme
4730!! described in Matsuno, 1966, J. Met. Soc. Japan, 44, 85-88, or on
4731!! a two-step quasi-2nd order Runge-Kutta scheme. These are much
4732!! slower than the split time-stepping scheme, but they are useful
4733!! for providing a more robust solution for debugging cases where the
4734!! more complicated split time-stepping scheme may be giving suspect
4735!! solutions.
4736!!
4737!! There are a range of closure options available. Horizontal
4738!! velocities are subject to a combination of horizontal biharmonic
4739!! and Laplacian friction (based on a stress tensor formalism) and a
4740!! vertical Fickian viscosity (perhaps using the kinematic viscosity
4741!! of water). The horizontal viscosities may be constant, spatially
4742!! varying or may be dynamically calculated using Smagorinsky's
4743!! approach. A diapycnal diffusion of density and thermodynamic
4744!! quantities is also allowed, but not required, as is horizontal
4745!! diffusion of interface heights (akin to the Gent-McWilliams
4746!! closure of geopotential coordinate models). The diapycnal mixing
4747!! may use a fixed diffusivity or it may use the shear Richardson
4748!! number dependent closure, like that described in Jackson et al.
4749!! (JPO, 2008). When there is diapycnal diffusion, it applies to
4750!! momentum as well. As this is in addition to the vertical viscosity,
4751!! the vertical Prandtl always exceeds 1. A refined bulk-mixed layer
4752!! is often used to describe the planetary boundary layer in realistic
4753!! ocean simulations.
4754!!
4755!! MOM has a number of noteworthy debugging capabilities.
4756!! Excessively large velocities are truncated and MOM will stop
4757!! itself after a number of such instances to keep the model from
4758!! crashing altogether. This is useful in diagnosing failures,
4759!! or (by accepting some truncations) it may be useful for getting
4760!! the model past the adjustment from an ill-balanced initial
4761!! condition. In addition, all of the accelerations in the columns
4762!! with excessively large velocities may be directed to a text file.
4763!! Parallelization errors may be diagnosed using the DEBUG option,
4764!! which causes extensive checksums to be written out along with
4765!! comments indicating where in the algorithm the sums originate and
4766!! what variable is being summed. The point where these checksums
4767!! differ between runs is usually a good indication of where in the
4768!! code the problem lies. All of the test cases provided with MOM
4769!! are routinely tested to ensure that they give bitwise identical
4770!! results regardless of the domain decomposition, or whether they
4771!! use static or dynamic memory allocation.
4772!!
4773!! \section section_structure Structure of MOM
4774!!
4775!! About 115 other files of source code and 4 header files comprise
4776!! the MOM code, although there are several hundred more files that
4777!! make up the FMS infrastructure upon which MOM is built. Each of
4778!! the MOM files contains comments documenting what it does, and
4779!! most of the file names are fairly self-evident. In addition, all
4780!! subroutines and data types are referenced via a module use, only
4781!! statement, and the module names are consistent with the file names,
4782!! so it is not too hard to find the source file for a subroutine.
4783!!
4784!! The typical MOM directory tree is as follows:
4785!!
4786!! \verbatim
4787!! ../MOM
4788!! |-- ac
4789!! |-- config_src
4790!! | |-- drivers
4791!! | ! |-- FMS_cap
4792!! | ! |-- ice_solo_driver
4793!! | ! |-- mct_cap
4794!! | ! |-- nuopc_cap
4795!! | ! |-- solo_driver
4796!! | ! `-- unit_drivers
4797!! | |-- external
4798!! | ! |-- drifters
4799!! | ! |-- GFDL_ocean_BGC
4800!! | ! `-- ODA_hooks
4801!! | |-- infra
4802!! | ! |-- FMS1
4803!! | ! `-- FMS2
4804!! | `-- memory
4805!! | ! |-- dynamic_nonsymmetric
4806!! | ! `-- dynamic_symmetric
4807!! |-- docs
4808!! |-- pkg
4809!! | |-- CVMix-src
4810!! | |-- ...
4811!! | `-- MOM6_DA_hooks
4812!! `-- src
4813!! |-- ALE
4814!! |-- core
4815!! |-- diagnostics
4816!! |-- equation_of_state
4817!! |-- framework
4818!! |-- ice_shelf
4819!! |-- initialization
4820!! |-- ocean_data_assim
4821!! |-- parameterizations
4822!! | |-- CVMix
4823!! | |-- lateral
4824!! | `-- vertical
4825!! |-- tracer
4826!! `-- user
4827!! \endverbatim
4828!!
4829!! Rather than describing each file here, selected directory contents
4830!! will be described to give a broad overview of the MOM code
4831!! structure.
4832!!
4833!! The directories under config_src contain files that are used for
4834!! configuring the code, for instance for coupled or ocean-only runs.
4835!! Only one or two of these directories are used in compiling any,
4836!! particular run.
4837!!
4838!! * config_src/drivers/FMS-cap:
4839!! The files here are used to couple MOM as a component in a larger
4840!! run driven by the FMS coupler. This includes code that converts
4841!! various forcing fields into the code structures and flux and unit
4842!! conventions used by MOM, and converts the MOM surface fields
4843!! back to the forms used by other FMS components.
4844!!
4845!! * config_src/drivers/nuopc-cap:
4846!! The files here are used to couple MOM as a component in a larger
4847!! run driven by the NUOPC coupler. This includes code that converts
4848!! various forcing fields into the code structures and flux and unit
4849!! conventions used by MOM, and converts the MOM surface fields
4850!! back to the forms used by other NUOPC components.
4851!!
4852!! * config_src/drivers/solo_driver:
4853!! The files here are include the _main driver that is used when
4854!! MOM is configured as an ocean-only model, as well as the files
4855!! that specify the surface forcing in this configuration.
4856!!
4857!! * config_src/external:
4858!! The files here are mostly just stubs, so that MOM6 can compile
4859!! with calls to the public interfaces external packages, but
4860!! without actually requiring those packages themselves. In more
4861!! elaborate configurations, would be linked to the actual code for
4862!! those external packages rather than these simple stubs.
4863!!
4864!! * config_src/memory/dynamic-symmetric:
4865!! The only file here is the version of MOM_memory.h that is used
4866!! for dynamic memory configurations of MOM.
4867!!
4868!! The directories under src contain most of the MOM files. These
4869!! files are used in every configuration using MOM.
4870!!
4871!! * src/core:
4872!! The files here constitute the MOM dynamic core. This directory
4873!! also includes files with the types that describe the model's
4874!! lateral grid and have defined types that are shared across
4875!! various MOM modules to allow for more succinct and flexible
4876!! subroutine argument lists.
4877!!
4878!! * src/diagnostics:
4879!! The files here calculate various diagnostics that are ancilliary
4880!! to the model itself. While most of these diagnostics do not
4881!! directly affect the model's solution, there are some, like the
4882!! calculation of the deformation radius, that are used in some
4883!! of the process parameterizations.
4884!!
4885!! * src/equation_of_state:
4886!! These files describe the physical properties of sea-water,
4887!! including both the equation of state and when it freezes.
4888!!
4889!! * src/framework:
4890!! These files provide infrastructure utilities for MOM. Many are
4891!! simply wrappers for capabilities provided by FMS, although others
4892!! provide capabilities (like the file_parser) that are unique to
4893!! MOM. When MOM is adapted to use a modeling infrastructure
4894!! distinct from FMS, most of the required changes are in this
4895!! directory.
4896!!
4897!! * src/initialization:
4898!! These are the files that are used to initialize the MOM grid
4899!! or provide the initial physical state for MOM. These files are
4900!! not intended to be modified, but provide a means for calling
4901!! user-specific initialization code like the examples in src/user.
4902!!
4903!! * src/parameterizations/lateral:
4904!! These files implement a number of quasi-lateral (along-layer)
4905!! process parameterizations, including lateral viscosities,
4906!! parameterizations of eddy effects, and the calculation of tidal
4907!! forcing.
4908!!
4909!! * src/parameterizations/vertical:
4910!! These files implement a number of vertical mixing or diabatic
4911!! processes, including the effects of vertical viscosity and
4912!! code to parameterize the planetary boundary layer. There is a
4913!! separate driver that orchestrates this portion of the algorithm,
4914!! and there is a diversity of parameterizations to be found here.
4915!!
4916!! * src/tracer:
4917!! These files handle the lateral transport and diffusion of
4918!! tracers, or are the code to implement various passive tracer
4919!! packages. Additional tracer packages are readily accommodated.
4920!!
4921!! * src/user:
4922!! These are either stub routines that a user could use to change
4923!! the model's initial conditions or forcing, or are examples that
4924!! implement specific test cases. These files can easily be hand
4925!! edited to create new analytically specified configurations.
4926!!
4927!!
4928!! Most simulations can be set up by modifying only the files
4929!! MOM_input, and possibly one or two of the files in src/user.
4930!! In addition, the diag_table (MOM_diag_table) will commonly be
4931!! modified to tailor the output to the needs of the question at
4932!! hand. The FMS utility mkmf works with a file called path_names
4933!! to build an appropriate makefile, and path_names should be edited
4934!! to reflect the actual location of the desired source code.
4935!!
4936!! The separate MOM-examples git repository provides a large number
4937!! of working configurations of MOM, along with reference solutions for several
4938!! different compilers on GFDL's latest large computer. The versions
4939!! of MOM_memory.h in these directories need not be used if dynamic
4940!! memory allocation is desired, and the answers should be unchanged.
4941!!
4942!!
4943!! There are 3 publicly visible subroutines in this file (MOM.F90).
4944!! * step_MOM steps MOM over a specified interval of time.
4945!! * MOM_initialize calls initialize and does other initialization
4946!! that does not warrant user modification.
4947!! * extract_surface_state determines the surface (bulk mixed layer
4948!! if traditional isopycnal vertical coordinate) properties of the
4949!! current model state and packages pointers to these fields into an
4950!! exported structure.
4951!!
4952!! The remaining subroutines in this file (src/core/MOM.F90) are:
4953!! * find_total_transport determines the barotropic mass transport.
4954!! * register_diags registers many diagnostic fields for the dynamic
4955!! solver, or of the main model variables.
4956!! * MOM_timing_init initializes various CPU time clocks.
4957!! * write_static_fields writes out various time-invariant fields.
4958!! * set_restart_fields is used to specify those fields that are
4959!! written to and read from the restart file.
4960!!
4961!! \section section_heat_budget Diagnosing MOM heat budget
4962!!
4963!! Here are some example heat budgets for the ALE version of MOM6.
4964!!
4965!! \subsection subsection_2d_heat_budget Depth integrated heat budget
4966!!
4967!! Depth integrated heat budget diagnostic for MOM.
4968!!
4969!! * OPOTTEMPTEND_2d = T_ADVECTION_XY_2d + OPOTTEMPPMDIFF_2d + HFDS + HFGEOU
4970!!
4971!! * T_ADVECTION_XY_2d = horizontal advection
4972!! * OPOTTEMPPMDIFF_2d = neutral diffusion
4973!! * HFDS = net surface boundary heat flux
4974!! * HFGEOU = geothermal heat flux
4975!!
4976!! * HFDS = net surface boundary heat flux entering the ocean
4977!! = rsntds + rlntds + hfls + hfss + heat_pme + hfsifrazil
4978!!
4979!! * More heat flux cross-checks
4980!! * hfds = net_heat_coupler + hfsifrazil + heat_pme
4981!! * heat_pme = heat_content_surfwater
4982!! = heat_content_massin + heat_content_massout
4983!! = heat_content_fprec + heat_content_cond + heat_content_vprec
4984!! + hfrunoffds + hfevapds + hfrainds
4985!!
4986!! \subsection subsection_3d_heat_budget Depth integrated heat budget
4987!!
4988!! Here is an example 3d heat budget diagnostic for MOM.
4989!!
4990!! * OPOTTEMPTEND = T_ADVECTION_XY + TH_TENDENCY_VERT_REMAP + OPOTTEMPDIFF + OPOTTEMPPMDIFF
4991!! + BOUNDARY_FORCING_HEAT_TENDENCY + FRAZIL_HEAT_TENDENCY
4992!!
4993!! * OPOTTEMPTEND = net tendency of heat as diagnosed in MOM.F90
4994!! * T_ADVECTION_XY = heating of a cell from lateral advection
4995!! * TH_TENDENCY_VERT_REMAP = heating of a cell from vertical remapping
4996!! * OPOTTEMPDIFF = heating of a cell from diabatic diffusion
4997!! * OPOTTEMPPMDIFF = heating of a cell from neutral diffusion
4998!! * BOUNDARY_FORCING_HEAT_TENDENCY = heating of cell from boundary fluxes
4999!! * FRAZIL_HEAT_TENDENCY = heating of cell from frazil
5000!!
5001!! * TH_TENDENCY_VERT_REMAP has zero vertical sum, as it redistributes heat in vertical.
5002!!
5003!! * OPOTTEMPDIFF has zero vertical sum, as it redistributes heat in the vertical.
5004!!
5005!! * BOUNDARY_FORCING_HEAT_TENDENCY generally has 3d structure, with k > 1 contributions from
5006!! penetrative shortwave, and from other fluxes for the case when layers are tiny, in which
5007!! case MOM6 partitions tendencies into k > 1 layers.
5008!!
5009!! * FRAZIL_HEAT_TENDENCY generally has 3d structure, since MOM6 frazil calculation checks the
5010!! full ocean column.
5011!!
5012!! * FRAZIL_HEAT_TENDENCY[k=\@sum] = HFSIFRAZIL = column integrated frazil heating.
5013!!
5014!! * HFDS = FRAZIL_HEAT_TENDENCY[k=\@sum] + BOUNDARY_FORCING_HEAT_TENDENCY[k=\@sum]
5015!!
5016!! Here is an example 2d heat budget (depth summed) diagnostic for MOM.
5017!!
5018!! * OPOTTEMPTEND_2d = T_ADVECTION_XY_2d + OPOTTEMPPMDIFF_2d + HFDS
5019!!
5020!!
5021!! Here is an example 3d salt budget diagnostic for MOM.
5022!!
5023!! * OSALTTEND = S_ADVECTION_XY + SH_TENDENCY_VERT_REMAP + OSALTDIFF + OSALTPMDIFF
5024!! + BOUNDARY_FORCING_SALT_TENDENCY
5025!!
5026!! * OSALTTEND = net tendency of salt as diagnosed in MOM.F90
5027!! * S_ADVECTION_XY = salt convergence to cell from lateral advection
5028!! * SH_TENDENCY_VERT_REMAP = salt convergence to cell from vertical remapping
5029!! * OSALTDIFF = salt convergence to cell from diabatic diffusion
5030!! * OSALTPMDIFF = salt convergence to cell from neutral diffusion
5031!! * BOUNDARY_FORCING_SALT_TENDENCY = salt convergence to cell from boundary fluxes
5032!!
5033!! * SH_TENDENCY_VERT_REMAP has zero vertical sum, as it redistributes salt in vertical.
5034!!
5035!! * OSALTDIFF has zero vertical sum, as it redistributes salt in the vertical.
5036!!
5037!! * BOUNDARY_FORCING_SALT_TENDENCY generally has 3d structure, with k > 1 contributions from
5038!! the case when layers are tiny, in which case MOM6 partitions tendencies into k > 1 layers.
5039!!
5040!! * SFDSI = BOUNDARY_FORCING_SALT_TENDENCY[k=\@sum]
5041!!
5042!! Here is an example 2d salt budget (depth summed) diagnostic for MOM.
5043!!
5044!! * OSALTTEND_2d = S_ADVECTION_XY_2d + OSALTPMDIFF_2d + SFDSI (+ SALT_FLUX_RESTORE)
5045!!
5046!!
5047!!
5048end module mom