MOM6
Data Types | Functions/Subroutines
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...
 

Functions/Subroutines

subroutine, public step_mom (forces_in, fluxes_in, 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_in)
 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...
 
subroutine rotate_initial_state (u_in, v_in, h_in, T_in, S_in, use_temperature, turns, u, v, h, T, S)
 Rotate initialization fields from input to rotated arrays.
 
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, C_p_scaled, 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_atmOcean surface pressure [R L2 T-2 ~> Pa]
[in]use_eosIf true, calculate the density for the SSH correction using the equation of state.

Definition at line 2964 of file MOM.F90.

2964  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
2965  type(ocean_grid_type), intent(in) :: G !< ocean grid structure
2966  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2967  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2968  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: ssh !< time mean surface height [m]
2969  real, dimension(:,:), optional, pointer :: p_atm !< Ocean surface pressure [R L2 T-2 ~> Pa]
2970  logical, optional, intent(in) :: use_EOS !< If true, calculate the density for
2971  !! the SSH correction using the equation of state.
2972 
2973  real :: Rho_conv(SZI_(G)) ! The density used to convert surface pressure to
2974  ! a corrected effective SSH [R ~> kg m-3].
2975  real :: IgR0 ! The SSH conversion factor from R L2 T-2 to m [m T2 R-1 L-2 ~> m Pa-1].
2976  logical :: calc_rho
2977  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
2978  integer :: i, j, is, ie, js, je
2979 
2980  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2981  eosdom(:) = eos_domain(g%HI)
2982  if (present(p_atm)) then ; if (associated(p_atm)) then
2983  calc_rho = associated(tv%eqn_of_state)
2984  if (present(use_eos) .and. calc_rho) calc_rho = use_eos
2985  ! Correct the output sea surface height for the contribution from the ice pressure.
2986  do j=js,je
2987  if (calc_rho) then
2988  call calculate_density(tv%T(:,j,1), tv%S(:,j,1), 0.5*p_atm(:,j), rho_conv, &
2989  tv%eqn_of_state, eosdom)
2990  do i=is,ie
2991  igr0 = us%Z_to_m / (rho_conv(i) * gv%g_Earth)
2992  ssh(i,j) = ssh(i,j) + p_atm(i,j) * igr0
2993  enddo
2994  else
2995  do i=is,ie
2996  ssh(i,j) = ssh(i,j) + p_atm(i,j) * (us%Z_to_m / (gv%Rho0 * gv%g_Earth))
2997  enddo
2998  endif
2999  enddo
3000  endif ; endif
3001 

◆ extract_surface_state()

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

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_state_intransparent ocean surface state structure shared with the calling routine data in this structure is intent out.

Definition at line 3008 of file MOM.F90.

3008  type(MOM_control_struct), pointer :: CS !< Master MOM control structure
3009  type(surface), target, intent(inout) :: sfc_state_in !< transparent ocean surface state
3010  !! structure shared with the calling routine
3011  !! data in this structure is intent out.
3012 
3013  ! Local variables
3014  real :: hu, hv ! Thicknesses interpolated to velocity points [H ~> m or kg m-2]
3015  type(ocean_grid_type), pointer :: G => null() !< pointer to a structure containing
3016  !! metrics and related information
3017  type(ocean_grid_type), pointer :: G_in => null() !< Input grid metric
3018  type(verticalGrid_type), pointer :: GV => null() !< structure containing vertical grid info
3019  type(unit_scale_type), pointer :: US => null() !< structure containing various unit conversion factors
3020  type(surface), pointer :: sfc_state => null() ! surface state on the model grid
3021  real, dimension(:,:,:), pointer :: &
3022  h => null() !< h : layer thickness [H ~> m or kg m-2]
3023  real :: depth(SZI_(CS%G)) !< Distance from the surface in depth units [Z ~> m] or [H ~> m or kg m-2]
3024  real :: depth_ml !< Depth over which to average to determine mixed
3025  !! layer properties [Z ~> m] or [H ~> m or kg m-2]
3026  real :: dh !< Thickness of a layer within the mixed layer [Z ~> m] or [H ~> m or kg m-2]
3027  real :: mass !< Mass per unit area of a layer [R Z ~> kg m-2]
3028  real :: T_freeze !< freezing temperature [degC]
3029  real :: I_depth !< The inverse of depth [Z-1 ~> m-1] or [H-1 ~> m-1 or m2 kg-1]
3030  real :: missing_depth !< The portion of depth_ml that can not be found in a column [H ~> m or kg m-2]
3031  real :: H_rescale !< A conversion factor from thickness units to the units used in the
3032  !! calculation of properties of the uppermost ocean [nondim] or [Z H-1 ~> 1 or m3 kg-1]
3033  ! After the ANSWERS_2018 flag has been obsoleted, H_rescale will be 1.
3034  real :: delT(SZI_(CS%G)) !< Depth integral of T-T_freeze [Z degC ~> m degC]
3035  logical :: use_temperature !< If true, temp and saln used as state variables.
3036  integer :: i, j, k, is, ie, js, je, nz, numberOfErrors, ig, jg
3037  integer :: isd, ied, jsd, jed
3038  integer :: iscB, iecB, jscB, jecB, isdB, iedB, jsdB, jedB
3039  logical :: localError
3040  character(240) :: msg
3041  integer :: turns ! Number of quarter turns
3042 
3043  call calltree_enter("extract_surface_state(), MOM.F90")
3044  g => cs%G ; g_in => cs%G_in ; gv => cs%GV ; us => cs%US
3045  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
3046  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3047  iscb = g%iscB ; iecb = g%iecB; jscb = g%jscB ; jecb = g%jecB
3048  isdb = g%isdB ; iedb = g%iedB; jsdb = g%jsdB ; jedb = g%jedB
3049  h => cs%h
3050 
3051  use_temperature = associated(cs%tv%T)
3052 
3053  turns = 0
3054  if (cs%rotate_index) &
3055  turns = g%HI%turns
3056 
3057  if (.not.sfc_state_in%arrays_allocated) &
3058  ! Consider using a run-time flag to determine whether to do the vertical
3059  ! integrals, since the 3-d sums are not negligible in cost.
3060  call allocate_surface_state(sfc_state_in, g_in, use_temperature, &
3061  do_integrals=.true., omit_frazil=.not.associated(cs%tv%frazil))
3062 
3063  if (cs%rotate_index) then
3064  allocate(sfc_state)
3065  call allocate_surface_state(sfc_state, g, use_temperature, &
3066  do_integrals=.true., omit_frazil=.not.associated(cs%tv%frazil))
3067  else
3068  sfc_state => sfc_state_in
3069  endif
3070 
3071  sfc_state%T_is_conT = cs%tv%T_is_conT
3072  sfc_state%S_is_absS = cs%tv%S_is_absS
3073 
3074  do j=js,je ; do i=is,ie
3075  sfc_state%sea_lev(i,j) = us%m_to_Z*cs%ave_ssh_ibc(i,j)
3076  enddo ; enddo
3077 
3078  if (allocated(sfc_state%frazil) .and. associated(cs%tv%frazil)) then ; do j=js,je ; do i=is,ie
3079  sfc_state%frazil(i,j) = cs%tv%frazil(i,j)
3080  enddo ; enddo ; endif
3081 
3082  ! copy Hml into sfc_state, so that caps can access it
3083  if (associated(cs%Hml)) then
3084  do j=js,je ; do i=is,ie
3085  sfc_state%Hml(i,j) = cs%Hml(i,j)
3086  enddo ; enddo
3087  endif
3088 
3089  if (cs%Hmix < 0.0) then ! A bulk mixed layer is in use, so layer 1 has the properties
3090  if (use_temperature) then ; do j=js,je ; do i=is,ie
3091  sfc_state%SST(i,j) = cs%tv%T(i,j,1)
3092  sfc_state%SSS(i,j) = cs%tv%S(i,j,1)
3093  enddo ; enddo ; endif
3094  do j=js,je ; do i=is-1,ie
3095  sfc_state%u(i,j) = cs%u(i,j,1)
3096  enddo ; enddo
3097  do j=js-1,je ; do i=is,ie
3098  sfc_state%v(i,j) = cs%v(i,j,1)
3099  enddo ; enddo
3100 
3101  else ! (CS%Hmix >= 0.0)
3102  h_rescale = 1.0 ; if (cs%answers_2018) h_rescale = gv%H_to_Z
3103  depth_ml = cs%Hmix
3104  if (.not.cs%answers_2018) depth_ml = cs%Hmix*gv%Z_to_H
3105  ! Determine the mean tracer properties of the uppermost depth_ml fluid.
3106 
3107  !$OMP parallel do default(shared) private(depth,dh)
3108  do j=js,je
3109  do i=is,ie
3110  depth(i) = 0.0
3111  if (use_temperature) then
3112  sfc_state%SST(i,j) = 0.0 ; sfc_state%SSS(i,j) = 0.0
3113  else
3114  sfc_state%sfc_density(i,j) = 0.0
3115  endif
3116  enddo
3117 
3118  do k=1,nz ; do i=is,ie
3119  if (depth(i) + h(i,j,k)*h_rescale < depth_ml) then
3120  dh = h(i,j,k)*h_rescale
3121  elseif (depth(i) < depth_ml) then
3122  dh = depth_ml - depth(i)
3123  else
3124  dh = 0.0
3125  endif
3126  if (use_temperature) then
3127  sfc_state%SST(i,j) = sfc_state%SST(i,j) + dh * cs%tv%T(i,j,k)
3128  sfc_state%SSS(i,j) = sfc_state%SSS(i,j) + dh * cs%tv%S(i,j,k)
3129  else
3130  sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) + dh * gv%Rlay(k)
3131  endif
3132  depth(i) = depth(i) + dh
3133  enddo ; enddo
3134  ! Calculate the average properties of the mixed layer depth.
3135  do i=is,ie
3136  if (cs%answers_2018) then
3137  if (depth(i) < gv%H_subroundoff*h_rescale) &
3138  depth(i) = gv%H_subroundoff*h_rescale
3139  if (use_temperature) then
3140  sfc_state%SST(i,j) = sfc_state%SST(i,j) / depth(i)
3141  sfc_state%SSS(i,j) = sfc_state%SSS(i,j) / depth(i)
3142  else
3143  sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) / depth(i)
3144  endif
3145  else
3146  if (depth(i) < gv%H_subroundoff*h_rescale) then
3147  i_depth = 1.0 / (gv%H_subroundoff*h_rescale)
3148  missing_depth = gv%H_subroundoff*h_rescale - depth(i)
3149  if (use_temperature) then
3150  sfc_state%SST(i,j) = (sfc_state%SST(i,j) + missing_depth*cs%tv%T(i,j,1)) * i_depth
3151  sfc_state%SSS(i,j) = (sfc_state%SSS(i,j) + missing_depth*cs%tv%S(i,j,1)) * i_depth
3152  else
3153  sfc_state%sfc_density(i,j) = (sfc_state%sfc_density(i,j) + &
3154  missing_depth*gv%Rlay(1)) * i_depth
3155  endif
3156  else
3157  i_depth = 1.0 / depth(i)
3158  if (use_temperature) then
3159  sfc_state%SST(i,j) = sfc_state%SST(i,j) * i_depth
3160  sfc_state%SSS(i,j) = sfc_state%SSS(i,j) * i_depth
3161  else
3162  sfc_state%sfc_density(i,j) = sfc_state%sfc_density(i,j) * i_depth
3163  endif
3164  endif
3165  endif
3166  enddo
3167  enddo ! end of j loop
3168 
3169 ! Determine the mean velocities in the uppermost depth_ml fluid.
3170  ! NOTE: Velocity loops start on `[ij]s-1` in order to update halo values
3171  ! required by the speed diagnostic on the non-symmetric grid.
3172  ! This assumes that u and v halos have already been updated.
3173  if (cs%Hmix_UV>0.) then
3174  depth_ml = cs%Hmix_UV
3175  if (.not.cs%answers_2018) depth_ml = cs%Hmix_UV*gv%Z_to_H
3176  !$OMP parallel do default(shared) private(depth,dh,hv)
3177  do j=js-1,ie
3178  do i=is,ie
3179  depth(i) = 0.0
3180  sfc_state%v(i,j) = 0.0
3181  enddo
3182  do k=1,nz ; do i=is,ie
3183  hv = 0.5 * (h(i,j,k) + h(i,j+1,k)) * h_rescale
3184  if (depth(i) + hv < depth_ml) then
3185  dh = hv
3186  elseif (depth(i) < depth_ml) then
3187  dh = depth_ml - depth(i)
3188  else
3189  dh = 0.0
3190  endif
3191  sfc_state%v(i,j) = sfc_state%v(i,j) + dh * cs%v(i,j,k)
3192  depth(i) = depth(i) + dh
3193  enddo ; enddo
3194  ! Calculate the average properties of the mixed layer depth.
3195  do i=is,ie
3196  sfc_state%v(i,j) = sfc_state%v(i,j) / max(depth(i), gv%H_subroundoff*h_rescale)
3197  enddo
3198  enddo ! end of j loop
3199 
3200  !$OMP parallel do default(shared) private(depth,dh,hu)
3201  do j=js,je
3202  do i=is-1,ie
3203  depth(i) = 0.0
3204  sfc_state%u(i,j) = 0.0
3205  enddo
3206  do k=1,nz ; do i=is-1,ie
3207  hu = 0.5 * (h(i,j,k) + h(i+1,j,k)) * h_rescale
3208  if (depth(i) + hu < depth_ml) then
3209  dh = hu
3210  elseif (depth(i) < depth_ml) then
3211  dh = depth_ml - depth(i)
3212  else
3213  dh = 0.0
3214  endif
3215  sfc_state%u(i,j) = sfc_state%u(i,j) + dh * cs%u(i,j,k)
3216  depth(i) = depth(i) + dh
3217  enddo ; enddo
3218  ! Calculate the average properties of the mixed layer depth.
3219  do i=is-1,ie
3220  sfc_state%u(i,j) = sfc_state%u(i,j) / max(depth(i), gv%H_subroundoff*h_rescale)
3221  enddo
3222  enddo ! end of j loop
3223  else ! Hmix_UV<=0.
3224  do j=js,je ; do i=is-1,ie
3225  sfc_state%u(i,j) = cs%u(i,j,1)
3226  enddo ; enddo
3227  do j=js-1,je ; do i=is,ie
3228  sfc_state%v(i,j) = cs%v(i,j,1)
3229  enddo ; enddo
3230  endif
3231  endif ! (CS%Hmix >= 0.0)
3232 
3233 
3234  if (allocated(sfc_state%melt_potential)) then
3235  !$OMP parallel do default(shared) private(depth_ml, dh, T_freeze, depth, delT)
3236  do j=js,je
3237  do i=is,ie
3238  depth(i) = 0.0
3239  delt(i) = 0.0
3240  enddo
3241 
3242  do k=1,nz ; do i=is,ie
3243  depth_ml = min(cs%HFrz, cs%visc%MLD(i,j))
3244  if (depth(i) + h(i,j,k)*gv%H_to_Z < depth_ml) then
3245  dh = h(i,j,k)*gv%H_to_Z
3246  elseif (depth(i) < depth_ml) then
3247  dh = depth_ml - depth(i)
3248  else
3249  dh = 0.0
3250  endif
3251 
3252  ! p=0 OK, HFrz ~ 10 to 20m
3253  call calculate_tfreeze(cs%tv%S(i,j,k), 0.0, t_freeze, cs%tv%eqn_of_state)
3254  depth(i) = depth(i) + dh
3255  delt(i) = delt(i) + dh * (cs%tv%T(i,j,k) - t_freeze)
3256  enddo ; enddo
3257 
3258  do i=is,ie
3259  ! set melt_potential to zero to avoid passing previous values
3260  sfc_state%melt_potential(i,j) = 0.0
3261 
3262  if (g%mask2dT(i,j)>0.) then
3263  ! instantaneous melt_potential [Q R Z ~> J m-2]
3264  sfc_state%melt_potential(i,j) = cs%tv%C_p * gv%Rho0 * delt(i)
3265  endif
3266  enddo
3267  enddo ! end of j loop
3268  endif ! melt_potential
3269 
3270  if (allocated(sfc_state%salt_deficit) .and. associated(cs%tv%salt_deficit)) then
3271  !$OMP parallel do default(shared)
3272  do j=js,je ; do i=is,ie
3273  ! Convert from gSalt to kgSalt
3274  sfc_state%salt_deficit(i,j) = 0.001 * cs%tv%salt_deficit(i,j)
3275  enddo ; enddo
3276  endif
3277  if (allocated(sfc_state%TempxPmE) .and. associated(cs%tv%TempxPmE)) then
3278  !$OMP parallel do default(shared)
3279  do j=js,je ; do i=is,ie
3280  sfc_state%TempxPmE(i,j) = cs%tv%TempxPmE(i,j)
3281  enddo ; enddo
3282  endif
3283  if (allocated(sfc_state%internal_heat) .and. associated(cs%tv%internal_heat)) then
3284  !$OMP parallel do default(shared)
3285  do j=js,je ; do i=is,ie
3286  sfc_state%internal_heat(i,j) = cs%tv%internal_heat(i,j)
3287  enddo ; enddo
3288  endif
3289  if (allocated(sfc_state%taux_shelf) .and. associated(cs%visc%taux_shelf)) then
3290  !$OMP parallel do default(shared)
3291  do j=js,je ; do i=is-1,ie
3292  sfc_state%taux_shelf(i,j) = cs%visc%taux_shelf(i,j)
3293  enddo ; enddo
3294  endif
3295  if (allocated(sfc_state%tauy_shelf) .and. associated(cs%visc%tauy_shelf)) then
3296  !$OMP parallel do default(shared)
3297  do j=js-1,je ; do i=is,ie
3298  sfc_state%tauy_shelf(i,j) = cs%visc%tauy_shelf(i,j)
3299  enddo ; enddo
3300  endif
3301 
3302  if (allocated(sfc_state%ocean_mass) .and. allocated(sfc_state%ocean_heat) .and. &
3303  allocated(sfc_state%ocean_salt)) then
3304  !$OMP parallel do default(shared)
3305  do j=js,je ; do i=is,ie
3306  sfc_state%ocean_mass(i,j) = 0.0
3307  sfc_state%ocean_heat(i,j) = 0.0 ; sfc_state%ocean_salt(i,j) = 0.0
3308  enddo ; enddo
3309  !$OMP parallel do default(shared) private(mass)
3310  do j=js,je ; do k=1,nz; do i=is,ie
3311  mass = gv%H_to_RZ*h(i,j,k)
3312  sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + mass
3313  sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass * cs%tv%T(i,j,k)
3314  sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + mass * (1.0e-3*cs%tv%S(i,j,k))
3315  enddo ; enddo ; enddo
3316  else
3317  if (allocated(sfc_state%ocean_mass)) then
3318  !$OMP parallel do default(shared)
3319  do j=js,je ; do i=is,ie ; sfc_state%ocean_mass(i,j) = 0.0 ; enddo ; enddo
3320  !$OMP parallel do default(shared)
3321  do j=js,je ; do k=1,nz ; do i=is,ie
3322  sfc_state%ocean_mass(i,j) = sfc_state%ocean_mass(i,j) + gv%H_to_RZ*h(i,j,k)
3323  enddo ; enddo ; enddo
3324  endif
3325  if (allocated(sfc_state%ocean_heat)) then
3326  !$OMP parallel do default(shared)
3327  do j=js,je ; do i=is,ie ; sfc_state%ocean_heat(i,j) = 0.0 ; enddo ; enddo
3328  !$OMP parallel do default(shared) private(mass)
3329  do j=js,je ; do k=1,nz ; do i=is,ie
3330  mass = gv%H_to_RZ*h(i,j,k)
3331  sfc_state%ocean_heat(i,j) = sfc_state%ocean_heat(i,j) + mass*cs%tv%T(i,j,k)
3332  enddo ; enddo ; enddo
3333  endif
3334  if (allocated(sfc_state%ocean_salt)) then
3335  !$OMP parallel do default(shared)
3336  do j=js,je ; do i=is,ie ; sfc_state%ocean_salt(i,j) = 0.0 ; enddo ; enddo
3337  !$OMP parallel do default(shared) private(mass)
3338  do j=js,je ; do k=1,nz ; do i=is,ie
3339  mass = gv%H_to_RZ*h(i,j,k)
3340  sfc_state%ocean_salt(i,j) = sfc_state%ocean_salt(i,j) + mass * (1.0e-3*cs%tv%S(i,j,k))
3341  enddo ; enddo ; enddo
3342  endif
3343  endif
3344 
3345  if (associated(cs%tracer_flow_CSp)) then
3346  call call_tracer_surface_state(sfc_state, h, g, cs%tracer_flow_CSp)
3347  endif
3348 
3349  if (cs%check_bad_sfc_vals) then
3350  numberoferrors=0 ! count number of errors
3351  do j=js,je; do i=is,ie
3352  if (g%mask2dT(i,j)>0.) then
3353  localerror = sfc_state%sea_lev(i,j) <= -g%bathyT(i,j) &
3354  .or. sfc_state%sea_lev(i,j) >= cs%bad_val_ssh_max &
3355  .or. sfc_state%sea_lev(i,j) <= -cs%bad_val_ssh_max &
3356  .or. sfc_state%sea_lev(i,j) + g%bathyT(i,j) < cs%bad_val_col_thick
3357  if (use_temperature) localerror = localerror &
3358  .or. sfc_state%SSS(i,j)<0. &
3359  .or. sfc_state%SSS(i,j)>=cs%bad_val_sss_max &
3360  .or. sfc_state%SST(i,j)< cs%bad_val_sst_min &
3361  .or. sfc_state%SST(i,j)>=cs%bad_val_sst_max
3362  if (localerror) then
3363  numberoferrors=numberoferrors+1
3364  if (numberoferrors<9) then ! Only report details for the first few errors
3365  ig = i + g%HI%idg_offset ! Global i-index
3366  jg = j + g%HI%jdg_offset ! Global j-index
3367  if (use_temperature) then
3368  write(msg(1:240),'(2(a,i4,x),4(a,f8.3,x),8(a,es11.4,x))') &
3369  'Extreme surface sfc_state detected: i=',ig,'j=',jg, &
3370  'lon=',g%geoLonT(i,j), 'lat=',g%geoLatT(i,j), &
3371  'x=',g%gridLonT(ig), 'y=',g%gridLatT(jg), &
3372  'D=',cs%US%Z_to_m*g%bathyT(i,j), 'SSH=',cs%US%Z_to_m*sfc_state%sea_lev(i,j), &
3373  'SST=',sfc_state%SST(i,j), 'SSS=',sfc_state%SSS(i,j), &
3374  'U-=',us%L_T_to_m_s*sfc_state%u(i-1,j), 'U+=',us%L_T_to_m_s*sfc_state%u(i,j), &
3375  'V-=',us%L_T_to_m_s*sfc_state%v(i,j-1), 'V+=',us%L_T_to_m_s*sfc_state%v(i,j)
3376  else
3377  write(msg(1:240),'(2(a,i4,x),4(a,f8.3,x),6(a,es11.4))') &
3378  'Extreme surface sfc_state detected: i=',ig,'j=',jg, &
3379  'lon=',g%geoLonT(i,j), 'lat=',g%geoLatT(i,j), &
3380  'x=',g%gridLonT(i), 'y=',g%gridLatT(j), &
3381  'D=',cs%US%Z_to_m*g%bathyT(i,j), 'SSH=',cs%US%Z_to_m*sfc_state%sea_lev(i,j), &
3382  'U-=',us%L_T_to_m_s*sfc_state%u(i-1,j), 'U+=',us%L_T_to_m_s*sfc_state%u(i,j), &
3383  'V-=',us%L_T_to_m_s*sfc_state%v(i,j-1), 'V+=',us%L_T_to_m_s*sfc_state%v(i,j)
3384  endif
3385  call mom_error(warning, trim(msg), all_print=.true.)
3386  elseif (numberoferrors==9) then ! Indicate once that there are more errors
3387  call mom_error(warning, 'There were more unreported extreme events!', all_print=.true.)
3388  endif ! numberOfErrors
3389  endif ! localError
3390  endif ! mask2dT
3391  enddo ; enddo
3392  call sum_across_pes(numberoferrors)
3393  if (numberoferrors>0) then
3394  write(msg(1:240),'(3(a,i9,x))') 'There were a total of ',numberoferrors, &
3395  'locations detected with extreme surface values!'
3396  call mom_error(fatal, trim(msg))
3397  endif
3398  endif
3399 
3400  if (cs%debug) call mom_surface_chksum("Post extract_sfc", sfc_state, g, us, haloshift=0)
3401 
3402  ! Rotate sfc_state back onto the input grid, sfc_state_in
3403  if (cs%rotate_index) then
3404  call rotate_surface_state(sfc_state, g, sfc_state_in, g_in, -turns)
3405  call deallocate_surface_state(sfc_state)
3406  endif
3407 
3408  call calltree_leave("extract_surface_sfc_state()")

