MOM6
Data Types | Functions/Subroutines | Variables
mom Module Reference

Detailed Description

This is the main routine for MOM.

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 accomodated.
  • 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.

  • calculate_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 this module. More...
 

Functions/Subroutines

subroutine, public step_mom (fluxes, state, Time_start, time_interval, CS)
 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_thermo (CS, G, GV, u, v, h, tv, fluxes, dtdia)
 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 (fluxes, 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, param_file, dirs, CS, Time_in, offline_tracer_mode)
 This subroutine initializes MOM. More...
 
subroutine, public finish_mom_initialization (Time, dirs, CS, fluxes)
 This subroutine finishes initializing MOM and writes out the initial conditions. More...
 
subroutine register_diags (Time, G, GV, CS, ADp)
 Register the diagnostics. More...
 
subroutine register_diags_ts_tendency (Time, G, CS)
 Initialize diagnostics for temp/heat and saln/salt tendencies. More...
 
subroutine register_diags_ts_vardec (Time, HI, GV, param_file, CS)
 Initialize diagnostics for the variance decay of temp/salt across regridding/remapping. More...
 
subroutine mom_timing_init (CS)
 This subroutine sets up clock IDs for timing various subroutines. More...
 
subroutine post_transport_diagnostics (G, GV, CS, diag, dt_trans, h, h_pre_dyn)
 This routine posts diagnostics of the transports, including the subgridscale contributions. More...
 
subroutine post_ts_diagnostics (CS, G, GV, tv, diag, dt)
 Post diagnostics of temperatures and salinities, their fluxes, and tendencies. More...
 
subroutine post_diags_ts_vardec (G, CS, dt)
 Calculate and post variance decay diagnostics for temp/salt. More...
 
subroutine post_integrated_diagnostics (CS, G, GV, diag, dt_int, tv, fluxes)
 This routine posts diagnostics of various integrated quantities. More...
 
subroutine post_surface_diagnostics (CS, G, diag, state)
 This routine posts diagnostics of various ocean surface quantities. More...
 
subroutine write_static_fields (G, diag)
 Offers the static fields in the ocean grid type for output via the diag_manager. More...
 
subroutine set_restart_fields (GV, param_file, CS)
 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 (CS, G, GV, ssh, p_atm)
 This subroutine applies a correction to the sea surface height to compensate for the atmospheric pressure (the inverse barometer). More...
 
subroutine, public calculate_surface_state (state, u, v, h, ssh, G, GV, CS)
 This subroutine sets the surface (return) properties of the ocean model by setting the appropriate fields in state. Unused fields are set to NULL or are unallocated. More...
 
subroutine, public mom_end (CS)
 End of model. More...
 

Variables

integer id_clock_ocean
 
integer id_clock_dynamics
 
integer id_clock_thermo
 
integer id_clock_tracer
 
integer id_clock_diabatic
 
integer id_clock_continuity
 
integer id_clock_thick_diff
 
integer id_clock_bbl_visc
 
integer id_clock_ml_restrat
 
integer id_clock_diagnostics
 
integer id_clock_z_diag
 
integer id_clock_init
 
integer id_clock_mom_init
 
integer id_clock_pass
 
integer id_clock_pass_init
 
integer id_clock_ale
 
integer id_clock_other
 
integer id_clock_offline_tracer
 

Function/Subroutine Documentation

◆ adjust_ssh_for_p_atm()

subroutine mom::adjust_ssh_for_p_atm ( type(mom_control_struct), intent(in)  CS,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
real, dimension(szi_(g),szj_(g)), intent(inout)  ssh,
real, dimension(:,:), optional, pointer  p_atm 
)
private

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

Parameters
[in]cscontrol structure
[in]gocean grid structure
[in]gvocean vertical grid structure
[in,out]sshtime mean surface height (m)
p_atmatmospheric pressure (Pascal)

Definition at line 3408 of file MOM.F90.

Referenced by step_mom(), and step_offline().

3408  type(mom_control_struct), intent(in) :: cs !< control structure
3409  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
3410  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
3411  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: ssh !< time mean surface height (m)
3412  real, dimension(:,:), optional, pointer :: p_atm !< atmospheric pressure (Pascal)
3413 
3414  real :: rho_conv ! The density used to convert surface pressure to
3415  ! a corrected effective SSH, in kg m-3.
3416  real :: igr0 ! The SSH conversion factor from Pa to m.
3417  integer :: i, j, is, ie, js, je
3418 
3419  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3420  if (ASSOCIATED(p_atm)) then
3421  ! Correct the output sea surface height for the contribution from the
3422  ! atmospheric pressure
3423  do j=js,je ; do i=is,ie
3424  if ((ASSOCIATED(cs%tv%eqn_of_state)) .and. (cs%calc_rho_for_sea_lev)) then
3425  call calculate_density(cs%tv%T(i,j,1), cs%tv%S(i,j,1), p_atm(i,j)/2.0, &
3426  rho_conv, cs%tv%eqn_of_state)
3427  else
3428  rho_conv=gv%Rho0
3429  endif
3430  igr0 = 1.0 / (rho_conv * gv%g_Earth)
3431  ssh(i,j) = ssh(i,j) + p_atm(i,j) * igr0
3432  enddo ; enddo
3433  endif
3434 
Here is the caller graph for this function:

◆ calculate_surface_state()

subroutine, public mom::calculate_surface_state ( type(surface), intent(inout)  state,
real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)  u,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szi_(g),szj_(g)), intent(in)  ssh,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(inout)  GV,
type(mom_control_struct), intent(inout)  CS 
)

This subroutine sets the surface (return) properties of the ocean model by setting the appropriate fields in state. Unused fields are set to NULL or are unallocated.

Parameters
[in,out]gocean grid structure
[in,out]gvocean vertical grid structure
[in,out]stateocean surface state
[in]uzonal velocity (m/s)
[in]vmeridional velocity (m/s)
[in]hlayer thickness (m or kg/m2)
[in]sshtime mean surface height (m)
[in,out]cscontrol structure

Definition at line 3441 of file MOM.F90.

Referenced by mom_main(), step_mom(), and step_offline().

3441  type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
3442  type(verticalgrid_type), intent(inout) :: gv !< ocean vertical grid structure
3443  type(surface), intent(inout) :: state !< ocean surface state
3444  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< zonal velocity (m/s)
3445  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< meridional velocity (m/s)
3446  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< layer thickness (m or kg/m2)
3447  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh !< time mean surface height (m)
3448  type(mom_control_struct), intent(inout) :: cs !< control structure
3449 
3450  ! local
3451  real :: depth(szi_(g)) ! distance from the surface (meter)
3452  real :: depth_ml ! depth over which to average to
3453  ! determine mixed layer properties (meter)
3454  real :: dh ! thickness of a layer within mixed layer (meter)
3455  real :: mass ! mass per unit area of a layer (kg/m2)
3456 
3457  real :: hu, hv
3458  integer :: i, j, k, is, ie, js, je, nz, numberoferrors
3459  integer :: isd, ied, jsd, jed
3460  integer :: iscb, iecb, jscb, jecb, isdb, iedb, jsdb, jedb
3461  logical :: localerror
3462  character(240) :: msg
3463 
3464  call calltree_enter("calculate_surface_state(), MOM.F90")
3465  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3466  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3467  iscb = g%iscB ; iecb = g%iecB; jscb = g%jscB ; jecb = g%jecB
3468  isdb = g%isdB ; iedb = g%iedB; jsdb = g%jsdB ; jedb = g%jedB
3469 
3470  if (.not.state%arrays_allocated) then
3471  if (cs%use_temperature) then
3472  allocate(state%SST(isd:ied,jsd:jed)) ; state%SST(:,:) = 0.0
3473  allocate(state%SSS(isd:ied,jsd:jed)) ; state%SSS(:,:) = 0.0
3474  else
3475  allocate(state%sfc_density(isd:ied,jsd:jed)) ; state%sfc_density(:,:) = 0.0
3476  endif
3477  allocate(state%sea_lev(isd:ied,jsd:jed)) ; state%sea_lev(:,:) = 0.0
3478  allocate(state%Hml(isd:ied,jsd:jed)) ; state%Hml(:,:) = 0.0
3479  allocate(state%u(isdb:iedb,jsd:jed)) ; state%u(:,:) = 0.0
3480  allocate(state%v(isd:ied,jsdb:jedb)) ; state%v(:,:) = 0.0
3481 
3482  ! Allocate structures for ocean_mass, ocean_heat, and ocean_salt. This could
3483  ! be wrapped in a run-time flag to disable it for economy, since the 3-d
3484  ! sums are not negligible.
3485  allocate(state%ocean_mass(isd:ied,jsd:jed)) ; state%ocean_mass(:,:) = 0.0
3486  if (cs%use_temperature) then
3487  allocate(state%ocean_heat(isd:ied,jsd:jed)) ; state%ocean_heat(:,:) = 0.0
3488  allocate(state%ocean_salt(isd:ied,jsd:jed)) ; state%ocean_salt(:,:) = 0.0
3489  endif
3490  allocate(state%salt_deficit(isd:ied,jsd:jed)) ; state%salt_deficit(:,:) = 0.0
3491 
3492  state%arrays_allocated = .true.
3493  endif
3494  state%frazil => cs%tv%frazil
3495  state%TempxPmE => cs%tv%TempxPmE
3496  state%internal_heat => cs%tv%internal_heat
3497  if (associated(cs%visc%taux_shelf)) state%taux_shelf => cs%visc%taux_shelf
3498  if (associated(cs%visc%tauy_shelf)) state%tauy_shelf => cs%visc%tauy_shelf
3499 
3500  do j=js,je ; do i=is,ie
3501  state%sea_lev(i,j) = ssh(i,j)
3502  enddo ; enddo
3503 
3504  if (cs%bulkmixedlayer) then
3505  if (cs%use_temperature) then ; do j=js,je ; do i=is,ie
3506  state%SST(i,j) = cs%tv%T(i,j,1)
3507  state%SSS(i,j) = cs%tv%S(i,j,1)
3508  enddo ; enddo ; endif
3509  do j=js,je ; do i=iscb,iecb
3510  state%u(i,j) = u(i,j,1)
3511  enddo ; enddo
3512  do j=jscb,jecb ; do i=is,ie
3513  state%v(i,j) = v(i,j,1)
3514  enddo ; enddo
3515 
3516  if (associated(cs%Hml)) then ; do j=js,je ; do i=is,ie
3517  state%Hml(i,j) = cs%Hml(i,j)
3518  enddo ; enddo ; endif
3519  else
3520 
3521  depth_ml = cs%Hmix
3522  ! Determine the mean tracer properties of the uppermost depth_ml fluid.
3523  !$OMP parallel do default(shared) private(depth,dh)
3524  do j=js,je
3525  do i=is,ie
3526  depth(i) = 0.0
3527  if (cs%use_temperature) then
3528  state%SST(i,j) = 0.0 ; state%SSS(i,j) = 0.0
3529  else
3530  state%sfc_density(i,j) = 0.0
3531  endif
3532  enddo
3533 
3534  do k=1,nz ; do i=is,ie
3535  if (depth(i) + h(i,j,k)*gv%H_to_m < depth_ml) then
3536  dh = h(i,j,k)*gv%H_to_m
3537  elseif (depth(i) < depth_ml) then
3538  dh = depth_ml - depth(i)
3539  else
3540  dh = 0.0
3541  endif
3542  if (cs%use_temperature) then
3543  state%SST(i,j) = state%SST(i,j) + dh * cs%tv%T(i,j,k)
3544  state%SSS(i,j) = state%SSS(i,j) + dh * cs%tv%S(i,j,k)
3545  else
3546  state%sfc_density(i,j) = state%sfc_density(i,j) + dh * gv%Rlay(k)
3547  endif
3548  depth(i) = depth(i) + dh
3549  enddo ; enddo
3550  ! Calculate the average properties of the mixed layer depth.
3551  do i=is,ie
3552  if (depth(i) < gv%H_subroundoff*gv%H_to_m) &
3553  depth(i) = gv%H_subroundoff*gv%H_to_m
3554  if (cs%use_temperature) then
3555  state%SST(i,j) = state%SST(i,j) / depth(i)
3556  state%SSS(i,j) = state%SSS(i,j) / depth(i)
3557  else
3558  state%sfc_density(i,j) = state%sfc_density(i,j) / depth(i)
3559  endif
3560  state%Hml(i,j) = depth(i)
3561  enddo
3562  enddo ! end of j loop
3563 
3564 ! Determine the mean velocities in the uppermost depth_ml fluid.
3565  if (cs%Hmix_UV>0.) then
3566  depth_ml = cs%Hmix_UV
3567  !$OMP parallel do default(shared) private(depth,dh,hv)
3568  do j=jscb,jecb
3569  do i=is,ie
3570  depth(i) = 0.0
3571  state%v(i,j) = 0.0
3572  enddo
3573  do k=1,nz ; do i=is,ie
3574  hv = 0.5 * (h(i,j,k) + h(i,j+1,k)) * gv%H_to_m
3575  if (depth(i) + hv < depth_ml) then
3576  dh = hv
3577  elseif (depth(i) < depth_ml) then
3578  dh = depth_ml - depth(i)
3579  else
3580  dh = 0.0
3581  endif
3582  state%v(i,j) = state%v(i,j) + dh * v(i,j,k)
3583  depth(i) = depth(i) + dh
3584  enddo ; enddo
3585  ! Calculate the average properties of the mixed layer depth.
3586  do i=is,ie
3587  if (depth(i) < gv%H_subroundoff*gv%H_to_m) &
3588  depth(i) = gv%H_subroundoff*gv%H_to_m
3589  state%v(i,j) = state%v(i,j) / depth(i)
3590  enddo
3591  enddo ! end of j loop
3592 
3593  !$OMP parallel do default(shared) private(depth,dh,hu)
3594  do j=js,je
3595  do i=iscb,iecb
3596  depth(i) = 0.0
3597  state%u(i,j) = 0.0
3598  enddo
3599  do k=1,nz ; do i=iscb,iecb
3600  hu = 0.5 * (h(i,j,k) + h(i+1,j,k)) * gv%H_to_m
3601  if (depth(i) + hu < depth_ml) then
3602  dh = hu
3603  elseif (depth(i) < depth_ml) then
3604  dh = depth_ml - depth(i)
3605  else
3606  dh = 0.0
3607  endif
3608  state%u(i,j) = state%u(i,j) + dh * u(i,j,k)
3609  depth(i) = depth(i) + dh
3610  enddo ; enddo
3611  ! Calculate the average properties of the mixed layer depth.
3612  do i=iscb,iecb
3613  if (depth(i) < gv%H_subroundoff*gv%H_to_m) &
3614  depth(i) = gv%H_subroundoff*gv%H_to_m
3615  state%u(i,j) = state%u(i,j) / depth(i)
3616  enddo
3617  enddo ! end of j loop
3618  else ! Hmix_UV<=0.
3619  do j=js,je ; do i=iscb,iecb
3620  state%u(i,j) = u(i,j,1)
3621  enddo ; enddo
3622  do j=jscb,jecb ; do i=is,ie
3623  state%v(i,j) = v(i,j,1)
3624  enddo ; enddo
3625  endif
3626  endif ! end BULKMIXEDLAYER
3627 
3628  if (allocated(state%salt_deficit) .and. associated(cs%tv%salt_deficit)) then
3629  !$OMP parallel do default(shared)
3630  do j=js,je ; do i=is,ie
3631  ! Convert from gSalt to kgSalt
3632  state%salt_deficit(i,j) = 1000.0 * cs%tv%salt_deficit(i,j)
3633  enddo ; enddo
3634  endif
3635 
3636  if (allocated(state%ocean_mass) .and. allocated(state%ocean_heat) .and. &
3637  allocated(state%ocean_salt)) then
3638  !$OMP parallel do default(shared)
3639  do j=js,je ; do i=is,ie
3640  state%ocean_mass(i,j) = 0.0
3641  state%ocean_heat(i,j) = 0.0 ; state%ocean_salt(i,j) = 0.0
3642  enddo ; enddo
3643  !$OMP parallel do default(shared) private(mass)
3644  do j=js,je ; do k=1,nz; do i=is,ie
3645  mass = gv%H_to_kg_m2*h(i,j,k)
3646  state%ocean_mass(i,j) = state%ocean_mass(i,j) + mass
3647  state%ocean_heat(i,j) = state%ocean_heat(i,j) + mass*cs%tv%T(i,j,k)
3648  state%ocean_salt(i,j) = state%ocean_salt(i,j) + &
3649  mass * (1.0e-3*cs%tv%S(i,j,k))
3650  enddo ; enddo ; enddo
3651  else
3652  if (allocated(state%ocean_mass)) then
3653  !$OMP parallel do default(shared)
3654  do j=js,je ; do i=is,ie ; state%ocean_mass(i,j) = 0.0 ; enddo ; enddo
3655  !$OMP parallel do default(shared)
3656  do j=js,je ; do k=1,nz ; do i=is,ie
3657  state%ocean_mass(i,j) = state%ocean_mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
3658  enddo ; enddo ; enddo
3659  endif
3660  if (allocated(state%ocean_heat)) then
3661  !$OMP parallel do default(shared)
3662  do j=js,je ; do i=is,ie ; state%ocean_heat(i,j) = 0.0 ; enddo ; enddo
3663  !$OMP parallel do default(shared) private(mass)
3664  do j=js,je ; do k=1,nz ; do i=is,ie
3665  mass = gv%H_to_kg_m2*h(i,j,k)
3666  state%ocean_heat(i,j) = state%ocean_heat(i,j) + mass*cs%tv%T(i,j,k)
3667  enddo ; enddo ; enddo
3668  endif
3669  if (allocated(state%ocean_salt)) then
3670  !$OMP parallel do default(shared)
3671  do j=js,je ; do i=is,ie ; state%ocean_salt(i,j) = 0.0 ; enddo ; enddo
3672  !$OMP parallel do default(shared) private(mass)
3673  do j=js,je ; do k=1,nz ; do i=is,ie
3674  mass = gv%H_to_kg_m2*h(i,j,k)
3675  state%ocean_salt(i,j) = state%ocean_salt(i,j) + &
3676  mass * (1.0e-3*cs%tv%S(i,j,k))
3677  enddo ; enddo ; enddo
3678  endif
3679  endif
3680 
3681  if (associated(cs%tracer_flow_CSp)) then
3682  if (.not.associated(state%tr_fields)) allocate(state%tr_fields)
3683  call call_tracer_surface_state(state, h, g, cs%tracer_flow_CSp)
3684  endif
3685 
3686  if (cs%check_bad_surface_vals) then
3687  numberoferrors=0 ! count number of errors
3688  do j=js,je; do i=is,ie
3689  if (g%mask2dT(i,j)>0.) then
3690  localerror = state%sea_lev(i,j)<=-g%bathyT(i,j) &
3691  .or. state%sea_lev(i,j)>= cs%bad_val_ssh_max &
3692  .or. state%sea_lev(i,j)<=-cs%bad_val_ssh_max &
3693  .or. state%sea_lev(i,j)+g%bathyT(i,j) < cs%bad_val_column_thickness
3694  if (cs%use_temperature) localerror = localerror &
3695  .or. state%SSS(i,j)<0. &
3696  .or. state%SSS(i,j)>=cs%bad_val_sss_max &
3697  .or. state%SST(i,j)< cs%bad_val_sst_min &
3698  .or. state%SST(i,j)>=cs%bad_val_sst_max
3699  if (localerror) then
3700  numberoferrors=numberoferrors+1
3701  if (numberoferrors<9) then ! Only report details for the first few errors
3702  if (cs%use_temperature) then
3703  write(msg(1:240),'(2(a,i4,x),2(a,f8.3,x),8(a,es11.4,x))') &
3704  'Extreme surface state detected: i=',i,'j=',j, &
3705  'x=',g%geoLonT(i,j),'y=',g%geoLatT(i,j), &
3706  'D=',g%bathyT(i,j), &
3707  'SSH=',state%sea_lev(i,j), &
3708  'SST=',state%SST(i,j), &
3709  'SSS=',state%SSS(i,j), &
3710  'U-=',state%u(i-1,j), &
3711  'U+=',state%u(i,j), &
3712  'V-=',state%v(i,j-1), &
3713  'V+=',state%v(i,j)
3714  else
3715  write(msg(1:240),'(2(a,i4,x),2(a,f8.3,x),6(a,es11.4))') &
3716  'Extreme surface state detected: i=',i,'j=',j, &
3717  'x=',g%geoLonT(i,j),'y=',g%geoLatT(i,j), &
3718  'D=',g%bathyT(i,j), &
3719  'SSH=',state%sea_lev(i,j), &
3720  'U-=',state%u(i-1,j), &
3721  'U+=',state%u(i,j), &
3722  'V-=',state%v(i,j-1), &
3723  'V+=',state%v(i,j)
3724  endif
3725  call mom_error(warning, trim(msg), all_print=.true.)
3726  elseif (numberoferrors==9) then ! Indicate once that there are more errors
3727  call mom_error(warning, 'There were more unreported extreme events!', all_print=.true.)
3728  endif ! numberOfErrors
3729  endif ! localError
3730  endif ! mask2dT
3731  enddo; enddo
3732  call sum_across_pes(numberoferrors)
3733  if (numberoferrors>0) then
3734  write(msg(1:240),'(3(a,i9,x))') 'There were a total of ',numberoferrors, &
3735  'locations detected with extreme surface values!'
3736  call mom_error(fatal, trim(msg))
3737  endif
3738  endif
3739 
3740  call calltree_leave("calculate_surface_state()")
Here is the caller graph for this function:

