ad_main3d.F

Go to the documentation of this file.

#include "cppdefs.h"
#if defined ADJOINT && defined SOLVE3D
      SUBROUTINE ad_main3d (RunInterval)
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group            Andrew M. Moore   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This subroutine is the main driver for adjoint  ROMS when           !
!  configurated as a full 3D baroclinic  ocean model. It advances      !
!  backwards the adjoint 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
# ifdef FOUR_DVAR
      USE mod_fourdvar
# endif
      USE mod_iounits
# ifdef SENSITIVITY_4DVAR
      USE mod_ncparam
# endif
# ifdef NESTING
      USE mod_nesting
# endif
      USE mod_scalars
      USE mod_stepping
# ifdef SO_SEMI
      USE mod_storage
# endif
!
# if defined AD_SENSITIVITY   || defined I4DVAR_ANA_SENSITIVITY || \
     defined opt_observations || defined sensitivity_4dvar
      USE adsen_force_mod,         ONLY : adsen_force
# endif
# ifdef ANA_VMIX
      USE analytical_mod,          ONLY : ana_vmix
# endif
# ifdef BIOLOGY
      USE ad_biology_mod,          ONLY : ad_biology
# endif
# ifdef BBL_MODEL_NOT_YET
!!    USE ad_bbl_mod,              ONLY : ad_bblm
# endif
# if defined BULK_FLUXES_NOT_YET && !defined PRIOR_BULK_FLUXES
!!    USE ad_bulk_flux_mod,        ONLY : ad_bulk_flux
# endif
# ifdef BVF_MIXING_NOT_YET
!!    USE ad_bvf_mix_mod,          ONLY : ad_bvf_mix
# endif
      USE ad_diag_mod, ONLY : ad_diag
# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
      USE ad_frc_adjust_mod,       ONLY : ad_frc_adjust
# endif
# ifdef AD_OUTPUT_STATE
      USE ad_ini_fields_mod,       ONLY : ad_out_fields, ad_out_zeta
# endif
# ifdef WEAK_CONSTRAINT
      USE ad_force_dual_mod,       ONLY : ad_force_dual
# endif
# ifdef FORCING_SV
      USE ad_forcing_mod,          ONLY : ad_forcing
# endif
# ifdef GLS_MIXING_NOT_YET
!!    USE ad_gls_corstep_mod,      ONLY : ad_gls_corstep
!!    USE ad_gls_prestep_mod,      ONLY : ad_gls_prestep
# endif
# ifdef LMD_MIXING_NOT_YET
!!    USE ad_lmd_vmix_mod,         ONLY : ad_lmd_vmix
# endif
# if defined FOUR_DVAR && defined OBSERVATIONS
#  ifdef WEAK_CONSTRAINT
      USE ad_htobs_mod,            ONLY : ad_htobs
#  else
      USE ad_misfit_mod,           ONLY : ad_misfit
#  endif
# endif
# ifdef MY25_MIXING
!!    USE ad_my25_corstep_mod,     ONLY : ad_my25_corstep
!!    USE ad_my25_prestep_mod,     ONLY : ad_my25_prestep
# endif
# ifdef NESTING
      USE nesting_mod,             ONLY : nesting
      USE ad_nesting_mod,          ONLY : ad_nesting
#  ifndef ONE_WAY
      USE nesting_mod,             ONLY : do_twoway
#  endif
# endif
# ifdef ADJUST_BOUNDARY
      USE ad_obc_adjust_mod,       ONLY : ad_obc_adjust
      USE ad_obc_adjust_mod,       ONLY : ad_obc2d_adjust
      USE ad_set_depth_mod,        ONLY : ad_set_depth_bry
# endif
      USE ad_omega_mod,            ONLY : ad_omega
# ifndef FORCING_SV
      USE ad_post_initial_mod,     ONLY : ad_post_initial
