tl_main3d.F

Go to the documentation of this file.

#include "cppdefs.h"
#if defined TANGENT && defined SOLVE3D
      SUBROUTINE tl_main3d (RunInterval)
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This routine is the main driver for tangent linear ROMS when        !
!  configurated as a full  3D baroclinic  ocean model.  It  advances   !
!  forward the tangent linear model equations  for all nested grids,   !
!  if any, by the specified time interval (seconds), RunInterval.      !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
# if defined MODEL_COUPLING && defined MCT_LIB
      USE mod_coupler
# endif
      USE mod_iounits
# ifdef NESTING
      USE mod_nesting
# endif
      USE mod_scalars
      USE mod_stepping
!
# ifdef ANA_VMIX
      USE analytical_mod,       ONLY : ana_vmix
# endif
      USE dateclock_mod,        ONLY : time_string
# ifdef TLM_CHECK
      USE dotproduct_mod,       ONLY : tl_dotproduct
# endif
# ifdef TIDE_GENERATING_FORCES
      USE equilibrium_tide_mod, ONLY : equilibrium_tide
# endif
# if defined ATM_COUPLING_NOT_YET && defined MCT_LIB
      USE mct_coupler_mod,      ONLY : ocn2atm_coupling
# endif
# if defined WAV_COUPLING_NOT_YET && defined MCT_LIB
      USE mct_coupler_mod,      ONLY : ocn2wav_coupling
# endif
# if defined FORWARD_READ || defined JEDI
      USE omega_mod,            ONLY : omega
      USE set_depth_mod,        ONLY : set_depth
      USE set_massflux_mod,     ONLY : set_massflux
# endif
      USE strings_mod,          ONLY : founderror
# ifdef BIOLOGY
      USE tl_biology_mod,       ONLY : tl_biology
# endif
# ifdef BBL_MODEL_NOT_YET
!!    USE tl_bbl_mod,           ONLY : tl_bblm
# endif
# if defined BULK_FLUXES_NOT_YET && !defined PRIOR_BULK_FLUXES
!!    USE tl_bulk_flux_mod,     ONLY : tl_bulk_flux
# endif
# ifdef BVF_MIXING_NOT_YET
!!    USE tl_bvf_mix_mod,       ONLY : tl_bvf_mix
# endif
      USE tl_diag_mod,          ONLY : tl_diag
# if defined WEAK_CONSTRAINT || defined FORCING_SV
      USE tl_forcing_mod,       ONLY : tl_forcing
# endif
# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
      USE tl_frc_adjust_mod,    ONLY : tl_frc_adjust
# endif
# ifdef GLS_MIXING_NOT_YET
!!    USE tl_gls_corstep_mod,   ONLY : tl_gls_corstep
!!    USE tl_gls_prestep_mod,   ONLY : tl_gls_prestep
# endif
# ifdef LMD_MIXING_NOT_YET
!!    USE tl_lmd_vmix_mod,      ONLY : tl_lmd_vmix
# endif
# ifdef MY25_MIXING_NOT_YET
!!    USE tl_my25_corstep_mod,  ONLY : tl_my25_corstep
!!    USE tl_my25_prestep_mod,  ONLY : tl_my25_prestep
# endif
# ifdef NESTING
      USE nesting_mod,          ONLY : nesting
      USE tl_nesting_mod,       ONLY : tl_nesting
#  ifndef ONE_WAY
      USE nesting_mod,          ONLY : do_twoway
#  endif
# endif
# ifdef ADJUST_BOUNDARY
      USE tl_obc_adjust_mod,    ONLY : tl_obc_adjust
      USE tl_obc_adjust_mod,    ONLY : tl_obc2d_adjust
      USE tl_set_depth_mod,     ONLY : tl_set_depth_bry
# endif
      USE tl_omega_mod,         ONLY : tl_omega
# if !(defined FORCING_SV || defined JEDI)
      USE tl_post_initial_mod,  ONLY : tl_post_initial
# endif
# ifdef NEARSHORE_MELLOR_NOT_YET
!!    USE tl_radiation_stress_mod, ONLY : tl_radiation_stress
# endif
# ifndef TS_FIXED
      USE tl_rho_eos_mod,       ONLY : tl_rho_eos
