mod_scalars.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE mod_scalars
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!
        USE mod_kinds
#ifdef NO_4BYTE_REALS
        USE netcdf
#endif
!
        PUBLIC :: allocate_scalars
        PUBLIC :: deallocate_scalars
        PUBLIC :: initialize_scalars
!
!-----------------------------------------------------------------------
!  Multiple grid structure.
!-----------------------------------------------------------------------
!
#if defined READ_WATER && defined MASKING && defined DISTRIBUTE
!    IJwater       IJ-indices of water points.
#endif
!    Fstate        Logical switches to control computations of the
!                    Forcing Singular Vectors or Stochastic Optimals.
!    Lstate        Logical switches to control computations of the
!                    model state.
#ifdef STATIONS
!    Sflag         Station extraction special flag:
!                    Sflag = 0  => locations in terms of (I,J) pairs.
!                    Sflag = 1  => locations in terms of (lon,lat) pairs.
!    SposX         Longitude or frational I-coordinate station location.
!    SposY         Latitude  or frational J-coordinate station location.
#endif
!    Cs_r          Set of S-curves used to stretch the vertical grid
!                    that follows the bathymetry at vertical RHO-points.
!    Cs_w          Set of S-curves used to stretch the vertical grid
!                    that follows the bathymetry at vertical W-points.
!    sc_r          S-coordinate independent variable, [-1 < sc < 0] at
!                    vertical RHO-points.
!    sc_w          S-coordinate independent variable, [-1 < sc < 0] at
!                    vertical W-points.
!
        TYPE t_scalars
          logical, pointer :: fstate(:)
          logical, pointer :: lstate(:)
#if defined READ_WATER && defined MASKING && defined DISTRIBUTE
          integer , pointer :: ijwater(:,:)
#endif
#ifdef STATIONS
          integer,  pointer :: sflag(:)
 
          real(r8), pointer :: sposx(:)
          real(r8), pointer :: sposy(:)
#endif
          real(dp), pointer :: cs_r(:)
          real(dp), pointer :: cs_w(:)
          real(dp), pointer :: sc_r(:)
          real(dp), pointer :: sc_w(:)
        END TYPE t_scalars
!
        TYPE (t_scalars), allocatable :: scalars(:)
!
!-----------------------------------------------------------------------
!  Time clock structure.
!-----------------------------------------------------------------------
!
!  Reference time (yyyymmdd.f) used to compute relative time. The
!  application date clock is measured ad elapsed time interval since
!  reference-time. This parameter also provides information about the
!  calendar used:
!
!    If TIME_REF = -2, the model time and DSTART are in modified Julian
!                  days units.  The time "units" attribute is:
!
!                  'time-units since 1968-05-23 00:00:00 GMT'
!
!    If TIME_REF = -1, the model time and DSTART are in a calendar
!                  with 360 days in every year (30 days each month).
!                  The time "units" attribute is:
!
!                  'time-units since 0001-01-01 00:00:00'
!
!    If TIME_REF = 0, the model time and DSTART are in a common year
!                  calendar with 365.2524 days.  The "units" attribute
!                  is:
!
!                  'time-units since 0001-01-01 00:00:00'
!
!    If TIME_REF > 0, the model time and DSTART are the elapsed time
!                  units since specified reference time.  For example,
!                  TIME_REF=20020115.5 will yield the following
!                  time "units" attribute:
!
!                  'time-units since 2002-01-15 12:00:00'
!
        real(dp) :: time_ref = 0.0_dp                    ! YYYYMMDD.dd
!
      TYPE t_clock
        integer :: yday                ! day of the year
        integer :: year                ! year including century (YYYY)
        integer :: month               ! month of the year (1,...,12)
        integer :: day                 ! day of the month
        integer :: hour                ! hour of the day (1,...,23)
        integer :: minutes             ! minutes of the hour
 
        real(dp) :: seconds            ! frational seconds of the minute
        real(dp) :: base               ! reference date (YYYYMMDD.dd)
        real(dp) :: datenumber(2)      ! date number, [1]: days
                                       !              [2]: seconds
        real(dp) :: tide_datenumber(2) ! tide reference date number,
                                       !   [1]: days  [2]: seconds
        character (len=22) :: string   ! YYYY-MM-DD hh:mm:ss.ss
        character (len=25) :: calendar ! date calendar
      END TYPE t_clock
!
        TYPE (t_clock) :: rclock       ! reference/base date
!
!-----------------------------------------------------------------------
!  Tracer identification indices.
!-----------------------------------------------------------------------
!
        integer :: itemp              ! Potential temperature
        integer :: isalt              ! Salinity
#ifdef T_PASSIVE
        integer, pointer :: inert(:)  ! inert tracers
#endif
#ifdef DIAGNOSTICS
!
!-----------------------------------------------------------------------
!  Diagnostic fields identification indices.
!-----------------------------------------------------------------------
!
# ifdef DIAGNOSTICS_TS
        integer :: itrate          ! Tracer, time rate of change
        integer :: itvadv          ! Tracer, vertical advection
        integer :: ithadv          ! Tracer, horizontal advection
        integer :: itxadv          ! Tracer, horizontal X-advection
        integer :: ityadv          ! Tracer, horizontal Y-advection
        integer :: itvdif          ! Tracer, vertical diffusion
        integer :: ithdif          ! Tracer, horizontal diffusion
        integer :: itxdif          ! Tracer, horizontal X-diffusion
        integer :: itydif          ! Tracer, horizontal Y-diffusion
        integer :: itsdif          ! Tracer, horizontal S-diffusion
# endif
# ifdef DIAGNOSTICS_UV
        integer :: m2fcor          ! 2D momentum, Coriolis
        integer :: m2hadv          ! 2D momentum, horizontal advection
        integer :: m2xadv          ! 2D momentum, horizontal X-advection
        integer :: m2yadv          ! 2D momentum, horizontal Y-advection
#  ifdef WEC_VF
        integer :: m2hjvf          ! 2D momentum, horizontal J vortex force
        integer :: m2kvrf          ! 2D momentum, K vortex force
        integer :: m2fsco          ! 2D momentum, coriolis-stokes
        integer :: m2bstm          ! 2D momentum, bottom streaming
        integer :: m2sstm          ! 2D momentum, surface streaming
        integer :: m2wrol          ! 2D momentum, wave roller accel
        integer :: m2wbrk          ! 2D momentum, wave breaking
        integer :: m2zeta          ! 2D momentum, Eulerian sea level adjustment
        integer :: m2zetw          ! 2D momentum, quasi-static sea level adjustment
        integer :: m2zqsp          ! 2D momentum, quasi-static pressure
        integer :: m2zbeh          ! 2D momentum, Bernoulli head
#  endif
        integer :: m2pgrd          ! 2D momentum, pressure gradient
        integer :: m2hvis          ! 2D momentum, horizontal viscosity
        integer :: m2xvis          ! 2D momentum, horizontal X-viscosity
        integer :: m2yvis          ! 2D momentum, horizontal Y-viscosity
        integer :: m2sstr          ! 2D momentum, surface stress
        integer :: m2bstr          ! 2D momentum, bottom stress
        integer :: m2rate          ! 2D momentum, time rate of change
#  ifdef SOLVE3D
        integer :: m3fcor          ! 3D momentum, Coriolis
        integer :: m3vadv          ! 3D momentum, vertical advection
        integer :: m3hadv          ! 3D momentum, horizontal advection
        integer :: m3xadv          ! 3D momentum, horizontal X-advection
        integer :: m3yadv          ! 3D momentum, horizontal Y-advection
#   ifdef WEC_VF
        integer :: m3vjvf          ! 3D momentum, vertical J vortex force
        integer :: m3hjvf          ! 3D momentum, horizontal J vortex force
        integer :: m3kvrf          ! 3D momentum, K vortex force
        integer :: m3fsco          ! 3D momentum, coriolis-stokes
        integer :: m3bstm          ! 3D momentum, bottom streaming
        integer :: m3sstm          ! 3D momentum, surface streaming
        integer :: m3wrol          ! 3D momentum, wave roller accel
        integer :: m3wbrk          ! 3D momentum, wave breaking
#   endif
        integer :: m3pgrd          ! 3D momentum, pressure gradient
        integer :: m3vvis          ! 3D momentum, vertical viscosity
        integer :: m3hvis          ! 3D momentum, horizontal viscosity
        integer :: m3xvis          ! 3D momentum, horizontal X-viscosity
        integer :: m3yvis          ! 3D momentum, horizontal Y-viscosity
        integer :: m3rate          ! 3D momentum, time rate of change
#  endif
# endif
#endif
!
!-----------------------------------------------------------------------
!  Time stepping indices, variables, and clocks.
!-----------------------------------------------------------------------
!
!    indx1         2D timestep rolling counter.
!    iic           Timestep counter for 3D primitive equations.
!    iif           Timestep counter for 2D primitive equations.
#ifdef JEDI
!    jic           Timestep counter for roms-jedi interface.
#endif
!    ndtfast       Number of barotropic timesteps between each
!                    baroclinic timestep.
!    nfast         Number of barotropic timesteps needed to compute
!                    time-averaged barotropic variables centered at
!                    time level n+1.
!    dt            Size baroclinic timestep (s).
!    dtfast        Size barotropic timestep (s).
!    run_time      Total run time for all nested grids (s), it is
!                    set in Masters/ocean.h
!    MyRunInterval Total run time for all nested grids (s), it is
!                    set in Drivers/nl_ocean.h (coupling window)
!    io_time       Current I/O time (s) processed in "get_state".
!    tdays         Model time clock (days).
!    time          Model time clock (s).
#ifdef JEDI
!    time4jedi     Model time clock (s) to pass to roms-jedi interface.
#endif
!    time_code     Model time clock (string, YYYY-MM-DD hh:mm:ss.ss)
!    AVGtime       Model time clock for averages output (s).
!    DIAtime       Model time clock for diagnostics output (s).
!    F_code        Final time string for simulation
!    I_code        Initial time string for simulation
!    INItime       Nonlinear model initial conditions time (s).
!    INItimeS      Saved nonlinear model initial conditions time (s).
!    IMPtime       Impulse forcing time (s) to process.
!    ObsTime       Observation time (s) to process.
!    FrcTime       Adjoint or tangent linear Impulse forcing time (s).
!    dstart        Time stamp assigned to model initialization (usually
!                    a Calendar day, like modified Julian Day).
#ifdef SP4DVAR
!    dstartS       Saved original dstart for saddle-point 4D-Var.
#endif
!    tide_start    Reference time for tidal forcing (days).
!
        logical, allocatable :: perfectrst(:)
        logical, allocatable :: predictor_2d_step(:)
!$OMP THREADPRIVATE (PREDICTOR_2D_STEP)
 
        integer, allocatable :: indx1(:)
        integer, allocatable :: iic(:)
        integer, allocatable :: iif(:)
        integer, allocatable :: next_kstp(:)
!$OMP THREADPRIVATE (indx1, iic, iif, next_kstp)
 
#ifdef JEDI
        integer, allocatable :: jic(:)
#endif
 
        integer, allocatable :: ndtfast(:)
        integer, allocatable :: nfast(:)
 
        real(dp), allocatable :: tdays(:)                ! days
        real(dp), allocatable :: time(:)                 ! seconds
!$OMP THREADPRIVATE (tdays, time)
 
#ifdef JEDI
        real(dp), allocatable :: time4jedi(:)            ! seconds
#endif
 
        real(dp), allocatable :: dt(:)                   ! seconds
        real(dp), allocatable :: dtfast(:)               ! seconds
 
        real(dp), allocatable :: timeend(:)              ! seconds
        real(dp), allocatable :: avgtime(:)              ! seconds
        real(dp), allocatable :: diatime(:)              ! seconds
        real(dp), allocatable :: imptime(:)              ! seconds
        real(dp), allocatable :: initime(:)              ! seconds
        real(dp), allocatable :: initimes(:)             ! seconds
        real(dp), allocatable :: obstime(:)              ! seconds
        real(dp), allocatable :: frctime(:)              ! seconds
 
        real(dp) :: dstart = 0.0_dp                      ! days
#ifdef SP4DVAR
        real(dp) :: dstarts = 0.0_dp                     ! days
#endif
        real(dp) :: io_time = 0.0_dp                     ! seconds
        real(dp) :: run_time = 0.0_dp                    ! seconds
        real(dp) :: myruninterval = 0.0_dp               ! seconds
        real(dp) :: tide_start = 0.0_dp                  ! days
 
        character (len=22) :: f_code, i_code
 
        character (len=22), allocatable :: time_code(:)  ! date string
!$OMP THREADPRIVATE (time_code)
 
#if defined POWER_LAW && defined SOLVE3D
!
!  Power-law shape filter parameters for time-averaging of barotropic
!  Fields.  The power-law shape filters are given by:
!
!     F(xi)=xi^Falpha*(1-xi^Fbeta)-Fgamma*xi
!
!  Possible settings of parameters to yield the second-order accuracy:
!
!     Falpha  Fbeta      Fgamma
!     ------------------------------
!      2.0     1.0    0.1181  0.169     The problem here is setting
!      2.0     2.0    0.1576  0.234     Fgamma. Its value here is
!      2.0     3.0    0.1772  0.266     understood as the MAXIMUM
!      2.0     4.0    0.1892  0.284     allowed. It is computed using
!      2.0     5.0    0.1976  0.296     a Newton iteration scheme.
!      2.0     6.0    0.2039  0.304
!      2.0     8.0    0.2129  0.314
!
!  NOTE: Theoretical values of Fgamma presented in the table above are
!  derived assuming "exact" barotropic mode stepping. Consequently, it
!  does not account for effects caused by Forward-Euler (FE) startup
!  of the barotropic mode at every 3D time step.  As the result, the
!  code may become unstable if the theoretical value of Fgamma is used
!  when mode splitting ratio "ndtfast" is small, thus yielding non-
!  negligible start up effects.  To compensate this, the accepted
!  value of Fgamma is reduced relatively to theoretical one, depending
!  on splitting ratio "ndtfast".  This measure is empirical. It is
!  shown to work with setting of "ndtfast" as low as 15, which is
!  more robust that the Hamming Window the squared cosine weights
!  options in "set_weights".
!
        real(dp) :: falpha = 2.0_dp
        real(dp) :: fbeta  = 4.0_dp
        real(dp) :: fgamma = 0.284_dp
#endif
!
!  Total number timesteps in current run. In 3D configurations, "ntimes"
!  is the total of baroclinic timesteps. In 2D configuration, "ntimes"
!  is the total of barotropic timesteps.
!
        integer, allocatable :: ntimes(:)
!
!  Time-step counter for current execution time-window.
!
        integer, allocatable :: step_counter(:)
 
!$OMP THREADPRIVATE (step_counter)
!
!  Number of time interval divisions for Stochastic Optimals.  It must
!  a multiple of "ntimes".
!
        integer :: nintervals = 1
!
!  Starting, current, and ending ensemble run parameters.
!
        integer :: erstr = 1                    ! Starting value
        integer :: erend = 1                    ! Ending value
        integer :: ninner = 1                   ! number of inner loops
        integer :: nouter = 1                   ! number of outer loops
        integer :: nrun = 1                     ! Current counter
#ifdef SENSITIVITY_4DVAR
        integer :: nrunsave = 0                 ! Loop counter
#endif
        integer :: outerloop = 0                ! split outer loop
        integer :: inner = 0                    ! inner loop counter
        integer :: outer = 0                    ! outer loop counter
!
! Set checksum algorithm: "adler32" or "crc32"
!
        character (len=*), parameter :: hashmethod = "bitsum"
!
! Split 4D-Var phase.
!
        character (len=20) :: phase4dvar
!
! Number of sadde point 4D-Var intervals.
!
        integer :: nsaddle = 1
!
!  First, starting, and ending timestepping parameters
!
        integer, allocatable :: ntfirst(:)      ! Forward-Euler step
        integer, allocatable :: ntstart(:)      ! Start step
        integer, allocatable :: ntend(:)        ! End step
!!$OMP THREADPRIVATE (ntfirst, ntstart, ntend)
!
!  Adjoint model or tangent linear model impulse forcing time record
!  counter and number of records available.
!
        integer, allocatable :: frcrec(:)
!$OMP THREADPRIVATE (FrcRec)
 
        integer, allocatable :: nrecfrc(:)
!
!  HSIMT tracer advection coefficients for the TVD limiter (Wu and Zhu,
!  2010).
!
        real(r8) :: cc1 = 0.25_r8
        real(r8) :: cc2 = 0.5_r8
        real(r8) :: cc3 = 1.0_r8/12.0_r8
!
!-----------------------------------------------------------------------
!  Control switches.
!-----------------------------------------------------------------------
!
!  Switch to use three ghost-points in the halo region.
!
        logical :: threeghostpoints = .false.
!
!  Switch to set-up application grid, metrics, and associated variables
!  and parameters.
!
        logical, allocatable :: setgridconfig(:)
!
!  Switch to procces (Get/Set) all input data in an application. It is
!  set to .TRUE. by default. It is manipulated by the ROMS-JEDI
!  interface during ROMS Phase 3 initialization that applies lateral
!  boundary condition to the initial state vector passed by OOPS. Some
!  applications requited input data for lateral boundary conditons.
!
        logical, allocatable :: processinputdata(:)
!
!  Switch to proccess nudging coefficients for radiation open boundary
!  conditions.
!
        logical, allocatable :: nudgingcoeff(:)
!
!  Switch to proccess input boundary data.
!
        logical, allocatable :: obcdata(:)
!
!  These switches are designed to control computational options within
!  nested and/or multiple connected grids.  They are .TRUE. by default.
!  They can turned off for a particular grind in input scripts.
!
        logical, allocatable :: lbiology(:)
        logical, allocatable :: lfloats(:)
        logical, allocatable :: lsediment(:)
        logical, allocatable :: lstations(:)
!
!  If equilibrium tides, switch to apply the 18.6-year lunar nodal
!  cycle correction.
!
        logical :: lnodal = .true.
!
!-----------------------------------------------------------------------
!  Physical constants.
!-----------------------------------------------------------------------
!
!    Cp            Specific heat for seawater (Joules/Kg/degC).
!    Csolar        Solar irradiantion constant (W/m2).
!    Eradius       Earth equatorial radius (m).
!    Infinity      Value resulting when dividing by zero.
!    StefBo        Stefan-Boltzmann constant (W/m2/K4).
!    emmiss        Infrared emissivity.
!    g             Acceleration due to gravity (m/s2).
!    gorho0        gravity divided by mean density anomaly.
!    rhow          fresh water density (kg/m3).
!    vonKar        von Karman constant.
!
        real(dp) :: cp = 3985.0_dp              ! Joules/kg/degC
        real(dp) :: csolar = 1353.0_dp          ! 1360-1380 W/m2
        real(dp) :: infinity                    ! Infinity = 1.0/0.0
        real(dp) :: eradius = 6371315.0_dp      ! m
        real(dp) :: stefbo = 5.67e-8_dp         ! Watts/m2/K4
        real(dp) :: emmiss = 0.97_dp            ! non_dimensional
        real(dp) :: rhow = 1000.0_dp            ! kg/m3
#ifdef SOLITON
        real(dp) :: g = 1.0_dp                  ! non-dimensional
#else
        real(dp) :: g = 9.81_dp                 ! m/s2
#endif
        real(dp) :: gorho0                      ! m4/s2/kg
        real(dp) :: vonkar = 0.41_dp            ! non-dimensional