# endif
# ifdef NEARSHORE_MELLOR_NOT_YET
!!    USE ad_radiation_stress_mod, ONLY : ad_radiation_stress
# endif
# ifndef TS_FIXED
      USE ad_rho_eos_mod,          ONLY : ad_rho_eos
# endif
      USE ad_rhs3d_mod,            ONLY : ad_rhs3d
# ifdef SEDIMENT_NOT_YET
!!    USE ad_sediment_mod,         ONLY : ad_sediment
# endif
# ifdef AD_AVERAGES
      USE ad_set_avg_mod,          ONLY : ad_set_avg
# endif
      USE ad_set_depth_mod,        ONLY : ad_set_depth
      USE ad_set_massflux_mod,     ONLY : ad_set_massflux
# if defined SSH_TIDES_NOT_YET || defined UV_TIDES_NOT_YET
!!    USE ad_set_tides_mod,        ONLY : ad_set_tides
# endif
      USE ad_set_vbc_mod,          ONLY : ad_set_vbc
      USE ad_set_zeta_mod,         ONLY : ad_set_zeta
      USE ad_step2d_mod,           ONLY : ad_step2d
# ifndef TS_FIXED
      USE ad_step3d_t_mod,         ONLY : ad_step3d_t
# endif
      USE ad_step3d_uv_mod,        ONLY : ad_step3d_uv
# ifdef FLOATS_NOT_YET
!!    USE ad_step_floats_mod,      ONLY : ad_step_floats
# endif
# if defined BULK_FLUXES && !defined PRIOR_BULK_FLUXES
      USE bulk_flux_mod,           ONLY : bulk_flux
# endif
      USE dateclock_mod,           ONLY : time_string
# ifdef TIDE_GENERATING_FORCES
      USE equilibrium_tide_mod,    ONLY : equilibrium_tide
# endif
# ifdef WEAK_CONSTRAINT
      USE frc_weak_mod,            ONLY : frc_adgather, frc_clear
# 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 FOUR_DVAR    && !defined I4DVAR_ANA_SENSITIVITY) && \
      defined observations
      USE obs_read_mod,            ONLY : obs_read
# endif
      USE omega_mod,               ONLY : omega
# ifdef SO_SEMI
      USE packing_mod,             ONLY : ad_pack
# endif
      USE rho_eos_mod,             ONLY : rho_eos
      USE set_depth_mod,           ONLY : set_depth
      USE set_massflux_mod
      USE strings_mod,             ONLY : founderror
# ifdef SENSITIVITY_4DVAR
      USE ad_wvelocity_mod,        ONLY : ad_wvelocity
# endif
!
      implicit none
!
!  Imported variable declarations.
!
      real(dp), intent(in) :: RunInterval
!
!  Local variable declarations.
!
      logical :: backward = .true.
      logical :: ad_advance
      logical :: DoNestLayer, Time_Step
      integer :: Nsteps, Rsteps
      integer :: ig, il, is, istep, ng, nl, tile
      integer :: my_iif, next_indx1
# ifdef FLOATS_NOT_YET
      integer :: Lend, Lstr, chunk_size
# endif
      integer :: ks, kt
# ifdef NESTING
      integer :: icount, itcount
# endif
!
      real(r8) :: HalfDT, MaxDT, my_StepTime
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__
!
!=======================================================================
!  Time-step adjoint 3D primitive equations backwards.
!=======================================================================
 
# ifdef NESTING
!
!  Compute nonlinear model timestepping sequence and load it into
!  StepInfo array.  The sequence of values are reversed to timestep
!  the adjoint model backwards.
!
      CALL nlm_step_sequence (runinterval, icount)
      IF (founderror(exit_flag, noerror, __line__,                      &
     &               __file__)) RETURN
# endif
!
!  Initialize.
!
      time_step=.true.
      donestlayer=.true.
# ifdef NESTING
      itcount=icount+1
# endif
!
      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
        DO ng=1,ngrids
          twowayinterval(ng)=dt(1)
        END DO
# endif
!
        nest_layer : DO WHILE (donestlayer)
 