# endif
      USE tl_rhs3d_mod,         ONLY : tl_rhs3d
# ifdef SEDIMENT_NOT_YET
!!    USE tl_sediment_mod,      ONLY : tl_sediment
# endif
# ifdef TL_AVERAGES
      USE tl_set_avg_mod,       ONLY : tl_set_avg
# endif
      USE tl_set_depth_mod,     ONLY : tl_set_depth
      USE tl_set_massflux_mod,  ONLY : tl_set_massflux
# if defined SSH_TIDES_NOT_YET || defined UV_TIDES_NOT_YET
!!    USE tl_set_tides_mod,     ONLY : tl_set_tides
# endif
      USE tl_set_vbc_mod,       ONLY : tl_set_vbc
      USE tl_set_zeta_mod,      ONLY : tl_set_zeta
      USE tl_step2d_mod,        ONLY : tl_step2d
# ifndef TS_FIXED
      USE tl_step3d_t_mod,      ONLY : tl_step3d_t
# endif
      USE tl_step3d_uv_mod,     ONLY : tl_step3d_uv
# ifdef FLOATS_NOT_YET
!!    USE tl_step_floats_mod,   ONLY : tl_step_floats
# endif
!!    USE wvelocity_mod,        ONLY : wvelocity
!
      implicit none
!
!  Imported variable declarations.
!
      real(dp), intent(in) :: RunInterval
!
!  Local variable declarations.
!
      logical :: DoNestLayer, Time_Step
 
      integer :: Nsteps, Rsteps
      integer :: ig, il, istep, ng, nl, tile
      integer :: my_iif, next_indx1
# ifdef FLOATS_NOT_YET
      integer :: Lend, Lstr, chunk_size
# endif
!
      real(r8) :: MaxDT, my_StepTime
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__
!
!=======================================================================
!  Time-step tangent linear 3D primitive equations.
!=======================================================================
!
!  Time-step the 3D kernel for the specified time interval (seconds),
!  RunInterval.
!
      time_step=.true.
      donestlayer=.true.
!
      kernel_loop : DO WHILE (time_step)
!
!  In nesting applications, the number of nesting layers (NestLayers) is
!  used to facilitate refinement grids and composite/refinament grids
!  combinations. Otherwise, the solution it is looped once for a single
!  grid application (NestLayers = 1).
!
        nl=0
#ifdef NESTING
        twowayinterval(1:ngrids)=0.0_r8
#endif
!
        nest_layer : DO WHILE (donestlayer)
!
!  Determine number of time steps to compute in each nested grid layer
!  based on the specified time interval (seconds), RunInterval. Non
!  nesting applications have NestLayers=1. Notice that RunInterval is
!  set in the calling driver. Its value may span the full period of the
!  simulation, a multi-model coupling interval (RunInterval > ifac*dt),
!  or just a single step (RunInterval=0).
!
          CALL ntimesteps (itlm, runinterval, nl, nsteps, rsteps)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          IF ((nl.le.0).or.(nl.gt.nestlayers)) EXIT
!
!  Time-step governing equations for Nsteps.
!
          step_loop : DO istep=1,nsteps
!
!  Set time indices and time clock.
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              nstp(ng)=1+mod(iic(ng)-ntstart(ng),2)
              nnew(ng)=3-nstp(ng)
              nrhs(ng)=nstp(ng)
# ifdef JEDI
              jic(ng)=jic(ng)+1
              time4jedi(ng)=time4jedi(ng)+dt(ng)
# endif
              tdays(ng)=time(ng)*sec2day
              IF (step_counter(ng).eq.rsteps) time_step=.false.
            END DO
!
!-----------------------------------------------------------------------
!  Read in required data, if any, from input NetCDF files.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
!$OMP MASTER
              IF (processinputdata(ng)) THEN
                CALL tl_get_data (ng)
              END IF
!$OMP END MASTER
!$OMP BARRIER
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
!
!-----------------------------------------------------------------------
!  If applicable, process input data: time interpolate between data
!  snapshots. Compute BASIC STATE depths and thickness.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=first_tile(ng),last_tile(ng),+1
                IF (processinputdata(ng)) THEN
                  CALL tl_set_data (ng, tile)
                END IF