◆ 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 2774 of file MOM.F90.

2774  type(time_type), intent(in) :: Time !< model time, used in this routine
2775  type(directories), intent(in) :: dirs !< structure with directory paths
2776  type(MOM_control_struct), pointer :: CS !< pointer to MOM control structure
2777  type(MOM_restart_CS), pointer :: restart_CSp !< pointer to the restart control
2778  !! structure that will be used for MOM.
2779  ! Local variables
2780  type(ocean_grid_type), pointer :: G => null() ! pointer to a structure containing
2781  ! metrics and related information
2782  type(verticalGrid_type), pointer :: GV => null() ! Pointer to the vertical grid structure
2783  type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing
2784  ! various unit conversion factors
2785  type(MOM_restart_CS), pointer :: restart_CSp_tmp => null()
2786  real, allocatable :: z_interface(:,:,:) ! Interface heights [m]
2787  type(vardesc) :: vd
2788 
2789  call cpu_clock_begin(id_clock_init)
2790  call calltree_enter("finish_MOM_initialization()")
2791 
2792  ! Pointers for convenience
2793  g => cs%G ; gv => cs%GV ; us => cs%US
2794 
2795  !### Move to initialize_MOM?
2796  call fix_restart_scaling(gv)
2797  call fix_restart_unit_scaling(us)
2798 
2799  ! Write initial conditions
2800  if (cs%write_IC) then
2801  allocate(restart_csp_tmp)
2802  restart_csp_tmp = restart_csp
2803  allocate(z_interface(szi_(g),szj_(g),szk_(g)+1))
2804  call find_eta(cs%h, cs%tv, g, gv, us, z_interface, eta_to_m=1.0)
2805  call register_restart_field(z_interface, "eta", .true., restart_csp_tmp, &
2806  "Interface heights", "meter", z_grid='i')
2807 
2808  call save_restart(dirs%output_directory, time, cs%G_in, &
2809  restart_csp_tmp, filename=cs%IC_file, gv=gv)
2810  deallocate(z_interface)
2811  deallocate(restart_csp_tmp)
2812  endif
2813 
2814  call write_energy(cs%u, cs%v, cs%h, cs%tv, time, 0, g, gv, us, &
2815  cs%sum_output_CSp, cs%tracer_flow_CSp)
2816 
2817  call calltree_leave("finish_MOM_initialization()")
2818  call cpu_clock_end(id_clock_init)
2819 

◆ 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,
real, intent(out), optional  C_p_scaled,
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 [J kg degC-1]
[out]c_p_scaledThe heat capacity in scaled units [Q degC-1 ~> J kg degC-1]
[out]use_tempTrue if temperature is a state variable

Definition at line 3450 of file MOM.F90.

3450  type(MOM_control_struct), pointer :: CS !< MOM control structure
3451  type(ocean_grid_type), optional, pointer :: G !< structure containing metrics and grid info
3452  type(verticalGrid_type), optional, pointer :: GV !< structure containing vertical grid info
3453  type(unit_scale_type), optional, pointer :: US !< A dimensional unit scaling type
3454  real, optional, intent(out) :: C_p !< The heat capacity [J kg degC-1]
3455  real, optional, intent(out) :: C_p_scaled !< The heat capacity in scaled
3456  !! units [Q degC-1 ~> J kg degC-1]
3457  logical, optional, intent(out) :: use_temp !< True if temperature is a state variable
3458 
3459  if (present(g)) g => cs%G_in
3460  if (present(gv)) gv => cs%GV
3461  if (present(us)) us => cs%US
3462  if (present(c_p)) c_p = cs%US%Q_to_J_kg * cs%tv%C_p
3463  if (present(c_p_scaled)) c_p_scaled = cs%tv%C_p
3464  if (present(use_temp)) use_temp = associated(cs%tv%T)

◆ 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 3469 of file MOM.F90.

3469  type(MOM_control_struct), pointer :: CS !< MOM control structure
3470  real, optional, intent(out) :: heat !< The globally integrated integrated ocean heat [J].
3471  real, optional, intent(out) :: salt !< The globally integrated integrated ocean salt [kg].
3472  real, optional, intent(out) :: mass !< The globally integrated integrated ocean mass [kg].
3473  logical, optional, intent(in) :: on_PE_only !< If present and true, only sum on the local PE.
3474 
3475  if (present(mass)) &
3476  mass = global_mass_integral(cs%h, cs%G, cs%GV, on_pe_only=on_pe_only)
3477  if (present(heat)) &
3478  heat = cs%US%Q_to_J_kg*cs%tv%C_p * global_mass_integral(cs%h, cs%G, cs%GV, cs%tv%T, on_pe_only=on_pe_only)
3479  if (present(salt)) &
3480  salt = 1.0e-3 * global_mass_integral(cs%h, cs%G, cs%GV, cs%tv%S, on_pe_only=on_pe_only)
3481 

◆ 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 1604 of file MOM.F90.