!
!-----------------------------------------------------------------------
!  Various model parameters.  Some of these parameters are overwritten
!  with the values provided from model standard input script.
!-----------------------------------------------------------------------
!
!  Switch for spherical grid (lon,lat) configurations.
!
        logical :: spherical = .false.
!
!  Switch to indicate if ROMS kernel arrays have been allocated.  It
!  used to skip memory reports during early termination.
!
        logical :: lallocatedmemory = .false.
!$OMP THREADPRIVATE (LallocatedMemory)
!
!  Switch to compute the grid stiffness.
!
        logical :: lstiffness = .true.
!$OMP THREADPRIVATE (Lstiffness)
!
!  Composite grid a refined grids switches. They are .FALSE. by default.
!
        logical, allocatable :: compositegrid(:,:)
        logical, allocatable :: refinedgrid(:)
!
!  Refinement grid scale factor from donor grid.
!
        integer, allocatable :: refinescale(:)
!
!  Switch to extract donor grid (coarse) data at the refinement grid
!  contact point locations. The coarse data is extracted at the first
!  sub-refined time step.  Recall that the finer grid time-step is
!  smaller than the coarser grid by a factor of RefineScale(:). This
!  switch is only relevant during refinement nesting.
!
        logical, allocatable :: getdonordata(:)
!
!  Periodic boundary swiches for distributed-memory exchanges.
!
        logical, allocatable :: ewperiodic(:)
        logical, allocatable :: nsperiodic(:)
!
!  Lateral open boundary edges volume conservation switches.
!
        logical, allocatable :: volcons(:,:)
#if defined ADJOINT || defined TANGENT || defined TL_IOMS
        logical, allocatable :: ad_volcons(:,:)
        logical, allocatable :: tl_volcons(:,:)
#endif
!
!  Switches to read and process climatology fields.
!
        logical, allocatable :: clm_file(:)          ! Process NetCDF
        logical, allocatable :: lclimatology(:)      ! any field
        logical, allocatable :: lsshclm(:)           ! free-surface
        logical, allocatable :: lm2clm(:)            ! 2D momentum
        logical, allocatable :: lm3clm(:)            ! 3D momentum
        logical, allocatable :: ltracerclm(:,:)      ! tracers
!
!  Switched to nudge to climatology fields.
!
        logical, allocatable :: lnudging(:)          ! any field
        logical, allocatable :: lnudgem2clm(:)       ! 2D momentum
        logical, allocatable :: lnudgem3clm(:)       ! 3D momentum
        logical, allocatable :: lnudgetclm(:,:)      ! tracers
!
!  Switches to activate point Source/Sinks in an application:
!    * Horizontal momentum transport (u or v)
!    * Vertical mass transport (w)
!    * Tracer transport
!
        logical, allocatable :: luvsrc(:)            ! momentum
        logical, allocatable :: lwsrc(:)             ! mass
        logical, allocatable :: ltracersrc(:,:)      ! tracers
!
!  Execution termination flag.
!
!    exit_flag = 0   No error
!    exit_flag = 1   Blows up
!    exit_flag = 2   Input error
!    exit_flag = 3   Output error
!    exit_flag = 4   IO error
!    exit_flag = 5   Configuration error
!    exit_flag = 6   Partition error
!    exit_flag = 7   Illegal input parameter
!    exit_flag = 8   Fatal algorithm result
!    exit_flag = 9   coupling error
!
        integer :: exit_flag = 0
        integer :: blowup = 0
        integer :: noerror = 0
!
!  Blow-up string.
!
        character (len=80) :: blowup_string
!
!  Set threshold maximum speed (m/s) and density anomaly (kg/m3) to
!  test if the model is blowing-up.
!
        real(dp), allocatable :: maxspeed(:)
        real(dp), allocatable :: maxrho(:)
!
        real(dp) :: max_speed = 20.0_dp         ! m/s
        real(dp) :: max_rho = 200.0_dp          ! kg/m3
!
!  Interpolation scheme.
!
        integer, parameter :: linear = 0        ! linear interpolation
        integer, parameter :: cubic  = 1        ! cubic  interpolation
!
        integer :: interpflag = linear          ! interpolation flag
!
!  Shallowest and Deepest levels to apply bottom momentum stresses as
!  a bodyforce
!
        integer, allocatable :: levsfrc(:)
        integer, allocatable :: levbfrc(:)
!
!  Vertical coordinates transform.  Currently, there are two vertical
!  transformation equations (see set_scoord.F for details):
!
!    Original transform (Vtransform=1):
!
!         z_r(x,y,s,t) = Zo_r + zeta(x,y,t) * [1.0 + Zo_r / h(x,y)]
!
!                 Zo_r = hc * [s(k) - C(k)] + C(k) * h(x,y)
!
!    New transform (Vtransform=2):
!
!         z_r(x,y,s,t) = zeta(x,y,t) + [zeta(x,y,t)+ h(x,y)] * Zo_r
!
!                 Zo_r = [hc * s(k) + C(k) * h(x,y)] / [hc + h(x,y)]
!
        integer, allocatable :: vtransform(:)
!
!  Vertical grid stretching function flag:
!
!    Vstretcing = 1   Original function (Song and Haidvogel, 1994)
!               = 2   A. Shchepetkin (ROMS-UCLA) function
!               = 3   R. Geyer BBL function
!
        integer, allocatable :: vstretching(:)
!
!  Vertical grid stretching parameters.
!
!    Tcline        Width (m) of surface or bottom boundary layer in
!                    which higher vertical resolution is required
!                    during stretching.
!    hc            S-coordinate critical depth, hc=MIN(hmin,Tcline).
!    theta_s       S-coordinate surface control parameter.
!    theta_b       S-coordinate bottom control parameter.
!
        real(dp), allocatable :: tcline(:)      ! m, positive
        real(dp), allocatable :: hc(:)          ! m, positive
        real(dp), allocatable :: theta_s(:)     ! 0 < theta_s < 20
        real(dp), allocatable :: theta_b(:)     ! 0 < theta_b < 1
!
!  Bathymetry range values.
!
        real(dp), allocatable :: hmin(:)        ! m, positive
        real(dp), allocatable :: hmax(:)        ! m, positive
!
!  Length (m) of domain box in the XI- and ETA-directions.
!
        real(r8), allocatable :: xl(:)          ! m
        real(r8), allocatable :: el(:)          ! m
!
!  Minimum and Maximum longitude and latitude at RHO-points
!
        real(r8), allocatable :: lonmin(:)      ! degrees east
        real(r8), allocatable :: lonmax(:)      ! degrees east
        real(r8), allocatable :: latmin(:)      ! degrees north
        real(r8), allocatable :: latmax(:)      ! degrees north
!
!  Constant used in the Shchepetkin boundary conditions for 2D momentum,
!  Co = 1.0_r8/(2.0_r8+SQRT(2.0_r8)).
!
        real(r8) :: co
!
!  Number of digits in grid size for format statements.
!
        integer, allocatable :: idigits(:)
        integer, allocatable :: jdigits(:)
#ifdef SOLVE3D
        integer, allocatable :: kdigits(:)
#endif
!
!  Diagnostic volume averaged variables.
!
        integer, allocatable :: first_time(:)
 
        real(dp) :: avgke = 0.0_dp              ! Kinetic energy
        real(dp) :: avgpe = 0.0_dp              ! Potential energy
        real(dp) :: avgkp = 0.0_dp              ! Total energy
        real(dp) :: volume = 0.0_dp             ! diagnostics volume
        real(dp) :: ad_volume = 0.0_dp          ! adjoint volume
 
        real(dp), allocatable :: totvolume(:)   ! Total volume
        real(dp), allocatable :: minvolume(:)   ! Minimum cell volume
        real(dp), allocatable :: maxvolume(:)   ! Maximum cell volume
!
!  Minimun and maximum grid spacing
!
        real(dp), allocatable :: dxmin(:)       ! all grid points
        real(dp), allocatable :: dxmax(:)
        real(dp), allocatable :: dymin(:)
        real(dp), allocatable :: dymax(:)
#ifdef MASKING
        real(dp), allocatable :: dxminw(:)      ! only water points
        real(dp), allocatable :: dxmaxw(:)
        real(dp), allocatable :: dyminw(:)
        real(dp), allocatable :: dymaxw(:)
#endif
#ifdef SOLVE3D
        real(dp), allocatable :: dzmin(:)       ! all grid points
        real(dp), allocatable :: dzmax(:)
# ifdef MASKING
        real(dp), allocatable :: dzminw(:)      ! only water points
        real(dp), allocatable :: dzmaxw(:)
# endif
#endif
!
!  Maximum size of a grid node (m) over the whole curvilinear grid
!  application. Used for scaling horizontal mixing by the grid size.
!
        real(dp), allocatable :: grdmax(:)
#ifdef DIFF_3DCOEF
        real(dp), allocatable :: diffmin(:)     ! Minimun diffusion
        real(dp), allocatable :: diffmax(:)     ! Maximum diffusion
#endif
#ifdef VISC_3DCOEF
        real(dp), allocatable :: viscmin(:)     ! Minimum viscosity
        real(dp), allocatable :: viscmax(:)     ! Maximum viscosity
#endif
!
!  Courant Numbers due to gravity wave speed limits.
!
        real(dp), allocatable :: cg_min(:)      ! Minimun barotropic
        real(dp), allocatable :: cg_max(:)      ! Maximun barotropic
        real(dp), allocatable :: cg_cor(:)      ! Maximun Coriolis
!
!  Time dependent Counrant Numbers due to velocity components and
!  indices location of maximum value.
!
        integer :: max_ci = 0                   ! maximum I-location
        integer :: max_cj = 0                   ! maximum J-location
        integer :: max_ck = 0                   ! maximum K-location
        real(r8) :: max_c = 0.0_r8              ! maximum total
        real(r8) :: max_cu = 0.0_r8             ! maximum I-component
        real(r8) :: max_cv = 0.0_r8             ! maximum J-component
#ifdef SOLVE3D
        real(r8) :: max_cw = 0.0_r8             ! maximum K-component
#endif
!
!  Linear equation of state parameters.
!
!    R0            Background constant density anomaly (kg/m3).
!    Tcoef         Thermal expansion coefficient (1/Celsius).
!    Scoef         Saline contraction coefficient (1/PSU).
!
        real(r8), allocatable :: r0(:)
        real(r8), allocatable :: tcoef(:)
        real(r8), allocatable :: scoef(:)
!
!  Background potential temperature (Celsius) and salinity (PSU) values
!  used in analytical initializations.
!
        real(r8), allocatable :: t0(:)
        real(r8), allocatable :: s0(:)
!
!  Slipperiness variable, either 1.0 (free slip) or -1.0 (no slip).
!
        real(r8), allocatable :: gamma2(:)
!
!  Weighting coefficient for the newest (implicit) time step derivatives
!  using either a Crack-Nicolson implicit scheme (lambda=0.5) or a
!  backward implicit scheme (lambda=1.0).
!
#if defined SPLINES_VDIFF || defined SPLINES_VVISC
        real(r8) :: lambda = 1.0_r8
#else
!!      real(r8) :: lambda = 0.5_r8
        real(r8) :: lambda = 1.0_r8
#endif
!
!  Jerlov water type to assign everywhere, range values: 1 - 5.
!
        integer, allocatable :: lmd_jwt(:)
!
!  Grid r-factor (non-dimensional).
!
        real(dp), allocatable :: rx0(:)         ! Beckmann and Haidvogel
        real(dp), allocatable :: rx1(:)         ! Haney
!
!  Linear (m/s) and quadratic (nondimensional) bottom drag coefficients.
!
        real(r8), allocatable :: rdrg(:)
        real(r8), allocatable :: rdrg2(:)
!
!  Minimum and maximum threshold for transfer coefficient of momentum.
!
        real(dp) :: cdb_min = 0.000001_dp
        real(dp) :: cdb_max = 0.5_dp
!
!  Surface and bottom roughness (m)
!
        real(r8), allocatable :: zos(:)
        real(r8), allocatable :: zob(:)
!
!  Minimum depth for wetting and drying (m).
!
        real(r8), allocatable :: dcrit(:)
!
!  Mean density (Kg/m3) used when the Boussinesq approximation is
!  inferred.
!
        real(dp) :: rho0 = 1025.0_dp
!
!  Background Brunt-Vaisala frequency (1/s2).
!
        real(dp) :: bvf_bak = 0.00001_dp
 
#ifdef PROPAGATOR
!
!  Number of converged Ritz values and relative accuracy of computed
!  Ritz values.
!
        integer, allocatable :: nconv(:)
        real(dp) :: ritz_tol = 1.0e-15_dp
#endif
!
!  Vector containing USER generic parameters.
!
        integer :: nuser
        real(r8), allocatable :: user(:)
!
!  Weights for the time average of 2D fields.
!
        real(dp), allocatable :: weight(:,:,:)
!
!  Constants.
!
        real(dp), parameter :: pi = 3.14159265358979323846_dp
        real(dp), parameter :: deg2rad = pi / 180.0_dp
        real(dp), parameter :: rad2deg = 180.0_dp / pi
        real(dp), parameter :: day2sec = 86400.0_dp
        real(dp), parameter :: sec2day = 1.0_dp / 86400.0_dp
#ifdef NO_4BYTE_REALS
        real(dp), parameter :: spval = nf90_fill_double
#else
        real(dp), parameter :: spval = 1.0e+37_dp
#endif
        real(dp), parameter :: large = 1.0e+20_dp
        real(dp), parameter :: jul_off = 2440000.0_dp
!
!  Set special check value.  Notice that a smaller value is assigned
!  to account for both NetCDF fill value and roundoff. There are
!  many Matlab scripts out there that do not inquire correctly
!  the spval from the _FillValue attribute in single/double
!  precision.
!
        real(dp), parameter :: spval_check = 1.0e+35_dp
!
!-----------------------------------------------------------------------
!  Horizontal and vertical constant mixing coefficients.
!-----------------------------------------------------------------------
!
!    Akk_bak       Background vertical mixing coefficient (m2/s) for
!                    turbulent energy.
!    Akp_bak       Background vertical mixing coefficient (m2/s) for
!                    generic statistical field "psi".
!    Akt_bak       Background vertical mixing coefficient (m2/s) for
!                    tracers.
!    Akv_bak       Background vertical mixing coefficient (m2/s) for
!                    momentum.
!    Akt_limit     Upper threshold vertical mixing coefficient (m2/s)
!                    for tracers.
!    Akv_limit     Upper threshold vertical mixing coefficient (m2/s)
!                    for momentum.
!    Kdiff         Isopycnal mixing thickness diffusivity (m2/s) for
!                    tracers.
!    ad_visc2      ADM lateral harmonic constant mixing coefficient
!                    (m2/s) for momentum.
!    nl_visc2      NLM lateral harmonic constant mixing coefficient
!                    (m2/s) for momentum.
!    tl_visc2      TLM lateral harmonic constant mixing coefficient
!                    (m2/s) for momentum.
!    visc2         Current lateral harmonic constant mixing coefficient
!                    (m2/s) for momentum.
!    ad_visc4      ADM lateral biharmonic (squared root) constant
!                     mixing coefficient (m2 s^-1/2) for momentum.
!    nl_visc4      NLM lateral biharmonic (squared root) constant
!                     mixing coefficient (m2 s^-1/2) for momentum.
!    tl_visc4      TLM lateral biharmonic (squared root) constant
!                     mixing coefficient (m2 s^-1/2) for momentum.
!    visc4         Current lateral biharmonic (squared root) constant
!                     mixing coefficient (m2 s^-1/2) for momentum.
!    ad_tnu2       ADM lateral harmonic constant mixing coefficient
!                    (m2/s) for tracer type variables.
!    nl_tnu2       NLM lateral harmonic constant mixing coefficient
!                    (m2/s) for tracer type variables.
!    tl_tnu2       TLM lateral harmonic constant mixing coefficient
!                    (m2/s) for tracer type variables.
!    tnu2          Current lateral harmonic constant mixing coefficient
!                    (m2/s) for tracer type variables.
!    ad_tnu4       ADM lateral biharmonic (squared root) constant
!                     mixing coefficient (m2 s^-1/2) for tracers.
!    nl_tnu4       NLM lateral biharmonic (squared root) constant
!                     mixing coefficient (m2 s^-1/2) for tracers.
!    tl_tnu4       TLM lateral biharmonic (squared root) constant
!                     mixing coefficient (m2 s^-1/2) for tracers.
!    tnu4          Current lateral biharmonic (squared root) constant
!                     mixing coefficient (m2 s^-1/2) for tracers.
!    tkenu2        Lateral harmonic constant mixing coefficient
!                    (m2/s) for turbulent energy.
!    tkenu4        Lateral biharmonic (squared root) constant mixing
!                    coefficient (m2 s^-1/2) for turbulent energy.
!
        real(r8), allocatable :: akk_bak(:)          ! m2/s
        real(r8), allocatable :: akp_bak(:)          ! m2/s
        real(r8), allocatable :: akv_bak(:)          ! m2/s
        real(r8), allocatable :: akv_limit(:)        ! m2/s
 
        real(r8), allocatable :: ad_visc2(:)         ! m2/s
        real(r8), allocatable :: nl_visc2(:)         ! m2/s
        real(r8), allocatable :: tl_visc2(:)         ! m2/s
        real(r8), allocatable :: visc2(:)            ! m2/s
 
        real(r8), allocatable :: ad_visc4(:)         ! m2 s-1/2
        real(r8), allocatable :: nl_visc4(:)         ! m2 s-1/2
        real(r8), allocatable :: tl_visc4(:)         ! m2 s-1/2
        real(r8), allocatable :: visc4(:)            ! m2 s-1/2
 
        real(r8), allocatable :: tkenu2(:)           ! m2/s
        real(r8), allocatable :: tkenu4(:)           ! m2 s-1/2
 
        real(r8), allocatable :: akt_bak(:,:)        ! m2/s
        real(r8), allocatable :: akt_limit(:,:)      ! m2/s
        real(r8), allocatable :: kdiff(:,:)          ! m2/s
 
        real(r8), allocatable :: ad_tnu2(:,:)        ! m2/s
        real(r8), allocatable :: nl_tnu2(:,:)        ! m2/s
        real(r8), allocatable :: tl_tnu2(:,:)        ! m2/s
        real(r8), allocatable :: tnu2(:,:)           ! m2/s
 
        real(r8), allocatable :: ad_tnu4(:,:)        ! m2 s-1/2
        real(r8), allocatable :: nl_tnu4(:,:)        ! m2 s-1/2
        real(r8), allocatable :: tl_tnu4(:,:)        ! m2 s-1/2
        real(r8), allocatable :: tnu4(:,:)           ! m2 s-1/2
!
!  Horizontal diffusive relaxation coefficients (m2/s) used to smooth
!  representer tangent linear solution during Picard iterations to
!  improve stability and convergence.
!
        real(r8), allocatable :: tl_m2diff(:)        ! 2D momentum
        real(r8), allocatable :: tl_m3diff(:)        ! 3D momentum
 
        real(r8), allocatable :: tl_tdiff(:,:)       ! tracers
!
!  Basic state vertical mixing coefficient scale factors for adjoint
!  based algorithms. In some applications, a smaller/larger values of
!  vertical mixing are necessary for stability.
!
        real(r8), allocatable :: ad_akv_fac(:)       ! ADM momentum
        real(r8), allocatable :: tl_akv_fac(:)       ! TLM momentum
 
        real(r8), allocatable :: ad_akt_fac(:,:)     ! ADM tracers
        real(r8), allocatable :: tl_akt_fac(:,:)     ! TLM tracers
 
