MOM6
Data Types
mom Module Reference

Detailed Description

The central module of the MOM6 ocean model.

Modular Ocean Model (MOM) Version 6.0 (MOM6)

Authors
Alistair Adcroft, Robert Hallberg, and Stephen Griffies

Additional contributions from:

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

MOM ice-shelf code was developed by

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

Overview of MOM

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

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

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

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

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

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

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

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

Structure of MOM

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

The typical MOM directory tree is as follows:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Diagnosing MOM heat budget

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

Depth integrated heat budget

Depth integrated heat budget diagnostic for MOM.

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

Depth integrated heat budget

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

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

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

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

  • OPOTTEMPTEND_2d = T_ADVECTION_XY_2d + OPOTTEMPPMDIFF_2d + HFDS

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

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

  • SFDSI = BOUNDARY_FORCING_SALT_TENDENCY[k=@sum]

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

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

Data Types

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

Function/Subroutine Documentation

◆ adjust_ssh_for_p_atm()

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

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

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

Definition at line 2642 of file MOM.F90.

2642  type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
2643  type(ocean_grid_type), intent(in) :: G !< ocean grid structure
2644  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2645  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2646  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: ssh !< time mean surface height [m]
2647  real, dimension(:,:), optional, pointer :: p_atm !< atmospheric pressure [Pa]
2648  logical, optional, intent(in) :: use_EOS !< If true, calculate the density for
2649  !! the SSH correction using the equation of state.
2650 
2651  real :: Rho_conv ! The density used to convert surface pressure to
2652  ! a corrected effective SSH [kg m-3].
2653  real :: IgR0 ! The SSH conversion factor from Pa to m [m Pa-1].
2654  logical :: calc_rho
2655  integer :: i, j, is, ie, js, je
2656 
2657  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2658  if (present(p_atm)) then ; if (associated(p_atm)) then
2659  calc_rho = associated(tv%eqn_of_state)
2660  if (present(use_eos) .and. calc_rho) calc_rho = use_eos
2661  ! Correct the output sea surface height for the contribution from the
2662  ! atmospheric pressure
2663  do j=js,je ; do i=is,ie
2664  if (calc_rho) then
2665  call calculate_density(tv%T(i,j,1), tv%S(i,j,1), p_atm(i,j)/2.0, &
2666  rho_conv, tv%eqn_of_state)
2667  else
2668  rho_conv=gv%Rho0
2669  endif
2670  igr0 = 1.0 / (rho_conv * gv%mks_g_Earth)
2671  ssh(i,j) = ssh(i,j) + p_atm(i,j) * igr0
2672  enddo ; enddo
2673  endif ; endif
2674 

◆ extract_surface_state()

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

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

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

Definition at line 2681 of file MOM.F90.

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

◆ finish_mom_initialization()

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

Finishes initializing MOM and writes out the initial conditions.

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

Definition at line 2459 of file MOM.F90.

2459  type(time_type), intent(in) :: Time !< model time, used in this routine
2460  type(directories), intent(in) :: dirs !< structure with directory paths
2461  type(MOM_control_struct), pointer :: CS !< pointer to MOM control structure
2462  type(MOM_restart_CS), pointer :: restart_CSp !< pointer to the restart control
2463  !! structure that will be used for MOM.
2464  ! Local variables
2465  type(ocean_grid_type), pointer :: G => null() ! pointer to a structure containing
2466  ! metrics and related information
2467  type(verticalGrid_type), pointer :: GV => null() ! Pointer to the vertical grid structure
2468  type(unit_scale_type), pointer :: US => null() ! Pointer to a structure containing
2469  ! various unit conversion factors
2470  type(MOM_restart_CS), pointer :: restart_CSp_tmp => null()
2471  real, allocatable :: z_interface(:,:,:) ! Interface heights [m]
2472  type(vardesc) :: vd
2473 
2474  call cpu_clock_begin(id_clock_init)
2475  call calltree_enter("finish_MOM_initialization()")
2476 
2477  ! Pointers for convenience
2478  g => cs%G ; gv => cs%GV ; us => cs%US
2479 
2480  !### Move to initialize_MOM?
2481  call fix_restart_scaling(gv)
2482  call fix_restart_unit_scaling(us)
2483 
2484  ! Write initial conditions
2485  if (cs%write_IC) then
2486  allocate(restart_csp_tmp)
2487  restart_csp_tmp = restart_csp
2488  allocate(z_interface(szi_(g),szj_(g),szk_(g)+1))
2489  call find_eta(cs%h, cs%tv, g, gv, us, z_interface, eta_to_m=1.0)
2490  call register_restart_field(z_interface, "eta", .true., restart_csp_tmp, &
2491  "Interface heights", "meter", z_grid='i')
2492 
2493  call save_restart(dirs%output_directory, time, g, &
2494  restart_csp_tmp, filename=cs%IC_file, gv=gv)
2495  deallocate(z_interface)
2496  deallocate(restart_csp_tmp)
2497  endif
2498 
2499  call write_energy(cs%u, cs%v, cs%h, cs%tv, time, 0, g, gv, us, &
2500  cs%sum_output_CSp, cs%tracer_flow_CSp)
2501 
2502  call calltree_leave("finish_MOM_initialization()")
2503  call cpu_clock_end(id_clock_init)
2504 