1604  type(time_type), target, intent(inout) :: Time !< model time, set in this routine
1605  type(time_type), intent(in) :: Time_init !< The start time for the coupled model's calendar
1606  type(param_file_type), intent(out) :: param_file !< structure indicating parameter file to parse
1607  type(directories), intent(out) :: dirs !< structure with directory paths
1608  type(MOM_control_struct), pointer :: CS !< pointer set in this routine to MOM control structure
1609  type(MOM_restart_CS), pointer :: restart_CSp !< pointer set in this routine to the
1610  !! restart control structure that will
1611  !! be used for MOM.
1612  type(time_type), optional, intent(in) :: Time_in !< time passed to MOM_initialize_state when
1613  !! model is not being started from a restart file
1614  logical, optional, intent(out) :: offline_tracer_mode !< True is returned if tracers are being run offline
1615  character(len=*),optional, intent(in) :: input_restart_file !< If present, name of restart file to read
1616  type(diag_ctrl), optional, pointer :: diag_ptr !< A pointer set in this routine to the diagnostic
1617  !! regulatory structure
1618  type(tracer_flow_control_CS), &
1619  optional, pointer :: tracer_flow_CSp !< A pointer set in this routine to
1620  !! the tracer flow control structure.
1621  logical, optional, intent(in) :: count_calls !< If true, nstep_tot counts the number of
1622  !! calls to step_MOM instead of the number of
1623  !! dynamics timesteps.
1624  ! local variables
1625  type(ocean_grid_type), pointer :: G => null() ! A pointer to the metric grid use for the run
1626  type(ocean_grid_type), pointer :: G_in => null() ! Pointer to the input grid
1627  type(hor_index_type), pointer :: HI => null() ! A hor_index_type for array extents
1628  type(hor_index_type), target :: HI_in ! HI on the input grid
1629  type(verticalGrid_type), pointer :: GV => null()
1630  type(dyn_horgrid_type), pointer :: dG => null()
1631  type(dyn_horgrid_type), pointer :: dG_in => null()
1632  type(diag_ctrl), pointer :: diag => null()
1633  type(unit_scale_type), pointer :: US => null()
1634  character(len=4), parameter :: vers_num = 'v2.0'
1635  integer :: turns ! Number of grid quarter-turns
1636 
1637  ! Initial state on the input index map
1638  real, allocatable, dimension(:,:,:) :: u_in, v_in, h_in
1639  real, allocatable, dimension(:,:,:), target :: T_in, S_in
1640  type(ocean_OBC_type), pointer :: OBC_in => null()
1641  type(sponge_CS), pointer :: sponge_in_CSp => null()
1642  type(ALE_sponge_CS), pointer :: ALE_sponge_in_CSp => null()
1643 
1644  ! This include declares and sets the variable "version".
1645 # include "version_variable.h"
1646 
1647  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
1648  integer :: IsdB, IedB, JsdB, JedB
1649  real :: dtbt ! The barotropic timestep [s]
1650 
1651  real, allocatable, dimension(:,:) :: eta ! free surface height or column mass [H ~> m or kg m-2]
1652  real, allocatable, dimension(:,:) :: area_shelf_h ! area occupied by ice shelf [L2 ~> m2]
1653  real, dimension(:,:), allocatable, target :: frac_shelf_h ! fraction of total area occupied by ice shelf [nondim]
1654  real, dimension(:,:), pointer :: shelf_area => null()
1655  type(MOM_restart_CS), pointer :: restart_CSp_tmp => null()
1656  type(group_pass_type) :: tmp_pass_uv_T_S_h, pass_uv_T_S_h
1657 
1658  real :: default_val ! default value for a parameter
1659  logical :: write_geom_files ! If true, write out the grid geometry files.
1660  logical :: ensemble_ocean ! If true, perform an ensemble gather at the end of step_MOM
1661  logical :: new_sim
1662  logical :: use_geothermal ! If true, apply geothermal heating.
1663  logical :: use_EOS ! If true, density calculated from T & S using an equation of state.
1664  logical :: symmetric ! If true, use symmetric memory allocation.
1665  logical :: save_IC ! If true, save the initial conditions.
1666  logical :: do_unit_tests ! If true, call unit tests.
1667  logical :: test_grid_copy = .false.
1668 
1669  logical :: bulkmixedlayer ! If true, a refined bulk mixed layer scheme is used
1670  ! with nkml sublayers and nkbl buffer layer.
1671  logical :: use_temperature ! If true, temp and saln used as state variables.
1672  logical :: use_frazil ! If true, liquid seawater freezes if temp below freezing,
1673  ! with accumulated heat deficit returned to surface ocean.
1674  logical :: bound_salinity ! If true, salt is added to keep salinity above
1675  ! a minimum value, and the deficit is reported.
1676  logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags.
1677  logical :: use_conT_absS ! If true, the prognostics T & S are conservative temperature
1678  ! and absolute salinity. Care should be taken to convert them
1679  ! to potential temperature and practical salinity before
1680  ! exchanging them with the coupler and/or reporting T&S diagnostics.
1681  logical :: advect_TS ! If false, then no horizontal advection of temperature
1682  ! and salnity is performed
1683  logical :: use_ice_shelf ! Needed for ALE
1684  logical :: global_indexing ! If true use global horizontal index values instead
1685  ! of having the data domain on each processor start at 1.
1686  logical :: bathy_at_vel ! If true, also define bathymetric fields at the
1687  ! the velocity points.
1688  logical :: calc_dtbt ! Indicates whether the dynamically adjusted barotropic
1689  ! time step needs to be updated before it is used.
1690  logical :: debug_truncations ! If true, turn on diagnostics useful for debugging truncations.
1691  integer :: first_direction ! An integer that indicates which direction is to be
1692  ! updated first in directionally split parts of the
1693  ! calculation. This can be altered during the course
1694  ! of the run via calls to set_first_direction.
1695  integer :: nkml, nkbl, verbosity, write_geom
1696  integer :: dynamics_stencil ! The computational stencil for the calculations
1697  ! in the dynamic core.
1698  real :: conv2watt, conv2salt
1699  real :: RL2_T2_rescale, Z_rescale, QRZ_rescale ! Unit conversion factors
1700  character(len=48) :: flux_units, S_flux_units
1701 
1702  type(vardesc) :: vd_T, vd_S ! Structures describing temperature and salinity variables.
1703  type(time_type) :: Start_time
1704  type(ocean_internal_state) :: MOM_internal_state
1705  character(len=200) :: area_varname, ice_shelf_file, inputdir, filename
1706 
1707  if (associated(cs)) then
1708  call mom_error(warning, "initialize_MOM called with an associated "// &
1709  "control structure.")
1710  return
1711  endif
1712  allocate(cs)
1713 
1714  cs%Time => time
1715 
1716  id_clock_init = cpu_clock_id('Ocean Initialization', grain=clock_subcomponent)
1717  call cpu_clock_begin(id_clock_init)
1718 
1719  start_time = time ; if (present(time_in)) start_time = time_in
1720 
1721  ! Read paths and filenames from namelist and store in "dirs".
1722  ! Also open the parsed input parameter file(s) and setup param_file.
1723  call get_mom_input(param_file, dirs, default_input_filename=input_restart_file)
1724 
1725  verbosity = 2 ; call read_param(param_file, "VERBOSITY", verbosity)
1726  call mom_set_verbosity(verbosity)
1727  call calltree_enter("initialize_MOM(), MOM.F90")
1728 
1729  call find_obsolete_params(param_file)
1730 
1731  ! Determining the internal unit scaling factors for this run.
1732  call unit_scaling_init(param_file, cs%US)
1733  us => cs%US
1734 
1735  ! Read relevant parameters and write them to the model log.
1736  call log_version(param_file, "MOM", version, "", log_to_all=.true., layout=.true., debugging=.true.)
1737  call get_param(param_file, "MOM", "VERBOSITY", verbosity, &
1738  "Integer controlling level of messaging\n" // &
1739  "\t0 = Only FATAL messages\n" // &
1740  "\t2 = Only FATAL, WARNING, NOTE [default]\n" // &
1741  "\t9 = All)", default=2, debuggingparam=.true.)
1742  call get_param(param_file, "MOM", "DO_UNIT_TESTS", do_unit_tests, &
1743  "If True, exercises unit tests at model start up.", &
1744  default=.false., debuggingparam=.true.)
1745  if (do_unit_tests) then
1746  id_clock_unit_tests = cpu_clock_id('(Ocean unit tests)', grain=clock_module)
1747  call cpu_clock_begin(id_clock_unit_tests)
1748  call unit_tests(verbosity)
1749  call cpu_clock_end(id_clock_unit_tests)
1750  endif
1751 
1752  call get_param(param_file, "MOM", "SPLIT", cs%split, &
1753  "Use the split time stepping if true.", default=.true.)
1754  if (cs%split) then
1755  cs%use_RK2 = .false.
1756  else
1757  call get_param(param_file, "MOM", "USE_RK2", cs%use_RK2, &
1758  "If true, use RK2 instead of RK3 in the unsplit time stepping.", &
1759  default=.false.)
1760  endif
1761 
1762  call get_param(param_file, "MOM", "CALC_RHO_FOR_SEA_LEVEL", cs%calc_rho_for_sea_lev, &
1763  "If true, the in-situ density is used to calculate the "//&
1764  "effective sea level that is returned to the coupler. If false, "//&
1765  "the Boussinesq parameter RHO_0 is used.", default=.false.)
1766  call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", use_temperature, &
1767  "If true, Temperature and salinity are used as state "//&
1768  "variables.", default=.true.)
1769  call get_param(param_file, "MOM", "USE_EOS", use_eos, &
1770  "If true, density is calculated from temperature and "//&
1771  "salinity with an equation of state. If USE_EOS is "//&
1772  "true, ENABLE_THERMODYNAMICS must be true as well.", &
1773  default=use_temperature)
1774  call get_param(param_file, "MOM", "DIABATIC_FIRST", cs%diabatic_first, &
1775  "If true, apply diabatic and thermodynamic processes, "//&
1776  "including buoyancy forcing and mass gain or loss, "//&
1777  "before stepping the dynamics forward.", default=.false.)
1778  call get_param(param_file, "MOM", "USE_CONTEMP_ABSSAL", use_cont_abss, &
1779  "If true, the prognostics T&S are the conservative temperature "//&
1780  "and absolute salinity. Care should be taken to convert them "//&
1781  "to potential temperature and practical salinity before "//&
1782  "exchanging them with the coupler and/or reporting T&S diagnostics.", &
1783  default=.false.)
1784  cs%tv%T_is_conT = use_cont_abss ; cs%tv%S_is_absS = use_cont_abss
1785  call get_param(param_file, "MOM", "ADIABATIC", cs%adiabatic, &
1786  "There are no diapycnal mass fluxes if ADIABATIC is "//&
1787  "true. This assumes that KD = KDML = 0.0 and that "//&
1788  "there is no buoyancy forcing, but makes the model "//&
1789  "faster by eliminating subroutine calls.", default=.false.)
1790  call get_param(param_file, "MOM", "DO_DYNAMICS", cs%do_dynamics, &
1791  "If False, skips the dynamics calls that update u & v, as well as "//&
1792  "the gravity wave adjustment to h. This may be a fragile feature, "//&
1793  "but can be useful during development", default=.true.)
1794  call get_param(param_file, "MOM", "ADVECT_TS", advect_ts, &
1795  "If True, advect temperature and salinity horizontally "//&
1796  "If False, T/S are registered for advection. "//&
1797  "This is intended only to be used in offline tracer mode "//&
1798  "and is by default false in that case.", &
1799  do_not_log = .true., default=.true. )
1800  if (present(offline_tracer_mode)) then ! Only read this parameter in enabled modes
1801  call get_param(param_file, "MOM", "OFFLINE_TRACER_MODE", cs%offline_tracer_mode, &
1802  "If true, barotropic and baroclinic dynamics, thermodynamics "//&
1803  "are all bypassed with all the fields necessary to integrate "//&
1804  "the tracer advection and diffusion equation are read in from "//&
1805  "files stored from a previous integration of the prognostic model. "//&
1806  "NOTE: This option only used in the ocean_solo_driver.", default=.false.)
1807  if (cs%offline_tracer_mode) then
1808  call get_param(param_file, "MOM", "ADVECT_TS", advect_ts, &
1809  "If True, advect temperature and salinity horizontally "//&
1810  "If False, T/S are registered for advection. "//&
1811  "This is intended only to be used in offline tracer mode."//&
1812  "and is by default false in that case", &
1813  default=.false. )
1814  endif
1815  endif
1816  call get_param(param_file, "MOM", "USE_REGRIDDING", cs%use_ALE_algorithm, &
1817  "If True, use the ALE algorithm (regridding/remapping). "//&
1818  "If False, use the layered isopycnal algorithm.", default=.false. )
1819  call get_param(param_file, "MOM", "BULKMIXEDLAYER", bulkmixedlayer, &
1820  "If true, use a Kraus-Turner-like bulk mixed layer "//&
1821  "with transitional buffer layers. Layers 1 through "//&
1822  "NKML+NKBL have variable densities. There must be at "//&
1823  "least NKML+NKBL+1 layers if BULKMIXEDLAYER is true. "//&
1824  "BULKMIXEDLAYER can not be used with USE_REGRIDDING. "//&
1825  "The default is influenced by ENABLE_THERMODYNAMICS.", &
1826  default=use_temperature .and. .not.cs%use_ALE_algorithm)
1827  call get_param(param_file, "MOM", "THICKNESSDIFFUSE", cs%thickness_diffuse, &
1828  "If true, interface heights are diffused with a "//&
1829  "coefficient of KHTH.", default=.false.)
1830  call get_param(param_file, "MOM", "THICKNESSDIFFUSE_FIRST", &
1831  cs%thickness_diffuse_first, &
1832  "If true, do thickness diffusion before dynamics. "//&
1833  "This is only used if THICKNESSDIFFUSE is true.", &
1834  default=.false.)
1835  if (.not.cs%thickness_diffuse) cs%thickness_diffuse_first = .false.
1836  call get_param(param_file, "MOM", "BATHYMETRY_AT_VEL", bathy_at_vel, &
1837  "If true, there are separate values for the basin depths "//&
1838  "at velocity points. Otherwise the effects of topography "//&
1839  "are entirely determined from thickness points.", &
1840  default=.false.)
1841  call get_param(param_file, "MOM", "USE_WAVES", cs%UseWaves, default=.false., &
1842  do_not_log=.true.)
1843 
1844  call get_param(param_file, "MOM", "DEBUG", cs%debug, &
1845  "If true, write out verbose debugging data.", &
1846  default=.false., debuggingparam=.true.)
1847  call get_param(param_file, "MOM", "DEBUG_TRUNCATIONS", debug_truncations, &
1848  "If true, calculate all diagnostics that are useful for "//&
1849  "debugging truncations.", default=.false., debuggingparam=.true.)
1850 
1851  call get_param(param_file, "MOM", "DT", cs%dt, &
1852  "The (baroclinic) dynamics time step. The time-step that "//&
1853  "is actually used will be an integer fraction of the "//&
1854  "forcing time-step (DT_FORCING in ocean-only mode or the "//&
1855  "coupling timestep in coupled mode.)", units="s", scale=us%s_to_T, &
1856  fail_if_missing=.true.)
1857  call get_param(param_file, "MOM", "DT_THERM", cs%dt_therm, &
1858  "The thermodynamic and tracer advection time step. "//&
1859  "Ideally DT_THERM should be an integer multiple of DT "//&
1860  "and less than the forcing or coupling time-step, unless "//&
1861  "THERMO_SPANS_COUPLING is true, in which case DT_THERM "//&
1862  "can be an integer multiple of the coupling timestep. By "//&
1863  "default DT_THERM is set to DT.", &
1864  units="s", scale=us%s_to_T, default=us%T_to_s*cs%dt)
1865  call get_param(param_file, "MOM", "THERMO_SPANS_COUPLING", cs%thermo_spans_coupling, &
1866  "If true, the MOM will take thermodynamic and tracer "//&
1867  "timesteps that can be longer than the coupling timestep. "//&
1868  "The actual thermodynamic timestep that is used in this "//&
1869  "case is the largest integer multiple of the coupling "//&
1870  "timestep that is less than or equal to DT_THERM.", default=.false.)
1871 
1872  if (bulkmixedlayer) then
1873  cs%Hmix = -1.0 ; cs%Hmix_UV = -1.0
1874  else
1875  call get_param(param_file, "MOM", "HMIX_SFC_PROP", cs%Hmix, &
1876  "If BULKMIXEDLAYER is false, HMIX_SFC_PROP is the depth "//&
1877  "over which to average to find surface properties like "//&
1878  "SST and SSS or density (but not surface velocities).", &
1879  units="m", default=1.0, scale=us%m_to_Z)
1880  call get_param(param_file, "MOM", "HMIX_UV_SFC_PROP", cs%Hmix_UV, &
1881  "If BULKMIXEDLAYER is false, HMIX_UV_SFC_PROP is the depth "//&
1882  "over which to average to find surface flow properties, "//&
1883  "SSU, SSV. A non-positive value indicates no averaging.", &
1884  units="m", default=0.0, scale=us%m_to_Z)
1885  endif
1886  call get_param(param_file, "MOM", "HFREEZE", cs%HFrz, &
1887  "If HFREEZE > 0, melt potential will be computed. The actual depth "//&
1888  "over which melt potential is computed will be min(HFREEZE, OBLD), "//&
1889  "where OBLD is the boundary layer depth. If HFREEZE <= 0 (default), "//&
1890  "melt potential will not be computed.", units="m", default=-1.0, scale=us%m_to_Z)
1891  call get_param(param_file, "MOM", "INTERPOLATE_P_SURF", cs%interp_p_surf, &
1892  "If true, linearly interpolate the surface pressure "//&
1893  "over the coupling time step, using the specified value "//&
1894  "at the end of the step.", default=.false.)
1895 
1896  if (cs%split) then
1897  call get_param(param_file, "MOM", "DTBT", dtbt, default=-0.98)
1898  default_val = us%T_to_s*cs%dt_therm ; if (dtbt > 0.0) default_val = -1.0
1899  cs%dtbt_reset_period = -1.0
1900  call get_param(param_file, "MOM", "DTBT_RESET_PERIOD", cs%dtbt_reset_period, &
1901  "The period between recalculations of DTBT (if DTBT <= 0). "//&
1902  "If DTBT_RESET_PERIOD is negative, DTBT is set based "//&
1903  "only on information available at initialization. If 0, "//&
1904  "DTBT will be set every dynamics time step. The default "//&
1905  "is set by DT_THERM. This is only used if SPLIT is true.", &
1906  units="s", default=default_val, do_not_read=(dtbt > 0.0))
1907  endif
1908 
1909  ! This is here in case these values are used inappropriately.
1910  use_frazil = .false. ; bound_salinity = .false.
1911  cs%tv%P_Ref = 2.0e7*us%kg_m3_to_R*us%m_s_to_L_T**2
1912  if (use_temperature) then
1913  call get_param(param_file, "MOM", "FRAZIL", use_frazil, &
1914  "If true, water freezes if it gets too cold, and the "//&
1915  "accumulated heat deficit is returned in the "//&
1916  "surface state. FRAZIL is only used if "//&
1917  "ENABLE_THERMODYNAMICS is true.", default=.false.)
1918  call get_param(param_file, "MOM", "DO_GEOTHERMAL", use_geothermal, &
1919  "If true, apply geothermal heating.", default=.false.)
1920  call get_param(param_file, "MOM", "BOUND_SALINITY", bound_salinity, &
1921  "If true, limit salinity to being positive. (The sea-ice "//&
1922  "model may ask for more salt than is available and "//&
1923  "drive the salinity negative otherwise.)", default=.false.)
1924  call get_param(param_file, "MOM", "MIN_SALINITY", cs%tv%min_salinity, &
1925  "The minimum value of salinity when BOUND_SALINITY=True.", &
1926  units="PPT", default=0.0, do_not_log=.not.bound_salinity)
1927  call get_param(param_file, "MOM", "C_P", cs%tv%C_p, &
1928  "The heat capacity of sea water, approximated as a "//&
1929  "constant. This is only used if ENABLE_THERMODYNAMICS is "//&
1930  "true. The default value is from the TEOS-10 definition "//&
1931  "of conservative temperature.", units="J kg-1 K-1", &
1932  default=3991.86795711963, scale=us%J_kg_to_Q)
1933  call get_param(param_file, "MOM", "USE_PSURF_IN_EOS", cs%use_p_surf_in_EOS, &
1934  "If true, always include the surface pressure contributions "//&
1935  "in equation of state calculations.", default=.true.)
1936  endif
1937  if (use_eos) call get_param(param_file, "MOM", "P_REF", cs%tv%P_Ref, &
1938  "The pressure that is used for calculating the coordinate "//&
1939  "density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
1940  "This is only used if USE_EOS and ENABLE_THERMODYNAMICS are true.", &
1941  units="Pa", default=2.0e7, scale=us%kg_m3_to_R*us%m_s_to_L_T**2)
1942 
1943  if (bulkmixedlayer) then
1944  call get_param(param_file, "MOM", "NKML", nkml, &
1945  "The number of sublayers within the mixed layer if "//&
1946  "BULKMIXEDLAYER is true.", units="nondim", default=2)
1947  call get_param(param_file, "MOM", "NKBL", nkbl, &
1948  "The number of layers that are used as variable density "//&
1949  "buffer layers if BULKMIXEDLAYER is true.", units="nondim", &
1950  default=2)
1951  endif
1952 
1953  call get_param(param_file, "MOM", "GLOBAL_INDEXING", global_indexing, &
1954  "If true, use a global lateral indexing convention, so "//&
1955  "that corresponding points on different processors have "//&
1956  "the same index. This does not work with static memory.", &
1957  default=.false., layoutparam=.true.)
1958 #ifdef STATIC_MEMORY_
1959  if (global_indexing) call mom_error(fatal, "initialize_MOM: "//&
1960  "GLOBAL_INDEXING can not be true with STATIC_MEMORY.")
1961 #endif
1962  call get_param(param_file, "MOM", "FIRST_DIRECTION", first_direction, &
1963  "An integer that indicates which direction goes first "//&
1964  "in parts of the code that use directionally split "//&
1965  "updates, with even numbers (or 0) used for x- first "//&
1966  "and odd numbers used for y-first.", default=0)
1967 
1968  call get_param(param_file, "MOM", "CHECK_BAD_SURFACE_VALS", cs%check_bad_sfc_vals, &
1969  "If true, check the surface state for ridiculous values.", &
1970  default=.false.)
1971  if (cs%check_bad_sfc_vals) then
1972  call get_param(param_file, "MOM", "BAD_VAL_SSH_MAX", cs%bad_val_ssh_max, &
1973  "The value of SSH above which a bad value message is "//&
1974  "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
1975  units="m", default=20.0, scale=us%m_to_Z)
1976  call get_param(param_file, "MOM", "BAD_VAL_SSS_MAX", cs%bad_val_sss_max, &
1977  "The value of SSS above which a bad value message is "//&
1978  "triggered, if CHECK_BAD_SURFACE_VALS is true.", units="PPT", &
1979  default=45.0)
1980  call get_param(param_file, "MOM", "BAD_VAL_SST_MAX", cs%bad_val_sst_max, &
1981  "The value of SST above which a bad value message is "//&
1982  "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
1983  units="deg C", default=45.0)
1984  call get_param(param_file, "MOM", "BAD_VAL_SST_MIN", cs%bad_val_sst_min, &
1985  "The value of SST below which a bad value message is "//&
1986  "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
1987  units="deg C", default=-2.1)
1988  call get_param(param_file, "MOM", "BAD_VAL_COLUMN_THICKNESS", cs%bad_val_col_thick, &
1989  "The value of column thickness below which a bad value message is "//&
1990  "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
1991  units="m", default=0.0, scale=us%m_to_Z)
1992  endif
1993  call get_param(param_file, "MOM", "DEFAULT_2018_ANSWERS", default_2018_answers, &
1994  "This sets the default value for the various _2018_ANSWERS parameters.", &
1995  default=.false.)
1996  call get_param(param_file, "MOM", "SURFACE_2018_ANSWERS", cs%answers_2018, &
1997  "If true, use expressions for the surface properties that recover the answers "//&
1998  "from the end of 2018. Otherwise, use more appropriate expressions that differ "//&
1999  "at roundoff for non-Boussinesq cases.", default=default_2018_answers)
2000 
2001  call get_param(param_file, "MOM", "SAVE_INITIAL_CONDS", save_ic, &
2002  "If true, write the initial conditions to a file given "//&
2003  "by IC_OUTPUT_FILE.", default=.false.)
2004  call get_param(param_file, "MOM", "IC_OUTPUT_FILE", cs%IC_file, &
2005  "The file into which to write the initial conditions.", &
2006  default="MOM_IC")
2007  call get_param(param_file, "MOM", "WRITE_GEOM", write_geom, &
2008  "If =0, never write the geometry and vertical grid files. "//&
2009  "If =1, write the geometry and vertical grid files only for "//&
2010  "a new simulation. If =2, always write the geometry and "//&
2011  "vertical grid files. Other values are invalid.", default=1)
2012  if (write_geom<0 .or. write_geom>2) call mom_error(fatal,"MOM: "//&
2013  "WRITE_GEOM must be equal to 0, 1 or 2.")
2014  write_geom_files = ((write_geom==2) .or. ((write_geom==1) .and. &
2015  ((dirs%input_filename(1:1)=='n') .and. (len_trim(dirs%input_filename)==1))))
2016 ! If the restart file type had been initialized, this could become:
2017 ! write_geom_files = ((write_geom==2) .or. &
2018 ! ((write_geom==1) .and. is_new_run(restart_CSp)))
2019 
2020  ! Check for inconsistent parameter settings.
2021  if (cs%use_ALE_algorithm .and. bulkmixedlayer) call mom_error(fatal, &
2022  "MOM: BULKMIXEDLAYER can not currently be used with the ALE algorithm.")
2023  if (cs%use_ALE_algorithm .and. .not.use_temperature) call mom_error(fatal, &
2024  "MOM: At this time, USE_EOS should be True when using the ALE algorithm.")
2025  if (cs%adiabatic .and. use_temperature) call mom_error(warning, &
2026  "MOM: ADIABATIC and ENABLE_THERMODYNAMICS both defined is usually unwise.")
2027  if (use_eos .and. .not.use_temperature) call mom_error(fatal, &
2028  "MOM: ENABLE_THERMODYNAMICS must be defined to use USE_EOS.")
2029  if (cs%adiabatic .and. bulkmixedlayer) call mom_error(fatal, &
2030  "MOM: ADIABATIC and BULKMIXEDLAYER can not both be defined.")
2031  if (bulkmixedlayer .and. .not.use_eos) call mom_error(fatal, &
2032  "initialize_MOM: A bulk mixed layer can only be used with T & S as "//&
2033  "state variables. Add USE_EOS = True to MOM_input.")
2034 
2035  call get_param(param_file, 'MOM', "ICE_SHELF", use_ice_shelf, default=.false., do_not_log=.true.)
2036  if (use_ice_shelf) then
2037  inputdir = "." ; call get_param(param_file, 'MOM', "INPUTDIR", inputdir)
2038  inputdir = slasher(inputdir)
2039  call get_param(param_file, 'MOM', "ICE_THICKNESS_FILE", ice_shelf_file, &
2040  "The file from which the ice bathymetry and area are read.", &
2041  fail_if_missing=.true.)
2042  call get_param(param_file, 'MOM', "ICE_AREA_VARNAME", area_varname, &
2043  "The name of the area variable in ICE_THICKNESS_FILE.", &
2044  fail_if_missing=.true.)
2045  endif
2046 
2047 
2048  cs%ensemble_ocean=.false.
2049  call get_param(param_file, "MOM", "ENSEMBLE_OCEAN", cs%ensemble_ocean, &
2050  "If False, The model is being run in serial mode as a single realization. "//&
2051  "If True, The current model realization is part of a larger ensemble "//&
2052  "and at the end of step MOM, we will perform a gather of the ensemble "//&
2053  "members for statistical evaluation and/or data assimilation.", default=.false.)
2054 
2055  call calltree_waypoint("MOM parameters read (initialize_MOM)")
2056 
2057  ! Grid rotation test
2058  call get_param(param_file, "MOM", "ROTATE_INDEX", cs%rotate_index, &
2059  "Enable rotation of the horizontal indices.", default=.false., &
2060  debuggingparam=.true.)
2061  if (cs%rotate_index) then
2062  ! TODO: Index rotation currently only works when index rotation does not
2063  ! change the MPI rank of each domain. Resolving this will require a
2064  ! modification to FMS PE assignment.
2065  ! For now, we only permit single-core runs.
2066 
2067  if (num_pes() /= 1) &
2068  call mom_error(fatal, "Index rotation is only supported on one PE.")
2069 
2070  call get_param(param_file, "MOM", "INDEX_TURNS", turns, &
2071  "Number of counterclockwise quarter-turn index rotations.", &
2072  default=1, debuggingparam=.true.)
2073  endif
2074 
2075  ! Set up the model domain and grids.
2076 #ifdef SYMMETRIC_MEMORY_
2077  symmetric = .true.
2078 #else
2079  symmetric = .false.
2080 #endif
2081  g_in => cs%G_in
2082 #ifdef STATIC_MEMORY_
2083  call mom_domains_init(g_in%domain, param_file, symmetric=symmetric, &
2084  static_memory=.true., nihalo=nihalo_, njhalo=njhalo_, &
2085  niglobal=niglobal_, njglobal=njglobal_, niproc=niproc_, &
2086  njproc=njproc_)
2087 #else
2088  call mom_domains_init(g_in%domain, param_file, symmetric=symmetric, &
2089  domain_name="MOM_in")
2090 #endif
2091 
2092  ! Copy input grid (G_in) domain to active grid G
2093  ! Swap axes for quarter and 3-quarter turns
2094  if (cs%rotate_index) then
2095  allocate(cs%G)
2096  call clone_mom_domain(g_in%Domain, cs%G%Domain, turns=turns)
2097  first_direction = modulo(first_direction + turns, 2)
2098  else
2099  cs%G => g_in
2100  endif
2101 
2102  ! TODO: It is unlikey that test_grid_copy and rotate_index would work at the
2103  ! same time. It may be possible to enable both but for now we prevent it.
2104  if (test_grid_copy .and. cs%rotate_index) &
2105  call mom_error(fatal, "Grid cannot be copied during index rotation.")
2106 
2107  if (test_grid_copy) then ; allocate(g)
2108  else ; g => cs%G ; endif
2109 
2110  call calltree_waypoint("domains initialized (initialize_MOM)")
2111 
2112  call mom_debugging_init(param_file)
2113  call diag_mediator_infrastructure_init()
2114  call mom_io_init(param_file)
2115 
2116  ! Create HI and dG on the input index map.
2117  call hor_index_init(g_in%Domain, hi_in, param_file, &
2118  local_indexing=.not.global_indexing)
2119  call create_dyn_horgrid(dg_in, hi_in, bathymetry_at_vel=bathy_at_vel)
2120  call clone_mom_domain(g_in%Domain, dg_in%Domain)
2121 
2122  ! Allocate initialize time-invariant MOM variables.
2123  call mom_initialize_fixed(dg_in, us, obc_in, param_file, write_geom_files, &
2124  dirs%output_directory)
2125 
2126  call calltree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)")
2127 
2128  ! Determine HI and dG for the model index map.
2129  if (cs%rotate_index) then
2130  allocate(hi)
2131  call rotate_hor_index(hi_in, turns, hi)
2132  call create_dyn_horgrid(dg, hi, bathymetry_at_vel=bathy_at_vel)
2133  call clone_mom_domain(g%Domain, dg%Domain)
2134  call rotate_dyngrid(dg_in, dg, us, turns)
2135  if (associated(obc_in)) then
2136  ! TODO: General OBC index rotations is not yet supported.
2137  if (modulo(turns, 4) /= 1) &
2138  call mom_error(fatal, "OBC index rotation of 180 and 270 degrees is " &
2139  // "not yet unsupported.")
2140  allocate(cs%OBC)
2141  call rotate_obc_config(obc_in, dg_in, cs%OBC, dg, turns)
2142  endif
2143  else
2144  hi => hi_in
2145  dg => dg_in
2146  cs%OBC => obc_in
2147  endif
2148 
2149  call verticalgridinit( param_file, cs%GV, us )
2150  gv => cs%GV
2151 
2152  ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
2153  if (cs%debug .or. dg%symmetric) &
2154  call clone_mom_domain(dg%Domain, dg%Domain_aux, symmetric=.false.)
2155 
2156  call calltree_waypoint("grids initialized (initialize_MOM)")
2157 
2158  call mom_timing_init(cs)
2159 
2160  if (associated(cs%OBC)) call call_obc_register(param_file, cs%update_OBC_CSp, cs%OBC)
2161 
2162  call tracer_registry_init(param_file, cs%tracer_Reg)
2163 
2164  ! Allocate and initialize space for the primary time-varying MOM variables.
2165  is = dg%isc ; ie = dg%iec ; js = dg%jsc ; je = dg%jec ; nz = gv%ke
2166  isd = dg%isd ; ied = dg%ied ; jsd = dg%jsd ; jed = dg%jed
2167  isdb = dg%IsdB ; iedb = dg%IedB ; jsdb = dg%JsdB ; jedb = dg%JedB
2168  alloc_(cs%u(isdb:iedb,jsd:jed,nz)) ; cs%u(:,:,:) = 0.0
2169  alloc_(cs%v(isd:ied,jsdb:jedb,nz)) ; cs%v(:,:,:) = 0.0
2170  alloc_(cs%h(isd:ied,jsd:jed,nz)) ; cs%h(:,:,:) = gv%Angstrom_H
2171  alloc_(cs%uh(isdb:iedb,jsd:jed,nz)) ; cs%uh(:,:,:) = 0.0
2172  alloc_(cs%vh(isd:ied,jsdb:jedb,nz)) ; cs%vh(:,:,:) = 0.0
2173  if (use_temperature) then
2174  alloc_(cs%T(isd:ied,jsd:jed,nz)) ; cs%T(:,:,:) = 0.0
2175  alloc_(cs%S(isd:ied,jsd:jed,nz)) ; cs%S(:,:,:) = 0.0
2176  cs%tv%T => cs%T ; cs%tv%S => cs%S
2177  if (cs%tv%T_is_conT) then
2178  vd_t = var_desc(name="contemp", units="Celsius", longname="Conservative Temperature", &
2179  cmor_field_name="thetao", cmor_longname="Sea Water Potential Temperature", &
2180  conversion=us%Q_to_J_kg*cs%tv%C_p)
2181  else
2182  vd_t = var_desc(name="temp", units="degC", longname="Potential Temperature", &
2183  cmor_field_name="thetao", cmor_longname="Sea Water Potential Temperature", &
2184  conversion=us%Q_to_J_kg*cs%tv%C_p)
2185  endif
2186  if (cs%tv%S_is_absS) then
2187  vd_s = var_desc(name="abssalt",units="g kg-1",longname="Absolute Salinity", &
2188  cmor_field_name="so", cmor_longname="Sea Water Salinity", &
2189  conversion=0.001)
2190  else
2191  vd_s = var_desc(name="salt",units="psu",longname="Salinity", &
2192  cmor_field_name="so", cmor_longname="Sea Water Salinity", &
2193  conversion=0.001)
2194  endif
2195 
2196  if (advect_ts) then
2197  s_flux_units = get_tr_flux_units(gv, "psu") ! Could change to "kg m-2 s-1"?
2198  conv2watt = gv%H_to_kg_m2 * us%Q_to_J_kg*cs%tv%C_p
2199  if (gv%Boussinesq) then
2200  conv2salt = gv%H_to_m ! Could change to GV%H_to_kg_m2 * 0.001?
2201  else
2202  conv2salt = gv%H_to_kg_m2
2203  endif
2204  call register_tracer(cs%tv%T, cs%tracer_Reg, param_file, dg%HI, gv, &
2205  tr_desc=vd_t, registry_diags=.true., flux_nameroot='T', &
2206  flux_units='W', flux_longname='Heat', &
2207  flux_scale=conv2watt, convergence_units='W m-2', &
2208  convergence_scale=conv2watt, cmor_tendprefix="opottemp", diag_form=2)
2209  call register_tracer(cs%tv%S, cs%tracer_Reg, param_file, dg%HI, gv, &
2210  tr_desc=vd_s, registry_diags=.true., flux_nameroot='S', &
2211  flux_units=s_flux_units, flux_longname='Salt', &
2212  flux_scale=conv2salt, convergence_units='kg m-2 s-1', &
2213  convergence_scale=0.001*gv%H_to_kg_m2, cmor_tendprefix="osalt", diag_form=2)
2214  endif
2215  ! NOTE: register_temp_salt_segments includes allocation of tracer fields
2216  ! along segments. Bit reproducibility requires that MOM_initialize_state
2217  ! be called on the input index map, so we must setup both OBC and OBC_in.
2218  !
2219  ! XXX: This call on OBC_in allocates the tracer fields on the unrotated
2220  ! grid, but also incorrectly stores a pointer to a tracer_type for the
2221  ! rotated registry (e.g. segment%tr_reg%Tr(n)%Tr) from CS%tracer_reg.
2222  !
2223  ! While incorrect and potentially dangerous, it does not seem that this
2224  ! pointer is used during initialization, so we leave it for now.
2225  if (cs%rotate_index .and. associated(obc_in)) &
2226  call register_temp_salt_segments(gv, obc_in, cs%tracer_Reg, param_file)
2227  if (associated(cs%OBC)) &
2228  call register_temp_salt_segments(gv, cs%OBC, cs%tracer_Reg, param_file)
2229  endif
2230  if (use_frazil) then
2231  allocate(cs%tv%frazil(isd:ied,jsd:jed)) ; cs%tv%frazil(:,:) = 0.0
2232  endif
2233  if (bound_salinity) then
2234  allocate(cs%tv%salt_deficit(isd:ied,jsd:jed)) ; cs%tv%salt_deficit(:,:) = 0.0
2235  endif
2236 
2237  if (bulkmixedlayer .or. use_temperature) then
2238  allocate(cs%Hml(isd:ied,jsd:jed)) ; cs%Hml(:,:) = 0.0
2239  endif
2240 
2241  if (bulkmixedlayer) then
2242  gv%nkml = nkml ; gv%nk_rho_varies = nkml + nkbl
2243  else
2244  gv%nkml = 0 ; gv%nk_rho_varies = 0
2245  endif
2246  if (cs%use_ALE_algorithm) then
2247  call get_param(param_file, "MOM", "NK_RHO_VARIES", gv%nk_rho_varies, default=0) ! Will default to nz later... -AJA
2248  endif
2249 
2250  alloc_(cs%uhtr(isdb:iedb,jsd:jed,nz)) ; cs%uhtr(:,:,:) = 0.0
2251  alloc_(cs%vhtr(isd:ied,jsdb:jedb,nz)) ; cs%vhtr(:,:,:) = 0.0
2252  cs%t_dyn_rel_adv = 0.0 ; cs%t_dyn_rel_thermo = 0.0 ; cs%t_dyn_rel_diag = 0.0
2253 
2254  if (debug_truncations) then
2255  allocate(cs%u_prev(isdb:iedb,jsd:jed,nz)) ; cs%u_prev(:,:,:) = 0.0
2256  allocate(cs%v_prev(isd:ied,jsdb:jedb,nz)) ; cs%v_prev(:,:,:) = 0.0
2257  mom_internal_state%u_prev => cs%u_prev
2258  mom_internal_state%v_prev => cs%v_prev
2259  call safe_alloc_ptr(cs%ADp%du_dt_visc,isdb,iedb,jsd,jed,nz)
2260  call safe_alloc_ptr(cs%ADp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
2261  if (.not.cs%adiabatic) then
2262  call safe_alloc_ptr(cs%ADp%du_dt_dia,isdb,iedb,jsd,jed,nz)
2263  call safe_alloc_ptr(cs%ADp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
2264  endif
2265  endif
2266 
2267  mom_internal_state%u => cs%u ; mom_internal_state%v => cs%v
2268  mom_internal_state%h => cs%h
2269  mom_internal_state%uh => cs%uh ; mom_internal_state%vh => cs%vh
2270  if (use_temperature) then
2271  mom_internal_state%T => cs%T ; mom_internal_state%S => cs%S
2272  endif
2273 
2274  cs%CDp%uh => cs%uh ; cs%CDp%vh => cs%vh
2275 
2276  if (cs%interp_p_surf) then
2277  allocate(cs%p_surf_prev(isd:ied,jsd:jed)) ; cs%p_surf_prev(:,:) = 0.0
2278  endif
2279 
2280  alloc_(cs%ssh_rint(isd:ied,jsd:jed)) ; cs%ssh_rint(:,:) = 0.0
2281  alloc_(cs%ave_ssh_ibc(isd:ied,jsd:jed)) ; cs%ave_ssh_ibc(:,:) = 0.0
2282  alloc_(cs%eta_av_bc(isd:ied,jsd:jed)) ; cs%eta_av_bc(:,:) = 0.0
2283  cs%time_in_cycle = 0.0 ; cs%time_in_thermo_cycle = 0.0
2284 
2285  ! Use the Wright equation of state by default, unless otherwise specified
2286  ! Note: this line and the following block ought to be in a separate
2287  ! initialization routine for tv.
2288  if (use_eos) call eos_init(param_file, cs%tv%eqn_of_state, us)
2289  if (use_temperature) then
2290  allocate(cs%tv%TempxPmE(isd:ied,jsd:jed)) ; cs%tv%TempxPmE(:,:) = 0.0
2291  if (use_geothermal) then
2292  allocate(cs%tv%internal_heat(isd:ied,jsd:jed)) ; cs%tv%internal_heat(:,:) = 0.0
2293  endif
2294  endif
2295  call calltree_waypoint("state variables allocated (initialize_MOM)")
2296 
2297  ! Set the fields that are needed for bitwise identical restarting
2298  ! the time stepping scheme.
2299  call restart_init(param_file, restart_csp)
2300  call set_restart_fields(gv, us, param_file, cs, restart_csp)
2301  if (cs%split) then
2302  call register_restarts_dyn_split_rk2(dg%HI, gv, param_file, &
2303  cs%dyn_split_RK2_CSp, restart_csp, cs%uh, cs%vh)
2304  elseif (cs%use_RK2) then
2305  call register_restarts_dyn_unsplit_rk2(dg%HI, gv, param_file, &
2306  cs%dyn_unsplit_RK2_CSp, restart_csp)
2307  else
2308  call register_restarts_dyn_unsplit(dg%HI, gv, param_file, &
2309  cs%dyn_unsplit_CSp, restart_csp)
2310  endif
2311 
2312  ! This subroutine calls user-specified tracer registration routines.
2313  ! Additional calls can be added to MOM_tracer_flow_control.F90.
2314  call call_tracer_register(dg%HI, gv, us, param_file, cs%tracer_flow_CSp, &
2315  cs%tracer_Reg, restart_csp)
2316 
2317  call meke_alloc_register_restart(dg%HI, param_file, cs%MEKE, restart_csp)
2318  call set_visc_register_restarts(dg%HI, gv, param_file, cs%visc, restart_csp)
2319  call mixedlayer_restrat_register_restarts(dg%HI, param_file, &
2320  cs%mixedlayer_restrat_CSp, restart_csp)
2321 
2322  if (associated(cs%OBC)) &
2323  call open_boundary_register_restarts(dg%HI, gv, cs%OBC, cs%tracer_Reg, &
2324  param_file, restart_csp, use_temperature)
2325 
2326  call calltree_waypoint("restart registration complete (initialize_MOM)")
2327 
2328  ! Initialize dynamically evolving fields, perhaps from restart files.
2329  call cpu_clock_begin(id_clock_mom_init)
2330  call mom_initialize_coord(gv, us, param_file, write_geom_files, &
2331  dirs%output_directory, cs%tv, dg%max_depth)
2332  call calltree_waypoint("returned from MOM_initialize_coord() (initialize_MOM)")
2333 
2334  if (cs%use_ALE_algorithm) then
2335  call ale_init(param_file, gv, us, dg%max_depth, cs%ALE_CSp)
2336  call calltree_waypoint("returned from ALE_init() (initialize_MOM)")
2337  endif
2338 
2339  ! Shift from using the temporary dynamic grid type to using the final
2340  ! (potentially static) ocean-specific grid type.
2341  ! The next line would be needed if G%Domain had not already been init'd above:
2342  ! call clone_MOM_domain(dG%Domain, G%Domain)
2343 
2344  ! NOTE: If indices are rotated, then G and G_in must both be initialized.
2345  ! If not rotated, then G_in and G are the same grid.
2346  if (cs%rotate_index) then
2347  call mom_grid_init(g, param_file, us, hi, bathymetry_at_vel=bathy_at_vel)
2348  call copy_dyngrid_to_mom_grid(dg, g, us)
2349  call destroy_dyn_horgrid(dg)
2350  endif
2351  call mom_grid_init(g_in, param_file, us, hi_in, bathymetry_at_vel=bathy_at_vel)
2352  call copy_dyngrid_to_mom_grid(dg_in, g_in, us)
2353  call destroy_dyn_horgrid(dg_in)
2354 
2355  ! Set a few remaining fields that are specific to the ocean grid type.
2356  call set_first_direction(g, first_direction)
2357  ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
2358  if (cs%debug .or. g%symmetric) then
2359  call clone_mom_domain(g%Domain, g%Domain_aux, symmetric=.false.)
2360  else ; g%Domain_aux => g%Domain ; endif
2361  ! Copy common variables from the vertical grid to the horizontal grid.
2362  ! Consider removing this later?
2363  g%ke = gv%ke
2364 
2365  if (cs%rotate_index) then
2366  g_in%ke = gv%ke
2367 
2368  allocate(u_in(g_in%IsdB:g_in%IedB, g_in%jsd:g_in%jed, nz))
2369  allocate(v_in(g_in%isd:g_in%ied, g_in%JsdB:g_in%JedB, nz))
2370  allocate(h_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed, nz))
2371  u_in(:,:,:) = 0.0
2372  v_in(:,:,:) = 0.0
2373  h_in(:,:,:) = gv%Angstrom_H
2374 
2375  if (use_temperature) then
2376  allocate(t_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed, nz))
2377  allocate(s_in(g_in%isd:g_in%ied, g_in%jsd:g_in%jed, nz))
2378  t_in(:,:,:) = 0.0
2379  s_in(:,:,:) = 0.0
2380 
2381  cs%tv%T => t_in
2382  cs%tv%S => s_in
2383  endif
2384 
2385  call mom_initialize_state(u_in, v_in, h_in, cs%tv, time, g_in, gv, us, &
2386  param_file, dirs, restart_csp, cs%ALE_CSp, cs%tracer_Reg, &
2387  sponge_in_csp, ale_sponge_in_csp, obc_in, time_in)
2388 
2389  if (use_temperature) then
2390  cs%tv%T => cs%T
2391  cs%tv%S => cs%S
2392  endif
2393 
2394  call rotate_initial_state(u_in, v_in, h_in, t_in, s_in, use_temperature, &
2395  turns, cs%u, cs%v, cs%h, cs%T, cs%S)
2396 
2397  if (associated(sponge_in_csp)) then
2398  ! TODO: Implementation and testing of non-ALE spong rotation
2399  call mom_error(fatal, "Index rotation of non-ALE sponge is not yet implemented.")
2400  endif
2401 
2402  if (associated(ale_sponge_in_csp)) then
2403  call rotate_ale_sponge(ale_sponge_in_csp, g_in, cs%ALE_sponge_CSp, g, turns, param_file)
2404  call update_ale_sponge_field(cs%ALE_sponge_CSp, t_in, g, gv, cs%T)
2405  call update_ale_sponge_field(cs%ALE_sponge_CSp, s_in, g, gv, cs%S)
2406  endif
2407 
2408  if (associated(obc_in)) &
2409  call rotate_obc_init(obc_in, g, gv, us, param_file, cs%tv, restart_csp, cs%OBC)
2410 
2411  deallocate(u_in)
2412  deallocate(v_in)
2413  deallocate(h_in)
2414  if (use_temperature) then
2415  deallocate(t_in)
2416  deallocate(s_in)
2417  endif
2418  else
2419  call mom_initialize_state(cs%u, cs%v, cs%h, cs%tv, time, g, gv, us, &
2420  param_file, dirs, restart_csp, cs%ALE_CSp, cs%tracer_Reg, &
2421  cs%sponge_CSp, cs%ALE_sponge_CSp, cs%OBC, time_in)
2422  endif
2423 
2424  call cpu_clock_end(id_clock_mom_init)
2425  call calltree_waypoint("returned from MOM_initialize_state() (initialize_MOM)")
2426 
2427 ! ! Need this after MOM_initialize_state for DOME OBC stuff.
2428 ! if (associated(CS%OBC)) &
2429 ! call open_boundary_register_restarts(G%HI, GV, CS%OBC, CS%tracer_Reg, &
2430 ! param_file, restart_CSp, use_temperature)
2431 
2432 ! call callTree_waypoint("restart registration complete (initialize_MOM)")
2433 
2434  ! From this point, there may be pointers being set, so the final grid type
2435  ! that will persist throughout the run has to be used.
2436 
2437  if (test_grid_copy) then
2438  ! Copy the data from the temporary grid to the dyn_hor_grid to CS%G.
2439  call create_dyn_horgrid(dg, g%HI)
2440  call clone_mom_domain(g%Domain, dg%Domain)
2441 
2442  call clone_mom_domain(g%Domain, cs%G%Domain)
2443  call mom_grid_init(cs%G, param_file, us)
2444 
2445  call copy_mom_grid_to_dyngrid(g, dg, us)
2446  call copy_dyngrid_to_mom_grid(dg, cs%G, us)
2447 
2448  call destroy_dyn_horgrid(dg)
2449  call mom_grid_end(g) ; deallocate(g)
2450 
2451  g => cs%G
2452  if (cs%debug .or. cs%G%symmetric) then
2453  call clone_mom_domain(cs%G%Domain, cs%G%Domain_aux, symmetric=.false.)
2454  else ; cs%G%Domain_aux => cs%G%Domain ; endif
2455  g%ke = gv%ke
2456  endif
2457 
2458  ! At this point, all user-modified initialization code has been called. The
2459  ! remainder of this subroutine is controlled by the parameters that have
2460  ! have already been set.
2461 
2462  if (ale_remap_init_conds(cs%ALE_CSp) .and. .not. query_initialized(cs%h,"h",restart_csp)) then
2463  ! This block is controlled by the ALE parameter REMAP_AFTER_INITIALIZATION.
2464  ! \todo This block exists for legacy reasons and we should phase it out of
2465  ! all examples. !###
2466  if (cs%debug) then
2467  call uvchksum("Pre ALE adjust init cond [uv]", &
2468  cs%u, cs%v, g%HI, haloshift=1)
2469  call hchksum(cs%h,"Pre ALE adjust init cond h", g%HI, haloshift=1, scale=gv%H_to_m)
2470  endif
2471  call calltree_waypoint("Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)")
2472  call adjustgridforintegrity(cs%ALE_CSp, g, gv, cs%h )
2473  call calltree_waypoint("Calling ALE_main() to remap initial conditions (initialize_MOM)")
2474  if (use_ice_shelf) then
2475  filename = trim(inputdir)//trim(ice_shelf_file)
2476  if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
2477  "MOM: Unable to open "//trim(filename))
2478 
2479  allocate(area_shelf_h(isd:ied,jsd:jed))
2480  allocate(frac_shelf_h(isd:ied,jsd:jed))
2481  call mom_read_data(filename, trim(area_varname), area_shelf_h, g%Domain, scale=us%m_to_L**2)
2482  ! initialize frac_shelf_h with zeros (open water everywhere)
2483  frac_shelf_h(:,:) = 0.0
2484  ! compute fractional ice shelf coverage of h
2485  do j=jsd,jed ; do i=isd,ied
2486  if (g%areaT(i,j) > 0.0) &
2487  frac_shelf_h(i,j) = area_shelf_h(i,j) / g%areaT(i,j)
2488  enddo ; enddo
2489  ! pass to the pointer
2490  shelf_area => frac_shelf_h
2491  call ale_main(g, gv, us, cs%h, cs%u, cs%v, cs%tv, cs%tracer_Reg, cs%ALE_CSp, &
2492  cs%OBC, frac_shelf_h=shelf_area)
2493  else
2494  call ale_main( g, gv, us, cs%h, cs%u, cs%v, cs%tv, cs%tracer_Reg, cs%ALE_CSp, cs%OBC)
2495  endif
2496 
2497  call cpu_clock_begin(id_clock_pass_init)
2498  call create_group_pass(tmp_pass_uv_t_s_h, cs%u, cs%v, g%Domain)
2499  if (use_temperature) then
2500  call create_group_pass(tmp_pass_uv_t_s_h, cs%tv%T, g%Domain, halo=1)
2501  call create_group_pass(tmp_pass_uv_t_s_h, cs%tv%S, g%Domain, halo=1)
2502  endif
2503  call create_group_pass(tmp_pass_uv_t_s_h, cs%h, g%Domain, halo=1)
2504  call do_group_pass(tmp_pass_uv_t_s_h, g%Domain)
2505  call cpu_clock_end(id_clock_pass_init)
2506 
2507  if (cs%debug) then
2508  call uvchksum("Post ALE adjust init cond [uv]", cs%u, cs%v, g%HI, haloshift=1)
2509  call hchksum(cs%h, "Post ALE adjust init cond h", g%HI, haloshift=1, scale=gv%H_to_m)
2510  endif
2511  endif
2512  if ( cs%use_ALE_algorithm ) call ale_updateverticalgridtype( cs%ALE_CSp, gv )
2513 
2514  diag => cs%diag
2515  ! Initialize the diag mediator.
2516  call diag_mediator_init(g, gv, us, gv%ke, param_file, diag, doc_file_dir=dirs%output_directory)
2517  if (present(diag_ptr)) diag_ptr => cs%diag
2518 
2519  ! Initialize the diagnostics masks for native arrays.
2520  ! This step has to be done after call to MOM_initialize_state
2521  ! and before MOM_diagnostics_init
2522  call diag_masks_set(g, gv%ke, diag)
2523 
2524  ! Set up pointers within diag mediator control structure,
2525  ! this needs to occur _after_ CS%h etc. have been allocated.
2526  call diag_set_state_ptrs(cs%h, cs%T, cs%S, cs%tv%eqn_of_state, diag)
2527 
2528  ! This call sets up the diagnostic axes. These are needed,
2529  ! e.g. to generate the target grids below.
2530  call set_axes_info(g, gv, us, param_file, diag)
2531 
2532  ! Whenever thickness/T/S changes let the diag manager know, target grids
2533  ! for vertical remapping may need to be regenerated.
2534  ! FIXME: are h, T, S updated at the same time? Review these for T, S updates.
2535  call diag_update_remap_grids(diag)
2536 
2537  ! Setup the diagnostic grid storage types
2538  call diag_grid_storage_init(cs%diag_pre_sync, g, diag)
2539  call diag_grid_storage_init(cs%diag_pre_dyn, g, diag)
2540 
2541  ! Calculate masks for diagnostics arrays in non-native coordinates
2542  ! This step has to be done after set_axes_info() because the axes needed
2543  ! to be configured, and after diag_update_remap_grids() because the grids
2544  ! must be defined.
2545  call set_masks_for_axes(g, diag)
2546 
2547  ! Diagnose static fields AND associate areas/volumes with axes
2548  call write_static_fields(g, gv, us, cs%tv, cs%diag)
2549  call calltree_waypoint("static fields written (initialize_MOM)")
2550 
2551  ! Register the volume cell measure (must be one of first diagnostics)
2552  call register_cell_measure(g, cs%diag, time)
2553 
2554  call cpu_clock_begin(id_clock_mom_init)
2555  if (cs%use_ALE_algorithm) then
2556  call ale_writecoordinatefile( cs%ALE_CSp, gv, dirs%output_directory )
2557  endif
2558  call cpu_clock_end(id_clock_mom_init)
2559  call calltree_waypoint("ALE initialized (initialize_MOM)")
2560 
2561  cs%useMEKE = meke_init(time, g, us, param_file, diag, cs%MEKE_CSp, cs%MEKE, restart_csp)
2562 
2563  call varmix_init(time, g, gv, us, param_file, diag, cs%VarMix)
2564  call set_visc_init(time, g, gv, us, param_file, diag, cs%visc, cs%set_visc_CSp, restart_csp, cs%OBC)
2565  call thickness_diffuse_init(time, g, gv, us, param_file, diag, cs%CDp, cs%thickness_diffuse_CSp)
2566 
2567  if (cs%split) then
2568  allocate(eta(szi_(g),szj_(g))) ; eta(:,:) = 0.0
2569  call initialize_dyn_split_rk2(cs%u, cs%v, cs%h, cs%uh, cs%vh, eta, time, &
2570  g, gv, us, param_file, diag, cs%dyn_split_RK2_CSp, restart_csp, &
2571  cs%dt, cs%ADp, cs%CDp, mom_internal_state, cs%VarMix, cs%MEKE, &
2572  cs%thickness_diffuse_CSp, &
2573  cs%OBC, cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, &
2574  cs%visc, dirs, cs%ntrunc, calc_dtbt=calc_dtbt, cont_stencil=cs%cont_stencil)
2575  if (cs%dtbt_reset_period > 0.0) then
2576  cs%dtbt_reset_interval = real_to_time(cs%dtbt_reset_period)
2577  ! Set dtbt_reset_time to be the next even multiple of dtbt_reset_interval.
2578  cs%dtbt_reset_time = time_init + cs%dtbt_reset_interval * &
2579  ((time - time_init) / cs%dtbt_reset_interval)
2580  if ((cs%dtbt_reset_time > time) .and. calc_dtbt) then
2581  ! Back up dtbt_reset_time one interval to force dtbt to be calculated,
2582  ! because the restart was not aligned with the interval to recalculate
2583  ! dtbt, and dtbt was not read from a restart file.
2584  cs%dtbt_reset_time = cs%dtbt_reset_time - cs%dtbt_reset_interval
2585  endif
2586  endif
2587  elseif (cs%use_RK2) then
2588  call initialize_dyn_unsplit_rk2(cs%u, cs%v, cs%h, time, g, gv, us, &
2589  param_file, diag, cs%dyn_unsplit_RK2_CSp, restart_csp, &
2590  cs%ADp, cs%CDp, mom_internal_state, cs%MEKE, cs%OBC, &
2591  cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, cs%visc, dirs, &
2592  cs%ntrunc, cont_stencil=cs%cont_stencil)
2593  else
2594  call initialize_dyn_unsplit(cs%u, cs%v, cs%h, time, g, gv, us, &
2595  param_file, diag, cs%dyn_unsplit_CSp, restart_csp, &
2596  cs%ADp, cs%CDp, mom_internal_state, cs%MEKE, cs%OBC, &
2597  cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, cs%visc, dirs, &
2598  cs%ntrunc, cont_stencil=cs%cont_stencil)
2599  endif
2600 
2601  call calltree_waypoint("dynamics initialized (initialize_MOM)")
2602 
2603  cs%mixedlayer_restrat = mixedlayer_restrat_init(time, g, gv, us, param_file, diag, &
2604  cs%mixedlayer_restrat_CSp, restart_csp)
2605  if (cs%mixedlayer_restrat) then
2606  if (.not.(bulkmixedlayer .or. cs%use_ALE_algorithm)) &
2607  call mom_error(fatal, "MOM: MIXEDLAYER_RESTRAT true requires a boundary layer scheme.")
2608  ! When DIABATIC_FIRST=False and using CS%visc%ML in mixedlayer_restrat we need to update after a restart
2609  if (.not. cs%diabatic_first .and. associated(cs%visc%MLD)) &
2610  call pass_var(cs%visc%MLD, g%domain, halo=1)
2611  endif
2612 
2613  call mom_diagnostics_init(mom_internal_state, cs%ADp, cs%CDp, time, g, gv, us, &
2614  param_file, diag, cs%diagnostics_CSp, cs%tv)
2615  call diag_copy_diag_to_storage(cs%diag_pre_sync, cs%h, cs%diag)
2616 
2617  if (associated(cs%sponge_CSp)) &
2618  call init_sponge_diags(time, g, gv, us, diag, cs%sponge_CSp)
2619 
2620  if (associated(cs%ALE_sponge_CSp)) &
2621  call init_ale_sponge_diags(time, g, diag, cs%ALE_sponge_CSp)
2622 
2623  if (cs%adiabatic) then
2624  call adiabatic_driver_init(time, g, param_file, diag, cs%diabatic_CSp, &
2625  cs%tracer_flow_CSp)
2626  else
2627  call diabatic_driver_init(time, g, gv, us, param_file, cs%use_ALE_algorithm, diag, &
2628  cs%ADp, cs%CDp, cs%diabatic_CSp, cs%tracer_flow_CSp, &
2629  cs%sponge_CSp, cs%ALE_sponge_CSp)
2630  endif
2631 
2632  call tracer_advect_init(time, g, us, param_file, diag, cs%tracer_adv_CSp)
2633  call tracer_hor_diff_init(time, g, us, param_file, diag, cs%tv%eqn_of_state, cs%diabatic_CSp, &
2634  cs%tracer_diff_CSp)
2635 
2636  call lock_tracer_registry(cs%tracer_Reg)
2637  call calltree_waypoint("tracer registry now locked (initialize_MOM)")
2638 
2639  ! now register some diagnostics since the tracer registry is now locked
2640  call register_surface_diags(time, g, us, cs%sfc_IDs, cs%diag, cs%tv)
2641  call register_diags(time, g, gv, us, cs%IDs, cs%diag)
2642  call register_transport_diags(time, g, gv, us, cs%transport_IDs, cs%diag)
2643  call register_tracer_diagnostics(cs%tracer_Reg, cs%h, time, diag, g, gv, us, &
2644  cs%use_ALE_algorithm)
2645  if (cs%use_ALE_algorithm) then
2646  call ale_register_diags(time, g, gv, us, diag, cs%ALE_CSp)
2647  endif
2648 
2649  ! This subroutine initializes any tracer packages.
2650  new_sim = is_new_run(restart_csp)
2651  call tracer_flow_control_init(.not.new_sim, time, g, gv, us, cs%h, param_file, &
2652  cs%diag, cs%OBC, cs%tracer_flow_CSp, cs%sponge_CSp, &
2653  cs%ALE_sponge_CSp, cs%tv)
2654  if (present(tracer_flow_csp)) tracer_flow_csp => cs%tracer_flow_CSp
2655 
2656  ! If running in offline tracer mode, initialize the necessary control structure and
2657  ! parameters
2658  if (present(offline_tracer_mode)) offline_tracer_mode=cs%offline_tracer_mode
2659 
2660  if (cs%offline_tracer_mode) then
2661  ! Setup some initial parameterizations and also assign some of the subtypes
2662  call offline_transport_init(param_file, cs%offline_CSp, cs%diabatic_CSp, g, gv, us)
2663  call insert_offline_main( cs=cs%offline_CSp, ale_csp=cs%ALE_CSp, diabatic_csp=cs%diabatic_CSp, &
2664  diag=cs%diag, obc=cs%OBC, tracer_adv_csp=cs%tracer_adv_CSp, &
2665  tracer_flow_csp=cs%tracer_flow_CSp, tracer_reg=cs%tracer_Reg, &
2666  tv=cs%tv, x_before_y = (mod(first_direction,2)==0), debug=cs%debug )
2667  call register_diags_offline_transport(time, cs%diag, cs%offline_CSp)
2668  endif
2669 
2670  !--- set up group pass for u,v,T,S and h. pass_uv_T_S_h also is used in step_MOM
2671  call cpu_clock_begin(id_clock_pass_init)
2672  dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
2673  call create_group_pass(pass_uv_t_s_h, cs%u, cs%v, g%Domain, halo=dynamics_stencil)
2674  if (use_temperature) then
2675  call create_group_pass(pass_uv_t_s_h, cs%tv%T, g%Domain, halo=dynamics_stencil)
2676  call create_group_pass(pass_uv_t_s_h, cs%tv%S, g%Domain, halo=dynamics_stencil)
2677  endif
2678  call create_group_pass(pass_uv_t_s_h, cs%h, g%Domain, halo=dynamics_stencil)
2679 
2680  call do_group_pass(pass_uv_t_s_h, g%Domain)
2681 
2682  if (associated(cs%visc%Kv_shear)) &
2683  call pass_var(cs%visc%Kv_shear, g%Domain, to_all+omit_corners, halo=1)
2684 
2685  if (associated(cs%visc%Kv_slow)) &
2686  call pass_var(cs%visc%Kv_slow, g%Domain, to_all+omit_corners, halo=1)
2687 
2688  call cpu_clock_end(id_clock_pass_init)
2689 
2690  call register_obsolete_diagnostics(param_file, cs%diag)
2691 
2692  if (use_frazil) then
2693  if (query_initialized(cs%tv%frazil, "frazil", restart_csp)) then
2694  ! Test whether the dimensional rescaling has changed for heat content.
2695  if ((us%kg_m3_to_R_restart*us%m_to_Z_restart*us%J_kg_to_Q_restart /= 0.0) .and. &
2696  ((us%J_kg_to_Q*us%kg_m3_to_R*us%m_to_Z) /= &
2697  (us%J_kg_to_Q_restart*us%kg_m3_to_R_restart*us%m_to_Z_restart)) ) then
2698  qrz_rescale = (us%J_kg_to_Q*us%kg_m3_to_R*us%m_to_Z) / &
2699  (us%J_kg_to_Q_restart*us%kg_m3_to_R_restart*us%m_to_Z_restart)
2700  do j=js,je ; do i=is,ie
2701  cs%tv%frazil(i,j) = qrz_rescale * cs%tv%frazil(i,j)
2702  enddo ; enddo
2703  endif
2704  else
2705  cs%tv%frazil(:,:) = 0.0
2706  endif
2707  endif
2708 
2709  if (cs%interp_p_surf) then
2710  cs%p_surf_prev_set = query_initialized(cs%p_surf_prev,"p_surf_prev",restart_csp)
2711 
2712  if (cs%p_surf_prev_set) then
2713  ! Test whether the dimensional rescaling has changed for pressure.
2714  if ((us%kg_m3_to_R_restart*us%s_to_T_restart*us%m_to_L_restart /= 0.0) .and. &
2715  ((us%kg_m3_to_R*(us%m_to_L*us%s_to_T_restart)**2) /= &
2716  (us%kg_m3_to_R_restart*(us%m_to_L_restart*us%s_to_T)**2)) ) then
2717  rl2_t2_rescale = (us%kg_m3_to_R*(us%m_to_L*us%s_to_T_restart)**2) / &
2718  (us%kg_m3_to_R_restart*(us%m_to_L_restart*us%s_to_T)**2)
2719  do j=js,je ; do i=is,ie
2720  cs%p_surf_prev(i,j) = rl2_t2_rescale * cs%p_surf_prev(i,j)
2721  enddo ; enddo
2722  endif
2723 
2724  call pass_var(cs%p_surf_prev, g%domain)
2725  endif
2726  endif
2727 
2728  if (use_ice_shelf .and. associated(cs%Hml)) then
2729  if (query_initialized(cs%Hml, "hML", restart_csp)) then
2730  ! Test whether the dimensional rescaling has changed for depths.
2731  if ((us%m_to_Z_restart /= 0.0) .and. (us%m_to_Z /= us%m_to_Z_restart) ) then
2732  z_rescale = us%m_to_Z / us%m_to_Z_restart
2733  do j=js,je ; do i=is,ie
2734  cs%Hml(i,j) = z_rescale * cs%Hml(i,j)
2735  enddo ; enddo
2736  endif
2737  endif
2738  endif
2739 
2740  if (.not.query_initialized(cs%ave_ssh_ibc,"ave_ssh",restart_csp)) then
2741  if (cs%split) then
2742  call find_eta(cs%h, cs%tv, g, gv, us, cs%ave_ssh_ibc, eta, eta_to_m=1.0)
2743  else
2744  call find_eta(cs%h, cs%tv, g, gv, us, cs%ave_ssh_ibc, eta_to_m=1.0)
2745  endif
2746  endif
2747  if (cs%split) deallocate(eta)
2748 
2749  cs%nstep_tot = 0
2750  if (present(count_calls)) cs%count_calls = count_calls
2751  call mom_sum_output_init(g_in, us, param_file, dirs%output_directory, &
2752  cs%ntrunc, time_init, cs%sum_output_CSp)
2753 
2754  ! Flag whether to save initial conditions in finish_MOM_initialization() or not.
2755  cs%write_IC = save_ic .and. &
2756  .not.((dirs%input_filename(1:1) == 'r') .and. &
2757  (len_trim(dirs%input_filename) == 1))
2758 
2759  if (cs%ensemble_ocean) then
2760  call init_oda(time, g, gv, cs%odaCS)
2761  endif
2762 
2763  !### This could perhaps go here instead of in finish_MOM_initialization?
2764  ! call fix_restart_scaling(GV)
2765  ! call fix_restart_unit_scaling(US)
2766 
2767  call calltree_leave("initialize_MOM()")
2768  call cpu_clock_end(id_clock_init)
2769 

◆ 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 3486 of file MOM.F90.

3486  type(MOM_control_struct), pointer :: CS !< MOM control structure
3487 
3488  if (cs%use_ALE_algorithm) call ale_end(cs%ALE_CSp)
3489 
3490  dealloc_(cs%u) ; dealloc_(cs%v) ; dealloc_(cs%h)
3491  dealloc_(cs%uh) ; dealloc_(cs%vh)
3492 
3493  if (associated(cs%tv%T)) then
3494  dealloc_(cs%T) ; cs%tv%T => null() ; dealloc_(cs%S) ; cs%tv%S => null()
3495  endif
3496  if (associated(cs%tv%frazil)) deallocate(cs%tv%frazil)
3497  if (associated(cs%tv%salt_deficit)) deallocate(cs%tv%salt_deficit)
3498  if (associated(cs%Hml)) deallocate(cs%Hml)
3499 
3500  call tracer_advect_end(cs%tracer_adv_CSp)
3501  call tracer_hor_diff_end(cs%tracer_diff_CSp)
3502  call tracer_registry_end(cs%tracer_Reg)
3503  call tracer_flow_control_end(cs%tracer_flow_CSp)
3504 
3505  call diabatic_driver_end(cs%diabatic_CSp)
3506 
3507  if (cs%offline_tracer_mode) call offline_transport_end(cs%offline_CSp)
3508 
3509  dealloc_(cs%uhtr) ; dealloc_(cs%vhtr)
3510  if (cs%split) then
3511  call end_dyn_split_rk2(cs%dyn_split_RK2_CSp)
3512  elseif (cs%use_RK2) then
3513  call end_dyn_unsplit_rk2(cs%dyn_unsplit_RK2_CSp)
3514  else
3515  call end_dyn_unsplit(cs%dyn_unsplit_CSp)
3516  endif
3517  dealloc_(cs%ave_ssh_ibc) ; dealloc_(cs%ssh_rint)
3518  if (associated(cs%update_OBC_CSp)) call obc_register_end(cs%update_OBC_CSp)
3519 
3520  call verticalgridend(cs%GV)
3521  call unit_scaling_end(cs%US)
3522  call mom_grid_end(cs%G)
3523 
3524  deallocate(cs)
3525 

◆ 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 3429 of file MOM.F90.

3429  type(MOM_control_struct), pointer :: CS !< MOM control structure
3430  logical, optional, intent(in) :: adv_dyn !< If present and true, only check
3431  !! whether the advection is up-to-date with
3432  !! the dynamics.
3433  logical :: in_synch !< True if all phases of the update are synchronized.
3434 
3435  logical :: adv_only
3436 
3437  adv_only = .false. ; if (present(adv_dyn)) adv_only = adv_dyn
3438 
3439  if (adv_only) then
3440  in_synch = (cs%t_dyn_rel_adv == 0.0)
3441  else
3442  in_synch = ((cs%t_dyn_rel_adv == 0.0) .and. (cs%t_dyn_rel_thermo == 0.0))
3443  endif
3444 

◆ 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 2855 of file MOM.F90.

2855  type(MOM_control_struct), intent(in) :: CS !< control structure set up by initialize_MOM.
2856 
2857  id_clock_ocean = cpu_clock_id('Ocean', grain=clock_component)
2858  id_clock_dynamics = cpu_clock_id('Ocean dynamics', grain=clock_subcomponent)
2859  id_clock_thermo = cpu_clock_id('Ocean thermodynamics and tracers', grain=clock_subcomponent)
2860  id_clock_other = cpu_clock_id('Ocean Other', grain=clock_subcomponent)
2861  id_clock_tracer = cpu_clock_id('(Ocean tracer advection)', grain=clock_module_driver)
2862  if (.not.cs%adiabatic) then
2863  id_clock_diabatic = cpu_clock_id('(Ocean diabatic driver)', grain=clock_module_driver)
2864  else
2865  id_clock_adiabatic = cpu_clock_id('(Ocean adiabatic driver)', grain=clock_module_driver)
2866  endif
2867 
2868  id_clock_continuity = cpu_clock_id('(Ocean continuity equation *)', grain=clock_module)
2869  id_clock_bbl_visc = cpu_clock_id('(Ocean set BBL viscosity)', grain=clock_module)
2870  id_clock_pass = cpu_clock_id('(Ocean message passing *)', grain=clock_module)
2871  id_clock_mom_init = cpu_clock_id('(Ocean MOM_initialize_state)', grain=clock_module)
2872  id_clock_pass_init = cpu_clock_id('(Ocean init message passing *)', grain=clock_routine)
2873  if (cs%thickness_diffuse) &
2874  id_clock_thick_diff = cpu_clock_id('(Ocean thickness diffusion *)', grain=clock_module)
2875 !if (CS%mixedlayer_restrat) &
2876  id_clock_ml_restrat = cpu_clock_id('(Ocean mixed layer restrat)', grain=clock_module)
2877  id_clock_diagnostics = cpu_clock_id('(Ocean collective diagnostics)', grain=clock_module)
2878  id_clock_z_diag = cpu_clock_id('(Ocean Z-space diagnostics)', grain=clock_module)
2879  id_clock_ale = cpu_clock_id('(Ocean ALE)', grain=clock_module)
2880  if (cs%offline_tracer_mode) then
2881  id_clock_offline_tracer = cpu_clock_id('Ocean offline tracers', grain=clock_subcomponent)
2882  endif
2883 

◆ 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 2824 of file MOM.F90.

2824  type(time_type), intent(in) :: Time !< current model time
2825  type(ocean_grid_type), intent(in) :: G !< ocean grid structure
2826  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2827  type(unit_scale_type), intent(inout) :: US !< A dimensional unit scaling type
2828  type(MOM_diag_IDs), intent(inout) :: IDs !< A structure with the diagnostic IDs.
2829  type(diag_ctrl), intent(inout) :: diag !< regulates diagnostic output
2830 
2831  real :: H_convert
2832  character(len=48) :: thickness_units
2833 
2834  thickness_units = get_thickness_units(gv)
2835  if (gv%Boussinesq) then
2836  h_convert = gv%H_to_m
2837  else
2838  h_convert = gv%H_to_kg_m2
2839  endif
2840 
2841  ! Diagnostics of the rapidly varying dynamic state
2842  ids%id_u = register_diag_field('ocean_model', 'u_dyn', diag%axesCuL, time, &
2843  'Zonal velocity after the dynamics update', 'm s-1', conversion=us%L_T_to_m_s)
2844  ids%id_v = register_diag_field('ocean_model', 'v_dyn', diag%axesCvL, time, &
2845  'Meridional velocity after the dynamics update', 'm s-1', conversion=us%L_T_to_m_s)
2846  ids%id_h = register_diag_field('ocean_model', 'h_dyn', diag%axesTL, time, &
2847  'Layer Thickness after the dynamics update', thickness_units, &
2848  v_extensive=.true., conversion=h_convert)
2849  ids%id_ssh_inst = register_diag_field('ocean_model', 'SSH_inst', diag%axesT1, &
2850  time, 'Instantaneous Sea Surface Height', 'm')

◆ 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 initialize_MOM
restart_csppointer to the restart control structure that will be used for MOM.

Definition at line 2896 of file MOM.F90.

2896  type(verticalGrid_type), intent(inout) :: GV !< ocean vertical grid structure
2897  type(unit_scale_type), intent(inout) :: US !< A dimensional unit scaling type
2898  type(param_file_type), intent(in) :: param_file !< opened file for parsing to get parameters
2899  type(MOM_control_struct), intent(in) :: CS !< control structure set up by initialize_MOM
2900  type(MOM_restart_CS), pointer :: restart_CSp !< pointer to the restart control
2901  !! structure that will be used for MOM.
2902  ! Local variables
2903  logical :: use_ice_shelf ! Needed to determine whether to add CS%Hml to restarts
2904  character(len=48) :: thickness_units, flux_units
2905  type(vardesc) :: u_desc, v_desc
2906 
2907  thickness_units = get_thickness_units(gv)
2908  flux_units = get_flux_units(gv)
2909 
2910  u_desc = var_desc("u", "m s-1", "Zonal velocity", hor_grid='Cu')
2911  v_desc = var_desc("v", "m s-1", "Meridional velocity", hor_grid='Cv')
2912 
2913  if (associated(cs%tv%T)) &
2914  call register_restart_field(cs%tv%T, "Temp", .true., restart_csp, &
2915  "Potential Temperature", "degC")
2916  if (associated(cs%tv%S)) &
2917  call register_restart_field(cs%tv%S, "Salt", .true., restart_csp, &
2918  "Salinity", "PPT")
2919 
2920  call register_restart_field(cs%h, "h", .true., restart_csp, &
2921  "Layer Thickness", thickness_units)
2922 
2923  call register_restart_pair(cs%u, cs%v, u_desc, v_desc, .true., restart_csp)
2924 
2925  if (associated(cs%tv%frazil)) &
2926  call register_restart_field(cs%tv%frazil, "frazil", .false., restart_csp, &
2927  "Frazil heat flux into ocean", "J m-2")
2928 
2929  if (cs%interp_p_surf) then
2930  call register_restart_field(cs%p_surf_prev, "p_surf_prev", .false., restart_csp, &
2931  "Previous ocean surface pressure", "Pa")
2932  endif
2933 
2934  call register_restart_field(cs%ave_ssh_ibc, "ave_ssh", .false., restart_csp, &
2935  "Time average sea surface height", "meter")
2936 
2937  ! hML is needed when using the ice shelf module
2938  call get_param(param_file, '', "ICE_SHELF", use_ice_shelf, default=.false., &
2939  do_not_log=.true.)
2940  if (use_ice_shelf .and. associated(cs%Hml)) then
2941  call register_restart_field(cs%Hml, "hML", .false., restart_csp, &
2942  "Mixed layer thickness", "meter")
2943  endif
2944 
2945  ! Register scalar unit conversion factors.
2946  call register_restart_field(us%m_to_Z_restart, "m_to_Z", .false., restart_csp, &
2947  "Height unit conversion factor", "Z meter-1")
2948  call register_restart_field(gv%m_to_H_restart, "m_to_H", .false., restart_csp, &
2949  "Thickness unit conversion factor", "H meter-1")
2950  call register_restart_field(us%m_to_L_restart, "m_to_L", .false., restart_csp, &
2951  "Length unit conversion factor", "L meter-1")
2952  call register_restart_field(us%s_to_T_restart, "s_to_T", .false., restart_csp, &
2953  "Time unit conversion factor", "T second-1")
2954  call register_restart_field(us%kg_m3_to_R_restart, "kg_m3_to_R", .false., restart_csp, &
2955  "Density unit conversion factor", "R m3 kg-1")
2956  call register_restart_field(us%J_kg_to_Q_restart, "J_kg_to_Q", .false., restart_csp, &
2957  "Heat content unit conversion factor.", units="Q kg J-1")
2958 

◆ step_mom()

subroutine, public mom::step_mom ( type(mech_forcing), intent(inout), target  forces_in,
type(forcing), intent(inout), target  fluxes_in,
type(surface), intent(inout), target  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]forces_inA structure with the driving mechanical forces
[in,out]fluxes_inA 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 422 of file MOM.F90.

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

◆ 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 [R L2 T-2 ~> Pa].
p_surf_endA pointer (perhaps NULL) to the surface pressure at the end of this dynamic step, intent in [R L2 T-2 ~> 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 946 of file MOM.F90.

946  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
947  real, dimension(:,:), pointer :: p_surf_begin !< A pointer (perhaps NULL) to the surface
948  !! pressure at the beginning of this dynamic
949  !! step, intent in [R L2 T-2 ~> Pa].
950  real, dimension(:,:), pointer :: p_surf_end !< A pointer (perhaps NULL) to the surface
951  !! pressure at the end of this dynamic step,
952  !! intent in [R L2 T-2 ~> Pa].
953  real, intent(in) :: dt !< time interval covered by this call [T ~> s].
954  real, intent(in) :: dt_thermo !< time interval covered by any updates that may
955  !! span multiple dynamics steps [T ~> s].
956  real, intent(in) :: bbl_time_int !< time interval over which updates to the
957  !! bottom boundary layer properties will apply [T ~> s],
958  !! or zero not to update the properties.
959  type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM
960  type(time_type), intent(in) :: Time_local !< End time of a segment, as a time type
961  type(wave_parameters_CS), &
962  optional, pointer :: Waves !< Container for wave related parameters; the
963  !! fields in Waves are intent in here.
964 
965  ! local variables
966  type(ocean_grid_type), pointer :: G => null() ! pointer to a structure containing
967  ! metrics and related information
968  type(verticalGrid_type), pointer :: GV => null() ! Pointer to the vertical grid structure
969  type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing
970  ! various unit conversion factors
971  type(MOM_diag_IDs), pointer :: IDs => null() ! A structure with the diagnostic IDs.
972  real, dimension(:,:,:), pointer :: &
973  u => null(), & ! u : zonal velocity component [L T-1 ~> m s-1]
974  v => null(), & ! v : meridional velocity component [L T-1 ~> m s-1]
975  h => null() ! h : layer thickness [H ~> m or kg m-2]
976 
977  logical :: calc_dtbt ! Indicates whether the dynamically adjusted
978  ! barotropic time step needs to be updated.
979  logical :: showCallTree
980 
981  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
982  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
983 
984  g => cs%G ; gv => cs%GV ; us => cs%US ; ids => cs%IDs
985  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
986  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
987  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
988  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
989  u => cs%u ; v => cs%v ; h => cs%h
990  showcalltree = calltree_showquery()
991 
992  call cpu_clock_begin(id_clock_dynamics)
993 
994  if ((cs%t_dyn_rel_adv == 0.0) .and. cs%thickness_diffuse .and. cs%thickness_diffuse_first) then
995 
996  call enable_averages(dt_thermo, time_local+real_to_time(us%T_to_s*(dt_thermo-dt)), cs%diag)
997  call cpu_clock_begin(id_clock_thick_diff)
998  if (associated(cs%VarMix)) &
999  call calc_slope_functions(h, cs%tv, dt, g, gv, us, cs%VarMix, obc=cs%OBC)
1000  call thickness_diffuse(h, cs%uhtr, cs%vhtr, cs%tv, dt_thermo, g, gv, us, &
1001  cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp)
1002  call cpu_clock_end(id_clock_thick_diff)
1003  call pass_var(h, g%Domain, clock=id_clock_pass, halo=max(2,cs%cont_stencil))
1004  call disable_averaging(cs%diag)
1005  if (showcalltree) call calltree_waypoint("finished thickness_diffuse_first (step_MOM)")
1006 
1007  ! Whenever thickness changes let the diag manager know, target grids
1008  ! for vertical remapping may need to be regenerated.
1009  call diag_update_remap_grids(cs%diag)
1010  endif
1011 
1012  ! The bottom boundary layer properties need to be recalculated.
1013  if (bbl_time_int > 0.0) then
1014  call enable_averages(bbl_time_int, &
1015  time_local + real_to_time(us%T_to_s*(bbl_time_int-dt)), cs%diag)
1016  ! Calculate the BBL properties and store them inside visc (u,h).
1017  call cpu_clock_begin(id_clock_bbl_visc)
1018  call set_viscous_bbl(cs%u, cs%v, cs%h, cs%tv, cs%visc, g, gv, us, &
1019  cs%set_visc_CSp, symmetrize=.true.)
1020  call cpu_clock_end(id_clock_bbl_visc)
1021  if (showcalltree) call calltree_waypoint("done with set_viscous_BBL (step_MOM)")
1022  call disable_averaging(cs%diag)
1023  endif
1024 
1025 
1026  if (cs%do_dynamics .and. cs%split) then !--------------------------- start SPLIT
1027  ! This section uses a split time stepping scheme for the dynamic equations,
1028  ! basically the stacked shallow water equations with viscosity.
1029 
1030  calc_dtbt = .false.
1031  if (cs%dtbt_reset_period == 0.0) calc_dtbt = .true.
1032  if (cs%dtbt_reset_period > 0.0) then
1033  if (time_local >= cs%dtbt_reset_time) then !### Change >= to > here.
1034  calc_dtbt = .true.
1035  cs%dtbt_reset_time = cs%dtbt_reset_time + cs%dtbt_reset_interval
1036  endif
1037  endif
1038 
1039  call step_mom_dyn_split_rk2(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
1040  p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
1041  cs%eta_av_bc, g, gv, us, cs%dyn_split_RK2_CSp, calc_dtbt, cs%VarMix, &
1042  cs%MEKE, cs%thickness_diffuse_CSp, waves=waves)
1043  if (showcalltree) call calltree_waypoint("finished step_MOM_dyn_split (step_MOM)")
1044 
1045  elseif (cs%do_dynamics) then ! ------------------------------------ not SPLIT
1046  ! This section uses an unsplit stepping scheme for the dynamic
1047  ! equations; basically the stacked shallow water equations with viscosity.
1048  ! Because the time step is limited by CFL restrictions on the external
1049  ! gravity waves, the unsplit is usually much less efficient that the split
1050  ! approaches. But because of its simplicity, the unsplit method is very
1051  ! useful for debugging purposes.
1052 
1053  if (cs%use_RK2) then
1054  call step_mom_dyn_unsplit_rk2(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
1055  p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
1056  cs%eta_av_bc, g, gv, us, cs%dyn_unsplit_RK2_CSp, cs%VarMix, cs%MEKE)
1057  else
1058  call step_mom_dyn_unsplit(u, v, h, cs%tv, cs%visc, time_local, dt, forces, &
1059  p_surf_begin, p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
1060  cs%eta_av_bc, g, gv, us, cs%dyn_unsplit_CSp, cs%VarMix, cs%MEKE, waves=waves)
1061  endif
1062  if (showcalltree) call calltree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)")
1063 
1064  endif ! -------------------------------------------------- end SPLIT
1065 
1066  if (cs%thickness_diffuse .and. .not.cs%thickness_diffuse_first) then
1067  call cpu_clock_begin(id_clock_thick_diff)
1068 
1069  if (cs%debug) call hchksum(h,"Pre-thickness_diffuse h", g%HI, haloshift=0, scale=gv%H_to_m)
1070 
1071  if (associated(cs%VarMix)) &
1072  call calc_slope_functions(h, cs%tv, dt, g, gv, us, cs%VarMix, obc=cs%OBC)
1073  call thickness_diffuse(h, cs%uhtr, cs%vhtr, cs%tv, dt, g, gv, us, &
1074  cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp)
1075 
1076  if (cs%debug) call hchksum(h,"Post-thickness_diffuse h", g%HI, haloshift=1, scale=gv%H_to_m)
1077  call cpu_clock_end(id_clock_thick_diff)
1078  call pass_var(h, g%Domain, clock=id_clock_pass, halo=max(2,cs%cont_stencil))
1079  if (showcalltree) call calltree_waypoint("finished thickness_diffuse (step_MOM)")
1080  endif
1081 
1082  ! apply the submesoscale mixed layer restratification parameterization
1083  if (cs%mixedlayer_restrat) then
1084  if (cs%debug) then
1085  call hchksum(h,"Pre-mixedlayer_restrat h", g%HI, haloshift=1, scale=gv%H_to_m)
1086  call uvchksum("Pre-mixedlayer_restrat uhtr", &
1087  cs%uhtr, cs%vhtr, g%HI, haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1088  endif
1089  call cpu_clock_begin(id_clock_ml_restrat)
1090  call mixedlayer_restrat(h, cs%uhtr, cs%vhtr, cs%tv, forces, dt, cs%visc%MLD, &
1091  cs%VarMix, g, gv, us, cs%mixedlayer_restrat_CSp)
1092  call cpu_clock_end(id_clock_ml_restrat)
1093  call pass_var(h, g%Domain, clock=id_clock_pass, halo=max(2,cs%cont_stencil))
1094  if (cs%debug) then
1095  call hchksum(h,"Post-mixedlayer_restrat h", g%HI, haloshift=1, scale=gv%H_to_m)
1096  call uvchksum("Post-mixedlayer_restrat [uv]htr", &
1097  cs%uhtr, cs%vhtr, g%HI, haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1098  endif
1099  endif
1100 
1101  ! Whenever thickness changes let the diag manager know, target grids
1102  ! for vertical remapping may need to be regenerated.
1103  call diag_update_remap_grids(cs%diag)
1104 
1105  if (cs%useMEKE) call step_forward_meke(cs%MEKE, h, cs%VarMix%SN_u, cs%VarMix%SN_v, &
1106  cs%visc, dt, g, gv, us, cs%MEKE_CSp, cs%uhtr, cs%vhtr)
1107  call disable_averaging(cs%diag)
1108 
1109  ! Advance the dynamics time by dt.
1110  cs%t_dyn_rel_adv = cs%t_dyn_rel_adv + dt
1111  cs%t_dyn_rel_thermo = cs%t_dyn_rel_thermo + dt
1112  if (abs(cs%t_dyn_rel_thermo) < 1e-6*dt) cs%t_dyn_rel_thermo = 0.0
1113  cs%t_dyn_rel_diag = cs%t_dyn_rel_diag + dt
1114 
1115  call cpu_clock_end(id_clock_dynamics)
1116 
1117  call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
1118  call enable_averages(dt, time_local, cs%diag)
1119  ! These diagnostics are available after every time dynamics step.
1120  if (ids%id_u > 0) call post_data(ids%id_u, u, cs%diag)
1121  if (ids%id_v > 0) call post_data(ids%id_v, v, cs%diag)
1122  if (ids%id_h > 0) call post_data(ids%id_h, h, cs%diag)
1123  call disable_averaging(cs%diag)
1124  call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
1125 

◆ 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 1213 of file MOM.F90.

1213  type(MOM_control_struct), intent(inout) :: CS !< Master MOM control structure
1214  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1215  type(verticalGrid_type), intent(inout) :: GV !< ocean vertical grid structure
1216  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1217  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1218  intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1219  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1220  intent(inout) :: v !< meridional velocity [L T-1 ~> m s-1]
1221  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1222  intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
1223  type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
1224  type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
1225  real, intent(in) :: dtdia !< The time interval over which to advance [T ~> s]
1226  type(time_type), intent(in) :: Time_end_thermo !< End of averaging interval for thermo diags
1227  logical, intent(in) :: update_BBL !< If true, calculate the bottom boundary layer properties.
1228  type(wave_parameters_CS), &
1229  optional, pointer :: Waves !< Container for wave related parameters
1230  !! the fields in Waves are intent in here.
1231 
1232  logical :: use_ice_shelf ! Needed for selecting the right ALE interface.
1233  logical :: showCallTree
1234  type(group_pass_type) :: pass_T_S, pass_T_S_h, pass_uv_T_S_h
1235  integer :: dynamics_stencil ! The computational stencil for the calculations
1236  ! in the dynamic core.
1237  integer :: halo_sz ! The size of a halo where data must be valid.
1238  integer :: i, j, k, is, ie, js, je, nz
1239 
1240  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1241  showcalltree = calltree_showquery()
1242  if (showcalltree) call calltree_enter("step_MOM_thermo(), MOM.F90")
1243 
1244  use_ice_shelf = .false.
1245  if (associated(fluxes%frac_shelf_h)) use_ice_shelf = .true.
1246 
1247  call enable_averages(dtdia, time_end_thermo, cs%diag)
1248 
1249  if (associated(cs%odaCS)) then
1250  call apply_oda_tracer_increments(us%T_to_s*dtdia,g,tv,h,cs%odaCS)
1251  endif
1252 
1253  if (associated(fluxes%p_surf) .or. associated(fluxes%p_surf_full)) then
1254  call extract_diabatic_member(cs%diabatic_CSp, diabatic_halo=halo_sz)
1255  if (halo_sz > 0) then
1256  if (associated(fluxes%p_surf_full)) &
1257  call pass_var(fluxes%p_surf_full, g%Domain, &
1258  clock=id_clock_pass, halo=halo_sz, complete=.not.associated(fluxes%p_surf))
1259  call pass_var(fluxes%p_surf, g%Domain, clock=id_clock_pass, halo=halo_sz, complete=.true.)
1260  endif
1261  endif
1262 
1263  if (update_bbl) then
1264  ! Calculate the BBL properties and store them inside visc (u,h).
1265  ! This is here so that CS%visc is updated before diabatic() when
1266  ! DIABATIC_FIRST=True. Otherwise diabatic() is called after the dynamics
1267  ! and set_viscous_BBL is called as a part of the dynamic stepping.
1268  call cpu_clock_begin(id_clock_bbl_visc)
1269  call set_viscous_bbl(u, v, h, tv, cs%visc, g, gv, us, cs%set_visc_CSp, symmetrize=.true.)
1270  call cpu_clock_end(id_clock_bbl_visc)
1271  if (showcalltree) call calltree_waypoint("done with set_viscous_BBL (step_MOM_thermo)")
1272  endif
1273 
1274  call cpu_clock_begin(id_clock_thermo)
1275  if (.not.cs%adiabatic) then
1276  if (cs%debug) then
1277  call uvchksum("Pre-diabatic [uv]", u, v, g%HI, haloshift=2, scale=us%L_T_to_m_s)
1278  call hchksum(h,"Pre-diabatic h", g%HI, haloshift=1, scale=gv%H_to_m)
1279  call uvchksum("Pre-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1280  haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1281  ! call MOM_state_chksum("Pre-diabatic ", u, v, h, CS%uhtr, CS%vhtr, G, GV, vel_scale=1.0)
1282  call mom_thermo_chksum("Pre-diabatic ", tv, g, us, haloshift=0)
1283  call check_redundant("Pre-diabatic ", u, v, g)
1284  call mom_forcing_chksum("Pre-diabatic", fluxes, g, us, haloshift=0)
1285  endif
1286 
1287  call cpu_clock_begin(id_clock_diabatic)
1288 
1289  call diabatic(u, v, h, tv, cs%Hml, fluxes, cs%visc, cs%ADp, cs%CDp, dtdia, &
1290  time_end_thermo, g, gv, us, cs%diabatic_CSp, obc=cs%OBC, waves=waves)
1291  fluxes%fluxes_used = .true.
1292 
1293  if (showcalltree) call calltree_waypoint("finished diabatic (step_MOM_thermo)")
1294 
1295  ! Regridding/remapping is done here, at end of thermodynamics time step
1296  ! (that may comprise several dynamical time steps)
1297  ! The routine 'ALE_main' can be found in 'MOM_ALE.F90'.
1298  if ( cs%use_ALE_algorithm ) then
1299  call enable_averages(dtdia, time_end_thermo, cs%diag)
1300 ! call pass_vector(u, v, G%Domain)
1301  call cpu_clock_begin(id_clock_pass)
1302  if (associated(tv%T)) &
1303  call create_group_pass(pass_t_s_h, tv%T, g%Domain, to_all+omit_corners, halo=1)
1304  if (associated(tv%S)) &
1305  call create_group_pass(pass_t_s_h, tv%S, g%Domain, to_all+omit_corners, halo=1)
1306  call create_group_pass(pass_t_s_h, h, g%Domain, to_all+omit_corners, halo=1)
1307  call do_group_pass(pass_t_s_h, g%Domain)
1308  call cpu_clock_end(id_clock_pass)
1309 
1310  call preale_tracer_diagnostics(cs%tracer_Reg, g, gv)
1311 
1312  if (cs%debug) then
1313  call mom_state_chksum("Pre-ALE ", u, v, h, cs%uh, cs%vh, g, gv, us)
1314  call hchksum(tv%T,"Pre-ALE T", g%HI, haloshift=1)
1315  call hchksum(tv%S,"Pre-ALE S", g%HI, haloshift=1)
1316  call check_redundant("Pre-ALE ", u, v, g)
1317  endif
1318  call cpu_clock_begin(id_clock_ale)
1319  if (use_ice_shelf) then
1320  call ale_main(g, gv, us, h, u, v, tv, cs%tracer_Reg, cs%ALE_CSp, cs%OBC, &
1321  dtdia, fluxes%frac_shelf_h)
1322  else
1323  call ale_main(g, gv, us, h, u, v, tv, cs%tracer_Reg, cs%ALE_CSp, cs%OBC, dtdia)
1324  endif
1325 
1326  if (showcalltree) call calltree_waypoint("finished ALE_main (step_MOM_thermo)")
1327  call cpu_clock_end(id_clock_ale)
1328  endif ! endif for the block "if ( CS%use_ALE_algorithm )"
1329 
1330  dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
1331  call create_group_pass(pass_uv_t_s_h, u, v, g%Domain, halo=dynamics_stencil)
1332  if (associated(tv%T)) &
1333  call create_group_pass(pass_uv_t_s_h, tv%T, g%Domain, halo=dynamics_stencil)
1334  if (associated(tv%S)) &
1335  call create_group_pass(pass_uv_t_s_h, tv%S, g%Domain, halo=dynamics_stencil)
1336  call create_group_pass(pass_uv_t_s_h, h, g%Domain, halo=dynamics_stencil)
1337  call do_group_pass(pass_uv_t_s_h, g%Domain, clock=id_clock_pass)
1338 
1339  if (cs%debug .and. cs%use_ALE_algorithm) then
1340  call mom_state_chksum("Post-ALE ", u, v, h, cs%uh, cs%vh, g, gv, us)
1341  call hchksum(tv%T, "Post-ALE T", g%HI, haloshift=1)
1342  call hchksum(tv%S, "Post-ALE S", g%HI, haloshift=1)
1343  call check_redundant("Post-ALE ", u, v, g)
1344  endif
1345 
1346  ! Whenever thickness changes let the diag manager know, target grids
1347  ! for vertical remapping may need to be regenerated. This needs to
1348  ! happen after the H update and before the next post_data.
1349  call diag_update_remap_grids(cs%diag)
1350 
1351  !### Consider moving this up into the if ALE block.
1352  call postale_tracer_diagnostics(cs%tracer_Reg, g, gv, cs%diag, dtdia)
1353 
1354  if (cs%debug) then
1355  call uvchksum("Post-diabatic u", u, v, g%HI, haloshift=2, scale=us%L_T_to_m_s)
1356  call hchksum(h, "Post-diabatic h", g%HI, haloshift=1, scale=gv%H_to_m)
1357  call uvchksum("Post-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1358  haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1359  ! call MOM_state_chksum("Post-diabatic ", u, v, &
1360  ! h, CS%uhtr, CS%vhtr, G, GV, haloshift=1)
1361  if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", g%HI, haloshift=1)
1362  if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", g%HI, haloshift=1)
1363  if (associated(tv%frazil)) call hchksum(tv%frazil, "Post-diabatic frazil", g%HI, haloshift=0, &
1364  scale=g%US%Q_to_J_kg*g%US%RZ_to_kg_m2)
1365  if (associated(tv%salt_deficit)) call hchksum(tv%salt_deficit, &
1366  "Post-diabatic salt deficit", g%HI, haloshift=0, scale=us%RZ_to_kg_m2)
1367  ! call MOM_thermo_chksum("Post-diabatic ", tv, G, US)
1368  call check_redundant("Post-diabatic ", u, v, g)
1369  endif
1370  call disable_averaging(cs%diag)
1371 
1372  call cpu_clock_end(id_clock_diabatic)
1373  else ! complement of "if (.not.CS%adiabatic)"
1374 
1375  call cpu_clock_begin(id_clock_adiabatic)
1376  call adiabatic(h, tv, fluxes, dtdia, g, gv, us, cs%diabatic_CSp)
1377  fluxes%fluxes_used = .true.
1378  call cpu_clock_end(id_clock_adiabatic)
1379 
1380  if (associated(tv%T)) then
1381  call create_group_pass(pass_t_s, tv%T, g%Domain, to_all+omit_corners, halo=1)
1382  call create_group_pass(pass_t_s, tv%S, g%Domain, to_all+omit_corners, halo=1)
1383  call do_group_pass(pass_t_s, g%Domain, clock=id_clock_pass)
1384  if (cs%debug) then
1385  if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", g%HI, haloshift=1)
1386  if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", g%HI, haloshift=1)
1387  endif
1388  endif
1389 
1390  endif ! endif for the block "if (.not.CS%adiabatic)"
1391  call cpu_clock_end(id_clock_thermo)
1392 
1393  call disable_averaging(cs%diag)
1394 
1395  if (showcalltree) call calltree_leave("step_MOM_thermo(), MOM.F90")
1396 

◆ 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 1132 of file MOM.F90.

1132  type(MOM_control_struct), intent(inout) :: CS !< control structure
1133  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1134  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1135  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1136  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1137  intent(in) :: h !< layer thicknesses after the transports [H ~> m or kg m-2]
1138  type(time_type), intent(in) :: Time_local !< The model time at the end
1139  !! of the time step.
1140  type(group_pass_type) :: pass_T_S
1141  integer :: halo_sz ! The size of a halo where data must be valid.
1142  logical :: showCallTree
1143  showcalltree = calltree_showquery()
1144 
1145  if (cs%debug) then
1146  call cpu_clock_begin(id_clock_other)
1147  call hchksum(h,"Pre-advection h", g%HI, haloshift=1, scale=gv%H_to_m)
1148  call uvchksum("Pre-advection uhtr", cs%uhtr, cs%vhtr, g%HI, &
1149  haloshift=0, scale=gv%H_to_m*us%L_to_m**2)
1150  if (associated(cs%tv%T)) call hchksum(cs%tv%T, "Pre-advection T", g%HI, haloshift=1)
1151  if (associated(cs%tv%S)) call hchksum(cs%tv%S, "Pre-advection S", g%HI, haloshift=1)
1152  if (associated(cs%tv%frazil)) call hchksum(cs%tv%frazil, "Pre-advection frazil", g%HI, haloshift=0, &
1153  scale=g%US%Q_to_J_kg*g%US%RZ_to_kg_m2)
1154  if (associated(cs%tv%salt_deficit)) call hchksum(cs%tv%salt_deficit, &
1155  "Pre-advection salt deficit", g%HI, haloshift=0, scale=us%RZ_to_kg_m2)
1156  ! call MOM_thermo_chksum("Pre-advection ", CS%tv, G, US)
1157  call cpu_clock_end(id_clock_other)
1158  endif
1159 
1160  call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer)
1161  call enable_averages(cs%t_dyn_rel_adv, time_local, cs%diag)
1162 
1163  call advect_tracer(h, cs%uhtr, cs%vhtr, cs%OBC, cs%t_dyn_rel_adv, g, gv, us, &
1164  cs%tracer_adv_CSp, cs%tracer_Reg)
1165  call tracer_hordiff(h, cs%t_dyn_rel_adv, cs%MEKE, cs%VarMix, g, gv, us, &
1166  cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1167  if (showcalltree) call calltree_waypoint("finished tracer advection/diffusion (step_MOM)")
1168  call update_segment_tracer_reservoirs(g, gv, cs%uhtr, cs%vhtr, h, cs%OBC, &
1169  cs%t_dyn_rel_adv, cs%tracer_Reg)
1170  call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo)
1171 
1172  call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
1173  call post_transport_diagnostics(g, gv, us, cs%uhtr, cs%vhtr, h, cs%transport_IDs, &
1174  cs%diag_pre_dyn, cs%diag, cs%t_dyn_rel_adv, cs%tracer_reg)
1175  ! Rebuild the remap grids now that we've posted the fields which rely on thicknesses
1176  ! from before the dynamics calls
1177  call diag_update_remap_grids(cs%diag)
1178 
1179  call disable_averaging(cs%diag)
1180  call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
1181 
1182  ! Reset the accumulated transports to 0 and record that the dynamics
1183  ! and advective times now agree.
1184  call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer)
1185  cs%uhtr(:,:,:) = 0.0
1186  cs%vhtr(:,:,:) = 0.0
1187  cs%t_dyn_rel_adv = 0.0
1188  call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo)
1189 
1190  if (associated(cs%tv%T)) then
1191  call extract_diabatic_member(cs%diabatic_CSp, diabatic_halo=halo_sz)
1192  if (halo_sz > 0) then
1193  call create_group_pass(pass_t_s, cs%tv%T, g%Domain, to_all, halo=halo_sz)
1194  call create_group_pass(pass_t_s, cs%tv%S, g%Domain, to_all, halo=halo_sz)
1195  call do_group_pass(pass_t_s, g%Domain, clock=id_clock_pass)
1196  elseif (cs%diabatic_first) then
1197  ! Temperature and salinity need halo updates because they will be used
1198  ! in the dynamics before they are changed again.
1199  call create_group_pass(pass_t_s, cs%tv%T, g%Domain, to_all+omit_corners, halo=1)
1200  call create_group_pass(pass_t_s, cs%tv%S, g%Domain, to_all+omit_corners, halo=1)
1201  call do_group_pass(pass_t_s, g%Domain, clock=id_clock_pass)
1202  endif
1203  endif
1204 
1205  cs%preadv_h_stored = .false.
1206 

◆ 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 1405 of file MOM.F90.

1405  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1406  type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
1407  type(surface), intent(inout) :: sfc_state !< surface ocean state
1408  type(time_type), intent(in) :: Time_start !< starting time of a segment, as a time type
1409  real, intent(in) :: time_interval !< time interval
1410  type(MOM_control_struct), pointer :: CS !< control structure from initialize_MOM
1411 
1412  ! Local pointers
1413  type(ocean_grid_type), pointer :: G => null() ! Pointer to a structure containing
1414  ! metrics and related information
1415  type(verticalGrid_type), pointer :: GV => null() ! Pointer to structure containing information
1416  ! about the vertical grid
1417  type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing
1418  ! various unit conversion factors
1419 
1420  logical :: first_iter !< True if this is the first time step_offline has been called in a given interval
1421  logical :: last_iter !< True if this is the last time step_tracer is to be called in an offline interval
1422  logical :: do_vertical !< If enough time has elapsed, do the diabatic tracer sources/sinks
1423  logical :: adv_converged !< True if all the horizontal fluxes have been used
1424 
1425  real :: dt_offline ! The offline timestep for advection [T ~> s]
1426  real :: dt_offline_vertical ! The offline timestep for vertical fluxes and remapping [T ~> s]
1427  logical :: skip_diffusion
1428  integer :: id_eta_diff_end
1429 
1430  type(time_type), pointer :: accumulated_time => null()
1431  type(time_type), pointer :: vertical_time => null()
1432  integer :: i,j,k
1433  integer :: is, ie, js, je, isd, ied, jsd, jed
1434 
1435  ! 3D pointers
1436  real, dimension(:,:,:), pointer :: &
1437  uhtr => null(), vhtr => null(), &
1438  eatr => null(), ebtr => null(), &
1439  h_end => null()
1440 
1441  ! 2D Array for diagnostics
1442  real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: eta_pre, eta_end
1443  type(time_type) :: Time_end ! End time of a segment, as a time type
1444 
1445  ! Grid-related pointer assignments
1446  g => cs%G ; gv => cs%GV ; us => cs%US
1447 
1448  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1449  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1450 
1451  call cpu_clock_begin(id_clock_offline_tracer)
1452  call extract_offline_main(cs%offline_CSp, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, &
1453  vertical_time, dt_offline, dt_offline_vertical, skip_diffusion)
1454  time_end = increment_date(time_start, seconds=floor(time_interval+0.001))
1455 
1456  call enable_averaging(time_interval, time_end, cs%diag)
1457 
1458  ! Check to see if this is the first iteration of the offline interval
1459  if (accumulated_time == real_to_time(0.0)) then
1460  first_iter = .true.
1461  else ! This is probably unnecessary but is used to guard against unwanted behavior
1462  first_iter = .false.
1463  endif
1464 
1465  ! Check to see if vertical tracer functions should be done
1466  if (first_iter .or. (accumulated_time >= vertical_time)) then
1467  do_vertical = .true.
1468  vertical_time = accumulated_time + real_to_time(us%T_to_s*dt_offline_vertical)
1469  else
1470  do_vertical = .false.
1471  endif
1472 
1473  ! Increment the amount of time elapsed since last read and check if it's time to roll around
1474  accumulated_time = accumulated_time + real_to_time(time_interval)
1475 
1476  last_iter = (accumulated_time >= real_to_time(us%T_to_s*dt_offline))
1477 
1478  if (cs%use_ALE_algorithm) then
1479  ! If this is the first iteration in the offline timestep, then we need to read in fields and
1480  ! perform the main advection.
1481  if (first_iter) then
1482  call mom_mesg("Reading in new offline fields")
1483  ! Read in new transport and other fields
1484  ! call update_transport_from_files(G, GV, CS%offline_CSp, h_end, eatr, ebtr, uhtr, vhtr, &
1485  ! CS%tv%T, CS%tv%S, fluxes, CS%use_ALE_algorithm)
1486  ! call update_transport_from_arrays(CS%offline_CSp)
1487  call update_offline_fields(cs%offline_CSp, cs%h, fluxes, cs%use_ALE_algorithm)
1488 
1489  ! Apply any fluxes into the ocean
1490  call offline_fw_fluxes_into_ocean(g, gv, cs%offline_CSp, fluxes, cs%h)
1491 
1492  if (.not.cs%diabatic_first) then
1493  call offline_advection_ale(fluxes, time_start, time_interval, cs%offline_CSp, id_clock_ale, &
1494  cs%h, uhtr, vhtr, converged=adv_converged)
1495 
1496  ! Redistribute any remaining transport
1497  call offline_redistribute_residual(cs%offline_CSp, cs%h, uhtr, vhtr, adv_converged)
1498 
1499  ! Perform offline diffusion if requested
1500  if (.not. skip_diffusion) then
1501  if (associated(cs%VarMix)) then
1502  call pass_var(cs%h, g%Domain)
1503  call calc_resoln_function(cs%h, cs%tv, g, gv, us, cs%VarMix)
1504  call calc_depth_function(g, cs%VarMix)
1505  call calc_slope_functions(cs%h, cs%tv, dt_offline, g, gv, us, cs%VarMix, obc=cs%OBC)
1506  endif
1507  call tracer_hordiff(cs%h, dt_offline, cs%MEKE, cs%VarMix, g, gv, us, &
1508  cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1509  endif
1510  endif
1511  endif
1512  ! The functions related to column physics of tracers is performed separately in ALE mode
1513  if (do_vertical) then
1514  call offline_diabatic_ale(fluxes, time_start, time_end, cs%offline_CSp, cs%h, eatr, ebtr)
1515  endif
1516 
1517  ! Last thing that needs to be done is the final ALE remapping
1518  if (last_iter) then
1519  if (cs%diabatic_first) then
1520  call offline_advection_ale(fluxes, time_start, time_interval, cs%offline_CSp, id_clock_ale, &
1521  cs%h, uhtr, vhtr, converged=adv_converged)
1522 
1523  ! Redistribute any remaining transport and perform the remaining advection
1524  call offline_redistribute_residual(cs%offline_CSp, cs%h, uhtr, vhtr, adv_converged)
1525  ! Perform offline diffusion if requested
1526  if (.not. skip_diffusion) then
1527  if (associated(cs%VarMix)) then
1528  call pass_var(cs%h, g%Domain)
1529  call calc_resoln_function(cs%h, cs%tv, g, gv, us, cs%VarMix)
1530  call calc_depth_function(g, cs%VarMix)
1531  call calc_slope_functions(cs%h, cs%tv, dt_offline, g, gv, us, cs%VarMix, obc=cs%OBC)
1532  endif
1533  call tracer_hordiff(cs%h, dt_offline, cs%MEKE, cs%VarMix, g, gv, us, &
1534  cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1535  endif
1536  endif
1537 
1538  call mom_mesg("Last iteration of offline interval")
1539 
1540  ! Apply freshwater fluxes out of the ocean
1541  call offline_fw_fluxes_out_ocean(g, gv, cs%offline_CSp, fluxes, cs%h)
1542  ! These diagnostic can be used to identify which grid points did not converge within
1543  ! the specified number of advection sub iterations
1544  call post_offline_convergence_diags(cs%offline_CSp, cs%h, h_end, uhtr, vhtr)
1545 
1546  ! Call ALE one last time to make sure that tracers are remapped onto the layer thicknesses
1547  ! stored from the forward run
1548  call cpu_clock_begin(id_clock_ale)
1549  call ale_offline_tracer_final( g, gv, cs%h, cs%tv, h_end, cs%tracer_Reg, cs%ALE_CSp, cs%OBC)
1550  call cpu_clock_end(id_clock_ale)
1551  call pass_var(cs%h, g%Domain)
1552  endif
1553  else ! NON-ALE MODE...NOT WELL TESTED
1554  call mom_error(warning, &
1555  "Offline tracer mode in non-ALE configuration has not been thoroughly tested")
1556  ! Note that for the layer mode case, the calls to tracer sources and sinks is embedded in
1557  ! main_offline_advection_layer. Warning: this may not be appropriate for tracers that
1558  ! exchange with the atmosphere
1559  if (abs(time_interval - us%T_to_s*dt_offline) > 1.0e-6) then
1560  call mom_error(fatal, &
1561  "For offline tracer mode in a non-ALE configuration, dt_offline must equal time_interval")
1562  endif
1563  call update_offline_fields(cs%offline_CSp, cs%h, fluxes, cs%use_ALE_algorithm)
1564  call offline_advection_layer(fluxes, time_start, time_interval, cs%offline_CSp, &
1565  cs%h, eatr, ebtr, uhtr, vhtr)
1566  ! Perform offline diffusion if requested
1567  if (.not. skip_diffusion) then
1568  call tracer_hordiff(h_end, dt_offline, cs%MEKE, cs%VarMix, g, gv, us, &
1569  cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
1570  endif
1571 
1572  cs%h = h_end
1573 
1574  call pass_var(cs%tv%T, g%Domain)
1575  call pass_var(cs%tv%S, g%Domain)
1576  call pass_var(cs%h, g%Domain)
1577 
1578  endif
1579 
1580  call adjust_ssh_for_p_atm(cs%tv, g, gv, us, cs%ave_ssh_ibc, forces%p_surf_SSH, &
1581  cs%calc_rho_for_sea_lev)
1582  call extract_surface_state(cs, sfc_state)
1583 
1584  call disable_averaging(cs%diag)
1585  call pass_var(cs%tv%T, g%Domain)
1586  call pass_var(cs%tv%S, g%Domain)
1587  call pass_var(cs%h, g%Domain)
1588 
1589  fluxes%fluxes_used = .true.
1590 
1591  if (last_iter) then
1592  accumulated_time = real_to_time(0.0)
1593  endif
1594 
1595  call cpu_clock_end(id_clock_offline_tracer)
1596