MOM6
Data Types
mom Module Reference

Detailed Description

The central module of the MOM6 ocean model.

Modular Ocean Model (MOM) Version 6.0 (MOM6)

Authors
Alistair Adcroft, Robert Hallberg, and Stephen Griffies

Additional contributions from:

  • Whit Anderson
  • Brian Arbic
  • Will Cooke
  • Anand Gnanadesikan
  • Matthew Harrison
  • Mehmet Ilicak
  • Laura Jackson
  • Jasmine John
  • John Krasting
  • Zhi Liang
  • Bonnie Samuels
  • Harper Simmons
  • Laurent White
  • Niki Zadeh

MOM ice-shelf code was developed by

  • Daniel Goldberg
  • Robert Hallberg
  • Chris Little
  • Olga Sergienko

Overview of MOM

This program (MOM) simulates the ocean by numerically solving the hydrostatic primitive equations in generalized Lagrangian vertical coordinates, typically tracking stretched pressure (p*) surfaces or following isopycnals in the ocean's interior, and general orthogonal horizontal coordinates. Unlike earlier versions of MOM, in MOM6 these equations are horizontally discretized on an Arakawa C-grid. (It remains to be seen whether a B-grid dynamic core will be revived in MOM6 at a later date; for now applications requiring a B-grid discretization should use MOM5.1.) MOM6 offers a range of options for the physical parameterizations, from those most appropriate to highly idealized models for geophysical fluid dynamics studies to a rich suite of processes appropriate for realistic ocean simulations. The thermodynamic options typically use conservative temperature and preformed salinity as conservative state variables and a full nonlinear equation of state, but there are also idealized adiabatic configurations of the model that use fixed density layers. Version 6.0 of MOM continues in the long tradition of a commitment to climate-quality ocean simulations embodied in previous versions of MOM, even as it draws extensively on the lessons learned in the development of the Generalized Ocean Layered Dynamics (GOLD) ocean model, which was also primarily developed at NOAA/GFDL. MOM has also benefited tremendously from the FMS infrastructure, which it utilizes and shares with other component models developed at NOAA/GFDL.

When run is isopycnal-coordinate mode, the uppermost few layers are often used to describe a bulk mixed layer, including the effects of penetrating shortwave radiation. Either a split- explicit time stepping scheme or a non-split scheme may be used for the dynamics, while the time stepping may be split (and use different numbers of steps to cover the same interval) for the forcing, the thermodynamics, and for the dynamics. Most of the numerics are second order accurate in space. MOM can run with an absurdly thin minimum layer thickness. A variety of non-isopycnal vertical coordinate options are under development, but all exploit the advantages of a Lagrangian vertical coordinate, as discussed in detail by Adcroft and Hallberg (Ocean Modelling, 2006).

Details of the numerics and physical parameterizations are provided in the appropriate source files. All of the available options are selected at run-time by parsing the input files, usually MOM_input and MOM_override, and the options choices are then documented for each run in MOM_param_docs.

MOM6 integrates the equations forward in time in three distinct phases. In one phase, the dynamic equations for the velocities and layer thicknesses are advanced, capturing the propagation of external and internal inertia-gravity waves, Rossby waves, and other strictly adiabatic processes, including lateral stresses, vertical viscosity and momentum forcing, and interface height diffusion (commonly called Gent-McWilliams diffusion in depth- coordinate models). In the second phase, all tracers are advected and diffused along the layers. The third phase applies diabatic processes, vertical mixing of water properties, and perhaps vertical remapping to cause the layers to track the desired vertical coordinate.

The present file (MOM.F90) orchestrates the main time stepping loops. One time integration option for the dynamics uses a split explicit time stepping scheme to rapidly step the barotropic pressure and velocity fields. The barotropic velocities are averaged over the baroclinic time step before they are used to advect thickness and determine the baroclinic accelerations. As described in Hallberg and Adcroft (2009), a barotropic correction is applied to the time-mean layer velocities to ensure that the sum of the layer transports agrees with the time-mean barotropic transport, thereby ensuring that the estimates of the free surface from the sum of the layer thicknesses agrees with the final free surface height as calculated by the barotropic solver. The barotropic and baroclinic velocities are kept consistent by recalculating the barotropic velocities from the baroclinic transports each time step. This scheme is described in Hallberg, 1997, J. Comp. Phys. 135, 54-65 and in Hallberg and Adcroft, 2009, Ocean Modelling, 29, 15-26.

The other time integration options use non-split time stepping schemes based on the 3-step third order Runge-Kutta scheme described in Matsuno, 1966, J. Met. Soc. Japan, 44, 85-88, or on a two-step quasi-2nd order Runge-Kutta scheme. These are much slower than the split time-stepping scheme, but they are useful for providing a more robust solution for debugging cases where the more complicated split time-stepping scheme may be giving suspect solutions.

There are a range of closure options available. Horizontal velocities are subject to a combination of horizontal biharmonic and Laplacian friction (based on a stress tensor formalism) and a vertical Fickian viscosity (perhaps using the kinematic viscosity of water). The horizontal viscosities may be constant, spatially varying or may be dynamically calculated using Smagorinsky's approach. A diapycnal diffusion of density and thermodynamic quantities is also allowed, but not required, as is horizontal diffusion of interface heights (akin to the Gent-McWilliams closure of geopotential coordinate models). The diapycnal mixing may use a fixed diffusivity or it may use the shear Richardson number dependent closure, like that described in Jackson et al. (JPO, 2008). When there is diapycnal diffusion, it applies to momentum as well. As this is in addition to the vertical viscosity, the vertical Prandtl always exceeds 1. A refined bulk-mixed layer is often used to describe the planetary boundary layer in realistic ocean simulations.

MOM has a number of noteworthy debugging capabilities. Excessively large velocities are truncated and MOM will stop itself after a number of such instances to keep the model from crashing altogether. This is useful in diagnosing failures, or (by accepting some truncations) it may be useful for getting the model past the adjustment from an ill-balanced initial condition. In addition, all of the accelerations in the columns with excessively large velocities may be directed to a text file. Parallelization errors may be diagnosed using the DEBUG option, which causes extensive checksums to be written out along with comments indicating where in the algorithm the sums originate and what variable is being summed. The point where these checksums differ between runs is usually a good indication of where in the code the problem lies. All of the test cases provided with MOM are routinely tested to ensure that they give bitwise identical results regardless of the domain decomposition, or whether they use static or dynamic memory allocation.

Structure of MOM

About 115 other files of source code and 4 header files comprise the MOM code, although there are several hundred more files that make up the FMS infrastructure upon which MOM is built. Each of the MOM files contains comments documenting what it does, and most of the file names are fairly self-evident. In addition, all subroutines and data types are referenced via a module use, only statement, and the module names are consistent with the file names, so it is not too hard to find the source file for a subroutine.