◆ get_mom_state_elements()

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

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

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

Definition at line 3046 of file MOM.F90.

3046  type(MOM_control_struct), pointer :: CS !< MOM control structure
3047  type(ocean_grid_type), &
3048  optional, pointer :: G !< structure containing metrics and grid info
3049  type(verticalGrid_type), &
3050  optional, pointer :: GV !< structure containing vertical grid info
3051  type(unit_scale_type), &
3052  optional, pointer :: US !< A dimensional unit scaling type
3053  real, optional, intent(out) :: C_p !< The heat capacity
3054  logical, optional, intent(out) :: use_temp !< Indicates whether temperature is a state variable
3055 
3056  if (present(g)) g => cs%G
3057  if (present(gv)) gv => cs%GV
3058  if (present(us)) us => cs%US
3059  if (present(c_p)) c_p = cs%tv%C_p
3060  if (present(use_temp)) use_temp = associated(cs%tv%T)

◆ get_ocean_stocks()

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

Find the global integrals of various quantities.

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

Definition at line 3065 of file MOM.F90.

3065  type(MOM_control_struct), pointer :: CS !< MOM control structure
3066  real, optional, intent(out) :: heat !< The globally integrated integrated ocean heat [J].
3067  real, optional, intent(out) :: salt !< The globally integrated integrated ocean salt [kg].
3068  real, optional, intent(out) :: mass !< The globally integrated integrated ocean mass [kg].
3069  logical, optional, intent(in) :: on_PE_only !< If present and true, only sum on the local PE.
3070 
3071  if (present(mass)) &
3072  mass = global_mass_integral(cs%h, cs%G, cs%GV, on_pe_only=on_pe_only)
3073  if (present(heat)) &
3074  heat = cs%tv%C_p * global_mass_integral(cs%h, cs%G, cs%GV, cs%tv%T, on_pe_only=on_pe_only)
3075  if (present(salt)) &
3076  salt = 1.0e-3 * global_mass_integral(cs%h, cs%G, cs%GV, cs%tv%S, on_pe_only=on_pe_only)
3077 

◆ initialize_mom()

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

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

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

Definition at line 1493 of file MOM.F90.

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

◆ mom_end()

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

End of ocean model, including memory deallocation.

Parameters
csMOM control structure

Definition at line 3082 of file MOM.F90.

3082  type(MOM_control_struct), pointer :: CS !< MOM control structure
3083 
3084  if (cs%use_ALE_algorithm) call ale_end(cs%ALE_CSp)
3085 
3086  dealloc_(cs%u) ; dealloc_(cs%v) ; dealloc_(cs%h)
3087  dealloc_(cs%uh) ; dealloc_(cs%vh)
3088 
3089  if (associated(cs%tv%T)) then
3090  dealloc_(cs%T) ; cs%tv%T => null() ; dealloc_(cs%S) ; cs%tv%S => null()
3091  endif
3092  if (associated(cs%tv%frazil)) deallocate(cs%tv%frazil)
3093  if (associated(cs%tv%salt_deficit)) deallocate(cs%tv%salt_deficit)
3094  if (associated(cs%Hml)) deallocate(cs%Hml)
3095 
3096  call tracer_advect_end(cs%tracer_adv_CSp)
3097  call tracer_hor_diff_end(cs%tracer_diff_CSp)
3098  call tracer_registry_end(cs%tracer_Reg)
3099  call tracer_flow_control_end(cs%tracer_flow_CSp)
3100 
3101  call diabatic_driver_end(cs%diabatic_CSp)
3102 
3103  if (cs%offline_tracer_mode) call offline_transport_end(cs%offline_CSp)
3104 
3105  dealloc_(cs%uhtr) ; dealloc_(cs%vhtr)
3106  if (cs%split) then
3107  call end_dyn_split_rk2(cs%dyn_split_RK2_CSp)
3108  elseif (cs%use_RK2) then
3109  call end_dyn_unsplit_rk2(cs%dyn_unsplit_RK2_CSp)
3110  else
3111  call end_dyn_unsplit(cs%dyn_unsplit_CSp)
3112  endif
3113  dealloc_(cs%ave_ssh_ibc) ; dealloc_(cs%ssh_rint)
3114  if (associated(cs%update_OBC_CSp)) call obc_register_end(cs%update_OBC_CSp)
3115 
3116  call verticalgridend(cs%GV)
3117  call unit_scaling_end(cs%US)
3118  call mom_grid_end(cs%G)
3119 
3120  deallocate(cs)
3121 

