obs_sen_r4dvar_analysis.h

Go to the documentation of this file.

      MODULE roms_kernel_mod
!
!git $Id$
!=================================================== Andrew M. Moore ===
!  Copyright (c) 2002-2025 The ROMS Group           Hernan G. Arango   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  ROMS Strong/Weak Constraint 4-Dimensional Variational Data          !
!         Assimilation and Observation Sensitivity Driver: Indirect    !
!         Representer Approach (R4D-Var).                              !
!         Dual formulation in observarion space.                       !
!                                                                      !
!  This driver is used for weak constraint 4D-Var where errors are     !
!  considered in both model and observations. It also computes the     !
!  the sensitivity of the assimilation system to each observation.     !
!  It measures the degree to which each observation contributes to     !
!  the uncertainty in the estimate.  This analysis  can be used to     !
!  determine the type of measurements that need to be made,  where     !
!  to observe, and when.                                               !
!                                                                      !
!  The routines in this driver control the initialization,  time-      !
!  stepping, and finalization of  ROMS  model following ESMF           !
!  conventions:                                                        !
!                                                                      !
!     ROMS_initialize                                                  !
!     ROMS_run                                                         !
!     ROMS_finalize                                                    !
!                                                                      !
!  References:                                                         !
!                                                                      !
!  Moore, A.M., H.G. Arango, G. Broquet, B.S. Powell, A.T. Weaver,     !
!    and J. Zavala-Garay, 2011: The Regional Ocean Modeling System     !
!    (ROMS)  4-dimensional variational data assimilations systems,     !
!    Part I - System overview and formulation, Prog. Oceanogr., 91,    !
!    34-49, doi:10.1016/j.pocean.2011.05.004.                          !
!                                                                      !
!  Moore, A.M., H.G. Arango, G. Broquet, C. Edward, M. Veneziani,      !
!    B. Powell, D. Foley, J.D. Doyle, D. Costa, and P. Robinson,       !
!    2011: The Regional Ocean Modeling System (ROMS) 4-dimensional     !
!    variational data assimilations systems, Part II - Performance     !
!    and application to the California Current System, Prog.           !
!    Oceanogr., 91, 50-73, doi:10.1016/j.pocean.2011.05.003.           !
!                                                                      !
!  Moore, A.M., H.G. Arango, G. Broquet, C. Edward, M. Veneziani,      !
!    B. Powell, D. Foley, J.D. Doyle, D. Costa, and P. Robinson,       !
!    2011: The Regional Ocean Modeling System (ROMS) 4-dimensional     !
!    variational data assimilations systems, Part III - Observation    !
!    impact and observation sensitivity in the California Current      !
!    System, Prog. Oceanogr., 91, 74-94,                               !
!    doi:10.1016/j.pocean.2011.05.005.                                 !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_arrays
      USE mod_fourdvar
      USE mod_iounits
      USE mod_ncparam
      USE mod_netcdf
#if defined PIO_LIB && defined DISTRIBUTE
      USE mod_pio_netcdf
#endif
      USE mod_scalars
      USE mod_stepping
!
#ifdef ADJUST_BOUNDARY
      USE mod_boundary,      ONLY : initialize_boundary
#endif
      USE mod_forces,        ONLY : initialize_forces
      USE mod_ocean,         ONLY : initialize_ocean
!
      USE ad_wrt_his_mod,    ONLY : ad_wrt_his
      USE close_io_mod,      ONLY : close_file, close_inp, close_out
      USE congrad_mod,       ONLY : congrad
      USE convolve_mod,      ONLY : error_covariance
      USE def_impulse_mod,   ONLY : def_impulse
      USE def_mod_mod,       ONLY : def_mod
      USE def_norm_mod,      ONLY : def_norm
      USE get_state_mod,     ONLY : get_state
      USE inp_par_mod,       ONLY : inp_par
#ifdef MCT_LIB
# ifdef ATM_COUPLING
      USE mct_coupler_mod,   ONLY : initialize_ocn2atm_coupling
# endif
# ifdef WAV_COUPLING
      USE mct_coupler_mod,   ONLY : initialize_ocn2wav_coupling
# endif
#endif
      USE normalization_mod, ONLY : normalization
      USE stats_modobs_mod,  ONLY : stats_modobs
      USE stdout_mod,        ONLY : set_stdoutunit, stdout_unit
      USE strings_mod,       ONLY : founderror, uppercase
      USE tl_def_ini_mod,    ONLY : tl_def_ini
      USE wrt_ini_mod,       ONLY : wrt_ini
      USE wrt_rst_mod,       ONLY : wrt_rst
#if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
      USE zeta_balance_mod,  ONLY : balance_ref, biconj
#endif
!
      implicit none
!
      PUBLIC :: roms_initialize
      PUBLIC :: roms_run
      PUBLIC :: roms_finalize
!
      CONTAINS
!
      SUBROUTINE roms_initialize (first, mpiCOMM)
!
!=======================================================================
!                                                                      !
!  This routine allocates and initializes ROMS state variables         !
!  and internal and external parameters.                               !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations.
!
      logical, intent(inout) :: first
!
      integer, intent(in), optional :: mpiCOMM
!
!  Local variable declarations.
!
      logical :: allocate_vars = .true.
!
#ifdef DISTRIBUTE
      integer :: MyError, MySize
#endif
      integer :: STDrec, Tindex
      integer :: chunk_size, ng, thread
#ifdef _OPENMP
      integer :: my_threadnum
#endif
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_initialize"
 
#ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Set distribute-memory (mpi) world communictor.
!-----------------------------------------------------------------------
!
      IF (PRESENT(mpicomm)) THEN
        ocn_comm_world=mpicomm
      ELSE
        ocn_comm_world=mpi_comm_world
      END IF
      CALL mpi_comm_rank (ocn_comm_world, myrank, myerror)
      CALL mpi_comm_size (ocn_comm_world, mysize, myerror)
#endif
!
!-----------------------------------------------------------------------
!  On first pass, initialize model parameters a variables for all
!  nested/composed grids.  Notice that the logical switch "first"
!  is used to allow multiple calls to this routine during ensemble
!  configurations.
!-----------------------------------------------------------------------
!
      IF (first) THEN
        first=.false.
!
!  Initialize parallel control switches. These scalars switches are
!  independent from standard input parameters.
!
        CALL initialize_parallel
!
!  Set the ROMS standard output unit to write verbose execution info.
!  Notice that the default standard out unit in Fortran is 6.
!
!  In some applications like coupling or disjointed mpi-communications,
!  it is advantageous to write standard output to a specific filename
!  instead of the default Fortran standard output unit 6. If that is
!  the case, it opens such formatted file for writing.
!
        IF (set_stdoutunit) THEN
          stdout=stdout_unit(master)
          set_stdoutunit=.false.
        END IF
!
!  Read in model tunable parameters from standard input. Allocate and
!  initialize variables in several modules after the number of nested
!  grids and dimension parameters are known.
!
        CALL inp_par (inlm)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Set domain decomposition tile partition range.  This range is
!  computed only once since the "first_tile" and "last_tile" values
!  are private for each parallel thread/node.
!
#if defined _OPENMP
      mythread=my_threadnum()
#elif defined DISTRIBUTE
      mythread=myrank
#else
      mythread=0
#endif
      DO ng=1,ngrids
        chunk_size=(ntilex(ng)*ntilee(ng)+numthreads-1)/numthreads
        first_tile(ng)=mythread*chunk_size
        last_tile(ng)=first_tile(ng)+chunk_size-1
      END DO
!
!  Initialize internal wall clocks. Notice that the timings does not
!  includes processing standard input because several parameters are
!  needed to allocate clock variables.
!
        IF (master) THEN
          WRITE (stdout,10)
 10       FORMAT (/,' Process Information:',/)
        END IF
!
        DO ng=1,ngrids
          DO thread=thread_range
            CALL wclock_on (ng, inlm, 0, __line__, myfile)
          END DO
        END DO
!
!  Allocate and initialize modules variables.
!
        CALL roms_allocate_arrays (allocate_vars)
        CALL roms_initialize_arrays
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      END IF
 
#if defined MCT_LIB && (defined ATM_COUPLING || defined WAV_COUPLING)
!
!-----------------------------------------------------------------------
!  Initialize coupling streams between model(s).
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
# ifdef ATM_COUPLING
        CALL initialize_ocn2atm_coupling (ng, myrank)
