rbl4dvar.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE rbl4dvar_mod
 
#ifdef RBL4DVAR
!
!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 module splits the RBL4D-Var data assimilation algorithm into   !
!  its logical components routines.                                    !
!                                                                      !
!    background_initialize:                                            !
!                                                                      !
!      Initializes the nonlinear model kernel. It is separated from    !
!      the 'background' phase to allow ESM coupling.                   !
!                                                                      !
!    background:                                                       !
!                                                                      !
!      Timesteps the nonlinear model to compute the basic state Xb(t)  !
!      used to linearize the tangent linear and adjoint models. It     !
!      interpolates the nonlinear model trajectory to the observations !
!      locations.                                                      !
!                                                                      !
!    increment:                                                        !
!                                                                      !
!      Minimizes of the 4D-Var cost function over Ninner inner loops   !
!      iterations to compute the data assimilation increment, dXa. It  !
!      also computes the nonlinear model initial conditions by adding  !
!      the increment to the background, Xa = Xb(t=0) + dXa. The new    !
!      NLM initial conditions is written inot INI(ng)%name NetCDF file.!
!                                                                      !
!    analysis_initialize:                                              !
!                                                                      !
!      Initializes the nonlinear model with the 4D-Var analysis, Xa.   !
!      It is separated from 'analysis' phase to allow ESM coupling.    !
!                                                                      !
!    analysis:                                                         !
!                                                                      !
!      Timesteps the nonlinear model to compute new state trajectory   !
!      that becomes the basic state for the next outer loop. It also   !
!      interpolates the solution at the observation locations.         !
!                                                                      !
!    prior_error:                                                      !
!                                                                      !
!      Processes the background and model prior error covariance and   !
!      its normalization coefficients.                                 !
!                                                                      !
!    posterior_error:                                                  !
!                                                                      !
!      Computes the posterior error analysis.                          !
!                                                                      !
!                                                                      !
!  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.         !
!                                                                      !
!    Gurol, S., A.T. Weaver, A.M. Moore, A. Piacentini, H.G. Arango,   !
!      S. Gratton, 2014: B-preconditioned minimization algorithms for  !
!      data assimilation with the dual formulation, QJRMS, 140,        !
!      539-556.                                                        !
!                                                                      !
!=======================================================================
!
      USE mod_kinds
      USE mod_param
      USE mod_parallel
      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_mixing,         ONLY : initialize_mixing
      USE mod_ocean,          ONLY : initialize_ocean
!
# ifdef SPLIT_4DVAR
      USE ad_def_his_mod,     ONLY : ad_def_his
# endif
      USE ad_wrt_his_mod,     ONLY : ad_wrt_his
# ifdef STD_MODEL
      USE background_std_mod, ONLY : background_std
# endif
# ifndef RPCG
      USE congrad_mod,        ONLY : congrad
#  ifdef SPLIT_4DVAR
      USE congrad_mod,        ONLY : cg_read_congrad
#  endif
# endif
# if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F || \
     defined posterior_eofs
      USE convolve_mod,       ONLY : convolve
# endif
      USE convolve_mod,       ONLY : error_covariance
# if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F
      USE def_error_mod,      ONLY : def_error
# endif
# if defined POSTERIOR_EOFS    || defined POSTERIOR_ERROR_I || \
     defined posterior_error_f
      USE def_hessian_mod,    ONLY : def_hessian
# endif
# ifdef SPLIT_4DVAR
      USE def_ini_mod,        ONLY : def_ini
# endif
      USE def_impulse_mod,    ONLY : def_impulse
      USE def_mod_mod,        ONLY : def_mod
      USE def_norm_mod,       ONLY : def_norm
# ifdef STD_MODEL
      USE def_std_mod,        ONLY : def_std
# endif
# ifdef RPCG
      USE frc_iau_mod,        ONLY : frc_iau_ini
# endif
      USE get_state_mod,      ONLY : get_state
# ifdef POSTERIOR_ERROR_F
      USE ini_adjust_mod,     ONLY : load_tltoad
# endif
      USE normalization_mod,  ONLY : normalization
# if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F || \
     defined posterior_eofs
      USE posterior_mod,      ONLY : posterior
# endif
#  if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F
      USE posterior_var_mod,  ONLY : posterior_var
#  endif
# ifdef RPCG
      USE rpcg_lanczos_mod,   ONLY : rpcg_lanczos
#  ifdef SPLIT_4DVAR
      USE rpcg_lanczos_mod,   ONLY : cg_read_rpcg
#  endif
# endif
      USE strings_mod,        ONLY : founderror
      USE tl_def_ini_mod,     ONLY : tl_def_ini
      USE tl_wrt_ini_mod,     ONLY : tl_wrt_ini
# ifdef RPCG
      USE wrt_aug_imp_mod,    ONLY : wrt_aug_imp
# endif
# if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F
      USE wrt_error_mod,      ONLY : wrt_error
# endif
# ifdef POSTERIOR_ERROR_F
      USE wrt_hessian_mod,    ONLY : wrt_hessian
# endif
      USE wrt_impulse_mod,    ONLY : wrt_impulse
      USE wrt_ini_mod,        ONLY : wrt_ini
# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
     defined adjust_boundary
      USE wrt_ini_mod,        ONLY : wrt_frc_ad
# endif
# ifdef STD_MODEL
      USE wrt_std_mod,        ONLY : wrt_std
# endif
# if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
      USE zeta_balance_mod,   ONLY : balance_ref, biconj
# endif
!
      implicit none
!
      PUBLIC  :: background_initialize
      PUBLIC  :: background
      PUBLIC  :: increment
      PUBLIC  :: analysis_initialize
      PUBLIC  :: analysis
      PUBLIC  :: prior_error
 
# if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F || \
     defined posterior_eofs
      PUBLIC  :: posterior_error
# endif
!
!  Set module internal parameters.
!
# ifdef SPLIT_4DVAR
      logical :: ldone
!
# endif
      integer :: lbck = 2
      integer :: lini = 1
      integer :: ltlm1 = 1
      integer :: ltlm2 = 2
      integer :: rec1 = 1
      integer :: rec2 = 2
      integer :: rec3 = 3
      integer :: rec4 = 4
      integer :: rec5 = 5
!
      CONTAINS
!
      SUBROUTINE background_initialize (my_outer)
!
!=======================================================================
!                                                                      !
!  This routine initializes ROMS nonlinear model trajectory Xb_n-1(0)  !
!  for each (n-1) outer loops. It is separated from the 'background'   !
!  4D-Var phase in ESM coupling applications that use generic methods  !
!  for 'initialize', 'run', and 'finalize'.                            !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     my_outer        Outer-loop counter (integer)                     !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      integer, intent(in) :: my_outer
!
!  Local variable declarations.
!
      integer :: lstr, ng, tile
# ifdef STD_MODEL
      integer :: lstd
# endif
 
      integer :: fcount, tindex
# ifdef PROFILE
      integer :: thread
# endif
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", background_initialize"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Initialize nonlinear model kernel.
!-----------------------------------------------------------------------
 
# ifdef PROFILE
!
!  Start profile clock.
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_on (ng, inlm, 86, __line__, myfile)
        END DO
      END DO
# endif
!
!  Initialize the switch to gather weak constraint forcing.
!
      DO ng=1,ngrids
        wrtforce(ng)=.false.
      END DO
!
!  Clear TLM mixing arrays.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL initialize_mixing (ng, tile, itlm)
        END DO
      END DO
!
!  Get nonlinear model initial conditions.
!
      DO ng=1,ngrids
        wrtnlmod(ng)=.true.
        wrttlmod(ng)=.false.
# ifdef FORWARD_FLUXES
        lreadblk(ng)=.false.
# endif
# ifdef FRC_FILES
        lreadfrc(ng)=.true.
# endif
        lreadfwd(ng)=.false.
        rst(ng)%Rindex=0
        fcount=rst(ng)%load
        rst(ng)%Nrec(fcount)=0
      END DO
!
      CALL initial
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Save nonlinear initial conditions (currently in time index 1,
!  background) into record "Lbck" of INI(ng)%name NetCDF file. The
!  record "Lbck" becomes the background state and "Lini" becomes
!  the current nonlinear initial conditions.
!
      tindex=1
      DO ng=1,ngrids
        ini(ng)%Rindex=1
        fcount=ini(ng)%load
        ini(ng)%Nrec(fcount)=1
# ifdef DISTRIBUTE
        CALL wrt_ini (ng, myrank, tindex, lbck)
# else
        CALL wrt_ini (ng, -1, tindex, lbck)
# endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
 
# ifdef STD_MODEL
!
!  Estimate initial conditions standard deviation fields for the
!  background error covariance using the method of Mogensen et al.
!  (2012), which assumes that the background errors are proportional
!  to the vertical derivatives of the state. The error have simmilar
!  shape as the prior profile but iwth a displacement over the water
!  column.
!
      IF (my_outer.eq.0) THEN
        lstd=1                          ! time index in error arrays
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL background_std (ng, tile, tindex, lstd)
          END DO
        END DO
!
!  Write out estimated background standard deviation fields.
!
        DO ng=1,ngrids
          IF (lwrtstd(ng)) THEN
            std(5,ng)%Rindex=0
            fcount=std(5,ng)%load
            std(5,ng)%Nrec(fcount)=1
#  ifdef DISTRIBUTE
            CALL wrt_std (ng, myrank, lstd)
#  else
            CALL wrt_std (ng, -1, lstd)
#  endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END IF
        END DO
      END IF
# endif
!
!  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), my_outer
        lstr=len_trim(his(ng)%name)
        his(ng)%base=his(ng)%name(1:lstr-3)
      END DO