# ifdef NESTING
!
!  Determine current nested layer (nl), number of timesteps (Nsteps)
!  to execute for grids in current nested layer, number of timesteps
!  (Rsteps) to complete the RunInterval time window, and timestep
!  counter (step_counter) for each grid. Their values are reversed
!  from the saved nonlinear model time stepping sequence.
!
          itcount=itcount-1             ! count backwards for adjoint
          IF (itcount.gt.0) THEN
            nl=stepinfo(itcount,1)
            nsteps=stepinfo(itcount,2)
            IF (itcount.eq.icount) nsteps=1
            rsteps=stepinfo(itcount,3)
            DO ng=1,ngrids
              step_counter(ng)=stepinfo(itcount,ng+3)
            END DO
          ELSE
            nl=0
          END IF
# else
!
!  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 (iadm, runinterval, nl, nsteps, rsteps)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
!
          IF ((nl.le.0).or.(nl.gt.nestlayers)) EXIT
!
!  Time-step governing equations for Nsteps.
!
          step_loop : DO istep=nsteps,1,-1
!
!  Set time indices and time clock.
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
# ifdef JEDI
              jic(ng)=jic(ng)-1
              time4jedi(ng)=time4jedi(ng)-dt(ng)
# endif
              tdays(ng)=time(ng)*sec2day
              IF (step_counter(ng).eq.1) time_step=.false.
              IF ((iic(ng).eq.ntstart(ng)).and.(ng.eq.ngrids)) THEN
                ad_advance=.false.
# ifdef NESTING
                DO ng=1,ngrids
                  twowayinterval(ng)=0.0_r8
                END DO
# endif
              ELSE
                ad_advance=.true.
              END IF
            END DO
!
!-----------------------------------------------------------------------
!  Read in required data, if any, from input NetCDF files.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              CALL ad_get_data (ng)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
!
!-----------------------------------------------------------------------
!  Process input data, if any: time interpolate between snapshots.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=first_tile(ng),last_tile(ng),+1
                CALL ad_set_data (ng, tile)
# if defined FORWARD_READ || defined JEDI
                CALL set_depth (ng, tile, iadm)
# endif
# ifdef TIDE_GENERATING_FORCES
                CALL equilibrium_tide (ng, tile, iadm)
# endif
              END DO
            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, iadm)
                CALL rho_eos (ng, tile, iadm)
#  if defined BULK_FLUXES && !defined PRIOR_BULK_FLUXES
!!              CALL bulk_flux (ng, tile)
#  endif
              END DO
            END DO
# endif
!
!  Avoid time-stepping if additional delayed IO time-step.
!
            advance : IF (ad_advance) THEN
 
# 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.
!-----------------------------------------------------------------------
!
!  Shift floats time indices.
!
              DO ig=1,gridsinlayer(nl)
                ng=gridnumber(ig,nl)
                IF (lfloats(ng)) THEN
                  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)
!
#  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 ad_step_floats (ng, lstr, lend)
                END IF
              END DO
# endif
 
# ifdef NESTING
!
!-----------------------------------------------------------------------
!  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 ad_nesting (ng, iadm, ngetd)
                END IF
              END DO
#  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=1,nestlayers
                DO ig=1,gridsinlayer(il)
                  ng=gridnumber(ig,il)
                  IF (do_twoway(iadm,nl,il,ng,istep)) THEN
                    CALL ad_nesting (ng, iadm, n2way)
                  END IF
                END DO
              END DO
#  endif
# endif
 
# ifndef TS_FIXED
!
!-----------------------------------------------------------------------
!  Time-step adjoint tracer equations.
!-----------------------------------------------------------------------
!
#  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 ad_nesting (ng, iadm, n3dtv)
                END IF
              END DO