# endif
# ifdef WAV_COUPLING
        CALL initialize_ocn2wav_coupling (ng, myrank)
# endif
      END DO
#endif
 
#if !defined RECOMPUTE_4DVAR
!
!-----------------------------------------------------------------------
!  If the required vectors and arrays from congrad from a previous
!  run of the assimilation cycle are available, read them here from
!  LCZ(ng)%name NetCDF file.
!-----------------------------------------------------------------------
!
      sourcefile=myfile
      DO ng=1,ngrids
        SELECT CASE (lcz(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,               &
     &                            'cg_beta', cg_beta)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,               &
     &                            'cg_delta', cg_delta)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,               &
     &                            'cg_Gnorm_v', cg_gnorm_v)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,               &
     &                            'cg_dla', cg_dla)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,               &
     &                            'cg_QG', cg_qg)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,               &
     &                            'zgrad0', zgrad0)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,               &
     &                            'zcglwk', zcglwk)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL netcdf_get_fvar (ng, itlm, lcz(ng)%name,               &
     &                            'TLmodVal_S', tlmodval_s,             &
     &                            broadcast = .false.) ! Master use only
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,           &
     &                                'cg_beta', cg_beta)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,           &
     &                                'cg_delta', cg_delta)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,           &
     &                                'cg_Gnorm_v', cg_gnorm_v)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,           &
     &                                'cg_dla', cg_dla)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,           &
     &                                'cg_QG', cg_qg)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,           &
     &                                'zgrad0', zgrad0)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,           &
     &                                'zcglwk', zcglwk)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            CALL pio_netcdf_get_fvar (ng, itlm, lcz(ng)%name,           &
     &                                'TLmodVal_S', tlmodval_s)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
        END SELECT
      END DO
#endif
!
!-----------------------------------------------------------------------
!  Read in standard deviation factors for error covariance.
!-----------------------------------------------------------------------
!
!  Initial conditions standard deviation. They are loaded in Tindex=1
!  of the e_var(...,Tindex) state variables.
!
      stdrec=1
      tindex=1
      DO ng=1,ngrids
        CALL get_state (ng, 10, 10, std(1,ng), stdrec, tindex)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Model error standard deviation. They are loaded in Tindex=2
!  of the e_var(...,Tindex) state variables.
!
      stdrec=1
      tindex=2
      IF (nsa.eq.2) THEN
        DO ng=1,ngrids
          CALL get_state (ng, 11, 11, std(2,ng), stdrec, tindex)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END IF
 
#ifdef ADJUST_BOUNDARY
!
!  Open boundary conditions standard deviation.
!
      stdrec=1
      tindex=1
      DO ng=1,ngrids
        CALL get_state (ng, 12, 12, std(3,ng), stdrec, tindex)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
#endif
#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
!
!  Surface forcing standard deviation.
!
      stdrec=1
      tindex=1
      DO ng=1,ngrids
        CALL get_state (ng, 13, 13, std(4,ng), stdrec, tindex)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
#endif
!
      RETURN
      END SUBROUTINE roms_initialize
!
      SUBROUTINE roms_run (RunInterval)
!
!=======================================================================
!                                                                      !
!  This routine time-steps ROMS nonlinear, tangent linear and          !
!  adjoint models.                                                     !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      real(dp), intent(in) :: RunInterval            ! seconds
!
!  Local variable declarations.
!
      logical :: Lcgini, Linner, Lposterior
!
      integer :: my_inner, my_outer
      integer :: Lbck, Lini, Rec, Rec1, Rec2
      integer :: i, lstr, ng, status, tile
      integer :: Fcount, NRMrec
 
      integer, dimension(Ngrids) :: Nrec
!
      real(r8) :: str_day, end_day
!
      character (len=14) :: driver
      character (len=20) :: string
 
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_run"
!
!=======================================================================
!  Run model for all nested grids, if any.
!=======================================================================
!
!  Initialize relevant parameters.
!
      DO ng=1,ngrids
#if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
        lfinp(ng)=1         ! forcing index for input
        lfout(ng)=1         ! forcing index for output history files
#endif
#ifdef ADJUST_BOUNDARY
        lbinp(ng)=1         ! boundary index for input
        lbout(ng)=1         ! boundary index for output history files
#endif
        lold(ng)=1          ! old minimization time index
        lnew(ng)=2          ! new minimization time index
      END DO
      lini=1                ! NLM initial conditions record in INI
      lbck=2                ! background record in INI
      rec1=1
      rec2=2
      nrun=1
      outer=0
      inner=0
      erstr=1
      erend=nouter
      driver='obs_sen_w4dvar'
!
!-----------------------------------------------------------------------
!  Configure weak constraint 4DVAR algorithm: Indirect Representer
!  Approach.
!-----------------------------------------------------------------------
!
!  Initialize the switch to gather weak constraint forcing.
!
      DO ng=1,ngrids
        wrtforce(ng)=.false.
      END DO
!
!  Initialize and set nonlinear model initial conditions.
!
      DO ng=1,ngrids
        wrtnlmod(ng)=.true.
        wrtrpmod(ng)=.false.
        wrttlmod(ng)=.false.
      END DO
!
      CALL initial
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Save nonlinear initial conditions (currently in time index 1,
!  background) into record "Lini" of INI(ng)%name NetCDF file. The
!  record "Lbck" becomes the background state record and the record
!  "Lini" becomes current nonlinear initial conditions.
!
      DO ng=1,ngrids
        ini(ng)%Rindex=1
        fcount=ini(ng)%load
        ini(ng)%Nrec(fcount)=1
#ifdef DISTRIBUTE
        CALL wrt_ini (ng, myrank, 1)
#else
        CALL wrt_ini (ng, -1, 1)
#endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Set nonlinear output history file as the initial basic state
!  trajectory.
!
      DO ng=1,ngrids
        ldefhis(ng)=.true.
        lwrthis(ng)=.true.
        WRITE (his(ng)%name,10) trim(fwd(ng)%head), outer
        lstr=len_trim(his(ng)%name)
        his(ng)%base=his(ng)%name(1:lstr-3)
      END DO
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Model-error covariance normalization and stardard deviation factors.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
!  Compute or read in the error covariance normalization factors.
!  If computing, write out factors to NetCDF. This is an expensive
!  computation that needs to be computed only once for a particular
!  application grid and decorrelation scales.
!
      DO ng=1,ngrids
        IF (any(lwrtnrm(:,ng))) THEN
          CALL def_norm (ng, inlm, 1)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          IF (nsa.eq.2) THEN
            CALL def_norm (ng, inlm, 2)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END IF
 
#ifdef ADJUST_BOUNDARY
          CALL def_norm (ng, inlm, 3)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
 
#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
          CALL def_norm (ng, inlm, 4)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
 
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL normalization (ng, tile, 2)
          END DO
          ldefnrm(1:4,ng)=.false.
          lwrtnrm(1:4,ng)=.false.
        ELSE
          nrmrec=1
          CALL get_state (ng, 14, 14, nrm(1,ng), nrmrec, 1)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          IF (nsa.eq.2) THEN
            CALL get_state (ng, 15, 15, nrm(2,ng), nrmrec, 2)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END IF
 
#ifdef ADJUST_BOUNDARY
          CALL get_state (ng, 16, 16, nrm(3,ng), nrmrec, 1)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
 
#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
          CALL get_state (ng, 17, 17, nrm(4,ng), nrmrec, 1)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#endif
        END IF
      END DO
 
#if !defined RECOMPUTE_4DVAR && defined BALANCE_OPERATOR && \
     defined zeta_elliptic
!
!  Compute the reference zeta and biconjugate gradient arrays
!  required for the balance of free surface.
!
      IF (balance(isfsur)) THEN
        DO ng=1,ngrids
          CALL get_state (ng, inlm, 2, ini(ng), lini, lini)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL balance_ref (ng, tile, lini)
            CALL biconj (ng, tile, inlm, lini)
          END DO
          wrtzetaref(ng)=.true.
        END DO
      END IF