◆ mom_state_is_synchronized()

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

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

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

Definition at line 3025 of file MOM.F90.

3025  type(MOM_control_struct), pointer :: CS !< MOM control structure
3026  logical, optional, intent(in) :: adv_dyn !< If present and true, only check
3027  !! whether the advection is up-to-date with
3028  !! the dynamics.
3029  logical :: in_synch !< True if all phases of the update are synchronized.
3030 
3031  logical :: adv_only
3032 
3033  adv_only = .false. ; if (present(adv_dyn)) adv_only = adv_dyn
3034 
3035  if (adv_only) then
3036  in_synch = (cs%t_dyn_rel_adv == 0.0)
3037  else
3038  in_synch = ((cs%t_dyn_rel_adv == 0.0) .and. (cs%t_dyn_rel_thermo == 0.0))
3039  endif
3040 

◆ mom_timing_init()

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

Set up CPU clock IDs for timing various subroutines.

Parameters
[in]cscontrol structure set up by initialize_MOM.

Definition at line 2539 of file MOM.F90.

2539  type(MOM_control_struct), intent(in) :: CS !< control structure set up by initialize_MOM.
2540 
2541  id_clock_ocean = cpu_clock_id('Ocean', grain=clock_component)
2542  id_clock_dynamics = cpu_clock_id('Ocean dynamics', grain=clock_subcomponent)
2543  id_clock_thermo = cpu_clock_id('Ocean thermodynamics and tracers', grain=clock_subcomponent)
2544  id_clock_other = cpu_clock_id('Ocean Other', grain=clock_subcomponent)
2545  id_clock_tracer = cpu_clock_id('(Ocean tracer advection)', grain=clock_module_driver)
2546  if (.not.cs%adiabatic) &
2547  id_clock_diabatic = cpu_clock_id('(Ocean diabatic driver)', grain=clock_module_driver)
2548 
2549  id_clock_continuity = cpu_clock_id('(Ocean continuity equation *)', grain=clock_module)
2550  id_clock_bbl_visc = cpu_clock_id('(Ocean set BBL viscosity)', grain=clock_module)
2551  id_clock_pass = cpu_clock_id('(Ocean message passing *)', grain=clock_module)
2552  id_clock_mom_init = cpu_clock_id('(Ocean MOM_initialize_state)', grain=clock_module)
2553  id_clock_pass_init = cpu_clock_id('(Ocean init message passing *)', grain=clock_routine)
2554  if (cs%thickness_diffuse) &
2555  id_clock_thick_diff = cpu_clock_id('(Ocean thickness diffusion *)', grain=clock_module)
2556 !if (CS%mixedlayer_restrat) &
2557  id_clock_ml_restrat = cpu_clock_id('(Ocean mixed layer restrat)', grain=clock_module)
2558  id_clock_diagnostics = cpu_clock_id('(Ocean collective diagnostics)', grain=clock_module)
2559  id_clock_z_diag = cpu_clock_id('(Ocean Z-space diagnostics)', grain=clock_module)
2560  id_clock_ale = cpu_clock_id('(Ocean ALE)', grain=clock_module)
2561  if (cs%offline_tracer_mode) then
2562  id_clock_offline_tracer = cpu_clock_id('Ocean offline tracers', grain=clock_subcomponent)
2563  endif
2564 

◆ register_diags()

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

Register certain diagnostics.

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

Definition at line 2509 of file MOM.F90.

2509  type(time_type), intent(in) :: Time !< current model time
2510  type(ocean_grid_type), intent(in) :: G !< ocean grid structure
2511  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2512  type(MOM_diag_IDs), intent(inout) :: IDs !< A structure with the diagnostic IDs.
2513  type(diag_ctrl), intent(inout) :: diag !< regulates diagnostic output
2514 
2515  real :: H_convert
2516  character(len=48) :: thickness_units
2517 
2518  thickness_units = get_thickness_units(gv)
2519  if (gv%Boussinesq) then
2520  h_convert = gv%H_to_m
2521  else
2522  h_convert = gv%H_to_kg_m2
2523  endif
2524 
2525  ! Diagnostics of the rapidly varying dynamic state
2526  ids%id_u = register_diag_field('ocean_model', 'u_dyn', diag%axesCuL, time, &
2527  'Zonal velocity after the dynamics update', 'm s-1')
2528  ids%id_v = register_diag_field('ocean_model', 'v_dyn', diag%axesCvL, time, &
2529  'Meridional velocity after the dynamics update', 'm s-1')
2530  ids%id_h = register_diag_field('ocean_model', 'h_dyn', diag%axesTL, time, &
2531  'Layer Thickness after the dynamics update', thickness_units, &
2532  v_extensive=.true., conversion=h_convert)
2533  ids%id_ssh_inst = register_diag_field('ocean_model', 'SSH_inst', diag%axesT1, &
2534  time, 'Instantaneous Sea Surface Height', 'm')