#  endif
!  Compute intermediate BASIC STATE mass fluxes (Huon,Hvom) for use in
!  the adjoint horizontal advection (ad_step3d_t) and adjoint vertical
!  velocity (ad_omega).
!
              DO ig=1,gridsinlayer(nl)
                ng=gridnumber(ig,nl)
                DO tile=first_tile(ng),last_tile(ng),+1
                  CALL reset_massflux (ng, tile, iadm)    ! intermediate
                END DO                                    ! fluxes
              END DO
!
!  Compute basic STATE omega vertical velocity with intermediate mass
!  fluxes. Time-step adjoint tracer equations.
!
              DO ig=1,gridsinlayer(nl)
                ng=gridnumber(ig,nl)
                DO tile=last_tile(ng),first_tile(ng),-1
                  CALL omega (ng, tile, iadm)   ! BASIC STATE w-velocity
                  CALL ad_step3d_t (ng, tile)
                END DO
              END DO
# endif
!
!-----------------------------------------------------------------------
!  Time-step adjoint vertical mixing turbulent equations and passive
!  tracer source and sink terms, if applicable. Reinstate original
!  BASIC STATE (Huon,Hvom).
!-----------------------------------------------------------------------
!
              DO ig=1,gridsinlayer(nl)
                ng=gridnumber(ig,nl)
                DO tile=first_tile(ng),last_tile(ng),+1
# ifdef SEDIMENT_NOT_YET
                  CALL ad_sediment (ng, tile)
# endif
# ifdef BIOLOGY
                  CALL ad_biology (ng, tile)
# endif
# ifdef MY25_MIXING_NOT_YET
                  CALL ad_my25_corstep (ng, tile)
# elif defined GLS_MIXING_NOT_YET
                  CALL ad_gls_corstep (ng, tile)
# endif
                  CALL ad_omega (ng, tile, iadm)
                  CALL set_massflux (ng, tile, iadm)     ! BASIC STATE
                END DO                                   ! mass fluxes
              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 ad_nesting (ng, iadm, n3duv)
                END IF
              END DO
# endif
!
!-----------------------------------------------------------------------
!  Time-step adjoint 3D equations.
!-----------------------------------------------------------------------
!
!  Reinstate original BASIC STATE omega vertical velocity. Time-step
!  3D adjoint 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 omega (ng, tile, iadm)   ! BASIC STATE w-velocity
                  CALL ad_step3d_uv (ng, tile)
                END DO
              END DO
!
!-----------------------------------------------------------------------
!  Adjoint of recompute depths and thicknesses using the new time
!  filtered free-surface. This call was moved from "ad_step2d" to here.
!-----------------------------------------------------------------------
!
# 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 ad_nesting (ng, iadm, nzwgt)
              END DO
# endif
              DO ig=1,gridsinlayer(nl)
                ng=gridnumber(ig,nl)
                DO tile=last_tile(ng),first_tile(ng),-1
                  CALL ad_set_depth (ng, tile, iadm)
                END DO
              END DO
 
# ifdef NESTING
!
!  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 ad_nesting (ng, iadm, n2dfx)
                END IF
              END DO
 
#  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 ad_nesting (ng, iadm, nmask)
              END DO
#  endif
# endif
!
!-----------------------------------------------------------------------
!  Solve adjoint vertically integrated primitive equations for the
!  free-surface and barotropic momentum components.
!-----------------------------------------------------------------------
!
              loop_2d : DO my_iif=maxval(nfast)+1,1,-1
 
# 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.
!
                DO ig=1,gridsinlayer(nl)
                  ng=gridnumber(ig,nl)
                  IF (any(compositegrid(:,ng))) THEN
                    CALL ad_nesting (ng, iadm, n2dcs)
                  END IF
                  IF (refinedgrid(ng).and.(refinescale(ng).gt.0)) THEN
                    CALL ad_nesting (ng, iadm, nmflx)
                  END IF
                END DO
# endif
!
!  Corrector step - Apply 2D adjoint time-step corrector scheme.  Notice
!  ==============    that there is not need for a corrector step during
!  the auxiliary (nfast+1) time-step.
!
                DO ig=1,gridsinlayer(nl)
                  ng=gridnumber(ig,nl)
                  iif(ng)=my_iif
                  IF (iif(ng).lt.(nfast(ng)+1)) THEN
                    DO tile=first_tile(ng),last_tile(ng),+1
                      CALL ad_step2d (ng, tile)
                    END DO
                  END IF