◆ finish_mom_initialization()

subroutine, public mom::finish_mom_initialization ( type(time_type), intent(in)  Time,
type(directories), intent(in)  dirs,
type(mom_control_struct), pointer  CS,
type(forcing), intent(inout)  fluxes 
)

This subroutine finishes initializing MOM and writes out the initial conditions.

Parameters
[in]timemodel time, used in this routine
[in]dirsstructure with directory paths
cspointer set in this routine to MOM control structure
[in,out]fluxespointers to forcing fields

Definition at line 2345 of file MOM.F90.

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

Referenced by mom_main().

2345  type(time_type), intent(in) :: time !< model time, used in this routine
2346  type(directories), intent(in) :: dirs !< structure with directory paths
2347  type(mom_control_struct), pointer :: cs !< pointer set in this routine to MOM control structure
2348  type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
2349  ! Local variables
2350  type(ocean_grid_type), pointer :: g => null()
2351  type(verticalgrid_type), pointer :: gv => null()
2352  type(mom_restart_cs), pointer :: restart_csp_tmp => null()
2353  real, allocatable :: z_interface(:,:,:) ! Interface heights (meter)
2354  real, allocatable :: eta(:,:) ! Interface heights (meter)
2355  type(vardesc) :: vd
2356 
2357  call cpu_clock_begin(id_clock_init)
2358  call calltree_enter("finish_MOM_initialization()")
2359 
2360  ! Pointers for convenience
2361  g => cs%G ; gv => cs%GV
2362 
2363  ! Write initial conditions
2364  if (cs%write_IC) then
2365  allocate(restart_csp_tmp)
2366  restart_csp_tmp = cs%restart_CSp
2367  allocate(z_interface(szi_(g),szj_(g),szk_(g)+1))
2368  call find_eta(cs%h, cs%tv, gv%g_Earth, g, gv, z_interface)
2369  vd = var_desc("eta","meter","Interface heights",z_grid='i')
2370  call register_restart_field(z_interface, vd, .true., restart_csp_tmp)
2371 
2372  call save_restart(dirs%output_directory, time, g, &
2373  restart_csp_tmp, filename=cs%IC_file, gv=gv)
2374  deallocate(z_interface)
2375  deallocate(restart_csp_tmp)
2376  endif
2377 
2378  call calltree_leave("finish_MOM_initialization()")
2379  call cpu_clock_end(id_clock_init)
2380 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ initialize_mom()

subroutine, public mom::initialize_mom ( type(time_type), intent(inout), target  Time,
type(param_file_type), intent(out)  param_file,
type(directories), intent(out)  dirs,
type(mom_control_struct), pointer  CS,
type(time_type), intent(in), optional  Time_in,
logical, intent(out), optional  offline_tracer_mode 
)

This subroutine initializes MOM.

Parameters
[in,out]timemodel time, set in this routine
[out]param_filestructure indicating paramater file to parse
[out]dirsstructure with directory paths
cspointer set in this routine to MOM control structure
[in]time_intime passed to MOM_initialize_state when model is not being started from a restart file
[out]offline_tracer_modeTrue if tracers are being run offline

Definition at line 1480 of file MOM.F90.

References mom_diabatic_driver::adiabatic_driver_init(), mom_boundary_update::call_obc_register(), mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), mom_error_handler::calltree_waypoint(), mom_transcribe_grid::copy_dyngrid_to_mom_grid(), mom_transcribe_grid::copy_mom_grid_to_dyngrid(), mom_dyn_horgrid::create_dyn_horgrid(), mom_dyn_horgrid::destroy_dyn_horgrid(), mom_diag_mediator::diag_masks_set(), mom_obsolete_params::find_obsolete_params(), mom_get_input::get_mom_input(), mom_hor_index::hor_index_init(), id_clock_init, id_clock_mom_init, id_clock_pass_init, mom_ale_sponge::init_ale_sponge_diags(), mom_sponge::init_sponge_diags(), mom_tracer_registry::lock_tracer_registry(), mom_meke::meke_alloc_register_restart(), mom_meke::meke_init(), mom_mixed_layer_restrat::mixedlayer_restrat_register_restarts(), mom_grid::mom_grid_end(), mom_grid::mom_grid_init(), mom_coord_initialization::mom_initialize_coord(), mom_fixed_initialization::mom_initialize_fixed(), mom_state_initialization::mom_initialize_state(), mom_timing_init(), mom_neutral_diffusion::neutral_diffusion_diag_init(), register_diags(), register_diags_ts_tendency(), register_diags_ts_vardec(), mom_obsolete_diagnostics::register_obsolete_diagnostics(), mom_restart::restart_init(), mom_diag_mediator::set_axes_info(), set_restart_fields(), mom_set_visc::set_visc_register_restarts(), mom_thickness_diffuse::thickness_diffuse_init(), mom_unit_tests::unit_tests(), and write_static_fields().

Referenced by mom_main().

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

◆ mom_end()

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

End of model.

Parameters
csMOM control structure

Definition at line 3746 of file MOM.F90.

Referenced by mom_main().

3746  type(mom_control_struct), pointer :: cs !< MOM control structure
3747 
3748  if (cs%use_ALE_algorithm) then
3749  call ale_end(cs%ALE_CSp)
3750  endif
3751 
3752  dealloc_(cs%u) ; dealloc_(cs%v) ; dealloc_(cs%h)
3753  dealloc_(cs%uh) ; dealloc_(cs%vh)
3754 
3755  if (cs%use_temperature) then
3756  dealloc_(cs%T) ; cs%tv%T => null() ; dealloc_(cs%S) ; cs%tv%S => null()
3757  endif
3758  if (associated(cs%tv%frazil)) deallocate(cs%tv%frazil)
3759  if (associated(cs%tv%salt_deficit)) deallocate(cs%tv%salt_deficit)
3760  if (associated(cs%Hml)) deallocate(cs%Hml)
3761 
3762  call tracer_advect_end(cs%tracer_adv_CSp)
3763  call tracer_hor_diff_end(cs%tracer_diff_CSp)
3764  call tracer_registry_end(cs%tracer_Reg)
3765  call tracer_flow_control_end(cs%tracer_flow_CSp)
3766 
3767  if(cs%offline_tracer_mode) then
3768  call offline_transport_end(cs%offline_CSp)
3769  endif
3770 
3771  dealloc_(cs%uhtr) ; dealloc_(cs%vhtr)
3772  if (cs%split) then ; if (cs%legacy_split) then
3773  call end_dyn_legacy_split(cs%dyn_legacy_split_CSp)
3774  else
3775  call end_dyn_split_rk2(cs%dyn_split_RK2_CSp)
3776  endif ; else ; if (cs%use_RK2) then
3777  call end_dyn_unsplit_rk2(cs%dyn_unsplit_RK2_CSp)
3778  else
3779  call end_dyn_unsplit(cs%dyn_unsplit_CSp)
3780  endif ; endif
3781  dealloc_(cs%ave_ssh)
3782  if (associated(cs%update_OBC_CSp)) call obc_register_end(cs%update_OBC_CSp)
3783 
3784  call verticalgridend(cs%GV)
3785  call mom_grid_end(cs%G)
3786 
3787  deallocate(cs)
3788 
Here is the caller graph for this function:

◆ mom_timing_init()

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

This subroutine sets up clock IDs for timing various subroutines.

Parameters
[in]cscontrol structure set up by initialize_MOM.

Definition at line 2763 of file MOM.F90.

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

Referenced by initialize_mom().

2763  type(mom_control_struct), intent(in) :: cs !< control structure set up by initialize_MOM.
2764 
2765  id_clock_ocean = cpu_clock_id('Ocean', grain=clock_component)
2766  id_clock_dynamics = cpu_clock_id('Ocean dynamics', grain=clock_subcomponent)
2767  id_clock_thermo = cpu_clock_id('Ocean thermodynamics and tracers', grain=clock_subcomponent)
2768  id_clock_other = cpu_clock_id('Ocean Other', grain=clock_subcomponent)
2769  id_clock_tracer = cpu_clock_id('(Ocean tracer advection)', grain=clock_module_driver)
2770  if (.not.cs%adiabatic) &
2771  id_clock_diabatic = cpu_clock_id('(Ocean diabatic driver)', grain=clock_module_driver)
2772 
2773  id_clock_continuity = cpu_clock_id('(Ocean continuity equation *)', grain=clock_module)
2774  id_clock_bbl_visc = cpu_clock_id('(Ocean set BBL viscosity)', grain=clock_module)
2775  id_clock_pass = cpu_clock_id('(Ocean message passing *)', grain=clock_module)
2776  id_clock_mom_init = cpu_clock_id('(Ocean MOM_initialize_state)', grain=clock_module)
2777  id_clock_pass_init = cpu_clock_id('(Ocean init message passing *)', grain=clock_routine)
2778  if (cs%thickness_diffuse) &
2779  id_clock_thick_diff = cpu_clock_id('(Ocean thickness diffusion *)', grain=clock_module)
2780 !if (CS%mixedlayer_restrat) &
2781  id_clock_ml_restrat = cpu_clock_id('(Ocean mixed layer restrat)', grain=clock_module)
2782  id_clock_diagnostics = cpu_clock_id('(Ocean collective diagnostics)', grain=clock_module)
2783  id_clock_z_diag = cpu_clock_id('(Ocean Z-space diagnostics)', grain=clock_module)
2784  id_clock_ale = cpu_clock_id('(Ocean ALE)', grain=clock_module)
2785  if(cs%offline_tracer_mode) then
2786  id_clock_offline_tracer = cpu_clock_id('Ocean offline tracers', grain=clock_subcomponent)
2787  endif
2788 
Here is the caller graph for this function:

◆ post_diags_ts_vardec()

subroutine mom::post_diags_ts_vardec ( type(ocean_grid_type), intent(in)  G,
type(mom_control_struct), intent(in)  CS,
real, intent(in)  dt 
)
private

Calculate and post variance decay diagnostics for temp/salt.

Parameters
[in]gocean grid structure
[in]cscontrol structure
[in]dttotal time step

Definition at line 3017 of file MOM.F90.

Referenced by step_mom_thermo().

3017  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
3018  type(mom_control_struct), intent(in) :: cs !< control structure
3019  real, intent(in) :: dt !< total time step
3020 
3021  real :: work(szi_(g),szj_(g),szk_(g))
3022  real :: idt
3023  integer :: i, j, k, is, ie, js, je, nz
3024  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3025 
3026  idt = 0.; if (dt/=0.) idt = 1.0 / dt ! The "if" is in case the diagnostic is called for a zero length interval
3027 
3028  if (cs%id_T_vardec > 0) then
3029  do k=1,nz ; do j=js,je ; do i=is,ie
3030  work(i,j,k) = (cs%T_squared(i,j,k) - cs%tv%T(i,j,k)**2) * idt
3031  enddo ; enddo ; enddo
3032  call post_data(cs%id_T_vardec, work, cs%diag)
3033  endif
3034 
3035  if (cs%id_S_vardec > 0) then
3036  do k=1,nz ; do j=js,je ; do i=is,ie
3037  work(i,j,k) = (cs%S_squared(i,j,k) - cs%tv%S(i,j,k)**2) * idt
3038  enddo ; enddo ; enddo
3039  call post_data(cs%id_S_vardec, work, cs%diag)
3040  endif
Here is the caller graph for this function:

◆ post_integrated_diagnostics()

subroutine mom::post_integrated_diagnostics ( type(mom_control_struct), intent(in)  CS,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(diag_ctrl), intent(in)  diag,
real, intent(in)  dt_int,
type(thermo_var_ptrs), intent(in)  tv,
type(forcing), intent(in)  fluxes 
)
private

This routine posts diagnostics of various integrated quantities.

Parameters
[in]cscontrol structure
[in]gocean grid structure
[in]gvocean vertical grid structure
[in]diagregulates diagnostic output
[in]dt_inttotal time step associated with these diagnostics, in s.
[in]tvA structure pointing to various thermodynamic variables type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
[in]fluxespointers to forcing fields

Definition at line 3045 of file MOM.F90.

Referenced by step_mom().

3045  type(mom_control_struct), intent(in) :: cs !< control structure
3046  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
3047  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
3048  type(diag_ctrl), intent(in) :: diag !< regulates diagnostic output
3049  real, intent(in) :: dt_int !< total time step associated with these diagnostics, in s.
3050  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
3051  type(forcing), intent(in) :: fluxes !< pointers to forcing fields
3052 
3053  real, allocatable, dimension(:,:) :: &
3054  tmp, & ! temporary 2d field
3055  zos, & ! dynamic sea lev (zero area mean) from inverse-barometer adjusted ssh (meter)
3056  zossq, & ! square of zos (m^2)
3057  sfc_speed, & ! sea surface speed at h-points (m/s)
3058  frazil_ave, & ! average frazil heat flux required to keep temp above freezing (W/m2)
3059  salt_deficit_ave, & ! average salt flux required to keep salinity above 0.01ppt (gSalt m-2 s-1)
3060  heat_pme_ave, & ! average effective heat flux into the ocean due to
3061  ! the exchange of water with other components, times the
3062  ! heat capacity of water, in W m-2.
3063  intern_heat_ave ! avg heat flux into ocean from geothermal or
3064  ! other internal heat sources (W/m2)
3065  real :: i_time_int ! The inverse of the time interval in s-1.
3066  real :: zos_area_mean, volo, ssh_ga
3067  integer :: i, j, k, is, ie, js, je, nz! , Isq, Ieq, Jsq, Jeq
3068 
3069  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3070 
3071  ! area mean SSH
3072  if (cs%id_ssh_ga > 0) then
3073  ssh_ga = global_area_mean(cs%ave_ssh, g)
3074  call post_data(cs%id_ssh_ga, ssh_ga, diag)
3075  endif
3076 
3077  i_time_int = 1.0 / dt_int
3078  if (cs%id_ssh > 0) &
3079  call post_data(cs%id_ssh, cs%ave_ssh, diag, mask=g%mask2dT)
3080 
3081  ! post the dynamic sea level, zos, and zossq.
3082  ! zos is ave_ssh with sea ice inverse barometer removed,
3083  ! and with zero global area mean.
3084  if(cs%id_zos > 0 .or. cs%id_zossq > 0) then
3085  allocate(zos(g%isd:g%ied,g%jsd:g%jed))
3086  zos(:,:) = 0.0
3087  do j=js,je ; do i=is,ie
3088  zos(i,j) = cs%ave_ssh(i,j)
3089  enddo ; enddo
3090  if (ASSOCIATED(fluxes%p_surf)) then
3091  do j=js,je ; do i=is,ie
3092  zos(i,j) = zos(i,j) + g%mask2dT(i,j)*fluxes%p_surf(i,j) / &
3093  (gv%Rho0 * gv%g_Earth)
3094  enddo ; enddo
3095  endif
3096  zos_area_mean = global_area_mean(zos, g)
3097  do j=js,je ; do i=is,ie
3098  zos(i,j) = zos(i,j) - g%mask2dT(i,j)*zos_area_mean
3099  enddo ; enddo
3100  if(cs%id_zos > 0) then
3101  call post_data(cs%id_zos, zos, diag, mask=g%mask2dT)
3102  endif
3103  if(cs%id_zossq > 0) then
3104  allocate(zossq(g%isd:g%ied,g%jsd:g%jed))
3105  zossq(:,:) = 0.0
3106  do j=js,je ; do i=is,ie
3107  zossq(i,j) = zos(i,j)*zos(i,j)
3108  enddo ; enddo
3109  call post_data(cs%id_zossq, zossq, diag, mask=g%mask2dT)
3110  deallocate(zossq)
3111  endif
3112  deallocate(zos)
3113  endif
3114 
3115  ! post total volume of the liquid ocean
3116  if(cs%id_volo > 0) then
3117  allocate(tmp(g%isd:g%ied,g%jsd:g%jed))
3118  do j=js,je ; do i=is,ie
3119  tmp(i,j) = g%mask2dT(i,j)*(cs%ave_ssh(i,j) + g%bathyT(i,j))
3120  enddo ; enddo
3121  volo = global_area_integral(tmp, g)
3122  call post_data(cs%id_volo, volo, diag)
3123  deallocate(tmp)
3124  endif
3125 
3126  ! post frazil
3127  if (ASSOCIATED(tv%frazil) .and. (cs%id_fraz > 0)) then
3128  allocate(frazil_ave(g%isd:g%ied,g%jsd:g%jed))
3129  do j=js,je ; do i=is,ie
3130  frazil_ave(i,j) = tv%frazil(i,j) * i_time_int
3131  enddo ; enddo
3132  call post_data(cs%id_fraz, frazil_ave, diag, mask=g%mask2dT)
3133  deallocate(frazil_ave)
3134  endif
3135 
3136  ! post the salt deficit
3137  if (ASSOCIATED(tv%salt_deficit) .and. (cs%id_salt_deficit > 0)) then
3138  allocate(salt_deficit_ave(g%isd:g%ied,g%jsd:g%jed))
3139  do j=js,je ; do i=is,ie
3140  salt_deficit_ave(i,j) = tv%salt_deficit(i,j) * i_time_int
3141  enddo ; enddo
3142  call post_data(cs%id_salt_deficit, salt_deficit_ave, diag, mask=g%mask2dT)
3143  deallocate(salt_deficit_ave)
3144  endif
3145 
3146  ! post temperature of P-E+R
3147  if (ASSOCIATED(tv%TempxPmE) .and. (cs%id_Heat_PmE > 0)) then
3148  allocate(heat_pme_ave(g%isd:g%ied,g%jsd:g%jed))
3149  do j=js,je ; do i=is,ie
3150  heat_pme_ave(i,j) = tv%TempxPmE(i,j) * (tv%C_p * i_time_int)
3151  enddo ; enddo
3152  call post_data(cs%id_Heat_PmE, heat_pme_ave, diag, mask=g%mask2dT)
3153  deallocate(heat_pme_ave)
3154  endif
3155 
3156  ! post geothermal heating or internal heat source/sinks
3157  if (ASSOCIATED(tv%internal_heat) .and. (cs%id_intern_heat > 0)) then
3158  allocate(intern_heat_ave(g%isd:g%ied,g%jsd:g%jed))
3159  do j=js,je ; do i=is,ie
3160  intern_heat_ave(i,j) = tv%internal_heat(i,j) * (tv%C_p * i_time_int)
3161  enddo ; enddo
3162  call post_data(cs%id_intern_heat, intern_heat_ave, diag, mask=g%mask2dT)
3163  deallocate(intern_heat_ave)
3164  endif
3165 
Here is the caller graph for this function:

◆ post_surface_diagnostics()

subroutine mom::post_surface_diagnostics ( type(mom_control_struct), intent(in)  CS,
type(ocean_grid_type), intent(in)  G,
type(diag_ctrl), intent(in)  diag,
type(surface), intent(in)  state 
)
private

This routine posts diagnostics of various ocean surface quantities.

Parameters
[in]cscontrol structure
[in]gocean grid structure
[in]diagregulates diagnostic output
[in]stateocean surface state

Definition at line 3170 of file MOM.F90.

Referenced by step_mom().

3170  type(mom_control_struct), intent(in) :: cs !< control structure
3171  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
3172  type(diag_ctrl), intent(in) :: diag !< regulates diagnostic output
3173  type(surface), intent(in) :: state !< ocean surface state
3174 
3175  real, dimension(SZI_(G),SZJ_(G)) :: &
3176  pottemp, & ! TEOS10 potential temperature (deg C)
3177  pracsal, & ! TEOS10 practical salinity
3178  sst_sq, & ! Surface temperature squared, in degC^2
3179  sss_sq, & ! Surface salinity squared, in salnity units^2
3180  sfc_speed ! sea surface speed at h-points (m/s)
3181 
3182  integer :: i, j, is, ie, js, je
3183 
3184  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3185 
3186  if (.NOT.cs%use_conT_absS) then
3187  !Internal T&S variables are assumed to be potential&practical
3188  if (cs%id_sst > 0) call post_data(cs%id_sst, state%SST, diag, mask=g%mask2dT)
3189  if (cs%id_sss > 0) call post_data(cs%id_sss, state%SSS, diag, mask=g%mask2dT)
3190  else
3191  !Internal T&S variables are assumed to be conservative&absolute
3192  if (cs%id_sstcon > 0) call post_data(cs%id_sstcon, state%SST, diag, mask=g%mask2dT)
3193  if (cs%id_sssabs > 0) call post_data(cs%id_sssabs, state%SSS, diag, mask=g%mask2dT)
3194  !Using TEOS-10 function calls convert T&S diagnostics
3195  !from conservative temp to potential temp and
3196  !from absolute salinity to practical salinity
3197  do j=js,je ; do i=is,ie
3198  pracsal(i,j) = gsw_sp_from_sr(state%SSS(i,j))
3199  pottemp(i,j) = gsw_pt_from_ct(state%SSS(i,j),state%SST(i,j))
3200  enddo ; enddo
3201  if (cs%id_sst > 0) call post_data(cs%id_sst, pottemp, diag, mask=g%mask2dT)
3202  if (cs%id_sss > 0) call post_data(cs%id_sss, pracsal, diag, mask=g%mask2dT)
3203  endif
3204 
3205  if (cs%id_sst_sq > 0) then
3206  do j=js,je ; do i=is,ie
3207  sst_sq(i,j) = state%SST(i,j)*state%SST(i,j)
3208  enddo ; enddo
3209  call post_data(cs%id_sst_sq, sst_sq, diag, mask=g%mask2dT)
3210  endif
3211  if (cs%id_sss_sq > 0) then
3212  do j=js,je ; do i=is,ie
3213  sss_sq(i,j) = state%SSS(i,j)*state%SSS(i,j)
3214  enddo ; enddo
3215  call post_data(cs%id_sss_sq, sss_sq, diag, mask=g%mask2dT)
3216  endif
3217 
3218  if (cs%id_ssu > 0) &
3219  call post_data(cs%id_ssu, state%u, diag, mask=g%mask2dCu)
3220  if (cs%id_ssv > 0) &
3221  call post_data(cs%id_ssv, state%v, diag, mask=g%mask2dCv)
3222 
3223  if (cs%id_speed > 0) then
3224  do j=js,je ; do i=is,ie
3225  sfc_speed(i,j) = sqrt(0.5*(state%u(i-1,j)**2 + state%u(i,j)**2) + &
3226  0.5*(state%v(i,j-1)**2 + state%v(i,j)**2))
3227  enddo ; enddo
3228  call post_data(cs%id_speed, sfc_speed, diag, mask=g%mask2dT)
3229  endif
3230 
Here is the caller graph for this function:

◆ post_transport_diagnostics()

subroutine mom::post_transport_diagnostics ( type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(mom_control_struct), intent(in)  CS,
type(diag_ctrl), intent(inout)  diag,
real, intent(in)  dt_trans,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h_pre_dyn 
)
private

This routine posts diagnostics of the transports, including the subgridscale contributions.

Parameters
[in,out]gocean grid structure
[in]gvocean vertical grid structure
[in]cscontrol structure
[in,out]diagregulates diagnostic output
[in]dt_transtotal time step associated with the transports, in s.
[in]hThe updated layer thicknesses, in H
[in]h_pre_dynThe thickness before the transports, in H.

Definition at line 2794 of file MOM.F90.

References id_clock_z_diag.

Referenced by step_mom().

2794  type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
2795  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
2796  type(mom_control_struct), intent(in) :: cs !< control structure
2797  type(diag_ctrl), intent(inout) :: diag !< regulates diagnostic output
2798  real , intent(in) :: dt_trans !< total time step associated with the transports, in s.
2799  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
2800  intent(in) :: h !< The updated layer thicknesses, in H
2801  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
2802  intent(in) :: h_pre_dyn !< The thickness before the transports, in H.
2803 
2804  real, dimension(SZIB_(G), SZJ_(G)) :: umo2d ! Diagnostics of integrated mass transport, in kg s-1
2805  real, dimension(SZI_(G), SZJB_(G)) :: vmo2d ! Diagnostics of integrated mass transport, in kg s-1
2806  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)) :: umo ! Diagnostics of layer mass transport, in kg s-1
2807  real, dimension(SZI_(G), SZJB_(G), SZK_(G)) :: vmo ! Diagnostics of layer mass transport, in kg s-1
2808  real :: h_to_kg_m2_dt ! A conversion factor from accumulated transports to fluxes, in kg m-2 H-1 s-1.
2809  integer :: i, j, k, is, ie, js, je, nz
2810  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2811 
2812  call cpu_clock_begin(id_clock_z_diag)
2813  call calculate_z_transport(cs%uhtr, cs%vhtr, h, dt_trans, g, gv, &
2814  cs%diag_to_Z_CSp)
2815  call cpu_clock_end(id_clock_z_diag)
2816 
2817  ! Post mass transports, including SGS
2818  ! Build the remap grids using the layer thicknesses from before the dynamics
2819  call diag_update_remap_grids(diag, alt_h = h_pre_dyn)
2820 
2821  h_to_kg_m2_dt = gv%H_to_kg_m2 / dt_trans
2822  if (cs%id_umo_2d > 0) then
2823  umo2d(:,:) = 0.0
2824  do k=1,nz ; do j=js,je ; do i=is-1,ie
2825  umo2d(i,j) = umo2d(i,j) + cs%uhtr(i,j,k) * h_to_kg_m2_dt
2826  enddo ; enddo ; enddo
2827  call post_data(cs%id_umo_2d, umo2d, diag)
2828  endif
2829  if (cs%id_umo > 0) then
2830  ! Convert to kg/s. Modifying the array for diagnostics is allowed here since it is set to zero immediately below
2831  do k=1,nz ; do j=js,je ; do i=is-1,ie
2832  umo(i,j,k) = cs%uhtr(i,j,k) * h_to_kg_m2_dt
2833  enddo ; enddo ; enddo
2834  call post_data(cs%id_umo, umo, diag, alt_h = h_pre_dyn)
2835  endif
2836  if (cs%id_vmo_2d > 0) then
2837  vmo2d(:,:) = 0.0
2838  do k=1,nz ; do j=js-1,je ; do i=is,ie
2839  vmo2d(i,j) = vmo2d(i,j) + cs%vhtr(i,j,k) * h_to_kg_m2_dt
2840  enddo ; enddo ; enddo
2841  call post_data(cs%id_vmo_2d, vmo2d, diag)
2842  endif
2843  if (cs%id_vmo > 0) then
2844  ! Convert to kg/s. Modifying the array for diagnostics is allowed here since it is set to zero immediately below
2845  do k=1,nz ; do j=js-1,je ; do i=is,ie
2846  vmo(i,j,k) = cs%vhtr(i,j,k) * h_to_kg_m2_dt
2847  enddo ; enddo ; enddo
2848  call post_data(cs%id_vmo, vmo, diag, alt_h = h_pre_dyn)
2849  endif
2850 
2851  if (cs%id_uhtr > 0) call post_data(cs%id_uhtr, cs%uhtr, diag, alt_h = h_pre_dyn)
2852  if (cs%id_vhtr > 0) call post_data(cs%id_vhtr, cs%vhtr, diag, alt_h = h_pre_dyn)
2853 
Here is the caller graph for this function:

◆ post_ts_diagnostics()

subroutine mom::post_ts_diagnostics ( type(mom_control_struct), intent(inout)  CS,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(thermo_var_ptrs), intent(in)  tv,
type(diag_ctrl), intent(in)  diag,
real, intent(in)  dt 
)
private

Post diagnostics of temperatures and salinities, their fluxes, and tendencies.

Parameters
[in,out]cscontrol structure
[in]gocean grid structure
[in]gvocean vertical grid structure
[in]tvA structure pointing to various thermodynamic variables type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
[in]diagregulates diagnostic output
[in]dttotal time step for T,S update

Definition at line 2858 of file MOM.F90.

Referenced by step_mom().

2858  type(mom_control_struct), intent(inout) :: cs !< control structure
2859  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
2860  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
2861  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
2862  type(diag_ctrl), intent(in) :: diag !< regulates diagnostic output
2863  real, intent(in) :: dt !< total time step for T,S update
2864 
2865  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: pottemp, pracsal !TEOS10 Diagnostics
2866  real :: work3d(szi_(g),szj_(g),szk_(g))
2867  real :: work2d(szi_(g),szj_(g))
2868  real :: idt, ppt2mks
2869  integer :: i, j, k, is, ie, js, je, nz
2870  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2871 
2872  if (.NOT. cs%use_conT_absS) then
2873  !Internal T&S variables are assumed to be potential&practical
2874  if (cs%id_T > 0) call post_data(cs%id_T, tv%T, diag)
2875  if (cs%id_S > 0) call post_data(cs%id_S, tv%S, diag)
2876 
2877  if (cs%id_tob > 0) call post_data(cs%id_tob, tv%T(:,:,g%ke), diag, mask=g%mask2dT)
2878  if (cs%id_sob > 0) call post_data(cs%id_sob, tv%S(:,:,g%ke), diag, mask=g%mask2dT)
2879  else
2880  !Internal T&S variables are assumed to be conservative&absolute
2881  if (cs%id_Tcon > 0) call post_data(cs%id_Tcon, tv%T, diag)
2882  if (cs%id_Sabs > 0) call post_data(cs%id_Sabs, tv%S, diag)
2883  !Using TEOS-10 function calls convert T&S diagnostics
2884  !from conservative temp to potential temp and
2885  !from absolute salinity to practical salinity
2886  do k=1,nz ; do j=js,je ; do i=is,ie
2887  pracsal(i,j,k) = gsw_sp_from_sr(tv%S(i,j,k))
2888  pottemp(i,j,k) = gsw_pt_from_ct(tv%S(i,j,k),tv%T(i,j,k))
2889  enddo; enddo ; enddo
2890  if (cs%id_T > 0) call post_data(cs%id_T, pottemp, diag)
2891  if (cs%id_S > 0) call post_data(cs%id_S, pracsal, diag)
2892  if (cs%id_tob > 0) call post_data(cs%id_tob, pottemp(:,:,g%ke), diag, mask=g%mask2dT)
2893  if (cs%id_sob > 0) call post_data(cs%id_sob, pracsal(:,:,g%ke), diag, mask=g%mask2dT)
2894  endif
2895 
2896  if (cs%id_Tadx > 0) call post_data(cs%id_Tadx, cs%T_adx, diag)
2897  if (cs%id_Tady > 0) call post_data(cs%id_Tady, cs%T_ady, diag)
2898  if (cs%id_Tdiffx > 0) call post_data(cs%id_Tdiffx, cs%T_diffx, diag)
2899  if (cs%id_Tdiffy > 0) call post_data(cs%id_Tdiffy, cs%T_diffy, diag)
2900 
2901  if (cs%id_Sadx > 0) call post_data(cs%id_Sadx, cs%S_adx, diag)
2902  if (cs%id_Sady > 0) call post_data(cs%id_Sady, cs%S_ady, diag)
2903  if (cs%id_Sdiffx > 0) call post_data(cs%id_Sdiffx, cs%S_diffx, diag)
2904  if (cs%id_Sdiffy > 0) call post_data(cs%id_Sdiffy, cs%S_diffy, diag)
2905 
2906  if (cs%id_Tadx_2d > 0) call post_data(cs%id_Tadx_2d, cs%T_adx_2d, diag)
2907  if (cs%id_Tady_2d > 0) call post_data(cs%id_Tady_2d, cs%T_ady_2d, diag)
2908  if (cs%id_Tdiffx_2d > 0) call post_data(cs%id_Tdiffx_2d, cs%T_diffx_2d, diag)
2909  if (cs%id_Tdiffy_2d > 0) call post_data(cs%id_Tdiffy_2d, cs%T_diffy_2d, diag)
2910 
2911  if (cs%id_Sadx_2d > 0) call post_data(cs%id_Sadx_2d, cs%S_adx_2d, diag)
2912  if (cs%id_Sady_2d > 0) call post_data(cs%id_Sady_2d, cs%S_ady_2d, diag)
2913  if (cs%id_Sdiffx_2d > 0) call post_data(cs%id_Sdiffx_2d, cs%S_diffx_2d, diag)
2914  if (cs%id_Sdiffy_2d > 0) call post_data(cs%id_Sdiffy_2d, cs%S_diffy_2d, diag)
2915 
2916  if(.not. cs%tendency_diagnostics) return
2917 
2918  idt = 0.; if (dt/=0.) idt = 1.0 / dt ! The "if" is in case the diagnostic is called for a zero length interval
2919  ppt2mks = 0.001
2920  work3d(:,:,:) = 0.0
2921  work2d(:,:) = 0.0
2922 
2923  ! Diagnose tendency of heat from convergence of lateral advective,
2924  ! fluxes, where advective transport arises from residual mean velocity.
2925  if (cs%id_T_advection_xy > 0 .or. cs%id_T_advection_xy_2d > 0) then
2926  do k=1,nz ; do j=js,je ; do i=is,ie
2927  work3d(i,j,k) = cs%T_advection_xy(i,j,k) * gv%H_to_kg_m2 * tv%C_p
2928  enddo ; enddo ; enddo
2929  if (cs%id_T_advection_xy > 0) call post_data(cs%id_T_advection_xy, work3d, diag)
2930  if (cs%id_T_advection_xy_2d > 0) then
2931  do j=js,je ; do i=is,ie
2932  work2d(i,j) = 0.0
2933  do k=1,nz
2934  work2d(i,j) = work2d(i,j) + work3d(i,j,k)
2935  enddo
2936  enddo ; enddo
2937  call post_data(cs%id_T_advection_xy_2d, work2d, diag)
2938  endif
2939  endif
2940 
2941  ! Diagnose tendency of salt from convergence of lateral advective
2942  ! fluxes, where advective transport arises from residual mean velocity.
2943  if (cs%id_S_advection_xy > 0 .or. cs%id_S_advection_xy_2d > 0) then
2944  do k=1,nz ; do j=js,je ; do i=is,ie
2945  work3d(i,j,k) = cs%S_advection_xy(i,j,k) * gv%H_to_kg_m2 * ppt2mks
2946  enddo ; enddo ; enddo
2947  if (cs%id_S_advection_xy > 0) call post_data(cs%id_S_advection_xy, work3d, diag)
2948  if (cs%id_S_advection_xy_2d > 0) then
2949  do j=js,je ; do i=is,ie
2950  work2d(i,j) = 0.0
2951  do k=1,nz
2952  work2d(i,j) = work2d(i,j) + work3d(i,j,k)
2953  enddo
2954  enddo ; enddo
2955  call post_data(cs%id_S_advection_xy_2d, work2d, diag)
2956  endif
2957  endif
2958 
2959  ! diagnose net tendency for temperature over a time step and update T_prev
2960  if (cs%id_T_tendency > 0) then
2961  do k=1,nz ; do j=js,je ; do i=is,ie
2962  work3d(i,j,k) = (tv%T(i,j,k) - cs%T_prev(i,j,k))*idt
2963  cs%T_prev(i,j,k) = tv%T(i,j,k)
2964  enddo ; enddo ; enddo
2965  call post_data(cs%id_T_tendency, work3d, diag)
2966  endif
2967 
2968  ! diagnose net tendency for salinity over a time step and update S_prev
2969  if (cs%id_S_tendency > 0) then
2970  do k=1,nz ; do j=js,je ; do i=is,ie
2971  work3d(i,j,k) = (tv%S(i,j,k) - cs%S_prev(i,j,k))*idt
2972  cs%S_prev(i,j,k) = tv%S(i,j,k)
2973  enddo ; enddo ; enddo
2974  call post_data(cs%id_S_tendency, work3d, diag)
2975  endif
2976 
2977  ! diagnose net tendency for heat content of a grid cell over a time step and update Th_prev
2978  if (cs%id_Th_tendency > 0 .or. cs%id_Th_tendency_2d > 0) then
2979  do k=1,nz ; do j=js,je ; do i=is,ie
2980  work3d(i,j,k) = (tv%T(i,j,k)*cs%h(i,j,k) - cs%Th_prev(i,j,k)) * idt * gv%H_to_kg_m2 * tv%C_p
2981  cs%Th_prev(i,j,k) = tv%T(i,j,k)*cs%h(i,j,k)
2982  enddo ; enddo ; enddo
2983  if (cs%id_Th_tendency > 0) call post_data(cs%id_Th_tendency, work3d, diag)
2984  if (cs%id_Th_tendency_2d > 0) then
2985  do j=js,je ; do i=is,ie
2986  work2d(i,j) = 0.0
2987  do k=1,nz
2988  work2d(i,j) = work2d(i,j) + work3d(i,j,k)
2989  enddo
2990  enddo ; enddo
2991  call post_data(cs%id_Th_tendency_2d, work2d, diag)
2992  endif
2993  endif
2994 
2995  ! diagnose net tendency for salt content of a grid cell over a time step and update Sh_prev
2996  if (cs%id_Sh_tendency > 0 .or. cs%id_Sh_tendency_2d > 0) then
2997  do k=1,nz ; do j=js,je ; do i=is,ie
2998  work3d(i,j,k) = (tv%S(i,j,k)*cs%h(i,j,k) - cs%Sh_prev(i,j,k)) * idt * gv%H_to_kg_m2 * ppt2mks
2999  cs%Sh_prev(i,j,k) = tv%S(i,j,k)*cs%h(i,j,k)
3000  enddo ; enddo ; enddo
3001  if (cs%id_Sh_tendency > 0) call post_data(cs%id_Sh_tendency, work3d, diag)
3002  if (cs%id_Sh_tendency_2d > 0) then
3003  do j=js,je ; do i=is,ie
3004  work2d(i,j) = 0.0
3005  do k=1,nz
3006  work2d(i,j) = work2d(i,j) + work3d(i,j,k)
3007  enddo
3008  enddo ; enddo
3009  call post_data(cs%id_Sh_tendency_2d, work2d, diag)
3010  endif
3011  endif
3012 
Here is the caller graph for this function:

◆ register_diags()

subroutine mom::register_diags ( type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(inout)  GV,
type(mom_control_struct), pointer  CS,
type(accel_diag_ptrs), intent(inout)  ADp 
)
private

Register the diagnostics.

Parameters
[in]timecurrent model time
[in,out]gocean grid structu
[in,out]gvocean vertical grid structure
cscontrol structure set up by initialize_MOM
[in,out]adpstructure pointing to accelerations in momentum equation

Definition at line 2385 of file MOM.F90.

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

Referenced by initialize_mom().

2385  type(time_type), intent(in) :: time !< current model time
2386  type(ocean_grid_type), intent(inout) :: g !< ocean grid structu
2387  type(verticalgrid_type), intent(inout) :: gv !< ocean vertical grid structure
2388  type(mom_control_struct), pointer :: cs !< control structure set up by initialize_MOM
2389  type(accel_diag_ptrs), intent(inout) :: adp !< structure pointing to accelerations in momentum equation
2390 
2391  character(len=48) :: thickness_units, flux_units, t_flux_units, s_flux_units
2392  type(diag_ctrl), pointer :: diag
2393  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz
2394  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
2395  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
2396 
2397  diag => cs%diag
2398 
2399  thickness_units = get_thickness_units(gv)
2400  flux_units = get_flux_units(gv)
2401  t_flux_units = get_tr_flux_units(gv, "Celsius")
2402  s_flux_units = get_tr_flux_units(gv, "PPT")
2403 
2404  !Initialize the diagnostics mask arrays.
2405  !This has to be done after MOM_initialize_state call.
2406  !call diag_masks_set(G, CS%missing)
2407 
2408  cs%id_u = register_diag_field('ocean_model', 'u', diag%axesCuL, time, &
2409  'Zonal velocity', 'meter second-1', cmor_field_name='uo', cmor_units='m s-1', &
2410  cmor_standard_name='sea_water_x_velocity', cmor_long_name='Sea Water X Velocity')
2411  cs%id_v = register_diag_field('ocean_model', 'v', diag%axesCvL, time, &
2412  'Meridional velocity', 'meter second-1', cmor_field_name='vo', cmor_units='m s-1', &
2413  cmor_standard_name='sea_water_y_velocity', cmor_long_name='Sea Water Y Velocity')
2414  cs%id_h = register_diag_field('ocean_model', 'h', diag%axesTL, time, &
2415  'Layer Thickness', thickness_units, v_extensive=.true.)
2416 
2417  cs%id_volo = register_scalar_field('ocean_model', 'volo', time, diag,&
2418  long_name='Total volume of liquid ocean', units='m3', &
2419  standard_name='sea_water_volume')
2420  cs%id_zos = register_diag_field('ocean_model', 'zos', diag%axesT1, time,&
2421  standard_name = 'sea_surface_height_above_geoid', &
2422  long_name= 'Sea surface height above geoid', units='meter', missing_value=cs%missing)
2423  cs%id_zossq = register_diag_field('ocean_model', 'zossq', diag%axesT1, time,&
2424  standard_name='square_of_sea_surface_height_above_geoid', &
2425  long_name='Square of sea surface height above geoid', units='m2', missing_value=cs%missing)
2426  cs%id_ssh = register_diag_field('ocean_model', 'SSH', diag%axesT1, time, &
2427  'Sea Surface Height', 'meter', cs%missing)
2428  cs%id_ssh_ga = register_scalar_field('ocean_model', 'ssh_ga', time, diag,&
2429  long_name='Area averaged sea surface height', units='m', &
2430  standard_name='area_averaged_sea_surface_height')
2431  cs%id_ssh_inst = register_diag_field('ocean_model', 'SSH_inst', diag%axesT1, time, &
2432  'Instantaneous Sea Surface Height', 'meter', cs%missing)
2433  cs%id_ssu = register_diag_field('ocean_model', 'SSU', diag%axesCu1, time, &
2434  'Sea Surface Zonal Velocity', 'meter second-1', cs%missing)
2435  cs%id_ssv = register_diag_field('ocean_model', 'SSV', diag%axesCv1, time, &
2436  'Sea Surface Meridional Velocity', 'meter second-1', cs%missing)
2437  cs%id_speed = register_diag_field('ocean_model', 'speed', diag%axesT1, time, &
2438  'Sea Surface Speed', 'meter second-1', cs%missing)
2439 
2440  if (cs%use_temperature) then
2441  cs%id_T = register_diag_field('ocean_model', 'temp', diag%axesTL, time, &
2442  'Potential Temperature', 'Celsius', &
2443  cmor_field_name="thetao", cmor_units="C", &
2444  cmor_standard_name="sea_water_potential_temperature", &
2445  cmor_long_name ="Sea Water Potential Temperature")
2446  cs%id_S = register_diag_field('ocean_model', 'salt', diag%axesTL, time, &
2447  long_name='Salinity', units='PPT', cmor_field_name='so', &
2448  cmor_long_name='Sea Water Salinity', cmor_units='ppt', &
2449  cmor_standard_name='sea_water_salinity')
2450  cs%id_tob = register_diag_field('ocean_model','tob', diag%axesT1, time, &
2451  long_name='Sea Water Potential Temperature at Sea Floor', &
2452  standard_name='sea_water_potential_temperature_at_sea_floor', units='degC')
2453  cs%id_sob = register_diag_field('ocean_model','sob',diag%axesT1, time, &
2454  long_name='Sea Water Salinity at Sea Floor', &
2455  standard_name='sea_water_salinity_at_sea_floor', units='ppt')
2456  cs%id_sst = register_diag_field('ocean_model', 'SST', diag%axesT1, time, &
2457  'Sea Surface Temperature', 'Celsius', cs%missing, cmor_field_name='tos', &
2458  cmor_long_name='Sea Surface Temperature', cmor_units='degC', &
2459  cmor_standard_name='sea_surface_temperature')
2460  cs%id_sst_sq = register_diag_field('ocean_model', 'SST_sq', diag%axesT1, time, &
2461  'Sea Surface Temperature Squared', 'Celsius**2', cs%missing, cmor_field_name='tossq', &
2462  cmor_long_name='Square of Sea Surface Temperature ', cmor_units='degC^2', &
2463  cmor_standard_name='square_of_sea_surface_temperature')
2464  cs%id_sss = register_diag_field('ocean_model', 'SSS', diag%axesT1, time, &
2465  'Sea Surface Salinity', 'PPT', cs%missing, cmor_field_name='sos', &
2466  cmor_long_name='Sea Surface Salinity', cmor_units='ppt', &
2467  cmor_standard_name='sea_surface_salinity')
2468  cs%id_sss_sq = register_diag_field('ocean_model', 'SSS_sq', diag%axesT1, time, &
2469  'Sea Surface Salinity Squared', 'ppt**2', cs%missing, cmor_field_name='sossq', &
2470  cmor_long_name='Square of Sea Surface Salinity ', cmor_units='ppt^2', &
2471  cmor_standard_name='square_of_sea_surface_salinity')
2472  cs%id_Tcon = register_diag_field('ocean_model', 'contemp', diag%axesTL, time, &
2473  'Conservative Temperature', 'Celsius')
2474  cs%id_Sabs = register_diag_field('ocean_model', 'abssalt', diag%axesTL, time, &
2475  long_name='Absolute Salinity', units='g/Kg')
2476  cs%id_sstcon = register_diag_field('ocean_model', 'conSST', diag%axesT1, time, &
2477  'Sea Surface Conservative Temperature', 'Celsius', cs%missing)
2478  cs%id_sssabs = register_diag_field('ocean_model', 'absSSS', diag%axesT1, time, &
2479  'Sea Surface Absolute Salinity', 'g/Kg', cs%missing)
2480 
2481  endif
2482 
2483  if (cs%use_temperature .and. cs%use_frazil) then
2484  cs%id_fraz = register_diag_field('ocean_model', 'frazil', diag%axesT1, time, &
2485  'Heat from frazil formation', 'Watt meter-2', cmor_field_name='hfsifrazil', &
2486  cmor_units='W m-2', cmor_standard_name='heat_flux_into_sea_water_due_to_frazil_ice_formation', &
2487  cmor_long_name='Heat Flux into Sea Water due to Frazil Ice Formation')
2488  endif
2489 
2490  cs%id_salt_deficit = register_diag_field('ocean_model', 'salt_deficit', diag%axesT1, time, &
2491  'Salt sink in ocean due to ice flux', 'g Salt meter-2 s-1')
2492  cs%id_Heat_PmE = register_diag_field('ocean_model', 'Heat_PmE', diag%axesT1, time, &
2493  'Heat flux into ocean from mass flux into ocean', 'Watt meter-2')
2494  cs%id_intern_heat = register_diag_field('ocean_model', 'internal_heat', diag%axesT1, time,&
2495  'Heat flux into ocean from geothermal or other internal sources', 'Watt meter-2')
2496 
2497 
2498  ! lateral heat advective and diffusive fluxes
2499  cs%id_Tadx = register_diag_field('ocean_model', 'T_adx', diag%axesCuL, time, &
2500  'Advective (by residual mean) Zonal Flux of Potential Temperature', t_flux_units)
2501  cs%id_Tady = register_diag_field('ocean_model', 'T_ady', diag%axesCvL, time, &
2502  'Advective (by residual mean) Meridional Flux of Potential Temperature', t_flux_units)
2503  cs%id_Tdiffx = register_diag_field('ocean_model', 'T_diffx', diag%axesCuL, time, &
2504  'Diffusive Zonal Flux of Potential Temperature', t_flux_units)
2505  cs%id_Tdiffy = register_diag_field('ocean_model', 'T_diffy', diag%axesCvL, time, &
2506  'Diffusive Meridional Flux of Potential Temperature', t_flux_units)
2507  if (cs%id_Tadx > 0) call safe_alloc_ptr(cs%T_adx,isdb,iedb,jsd,jed,nz)
2508  if (cs%id_Tady > 0) call safe_alloc_ptr(cs%T_ady,isd,ied,jsdb,jedb,nz)
2509  if (cs%id_Tdiffx > 0) call safe_alloc_ptr(cs%T_diffx,isdb,iedb,jsd,jed,nz)
2510  if (cs%id_Tdiffy > 0) call safe_alloc_ptr(cs%T_diffy,isd,ied,jsdb,jedb,nz)
2511 
2512 
2513  ! lateral salt advective and diffusive fluxes
2514  cs%id_Sadx = register_diag_field('ocean_model', 'S_adx', diag%axesCuL, time, &
2515  'Advective (by residual mean) Zonal Flux of Salinity', s_flux_units)
2516  cs%id_Sady = register_diag_field('ocean_model', 'S_ady', diag%axesCvL, time, &
2517  'Advective (by residual mean) Meridional Flux of Salinity', s_flux_units)
2518  cs%id_Sdiffx = register_diag_field('ocean_model', 'S_diffx', diag%axesCuL, time, &
2519  'Diffusive Zonal Flux of Salinity', s_flux_units)
2520  cs%id_Sdiffy = register_diag_field('ocean_model', 'S_diffy', diag%axesCvL, time, &
2521  'Diffusive Meridional Flux of Salinity', s_flux_units)
2522  if (cs%id_Sadx > 0) call safe_alloc_ptr(cs%S_adx,isdb,iedb,jsd,jed,nz)
2523  if (cs%id_Sady > 0) call safe_alloc_ptr(cs%S_ady,isd,ied,jsdb,jedb,nz)
2524  if (cs%id_Sdiffx > 0) call safe_alloc_ptr(cs%S_diffx,isdb,iedb,jsd,jed,nz)
2525  if (cs%id_Sdiffy > 0) call safe_alloc_ptr(cs%S_diffy,isd,ied,jsdb,jedb,nz)
2526 
2527 
2528  ! vertically integrated lateral heat advective and diffusive fluxes
2529  cs%id_Tadx_2d = register_diag_field('ocean_model', 'T_adx_2d', diag%axesCu1, time, &
2530  'Vertically Integrated Advective Zonal Flux of Potential Temperature', t_flux_units)
2531  cs%id_Tady_2d = register_diag_field('ocean_model', 'T_ady_2d', diag%axesCv1, time, &
2532  'Vertically Integrated Advective Meridional Flux of Potential Temperature', t_flux_units)
2533  cs%id_Tdiffx_2d = register_diag_field('ocean_model', 'T_diffx_2d', diag%axesCu1, time, &
2534  'Vertically Integrated Diffusive Zonal Flux of Potential Temperature', t_flux_units)
2535  cs%id_Tdiffy_2d = register_diag_field('ocean_model', 'T_diffy_2d', diag%axesCv1, time, &
2536  'Vertically Integrated Diffusive Meridional Flux of Potential Temperature', t_flux_units)
2537  if (cs%id_Tadx_2d > 0) call safe_alloc_ptr(cs%T_adx_2d,isdb,iedb,jsd,jed)
2538  if (cs%id_Tady_2d > 0) call safe_alloc_ptr(cs%T_ady_2d,isd,ied,jsdb,jedb)
2539  if (cs%id_Tdiffx_2d > 0) call safe_alloc_ptr(cs%T_diffx_2d,isdb,iedb,jsd,jed)
2540  if (cs%id_Tdiffy_2d > 0) call safe_alloc_ptr(cs%T_diffy_2d,isd,ied,jsdb,jedb)
2541 
2542  ! vertically integrated lateral salt advective and diffusive fluxes
2543  cs%id_Sadx_2d = register_diag_field('ocean_model', 'S_adx_2d', diag%axesCu1, time, &
2544  'Vertically Integrated Advective Zonal Flux of Salinity', s_flux_units)
2545  cs%id_Sady_2d = register_diag_field('ocean_model', 'S_ady_2d', diag%axesCv1, time, &
2546  'Vertically Integrated Advective Meridional Flux of Salinity', s_flux_units)
2547  cs%id_Sdiffx_2d = register_diag_field('ocean_model', 'S_diffx_2d', diag%axesCu1, time, &
2548  'Vertically Integrated Diffusive Zonal Flux of Salinity', s_flux_units)
2549  cs%id_Sdiffy_2d = register_diag_field('ocean_model', 'S_diffy_2d', diag%axesCv1, time, &
2550  'Vertically Integrated Diffusive Meridional Flux of Salinity', s_flux_units)
2551  if (cs%id_Sadx_2d > 0) call safe_alloc_ptr(cs%S_adx_2d,isdb,iedb,jsd,jed)
2552  if (cs%id_Sady_2d > 0) call safe_alloc_ptr(cs%S_ady_2d,isd,ied,jsdb,jedb)
2553  if (cs%id_Sdiffx_2d > 0) call safe_alloc_ptr(cs%S_diffx_2d,isdb,iedb,jsd,jed)
2554  if (cs%id_Sdiffy_2d > 0) call safe_alloc_ptr(cs%S_diffy_2d,isd,ied,jsdb,jedb)
2555 
2556 
2557  if (cs%debug_truncations) then
2558  call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)
2559  call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
2560  if (.not.cs%adiabatic) then
2561  call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)
2562  call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
2563  endif
2564  endif
2565 
2566  ! fields posted prior to dynamics step
2567  cs%id_h_pre_dyn = register_diag_field('ocean_model', 'h_pre_dyn', diag%axesTL, time, &
2568  'Layer Thickness before dynamics step', thickness_units)
2569 
2570  ! diagnostics for values prior to diabatic and prior to ALE
2571  cs%id_u_predia = register_diag_field('ocean_model', 'u_predia', diag%axesCuL, time, &
2572  'Zonal velocity before diabatic forcing', 'meter second-1')
2573  cs%id_v_predia = register_diag_field('ocean_model', 'v_predia', diag%axesCvL, time, &
2574  'Meridional velocity before diabatic forcing', 'meter second-1')
2575  cs%id_h_predia = register_diag_field('ocean_model', 'h_predia', diag%axesTL, time, &
2576  'Layer Thickness before diabatic forcing', thickness_units, v_extensive=.true.)
2577  cs%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, time, &
2578  'Interface Heights before diabatic forcing', 'meter')
2579  if (.not. cs%adiabatic) then
2580  cs%id_u_preale = register_diag_field('ocean_model', 'u_preale', diag%axesCuL, time, &
2581  'Zonal velocity before remapping', 'meter second-1')
2582  cs%id_v_preale = register_diag_field('ocean_model', 'v_preale', diag%axesCvL, time, &
2583  'Meridional velocity before remapping', 'meter second-1')
2584  cs%id_h_preale = register_diag_field('ocean_model', 'h_preale', diag%axesTL, time, &
2585  'Layer Thickness before remapping', thickness_units, v_extensive=.true.)
2586  cs%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, time, &
2587  'Temperature before remapping', 'degC')
2588  cs%id_S_preale = register_diag_field('ocean_model', 'S_preale', diag%axesTL, time, &
2589  'Salinity before remapping', 'ppt')
2590  cs%id_e_preale = register_diag_field('ocean_model', 'e_preale', diag%axesTi, time, &
2591  'Interface Heights before remapping', 'meter')
2592  endif
2593 
2594  if (cs%use_temperature) then
2595  cs%id_T_predia = register_diag_field('ocean_model', 'temp_predia', diag%axesTL, time, &
2596  'Potential Temperature', 'Celsius')
2597  cs%id_S_predia = register_diag_field('ocean_model', 'salt_predia', diag%axesTL, time, &
2598  'Salinity', 'PPT')
2599  endif
2600 
2601  ! Diagnostics related to tracer transport
2602  cs%id_uhtr = register_diag_field('ocean_model', 'uhtr', diag%axesCuL, time, &
2603  'Accumulated zonal thickness fluxes to advect tracers', 'kg', &
2604  y_cell_method='sum', v_extensive=.true.)
2605  cs%id_vhtr = register_diag_field('ocean_model', 'vhtr', diag%axesCvL, time, &
2606  'Accumulated meridional thickness fluxes to advect tracers', 'kg', &
2607  x_cell_method='sum', v_extensive=.true.)
2608  cs%id_umo = register_diag_field('ocean_model', 'umo', &
2609  diag%axesCuL, time, 'Ocean Mass X Transport', 'kg/s', &
2610  standard_name='ocean_mass_x_transport', y_cell_method='sum', v_extensive=.true.)
2611  cs%id_vmo = register_diag_field('ocean_model', 'vmo', &
2612  diag%axesCvL, time, 'Ocean Mass Y Transport', 'kg/s', &
2613  standard_name='ocean_mass_y_transport', x_cell_method='sum', v_extensive=.true.)
2614  cs%id_umo_2d = register_diag_field('ocean_model', 'umo_2d', &
2615  diag%axesCu1, time, 'Ocean Mass X Transport Vertical Sum', 'kg/s', &
2616  standard_name='ocean_mass_x_transport_vertical_sum', y_cell_method='sum')
2617  cs%id_vmo_2d = register_diag_field('ocean_model', 'vmo_2d', &
2618  diag%axesCv1, time, 'Ocean Mass Y Transport Vertical Sum', 'kg/s', &
2619  standard_name='ocean_mass_y_transport_vertical_sum', x_cell_method='sum')
2620 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ register_diags_ts_tendency()