!
!  Set the nonlinear model to output the quicksave history file as a
!  function of the outer loop. It may be used as the basic state
!  trajectory for the surface fluxes (wind stress, shortwave, heat
!  flux, and E-P) because they can be saved frequently to resolve
!  the daily cycle while avoiding large files.
!
      DO ng=1,ngrids
        ldefqck(ng)=.true.
        lwrtqck(ng)=.true.
        WRITE (qck(ng)%name,10) trim(qck(ng)%head), my_outer
        lstr=len_trim(qck(ng)%name)
        qck(ng)%base=qck(ng)%name(1:lstr-3)
      END DO
!
!  Define output 4D-Var NetCDF file (DAV struc) 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
 
# ifdef PROFILE
!
!  Stop profile clock
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_off (ng, inlm, 86, __line__, myfile)
        END DO
      END DO
# endif
!
  10  FORMAT (a,'_outer',i0,'.nc')
!
      RETURN
      END SUBROUTINE background_initialize
!
      SUBROUTINE background (my_outer, RunInterval)
!
!=======================================================================
!                                                                      !
!  This routine computes the backgound state trajectory, Xb_n-1(t),    !
!  used to linearize the tangent linear and adjoint models in the      !
!  inner loops. It interpolates the background at the observations     !
!  locations, and computes the accept/reject quality control flag,     !
!  ObsScale.                                                           !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     my_outer        Outer-loop counter (integer)                     !
!     RunInterval     NLM kernel time stepping window (seconds)        !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      integer, intent(in) :: my_outer
!
      real(dp), intent(in) :: runinterval
!
!  Local variable declarations.
!
      logical :: donestepping
!
      integer :: i, lstr, ng, tile
      integer :: fcount, tindex
# ifdef PROFILE
      integer :: thread
# endif
# if defined MODEL_COUPLING && !defined MCT_LIB
      integer :: nstrstep, nendstep, extra
!
      real(dp) :: endtime, nexttime
# endif
!
      character (len=20) :: string
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", background"
!
      sourcefile=myfile
 
# if !(defined MODEL_COUPLING && defined ESMF_LIB)
!
!-----------------------------------------------------------------------
!  Set nonlinear model initial conditions.
!-----------------------------------------------------------------------
!
      CALL background_initialize (my_outer)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
!
!-----------------------------------------------------------------------
!  Run nonlinear model and compute background state trajectory,
!  Xb_n-1(t), and the background values at the observation points
!  and times. It processes and writes the observations accept/reject
!  flag (ObsScale) once to allow background quality control, if any.
# if defined MODEL_COUPLING && !defined MCT_LIB
!  Since the ROMS kernel has a delayed output and line diagnostics by
!  one timestep, subtact an extra value to the report of starting and
!  ending timestep for clarity. Usually, the model coupling interval
!  is of the same size as ROMS timestep.
# endif
!-----------------------------------------------------------------------
 
# ifdef PROFILE
!
!  Start profile clock.
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_on (ng, inlm, 86, __line__, myfile)
        END DO
      END DO
# endif
!
!  Initialize various parameters and switches.
!
      myruninterval=runinterval
!
      DO ng=1,ngrids
# ifdef AVERAGES
        ldefavg(ng)=.true.
        lwrtavg(ng)=.true.
        WRITE (avg(ng)%name,10) trim(avg(ng)%head), my_outer
        lstr=len_trim(avg(ng)%name)
        avg(ng)%base=avg(ng)%name(1:lstr-3)
# endif
# ifdef DIAGNOSTICS
        ldefdia(ng)=.true.
        lwrtdia(ng)=.true.
        WRITE (dia(ng)%name,10) trim(dia(ng)%head), my_outer
        lstr=len_trim(dia(ng)%name)
        dia(ng)%base=dia(ng)%name(1:lstr-3)
# endif
        wrtmisfit(ng)=.true.
        wrtobsscale(ng)=.true.
        sporadicimpulse(ng)=.false.
        frequentimpulse(ng)=.false.
# if defined MODEL_COUPLING && !defined MCT_LIB
!
        nexttime=time(ng)+runinterval
        endtime=initime(ng)+(ntimes(ng)-1)*dt(ng)
        IF ((nexttime.eq.endtime).and.(ng.eq.1)) THEN
          extra=0                                   ! last time interval
        ELSE
          extra=1
        END IF
        step_counter(ng)=0
        nstrstep=iic(ng)
        nendstep=nstrstep+int((myruninterval)/dt(ng))-extra
        IF (master) WRITE (stdout,20) 'NL', ng, my_outer, 0,            &
     &                                nstrstep, nendstep
# else
 
        IF (master) WRITE (stdout,20) 'NL', ng, my_outer, 0,            &
     &                                ntstart(ng), ntend(ng)
# endif
      END DO
!
!  Run nonlinear model kernel.
!
# ifdef SOLVE3D
      CALL main3d (runinterval)
# else
      CALL main2d (runinterval)
# endif
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!-----------------------------------------------------------------------
!  If completed stepping, reset switches and write out cost function.
# if defined MODEL_COUPLING && !defined MCT_LIB
!  In coupled applications, RunInterval is much less than ntimes*dt,
!  so we need to wait until the last coupling interval is finished.
!  Otherwise, the control switches will be turned off prematurely.
# endif
!-----------------------------------------------------------------------
!
# if defined MODEL_COUPLING && !defined MCT_LIB
      IF (nendstep.eq.ntend(1)) THEN
        donestepping=.true.
      ELSE
        donestepping=.false.
      END IF
# else
      donestepping=.true.
# endif
!
      IF (donestepping) THEN
        DO ng=1,ngrids
          ldefqck(ng)=.false.
          lwrtqck(ng)=.false.
# ifdef AVERAGES
          ldefavg(ng)=.false.
          lwrtavg(ng)=.false.
# endif
# ifdef DIAGNOSTICS
          ldefdia(ng)=.false.
          lwrtdia(ng)=.false.
# endif
          wrtnlmod(ng)=.false.
          wrtobsscale(ng)=.false.
        END DO
!
!  Report data penalty function.
!
        DO ng=1,ngrids
          IF (master) THEN
            DO i=0,nobsvar(ng)
              IF (i.eq.0) THEN
                string='Total'
              ELSE
                string=obsname(i)
              END IF
              IF (fourdvar(ng)%NLPenalty(i).ne.0.0_r8) THEN
                WRITE (stdout,30) my_outer, 0, 'NLM',                   &
     &                            fourdvar(ng)%NLPenalty(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, inlm, dav(ng)%name,             &
     &                              'NL_iDataPenalty',                  &
     &                              fourdvar(ng)%NLPenalty(0:),         &
     &                              (/1/), (/nobsvar(ng)+1/),           &
     &                              ncid = dav(ng)%ncid)
 
# if defined PIO_LIB && defined DISTRIBUTE
            CASE (io_pio)
              CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name,         &
     &                                  'NL_iDataPenalty',              &
     &                                  fourdvar(ng)%NLPenalty(0:),     &
     &                                  (/1/), (/nobsvar(ng)+1/),       &
     &                                  piofile = dav(ng)%pioFile)
# endif
          END SELECT
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Clean penalty array before next run of NL model.
!
          fourdvar(ng)%NLPenalty=0.0_r8
        END DO
      END IF
 
# ifdef PROFILE
!
!  Stop profile clock
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_off (ng, inlm, 86, __line__, myfile)
        END DO
      END DO
# endif
!
 10   FORMAT (a,'_outer',i0,'.nc')
 20   FORMAT (/,1x,a,1x,'ROMS: started time-stepping:',                 &
     &        ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3,          &
              ', TimeSteps: ',i0,' - ',i0,')',/)
 30   FORMAT (' (',i3.3,',',i3.3,'): ',a,' data penalty, Jdata = ',     &
     &        1p,e17.10,0p,t68,a)
 
      RETURN
      END SUBROUTINE background
!
      SUBROUTINE increment (my_outer, RunInterval)
!
!=======================================================================
!                                                                      !
!  This routine computes the 4D-Var data assimilation state increment, !
!  dXa, by iterating the inner loops and minimizing the cost function. !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     my_outer        Outer-loop counter (integer)                     !
!     RunInterval     TLM/ADM kernels time stepping window (seconds)   !
!                                                                      !
!=======================================================================
 
# ifdef RPCG
!
      USE comp_jb0_mod,   ONLY : aug_oper, comp_jb0
      USE ini_adjust_mod, ONLY : ini_adjust
      USE sum_grad_mod,   ONLY : sum_grad
      USE sum_imp_mod,    ONLY : sum_imp
# endif
!
!  Imported variable declarations
!
      integer, intent(in) :: my_outer
!
      real(dp), intent(in) :: runinterval
!
!  Local variable declarations.
!
      logical :: lcgini, linner, lposterior
!
      integer :: i, ifile, lstr, my_inner, ng, tile
      integer :: fcount, inprec
# ifdef RPCG
      integer :: adrec, nlast
      integer :: irec, jrec, jrec1, jrec2
      integer :: linl
# endif
# ifdef PROFILE
      integer :: thread
# endif
!
      integer, dimension(Ngrids) :: nrec
# ifdef RPCG
      integer, dimension(Ngrids) :: nadrec
# endif
# ifdef SPLIT_4DVAR
!
      real(dp) :: stime
# endif
!
      character (len=8 ) :: driver = 'rbl4dvar'
      character (len=10) :: suffix
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", increment"
!
      sourcefile=myfile
!
!=======================================================================
!  Compute 4D-Var increment.
!=======================================================================
 
# ifdef PROFILE
!
!  Start profile clock.
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_on (ng, itlm, 87, __line__, myfile)
        END DO
      END DO
# endif
 
# ifdef SPLIT_4DVAR
!
!-----------------------------------------------------------------------
!  If split 4D-Var algorithm, set several variables that are computed
!  or assigned in other 4D-Var phase executable.
!-----------------------------------------------------------------------
!
!  Reset Nrun counter to a value greater than one.
!
      IF (my_outer.gt.1) THEN
        nrun=1+(my_outer-1)*ninner
      END IF