!
!  Set time indices for adjoint predictor step.
!
                  next_indx1=3-indx1(ng)
                  IF (.not.predictor_2d_step(ng)) THEN
                    predictor_2d_step(ng)=.true.
!^                  knew(ng)=next_indx1
!^                  kstp(ng)=3-knew(ng)
!^                  krhs(ng)=3
!^
                    kt=knew(ng)
                    ks=kstp(ng)
                    knew(ng)=krhs(ng)
                    kstp(ng)=kt
                    krhs(ng)=ks
!^                  IF (my_iif.lt.(nfast(ng)+1)) indx1(ng)=next_indx1
                  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.
!
                DO ig=1,gridsinlayer(nl)
                  ng=gridnumber(ig,nl)
                  IF (any(compositegrid(:,ng))) THEN
                    CALL ad_nesting (ng, iadm, n2dps)
                  END IF
                  IF (refinedgrid(ng).and.(refinescale(ng).gt.0)) THEN
                    CALL ad_nesting (ng, iadm, nmflx)
                  END IF
                END DO
# endif
!
!  Predictor step - Advance adjoint 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.
!
                DO ig=1,gridsinlayer(nl)
                  ng=gridnumber(ig,nl)
                  IF (my_iif.le.(nfast(ng)+1)) THEN
                    DO tile=last_tile(ng),first_tile(ng),-1
                      CALL ad_step2d (ng, tile)
                    END DO
                  END IF
!
!  Set time indices for next adjoint corrector step. The
!  PREDICTOR_2D_STEP switch it is assumed to be false before the
!  first time-step.
!
                  IF (predictor_2d_step(ng).and.                        &
     &                my_iif.le.(nfast(ng)+1)) THEN
                    predictor_2d_step(ng)=.false.
!^                  IF (FIRST_2D_STEP) THEN
!^                    kstp(ng)=indx1(ng)
!^                  ELSE
!^                    kstp(ng)=3-indx1(ng)
!^                  END IF
!^                  knew(ng)=3
!^                  krhs(ng)=indx1(ng)
!^
                    ks=knew(ng)
                    knew(ng)=krhs(ng)
                    krhs(ng)=ks
                  END IF
                END DO
 
              END DO loop_2d
 
            END IF advance
 
# if (defined FOUR_DVAR    && !defined I4DVAR_ANA_SENSITIVITY) && \
      defined observations
!
!-----------------------------------------------------------------------
!  If appropriate, read observation and model state at observation
!  locations.  Then, compute adjoint forcing terms due to observations.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
#  ifdef SENSITIVITY_4DVAR
#   ifdef RBL4DVAR_FCT_SENSITIVITY
              IF (.not.lsen4dvar(ng).and.lsenfct(ng)) THEN
#   else
              IF (.not.lsen4dvar(ng)) THEN
#   endif
#  endif
              halfdt=0.5_r8*dt(ng)
              IF (((time(ng)-halfdt).le.obstime(ng)).and.               &
     &             (obstime(ng).lt.(time(ng)+halfdt))) THEN
                processobs(ng)=.true.
                CALL obs_read (ng, iadm, backward)
              ELSE
                processobs(ng)=.false.
              END IF
!
#  ifdef SP4DVAR
!
!  Skip assimilation of obs on first timestep unless inter=Nsaddle.
!
              IF ((iic(ng).ne.ntstart(ng)).or.lsadd(ng)) THEN
#  endif
 
              DO tile=first_tile(ng),last_tile(ng),+1
#  ifdef WEAK_CONSTRAINT
                CALL ad_htobs (ng, tile, iadm)
#  else
                CALL ad_misfit (ng, tile, iadm)