!
!  Switches to increase/decrease horizontal viscosity and/or diffusion
!  in specific areas of the application domain (like sponge areas).
!
        logical, allocatable :: lsponge(:)
        logical, allocatable :: luvsponge(:)         ! viscosity
        logical, allocatable :: ltracersponge(:,:)   ! diffusion
!
!-----------------------------------------------------------------------
!  IO parameters.
!-----------------------------------------------------------------------
!
!  Switches to activate creation and writing of output NetCDF files.
!
        logical, allocatable :: ldefadj(:)       ! Adjoint file
        logical, allocatable :: ldefavg(:)       ! Average file
        logical, allocatable :: ldefdai(:)       ! DA initial/restart
        logical, allocatable :: ldefdia(:)       ! Diagnostics file
        logical, allocatable :: ldeferr(:)       ! 4DVar error file
        logical, allocatable :: ldefflt(:)       ! Floats file
        logical, allocatable :: ldefhis(:)       ! History file
        logical, allocatable :: ldefhss(:)       ! Hessian file
        logical, allocatable :: ldefini(:)       ! Initial file
        logical, allocatable :: ldefirp(:)       ! Initial RPM file
        logical, allocatable :: ldefitl(:)       ! Initial TLM file
        logical, allocatable :: ldeflcz(:)       ! Lanczos Vectors file
        logical, allocatable :: ldeflze(:)       ! Evolved Lanczos file
        logical, allocatable :: ldefmod(:)       ! 4DVAR file
        logical, allocatable :: ldefqck(:)       ! Quicksave file
        logical, allocatable :: ldefrst(:)       ! Restart file
#ifdef SP4DVAR
        logical, allocatable :: ldefsca(:)       ! AD scratch state file
        logical, allocatable :: ldefsct(:)       ! TL scratch state file
        logical, allocatable :: ldefspa(:)       ! AD Arnoldi state file
        logical, allocatable :: ldefspt(:)       ! TL Arnoldi state file
#endif
        logical, allocatable :: ldefsta(:)       ! Stations file
#ifdef STD_MODEL
        logical, allocatable :: ldefstd(:)       ! STD file
#endif
        logical, allocatable :: ldeftide(:)      ! tide forcing file
        logical, allocatable :: ldeftlm(:)       ! Tangent linear file
        logical, allocatable :: ldeftlf(:)       ! TLM/RPM impulse file
        logical, allocatable :: ldefxtr(:)       ! TLM/RPM impulse file
 
        logical, allocatable :: lreadadm(:)      ! Read ADM multi-files
        logical, allocatable :: lreadblk(:)      ! Read NLM bulk fluxes
        logical, allocatable :: lreadfrc(:)      ! Read FRC files
        logical, allocatable :: lreadfwd(:)      ! Read FWD trajectory
        logical, allocatable :: lreadqck(:)      ! Read QCK trajectory
        logical, allocatable :: lreadstd(:)      ! Read STD file
        logical, allocatable :: lreadtlm(:)      ! Read TLM multi-files
 
        logical, allocatable :: lwrtadj(:)       ! Write adjoint file
        logical, allocatable :: lwrtavg(:)       ! Write average file
        logical, allocatable :: lwrtdia(:)       ! Write diagnostic file
        logical, allocatable :: lwrthis(:)       ! Write history file
        logical, allocatable :: lwrtper(:)       ! Write during ensemble
        logical, allocatable :: lwrtqck(:)       ! write quicksave file
        logical, allocatable :: lwrtrst(:)       ! Write restart file
#ifdef STD_MODEL
        logical, allocatable :: lwrtstd(:)       ! Write STD file
#endif
        logical, allocatable :: lwrttlf(:)       ! Write impulse file
        logical, allocatable :: lwrttlm(:)       ! Write tangent file
        logical, allocatable :: lwrtxtr(:)       ! Write extraction file
 
        logical, allocatable :: ldefnrm(:,:)     ! Norm file
        logical, allocatable :: lwrtnrm(:,:)     ! Write norm file
!
!  Switch to write out adjoint 2D state arrays instead of IO solution
!  arrays and adjoint ocean time. This is used in 4DVAR for IO
!  maniputations.
!
#ifdef AD_OUTPUT_STATE
        logical, allocatable :: lwrtstate3d(:)
#endif
#if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
        logical, allocatable :: soinitial(:)
#endif
        logical, allocatable :: lwrtstate2d(:)
        logical, allocatable :: lwrttime(:)
!
!  Switch to write 4D-Var cost functions to DAV NetCDF.
!
        logical, allocatable :: lwrtcost(:)
!
!  Switch to write out adjoint surface forcing fields adjusted by the
!  4DVAR algorithms.
!
        logical, allocatable :: ladjusted(:)
!
!  Switch to append information to an existing ROMS standard output
!  log file.
!
        logical :: lappend = .false.
!
!  Switch to read input open boundary conditions data.
!
        logical, allocatable :: lprocessobc(:)
!
!  Switch to read input tidal forcing data.
!
        logical, allocatable :: lprocesstides(:)
!
!  Switch to write application set-up information to standard output.
!
        logical, allocatable :: lwrtinfo(:)
!
!  Switch used to create new output NetCDF files. If TRUE, new output
!  files are created. If FALSE, data is appended to an existing output
!  files.  Used only for history, average and station files.
!
        logical, allocatable :: ldefout(:)       ! New output files
!
!  Number of timesteps between creation of new output files.
!
        integer, allocatable :: ndefadj(:)       ! Adjoint file
        integer, allocatable :: ndefavg(:)       ! Average file
        integer, allocatable :: ndefdia(:)       ! Diagnostics file
        integer, allocatable :: ndefhis(:)       ! History file
        integer, allocatable :: ndefqck(:)       ! Quicksave file
        integer, allocatable :: ndeftlm(:)       ! Tangent linear file
        integer, allocatable :: ndefxtr(:)       ! extraction file
!
!  Starting timestep for accumulation of output.
!
        integer, allocatable :: ntsavg(:)        ! Average file
        integer, allocatable :: ntsdia(:)        ! Diagnostics file
!
!  Number of timesteps between writing of output data.
!
        integer, allocatable :: nadj(:)          ! Adjoint file
        integer, allocatable :: navg(:)          ! Average file
        integer, allocatable :: ndia(:)          ! Diagnostics file
        integer, allocatable :: nflt(:)          ! Floats file
        integer, allocatable :: nhis(:)          ! History file
        integer, allocatable :: nqck(:)          ! Quicksave file
        integer, allocatable :: nrst(:)          ! Restart file
        integer, allocatable :: nsta(:)          ! Stations file
        integer, allocatable :: ntlm(:)          ! Tangent linear file
        integer, allocatable :: nxtr(:)          ! extraction file
#ifdef SP4DVAR
        integer, allocatable :: nsca(:)          ! AD scratch state file
        integer, allocatable :: nsct(:)          ! TL scratch state file
        integer, allocatable :: nspa(:)          ! AD Arnoldi state file
        integer, allocatable :: nspt(:)          ! TL Arnoldi state file
        integer, allocatable :: crec(:)          ! Convolution records
#endif
!
!  Field extraction flag to interpolate or decimate solution to the
!  provided grid geometry.
!  (For 4D-Var coarser inner loops use ExtractFlag = 2)
!
!    ExtractFlag = 0    no extraction
!    ExtractFlag = 1    extraction by interpolation
!    ExtractFlag > 1    extraction by decimation
!
        integer, allocatable :: extractflag(:)
!
!  Number of timesteps between print of single line information to
!  standard output.
!
        integer, allocatable :: ninfo(:)
!
!  Number of timesteps between 4DVAR adjustment of open boundaries.
!  In strong constraint 4DVAR, it is possible to open bounadies at
!  other intervals in addition to initial time. These parameters are
!  used to store the appropriate number of open boundary records in
!  output history NetCDF files.
!
!    Nbrec(:) = 1 + ntimes(:) / nOBC(:)
!
!  Here, it is assumed that nOBC is a multiple of NTIMES or greater
!  than NTIMES. If nOBC > NTIMES, only one record is stored in the
!  output history NetCDF files and the adjustment is for constant
!  open boundaries with constant correction.
!
        integer, allocatable :: nobc(:)          ! number of timesteps
        integer, allocatable :: nbrec(:)         ! number of records
        integer, allocatable :: obccount(:)      ! record counter
 
#ifdef ADJUST_BOUNDARY
!
!  Logical switches to process open boundary arrays during 4DVar
!  adjustments.
!
        logical, allocatable :: lobc(:,:,:)
!
!  Time (s) of surface forcing adjustment.
!
        real(dp), allocatable :: obc_time(:,:)
#endif
!
!  Number of timesteps between adjustment of 4DVAR surface forcing
!  fields. In strong constraint 4DVAR, it is possible to adjust surface
!  forcing fields at other intervals in addition to initial time.
!  These parameters are used to store the appropriate number of
!  surface forcing records in output history NetCDF files.
!
!    Nfrec(:) = 1 + ntimes(:) / nSFF(:)
!
!  Here, it is assumed that nSFF is a multiple of NTIMES or greater
!  than NTIMES. If nSFF > NTIMES, only one record is stored in the
!  output history NetCDF files and the adjustment is for constant
!  forcing with constant correction.
!
        integer, allocatable :: nsff(:)          ! number of timesteps
        integer, allocatable :: nfrec(:)         ! number of records
        integer, allocatable :: sfcount(:)       ! record counter
 
#ifdef ADJUST_STFLUX
!
!  Logical switches to process surface tracer fluxes during 4DVar
!  adjustments.
!
        logical, allocatable :: lstflux(:,:)
#endif
#if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
!
!  Time (s) of surface forcing adjustment.
!
        real(dp), allocatable :: sf_time(:,:)
#endif
!
!  Restart time record to read from disk and use as the initial
!  conditions. Use nrrec=0 for new solutions. If nrrec is negative
!  (say, nrrec=-1), the model will restart from the most recent
!  time record. That is, the initialization record is assigned
!  internally.
!
        integer, allocatable :: nrrec(:)
!
!  Switch to activate processing of input data.  This switch becomes
!  very useful when reading input data serially in parallel
!  applications.
!
        logical, allocatable :: synchro_flag(:)
!$OMP THREADPRIVATE (synchro_flag)
!
!  Switch to inialize model with latest time record from initial
!  (restart/history) NetCDF file.
!
        logical, allocatable :: lastrec(:)
!
!  Generalized Statbility Theory (GST) parameters.
!
        logical :: lmultigst          ! multiple eigenvector file switch
        logical :: lrstgst            ! restart switch
        integer :: maxitergst         ! Number of iterations
        integer :: ngst               ! check pointing interval
!
!  Switches used to recycle time records in some output file. If TRUE,
!  only the latest two time records are maintained.  If FALSE, all
!  field records are saved.
!
        logical, allocatable :: lcycleadj(:)
        logical, allocatable :: lcyclerst(:)
        logical, allocatable :: lcycletlm(:)
 
#if defined AVERAGES  && defined AVERAGES_DETIDE && \
   (defined ssh_tides || defined uv_tides)
!
!  Counter storing the number of accumulated harmonic records used
!  for detiding.
!
        integer, allocatable :: hcount(:)
#endif
!
!-----------------------------------------------------------------------
!  Adjoint sensitivity parameters.
!-----------------------------------------------------------------------
!
!  Starting and ending vertical levels of the 3D adjoint state whose
!  sensitivity is required.
!
        integer, allocatable :: kstrs(:)           ! starting level
        integer, allocatable :: kends(:)           ! ending level
!
!  Starting and ending day for adjoint sensitivity forcing.
!
        real(r8), allocatable :: dstrs(:)          ! starting day
        real(r8), allocatable :: dends(:)          ! ending day
!
!-----------------------------------------------------------------------
!  Stochastic optimals parameters.
!-----------------------------------------------------------------------
!
!  Stochastic optimals forcing records counter.
!
        integer, allocatable :: sorec(:)
!$OMP THREADPRIVATE (SOrec)
!
!  Trace of stochastic optimals matrix.
!
        real(r8), allocatable :: trnorm(:)
!
!  Stochastic optimals time decorrelation scale (days) assumed for
!  red noise processes.
!
        real(r8), allocatable :: so_decay(:)
!
!  Stochastic optimals surface forcing standard deviation for
!  dimensionalization.
!
        real(r8), allocatable :: so_sdev(:,:)
 
#if defined FOUR_DVAR || defined VERIFICATION
!
!------------------------------------------------------------------------
!  Background/model error covariance parameters.
!------------------------------------------------------------------------
!
!  Maximum number of model state variables to process.
!
        integer :: mstatevar
!
!  Logical switch to compute initial conditions, model and surface
!  forcing error covariance normalization factors.
!
        logical,  allocatable :: cnorm(:,:)
!
!  Logical switch to compute boundary conditions error covariance
!  normalization factors.
!
        logical,  allocatable :: cnormb(:,:)
!
!  Logical switches to process weak constraint forcing as intermittent
!  or continuous impulses.
!
        logical, allocatable :: sporadicimpulse(:)     ! intermittent
        logical, allocatable :: frequentimpulse(:)     ! continuous
!
!  Logical switch to control 4D-Var incremental analysis update forcing.
!
        logical, allocatable :: iauswitch(:)
!
!  Duration of the 4D-Var incremental inalysis update (seconds).
!
        real(dp), allocatable :: timeiau(:)
!
!  Stability and accuracy factor used to scale the time-step of the
!  horizontal and vertical convolution operator below its theoretical
!  (CFL) limit. Notice that four values are needed for these factors to
!  facilitate the error covariance modeling for initial conditions (1),
!  model (2), boundary conditions (3), and surface forcing (4).
!
        real(r8), dimension(4) :: hgamma
        real(r8), dimension(4) :: vgamma
!
!  Parameters used to compute balanced salinity in terms of temperature
!  using empirical T-S relationships in error covariance balance
!  operator.
!
        real(r8), allocatable :: dtdz_min(:)     ! minimum dT/ds (C/m)
        real(r8), allocatable :: ml_depth(:)     ! mixed-layer depth (m)
!
!  Balance operator level of no motion depth (m; positive) used when
!  computing balanced free-surface contribution.
!
        real(r8), allocatable :: lnm_depth(:)
!
! Balance operator level of no motion flag used to compute balanced
! free-surface contribution:
!
!      [0] Integrate from local bottom to the surface
!      [1] Integrate from LNM_depth to surface or integrate from local bottom
!          if shallower than LNM_depth
!
        integer :: lnm_flag
!
! Balance operator logical switches for state variables to consider in the
! error covariance multivariate constraints.
!
        logical, allocatable :: balance(:)
!
!  Initial conditions, model and surface forcing error covariance
!  horizontal decorrelation scales (m).
!
        real(r8), allocatable :: hdecay(:,:,:)
!
!  Initial conditions, model and surface forcing error covariance
!  temporal decorrelation scales (second).
!
        real(r8), allocatable :: tdecay(:,:)
!
!  Initial conditions, model and surface forcing error covariance
!  vertical decorrelation scales (m).
!
        real(r8), allocatable :: vdecay(:,:,:)
!
!  Boundary conditions error covariance horizontal decorrelation
!  scales (m).
!
        real(r8), allocatable :: hdecayb(:,:,:)
!
!  Boundary conditions error covariance vertical decorrelation
!  scales (m).
!
        real(r8), allocatable :: vdecayb(:,:,:)
!
!  Method for background quality control of observations, [Ngrids].
!
!      [1] Quality control in terms of state variable indices
!      [2] Quality control in terms of observation provenance
!
        integer, allocatable :: bgqc_type(:)
!
!  Background quality control standard deviation value for not
!  rejection of observations.
!
        real(r8) :: bgqc_large = 1.0e+5_r8
!
!  Number of observation provenances used in background quality control,
!  [Ngrids]. Only used when bgqc_type(ng)=2.
!
        integer, allocatable :: nprovenance(:)
!
!  Observation provenance indices to process during background quality
!  control of observations, [MAXVAL(Nprovenance),Ngrids].
!
        integer, allocatable :: iprovenance(:,:)
!
!  Background quality control threshold standard deviations in terms
!  of observation provenance indices, [MAXVAL(Nprovenance),Ngrids].
!
        real(r8), allocatable :: p_bgqc(:,:)
!
!  Background quality control threshold standard deviations in terms
!  of state variable indices, [MstateVar,Ngrids]
!
        real(r8), allocatable :: s_bgqc(:,:)
 
# ifdef STD_MODEL
!
!  Background error covariance, standard deviation modeling parameters.
!
!  This approach computes the standard deviation (STD) directly from
!  the background (prior) field as an alternative to climatological
!  values read from the input NetCDF files. It follows the work of
!  Mogensen et al. (2012) by assuming the background errors are
!  proportional to the vertical derivatives of the background field.
!  The field error has a similar profile shape, but the difference
!  with the actual error value is due to due to a vertical displacemnt.
!  Thus, the parameters below, per state variable, are used to process
!  the STD profile in terms od the location of the mixed layer where
!  the error is larger.
!
        real(r8), allocatable :: sigma_max(:,:) ! Maximum value
        real(r8), allocatable :: sigma_ml(:,:)  ! Minimum at mixed layer
        real(r8), allocatable :: sigma_do(:,:)  ! Minimum in deep ocean
        real(r8), allocatable :: sigma_dz(:,:)  ! vetical displacement
!
!  The computation of the mixed-layer depth is computed using the
!  approach of Kara et al. (2000), if activated. Otherwiae, a constant
!  value is used.
!
        real(r8), allocatable :: mld_uniform(:) ! mixed layer depth
# endif
#endif
!
!-----------------------------------------------------------------------
!  Nudging variables for passive (outflow) and active (inflow) oepn
!  boundary conditions.
!-----------------------------------------------------------------------
!
!    iwest         West  identification index in boundary arrays.
!    isouth        South identification index in boundary arrays.
!    ieast         East  identification index in boundary arrays.
!    inorth        North identification index in boundary arrays.
!    obcfac        Factor between passive and active open boundary
!                    conditions (nondimensional and greater than one).
!                    The nudging time scales for the active conditions
!                    are obtained by multiplying the passive values by
!                    factor.
!    FSobc_in      Active and strong time-scale (1/sec) coefficients
!                    for nudging towards free-surface data at  inflow.
!    FSobc_out     Passive and weak  time-scale (1/sec) coefficients
!                    for nudging towards free-surface data at outflow.
!    M2obc_in      Active and strong time-scale (1/sec) coefficients
!                    for nudging towards 2D momentum data at  inflow.
!    M2obc_out     Passive and weak  time-scale (1/sec) coefficients
!                    for nudging towards 2D momentum data at outflow.
!    M3obc_in      Active and strong time-scale (1/sec) coefficients
!                    for nudging towards 3D momentum data at  inflow.
!    M3obc_out     Passive and weak  time-scale (1/sec) coefficients
!                    for nudging towards 3D momentum data at outflow.
!    Tobc_in       Active and strong time-scale (1/sec) coefficients
!                    for nudging towards tracer data at  inflow.
!    Tobc_out      Passive and weak  time-scale (1/sec) coefficients
!                    for nudging towards tracer data at outflow.
!
        integer, parameter :: iwest = 1
        integer, parameter :: isouth = 2
        integer, parameter :: ieast = 3
        integer, parameter :: inorth = 4
 
        real(dp), allocatable :: obcfac(:)
        real(dp), allocatable :: fsobc_in(:,:)
        real(dp), allocatable :: fsobc_out(:,:)
        real(dp), allocatable :: m2obc_in(:,:)
        real(dp), allocatable :: m2obc_out(:,:)