subroutine mom::register_diags_ts_tendency ( type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(inout)  G,
type(mom_control_struct), pointer  CS 
)
private

Initialize diagnostics for temp/heat and saln/salt tendencies.

Parameters
[in]timecurrent model time
[in,out]gocean grid structure
cscontrol structure set up by initialize_MOM

Definition at line 2626 of file MOM.F90.

Referenced by initialize_mom().

2626  type(time_type), intent(in) :: time !< current model time
2627  type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
2628  type(mom_control_struct), pointer :: cs !< control structure set up by initialize_MOM
2629 
2630  type(diag_ctrl), pointer :: diag
2631  integer :: i, j, k
2632  integer :: isd, ied, jsd, jed, nz
2633  integer :: is, ie, js, je
2634 
2635  diag => cs%diag
2636  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
2637  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2638 
2639 
2640  ! heat tendencies from lateral advection
2641  cs%id_T_advection_xy = register_diag_field('ocean_model', 'T_advection_xy', diag%axesTL, time, &
2642  'Horizontal convergence of residual mean heat advective fluxes', 'W/m2',v_extensive=.true.)
2643  cs%id_T_advection_xy_2d = register_diag_field('ocean_model', 'T_advection_xy_2d', diag%axesT1, time,&
2644  'Vertical sum of horizontal convergence of residual mean heat advective fluxes', 'W/m2')
2645  if (cs%id_T_advection_xy > 0 .or. cs%id_T_advection_xy_2d > 0) then
2646  call safe_alloc_ptr(cs%T_advection_xy,isd,ied,jsd,jed,nz)
2647  cs%tendency_diagnostics = .true.
2648  endif
2649 
2650  ! net temperature and heat tendencies
2651  cs%id_T_tendency = register_diag_field('ocean_model', 'T_tendency', diag%axesTL, time, &
2652  'Net time tendency for temperature', 'degC/s')
2653  cs%id_Th_tendency = register_diag_field('ocean_model', 'Th_tendency', diag%axesTL, time, &
2654  'Net time tendency for heat', 'W/m2', &
2655  cmor_field_name="opottemptend", cmor_units="W m-2", &
2656  cmor_standard_name="tendency_of_sea_water_potential_temperature_expressed_as_heat_content", &
2657  cmor_long_name ="Tendency of Sea Water Potential Temperature Expressed as Heat Content", &
2658  v_extensive=.true.)
2659  cs%id_Th_tendency_2d = register_diag_field('ocean_model', 'Th_tendency_2d', diag%axesT1, time, &
2660  'Vertical sum of net time tendency for heat', 'W/m2', &
2661  cmor_field_name="opottemptend_2d", cmor_units="W m-2", &
2662  cmor_standard_name="tendency_of_sea_water_potential_temperature_expressed_as_heat_content_vertical_sum",&
2663  cmor_long_name ="Tendency of Sea Water Potential Temperature Expressed as Heat Content Vertical Sum")
2664  if (cs%id_T_tendency > 0) then
2665  cs%tendency_diagnostics = .true.
2666  call safe_alloc_ptr(cs%T_prev,isd,ied,jsd,jed,nz)
2667  do k=1,nz ; do j=js,je ; do i=is,ie
2668  cs%T_prev(i,j,k) = cs%tv%T(i,j,k)
2669  enddo ; enddo ; enddo
2670  endif
2671  if (cs%id_Th_tendency > 0 .or. cs%id_Th_tendency_2d > 0) then
2672  cs%tendency_diagnostics = .true.
2673  call safe_alloc_ptr(cs%Th_prev,isd,ied,jsd,jed,nz)
2674  do k=1,nz ; do j=js,je ; do i=is,ie
2675  cs%Th_prev(i,j,k) = cs%tv%T(i,j,k) * cs%h(i,j,k)
2676  enddo ; enddo ; enddo
2677  endif
2678 
2679 
2680  ! salt tendencies from lateral advection
2681  cs%id_S_advection_xy = register_diag_field('ocean_model', 'S_advection_xy', diag%axesTL, time, &
2682  'Horizontal convergence of residual mean salt advective fluxes', 'kg/(m2 * s)', v_extensive=.true.)
2683  cs%id_S_advection_xy_2d = register_diag_field('ocean_model', 'S_advection_xy_2d', diag%axesT1, time,&
2684  'Vertical sum of horizontal convergence of residual mean salt advective fluxes', 'kg/(m2 * s)')
2685  if (cs%id_S_advection_xy > 0 .or. cs%id_S_advection_xy_2d > 0) then
2686  call safe_alloc_ptr(cs%S_advection_xy,isd,ied,jsd,jed,nz)
2687  cs%tendency_diagnostics = .true.
2688  endif
2689 
2690  ! net salinity and salt tendencies
2691  cs%id_S_tendency = register_diag_field('ocean_model', 'S_tendency', diag%axesTL, time, &
2692  'Net time tendency for salinity', 'PPT/s')
2693  cs%id_Sh_tendency = register_diag_field('ocean_model', 'Sh_tendency', diag%axesTL, time,&
2694  'Net time tendency for salt', 'kg/(m2 * s)', &
2695  cmor_field_name="osalttend", cmor_units="kg m-2 s-1", &
2696  cmor_standard_name="tendency_of_sea_water_salinity_expressed_as_salt_content", &
2697  cmor_long_name ="Tendency of Sea Water Salinity Expressed as Salt Content", &
2698  v_extensive=.true.)
2699  cs%id_Sh_tendency_2d = register_diag_field('ocean_model', 'Sh_tendency_2d', diag%axesT1, time, &
2700  'Vertical sum of net time tendency for salt', 'kg/(m2 * s)', &
2701  cmor_field_name="osalttend_2d", cmor_units="kg m-2 s-1", &
2702  cmor_standard_name="tendency_of_sea_water_salinity_expressed_as_salt_content_vertical_sum",&
2703  cmor_long_name ="Tendency of Sea Water Salinity Expressed as Salt Content Vertical Sum")
2704  if (cs%id_S_tendency > 0) then
2705  cs%tendency_diagnostics = .true.
2706  call safe_alloc_ptr(cs%S_prev,isd,ied,jsd,jed,nz)
2707  do k=1,nz ; do j=js,je ; do i=is,ie
2708  cs%S_prev(i,j,k) = cs%tv%S(i,j,k)
2709  enddo ; enddo ; enddo
2710  endif
2711  if (cs%id_Sh_tendency > 0 .or. cs%id_Sh_tendency_2d > 0) then
2712  cs%tendency_diagnostics = .true.
2713  call safe_alloc_ptr(cs%Sh_prev,isd,ied,jsd,jed,nz)
2714  do k=1,nz ; do j=js,je ; do i=is,ie
2715  cs%Sh_prev(i,j,k) = cs%tv%S(i,j,k) * cs%h(i,j,k)
2716  enddo ; enddo ; enddo
2717  endif
2718 
Here is the caller graph for this function:

◆ register_diags_ts_vardec()

subroutine mom::register_diags_ts_vardec ( type(time_type), intent(in)  Time,
type(hor_index_type), intent(in)  HI,
type(verticalgrid_type), intent(in)  GV,
type(param_file_type), intent(in)  param_file,
type(mom_control_struct), pointer  CS 
)
private

Initialize diagnostics for the variance decay of temp/salt across regridding/remapping.

Parameters
[in]timecurrent model time
[in]hihorizontal index type
[in]gvocean vertical grid structure
[in]param_fileparameter file
cscontrol structure for MOM

Definition at line 2725 of file MOM.F90.

Referenced by initialize_mom().

2725  type(time_type), intent(in) :: time !< current model time
2726  type(hor_index_type), intent(in) :: hi !< horizontal index type
2727  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
2728  type(param_file_type), intent(in) :: param_file !< parameter file
2729  type(mom_control_struct), pointer :: cs !< control structure for MOM
2730 
2731  integer :: isd, ied, jsd, jed, nz
2732  type(vardesc) :: vd_tmp
2733  type(diag_ctrl), pointer :: diag
2734 
2735  diag => cs%diag
2736  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
2737 
2738  ! variancy decay through ALE operation
2739  cs%id_T_vardec = register_diag_field('ocean_model', 'T_vardec', diag%axesTL, time, &
2740  'ALE variance decay for temperature', 'degC2 s-1')
2741  if (cs%id_T_vardec > 0) then
2742  call safe_alloc_ptr(cs%T_squared,isd,ied,jsd,jed,nz)
2743  cs%T_squared(:,:,:) = 0.
2744 
2745  vd_tmp = var_desc(name="T2", units="degC2", longname="Squared Potential Temperature")
2746  call register_tracer(cs%T_squared, vd_tmp, param_file, hi, gv, cs%tracer_reg)
2747  endif
2748 
2749  cs%id_S_vardec = register_diag_field('ocean_model', 'S_vardec', diag%axesTL, time, &
2750  'ALE variance decay for salinity', 'PPT2 s-1')
2751  if (cs%id_S_vardec > 0) then
2752  call safe_alloc_ptr(cs%S_squared,isd,ied,jsd,jed,nz)
2753  cs%S_squared(:,:,:) = 0.
2754 
2755  vd_tmp = var_desc(name="S2", units="PPT2", longname="Squared Salinity")
2756  call register_tracer(cs%S_squared, vd_tmp, param_file, hi, gv, cs%tracer_reg)
2757  endif
2758 
Here is the caller graph for this function:

◆ set_restart_fields()

subroutine mom::set_restart_fields ( type(verticalgrid_type), intent(inout)  GV,
type(param_file_type), intent(in)  param_file,
type(mom_control_struct), intent(in)  CS 
)
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]param_fileopened file for parsing to get parameters
[in]cscontrol structure set up by inialize_MOM

Definition at line 3354 of file MOM.F90.

Referenced by initialize_mom().

3354  type(verticalgrid_type), intent(inout) :: gv !< ocean vertical grid structure
3355  type(param_file_type), intent(in) :: param_file !< opened file for parsing to get parameters
3356  type(mom_control_struct), intent(in) :: cs !< control structure set up by inialize_MOM
3357  ! Local variables
3358  logical :: use_ice_shelf ! Needed to determine whether to add CS%Hml to restarts
3359  type(vardesc) :: vd
3360  character(len=48) :: thickness_units, flux_units
3361 
3362  call get_param(param_file, '', "ICE_SHELF", use_ice_shelf, default=.false., do_not_log=.true.)
3363 
3364  thickness_units = get_thickness_units(gv)
3365  flux_units = get_flux_units(gv)
3366 
3367  if (cs%use_temperature) then
3368  vd = var_desc("Temp","degC","Potential Temperature")
3369  call register_restart_field(cs%tv%T, vd, .true., cs%restart_CSp)
3370 
3371  vd = var_desc("Salt","PPT","Salinity")
3372  call register_restart_field(cs%tv%S, vd, .true., cs%restart_CSp)
3373  endif
3374 
3375  vd = var_desc("h",thickness_units,"Layer Thickness")
3376  call register_restart_field(cs%h, vd, .true., cs%restart_CSp)
3377 
3378  vd = var_desc("u","meter second-1","Zonal velocity",'u','L')
3379  call register_restart_field(cs%u, vd, .true., cs%restart_CSp)
3380 
3381  vd = var_desc("v","meter second-1","Meridional velocity",'v','L')
3382  call register_restart_field(cs%v, vd, .true., cs%restart_CSp)
3383 
3384  if (cs%use_frazil) then
3385  vd = var_desc("frazil","J m-2","Frazil heat flux into ocean",'h','1')
3386  call register_restart_field(cs%tv%frazil, vd, .false., cs%restart_CSp)
3387  endif
3388 
3389  if (cs%interp_p_surf) then
3390  vd = var_desc("p_surf_prev","Pa","Previous ocean surface pressure",'h','1')
3391  call register_restart_field(cs%p_surf_prev, vd, .false., cs%restart_CSp)
3392  endif
3393 
3394  vd = var_desc("ave_ssh","meter","Time average sea surface height",'h','1')
3395  call register_restart_field(cs%ave_ssh, vd, .false., cs%restart_CSp)
3396 
3397  ! hML is needed when using the ice shelf module
3398  if (use_ice_shelf .and. associated(cs%Hml)) then
3399  vd = var_desc("hML","meter","Mixed layer thickness",'h','1')
3400  call register_restart_field(cs%Hml, vd, .false., cs%restart_CSp)
3401  endif
3402 
Here is the caller graph for this function:

◆ step_mom()

subroutine, public mom::step_mom ( type(forcing), intent(inout)  fluxes,
type(surface), intent(inout)  state,
type(time_type), intent(in)  Time_start,
real, intent(in)  time_interval,
type(mom_control_struct), pointer  CS 
)

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]fluxespointers to forcing fields
[in,out]statesurface ocean state
[in]time_startstarting time of a segment, as a time type
[in]time_intervaltime interval covered by this run segment, in s.
cscontrol structure from initialize_MOM

Definition at line 466 of file MOM.F90.

References adjust_ssh_for_p_atm(), mom_dynamics_legacy_split::adjustments_dyn_legacy_split(), mom_lateral_mixing_coeffs::calc_resoln_function(), calculate_surface_state(), mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), mom_error_handler::calltree_waypoint(), mom_domains::complete_group_pass(), id_clock_bbl_visc, id_clock_diagnostics, id_clock_dynamics, id_clock_ml_restrat, id_clock_ocean, id_clock_other, id_clock_pass, id_clock_thermo, id_clock_thick_diff, id_clock_tracer, id_clock_z_diag, mom_forcing_type::mom_forcing_chksum(), post_integrated_diagnostics(), post_surface_diagnostics(), post_transport_diagnostics(), post_ts_diagnostics(), mom_domains::start_group_pass(), mom_meke::step_forward_meke(), step_mom_thermo(), and mom_thickness_diffuse::thickness_diffuse().

Referenced by mom_main().