!
!  Open Nonlinear model initial conditions NetCDF file (INI struc) and
!  inquire about its variables IDs.
!
      DO ng=1,ngrids
        ldefini(ng)=.false.
        CALL def_ini (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  If outer>1, open tangent linear initial conditions NetCDF file
!  (ITL struc) and inquire about its variables IDs.
!
      IF (my_outer.gt.1) THEN
        DO ng=1,ngrids
          ldefitl(ng)=.false.
          CALL tl_def_ini (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END IF
!
!  If outer>1, open impulse forcing NetCDF file (TLF struc) and inquire
!  about its variables IDs.
!
      IF (my_outer.gt.1) THEN
        DO ng=1,ngrids
          ldeftlf(ng)=.false.
          CALL def_impulse (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END IF
!
!  If outer>1, open adjoint history NetCDF file (ADM struc) and inquire
!  about its variables IDs and set "FrcRec".
!
      IF (my_outer.gt.1) THEN
        DO ng=1,ngrids
          ldefadj(ng)=.false.
          CALL ad_def_his (ng, ldefadj(ng))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          fcount=adm(ng)%Fcount
          frcrec(ng)=adm(ng)%Nrec(fcount)
        END DO
      END IF
!
!  Open 4D-Var NetCDF file (DAV struc) and inquire about its variables.
!  Deactivate "haveNLmod" since its values are read below. Otherwise,
!  the TLM will read zero values when calling "obs_read" for outer>1.
!  The DAV file does not have values for NLmodel_value(:,outer) yet.
!
      DO ng=1,ngrids
        ldefmod(ng)=.false.
        CALL def_mod (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        havenlmod(ng)=.false.        ! because is activated in "def_mod"
      END DO
!
!  If split 4D-Var and outer>1, read several variables from 4D-VAR
!  NetCDF (DAV struc) file needed in the conjugate gradient algorithm.
!  In the unsplit case, such values are available in memory.
!
      IF (my_outer.gt.1) THEN
        DO ng=1,ngrids
#  ifdef RPCG
          CALL cg_read_rpcg (ng, itlm, my_outer)
#  else
          CALL cg_read_congrad (ng, itlm, my_outer)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END IF
!
!  In the split 4D-Var algorithm, we need to read in initial NLmodVal
!  ("NLmodel_initial") and ObsScale ("obs_scale"), computed in the
!  background phase, from the DAV NetCDF file. Such data is in memory
!  in the unsplit algorithm.
!
!  HGA: What to do in 4D-Var with nested grids?
!
      IF (my_outer.eq.1) THEN        ! initial values from "background"
        ng=1
        SELECT CASE (dav(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_get_fvar (ng, itlm, dav(ng)%name,               &
     &                            vname(1,idnlmi), nlmodval,            &
     &                            ncid = dav(ng)%ncid)
 
#  if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name,           &
     &                                vname(1,idnlmi), nlmodval,        &
     &                                piofile = dav(ng)%pioFile)
 
#  endif
        END SELECT
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
      ng=1
      SELECT CASE (dav(ng)%IOtype)
        CASE (io_nf90)
          CALL netcdf_get_fvar (ng, itlm, dav(ng)%name,                 &
     &                          vname(1,idobss), obsscale,              &
     &                          ncid = dav(ng)%ncid)
 
#  if defined PIO_LIB && defined DISTRIBUTE
        CASE (io_pio)
          CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name,             &
     &                              vname(1,idobss), obsscale,          &
     &                              piofile = dav(ng)%pioFile)
#  endif
      END SELECT
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  In the split 4D-Var alorithm, we also need to read in ObsVal
!  ("obs_value") and ObsErr ("obs_error") from the observation file.
!  They are used in the initialization of the conjugate gradient
!  algorithm.
!
      ng=1
      SELECT CASE (obs(ng)%IOtype)
        CASE (io_nf90)
          CALL netcdf_get_fvar (ng, itlm, obs(ng)%name,                 &
     &                          vname(1,idoval), obsval)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          CALL netcdf_get_fvar (ng, itlm, obs(ng)%name,                 &
     &                          vname(1,idoerr), obserr)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#  if defined PIO_LIB && defined DISTRIBUTE
        CASE (io_pio)
          CALL pio_netcdf_get_fvar (ng, itlm, obs(ng)%name,             &
     &                              vname(1,idoval), obsval)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
          CALL pio_netcdf_get_fvar (ng, itlm, obs(ng)%name,             &
     &                              vname(1,idoerr), obserr)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
      END SELECT
!
!  Set nonlinear output history file to be used as the basic state
!  trajectory. The 4D-Var increment phase is computed by a different
!  executable and needs to know some of the HIS structure information.
!
      DO ng=1,ngrids
        WRITE (his(ng)%name,10) trim(fwd(ng)%head), my_outer-1
        lstr=len_trim(his(ng)%name)
        his(ng)%base=his(ng)%name(1:lstr-3)
        IF (his(ng)%Nfiles.gt.1) THEN
          DO ifile=1,his(ng)%Nfiles
            WRITE (suffix,"('_',i4.4,'.nc')") ifile
            his(ng)%files(ifile)=trim(his(ng)%base)//trim(suffix)
          END DO
          his(ng)%name=trim(his(ng)%files(1))
        ELSE
          his(ng)%files(1)=trim(his(ng)%name)
        END IF
      END DO
!
!  Set the nonlinear model quicksave-history file as the basic state for
!  the surface fluxes computed in "bulk_flux", which may be available at
!  more frequent intervals while avoiding large files. Since the 4D-Var
!  increment phase is calculated by a different executable and needs to
!  know some of the QCK structure information.
!
      DO ng=1,ngrids
        WRITE (qck(ng)%name,10) trim(qck(ng)%head), my_outer-1
        lstr=len_trim(qck(ng)%name)
        qck(ng)%base=qck(ng)%name(1:lstr-3)
        IF (qck(ng)%Nfiles.gt.1) THEN
          DO ifile=1,qck(ng)%Nfiles
            WRITE (suffix,"('_',i4.4,'.nc')") ifile
            qck(ng)%files(ifile)=trim(qck(ng)%base)//trim(suffix)
          END DO
          qck(ng)%name=trim(qck(ng)%files(1))
        ELSE
          qck(ng)%files(1)=trim(qck(ng)%name)
        END IF
      END DO
!
!  Read in 4D-Var starting time (sec) from nonlinear trajectory.
!  Initialize "tday" which are needed to write the correct time in
!  the ITL NetCDF file.  It is alse need for boundary and surface
!  forcing adjustments, if any.
!
      inprec=1
      DO ng=1,ngrids
        SELECT CASE (his(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_get_fvar (ng, itlm, his(ng)%name,               &
     &                            vname(1,idtime), stime,               &
     &                            start = (/inprec/),                   &
     &                            total = (/1/))
 
#  if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL pio_netcdf_get_fvar (ng, itlm, his(ng)%name,           &
     &                                vname(1,idtime), stime,           &
     &                                start = (/inprec/),               &
     &                                total = (/1/))
#  endif
        END SELECT
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        initime(ng)=stime
 
#  ifdef ADJUST_BOUNDARY
!
!  Set time (sec) for the open boundary adjustment.
!
        obc_time(1,ng)=stime
        DO i=2,nbrec(ng)
          obc_time(i,ng)=obc_time(i-1,ng)+nobc(ng)*dt(ng)
        END DO
#  endif
 
#  if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
!
!  Set time (sec) for the surface forcing adjustment.
!
        sf_time(1,ng)=stime
        DO i=2,nfrec(ng)
          sf_time(i,ng)=sf_time(i-1,ng)+nsff(ng)*dt(ng)
        END DO
#  endif
      END DO
# endif
!
!-----------------------------------------------------------------------
!  Set few variables needed in the "increment" phase.
!-----------------------------------------------------------------------
!
!  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.
!
      CALL edit_multifile ('HIS2FWD')
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      DO ng=1,ngrids
        lreadfwd(ng)=.true.
      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 RPCG
!
!  If weak constraint, set the number of TLF records used in the
!  augmented solver.
!
      DO ng=1,ngrids
        IF (nadj(ng).lt.ntimes(ng)) THEN
          nadrec(ng)=2*(1+ntimes(ng)/nadj(ng))
        ELSE
          nadrec(ng)=0
        END IF
      END DO
# endif
!
!-----------------------------------------------------------------------
!  On first pass (outer=1), create NetCDF files and initialize all the
!  records in the ITL file to zero.
!-----------------------------------------------------------------------
!
      check_outer1 : IF (my_outer.eq.1) THEN
!
!  If outer=1, define tangent linear initial conditions file.
!
        DO ng=1,ngrids
          ldefitl(ng)=.true.
          CALL tl_def_ini (ng)
          ldefitl(ng)=.false.
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
 
# ifdef RPCG
!
!  If outer=1, Initialize all records of the ITL file to zero. The TLM
!  state variables are zero when outer=1.
!
        DO ng=1,ngrids
          tdays(ng)=initime(ng)*sec2day
#  ifdef DISTRIBUTE
          CALL tl_wrt_ini (ng, myrank, rec1, rec1)
#  else
          CALL tl_wrt_ini (ng, -1, rec1, rec1)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  ifdef DISTRIBUTE
          CALL tl_wrt_ini (ng, myrank, rec1, rec2)
#  else
          CALL tl_wrt_ini (ng, -1, rec1, rec2)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  ifdef DISTRIBUTE
          CALL tl_wrt_ini (ng, myrank, rec1, rec3)
#  else
          CALL tl_wrt_ini (ng, -1, rec1, rec3)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  ifdef DISTRIBUTE
          CALL tl_wrt_ini (ng, myrank, rec1, rec4)
#  else
          CALL tl_wrt_ini (ng, -1, rec1, rec4)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  ifdef DISTRIBUTE
          CALL tl_wrt_ini (ng, myrank, rec1, rec5)
#  else
          CALL tl_wrt_ini (ng, -1, rec1, rec5)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          IF (nadj(ng).lt.ntimes(ng)) THEN
            nlast=rec5
            DO irec=1,nadrec(ng)
#  ifdef DISTRIBUTE
              CALL tl_wrt_ini (ng, myrank, rec1, nlast+irec)
#  else
              CALL tl_wrt_ini (ng, -1, rec1, nlast+irec)
#  endif
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
          END IF
        END DO
# endif
!
!  If outer=1, define 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
 
# if defined POSTERIOR_EOFS    || defined POSTERIOR_ERROR_I || \
     defined posterior_error_f
!
!  If outer=1, define output Hessian NetCDF file that will eventually
!  contain the intermediate posterior analysis error covariance matrix
!  fields or its EOFs.
!
        DO ng=1,ngrids
          ldefhss(ng)=.true.
          CALL def_hessian (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
# endif
 
# if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F
!
!  If outer=1, define output initial or final full posterior error
!  covariance (diagonal) matrix NetCDF.
!
        DO ng=1,ngrids
          ldeferr(ng)=.true.
          CALL def_error (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
# endif
 
      END IF check_outer1
 
# ifdef RPCG
!
!-----------------------------------------------------------------------
!  If outer>1, integrate the tangent linear model to evaluate
!  G[x(k)-x(k-1)], which is written to the DAV NetCDF in the
!  variable "TLmodel_value". It is needed in routine "rpcg_lanczos".
!  In the old unsplit algorithm, the TLM was run in the "analysis"
!  phase. It is done here to allow lower resolution in the inner
!  loops in the future.
!-----------------------------------------------------------------------
!
      check_outer2 : IF (my_outer.gt.1) THEN
!
!  Clear tangent linear forcing arrays. This is very important since
!  these arrays are non-zero and must be zero when running the tangent
!  linear model. We also need to clear the nonlinear state arrays to
!  make sure that the RHS terms rzeta, rubar, rvbar, ru, and rv are
!  zero.  Otherwise, those array will have the last computed values
!  when running the nonlinear model if not processing the forward
!  trajectory RHS terms.  This needs to be done to get identical
!  solutions with the split schemes.
!
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL initialize_ocean (ng, tile, inlm)
            CALL initialize_forces (ng, tile, itlm)
#  ifdef ADJUST_BOUNDARY
            CALL initialize_boundary (ng, tile, itlm)
#  endif
          END DO
        END DO
!
!  Initialize tangent linear model from initial impulse which is now
!  stored in file ITL(ng)%name.
!
        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 ITLname, record Rec3.
!  The sum of the initial condition increments is in record 3.
!
        DO ng=1,ngrids
          itl(ng)%Rindex=rec3
          CALL tl_initial (ng)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
!  Run tangent linear model.
!
        DO ng=1,ngrids
          IF (master) THEN
            WRITE (stdout,20) 'TL', ng, my_outer, 0,                    &
     &                        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 IF check_outer2
# endif
!
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!  4D-Var inner loops.
!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
!  Clear tangent linear forcing arrays before entering inner-loop.
!  This is very important since these arrays are non-zero 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_ocean (ng, tile, inlm)
            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 loops iteration.
!
      inner_loop : DO my_inner=0,ninner
        inner=my_inner
 
# ifdef RPCG
        IF (inner.ne.ninner) THEN
          linner=.true.
        ELSE
          linner=.false.
        END IF
!
!  Retrieve TLmodVal and NLmodVal when inner=0 and outer>1 for use as
!  Hbk and BCKmodVal, respectively.
!
        IF ((inner.eq.0).and.(my_outer.gt.1)) THEN
          DO ng=1,ngrids
            SELECT CASE (dav(ng)%IOtype)
              CASE (io_nf90)
                CALL netcdf_get_fvar (ng, itlm, dav(ng)%name,           &
     &                                'TLmodel_value', tlmodval)
                IF (founderror(exit_flag, noerror, __line__,            &
     &                         myfile)) RETURN
 
                CALL netcdf_get_fvar (ng, itlm, dav(ng)%name,           &
     &                                'NLmodel_value', nlmodval,        &
     &                                ncid = dav(ng)%ncid,              &
     &                                start = (/1,my_outer-1/),         &
     &                                total = (/ndatum(ng),1/))
                IF (founderror(exit_flag, noerror, __line__,            &
     &                         myfile)) RETURN
 
#  if defined PIO_LIB && defined DISTRIBUTE
              CASE (io_pio)
                CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name,       &
     &                                    'TLmodel_value', tlmodval)
                IF (founderror(exit_flag, noerror, __line__,            &
     &                         myfile)) RETURN
 
                CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name,       &
     &                                    'NLmodel_value', nlmodval,    &
     &                                    piofile = dav(ng)%pioFile,    &
     &                                    start = (/1,my_outer-1/),     &
     &                                    total = (/ndatum(ng),1/))
                IF (founderror(exit_flag, noerror, __line__,            &
     &                         myfile)) RETURN
#  endif
            END SELECT
          END DO
        END IF
!
        IF (inner.eq.0) THEN
          lcgini=.true.
        END IF
        DO ng=1,ngrids
          CALL rpcg_lanczos (ng, itlm, my_outer, inner, ninner, lcgini)
        END DO
# else
!
!  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, my_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
# endif
!
!  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
            CALL ad_initial (ng)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
            wrtmisfit(ng)=.false.
          END DO
!
!  Set adjoint history NetCDF parameters.  Define adjoint history
!  file only once to avoid opening too many files.
!
          DO ng=1,ngrids
# ifdef WEAK_NOINTERP
            wrtforce(ng)=.false.
# else
            wrtforce(ng)=.true.
# endif
            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, my_outer, my_inner,           &
     &                          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
!
!  Activate "LwrtState2d" switch to write the correct ad_ubar, ad_vbar,
!  ad_zeta here instead of ad_ubar_sol, ad_vbar_sol, and ad_zeta_sol in
!  the calls to "ad_wrt_his" below (HGA).
!  Here, ad_zeta(:,:,kout)=ad_zeta_sol. However, ad_ubar and ad_vbar
!  not equal to ad_ubar_sol and ad_vbar_sol, respectively. It is
!  irrelevant because ubar and vbar are not part of the state in
!  3D application.  Notice that "LwrtState2d" will be turned off at
!  the bottom of "error_covariance".
!
          DO ng=1,ngrids
            lwrtstate2d(ng)=.true.
          END DO
 
# ifndef WEAK_NOINTERP
!
!  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
# endif
!
!  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.
!
# ifdef POSTERIOR_ERROR_I
          lposterior=.true.
# else
          lposterior=.false.
# endif
# ifdef RPCG
!
!  Set the flag that controls the augmentation of the model error
!  forcing terms. This is ONLY done in the outer-loop so
!  LaugWeak=.FALSE. here.
!
          laugweak=.false.
# endif
          CALL error_covariance (itlm, driver, my_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,30) my_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.
!-----------------------------------------------------------------------
!
!  Initialize tangent linear model from initial impulse which is now
!  stored in file ITL(ng)%name.
!
          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 ((my_outer.eq.1).and.(inner.eq.1)) THEN
            DO ng=1,ngrids
              wrtmisfit(ng)=.true.
            END DO
          END IF
!
!  Run tangent linear model forward and force with convolved adjoint
!  trajectory impulses. Compute (H M B M' 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, my_outer, my_inner,           &
     &                          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 POSTERIOR_ERROR_F
!
!  Copy the final time tl_var(Lold) into ad_var(Lold) so that it can be
!  written to the Hessian NetCDF file.
!
          add=.false.
          DO ng=1,ngrids
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL load_tltoad (ng, tile, lold(ng), lold(ng), add)
            END DO
          END DO
!
!  Write evolved tangent solution into hessian netcdf file for use
!  later.
!
          IF (inner.ne.0) THEN
            DO ng=1,ngrids
#  ifdef DISTRIBUTE
              CALL wrt_hessian (ng, tile, lold(ng), lold(ng))
#  else
              CALL wrt_hessian (ng, -1, lold(ng), lold(ng))
#  endif
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
          END IF
# endif
!
          nrun=nrun+1
 
# ifndef RPCG
!
!  Use conjugate gradient algorithm to find a better approximation
!  PSI to coefficients Beta_n.
!
          DO ng=1,ngrids
            lcgini=.false.
            CALL congrad (ng, itlm, my_outer, inner, ninner, lcgini)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
# else
          lcgini=.false.
# endif
 
        END IF inner_compute
 
      END DO inner_loop
!
!-----------------------------------------------------------------------
!  Once the w_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 nonlinear model.
!-----------------------------------------------------------------------
!
!  Initialize the adjoint model always from rest.
!
      DO ng=1,ngrids
        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
# ifdef WEAK_NOINTERP
        wrtforce(ng)=.false.
# else
        wrtforce(ng)=.true.
# endif
        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 coefficients,
!  Beta_n.
!
      DO ng=1,ngrids
        IF (master) THEN
          WRITE (stdout,20) 'AD', ng, my_outer, ninner,                 &
     &                      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
 
# ifndef WEAK_NOINTERP
!
!  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
# endif
!
!  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
 
# ifdef RPCG
!
!  Get number of records in adjoint NetCDF.
!  We need to do this here since ADM(ng)%Nrec is reset to zero in
!  error_covariance.
!
      DO ng=1,ngrids
        fcount=adm(ng)%load
        nrec(ng)=adm(ng)%Nrec(fcount)
      END DO
# endif
!
!  Convolve adjoint trajectory with error covariances.
!
      lposterior=.false.
 
# ifdef RPCG
!
!  Set the flag that controls the augmentation of the model error
!  forcing terms. This is ONLY done in the outer-loop so
!  LaugWeak=.TRUE. here.
!
      laugweak=.true.
      linl=my_outer+1
      CALL error_covariance (inlm, driver, my_outer, inner,             &
     &                       linl, lini, lold, lnew,                    &
     &                       rec1, rec2, lposterior)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# else
      CALL error_covariance (inlm, driver, my_outer, inner,             &
     &                       lbck, lini, lold, lnew,                    &
     &                       rec1, rec2, lposterior)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
 
# ifdef RPCG
!
!  Augmented solver:
!
!  NOTES: The ITL file contains 5 records -
!         Rec2 = the new TL initial condition
!         Rec3 = the sum of the TL initial conditions
!         Rec4 = B^-1(xb-xk)=sum_j^k-1 G_j lambda_j
!         Rec5 = the augmented correction to Rec2.
!         Rec5+1 to Rec5+nADrec/2 = sum of the TLF forcing increments
!         Rec5+nADrec/2+1 to Rec5+nADrec = B^-1 of the sum of the TLF
!                                         forcing increments.
!  Reset the flag LaugWeak flag.
!
      laugweak=.false.
 
# else
!
!  Write out nonlinear model initial conditions into INIname, record
!  INI(ng)%Rindex.
!
      DO ng=1,ngrids
#  ifdef DISTRIBUTE
        CALL wrt_ini (ng, myrank, lnew(ng))
#  else
        CALL wrt_ini (ng, -1, lnew(ng))
#  endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#  if defined ADJUST_STFLUX   || defined ADJUST_WSTRESS || \
      defined adjust_boundary
#   ifdef DISTRIBUTE
        CALL wrt_frc_ad (ng, myrank, lold(ng), ini(ng)%Rindex)
#   else
        CALL wrt_frc_ad (ng, -1, lold(ng), ini(ng)%Rindex)
#   endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
      END DO
# endif
 
# ifdef RPCG
!
! Compute the augmented correction to the adjoint propagator.
! We need to use sum (x(k)-x(k-1)) before it is updated. This is
! in record 3 of the ITL file.
!
      DO ng=1,ngrids
        CALL get_state (ng, itlm, 4, itl(ng), rec3, ltlm1)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL aug_oper (ng, tile, ltlm1, ltlm2)
        END DO
      END DO
!
! Save this augmented piece in record 5 of the ITL file.
!
      DO ng=1,ngrids
#  ifdef DISTRIBUTE
        CALL tl_wrt_ini (ng, myrank, ltlm2, rec5)
#  else
        CALL tl_wrt_ini (ng, -1, ltlm2, rec5)
#  endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
! Complete the computation of the TL initial condition by adding the
! contribution from the augmented adjoint propagator.
!
      DO ng=1,ngrids
        CALL get_state (ng, itlm, 4, itl(ng), rec2, ltlm1)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        CALL get_state (ng, itlm, 4, itl(ng), rec5, ltlm2)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL sum_grad (ng, tile, ltlm1, ltlm2)
        END DO
      END DO
!
! Write the final TL increment to Rec2 of the ITL file.
!
      DO ng=1,ngrids
#  ifdef DISTRIBUTE
        CALL tl_wrt_ini (ng, myrank, ltlm2, rec2)
#  else
        CALL tl_wrt_ini (ng, -1, ltlm2, rec2)
#  endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
! Now update the non-linear model initial conditions.
!
      linl=my_outer+1
      DO ng=1,ngrids
        CALL get_state (ng, inlm, 9, ini(ng), linl, lnew(ng))
        IF (timeiau(ng).eq.0.0_dp) THEN
          CALL get_state (ng, iadm, 4, itl(ng), rec2, ltlm1)
        END IF
      END DO
      DO ng=1,ngrids
        IF (timeiau(ng).eq.0.0_dp) THEN
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL ini_adjust (ng, tile, ltlm1, lnew(ng))
          END DO
        END IF
      END DO
!
!  Write out nonlinear model initial conditions into INIname, record
!  INI(ng)%Rindex.
!
      DO ng=1,ngrids
#  ifdef DISTRIBUTE
        CALL wrt_ini (ng, myrank, lnew(ng))
#  else
        CALL wrt_ini (ng, -1, lnew(ng))
#  endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
#  if defined ADJUST_STFLUX   || defined ADJUST_WSTRESS || \
      defined adjust_boundary
#   ifdef DISTRIBUTE
        CALL wrt_frc_ad (ng, myrank, ltlm1, ini(ng)%Rindex)
#   else
        CALL wrt_frc_ad (ng, -1, ltlm1, ini(ng)%Rindex)
#   endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
      END DO
!
! Compute the new B^-1(x(k)-x(k-1)) term.
! Gather the final adjoint solutions increments, sum and
! save in record 4 of the ITL file. Use the tl arrays as temporary
! storage.
!
! First add the augmented piece which is computed from the previous
! sum.
!
      DO ng=1,ngrids
        CALL get_state (ng, itlm, 4, itl(ng), rec4, ltlm1)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL aug_oper (ng, tile, ltlm1, ltlm2)
        END DO
      END DO
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL sum_grad (ng, tile, ltlm1, ltlm2)
        END DO
      END DO
!
      DO ng=1,ngrids
        adrec=nrec(ng)
        CALL get_state (ng, itlm, 4, adm(ng), adrec, ltlm1)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL sum_grad (ng, tile, ltlm1, ltlm2)
        END DO
      END DO
!
!  Write the current sum of adjoint solutions into record 4 of the ITL
!  file.
!
      DO ng=1,ngrids
#  ifdef DISTRIBUTE
        CALL tl_wrt_ini (ng, myrank, ltlm2, rec4)
#  else
        CALL tl_wrt_ini (ng, -1, ltlm2, rec4)
#  endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
# endif
!
!  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,30) my_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
 
# ifdef RPCG
!
!  Now Compute the augmented corrections to the weak constraint
!  forcing terms. The sums so far are in records 6 to
!  6+nADrec/2.
!
      DO ng=1,ngrids
        DO i=1,nadrec(ng)/2
          irec=i
          jrec=rec5+i
          CALL get_state (ng, itlm, 4, itl(ng), jrec, ltlm1)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL aug_oper (ng, tile, ltlm1, ltlm2)
          END DO
!
!  Complete the computation of the TLF forcing term by adding the
!  contribution from the augmented adjoint propagator. Specify
!  the value of 7 for the model variable since this the special
!  case in get_state for reading the impulse forcing.
!
!!        CALL get_state (ng, iTLM, 4, TLF(ng), irec, LTLM1)      ! TEST
          CALL get_state (ng, 7, 4, tlf(ng), irec, ltlm1)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          DO tile=first_tile(ng),last_tile(ng),+1
!!          CALL sum_grad (ng, tile, LTLM1, LTLM2)                ! TEST
            CALL sum_imp (ng, tile, ltlm2)
          END DO
!
!  Write the final forcing increment to he TLF file.
!  Note the original TLF file is overwritten at this point.
!
          tlf(ng)%Rindex=0
# ifdef DISTRIBUTE
          CALL wrt_aug_imp (ng, myrank, itlm, ltlm2, i)
# else
          CALL wrt_aug_imp (ng, -1, itlm, ltlm2, i)
# endif
 
        END DO
      END DO
!
!  Gather the increments from the final inner-loop and
!  save in record 3 of the ITL file. The current increment
!  is in record 2 and the sum so far is in record 3.
!
      DO ng=1,ngrids
        CALL get_state (ng, itlm, 4, itl(ng), rec2, ltlm1)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        CALL get_state (ng, itlm, 4, itl(ng), rec3, ltlm2)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL sum_grad (ng, tile, ltlm1, ltlm2)
        END DO
      END DO
!
!  Write the current sum into record 3 of the ITL file.
!
      DO ng=1,ngrids
# ifdef DISTRIBUTE
        CALL tl_wrt_ini (ng, myrank, ltlm2, rec3)
# else
        CALL tl_wrt_ini (ng, -1, ltlm2, rec3)
# endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
 
# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
     defined adjust_boundary
!
!  Write the current sum of the forcing and boundary increments
!  into the INI file into index INI(ng)%Rindex-1 (i.e. we are
!  overwriting the previous fields. Read the current sum into the
!  adjoint arrays since this is what wrt_frc_AD uses.
!
      DO ng=1,ngrids
        CALL get_state (ng, iadm, 4, itl(ng), rec3, ltlm1)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
      DO ng=1,ngrids
#  ifdef DISTRIBUTE
        CALL wrt_frc_ad (ng, myrank, ltlm1, ini(ng)%Rindex)
#  else
        CALL wrt_frc_ad (ng, -1, ltlm1, ini(ng)%Rindex)
#  endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
# endif
!
!  Gather the model error forcing increments and update the
!  sum in records 6 to 6+nADrec/2.
!
      DO ng=1,ngrids
        DO i=1,nadrec(ng)/2
          irec=i
          jrec=rec5+i
          CALL get_state (ng, itlm, 4, itl(ng), jrec, ltlm1)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!!        CALL get_state (ng, iTLM, 4, TLF(ng)%, irec, LTLM2)     ! TEST
          CALL get_state (ng, 7, 4, tlf(ng), irec, ltlm2)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          DO tile=first_tile(ng),last_tile(ng),+1
!!          CALL sum_grad (ng, tile, LTLM1, LTLM2)                ! TEST
            CALL sum_imp (ng, tile, ltlm1)
          END DO
!
!  Write the current sum into record jrec of the ITL file.
!
# ifdef DISTRIBUTE
!!        CALL tl_wrt_ini (ng, MyRank, LTLM2, jrec)               ! TEST
          CALL tl_wrt_ini (ng, myrank, ltlm1, jrec)
# else
!!        CALL tl_wrt_ini (ng, -1, LTLM2, jrec)                   ! TEST
          CALL tl_wrt_ini (ng, -1, ltlm1, jrec)
# endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END DO
!
!  Now compute the background cost function Jb0.
!  First compute the contribution from the increments in the
!  initial conditions, surface forcing, and boundary conditions.
!
      jb0(my_outer)=0.0_r8
      DO ng=1,ngrids
        CALL get_state (ng, iadm, 8, itl(ng), rec4, ltlm1)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        CALL get_state (ng, itlm, 8, itl(ng), rec3, ltlm1)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL comp_jb0 (ng, tile, itlm, my_outer, ltlm1, ltlm1)
        END DO
      END DO
!
!  Now compute the contribution of model error terms to the
!  background cost function Jb0.
!
!  NOTE: I THINK WE NEED A BETTER WAY OF COMPUTING THE OBS ERROR
!  CONTRIBUTION TO Jb. THE FOLLOWING CALCULATION IS DANGEROUS BECAUSE
!  IT RELIES ON THE FACT THAT THE SURFACING FORCING AND OBC TL ARRAYS
!  READ FROM jrec2 ARE ZERO. THIS MEANS THAT ONLY THE MODEL STATE
!  CONTRIBUTES TO Jb FOR THE MODEL ERROR TERM, WHICH OF COURSE IS
!  WHAT SHOULD BE THE CASE. THE SURFACE FORCING AND OBC CONTRIBUTIONS
!  ARE COMPUTED IN THE PREVIOUS CALL TO COMP_JB0.
!
      DO ng=1,ngrids
        DO irec=1,nadrec(ng)/2
          jrec1=rec5+irec
          jrec2=rec5+nadrec(ng)/2+irec
          CALL get_state (ng, iadm, 8, itl(ng), jrec1, ltlm1)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          CALL get_state (ng, itlm, 8, itl(ng), jrec2, ltlm1)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL comp_jb0 (ng, tile, itlm, my_outer, ltlm1, ltlm1)
          END DO
        END DO
      END DO
!
!  Overwrite the TFL netcdf file with the sum of the model error
!  forcing increments - required for the next run of the NLM and
!  TLM.
!
      DO ng=1,ngrids
        DO irec=1,nadrec(ng)/2
          jrec=rec5+irec
          CALL get_state (ng, itlm, 8, itl(ng), jrec, ltlm1)
!
          tlf(ng)%Rindex=0
# ifdef DISTRIBUTE
          CALL wrt_aug_imp (ng, myrank, itlm, ltlm1, irec)
# else
          CALL wrt_aug_imp (ng, -1, itlm, ltlm1, irec)
# endif
        END DO
      END DO
# endif
 
# ifdef PROFILE
!
!  Stop profile clock
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_off (ng, itlm, 87, __line__, myfile)
        END DO
      END DO
# endif
!
 10   FORMAT (a,'_outer',i0,'.nc')
 20   FORMAT (/,1x,a,1x,'ROMS: started time-stepping:',                 &
     &        ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3,          &
              ', TimeSteps: ',i0,' - ',i0,')',/)
 30   FORMAT (/,' Converting Convolved Adjoint Trajectory to',          &
     &          ' Impulses: Outer = ',i3.3,' Inner = ',i3.3,/)
!
      RETURN
      END SUBROUTINE increment
!
      SUBROUTINE analysis_initialize (my_outer)
!
!=======================================================================
!                                                                      !
!  This routine initializes the nonlinear kernel with the 4D-Var new   !
!  state estimate Xa = Xb + dXa.  It is separated from the 'analysis'  !
!  4D-Var phase in ESM coupling applications that use generic methods  !
!  for 'initialize', 'run', and 'finalize'.                            !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     my_outer        Outer-loop counter (integer)                     !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      integer, intent(in) :: my_outer
!
!  Local variable declarations.
!
      integer :: ifile, lstr, ng, tile
      integer :: fcount
# ifdef PROFILE
      integer :: thread
# endif
!
      integer, dimension(Ngrids) :: indxsave
!
      character (len=10) :: suffix
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", analysis_initialize"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Initialize nonlinear model kernel with new 4D-Var state estimate.
!-----------------------------------------------------------------------
 
# ifdef PROFILE
!
!  Start profile clock.
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_on (ng, inlm, 88, __line__, myfile)
        END DO
      END DO
# endif
!
!  Clear tangent arrays and the nonlinear model mixing arrays
!  before running nonlinear model (important).
!
      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)
# ifdef ADJUST_BOUNDARY
          CALL initialize_boundary (ng, tile, itlm)
# endif
          CALL initialize_mixing (ng, tile, inlm)
        END DO
      END DO
 
# ifdef SPLIT_4DVAR
!
!-----------------------------------------------------------------------
!  If split 4D-Var algorithm, set several variables that computed or
!  assigned in other 4D-Var phase executable.
!-----------------------------------------------------------------------
!
!  Set Nrun>1, to read in surface forcing and open boundary conditions
!  increments in "initial", if appropriate.
!
      nrun=ninner+1
!
!  Set ERstr=Nrun, to set the open boundary condition (OBC_time) and
!  surface forcing (SF_time) adjustment times, if needed.
!
      erstr=nrun
#  if defined MODEL_COUPLING && defined ESMF_LIB
      erend=nouter
#  endif
!
!  Open Nonlinear model initial conditions NetCDF file (INI struc) and
!  inquire about its variable IDs. In particular, set the current value
!  of INI(ng)%Rindex, which is needed in "initial".
!
      DO ng=1,ngrids
        ldefini(ng)=.false.
        CALL def_ini (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Open 4D-Var NetCDF file (DAV struc) and inquire about its variables.
!
      DO ng=1,ngrids
        ldefmod(ng)=.false.
        CALL def_mod (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  In the split 4D-Var algorithm, we need to read the values of ObsScale
!  computed in the background phase from the DAV NetCDF file. Such data
!  is in memory in the unsplit algorithm.
!
!  HGA: What to do in 4D-Var with nested grids?
!
      ng=1
      SELECT CASE (dav(ng)%IOtype)
        CASE (io_nf90)
          CALL netcdf_get_fvar (ng, itlm, dav(ng)%name,                 &
     &                          vname(1,idobss), obsscale,              &
     &                          ncid = dav(ng)%ncid)
 
#  if defined PIO_LIB && defined DISTRIBUTE
        CASE (io_pio)
          CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name,             &
     &                              vname(1,idobss), obsscale,          &
     &                              piofile = dav(ng)%pioFile)
#  endif
      END SELECT
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  If outer=Nouter, read initial NLmodVal ("NLmodel_initial") and load
!  it in the scratch "NLincrement" array.
!
      IF (my_outer.eq.nouter) THEN
        ng=1
        SELECT CASE (dav(ng)%IOtype)
          CASE (io_nf90)
            CALL netcdf_get_fvar (ng, itlm, dav(ng)%name,               &
     &                            vname(1,idnlmi), nlincrement,         &
     &                            ncid = dav(ng)%ncid)
 
#  if defined PIO_LIB && defined DISTRIBUTE
          CASE (io_pio)
            CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name,           &
     &                                vname(1,idnlmi), nlincrement,     &
     &                                piofile = dav(ng)%pioFile)
#  endif
        END SELECT
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
 
#  ifdef FORWARD_FLUXES
!
!  If not first NLM run, set BLK structure to the previous outer loop
!  quicksave trajectory.  There is a logic in "get_data" that reads
!  from BLK.
!  (HGA: This probably legacy code that it is no longer needed)
!
      DO ng=1,ngrids
        WRITE (qck(ng)%name,10) trim(qck(ng)%head), my_outer-1
        lstr=len_trim(qck(ng)%name)
        qck(ng)%base=qck(ng)%name(1:lstr-3)
      END DO
!
      CALL edit_multifile ('QCK2BLK')
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
# endif
!
!-----------------------------------------------------------------------
!  Initialize nonlinear model.
!-----------------------------------------------------------------------
!
!  Set new basic state trajectory for next outer loop.  Notice that the
!  LreadFWD switch is turned off to suppress processing of the structure
!  when "check_multifile" during nonlinear model execution.
!
      DO ng=1,ngrids
        idefhis(ng)=-1
        ldefhis(ng)=.true.
        lwrthis(ng)=.true.
        wrtnlmod(ng)=.true.
        wrttlmod(ng)=.false.
# ifdef FRC_FILES
        lreadfrc(ng)=.true.
# endif
        lreadfwd(ng)=.false.
        WRITE (his(ng)%name,10) trim(fwd(ng)%head), my_outer
        lstr=len_trim(his(ng)%name)
        his(ng)%base=his(ng)%name(1:lstr-3)
        IF (his(ng)%Nfiles.gt.1) THEN
          DO ifile=1,his(ng)%Nfiles
            WRITE (suffix,"('_',i4.4,'.nc')") ifile
            his(ng)%files(ifile)=trim(his(ng)%base)//trim(suffix)
          END DO
          his(ng)%name=trim(his(ng)%files(1))
        ELSE
          his(ng)%files(1)=trim(his(ng)%name)
        END IF
      END DO
!
!  Set nonlinear output quicksave history file for the next outer loop.
!  Notice that the LreadBLK switch is turned off to suppress processing
!  of the structure when "check_multifile" during nonlinear model
!  execution.
!
      DO ng=1,ngrids
        idefqck(ng)=-1
        ldefqck(ng)=.true.
        lwrtqck(ng)=.true.
# ifdef FORWARD_FLUXES
        lreadblk(ng)=.false.
# endif
        WRITE (qck(ng)%name,10) trim(qck(ng)%head), my_outer
        lstr=len_trim(qck(ng)%name)
        qck(ng)%base=qck(ng)%name(1:lstr-3)
        IF (qck(ng)%Nfiles.gt.1) THEN
          DO ifile=1,qck(ng)%Nfiles
            WRITE (suffix,"('_',i4.4,'.nc')") ifile
            qck(ng)%files(ifile)=trim(qck(ng)%base)//trim(suffix)
          END DO
          qck(ng)%name=trim(qck(ng)%files(1))
        ELSE
          qck(ng)%files(1)=trim(qck(ng)%name)
        END IF
      END DO
!
!  If weak constraint, the impulses are time-interpolated at each
!  time-steps.
!
      DO ng=1,ngrids
# ifdef SPLIT_4DVAR
        IF (nadj(ng).lt.ntimes(ng)) THEN
          SELECT CASE (adm(ng)%IOtype)
            CASE (io_nf90)
              CALL netcdf_get_dim (ng, inlm, adm(ng)%name,              &
     &                             dimname = 'ocean_time',              &
     &                             dimsize = frcrec(ng))
#  if defined PIO_LIB && defined DISTRIBUTE
            CASE (io_pio)
              CALL pio_netcdf_get_dim (ng, inlm, adm(ng)%name,          &
     &                                 dimname = 'ocean_time',          &
     &                                 dimsize = frcrec(ng))
#  endif
          END SELECT
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END IF
# endif
        IF (frcrec(ng).gt.3) THEN
          frequentimpulse(ng)=.true.
        END IF
      END DO
!
!  Initialize nonlinear model INI(ng)%name file, record outer+2.
!  Notice that NetCDF record index counter is saved because this
!  counter is used to write initial conditions.
!
      DO ng=1,ngrids
        indxsave(ng)=ini(ng)%Rindex
        ini(ng)%Rindex=my_outer+2
        rst(ng)%Rindex=0
        fcount=rst(ng)%load
        rst(ng)%Nrec(fcount)=0
      END DO
 
# ifdef RPCG
!
! Compute the incremental analysis update forcing for the NLM.
!
!    Read the accumulated increment from Rec3 of the ITL file.
!
        DO ng=1,ngrids
          IF (timeiau(ng).gt.0.0_dp) THEN
            iauswitch(ng)=.true.
            CALL get_state (ng, iadm, 4, itl(ng), rec3, ltlm1)
          END IF
        END DO
!
!    Set up the NL model forcing contribution for the incremental
!    analysis update.
!
        DO ng=1,ngrids
          IF (timeiau(ng).gt.0.0_dp) THEN
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL frc_iau_ini (ng, tile, ltlm1)
            END DO
          END IF
        END DO
# endif
!
      CALL initial
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      DO ng=1,ngrids
        ini(ng)%Rindex=indxsave(ng)
      END DO
!
!  Activate switch to write out final misfit between model and
!  observations.
!
      IF (my_outer.eq.nouter) THEN
        DO ng=1,ngrids
          wrtmisfit(ng)=.true.
        END DO
      END IF
!
 10   FORMAT (a,'_outer',i0,'.nc')
!
      RETURN
      END SUBROUTINE analysis_initialize
!
      SUBROUTINE analysis (my_outer, RunInterval)
!
!=======================================================================
!                                                                      !
!  This routine computes 4D-Var data assimilation analysis, Xa. The    !
!  nonlinear model initial conditions are computed by adding the       !
!  4D-Var increments to the current background:  Xa = Xb + dXa.        !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     my_outer        Outer-loop counter (integer)                     !
!     RunInterval     NLM kernel time stepping window (seconds)        !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      integer, intent(in) :: my_outer
!
      real(dp), intent(in) :: runinterval
!
!  Local variable declarations.
!
      logical :: donestepping
!
      integer :: i, lstr, ng
      integer :: fcount
# ifdef PROFILE
      integer :: thread
# endif
# if defined MODEL_COUPLING && !defined MCT_LIB
      integer :: nstrstep, nendstep, extra
!
      real(dp) :: endtime, nexttime
# endif
!
      character (len=20) :: string
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", analysis"
!
      sourcefile=myfile
!
# if !(defined MODEL_COUPLING && defined ESMF_LIB)
!
!-----------------------------------------------------------------------
!  Set nonlinear model initial conditions.
!-----------------------------------------------------------------------
!
      CALL analysis_initialize (my_outer)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
!
!-----------------------------------------------------------------------
!  Run nonlinear model and compute a "new estimate" of the state
!  trajectory, Xb(t)|n.
# if defined MODEL_COUPLING && !defined MCT_LIB
!  Since the ROMS kernel has a delayed output and line diagnostics by
!  one timestep, subtact an extra value to the report of starting and
!  ending timestep for clarity. Usually, the model coupling interval
!  is of the same size as ROMS timestep.
# endif
!-----------------------------------------------------------------------
 
# ifdef PROFILE
!
!  Start profile clock.
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_on (ng, inlm, 88, __line__, myfile)
        END DO
      END DO
# endif
!
!  Initialize various parameters and switches.
!
      myruninterval=runinterval
!
      DO ng=1,ngrids
# ifdef AVERAGES
        idefavg(ng)=-1
        ldefavg(ng)=.true.
        lwrtavg(ng)=.true.
        WRITE (avg(ng)%name,10) trim(avg(ng)%head), my_outer
        lstr=len_trim(avg(ng)%name)
        avg(ng)%base=avg(ng)%name(1:lstr-3)
# endif
# ifdef DIAGNOSTICS
        idefdia(ng)=-1
        ldefdia(ng)=.true.
        lwrtdia(ng)=.true.
        WRITE (dia(ng)%name,10) trim(dia(ng)%head), my_outer
        lstr=len_trim(dia(ng)%name)
        dia(ng)%base=dia(ng)%name(1:lstr-3)
# endif
# if defined MODEL_COUPLING && !defined MCT_LIB
!
        nexttime=time(ng)+runinterval
        endtime=initime(ng)+(ntimes(ng)-1)*dt(ng)
        IF ((nexttime.eq.endtime).and.(ng.eq.1)) THEN
          extra=0                                   ! last time interval
        ELSE
          extra=1
        END IF
        step_counter(ng)=0
        nstrstep=iic(ng)
        nendstep=nstrstep+int((myruninterval)/dt(ng))-extra
        IF (master) WRITE (stdout,20) 'NL', ng, my_outer, ninner,       &
     &                                nstrstep, nendstep
# else
        IF (master) WRITE (stdout,20) 'NL', ng, my_outer, ninner,       &
     &                                ntstart(ng), ntend(ng)
# endif
      END DO
!
!  Run nonlinear model kernel.
!
# ifdef SOLVE3D
      CALL main3d (runinterval)
# else
      CALL main2d (runinterval)
# endif
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!-----------------------------------------------------------------------
!  If completed stepping, reset switches and write out cost function.
# if defined MODEL_COUPLING && !defined MCT_LIB
!  In coupled applications, RunInterval is much less than ntimes*dt,
!  so we need to wait until the last coupling interval is finished.
!  Otherwise, the control switches will be turned off prematurely.
# endif
!-----------------------------------------------------------------------
!
# if defined MODEL_COUPLING && !defined MCT_LIB
      IF (nendstep.eq.ntend(1)) THEN
        donestepping=.true.
      ELSE
        donestepping=.false.
      END IF
# else
      donestepping=.true.
# endif
!
      IF (donestepping) THEN
        DO ng=1,ngrids
# ifdef AVERAGES
          ldefavg(ng)=.false.
          lwrtavg(ng)=.false.
# endif
# ifdef DIAGNOSTICS
          ldefdia(ng)=.false.
          lwrtdia(ng)=.false.
# endif
          wrtnlmod(ng)=.false.
          wrttlmod(ng)=.false.
        END DO
 
# ifdef RPCG
!
!  Deactivate the NL forcing term for iau.
!
        DO ng=1,ngrids
          iauswitch(ng)=.false.
        END DO
# endif
!
!  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.
!
        CALL edit_multifile ('HIS2FWD')
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        DO ng=1,ngrids
          lreadfwd(ng)=.true.
        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.
!
        CALL edit_multifile ('QCK2BLK')
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        DO ng=1,ngrids
          lreadblk(ng)=.true.
        END DO
# endif
!
!  Report data penalty function.
!
        DO ng=1,ngrids
          IF (master) THEN
            DO i=0,nobsvar(ng)
              IF (i.eq.0) THEN
                string='Total'
              ELSE
                string=obsname(i)
              END IF
              IF (fourdvar(ng)%NLPenalty(i).ne.0.0_r8) THEN
                WRITE (stdout,30) my_outer, ninner, 'NLM',              &
     &                            fourdvar(ng)%NLPenalty(i),            &
     &                            trim(string)
              END IF
            END DO
          END IF
!
!  Write out final data penalty function to NetCDF file.
!
          sourcefile=myfile
          SELECT CASE (dav(ng)%IOtype)
            CASE (io_nf90)
              CALL netcdf_put_fvar (ng, inlm, dav(ng)%name,             &
     &                              'NL_fDataPenalty',                  &
     &                              fourdvar(ng)%NLPenalty(0:),         &
     &                              (/1,my_outer/),                     &
     &                              (/nobsvar(ng)+1,1/),                &
     &                              ncid = dav(ng)%ncid)
 
# if defined PIO_LIB && defined DISTRIBUTE
            CASE (io_pio)
              CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name,         &
     &                                  'NL_fDataPenalty',              &
     &                                  fourdvar(ng)%NLPenalty(0:),     &
     &                                  (/1,my_outer/),                 &
     &                                  (/nobsvar(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 NL model.
!
          fourdvar(ng)%NLPenalty=0.0_r8
        END DO
      END IF
 
# ifdef PROFILE
!
!  Stop profile clock
!
      DO ng=1,ngrids
        DO thread=thread_range
          CALL wclock_off (ng, inlm, 88, __line__, myfile)
        END DO
      END DO
# endif
!
 10   FORMAT (a,'_outer',i0,'.nc')
 20   FORMAT (/,1x,a,1x,'ROMS: started time-stepping:',                 &
     &        ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3,          &
              ', TimeSteps: ',i0,' - ',i0,')',/)
 30   FORMAT (' (',i3.3,',',i3.3,'): ',a,' data penalty, Jdata = ',     &
     &        1p,e17.10,0p,t68,a)
!
      RETURN
      END SUBROUTINE analysis
!
      SUBROUTINE prior_error (ng)
!
!=======================================================================
!                                                                      !
!  This routine processes background prior error covariance standard   !
!  deviations and normalization coefficients.                          !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng              Nested grid number                               !
!                                                                      !
!=======================================================================
!
      USE def_norm_mod,      ONLY : def_norm
      USE normalization_mod, ONLY : normalization
      USE strings_mod,       ONLY : founderror
!
!  Imported variable declarations
!
      integer, intent(in) :: ng
!
!  Local variable declarations.
!
      integer :: tile
      integer :: nrmrec, stdrec, tindex
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", prior_error"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Set application grid, metrics, and associated variables and
!  parameters.
!-----------------------------------------------------------------------
!
!  The ROMS application grid configuration is done once. It is usually
!  done in the "initial" kernel routine. However, since we are calling
!  the "normalization" routine here, we need several grid variables and
!  parameter.  Also, if reading only water points, we need to know the
!  land/sea mask arrays to unpack.
!
      IF (setgridconfig(ng)) THEN
        CALL set_grid (ng, inlm)
      END IF
!
!-----------------------------------------------------------------------
!  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.
!
      IF (lreadstd(ng)) THEN
        stdrec=1
        tindex=1
# ifdef STD_MODEL
        CALL get_state (ng, 10, 10, std(5,ng), stdrec, tindex)
# else
        CALL get_state (ng, 10, 10, std(1,ng), stdrec, tindex)
# endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  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
        CALL get_state (ng, 11, 11, std(2,ng), stdrec, tindex)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
 
# ifdef ADJUST_BOUNDARY
!
!  Open boundary conditions standard deviation.
!
      stdrec=1
      tindex=1
      CALL get_state (ng, 12, 12, std(3,ng), stdrec, tindex)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
 
# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
!
!  Surface forcing standard deviation.
!
      stdrec=1
      tindex=1
      CALL get_state (ng, 13, 13, std(4,ng), stdrec, tindex)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
# endif
 
# ifdef STD_MODEL
!
!  If computing the standard deviation (STD) directly from the
!  background (prior) field as an alternative to climatological
!  read from the input NetCDF file, define the output STD NetCDF
!  file.
!
      IF (ldefstd(ng)) THEN
        CALL def_std (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
# endif
!
!-----------------------------------------------------------------------
!  Error covariance normalization coefficients.
!-----------------------------------------------------------------------
!
!  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.
!
      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
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        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
!
      RETURN
      END SUBROUTINE prior_error
 
# if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F || \
     defined posterior_eofs
!
      SUBROUTINE posterior_error (RunInterval)
!
!=======================================================================
!                                                                      !
!  This routine computes posterior analysis error covariance matrix    !
!  using the Lanczos vectors.                                          !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     RunInterval     Kernel time stepping window (seconds)            !
!                                                                      !
!  Warning:                                                            !
!                                                                      !
!  Currently, this code only works for a single outer-loop.            !
!                                                                      !
!=======================================================================
!
      USE random_ic_mod,     ONLY : random_ic
!
!  Imported variable declarations
!
      real(dp), intent(in) :: runinterval
!
!  Local variable declarations.
!
      logical :: ltrace
!
      integer :: my_inner, my_outer, ng, tile
      integer :: fcount, rec
!
      character (len=8) :: driver = 'rbl4dvar'
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", posterior_error"
!
      sourcefile=myfile
 
#  if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F
!
!-----------------------------------------------------------------------
!  Compute full (diagonal) posterior analysis error covariance matrix.
!-----------------------------------------------------------------------
!
!  Clear tangent and adjoint arrays because they are used as
!  work arrays below.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL initialize_ocean (ng, tile, iadm)
          CALL initialize_ocean (ng, tile, itlm)
          CALL initialize_forces (ng, tile, iadm)
          CALL initialize_forces (ng, tile, itlm)
#   ifdef ADJUST_BOUNDARY
          CALL initialize_boundary (ng, tile, iadm)
          CALL initialize_boundary (ng, tile, itlm)
#   endif
        END DO
      END DO
!
!  Compute the diagonal of the posterior/analysis error covariance
!  matrix. The result is written to record 2 of the ITL netcdf file.
!
      var_oloop : DO my_outer=1,nouter
        outer=my_outer
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL posterior_var (ng, tile, itlm, outer)
          END DO
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END DO var_oloop
!
!  Write out the diagonal of the posterior/analysis covariance matrix
!  which is in tl_var(Rec1) to 4DVar error NetCDF file.
!
      DO ng=1,ngrids
#   ifdef DISTRIBUTE
        CALL wrt_error (ng, myrank, rec1, rec1)
#   else
        CALL wrt_error (ng, -1, rec1, rec1)
#   endif
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!  Clear tangent and adjoint arrays because they are used as
!  work arrays below.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL initialize_ocean (ng, tile, iadm)
          CALL initialize_ocean (ng, tile, itlm)
          CALL initialize_forces (ng, tile, iadm)
          CALL initialize_forces (ng, tile, itlm)
#   ifdef ADJUST_BOUNDARY
          CALL initialize_boundary (ng, tile, iadm)
          CALL initialize_boundary (ng, tile, itlm)
#   endif
        END DO
      END DO
#  endif
 
#  ifdef POSTERIOR_EOFS
!
!-----------------------------------------------------------------------
!  Compute the posterior analysis error covariance matrix EOFs using a
!  Lanczos algorithm.
!
!  NOTE: Currently, this code only works for a single outer-loop.
!-----------------------------------------------------------------------
!
      IF (master) WRITE (stdout,10)
!
!  Estimate first the trace of the posterior analysis error
!  covariance matrix since the evolved and convolved Lanczos
!  vectors stored in the Hessian NetCDF file will be destroyed
!  later.
!
      ltrace=.true.
 
      trace_oloop : DO my_outer=1,nouter
        outer=my_outer
        inner=0
 
        trace_iloop : DO my_inner=1,nposti
          inner=my_inner
!
!  Initialize the tangent linear variables with a random vector
!  comprised of +1 and -1 elements randomly chosen.
!
          DO ng=1,ngrids
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL random_ic (ng, tile, itlm, inner, outer,             &
     &                        lold(ng), ltrace)
            END DO
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
!
!  Apply horizontal convolution.
!
          CALL convolve (driver, lini, lold, lnew)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Compute Lanczos vector and eigenvectors of the posterior analysis
!  error covariance matrix.
!
          DO ng=1,ngrids
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL posterior (ng, tile, itlm, inner, outer, ltrace)
            END DO
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
 
        END DO trace_iloop
 
      END DO trace_oloop
!
!  Estimate posterior analysis error covariance matrix.
!
      ltrace=.false.
 
      post_oloop : DO my_outer=1,nouter
        outer=my_outer
        inner=0
!
!  The Lanczos algorithm requires to save all the Lanczos vectors.
!  They are used to compute the posterior EOFs.
!
        DO ng=1,ngrids
          adm(ng)%Rindex=0
          fcount=adm(ng)%load
          adm(ng)%Nrec(fcount)=0
        END DO
 
        post_iloop : DO my_inner=0,nposti
          inner=my_inner
!
!  Read first record of ITL file and apply convolutions.
!
!  NOTE: If inner=0, we would like to use a random starting vector.
!        For now we can use what ever is in record 1.
!
          IF (inner.ne.0) THEN
            DO ng=1,ngrids
              rec=1
              CALL get_state (ng, itlm, 1, itl(ng), rec, lold(ng))
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
          ELSE
 
            DO ng=1,ngrids
              DO tile=first_tile(ng),last_tile(ng),+1
                CALL random_ic (ng, tile, itlm, inner, outer,           &
     &                          lold(ng), ltrace)
              END DO
              IF (founderror(exit_flag, noerror,                        &
     &                       __line__, myfile)) RETURN
            END DO
          END IF
!
!  Apply horizontal convolution.
!
          CALL convolve (driver, lini, lold, lnew)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Compute Lanczos vector and eigenvectors of the posterior analysis
!  error covariance matrix.
!
          DO ng=1,ngrids
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL posterior (ng, tile, itlm, inner, outer, ltrace)
            END DO
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
!
!   Write the Lanczos vectors of the posterior error covariance
!   to the adjoint NetCDF file.
!
          DO ng=1,ngrids
#   if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
            lfout(ng)=lnew(ng)
#   endif
#   ifdef ADJUST_BOUNDARY
            lbout(ng)=lnew(ng)
#   endif
            kstp(ng)=lnew(ng)
#   ifdef SOLVE3D
            nstp(ng)=lnew(ng)
#   endif
            lwrtstate2d(ng)=.true.
#   ifdef DISTRIBUTE
            CALL ad_wrt_his (ng, myrank)
#   else
            CALL ad_wrt_his (ng, -1)
#   endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
            lwrtstate2d(ng)=.false.
          END DO
!
!  Write out tangent linear model initial conditions and tangent
!  linear surface forcing adjustments for next inner
!  loop into ITL(ng)%name (record Rec1).
!
          DO ng=1,ngrids
#   ifdef DISTRIBUTE
            CALL tl_wrt_ini (ng, myrank, lnew(ng), rec1)
#   else
            CALL tl_wrt_ini (ng, -1, lnew(ng), rec1)
#   endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
 
        END DO post_iloop
 
      END DO post_oloop
 
#  endif
!
!  Done.  Set history file ID to closed state since we manipulated
!  its indices with the forward file ID which was closed above.
!
      DO ng=1,ngrids
        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)
          his(ng)%pioFile%fh=-1
#  endif
        END IF
      END DO
 
#  ifdef POSTERIOR_EOFS
!
 10   FORMAT (/,' <<<< Posterior Analysis Error Covariance Matrix',     &
     &          ' Estimation >>>>',/)
#  endif
!
      RETURN
      END SUBROUTINE posterior_error
# endif
 
#endif
      END MODULE rbl4dvar_mod