#ifdef SOLVE3D
        real(dp), allocatable :: m3obc_in(:,:)
        real(dp), allocatable :: m3obc_out(:,:)
        real(dp), allocatable :: tobc_in(:,:,:)
        real(dp), allocatable :: tobc_out(:,:,:)
#endif
!
!  Inverse time-scales (1/s) for nudging at open boundaries and sponge
!  areas.
!
        real(dp), allocatable :: znudg(:)          ! Free-surface
        real(dp), allocatable :: m2nudg(:)         ! 2D momentum
        real(dp), allocatable :: m3nudg(:)         ! 3D momentum
        real(dp), allocatable :: tnudg(:,:)        ! Tracers
!
!  Variables used to impose mass flux conservation in open boundary
!  configurations.
!
        real(dp) :: bc_area = 0.0_dp
        real(dp) :: bc_flux = 0.0_dp
        real(dp) :: ubar_xs = 0.0_dp
#if defined TANGENT || defined TL_IOMS
        real(dp) :: tl_bc_area = 0.0_dp
        real(dp) :: tl_bc_flux = 0.0_dp
        real(dp) :: tl_ubar_xs = 0.0_dp
#endif
#ifdef ADJOINT
        real(dp) :: ad_bc_area = 0.0_dp
        real(dp) :: ad_bc_flux = 0.0_dp
        real(dp) :: ad_ubar_xs = 0.0_dp
#endif
#ifdef SP4DVAR
        real(r8) :: sp_sumdot = 0.0_r8
#endif
 
#ifdef BULK_FLUXES
!
!-----------------------------------------------------------------------
! Constants used in surface fluxes bulk parameterization.
!-----------------------------------------------------------------------
!
!    blk_Cpa       Specific heat capacity for dry air (J/kg/K).
!    blk_Cpw       Specific heat capacity for seawater (J/kg/K).
!    blk_Rgas      Gas constant for dry air (J/kg/K).
!    blk_Zabl      Height (m) of atmospheric boundary layer.
!    blk_ZQ        Height (m) of surface air humidity measurement.
!    blk_ZT        Height (m) of surface air temperature measurement.
!    blk_ZW        Height (m) of surface winds measurement.
!    blk_beta      Beta parameter evaluated from Fairall low windspeed
!                    turbulence data.
!    blk_dter      Temperature change.
!    blk_tcw       Thermal conductivity of water (W/m/K).
!    blk_visw      Kinematic viscosity water (m2/s).
!
        real(dp) :: blk_cpa = 1004.67_dp      ! (J/kg/K), Businger 1982
        real(dp) :: blk_cpw = 4000.0_dp       ! (J/kg/K)
        real(dp) :: blk_rgas = 287.1_dp       ! (J/kg/K)
        real(dp) :: blk_zabl = 600.0_dp       ! (m)
        real(dp) :: blk_beta = 1.2_dp         ! non-dimensional
        real(dp) :: blk_dter = 0.3_dp         ! (K)
        real(dp) :: blk_tcw = 0.6_dp          ! (W/m/K)
        real(dp) :: blk_visw = 0.000001_dp    ! (m2/s)
 
        real(r8), allocatable :: blk_zq(:)        ! (m)
        real(r8), allocatable :: blk_zt(:)        ! (m)
        real(r8), allocatable :: blk_zw(:)        ! (m)
#endif
# if defined SG_BBL || defined SSW_BBL
!
!-----------------------------------------------------------------------
!  Closure parameters associated with Styles and Glenn (1999) bottom
!  currents and waves boundary layer.
!-----------------------------------------------------------------------
!
!    sg_Cdmax      Upper limit on bottom darg coefficient.
!    sg_alpha      Free parameter indicating the constant stress
!                    region of the wave boundary layer.
!    sg_g          Acceleration of gravity (m/s2).
!    sg_kappa      Von Karman constant.
!    sg_mp         Nondimensional closure constant.
!    sg_n          Maximum number of iterations for bisection method.
!    sg_nu         Kinematic viscosity of seawater (m2/s).
!    sg_pi         Ratio of circumference to diameter.
!    sg_tol        Convergence criterion.
!    sg_ustarcdef  Default bottom stress (m/s).
!    sg_z100       Depth (m), 100 cm above bottom.
!    sg_z1p        Nondimensional closure constant.
!    sg_zrmin      Minimum allowed height (m) of current above bed.
!                    Otherwise, logarithmic interpolation is used.
!    sg_znotcdef   Default apparent hydraulic roughness (m).
!    sg_znotdef    Default hydraulic roughness (m).
!
        integer, parameter :: sg_n = 20
 
        real(dp), parameter :: sg_pi = pi
 
        real(dp) :: sg_cdmax = 0.01_dp         ! non-dimensional
        real(dp) :: sg_alpha = 1.0_dp          ! non-dimensional
        real(dp) :: sg_g = 9.81_dp             ! (m/s2)
        real(dp) :: sg_kappa = 0.41_dp         ! non-dimensional
        real(dp) :: sg_nu = 0.00000119_dp      ! (m2/s)
        real(dp) :: sg_tol = 0.0001_dp         ! non-dimensional
        real(dp) :: sg_ustarcdef = 0.01_dp     ! (m/s)
        real(dp) :: sg_z100 = 1.0_dp           ! (m)
        real(dp) :: sg_z1min = 0.20_dp         ! (m)
        real(dp) :: sg_z1p                     ! non-dimensional
        real(dp) :: sg_znotcdef = 0.01_dp      ! (m)
        real(dp) :: sg_znotdef                 ! (m)
 
        complex(c8) :: sg_mp
# endif
#if defined LMD_SKPP || defined SOLAR_SOURCE
!
!-----------------------------------------------------------------------
!  Water clarity parameters.
!-----------------------------------------------------------------------
!
!  The water type classification is based on Jerlov water type using
!  a double exponential function for light absorption:
!
!    Array
!    Index   WaterType   Examples
!    -----   ---------   --------
!
!      1         I       Open Pacific
!      2         IA      Eastern Mediterranean, Indian Ocean
!      3         IB      Western Mediterranean, Open Atlantic
!      4         II      Coastal waters, Azores
!      5         III     Coastal waters, North Sea
!      6         1       Skagerrak Strait
!      7         3       Baltic
!      8         5       Black Sea
!      9         7       Dark coastal water
!
!    lmd_mu1       Reciprocal of the absorption coefficient for solar
!                    wavelength band 1 as a function of the Jerlov
!                    water type.
!    lmd_mu2       Reciprocal of the absorption coefficient for solar
!                    wavelength band 2 as a function of the Jerlov
!                    water type.
!    lmd_r1        Fraction of total radiance for wavelength band 1 as
!                    a function of the Jerlov water type.
!
        real(r8), dimension(9) :: lmd_mu1 =                             &
     &            (/ 0.35_r8, 0.6_r8,  1.0_r8,  1.5_r8, 1.4_r8,         &
     &               0.42_r8, 0.37_r8, 0.33_r8, 0.00468592_r8 /)
 
        real(r8), dimension(9) :: lmd_mu2 =                             &
     &            (/ 23.0_r8,  20.0_r8,  17.0_r8,  14.0_r8, 7.9_r8,     &
     &                5.13_r8,  3.54_r8,  2.34_r8,  1.51_r8 /)
 
        real(r8), dimension(9) :: lmd_r1 =                              &
     &            (/ 0.58_r8, 0.62_r8, 0.67_r8, 0.77_r8, 0.78_r8,       &
     &               0.57_r8, 0.57_r8, 0.57_r8, 0.55_r8 /)
#endif
 
#ifdef LMD_MIXING
!
!-----------------------------------------------------------------------
!  Large et al. (1994) K-profile parameterization.
!-----------------------------------------------------------------------
!
!    lmd_Ri0       Critical gradient Richardson number below which
!                    turbulent mixing occurs.
!    lmd_Rrho0     Value of double-diffusive density ratio where
!                    mixing goes to zero in salt fingering.
!    lmd_bvfcon    Brunt-Vaisala frequency (1/s2) limit for convection.
!    lmd_fdd       Scaling factor for double diffusion of temperature
!                    in salt fingering case (lmd_fdd=0.7).
!    lmd_nu        Molecular viscosity (m2/s).
!    lmd_nu0c      Maximum interior convective viscosity and diffusivity
!                    due to shear instability.
!    lmd_nu0m      Maximum interior viscosity (m2/s) due shear
!                    instability.
!    lmd_nu0s      Maximum interior diffusivity (m2/s) due shear
!                    instability.
!    lmd_nuf       Scaling factor for double diffusion in salt
!                    fingering.
!    lmd_nuwm      Interior viscosity (m2/s) due to wave breaking.
!    lmd_nuws      Interior diffusivity (m2/s) due to wave breaking.
!    lmd_sdd1      Double diffusion constant for salinity in diffusive
!                    convection case (lmd_sdd1=0.15).
!    lmd_sdd2      Double diffusion constant for salinity in diffusive
!                    convection case (lmd_sdd2=1.85).
!    lmd_sdd3      Double diffusion constant for salinity in diffusive
!                    convection case (lmd_sdd3=0.85).
!    lmd_tdd1      Double diffusion constant for temperature
!                    in diffusive convection case (lmd_tdd1=0.909).
!    lmd_tdd2      Double diffusion constant for temperature in
!                    diffusive convection case (lmd_tdd2=4.6).
!    lmd_tdd3      Double diffusion constant for temperature in
!                    diffusive convection case (lmd_tdd3=0.54).
!
        real(r8) :: lmd_ri0 = 0.7_r8          ! non-dimensional
        real(r8) :: lmd_rrho0 = 1.9_r8        ! m2/s
        real(r8) :: lmd_bvfcon = -2.0e-5_r8   ! 1/s2
        real(r8) :: lmd_fdd = 0.7_r8          ! non-dimensional
        real(r8) :: lmd_nu = 1.5e-6_r8        ! m2/s
        real(r8) :: lmd_nu0c = 0.01_r8        ! m2/s
!!      real(r8) :: lmd_nu0c = 0.1_r8         ! m2/s
!!      real(r8) :: lmd_nu0c = 0.05_r8        ! m2/s
        real(r8) :: lmd_nu0m = 10.0e-4_r8     ! m2/s
        real(r8) :: lmd_nu0s = 10.0e-4_r8     ! m2/s
!!      real(r8) :: lmd_nu0m = 50.0E-4_r8     ! m2/s
!!      real(r8) :: lmd_nu0s = 50.0E-4_r8     ! m2/s
        real(r8) :: lmd_nuf = 10.0e-4_r8      ! m2/s
# ifdef DAMEE_4
        real(r8) :: lmd_nuwm = 1.0e-4_r8      ! m2/s
        real(r8) :: lmd_nuws = 1.0e-5_r8      ! m2/s
# else
        real(r8) :: lmd_nuwm = 1.0e-5_r8      ! m2/s
        real(r8) :: lmd_nuws = 1.0e-6_r8      ! m2/s
# endif
        real(r8) :: lmd_sdd1 = 0.15_r8        ! non-dimensional
        real(r8) :: lmd_sdd2 = 1.85_r8        ! non-dimensional
        real(r8) :: lmd_sdd3 = 0.85_r8        ! non-dimensional
        real(r8) :: lmd_tdd1 = 0.909_r8       ! non-dimensional
        real(r8) :: lmd_tdd2 = 4.6_r8         ! non-dimensional
        real(r8) :: lmd_tdd3 = 0.54_r8        ! non-dimensional
 
# if defined LMD_SKPP || defined LMD_BKPP
!
!-----------------------------------------------------------------------
!  Large et al. (1994) oceanic boundary layer parameters.
!-----------------------------------------------------------------------
!
!    lmd_Cg        Proportionality coefficient parameterizing nonlocal
!                    transport.
!    lmd_Cstar     Proportionality coefficient parameterizing nonlocal
!                    transport.
!    lmd_Cv        Ratio of interior Brunt-Vaisala frequency to that
!                    at entrainment depth "he".
!    lmd_Ric       Critical bulk Richardson number.
!    lmd_am        Coefficient of flux profile for momentum in their
!                    1/3 power law regimes.
!    lmd_as        Coefficient of flux profile for tracers in their
!                    1/3 power law regimes.
!    lmd_betaT     Ratio of entrainment flux to surface buoyancy flux.
!    lmd_cekman    Constant used in the computation of Ekman depth.
!    lmd_cmonob    Constant used in the computation of Monin-Obukhov
!                    depth.
!    lmd_cm        Coefficient of flux profile for momentum in their
!                    1/3 power law regimes.
!    lmd_cs        Coefficient of flux profile for tracers in their
!                    1/3 power law regimes.
!    lmd_epsilon   Non-dimensional extent of the surface layer.
!    lmd_zetam     Maximum stability parameter "zeta" value of the 1/3
!                    power law regime of flux profile for momentum.
!    lmd_zetas     Maximum stability parameter "zeta" value of the 1/3
!                    power law regime of flux profile for tracers.
!
        real(r8) :: lmd_cg
        real(r8) :: lmd_cstar = 10.0_r8
        real(r8) :: lmd_cv = 1.25_r8
!!      real(r8) :: lmd_Cv = 1.4_r8
!!      real(r8) :: lmd_Cv = 1.6_r8
!!      real(r8) :: lmd_Cv = 1.8_r8
!!      real(r8) :: lmd_Ric = 0.25_r8
        real(r8) :: lmd_ric = 0.3_r8
!!      real(r8) :: lmd_Ric = 0.5_r8
!!      real(r8) :: lmd_Ric = 0.75_r8
        real(r8) :: lmd_am = 1.257_r8
        real(r8) :: lmd_as = -28.86_r8
        real(r8) :: lmd_betat = -0.2_r8
        real(r8) :: lmd_cekman = 0.7_r8
        real(r8) :: lmd_cmonob = 1.0_r8
        real(r8) :: lmd_cm = 8.36_r8
        real(r8) :: lmd_cs = 98.96_r8
        real(r8) :: lmd_epsilon = 0.1_r8
        real(r8) :: lmd_zetam = -0.2_r8
        real(r8) :: lmd_zetas = -1.0_r8
# endif
#endif
!
!-----------------------------------------------------------------------
!  Generic Length Scale parameters.
!-----------------------------------------------------------------------
!
!    gls_Gh0
!    gls_Ghcri
!    gls_Ghmin
!    gls_Kmin      Minimum value of specific turbulent kinetic energy.
!    gls_Pmin      Minimum Value of dissipation.
!    gls_cmu0      Stability coefficient (non-dimensional).
!    gls_c1        Shear production coefficient (non-dimensional).
!    gls_c2        Dissipation coefficient (non-dimensional).
!    gls_c3m       Buoyancy production coefficient (minus).
!    gls_c3p       Buoyancy production coefficient (plus).
!    gls_E2
!    gls_m         Turbulent kinetic energy exponent (non-dimensional).
!    gls_n         Turbulent length scale exponent (non-dimensional).
!    gls_p         Stability exponent (non-dimensional).
!    gls_sigk      Constant Schmidt number (non-dimensional) for
!                    turbulent kinetic energy diffusivity.
!    gls_sigp      Constant Schmidt number (non-dimensional) for
!                    turbulent generic statistical field, "psi".
!
        real(r8), allocatable :: gls_m(:)
        real(r8), allocatable :: gls_n(:)
        real(r8), allocatable :: gls_p(:)
        real(r8), allocatable :: gls_sigk(:)
        real(r8), allocatable :: gls_sigp(:)
        real(r8), allocatable :: gls_cmu0(:)
        real(r8), allocatable :: gls_cmupr(:)
        real(r8), allocatable :: gls_c1(:)
        real(r8), allocatable :: gls_c2(:)
        real(r8), allocatable :: gls_c3m(:)
        real(r8), allocatable :: gls_c3p(:)
        real(r8), allocatable :: gls_kmin(:)
        real(r8), allocatable :: gls_pmin(:)
#ifdef GLS_MIXING
# if defined CANUTO_A || defined CANUTO_B
        real(r8) :: gls_s0
        real(r8) :: gls_s1
        real(r8) :: gls_s2
        real(r8) :: gls_s3
        real(r8) :: gls_s4
        real(r8) :: gls_s5
        real(r8) :: gls_s6
        real(r8) :: gls_b0
        real(r8) :: gls_b1
        real(r8) :: gls_b2
        real(r8) :: gls_b3
        real(r8) :: gls_b4
        real(r8) :: gls_b5
# endif
# ifdef CANUTO_A
        real(r8), parameter :: gls_gh0 = 0.0329_r8 ! 0.0329 GOTM, 0.0673 Burchard
        real(r8), parameter :: gls_ghcri = 0.03_r8
        real(r8), parameter :: gls_l1 = 0.107_r8
        real(r8), parameter :: gls_l2 = 0.0032_r8
        real(r8), parameter :: gls_l3 = 0.0864_r8
        real(r8), parameter :: gls_l4 = 0.12_r8
        real(r8), parameter :: gls_l5 = 11.9_r8
        real(r8), parameter :: gls_l6 = 0.4_r8
        real(r8), parameter :: gls_l7 = 0.0_r8
        real(r8), parameter :: gls_l8 = 0.48_r8
# elif defined CANUTO_B
        real(r8), parameter :: gls_gh0 = 0.0444_r8 ! 0.044 GOTM, 0.0673 Burchard
        real(r8), parameter :: gls_ghcri = 0.0414_r8
        real(r8), parameter :: gls_l1 = 0.127_r8
        real(r8), parameter :: gls_l2 = 0.00336_r8
        real(r8), parameter :: gls_l3 = 0.0906_r8
        real(r8), parameter :: gls_l4 = 0.101_r8
        real(r8), parameter :: gls_l5 = 11.2_r8
        real(r8), parameter :: gls_l6 = 0.4_r8
        real(r8), parameter :: gls_l7 = 0.0_r8
        real(r8), parameter :: gls_l8 = 0.318_r8
# else
        real(r8), parameter :: gls_gh0 = 0.028_r8
        real(r8), parameter :: gls_ghcri = 0.02_r8
# endif
        real(r8), parameter :: gls_ghmin = -0.28_r8
        real(r8), parameter :: gls_e2 = 1.33_r8
#endif
!
! Constants used in the various formulation of surface flux boundary
! conditions for the GLS vertical turbulence closure in terms of
! Charnock surface roughness (CHARNOK_ALPHA), roughness from wave
! amplitude (zos_hsig_alpha), wave dissipation (SZ_ALPHA), and
! Craig and Banner wave breaking (CRGBAN_CW).
!
        real(r8), allocatable :: charnok_alpha(:)
        real(r8), allocatable :: zos_hsig_alpha(:)
        real(r8), allocatable :: sz_alpha(:)
        real(r8), allocatable :: crgban_cw(:)
!
! Waves Effect on Currents dissipation energy partion scale:
!
!   [0.0] All wave dissipation goes to breaking and none to roller
!   [1.0] All wave dissipation goes to roller and none to breaking
!
        real(r8), allocatable :: wec_alpha(:)
 