#endif
!
!  Define tangent linear initial conditions file.
!
      DO ng=1,ngrids
        ldefitl(ng)=.true.
        CALL tl_def_ini (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Define TLM/RPM impulse forcing NetCDF file.
!
      DO ng=1,ngrids
        ldeftlf(ng)=.true.
        CALL def_impulse (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Define output 4DVAR NetCDF file containing all processed data
!  at observation locations.
!
      DO ng=1,ngrids
        ldefmod(ng)=.true.
        CALL def_mod (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Write out outer loop beeing processed.
!
      sourcefile=myfile
      DO ng=1,ngrids
        SELECT CASE (dav(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_put_ivar (ng, inlm, dav(ng)%name,               &
     &                            'Nimpact', nimpact,                   &
     &                            (/0/), (/0/),                         &
     &                            ncid = dav(ng)%ncid)
 
#if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL pio_netcdf_put_ivar (ng, inlm, dav(ng)%name,           &
     &                                'Nimpact', nimpact,               &
     &                                (/0/), (/0/),                     &
     &                                piofile = dav(ng)%pioFile)
#endif
        END SELECT
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Run nonlinear model and compute basic state trajectory. It processes
!  and writes the observations accept/reject flag (ObsScale) once to
!  allow background quality control, if any.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
      DO ng=1,ngrids
        wrtobsscale(ng)=.true.
        IF (master) THEN
          WRITE (stdout,20) 'NL', ng, ntstart(ng), ntend(ng)
        END IF
      END DO
!
#ifdef SOLVE3D
      CALL main3d (runinterval)
#else
      CALL main2d (runinterval)
#endif
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      DO ng=1,ngrids
        wrtnlmod(ng)=.false.
        wrtobsscale(ng)=.false.
      END DO
 
#ifdef FORWARD_FLUXES
!
!  Set the BLK structure to contain the nonlinear model surface fluxes
!  needed by the tangent linear and adjoint models. Also, set switches
!  to process that structure in routine "check_multifile". Notice that
!  it is possible to split the solution into multiple NetCDF files to
!  reduce their size.
!
!  The switch LreadFRC is deactivated because all the atmospheric
!  forcing, including shortwave radiation, is read from the NLM
!  surface fluxes or is assigned during ESM coupling.  Such fluxes
!  are available from the QCK structure. There is no need for reading
!  and processing from the FRC structure input forcing-files.
!
      CALL edit_multifile ('QCK2BLK')
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      DO ng=1,ngrids
        lreadblk(ng)=.true.
        lreadfrc(ng)=.false.
        lreadqck(ng)=.false.
      END DO
#endif
 
#ifdef RECOMPUTE_4DVAR
!
!-----------------------------------------------------------------------
!  Solve the system:
!
!              (R_n + Cobs) * Beta_n = h_n
!
!              h_n = Xo - H * X_n
!
!  where R_n is the representer matrix, Cobs is the observation-error
!  covariance, Beta_n are the representer coefficients, h_n is the
!  misfit between observations (Xo) and model (H * X_n), and H is
!  the linearized observation operator. Here, _n denotes a sequence
!  of estimates.
!
!  The system does not need to be solved explicitly by inverting the
!  symmetric stabilized representer matrix, P_n:
!
!              P_n = R_n + Cobs
!
!  but by computing the action of P_n on any vector PSI, such that
!
!              P_n * PSI = R_n * PSI + Cobs * PSI
!
!  The representer matrix is not explicitly computed but evaluated by
!  one integration backward of the adjoint model and one integration
!  forward of the tangent linear model for any forcing vector PSI.
!
!  A preconditioned conjugate gradient algorithm is used to compute
!  an approximation PSI for Beta_n.
!
!-----------------------------------------------------------------------
!
!  If the required vectors and arrays from congrad from a previous run
!  of the assimilation cycle are not available, rerun the 4D-Var cycle.
!
      outer_loop : DO my_outer=1,1
        outer=my_outer
        inner=0
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Run representer model and compute a "prior estimate" state
!  trajectory, X_n(t). Use linearized state trajectory (X_n-1) as
!  basic state.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
!  Set representer model basic state trajectory file to previous outer
!  loop file (outer-1). If outer=1, the basic state trajectory is the
!  nonlinear model.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(fwd(ng)%head), outer-1
          lstr=len_trim(fwd(ng)%name)
          fwd(ng)%base=fwd(ng)%name(1:lstr-3)
        END DO
!
!  Set structure for the nonlinear forward trajectory to be processed
!  by the tangent linear and adjoint models. Also, set switches to
!  process the FWD structure in routine "check_multifile". Notice that
!  it is possible to split solution into multiple NetCDF files to reduce
!  their size.
!
        IF (outer.eq.1) THEN
          CALL edit_multifile ('HIS2FWD')
        ELSE
          CALL edit_multifile ('TLM2FWD')
        END IF
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        DO ng=1,ngrids
          lreadfwd(ng)=.true.
        END DO
!
!  Set representer model output file name.  The strategy is to write
!  the representer solution at the beginning of each outer loop.
!
        DO ng=1,ngrids
          ideftlm(ng)=-1
          ldeftlm(ng)=.true.
          lwrttlm(ng)=.true.
          WRITE (tlm(ng)%name,10) trim(tlm(ng)%head), outer
          lstr=len_trim(tlm(ng)%name)
          tlm(ng)%base=tlm(ng)%name(1:lstr-3)
        END DO
!
!  Activate switch to write the representer model at observation points.
!  Turn off writing into history file and turn off impulse forcing.
!
        DO ng=1,ngrids
          wrtrpmod(ng)=.true.
          sporadicimpulse(ng)=.false.
          frequentimpulse(ng)=.false.
        END DO
 
# ifndef DATALESS_LOOPS
!
!  As in the nonlinear model, initialize always the representer model
!  here with the background or reference state (IRP(ng)%name, record
!  Rec1).
!
        DO ng=1,ngrids
          irp(ng)%Rindex=rec1
          CALL rp_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Run representer model using the nonlinear trajectory as a basic
!  state.  Compute model solution at observation points, H * X_n.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) 'RP', ng, ntstart(ng), ntend(ng)
          END IF
        END DO
!
#  ifdef SOLVE3D
        CALL rp_main3d (runinterval)
#  else
        CALL rp_main2d (runinterval)
#  endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Report data penalty function.
!
        DO ng=1,ngrids
          IF (master) THEN
            DO i=0,nstatevar(ng)
              IF (i.eq.0) THEN
                string='Total'
              ELSE
                string=vname(1,idsvar(i))
              END IF
              IF (fourdvar(ng)%DataPenalty(i).ne.0.0_r8) THEN
                WRITE (stdout,30) outer, inner, 'RPM',                  &
     &                            fourdvar(ng)%DataPenalty(i),          &
     &                            trim(string)
              END IF
            END DO
          END IF
!
!  Write out initial data penalty function to NetCDF file.
!
          sourcefile=myfile
          SELECT CASE (dav(ng)%IOtype)
            CASE (io_nf90)
              CALL netcdf_put_fvar (ng, irpm, dav(ng)%name,             &
     &                              'RP_iDataPenalty',                  &
     &                              fourdvar(ng)%DataPenalty(0:),       &
     &                              (/1,outer/),                        &
     &                              (/nstatevar(ng)+1,1/),              &
     &                              ncid = dav(ng)%ncid)
#  if defined PIO_LIB && defined DISTRIBUTE
            CASE (io_pio)
              CALL pio_netcdf_put_fvar (ng, irpm, dav(ng)%name,         &
     &                                  'RP_iDataPenalty',              &
     &                                  fourdvar(ng)%DataPenalty(0:),   &
     &                                  (/1,outer/),                    &
     &                                  (/nstatevar(ng)+1,1/),          &
     &                                  piofile = dav(ng)%pioFile)
#  endif
          END SELECT
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Clean penalty array before next run of RP model.
!
          fourdvar(ng)%DataPenalty=0.0_r8
        END DO
!
!  Turn off IO switches.
!
        DO ng=1,ngrids
          ldeftlm(ng)=.false.
          lwrttlm(ng)=.false.
          wrtrpmod(ng)=.false.
        END DO
!
!  Clear tangent linear forcing arrays before entering inner-loop.
!  This is very important since these arrays are non-zero after
!  running the representer model and must be zero when running the
!  tangent linear model.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
#  ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
#  endif
          END DO
        END DO
 
#  if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
!
!  Compute the reference zeta and biconjugate gradient arrays
!  required for the balance of free surface.
!
        IF (balance(isfsur)) THEN
          DO ng=1,ngrids
            CALL get_state (ng, inlm, 2, ini(ng), lini, lini)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL balance_ref (ng, tile, lini)
              CALL biconj (ng, tile, inlm, lini)
            END DO
            wrtzetaref(ng)=.true.
          END DO
        END IF
#  endif
!
        inner_loop : DO my_inner=0,ninner
          inner=my_inner
!
! Initialize conjugate gradient algorithm depending on hot start or
! outer loop index.
!
          IF (inner.eq.0) THEN
            lcgini=.true.
            DO ng=1,ngrids
              CALL congrad (ng, irpm, outer, inner, ninner, lcgini)
            END DO
          END IF
!
!  If initialization step, skip the inner-loop computations.
!
          linner=.false.
          IF ((inner.ne.0).or.(nrun.ne.1)) THEN
            IF (((inner.eq.0).and.lhotstart).or.(inner.ne.0)) THEN
              linner=.true.
            END IF
          END IF
!
!  Start inner loop computations.
!
          inner_compute : IF (linner) THEN
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate adjoint model forced with any vector PSI at the observation
!  locations and generate adjoint trajectory, Lambda_n(t).
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
!  Initialize the adjoint model from rest.
!
            DO ng=1,ngrids
              lsen4dvar(ng)=.false.
              CALL ad_initial (ng)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
              wrtmisfit(ng)=.false.
            END DO
 
#  ifdef RPM_RELAXATION
!
!  Adjoint of representer relaxation is not applied during the
!  inner-loops.
!
            DO ng=1,ngrids
              lweakrelax(ng)=.false.
            END DO
#  endif
!
!  Set adjoint history NetCDF parameters.  Define adjoint history
!  file only once to avoid opening too many files.
!
            DO ng=1,ngrids
              wrtforce(ng)=.true.
              IF (nrun.gt.1) ldefadj(ng)=.false.
              fcount=adm(ng)%load
              adm(ng)%Nrec(fcount)=0
              adm(ng)%Rindex=0
            END DO
!
!  Time-step adjoint model backwards forced with current PSI vector.
!
            DO ng=1,ngrids
              IF (master) THEN
                WRITE (stdout,20) 'AD', ng, ntstart(ng), ntend(ng)
              END IF
            END DO
!
#  ifdef SOLVE3D
            CALL ad_main3d (runinterval)
#  else
            CALL ad_main2d (runinterval)
#  endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Write out last weak-constraint forcing (WRTforce is still .TRUE.)
!  record into the adjoint history file.  Note that the weak-constraint
!  forcing is delayed by nADJ time-steps.
!
            DO ng=1,ngrids
#  ifdef DISTRIBUTE
              CALL ad_wrt_his (ng, myrank)
#  else
              CALL ad_wrt_his (ng, -1)
#  endif
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
!
!  Write out adjoint initial condition record into the adjoint
!  history file.
!
            DO ng=1,ngrids
              wrtforce(ng)=.false.
#  ifdef DISTRIBUTE
              CALL ad_wrt_his (ng, myrank)
#  else
              CALL ad_wrt_his (ng, -1)
#  endif
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
!
!  Convolve adjoint trajectory with error covariances.
!
            lposterior=.false.
            CALL error_covariance (itlm, driver, outer, inner,          &
     &                             lbck, lini, lold, lnew,              &
     &                             rec1, rec2, lposterior)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Convert the current adjoint solution in ADM(ng)%name to impulse
!  forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
!  To facilitate the forcing to the TLM and RPM, the forcing is
!  processed and written in increasing time coordinates (recall that
!  the adjoint solution in ADM(ng)%name is backwards in time).
!
            IF (master) THEN
              WRITE (stdout,40) outer, inner
            END IF
            DO ng=1,ngrids
              tlf(ng)%Rindex=0
#  ifdef DISTRIBUTE
              CALL wrt_impulse (ng, myrank, iadm, adm(ng)%name)
#  else
              CALL wrt_impulse (ng, -1, iadm, adm(ng)%name)
#  endif
 
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate tangent linear model forced by the convolved adjoint
!  trajectory (impulse forcing) to compute R_n * PSI at observation
!  points.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
            DO ng=1,ngrids
              wrtnlmod(ng)=.false.
              wrttlmod(ng)=.true.
            END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
            DO ng=1,ngrids
              IF (frcrec(ng).gt.3) THEN
                frequentimpulse(ng)=.true.
              END IF
            END DO
!
!  Initialize tangent linear model from ITL(ng)%name, record Rec1.
!
            DO ng=1,ngrids
              itl(ng)%Rindex=rec1
              CALL tl_initial (ng)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
!
!  Activate switch to write out initial misfit between model and
!  observations.
!
            IF ((outer.eq.1).and.(inner.eq.1)) THEN
              DO ng=1,ngrids
                wrtmisfit(ng)=.true.
              END DO
            END IF
!
!  Set tangent linear history NetCDF parameters.  Define tangent linear
!  history file at the beggining of each inner loop  to avoid opening
!  too many NetCDF files.
!
            DO ng=1,ngrids
              IF (inner.gt.1) ldeftlm(ng)=.false.
              fcount=tlm(ng)%load
              tlm(ng)%Nrec(fcount)=0
              tlm(ng)%Rindex=0
            END DO
!
!  Run tangent linear model forward and force with convolved adjoint
!  trajectory impulses. Compute R_n * PSI at observation points which
!  are used in the conjugate gradient algorithm.
!
            DO ng=1,ngrids
              IF (master) THEN
                WRITE (stdout,20) 'TL', ng, ntstart(ng), ntend(ng)
              END IF
            END DO
!
#  ifdef SOLVE3D
            CALL tl_main3d (runinterval)
#  else
            CALL tl_main2d (runinterval)
#  endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
            DO ng=1,ngrids
              wrtnlmod(ng)=.false.
              wrttlmod(ng)=.false.
            END DO
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Use conjugate gradient algorithm to find a better approximation
!  PSI to representer coefficients Beta_n.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
            nrun=nrun+1
            DO ng=1,ngrids
              lcgini=.false.
              CALL congrad (ng, irpm, outer, inner, ninner, lcgini)
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
 
          END IF inner_compute
 
        END DO inner_loop
!
!  Close tangent linear NetCDF file.
!
        sourcefile=myfile
        DO ng=1,ngrids
          CALL close_file (ng, itlm, tlm(ng))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!-----------------------------------------------------------------------
!  Once that the representer coefficients, Beta_n, have been
!  approximated with sufficient accuracy, compute estimates of
!  Lambda_n and Xhat_n by carrying out one backward intergration
!  of the adjoint model and one forward itegration of the representer
!  model.
!-----------------------------------------------------------------------
!
!  Initialize the adjoint model always from rest.
!
        DO ng=1,ngrids
          lsen4dvar(ng)=.false.
          CALL ad_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
 
#  ifdef RPM_RELAXATION
!
!  Adjoint of representer relaxation is applied during the
!  outer-loops.
!
        DO ng=1,ngrids
          lweakrelax(ng)=.true.
        END DO
#  endif
!
!  Set adjoint history NetCDF parameters.  Define adjoint history
!  file one to avoid opening to many files.
!
        DO ng=1,ngrids
          wrtforce(ng)=.true.
          IF (nrun.gt.1) ldefadj(ng)=.false.
          fcount=adm(ng)%load
          adm(ng)%Nrec(fcount)=0
          adm(ng)%Rindex=0
        END DO
!
!  Time-step adjoint model backwards forced with estimated representer
!  coefficients, Beta_n.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) 'AD', ng, ntstart(ng), ntend(ng)
          END IF
        END DO
!
#  ifdef SOLVE3D
        CALL ad_main3d (runinterval)
#  else
        CALL ad_main2d (runinterval)
#  endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Write out last weak-constraint forcing (WRTforce is still .TRUE.)
!  record into the adjoint history file.  Note that the weak-constraint
!  forcing is delayed by nADJ time-steps.
!
        DO ng=1,ngrids
#  ifdef DISTRIBUTE
          CALL ad_wrt_his (ng, myrank)
#  else
          CALL ad_wrt_his (ng, -1)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Write out adjoint initial condition record into the adjoint
!  history file.
!
        DO ng=1,ngrids
          wrtforce(ng)=.false.
#  ifdef DISTRIBUTE
          CALL ad_wrt_his (ng, myrank)
#  else
          CALL ad_wrt_his (ng, -1)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Convolve adjoint trajectory with error covariances.
!
        lposterior=.false.
        CALL error_covariance (irpm, driver, outer, inner,              &
     &                         lbck, lini, lold, lnew,                  &
     &                         rec1, rec2, lposterior)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Convert the current adjoint solution in ADM(ng)%name to impulse
!  forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
!  To facilitate the forcing to the TLM and RPM, the forcing is
!  processed and written in increasing time coordinates (recall that
!  the adjoint solution in ADM(ng)%name is backwards in time).
!
        IF (master) THEN
          WRITE (stdout,40) outer, inner
        END IF
        DO ng=1,ngrids
          tlf(ng)%Rindex=0
#  ifdef DISTRIBUTE
          CALL wrt_impulse (ng, myrank, iadm, adm(ng)%name)
#  else
          CALL wrt_impulse (ng, -1, iadm, adm(ng)%name)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
 
# endif /* !DATALESS_LOOPS */
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Run representer model and compute a "new estimate" of the state
!  trajectory, X_n(t).
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
!  Set new basic state trajectory for next outer loop.
!
        DO ng=1,ngrids
          ideftlm(ng)=-1
          ldeftlm(ng)=.true.
          lwrttlm(ng)=.true.
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.true.
          wrtrpmod(ng)=.true.
          WRITE (tlm(ng)%name,10) trim(fwd(ng)%head), outer
          lstr=len_trim(tlm(ng)%name)
          tlm(ng)%base=tlm(ng)%name(1:lstr-3)
        END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
        DO ng=1,ngrids
          IF (frcrec(ng).gt.3) THEN
            frequentimpulse(ng)=.true.
          END IF
        END DO
!
!  Initialize representer model IRP(ng)%name file, record Rec2.
!
        DO ng=1,ngrids
# ifdef DATALESS_LOOPS
          irp(ng)%Rindex=rec1
# else
          irp(ng)%Rindex=rec2
# endif
          CALL rp_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Activate switch to write out final misfit between model and
!  observations.
!
        IF (outer.eq.nouter) THEN
          DO ng=1,ngrids
            wrtmisfit(ng)=.true.
          END DO
        END IF
!
!  Run representer model using previous linearized trajectory, X_n-1, as
!  basic state and forced with convolved adjoint trajectory impulses.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) 'RP', ng, ntstart(ng), ntend(ng)
          END IF
        END DO
!
# ifdef SOLVE3D
        CALL rp_main3d (runinterval)
# else
        CALL rp_main2d (runinterval)
# endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Report data penalty function.
!
        DO ng=1,ngrids
          IF (master) THEN
            DO i=0,nstatevar(ng)
              IF (i.eq.0) THEN
                string='Total'
              ELSE
                string=vname(1,idsvar(i))
              END IF
              IF (fourdvar(ng)%DataPenalty(i).ne.0.0_r8) THEN
                WRITE (stdout,30) outer, inner, 'RPM',                  &
     &                            fourdvar(ng)%DataPenalty(i),          &
     &                            trim(string)
# ifdef DATALESS_LOOPS
                WRITE (stdout,30) outer, inner, 'NLM',                  &
     &                            fourdvar(ng)%NLPenalty(i),            &
     &                            trim(string)
# endif
              END IF
            END DO
          END IF
        END DO
!
!  Write out final data penalty function to NetCDF file.
!
        sourcefile=myfile
        DO ng=1,ngrids
          SELECT CASE (dav(ng)%IOtype)
            CASE (io_nf90)
              CALL netcdf_put_fvar (ng, irpm, dav(ng)%name,             &
     &                              'RP_fDataPenalty',                  &
     &                              fourdvar(ng)%DataPenalty(0:),       &
     &                              (/1,outer/),                        &
     &                              (/nstatevar(ng)+1,1/),              &
     &                              ncid = dav(ng)%ncid)
 
# if defined PIO_LIB && defined DISTRIBUTE
            CASE (io_pio)
              CALL pio_netcdf_put_fvar (ng, irpm, dav(ng)%name,         &
     &                                  'RP_fDataPenalty',              &
     &                                  fourdvar(ng)%DataPenalty(0:),   &
     &                                  (/1,outer/),                    &
     &                                  (/nstatevar(ng)+1,1/),          &
     &                                  piofile = dav(ng)%pioFile)
# endif
          END SELECT
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Clean array before next run of RP model.
!
        DO ng=1,ngrids
          fourdvar(ng)%DataPenalty=0.0_r8
# ifdef DATALESS_LOOPS
          fourdvar(ng)%NLPenalty=0.0_r8
# endif
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.false.
        END DO
!
!  Close current forward NetCDF file.
!
        DO ng=1,ngrids
          CALL close_file (ng, irpm, fwd(ng))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          IF (his(ng)%IOtype.eq.io_nf90) THEN
            his(ng)%ncid=-1
#if defined PIO_LIB && defined DISTRIBUTE
          ELSE IF (his(ng)%IOtype.eq.io_pio) THEN
            his(ng)%pioFile%fh=-1
#endif
          END IF
        END DO
 
      END DO outer_loop
 
#endif /* RECOMPUTE_4DVAR */
!!
!! Compute and report model-observation comparison statistics.
!!
!!    DO ng=1,Ngrids
#ifdef DISTRIBUTE
!!      CALL stats_modobs (ng, MyRank)
#else
!!      CALL stats_modobs (ng, -1)
#endif
!!    END DO
!
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!  Adjoint of W4DVar to compute the observation sensitivity.
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
!
!  Reset the start and end times for the adjoint forcing.
!
      DO ng=1,ngrids
        str_day=tdays(ng)
        end_day=str_day-ntimes(ng)*dt(ng)*sec2day
        IF ((dstrs(ng).eq.0.0_r8).and.(dends(ng).eq.0.0_r8)) THEN
          dstrs(ng)=end_day
          dends(ng)=str_day
        END IF
        IF (master) THEN
          WRITE (stdout,70) 'AD', dends(ng), dstrs(ng)
        END IF
      END DO
!
!  WARNING: ONLY 1 outer loop can be used for this application.
!  =======  For more than 1 outer-loop, we require the second
!  derivative of each model operator (i.e. the tangent linear
!  of the tangent linear operator).
!
      ad_outer_loop : DO my_outer=1,1,-1
        outer=my_outer
        inner=0
!
!-----------------------------------------------------------------------
!  Run the adjoint model initialized and forced by dI/dx where I is the
!  chosen function of the analysis/forecast state x.
!-----------------------------------------------------------------------
!
!  Set basic state trajectory.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(fwd(ng)%head), outer-1
          lstr=len_trim(fwd(ng)%name)
          fwd(ng)%base=fwd(ng)%name(1:lstr-3)
        END DO
        IF ((outer.eq.1).and.master) THEN
          WRITE (stdout,50)
        END IF
!
!  Initialize the adjoint model: initialize using dI/dxf is
!  appropriate.
!
        lstiffness=.false.
        DO ng=1,ngrids
          lsen4dvar(ng)=.true.
          CALL ad_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Set adjoint history NetCDF parameters.  Define adjoint history
!  file one to avoid opening to many files.
!
        DO ng=1,ngrids
          wrtforce(ng)=.true.
          IF (nrun.gt.1) ldefadj(ng)=.false.
          fcount=adm(ng)%load
          adm(ng)%Nrec(fcount)=0
          adm(ng)%Rindex=0
        END DO
!
!  NOTE: THE ADM IS FORCED BY dI/dx ONLY when outer=Nouter.
!
!  Time-step adjoint model backwards.
!  ??? What do we do in the case of model error? Save forcing for TLM?
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) 'AD', ng, ntstart(ng), ntend(ng)
          END IF
        END DO
