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