#if defined MY25_MIXING || defined GLS_MIXING
!
!-----------------------------------------------------------------------
!  Mellor-Yamada (1982) Level 2.5 vertical mixing variables.
!-----------------------------------------------------------------------
!
!    my_A1         Turbulent closure A1 constant.
!    my_A2         Turbulent closure A2 constant.
!    my_B1         Turbulent closure B1 constant.
!    my_B1p2o3     B1**(2/3).
!    my_B1pm1o3    B1**(-1/3).
!    my_B2         Turbulent closure B2 constant.
!    my_C1         Turbulent closure C1 constant.
!    my_C2         Turbulent closure C2 constant.
!    my_C3         Turbulent closure C3 constant.
!    my_E1         Turbulent closure E1 constant.
!    my_E1o2       0.5*E1
!    my_E2         Turbulent closure E2 constant.
!    my_Gh0        Lower bound on Galperin et al. stability function.
!    my_Sh1        Tracers stability function constant factor.
!    my_Sh2        Tracers stability function constant factor.
!    my_Sm1        Momentum stability function constant factor.
!    my_Sm2        Momentum stability function constant factor.
!    my_Sm3        Momentum stability function constant factor.
!    my_Sm4        Momentum stability function constant factor.
!    my_Sq         Scale for vertical mixing of turbulent energy.
!    my_dtfac      Asselin time filter coefficient.
!    my_lmax       Upper bound on the turbulent length scale.
!    my_qmin       Lower bound on turbulent energy "tke" and "gls".
!
        real(r8), parameter :: my_a1 = 0.92_r8
        real(r8), parameter :: my_a2 = 0.74_r8
        real(r8), parameter :: my_b1 = 16.6_r8
        real(r8), parameter :: my_b2 = 10.1_r8
        real(r8), parameter :: my_c1 = 0.08_r8
        real(r8), parameter :: my_c2 = 0.7_r8
        real(r8), parameter :: my_c3 = 0.2_r8
        real(r8), parameter :: my_e1 = 1.8_r8
        real(r8), parameter :: my_e2 = 1.33_r8
        real(r8), parameter :: my_gh0 = 0.0233_r8
        real(r8), parameter :: my_sq = 0.2_r8
        real(r8), parameter :: my_dtfac = 0.05_r8
        real(r8), parameter :: my_lmax = 0.53_r8
        real(r8), parameter :: my_qmin = 1.0e-8_r8
 
        real(r8) :: my_b1p2o3
        real(r8) :: my_b1pm1o3
        real(r8) :: my_e1o2
        real(r8) :: my_sh1
        real(r8) :: my_sh2
        real(r8) :: my_sm1
        real(r8) :: my_sm2
        real(r8) :: my_sm3
        real(r8) :: my_sm4
#endif
#ifdef BVF_MIXING
!
!-----------------------------------------------------------------------
!  Brunt-Vaisala frequency based vertical mixing variables.
!-----------------------------------------------------------------------
!
!    bvf_numax     Upper bound vertical diffusion (m2/s).
!    bvf_numin     Lower bound vertical diffusion (m2/s).
!    bvf_nu0       Proportionality constant (m2/s2).
!    bvf_nu0c      Convective diffusion (m2/s) in static unstable
!                    regime.
!
        real(r8) :: bvf_numax = 4.0e-4_r8      ! m2/s
        real(r8) :: bvf_numin = 3.0e-5_r8      ! m2/s
        real(r8) :: bvf_nu0 = 1.0e-7_r8        ! m2/s2
        real(r8) :: bvf_nu0c = 1.0_r8          ! m2/s
#endif
!
!-----------------------------------------------------------------------
!  Tangent linear and adjoint model parameters.
!-----------------------------------------------------------------------
!
!  Tangent linear and adjoint model control switches.
!
        logical :: tlmodel = .false.
        logical :: admodel = .false.
!
      CONTAINS
!
      SUBROUTINE allocate_scalars
!
!=======================================================================
!                                                                      !
!  This routine allocates structure and  several variables in module   !
!  that depend on the number of nested grids.                          !
!                                                                      !
!=======================================================================
!
      USE mod_param
!
!  Local variable declarations.
!
      integer :: ng
 
      real(r8), parameter :: inival = 0.0_r8
!
!-----------------------------------------------------------------------
!  Allocate and initialize variables in module structure.
!-----------------------------------------------------------------------
!
      IF (.not.allocated(scalars)) THEN
        allocate ( scalars(ngrids) )
 
        DO ng=1,ngrids
 
#ifdef SOLVE3D
          allocate ( scalars(ng) % Fstate(9+2*mt) )
          dmem(ng)=dmem(ng)+real(9+2*mt,r8)
          scalars(ng) % Fstate = .false.
 
          allocate ( scalars(ng) % Lstate(8+2*mt) )
          dmem(ng)=dmem(ng)+real(9+2*mt,r8)
          scalars(ng) % Lstate = .false.
#else
          allocate ( scalars(ng) % Fstate(5) )
          dmem(ng)=dmem(ng)+5.0_r8
          scalars(ng) % Fstate = .false.
 
          allocate ( scalars(ng) % Lstate(3) )
          dmem(ng)=dmem(ng)+3.0_r8
          scalars(ng) % Lstate = .false.
#endif
 
#if defined READ_WATER && defined MASKING && defined DISTRIBUTE
          allocate ( scalars(ng) % IJwater((lm(ng)+2)*(mm(ng)+2),4) )
          dmem(ng)=dmem(ng)+4.0_r8*real((lm(ng)+2)*(mm(ng)+2),r8)
          scalars(ng) % IJwater = 0
#endif
 
          allocate ( scalars(ng) % Cs_r(n(ng)) )
          dmem(ng)=dmem(ng)+real(n(ng),r8)
          scalars(ng) % Cs_r = inival
 
          allocate ( scalars(ng) % Cs_w(0:n(ng)) )
          dmem(ng)=dmem(ng)+real(n(ng)+1,r8)
          scalars(ng) % Cs_w = inival
 
          allocate ( scalars(ng) % sc_r(n(ng)) )
          dmem(ng)=dmem(ng)+real(n(ng),r8)
          scalars(ng) % sc_r = inival
 
          allocate ( scalars(ng) % sc_w(0:n(ng)) )
          dmem(ng)=dmem(ng)+real(n(ng)+1,r8)
          scalars(ng) % sc_w = inival
 
        END DO
      END IF