!
#ifdef SOLVE3D
        CALL ad_main3d (runinterval)
#else
        CALL ad_main2d (runinterval)
#endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Write out last weak-constraint forcing (WRTforce is still .TRUE.)
!  record into the adjoint history file.  Note that the weak-constraint
!  forcing is delayed by nADJ time-steps.
!
        DO ng=1,ngrids
#ifdef DISTRIBUTE
          CALL ad_wrt_his (ng, myrank)
#else
          CALL ad_wrt_his (ng, -1)
#endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Write out adjoint initial condition record into the adjoint
!  history file.
!
        DO ng=1,ngrids
          wrtforce(ng)=.false.
#ifdef DISTRIBUTE
          CALL ad_wrt_his (ng, myrank)
#else
          CALL ad_wrt_his (ng, -1)
#endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Convolve adjoint trajectory with error covariances.
!
        lposterior=.false.
        CALL error_covariance (itlm, driver, outer, inner,              &
     &                         lbck, lini, lold, lnew,                  &
     &                         rec1, rec2, lposterior)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Convert the current adjoint solution in ADM(ng)%name to impulse
!  forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
!  To facilitate the forcing to the TLM and RPM, the forcing is
!  processed and written in increasing time coordinates (recall that
!  the adjoint solution in ADM(ng)%name is backwards in time).
!
!  AMM: Do not know what to do in the weak constraint case yet.
!
!!      IF (Master) THEN
!!        WRITE (stdout,40) outer, inner
!!      END IF
!!      DO ng=1,Ngrids
!!        TLF(ng)%Rindex=0
#ifdef DISTRIBUTE
!!        CALL wrt_impulse (ng, MyRank, iADM, ADM(ng)%name)
#else
!!        CALL wrt_impulse (ng, -1, iADM, ADM(ng)%name)
#endif
!!        IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN
!!      END DO
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate tangent linear model forced by the convolved adjoint
!  trajectory.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.true.
          lwrttlm(ng)=.false.
        END DO
