MOM.F90

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