466  type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
467  type(surface), intent(inout) :: state !< surface ocean state
468  type(time_type), intent(in) :: time_start !< starting time of a segment, as a time type
469  real, intent(in) :: time_interval !< time interval covered by this run segment, in s.
470  type(mom_control_struct), pointer :: cs !< control structure from initialize_MOM
471 
472  ! local
473  type(ocean_grid_type), pointer :: g ! pointer to a structure containing
474  ! metrics and related information
475  type(verticalgrid_type), pointer :: gv => null()
476  integer, save :: nt_debug = 1 ! running number of iterations, for debugging only.
477  integer :: ntstep ! time steps between tracer updates or diabatic forcing
478  integer :: n_max ! number of steps to take in this call
479 
480  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, n
481  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
482 
483  real :: dt ! baroclinic time step (sec)
484  real :: dtth ! time step for thickness diffusion (sec)
485  real :: dtdia ! time step for diabatic processes (sec)
486  real :: dt_therm ! a limited and quantized version of CS%dt_therm (sec)
487  real :: dtbt_reset_time ! value of CS%rel_time when DTBT was last calculated (sec)
488 
489  real :: mass_src_time ! The amount of time for the surface mass source from
490  ! precipitation-evaporation, rivers, etc., that should
491  ! be applied to the start of the barotropic solver to
492  ! avoid generating tsunamis, in s. This is negative
493  ! if the precipation has already been applied to the
494  ! layers, and positive if it will be applied later.
495 
496  real :: wt_end, wt_beg
497  real :: bbl_time_int ! The amount of time over which the calculated BBL
498  ! properties will apply, for use in diagnostics.
499 
500  logical :: calc_dtbt ! Indicates whether the dynamically adjusted
501  ! barotropic time step needs to be updated.
502  logical :: do_advection ! If true, it is time to advect tracers.
503  logical :: do_calc_bbl ! If true, calculate the boundary layer properties.
504  logical :: thermo_does_span_coupling ! If true, thermodynamic forcing spans
505  ! multiple dynamic timesteps.
506  real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
507  eta_av, & ! average sea surface height or column mass over a timestep (meter or kg/m2)
508  ssh ! sea surface height based on eta_av (meter or kg/m2)
509 
510  real, pointer, dimension(:,:,:) :: &
511  u, & ! u : zonal velocity component (m/s)
512  v, & ! v : meridional velocity component (m/s)
513  h ! h : layer thickness (meter (Bouss) or kg/m2 (non-Bouss))
514 
515  ! Store the layer thicknesses before any changes by the dynamics. This is necessary for remapped
516  ! mass transport diagnostics
517  real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: h_pre_dyn
518  real :: tot_wt_ssh, itot_wt_ssh
519 
520  type(time_type) :: time_local
521  logical :: showcalltree
522  ! These are used for group halo passes.
523  logical :: do_pass_kv_turb, do_pass_ray, do_pass_kv_bbl_thick
524 
525  g => cs%G ; gv => cs%GV
526  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
527  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
528  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
529  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
530  u => cs%u ; v => cs%v ; h => cs%h
531 
532  call cpu_clock_begin(id_clock_ocean)
533  call cpu_clock_begin(id_clock_other)
534 
535  if (cs%debug) then
536  call mom_state_chksum("Beginning of step_MOM ", u, v, h, cs%uh, cs%vh, g, gv)
537  call hchksum(cs%h,"CS%h beginning of step_MOM",g%HI, scale=gv%H_to_m)
538  endif
539 
540  showcalltree = calltree_showquery()
541  if (showcalltree) call calltree_enter("step_MOM(), MOM.F90")
542 
543  ! First determine the time step that is consistent with this call.
544  ! It is anticipated that the time step will almost always coincide
545  ! with dt. In addition, ntstep is determined, subject to the constraint
546  ! that ntstep cannot exceed n_max.
547  if (time_interval <= cs%dt) then
548  n_max = 1
549  else
550  n_max = ceiling(time_interval/cs%dt - 0.001)
551  endif
552 
553  dt = time_interval / real(n_max)
554  dtdia = 0.0
555  thermo_does_span_coupling = (cs%thermo_spans_coupling .and. &
556  (cs%dt_therm > 1.5*time_interval))
557  if (thermo_does_span_coupling) then
558  ! Set dt_therm to be an integer multiple of the coupling time step.
559  dt_therm = time_interval * floor(cs%dt_therm / time_interval + 0.001)
560  ntstep = floor(dt_therm/dt + 0.001)
561  else
562  ntstep = max(1,min(n_max,floor(cs%dt_therm/dt + 0.001)))
563  dt_therm = dt*ntstep
564  endif
565 
566  if (.not.ASSOCIATED(fluxes%p_surf)) cs%interp_p_surf = .false.
567 
568  !---------- Begin setup for group halo pass
569 
570  call cpu_clock_begin(id_clock_pass)
571  call create_group_pass(cs%pass_tau_ustar_psurf, fluxes%taux, fluxes%tauy, g%Domain)
572  if (ASSOCIATED(fluxes%ustar)) &
573  call create_group_pass(cs%pass_tau_ustar_psurf, fluxes%ustar(:,:), g%Domain)
574  if (ASSOCIATED(fluxes%p_surf)) &
575  call create_group_pass(cs%pass_tau_ustar_psurf, fluxes%p_surf(:,:), g%Domain)
576 
577  do_pass_ray = .false.
578  if ((.not.g%Domain%symmetric) .and. &
579  associated(cs%visc%Ray_u) .and. associated(cs%visc%Ray_v)) then
580  call create_group_pass(cs%pass_ray, cs%visc%Ray_u, cs%visc%Ray_v, g%Domain, &
581  to_north+to_east+scalar_pair+omit_corners, cgrid_ne, halo=1)
582  do_pass_ray = .true.
583  endif
584  do_pass_kv_bbl_thick = .false.
585  if (associated(cs%visc%bbl_thick_u) .and. associated(cs%visc%bbl_thick_v)) then
586  call create_group_pass(cs%pass_bbl_thick_kv_bbl, cs%visc%bbl_thick_u, &
587  cs%visc%bbl_thick_v, g%Domain, &
588  to_north+to_east+scalar_pair+omit_corners, cgrid_ne, halo=1)
589  do_pass_kv_bbl_thick = .true.
590  endif
591  if (associated(cs%visc%kv_bbl_u) .and. associated(cs%visc%kv_bbl_v)) then
592  call create_group_pass(cs%pass_bbl_thick_kv_bbl, cs%visc%kv_bbl_u, &
593  cs%visc%kv_bbl_v, g%Domain, &
594  to_north+to_east+scalar_pair+omit_corners, cgrid_ne, halo=1)
595  do_pass_kv_bbl_thick = .true.
596  endif
597  do_pass_kv_turb = associated(cs%visc%Kv_turb)
598  if (associated(cs%visc%Kv_turb)) &
599  call create_group_pass(cs%pass_kv_turb, cs%visc%Kv_turb, g%Domain, to_all+omit_corners, halo=1)
600 
601  if (.not.cs%adiabatic .AND. cs%use_ALE_algorithm ) then
602  if (cs%use_temperature) then
603  call create_group_pass(cs%pass_T_S_h, cs%tv%T, g%Domain, to_all+omit_corners, halo=1)
604  call create_group_pass(cs%pass_T_S_h, cs%tv%S, g%Domain, to_all+omit_corners, halo=1)
605  endif
606  call create_group_pass(cs%pass_T_S_h, h, g%Domain, to_all+omit_corners, halo=1)
607  endif
608 
609  if ((cs%adiabatic .OR. cs%diabatic_first) .AND. cs%use_temperature) then
610  call create_group_pass(cs%pass_T_S, cs%tv%T, g%Domain, halo=1) ! Could also omit corners?
611  call create_group_pass(cs%pass_T_S, cs%tv%S, g%Domain, halo=1)
612  endif
613 
614  !---------- End setup for group halo pass
615 
616 
617  if (g%nonblocking_updates) then
618  call start_group_pass(cs%pass_tau_ustar_psurf, g%Domain)
619  else
620  call do_group_pass(cs%pass_tau_ustar_psurf, g%Domain)
621  endif
622  call cpu_clock_end(id_clock_pass)
623 
624  if (ASSOCIATED(cs%tv%frazil)) cs%tv%frazil(:,:) = 0.0
625  if (ASSOCIATED(cs%tv%salt_deficit)) cs%tv%salt_deficit(:,:) = 0.0
626  if (ASSOCIATED(cs%tv%TempxPmE)) cs%tv%TempxPmE(:,:) = 0.0
627  if (ASSOCIATED(cs%tv%internal_heat)) cs%tv%internal_heat(:,:) = 0.0
628 
629  cs%rel_time = 0.0
630 
631  tot_wt_ssh = 0.0
632  do j=js,je ; do i=is,ie ; cs%ave_ssh(i,j) = 0.0 ; ssh(i,j) = cs%missing; enddo ; enddo
633 
634  if (associated(cs%VarMix)) then
635  call enable_averaging(time_interval, time_start+set_time(int(time_interval)), &
636  cs%diag)
637  call calc_resoln_function(h, cs%tv, g, gv, cs%VarMix)
638  call disable_averaging(cs%diag)
639  endif
640 
641  if (g%nonblocking_updates) then
642  call cpu_clock_begin(id_clock_pass)
643  call complete_group_pass(cs%pass_tau_ustar_psurf, g%Domain)
644  call cpu_clock_end(id_clock_pass)
645  endif
646 
647  if (cs%interp_p_surf) then
648  if (.not.ASSOCIATED(cs%p_surf_end)) allocate(cs%p_surf_end(isd:ied,jsd:jed))
649  if (.not.ASSOCIATED(cs%p_surf_begin)) allocate(cs%p_surf_begin(isd:ied,jsd:jed))
650  if (.not.cs%p_surf_prev_set) then
651  do j=jsd,jed ; do i=isd,ied
652  cs%p_surf_prev(i,j) = fluxes%p_surf(i,j)
653  enddo ; enddo
654  cs%p_surf_prev_set = .true.
655  endif
656  else
657  cs%p_surf_end => fluxes%p_surf
658  endif
659 
660  if (cs%debug) then
661  call mom_state_chksum("Before steps ", u, v, h, cs%uh, cs%vh, g, gv)
662  call mom_forcing_chksum("Before steps", fluxes, g, haloshift=0)
663  call check_redundant("Before steps ", u, v, g)
664  call check_redundant("Before steps ", fluxes%taux, fluxes%tauy, g)
665  endif
666  call cpu_clock_end(id_clock_other)
667 
668  do n=1,n_max
669 
670  nt_debug = nt_debug + 1
671 
672  ! Set the universally visible time to the middle of the time step
673  cs%Time = time_start + set_time(int(floor(cs%rel_time+0.5*dt+0.5)))
674  cs%rel_time = cs%rel_time + dt
675 
676  ! Set the local time to the end of the time step.
677  time_local = time_start + set_time(int(floor(cs%rel_time+0.5)))
678  if (showcalltree) call calltree_enter("DT cycles (step_MOM) n=",n)
679 
680  !===========================================================================
681  ! This is the first place where the diabatic processes and remapping could occur.
682  if (cs%diabatic_first .and. (cs%t_dyn_rel_adv==0.0)) then ! do thermodynamics.
683 
684  if (thermo_does_span_coupling) then
685  dtdia = dt_therm
686  if ((fluxes%dt_buoy_accum > 0.0) .and. (dtdia > time_interval) .and. &
687  (abs(fluxes%dt_buoy_accum - dtdia) > 1e-6*dtdia)) then
688  call mom_error(fatal, "step_MOM: Mismatch between long thermodynamic "//&
689  "timestep and time over which buoyancy fluxes have been accumulated.")
690  endif
691  call mom_error(fatal, "MOM is not yet set up to have restarts that work "//&
692  "with THERMO_SPANS_COUPLING and DIABATIC_FIRST.")
693  else
694  dtdia = dt*min(ntstep,n_max-(n-1))
695  endif
696 
697  ! The end-time of the diagnostic interval needs to be set ahead if there
698  ! are multiple dynamic time steps worth of thermodynamics applied here.
699  call enable_averaging(dtdia, time_local + &
700  set_time(int(floor(dtdia-dt+0.5))), cs%diag)
701 
702  if (cs%debug) then
703  call uvchksum("Pre set_viscous_BBL [uv]", u, v, g%HI, haloshift=1)
704  call hchksum(h,"Pre set_viscous_BBL h", g%HI, haloshift=1, scale=gv%H_to_m)
705  if (associated(cs%tv%T)) call hchksum(cs%tv%T, "Pre set_viscous_BBL T", g%HI, haloshift=1)
706  if (associated(cs%tv%S)) call hchksum(cs%tv%S, "Pre set_viscous_BBL S", g%HI, haloshift=1)
707  endif
708 
709  ! Calculate the BBL properties and store them inside visc (u,h).
710  ! This is here so that CS%visc is updated before diabatic() when
711  ! DIABATIC_FIRST=True. Otherwise diabatic() is called after the dynamics
712  ! and set_viscous_BBL is called as a part of the dynamic stepping.
713  call cpu_clock_begin(id_clock_bbl_visc)
714  call set_viscous_bbl(u, v, h, cs%tv, cs%visc, g, gv, cs%set_visc_CSp)
715  call cpu_clock_end(id_clock_bbl_visc)
716 
717  call cpu_clock_begin(id_clock_pass)
718  if (do_pass_ray) call do_group_pass(cs%pass_ray, g%Domain )
719  if (do_pass_kv_bbl_thick) call do_group_pass(cs%pass_bbl_thick_kv_bbl, g%Domain)
720  call cpu_clock_end(id_clock_pass)
721  if (showcalltree) call calltree_waypoint("done with set_viscous_BBL (diabatic_first)")
722 
723  call cpu_clock_begin(id_clock_thermo)
724 
725  ! These adjustments are a part of an archaic time stepping algorithm.
726  if (cs%split .and. cs%legacy_split .and. .not.cs%adiabatic) then
727  call adjustments_dyn_legacy_split(u, v, h, dt, g, gv, cs%dyn_legacy_split_CSp)
728  endif
729 
730  ! Apply diabatic forcing, do mixing, and regrid.
731  call step_mom_thermo(cs, g, gv, u, v, h, cs%tv, fluxes, dtdia)
732 
733  ! The diabatic processes are now ahead of the dynamics by dtdia.
734  cs%t_dyn_rel_thermo = -dtdia
735  if (showcalltree) call calltree_waypoint("finished diabatic_first (step_MOM)")
736 
737  call disable_averaging(cs%diag)
738  call cpu_clock_end(id_clock_thermo)
739 
740  endif ! end of block "(CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0))"
741 
742  !===========================================================================
743  ! This is the start of the dynamics stepping part of the algorithm.
744 
745  call cpu_clock_begin(id_clock_dynamics)
746  call disable_averaging(cs%diag)
747 
748  if (cs%thickness_diffuse .and. cs%thickness_diffuse_first) then
749  if (cs%t_dyn_rel_adv == 0.0) then
750  if (thermo_does_span_coupling) then
751  dtth = dt_therm
752  else
753  dtth = dt*min(ntstep,n_max-n+1)
754  endif
755 
756  call enable_averaging(dtth,time_local+set_time(int(floor(dtth-dt+0.5))), cs%diag)
757  call cpu_clock_begin(id_clock_thick_diff)
758  if (associated(cs%VarMix)) &
759  call calc_slope_functions(h, cs%tv, dt, g, gv, cs%VarMix)
760  call thickness_diffuse(h, cs%uhtr, cs%vhtr, cs%tv, dtth, g, gv, &
761  cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp)
762  call cpu_clock_end(id_clock_thick_diff)
763  call cpu_clock_begin(id_clock_pass)
764  call pass_var(h, g%Domain) !###, halo=max(2,cont_stensil))
765  call cpu_clock_end(id_clock_pass)
766  call disable_averaging(cs%diag)
767  if (showcalltree) call calltree_waypoint("finished thickness_diffuse_first (step_MOM)")
768 
769  ! Whenever thickness changes let the diag manager know, target grids
770  ! for vertical remapping may need to be regenerated.
771  call diag_update_remap_grids(cs%diag)
772 
773  endif
774  endif
775 
776  ! The bottom boundary layer properties are out-of-date and need to be
777  ! recalculated. This always occurs at the start of a coupling time
778  ! step because the externally prescribed stresses may have changed.
779  do_calc_bbl = ((cs%t_dyn_rel_adv == 0.0) .or. (n==1))
780  if (do_calc_bbl) then
781  ! Calculate the BBL properties and store them inside visc (u,h).
782  call cpu_clock_begin(id_clock_bbl_visc)
783  bbl_time_int = max(dt, min(dt_therm - cs%t_dyn_rel_adv, dt*(1+n_max-n)) )
784  call enable_averaging(bbl_time_int, &
785  time_local+set_time(int(bbl_time_int-dt+0.5)), cs%diag)
786  call set_viscous_bbl(u, v, h, cs%tv, cs%visc, g, gv, cs%set_visc_CSp)
787  call disable_averaging(cs%diag)
788  call cpu_clock_end(id_clock_bbl_visc)
789  if (showcalltree) call calltree_waypoint("done with set_viscous_BBL (step_MOM)")
790  endif
791 
792  call cpu_clock_begin(id_clock_pass)
793  if (do_pass_kv_turb) call do_group_pass(cs%pass_kv_turb, g%Domain)
794  call cpu_clock_end(id_clock_pass)
795 
796  if (do_calc_bbl) then
797  call cpu_clock_begin(id_clock_pass)
798  if (g%nonblocking_updates) then
799  if (do_pass_ray) call start_group_pass(cs%pass_Ray, g%Domain)
800  if (do_pass_kv_bbl_thick) call start_group_pass(cs%pass_bbl_thick_kv_bbl, g%Domain)
801  ! do_calc_bbl will be set to .false. when the message passing is complete.
802  else
803  if (do_pass_ray) call do_group_pass(cs%pass_Ray, g%Domain)
804  if (do_pass_kv_bbl_thick) call do_group_pass(cs%pass_bbl_thick_kv_bbl, g%Domain)
805  endif
806  call cpu_clock_end(id_clock_pass)
807  endif
808 
809  if (cs%interp_p_surf) then
810  wt_end = real(n) / real(n_max)
811  wt_beg = real(n-1) / real(n_max)
812  do j=jsd,jed ; do i=isd,ied
813  cs%p_surf_end(i,j) = wt_end * fluxes%p_surf(i,j) + &
814  (1.0-wt_end) * cs%p_surf_prev(i,j)
815  cs%p_surf_begin(i,j) = wt_beg * fluxes%p_surf(i,j) + &
816  (1.0-wt_beg) * cs%p_surf_prev(i,j)
817  enddo ; enddo
818  endif
819 
820  if (associated(cs%u_prev) .and. associated(cs%v_prev)) then
821  do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb
822  cs%u_prev(i,j,k) = u(i,j,k)
823  enddo ; enddo ; enddo
824  do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
825  cs%v_prev(i,j,k) = u(i,j,k)
826  enddo ; enddo ; enddo
827  endif
828 
829  ! Store pre-dynamics layer thicknesses so that mass fluxes are remapped correctly
830  do k=1,nz ; do j=jsd,jed ; do i=isd,ied
831  h_pre_dyn(i,j,k) = h(i,j,k)
832  enddo ; enddo ; enddo
833 
834  if (g%nonblocking_updates) then ; if (do_calc_bbl) then
835  call cpu_clock_begin(id_clock_pass)
836  if (do_pass_ray) call complete_group_pass(cs%pass_Ray, g%Domain)
837  if (do_pass_kv_bbl_thick) call complete_group_pass(cs%pass_bbl_thick_kv_bbl, g%Domain)
838  call cpu_clock_end(id_clock_pass)
839  endif ; endif
840 
841  if (cs%do_dynamics .and. cs%split) then !--------------------------- start SPLIT
842  ! This section uses a split time stepping scheme for the dynamic equations,
843  ! basically the stacked shallow water equations with viscosity.
844 
845  calc_dtbt = .false.
846  if ((cs%dtbt_reset_period >= 0.0) .and. &
847  ((n==1) .or. (cs%dtbt_reset_period == 0.0) .or. &
848  (cs%rel_time >= dtbt_reset_time + 0.999*cs%dtbt_reset_period))) then
849  calc_dtbt = .true.
850  dtbt_reset_time = cs%rel_time
851  endif
852 
853  mass_src_time = cs%t_dyn_rel_thermo
854  if (cs%legacy_split) then
855  call step_mom_dyn_legacy_split(u, v, h, cs%tv, cs%visc, &
856  time_local, dt, fluxes, cs%p_surf_begin, cs%p_surf_end, &
857  mass_src_time, dt_therm, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
858  eta_av, g, gv, cs%dyn_legacy_split_CSp, calc_dtbt, cs%VarMix, cs%MEKE)
859  else
860  call step_mom_dyn_split_rk2(u, v, h, cs%tv, cs%visc, &
861  time_local, dt, fluxes, cs%p_surf_begin, cs%p_surf_end, &
862  mass_src_time, dt_therm, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
863  eta_av, g, gv, cs%dyn_split_RK2_CSp, calc_dtbt, cs%VarMix, cs%MEKE)
864  endif
865  if (showcalltree) call calltree_waypoint("finished step_MOM_dyn_split (step_MOM)")
866 
867 
868  elseif (cs%do_dynamics) then ! --------------------------------------------------- not SPLIT
869  ! This section uses an unsplit stepping scheme for the dynamic
870  ! equations; basically the stacked shallow water equations with viscosity.
871  ! Because the time step is limited by CFL restrictions on the external
872  ! gravity waves, the unsplit is usually much less efficient that the split
873  ! approaches. But because of its simplicity, the unsplit method is very
874  ! useful for debugging purposes.
875 
876  if (cs%use_RK2) then
877  call step_mom_dyn_unsplit_rk2(u, v, h, cs%tv, cs%visc, time_local, dt, fluxes, &
878  cs%p_surf_begin, cs%p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
879  eta_av, g, gv, cs%dyn_unsplit_RK2_CSp, cs%VarMix, cs%MEKE)
880  else
881  call step_mom_dyn_unsplit(u, v, h, cs%tv, cs%visc, time_local, dt, fluxes, &
882  cs%p_surf_begin, cs%p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
883  eta_av, g, gv, cs%dyn_unsplit_CSp, cs%VarMix, cs%MEKE)
884  endif
885  if (showcalltree) call calltree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)")
886 
887  endif ! -------------------------------------------------- end SPLIT
888 
889 
890 
891  if (cs%thickness_diffuse .and. .not.cs%thickness_diffuse_first) then
892  call cpu_clock_begin(id_clock_thick_diff)
893 
894  if (cs%debug) call hchksum(h,"Pre-thickness_diffuse h", g%HI, haloshift=0, scale=gv%H_to_m)
895 
896  if (associated(cs%VarMix)) &
897  call calc_slope_functions(h, cs%tv, dt, g, gv, cs%VarMix)
898  call thickness_diffuse(h, cs%uhtr, cs%vhtr, cs%tv, dt, g, gv, &
899  cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp)
900 
901  if (cs%debug) call hchksum(h,"Post-thickness_diffuse h", g%HI, haloshift=1, scale=gv%H_to_m)
902  call cpu_clock_end(id_clock_thick_diff)
903  call cpu_clock_begin(id_clock_pass)
904  call pass_var(h, g%Domain) !###, halo=max(2,cont_stensil))
905  call cpu_clock_end(id_clock_pass)
906  if (showcalltree) call calltree_waypoint("finished thickness_diffuse (step_MOM)")
907  endif
908 
909  ! apply the submesoscale mixed layer restratification parameterization
910  if (cs%mixedlayer_restrat) then
911  if (cs%debug) then
912  call hchksum(h,"Pre-mixedlayer_restrat h", g%HI, haloshift=1, scale=gv%H_to_m)
913  call uvchksum("Pre-mixedlayer_restrat uhtr", &
914  cs%uhtr, cs%vhtr, g%HI, haloshift=0)
915  endif
916  call cpu_clock_begin(id_clock_ml_restrat)
917  call mixedlayer_restrat(h, cs%uhtr ,cs%vhtr, cs%tv, fluxes, dt, cs%visc%MLD, &
918  cs%VarMix, g, gv, cs%mixedlayer_restrat_CSp)
919  call cpu_clock_end(id_clock_ml_restrat)
920  call cpu_clock_begin(id_clock_pass)
921  call pass_var(h, g%Domain) !###, halo=max(2,cont_stensil))
922  call cpu_clock_end(id_clock_pass)
923  if (cs%debug) then
924  call hchksum(h,"Post-mixedlayer_restrat h", g%HI, haloshift=1, scale=gv%H_to_m)
925  call uvchksum("Post-mixedlayer_restrat [uv]htr", &
926  cs%uhtr, cs%vhtr, g%HI, haloshift=0)
927  endif
928  endif
929 
930  ! Whenever thickness changes let the diag manager know, target grids
931  ! for vertical remapping may need to be regenerated.
932  call diag_update_remap_grids(cs%diag)
933 
934  if (cs%useMEKE) call step_forward_meke(cs%MEKE, h, cs%VarMix%SN_u, cs%VarMix%SN_v, &
935  cs%visc, dt, g, gv, cs%MEKE_CSp, cs%uhtr, cs%vhtr)
936  call disable_averaging(cs%diag)
937 
938  ! Advance the dynamics time by dt.
939  cs%t_dyn_rel_adv = cs%t_dyn_rel_adv + dt
940  cs%t_dyn_rel_thermo = cs%t_dyn_rel_thermo + dt
941  cs%t_dyn_rel_diag = cs%t_dyn_rel_diag + dt
942 
943  call cpu_clock_end(id_clock_dynamics)
944 
945  !===========================================================================
946  ! This is the start of the tracer advection part of the algorithm.
947 
948  if (thermo_does_span_coupling) then
949  do_advection = (cs%t_dyn_rel_adv + 0.5*dt > dt_therm)
950  else
951  do_advection = ((mod(n,ntstep) == 0) .or. (n==n_max))
952  endif
953 
954  if (do_advection) then ! Do advective transport and lateral tracer mixing.
955 
956  if (cs%debug) then
957  call cpu_clock_begin(id_clock_other)
958  call uvchksum("Pre-advection [uv]", u, v, g%HI, haloshift=2)
959  call hchksum(h,"Pre-advection h", g%HI, haloshift=1, scale=gv%H_to_m)
960  call uvchksum("Pre-advection uhtr", cs%uhtr, cs%vhtr, g%HI, &
961  haloshift=0, scale=gv%H_to_m)
962  ! call MOM_state_chksum("Pre-advection ", u, v, &
963  ! h, CS%uhtr, CS%vhtr, G, GV, haloshift=1)
964  if (associated(cs%tv%T)) call hchksum(cs%tv%T, "Pre-advection T", g%HI, haloshift=1)
965  if (associated(cs%tv%S)) call hchksum(cs%tv%S, "Pre-advection S", g%HI, haloshift=1)
966  if (associated(cs%tv%frazil)) call hchksum(cs%tv%frazil, &
967  "Pre-advection frazil", g%HI, haloshift=0)
968  if (associated(cs%tv%salt_deficit)) call hchksum(cs%tv%salt_deficit, &
969  "Pre-advection salt deficit", g%HI, haloshift=0)
970  ! call MOM_thermo_chksum("Pre-advection ", CS%tv, G)
971  call check_redundant("Pre-advection ", u, v, g)
972  call cpu_clock_end(id_clock_other)
973  endif
974 
975  call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer)
976  call enable_averaging(cs%t_dyn_rel_adv, time_local, cs%diag)
977 
978  call advect_tracer(h, cs%uhtr, cs%vhtr, cs%OBC, cs%t_dyn_rel_adv, g, gv, &
979  cs%tracer_adv_CSp, cs%tracer_Reg)
980  call tracer_hordiff(h, cs%t_dyn_rel_adv, cs%MEKE, cs%VarMix, g, gv, &
981  cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
982  if (showcalltree) call calltree_waypoint("finished tracer advection/diffusion (step_MOM)")
983  call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo)
984 
985  call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
986  call post_transport_diagnostics(g, gv, cs, cs%diag, cs%t_dyn_rel_adv, h, h_pre_dyn)
987  ! Rebuild the remap grids now that we've posted the fields which rely on thicknesses
988  ! from before the dynamics calls
989  call diag_update_remap_grids(cs%diag)
990 
991  call disable_averaging(cs%diag)
992  call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
993 
994  ! Reset the accumulated transports to 0 and record that the dynamics
995  ! and advective times now agree.
996  call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer)
997  cs%uhtr(:,:,:) = 0.0
998  cs%vhtr(:,:,:) = 0.0
999  cs%t_dyn_rel_adv = 0.0
1000  call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo)
1001 
1002  endif
1003 
1004  !===========================================================================
1005  ! This is the second place where the diabatic processes and remapping could occur.
1006  if (cs%t_dyn_rel_adv == 0.0) then
1007  call cpu_clock_begin(id_clock_thermo)
1008 
1009  if (.not.cs%diabatic_first) then
1010  dtdia = cs%t_dyn_rel_thermo
1011  if (thermo_does_span_coupling .and. (abs(dt_therm - dtdia) > 1e-6*dt_therm)) then
1012  call mom_error(fatal, "step_MOM: Mismatch between dt_therm and dtdia "//&
1013  "before call to diabatic.")
1014  endif
1015 
1016  call enable_averaging(cs%t_dyn_rel_thermo, time_local, cs%diag)
1017 
1018  ! These adjustments are a part of an archaic time stepping algorithm.
1019  if (cs%split .and. cs%legacy_split .and. .not.cs%adiabatic) then
1020  call adjustments_dyn_legacy_split(u, v, h, dt, g, gv, cs%dyn_legacy_split_CSp)
1021  endif
1022 
1023  ! Apply diabatic forcing, do mixing, and regrid.
1024  call step_mom_thermo(cs, g, gv, u, v, h, cs%tv, fluxes, dtdia)
1025 
1026  call disable_averaging(cs%diag)
1027 
1028  else ! "else branch for if (.not.CS%diabatic_first) then"
1029  if (abs(cs%t_dyn_rel_thermo) > 1e-6*dt) call mom_error(fatal, &
1030  "step_MOM: Mismatch between the dynamics and diabatic times "//&
1031  "with DIABATIC_FIRST.")
1032 
1033  ! Tracers have been advected and diffused, and need halo updates.
1034  if (cs%use_temperature) then
1035  call cpu_clock_begin(id_clock_pass)
1036  call do_group_pass(cs%pass_T_S, g%Domain)
1037  call cpu_clock_end(id_clock_pass)
1038  endif
1039  endif ! close of "if (.not.CS%diabatic_first) then ; if (.not.CS%adiabatic)"
1040 
1041  ! Record that the dynamics and diabatic processes are synchronized.
1042  cs%t_dyn_rel_thermo = 0.0
1043  call cpu_clock_end(id_clock_thermo)
1044  endif
1045 
1046  call cpu_clock_begin(id_clock_dynamics)
1047 
1048  ! Determining the time-average sea surface height is part of the algorithm.
1049  ! This may be eta_av if Boussinesq, or need to be diagnosed if not.
1050  tot_wt_ssh = tot_wt_ssh + dt
1051  call find_eta(h, cs%tv, gv%g_Earth, g, gv, ssh, eta_av)
1052  do j=js,je ; do i=is,ie
1053  cs%ave_ssh(i,j) = cs%ave_ssh(i,j) + dt*ssh(i,j)
1054  enddo ; enddo
1055  call cpu_clock_end(id_clock_dynamics)
1056 
1057  !===========================================================================
1058  ! Calculate diagnostics at the end of the time step.
1059  call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
1060 
1061  call enable_averaging(dt,time_local, cs%diag)
1062  ! These diagnostics are available every time step.
1063  if (cs%id_u > 0) call post_data(cs%id_u, u, cs%diag)
1064  if (cs%id_v > 0) call post_data(cs%id_v, v, cs%diag)
1065  if (cs%id_h > 0) call post_data(cs%id_h, h, cs%diag)
1066  if (cs%id_ssh_inst > 0) call post_data(cs%id_ssh_inst, ssh, cs%diag)
1067  call disable_averaging(cs%diag)
1068 
1069  if (cs%t_dyn_rel_adv == 0.0) then
1070  ! Diagnostics that require the complete state to be up-to-date can be calculated.
1071 
1072  call enable_averaging(cs%t_dyn_rel_diag, time_local, cs%diag)
1073  call calculate_diagnostic_fields(u, v, h, cs%uh, cs%vh, cs%tv, cs%ADp, &
1074  cs%CDp, fluxes, cs%t_dyn_rel_diag, g, gv, cs%diagnostics_CSp)
1075  call post_ts_diagnostics(cs, g, gv, cs%tv, cs%diag, cs%t_dyn_rel_diag)
1076  if (showcalltree) call calltree_waypoint("finished calculate_diagnostic_fields (step_MOM)")
1077  call disable_averaging(cs%diag)
1078  cs%t_dyn_rel_diag = 0.0
1079 
1080  call cpu_clock_begin(id_clock_z_diag)
1081  if (time_local + set_time(int(0.5*dt_therm)) > cs%Z_diag_time) then
1082  call enable_averaging(real(time_type_to_real(CS%Z_diag_interval)), &
1083  cs%z_diag_time, cs%diag)
1084  call calculate_z_diag_fields(u, v, h, ssh, fluxes%frac_shelf_h, &
1085  g, gv, cs%diag_to_Z_CSp)
1086  cs%Z_diag_time = cs%Z_diag_time + cs%Z_diag_interval
1087  call disable_averaging(cs%diag)
1088  if (showcalltree) call calltree_waypoint("finished calculate_Z_diag_fields (step_MOM)")
1089  endif
1090  call cpu_clock_end(id_clock_z_diag)
1091  endif
1092  call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)
1093 
1094  if (showcalltree) call calltree_leave("DT cycles (step_MOM)")
1095 
1096  enddo ! complete the n loop
1097 
1098  call cpu_clock_begin(id_clock_other)
1099 
1100  itot_wt_ssh = 1.0/tot_wt_ssh
1101  do j=js,je ; do i=is,ie
1102  cs%ave_ssh(i,j) = cs%ave_ssh(i,j)*itot_wt_ssh
1103  enddo ; enddo
1104 
1105  call enable_averaging(dt*n_max, time_local, cs%diag)
1106  call post_integrated_diagnostics(cs, g, gv, cs%diag, dt*n_max, cs%tv, fluxes)
1107  call disable_averaging(cs%diag)
1108 
1109  if (showcalltree) call calltree_waypoint("calling calculate_surface_state (step_MOM)")
1110  call adjust_ssh_for_p_atm(cs, g, gv, cs%ave_ssh, fluxes%p_surf_SSH)
1111  call calculate_surface_state(state, u, v, h, cs%ave_ssh, g, gv, cs)
1112 
1113  call enable_averaging(dt*n_max, time_local, cs%diag)
1114  call post_surface_diagnostics(cs, g, cs%diag, state)
1115  call disable_averaging(cs%diag)
1116 
1117  if (cs%interp_p_surf) then ; do j=jsd,jed ; do i=isd,ied
1118  cs%p_surf_prev(i,j) = fluxes%p_surf(i,j)
1119  enddo ; enddo ; endif
1120 
1121  call cpu_clock_end(id_clock_other)
1122 
1123  if (showcalltree) call calltree_leave("step_MOM()")
1124  call cpu_clock_end(id_clock_ocean)
1125 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ step_mom_thermo()