!
!  Clear tangent linear forcing arrays before entering inner-loop.
!  This is very important since these arrays are non-zero after
!  running the representer model and must be zero when running the
!  tangent linear model.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
#ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
#endif
          END DO
        END DO
!
!  Set basic state trajectory.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(fwd(ng)%head), outer-1
          lstr=len_trim(fwd(ng)%name)
          fwd(ng)%base=fwd(ng)%name(1:lstr-3)
        END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
        DO ng=1,ngrids
          IF (frcrec(ng).gt.3) THEN
            frequentimpulse(ng)=.true.
          END IF
        END DO
!
!  Initialize tangent linear model from ITL(ng)%name, record 1.
!
        DO ng=1,ngrids
          itl(ng)%Rindex=rec1
          CALL tl_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
 
#ifdef RPM_RELAXATION
!
!  Activate the tangent linear relaxation terms if used.
!
        DO ng=1,ngrids
          lweakrelax(ng)=.true.
        END DO
#endif
!
!  Run tangent linear model forward and force with convolved adjoint
!  trajectory impulses. Compute (HMBM'H')_n * PSI at observation points
!  which are used in the conjugate gradient algorithm.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) 'TL', ng, ntstart(ng), ntend(ng)
          END IF
        END DO
!
#ifdef SOLVE3D
        CALL tl_main3d (runinterval)