# if defined FORWARD_READ || defined JEDI
                CALL set_depth (ng, tile, itlm)
# endif
              END DO
!$OMP BARRIER
            END DO
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# if defined FORWARD_READ || defined JEDI
!
!-----------------------------------------------------------------------
!  Compute BASIC STATE horizontal mass fluxes (Hz*u/n and Hz*v/m).
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=last_tile(ng),first_tile(ng),-1
                CALL set_massflux (ng, tile, itlm)
              END DO
!$OMP BARRIER
            END DO
# endif
 
# if (defined WEAK_CONSTRAINT || defined FORCING_SV) && \
     !defined SP4DVAR
!
!-----------------------------------------------------------------------
!  If appropriate, add convolved adjoint solution impulse forcing to
!  the representer model solution. Notice that the forcing is only
!  needed after finishing all inner loops. The forcing is continuous.
!  That is, it is time interpolated at every time-step from available
!  snapshots (FrequentImpulse=TRUE).
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
#  ifdef WEAK_CONSTRAINT
#   ifdef WEAK_NOINTERP
              IF ((iic(ng).gt.1).and.(iic(ng).ne.ntend(ng)+1).and.      &
     &            (mod(iic(ng)-1,nadj(ng)).eq.0)) THEN
                IF (master) THEN
                  WRITE (stdout,*) ' FORCING TLM at iic = ',iic(ng)
                END IF
#   endif
                IF (frequentimpulse(ng)) THEN
                  DO tile=first_tile(ng),last_tile(ng),+1
                    CALL tl_forcing (ng, tile, kstp(ng), nstp(ng))
                    CALL tl_set_depth (ng, tile, itlm)
                  END DO
!$OMP BARRIER
                END IF
#   ifdef WEAK_NOINTERP
              END IF
#   endif
#  else
              DO tile=first_tile(ng),last_tile(ng),+1
                CALL tl_forcing (ng, tile, kstp(ng), nstp(ng))
                CALL tl_set_depth (ng, tile, itlm)
              END DO
!$OMP BARRIER
#  endif
            END DO
# endif
 
# if !(defined FORCING_SV || defined JEDI)
!
!-----------------------------------------------------------------------
!  On the first timestep, it computes the initial depths and level
!  thicknesses from the initial free-surface field. Additionally, it
!  initializes the tangent linear state variables for all time levels
!  and applies lateral boundary conditions.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (iic(ng).eq.ntstart(ng)) THEN
                CALL tl_post_initial (ng, itlm)
              END IF
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Compute horizontal mass fluxes (Hz*u/n and Hz*v/m), density related
!  quatities and report global diagnostics. Compute BASIC STATE omega
!  vertical velocity.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=first_tile(ng),last_tile(ng),+1
                CALL tl_set_massflux (ng, tile, itlm)
# ifndef TS_FIXED
                CALL tl_rho_eos (ng, tile, itlm)
# endif
# ifdef TIDE_GENERATING_FORCES
                CALL equilibrium_tide (ng, tile, itlm)
# endif
                CALL tl_diag (ng, tile)
# if defined FORWARD_READ || defined JEDI
                CALL omega (ng, tile, itlm)
# endif
# ifdef TLM_CHECK
                CALL tl_dotproduct (ng, tile, lnew(ng))
# endif
              END DO
!$OMP BARRIER
            END DO
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# if defined ATM_COUPLING_NOT_YET && defined MCT_LIB
!
!-----------------------------------------------------------------------
!  Couple to atmospheric model every CoupleSteps(Iatmos) timesteps: get
!  air/sea fluxes.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF ((iic(ng).ne.ntstart(ng)).and.                         &
     &            mod(iic(ng)-1,couplesteps(iatmos,ng)).eq.0) THEN
                DO tile=last_tile(ng),first_tile(ng),-1
                  CALL ocn2atm_coupling (ng, tile)
                END DO
!$OMP BARRIER
              END IF
            END DO
# endif
 