!
!  Allocate variables that require special treatment in shared-memory.
!  These variables are private for each thread to avoid collisions.
!
!$OMP PARALLEL
      IF (.not.allocated(predictor_2d_step)) THEN
        allocate ( predictor_2d_step(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(indx1)) THEN
        allocate ( indx1(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(iic)) THEN
        allocate ( iic(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(iif)) THEN
        allocate ( iif(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(next_kstp)) THEN
        allocate ( next_kstp(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(frcrec)) THEN
        allocate ( frcrec(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(sorec)) THEN
        allocate ( sorec(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
!$OMP END PARALLEL
 
#ifdef JEDI
      IF (.not.allocated(jic)) THEN
        allocate ( jic(ngrids) )
        jic=0
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
!
!-----------------------------------------------------------------------
!  Allocate variables.
!-----------------------------------------------------------------------
!
#if defined FOUR_DVAR || defined VERIFICATION
      mstatevar=6+mt
# if defined GLS_MIXING || defined MY25_MIXING
      mstatevar=mstatevar+1
# endif
# ifdef WEC
      mstatevar=mstatevar+2
#  if defined SOLVE3D
      mstatevar=mstatevar+2
#  endif
# endif
# ifdef ADJUST_WSTRESS
      mstatevar=mstatevar+2
# endif
# ifdef ADJUST_STFLUX
      mstatevar=mstatevar+mt
# endif
!
#endif
 
#ifdef T_PASSIVE
      IF (.not.associated(inert)) THEN
        allocate ( inert(npt) )
        dmem(1)=dmem(1)+real(npt,r8)
      END IF
#endif
      IF (.not.allocated(perfectrst)) THEN
        allocate ( perfectrst(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ndtfast)) THEN
        allocate ( ndtfast(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nfast)) THEN
        allocate ( nfast(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(dt)) THEN
        allocate ( dt(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(dtfast)) THEN
        allocate ( dtfast(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(timeend)) THEN
        allocate ( timeend(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(avgtime)) THEN
        allocate ( avgtime(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(diatime)) THEN
        allocate ( diatime(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(imptime)) THEN
        allocate ( imptime(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(initime)) THEN
        allocate ( initime(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(initimes)) THEN
        allocate ( initimes(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(obstime)) THEN
        allocate ( obstime(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(frctime)) THEN
        allocate ( frctime(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ntimes)) THEN
        allocate ( ntimes(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(first_time)) THEN
        allocate ( first_time(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ntfirst)) THEN
        allocate ( ntfirst(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ntstart)) THEN
        allocate ( ntstart(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ntend)) THEN
        allocate ( ntend(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
!$OMP PARALLEL
      IF (.not.allocated(synchro_flag)) THEN
        allocate ( synchro_flag(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(step_counter)) THEN
        allocate ( step_counter(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(tdays)) THEN
        allocate ( tdays(ngrids) )
        tdays=0.0_dp
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(time)) THEN
        allocate ( time(ngrids) )
        time=0.0_dp
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(time_code)) THEN
        allocate ( time_code(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
!$OMP END PARALLEL
 
#ifdef JEDI
      IF (.not.allocated(time4jedi)) THEN
        allocate ( time4jedi(ngrids) )
        time4jedi=0.0_dp
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
 
      IF (.not.allocated(nrecfrc)) THEN
        allocate ( nrecfrc(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(setgridconfig)) THEN
        allocate ( setgridconfig(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(processinputdata)) THEN
        allocate ( processinputdata(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(nudgingcoeff)) THEN
        allocate ( nudgingcoeff(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(obcdata)) THEN
        allocate ( obcdata(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(lbiology)) THEN
        allocate ( lbiology(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lfloats)) THEN
        allocate ( lfloats(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lsediment)) THEN
        allocate ( lsediment(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lstations)) THEN
        allocate ( lstations(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(compositegrid)) THEN
        allocate ( compositegrid(4,ngrids) )
        dmem(1)=dmem(1)+4.0_r8*real(ngrids,r8)
      END IF
      IF (.not.allocated(refinedgrid)) THEN
        allocate ( refinedgrid(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(refinescale)) THEN
        allocate ( refinescale(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(getdonordata)) THEN
        allocate ( getdonordata(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ewperiodic)) THEN
        allocate ( ewperiodic(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nsperiodic)) THEN
        allocate ( nsperiodic(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(volcons)) THEN
        allocate ( volcons(4,ngrids) )
        dmem(1)=dmem(1)+4.0_r8*real(ngrids,r8)
      END IF
#if defined ADJOINT || defined TANGENT || defined TL_IOMS
      IF (.not.allocated(ad_volcons)) THEN
        allocate ( ad_volcons(4,ngrids) )
        dmem(1)=dmem(1)+4.0_r8*real(ngrids,r8)
      END IF
      IF (.not.allocated(tl_volcons)) THEN
        allocate ( tl_volcons(4,ngrids) )
        dmem(1)=dmem(1)+4.0_r8*real(ngrids,r8)
      END IF
#endif
 
      IF (.not.allocated(lsponge)) THEN
        allocate ( lsponge(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(luvsponge)) THEN
        allocate ( luvsponge(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ltracersponge)) THEN
        allocate ( ltracersponge(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
 
      IF (.not.allocated(clm_file)) THEN
        allocate ( clm_file(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lclimatology)) THEN
        allocate ( lclimatology(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lsshclm)) THEN
        allocate ( lsshclm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lm2clm)) THEN
        allocate ( lm2clm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lm3clm)) THEN
        allocate ( lm3clm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ltracerclm)) THEN
        allocate ( ltracerclm(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
 
      IF (.not.allocated(lnudging)) THEN
        allocate ( lnudging(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lnudgem2clm)) THEN
        allocate ( lnudgem2clm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lnudgem3clm)) THEN
        allocate ( lnudgem3clm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lnudgetclm)) THEN
        allocate ( lnudgetclm(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
 
      IF (.not.allocated(luvsrc)) THEN
        allocate ( luvsrc(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwsrc)) THEN
        allocate ( lwsrc(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ltracersrc)) THEN
        allocate ( ltracersrc(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
 
      IF (.not.allocated(maxspeed)) THEN
        allocate ( maxspeed(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(maxrho)) THEN
        allocate ( maxrho(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(levsfrc)) THEN
        allocate ( levsfrc(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(levbfrc)) THEN
        allocate ( levbfrc(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(vtransform)) THEN
        allocate ( vtransform(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(vstretching)) THEN
        allocate ( vstretching(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(tcline)) THEN
        allocate ( tcline(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(hc)) THEN
        allocate ( hc(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(theta_s)) THEN
        allocate ( theta_s(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(theta_b)) THEN
        allocate ( theta_b(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(hmin)) THEN
        allocate ( hmin(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(hmax)) THEN
        allocate ( hmax(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(xl)) THEN
        allocate ( xl(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(el)) THEN
        allocate ( el(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(lonmin)) THEN
        allocate ( lonmin(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lonmax)) THEN
        allocate ( lonmax(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(latmin)) THEN
        allocate ( latmin(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(latmax)) THEN
        allocate ( latmax(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(idigits)) THEN
        allocate ( idigits(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(jdigits)) THEN
        allocate ( jdigits(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#ifdef SOLVE3D
      IF (.not.allocated(kdigits)) THEN
        allocate ( kdigits(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
      IF (.not.allocated(totvolume)) THEN
        allocate ( totvolume(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(minvolume)) THEN
        allocate ( minvolume(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(maxvolume)) THEN
        allocate ( maxvolume(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(dxmin)) THEN
        allocate ( dxmin(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(dxmax)) THEN
        allocate ( dxmax(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(dymin)) THEN
        allocate ( dymin(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(dymax)) THEN
        allocate ( dymax(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
#ifdef MASKING
      IF (.not.allocated(dxminw)) THEN
        allocate ( dxminw(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(dxmaxw)) THEN
        allocate ( dxmaxw(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(dyminw)) THEN
        allocate ( dyminw(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(dymaxw)) THEN
        allocate ( dymaxw(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
 
 
#ifdef SOLVE3D
      IF (.not.allocated(dzmin)) THEN
        allocate ( dzmin(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(dzmax)) THEN
        allocate ( dzmax(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
# ifdef MASKING
      IF (.not.allocated(dzminw)) THEN
        allocate ( dzminw(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(dzmaxw)) THEN
        allocate ( dzmaxw(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
# endif
#endif
 
      IF (.not.allocated(grdmax)) THEN
        allocate ( grdmax(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
#ifdef DIFF_3DCOEF
      IF (.not.allocated(diffmin)) THEN
        allocate ( diffmin(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(diffmax)) THEN
        allocate ( diffmax(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
 
      IF (.not.allocated(cg_min)) THEN
        allocate ( cg_min(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(cg_max)) THEN
        allocate ( cg_max(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(cg_cor)) THEN
        allocate ( cg_cor(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
#ifdef VISC_3DCOEF
      IF (.not.allocated(viscmin)) THEN
        allocate ( viscmin(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(viscmax)) THEN
        allocate ( viscmax(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
 
      IF (.not.allocated(r0)) THEN
        allocate ( r0(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(tcoef)) THEN
        allocate ( tcoef(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(scoef)) THEN
        allocate ( scoef(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(t0)) THEN
        allocate ( t0(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(s0)) THEN
        allocate ( s0(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(gamma2)) THEN
        allocate ( gamma2(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(lmd_jwt)) THEN
        allocate ( lmd_jwt(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(rx0)) THEN
        allocate ( rx0(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(rx1)) THEN
        allocate ( rx1(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(rdrg)) THEN
        allocate ( rdrg(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(rdrg2)) THEN
        allocate ( rdrg2(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(zos)) THEN
        allocate ( zos(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(zob)) THEN
        allocate ( zob(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(dcrit)) THEN
        allocate ( dcrit(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
#ifdef PROPAGATOR
      IF (.not.allocated(nconv)) THEN
        allocate ( nconv(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
 
      IF (.not.allocated(weight)) THEN
        allocate ( weight(2,0:256,ngrids) )
        dmem(1)=dmem(1)+514.0_r8*real(ngrids,r8)
      END IF
 
      IF (.not.allocated(akk_bak)) THEN
        allocate ( akk_bak(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(akp_bak)) THEN
        allocate ( akp_bak(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(akv_bak)) THEN
        allocate ( akv_bak(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(akv_limit)) THEN
        allocate ( akv_limit(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ad_visc2)) THEN
        allocate ( ad_visc2(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nl_visc2)) THEN
        allocate ( nl_visc2(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(tl_visc2)) THEN
        allocate ( tl_visc2(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(visc2)) THEN
        allocate ( visc2(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ad_visc4)) THEN
        allocate ( ad_visc4(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nl_visc4)) THEN
        allocate ( nl_visc4(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(tl_visc4)) THEN
        allocate ( tl_visc4(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(visc4)) THEN
        allocate ( visc4(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(tkenu2)) THEN
        allocate ( tkenu2(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(tkenu4)) THEN
        allocate ( tkenu4(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(akt_bak)) THEN
        allocate ( akt_bak(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
      IF (.not.allocated(akt_limit)) THEN
        allocate ( akt_limit(nat,ngrids) )
        dmem(1)=dmem(1)+real(nat*ngrids,r8)
      END IF
      IF (.not.allocated(kdiff)) THEN
        allocate ( kdiff(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
 
      IF (.not.allocated(ad_tnu2)) THEN
        allocate ( ad_tnu2(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
      IF (.not.allocated(nl_tnu2)) THEN
        allocate ( nl_tnu2(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
      IF (.not.allocated(tl_tnu2)) THEN
        allocate ( tl_tnu2(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
      IF (.not.allocated(tnu2)) THEN
        allocate ( tnu2(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
 
      IF (.not.allocated(ad_tnu4)) THEN
        allocate ( ad_tnu4(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
      IF (.not.allocated(nl_tnu4)) THEN
        allocate ( nl_tnu4(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
      IF (.not.allocated(tl_tnu4)) THEN
        allocate ( tl_tnu4(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
      IF (.not.allocated(tnu4)) THEN
        allocate ( tnu4(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
 
      IF (.not.allocated(tl_m2diff)) THEN
        allocate ( tl_m2diff(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(tl_m3diff)) THEN
        allocate ( tl_m3diff(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(tl_tdiff)) THEN
        allocate ( tl_tdiff(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
 
      IF (.not.allocated(ad_akv_fac)) THEN
        allocate ( ad_akv_fac(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(tl_akv_fac)) THEN
        allocate ( tl_akv_fac(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ad_akt_fac)) THEN
        allocate ( ad_akt_fac(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
      IF (.not.allocated(tl_akt_fac)) THEN
        allocate ( tl_akt_fac(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
 
      IF (.not.allocated(ldefadj)) THEN
        allocate ( ldefadj(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefavg)) THEN
        allocate ( ldefavg(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefdai)) THEN
        allocate ( ldefdai(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefdia)) THEN
        allocate ( ldefdia(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldeferr)) THEN
        allocate ( ldeferr(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefflt)) THEN
        allocate ( ldefflt(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefhis)) THEN
        allocate ( ldefhis(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefhss)) THEN
        allocate ( ldefhss(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefini)) THEN
        allocate ( ldefini(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefirp)) THEN
        allocate ( ldefirp(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefitl)) THEN
        allocate ( ldefitl(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldeflcz)) THEN
        allocate ( ldeflcz(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldeflze)) THEN
        allocate ( ldeflze(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefmod)) THEN
        allocate ( ldefmod(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefqck)) THEN
        allocate ( ldefqck(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefrst)) THEN
        allocate ( ldefrst(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefxtr)) THEN
        allocate ( ldefxtr(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#ifdef SP4DVAR
      IF (.not.allocated(ldefsca)) THEN
        allocate ( ldefsca(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefsct)) THEN
        allocate ( ldefsct(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefspa)) THEN
        allocate ( ldefspa(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldefspt)) THEN
        allocate ( ldefspt(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
      IF (.not.allocated(ldefsta)) THEN
        allocate ( ldefsta(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldeftide)) THEN
        allocate ( ldeftide(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldeftlm)) THEN
        allocate ( ldeftlm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ldeftlf)) THEN
        allocate ( ldeftlf(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(lreadadm)) THEN
        allocate ( lreadadm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lreadblk)) THEN
        allocate ( lreadblk(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lreadfrc)) THEN
        allocate ( lreadfrc(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lreadfwd)) THEN
        allocate ( lreadfwd(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lreadqck)) THEN
        allocate ( lreadqck(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lreadstd)) THEN
        allocate ( lreadstd(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lreadtlm)) THEN
        allocate ( lreadtlm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(lwrtadj)) THEN
        allocate ( lwrtadj(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrtavg)) THEN
        allocate ( lwrtavg(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrtdia)) THEN
        allocate ( lwrtdia(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrthis)) THEN
        allocate ( lwrthis(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrtper)) THEN
        allocate ( lwrtper(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrtqck)) THEN
        allocate ( lwrtqck(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrtrst)) THEN
        allocate ( lwrtrst(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrttlf)) THEN
        allocate ( lwrttlf(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrttlm)) THEN
        allocate ( lwrttlm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrtxtr)) THEN
        allocate ( lwrtxtr(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ldefnrm)) THEN
        allocate ( ldefnrm(4,ngrids) )
        dmem(1)=dmem(1)+4.0_r8*real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrtnrm)) THEN
        allocate ( lwrtnrm(4,ngrids) )
        dmem(1)=dmem(1)+4.0_r8*real(ngrids,r8)
      END IF
 
#ifdef STD_MODEL
      IF (.not.allocated(ldefstd)) THEN
        allocate ( ldefstd(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrtstd)) THEN
        allocate ( lwrtstd(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
 
#ifdef AD_OUTPUT_STATE
      IF (.not.allocated(lwrtstate3d)) THEN
        allocate ( lwrtstate3d(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
#if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
      IF (.not.allocated(soinitial)) THEN
        allocate ( soinitial(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
 
      IF (.not.allocated(lwrtstate2d)) THEN
        allocate ( lwrtstate2d(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrttime)) THEN
        allocate ( lwrttime(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lwrtcost)) THEN
        allocate ( lwrtcost(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ladjusted)) THEN
        allocate ( ladjusted(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(lprocessobc)) THEN
        allocate ( lprocessobc(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lprocesstides)) THEN
        allocate ( lprocesstides(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(lwrtinfo)) THEN
        allocate ( lwrtinfo(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ldefout)) THEN
        allocate ( ldefout(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ndefadj)) THEN
        allocate ( ndefadj(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ndefavg)) THEN
        allocate ( ndefavg(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ndefdia)) THEN
        allocate ( ndefdia(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ndefhis)) THEN
        allocate ( ndefhis(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ndefqck)) THEN
        allocate ( ndefqck(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ndeftlm)) THEN
        allocate ( ndeftlm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ndefxtr)) THEN
        allocate ( ndefxtr(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ntsavg)) THEN
        allocate ( ntsavg(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ntsdia)) THEN
        allocate ( ntsdia(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(nadj)) THEN
        allocate ( nadj(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(navg)) THEN
        allocate ( navg(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ndia)) THEN
        allocate ( ndia(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nflt)) THEN
        allocate ( nflt(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nhis)) THEN
        allocate ( nhis(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nqck)) THEN
        allocate ( nqck(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nrst)) THEN
        allocate ( nrst(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nsta)) THEN
        allocate ( nsta(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ntlm)) THEN
        allocate ( ntlm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nxtr)) THEN
        allocate ( nxtr(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#ifdef SP4DVAR
      IF (.not.allocated(nsca)) THEN
        allocate ( nsca(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nsct)) THEN
        allocate ( nsct(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nspa)) THEN
        allocate ( nspa(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nspt)) THEN
        allocate ( nspt(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(crec)) THEN
        allocate ( crec(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
 
      IF (.not.allocated(extractflag)) THEN
        allocate ( extractflag(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(ninfo)) THEN
        allocate ( ninfo(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(nobc)) THEN
        allocate ( nobc(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nbrec)) THEN
        allocate ( nbrec(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(obccount)) THEN
        allocate ( obccount(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
#ifdef ADJUST_BOUNDARY
      IF (.not.allocated(lobc)) THEN
        allocate ( lobc(4,mstatevar,ngrids) )
        dmem(1)=dmem(1)+4.0_r8*real(mstatevar*ngrids,r8)
      END IF
#endif
 
      IF (.not.allocated(nsff)) THEN
        allocate ( nsff(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nfrec)) THEN
        allocate ( nfrec(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(sfcount)) THEN
        allocate ( sfcount(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
#ifdef ADJUST_STFLUX
      IF (.not.allocated(lstflux)) THEN
        allocate ( lstflux(mt,ngrids) )
        dmem(1)=dmem(1)+real(mt*ngrids,r8)
      END IF
#endif
 
      IF (.not.allocated(nrrec)) THEN
        allocate ( nrrec(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(lastrec)) THEN
        allocate ( lastrec(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(lcycleadj)) THEN
        allocate ( lcycleadj(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lcyclerst)) THEN
        allocate ( lcyclerst(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(lcycletlm)) THEN
        allocate ( lcycletlm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
#if defined AVERAGES  && defined AVERAGES_DETIDE && \
   (defined ssh_tides || defined uv_tides)
      IF (.not.allocated(hcount)) THEN
        allocate ( hcount(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
 
      IF (.not.allocated(kstrs)) THEN
        allocate ( kstrs(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(kends)) THEN
        allocate ( kends(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(dstrs)) THEN
        allocate ( dstrs(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(dends)) THEN
        allocate ( dends(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(trnorm)) THEN
        allocate ( trnorm(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(so_decay)) THEN
        allocate ( so_decay(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#ifdef SOLVE3D
      IF (.not.allocated(so_sdev)) THEN
        allocate ( so_sdev(7+2*mt,ngrids) )
        dmem(1)=dmem(1)+real((7+2*mt)*ngrids,r8)
      END IF
#else
      IF (.not.allocated(so_sdev)) THEN
        allocate ( so_sdev(5,ngrids) )
        dmem(1)=dmem(1)+5.0_r8*real(ngrids,r8)
      END IF
#endif
 
#if defined FOUR_DVAR || defined VERIFICATION
      IF (.not.allocated(cnorm)) THEN
        allocate ( cnorm(2,mstatevar) )
        dmem(1)=dmem(1)+2.0_r8*real(mstatevar,r8)
      END IF
      IF (.not.allocated(cnormb)) THEN
        allocate ( cnormb(mstatevar,4) )
        dmem(1)=dmem(1)+4.0_r8*real(mstatevar,r8)
      END IF
 
      IF (.not.allocated(sporadicimpulse)) THEN
        allocate ( sporadicimpulse(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(frequentimpulse)) THEN
        allocate ( frequentimpulse(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(iauswitch)) THEN
        allocate ( iauswitch(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(timeiau)) THEN
        allocate ( timeiau(ngrids) )
        timeiau=0.0_dp
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(dtdz_min)) THEN
        allocate ( dtdz_min(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(ml_depth)) THEN
        allocate ( ml_depth(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(lnm_depth)) THEN
        allocate ( lnm_depth(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(balance)) THEN
        allocate ( balance(mstatevar) )
        dmem(1)=dmem(1)+real(mstatevar,r8)
      END IF
 
      IF (.not.allocated(hdecay)) THEN
        allocate ( hdecay(2,mstatevar,ngrids) )
        dmem(1)=dmem(1)+2.0_r8*real(mstatevar*ngrids,r8)
      END IF
      IF (.not.allocated(vdecay)) THEN
        allocate ( vdecay(2,mstatevar,ngrids) )
        dmem(1)=dmem(1)+2.0_r8*real(mstatevar*ngrids,r8)
      END IF
 
      IF (.not.allocated(tdecay)) THEN
        allocate ( tdecay(mstatevar,ngrids) )
        dmem(1)=dmem(1)+real(mstatevar*ngrids,r8)
      END IF
 
      IF (.not.allocated(hdecayb)) THEN
        allocate ( hdecayb(mstatevar,4,ngrids) )
        dmem(1)=dmem(1)+4.0_r8*real(mstatevar*ngrids,r8)
      END IF
      IF (.not.allocated(vdecayb)) THEN
        allocate ( vdecayb(mstatevar,4,ngrids) )
        dmem(1)=dmem(1)+4.0_r8*real(mstatevar*ngrids,r8)
      END IF
 
      IF (.not.allocated(bgqc_type)) THEN
        allocate ( bgqc_type(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(nprovenance)) THEN
        allocate ( nprovenance(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(s_bgqc)) THEN
        allocate ( s_bgqc(mstatevar,ngrids) )
        dmem(1)=dmem(1)+real(mstatevar*ngrids,r8)
      END IF
 
# ifdef STD_MODEL
      IF (.not.allocated(sigma_max)) THEN
        allocate ( sigma_max(mstatevar,ngrids) )
        dmem(1)=dmem(1)+real(mstatevar*ngrids,r8)
      END IF
      IF (.not.allocated(sigma_ml)) THEN
        allocate ( sigma_ml(mstatevar,ngrids) )
        dmem(1)=dmem(1)+real(mstatevar*ngrids,r8)
      END IF
      IF (.not.allocated(sigma_do)) THEN
        allocate ( sigma_do(mstatevar,ngrids) )
        dmem(1)=dmem(1)+real(mstatevar*ngrids,r8)
      END IF
      IF (.not.allocated(sigma_dz)) THEN
        allocate ( sigma_dz(mstatevar,ngrids) )
        dmem(1)=dmem(1)+real(mstatevar*ngrids,r8)
      END IF
 
      IF (.not.allocated(mld_uniform)) THEN
        allocate ( mld_uniform(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
# endif
#endif
 
      IF (.not.allocated(obcfac)) THEN
        allocate ( obcfac(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(fsobc_in)) THEN
        allocate ( fsobc_in(ngrids,4) )
        dmem(1)=dmem(1)+4.0_r8*real(ngrids,r8)
      END IF
      IF (.not.allocated(fsobc_out)) THEN
        allocate ( fsobc_out(ngrids,4) )
        dmem(1)=dmem(1)+4.0_r8*real(ngrids,r8)
      END IF
      IF (.not.allocated(m2obc_in)) THEN
        allocate ( m2obc_in(ngrids,4) )
        dmem(1)=dmem(1)+4.0_r8*real(ngrids,r8)
      END IF
      IF (.not.allocated(m2obc_out)) THEN
        allocate ( m2obc_out(ngrids,4) )
        dmem(1)=dmem(1)+4.0_r8*real(ngrids,r8)
      END IF
 
#ifdef SOLVE3D
      IF (.not.allocated(m3obc_in)) THEN
        allocate ( m3obc_in(ngrids,4) )
        dmem(1)=dmem(1)+4.0_r8*real(ngrids,r8)
      END IF
      IF (.not.allocated(m3obc_out)) THEN
        allocate ( m3obc_out(ngrids,4) )
        dmem(1)=dmem(1)+4.0_r8*real(ngrids,r8)
      END IF
      IF (.not.allocated(tobc_in)) THEN
        allocate ( tobc_in(mt,ngrids,4) )
        dmem(1)=dmem(1)+4.0_r8*real(mt*ngrids,r8)
      END IF
      IF (.not.allocated(tobc_out)) THEN
        allocate ( tobc_out(mt,ngrids,4) )
        dmem(1)=dmem(1)+4.0_r8*real(mt*ngrids,r8)
      END IF
#endif
 
      IF (.not.allocated(znudg)) THEN
        allocate ( znudg(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(m2nudg)) THEN
        allocate ( m2nudg(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(m3nudg)) THEN
        allocate ( m3nudg(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(tnudg)) THEN
        allocate ( tnudg(mt,ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
#ifdef BULK_FLUXES
      IF (.not.allocated(blk_zq)) THEN
        allocate ( blk_zq(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(blk_zt)) THEN
        allocate ( blk_zt(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(blk_zw)) THEN
        allocate ( blk_zw(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
#endif
 
      IF (.not.allocated(gls_m)) THEN
        allocate ( gls_m(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(gls_n)) THEN
        allocate ( gls_n(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(gls_p)) THEN
        allocate ( gls_p(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(gls_sigk)) THEN
        allocate ( gls_sigk(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(gls_sigp)) THEN
        allocate ( gls_sigp(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(gls_cmu0)) THEN
        allocate ( gls_cmu0(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(gls_cmupr)) THEN
        allocate ( gls_cmupr(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(gls_c1)) THEN
        allocate ( gls_c1(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(gls_c2)) THEN
        allocate ( gls_c2(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(gls_c3m)) THEN
        allocate ( gls_c3m(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(gls_c3p)) THEN
        allocate ( gls_c3p(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(gls_kmin)) THEN
        allocate ( gls_kmin(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(gls_pmin)) THEN
        allocate ( gls_pmin(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
 
      IF (.not.allocated(charnok_alpha)) THEN
        allocate ( charnok_alpha(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(zos_hsig_alpha)) THEN
        allocate ( zos_hsig_alpha(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(sz_alpha)) THEN
        allocate ( sz_alpha(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(crgban_cw)) THEN
        allocate ( crgban_cw(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
      IF (.not.allocated(wec_alpha)) THEN
        allocate ( wec_alpha(ngrids) )
        dmem(1)=dmem(1)+real(ngrids,r8)
      END IF
!
      RETURN
      END SUBROUTINE allocate_scalars
!
      SUBROUTINE deallocate_scalars (ng)
!
!=======================================================================
!                                                                      !
!  This routine deallocates structures variables and module variables. !
!  Notice that "destroy" cannot be use to deallocate pointer variables !
!  because of cyclic dependencies.                                     !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      integer :: ng
 
#ifdef SUBOBJECT_DEALLOCATION
!
!-----------------------------------------------------------------------
!  Deallocate each variable in the derived-type T_SCALARS structure
!  separately.
!-----------------------------------------------------------------------
!
# ifdef SOLVE3D
      IF (associated(scalars(ng)%Fstate)) THEN
        deallocate ( scalars(ng)%Fstate )
      END IF
 
      IF (associated(scalars(ng)%Lstate)) THEN
        deallocate ( scalars(ng)%Lstate )
      END IF
# else
      IF (associated(scalars(ng)%Fstate)) THEN
        deallocate ( scalars(ng)%Fstate )
      END IF
 
      IF (associated(scalars(ng)%Lstate)) THEN
        deallocate ( scalars(ng)%Lstate )
      END IF
# endif
 
# if defined READ_WATER && defined MASKING && defined DISTRIBUTE
      IF (associated(scalars(ng)%IJwater)) THEN
        deallocate ( scalars(ng)%IJwater )
      END IF
# endif
      IF (associated(scalars(ng)%Cs_r)) THEN
        deallocate ( scalars(ng)%Cs_r )
      END IF
 
      IF (associated(scalars(ng)%Cs_w)) THEN
        deallocate ( scalars(ng)%Cs_w )
      END IF
 
      IF (associated(scalars(ng)%sc_r)) THEN
        deallocate ( scalars(ng)%sc_r )
      END IF
 
      IF (associated(scalars(ng)%sc_w)) THEN
        deallocate ( scalars(ng)%sc_w )
      END IF
#endif
!
!-----------------------------------------------------------------------
!  Deallocate derived-type SCALARS structure.
!-----------------------------------------------------------------------
!
      IF (ng.eq.ngrids) THEN
        IF (allocated(scalars)) deallocate ( scalars )
      END IF
!
!-----------------------------------------------------------------------
!  Deallocate variables in modules.
!-----------------------------------------------------------------------
!
!  Deallocate variables that require special treatment in shared-memory.
!  These variables are private for each thread to avoid collisions.
!
!$OMP PARALLEL
      IF (allocated(predictor_2d_step)) deallocate ( predictor_2d_step )
 
      IF (allocated(indx1))             deallocate ( indx1 )
 
      IF (allocated(iic))               deallocate ( iic )
 
      IF (allocated(iif))               deallocate ( iif )
 
      IF (allocated(next_kstp))         deallocate ( next_kstp )
 
      IF (allocated(frcrec))            deallocate ( frcrec )
 
      IF (allocated(sorec))             deallocate ( sorec )
 
      IF (allocated(synchro_flag))      deallocate ( synchro_flag )
 
      IF (allocated(step_counter))      deallocate ( step_counter )
 
      IF (allocated(tdays))             deallocate ( tdays )
 
      IF (allocated(time))              deallocate ( time )
 
      IF (allocated(time_code))         deallocate ( time_code )
!$OMP END PARALLEL
!
!  Deallocate regular variables
!
#ifdef T_PASSIVE
      IF (associated(inert))            deallocate ( inert )
#endif
      IF (allocated(perfectrst))        deallocate ( perfectrst )
 
      IF (allocated(ndtfast))           deallocate ( ndtfast )
 
      IF (allocated(nfast))             deallocate ( nfast )
 
      IF (allocated(dt))                deallocate ( dt )
 
      IF (allocated(dtfast))            deallocate ( dtfast )
 
      IF (allocated(timeend))           deallocate ( timeend )
 
      IF (allocated(avgtime))           deallocate ( avgtime )
 
      IF (allocated(diatime))           deallocate ( diatime )
 
      IF (allocated(imptime))           deallocate ( imptime )
 
      IF (allocated(initime))           deallocate ( initime )
 
      IF (allocated(initimes))          deallocate ( initimes )
 
      IF (allocated(obstime))           deallocate ( obstime )
 
      IF (allocated(frctime))           deallocate ( frctime )
 
      IF (allocated(ntimes))            deallocate ( ntimes )
 
      IF (allocated(first_time))        deallocate ( first_time )
 
      IF (allocated(ntfirst))           deallocate ( ntfirst )
 
      IF (allocated(ntstart))           deallocate ( ntstart )
 
      IF (allocated(ntend))             deallocate ( ntend )
 
      IF (allocated(nrecfrc))           deallocate ( nrecfrc )
 
      IF (allocated(setgridconfig))     deallocate ( setgridconfig )
 
      IF (allocated(processinputdata))  deallocate ( processinputdata )
 
      IF (allocated(nudgingcoeff))      deallocate ( nudgingcoeff )
 
      IF (allocated(obcdata))           deallocate ( obcdata )
 
      IF (allocated(lbiology))          deallocate ( lbiology )
 
      IF (allocated(lfloats))           deallocate ( lfloats )
 
      IF (allocated(lsediment))         deallocate ( lsediment )
 
      IF (allocated(lstations))         deallocate ( lstations )
 
      IF (allocated(compositegrid))     deallocate ( compositegrid )
 
      IF (allocated(refinedgrid))       deallocate ( refinedgrid )
 
      IF (allocated(refinescale))       deallocate ( refinescale )
 
      IF (allocated(getdonordata))      deallocate ( getdonordata )
 
      IF (allocated(ewperiodic))        deallocate ( ewperiodic )
 
      IF (allocated(nsperiodic))        deallocate ( nsperiodic )
 
      IF (allocated(volcons))           deallocate ( volcons )
 
#if defined ADJOINT || defined TANGENT || defined TL_IOMS
      IF (allocated(ad_volcons))        deallocate ( ad_volcons )
 
      IF (allocated(tl_volcons))        deallocate ( tl_volcons )
#endif
 
      IF (allocated(lsponge))           deallocate ( lsponge )
 
      IF (allocated(luvsponge))         deallocate ( luvsponge )
 
      IF (allocated(ltracersponge))     deallocate ( ltracersponge )
 
      IF (allocated(clm_file))          deallocate ( clm_file )
 
      IF (allocated(lclimatology))      deallocate ( lclimatology )
 
      IF (allocated(lsshclm))           deallocate ( lsshclm )
 
      IF (allocated(lm2clm))            deallocate ( lm2clm )
 
      IF (allocated(lm3clm))            deallocate ( lm3clm )
 
      IF (allocated(ltracerclm))        deallocate ( ltracerclm )
 
      IF (allocated(lnudging))          deallocate ( lnudging )
 
      IF (allocated(lnudgem2clm))       deallocate ( lnudgem2clm )
 
      IF (allocated(lnudgem3clm))       deallocate ( lnudgem3clm )
 
      IF (allocated(lnudgetclm))        deallocate ( lnudgetclm )
 
      IF (allocated(luvsrc))            deallocate ( luvsrc )
 
      IF (allocated(lwsrc))             deallocate ( lwsrc )
 
      IF (allocated(ltracersrc))        deallocate ( ltracersrc )
 
      IF (allocated(maxspeed))          deallocate ( maxspeed )
 
      IF (allocated(maxrho))            deallocate ( maxrho )
 
      IF (allocated(levsfrc))           deallocate ( levsfrc )
 
      IF (allocated(levbfrc))           deallocate ( levbfrc )
 
      IF (allocated(vtransform))        deallocate ( vtransform )
 
      IF (allocated(vstretching))       deallocate ( vstretching )
 
      IF (allocated(tcline))            deallocate ( tcline )
 
      IF (allocated(hc))                deallocate ( hc )
 
      IF (allocated(theta_s))           deallocate ( theta_s )
 
      IF (allocated(theta_b))           deallocate ( theta_b )
 
      IF (allocated(hmin))              deallocate ( hmin )
 
      IF (allocated(hmax))              deallocate ( hmax )
 
      IF (allocated(xl))                deallocate ( xl )
 
      IF (allocated(el))                deallocate ( el )
 
      IF (allocated(lonmin))            deallocate ( lonmin )
 
      IF (allocated(lonmax))            deallocate ( lonmax )
 
      IF (allocated(latmin))            deallocate ( latmin )
 
      IF (allocated(latmax))            deallocate ( latmax )
 
      IF (allocated(idigits))           deallocate ( idigits )
 
      IF (allocated(jdigits))           deallocate ( jdigits )
 
#ifdef SOLVE3D
      IF (allocated(kdigits))           deallocate ( kdigits )
#endif
 
      IF (allocated(totvolume))         deallocate ( totvolume )
 
      IF (allocated(minvolume))         deallocate ( minvolume )
 
      IF (allocated(maxvolume))         deallocate ( maxvolume )
 
      IF (allocated(dxmin))             deallocate ( dxmin )
 
      IF (allocated(dxmax))             deallocate ( dxmax )
 
      IF (allocated(dymin))             deallocate ( dymin )
 
      IF (allocated(dymax))             deallocate ( dymax )
 
#ifdef MASKING
      IF (allocated(dxminw))            deallocate ( dxminw )
 
      IF (allocated(dxmaxw))            deallocate ( dxmaxw )
 
      IF (allocated(dyminw))            deallocate ( dyminw )
 
      IF (allocated(dymaxw))            deallocate ( dymaxw )
 
#endif
 
#ifdef SOLVE3D
      IF (allocated(dzmin))             deallocate ( dzmin )
 
      IF (allocated(dzmax))             deallocate ( dzmax )
 
# ifdef MASKING
      IF (allocated(dzminw))            deallocate ( dzminw )
 
      IF (allocated(dzmaxw))            deallocate ( dzmaxw )
# endif
#endif
 
      IF (allocated(grdmax))            deallocate ( grdmax )
 
#ifdef DIFF_3DCOEF
      IF (allocated(diffmin))           deallocate ( diffmin )
 
      IF (allocated(diffmax))           deallocate ( diffmax )
#endif
 
      IF (allocated(cg_min))            deallocate ( cg_min )
 
      IF (allocated(cg_max))            deallocate ( cg_max )
 
      IF (allocated(cg_cor))            deallocate ( cg_cor )
 
#ifdef VISC_3DCOEF
      IF (allocated(viscmin))           deallocate ( viscmin )
 
      IF (allocated(viscmax))           deallocate ( viscmax )
#endif
 
      IF (allocated(r0))                deallocate ( r0 )
 
      IF (allocated(tcoef))             deallocate ( tcoef )
 
      IF (allocated(scoef))             deallocate ( scoef )
 
      IF (allocated(t0))                deallocate ( t0 )
 
      IF (allocated(s0))                deallocate ( s0 )
 
      IF (allocated(gamma2))            deallocate ( gamma2 )
 
      IF (allocated(lmd_jwt))           deallocate ( lmd_jwt )
 
      IF (allocated(rx0))               deallocate ( rx0 )
 
      IF (allocated(rx1))               deallocate ( rx1 )
 
      IF (allocated(rdrg))              deallocate ( rdrg )
 
      IF (allocated(rdrg2))             deallocate ( rdrg2 )
 
      IF (allocated(zos))               deallocate ( zos )
 
      IF (allocated(zob))               deallocate ( zob )
 
      IF (allocated(dcrit))             deallocate ( dcrit )
 
#ifdef PROPAGATOR
      IF (allocated(nconv))             deallocate ( nconv )
#endif
 
      IF (allocated(weight))            deallocate ( weight )
 
      IF (allocated(akk_bak))           deallocate ( akk_bak )
 
      IF (allocated(akp_bak))           deallocate ( akp_bak )
 
      IF (allocated(akv_bak))           deallocate ( akv_bak )
 
      IF (allocated(akv_limit))         deallocate ( akv_limit )
 
      IF (allocated(ad_visc2))          deallocate ( ad_visc2 )
 
      IF (allocated(nl_visc2))          deallocate ( nl_visc2 )
 
      IF (allocated(tl_visc2))          deallocate ( tl_visc2 )
 
      IF (allocated(visc2))             deallocate ( visc2 )
 
      IF (allocated(ad_visc4))          deallocate ( ad_visc4 )
 
      IF (allocated(nl_visc4))          deallocate ( nl_visc4 )
 
      IF (allocated(tl_visc4))          deallocate ( tl_visc4 )
 
      IF (allocated(visc4))             deallocate ( visc4 )
 
      IF (allocated(tkenu2))            deallocate ( tkenu2 )
 
      IF (allocated(tkenu4))            deallocate ( tkenu4 )
 
      IF (allocated(akt_bak))           deallocate ( akt_bak )
 
      IF (allocated(akt_limit))         deallocate ( akt_limit )
 
      IF (allocated(kdiff))             deallocate ( kdiff )
 
      IF (allocated(ad_tnu2))           deallocate ( ad_tnu2 )
 
      IF (allocated(nl_tnu2))           deallocate ( nl_tnu2 )
 
      IF (allocated(tl_tnu2))           deallocate ( tl_tnu2 )
 
      IF (allocated(tnu2))              deallocate ( tnu2 )
 
      IF (allocated(ad_tnu4))           deallocate ( ad_tnu4 )
 
      IF (allocated(nl_tnu4))           deallocate ( nl_tnu4 )
 
      IF (allocated(tl_tnu4))           deallocate ( tl_tnu4 )
 
      IF (allocated(tnu4))              deallocate ( tnu4 )
 
      IF (allocated(tl_m2diff))         deallocate ( tl_m2diff )
 
      IF (allocated(tl_m3diff))         deallocate ( tl_m3diff )
 
      IF (allocated(tl_tdiff))          deallocate ( tl_tdiff )
 
      IF (allocated(ad_akv_fac))        deallocate ( ad_akv_fac )
 
      IF (allocated(tl_akv_fac))        deallocate ( tl_akv_fac )
 
      IF (allocated(ad_akt_fac))        deallocate ( ad_akt_fac )
 
      IF (allocated(tl_akt_fac))        deallocate ( tl_akt_fac )
 
      IF (allocated(ldefadj))           deallocate ( ldefadj )
 
      IF (allocated(ldefavg))           deallocate ( ldefavg )
 
      IF (allocated(ldefdai))           deallocate ( ldefdai )
 
      IF (allocated(ldefdia))           deallocate ( ldefdia )
 
      IF (allocated(ldeferr))           deallocate ( ldeferr )
 
      IF (allocated(ldefflt))           deallocate ( ldefflt )
 
      IF (allocated(ldefhis))           deallocate ( ldefhis )
 
      IF (allocated(ldefhss))           deallocate ( ldefhss )
 
      IF (allocated(ldefini))           deallocate ( ldefini )
 
      IF (allocated(ldefirp))           deallocate ( ldefirp )
 
      IF (allocated(ldefitl))           deallocate ( ldefitl )
 
      IF (allocated(ldeflcz))           deallocate ( ldeflcz )
 
      IF (allocated(ldeflze))           deallocate ( ldeflze )
 
      IF (allocated(ldefmod))           deallocate ( ldefmod )
 
      IF (allocated(ldefqck))           deallocate ( ldefqck )
 
      IF (allocated(ldefrst))           deallocate ( ldefrst )
 
#ifdef STD_MODEL
      IF (allocated(ldefstd))           deallocate ( ldefstd )
#endif
 
#ifdef SP4DVAR
      IF (allocated(ldefsca))           deallocate ( ldefsca )
 
      IF (allocated(ldefsct))           deallocate ( ldefsct )
 
      IF (allocated(ldefspa))           deallocate ( ldefspa )
 
      IF (allocated(ldefspt))           deallocate ( ldefspt )
#endif
 
      IF (allocated(ldefsta))           deallocate ( ldefsta )
 
      IF (allocated(ldeftide))          deallocate ( ldeftide )
 
      IF (allocated(ldeftlf))           deallocate ( ldeftlf )
 
      IF (allocated(ldeftlm))           deallocate ( ldeftlm )
 
      IF (allocated(ldefxtr))           deallocate ( ldefxtr )
 
      IF (allocated(lreadadm))          deallocate ( lreadadm )
 
      IF (allocated(lreadblk))          deallocate ( lreadblk )
 
      IF (allocated(lreadfrc))          deallocate ( lreadfrc )
 
      IF (allocated(lreadfwd))          deallocate ( lreadfwd )
 
      IF (allocated(lreadqck))          deallocate ( lreadqck )
 
      IF (allocated(lreadstd))          deallocate ( lreadstd )
 
      IF (allocated(lreadtlm))          deallocate ( lreadtlm )
 
      IF (allocated(lwrtadj))           deallocate ( lwrtadj )
 
      IF (allocated(lwrtavg))           deallocate ( lwrtavg )
 
      IF (allocated(lwrtdia))           deallocate ( lwrtdia )
 
      IF (allocated(lwrthis))           deallocate ( lwrthis )
 
      IF (allocated(lwrtper))           deallocate ( lwrtper )
 
      IF (allocated(lwrtqck))           deallocate ( lwrtqck )
 
      IF (allocated(lwrtrst))           deallocate ( lwrtrst )
 
#ifdef STD_MODEL
      IF (allocated(lwrtstd))           deallocate ( lwrtstd )
#endif
 
      IF (allocated(lwrttlf))           deallocate ( lwrttlf )
 
      IF (allocated(lwrttlm))           deallocate ( lwrttlm )
 
      IF (allocated(lwrtxtr))           deallocate ( lwrtxtr )
 
 
      IF (allocated(ldefnrm))           deallocate ( ldefnrm )
 
      IF (allocated(lwrtnrm))           deallocate ( lwrtnrm )
 
#ifdef AD_OUTPUT_STATE
      IF (allocated(lwrtstate3d))       deallocate ( lwrtstate3d )
#endif
 
#if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
      IF (allocated(soinitial))         deallocate ( soinitial )
#endif
 
      IF (allocated(lwrtstate2d))       deallocate ( lwrtstate2d )
 
      IF (allocated(lwrttime))          deallocate ( lwrttime )
 
      IF (allocated(lwrtcost))          deallocate ( lwrtcost )
 
      IF (allocated(ladjusted))         deallocate ( ladjusted )
 
      IF (allocated(lprocessobc))       deallocate ( lprocessobc )
 
      IF (allocated(lprocesstides))     deallocate ( lprocesstides )
 
      IF (allocated(lwrtinfo))          deallocate ( lwrtinfo )
 
      IF (allocated(ldefout))           deallocate ( ldefout )
 
      IF (allocated(ndefadj))           deallocate ( ndefadj )
 
      IF (allocated(ndefavg))           deallocate ( ndefavg )
 
      IF (allocated(ndefdia))           deallocate ( ndefdia )
 
      IF (allocated(ndefhis))           deallocate ( ndefhis )
 
      IF (allocated(ndefqck))           deallocate ( ndefqck )
 
      IF (allocated(ndeftlm))           deallocate ( ndeftlm )
 
      IF (allocated(ndeftlm))           deallocate ( ndefxtr )
 
      IF (allocated(ntsavg))            deallocate ( ntsavg )
 
      IF (allocated(ntsdia))            deallocate ( ntsdia )
 
      IF (allocated(nadj))              deallocate ( nadj )
 
      IF (allocated(navg))              deallocate ( navg )
 
      IF (allocated(ndia))              deallocate ( ndia )
 
      IF (allocated(nflt))              deallocate ( nflt )
 
      IF (allocated(nhis))              deallocate ( nhis )
 
      IF (allocated(nqck))              deallocate ( nqck )
 
      IF (allocated(nrst))              deallocate ( nrst )
 
      IF (allocated(nsta))              deallocate ( nsta )
 
      IF (allocated(ntlm))              deallocate ( ntlm )
 
#ifdef SP4DVAR
      IF (allocated(nsca))              deallocate ( nsca )
 
      IF (allocated(nsct))              deallocate ( nsct )
 
      IF (allocated(nspa))              deallocate ( nspa )
 
      IF (allocated(nspt))              deallocate ( nspt )
 
      IF (allocated(crec))              deallocate ( crec )
#endif
 
      IF (allocated(extractflag))       deallocate ( extractflag )
 
      IF (allocated(ninfo))             deallocate ( ninfo )
 
      IF (allocated(nobc))              deallocate ( nobc )
 
      IF (allocated(nbrec))             deallocate ( nbrec )
 
      IF (allocated(obccount))          deallocate ( obccount )
 
#ifdef ADJUST_BOUNDARY
      IF (allocated(lobc))              deallocate ( lobc )
#endif
 
      IF (allocated(nsff))              deallocate ( nsff )
 
      IF (allocated(nfrec))             deallocate ( nfrec )
 
      IF (allocated(sfcount))           deallocate ( sfcount )
 
#ifdef ADJUST_STFLUX
      IF (allocated(lstflux))           deallocate ( lstflux )
#endif
 
      IF (allocated(nrrec))             deallocate ( nrrec )
 
      IF (allocated(lastrec))           deallocate ( lastrec )
 
      IF (allocated(lcycleadj))         deallocate ( lcycleadj )
 
      IF (allocated(lcyclerst))         deallocate ( lcyclerst )
 
      IF (allocated(lcycletlm))         deallocate ( lcycletlm )
 
#if defined AVERAGES  && defined AVERAGES_DETIDE && \
   (defined ssh_tides || defined uv_tides)
      IF (allocated(hcount))            deallocate ( hcount )
#endif
 
      IF (allocated(kstrs))             deallocate ( kstrs )
 
      IF (allocated(kends))             deallocate ( kends )
 
      IF (allocated(dstrs))             deallocate ( dstrs )
 
      IF (allocated(dends))             deallocate ( dends )
 
      IF (allocated(trnorm))            deallocate ( trnorm )
 
      IF (allocated(so_decay))          deallocate ( so_decay )
 
      IF (allocated(so_sdev))           deallocate ( so_sdev )
 
#if defined FOUR_DVAR || defined VERIFICATION
      IF (allocated(cnorm))             deallocate ( cnorm )
 
      IF (allocated(cnormb))            deallocate ( cnormb )
 
      IF (allocated(sporadicimpulse))   deallocate ( sporadicimpulse )
 
      IF (allocated(frequentimpulse))   deallocate ( frequentimpulse )
 
      IF (allocated(dtdz_min))          deallocate ( dtdz_min )
 
      IF (allocated(ml_depth))          deallocate ( ml_depth )
 
      IF (allocated(lnm_depth))         deallocate ( lnm_depth )
 
      IF (allocated(balance))           deallocate ( balance )
 
      IF (allocated(hdecay))            deallocate ( hdecay )
 
      IF (allocated(vdecay))            deallocate ( vdecay )
 
      IF (allocated(tdecay))            deallocate ( tdecay )
 
      IF (allocated(hdecayb))           deallocate ( hdecayb )
 
      IF (allocated(vdecayb))           deallocate ( vdecayb )
 
      IF (allocated(bgqc_type))         deallocate ( bgqc_type )
 
      IF (allocated(nprovenance))       deallocate ( nprovenance )
 
      IF (allocated(s_bgqc))            deallocate ( s_bgqc )
 
# ifdef STD_MODEL
      IF (allocated(sigma_max))         deallocate ( sigma_max )
 
      IF (allocated(sigma_ml))          deallocate ( sigma_ml )
 
      IF (allocated(sigma_do))          deallocate ( sigma_do )
 
      IF (allocated(sigma_dz))          deallocate ( sigma_dz )
 
      IF (allocated(mld_uniform))       deallocate ( mld_uniform )
# endif
#endif
 
      IF (allocated(obcfac))            deallocate ( obcfac )
 
      IF (allocated(fsobc_in))          deallocate ( fsobc_in )
 
      IF (allocated(fsobc_out))         deallocate ( fsobc_out )
 
      IF (allocated(m2obc_in))          deallocate ( m2obc_in )
 
      IF (allocated(m2obc_out))         deallocate ( m2obc_out )
 
#ifdef SOLVE3D
      IF (allocated(m3obc_in))          deallocate ( m3obc_in )
 
      IF (allocated(m3obc_out))         deallocate ( m3obc_out )
 
      IF (allocated(tobc_in))           deallocate ( tobc_in )
 
      IF (allocated(tobc_out))          deallocate ( tobc_out )
#endif
 
      IF (allocated(znudg))             deallocate ( znudg )
 
      IF (allocated(m2nudg))            deallocate ( m2nudg )
 
      IF (allocated(m3nudg))            deallocate ( m3nudg )
 
      IF (allocated(tnudg))             deallocate ( tnudg )
 
#ifdef BULK_FLUXES
      IF (allocated(blk_zq))            deallocate ( blk_zq )
 
      IF (allocated(blk_zt))            deallocate ( blk_zt )
 
      IF (allocated(blk_zw))            deallocate ( blk_zw )
#endif
 
      IF (allocated(gls_m))             deallocate ( gls_m )
 
      IF (allocated(gls_n))             deallocate ( gls_n )
 
      IF (allocated(gls_p))             deallocate ( gls_p )
 
      IF (allocated(gls_sigk))          deallocate ( gls_sigk )
 
      IF (allocated(gls_sigp))          deallocate ( gls_sigp )
 
      IF (allocated(gls_cmu0))          deallocate ( gls_cmu0 )
 
      IF (allocated(gls_cmupr))         deallocate ( gls_cmupr )
 
      IF (allocated(gls_c1))            deallocate ( gls_c1 )
 
      IF (allocated(gls_c2))            deallocate ( gls_c2 )
 
      IF (allocated(gls_c3m))           deallocate ( gls_c3m )
 
      IF (allocated(gls_c3p))           deallocate ( gls_c3p )
 
      IF (allocated(gls_kmin))          deallocate ( gls_kmin )
 
      IF (allocated(gls_pmin))          deallocate ( gls_pmin )
 
      IF (allocated(charnok_alpha))     deallocate ( charnok_alpha )
 
      IF (allocated(zos_hsig_alpha))    deallocate ( zos_hsig_alpha )
 
      IF (allocated(sz_alpha))          deallocate ( sz_alpha )
 
      IF (allocated(crgban_cw))         deallocate ( crgban_cw )
!
      RETURN
      END SUBROUTINE deallocate_scalars
!
      SUBROUTINE initialize_scalars
!
!=======================================================================
!                                                                      !
!  This routine initializes several variables in module for all nested !
!  grids.                                                              !
!                                                                      !
!=======================================================================
!
      USE mod_param
!
!  Local variable declarations.
!
      integer :: i, ic, j, ng, itrc
 
      real(r8) :: one, zero
      real(r8), parameter :: inival = 0.0_r8
!
!---------------------------------------------------------------------
!  Set tracer identification indices.
!---------------------------------------------------------------------
!
      itemp=1
      isalt=2
      ic=nat
 
#ifdef T_PASSIVE
!
!  Indices for inert passive tracers to advect and diffuse.
!
      DO i=1,npt
        ic=ic+1
        inert(i)=ic
      END DO
#endif
 
#ifdef DIAGNOSTICS
!
!---------------------------------------------------------------------
!  Set diagnostic fields identification indices.
!---------------------------------------------------------------------
 
# ifdef DIAGNOSTICS_TS
!
!  Indices for tracer diagnostic variables.
!
      ithadv=1
      itxadv=2
      ityadv=3
      itvadv=4
      ic=4
#  if defined TS_DIF2 || defined TS_DIF4
      ithdif=ic+1
      itxdif=ic+2
      itydif=ic+3
      ic=ic+3
#   if defined MIX_GEO_TS || defined MIX_ISO_TS
      itsdif=ic+1
      ic=ic+1
#   endif
#  endif
      itvdif=ic+1
      itrate=ic+2
# endif
# ifdef DIAGNOSTICS_UV
!
!  Indices for 2D momentum diagnostic variables.  In some places in
!  the code a compact DO-loop (idiag=1:M2pgrd) is used to improve
!  flexibility. Therefore, the order of indices is very important.
!  Only those fields that require special treatment are set below
!  the M2pgrd index.
!
      ic=0
#  if defined UV_COR
      m2fcor=ic+1
      ic=ic+1
#  endif
#  if defined UV_ADV
      m2hadv=ic+1
      m2xadv=ic+2
      m2yadv=ic+3
      ic=ic+3
#  endif
#  if defined WEC_VF
#   if defined UV_COR
      m2fsco=ic+1
      ic=ic+1
#   endif
#   ifdef BOTTOM_STREAMING
      m2bstm=ic+1
      ic=ic+1
#   endif
#   ifdef SURFACE_STREAMING
      m2sstm=ic+1
      ic=ic+1
#   endif
      m2hjvf=ic+1
      m2kvrf=ic+2
      m2wrol=ic+3
      m2wbrk=ic+4
      ic=ic+4
#  endif
#  if defined VEGETATION && defined VEG_DRAG
      m2fveg=ic+1
      ic=ic+1
#  endif
#  if defined UV_VIS2 || defined UV_VIS4
      m2hvis=ic+1
      m2xvis=ic+2
      m2yvis=ic+3
      ic=ic+3
#  endif
      m2pgrd=ic+1
      m2sstr=ic+2             ! These indices need to be
      m2bstr=ic+3             ! specified last to allow a
#  if defined WEC_VF
      m2zeta=ic+4
      m2zetw=ic+5
      m2zqsp=ic+6
      m2zbeh=ic+7
#  endif
      m2rate=ndm2d            ! compact DO-loop structure
#  ifdef SOLVE3D
!
!  Indices for 3D momentum diagnostic variables.  In some places in
!  the code a compact DO-loop (idiag=1:M3pgrd) is used to improve
!  flexibility. Therefore, the order of indices is very important.
!  Only those fields that require special treatment are set below
!  the M3pgrd index.
!
      ic=0
#   if defined UV_COR
      m3fcor=ic+1
      ic=ic+1
#   endif
#   if defined UV_ADV
      m3vadv=ic+1
      m3hadv=ic+2
      m3xadv=ic+3
      m3yadv=ic+4
      ic=ic+4
#   endif
#   if defined WEC_VF
#    if defined UV_COR
      m3fsco=ic+1
      ic=ic+1
#    endif
#    ifdef BOTTOM_STREAMING
      m3bstm=ic+1
      ic=ic+1
#    endif
#    ifdef SURFACE_STREAMING
      m3sstm=ic+1
      ic=ic+1
#    endif
      m3vjvf=ic+1
      m3hjvf=ic+2
      m3kvrf=ic+3
      m3wrol=ic+4
      m3wbrk=ic+5
      ic=ic+5
#   endif
#   if defined VEGETATION && defined VEG_DRAG
      m3fveg=ic+1
      ic=ic+1
#   endif
      m3pgrd=ic+1             ! needs to be here, indices below
      m3vvis=ic+2             ! require special treatment
#   if defined UV_VIS2 || defined UV_VIS4
      m3hvis=ic+3
      m3xvis=ic+4
      m3yvis=ic+5
#   endif
      m3rate=ndm3d
#  endif
# endif
#endif
!
!-----------------------------------------------------------------------
!  Activate all computation control switches.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        lastrec(ng)=.false.
        compositegrid(1:4,ng)=.false.
        refinedgrid(ng)=.false.
        getdonordata(ng)=.false.
        lbiology(ng)=.true.
        lcycleadj(ng)=.false.
        lcyclerst(ng)=.false.
        lcycletlm(ng)=.false.
        lfloats(ng)=.true.
        lsediment(ng)=.true.
        lstations(ng)=.true.
#if defined FOUR_DVAR || defined VERIFICATION
        frequentimpulse(ng)=.false.
        sporadicimpulse(ng)=.false.
#endif
      END DO
 
#ifdef ADJUST_BOUNDARY
!
!-----------------------------------------------------------------------
!  Initilize switches to process open boundary arrays for 4DVar
!  adjustments.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO j=1,mstatevar
          DO i=1,4
            lobc(i,j,ng)=.false.
          END DO
        END DO
      END DO
#endif
#ifdef ADJUST_STFLUX
!
!-----------------------------------------------------------------------
!  Initilize switches to process surface tracer fluexes for 4DVar
!  adjustments.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO i=1,mt
          lstflux(i,ng)=.false.
        END DO
      END DO
#endif
!
!-----------------------------------------------------------------------
!  Initialize several scalar variables.
!-----------------------------------------------------------------------
!
      one=1.0_r8
      zero=0.0_r8
      co=1.0_r8/(2.0_r8+sqrt(2.0_r8))
      gorho0=g/rho0
      DO ng=1,ngrids
        ewperiodic(ng)=.false.
        nsperiodic(ng)=.false.
        nudgingcoeff(ng)=.false.
        obcdata(ng)=.false.
        setgridconfig(ng)=.true.
        processinputdata(ng)=.true.
        extractflag(ng)=0
        refinescale(ng)=0
        gamma2(ng)=-1.0_r8
        vtransform(ng)=1
        vstretching(ng)=1
#if defined AVERAGES  && defined AVERAGES_DETIDE && \
   (defined ssh_tides || defined uv_tides)
        hcount(ng)=0
#endif
#ifdef ADJUST_BOUNDARY
        obccount(ng)=0
#endif
#if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
        sfcount(ng)=0
#endif
        first_time(ng)=0
        idigits(ng)=int(log10(real(lm(ng),r8)))+1
        jdigits(ng)=int(log10(real(mm(ng),r8)))+1
#ifdef SOLVE3D
        kdigits(ng)=int(log10(real(n(ng),r8)))+1
#endif
        maxspeed(ng)=-large
        maxrho(ng)=-large
        totvolume(ng)=0.0_dp
        minvolume(ng)= large
        maxvolume(ng)=-large
        dxmin(ng)= large
        dxmax(ng)=-large
        dymin(ng)= large
        dymax(ng)=-large
#ifdef MASKING
        dxminw(ng)= large
        dxmaxw(ng)=-large
        dyminw(ng)= large
        dymaxw(ng)=-large
#endif
#ifdef SOLVE3D
        dzmin(ng)= large
        dzmax(ng)=-large
# ifdef MASKING
        dzminw(ng)= large
        dzmaxw(ng)=-large
# endif
#endif
        grdmax(ng)=-large
#ifdef DIFF_3DCOEF
        diffmin(ng)= large
        diffmax(ng)=-large
#endif
        cg_min(ng)= large
        cg_max(ng)=-large
        cg_cor(ng)=-large
#ifdef VISC_3DCOEF
        viscmin(ng)= large
        viscmax(ng)=-large
#endif
        rx0(ng)=-large
        rx1(ng)=-large
        clm_file(ng)=.false.
        lnudging(ng)=.false.
        lnudgem2clm(ng)=.false.
        lnudgem3clm(ng)=.false.
        lclimatology(ng)=.false.
        lm2clm(ng)=.false.
        lm3clm(ng)=.false.
        lsshclm(ng)=.false.
        lsponge(ng)=.false.
        luvsponge(ng)=.false.
        luvsrc(ng)=.false.
        lwsrc(ng)=.false.
        DO itrc=1,mt
          lnudgetclm(itrc,ng)=.false.
          ltracerclm(itrc,ng)=.false.
          ltracersrc(itrc,ng)=.false.
          ltracersponge(itrc,ng)=.false.
          ad_akt_fac(itrc,ng)=1.0_r8
          tl_akt_fac(itrc,ng)=1.0_r8
          ad_tnu2(itrc,ng)=inival
          nl_tnu2(itrc,ng)=inival
          tl_tnu2(itrc,ng)=inival
          tnu2(itrc,ng)=inival
          ad_tnu4(itrc,ng)=inival
          nl_tnu4(itrc,ng)=inival
          tl_tnu4(itrc,ng)=inival
          tnu4(itrc,ng)=inival
        END DO
        DO itrc=1,nat
          akt_limit(itrc,ng)=1.0e-3_r8
        END DO
        akv_limit(ng)=1.0e-3_r8
        ad_akv_fac(ng)=1.0_r8
        tl_akv_fac(ng)=1.0_r8
        ad_visc2(ng)=inival
        nl_visc2(ng)=inival
        tl_visc2(ng)=inival
        visc2(ng)=inival
        ad_visc4(ng)=inival
        nl_visc4(ng)=inival
        tl_visc4(ng)=inival
        visc4(ng)=inival
#ifdef BULK_FLUXES
        blk_zq(ng)=10.0_r8
        blk_zt(ng)=10.0_r8
        blk_zw(ng)=10.0_r8
#endif
        DO i=1,4
          volcons(i,ng)=.false.
#if defined ADJOINT || defined TANGENT || defined TL_IOMS
          ad_volcons(i,ng)=.false.
          tl_volcons(i,ng)=.false.
#endif
          fsobc_in(ng,i)=0.0_dp
          fsobc_out(ng,i)=0.0_dp
          m2obc_in(ng,i)=0.0_dp
          m2obc_out(ng,i)=0.0_dp
#ifdef SOLVE3D
          m3obc_in(ng,i)=0.0_dp
          m3obc_out(ng,i)=0.0_dp
#endif
        END DO
      END DO
#ifdef SOLVE3D
      tobc_in = 0.0_dp
      tobc_out = 0.0_dp
#endif
!
!  Initialize blowup string.
!
      DO i=1,len(blowup_string)
        blowup_string(i:i)=' '
      END DO
!
!  Initialize thread private variables.
!
!$OMP PARALLEL
      synchro_flag=.false.
 
      ntfirst=1
      ntstart=1
      ntend=0
      step_counter=0
!$OMP END PARALLEL
 
#if defined LMD_SKPP || defined LMD_BKPP
!
!  Proportionality coefficient parameterizing boundary layer
!  nonlocal transport.
!
      lmd_cg=lmd_cstar*                                                 &
     &       vonkar*(lmd_cs*vonkar*lmd_epsilon)**(1.0_r8/3.0_r8)
#endif
 
#if defined FOUR_DVAR || defined VERIFICATION
!
!  Initialize error covariace variables.
!
      balance=.false.
 
      cnorm=.false.
      hdecay=inival
      tdecay=inival
      vdecay=inival
 
      cnormb=.false.
      hdecayb=inival
      vdecayb=inival
 
# ifdef STD_MODEL
      sigma_max=inival
      sigma_ml=inival
      sigma_do=inival
      sigma_dz=inival
      mld_uniform=inival
# endif
#endif
!
!  Initialize several IO flags.
!
      lmultigst=.false.
      lrstgst=.false.
 
      DO ng=1,ngrids
 
        perfectrst(ng)=.false.
 
        ladjusted(ng)=.false.
 
        lprocessobc(ng)=.false.
        lprocesstides(ng)=.false.
 
        ldefadj(ng)=.false.
        ldefavg(ng)=.true.
        ldefdai(ng)=.false.
        ldefdia(ng)=.true.
        ldeferr(ng)=.false.
#ifdef FLOATS
        ldefflt(ng)=.true.
#endif
        ldefhis(ng)=.true.
#ifdef FOUR_DVAR
        ldefini(ng)=.true.
#else
        ldefini(ng)=.false.
#endif
        ldefirp(ng)=.false.
        ldefitl(ng)=.false.
        ldefmod(ng)=.false.
        ldefqck(ng)=.false.
        ldefrst(ng)=.true.
        ldefsta(ng)=.true.
#ifdef STD_MODEL
        ldefstd(ng)=.false.
#endif
        ldeftlm(ng)=.false.
#if defined AVERAGES  && defined AVERAGES_DETIDE && \
   (defined ssh_tides || defined uv_tides)
        ldeftide(ng)=.true.
#else
        ldeftide(ng)=.false.
#endif
        lreadadm(ng)=.false.
        lreadblk(ng)=.false.
#ifdef FRC_FILE
        lreadfrc(ng)=.true.
#else
        lreadfrc(ng)=.false.
#endif
        lreadfwd(ng)=.false.
        lreadqck(ng)=.false.
        lreadstd(ng)=.false.
        lreadtlm(ng)=.false.
        lwrtadj(ng)=.false.
        lwrtavg(ng)=.false.
        lwrtdia(ng)=.false.
        lwrthis(ng)=.false.
        lwrtper(ng)=.false.
        lwrtqck(ng)=.false.
        lwrtrst(ng)=.false.
#ifdef STD_MODEL
        lwrtstd(ng)=.false.
#endif
        lwrttlm(ng)=.false.
        lwrtxtr(ng)=.false.
        lwrtinfo(ng)=.true.
        lwrtstate2d(ng)=.false.
#ifdef AD_OUTPUT_STATE
        lwrtstate3d(ng)=.false.
#endif
#if defined STOCHASTIC_OPT && !defined STOCH_OPT_WHITE
        soinitial(ng)=.false.
#endif
        lwrttime(ng)=.true.
        lwrtcost(ng)=.false.
        ldefout(ng)=.false.
 
      END DO
!
!  Initialize the NLM initial conditions time to a negative number to
!  check if its value was assigned elsewhere.  It can be used during
!  the initialization of the adjoint model when DSTART is not the
!  same as the start of the simulation.
!
      DO ng=1,ngrids
        initime(ng)=-1.0_dp
        initimes(ng)=-1.0_dp
      END DO
 
# if defined SG_BBL || defined SSW_BBL
!
!  Nondimensional closure parameters associated with Styles and Glenn
!  (1999) bottom currents and waves boundary layer.
!
      sg_z1p=sg_alpha
      sg_mp=cmplx(sqrt(1.0_r8/(2.0_r8*sg_z1p)),                         &
     &            sqrt(1.0_r8/(2.0_r8*sg_z1p)))
#endif
#if defined GLS_MIXING || defined MY25_MIXING
# if defined CANUTO_A || defined CANUTO_B
!
! Compute parameters for Canuto et al. (2001) stability functions.
!   (Canuto, V.M., Cheng, H.Y., and Dubovikov, M.S., 2001: Ocean
!    turbulence. Part I: One-point closure model - momentum and
!    heat vertical diffusivities, JPO, 1413-1426).
!
      gls_s0=3.0_r8/2.0_r8*gls_l1*gls_l5**2
      gls_s1=-gls_l4*(gls_l6+gls_l7)+2.0_r8*gls_l4*gls_l5*              &
     &        (gls_l1-1.0_r8/3.0_r8*gls_l2-gls_l3)+3.0_r8/2.0_r8*       &
     &        gls_l1*gls_l5*gls_l8
      gls_s2=-3.0_r8/8.0_r8*gls_l1*(gls_l6**2-gls_l7**2)
      gls_s4=2.0_r8*gls_l5
      gls_s5=2.0_r8*gls_l4
      gls_s6=2.0_r8/3.0_r8*gls_l5*(3.0_r8*gls_l3**2-gls_l2**2)-         &
     &       1.0_r8/2.0_r8*gls_l5*gls_l1*(3.0_r8*gls_l3-gls_l2)+        &
     &       3.0_r8/4.0_r8*gls_l1*(gls_l6-gls_l7)
      gls_b0=3.0_r8*gls_l5**2
      gls_b1=gls_l5*(7.0_r8*gls_l4+3.0_r8*gls_l8)
      gls_b2=gls_l5**2*(3.0_r8*gls_l3**2-gls_l2**2)-                    &
     &       3.0_r8/4.0_r8*(gls_l6**2-gls_l7**2)
      gls_b3=gls_l4*(4.0_r8*gls_l4+3.0_r8*gls_l8)
      gls_b5=1.0_r8/4.0_r8*(gls_l2**2-3.0_r8*gls_l3**2)*                &
     &       (gls_l6**2-gls_l7**2)
      gls_b4=gls_l4*(gls_l2*gls_l6-3.0_r8*gls_l3*gls_l7-                &
     &       gls_l5*(gls_l2**2-gls_l3**2))+gls_l5*gls_l8*               &
     &       (3.0_r8*gls_l3**2-gls_l2**2)
# endif
!
!  Coefficients used to compute stability functions for tracer and
!  momentum.
!
      my_b1p2o3=my_b1**(2.0_r8/3.0_r8)
      my_b1pm1o3=1.0_r8/(my_b1**(1.0_r8/3.0_r8))
      my_e1o2=0.5_r8*my_e1
      my_sm1=my_a1*my_a2*((my_b2-3.0_r8*my_a2)*                         &
     &                    (1.0_r8-6.0_r8*my_a1/my_b1)-                  &
     &                    3.0_r8*my_c1*(my_b2+6.0_r8*my_a1))
      my_sm2=9.0_r8*my_a1*my_a2
      my_sh1=my_a2*(1.0_r8-6.0_r8*my_a1/my_b1)
# ifdef KANTHA_CLAYSON
      my_sh2=3.0_r8*my_a2*(6.0_r8*my_a1+my_b2*(1.0_r8-my_c3))
      my_sm4=18.0_r8*my_a1*my_a1+9.0_r8*my_a1*my_a2*(1.0_r8-my_c2)
# else
      my_sh2=3.0_r8*my_a2*(6.0_r8*my_a1+my_b2)
      my_sm3=my_a1*(1.0_r8-3.0_r8*my_c1-6.0_r8*my_a1/my_b1)
      my_sm4=18.0_r8*my_a1*my_a1+9.0_r8*my_a1*my_a2
# endif
#endif
 
      RETURN
      END SUBROUTINE initialize_scalars
!
      END MODULE mod_scalars