#else
        CALL tl_main2d (runinterval)
#endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.false.
        END DO
 
#ifdef OBS_IMPACT
!
!  Compute observation impact to the data assimilation system.
!
        DO ng=1,ngrids
          CALL rep_matrix (ng, itlm, outer, ninner)
        END DO
#else
!
!  Set basic state trajectory for adjoint inner-loops.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(fwd(ng)%head), outer-1
          lstr=len_trim(fwd(ng)%name)
          fwd(ng)%base=fwd(ng)%name(1:lstr-3)
        END DO
!
!  Clear tangent linear forcing arrays before entering inner-loop.
!  This is very important since these arrays are non-zero after
!  running the representer model and must be zero when running the
!  tangent linear model.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
# ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
# endif
          END DO
        END DO
!
        ad_inner_loop : DO my_inner=ninner,1,-1
          inner=my_inner
 
          IF (master) THEN
            WRITE (stdout,60) 'Adjoint of', uppercase('w4dvar'),        &
     &                        outer, inner
          END IF
!
!  Call adjoint conjugate gradient algorithm.
!
          lcgini=.false.
          DO ng=1,ngrids
            CALL ad_congrad (ng, iadm, outer, inner, ninner, lcgini)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
!
!  Initialize the adjoint model from rest.
!
          DO ng=1,ngrids
            lsen4dvar(ng)=.false.
            CALL ad_initial (ng)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
            wrtmisfit(ng)=.false.
          END DO
 
# ifdef RPM_RELAXATION
!
!  Adjoint of representer relaxation is not applied during the
!  inner-loops.
!
          DO ng=1,ngrids
            lweakrelax(ng)=.false.
          END DO
# endif
!
!  Set adjoint history NetCDF parameters.  Define adjoint history
!  file only once to avoid opening too many files.
!
          DO ng=1,ngrids
            wrtforce(ng)=.true.
            IF (nrun.gt.1) ldefadj(ng)=.false.
            fcount=adm(ng)%load
            adm(ng)%Nrec(fcount)=0
            adm(ng)%Rindex=0
          END DO
!
!  Time-step adjoint model backwards forced with current PSI vector.
!
          DO ng=1,ngrids
            IF (master) THEN
              WRITE (stdout,20) 'AD', ng, ntstart(ng), ntend(ng)
            END IF
          END DO
!
# ifdef SOLVE3D
          CALL ad_main3d (runinterval)
# else
          CALL ad_main2d (runinterval)
# endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Write out last weak-constraint forcing (WRTforce is still .TRUE.)
!  record into the adjoint history file.  Note that the weak-constraint
!  forcing is delayed by nADJ time-steps
!
          DO ng=1,ngrids
# ifdef DISTRIBUTE
            CALL ad_wrt_his (ng, myrank)
# else
            CALL ad_wrt_his (ng, -1)
# endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
!
!  Write out adjoint initial condition record into the adjoint
!  history file.
!
          DO ng=1,ngrids
            wrtforce(ng)=.false.
# ifdef DISTRIBUTE
            CALL ad_wrt_his (ng, myrank)
# else
            CALL ad_wrt_his (ng, -1)
# endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
!
!  Convolve adjoint trajectory with error covariances.
!
          lposterior=.false.
          CALL error_covariance (itlm, driver, outer, inner,            &
     &                           lbck, lini, lold, lnew,                &
     &                           rec1, rec2, lposterior)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  ??? Do we need the adjoint of impulse here???
!
!  Convert the current adjoint solution in ADM(ng)%name to impulse
!  forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
!  To facilitate the forcing to the TLM and RPM, the forcing is
!  processed and written in increasing time coordinates (recall that
!  the adjoint solution in ADM(ng)%name is backwards in time).
!
!!        IF (Master) THEN
!!          WRITE (stdout,40) outer, inner
!!        END IF
!!        DO ng=1,Ngrids
!!          TLF(ng)%Rindex=0
# ifdef DISTRIBUTE
!!          CALL wrt_impulse (ng, MyRank, iADM, ADM(ng)%name)
# else
!!          CALL wrt_impulse (ng, -1, iADM, ADM(ng)%name)
# endif
!!          IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN
!!
!! AMM: Don't know what to do yet in the weak constraint case.
!!
!!        END DO
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate tangent linear model.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
          DO ng=1,ngrids
            tlm(ng)%name=trim(tlm(ng)%base)//'.nc'
            wrtnlmod(ng)=.false.
            wrttlmod(ng)=.true.
          END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
          DO ng=1,ngrids
            IF (frcrec(ng).gt.3) THEN
              frequentimpulse(ng)=.true.
            END IF
          END DO
!
!  Initialize tangent linear model from ITL(ng)%name, record Rec1.
!
          DO ng=1,ngrids
            itl(ng)%Rindex=rec1
            CALL tl_initial (ng)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
 
# ifdef RPM_RELAXATION
!
!  Deactivate the tangent linear relaxation terms if used.
!
          DO ng=1,ngrids
            lweakrelax(ng)=.false.
          END DO
# endif
!
!  Activate switch to write out initial misfit between model and
!  observations.
!
          IF ((outer.eq.1).and.(inner.eq.1)) THEN
            DO ng=1,ngrids
              wrtmisfit(ng)=.true.
            END DO
          END IF
!
!  Set tangent linear history NetCDF parameters.  Define tangent linear
!  history file at the beggining of each inner loop  to avoid opening
!  too many NetCDF files.
!
          DO ng=1,ngrids
            IF (inner.lt.ninner) ldeftlm(ng)=.false.
            fcount=tlm(ng)%load
            tlm(ng)%Nrec(fcount)=0
            tlm(ng)%Rindex=0
          END DO
!
!  Run tangent linear model forward and force with convolved adjoint
!  trajectory impulses. Compute R_n * PSI at observation points which
!  are used in the conjugate gradient algorithm.
!
          DO ng=1,ngrids
            IF (master) THEN
              WRITE (stdout,20) 'TL', ng, ntstart(ng), ntend(ng)
            END IF
          END DO