The typical MOM directory tree is as follows:

        ../MOM
        |-- config_src
        |   |-- coupled_driver
        |   |-- dynamic
        |   `-- solo_driver
        |-- examples
        |   |-- CM2G
        |   |-- ...
        |   `-- torus_advection_test
        `-- src
            |-- core
            |-- diagnostics
            |-- equation_of_state
            |-- framework
            |-- ice_shelf
            |-- initialization
            |-- parameterizations
            |   |-- lateral
            |   `-- vertical
            |-- tracer
            `-- user

Rather than describing each file here, each directory contents will be described to give a broad overview of the MOM code structure.

The directories under config_src contain files that are used for configuring the code, for instance for coupled or ocean-only runs. Only one or two of these directories are used in compiling any, particular run.

  • config_src/coupled_driver: The files here are used to couple MOM as a component in a larger run driven by the FMS coupler. This includes code that converts various forcing fields into the code structures and flux and unit conventions used by MOM, and converts the MOM surface fields back to the forms used by other FMS components.
  • config_src/dynamic: The only file here is the version of MOM_memory.h that is used for dynamic memory configurations of MOM.

  • config_src/solo_driver: The files here are include the _main driver that is used when MOM is configured as an ocean-only model, as well as the files that specify the surface forcing in this configuration.

    The directories under examples provide a large number of working configurations of MOM, along with reference solutions for several different compilers on GFDL's latest large computer. The versions of MOM_memory.h in these directories need not be used if dynamic memory allocation is desired, and the answers should be unchanged.

    The directories under src contain most of the MOM files. These files are used in every configuration using MOM.

  • src/core: The files here constitute the MOM dynamic core. This directory also includes files with the types that describe the model's lateral grid and have defined types that are shared across various MOM modules to allow for more succinct and flexible subroutine argument lists.
  • src/diagnostics: The files here calculate various diagnostics that are anciliary to the model itself. While most of these diagnostics do not directly affect the model's solution, there are some, like the calculation of the deformation radius, that are used in some of the process parameterizations.
  • src/equation_of_state: These files describe the physical properties of sea-water, including both the equation of state and when it freezes.
  • src/framework: These files provide infrastructure utilities for MOM. Many are simply wrappers for capabilities provided by FMS, although others provide capabilities (like the file_parser) that are unique to MOM. When MOM is adapted to use a modeling infrastructure distinct from FMS, most of the required changes are in this directory.
  • src/initialization: These are the files that are used to initialize the MOM grid or provide the initial physical state for MOM. These files are not intended to be modified, but provide a means for calling user-specific initialization code like the examples in src/user.
  • src/parameterizations/lateral: These files implement a number of quasi-lateral (along-layer) process parameterizations, including lateral viscosities, parameterizations of eddy effects, and the calculation of tidal forcing.
  • src/parameterizations/vertical: These files implement a number of vertical mixing or diabatic processes, including the effects of vertical viscosity and code to parameterize the planetary boundary layer. There is a separate driver that orchestrates this portion of the algorithm, and there is a diversity of parameterizations to be found here.
  • src/tracer: These files handle the lateral transport and diffusion of tracers, or are the code to implement various passive tracer packages. Additional tracer packages are readily accommodated.
  • src/user: These are either stub routines that a user could use to change the model's initial conditions or forcing, or are examples that implement specific test cases. These files can easily be hand edited to create new analytically specified configurations.

Most simulations can be set up by modifying only the files MOM_input, and possibly one or two of the files in src/user. In addition, the diag_table (MOM_diag_table) will commonly be modified to tailor the output to the needs of the question at hand. The FMS utility mkmf works with a file called path_names to build an appropriate makefile, and path_names should be edited to reflect the actual location of the desired source code.

There are 3 publicly visible subroutines in this file (MOM.F90).

  • step_MOM steps MOM over a specified interval of time.
  • MOM_initialize calls initialize and does other initialization that does not warrant user modification.

  • extract_surface_state determines the surface (bulk mixed layer if traditional isoycnal vertical coordinate) properties of the current model state and packages pointers to these fields into an exported structure.

    The remaining subroutines in this file (src/core/MOM.F90) are:

  • find_total_transport determines the barotropic mass transport.
  • register_diags registers many diagnostic fields for the dynamic solver, or of the main model variables.
  • MOM_timing_init initializes various CPU time clocks.
  • write_static_fields writes out various time-invariant fields.
  • set_restart_fields is used to specify those fields that are written to and read from the restart file.

Diagnosing MOM heat budget

Here are some example heat budgets for the ALE version of MOM6.

Depth integrated heat budget

Depth integrated heat budget diagnostic for MOM.

  • OPOTTEMPTEND_2d = T_ADVECTION_XY_2d + OPOTTEMPPMDIFF_2d + HFDS + HFGEOU
  • T_ADVECTION_XY_2d = horizontal advection
  • OPOTTEMPPMDIFF_2d = neutral diffusion
  • HFDS = net surface boundary heat flux
  • HFGEOU = geothermal heat flux
  • HFDS = net surface boundary heat flux entering the ocean = rsntds + rlntds + hfls + hfss + heat_pme + hfsifrazil
  • More heat flux cross-checks
    • hfds = net_heat_coupler + hfsifrazil + heat_pme
    • heat_pme = heat_content_surfwater = heat_content_massin + heat_content_massout = heat_content_fprec + heat_content_cond + heat_content_vprec
      • hfrunoffds + hfevapds + hfrainds

Depth integrated heat budget

Here is an example 3d heat budget diagnostic for MOM.

  • OPOTTEMPTEND = T_ADVECTION_XY + TH_TENDENCY_VERT_REMAP + OPOTTEMPDIFF + OPOTTEMPPMDIFF
    • BOUNDARY_FORCING_HEAT_TENDENCY + FRAZIL_HEAT_TENDENCY
  • OPOTTEMPTEND = net tendency of heat as diagnosed in MOM.F90
  • T_ADVECTION_XY = heating of a cell from lateral advection
  • TH_TENDENCY_VERT_REMAP = heating of a cell from vertical remapping
  • OPOTTEMPDIFF = heating of a cell from diabatic diffusion
  • OPOTTEMPPMDIFF = heating of a cell from neutral diffusion
  • BOUNDARY_FORCING_HEAT_TENDENCY = heating of cell from boundary fluxes
  • FRAZIL_HEAT_TENDENCY = heating of cell from frazil
  • TH_TENDENCY_VERT_REMAP has zero vertical sum, as it redistributes heat in vertical.
  • OPOTTEMPDIFF has zero vertical sum, as it redistributes heat in the vertical.
  • BOUNDARY_FORCING_HEAT_TENDENCY generally has 3d structure, with k > 1 contributions from penetrative shortwave, and from other fluxes for the case when layers are tiny, in which case MOM6 partitions tendencies into k > 1 layers.
  • FRAZIL_HEAT_TENDENCY generally has 3d structure, since MOM6 frazil calculation checks the full ocean column.
  • FRAZIL_HEAT_TENDENCY[k=@sum] = HFSIFRAZIL = column integrated frazil heating.

  • HFDS = FRAZIL_HEAT_TENDENCY[k=@sum] + BOUNDARY_FORCING_HEAT_TENDENCY[k=@sum]

    Here is an example 2d heat budget (depth summed) diagnostic for MOM.

  • OPOTTEMPTEND_2d = T_ADVECTION_XY_2d + OPOTTEMPPMDIFF_2d + HFDS

    Here is an example 3d salt budget diagnostic for MOM.

  • OSALTTEND = S_ADVECTION_XY + SH_TENDENCY_VERT_REMAP + OSALTDIFF + OSALTPMDIFF
    • BOUNDARY_FORCING_SALT_TENDENCY
  • OSALTTEND = net tendency of salt as diagnosed in MOM.F90
  • S_ADVECTION_XY = salt convergence to cell from lateral advection
  • SH_TENDENCY_VERT_REMAP = salt convergence to cell from vertical remapping
  • OSALTDIFF = salt convergence to cell from diabatic diffusion
  • OSALTPMDIFF = salt convergence to cell from neutral diffusion
  • BOUNDARY_FORCING_SALT_TENDENCY = salt convergence to cell from boundary fluxes
  • SH_TENDENCY_VERT_REMAP has zero vertical sum, as it redistributes salt in vertical.
  • OSALTDIFF has zero vertical sum, as it redistributes salt in the vertical.
  • BOUNDARY_FORCING_SALT_TENDENCY generally has 3d structure, with k > 1 contributions from the case when layers are tiny, in which case MOM6 partitions tendencies into k > 1 layers.

  • SFDSI = BOUNDARY_FORCING_SALT_TENDENCY[k=@sum]

    Here is an example 2d salt budget (depth summed) diagnostic for MOM.

  • OSALTTEND_2d = S_ADVECTION_XY_2d + OSALTPMDIFF_2d + SFDSI (+ SALT_FLUX_RESTORE)

Data Types

type  mom_control_struct
 Control structure for the MOM module, including the variables that describe the state of the ocean. More...
 
type  mom_diag_ids
 A structure with diagnostic IDs of the state variables. More...
 
integer id_clock_ocean
 CPU time clock IDs. More...
 
integer id_clock_dynamics
 CPU time clock IDs. More...
 
integer id_clock_thermo
 CPU time clock IDs. More...
 
integer id_clock_tracer
 CPU time clock IDs. More...
 
integer id_clock_diabatic
 CPU time clock IDs. More...
 
integer id_clock_adiabatic
 CPU time clock IDs. More...
 
integer id_clock_continuity
 CPU time clock IDs. More...
 
integer id_clock_thick_diff
 CPU time clock IDs. More...
 
integer id_clock_bbl_visc
 CPU time clock IDs. More...
 
integer id_clock_ml_restrat
 CPU time clock IDs. More...
 
integer id_clock_diagnostics
 CPU time clock IDs. More...
 
integer id_clock_z_diag
 CPU time clock IDs. More...
 
integer id_clock_init
 CPU time clock IDs. More...
 
integer id_clock_mom_init
 CPU time clock IDs. More...
 
integer id_clock_pass
 CPU time clock IDs. More...
 
integer id_clock_pass_init
 CPU time clock IDs. More...
 
integer id_clock_ale
 CPU time clock IDs. More...
 
integer id_clock_other
 CPU time clock IDs. More...
 
integer id_clock_offline_tracer
 CPU time clock IDs. More...
 
subroutine, public step_mom (forces, fluxes, sfc_state, Time_start, time_int_in, CS, Waves, do_dynamics, do_thermodynamics, start_cycle, end_cycle, cycle_length, reset_therm)
 This subroutine orchestrates the time stepping of MOM. The adiabatic dynamics are stepped by calls to one of the step_MOM_dyn_...routines. The action of lateral processes on tracers occur in calls to advect_tracer and tracer_hordiff. Vertical mixing and possibly remapping occur inside of diabatic. More...
 
subroutine step_mom_dynamics (forces, p_surf_begin, p_surf_end, dt, dt_thermo, bbl_time_int, CS, Time_local, Waves)
 Time step the ocean dynamics, including the momentum and continuity equations. More...
 
subroutine step_mom_tracer_dyn (CS, G, GV, US, h, Time_local)
 step_MOM_tracer_dyn does tracer advection and lateral diffusion, bringing the tracers up to date with the changes in state due to the dynamics. Surface sources and sinks and remapping are handled via step_MOM_thermo. More...
 
subroutine step_mom_thermo (CS, G, GV, US, u, v, h, tv, fluxes, dtdia, Time_end_thermo, update_BBL, Waves)
 MOM_step_thermo orchestrates the thermodynamic time stepping and vertical remapping, via calls to diabatic (or adiabatic) and ALE_main. More...
 
subroutine, public step_offline (forces, fluxes, sfc_state, Time_start, time_interval, CS)
 step_offline is the main driver for running tracers offline in MOM6. This has been primarily developed with ALE configurations in mind. Some work has been done in isopycnal configuration, but the work is very preliminary. Some more detail about this capability along with some of the subroutines called here can be found in tracers/MOM_offline_control.F90 More...
 
subroutine, public initialize_mom (Time, Time_init, param_file, dirs, CS, restart_CSp, Time_in, offline_tracer_mode, input_restart_file, diag_ptr, count_calls, tracer_flow_CSp)
 Initialize MOM, including memory allocation, setting up parameters and diagnostics, initializing the ocean state variables, and initializing subsidiary modules. More...
 
subroutine, public finish_mom_initialization (Time, dirs, CS, restart_CSp)
 Finishes initializing MOM and writes out the initial conditions. More...
 
subroutine register_diags (Time, G, GV, US, IDs, diag)
 Register certain diagnostics. More...
 
subroutine mom_timing_init (CS)
 Set up CPU clock IDs for timing various subroutines. More...
 
subroutine set_restart_fields (GV, US, param_file, CS, restart_CSp)
 Set the fields that are needed for bitwise identical restarting the time stepping scheme. In addition to those specified here directly, there may be fields related to the forcing or to the barotropic solver that are needed; these are specified in sub- routines that are called from this one. More...
 
subroutine adjust_ssh_for_p_atm (tv, G, GV, US, ssh, p_atm, use_EOS)
 Apply a correction to the sea surface height to compensate for the atmospheric pressure (the inverse barometer). More...
 
subroutine, public extract_surface_state (CS, sfc_state)
 Set the surface (return) properties of the ocean model by setting the appropriate fields in sfc_state. Unused fields are set to NULL or are unallocated. More...
 
logical function, public mom_state_is_synchronized (CS, adv_dyn)
 Return true if all phases of step_MOM are at the same point in time. More...
 
subroutine, public get_mom_state_elements (CS, G, GV, US, C_p, use_temp)
 This subroutine offers access to values or pointers to other types from within the MOM_control_struct, allowing the MOM_control_struct to be opaque. More...
 
subroutine, public get_ocean_stocks (CS, mass, heat, salt, on_PE_only)
 Find the global integrals of various quantities. More...
 
subroutine, public mom_end (CS)
 End of ocean model, including memory deallocation. More...
 

Function/Subroutine Documentation

◆ adjust_ssh_for_p_atm()

subroutine mom::adjust_ssh_for_p_atm ( type(thermo_var_ptrs), intent(in)  tv,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(inout)  ssh,
real, dimension(:,:), optional, pointer  p_atm,
logical, intent(in), optional  use_EOS 
)
private

Apply a correction to the sea surface height to compensate for the atmospheric pressure (the inverse barometer).

Parameters
[in]tvA structure pointing to various thermodynamic variables
[in]gocean grid structure
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
[in,out]sshtime mean surface height [m]
p_atmatmospheric pressure [Pa]
[in]use_eosIf true, calculate the density for the SSH correction using the equation of state.

Definition at line 2673 of file MOM.F90.

2673  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
2674  type(ocean_grid_type), intent(in) :: G !< ocean grid structure
2675  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2676  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2677  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: ssh !< time mean surface height [m]
2678  real, dimension(:,:), optional, pointer :: p_atm !< atmospheric pressure [Pa]
2679  logical, optional, intent(in) :: use_EOS !< If true, calculate the density for
2680  !! the SSH correction using the equation of state.
2681 
2682  real :: Rho_conv ! The density used to convert surface pressure to
2683  ! a corrected effective SSH [R ~> kg m-3].
2684  real :: IgR0 ! The SSH conversion factor from Pa to m [m Pa-1].
2685  logical :: calc_rho
2686  integer :: i, j, is, ie, js, je
2687 
2688  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2689  if (present(p_atm)) then ; if (associated(p_atm)) then
2690  calc_rho = associated(tv%eqn_of_state)
2691  if (present(use_eos) .and. calc_rho) calc_rho = use_eos
2692  ! Correct the output sea surface height for the contribution from the
2693  ! atmospheric pressure
2694  do j=js,je ; do i=is,ie
2695  if (calc_rho) then
2696  call calculate_density(tv%T(i,j,1), tv%S(i,j,1), p_atm(i,j)/2.0, &
2697  rho_conv, tv%eqn_of_state, scale=us%kg_m3_to_R)
2698  else
2699  rho_conv = gv%Rho0
2700  endif
2701  igr0 = 1.0 / (rho_conv * us%R_to_kg_m3*gv%mks_g_Earth)
2702  ssh(i,j) = ssh(i,j) + p_atm(i,j) * igr0
2703  enddo ; enddo
2704  endif ; endif
2705 

Referenced by step_mom(), and step_offline().

Here is the caller graph for this function:

◆ extract_surface_state()

subroutine, public mom::extract_surface_state ( type(mom_control_struct), pointer  CS,
type(surface), intent(inout)  sfc_state 
)

Set the surface (return) properties of the ocean model by setting the appropriate fields in sfc_state. Unused fields are set to NULL or are unallocated.

Parameters
csMaster MOM control structure
[in,out]sfc_statetransparent ocean surface state structure shared with the calling routine data in this structure is intent out.

Definition at line 2712 of file MOM.F90.

2712  type(MOM_control_struct), pointer :: CS !< Master MOM control structure
2713  type(surface), intent(inout) :: sfc_state !< transparent ocean surface state
2714  !! structure shared with the calling routine
2715  !! data in this structure is intent out.
2716 
2717  ! local
2718  real :: hu, hv ! Thicknesses interpolated to velocity points [H ~> m or kg m-2]
2719  type(ocean_grid_type), pointer :: G => null() !< pointer to a structure containing
2720  !! metrics and related information
2721  type(verticalGrid_type), pointer :: GV => null() !< structure containing vertical grid info
2722  type(unit_scale_type), pointer :: US => null() !< structure containing various unit conversion factors
2723  real, dimension(:,:,:), pointer :: &
2724  h => null() !< h : layer thickness [H ~> m or kg m-2]
2725  real :: depth(SZI_(CS%G)) !< Distance from the surface in depth units [Z ~> m]
2726  real :: depth_ml !< Depth over which to average to determine mixed
2727  !! layer properties [Z ~> m]
2728  real :: dh !< Thickness of a layer within the mixed layer [Z ~> m]
2729  real :: mass !< Mass per unit area of a layer [kg m-2]
2730  real :: bathy_m !< The depth of bathymetry [m] (not Z), used for error checking.
2731  real :: T_freeze !< freezing temperature [degC]
2732  real :: delT(SZI_(CS%G)) !< T-T_freeze [degC]
2733  logical :: use_temperature !< If true, temp and saln used as state variables.
2734  integer :: i, j, k, is, ie, js, je, nz, numberOfErrors, ig, jg
2735  integer :: isd, ied, jsd, jed
2736  integer :: iscB, iecB, jscB, jecB, isdB, iedB, jsdB, jedB
2737  logical :: localError
2738  character(240) :: msg
2739 
2740  call calltree_enter("extract_surface_state(), MOM.F90")
2741  g => cs%G ; gv => cs%GV ; us => cs%US
2742  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
2743  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2744  iscb = g%iscB ; iecb = g%iecB; jscb = g%jscB ; jecb = g%jecB
2745  isdb = g%isdB ; iedb = g%iedB; jsdb = g%jsdB ; jedb = g%jedB
2746  h => cs%h
2747 
2748  use_temperature = associated(cs%tv%T)
2749 
2750  if (.not.sfc_state%arrays_allocated) then
2751  ! Consider using a run-time flag to determine whether to do the vertical
2752  ! integrals, since the 3-d sums are not negligible in cost.
2753  call allocate_surface_state(sfc_state, g, use_temperature, do_integrals=.true.)
2754  endif
2755  sfc_state%frazil => cs%tv%frazil
2756  sfc_state%T_is_conT = cs%tv%T_is_conT
2757  sfc_state%S_is_absS = cs%tv%S_is_absS
2758 
2759  do j=js,je ; do i=is,ie
2760  sfc_state%sea_lev(i,j) = cs%ave_ssh_ibc(i,j)
2761  enddo ; enddo
2762 
2763  ! copy Hml into sfc_state, so that caps can access it
2764  if (associated(cs%Hml)) then
2765  do j=js,je ; do i=is,ie
2766  sfc_state%Hml(i,j) = cs%Hml(i,j)
2767  enddo ; enddo
2768  endif
2769 
2770  if (cs%Hmix < 0.0) then ! A bulk mixed layer is in use, so layer 1 has the properties
2771  if (use_temperature) then ; do j=js,je ; do i=is,ie
2772  sfc_state%SST(i,j) = cs%tv%T(i,j,1)
2773  sfc_state%SSS(i,j) = cs%tv%S(i,j,1)
2774  enddo ; enddo ; endif
2775  do j=js,je ; do i=is-1,ie
2776  sfc_state%u(i,j) = us%L_T_to_m_s * cs%u(i,j,1)
2777  enddo ; enddo
2778  do j=js-1,je ; do i=is,ie
2779  sfc_state%v(i,j) = us%L_T_to_m_s * cs%v(i,j,1)
2780  enddo ; enddo
2781 
2782  else ! (CS%Hmix >= 0.0)
2783  !### This calculation should work in thickness (H) units instead of Z, but that
2784  !### would change answers at roundoff in non-Boussinesq cases.
2785  depth_ml = cs%Hmix
2786  ! Determine the mean tracer properties of the uppermost depth_ml fluid.
2787  !$OMP parallel do default(shared) private(depth,dh)
2788  do j=js,je
2789  do i=is,ie
2790  depth(i) = 0.0
2791  if (use_temperature) then
2792  sfc_state%SST(i,j) = 0.0 ; sfc_state%SSS(i,j) = 0.0
2793  else
2794  sfc_state%sfc_density(i,j) = 0.0
2795  endif
2796  enddo
2797 
2798  do k=1,nz ; do i=is,ie
2799  if (depth(i) + h(i,j,k)*gv%H_to_Z < depth_ml) then
2800  dh = h(i,j,k)*gv%H_to_Z
2801  elseif (depth(i) < depth_ml) then
2802  dh = depth_ml - depth(i)
2803  else
2804  dh = 0.0
2805  endif
2806  if (use_temperature) then
2807  sfc_state%SST(i,j) = sfc_state%SST(i,j) + dh * cs%tv%T(i,j,k)
2808  sfc_state%SSS(i,j) = sfc_state%SSS(i,j) + dh * cs%tv%S(i,j,k)
2809  else
2810  sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) + dh * us%R_to_kg_m3*gv%Rlay(k)
2811  endif
2812  depth(i) = depth(i) + dh
2813  enddo ; enddo
2814  ! Calculate the average properties of the mixed layer depth.
2815  do i=is,ie
2816  if (depth(i) < gv%H_subroundoff*gv%H_to_Z) &
2817  depth(i) = gv%H_subroundoff*gv%H_to_Z
2818  if (use_temperature) then
2819  sfc_state%SST(i,j) = sfc_state%SST(i,j) / depth(i)
2820  sfc_state%SSS(i,j) = sfc_state%SSS(i,j) / depth(i)
2821  else
2822  sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) / depth(i)
2823  endif
2824  !### Verify that this is no longer needed.
2825  ! sfc_state%Hml(i,j) = US%Z_to_m * depth(i)
2826  enddo
2827  enddo ! end of j loop
2828 
2829 ! Determine the mean velocities in the uppermost depth_ml fluid.
2830  ! NOTE: Velocity loops start on `[ij]s-1` in order to update halo values
2831  ! required by the speed diagnostic on the non-symmetric grid.
2832  ! This assumes that u and v halos have already been updated.
2833  if (cs%Hmix_UV>0.) then
2834  !### This calculation should work in thickness (H) units instead of Z, but that
2835  !### would change answers at roundoff in non-Boussinesq cases.
2836  depth_ml = cs%Hmix_UV
2837  !$OMP parallel do default(shared) private(depth,dh,hv)
2838  do j=js-1,ie
2839  do i=is,ie
2840  depth(i) = 0.0
2841  sfc_state%v(i,j) = 0.0
2842  enddo
2843  do k=1,nz ; do i=is,ie
2844  hv = 0.5 * (h(i,j,k) + h(i,j+1,k)) * gv%H_to_Z
2845  if (depth(i) + hv < depth_ml) then
2846  dh = hv
2847  elseif (depth(i) < depth_ml) then
2848  dh = depth_ml - depth(i)
2849  else
2850  dh = 0.0
2851  endif
2852  sfc_state%v(i,j) = sfc_state%v(i,j) + dh * us%L_T_to_m_s * cs%v(i,j,k)
2853  depth(i) = depth(i) + dh
2854  enddo ; enddo
2855  ! Calculate the average properties of the mixed layer depth.
2856  do i=is,ie
2857  if (depth(i) < gv%H_subroundoff*gv%H_to_Z) &
2858  depth(i) = gv%H_subroundoff*gv%H_to_Z
2859  sfc_state%v(i,j) = sfc_state%v(i,j) / depth(i)
2860  enddo
2861  enddo ! end of j loop
2862 
2863  !$OMP parallel do default(shared) private(depth,dh,hu)
2864  do j=js,je
2865  do i=is-1,ie
2866  depth(i) = 0.0
2867  sfc_state%u(i,j) = 0.0
2868  enddo
2869  do k=1,nz ; do i=is-1,ie
2870  hu = 0.5 * (h(i,j,k) + h(i+1,j,k)) * gv%H_to_Z
2871  if (depth(i) + hu < depth_ml) then
2872  dh = hu
2873  elseif (depth(i) < depth_ml) then
2874  dh = depth_ml - depth(i)
2875  else
2876  dh = 0.0
2877  endif
2878  sfc_state%u(i,j) = sfc_state%u(i,j) + dh * us%L_T_to_m_s * cs%u(i,j,k)
2879  depth(i) = depth(i) + dh
2880  enddo ; enddo
2881  ! Calculate the average properties of the mixed layer depth.
2882  do i=is-1,ie
2883  if (depth(i) < gv%H_subroundoff*gv%H_to_Z) &
2884  depth(i) = gv%H_subroundoff*gv%H_to_Z
2885  sfc_state%u(i,j) = sfc_state%u(i,j) / depth(i)
2886  enddo
2887  enddo ! end of j loop
2888  else ! Hmix_UV<=0.
2889  do j=js,je ; do i=is-1,ie
2890  sfc_state%u(i,j) = us%L_T_to_m_s * cs%u(i,j,1)
2891  enddo ; enddo
2892  do j=js-1,je ; do i=is,ie
2893  sfc_state%v(i,j) = us%L_T_to_m_s * cs%v(i,j,1)
2894  enddo ; enddo
2895  endif
2896  endif ! (CS%Hmix >= 0.0)
2897 
2898 
2899  if (allocated(sfc_state%melt_potential)) then
2900  !$OMP parallel do default(shared)
2901  do j=js,je
2902  do i=is,ie
2903  depth(i) = 0.0
2904  delt(i) = 0.0
2905  enddo
2906 
2907  do k=1,nz ; do i=is,ie
2908  depth_ml = min(cs%HFrz,cs%visc%MLD(i,j))
2909  if (depth(i) + h(i,j,k)*gv%H_to_m < depth_ml) then
2910  dh = h(i,j,k)*gv%H_to_m
2911  elseif (depth(i) < depth_ml) then
2912  dh = depth_ml - depth(i)
2913  else
2914  dh = 0.0
2915  endif
2916 
2917  ! p=0 OK, HFrz ~ 10 to 20m
2918  call calculate_tfreeze(cs%tv%S(i,j,k), 0.0, t_freeze, cs%tv%eqn_of_state)
2919  depth(i) = depth(i) + dh
2920  delt(i) = delt(i) + dh * (cs%tv%T(i,j,k) - t_freeze)
2921  enddo ; enddo
2922 
2923  do i=is,ie
2924  ! set melt_potential to zero to avoid passing previous values
2925  sfc_state%melt_potential(i,j) = 0.0
2926 
2927  if (g%mask2dT(i,j)>0.) then
2928  ! instantaneous melt_potential [J m-2]
2929  sfc_state%melt_potential(i,j) = cs%tv%C_p * us%R_to_kg_m3*gv%Rho0 * delt(i)
2930  endif
2931  enddo
2932  enddo ! end of j loop
2933  endif ! melt_potential
2934 
2935  if (allocated(sfc_state%salt_deficit) .and. associated(cs%tv%salt_deficit)) then
2936  !$OMP parallel do default(shared)
2937  do j=js,je ; do i=is,ie
2938  ! Convert from gSalt to kgSalt
2939  sfc_state%salt_deficit(i,j) = 1000.0 * us%R_to_kg_m3*us%Z_to_m*cs%tv%salt_deficit(i,j)
2940  enddo ; enddo
2941  endif
2942  if (allocated(sfc_state%TempxPmE) .and. associated(cs%tv%TempxPmE)) then
2943  !$OMP parallel do default(shared)
2944  do j=js,je ; do i=is,ie
2945  sfc_state%TempxPmE(i,j) = us%R_to_kg_m3*us%Z_to_m*cs%tv%TempxPmE(i,j)
2946  enddo ; enddo
2947  endif
2948  if (allocated(sfc_state%internal_heat) .and. associated(cs%tv%internal_heat)) then
2949  !$OMP parallel do default(shared)
2950  do j=js,je ; do i=is,ie
2951  sfc_state%internal_heat(i,j) = cs%tv%internal_heat(i,j)
2952  enddo ; enddo
2953  endif
2954  if (allocated(sfc_state%taux_shelf) .and. associated(cs%visc%taux_shelf)) then
2955  !$OMP parallel do default(shared)
2956  do j=js,je ; do i=is-1,ie
2957  sfc_state%taux_shelf(i,j) = us%R_to_kg_m3*us%L_T_to_m_s**2*us%Z_to_L*cs%visc%taux_shelf(i,j)
2958  enddo ; enddo
2959  endif
2960  if (allocated(sfc_state%tauy_shelf) .and. associated(cs%visc%tauy_shelf)) then
2961  !$OMP parallel do default(shared)
2962  do j=js-1,je ; do i=is,ie
2963  sfc_state%tauy_shelf(i,j) = us%R_to_kg_m3*us%L_T_to_m_s**2*us%Z_to_L*cs%visc%tauy_shelf(i,j)
2964  enddo ; enddo
2965  endif
2966 
2967  if (allocated(sfc_state%ocean_mass) .and. allocated(sfc_state%ocean_heat) .and. &
2968  allocated(sfc_state%ocean_salt)) then
2969  !$OMP parallel do default(shared)
2970  do j=js,je ; do i=is,ie
2971  sfc_state%ocean_mass(i,j) = 0.0
2972  sfc_state%ocean_heat(i,j) = 0.0 ; sfc_state%ocean_salt(i,j) = 0.0
2973  enddo ; enddo
2974  !$OMP parallel do default(shared) private(mass)
2975  do j=js,je ; do k=1,nz; do i=is,ie
2976  mass = gv%H_to_kg_m2*h(i,j,k)
2977  sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + mass
2978  sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*cs%tv%T(i,j,k)
2979  sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + &
2980  mass * (1.0e-3*cs%tv%S(i,j,k))
2981  enddo ; enddo ; enddo
2982  else
2983  if (allocated(sfc_state%ocean_mass)) then
2984  !$OMP parallel do default(shared)
2985  do j=js,je ; do i=is,ie ; sfc_state%ocean_mass(i,j) = 0.0 ; enddo ; enddo
2986  !$OMP parallel do default(shared)
2987  do j=js,je ; do k=1,nz ; do i=is,ie
2988  sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
2989  enddo ; enddo ; enddo
2990  endif
2991  if (allocated(sfc_state%ocean_heat)) then
2992  !$OMP parallel do default(shared)
2993  do j=js,je ; do i=is,ie ; sfc_state%ocean_heat(i,j) = 0.0 ; enddo ; enddo
2994  !$OMP parallel do default(shared) private(mass)
2995  do j=js,je ; do k=1,nz ; do i=is,ie
2996  mass = gv%H_to_kg_m2*h(i,j,k)
2997  sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*cs%tv%T(i,j,k)
2998  enddo ; enddo ; enddo
2999  endif
3000  if (allocated(sfc_state%ocean_salt)) then
3001  !$OMP parallel do default(shared)
3002  do j=js,je ; do i=is,ie ; sfc_state%ocean_salt(i,j) = 0.0 ; enddo ; enddo
3003  !$OMP parallel do default(shared) private(mass)
3004  do j=js,je ; do k=1,nz ; do i=is,ie
3005  mass = gv%H_to_kg_m2*h(i,j,k)
3006  sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + &
3007  mass * (1.0e-3*cs%tv%S(i,j,k))
3008  enddo ; enddo ; enddo
3009  endif
3010  endif
3011 
3012  if (associated(cs%tracer_flow_CSp)) then
3013  call call_tracer_surface_state(sfc_state, h, g, cs%tracer_flow_CSp)
3014  endif
3015 
3016  if (cs%check_bad_sfc_vals) then
3017  numberoferrors=0 ! count number of errors
3018  do j=js,je; do i=is,ie
3019  if (g%mask2dT(i,j)>0.) then
3020  bathy_m = cs%US%Z_to_m * g%bathyT(i,j)
3021  localerror = sfc_state%sea_lev(i,j)<=-bathy_m &
3022  .or. sfc_state%sea_lev(i,j)>= cs%bad_val_ssh_max &
3023  .or. sfc_state%sea_lev(i,j)<=-cs%bad_val_ssh_max &
3024  .or. sfc_state%sea_lev(i,j) + bathy_m < cs%bad_val_col_thick
3025  if (use_temperature) localerror = localerror &
3026  .or. sfc_state%SSS(i,j)<0. &
3027  .or. sfc_state%SSS(i,j)>=cs%bad_val_sss_max &
3028  .or. sfc_state%SST(i,j)< cs%bad_val_sst_min &
3029  .or. sfc_state%SST(i,j)>=cs%bad_val_sst_max
3030  if (localerror) then
3031  numberoferrors=numberoferrors+1
3032  if (numberoferrors<9) then ! Only report details for the first few errors
3033  ig = i + g%HI%idg_offset ! Global i-index
3034  jg = j + g%HI%jdg_offset ! Global j-index
3035  if (use_temperature) then
3036  write(msg(1:240),'(2(a,i4,x),4(a,f8.3,x),8(a,es11.4,x))') &
3037  'Extreme surface sfc_state detected: i=',ig,'j=',jg, &
3038  'lon=',g%geoLonT(i,j), 'lat=',g%geoLatT(i,j), &
3039  'x=',g%gridLonT(ig), 'y=',g%gridLatT(jg), &
3040  'D=',bathy_m, 'SSH=',sfc_state%sea_lev(i,j), &
3041  'SST=',sfc_state%SST(i,j), 'SSS=',sfc_state%SSS(i,j), &
3042  'U-=',sfc_state%u(i-1,j), 'U+=',sfc_state%u(i,j), &
3043  'V-=',sfc_state%v(i,j-1), 'V+=',sfc_state%v(i,j)
3044  else
3045  write(msg(1:240),'(2(a,i4,x),4(a,f8.3,x),6(a,es11.4))') &
3046  'Extreme surface sfc_state detected: i=',ig,'j=',jg, &
3047  'lon=',g%geoLonT(i,j), 'lat=',g%geoLatT(i,j), &
3048  'x=',g%gridLonT(i), 'y=',g%gridLatT(j), &
3049  'D=',bathy_m, 'SSH=',sfc_state%sea_lev(i,j), &
3050  'U-=',sfc_state%u(i-1,j), 'U+=',sfc_state%u(i,j), &
3051  'V-=',sfc_state%v(i,j-1), 'V+=',sfc_state%v(i,j)
3052  endif
3053  call mom_error(warning, trim(msg), all_print=.true.)
3054  elseif (numberoferrors==9) then ! Indicate once that there are more errors
3055  call mom_error(warning, 'There were more unreported extreme events!', all_print=.true.)
3056  endif ! numberOfErrors
3057  endif ! localError
3058  endif ! mask2dT
3059  enddo ; enddo
3060  call sum_across_pes(numberoferrors)
3061  if (numberoferrors>0) then
3062  write(msg(1:240),'(3(a,i9,x))') 'There were a total of ',numberoferrors, &
3063  'locations detected with extreme surface values!'
3064  call mom_error(fatal, trim(msg))
3065  endif
3066  endif
3067 
3068  if (cs%debug) call mom_surface_chksum("Post extract_sfc", sfc_state, g)
3069 
3070  call calltree_leave("extract_surface_sfc_state()")

References mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), and mom_checksum_packages::mom_surface_chksum().

Referenced by mom_main(), ocn_comp_mct::ocean_model_init_sfc(), step_mom(), and step_offline().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ finish_mom_initialization()

subroutine, public mom::finish_mom_initialization ( type(time_type), intent(in)  Time,
type(directories), intent(in)  dirs,
type(mom_control_struct), pointer  CS,
type(mom_restart_cs), pointer  restart_CSp 
)

Finishes initializing MOM and writes out the initial conditions.

Parameters
[in]timemodel time, used in this routine
[in]dirsstructure with directory paths
cspointer to MOM control structure
restart_csppointer to the restart control structure that will be used for MOM.

Definition at line 2484 of file MOM.F90.

2484  type(time_type), intent(in) :: Time !< model time, used in this routine
2485  type(directories), intent(in) :: dirs !< structure with directory paths
2486  type(MOM_control_struct), pointer :: CS !< pointer to MOM control structure
2487  type(MOM_restart_CS), pointer :: restart_CSp !< pointer to the restart control
2488  !! structure that will be used for MOM.
2489  ! Local variables
2490  type(ocean_grid_type), pointer :: G => null() ! pointer to a structure containing
2491  ! metrics and related information
2492  type(verticalGrid_type), pointer :: GV => null() ! Pointer to the vertical grid structure
2493  type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing
2494  ! various unit conversion factors
2495  type(MOM_restart_CS), pointer :: restart_CSp_tmp => null()
2496  real, allocatable :: z_interface(:,:,:) ! Interface heights [m]
2497  type(vardesc) :: vd
2498 
2499  call cpu_clock_begin(id_clock_init)
2500  call calltree_enter("finish_MOM_initialization()")
2501 
2502  ! Pointers for convenience
2503  g => cs%G ; gv => cs%GV ; us => cs%US
2504 
2505  !### Move to initialize_MOM?
2506  call fix_restart_scaling(gv)
2507  call fix_restart_unit_scaling(us)
2508 
2509  ! Write initial conditions
2510  if (cs%write_IC) then
2511  allocate(restart_csp_tmp)
2512  restart_csp_tmp = restart_csp
2513  allocate(z_interface(szi_(g),szj_(g),szk_(g)+1))
2514  call find_eta(cs%h, cs%tv, g, gv, us, z_interface, eta_to_m=1.0)
2515  call register_restart_field(z_interface, "eta", .true., restart_csp_tmp, &
2516  "Interface heights", "meter", z_grid='i')
2517 
2518  call save_restart(dirs%output_directory, time, g, &
2519  restart_csp_tmp, filename=cs%IC_file, gv=gv)
2520  deallocate(z_interface)
2521  deallocate(restart_csp_tmp)
2522  endif
2523 
2524  call write_energy(cs%u, cs%v, cs%h, cs%tv, time, 0, g, gv, us, &
2525  cs%sum_output_CSp, cs%tracer_flow_CSp)
2526 
2527  call calltree_leave("finish_MOM_initialization()")
2528  call cpu_clock_end(id_clock_init)
2529 

References mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), mom_unit_scaling::fix_restart_unit_scaling(), and id_clock_init.

Referenced by mom_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_mom_state_elements()

subroutine, public mom::get_mom_state_elements ( type(mom_control_struct), pointer  CS,
type(ocean_grid_type), optional, pointer  G,
type(verticalgrid_type), optional, pointer  GV,
type(unit_scale_type), optional, pointer  US,
real, intent(out), optional  C_p,
logical, intent(out), optional  use_temp 
)

This subroutine offers access to values or pointers to other types from within the MOM_control_struct, allowing the MOM_control_struct to be opaque.

Parameters
csMOM control structure
gstructure containing metrics and grid info
gvstructure containing vertical grid info
usA dimensional unit scaling type
[out]c_pThe heat capacity
[out]use_tempIndicates whether temperature is a state variable

Definition at line 3096 of file MOM.F90.

3096  type(MOM_control_struct), pointer :: CS !< MOM control structure
3097  type(ocean_grid_type), &
3098  optional, pointer :: G !< structure containing metrics and grid info
3099  type(verticalGrid_type), &
3100  optional, pointer :: GV !< structure containing vertical grid info
3101  type(unit_scale_type), &
3102  optional, pointer :: US !< A dimensional unit scaling type
3103  real, optional, intent(out) :: C_p !< The heat capacity
3104  logical, optional, intent(out) :: use_temp !< Indicates whether temperature is a state variable
3105 
3106  if (present(g)) g => cs%G
3107  if (present(gv)) gv => cs%GV
3108  if (present(us)) us => cs%US
3109  if (present(c_p)) c_p = cs%tv%C_p
3110  if (present(use_temp)) use_temp = associated(cs%tv%T)

Referenced by mom_main().

Here is the caller graph for this function:

◆ get_ocean_stocks()

subroutine, public mom::get_ocean_stocks ( type(mom_control_struct), pointer  CS,
real, intent(out), optional  mass,
real, intent(out), optional  heat,
real, intent(out), optional  salt,
logical, intent(in), optional  on_PE_only 
)

Find the global integrals of various quantities.

Parameters
csMOM control structure
[out]heatThe globally integrated integrated ocean heat [J].
[out]saltThe globally integrated integrated ocean salt [kg].
[out]massThe globally integrated integrated ocean mass [kg].
[in]on_pe_onlyIf present and true, only sum on the local PE.

Definition at line 3115 of file MOM.F90.

3115  type(MOM_control_struct), pointer :: CS !< MOM control structure
3116  real, optional, intent(out) :: heat !< The globally integrated integrated ocean heat [J].
3117  real, optional, intent(out) :: salt !< The globally integrated integrated ocean salt [kg].
3118  real, optional, intent(out) :: mass !< The globally integrated integrated ocean mass [kg].
3119  logical, optional, intent(in) :: on_PE_only !< If present and true, only sum on the local PE.
3120 
3121  if (present(mass)) &
3122  mass = global_mass_integral(cs%h, cs%G, cs%GV, on_pe_only=on_pe_only)
3123  if (present(heat)) &
3124  heat = cs%tv%C_p * global_mass_integral(cs%h, cs%G, cs%GV, cs%tv%T, on_pe_only=on_pe_only)
3125  if (present(salt)) &
3126  salt = 1.0e-3 * global_mass_integral(cs%h, cs%G, cs%GV, cs%tv%S, on_pe_only=on_pe_only)
3127 

References mom_spatial_means::global_mass_integral().

Referenced by mom_ocean_model_nuopc::ocean_stock_pe(), and mom_ocean_model_mct::ocean_stock_pe().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ initialize_mom()

subroutine, public mom::initialize_mom ( type(time_type), intent(inout), target  Time,
type(time_type), intent(in)  Time_init,
type(param_file_type), intent(out)  param_file,
type(directories), intent(out)  dirs,
type(mom_control_struct), pointer  CS,
type(mom_restart_cs), pointer  restart_CSp,
type(time_type), intent(in), optional  Time_in,
logical, intent(out), optional  offline_tracer_mode,
character(len=*), intent(in), optional  input_restart_file,
type(diag_ctrl), optional, pointer  diag_ptr,
logical, intent(in), optional  count_calls,
type(tracer_flow_control_cs), optional, pointer  tracer_flow_CSp 
)

Initialize MOM, including memory allocation, setting up parameters and diagnostics, initializing the ocean state variables, and initializing subsidiary modules.

Parameters
[in,out]timemodel time, set in this routine
[in]time_initThe start time for the coupled model's calendar
[out]param_filestructure indicating parameter file to parse
[out]dirsstructure with directory paths
cspointer set in this routine to MOM control structure
restart_csppointer set in this routine to the restart control structure that will be used for MOM.
[in]time_intime passed to MOM_initialize_state when model is not being started from a restart file
[out]offline_tracer_modeTrue is returned if tracers are being run offline
[in]input_restart_fileIf present, name of restart file to read
diag_ptrA pointer set in this routine to the diagnostic regulatory structure
tracer_flow_cspA pointer set in this routine to
[in]count_callsIf true, nstep_tot counts the number of calls to step_MOM instead of the number of dynamics timesteps.

Definition at line 1506 of file MOM.F90.

1506  type(time_type), target, intent(inout) :: Time !< model time, set in this routine
1507  type(time_type), intent(in) :: Time_init !< The start time for the coupled model's calendar
1508  type(param_file_type), intent(out) :: param_file !< structure indicating parameter file to parse
1509  type(directories), intent(out) :: dirs !< structure with directory paths
1510  type(MOM_control_struct), pointer :: CS !< pointer set in this routine to MOM control structure
1511  type(MOM_restart_CS), pointer :: restart_CSp !< pointer set in this routine to the
1512  !! restart control structure that will
1513  !! be used for MOM.
1514  type(time_type), optional, intent(in) :: Time_in !< time passed to MOM_initialize_state when
1515  !! model is not being started from a restart file
1516  logical, optional, intent(out) :: offline_tracer_mode !< True is returned if tracers are being run offline
1517  character(len=*),optional, intent(in) :: input_restart_file !< If present, name of restart file to read
1518  type(diag_ctrl), optional, pointer :: diag_ptr !< A pointer set in this routine to the diagnostic
1519  !! regulatory structure
1520  type(tracer_flow_control_CS), &
1521  optional, pointer :: tracer_flow_CSp !< A pointer set in this routine to
1522  !! the tracer flow control structure.
1523  logical, optional, intent(in) :: count_calls !< If true, nstep_tot counts the number of
1524  !! calls to step_MOM instead of the number of
1525  !! dynamics timesteps.
1526  ! local variables
1527  type(ocean_grid_type), pointer :: G => null() ! A pointer to a structure with metrics and related
1528  type(hor_index_type) :: HI ! A hor_index_type for array extents
1529  type(verticalGrid_type), pointer :: GV => null()
1530  type(dyn_horgrid_type), pointer :: dG => null()
1531  type(diag_ctrl), pointer :: diag => null()
1532  type(unit_scale_type), pointer :: US => null()
1533  character(len=4), parameter :: vers_num = 'v2.0'
1534 
1535  ! This include declares and sets the variable "version".
1536 # include "version_variable.h"
1537 
1538  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
1539  integer :: IsdB, IedB, JsdB, JedB
1540  real :: dtbt ! The barotropic timestep [s]
1541  real :: Z_diag_int ! minimum interval between calc depth-space diagnosetics [s]
1542 
1543  real, allocatable, dimension(:,:) :: eta ! free surface height or column mass [H ~> m or kg m-2]
1544  real, allocatable, dimension(:,:) :: area_shelf_h ! area occupied by ice shelf [m2]
1545  real, dimension(:,:), allocatable, target :: frac_shelf_h ! fraction of total area occupied by ice shelf [nondim]
1546  real, dimension(:,:), pointer :: shelf_area => null()
1547  type(MOM_restart_CS), pointer :: restart_CSp_tmp => null()
1548  type(group_pass_type) :: tmp_pass_uv_T_S_h, pass_uv_T_S_h
1549 
1550  real :: default_val ! default value for a parameter
1551  logical :: write_geom_files ! If true, write out the grid geometry files.
1552  logical :: ensemble_ocean ! If true, perform an ensemble gather at the end of step_MOM
1553  logical :: new_sim
1554  logical :: use_geothermal ! If true, apply geothermal heating.
1555  logical :: use_EOS ! If true, density calculated from T & S using an equation of state.
1556  logical :: symmetric ! If true, use symmetric memory allocation.
1557  logical :: save_IC ! If true, save the initial conditions.
1558  logical :: do_unit_tests ! If true, call unit tests.
1559  logical :: test_grid_copy = .false.
1560 
1561  logical :: bulkmixedlayer ! If true, a refined bulk mixed layer scheme is used
1562  ! with nkml sublayers and nkbl buffer layer.
1563  logical :: use_temperature ! If true, temp and saln used as state variables.
1564  logical :: use_frazil ! If true, liquid seawater freezes if temp below freezing,
1565  ! with accumulated heat deficit returned to surface ocean.
1566  logical :: bound_salinity ! If true, salt is added to keep salinity above
1567  ! a minimum value, and the deficit is reported.
1568  logical :: use_conT_absS ! If true, the prognostics T & S are conservative temperature
1569  ! and absolute salinity. Care should be taken to convert them
1570  ! to potential temperature and practical salinity before
1571  ! exchanging them with the coupler and/or reporting T&S diagnostics.
1572  logical :: advect_TS ! If false, then no horizontal advection of temperature
1573  ! and salnity is performed
1574  logical :: use_ice_shelf ! Needed for ALE
1575  logical :: global_indexing ! If true use global horizontal index values instead
1576  ! of having the data domain on each processor start at 1.
1577  logical :: bathy_at_vel ! If true, also define bathymetric fields at the
1578  ! the velocity points.
1579  logical :: calc_dtbt ! Indicates whether the dynamically adjusted barotropic
1580  ! time step needs to be updated before it is used.
1581  logical :: debug_truncations ! If true, turn on diagnostics useful for debugging truncations.
1582  integer :: first_direction ! An integer that indicates which direction is to be
1583  ! updated first in directionally split parts of the
1584  ! calculation. This can be altered during the course
1585  ! of the run via calls to set_first_direction.
1586  integer :: nkml, nkbl, verbosity, write_geom
1587  integer :: dynamics_stencil ! The computational stencil for the calculations
1588  ! in the dynamic core.
1589  real :: conv2watt, conv2salt
1590  character(len=48) :: flux_units, S_flux_units
1591 
1592  type(vardesc) :: vd_T, vd_S ! Structures describing temperature and salinity variables.
1593  type(time_type) :: Start_time
1594  type(ocean_internal_state) :: MOM_internal_state
1595  character(len=200) :: area_varname, ice_shelf_file, inputdir, filename
1596 
1597  if (associated(cs)) then
1598  call mom_error(warning, "initialize_MOM called with an associated "// &
1599  "control structure.")
1600  return
1601  endif
1602  allocate(cs)
1603 
1604  if (test_grid_copy) then ; allocate(g)
1605  else ; g => cs%G ; endif
1606 
1607  cs%Time => time
1608 
1609  id_clock_init = cpu_clock_id('Ocean Initialization', grain=clock_subcomponent)
1610  call cpu_clock_begin(id_clock_init)
1611 
1612  start_time = time ; if (present(time_in)) start_time = time_in
1613 
1614  ! Read paths and filenames from namelist and store in "dirs".
1615  ! Also open the parsed input parameter file(s) and setup param_file.
1616  call get_mom_input(param_file, dirs, default_input_filename=input_restart_file)
1617 
1618  verbosity = 2 ; call read_param(param_file, "VERBOSITY", verbosity)
1619  call mom_set_verbosity(verbosity)
1620  call calltree_enter("initialize_MOM(), MOM.F90")
1621 
1622  call find_obsolete_params(param_file)
1623 
1624  ! Read relevant parameters and write them to the model log.
1625  call log_version(param_file, "MOM", version, "")
1626  call get_param(param_file, "MOM", "VERBOSITY", verbosity, &
1627  "Integer controlling level of messaging\n" // &
1628  "\t0 = Only FATAL messages\n" // &
1629  "\t2 = Only FATAL, WARNING, NOTE [default]\n" // &
1630  "\t9 = All)", default=2, debuggingparam=.true.)
1631  call get_param(param_file, "MOM", "DO_UNIT_TESTS", do_unit_tests, &
1632  "If True, exercises unit tests at model start up.", &
1633  default=.false., debuggingparam=.true.)
1634  if (do_unit_tests) then
1635  call unit_tests(verbosity)
1636  endif
1637 
1638  ! Determining the internal unit scaling factors for this run.
1639  call unit_scaling_init(param_file, cs%US)
1640 
1641  us => cs%US
1642 
1643  call get_param(param_file, "MOM", "SPLIT", cs%split, &
1644  "Use the split time stepping if true.", default=.true.)
1645  if (cs%split) then
1646  cs%use_RK2 = .false.
1647  else
1648  call get_param(param_file, "MOM", "USE_RK2", cs%use_RK2, &
1649  "If true, use RK2 instead of RK3 in the unsplit time stepping.", &
1650  default=.false.)
1651  endif
1652 
1653  call get_param(param_file, "MOM", "CALC_RHO_FOR_SEA_LEVEL", cs%calc_rho_for_sea_lev, &
1654  "If true, the in-situ density is used to calculate the "//&
1655  "effective sea level that is returned to the coupler. If false, "//&
1656  "the Boussinesq parameter RHO_0 is used.", default=.false.)
1657  call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, &
1658  "If true, Temperature and salinity are used as state "//&
1659  "variables.", default=.true.)
1660  call get_param(param_file, "MOM", "USE_EOS", use_eos, &
1661  "If true, density is calculated from temperature and "//&
1662  "salinity with an equation of state. If USE_EOS is "//&
1663  "true, ENABLE_THERMODYNAMICS must be true as well.", &
1664  default=use_temperature)
1665  call get_param(param_file, "MOM", "DIABATIC_FIRST", cs%diabatic_first, &
1666  "If true, apply diabatic and thermodynamic processes, "//&
1667  "including buoyancy forcing and mass gain or loss, "//&
1668  "before stepping the dynamics forward.", default=.false.)
1669  call get_param(param_file, "MOM", "USE_CONTEMP_ABSSAL", use_cont_abss, &
1670  "If true, the prognostics T&S are the conservative temperature "//&
1671  "and absolute salinity. Care should be taken to convert them "//&
1672  "to potential temperature and practical salinity before "//&
1673  "exchanging them with the coupler and/or reporting T&S diagnostics.", &
1674  default=.false.)
1675  cs%tv%T_is_conT = use_cont_abss ; cs%tv%S_is_absS = use_cont_abss
1676  call get_param(param_file, "MOM", "ADIABATIC", cs%adiabatic, &
1677  "There are no diapycnal mass fluxes if ADIABATIC is "//&
1678  "true. This assumes that KD = KDML = 0.0 and that "//&
1679  "there is no buoyancy forcing, but makes the model "//&
1680  "faster by eliminating subroutine calls.", default=.false.)
1681  call get_param(param_file, "MOM", "DO_DYNAMICS", cs%do_dynamics, &
1682  "If False, skips the dynamics calls that update u & v, as well as "//&
1683  "the gravity wave adjustment to h. This may be a fragile feature, "//&
1684  "but can be useful during development", default=.true.)
1685  call get_param(param_file, "MOM", "ADVECT_TS", advect_ts, &
1686  "If True, advect temperature and salinity horizontally "//&
1687  "If False, T/S are registered for advection. "//&
1688  "This is intended only to be used in offline tracer mode "//&
1689  "and is by default false in that case.", &
1690  do_not_log = .true., default=.true. )
1691  if (present(offline_tracer_mode)) then ! Only read this parameter in enabled modes
1692  call get_param(param_file, "MOM", "OFFLINE_TRACER_MODE", cs%offline_tracer_mode, &
1693  "If true, barotropic and baroclinic dynamics, thermodynamics "//&
1694  "are all bypassed with all the fields necessary to integrate "//&
1695  "the tracer advection and diffusion equation are read in from "//&
1696  "files stored from a previous integration of the prognostic model. "//&
1697  "NOTE: This option only used in the ocean_solo_driver.", default=.false.)
1698  if (cs%offline_tracer_mode) then
1699  call get_param(param_file, "MOM", "ADVECT_TS", advect_ts, &
1700  "If True, advect temperature and salinity horizontally "//&
1701  "If False, T/S are registered for advection. "//&
1702  "This is intended only to be used in offline tracer mode."//&
1703  "and is by default false in that case", &
1704  default=.false. )
1705  endif
1706  endif
1707  call get_param(param_file, "MOM", "USE_REGRIDDING", cs%use_ALE_algorithm, &
1708  "If True, use the ALE algorithm (regridding/remapping). "//&
1709  "If False, use the layered isopycnal algorithm.", default=.false. )
1710  call get_param(param_file, "MOM", "BULKMIXEDLAYER", bulkmixedlayer, &
1711  "If true, use a Kraus-Turner-like bulk mixed layer "//&
1712  "with transitional buffer layers. Layers 1 through "//&
1713  "NKML+NKBL have variable densities. There must be at "//&
1714  "least NKML+NKBL+1 layers if BULKMIXEDLAYER is true. "//&
1715  "BULKMIXEDLAYER can not be used with USE_REGRIDDING. "//&
1716  "The default is influenced by ENABLE_THERMODYNAMICS.", &
1717  default=use_temperature .and. .not.cs%use_ALE_algorithm)
1718  call get_param(param_file, "MOM", "THICKNESSDIFFUSE", cs%thickness_diffuse, &
1719  "If true, interface heights are diffused with a "//&
1720  "coefficient of KHTH.", default=.false.)
1721  call get_param(param_file, "MOM", "THICKNESSDIFFUSE_FIRST", &
1722  cs%thickness_diffuse_first, &
1723  "If true, do thickness diffusion before dynamics. "//&
1724  "This is only used if THICKNESSDIFFUSE is true.", &
1725  default=.false.)
1726  if (.not.cs%thickness_diffuse) cs%thickness_diffuse_first = .false.
1727  call get_param(param_file, "MOM", "BATHYMETRY_AT_VEL", bathy_at_vel, &
1728  "If true, there are separate values for the basin depths "//&
1729  "at velocity points. Otherwise the effects of topography "//&
1730  "are entirely determined from thickness points.", &
1731  default=.false.)
1732  call get_param(param_file, "MOM", "USE_WAVES", cs%UseWaves, default=.false., &
1733  do_not_log=.true.)
1734 
1735  call get_param(param_file, "MOM", "DEBUG", cs%debug, &
1736  "If true, write out verbose debugging data.", &
1737  default=.false., debuggingparam=.true.)
1738  call get_param(param_file, "MOM", "DEBUG_TRUNCATIONS", debug_truncations, &
1739  "If true, calculate all diagnostics that are useful for "//&
1740  "debugging truncations.", default=.false., debuggingparam=.true.)
1741 
1742  call get_param(param_file, "MOM", "DT", cs%dt, &
1743  "The (baroclinic) dynamics time step. The time-step that "//&
1744  "is actually used will be an integer fraction of the "//&
1745  "forcing time-step (DT_FORCING in ocean-only mode or the "//&
1746  "coupling timestep in coupled mode.)", units="s", scale=us%s_to_T, &
1747  fail_if_missing=.true.)
1748  call get_param(param_file, "MOM", "DT_THERM", cs%dt_therm, &
1749  "The thermodynamic and tracer advection time step. "//&
1750  "Ideally DT_THERM should be an integer multiple of DT "//&
1751  "and less than the forcing or coupling time-step, unless "//&
1752  "THERMO_SPANS_COUPLING is true, in which case DT_THERM "//&
1753  "can be an integer multiple of the coupling timestep. By "//&
1754  "default DT_THERM is set to DT.", &
1755  units="s", scale=us%s_to_T, default=us%T_to_s*cs%dt)
1756  call get_param(param_file, "MOM", "THERMO_SPANS_COUPLING", cs%thermo_spans_coupling, &
1757  "If true, the MOM will take thermodynamic and tracer "//&
1758  "timesteps that can be longer than the coupling timestep. "//&
1759  "The actual thermodynamic timestep that is used in this "//&
1760  "case is the largest integer multiple of the coupling "//&
1761  "timestep that is less than or equal to DT_THERM.", default=.false.)
1762 
1763  if (bulkmixedlayer) then
1764  cs%Hmix = -1.0 ; cs%Hmix_UV = -1.0
1765  else
1766  call get_param(param_file, "MOM", "HMIX_SFC_PROP", cs%Hmix, &
1767  "If BULKMIXEDLAYER is false, HMIX_SFC_PROP is the depth "//&
1768  "over which to average to find surface properties like "//&
1769  "SST and SSS or density (but not surface velocities).", &
1770  units="m", default=1.0, scale=us%m_to_Z)
1771  call get_param(param_file, "MOM", "HMIX_UV_SFC_PROP", cs%Hmix_UV, &
1772  "If BULKMIXEDLAYER is false, HMIX_UV_SFC_PROP is the depth "//&
1773  "over which to average to find surface flow properties, "//&
1774  "SSU, SSV. A non-positive value indicates no averaging.", &
1775  units="m", default=0.0, scale=us%m_to_Z)
1776  endif
1777  call get_param(param_file, "MOM", "HFREEZE", cs%HFrz, &
1778  "If HFREEZE > 0, melt potential will be computed. The actual depth "//&
1779  "over which melt potential is computed will be min(HFREEZE, OBLD), "//&
1780  "where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), "//&
1781  "melt potential will not be computed.", units="m", default=-1.0)
1782  call get_param(param_file, "MOM", "INTERPOLATE_P_SURF", cs%interp_p_surf, &
1783  "If true, linearly interpolate the surface pressure "//&
1784  "over the coupling time step, using the specified value "//&
1785  "at the end of the step.", default=.false.)
1786 
1787  if (cs%split) then
1788  call get_param(param_file, "MOM", "DTBT", dtbt, default=-0.98)
1789  default_val = us%T_to_s*cs%dt_therm ; if (dtbt > 0.0) default_val = -1.0
1790  cs%dtbt_reset_period = -1.0
1791  call get_param(param_file, "MOM", "DTBT_RESET_PERIOD", cs%dtbt_reset_period, &
1792  "The period between recalculations of DTBT (if DTBT <= 0). "//&
1793  "If DTBT_RESET_PERIOD is negative, DTBT is set based "//&
1794  "only on information available at initialization. If 0, "//&
1795  "DTBT will be set every dynamics time step. The default "//&
1796  "is set by DT_THERM. This is only used if SPLIT is true.", &
1797  units="s", default=default_val, do_not_read=(dtbt > 0.0))
1798  endif
1799 
1800  ! This is here in case these values are used inappropriately.
1801  use_frazil = .false. ; bound_salinity = .false. ; cs%tv%P_Ref = 2.0e7
1802  if (use_temperature) then
1803  call get_param(param_file, "MOM", "FRAZIL", use_frazil, &
1804  "If true, water freezes if it gets too cold, and the "//&
1805  "the accumulated heat deficit is returned in the "//&
1806  "surface state. FRAZIL is only used if "//&
1807  "ENABLE_THERMODYNAMICS is true.", default=.false.)
1808  call get_param(param_file, "MOM", "DO_GEOTHERMAL", use_geothermal, &
1809  "If true, apply geothermal heating.", default=.false.)
1810  call get_param(param_file, "MOM", "BOUND_SALINITY", bound_salinity, &
1811  "If true, limit salinity to being positive. (The sea-ice "//&
1812  "model may ask for more salt than is available and "//&
1813  "drive the salinity negative otherwise.)", default=.false.)
1814  call get_param(param_file, "MOM", "MIN_SALINITY", cs%tv%min_salinity, &
1815  "The minimum value of salinity when BOUND_SALINITY=True. "//&
1816  "The default is 0.01 for backward compatibility but ideally "//&
1817  "should be 0.", units="PPT", default=0.01, do_not_log=.not.bound_salinity)
1818  call get_param(param_file, "MOM", "C_P", cs%tv%C_p, &
1819  "The heat capacity of sea water, approximated as a "//&
1820  "constant. This is only used if ENABLE_THERMODYNAMICS is "//&
1821  "true. The default value is from the TEOS-10 definition "//&
1822  "of conservative temperature.", units="J kg-1 K-1", &
1823  default=3991.86795711963)
1824  endif
1825  if (use_eos) call get_param(param_file, "MOM", "P_REF", cs%tv%P_Ref, &
1826  "The pressure that is used for calculating the coordinate "//&
1827  "density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
1828  "This is only used if USE_EOS and ENABLE_THERMODYNAMICS "//&
1829  "are true.", units="Pa", default=2.0e7)
1830 
1831  if (bulkmixedlayer) then
1832  call get_param(param_file, "MOM", "NKML", nkml, &
1833  "The number of sublayers within the mixed layer if "//&
1834  "BULKMIXEDLAYER is true.", units="nondim", default=2)
1835  call get_param(param_file, "MOM", "NKBL", nkbl, &
1836  "The number of layers that are used as variable density "//&
1837  "buffer layers if BULKMIXEDLAYER is true.", units="nondim", &
1838  default=2)
1839  endif
1840 
1841  call get_param(param_file, "MOM", "GLOBAL_INDEXING", global_indexing, &
1842  "If true, use a global lateral indexing convention, so "//&
1843  "that corresponding points on different processors have "//&
1844  "the same index. This does not work with static memory.", &
1845  default=.false., layoutparam=.true.)
1846 #ifdef STATIC_MEMORY_
1847  if (global_indexing) call mom_error(fatal, "initialize_MOM: "//&
1848  "GLOBAL_INDEXING can not be true with STATIC_MEMORY.")
1849 #endif
1850  call get_param(param_file, "MOM", "FIRST_DIRECTION", first_direction, &
1851  "An integer that indicates which direction goes first "//&
1852  "in parts of the code that use directionally split "//&
1853  "updates, with even numbers (or 0) used for x- first "//&
1854  "and odd numbers used for y-first.", default=0)
1855 
1856  call get_param(param_file, "MOM", "CHECK_BAD_SURFACE_VALS", cs%check_bad_sfc_vals, &
1857  "If true, check the surface state for ridiculous values.", &
1858  default=.false.)
1859  if (cs%check_bad_sfc_vals) then
1860  call get_param(param_file, "MOM", "BAD_VAL_SSH_MAX", cs%bad_val_ssh_max, &
1861  "The value of SSH above which a bad value message is "//&
1862  "triggered, if CHECK_BAD_SURFACE_VALS is true.", units="m", &
1863  default=20.0)
1864  call get_param(param_file, "MOM", "BAD_VAL_SSS_MAX", cs%bad_val_sss_max, &
1865  "The value of SSS above which a bad value message is "//&
1866  "triggered, if CHECK_BAD_SURFACE_VALS is true.", units="PPT", &
1867  default=45.0)
1868  call get_param(param_file, "MOM", "BAD_VAL_SST_MAX", cs%bad_val_sst_max, &
1869  "The value of SST above which a bad value message is "//&
1870  "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
1871  units="deg C", default=45.0)
1872  call get_param(param_file, "MOM", "BAD_VAL_SST_MIN", cs%bad_val_sst_min, &
1873  "The value of SST below which a bad value message is "//&
1874  "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
1875  units="deg C", default=-2.1)
1876  call get_param(param_file, "MOM", "BAD_VAL_COLUMN_THICKNESS", cs%bad_val_col_thick, &
1877  "The value of column thickness below which a bad value message is "//&
1878  "triggered, if CHECK_BAD_SURFACE_VALS is true.", units="m", &
1879  default=0.0)
1880  endif
1881 
1882  call get_param(param_file, "MOM", "SAVE_INITIAL_CONDS", save_ic, &
1883  "If true, write the initial conditions to a file given "//&
1884  "by IC_OUTPUT_FILE.", default=.false.)
1885  call get_param(param_file, "MOM", "IC_OUTPUT_FILE", cs%IC_file, &
1886  "The file into which to write the initial conditions.", &
1887  default="MOM_IC")
1888  call get_param(param_file, "MOM", "WRITE_GEOM", write_geom, &
1889  "If =0, never write the geometry and vertical grid files. "//&
1890  "If =1, write the geometry and vertical grid files only for "//&
1891  "a new simulation. If =2, always write the geometry and "//&
1892  "vertical grid files. Other values are invalid.", default=1)
1893  if (write_geom<0 .or. write_geom>2) call mom_error(fatal,"MOM: "//&
1894  "WRITE_GEOM must be equal to 0, 1 or 2.")
1895  write_geom_files = ((write_geom==2) .or. ((write_geom==1) .and. &
1896  ((dirs%input_filename(1:1)=='n') .and. (len_trim(dirs%input_filename)==1))))
1897 ! If the restart file type had been initialized, this could become:
1898 ! write_geom_files = ((write_geom==2) .or. &
1899 ! ((write_geom==1) .and. is_new_run(restart_CSp)))
1900 
1901  ! Check for inconsistent parameter settings.
1902  if (cs%use_ALE_algorithm .and. bulkmixedlayer) call mom_error(fatal, &
1903  "MOM: BULKMIXEDLAYER can not currently be used with the ALE algorithm.")
1904  if (cs%use_ALE_algorithm .and. .not.use_temperature) call mom_error(fatal, &
1905  "MOM: At this time, USE_EOS should be True when using the ALE algorithm.")
1906  if (cs%adiabatic .and. use_temperature) call mom_error(warning, &
1907  "MOM: ADIABATIC and ENABLE_THERMODYNAMICS both defined is usually unwise.")
1908  if (use_eos .and. .not.use_temperature) call mom_error(fatal, &
1909  "MOM: ENABLE_THERMODYNAMICS must be defined to use USE_EOS.")
1910  if (cs%adiabatic .and. bulkmixedlayer) call mom_error(fatal, &
1911  "MOM: ADIABATIC and BULKMIXEDLAYER can not both be defined.")
1912  if (bulkmixedlayer .and. .not.use_eos) call mom_error(fatal, &
1913  "initialize_MOM: A bulk mixed layer can only be used with T & S as "//&
1914  "state variables. Add USE_EOS = True to MOM_input.")
1915 
1916  call get_param(param_file, 'MOM', "ICE_SHELF", use_ice_shelf, default=.false., do_not_log=.true.)
1917  if (use_ice_shelf) then
1918  inputdir = "." ; call get_param(param_file, 'MOM', "INPUTDIR", inputdir)
1919  inputdir = slasher(inputdir)
1920  call get_param(param_file, 'MOM', "ICE_THICKNESS_FILE", ice_shelf_file, &
1921  "The file from which the ice bathymetry and area are read.", &
1922  fail_if_missing=.true.)
1923  call get_param(param_file, 'MOM', "ICE_AREA_VARNAME", area_varname, &
1924  "The name of the area variable in ICE_THICKNESS_FILE.", &
1925  fail_if_missing=.true.)
1926  endif
1927 
1928 
1929  cs%ensemble_ocean=.false.
1930  call get_param(param_file, "MOM", "ENSEMBLE_OCEAN", cs%ensemble_ocean, &
1931  "If False, The model is being run in serial mode as a single realization. "//&
1932  "If True, The current model realization is part of a larger ensemble "//&
1933  "and at the end of step MOM, we will perform a gather of the ensemble "//&
1934  "members for statistical evaluation and/or data assimilation.", default=.false.)
1935 
1936  call calltree_waypoint("MOM parameters read (initialize_MOM)")
1937 
1938  ! Set up the model domain and grids.
1939 #ifdef SYMMETRIC_MEMORY_
1940  symmetric = .true.
1941 #else
1942  symmetric = .false.
1943 #endif
1944 #ifdef STATIC_MEMORY_
1945  call mom_domains_init(g%domain, param_file, symmetric=symmetric, &
1946  static_memory=.true., nihalo=nihalo_, njhalo=njhalo_, &
1947  niglobal=niglobal_, njglobal=njglobal_, niproc=niproc_, &
1948  njproc=njproc_)
1949 #else
1950  call mom_domains_init(g%domain, param_file, symmetric=symmetric)
1951 #endif
1952  call calltree_waypoint("domains initialized (initialize_MOM)")
1953 
1954  call mom_debugging_init(param_file)
1955  call diag_mediator_infrastructure_init()
1956  call mom_io_init(param_file)
1957 
1958  call hor_index_init(g%Domain, hi, param_file, &
1959  local_indexing=.not.global_indexing)
1960 
1961  call create_dyn_horgrid(dg, hi, bathymetry_at_vel=bathy_at_vel)
1962  call clone_mom_domain(g%Domain, dg%Domain)
1963 
1964  call verticalgridinit( param_file, cs%GV, us )
1965  gv => cs%GV
1966 ! dG%g_Earth = GV%mks_g_Earth
1967 
1968  ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
1969  if (cs%debug .or. dg%symmetric) &
1970  call clone_mom_domain(dg%Domain, dg%Domain_aux, symmetric=.false.)
1971 
1972  call calltree_waypoint("grids initialized (initialize_MOM)")
1973 
1974  call mom_timing_init(cs)
1975 
1976  ! Allocate initialize time-invariant MOM variables.
1977  call mom_initialize_fixed(dg, us, cs%OBC, param_file, write_geom_files, dirs%output_directory)
1978  call calltree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)")
1979 
1980  if (associated(cs%OBC)) call call_obc_register(param_file, cs%update_OBC_CSp, cs%OBC)
1981 
1982  call tracer_registry_init(param_file, cs%tracer_Reg)
1983 
1984  ! Allocate and initialize space for the primary time-varying MOM variables.
1985  is = dg%isc ; ie = dg%iec ; js = dg%jsc ; je = dg%jec ; nz = gv%ke
1986  isd = dg%isd ; ied = dg%ied ; jsd = dg%jsd ; jed = dg%jed
1987  isdb = dg%IsdB ; iedb = dg%IedB ; jsdb = dg%JsdB ; jedb = dg%JedB
1988  alloc_(cs%u(isdb:iedb,jsd:jed,nz)) ; cs%u(:,:,:) = 0.0
1989  alloc_(cs%v(isd:ied,jsdb:jedb,nz)) ; cs%v(:,:,:) = 0.0
1990  alloc_(cs%h(isd:ied,jsd:jed,nz)) ; cs%h(:,:,:) = gv%Angstrom_H
1991  alloc_(cs%uh(isdb:iedb,jsd:jed,nz)) ; cs%uh(:,:,:) = 0.0
1992  alloc_(cs%vh(isd:ied,jsdb:jedb,nz)) ; cs%vh(:,:,:) = 0.0
1993  if (use_temperature) then
1994  alloc_(cs%T(isd:ied,jsd:jed,nz)) ; cs%T(:,:,:) = 0.0
1995  alloc_(cs%S(isd:ied,jsd:jed,nz)) ; cs%S(:,:,:) = 0.0
1996  cs%tv%T => cs%T ; cs%tv%S => cs%S
1997  if (cs%tv%T_is_conT) then
1998  vd_t = var_desc(name="contemp", units="Celsius", longname="Conservative Temperature", &
1999  cmor_field_name="thetao", cmor_longname="Sea Water Potential Temperature", &
2000  conversion=cs%tv%C_p)
2001  else
2002  vd_t = var_desc(name="temp", units="degC", longname="Potential Temperature", &
2003  cmor_field_name="thetao", cmor_longname="Sea Water Potential Temperature", &
2004  conversion=cs%tv%C_p)
2005  endif
2006  if (cs%tv%S_is_absS) then
2007  vd_s = var_desc(name="abssalt",units="g kg-1",longname="Absolute Salinity", &
2008  cmor_field_name="so", cmor_longname="Sea Water Salinity", &
2009  conversion=0.001)
2010  else
2011  vd_s = var_desc(name="salt",units="psu",longname="Salinity", &
2012  cmor_field_name="so", cmor_longname="Sea Water Salinity", &
2013  conversion=0.001)
2014  endif
2015 
2016  if (advect_ts) then
2017  s_flux_units = get_tr_flux_units(gv, "psu") ! Could change to "kg m-2 s-1"?
2018  conv2watt = gv%H_to_kg_m2 * cs%tv%C_p
2019  if (gv%Boussinesq) then
2020  conv2salt = gv%H_to_m ! Could change to GV%H_to_kg_m2 * 0.001?
2021  else
2022  conv2salt = gv%H_to_kg_m2
2023  endif
2024  call register_tracer(cs%tv%T, cs%tracer_Reg, param_file, dg%HI, gv, &
2025  tr_desc=vd_t, registry_diags=.true., flux_nameroot='T', &
2026  flux_units='W', flux_longname='Heat', &
2027  flux_scale=conv2watt, convergence_units='W m-2', &
2028  convergence_scale=conv2watt, cmor_tendprefix="opottemp", diag_form=2)
2029  call register_tracer(cs%tv%S, cs%tracer_Reg, param_file, dg%HI, gv, &
2030  tr_desc=vd_s, registry_diags=.true., flux_nameroot='S', &
2031  flux_units=s_flux_units, flux_longname='Salt', &
2032  flux_scale=conv2salt, convergence_units='kg m-2 s-1', &
2033  convergence_scale=0.001*gv%H_to_kg_m2, cmor_tendprefix="osalt", diag_form=2)
2034  endif
2035  if (associated(cs%OBC)) &
2036  call register_temp_salt_segments(gv, cs%OBC, cs%tracer_Reg, param_file)
2037  endif
2038  if (use_frazil) then
2039  allocate(cs%tv%frazil(isd:ied,jsd:jed)) ; cs%tv%frazil(:,:) = 0.0
2040  endif
2041  if (bound_salinity) then
2042  allocate(cs%tv%salt_deficit(isd:ied,jsd:jed)) ; cs%tv%salt_deficit(:,:)=0.0
2043  endif
2044 
2045  if (bulkmixedlayer .or. use_temperature) then
2046  allocate(cs%Hml(isd:ied,jsd:jed)) ; cs%Hml(:,:) = 0.0
2047  endif
2048 
2049  if (bulkmixedlayer) then
2050  gv%nkml = nkml ; gv%nk_rho_varies = nkml + nkbl
2051  else
2052  gv%nkml = 0 ; gv%nk_rho_varies = 0
2053  endif
2054  if (cs%use_ALE_algorithm) then
2055  call get_param(param_file, "MOM", "NK_RHO_VARIES", gv%nk_rho_varies, default=0) ! Will default to nz later... -AJA
2056  endif
2057 
2058  alloc_(cs%uhtr(isdb:iedb,jsd:jed,nz)) ; cs%uhtr(:,:,:) = 0.0
2059  alloc_(cs%vhtr(isd:ied,jsdb:jedb,nz)) ; cs%vhtr(:,:,:) = 0.0
2060  cs%t_dyn_rel_adv = 0.0 ; cs%t_dyn_rel_thermo = 0.0 ; cs%t_dyn_rel_diag = 0.0
2061 
2062  if (debug_truncations) then
2063  allocate(cs%u_prev(isdb:iedb,jsd:jed,nz)) ; cs%u_prev(:,:,:) = 0.0
2064  allocate(cs%v_prev(isd:ied,jsdb:jedb,nz)) ; cs%v_prev(:,:,:) = 0.0
2065  mom_internal_state%u_prev => cs%u_prev
2066  mom_internal_state%v_prev => cs%v_prev
2067  call safe_alloc_ptr(cs%ADp%du_dt_visc,isdb,iedb,jsd,jed,nz)
2068  call safe_alloc_ptr(cs%ADp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
2069  if (.not.cs%adiabatic) then
2070  call safe_alloc_ptr(cs%ADp%du_dt_dia,isdb,iedb,jsd,jed,nz)
2071  call safe_alloc_ptr(cs%ADp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
2072  endif
2073  endif
2074 
2075  mom_internal_state%u => cs%u ; mom_internal_state%v => cs%v
2076  mom_internal_state%h => cs%h
2077  mom_internal_state%uh => cs%uh ; mom_internal_state%vh => cs%vh
2078  if (use_temperature) then
2079  mom_internal_state%T => cs%T ; mom_internal_state%S => cs%S
2080  endif
2081 
2082  cs%CDp%uh => cs%uh ; cs%CDp%vh => cs%vh
2083 
2084  if (cs%interp_p_surf) then
2085  allocate(cs%p_surf_prev(isd:ied,jsd:jed)) ; cs%p_surf_prev(:,:) = 0.0
2086  endif
2087 
2088  alloc_(cs%ssh_rint(isd:ied,jsd:jed)) ; cs%ssh_rint(:,:) = 0.0
2089  alloc_(cs%ave_ssh_ibc(isd:ied,jsd:jed)) ; cs%ave_ssh_ibc(:,:) = 0.0
2090  alloc_(cs%eta_av_bc(isd:ied,jsd:jed)) ; cs%eta_av_bc(:,:) = 0.0
2091  cs%time_in_cycle = 0.0 ; cs%time_in_thermo_cycle = 0.0
2092 
2093  ! Use the Wright equation of state by default, unless otherwise specified
2094  ! Note: this line and the following block ought to be in a separate
2095  ! initialization routine for tv.
2096  if (use_eos) call eos_init(param_file, cs%tv%eqn_of_state)
2097  if (use_temperature) then
2098  allocate(cs%tv%TempxPmE(isd:ied,jsd:jed))
2099  cs%tv%TempxPmE(:,:) = 0.0
2100  if (use_geothermal) then
2101  allocate(cs%tv%internal_heat(isd:ied,jsd:jed))
2102  cs%tv%internal_heat(:,:) = 0.0
2103  endif
2104  endif
2105  call calltree_waypoint("state variables allocated (initialize_MOM)")
2106 
2107  ! Set the fields that are needed for bitwise identical restarting
2108  ! the time stepping scheme.
2109  call restart_init(param_file, restart_csp)
2110  call set_restart_fields(gv, us, param_file, cs, restart_csp)
2111  if (cs%split) then
2112  call register_restarts_dyn_split_rk2(dg%HI, gv, param_file, &
2113  cs%dyn_split_RK2_CSp, restart_csp, cs%uh, cs%vh)
2114  elseif (cs%use_RK2) then
2115  call register_restarts_dyn_unsplit_rk2(dg%HI, gv, param_file, &
2116  cs%dyn_unsplit_RK2_CSp, restart_csp)
2117  else
2118  call register_restarts_dyn_unsplit(dg%HI, gv, param_file, &
2119  cs%dyn_unsplit_CSp, restart_csp)
2120  endif
2121 
2122  ! This subroutine calls user-specified tracer registration routines.
2123  ! Additional calls can be added to MOM_tracer_flow_control.F90.
2124  call call_tracer_register(dg%HI, gv, us, param_file, cs%tracer_flow_CSp, &
2125  cs%tracer_Reg, restart_csp)
2126 
2127  call meke_alloc_register_restart(dg%HI, param_file, cs%MEKE, restart_csp)
2128  call set_visc_register_restarts(dg%HI, gv, param_file, cs%visc, restart_csp)
2129  call mixedlayer_restrat_register_restarts(dg%HI, param_file, &
2130  cs%mixedlayer_restrat_CSp, restart_csp)
2131 
2132  if (associated(cs%OBC)) &
2133  call open_boundary_register_restarts(dg%HI, gv, cs%OBC, cs%tracer_Reg, &
2134  param_file, restart_csp, use_temperature)
2135 
2136  call calltree_waypoint("restart registration complete (initialize_MOM)")
2137 
2138  ! Initialize dynamically evolving fields, perhaps from restart files.
2139  call cpu_clock_begin(id_clock_mom_init)
2140  call mom_initialize_coord(gv, us, param_file, write_geom_files, &
2141  dirs%output_directory, cs%tv, dg%max_depth)
2142  call calltree_waypoint("returned from MOM_initialize_coord() (initialize_MOM)")
2143 
2144  if (cs%use_ALE_algorithm) then
2145  call ale_init(param_file, gv, us, dg%max_depth, cs%ALE_CSp)
2146  call calltree_waypoint("returned from ALE_init() (initialize_MOM)")
2147  endif
2148 
2149  ! Shift from using the temporary dynamic grid type to using the final
2150  ! (potentially static) ocean-specific grid type.
2151  ! The next line would be needed if G%Domain had not already been init'd above:
2152  ! call clone_MOM_domain(dG%Domain, G%Domain)
2153  call mom_grid_init(g, param_file, us, hi, bathymetry_at_vel=bathy_at_vel)
2154  call copy_dyngrid_to_mom_grid(dg, g, us)
2155  call destroy_dyn_horgrid(dg)
2156 
2157  ! Set a few remaining fields that are specific to the ocean grid type.
2158  call set_first_direction(g, first_direction)
2159  ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
2160  if (cs%debug .or. g%symmetric) then
2161  call clone_mom_domain(g%Domain, g%Domain_aux, symmetric=.false.)
2162  else ; g%Domain_aux => g%Domain ; endif
2163  ! Copy common variables from the vertical grid to the horizontal grid.
2164  ! Consider removing this later?
2165  g%ke = gv%ke ; g%g_Earth = gv%mks_g_Earth
2166 
2167  call mom_initialize_state(cs%u, cs%v, cs%h, cs%tv, time, g, gv, us, param_file, &
2168  dirs, restart_csp, cs%ALE_CSp, cs%tracer_Reg, &
2169  cs%sponge_CSp, cs%ALE_sponge_CSp, cs%OBC, time_in)
2170  call cpu_clock_end(id_clock_mom_init)
2171  call calltree_waypoint("returned from MOM_initialize_state() (initialize_MOM)")
2172 
2173 ! ! Need this after MOM_initialize_state for DOME OBC stuff.
2174 ! if (associated(CS%OBC)) &
2175 ! call open_boundary_register_restarts(G%HI, GV, CS%OBC, CS%tracer_Reg, &
2176 ! param_file, restart_CSp, use_temperature)
2177 
2178 ! call callTree_waypoint("restart registration complete (initialize_MOM)")
2179 
2180  ! From this point, there may be pointers being set, so the final grid type
2181  ! that will persist throughout the run has to be used.
2182 
2183  if (test_grid_copy) then
2184  ! Copy the data from the temporary grid to the dyn_hor_grid to CS%G.
2185  call create_dyn_horgrid(dg, g%HI)
2186  call clone_mom_domain(g%Domain, dg%Domain)
2187 
2188  call clone_mom_domain(g%Domain, cs%G%Domain)
2189  call mom_grid_init(cs%G, param_file, us)
2190 
2191  call copy_mom_grid_to_dyngrid(g, dg, us)
2192  call copy_dyngrid_to_mom_grid(dg, cs%G, us)
2193 
2194  call destroy_dyn_horgrid(dg)
2195  call mom_grid_end(g) ; deallocate(g)
2196 
2197  g => cs%G
2198  if (cs%debug .or. cs%G%symmetric) then
2199  call clone_mom_domain(cs%G%Domain, cs%G%Domain_aux, symmetric=.false.)
2200  else ; cs%G%Domain_aux => cs%G%Domain ;endif
2201  g%ke = gv%ke ; g%g_Earth = gv%mks_g_Earth
2202  endif
2203 
2204  ! At this point, all user-modified initialization code has been called. The
2205  ! remainder of this subroutine is controlled by the parameters that have
2206  ! have already been set.
2207 
2208  if (ale_remap_init_conds(cs%ALE_CSp) .and. .not. query_initialized(cs%h,"h",restart_csp)) then
2209  ! This block is controlled by the ALE parameter REMAP_AFTER_INITIALIZATION.
2210  ! \todo This block exists for legacy reasons and we should phase it out of
2211  ! all examples. !###
2212  if (cs%debug) then
2213  call uvchksum("Pre ALE adjust init cond [uv]", &
2214  cs%u, cs%v, g%HI, haloshift=1)
2215  call hchksum(cs%h,"Pre ALE adjust init cond h", g%HI, haloshift=1, scale=gv%H_to_m)
2216  endif
2217  call calltree_waypoint("Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)")
2218  call adjustgridforintegrity(cs%ALE_CSp, g, gv, cs%h )
2219  call calltree_waypoint("Calling ALE_main() to remap initial conditions (initialize_MOM)")
2220  if (use_ice_shelf) then
2221  filename = trim(inputdir)//trim(ice_shelf_file)
2222  if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
2223  "MOM: Unable to open "//trim(filename))
2224 
2225  allocate(area_shelf_h(isd:ied,jsd:jed))
2226  allocate(frac_shelf_h(isd:ied,jsd:jed))
2227  call mom_read_data(filename, trim(area_varname), area_shelf_h, g%Domain)
2228  ! initialize frac_shelf_h with zeros (open water everywhere)
2229  frac_shelf_h(:,:) = 0.0
2230  ! compute fractional ice shelf coverage of h
2231  do j=jsd,jed ; do i=isd,ied
2232  if (g%areaT(i,j) > 0.0) &
2233  frac_shelf_h(i,j) = area_shelf_h(i,j) / (us%L_to_m**2*g%areaT(i,j))
2234  enddo ; enddo
2235  ! pass to the pointer
2236  shelf_area => frac_shelf_h
2237  call ale_main(g, gv, us, cs%h, cs%u, cs%v, cs%tv, cs%tracer_Reg, cs%ALE_CSp, &
2238  cs%OBC, frac_shelf_h=shelf_area)
2239  else
2240  call ale_main( g, gv, us, cs%h, cs%u, cs%v, cs%tv, cs%tracer_Reg, cs%ALE_CSp, cs%OBC)
2241  endif
2242 
2243  call cpu_clock_begin(id_clock_pass_init)
2244  call create_group_pass(tmp_pass_uv_t_s_h, cs%u, cs%v, g%Domain)
2245  if (use_temperature) then
2246  call create_group_pass(tmp_pass_uv_t_s_h, cs%tv%T, g%Domain, halo=1)
2247  call create_group_pass(tmp_pass_uv_t_s_h, cs%tv%S, g%Domain, halo=1)
2248  endif
2249  call create_group_pass(tmp_pass_uv_t_s_h, cs%h, g%Domain, halo=1)
2250  call do_group_pass(tmp_pass_uv_t_s_h, g%Domain)
2251  call cpu_clock_end(id_clock_pass_init)
2252 
2253  if (cs%debug) then
2254  call uvchksum("Post ALE adjust init cond [uv]", cs%u, cs%v, g%HI, haloshift=1)
2255  call hchksum(cs%h, "Post ALE adjust init cond h", g%HI, haloshift=1, scale=gv%H_to_m)
2256  endif
2257  endif
2258  if ( cs%use_ALE_algorithm ) call ale_updateverticalgridtype( cs%ALE_CSp, gv )
2259 
2260  diag => cs%diag
2261  ! Initialize the diag mediator.
2262  call diag_mediator_init(g, gv, us, gv%ke, param_file, diag, doc_file_dir=dirs%output_directory)
2263  if (present(diag_ptr)) diag_ptr => cs%diag
2264 
2265  ! Initialize the diagnostics masks for native arrays.
2266  ! This step has to be done after call to MOM_initialize_state
2267  ! and before MOM_diagnostics_init
2268  call diag_masks_set(g, gv%ke, diag)
2269 
2270  ! Set up pointers within diag mediator control structure,
2271  ! this needs to occur _after_ CS%h etc. have been allocated.
2272  call diag_set_state_ptrs(cs%h, cs%T, cs%S, cs%tv%eqn_of_state, diag)
2273 
2274  ! This call sets up the diagnostic axes. These are needed,
2275  ! e.g. to generate the target grids below.
2276  call set_axes_info(g, gv, us, param_file, diag)
2277 
2278  ! Whenever thickness/T/S changes let the diag manager know, target grids
2279  ! for vertical remapping may need to be regenerated.
2280  ! FIXME: are h, T, S updated at the same time? Review these for T, S updates.
2281  call diag_update_remap_grids(diag)
2282 
2283  ! Setup the diagnostic grid storage types
2284  call diag_grid_storage_init(cs%diag_pre_sync, g, diag)
2285  call diag_grid_storage_init(cs%diag_pre_dyn, g, diag)
2286 
2287  ! Calculate masks for diagnostics arrays in non-native coordinates
2288  ! This step has to be done after set_axes_info() because the axes needed
2289  ! to be configured, and after diag_update_remap_grids() because the grids
2290  ! must be defined.
2291  call set_masks_for_axes(g, diag)
2292 
2293  ! Diagnose static fields AND associate areas/volumes with axes
2294  call write_static_fields(g, gv, us, cs%tv, cs%diag)
2295  call calltree_waypoint("static fields written (initialize_MOM)")
2296 
2297  ! Register the volume cell measure (must be one of first diagnostics)
2298  call register_cell_measure(g, cs%diag, time)
2299 
2300  call cpu_clock_begin(id_clock_mom_init)
2301  if (cs%use_ALE_algorithm) then
2302  call ale_writecoordinatefile( cs%ALE_CSp, gv, dirs%output_directory )
2303  endif
2304  call cpu_clock_end(id_clock_mom_init)
2305  call calltree_waypoint("ALE initialized (initialize_MOM)")
2306 
2307  cs%useMEKE = meke_init(time, g, us, param_file, diag, cs%MEKE_CSp, cs%MEKE, restart_csp)
2308 
2309  call varmix_init(time, g, gv, us, param_file, diag, cs%VarMix)
2310  call set_visc_init(time, g, gv, us, param_file, diag, cs%visc, cs%set_visc_CSp, restart_csp, cs%OBC)
2311  call thickness_diffuse_init(time, g, gv, us, param_file, diag, cs%CDp, cs%thickness_diffuse_CSp)
2312 
2313  if (cs%split) then
2314  allocate(eta(szi_(g),szj_(g))) ; eta(:,:) = 0.0
2315  call initialize_dyn_split_rk2(cs%u, cs%v, cs%h, cs%uh, cs%vh, eta, time, &
2316  g, gv, us, param_file, diag, cs%dyn_split_RK2_CSp, restart_csp, &
2317  cs%dt, cs%ADp, cs%CDp, mom_internal_state, cs%VarMix, cs%MEKE, &
2318  cs%thickness_diffuse_CSp, &
2319  cs%OBC, cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, &
2320  cs%visc, dirs, cs%ntrunc, calc_dtbt=calc_dtbt)
2321  if (cs%dtbt_reset_period > 0.0) then
2322  cs%dtbt_reset_interval = real_to_time(cs%dtbt_reset_period)
2323  ! Set dtbt_reset_time to be the next even multiple of dtbt_reset_interval.
2324  cs%dtbt_reset_time = time_init + cs%dtbt_reset_interval * &
2325  ((time - time_init) / cs%dtbt_reset_interval)
2326  if ((cs%dtbt_reset_time > time) .and. calc_dtbt) then
2327  ! Back up dtbt_reset_time one interval to force dtbt to be calculated,
2328  ! because the restart was not aligned with the interval to recalculate
2329  ! dtbt, and dtbt was not read from a restart file.
2330  cs%dtbt_reset_time = cs%dtbt_reset_time - cs%dtbt_reset_interval
2331  endif
2332  endif
2333  elseif (cs%use_RK2) then
2334  call initialize_dyn_unsplit_rk2(cs%u, cs%v, cs%h, time, g, gv, us, &
2335  param_file, diag, cs%dyn_unsplit_RK2_CSp, restart_csp, &
2336  cs%ADp, cs%CDp, mom_internal_state, cs%MEKE, cs%OBC, &
2337  cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, cs%visc, dirs, &
2338  cs%ntrunc)
2339  else
2340  call initialize_dyn_unsplit(cs%u, cs%v, cs%h, time, g, gv, us, &
2341  param_file, diag, cs%dyn_unsplit_CSp, restart_csp, &
2342  cs%ADp, cs%CDp, mom_internal_state, cs%MEKE, cs%OBC, &
2343  cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, cs%visc, dirs, &
2344  cs%ntrunc)
2345  endif
2346 
2347  call calltree_waypoint("dynamics initialized (initialize_MOM)")
2348 
2349  cs%mixedlayer_restrat = mixedlayer_restrat_init(time, g, gv, us, param_file, diag, &
2350  cs%mixedlayer_restrat_CSp, restart_csp)
2351  if (cs%mixedlayer_restrat) then
2352  if (.not.(bulkmixedlayer .or. cs%use_ALE_algorithm)) &
2353  call mom_error(fatal, "MOM: MIXEDLAYER_RESTRAT true requires a boundary layer scheme.")
2354  ! When DIABATIC_FIRST=False and using CS%visc%ML in mixedlayer_restrat we need to update after a restart
2355  if (.not. cs%diabatic_first .and. associated(cs%visc%MLD)) &
2356  call pass_var(cs%visc%MLD, g%domain, halo=1)
2357  endif
2358 
2359  call mom_diagnostics_init(mom_internal_state, cs%ADp, cs%CDp, time, g, gv, us, &
2360  param_file, diag, cs%diagnostics_CSp, cs%tv)
2361  call diag_copy_diag_to_storage(cs%diag_pre_sync, cs%h, cs%diag)
2362 
2363  if (associated(cs%sponge_CSp)) &
2364  call init_sponge_diags(time, g, gv, us, diag, cs%sponge_CSp)
2365 
2366  if (associated(cs%ALE_sponge_CSp)) &
2367  call init_ale_sponge_diags(time, g, diag, cs%ALE_sponge_CSp)
2368 
2369  if (cs%adiabatic) then
2370  call adiabatic_driver_init(time, g, param_file, diag, cs%diabatic_CSp, &
2371  cs%tracer_flow_CSp)
2372  else
2373  call diabatic_driver_init(time, g, gv, us, param_file, cs%use_ALE_algorithm, diag, &
2374  cs%ADp, cs%CDp, cs%diabatic_CSp, cs%tracer_flow_CSp, &
2375  cs%sponge_CSp, cs%ALE_sponge_CSp)
2376  endif
2377 
2378  call tracer_advect_init(time, g, us, param_file, diag, cs%tracer_adv_CSp)
2379  call tracer_hor_diff_init(time, g, us, param_file, diag, cs%tv%eqn_of_state, cs%diabatic_CSp, &
2380  cs%tracer_diff_CSp)
2381 
2382  call lock_tracer_registry(cs%tracer_Reg)
2383  call calltree_waypoint("tracer registry now locked (initialize_MOM)")
2384 
2385  ! now register some diagnostics since the tracer registry is now locked
2386  call register_surface_diags(time, g, us, cs%sfc_IDs, cs%diag, cs%tv)
2387  call register_diags(time, g, gv, us, cs%IDs, cs%diag)
2388  call register_transport_diags(time, g, gv, us, cs%transport_IDs, cs%diag)
2389  call register_tracer_diagnostics(cs%tracer_Reg, cs%h, time, diag, g, gv, us, &
2390  cs%use_ALE_algorithm)
2391  if (cs%use_ALE_algorithm) then
2392  call ale_register_diags(time, g, gv, us, diag, cs%ALE_CSp)
2393  endif
2394 
2395  ! This subroutine initializes any tracer packages.
2396  new_sim = is_new_run(restart_csp)
2397  call tracer_flow_control_init(.not.new_sim, time, g, gv, us, cs%h, param_file, &
2398  cs%diag, cs%OBC, cs%tracer_flow_CSp, cs%sponge_CSp, &
2399  cs%ALE_sponge_CSp, cs%tv)
2400  if (present(tracer_flow_csp)) tracer_flow_csp => cs%tracer_flow_CSp
2401 
2402  ! If running in offline tracer mode, initialize the necessary control structure and
2403  ! parameters
2404  if (present(offline_tracer_mode)) offline_tracer_mode=cs%offline_tracer_mode
2405 
2406  if (cs%offline_tracer_mode) then
2407  ! Setup some initial parameterizations and also assign some of the subtypes
2408  call offline_transport_init(param_file, cs%offline_CSp, cs%diabatic_CSp, g, gv, us)
2409  call insert_offline_main( cs=cs%offline_CSp, ale_csp=cs%ALE_CSp, diabatic_csp=cs%diabatic_CSp, &
2410  diag=cs%diag, obc=cs%OBC, tracer_adv_csp=cs%tracer_adv_CSp, &
2411  tracer_flow_csp=cs%tracer_flow_CSp, tracer_reg=cs%tracer_Reg, &
2412  tv=cs%tv, x_before_y = (mod(first_direction,2)==0), debug=cs%debug )
2413  call register_diags_offline_transport(time, cs%diag, cs%offline_CSp)
2414  endif
2415 
2416  !--- set up group pass for u,v,T,S and h. pass_uv_T_S_h also is used in step_MOM
2417  call cpu_clock_begin(id_clock_pass_init)
2418  dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
2419  call create_group_pass(pass_uv_t_s_h, cs%u, cs%v, g%Domain, halo=dynamics_stencil)
2420  if (use_temperature) then
2421  call create_group_pass(pass_uv_t_s_h, cs%tv%T, g%Domain, halo=dynamics_stencil)
2422  call create_group_pass(pass_uv_t_s_h, cs%tv%S, g%Domain, halo=dynamics_stencil)
2423  endif
2424  call create_group_pass(pass_uv_t_s_h, cs%h, g%Domain, halo=dynamics_stencil)
2425 
2426  call do_group_pass(pass_uv_t_s_h, g%Domain)
2427 
2428  if (associated(cs%visc%Kv_shear)) &
2429  call pass_var(cs%visc%Kv_shear, g%Domain, to_all+omit_corners, halo=1)
2430 
2431  if (associated(cs%visc%Kv_slow)) &
2432  call pass_var(cs%visc%Kv_slow, g%Domain, to_all+omit_corners, halo=1)
2433 
2434  call cpu_clock_end(id_clock_pass_init)
2435 
2436  call register_obsolete_diagnostics(param_file, cs%diag)
2437 
2438  if (use_frazil) then
2439  if (.not.query_initialized(cs%tv%frazil,"frazil",restart_csp)) &
2440  cs%tv%frazil(:,:) = 0.0
2441  endif
2442 
2443  if (cs%interp_p_surf) then
2444  cs%p_surf_prev_set = &
2445  query_initialized(cs%p_surf_prev,"p_surf_prev",restart_csp)
2446 
2447  if (cs%p_surf_prev_set) call pass_var(cs%p_surf_prev, g%domain)
2448  endif
2449 
2450  if (.not.query_initialized(cs%ave_ssh_ibc,"ave_ssh",restart_csp)) then
2451  if (cs%split) then
2452  call find_eta(cs%h, cs%tv, g, gv, us, cs%ave_ssh_ibc, eta, eta_to_m=1.0)
2453  else
2454  call find_eta(cs%h, cs%tv, g, gv, us, cs%ave_ssh_ibc, eta_to_m=1.0)
2455  endif
2456  endif
2457  if (cs%split) deallocate(eta)
2458 
2459  cs%nstep_tot = 0
2460  if (present(count_calls)) cs%count_calls = count_calls
2461  call mom_sum_output_init(g, us, param_file, dirs%output_directory, &
2462  cs%ntrunc, time_init, cs%sum_output_CSp)
2463 
2464  ! Flag whether to save initial conditions in finish_MOM_initialization() or not.
2465  cs%write_IC = save_ic .and. &
2466  .not.((dirs%input_filename(1:1) == 'r') .and. &
2467  (len_trim(dirs%input_filename) == 1))
2468 
2469  if (cs%ensemble_ocean) then
2470  call init_oda(time, g, gv, cs%odaCS)
2471  endif
2472 
2473  !### This could perhaps go here instead of in finish_MOM_initialization?
2474  ! call fix_restart_scaling(GV)
2475  ! call fix_restart_unit_scaling(US)
2476 
2477  call calltree_leave("initialize_MOM()")
2478  call cpu_clock_end(id_clock_init)
2479 

References mom_diabatic_driver::adiabatic_driver_init(), mom_boundary_update::call_obc_register(), mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), mom_error_handler::calltree_waypoint(), mom_transcribe_grid::copy_dyngrid_to_mom_grid(), mom_transcribe_grid::copy_mom_grid_to_dyngrid(), mom_dyn_horgrid::create_dyn_horgrid(), mom_dyn_horgrid::destroy_dyn_horgrid(), mom_diag_mediator::diag_copy_diag_to_storage(), mom_eos::eos_init(), mom_obsolete_params::find_obsolete_params(), mom_get_input::get_mom_input(), mom_verticalgrid::get_tr_flux_units(), mom_hor_index::hor_index_init(), id_clock_init, id_clock_mom_init, id_clock_pass_init, mom_ale_sponge::init_ale_sponge_diags(), mom_sponge::init_sponge_diags(), mom_restart::is_new_run(), mom_tracer_registry::lock_tracer_registry(), mom_meke::meke_alloc_register_restart(), mom_meke::meke_init(), mom_mixed_layer_restrat::mixedlayer_restrat_register_restarts(), mom_coord_initialization::mom_initialize_coord(), mom_fixed_initialization::mom_initialize_fixed(), mom_state_initialization::mom_initialize_state(), mom_sum_output::mom_sum_output_init(), mom_timing_init(), register_diags(), mom_obsolete_diagnostics::register_obsolete_diagnostics(), mom_restart::restart_init(), mom_grid::set_first_direction(), set_restart_fields(), mom_set_visc::set_visc_register_restarts(), mom_thickness_diffuse::thickness_diffuse_init(), and mom_unit_tests::unit_tests().

Referenced by mom_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ mom_end()

subroutine, public mom::mom_end ( type(mom_control_struct), pointer  CS)

End of ocean model, including memory deallocation.

Parameters
csMOM control structure

Definition at line 3132 of file MOM.F90.

3132  type(MOM_control_struct), pointer :: CS !< MOM control structure
3133 
3134  if (cs%use_ALE_algorithm) call ale_end(cs%ALE_CSp)
3135 
3136  dealloc_(cs%u) ; dealloc_(cs%v) ; dealloc_(cs%h)
3137  dealloc_(cs%uh) ; dealloc_(cs%vh)
3138 
3139  if (associated(cs%tv%T)) then
3140  dealloc_(cs%T) ; cs%tv%T => null() ; dealloc_(cs%S) ; cs%tv%S => null()
3141  endif
3142  if (associated(cs%tv%frazil)) deallocate(cs%tv%frazil)
3143  if (associated(cs%tv%salt_deficit)) deallocate(cs%tv%salt_deficit)
3144  if (associated(cs%Hml)) deallocate(cs%Hml)
3145 
3146  call tracer_advect_end(cs%tracer_adv_CSp)
3147  call tracer_hor_diff_end(cs%tracer_diff_CSp)
3148  call tracer_registry_end(cs%tracer_Reg)
3149  call tracer_flow_control_end(cs%tracer_flow_CSp)
3150 
3151  call diabatic_driver_end(cs%diabatic_CSp)
3152 
3153  if (cs%offline_tracer_mode) call offline_transport_end(cs%offline_CSp)
3154 
3155  dealloc_(cs%uhtr) ; dealloc_(cs%vhtr)
3156  if (cs%split) then
3157  call end_dyn_split_rk2(cs%dyn_split_RK2_CSp)
3158  elseif (cs%use_RK2) then
3159  call end_dyn_unsplit_rk2(cs%dyn_unsplit_RK2_CSp)
3160  else
3161  call end_dyn_unsplit(cs%dyn_unsplit_CSp)
3162  endif
3163  dealloc_(cs%ave_ssh_ibc) ; dealloc_(cs%ssh_rint)
3164  if (associated(cs%update_OBC_CSp)) call obc_register_end(cs%update_OBC_CSp)
3165 
3166  call verticalgridend(cs%GV)
3167  call unit_scaling_end(cs%US)
3168  call mom_grid_end(cs%G)
3169 
3170  deallocate(cs)
3171 

References mom_diabatic_driver::diabatic_driver_end(), mom_boundary_update::obc_register_end(), mom_offline_main::offline_transport_end(), mom_tracer_advect::tracer_advect_end(), mom_tracer_flow_control::tracer_flow_control_end(), mom_tracer_hor_diff::tracer_hor_diff_end(), mom_tracer_registry::tracer_registry_end(), and mom_unit_scaling::unit_scaling_end().

Referenced by mom_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ mom_state_is_synchronized()

logical function, public mom::mom_state_is_synchronized ( type(mom_control_struct), pointer  CS,
logical, intent(in), optional  adv_dyn 
)

Return true if all phases of step_MOM are at the same point in time.

Parameters
csMOM control structure
[in]adv_dynIf present and true, only check whether the advection is up-to-date with the dynamics.
Returns
True if all phases of the update are synchronized.

Definition at line 3075 of file MOM.F90.

3075  type(MOM_control_struct), pointer :: CS !< MOM control structure
3076  logical, optional, intent(in) :: adv_dyn !< If present and true, only check
3077  !! whether the advection is up-to-date with
3078  !! the dynamics.
3079  logical :: in_synch !< True if all phases of the update are synchronized.
3080 
3081  logical :: adv_only
3082 
3083  adv_only = .false. ; if (present(adv_dyn)) adv_only = adv_dyn
3084 
3085  if (adv_only) then
3086  in_synch = (cs%t_dyn_rel_adv == 0.0)
3087  else
3088  in_synch = ((cs%t_dyn_rel_adv == 0.0) .and. (cs%t_dyn_rel_thermo == 0.0))
3089  endif
3090 

Referenced by mom_main(), and step_mom().

Here is the caller graph for this function:

◆ mom_timing_init()

subroutine mom::mom_timing_init ( type(mom_control_struct), intent(in)  CS)
private

Set up CPU clock IDs for timing various subroutines.

Parameters
[in]cscontrol structure set up by initialize_MOM.

Definition at line 2565 of file MOM.F90.

2565  type(MOM_control_struct), intent(in) :: CS !< control structure set up by initialize_MOM.
2566 
2567  id_clock_ocean = cpu_clock_id('Ocean', grain=clock_component)
2568  id_clock_dynamics = cpu_clock_id('Ocean dynamics', grain=clock_subcomponent)
2569  id_clock_thermo = cpu_clock_id('Ocean thermodynamics and tracers', grain=clock_subcomponent)
2570  id_clock_other = cpu_clock_id('Ocean Other', grain=clock_subcomponent)
2571  id_clock_tracer = cpu_clock_id('(Ocean tracer advection)', grain=clock_module_driver)
2572  if (.not.cs%adiabatic) then
2573  id_clock_diabatic = cpu_clock_id('(Ocean diabatic driver)', grain=clock_module_driver)
2574  else
2575  id_clock_adiabatic = cpu_clock_id('(Ocean adiabatic driver)', grain=clock_module_driver)
2576  endif
2577 
2578  id_clock_continuity = cpu_clock_id('(Ocean continuity equation *)', grain=clock_module)
2579  id_clock_bbl_visc = cpu_clock_id('(Ocean set BBL viscosity)', grain=clock_module)
2580  id_clock_pass = cpu_clock_id('(Ocean message passing *)', grain=clock_module)
2581  id_clock_mom_init = cpu_clock_id('(Ocean MOM_initialize_state)', grain=clock_module)
2582  id_clock_pass_init = cpu_clock_id('(Ocean init message passing *)', grain=clock_routine)
2583  if (cs%thickness_diffuse) &
2584  id_clock_thick_diff = cpu_clock_id('(Ocean thickness diffusion *)', grain=clock_module)
2585 !if (CS%mixedlayer_restrat) &
2586  id_clock_ml_restrat = cpu_clock_id('(Ocean mixed layer restrat)', grain=clock_module)
2587  id_clock_diagnostics = cpu_clock_id('(Ocean collective diagnostics)', grain=clock_module)
2588  id_clock_z_diag = cpu_clock_id('(Ocean Z-space diagnostics)', grain=clock_module)
2589  id_clock_ale = cpu_clock_id('(Ocean ALE)', grain=clock_module)
2590  if (cs%offline_tracer_mode) then
2591  id_clock_offline_tracer = cpu_clock_id('Ocean offline tracers', grain=clock_subcomponent)
2592  endif
2593 

References id_clock_adiabatic, id_clock_ale, id_clock_bbl_visc, id_clock_continuity, id_clock_diabatic, id_clock_diagnostics, id_clock_dynamics, id_clock_ml_restrat, id_clock_mom_init, id_clock_ocean, id_clock_offline_tracer, id_clock_other, id_clock_pass, id_clock_pass_init, id_clock_thermo, id_clock_thick_diff, id_clock_tracer, and id_clock_z_diag.

Referenced by initialize_mom().

Here is the caller graph for this function:

◆ register_diags()

subroutine mom::register_diags ( type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(inout)  US,
type(mom_diag_ids), intent(inout)  IDs,
type(diag_ctrl), intent(inout)  diag 
)
private

Register certain diagnostics.

Parameters
[in]timecurrent model time
[in]gocean grid structure
[in]gvocean vertical grid structure
[in,out]usA dimensional unit scaling type
[in,out]idsA structure with the diagnostic IDs.
[in,out]diagregulates diagnostic output

Definition at line 2534 of file MOM.F90.

2534  type(time_type), intent(in) :: Time !< current model time
2535  type(ocean_grid_type), intent(in) :: G !< ocean grid structure
2536  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2537  type(unit_scale_type), intent(inout) :: US !< A dimensional unit scaling type
2538  type(MOM_diag_IDs), intent(inout) :: IDs !< A structure with the diagnostic IDs.
2539  type(diag_ctrl), intent(inout) :: diag !< regulates diagnostic output
2540 
2541  real :: H_convert
2542  character(len=48) :: thickness_units
2543 
2544  thickness_units = get_thickness_units(gv)
2545  if (gv%Boussinesq) then
2546  h_convert = gv%H_to_m
2547  else
2548  h_convert = gv%H_to_kg_m2
2549  endif
2550 
2551  ! Diagnostics of the rapidly varying dynamic state
2552  ids%id_u = register_diag_field('ocean_model', 'u_dyn', diag%axesCuL, time, &
2553  'Zonal velocity after the dynamics update', 'm s-1', conversion=us%L_T_to_m_s)
2554  ids%id_v = register_diag_field('ocean_model', 'v_dyn', diag%axesCvL, time, &
2555  'Meridional velocity after the dynamics update', 'm s-1', conversion=us%L_T_to_m_s)
2556  ids%id_h = register_diag_field('ocean_model', 'h_dyn', diag%axesTL, time, &
2557  'Layer Thickness after the dynamics update', thickness_units, &
2558  v_extensive=.true., conversion=h_convert)
2559  ids%id_ssh_inst = register_diag_field('ocean_model', 'SSH_inst', diag%axesT1, &
2560  time, 'Instantaneous Sea Surface Height', 'm')

References mom_verticalgrid::get_thickness_units().

Referenced by initialize_mom().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_restart_fields()

subroutine mom::set_restart_fields ( type(verticalgrid_type), intent(inout)  GV,
type(unit_scale_type), intent(inout)  US,
type(param_file_type), intent(in)  param_file,
type(mom_control_struct), intent(in)  CS,
type(mom_restart_cs), pointer  restart_CSp 
)
private

Set the fields that are needed for bitwise identical restarting the time stepping scheme. In addition to those specified here directly, there may be fields related to the forcing or to the barotropic solver that are needed; these are specified in sub- routines that are called from this one.

This routine should be altered if there are any changes to the time stepping scheme. The CHECK_RESTART facility may be used to confirm that all needed restart fields have been included.

Parameters
[in,out]gvocean vertical grid structure
[in,out]usA dimensional unit scaling type
[in]param_fileopened file for parsing to get parameters
[in]cscontrol structure set up by inialize_MOM
restart_csppointer to the restart control structure that will be used for MOM.

Definition at line 2606 of file MOM.F90.

2606  type(verticalGrid_type), intent(inout) :: GV !< ocean vertical grid structure
2607  type(unit_scale_type), intent(inout) :: US !< A dimensional unit scaling type
2608  type(param_file_type), intent(in) :: param_file !< opened file for parsing to get parameters
2609  type(MOM_control_struct), intent(in) :: CS !< control structure set up by inialize_MOM
2610  type(MOM_restart_CS), pointer :: restart_CSp !< pointer to the restart control
2611  !! structure that will be used for MOM.
2612  ! Local variables
2613  logical :: use_ice_shelf ! Needed to determine whether to add CS%Hml to restarts
2614  character(len=48) :: thickness_units, flux_units
2615 
2616 
2617  thickness_units = get_thickness_units(gv)
2618  flux_units = get_flux_units(gv)
2619 
2620  if (associated(cs%tv%T)) &
2621  call register_restart_field(cs%tv%T, "Temp", .true., restart_csp, &
2622  "Potential Temperature", "degC")
2623  if (associated(cs%tv%S)) &
2624  call register_restart_field(cs%tv%S, "Salt", .true., restart_csp, &
2625  "Salinity", "PPT")
2626 
2627  call register_restart_field(cs%h, "h", .true., restart_csp, &
2628  "Layer Thickness", thickness_units)
2629 
2630  call register_restart_field(cs%u, "u", .true., restart_csp, &
2631  "Zonal velocity", "m s-1", hor_grid='Cu')
2632 
2633  call register_restart_field(cs%v, "v", .true., restart_csp, &
2634  "Meridional velocity", "m s-1", hor_grid='Cv')
2635 
2636  if (associated(cs%tv%frazil)) &
2637  call register_restart_field(cs%tv%frazil, "frazil", .false., restart_csp, &
2638  "Frazil heat flux into ocean", "J m-2")
2639 
2640  if (cs%interp_p_surf) then
2641  call register_restart_field(cs%p_surf_prev, "p_surf_prev", .false., restart_csp, &
2642  "Previous ocean surface pressure", "Pa")
2643  endif
2644 
2645  call register_restart_field(cs%ave_ssh_ibc, "ave_ssh", .false., restart_csp, &
2646  "Time average sea surface height", "meter")
2647 
2648  ! hML is needed when using the ice shelf module
2649  call get_param(param_file, '', "ICE_SHELF", use_ice_shelf, default=.false., &
2650  do_not_log=.true.)
2651  if (use_ice_shelf .and. associated(cs%Hml)) then
2652  call register_restart_field(cs%Hml, "hML", .false., restart_csp, &
2653  "Mixed layer thickness", "meter")
2654  endif
2655 
2656  ! Register scalar unit conversion factors.
2657  call register_restart_field(us%m_to_Z_restart, "m_to_Z", .false., restart_csp, &
2658  "Height unit conversion factor", "Z meter-1")
2659  call register_restart_field(gv%m_to_H_restart, "m_to_H", .false., restart_csp, &
2660  "Thickness unit conversion factor", "H meter-1")
2661  call register_restart_field(us%m_to_L_restart, "m_to_L", .false., restart_csp, &
2662  "Length unit conversion factor", "L meter-1")
2663  call register_restart_field(us%s_to_T_restart, "s_to_T", .false., restart_csp, &
2664  "Time unit conversion factor", "T second-1")
2665  call register_restart_field(us%kg_m3_to_R_restart, "kg_m3_to_R", .false., restart_csp, &
2666  "Density unit conversion factor", "R m3 kg-1")
2667 

References mom_verticalgrid::get_flux_units(), and mom_verticalgrid::get_thickness_units().

Referenced by initialize_mom().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ step_mom()

subroutine, public mom::step_mom ( type(mech_forcing), intent(inout)  forces,
type(forcing), intent(inout)  fluxes,
type(surface), intent(inout)  sfc_state,
type(time_type), intent(in)  Time_start,
real, intent(in)  time_int_in,
type(mom_control_struct), pointer  CS,
type(wave_parameters_cs), optional, pointer  Waves,
logical, intent(in), optional  do_dynamics,
logical, intent(in), optional  do_thermodynamics,
logical, intent(in), optional  start_cycle,
logical, intent(in), optional  end_cycle,
real, intent(in), optional  cycle_length,
logical, intent(in), optional  reset_therm 
)

This subroutine orchestrates the time stepping of MOM. The adiabatic dynamics are stepped by calls to one of the step_MOM_dyn_...routines. The action of lateral processes on tracers occur in calls to advect_tracer and tracer_hordiff. Vertical mixing and possibly remapping occur inside of diabatic.

Parameters
[in,out]forcesA structure with the driving mechanical forces
[in,out]fluxesA structure with pointers to themodynamic, tracer and mass exchange forcing fields
[in,out]sfc_statesurface ocean state
[in]time_startstarting time of a segment, as a time type
[in]time_int_intime interval covered by this run segment [s].
cscontrol structure from initialize_MOM
wavesAn optional pointer to a wave property CS
[in]do_dynamicsPresent and false, do not do updates due to the dynamics.
[in]do_thermodynamicsPresent and false, do not do updates due to the thermodynamics or remapping.
[in]start_cycleThis indicates whether this call is to be treated as the first call to step_MOM in a time-stepping cycle; missing is like true.
[in]end_cycleThis indicates whether this call is to be treated as the last call to step_MOM in a time-stepping cycle; missing is like true.
[in]cycle_lengthThe amount of time in a coupled time stepping cycle [s].
[in]reset_thermThis indicates whether the running sums of thermodynamic quantities should be reset. If missing, this is like start_cycle.

Definition at line 399 of file MOM.F90.

399  type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
400  type(forcing), intent(inout) :: fluxes !< A structure with pointers to themodynamic,
401  !! tracer and mass exchange forcing fields
402  type(surface), intent(inout) :: sfc_state !< surface ocean state
403  type(time_type), intent(in) :: Time_start !< starting time of a segment, as a time type
404  real, intent(in) :: time_int_in !< time interval covered by this run segment [s].
405  type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM
406  type(Wave_parameters_CS), &
407  optional, pointer :: Waves !< An optional pointer to a wave property CS
408  logical, optional, intent(in) :: do_dynamics !< Present and false, do not do updates due
409  !! to the dynamics.
410  logical, optional, intent(in) :: do_thermodynamics !< Present and false, do not do updates due
411  !! to the thermodynamics or remapping.
412  logical, optional, intent(in) :: start_cycle !< This indicates whether this call is to be
413  !! treated as the first call to step_MOM in a
414  !! time-stepping cycle; missing is like true.
415  logical, optional, intent(in) :: end_cycle !< This indicates whether this call is to be
416  !! treated as the last call to step_MOM in a
417  !! time-stepping cycle; missing is like true.
418  real, optional, intent(in) :: cycle_length !< The amount of time in a coupled time
419  !! stepping cycle [s].
420  logical, optional, intent(in) :: reset_therm !< This indicates whether the running sums of
421  !! thermodynamic quantities should be reset.
422  !! If missing, this is like start_cycle.
423 
424  ! local variables
425  type(ocean_grid_type), pointer :: G => null() ! pointer to a structure containing
426  ! metrics and related information
427  type(verticalGrid_type), pointer :: GV => null() ! Pointer to the vertical grid structure
428  type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing
429  ! various unit conversion factors
430  integer :: ntstep ! time steps between tracer updates or diabatic forcing
431  integer :: n_max ! number of steps to take in this call
432 
433  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, n
434  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
435 
436  real :: time_interval ! time interval covered by this run segment [T ~> s].
437  real :: dt ! baroclinic time step [T ~> s]
438  real :: dtdia ! time step for diabatic processes [T ~> s]
439  real :: dt_therm ! a limited and quantized version of CS%dt_therm [T ~> s]
440  real :: dt_therm_here ! a further limited value of dt_therm [T ~> s]
441 
442  real :: wt_end, wt_beg
443  real :: bbl_time_int ! The amount of time over which the calculated BBL
444  ! properties will apply, for use in diagnostics, or 0
445  ! if it is not to be calculated anew [T ~> s].
446  real :: rel_time = 0.0 ! relative time since start of this call [T ~> s].
447 
448  logical :: calc_dtbt ! Indicates whether the dynamically adjusted
449  ! barotropic time step needs to be updated.
450  logical :: do_advection ! If true, it is time to advect tracers.
451  logical :: do_calc_bbl ! If true, calculate the boundary layer properties.
452  logical :: thermo_does_span_coupling ! If true, thermodynamic forcing spans
453  ! multiple dynamic timesteps.
454  logical :: do_dyn ! If true, dynamics are updated with this call.
455  logical :: do_thermo ! If true, thermodynamics and remapping may be applied with this call.
456  logical :: cycle_start ! If true, do calculations that are only done at the start of
457  ! a stepping cycle (whatever that may mean).
458  logical :: cycle_end ! If true, do calculations and diagnostics that are only done at
459  ! the end of a stepping cycle (whatever that may mean).
460  logical :: therm_reset ! If true, reset running sums of thermodynamic quantities.
461  real :: cycle_time ! The length of the coupled time-stepping cycle [T ~> s].
462  real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
463  ssh ! sea surface height, which may be based on eta_av [m]
464 
465  real, dimension(:,:,:), pointer :: &
466  u => null(), & ! u : zonal velocity component [L T-1 ~> m s-1]
467  v => null(), & ! v : meridional velocity component [L T-1 ~> m s-1]
468  h => null() ! h : layer thickness [H ~> m or kg m-2]
469  real, dimension(:,:), pointer :: &
470  p_surf => null() ! A pointer to the ocean surface pressure [Pa].
471  real :: I_wt_ssh ! The inverse of the time weights [T-1 ~> s-1]
472 
473  type(time_type) :: Time_local, end_time_thermo, Time_temp
474  type(group_pass_type) :: pass_tau_ustar_psurf
475  logical :: showCallTree
476 
477  g => cs%G ; gv => cs%GV ; us => cs%US
478  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
479  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
480  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
481  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
482  u => cs%u ; v => cs%v ; h => cs%h
483 
484  time_interval = us%s_to_T*time_int_in
485  do_dyn = .true. ; if (present(do_dynamics)) do_dyn = do_dynamics
486  do_thermo = .true. ; if (present(do_thermodynamics)) do_thermo = do_thermodynamics
487  if (.not.(do_dyn .or. do_thermo)) call mom_error(fatal,"Step_MOM: "//&
488  "Both do_dynamics and do_thermodynamics are false, which makes no sense.")
489  cycle_start = .true. ; if (present(start_cycle)) cycle_start = start_cycle
490  cycle_end = .true. ; if (present(end_cycle)) cycle_end = end_cycle
491  cycle_time = time_interval ; if (present(cycle_length)) cycle_time = us%s_to_T*cycle_length
492  therm_reset = cycle_start ; if (present(reset_therm)) therm_reset = reset_therm
493 
494  call cpu_clock_begin(id_clock_ocean)
495  call cpu_clock_begin(id_clock_other)
496 
497  if (cs%debug) then
498  call mom_state_chksum("Beginning of step_MOM ", u, v, h, cs%uh, cs%vh, g, gv, us)
499  endif
500 
501  showcalltree = calltree_showquery()
502  if (showcalltree) call calltree_enter("step_MOM(), MOM.F90")
503 
504  ! First determine the time step that is consistent with this call and an
505  ! integer fraction of time_interval.
506  if (do_dyn) then
507  n_max = 1
508  if (time_interval > cs%dt) n_max = ceiling(time_interval/cs%dt - 0.001)
509  dt = time_interval / real(n_max)
510  thermo_does_span_coupling = (cs%thermo_spans_coupling .and. &
511  (cs%dt_therm > 1.5*cycle_time))
512  if (thermo_does_span_coupling) then
513  ! Set dt_therm to be an integer multiple of the coupling time step.
514  dt_therm = cycle_time * floor(cs%dt_therm / cycle_time + 0.001)
515  ntstep = floor(dt_therm/dt + 0.001)
516  elseif (.not.do_thermo) then
517  dt_therm = cs%dt_therm
518  if (present(cycle_length)) dt_therm = min(cs%dt_therm, us%s_to_T*cycle_length)
519  ! ntstep is not used.
520  else
521  ntstep = max(1, min(n_max, floor(cs%dt_therm/dt + 0.001)))
522  dt_therm = dt*ntstep
523  endif
524 
525  if (associated(forces%p_surf)) p_surf => forces%p_surf
526  if (.not.associated(forces%p_surf)) cs%interp_p_surf = .false.
527 
528  !---------- Initiate group halo pass of the forcing fields
529  call cpu_clock_begin(id_clock_pass)
530  call create_group_pass(pass_tau_ustar_psurf, forces%taux, forces%tauy, g%Domain)
531  if (associated(forces%ustar)) &
532  call create_group_pass(pass_tau_ustar_psurf, forces%ustar, g%Domain)
533  if (associated(forces%p_surf)) &
534  call create_group_pass(pass_tau_ustar_psurf, forces%p_surf, g%Domain)
535  if (g%nonblocking_updates) then
536  call start_group_pass(pass_tau_ustar_psurf, g%Domain)
537  else
538  call do_group_pass(pass_tau_ustar_psurf, g%Domain)
539  endif
540  call cpu_clock_end(id_clock_pass)
541  else
542  ! This step only updates the thermodynamics so setting timesteps is simpler.
543  n_max = 1
544  if ((time_interval > cs%dt_therm) .and. (cs%dt_therm > 0.0)) &
545  n_max = ceiling(time_interval/cs%dt_therm - 0.001)
546 
547  dt = time_interval / real(n_max)
548  dt_therm = dt ; ntstep = 1
549  if (associated(fluxes%p_surf)) p_surf => fluxes%p_surf
550 
551  if (cs%UseWaves) call pass_var(fluxes%ustar, g%Domain, clock=id_clock_pass)
552  endif
553 
554  if (therm_reset) then
555  cs%time_in_thermo_cycle = 0.0
556  if (associated(cs%tv%frazil)) cs%tv%frazil(:,:) = 0.0
557  if (associated(cs%tv%salt_deficit)) cs%tv%salt_deficit(:,:) = 0.0
558  if (associated(cs%tv%TempxPmE)) cs%tv%TempxPmE(:,:) = 0.0
559  if (associated(cs%tv%internal_heat)) cs%tv%internal_heat(:,:) = 0.0
560  endif
561 
562  if (cycle_start) then
563  cs%time_in_cycle = 0.0
564  do j=js,je ; do i=is,ie ; cs%ssh_rint(i,j) = 0.0 ; enddo ; enddo
565 
566  if (associated(cs%VarMix)) then
567  call enable_averages(cycle_time, time_start + real_to_time(us%T_to_s*cycle_time), cs%diag)
568  call calc_resoln_function(h, cs%tv, g, gv, us, cs%VarMix)
569  call calc_depth_function(g, cs%VarMix)
570  call disable_averaging(cs%diag)
571  endif
572  endif
573 
574  if (do_dyn) then
575  if (g%nonblocking_updates) &
576  call complete_group_pass(pass_tau_ustar_psurf, g%Domain, clock=id_clock_pass)
577 
578  if (cs%interp_p_surf) then
579  if (.not.associated(cs%p_surf_end)) allocate(cs%p_surf_end(isd:ied,jsd:jed))
580  if (.not.associated(cs%p_surf_begin)) allocate(cs%p_surf_begin(isd:ied,jsd:jed))
581  if (.not.cs%p_surf_prev_set) then
582  do j=jsd,jed ; do i=isd,ied
583  cs%p_surf_prev(i,j) = forces%p_surf(i,j)
584  enddo ; enddo
585  cs%p_surf_prev_set = .true.
586  endif
587  else
588  cs%p_surf_end => forces%p_surf
589  endif
590 
591  if (cs%UseWaves) then
592  ! Update wave information, which is presently kept static over each call to step_mom
593  call enable_averages(time_interval, time_start + real_to_time(us%T_to_s*time_interval), cs%diag)
594  call update_stokes_drift(g, gv, us, waves, h, forces%ustar)
595  call disable_averaging(cs%diag)
596  endif
597  else ! not do_dyn.
598  if (cs%UseWaves) & ! Diagnostics are not enabled in this call.
599  call update_stokes_drift(g, gv, us, waves, h, fluxes%ustar)
600  endif
601 
602  if (cs%debug) then
603  if (cycle_start) &
604  call mom_state_chksum("Before steps ", u, v, h, cs%uh, cs%vh, g, gv, us)
605  if (cycle_start) call check_redundant("Before steps ", u, v, g)
606  if (do_dyn) call mom_mech_forcing_chksum("Before steps", forces, g, us, haloshift=0)
607  if (do_dyn) call check_redundant("Before steps ", forces%taux, forces%tauy, g)
608  endif
609  call cpu_clock_end(id_clock_other)
610 
611  rel_time = 0.0
612  do n=1,n_max
613  rel_time = rel_time + dt ! The relative time at the end of the step.
614  ! Set the universally visible time to the middle of the time step.
615  cs%Time = time_start + real_to_time(us%T_to_s*(rel_time - 0.5*dt))
616  ! Set the local time to the end of the time step.
617  time_local = time_start + real_to_time(us%T_to_s*rel_time)
618 
619  if (showcalltree) call calltree_enter("DT cycles (step_MOM) n=",n)
620 
621  !===========================================================================
622  ! This is the first place where the diabatic processes and remapping could occur.
623  if (cs%diabatic_first .and. (cs%t_dyn_rel_adv==0.0) .and. do_thermo) then ! do thermodynamics.
624 
625  if (.not.do_dyn) then
626  dtdia = dt
627  elseif (thermo_does_span_coupling) then
628  dtdia = dt_therm
629  if ((fluxes%dt_buoy_accum > 0.0) .and. (dtdia > time_interval) .and. &
630  (abs(fluxes%dt_buoy_accum - dtdia) > 1e-6*dtdia)) then
631  call mom_error(fatal, "step_MOM: Mismatch between long thermodynamic "//&
632  "timestep and time over which buoyancy fluxes have been accumulated.")
633  endif
634  call mom_error(fatal, "MOM is not yet set up to have restarts that work "//&
635  "with THERMO_SPANS_COUPLING and DIABATIC_FIRST.")
636  else
637  dtdia = dt*min(ntstep,n_max-(n-1))
638  endif
639 
640  end_time_thermo = time_local
641  if (dtdia > dt) then
642  ! If necessary, temporarily reset CS%Time to the center of the period covered
643  ! by the call to step_MOM_thermo, noting that they begin at the same time.
644  cs%Time = cs%Time + real_to_time(0.5*us%T_to_s*(dtdia-dt))
645  ! The end-time of the diagnostic interval needs to be set ahead if there
646  ! are multiple dynamic time steps worth of thermodynamics applied here.
647  end_time_thermo = time_local + real_to_time(us%T_to_s*(dtdia-dt))
648  endif
649 
650  ! Apply diabatic forcing, do mixing, and regrid.
651  call step_mom_thermo(cs, g, gv, us, u, v, h, cs%tv, fluxes, dtdia, &
652  end_time_thermo, .true., waves=waves)
653  cs%time_in_thermo_cycle = cs%time_in_thermo_cycle + dtdia
654 
655  ! The diabatic processes are now ahead of the dynamics by dtdia.
656  cs%t_dyn_rel_thermo = -dtdia
657  if (showcalltree) call calltree_waypoint("finished diabatic_first (step_MOM)")
658 
659  if (dtdia > dt) & ! Reset CS%Time to its previous value.
660  cs%Time = time_start + real_to_time(us%T_to_s*(rel_time - 0.5*dt))
661  endif ! end of block "(CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0))"
662 
663  if (do_dyn) then
664  ! Store pre-dynamics grids for proper diagnostic remapping for transports
665  ! or advective tendencies. If there are more dynamics steps per advective
666  ! steps (i.e DT_THERM /= DT), this needs to be stored at the first call.
667  if (cs%ndyn_per_adv == 0 .and. cs%t_dyn_rel_adv == 0.) then
668  call diag_copy_diag_to_storage(cs%diag_pre_dyn, h, cs%diag)
669  cs%ndyn_per_adv = cs%ndyn_per_adv + 1
670  endif
671 
672  ! The pre-dynamics velocities might be stored for debugging truncations.
673  if (associated(cs%u_prev) .and. associated(cs%v_prev)) then
674  do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb
675  cs%u_prev(i,j,k) = u(i,j,k)
676  enddo ; enddo ; enddo
677  do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
678  cs%v_prev(i,j,k) = v(i,j,k)
679  enddo ; enddo ; enddo
680  endif
681 
682  dt_therm_here = dt_therm
683  if (do_thermo .and. do_dyn .and. .not.thermo_does_span_coupling) &
684  dt_therm_here = dt*min(ntstep, n_max-n+1)
685 
686  ! Indicate whether the bottom boundary layer properties need to be
687  ! recalculated, and if so for how long an interval they are valid.
688  bbl_time_int = 0.0
689  if (do_thermo) then
690  if ((cs%t_dyn_rel_adv == 0.0) .or. (n==1)) &
691  bbl_time_int = max(dt, min(dt_therm - cs%t_dyn_rel_adv, dt*(1+n_max-n)) )
692  else
693  if ((cs%t_dyn_rel_adv == 0.0) .or. ((n==1) .and. cycle_start)) &
694  bbl_time_int = min(dt_therm, cycle_time)
695  endif
696 
697  if (cs%interp_p_surf) then
698  wt_end = real(n) / real(n_max)
699  wt_beg = real(n-1) / real(n_max)
700  do j=jsd,jed ; do i=isd,ied
701  cs%p_surf_end(i,j) = wt_end * forces%p_surf(i,j) + &
702  (1.0-wt_end) * cs%p_surf_prev(i,j)
703  cs%p_surf_begin(i,j) = wt_beg * forces%p_surf(i,j) + &
704  (1.0-wt_beg) * cs%p_surf_prev(i,j)
705  enddo ; enddo
706  endif
707 
708  call step_mom_dynamics(forces, cs%p_surf_begin, cs%p_surf_end, dt, &
709  dt_therm_here, bbl_time_int, cs, &
710  time_local, waves=waves)
711 
712  !===========================================================================
713  ! This is the start of the tracer advection part of the algorithm.
714 
715  if (thermo_does_span_coupling .or. .not.do_thermo) then
716  do_advection = (cs%t_dyn_rel_adv + 0.5*dt > dt_therm)
717  else
718  do_advection = ((mod(n,ntstep) == 0) .or. (n==n_max))
719  endif
720 
721  if (do_advection) then ! Do advective transport and lateral tracer mixing.
722  call step_mom_tracer_dyn(cs, g, gv, us, h, time_local)
723  cs%ndyn_per_adv = 0
724  if (cs%diabatic_first .and. abs(cs%t_dyn_rel_thermo) > 1e-6*dt) call mom_error(fatal, &
725  "step_MOM: Mismatch between the dynamics and diabatic times "//&
726  "with DIABATIC_FIRST.")
727  endif
728  endif ! end of (do_dyn)
729 
730  !===========================================================================
731  ! This is the second place where the diabatic processes and remapping could occur.
732  if ((cs%t_dyn_rel_adv==0.0) .and. do_thermo .and. (.not.cs%diabatic_first)) then
733 
734  dtdia = cs%t_dyn_rel_thermo
735  ! If the MOM6 dynamic and thermodynamic time stepping is being orchestrated
736  ! by the coupler, the value of diabatic_first does not matter.
737  if ((cs%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) dtdia = dt
738 
739  if (cs%thermo_spans_coupling .and. (cs%dt_therm > 1.5*cycle_time) .and. &
740  (abs(dt_therm - dtdia) > 1e-6*dt_therm)) then
741  call mom_error(fatal, "step_MOM: Mismatch between dt_therm and dtdia "//&
742  "before call to diabatic.")
743  endif
744 
745  ! If necessary, temporarily reset CS%Time to the center of the period covered
746  ! by the call to step_MOM_thermo, noting that they end at the same time.
747  if (dtdia > dt) cs%Time = cs%Time - real_to_time(0.5*us%T_to_s*(dtdia-dt))
748 
749  ! Apply diabatic forcing, do mixing, and regrid.
750  call step_mom_thermo(cs, g, gv, us, u, v, h, cs%tv, fluxes, dtdia, &
751  time_local, .false., waves=waves)
752  cs%time_in_thermo_cycle = cs%time_in_thermo_cycle + dtdia
753 
754  if ((cs%t_dyn_rel_thermo==0.0) .and. .not.do_dyn) then
755  ! The diabatic processes are now ahead of the dynamics by dtdia.
756  cs%t_dyn_rel_thermo = -dtdia
757  else ! The diabatic processes and the dynamics are synchronized.
758  cs%t_dyn_rel_thermo = 0.0
759  endif
760 
761  if (dtdia > dt) & ! Reset CS%Time to its previous value.
762  cs%Time = time_start + real_to_time(us%T_to_s*(rel_time - 0.5*dt))
763  endif
764 
765  if (do_dyn) then
766  call cpu_clock_begin(id_clock_dynamics)
767  ! Determining the time-average sea surface height is part of the algorithm.
768  ! This may be eta_av if Boussinesq, or need to be diagnosed if not.
769  cs%time_in_cycle = cs%time_in_cycle + dt
770  call find_eta(h, cs%tv, g, gv, us, ssh, cs%eta_av_bc, eta_to_m=1.0)
771  do j=js,je ; do i=is,ie
772  cs%ssh_rint(i,j) = cs%ssh_rint(i,j) + dt*ssh(i,j)
773  enddo ; enddo
774  if (cs%IDs%id_ssh_inst > 0) call post_data(cs%IDs%id_ssh_inst, ssh, cs%diag)
775  call cpu_clock_end(id_clock_dynamics)
776  endif
777 
778  !===========================================================================
779  ! Calculate diagnostics at the end of the time step if the state is self-consistent.
780  if (mom_state_is_synchronized(cs)) then
781  !### Perhaps this should be if (CS%t_dyn_rel_thermo == 0.0)
782  call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
783  ! Diagnostics that require the complete state to be up-to-date can be calculated.
784 
785  call enable_averages(cs%t_dyn_rel_diag, time_local, cs%diag)
786  call calculate_diagnostic_fields(u, v, h, cs%uh, cs%vh, cs%tv, cs%ADp, &
787  cs%CDp, p_surf, cs%t_dyn_rel_diag, cs%diag_pre_sync,&
788  g, gv, us, cs%diagnostics_CSp)
789  call post_tracer_diagnostics(cs%Tracer_reg, h, cs%diag_pre_sync, cs%diag, g, gv, cs%t_dyn_rel_diag)
790  call diag_copy_diag_to_storage(cs%diag_pre_sync, h, cs%diag)
791  if (showcalltree) call calltree_waypoint("finished calculate_diagnostic_fields (step_MOM)")
792  call disable_averaging(cs%diag)
793  cs%t_dyn_rel_diag = 0.0
794 
795  call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
796  endif
797 
798  if (do_dyn .and. .not.cs%count_calls) cs%nstep_tot = cs%nstep_tot + 1
799  if (showcalltree) call calltree_leave("DT cycles (step_MOM)")
800 
801  enddo ! complete the n loop
802 
803  if (cs%count_calls .and. cycle_start) cs%nstep_tot = cs%nstep_tot + 1
804 
805  call cpu_clock_begin(id_clock_other)
806 
807  if (cs%time_in_cycle > 0.0) then
808  i_wt_ssh = 1.0/cs%time_in_cycle
809  do j=js,je ; do i=is,ie
810  ssh(i,j) = cs%ssh_rint(i,j)*i_wt_ssh
811  cs%ave_ssh_ibc(i,j) = ssh(i,j)
812  enddo ; enddo
813  if (do_dyn) then
814  call adjust_ssh_for_p_atm(cs%tv, g, gv, us, cs%ave_ssh_ibc, forces%p_surf_SSH, &
815  cs%calc_rho_for_sea_lev)
816  elseif (do_thermo) then
817  call adjust_ssh_for_p_atm(cs%tv, g, gv, us, cs%ave_ssh_ibc, fluxes%p_surf_SSH, &
818  cs%calc_rho_for_sea_lev)
819  endif
820  endif
821 
822  if (do_dyn .and. cs%interp_p_surf) then ; do j=jsd,jed ; do i=isd,ied
823  cs%p_surf_prev(i,j) = forces%p_surf(i,j)
824  enddo ; enddo ; endif
825 
826  if (cs%ensemble_ocean) then
827  ! update the time for the next analysis step if needed
828  call set_analysis_time(cs%Time,cs%odaCS)
829  ! store ensemble vector in odaCS
830  call set_prior_tracer(cs%Time, g, gv, cs%h, cs%tv, cs%odaCS)
831  ! call DA interface
832  call oda(cs%Time,cs%odaCS)
833  endif
834 
835  if (showcalltree) call calltree_waypoint("calling extract_surface_state (step_MOM)")
836  call extract_surface_state(cs, sfc_state)
837 
838  ! Do diagnostics that only occur at the end of a complete forcing step.
839  if (cycle_end) then
840  call cpu_clock_begin(id_clock_diagnostics)
841  if (cs%time_in_cycle > 0.0) then
842  call enable_averages(cs%time_in_cycle, time_local, cs%diag)
843  call post_surface_dyn_diags(cs%sfc_IDs, g, cs%diag, sfc_state, ssh)
844  endif
845  if (cs%time_in_thermo_cycle > 0.0) then
846  call enable_averages(cs%time_in_thermo_cycle, time_local, cs%diag)
847  call post_surface_thermo_diags(cs%sfc_IDs, g, gv, us, cs%diag, cs%time_in_thermo_cycle, &
848  sfc_state, cs%tv, ssh, cs%ave_ssh_ibc)
849  endif
850  call disable_averaging(cs%diag)
851  call cpu_clock_end(id_clock_diagnostics)
852  endif
853 
854  ! Accumulate the surface fluxes for assessing conservation
855  if (do_thermo .and. fluxes%fluxes_used) &
856  call accumulate_net_input(fluxes, sfc_state, cs%tv, fluxes%dt_buoy_accum, &
857  g, us, cs%sum_output_CSp)
858 
859  if (mom_state_is_synchronized(cs)) &
860  call write_energy(cs%u, cs%v, cs%h, cs%tv, time_local, cs%nstep_tot, &
861  g, gv, us, cs%sum_output_CSp, cs%tracer_flow_CSp, &
862  dt_forcing=real_to_time(us%T_to_s*time_interval) )
863 
864  call cpu_clock_end(id_clock_other)
865 
866  if (showcalltree) call calltree_leave("step_MOM()")
867  call cpu_clock_end(id_clock_ocean)
868 

References adjust_ssh_for_p_atm(), mom_lateral_mixing_coeffs::calc_depth_function(), mom_lateral_mixing_coeffs::calc_resoln_function(), mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), mom_error_handler::calltree_waypoint(), mom_domains::complete_group_pass(), mom_diag_mediator::diag_copy_diag_to_storage(), extract_surface_state(), id_clock_diagnostics, id_clock_dynamics, id_clock_ocean, id_clock_other, id_clock_pass, mom_forcing_type::mom_mech_forcing_chksum(), mom_state_is_synchronized(), mom_oda_driver_mod::set_analysis_time(), mom_oda_driver_mod::set_prior_tracer(), mom_domains::start_group_pass(), step_mom_dynamics(), step_mom_thermo(), step_mom_tracer_dyn(), and mom_wave_interface::update_stokes_drift().

Referenced by mom_main().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ step_mom_dynamics()

subroutine mom::step_mom_dynamics ( type(mech_forcing), intent(in)  forces,
real, dimension(:,:), pointer  p_surf_begin,
real, dimension(:,:), pointer  p_surf_end,
real, intent(in)  dt,
real, intent(in)  dt_thermo,
real, intent(in)  bbl_time_int,
type(mom_control_struct), pointer  CS,
type(time_type), intent(in)  Time_local,
type(wave_parameters_cs), optional, pointer  Waves 
)
private

Time step the ocean dynamics, including the momentum and continuity equations.

Parameters
[in]forcesA structure with the driving mechanical forces
p_surf_beginA pointer (perhaps NULL) to the surface pressure at the beginning of this dynamic step, intent in [Pa].
p_surf_endA pointer (perhaps NULL) to the surface pressure at the end of this dynamic step, intent in [Pa].
[in]dttime interval covered by this call [T ~> s].
[in]dt_thermotime interval covered by any updates that may span multiple dynamics steps [T ~> s].
[in]bbl_time_inttime interval over which updates to the bottom boundary layer properties will apply [T ~> s], or zero not to update the properties.
cscontrol structure from initialize_MOM
[in]time_localEnd time of a segment, as a time type
wavesContainer for wave related parameters; the

Definition at line 874 of file MOM.F90.

874  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
875  real, dimension(:,:), pointer :: p_surf_begin !< A pointer (perhaps NULL) to the surface
876  !! pressure at the beginning of this dynamic
877  !! step, intent in [Pa].
878  real, dimension(:,:), pointer :: p_surf_end !< A pointer (perhaps NULL) to the surface
879  !! pressure at the end of this dynamic step,
880  !! intent in [Pa].
881  real, intent(in) :: dt !< time interval covered by this call [T ~> s].
882  real, intent(in) :: dt_thermo !< time interval covered by any updates that may
883  !! span multiple dynamics steps [T ~> s].
884  real, intent(in) :: bbl_time_int !< time interval over which updates to the
885  !! bottom boundary layer properties will apply [T ~> s],
886  !! or zero not to update the properties.
887  type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM
888  type(time_type), intent(in) :: Time_local !< End time of a segment, as a time type
889  type(wave_parameters_CS), &
890  optional, pointer :: Waves !< Container for wave related parameters; the
891  !! fields in Waves are intent in here.
892 
893  ! local variables
894  type(ocean_grid_type), pointer :: G => null() ! pointer to a structure containing
895  ! metrics and related information
896  type(verticalGrid_type), pointer :: GV => null() ! Pointer to the vertical grid structure
897  type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing
898  ! various unit conversion factors
899  type(MOM_diag_IDs), pointer :: IDs => null() ! A structure with the diagnostic IDs.
900  real, dimension(:,:,:), pointer :: &
901  u => null(), & ! u : zonal velocity component [L T-1 ~> m s-1]
902  v => null(), & ! v : meridional velocity component [L T-1 ~> m s-1]
903  h => null() ! h : layer thickness [H ~> m or kg m-2]
904 
905  logical :: calc_dtbt ! Indicates whether the dynamically adjusted
906  ! barotropic time step needs to be updated.
907  logical :: showCallTree
908 
909  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
910  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
911 
912  g => cs%G ; gv => cs%GV ; us => cs%US ; ids => cs%IDs
913  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
914  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
915  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
916  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
917  u => cs%u ; v => cs%v ; h => cs%h
918  showcalltree = calltree_showquery()
919 
920  call cpu_clock_begin(id_clock_dynamics)
921 
922  if ((cs%t_dyn_rel_adv == 0.0) .and. cs%thickness_diffuse .and. cs%thickness_diffuse_first) then
923 
924  call enable_averages(dt_thermo, time_local+real_to_time(us%T_to_s*(dt_thermo-dt)), cs%diag)
925  call cpu_clock_begin(id_clock_thick_diff)
926  if (associated(cs%VarMix)) &
927  call calc_slope_functions(h, cs%tv, dt, g, gv, us, cs%VarMix)
928  call thickness_diffuse(h, cs%uhtr, cs%vhtr, cs%tv, dt_thermo, g, gv, us, &
929  cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp)
930  call cpu_clock_end(id_clock_thick_diff)
931  call pass_var(h, g%Domain, clock=id_clock_pass) !###, halo=max(2,cont_stensil))
932  call disable_averaging(cs%diag)
933  if (showcalltree) call calltree_waypoint("finished thickness_diffuse_first (step_MOM)")
934 
935  ! Whenever thickness changes let the diag manager know, target grids
936  ! for vertical remapping may need to be regenerated.
937  call diag_update_remap_grids(cs%diag)
938  endif
939 
940  ! The bottom boundary layer properties need to be recalculated.
941  if (bbl_time_int > 0.0) then
942  call enable_averages(bbl_time_int, &
943  time_local + real_to_time(us%T_to_s*(bbl_time_int-dt)), cs%diag)
944  ! Calculate the BBL properties and store them inside visc (u,h).
945  call cpu_clock_begin(id_clock_bbl_visc)
946  call set_viscous_bbl(cs%u(:,:,:), cs%v(:,:,:), cs%h, cs%tv, cs%visc, g, gv, us, &
947  cs%set_visc_CSp, symmetrize=.true.)
948  call cpu_clock_end(id_clock_bbl_visc)
949  if (showcalltree) call calltree_waypoint("done with set_viscous_BBL (step_MOM)")
950  call disable_averaging(cs%diag)
951  endif
952 
953 
954  if (cs%do_dynamics .and. cs%split) then !--------------------------- start SPLIT
955  ! This section uses a split time stepping scheme for the dynamic equations,
956  ! basically the stacked shallow water equations with viscosity.
957 
958  calc_dtbt = .false.
959  if (cs%dtbt_reset_period == 0.0) calc_dtbt = .true.
960  if (cs%dtbt_reset_period > 0.0) then
961  if (time_local >= cs%dtbt_reset_time) then !### Change >= to > here.
962  calc_dtbt = .true.
963  cs%dtbt_reset_time = cs%dtbt_reset_time + cs%dtbt_reset_interval
964  endif
965  endif
966 
967  call step_mom_dyn_split_rk2(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
968  p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
969  cs%eta_av_bc, g, gv, us, cs%dyn_split_RK2_CSp, calc_dtbt, cs%VarMix, &
970  cs%MEKE, cs%thickness_diffuse_CSp, waves=waves)
971  if (showcalltree) call calltree_waypoint("finished step_MOM_dyn_split (step_MOM)")
972 
973  elseif (cs%do_dynamics) then ! ------------------------------------ not SPLIT
974  ! This section uses an unsplit stepping scheme for the dynamic
975  ! equations; basically the stacked shallow water equations with viscosity.
976  ! Because the time step is limited by CFL restrictions on the external
977  ! gravity waves, the unsplit is usually much less efficient that the split
978  ! approaches. But because of its simplicity, the unsplit method is very
979  ! useful for debugging purposes.
980 
981  if (cs%use_RK2) then
982  call step_mom_dyn_unsplit_rk2(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
983  p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
984  cs%eta_av_bc, g, gv, us, cs%dyn_unsplit_RK2_CSp, cs%VarMix, cs%MEKE)
985  else
986  call step_mom_dyn_unsplit(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
987  p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
988  cs%eta_av_bc, g, gv, us, cs%dyn_unsplit_CSp, cs%VarMix, cs%MEKE, waves=waves)
989  endif
990  if (showcalltree) call calltree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)")
991 
992  endif ! -------------------------------------------------- end SPLIT
993 
994  if (cs%thickness_diffuse .and. .not.cs%thickness_diffuse_first) then
995  call cpu_clock_begin(id_clock_thick_diff)
996 
997  if (cs%debug) call hchksum(h,"Pre-thickness_diffuse h", g%HI, haloshift=0, scale=gv%H_to_m)
998 
999  if (associated(cs%VarMix)) &
1000  call calc_slope_functions(h, cs%tv, dt, g, gv, us, cs%VarMix)
1001  call thickness_diffuse(h, cs%uhtr, cs%vhtr, cs%tv, dt, g, gv, us, &
1002  cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp)
1003 
1004  if (cs%debug) call hchksum(h,"Post-thickness_diffuse h", g%HI, haloshift=1, scale=gv%H_to_m)
1005  call cpu_clock_end(id_clock_thick_diff)
1006  call pass_var(h, g%Domain, clock=id_clock_pass) !###, halo=max(2,cont_stensil))
1007  if (showcalltree) call calltree_waypoint("finished thickness_diffuse (step_MOM)")
1008  endif
1009 
1010  ! apply the submesoscale mixed layer restratification parameterization
1011  if (cs%mixedlayer_restrat) then
1012  if (cs%debug) then
1013  call hchksum(h,"Pre-mixedlayer_restrat h", g%HI, haloshift=1, scale=gv%H_to_m)
1014  call uvchksum("Pre-mixedlayer_restrat uhtr", &
1015  cs%uhtr, cs%vhtr, g%HI, haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1016  endif
1017  call cpu_clock_begin(id_clock_ml_restrat)
1018  call mixedlayer_restrat(h, cs%uhtr, cs%vhtr, cs%tv, forces, dt, cs%visc%MLD, &
1019  cs%VarMix, g, gv, us, cs%mixedlayer_restrat_CSp)
1020  call cpu_clock_end(id_clock_ml_restrat)
1021  call pass_var(h, g%Domain, clock=id_clock_pass) !###, halo=max(2,cont_stensil))
1022  if (cs%debug) then
1023  call hchksum(h,"Post-mixedlayer_restrat h", g%HI, haloshift=1, scale=gv%H_to_m)
1024  call uvchksum("Post-mixedlayer_restrat [uv]htr", &
1025  cs%uhtr, cs%vhtr, g%HI, haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1026  endif
1027  endif
1028 
1029  ! Whenever thickness changes let the diag manager know, target grids
1030  ! for vertical remapping may need to be regenerated.
1031  call diag_update_remap_grids(cs%diag)
1032 
1033  if (cs%useMEKE) call step_forward_meke(cs%MEKE, h, cs%VarMix%SN_u, cs%VarMix%SN_v, &
1034  cs%visc, dt, g, gv, us, cs%MEKE_CSp, cs%uhtr, cs%vhtr)
1035  call disable_averaging(cs%diag)
1036 
1037  ! Advance the dynamics time by dt.
1038  cs%t_dyn_rel_adv = cs%t_dyn_rel_adv + dt
1039  cs%t_dyn_rel_thermo = cs%t_dyn_rel_thermo + dt
1040  if (abs(cs%t_dyn_rel_thermo) < 1e-6*dt) cs%t_dyn_rel_thermo = 0.0
1041  cs%t_dyn_rel_diag = cs%t_dyn_rel_diag + dt
1042 
1043  call cpu_clock_end(id_clock_dynamics)
1044 
1045  call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
1046  call enable_averages(dt, time_local, cs%diag)
1047  ! These diagnostics are available after every time dynamics step.
1048  if (ids%id_u > 0) call post_data(ids%id_u, u, cs%diag)
1049  if (ids%id_v > 0) call post_data(ids%id_v, v, cs%diag)
1050  if (ids%id_h > 0) call post_data(ids%id_h, h, cs%diag)
1051  call disable_averaging(cs%diag)
1052  call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
1053 

References mom_error_handler::calltree_waypoint(), id_clock_bbl_visc, id_clock_diagnostics, id_clock_dynamics, id_clock_ml_restrat, id_clock_other, id_clock_pass, id_clock_thick_diff, mom_meke::step_forward_meke(), and mom_thickness_diffuse::thickness_diffuse().

Referenced by step_mom().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ step_mom_thermo()

subroutine mom::step_mom_thermo ( type(mom_control_struct), intent(inout)  CS,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(inout)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(inout)  u,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)  h,
type(thermo_var_ptrs), intent(inout)  tv,
type(forcing), intent(inout)  fluxes,
real, intent(in)  dtdia,
type(time_type), intent(in)  Time_end_thermo,
logical, intent(in)  update_BBL,
type(wave_parameters_cs), optional, pointer  Waves 
)
private

MOM_step_thermo orchestrates the thermodynamic time stepping and vertical remapping, via calls to diabatic (or adiabatic) and ALE_main.

Parameters
[in,out]csMaster MOM control structure
[in,out]gocean grid structure
[in,out]gvocean vertical grid structure
[in]usA dimensional unit scaling type
[in,out]uzonal velocity [L T-1 ~> m s-1]
[in,out]vmeridional velocity [L T-1 ~> m s-1]
[in,out]hlayer thickness [H ~> m or kg m-2]
[in,out]tvA structure pointing to various thermodynamic variables
[in,out]fluxespointers to forcing fields
[in]dtdiaThe time interval over which to advance [T ~> s]
[in]time_end_thermoEnd of averaging interval for thermo diags
[in]update_bblIf true, calculate the bottom boundary layer properties.
wavesContainer for wave related parameters

Definition at line 1131 of file MOM.F90.

1131  type(MOM_control_struct), intent(inout) :: CS !< Master MOM control structure
1132  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1133  type(verticalGrid_type), intent(inout) :: GV !< ocean vertical grid structure
1134  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1135  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1136  intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1137  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1138  intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
1139  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1140  intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
1141  type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
1142  type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
1143  real, intent(in) :: dtdia !< The time interval over which to advance [T ~> s]
1144  type(time_type), intent(in) :: Time_end_thermo !< End of averaging interval for thermo diags
1145  logical, intent(in) :: update_BBL !< If true, calculate the bottom boundary layer properties.
1146  type(wave_parameters_CS), &
1147  optional, pointer :: Waves !< Container for wave related parameters
1148  !! the fields in Waves are intent in here.
1149 
1150  logical :: use_ice_shelf ! Needed for selecting the right ALE interface.
1151  logical :: showCallTree
1152  type(group_pass_type) :: pass_T_S, pass_T_S_h, pass_uv_T_S_h
1153  integer :: dynamics_stencil ! The computational stencil for the calculations
1154  ! in the dynamic core.
1155  integer :: i, j, k, is, ie, js, je, nz! , Isq, Ieq, Jsq, Jeq, n
1156 
1157  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1158  showcalltree = calltree_showquery()
1159  if (showcalltree) call calltree_enter("step_MOM_thermo(), MOM.F90")
1160 
1161  use_ice_shelf = .false.
1162  if (associated(fluxes%frac_shelf_h)) use_ice_shelf = .true.
1163 
1164  call enable_averages(dtdia, time_end_thermo, cs%diag)
1165 
1166  if (associated(cs%odaCS)) then
1167  call apply_oda_tracer_increments(us%T_to_s*dtdia,g,tv,h,cs%odaCS)
1168  endif
1169 
1170  if (update_bbl) then
1171  ! Calculate the BBL properties and store them inside visc (u,h).
1172  ! This is here so that CS%visc is updated before diabatic() when
1173  ! DIABATIC_FIRST=True. Otherwise diabatic() is called after the dynamics
1174  ! and set_viscous_BBL is called as a part of the dynamic stepping.
1175  call cpu_clock_begin(id_clock_bbl_visc)
1176  call set_viscous_bbl(u, v, h, tv, cs%visc, g, gv, us, cs%set_visc_CSp, symmetrize=.true.)
1177  call cpu_clock_end(id_clock_bbl_visc)
1178  if (showcalltree) call calltree_waypoint("done with set_viscous_BBL (step_MOM_thermo)")
1179  endif
1180 
1181  call cpu_clock_begin(id_clock_thermo)
1182  if (.not.cs%adiabatic) then
1183  if (cs%debug) then
1184  call uvchksum("Pre-diabatic [uv]", u, v, g%HI, haloshift=2, scale=us%L_T_to_m_s)
1185  call hchksum(h,"Pre-diabatic h", g%HI, haloshift=1, scale=gv%H_to_m)
1186  call uvchksum("Pre-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1187  haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1188  ! call MOM_state_chksum("Pre-diabatic ", u, v, h, CS%uhtr, CS%vhtr, G, GV, vel_scale=1.0)
1189  call mom_thermo_chksum("Pre-diabatic ", tv, g, us, haloshift=0)
1190  call check_redundant("Pre-diabatic ", u, v, g)
1191  call mom_forcing_chksum("Pre-diabatic", fluxes, g, us, haloshift=0)
1192  endif
1193 
1194  call cpu_clock_begin(id_clock_diabatic)
1195 
1196  call diabatic(u, v, h, tv, cs%Hml, fluxes, cs%visc, cs%ADp, cs%CDp, &
1197  dtdia, time_end_thermo, g, gv, us, cs%diabatic_CSp, waves=waves)
1198  fluxes%fluxes_used = .true.
1199 
1200  if (showcalltree) call calltree_waypoint("finished diabatic (step_MOM_thermo)")
1201 
1202  ! Regridding/remapping is done here, at end of thermodynamics time step
1203  ! (that may comprise several dynamical time steps)
1204  ! The routine 'ALE_main' can be found in 'MOM_ALE.F90'.
1205  if ( cs%use_ALE_algorithm ) then
1206  call enable_averages(dtdia, time_end_thermo, cs%diag)
1207 ! call pass_vector(u, v, G%Domain)
1208  if (associated(tv%T)) &
1209  call create_group_pass(pass_t_s_h, tv%T, g%Domain, to_all+omit_corners, halo=1)
1210  if (associated(tv%S)) &
1211  call create_group_pass(pass_t_s_h, tv%S, g%Domain, to_all+omit_corners, halo=1)
1212  call create_group_pass(pass_t_s_h, h, g%Domain, to_all+omit_corners, halo=1)
1213  call do_group_pass(pass_t_s_h, g%Domain)
1214 
1215  call preale_tracer_diagnostics(cs%tracer_Reg, g, gv)
1216 
1217  if (cs%debug) then
1218  call mom_state_chksum("Pre-ALE ", u, v, h, cs%uh, cs%vh, g, gv, us)
1219  call hchksum(tv%T,"Pre-ALE T", g%HI, haloshift=1)
1220  call hchksum(tv%S,"Pre-ALE S", g%HI, haloshift=1)
1221  call check_redundant("Pre-ALE ", u, v, g)
1222  endif
1223  call cpu_clock_begin(id_clock_ale)
1224  if (use_ice_shelf) then
1225  call ale_main(g, gv, us, h, u, v, tv, cs%tracer_Reg, cs%ALE_CSp, cs%OBC, &
1226  dtdia, fluxes%frac_shelf_h)
1227  else
1228  call ale_main(g, gv, us, h, u, v, tv, cs%tracer_Reg, cs%ALE_CSp, cs%OBC, dtdia)
1229  endif
1230 
1231  if (showcalltree) call calltree_waypoint("finished ALE_main (step_MOM_thermo)")
1232  call cpu_clock_end(id_clock_ale)
1233  endif ! endif for the block "if ( CS%use_ALE_algorithm )"
1234 
1235  dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
1236  call create_group_pass(pass_uv_t_s_h, u, v, g%Domain, halo=dynamics_stencil)
1237  if (associated(tv%T)) &
1238  call create_group_pass(pass_uv_t_s_h, tv%T, g%Domain, halo=dynamics_stencil)
1239  if (associated(tv%S)) &
1240  call create_group_pass(pass_uv_t_s_h, tv%S, g%Domain, halo=dynamics_stencil)
1241  call create_group_pass(pass_uv_t_s_h, h, g%Domain, halo=dynamics_stencil)
1242  call do_group_pass(pass_uv_t_s_h, g%Domain, clock=id_clock_pass)
1243 
1244  if (cs%debug .and. cs%use_ALE_algorithm) then
1245  call mom_state_chksum("Post-ALE ", u, v, h, cs%uh, cs%vh, g, gv, us)
1246  call hchksum(tv%T, "Post-ALE T", g%HI, haloshift=1)
1247  call hchksum(tv%S, "Post-ALE S", g%HI, haloshift=1)
1248  call check_redundant("Post-ALE ", u, v, g)
1249  endif
1250 
1251  ! Whenever thickness changes let the diag manager know, target grids
1252  ! for vertical remapping may need to be regenerated. This needs to
1253  ! happen after the H update and before the next post_data.
1254  call diag_update_remap_grids(cs%diag)
1255 
1256  !### Consider moving this up into the if ALE block.
1257  call postale_tracer_diagnostics(cs%tracer_Reg, g, gv, cs%diag, dtdia)
1258 
1259  if (cs%debug) then
1260  call uvchksum("Post-diabatic u", u, v, g%HI, haloshift=2, scale=us%L_T_to_m_s)
1261  call hchksum(h, "Post-diabatic h", g%HI, haloshift=1, scale=gv%H_to_m)
1262  call uvchksum("Post-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1263  haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1264  ! call MOM_state_chksum("Post-diabatic ", u, v, &
1265  ! h, CS%uhtr, CS%vhtr, G, GV, haloshift=1)
1266  if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", g%HI, haloshift=1)
1267  if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", g%HI, haloshift=1)
1268  if (associated(tv%frazil)) call hchksum(tv%frazil, &
1269  "Post-diabatic frazil", g%HI, haloshift=0)
1270  if (associated(tv%salt_deficit)) call hchksum(tv%salt_deficit, &
1271  "Post-diabatic salt deficit", g%HI, haloshift=0, scale=us%R_to_kg_m3*us%Z_to_m)
1272  ! call MOM_thermo_chksum("Post-diabatic ", tv, G, US)
1273  call check_redundant("Post-diabatic ", u, v, g)
1274  endif
1275  call disable_averaging(cs%diag)
1276 
1277  call cpu_clock_end(id_clock_diabatic)
1278  else ! complement of "if (.not.CS%adiabatic)"
1279 
1280  call cpu_clock_begin(id_clock_adiabatic)
1281  call adiabatic(h, tv, fluxes, dtdia, g, gv, us, cs%diabatic_CSp)
1282  fluxes%fluxes_used = .true.
1283  call cpu_clock_end(id_clock_adiabatic)
1284 
1285  if (associated(tv%T)) then
1286  call create_group_pass(pass_t_s, tv%T, g%Domain, to_all+omit_corners, halo=1)
1287  call create_group_pass(pass_t_s, tv%S, g%Domain, to_all+omit_corners, halo=1)
1288  call do_group_pass(pass_t_s, g%Domain, clock=id_clock_pass)
1289  if (cs%debug) then
1290  if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", g%HI, haloshift=1)
1291  if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", g%HI, haloshift=1)
1292  endif
1293  endif
1294 
1295  endif ! endif for the block "if (.not.CS%adiabatic)"
1296  call cpu_clock_end(id_clock_thermo)
1297 
1298  call disable_averaging(cs%diag)
1299 
1300  if (showcalltree) call calltree_leave("step_MOM_thermo(), MOM.F90")
1301 

References mom_diabatic_driver::adiabatic(), mom_oda_driver_mod::apply_oda_tracer_increments(), mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), mom_error_handler::calltree_waypoint(), id_clock_adiabatic, id_clock_ale, id_clock_bbl_visc, id_clock_diabatic, id_clock_pass, id_clock_thermo, and mom_forcing_type::mom_forcing_chksum().

Referenced by step_mom().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ step_mom_tracer_dyn()

subroutine mom::step_mom_tracer_dyn ( type(mom_control_struct), intent(inout)  CS,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
type(time_type), intent(in)  Time_local 
)
private

step_MOM_tracer_dyn does tracer advection and lateral diffusion, bringing the tracers up to date with the changes in state due to the dynamics. Surface sources and sinks and remapping are handled via step_MOM_thermo.

Parameters
[in,out]cscontrol structure
[in,out]gocean grid structure
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
[in]hlayer thicknesses after the transports [H ~> m or kg m-2]
[in]time_localThe model time at the end of the time step.

Definition at line 1060 of file MOM.F90.

1060  type(MOM_control_struct), intent(inout) :: CS !< control structure
1061  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1062  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1063  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1064  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1065  intent(in) :: h !< layer thicknesses after the transports [H ~> m or kg m-2]
1066  type(time_type), intent(in) :: Time_local !< The model time at the end
1067  !! of the time step.
1068  type(group_pass_type) :: pass_T_S
1069  logical :: showCallTree
1070  showcalltree = calltree_showquery()
1071 
1072  if (cs%debug) then
1073  call cpu_clock_begin(id_clock_other)
1074  call hchksum(h,"Pre-advection h", g%HI, haloshift=1, scale=gv%H_to_m)
1075  call uvchksum("Pre-advection uhtr", cs%uhtr, cs%vhtr, g%HI, &
1076  haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1077  if (associated(cs%tv%T)) call hchksum(cs%tv%T, "Pre-advection T", g%HI, haloshift=1)
1078  if (associated(cs%tv%S)) call hchksum(cs%tv%S, "Pre-advection S", g%HI, haloshift=1)
1079  if (associated(cs%tv%frazil)) call hchksum(cs%tv%frazil, &
1080  "Pre-advection frazil", g%HI, haloshift=0)
1081  if (associated(cs%tv%salt_deficit)) call hchksum(cs%tv%salt_deficit, &
1082  "Pre-advection salt deficit", g%HI, haloshift=0, scale=us%R_to_kg_m3*us%Z_to_m)
1083  ! call MOM_thermo_chksum("Pre-advection ", CS%tv, G, US)
1084  call cpu_clock_end(id_clock_other)
1085  endif
1086 
1087  call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer)
1088  call enable_averages(cs%t_dyn_rel_adv, time_local, cs%diag)
1089 
1090  call advect_tracer(h, cs%uhtr, cs%vhtr, cs%OBC, cs%t_dyn_rel_adv, g, gv, us, &
1091  cs%tracer_adv_CSp, cs%tracer_Reg)
1092  call tracer_hordiff(h, cs%t_dyn_rel_adv, cs%MEKE, cs%VarMix, g, gv, us, &
1093  cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1094  if (showcalltree) call calltree_waypoint("finished tracer advection/diffusion (step_MOM)")
1095  call update_segment_tracer_reservoirs(g, gv, cs%uhtr, cs%vhtr, h, cs%OBC, &
1096  cs%t_dyn_rel_adv, cs%tracer_Reg)
1097  call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo)
1098 
1099  call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
1100  call post_transport_diagnostics(g, gv, us, cs%uhtr, cs%vhtr, h, cs%transport_IDs, &
1101  cs%diag_pre_dyn, cs%diag, cs%t_dyn_rel_adv, cs%tracer_reg)
1102  ! Rebuild the remap grids now that we've posted the fields which rely on thicknesses
1103  ! from before the dynamics calls
1104  call diag_update_remap_grids(cs%diag)
1105 
1106  call disable_averaging(cs%diag)
1107  call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
1108 
1109  ! Reset the accumulated transports to 0 and record that the dynamics
1110  ! and advective times now agree.
1111  call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer)
1112  cs%uhtr(:,:,:) = 0.0
1113  cs%vhtr(:,:,:) = 0.0
1114  cs%t_dyn_rel_adv = 0.0
1115  call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo)
1116 
1117  if (cs%diabatic_first .and. associated(cs%tv%T)) then
1118  ! Temperature and salinity need halo updates because they will be used
1119  ! in the dynamics before they are changed again.
1120  call create_group_pass(pass_t_s, cs%tv%T, g%Domain, to_all+omit_corners, halo=1)
1121  call create_group_pass(pass_t_s, cs%tv%S, g%Domain, to_all+omit_corners, halo=1)
1122  call do_group_pass(pass_t_s, g%Domain, clock=id_clock_pass)
1123  endif
1124 

References mom_error_handler::calltree_waypoint(), id_clock_diagnostics, id_clock_other, id_clock_pass, id_clock_thermo, id_clock_tracer, and mom_open_boundary::update_segment_tracer_reservoirs().

Referenced by step_mom().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ step_offline()

subroutine, public mom::step_offline ( type(mech_forcing), intent(in)  forces,
type(forcing), intent(inout)  fluxes,
type(surface), intent(inout)  sfc_state,
type(time_type), intent(in)  Time_start,
real, intent(in)  time_interval,
type(mom_control_struct), pointer  CS 
)

step_offline is the main driver for running tracers offline in MOM6. This has been primarily developed with ALE configurations in mind. Some work has been done in isopycnal configuration, but the work is very preliminary. Some more detail about this capability along with some of the subroutines called here can be found in tracers/MOM_offline_control.F90

Parameters
[in]forcesA structure with the driving mechanical forces
[in,out]fluxespointers to forcing fields
[in,out]sfc_statesurface ocean state
[in]time_startstarting time of a segment, as a time type
[in]time_intervaltime interval
cscontrol structure from initialize_MOM

Definition at line 1310 of file MOM.F90.

1310  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1311  type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
1312  type(surface), intent(inout) :: sfc_state !< surface ocean state
1313  type(time_type), intent(in) :: Time_start !< starting time of a segment, as a time type
1314  real, intent(in) :: time_interval !< time interval
1315  type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM
1316 
1317  ! Local pointers
1318  type(ocean_grid_type), pointer :: G => null() ! Pointer to a structure containing
1319  ! metrics and related information
1320  type(verticalGrid_type), pointer :: GV => null() ! Pointer to structure containing information
1321  ! about the vertical grid
1322  type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing
1323  ! various unit conversion factors
1324 
1325  logical :: first_iter !< True if this is the first time step_offline has been called in a given interval
1326  logical :: last_iter !< True if this is the last time step_tracer is to be called in an offline interval
1327  logical :: do_vertical !< If enough time has elapsed, do the diabatic tracer sources/sinks
1328  logical :: adv_converged !< True if all the horizontal fluxes have been used
1329 
1330  real :: dt_offline ! The offline timestep for advection [T ~> s]
1331  real :: dt_offline_vertical ! The offline timestep for vertical fluxes and remapping [T ~> s]
1332  logical :: skip_diffusion
1333  integer :: id_eta_diff_end
1334 
1335  integer, pointer :: accumulated_time => null()
1336  integer :: i,j,k
1337  integer :: is, ie, js, je, isd, ied, jsd, jed
1338 
1339  ! 3D pointers
1340  real, dimension(:,:,:), pointer :: &
1341  uhtr => null(), vhtr => null(), &
1342  eatr => null(), ebtr => null(), &
1343  h_end => null()
1344 
1345  ! 2D Array for diagnostics
1346  real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: eta_pre, eta_end
1347  type(time_type) :: Time_end ! End time of a segment, as a time type
1348 
1349  ! Grid-related pointer assignments
1350  g => cs%G ; gv => cs%GV ; us => cs%US
1351 
1352  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1353  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1354 
1355  call cpu_clock_begin(id_clock_offline_tracer)
1356  call extract_offline_main(cs%offline_CSp, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, &
1357  dt_offline, dt_offline_vertical, skip_diffusion)
1358  time_end = increment_date(time_start, seconds=floor(time_interval+0.001))
1359 
1360  call enable_averaging(time_interval, time_end, cs%diag)
1361 
1362  ! Check to see if this is the first iteration of the offline interval
1363  if (accumulated_time==0) then
1364  first_iter = .true.
1365  else ! This is probably unnecessary but is used to guard against unwanted behavior
1366  first_iter = .false.
1367  endif
1368 
1369  ! Check to see if vertical tracer functions should be done
1370  if ( mod(accumulated_time, floor(us%T_to_s*dt_offline_vertical + 1e-6)) == 0 ) then
1371  do_vertical = .true.
1372  else
1373  do_vertical = .false.
1374  endif
1375 
1376  ! Increment the amount of time elapsed since last read and check if it's time to roll around
1377  accumulated_time = mod(accumulated_time + int(time_interval), floor(us%T_to_s*dt_offline+1e-6))
1378  if (accumulated_time==0) then
1379  last_iter = .true.
1380  else
1381  last_iter = .false.
1382  endif
1383 
1384  if (cs%use_ALE_algorithm) then
1385  ! If this is the first iteration in the offline timestep, then we need to read in fields and
1386  ! perform the main advection.
1387  if (first_iter) then
1388  call mom_mesg("Reading in new offline fields")
1389  ! Read in new transport and other fields
1390  ! call update_transport_from_files(G, GV, CS%offline_CSp, h_end, eatr, ebtr, uhtr, vhtr, &
1391  ! CS%tv%T, CS%tv%S, fluxes, CS%use_ALE_algorithm)
1392  ! call update_transport_from_arrays(CS%offline_CSp)
1393  call update_offline_fields(cs%offline_CSp, cs%h, fluxes, cs%use_ALE_algorithm)
1394 
1395  ! Apply any fluxes into the ocean
1396  call offline_fw_fluxes_into_ocean(g, gv, cs%offline_CSp, fluxes, cs%h)
1397 
1398  if (.not.cs%diabatic_first) then
1399  call offline_advection_ale(fluxes, time_start, time_interval, cs%offline_CSp, id_clock_ale, &
1400  cs%h, uhtr, vhtr, converged=adv_converged)
1401 
1402  ! Redistribute any remaining transport
1403  call offline_redistribute_residual(cs%offline_CSp, cs%h, uhtr, vhtr, adv_converged)
1404 
1405  ! Perform offline diffusion if requested
1406  if (.not. skip_diffusion) then
1407  if (associated(cs%VarMix)) then
1408  call pass_var(cs%h, g%Domain)
1409  call calc_resoln_function(cs%h, cs%tv, g, gv, us, cs%VarMix)
1410  call calc_depth_function(g, cs%VarMix)
1411  call calc_slope_functions(cs%h, cs%tv, dt_offline, g, gv, us, cs%VarMix)
1412  endif
1413  call tracer_hordiff(cs%h, dt_offline, cs%MEKE, cs%VarMix, g, gv, us, &
1414  cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1415  endif
1416  endif
1417  endif
1418  ! The functions related to column physics of tracers is performed separately in ALE mode
1419  if (do_vertical) then
1420  call offline_diabatic_ale(fluxes, time_start, time_end, cs%offline_CSp, cs%h, eatr, ebtr)
1421  endif
1422 
1423  ! Last thing that needs to be done is the final ALE remapping
1424  if (last_iter) then
1425  if (cs%diabatic_first) then
1426  call offline_advection_ale(fluxes, time_start, time_interval, cs%offline_CSp, id_clock_ale, &
1427  cs%h, uhtr, vhtr, converged=adv_converged)
1428 
1429  ! Redistribute any remaining transport and perform the remaining advection
1430  call offline_redistribute_residual(cs%offline_CSp, cs%h, uhtr, vhtr, adv_converged)
1431  ! Perform offline diffusion if requested
1432  if (.not. skip_diffusion) then
1433  if (associated(cs%VarMix)) then
1434  call pass_var(cs%h, g%Domain)
1435  call calc_resoln_function(cs%h, cs%tv, g, gv, us, cs%VarMix)
1436  call calc_depth_function(g, cs%VarMix)
1437  call calc_slope_functions(cs%h, cs%tv, dt_offline, g, gv, us, cs%VarMix)
1438  endif
1439  call tracer_hordiff(cs%h, dt_offline, cs%MEKE, cs%VarMix, g, gv, us, &
1440  cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1441  endif
1442  endif
1443 
1444  call mom_mesg("Last iteration of offline interval")
1445 
1446  ! Apply freshwater fluxes out of the ocean
1447  call offline_fw_fluxes_out_ocean(g, gv, cs%offline_CSp, fluxes, cs%h)
1448  ! These diagnostic can be used to identify which grid points did not converge within
1449  ! the specified number of advection sub iterations
1450  call post_offline_convergence_diags(cs%offline_CSp, cs%h, h_end, uhtr, vhtr)
1451 
1452  ! Call ALE one last time to make sure that tracers are remapped onto the layer thicknesses
1453  ! stored from the forward run
1454  call cpu_clock_begin(id_clock_ale)
1455  call ale_offline_tracer_final( g, gv, cs%h, cs%tv, h_end, cs%tracer_Reg, cs%ALE_CSp, cs%OBC)
1456  call cpu_clock_end(id_clock_ale)
1457  call pass_var(cs%h, g%Domain)
1458  endif
1459  else ! NON-ALE MODE...NOT WELL TESTED
1460  call mom_error(warning, &
1461  "Offline tracer mode in non-ALE configuration has not been thoroughly tested")
1462  ! Note that for the layer mode case, the calls to tracer sources and sinks is embedded in
1463  ! main_offline_advection_layer. Warning: this may not be appropriate for tracers that
1464  ! exchange with the atmosphere
1465  if (abs(time_interval - us%T_to_s*dt_offline) > 1.0e-6) then
1466  call mom_error(fatal, &
1467  "For offline tracer mode in a non-ALE configuration, dt_offline must equal time_interval")
1468  endif
1469  call update_offline_fields(cs%offline_CSp, cs%h, fluxes, cs%use_ALE_algorithm)
1470  call offline_advection_layer(fluxes, time_start, time_interval, cs%offline_CSp, &
1471  cs%h, eatr, ebtr, uhtr, vhtr)
1472  ! Perform offline diffusion if requested
1473  if (.not. skip_diffusion) then
1474  call tracer_hordiff(h_end, dt_offline, cs%MEKE, cs%VarMix, g, gv, us, &
1475  cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1476  endif
1477 
1478  cs%h = h_end
1479 
1480  call pass_var(cs%tv%T, g%Domain)
1481  call pass_var(cs%tv%S, g%Domain)
1482  call pass_var(cs%h, g%Domain)
1483 
1484  endif
1485 
1486  call adjust_ssh_for_p_atm(cs%tv, g, gv, us, cs%ave_ssh_ibc, forces%p_surf_SSH, &
1487  cs%calc_rho_for_sea_lev)
1488  call extract_surface_state(cs, sfc_state)
1489 
1490  call disable_averaging(cs%diag)
1491  call pass_var(cs%tv%T, g%Domain)
1492  call pass_var(cs%tv%S, g%Domain)
1493  call pass_var(cs%h, g%Domain)
1494 
1495  fluxes%fluxes_used = .true.
1496 
1497  call cpu_clock_end(id_clock_offline_tracer)
1498 

References adjust_ssh_for_p_atm(), mom_ale::ale_offline_tracer_final(), mom_lateral_mixing_coeffs::calc_depth_function(), mom_lateral_mixing_coeffs::calc_resoln_function(), extract_surface_state(), id_clock_ale, id_clock_offline_tracer, and mom_offline_main::offline_advection_layer().

Referenced by mom_main(), mom_ocean_model_mct::update_ocean_model(), and mom_ocean_model_nuopc::update_ocean_model().

Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ id_clock_adiabatic

integer mom::id_clock_adiabatic
private

CPU time clock IDs.

Definition at line 373 of file MOM.F90.

373 integer :: id_clock_adiabatic

Referenced by mom_timing_init(), and step_mom_thermo().

◆ id_clock_ale

integer mom::id_clock_ale
private

CPU time clock IDs.

Definition at line 384 of file MOM.F90.

384 integer :: id_clock_ALE

Referenced by mom_timing_init(), step_mom_thermo(), and step_offline().

◆ id_clock_bbl_visc

integer mom::id_clock_bbl_visc
private

CPU time clock IDs.

Definition at line 376 of file MOM.F90.

376 integer :: id_clock_BBL_visc

Referenced by mom_timing_init(), step_mom_dynamics(), and step_mom_thermo().

◆ id_clock_continuity

integer mom::id_clock_continuity
private

CPU time clock IDs.

Definition at line 374 of file MOM.F90.

374 integer :: id_clock_continuity ! also in dynamics s/r

Referenced by mom_timing_init().

◆ id_clock_diabatic

integer mom::id_clock_diabatic
private

CPU time clock IDs.

Definition at line 372 of file MOM.F90.

372 integer :: id_clock_diabatic

Referenced by mom_timing_init(), and step_mom_thermo().

◆ id_clock_diagnostics

integer mom::id_clock_diagnostics
private

CPU time clock IDs.

Definition at line 378 of file MOM.F90.

378 integer :: id_clock_diagnostics

Referenced by mom_timing_init(), step_mom(), step_mom_dynamics(), and step_mom_tracer_dyn().

◆ id_clock_dynamics

integer mom::id_clock_dynamics
private

CPU time clock IDs.

Definition at line 369 of file MOM.F90.

369 integer :: id_clock_dynamics

Referenced by mom_timing_init(), step_mom(), and step_mom_dynamics().

◆ id_clock_init

integer mom::id_clock_init
private

CPU time clock IDs.

Definition at line 380 of file MOM.F90.

380 integer :: id_clock_init

Referenced by finish_mom_initialization(), and initialize_mom().

◆ id_clock_ml_restrat

integer mom::id_clock_ml_restrat
private

CPU time clock IDs.

Definition at line 377 of file MOM.F90.

377 integer :: id_clock_ml_restrat

Referenced by mom_timing_init(), and step_mom_dynamics().

◆ id_clock_mom_init

integer mom::id_clock_mom_init
private

CPU time clock IDs.

Definition at line 381 of file MOM.F90.

381 integer :: id_clock_MOM_init

Referenced by initialize_mom(), and mom_timing_init().

◆ id_clock_ocean

integer mom::id_clock_ocean

CPU time clock IDs.

Definition at line 368 of file MOM.F90.

368 integer :: id_clock_ocean

Referenced by mom_timing_init(), and step_mom().

◆ id_clock_offline_tracer

integer mom::id_clock_offline_tracer
private

CPU time clock IDs.

Definition at line 386 of file MOM.F90.

386 integer :: id_clock_offline_tracer

Referenced by mom_timing_init(), and step_offline().

◆ id_clock_other

integer mom::id_clock_other
private

CPU time clock IDs.

Definition at line 385 of file MOM.F90.

385 integer :: id_clock_other

Referenced by mom_timing_init(), step_mom(), step_mom_dynamics(), and step_mom_tracer_dyn().

◆ id_clock_pass

integer mom::id_clock_pass
private

CPU time clock IDs.

Definition at line 382 of file MOM.F90.

382 integer :: id_clock_pass ! also in dynamics d/r

Referenced by mom_timing_init(), step_mom(), step_mom_dynamics(), step_mom_thermo(), and step_mom_tracer_dyn().

◆ id_clock_pass_init

integer mom::id_clock_pass_init
private

CPU time clock IDs.

Definition at line 383 of file MOM.F90.

383 integer :: id_clock_pass_init ! also in dynamics d/r

Referenced by initialize_mom(), and mom_timing_init().

◆ id_clock_thermo

integer mom::id_clock_thermo
private

CPU time clock IDs.

Definition at line 370 of file MOM.F90.

370 integer :: id_clock_thermo

Referenced by mom_timing_init(), step_mom_thermo(), and step_mom_tracer_dyn().

◆ id_clock_thick_diff

integer mom::id_clock_thick_diff
private

CPU time clock IDs.

Definition at line 375 of file MOM.F90.

375 integer :: id_clock_thick_diff

Referenced by mom_timing_init(), and step_mom_dynamics().

◆ id_clock_tracer

integer mom::id_clock_tracer
private

CPU time clock IDs.

Definition at line 371 of file MOM.F90.

371 integer :: id_clock_tracer

Referenced by mom_timing_init(), and step_mom_tracer_dyn().

◆ id_clock_z_diag

integer mom::id_clock_z_diag
private

CPU time clock IDs.

Definition at line 379 of file MOM.F90.

379 integer :: id_clock_Z_diag

Referenced by mom_timing_init().