#  endif
              END DO
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
#  ifdef SENSITIVITY_4DVAR
              END IF
#  endif
#  ifdef SP4DVAR
              END IF
#  endif
            END DO
# endif
 
# ifdef WEAK_CONSTRAINT
!
!-----------------------------------------------------------------------
!  If appropriate, add representer coefficients (Beta hat) impulse
!  forcing to adjoint solution. Read next impulse record, if available.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (processobs(ng)) THEN
                DO tile=first_tile(ng),last_tile(ng),+1
                  CALL ad_force_dual (ng, tile, kstp(ng), nstp(ng))
                END DO
              END IF
            END DO
# endif
 
# 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 ad_nesting (ng, iadm, nrhst)
              END IF
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Compute adjoint right-hand-side terms for 3D equations. If
!  applicable, time-step turbulent mixing schemes.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=first_tile(ng),last_tile(ng),+1
# ifdef MY25_MIXING_NOT_YET
                CALL ad_my25_prestep (ng, tile)
# elif defined GLS_MIXING_NOT_YET
                CALL ad_gls_prestep (ng, tile)
# endif
                CALL ad_rhs3d (ng, tile)
              END DO
            END DO
 
# ifdef NESTING
!
!-----------------------------------------------------------------------
!  If refinement grid, interpolate (space, time) state variables
!  contact points from donor grid extracted data.
!
!  Also, fill BRY_CONTACT(:,:)%Mflux to check for mass conservation
!  between coarse and fine grids.  This is only done for diagnostic
!  purposes.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (refinedgrid(ng).and.(refinescale(ng).gt.0)) THEN
                CALL ad_nesting (ng, iadm, nmflx)
                CALL ad_nesting (ng, iadm, nputd)
              END IF
            END DO
# endif
!!          END IF ADVANCE
 
# 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 ad_nesting (ng, iadm, nzeta)
              END IF
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Set adjoint free-surface to it time-averaged value.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=first_tile(ng),last_tile(ng),+1
# ifdef DIAGNOSTICS
!!              CALL set_diags (ng, tile)
# endif
                CALL ad_set_zeta (ng, tile)
              END DO
            END DO
!
!-----------------------------------------------------------------------
!  Compute adjoint vertical mixing coefficients for momentum and
!  tracers. Compute adjoint 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
# ifdef SENSITIVITY_4DVAR
                IF (scalars(ng)%Lstate(iswvel)) THEN
                  CALL ad_wvelocity (ng, tile, nstp(ng))
                END IF
# endif
                CALL ad_omega (ng, tile, iadm)
# if defined ANA_VMIX_NOT_YET
                CALL ad_ana_vmix (ng, tile, iadm)
# elif defined LMD_MIXING_NOT_YET
                CALL ad_lmd_vmix (ng, tile)
# elif defined BVF_MIXING
                CALL ad_bvf_mix (ng, tile)
# endif
              END DO
            END DO
 
# ifdef SO_SEMI
!
!-----------------------------------------------------------------------
!  If stochastic optimals with respect the seminorm of chosen
!  functional, pack adjoint state surface forcing needed by the
!  dynamical propagator.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (mod(iic(ng)-1,nadj(ng)).eq.0) THEN
                sorec(ng)=sorec(ng)+1
                DO tile=first_tile(ng),last_tile(ng),+1
                  CALL ad_pack (ng, tile, nstr(ng), nend(ng),           &
     &                          storage(ng)%so_state(:,sorec(ng)))
                END DO
              END IF
            END DO
# endif
 
# 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 ad_nesting (ng, iadm, nbstr)
              END IF
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Set adjoint 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 SSH_TIDES_NOT_YET || defined UV_TIDES_NOT_YET
                CALL ad_set_tides (ng, tile)
# endif
                CALL ad_set_vbc (ng, tile)
# ifdef BBL_MODEL_NOT_YET
                CALL ad_bblm (ng, tile)
# endif
# if defined BULK_FLUXES_NOT_YET && !defined PRIOR_BULK_FLUXES
                CALL ad_bulk_flux (ng, tile)