!
# ifdef SOLVE3D
          CALL tl_main3d (runinterval)
# else
          CALL tl_main2d (runinterval)
# endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          DO ng=1,ngrids
            wrtnlmod(ng)=.false.
            wrttlmod(ng)=.false.
          END DO
 
        END DO ad_inner_loop
!
!  Call adjoint conjugate gradient algorithm.
!
        inner=0
        lcgini=.true.
        DO ng=1,ngrids
          CALL ad_congrad (ng, itlm, outer, inner, ninner, lcgini)
        END DO
 
#endif /* !OBS_IMPACT */
 
#ifdef OBS_IMPACT
!
!  Write out total observation impact.
!
        IF (outer.eq.1) THEN
          sourcefile=myfile
          DO ng=1,ngrids
            SELECT CASE (dav(ng)%IOtype)
              CASE (io_nf90)
                CALL netcdf_put_fvar (ng, itlm, dav(ng)%name,           &
     &                                'ObsImpact_total', ad_obsval,     &
     &                                (/1/), (/mobs/),                  &
     &                                ncid = dav(ng)%ncid)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
                CALL netcdf_sync (ng, inlm, dav(ng)%name, dav(ng)%ncid)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
# if defined PIO_LIB && defined DISTRIBUTE
              CASE (io_pio)
                CALL pio_netcdf_put_fvar (ng, itlm, dav(ng)%name,       &
     &                                    'ObsImpact_total', ad_obsval, &
     &                                    (/1/), (/mobs/),              &
     &                                    piofile = dav(ng)%pioFile)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
                CALL pio_netcdf_sync (ng, inlm, dav(ng)%name,           &
     &                                dav(ng)%pioFile)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
# endif
            END SELECT
          END DO
        END IF
#else
!
!  Write out observation sensitivity.
!
        IF (outer.eq.1) THEN
          sourcefile=myfile
          DO ng=1,ngrids
            SELECT CASE (dav(ng)%IOtype)
              CASE (io_nf90)
                CALL netcdf_put_fvar (ng, itlm, dav(ng)%name,           &
     &                                'ObsSens_total', ad_obsval,       &
     &                                (/1/), (/mobs/),                  &
     &                                ncid = dav(ng)%ncid)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
                CALL netcdf_sync (ng, inlm, dav(ng)%name, dav(ng)%ncid)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
# if defined PIO_LIB && defined DISTRIBUTE
              CASE (io_pio)
                CALL pio_netcdf_put_fvar (ng, itlm, dav(ng)%name,       &
     &                                    'ObsSens_total', ad_obsval,   &
     &                                    (/1/), (/mobs/),              &
     &                                    piofile = dav(ng)%pioFile)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
                CALL pio_netcdf_sync (ng, inlm, dav(ng)%name,           &
     &                                dav(ng)%pioFile)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
# endif
            END SELECT
          END DO
        END IF
#endif
!
!  Close tangent linear NetCDF file.
!
        sourcefile=myfile
        DO ng=1,ngrids
          CALL close_file (ng, itlm, tlm(ng))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
 
#if defined OBS_IMPACT && defined OBS_IMPACT_SPLIT
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate tangent linear model with initial condition increments
!  only to compute the observation impact associated with the initial
!  conditions.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.true.
          lwrttlm(ng)=.false.
        END DO
!
!  Clear tangent linear forcing arrays before entering inner-loop.
!  This is very important since these arrays are non-zero after
!  running the representer model and must be zero when running the
!  tangent linear model.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
# ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
# endif
          END DO
        END DO
!
!  Set basic state trajectory.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(fwd(ng)%head), outer-1
          lstr=len_trim(fwd(ng)%name)
          fwd(ng)%base=fwd(ng)%name(1:lstr-3)
        END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
        DO ng=1,ngrids
          IF (frcrec(ng).gt.3) THEN
            frequentimpulse(ng)=.true.
          END IF
        END DO
!
!  Initialize tangent linear model from ITL(ng)%name, record 1.
!
        DO ng=1,ngrids
          itl(ng)%Rindex=rec1
          CALL tl_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Clear tangent linear forcing arrays and boundary arrays
!  before the obs impact initial condition calculation.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
# ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
# endif
          END DO
        END DO
 
# ifdef RPM_RELAXATION
!
!  Activate the tangent linear relaxation terms if used.
!
        DO ng=1,ngrids
          lweakrelax(ng)=.true.
        END DO
# endif
!
!  Run tangent linear model forward and force with convolved adjoint
!  trajectory impulses. Compute (HMBM'H')_n * PSI at observation points
!  which are used in the conjugate gradient algorithm.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) 'TL', ng, ntstart(ng), ntend(ng)
          END IF
        END DO
!
# ifdef SOLVE3D
        CALL tl_main3d (runinterval)
# else
        CALL tl_main2d (runinterval)
# endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.false.
        END DO
!
!  Compute observation impact to the data assimilation system.
!
        DO ng=1,ngrids
          CALL rep_matrix (ng, itlm, outer, ninner)
        END DO
!
!  Write out observation sentivity.
!
        IF (outer.eq.1) THEN
          sourcefile=myfile
          DO ng=1,ngrids
            SELECT CASE (dav(ng)%IOtype)
              CASE (io_nf90)
                CALL netcdf_put_fvar (ng, itlm, dav(ng)%name,           &
     &                                'ObsImpact_IC', ad_obsval,        &
     &                                (/1/), (/mobs/),                  &
     &                                ncid = dav(ng)%ncid)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
                CALL netcdf_sync (ng, inlm, dav(ng)%name, dav(ng)%ncid)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
# if defined PIO_LIB && defined DISTRIBUTE
              CASE (io_pio)
                CALL pio_netcdf_put_fvar (ng, itlm, dav(ng)%name,       &
     &                                    'ObsImpact_IC', ad_obsval,    &
     &                                    (/1/), (/mobs/),              &
     &                                    piofile = dav(ng)%pioFile)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
                CALL pio_netcdf_sync (ng, inlm, dav(ng)%name,           &
     &                                dav(ng)%pioFile)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
# endif
            END SELECT
          END DO
        END IF
 
# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate tangent linear model with surface forcing increments
!  only to compute the observations impact associated with the
!  surface forcing.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.true.
          lwrttlm(ng)=.false.
        END DO
!
!  Clear tangent linear forcing arrays before entering inner-loop.
!  This is very important since these arrays are non-zero after
!  running the representer model and must be zero when running the
!  tangent linear model.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
#  ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
#  endif
          END DO
        END DO
!
!  Set basic state trajectory.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(fwd(ng)%head), outer-1
          lstr=len_trim(fwd(ng)%name)
          fwd(ng)%base=fwd(ng)%name(1:lstr-3)
        END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
        DO ng=1,ngrids
          IF (frcrec(ng).gt.3) THEN
            frequentimpulse(ng)=.true.
          END IF
        END DO
!
!  Initialize tangent linear model from ITL(ng)%name, record 1.
!
        DO ng=1,ngrids
          itl(ng)%Rindex=rec1
          CALL tl_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Clear tangent initial condition arrays and boundary arrays
!  before the obs impact initial condition calculation.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_ocean (ng, tile, itlm)
#  ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
#  endif
          END DO
        END DO
 
#  ifdef RPM_RELAXATION
!
!  Activate the tangent linear relaxation terms if used.
!
        DO ng=1,ngrids
          lweakrelax(ng)=.true.
        END DO
#  endif
!
!  Run tangent linear model forward and force with convolved adjoint
!  trajectory impulses. Compute (HMBM'H')_n * PSI at observation points
!  which are used in the conjugate gradient algorithm.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) 'TL', ng, ntstart(ng), ntend(ng)
          END IF
        END DO
!
#  ifdef SOLVE3D
        CALL tl_main3d (runinterval)
#  else
        CALL tl_main2d (runinterval)
#  endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.false.
        END DO
!
!  Compute observation impact to the data assimilation system.
!
        DO ng=1,ngrids
          CALL rep_matrix (ng, itlm, outer, ninner)
        END DO