subroutine mom::step_mom_thermo ( type(mom_control_struct), intent(inout)  CS,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(inout)  GV,
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 
)
private

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

Parameters
[in,out]cscontrol structure
[in,out]gocean grid structure
[in,out]gvocean vertical grid structure
[in,out]uzonal velocity (m/s)
[in,out]vmeridional velocity (m/s)
[in,out]hlayer thickness (m or kg/m2)
[in,out]tvA structure pointing to various thermodynamic variables
[in,out]fluxespointers to forcing fields
[in]dtdiaThe time interval over which to advance, in s

Definition at line 1131 of file MOM.F90.

References mom_diabatic_driver::adiabatic(), mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), mom_error_handler::calltree_waypoint(), id_clock_ale, id_clock_diabatic, id_clock_pass, mom_forcing_type::mom_forcing_chksum(), mom_checksum_packages::mom_thermo_chksum(), and post_diags_ts_vardec().

Referenced by step_mom().

1131  type(mom_control_struct), intent(inout) :: cs !< control structure
1132  type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
1133  type(verticalgrid_type), intent(inout) :: gv !< ocean vertical grid structure
1134  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1135  intent(inout) :: u !< zonal velocity (m/s)
1136  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1137  intent(inout) :: v !< meridional velocity (m/s)
1138  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1139  intent(inout) :: h !< layer thickness (m or kg/m2)
1140  type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
1141  type(forcing), intent(inout) :: fluxes !< pointers to forcing fields
1142  real, intent(in) :: dtdia !< The time interval over which to advance, in s
1143 
1144  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta_predia, eta_preale
1145  integer :: i, j, k, is, ie, js, je, nz! , Isq, Ieq, Jsq, Jeq, n
1146  logical :: use_ice_shelf ! Needed for selecting the right ALE interface.
1147  logical :: showcalltree
1148 
1149  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1150  showcalltree = calltree_showquery()
1151  if (showcalltree) call calltree_enter("step_MOM_thermo(), MOM.F90")
1152 
1153  use_ice_shelf = .false.
1154  if (associated(fluxes%frac_shelf_h)) use_ice_shelf = .true.
1155 
1156  if (.not.cs%adiabatic) then
1157 
1158  if (cs%debug) then
1159  call uvchksum("Pre-diabatic [uv]", u, v, g%HI, haloshift=2)
1160  call hchksum(h,"Pre-diabatic h", g%HI, haloshift=1, scale=gv%H_to_m)
1161  call uvchksum("Pre-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1162  haloshift=0, scale=gv%H_to_m)
1163  ! call MOM_state_chksum("Pre-diabatic ",u, v, h, CS%uhtr, CS%vhtr, G, GV)
1164  call mom_thermo_chksum("Pre-diabatic ", cs%tv, g,haloshift=0)
1165  call check_redundant("Pre-diabatic ", u, v, g)
1166  call mom_forcing_chksum("Pre-diabatic", fluxes, g, haloshift=0)
1167  endif
1168 
1169  if (cs%id_u_predia > 0) call post_data(cs%id_u_predia, u, cs%diag)
1170  if (cs%id_v_predia > 0) call post_data(cs%id_v_predia, v, cs%diag)
1171  if (cs%id_h_predia > 0) call post_data(cs%id_h_predia, h, cs%diag)
1172  if (cs%id_T_predia > 0) call post_data(cs%id_T_predia, cs%tv%T, cs%diag)
1173  if (cs%id_S_predia > 0) call post_data(cs%id_S_predia, cs%tv%S, cs%diag)
1174  if (cs%id_e_predia > 0) then
1175  call find_eta(h, cs%tv, gv%g_Earth, g, gv, eta_predia)
1176  call post_data(cs%id_e_predia, eta_predia, cs%diag)
1177  endif
1178 
1179  call cpu_clock_begin(id_clock_diabatic)
1180  call diabatic(u, v, h, tv, cs%Hml, fluxes, cs%visc, cs%ADp, cs%CDp, &
1181  dtdia, g, gv, cs%diabatic_CSp)
1182  fluxes%fluxes_used = .true.
1183  call cpu_clock_end(id_clock_diabatic)
1184 
1185  if (cs%id_u_preale > 0) call post_data(cs%id_u_preale, u, cs%diag)
1186  if (cs%id_v_preale > 0) call post_data(cs%id_v_preale, v, cs%diag)
1187  if (cs%id_h_preale > 0) call post_data(cs%id_h_preale, h, cs%diag)
1188  if (cs%id_T_preale > 0) call post_data(cs%id_T_preale, tv%T, cs%diag)
1189  if (cs%id_S_preale > 0) call post_data(cs%id_S_preale, tv%S, cs%diag)
1190  if (cs%id_e_preale > 0) then
1191  call find_eta(h, tv, gv%g_Earth, g, gv, eta_preale)
1192  call post_data(cs%id_e_preale, eta_preale, cs%diag)
1193  endif
1194 
1195  if (showcalltree) call calltree_waypoint("finished diabatic (step_MOM_thermo)")
1196 
1197  ! Regridding/remapping is done here, at end of thermodynamics time step
1198  ! (that may comprise several dynamical time steps)
1199  ! The routine 'ALE_main' can be found in 'MOM_ALE.F90'.
1200  if ( cs%use_ALE_algorithm ) then
1201 ! call pass_vector(u, v, G%Domain)
1202  call do_group_pass(cs%pass_T_S_h, g%Domain)
1203 
1204  ! update squared quantities
1205  if (associated(cs%S_squared)) then ; do k=1,nz ; do j=js,je ; do i=is,ie
1206  cs%S_squared(i,j,k) = tv%S(i,j,k)**2
1207  enddo ; enddo ; enddo ; endif
1208  if (associated(cs%T_squared)) then ; do k=1,nz ; do j=js,je ; do i=is,ie
1209  cs%T_squared(i,j,k) = tv%T(i,j,k)**2
1210  enddo ; enddo ; enddo ; endif
1211 
1212  if (cs%debug) then
1213  call mom_state_chksum("Pre-ALE ", u, v, h, cs%uh, cs%vh, g, gv)
1214  call hchksum(tv%T,"Pre-ALE T", g%HI, haloshift=1)
1215  call hchksum(tv%S,"Pre-ALE S", g%HI, haloshift=1)
1216  call check_redundant("Pre-ALE ", u, v, g)
1217  endif
1218  call cpu_clock_begin(id_clock_ale)
1219  if (use_ice_shelf) then
1220  call ale_main(g, gv, h, u, v, tv, cs%tracer_Reg, cs%ALE_CSp, dtdia, &
1221  fluxes%frac_shelf_h)
1222  else
1223  call ale_main(g, gv, h, u, v, tv, cs%tracer_Reg, cs%ALE_CSp, dtdia)
1224  endif
1225 
1226  if (showcalltree) call calltree_waypoint("finished ALE_main (step_MOM_thermo)")
1227  call cpu_clock_end(id_clock_ale)
1228  endif ! endif for the block "if ( CS%use_ALE_algorithm )"
1229 
1230  call cpu_clock_begin(id_clock_pass)
1231  call do_group_pass(cs%pass_uv_T_S_h, g%Domain)
1232  call cpu_clock_end(id_clock_pass)
1233 
1234  if (cs%debug .and. cs%use_ALE_algorithm) then
1235  call mom_state_chksum("Post-ALE ", u, v, h, cs%uh, cs%vh, g, gv)
1236  call hchksum(tv%T, "Post-ALE T", g%HI, haloshift=1)
1237  call hchksum(tv%S, "Post-ALE S", g%HI, haloshift=1)
1238  call check_redundant("Post-ALE ", u, v, g)
1239  endif
1240 
1241  ! Whenever thickness changes let the diag manager know, target grids
1242  ! for vertical remapping may need to be regenerated. This needs to
1243  ! happen after the H update and before the next post_data.
1244  call diag_update_remap_grids(cs%diag)
1245 
1246  call post_diags_ts_vardec(g, cs, dtdia)
1247 
1248  if (cs%debug) then
1249  call uvchksum("Post-diabatic u", u, v, g%HI, haloshift=2)
1250  call hchksum(h, "Post-diabatic h", g%HI, haloshift=1, scale=gv%H_to_m)
1251  call uvchksum("Post-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
1252  haloshift=0, scale=gv%H_to_m)
1253  ! call MOM_state_chksum("Post-diabatic ", u, v, &
1254  ! h, CS%uhtr, CS%vhtr, G, GV, haloshift=1)
1255  if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", g%HI, haloshift=1)
1256  if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", g%HI, haloshift=1)
1257  if (associated(tv%frazil)) call hchksum(tv%frazil, &
1258  "Post-diabatic frazil", g%HI, haloshift=0)
1259  if (associated(tv%salt_deficit)) call hchksum(tv%salt_deficit, &
1260  "Post-diabatic salt deficit", g%HI, haloshift=0)
1261  ! call MOM_thermo_chksum("Post-diabatic ", tv, G)
1262  call check_redundant("Post-diabatic ", u, v, g)
1263  endif
1264 
1265  else ! complement of "if (.not.CS%adiabatic)"
1266 
1267  call cpu_clock_begin(id_clock_diabatic)
1268  call adiabatic(h, tv, fluxes, dtdia, g, gv, cs%diabatic_CSp)
1269  fluxes%fluxes_used = .true.
1270  call cpu_clock_end(id_clock_diabatic)
1271 
1272  if (cs%use_temperature) then
1273  call cpu_clock_begin(id_clock_pass)
1274  call do_group_pass(cs%pass_T_S, g%Domain)
1275  call cpu_clock_end(id_clock_pass)
1276  if (cs%debug) then
1277  if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", g%HI, haloshift=1)
1278  if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", g%HI, haloshift=1)
1279  endif
1280  endif
1281 
1282  endif ! endif for the block "if (.not.CS%adiabatic)"
1283 
1284  if (showcalltree) call calltree_leave("step_MOM_thermo(), MOM.F90")
1285 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ step_offline()

subroutine, public mom::step_offline ( type(forcing), intent(inout)  fluxes,
type(surface), intent(inout)  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,out]fluxespointers to forcing fields
[in,out]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 1294 of file MOM.F90.

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

Referenced by mom_main(), and ocean_model_mod::update_ocean_model().

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

◆ write_static_fields()

subroutine mom::write_static_fields ( type(ocean_grid_type), intent(in)  G,
type(diag_ctrl), intent(inout), target  diag 
)
private

Offers the static fields in the ocean grid type for output via the diag_manager.

Parameters
[in]gocean grid structure
[in,out]diagregulates diagnostic output

Definition at line 3236 of file MOM.F90.

Referenced by initialize_mom().

3236  type(ocean_grid_type), intent(in) :: g !< ocean grid structure
3237  type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output
3238 
3239  ! The out_X arrays are needed because some of the elements of the grid
3240  ! type may be reduced rank macros.
3241  real :: out_h(szi_(g),szj_(g))
3242  integer :: id, i, j, is, ie, js, je
3243  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
3244 
3245  out_h(:,:) = 0.0
3246 
3247  id = register_static_field('ocean_model', 'geolat', diag%axesT1, &
3248  'Latitude of tracer (T) points', 'degrees_N')
3249  if (id > 0) call post_data(id, g%geoLatT, diag, .true.)
3250 
3251  id = register_static_field('ocean_model', 'geolon', diag%axesT1, &
3252  'Longitude of tracer (T) points', 'degrees_E')
3253  if (id > 0) call post_data(id, g%geoLonT, diag, .true.)
3254 
3255  id = register_static_field('ocean_model', 'geolat_c', diag%axesB1, &
3256  'Latitude of corner (Bu) points', 'degrees_N', interp_method='none')
3257  if (id > 0) call post_data(id, g%geoLatBu, diag, .true.)
3258 
3259  id = register_static_field('ocean_model', 'geolon_c', diag%axesB1, &
3260  'Longitude of corner (Bu) points', 'degrees_E', interp_method='none')
3261  if (id > 0) call post_data(id, g%geoLonBu, diag, .true.)
3262 
3263  id = register_static_field('ocean_model', 'geolat_v', diag%axesCv1, &
3264  'Latitude of meridional velocity (Cv) points', 'degrees_N', interp_method='none')
3265  if (id > 0) call post_data(id, g%geoLatCv, diag, .true.)
3266 
3267  id = register_static_field('ocean_model', 'geolon_v', diag%axesCv1, &
3268  'Longitude of meridional velocity (Cv) points', 'degrees_E', interp_method='none')
3269  if (id > 0) call post_data(id, g%geoLonCv, diag, .true.)
3270 
3271  id = register_static_field('ocean_model', 'geolat_u', diag%axesCu1, &
3272  'Latitude of zonal velocity (Cu) points', 'degrees_N', interp_method='none')
3273  if (id > 0) call post_data(id, g%geoLatCu, diag, .true.)
3274 
3275  id = register_static_field('ocean_model', 'geolon_u', diag%axesCu1, &
3276  'Longitude of zonal velocity (Cu) points', 'degrees_E', interp_method='none')
3277  if (id > 0) call post_data(id, g%geoLonCu, diag, .true.)
3278 
3279  id = register_static_field('ocean_model', 'area_t', diag%axesT1, &
3280  'Surface area of tracer (T) cells', 'm2', &
3281  cmor_field_name='areacello', cmor_standard_name='cell_area', &
3282  cmor_units='m2', cmor_long_name='Ocean Grid-Cell Area')
3283  if (id > 0) then
3284  do j=js,je ; do i=is,ie ; out_h(i,j) = g%areaT(i,j) ; enddo ; enddo
3285  call post_data(id, out_h, diag, .true.)
3286  call diag_register_area_ids(diag, id_area_t=id)
3287  endif
3288 
3289  id = register_static_field('ocean_model', 'depth_ocean', diag%axesT1, &
3290  'Depth of the ocean at tracer points', 'm', &
3291  standard_name='sea_floor_depth_below_geoid', &
3292  cmor_field_name='deptho', cmor_long_name='Sea Floor Depth', &
3293  cmor_units='m', cmor_standard_name='sea_floor_depth_below_geoid',&
3294  area=diag%axesT1%id_area)
3295  if (id > 0) call post_data(id, g%bathyT, diag, .true., mask=g%mask2dT)
3296 
3297  id = register_static_field('ocean_model', 'wet', diag%axesT1, &
3298  '0 if land, 1 if ocean at tracer points', 'none', area=diag%axesT1%id_area)
3299  if (id > 0) call post_data(id, g%mask2dT, diag, .true.)
3300 
3301  id = register_static_field('ocean_model', 'wet_c', diag%axesB1, &
3302  '0 if land, 1 if ocean at corner (Bu) points', 'none', interp_method='none')
3303  if (id > 0) call post_data(id, g%mask2dBu, diag, .true.)
3304 
3305  id = register_static_field('ocean_model', 'wet_u', diag%axesCu1, &
3306  '0 if land, 1 if ocean at zonal velocity (Cu) points', 'none', interp_method='none')
3307  if (id > 0) call post_data(id, g%mask2dCu, diag, .true.)
3308 
3309  id = register_static_field('ocean_model', 'wet_v', diag%axesCv1, &
3310  '0 if land, 1 if ocean at meridional velocity (Cv) points', 'none', interp_method='none')
3311  if (id > 0) call post_data(id, g%mask2dCv, diag, .true.)
3312 
3313  id = register_static_field('ocean_model', 'Coriolis', diag%axesB1, &
3314  'Coriolis parameter at corner (Bu) points', 's-1', interp_method='none')
3315  if (id > 0) call post_data(id, g%CoriolisBu, diag, .true.)
3316 
3317  id = register_static_field('ocean_model', 'dxt', diag%axesT1, &
3318  'Delta(x) at thickness/tracer points (meter)', 'm', interp_method='none')
3319  if (id > 0) call post_data(id, g%dxt, diag, .true.)
3320 
3321  id = register_static_field('ocean_model', 'dyt', diag%axesT1, &
3322  'Delta(y) at thickness/tracer points (meter)', 'm', interp_method='none')
3323  if (id > 0) call post_data(id, g%dyt, diag, .true.)
3324 
3325  id = register_static_field('ocean_model', 'dxCu', diag%axesCu1, &
3326  'Delta(x) at u points (meter)', 'm', interp_method='none')
3327  if (id > 0) call post_data(id, g%dxCu, diag, .true.)
3328 
3329  id = register_static_field('ocean_model', 'dyCu', diag%axesCu1, &
3330  'Delta(y) at u points (meter)', 'm', interp_method='none')
3331  if (id > 0) call post_data(id, g%dyCu, diag, .true.)
3332 
3333  id = register_static_field('ocean_model', 'dxCv', diag%axesCv1, &
3334  'Delta(x) at v points (meter)', 'm', interp_method='none')
3335  if (id > 0) call post_data(id, g%dxCv, diag, .true.)
3336 
3337  id = register_static_field('ocean_model', 'dyCv', diag%axesCv1, &
3338  'Delta(y) at v points (meter)', 'm', interp_method='none')
3339  if (id > 0) call post_data(id, g%dyCv, diag, .true.)
3340 
Here is the caller graph for this function:

Variable Documentation

◆ id_clock_ale

integer mom::id_clock_ale
private

Definition at line 453 of file MOM.F90.

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

453 integer :: id_clock_ale

◆ id_clock_bbl_visc

integer mom::id_clock_bbl_visc
private

Definition at line 445 of file MOM.F90.

Referenced by mom_timing_init(), and step_mom().

445 integer :: id_clock_bbl_visc

◆ id_clock_continuity

integer mom::id_clock_continuity
private

Definition at line 443 of file MOM.F90.

Referenced by mom_timing_init().

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

◆ id_clock_diabatic

integer mom::id_clock_diabatic
private

Definition at line 442 of file MOM.F90.

Referenced by mom_timing_init(), and step_mom_thermo().

442 integer :: id_clock_diabatic

◆ id_clock_diagnostics

integer mom::id_clock_diagnostics
private

Definition at line 447 of file MOM.F90.

Referenced by mom_timing_init(), and step_mom().

447 integer :: id_clock_diagnostics

◆ id_clock_dynamics

integer mom::id_clock_dynamics
private

Definition at line 439 of file MOM.F90.

Referenced by mom_timing_init(), and step_mom().

439 integer :: id_clock_dynamics

◆ id_clock_init

integer mom::id_clock_init
private

Definition at line 449 of file MOM.F90.

Referenced by finish_mom_initialization(), and initialize_mom().

449 integer :: id_clock_init

◆ id_clock_ml_restrat

integer mom::id_clock_ml_restrat
private

Definition at line 446 of file MOM.F90.

Referenced by mom_timing_init(), and step_mom().

446 integer :: id_clock_ml_restrat

◆ id_clock_mom_init

integer mom::id_clock_mom_init
private

Definition at line 450 of file MOM.F90.

Referenced by initialize_mom(), and mom_timing_init().

450 integer :: id_clock_mom_init

◆ id_clock_ocean

integer mom::id_clock_ocean

Definition at line 438 of file MOM.F90.

Referenced by mom_timing_init(), and step_mom().

438 integer :: id_clock_ocean

◆ id_clock_offline_tracer

integer mom::id_clock_offline_tracer
private

Definition at line 455 of file MOM.F90.

Referenced by mom_timing_init(), and step_offline().

455 integer :: id_clock_offline_tracer

◆ id_clock_other

integer mom::id_clock_other
private

Definition at line 454 of file MOM.F90.

Referenced by mom_timing_init(), and step_mom().

454 integer :: id_clock_other

◆ id_clock_pass

integer mom::id_clock_pass
private

Definition at line 451 of file MOM.F90.

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

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

◆ id_clock_pass_init

integer mom::id_clock_pass_init
private

Definition at line 452 of file MOM.F90.

Referenced by initialize_mom(), and mom_timing_init().

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

◆ id_clock_thermo

integer mom::id_clock_thermo
private

Definition at line 440 of file MOM.F90.

Referenced by mom_timing_init(), and step_mom().

440 integer :: id_clock_thermo

◆ id_clock_thick_diff

integer mom::id_clock_thick_diff
private

Definition at line 444 of file MOM.F90.

Referenced by mom_timing_init(), and step_mom().

444 integer :: id_clock_thick_diff

◆ id_clock_tracer

integer mom::id_clock_tracer
private

Definition at line 441 of file MOM.F90.

Referenced by mom_timing_init(), and step_mom().

441 integer :: id_clock_tracer

◆ id_clock_z_diag

integer mom::id_clock_z_diag
private

Definition at line 448 of file MOM.F90.

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

448 integer :: id_clock_z_diag