# if defined WAV_COUPLING_NOT_YET && defined MCT_LIB
!
!-----------------------------------------------------------------------
!  Couple to waves model every CoupleSteps(Iwaves) timesteps: get
!  waves/sea fluxes.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF ((iic(ng).ne.ntstart(ng)).and.                         &
     &            mod(iic(ng)-1,couplesteps(iwaves,ng)).eq.0) THEN
                DO tile=first_tile(ng),last_tile(ng),+1
                  CALL ocn2wav_coupling (ng, tile)
                END DO
!$OMP BARRIER
              END IF
            END DO
# endif
 
# ifdef NEARSHORE_MELLOR_NOT_YET
!
!-----------------------------------------------------------------------
!  Compute radiation stress terms.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=last_tile(ng),first_tile(ng),-1
                CALL tl_radiation_stress (ng, tile)
              END DO
!$OMP BARRIER
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Set fields for vertical boundary conditions. Process tidal forcing,
!  if any.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=first_tile(ng),last_tile(ng),+1
# if defined BULK_FLUXES_NOT_YET && !defined PRIOR_BULK_FLUXES
                CALL tl_bulk_flux (ng, tile)
# endif
# ifdef BBL_MODEL_NOT_YET
                CALL tl_bblm (ng, tile)
# endif
                CALL tl_set_vbc (ng, tile)
# if defined SSH_TIDES_NOT_YET || defined UV_TIDES_NOT_YET
                CALL tl_set_tides (ng, tile)
# endif
              END DO
!$OMP BARRIER
            END DO
 
# ifdef NESTING
!
!  If composite or mosaic grids, process additional points in the
!  contact zone between connected grids for bottom stress variables.
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (any(compositegrid(:,ng))) THEN
                CALL tl_nesting (ng, itlm, nbstr)
              END IF
            END DO
# endif
 
# ifdef ADJUST_BOUNDARY
!
!-----------------------------------------------------------------------
!  Interpolate open boundary increments and adjust open boundaries.
!  Skip the last output timestep.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (iic(ng).lt.(ntend(ng)+1)) THEN
                DO tile=first_tile(ng),last_tile(ng),+1
                  CALL tl_obc_adjust (ng, tile, lbinp(ng))
                  CALL tl_set_depth_bry (ng, tile, itlm)
                  CALL tl_obc2d_adjust (ng, tile, lbinp(ng))
                END DO
!$OMP BARRIER
              END IF
            END DO
# endif
 
# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
!
!-----------------------------------------------------------------------
!  Interpolate surface forcing increments and adjust surface forcing.
!  Skip the last output timestep.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (iic(ng).lt.(ntend(ng)+1)) THEN
                DO tile=first_tile(ng),last_tile(ng),+1
                  CALL tl_frc_adjust (ng, tile, lfinp(ng))
                END DO
!$OMP BARRIER
              END IF
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Compute tangent linear vertical mixing coefficients for momentum and
!  tracers. Compute S-coordinate vertical velocity, diagnostically from
!  horizontal mass divergence.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=last_tile(ng),first_tile(ng),-1
# if defined ANA_VMIX_NOT_YET
                CALL tl_ana_vmix (ng, tile)
# elif defined LMD_MIXING_NOT_YET
                CALL tl_lmd_vmix (ng, tile)
# elif defined BVF_MIXING_NOT_YET
                CALL tl_bvf_mix (ng, tile)
# endif
                CALL tl_omega (ng, tile, itlm)
!!              CALL wvelocity (ng, tile, nstp(ng))
              END DO
!$OMP BARRIER
            END DO
!
!-----------------------------------------------------------------------
!  Set free-surface to it time-averaged value.  If applicable,
!  accumulate time-averaged output data which needs a irreversible
!  loop in shared-memory jobs.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=first_tile(ng),last_tile(ng),+1     ! irreversible
                CALL tl_set_zeta (ng, tile)
#  ifdef DIAGNOSTICS
!!              CALL set_diags (ng, tile)
#  endif
#  ifdef TL_AVERAGES
                CALL tl_set_avg (ng, tile)
#  endif
              END DO
!$OMP BARRIER
            END DO
 