!
!  Write out observation sentivity.
!
        IF (outer.eq.1) THEN
          sourcefile=myfile
          DO ng=1,ngrids
            SELECT CASE (dav(ng)%IOtype)
              CASE (io_nf90)
                CALL netcdf_put_fvar (ng, itlm, dav(ng)%name,           &
     &                                'ObsImpact_FC', ad_obsval,        &
     &                                (/1/), (/mobs/),                  &
     &                                ncid = dav(ng)%ncid)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
                CALL netcdf_sync (ng, inlm, dav(ng)%name, dav(ng)%ncid)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
#  if defined PIO_LIB && defined DISTRIBUTE
              CASE (io_pio)
                CALL pio_netcdf_put_fvar (ng, itlm, dav(ng)%name,       &
     &                                    'ObsImpact_FC', ad_obsval,    &
     &                                    (/1/), (/mobs/),              &
     &                                    piofile = dav(ng)%pioFile)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
                CALL pio_netcdf_sync (ng, inlm, dav(ng)%name,           &
     &                                dav(ng)%pioFile)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
#  endif
            END SELECT
          END DO
        END IF
# endif
 
# if defined ADJUST_BOUNDARY
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  Integrate tangent linear model with boundary condition increments
!  only to compute the observation impact associated with the boundary
!  conditions.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.true.
          lwrttlm(ng)=.false.
        END DO
!
!  Clear tangent linear forcing arrays before entering inner-loop.
!  This is very important since these arrays are non-zero after
!  running the representer model and must be zero when running the
!  tangent linear model.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_forces (ng, tile, itlm)
#  ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
#  endif
          END DO
        END DO
!
!  Set basic state trajectory.
!
        DO ng=1,ngrids
          WRITE (fwd(ng)%name,10) trim(fwd(ng)%head), outer-1
          lstr=len_trim(fwd(ng)%name)
          fwd(ng)%base=fwd(ng)%name(1:lstr-3)
        END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
        DO ng=1,ngrids
          IF (frcrec(ng).gt.3) THEN
            frequentimpulse(ng)=.true.
          END IF
        END DO
!
!  Initialize tangent linear model from ITL(ng)%name, record Rec1.
!
        DO ng=1,ngrids
          itl(ng)%Rindex=rec1
          CALL tl_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Clear tangent linear initial condition and forcing arrays
!  before the obs impact initial condition calculation.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_ocean (ng, tile, itlm)
            CALL initialize_forces (ng, tile, itlm)
          END DO
        END DO
 
#  ifdef RPM_RELAXATION
!
!  Activate the tangent linear relaxation terms if used.
!
        DO ng=1,ngrids
          lweakrelax(ng)=.true.
        END DO
#  endif
!
!  Run tangent linear model forward and force with convolved adjoint
!  trajectory impulses. Compute (HMBM'H')_n * PSI at observation points
!  which are used in the conjugate gradient algorithm.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) 'TL', ng, ntstart(ng), ntend(ng)
          END IF
        END DO
!
#  ifdef SOLVE3D
        CALL tl_main3d (runinterval)
#  else
        CALL tl_main2d (runinterval)
#  endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        DO ng=1,ngrids
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.false.
        END DO
!
!  Compute observation impact to the data assimilation system.
!
        DO ng=1,ngrids
          CALL rep_matrix (ng, itlm, outer, ninner)
        END DO
!
!  Write out observation sentivity.
!
        IF (outer.eq.1) THEN
          sourcefile=myfile
          DO ng=1,ngrids
            SELECT CASE (dav(ng)%IOtype)
              CASE (io_nf90)
                CALL netcdf_put_fvar (ng, itlm, dav(ng)%name,           &
     &                                'ObsImpact_BC', ad_obsval,        &
     &                                (/1/), (/mobs/),                  &
     &                                ncid = dav(ng)%ncid)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
                CALL netcdf_sync (ng, inlm, dav(ng)%name, dav(ng)%ncid)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
#  if defined PIO_LIB && defined DISTRIBUTE
              CASE (io_pio)
                CALL pio_netcdf_put_fvar (ng, itlm, dav(ng)%name,       &
     &                                    'ObsImpact_BC', ad_obsval,    &
     &                                    (/1/), (/mobs/),              &
     &                                    piofile = dav(ng)%pioFile)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
 
                CALL pio_netcdf_sync (ng, inlm, dav(ng)%name,           &
     &                                dav(ng)%pioFile)
                IF (founderror(exit_flag, noerror,                      &
     &                         __line__, myfile)) RETURN
#  endif
            END SELECT
          END DO
        END IF
# endif
#endif /* OBS_IMPACT_SPLIT */
!
!  Close current forward NetCDF file.
!
        sourcefile=myfile
        DO ng=1,ngrids
          CALL close_file (ng, inlm, fwd(ng))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          IF (his(ng)%IOtype.eq.io_nf90) THEN
            his(ng)%ncid=-1
#if defined PIO_LIB && defined DISTRIBUTE
          ELSE IF (his(ng)%IOtype.eq.io_pio) THEN
            his(ng)%pioFile%fh=-1
#endif
          END IF
        END DO
 
      END DO ad_outer_loop
!
 10   FORMAT (a,'_outer',i0,'.nc')
 20   FORMAT (/,1x,a,1x,'ROMS: started time-stepping:',                 &
     &        ' (Grid: ',i2.2,' TimeSteps: ',i8.8,' - ',i8.8,')',/)
 30   FORMAT (' (',i3.3,',',i3.3,'): ',a,' data penalty, Jdata = ',     &
     &        1p,e17.10,0p,t68,a)
 40   FORMAT (/,' Converting Convolved Adjoint Trajectory to',          &
     &          ' Impulses: Outer = ',i3.3,' Inner = ',i3.3,/)
 50   FORMAT (/,'ROMS: Started adjoint Sensitivity calculation',        &
     &          ' ...',/)
 60   FORMAT (/,'ROMS: ',a,1x,a,', Outer = ',i3.3,                      &
     &          ' Inner = ',i3.3,/)
 70   FORMAT (/,1x,a,1x,'ROMS: adjoint forcing time range: ',           &
     &        f12.4,' - ',f12.4 ,/)
!
      RETURN
      END SUBROUTINE roms_run
!
      SUBROUTINE roms_finalize
!
!=======================================================================
!                                                                      !
!  This routine terminates ROMS nonlinear, tangent linear, and         !
!  adjoint models execution.                                           !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_iounits
      USE mod_ncparam
      USE mod_scalars
!
!  Local variable declarations.
!
      integer :: Fcount, ng, thread
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", ROMS_finalize"
!
!-----------------------------------------------------------------------
!  Read and write observation variables for completeness.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
#ifdef DISTRIBUTE
        CALL stats_modobs (ng, myrank)
#else
        CALL stats_modobs (ng, -1)
#endif
      END DO
!
!-----------------------------------------------------------------------
!  If blowing-up, save latest model state into RESTART NetCDF file.
!-----------------------------------------------------------------------
!
!  If cycling restart records, write solution into record 3.
!
      IF (exit_flag.eq.1) THEN
        DO ng=1,ngrids
          IF (lwrtrst(ng)) THEN
            IF (master) WRITE (stdout,10)
 10         FORMAT (/,' Blowing-up: Saving latest model state into ',   &
     &                ' RESTART file',/)
            fcount=rst(ng)%load
            IF (lcyclerst(ng).and.(rst(ng)%Nrec(fcount).ge.2)) THEN
              rst(ng)%Rindex=2
              lcyclerst(ng)=.false.
            END IF
            blowup=exit_flag
            exit_flag=noerror
#ifdef DISTRIBUTE
            CALL wrt_rst (ng, myrank)
#else
            CALL wrt_rst (ng, -1)
#endif
          END IF
        END DO
      END IF
!
!-----------------------------------------------------------------------
!  Stop model and time profiling clocks, report memory requirements, and
!  close output NetCDF files.
!-----------------------------------------------------------------------
!
!  Stop time clocks.
!
      IF (master) THEN
        WRITE (stdout,20)
 20     FORMAT (/,'Elapsed wall CPU time for each process (seconds):',/)
      END IF
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_off (ng, inlm, 0, __line__, myfile)
        END DO
      END DO
!
!  Report dynamic memory and automatic memory requirements.
!
      CALL memory
!
!  Close IO files.
!
      DO ng=1,ngrids
        CALL close_inp (ng, inlm)
      END DO
      CALL close_out
!
      RETURN
      END SUBROUTINE roms_finalize
 
      END MODULE roms_kernel_mod