◆ set_restart_fields()

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

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

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

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

Definition at line 2577 of file MOM.F90.

2577  type(verticalGrid_type), intent(inout) :: GV !< ocean vertical grid structure
2578  type(unit_scale_type), intent(inout) :: US !< A dimensional unit scaling type
2579  type(param_file_type), intent(in) :: param_file !< opened file for parsing to get parameters
2580  type(MOM_control_struct), intent(in) :: CS !< control structure set up by inialize_MOM
2581  type(MOM_restart_CS), pointer :: restart_CSp !< pointer to the restart control
2582  !! structure that will be used for MOM.
2583  ! Local variables
2584  logical :: use_ice_shelf ! Needed to determine whether to add CS%Hml to restarts
2585  character(len=48) :: thickness_units, flux_units
2586 
2587 
2588  thickness_units = get_thickness_units(gv)
2589  flux_units = get_flux_units(gv)
2590 
2591  if (associated(cs%tv%T)) &
2592  call register_restart_field(cs%tv%T, "Temp", .true., restart_csp, &
2593  "Potential Temperature", "degC")
2594  if (associated(cs%tv%S)) &
2595  call register_restart_field(cs%tv%S, "Salt", .true., restart_csp, &
2596  "Salinity", "PPT")
2597 
2598  call register_restart_field(cs%h, "h", .true., restart_csp, &
2599  "Layer Thickness", thickness_units)
2600 
2601  call register_restart_field(cs%u, "u", .true., restart_csp, &
2602  "Zonal velocity", "m s-1", hor_grid='Cu')
2603 
2604  call register_restart_field(cs%v, "v", .true., restart_csp, &
2605  "Meridional velocity", "m s-1", hor_grid='Cv')
2606 
2607  if (associated(cs%tv%frazil)) &
2608  call register_restart_field(cs%tv%frazil, "frazil", .false., restart_csp, &
2609  "Frazil heat flux into ocean", "J m-2")
2610 
2611  if (cs%interp_p_surf) then
2612  call register_restart_field(cs%p_surf_prev, "p_surf_prev", .false., restart_csp, &
2613  "Previous ocean surface pressure", "Pa")
2614  endif
2615 
2616  call register_restart_field(cs%ave_ssh_ibc, "ave_ssh", .false., restart_csp, &
2617  "Time average sea surface height", "meter")
2618 
2619  ! hML is needed when using the ice shelf module
2620  call get_param(param_file, '', "ICE_SHELF", use_ice_shelf, default=.false., &
2621  do_not_log=.true.)
2622  if (use_ice_shelf .and. associated(cs%Hml)) then
2623  call register_restart_field(cs%Hml, "hML", .false., restart_csp, &
2624  "Mixed layer thickness", "meter")
2625  endif
2626 
2627  ! Register scalar unit conversion factors.
2628  call register_restart_field(us%m_to_Z_restart, "m_to_Z", .false., restart_csp, &
2629  "Height unit conversion factor", "Z meter-1")
2630  call register_restart_field(gv%m_to_H_restart, "m_to_H", .false., restart_csp, &
2631  "Thickness unit conversion factor", "H meter-1")
2632  call register_restart_field(us%m_to_Z_restart, "m_to_L", .false., restart_csp, &
2633  "Length unit conversion factor", "L meter-1")
2634  call register_restart_field(us%s_to_T_restart, "s_to_T", .false., restart_csp, &
2635  "Time unit conversion factor", "T second-1")
2636 

◆ step_mom()

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

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

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

Definition at line 397 of file MOM.F90.

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

◆ step_mom_dynamics()

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

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

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

Definition at line 871 of file MOM.F90.

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

◆ step_mom_thermo()

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

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

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

Definition at line 1124 of file MOM.F90.

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

◆ step_mom_tracer_dyn()

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

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

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

Definition at line 1056 of file MOM.F90.

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

◆ step_offline()

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

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

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

Definition at line 1301 of file MOM.F90.

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