# ifdef NESTING
!
!  If composite or mosaic grids, process additional points in the
!  contact zone between connected grids for 3D kernel free-surface.
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (any(compositegrid(:,ng))) THEN
                CALL tl_nesting (ng, itlm, nzeta)
              END IF
            END DO
# endif
!
!-----------------------------------------------------------------------
!  If appropriate, write out fields into output NetCDF files.  Notice
!  that IO data is written in delayed and serial mode.  Exit if last
!  time step.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
!$OMP MASTER
              CALL tl_output (ng)
!$OMP END MASTER
!$OMP BARRIER
              IF ((founderror(exit_flag, noerror, __line__, myfile)).or.&
     &            ((iic(ng).eq.(ntend(ng)+1)).and.(ng.eq.ngrids))) THEN
                RETURN
              END IF
            END DO
 
# ifdef NESTING
!
!-----------------------------------------------------------------------
!  If refinement grid, interpolate (space, time) state variables
!  contact points from donor grid extracted data.
#  ifdef NESTING_DEBUG
!
!  Also, fill BRY_CONTACT(:,:)%Mflux to check for mass conservation
!  between coarse and fine grids.  This is only done for diagnostic
!  purposes. Also, debugging is possible with very verbose output
!  to fort.300 is allowed by activating uppercase(nesting_debug).
#  endif
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (refinedgrid(ng).and.(refinescale(ng).gt.0)) THEN
                CALL tl_nesting (ng, itlm, nputd)
#  ifdef NESTING_DEBUG
                CALL tl_nesting (ng, itlm, nmflx)
#  endif
              END IF
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Compute right-hand-side terms for 3D equations.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=last_tile(ng),first_tile(ng),-1
                CALL tl_rhs3d (ng, tile)
#  ifdef MY25_MIXING_NOT_YET
                CALL tl_my25_prestep (ng, tile)
#  elif defined GLS_MIXING_NOT_YET
                CALL tl_gls_prestep (ng, tile)
#  endif
              END DO
!$OMP BARRIER
            END DO
 
# ifdef NESTING
!
!  If composite or mosaic grids, process additional points in the
!  contact zone between connected grids for right-hand-side terms
!  (tracers).
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (any(compositegrid(:,ng))) THEN
                CALL tl_nesting (ng, itlm, nrhst)
              END IF
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Solve the vertically integrated primitive equations for the
!  free-surface and barotropic momentum components.
!-----------------------------------------------------------------------
!
            loop_2d : DO my_iif=1,maxval(nfast)+1
!
!  Set time indices for predictor step. The PREDICTOR_2D_STEP switch
!  it is assumed to be false before the first time-step.
!
              DO ig=1,gridsinlayer(nl)
                ng=gridnumber(ig,nl)
                next_indx1=3-indx1(ng)
                IF (.not.predictor_2d_step(ng).and.                     &
     &              my_iif.le.(nfast(ng)+1)) THEN
                  predictor_2d_step(ng)=.true.
                  iif(ng)=my_iif
                  IF (first_2d_step) THEN
                    kstp(ng)=indx1(ng)
                  ELSE
                    kstp(ng)=3-indx1(ng)
                  END IF
                  knew(ng)=3
                  krhs(ng)=indx1(ng)
                END IF
!
!  Predictor step - Advance barotropic equations using 2D time-step
!  ==============   predictor scheme.  No actual time-stepping is
!  performed during the auxiliary (nfast+1) time-step. It is needed
!  to finalize the fast-time averaging of 2D fields, if any, and
!  compute the new time-evolving depths.
!
                IF (my_iif.le.(nfast(ng)+1)) THEN
                  DO tile=last_tile(ng),first_tile(ng),-1
                    CALL tl_step2d (ng, tile)
                  END DO
!$OMP BARRIER
                END IF
              END DO
 