# endif
              END DO
            END DO
 
# 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 ad_radiation_stress (ng, tile)
              END DO
            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
              END IF
            END DO
# endif
 
# 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
              END IF
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Compute adjoint density related fields and horizontal mass fluxes
!  (Hz*u/n and Hz*v/m). If applicable, compute and report diagnostics
!  and accumulate time-averaged adjoint fields 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
# ifndef TS_FIXED
                CALL ad_rho_eos (ng, tile, iadm)
# endif
                CALL ad_set_massflux (ng, tile, iadm)
                CALL ad_diag (ng, tile)
# ifdef AD_AVERAGES
                CALL ad_set_avg (ng, tile)
# endif
              END DO
            END DO
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# ifdef FORCING_SV
!
!-----------------------------------------------------------------------
!  Compute the adjoint forcing for the forcing singular vectors.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              DO tile=first_tile(ng),last_tile(ng),+1
                CALL ad_set_depth (ng, tile, iadm)
                CALL ad_forcing (ng, tile, kstp(ng), nstp(ng))
              END DO
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Update 3D time-level indices.
!-----------------------------------------------------------------------
!
!  The original forward time-stepping indices are advanced as follows:
!
!     nstp(ng)=1+MOD(iic(ng)-ntstart(ng),2)
!     nnew(ng)=3-nstp(ng)
!     nrhs(ng)=nstp(ng)
!
!  This yields the only 2 permutations:
!
!     nstp  nnew  nrhs
!      1     2     1
!      2     1     2
!
!   With nstp=1, nnew=1 and nrhs=2 at time zero, this is equivalent to
!   the following:
!
!     nstp(ng)=nnew(ng)
!     nnew(ng)=nrhs(ng)
!     nrhs(ng)=nstp(ng)
!
!   The adjoint of this is as follows:
!
!     nstp(ng)=nrhs(ng)
!     nrhs(ng)=nnew(ng)
!     nnew(ng)=nstp(ng)
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (iic(ng).ne.ntend(ng)) THEN
                nrhs(ng)=nnew(ng)
                nnew(ng)=nstp(ng)
                nstp(ng)=nrhs(ng)
              END IF
            END DO
 
# ifdef AD_OUTPUT_STATE
!
!-----------------------------------------------------------------------
!  Set full adjoint output arrays. Due to the exact discrete adjoint,
!  the predictor/corrector time-stepping scheme with multiple time
!  levels, pieces of the adjoint solution are in two-time levels and
!  need to be added in the "_sol" arrays for output purposes.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (iic(ng).ne.ntend(ng)) THEN
                DO tile=last_tile(ng),first_tile(ng),-1
                  CALL ad_out_fields (ng, tile, iadm)
                END DO
                DO tile=first_tile(ng),last_tile(ng),+1
                  CALL ad_set_depth (ng, tile, iadm)
                  CALL ad_out_zeta (ng, tile, iadm)
                END DO
              END IF
            END DO
# endif
# ifndef FORCING_SV
!
!-----------------------------------------------------------------------
!  On the last timestep, it computes the adjoint of initial depths and
!  level thicknesses from the initial free-surface field. Additionally,
!  it initializes the adjoint state variables for all time levels and
!  and applies lateral boundary conditions.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (iic(ng).eq.ntend(ng)) THEN
                CALL ad_post_initial (ng, iadm)
              END IF
            END DO
# endif
 
# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
!
!-----------------------------------------------------------------------
!  Interpolate surface forcing increments and adjust surface forcing.
!  Skip first timestep.
!-----------------------------------------------------------------------
!
#  ifdef RBL4DVAR_FCT_SENSITIVITY
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (.not.lsen4dvar(ng)) THEN          ! ignore in forecast
                IF (iic(ng).ne.ntstart(ng)) THEN    !           interval
                  DO tile=first_tile(ng),last_tile(ng),+1
                    CALL ad_frc_adjust (ng, tile, lfout(ng))
                  END DO
                END IF
              END IF
            END DO