# ifdef NESTING
!
!  If composite or mosaic grids, process additional points in the
!  contact zone between connected grids for the state variables
!  associated with the 2D engine Predictor Step section.
!  If refinement, check mass flux conservation between coarse and
!  fine grids during debugging.
#  ifdef NESTING_DEBUG
!  Warning: very verbose output to fort.300 ascii file to check
!           mass flux conservation.
#  endif
!
              DO ig=1,gridsinlayer(nl)
                ng=gridnumber(ig,nl)
                IF (any(compositegrid(:,ng))) THEN
                  CALL tl_nesting (ng, itlm, n2dps)
                END IF
                IF (refinedgrid(ng).and.(refinescale(ng).gt.0)) THEN
                  CALL tl_nesting (ng, itlm, nmflx)
                END IF
              END DO
# endif
!
!  Set time indices for corrector step.
!
              DO ig=1,gridsinlayer(nl)
                ng=gridnumber(ig,nl)
                IF (predictor_2d_step(ng)) THEN
                  predictor_2d_step(ng)=.false.
                  knew(ng)=next_indx1
                  kstp(ng)=3-knew(ng)
                  krhs(ng)=3
                  IF (iif(ng).lt.(nfast(ng)+1)) indx1(ng)=next_indx1
                END IF
!
!  Corrector step - Apply 2D time-step corrector scheme.  Notice that
!  ==============   there is not need for a corrector step during the
!  auxiliary (nfast+1) time-step.
!
                IF (iif(ng).lt.(nfast(ng)+1)) THEN
                  DO tile=first_tile(ng),last_tile(ng),+1
                    CALL tl_step2d (ng, tile)
                  END DO
!$OMP BARRIER
                END IF
              END DO
 
# ifdef NESTING
!
!  If composite or mosaic grids, process additional points in the
!  contact zone between connected grids for the state variables
!  associated with the 2D engine Corrector step section.
!  If refinement, check mass flux conservation between coarse and
!  fine grids during debugging.
#  ifdef NESTING_DEBUG
!  Warning: very verbose output to fort.300 ascii file to check
!           mass flux conservation.
#  endif
!
              DO ig=1,gridsinlayer(nl)
                ng=gridnumber(ig,nl)
                IF (any(compositegrid(:,ng))) THEN
                  CALL tl_nesting (ng, itlm, n2dcs)
                END IF
                IF (refinedgrid(ng).and.(refinescale(ng).gt.0)) THEN
                  CALL tl_nesting (ng, itlm, nmflx)
                END IF
              END DO
# endif
            END DO loop_2d
 
# ifdef NESTING
#  if defined MASKING && defined WET_DRY
!
!  If nesting and wetting and drying, scale horizontal interpolation
!  weights to account for land/sea masking in contact areas. This needs
!  to be done at very time-step since the Land/Sea masking is time
!  dependent.
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              CALL tl_nesting (ng, itlm, nmask)
            END DO
#  endif
!
!  If composite or mosaic grids, process additional points in the
!  contact zone between connected grids for the time-averaged
!  momentum fluxes (DU_avg1, DV_avg1) and free-surface (Zt_avg).
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (any(compositegrid(:,ng))) THEN
                CALL tl_nesting (ng, itlm, n2dfx)
              END IF
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Recompute depths and thicknesses using the new time filtered
!  free-surface.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=last_tile(ng),first_tile(ng),-1
                CALL tl_set_depth (ng, tile, itlm)
              END DO
!$OMP BARRIER
            END DO
 
# ifdef NESTING
!
!  If nesting, determine vertical indices and vertical interpolation
!  weights in the contact zone using new depth arrays.
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              CALL tl_nesting (ng, itlm, nzwgt)
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Time-step 3D momentum equations.
!-----------------------------------------------------------------------
!
!  Time-step 3D momentum equations and couple with vertically
!  integrated equations.
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=last_tile(ng),first_tile(ng),-1
                CALL tl_step3d_uv (ng, tile)
              END DO
!$OMP BARRIER
            END DO
 
# ifdef NESTING
!
!  If composite or mosaic grids, process additional points in the
!  contact zone between connected grids for 3D momentum (u,v),
!  adjusted 2D momentum (ubar,vbar), and fluxes (Huon, Hvom).
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (any(compositegrid(:,ng))) THEN
                CALL tl_nesting (ng, itlm, n3duv)
              END IF
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Time-step vertical mixing turbulent equations and passive tracer
!  source and sink terms, if applicable.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=first_tile(ng),last_tile(ng),+1
                CALL tl_omega (ng, tile, itlm)