#  else
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (iic(ng).ne.ntstart(ng)) THEN
                DO tile=first_tile(ng),last_tile(ng),+1
                  CALL ad_frc_adjust (ng, tile, lfout(ng))
                END DO
              END IF
            END DO
#  endif
# endif
 
# ifdef ADJUST_BOUNDARY
!
!-----------------------------------------------------------------------
!  Interpolate open boundary increments and adjust open boundaries.
!  Skip first timestep.
!-----------------------------------------------------------------------
!
#  ifdef RBL4DVAR_FCT_SENSITIVITY
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (.not.lsen4dvar(ng)) THEN          ! ignore in forecast
                IF (iic(ng).ne.ntstart(ng)) THEN    !           interval
                  DO tile=first_tile(ng),last_tile(ng),+1
                    CALL ad_obc2d_adjust (ng, tile, lbout(ng))
                    CALL ad_set_depth_bry (ng, tile, iadm)
                    CALL ad_obc_adjust (ng, tile, lbout(ng))
                  END DO
                END IF
              END IF
            END DO
#  else
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (iic(ng).ne.ntstart(ng)) THEN
                DO tile=first_tile(ng),last_tile(ng),+1
                  CALL ad_obc2d_adjust (ng, tile, lbout(ng))
                  CALL ad_set_depth_bry (ng, tile, iadm)
                  CALL ad_obc_adjust (ng, tile, lbout(ng))
                END DO
              END IF
            END DO
#  endif
# endif
 
# if defined WEAK_CONSTRAINT && !defined SP4DVAR
!
!-----------------------------------------------------------------------
!  Gather weak constraint forcing to storage arrays.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (iic(ng).ne.ntstart(ng)) THEN
                DO tile=first_tile(ng),last_tile(ng),+1
                  CALL ad_set_depth (ng, tile, iadm)
                  CALL frc_adgather (ng, tile)
                END DO
              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.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              CALL ad_output (ng)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
 
# if defined WEAK_CONSTRAINT && !defined SP4DVAR
!
!-----------------------------------------------------------------------
!  Copy storage arrays index 1 into index 2, and clear index 1.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              IF (mod(iic(ng)-1,nadj(ng)).eq.0) THEN
                DO tile=first_tile(ng),last_tile(ng),+1
                  CALL frc_clear (ng, tile)
                END DO
              END IF
            END DO
# endif
 
# if (defined AD_SENSITIVITY   || defined I4DVAR_ANA_SENSITIVITY || \
      defined opt_observations || defined sensitivity_4dvar)     && \
     !defined OBS_SPACE
!
!-----------------------------------------------------------------------
!  Add appropriate forcing terms to the adjoint model. The form of the
!  forcing depends on the functional whose sensitivity is required.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
#  ifdef SENSITIVITY_4DVAR
              IF (lsen4dvar(ng)) THEN
#  endif
#  if !defined AD_IMPULSE
                IF ((dends(ng).ge.tdays(ng)).and.                       &
     &              (tdays(ng).ge.dstrs(ng))) THEN
#  endif
                  DO tile=first_tile(ng),last_tile(ng),+1
                    CALL adsen_force (ng, tile)
                  END DO
#  if !defined AD_IMPULSE
                END IF
#  endif
#  ifdef SENSITIVITY_4DVAR
              END IF
#  endif
            END DO
# endif
!
!-----------------------------------------------------------------------
!  Set adjoint time indices and time clock.
!-----------------------------------------------------------------------
!
            DO ig=1,gridsinlayer(nl)
              ng=gridnumber(ig,nl)
              iic(ng)=iic(ng)-1
# ifndef NESTING
              step_counter(ng)=step_counter(ng)-1
# endif
              time(ng)=time(ng)-dt(ng)
              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 ad_main3d
#else
      SUBROUTINE ad_main3d
      RETURN
      END SUBROUTINE ad_main3d
#endif