# ifdef MY25_MIXING_NOT_YET
                CALL tl_my25_corstep (ng, tile)
# elif defined GLS_MIXING_NOT_YET
                CALL tl_gls_corstep (ng, tile)
# endif
# ifdef BIOLOGY
                CALL tl_biology (ng, tile)
# endif
# ifdef SEDIMENT_NOT_YET
                CALL tl_sediment (ng, tile)
# endif
              END DO
!$OMP BARRIER
            END DO
 
# ifndef TS_FIXED
!
!-----------------------------------------------------------------------
!  Time-step tracer equations.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=last_tile(ng),first_tile(ng),-1
                CALL tl_step3d_t (ng, tile)
              END DO
!$OMP BARRIER
            END DO
 
#  ifdef NESTING
!
!  If composite or mosaic grids, process additional points in the
!  contact zone between connected grids for Tracer Variables.
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (any(compositegrid(:,ng))) THEN
                CALL tl_nesting (ng, itlm, n3dtv)
              END IF
            END DO
#  endif
# endif
 
# ifdef NESTING
#  ifndef ONE_WAY
!
!-----------------------------------------------------------------------
!  If refinement grids, perform two-way coupling between fine and
!  coarse grids. Correct coarse grid tracers values at the refinement
!  grid with refined accumulated fluxes.  Then, replace coarse grid
!  state variable with averaged refined grid values (two-way nesting).
!  Update coarse grid depth variables.
!
!  The two-way exchange of infomation between nested grids needs to be
!  done at the correct time-step and in the right sequence.
!-----------------------------------------------------------------------
!
            DO il=nestlayers,1,-1
              DO ig=1,gridsinlayer(il)
                ng=gridnumber(ig,il)
                IF (do_twoway(itlm, nl, il, ng, istep)) THEN
                  CALL tl_nesting (ng, itlm, n2way)
                END IF
              END DO
            END DO
#  endif
!
!-----------------------------------------------------------------------
!  If donor to a finer grid, extract data for the external contact
!  points. This is the latest solution for the coarser grid.
!
!  It is stored in the REFINED structure so it can be used for the
!  space-time interpolation when "nputD" argument is used above.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (donortofiner(ng)) THEN
                CALL tl_nesting (ng, itlm, ngetd)
              END IF
            END DO
# endif
 
# ifdef FLOATS_NOT_YET
!
!-----------------------------------------------------------------------
!  Compute Lagrangian drifters trajectories: Split all the drifters
!  between all the computational threads, except in distributed-memory
!  and serial configurations. In distributed-memory, the parallel node
!  containing the drifter is selected internally since the state
!  variables do not have a global scope.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (lfloats(ng)) THEN
#  ifdef _OPENMP
                chunk_size=(nfloats(ng)+numthreads-1)/numthreads
                lstr=1+mythread*chunk_size
                lend=min(nfloats(ng),lstr+chunk_size-1)
#  else
                lstr=1
                lend=nfloats(ng)
#  endif
                CALL tl_step_floats (ng, lstr, lend)
!$OMP BARRIER
!
!  Shift floats time indices.
!
                nfp1(ng)=mod(nfp1(ng)+1,nft+1)
                nf(ng)=mod(nf(ng)+1,nft+1)
                nfm1(ng)=mod(nfm1(ng)+1,nft+1)
                nfm2(ng)=mod(nfm2(ng)+1,nft+1)
                nfm3(ng)=mod(nfm3(ng)+1,nft+1)
              END IF
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Advance time index and time clock.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              iic(ng)=iic(ng)+1
              time(ng)=time(ng)+dt(ng)
              step_counter(ng)=step_counter(ng)-1
              CALL time_string (time(ng), time_code(ng))
            END DO
 
          END DO step_loop
 
        END DO nest_layer
 
      END DO kernel_loop
 
      RETURN
      END SUBROUTINE tl_main3d
#else
      SUBROUTINE tl_main3d
      RETURN
      END SUBROUTINE tl_main3d
#endif