convolve.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE convolve_mod
 
#ifdef FOUR_DVAR
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group            Andrew M. Moore   !
!    Licensed under a MIT/X style license                              !
!=======================================================================
!                                                                      !
!  These routines impose the error covariance matrices for the 4D-Var  !
!  control variables. The  error covariances matrices  are modeled by  !
!  the spatial convolutions of pseudo-diffusion equations.             !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     driver       Calling driver (string)                             !
!     model        Initial conditions model to process (iTLM or iRPM)  !
!     outLoop      Current 4D-Var outer loop                           !
!     innLoop      Current 4D-Var inner loop                           !
!     Rbck         NLM background record to process                    !
!     Rini         NLM initial condition record to process             !
!     Rold         State vector old minimization time index            !
!     Rnew         State vector new minimization time index            !
!     Rec1         TLM initial conditions record to write              !
!     Rec2         RPM initial conditions record to write              !
!     Lposterior   Switch to process posterior error covariance        !
!                                                                      !
!  Reference:                                                          !
!                                                                      !
!    Weaver, A. and P. Courtier, 2001: Correlation modeling on the     !
!      sphere using a generalized diffusion equation, Q.J.R. Meteo.    !
!      Soc, 127, 1815-1846.                                            !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
# ifdef RPCG
      USE mod_fourdvar
# endif
      USE mod_iounits
      USE mod_scalars
      USE mod_stepping
!
# ifdef ADJUST_BOUNDARY
      USE mod_boundary,       ONLY : initialize_boundary
# endif
      USE mod_forces,         ONLY : initialize_forces
      USE mod_ocean,          ONLY : initialize_ocean
!
# ifdef BALANCE_OPERATOR
      USE ad_balance_mod,     ONLY : ad_balance
# endif
      USE ad_convolution_mod, ONLY : ad_convolution
      USE ad_variability_mod, ONLY : ad_variability
      USE ad_wrt_his_mod,     ONLY : ad_wrt_his
# ifdef RPCG
      USE comp_jb0_mod,       ONLY : aug_oper
# endif
      USE get_state_mod,      ONLY : get_state
      USE ini_adjust_mod,     ONLY : load_adtotl
      USE ini_adjust_mod,     ONLY : load_tltoad
# if defined ARRAY_MODES            || defined R4DVAR || \
     defined r4dvar_ana_sensitivity
      USE ini_adjust_mod,     ONLY : rp_ini_adjust
# elif defined RBL4DVAR  || defined RBL4DVAR_ANA_SENSITIVITY
      USE ini_adjust_mod,     ONLY : ini_adjust
# endif
# if defined ARRAY_MODES            || defined R4DVAR || \
     defined r4dvar_ana_sensitivity
      USE rp_wrt_ini_mod,     ONLY : rp_wrt_ini
# endif
      USE strings_mod,        ONLY : founderror
# ifdef RPCG
      USE sum_grad_mod,       ONLY : sum_grad
# endif
# ifdef TIME_CONV
      USE time_corr_mod,      ONLY : time_corr
# endif
# ifdef BALANCE_OPERATOR
      USE tl_balance_mod,     ONLY : tl_balance
# endif
      USE tl_convolution_mod, ONLY : tl_convolution
      USE tl_variability_mod, ONLY : tl_variability
      USE tl_wrt_ini_mod,     ONLY : tl_wrt_ini
# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
     defined adjust_boundary
      USE wrt_ini_mod,        ONLY : wrt_frc_ad
# endif
      USE wrt_ini_mod,        ONLY : wrt_ini
# ifdef POSTERIOR_ERROR_I
      USE wrt_hessian_mod,    ONLY : wrt_hessian
# endif
# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
     defined adjust_boundary
      USE wrt_ini_mod,        ONLY : wrt_frc_ad
# endif
      USE wrt_ini_mod,        ONLY : wrt_ini
 
!
      implicit none
!
      PRIVATE
      PUBLIC  :: convolve
      PUBLIC  :: error_covariance
# ifdef SP4DVAR
      PUBLIC  :: saddlec
# endif
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE convolve (driver, Rini, Rold, Rnew)
!***********************************************************************
!
!  Imported variable declarations.
!
      integer, intent(in) :: rini
      integer, intent(in) :: rold(ngrids)
      integer, intent(in) :: rnew(ngrids)
!
      character (len=*), intent(in) :: driver
!
!  Local variable declarations.
!
      logical :: lweak, add
!
      integer :: inirec
      integer :: ng, tile
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Specify spatial error covariance on the TLM control vector, at time
!  index "Rold", via convolutions of a pseudo-diffusion operator.  The
!  convolved control vector is loaded into time index "Rnew" and index
!  "Rold" is preserved.
!-----------------------------------------------------------------------
 
# ifdef BALANCE_OPERATOR
!
!  Read NLM initial condition in readiness for the balance operator.
!
      inirec=rini
      DO ng=1,ngrids
        CALL get_state (ng, inlm, 2, ini(ng), inirec, inirec)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        nrhs(ng)=rini
      END DO
# endif
!
!  Load "Rold" time index TLM control vector into ADM state arrays.
!  Apply balance operator, if activated. Then, scale control vector
!  with error covariance standard deviations. Next, convolve resulting
!  vector with the squared-root adjoint diffusion operator. Notice
!  that the spatial convolution is only done for half of the diffusion
!  steps (squared-root filter).
!
      lweak=.false.
      add=.false.
      DO ng=1,ngrids
# ifdef PROFILE
        CALL wclock_on (ng, iadm, 82, __line__, myfile)
# endif
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
# ifdef BALANCE_OPERATOR
          CALL ad_balance (ng, tile, rini, rold(ng))
# endif
          CALL ad_variability (ng, tile, rold(ng), lweak)
          CALL ad_convolution (ng, tile, rold(ng), lweak, 2)
        END DO
# ifdef PROFILE
        CALL wclock_off (ng, iadm, 82, __line__, myfile)
# endif
      END DO
!
!  Since we wish to preserve what is in tl_var(Rold), load resulting
!  filtered solution from above, ad_var(Rold), into tl_var(Rnew).
!  Then, convolve with the squared-root (half of steps) tangent
!  linear diffusion operator. Next, scale results with error
!  covariance standard deviations. Apply balance operator, if
!  activated.
!
      add=.false.
      DO ng=1,ngrids
# ifdef PROFILE
        CALL wclock_on (ng, itlm, 82, __line__, myfile)
# endif
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL load_adtotl (ng, tile, rold(ng), rnew(ng), add)
          CALL tl_convolution (ng, tile, rnew(ng), lweak, 2)
          CALL tl_variability (ng, tile, rnew(ng), lweak)
# ifdef BALANCE_OPERATOR
          CALL tl_balance (ng, tile, rini, rnew(ng))
# endif
        END DO
# ifdef PROFILE
        CALL wclock_off (ng, itlm, 82, __line__, myfile)
# endif
      END DO
!
      RETURN
      END SUBROUTINE convolve
!
!***********************************************************************
      SUBROUTINE error_covariance (model, driver, outLoop, innLoop,     &
     &                             Rbck, Rini, Rold, Rnew,              &
     &                             Rec1, Rec2, Lposterior)
!***********************************************************************
!
!  Imported variable declarations.
!
      logical, intent(in) :: lposterior
!
      integer, intent(in) :: model, outloop, innloop
      integer, intent(in) :: rbck, rini, rec1, rec2
      integer, intent(in) :: rold(ngrids)
      integer, intent(in) :: rnew(ngrids)
!
      character (len=*), intent(in) :: driver
!
!  Local variable declarations.
!
      logical :: lweak, add
!
      integer :: adrec, bckrec, fcount, inirec
      integer :: i, irec, ng, tile
# ifdef RPCG
      integer :: ltlm1, ltlm2, rec5, jrec, nadrec
# endif
      integer, dimension(Ngrids) :: nrec
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", error_covariance"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Convolve adjoint trajectory.
!-----------------------------------------------------------------------
!
      IF (master) THEN
        IF ((outloop.lt.0).and.(innloop.lt.0)) THEN
          WRITE (stdout,10)
 10       FORMAT (/,' Convolving Adjoint Trajectory...')
        ELSE
          WRITE (stdout,20) outloop, innloop
 20       FORMAT (/,' Convolving Adjoint Trajectory: Outer = ',i3.3,    &
     &            ' Inner = ',i3.3)
        END IF
      END IF
!
!  Clear adjoint state arrays.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL initialize_ocean (ng, tile, iadm)
        END DO
      END DO
!
!  Get number of records in adjoint NetCDF and initialize indices.
!
      DO ng=1,ngrids
        fcount=adm(ng)%Fcount
        nrec(ng)=adm(ng)%Nrec(fcount)
        adm(ng)%Nrec(fcount)=0
        adm(ng)%Rindex=0
        lwrtstate2d(ng)=.true.
        lwrttime(ng)=.false.
      END DO
!
!  Read in adjoint control vector to convolve from NetCDF file
!  ADM(ng)%name for time record Nrec, which correspond to the
!  initial conditions time. If adjusting open boundaries and/or
!  surface forcing, only record Nrec is read since it is the
!  only record for which adjoint forcing arrays are complete.
!
!  Notice that since routine "get_state" loads data into the
!  ghost points, the adjoint contol vector is read into the
!  tangent linear state arrays by using iTLM instead of iADM
!  in the calling arguments.
!
      DO ng=1,ngrids
        adrec=nrec(ng)
        frcrec(ng)=nrec(ng)
        CALL get_state (ng, itlm, 4, adm(ng), adrec, rold(ng))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
 
# ifdef BALANCE_OPERATOR
!
!  Read NLM initial condition in readiness for the balance
!  operator.
!
      inirec=rini
      DO ng=1,ngrids
        CALL get_state (ng, inlm, 2, ini(ng), inirec, inirec)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        nrhs(ng)=rini
      END DO
# endif
!
!  Load ADM control read above from TLM into ADM state arrays.
!  Then, multiply control vector by the error covariance standard
!  deviations. Next, convolve resulting adjoint solution with the
!  squared-root adjoint diffusion operator to impose apriori error
!  hypothesis. Notice that the spatial convolution are only done
!  for half of the diffusion steps (squared-root filter). Clear
!  tangent linear state arrays when done.
!
      lweak=.false.
      add=.false.
      DO ng=1,ngrids
# ifdef PROFILE
        CALL wclock_on (ng, model, 82, __line__, myfile)
# endif
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
# ifdef BALANCE_OPERATOR
          CALL ad_balance (ng, tile, rini, rold(ng))
# endif
          CALL ad_variability (ng, tile, rold(ng), lweak)
          CALL ad_convolution (ng, tile, rold(ng), lweak, 2)
          CALL initialize_ocean (ng, tile, itlm)
          CALL initialize_forces (ng, tile, itlm)
# ifdef ADJUST_BOUNDARY
          CALL initialize_boundary (ng, tile, itlm)
# endif
        END DO
# ifdef PROFILE
        CALL wclock_off (ng, model, 82, __line__, myfile)
# endif
      END DO
!
!  To insure symmetry, convolve resulting filtered adjoint solution
!  from above with the squared-root (half of steps) tangent linear
!  diffusion operator. Then, multiply result with its corresponding
!  error covariance standard deviations. Since the convolved solution
!  is in the adjoint state arrays, first copy to tangent linear state
!  arrays including the ghosts points.
!
      lweak=.false.
      add=.false.
      DO ng=1,ngrids
# ifdef PROFILE
        CALL wclock_on (ng, model, 82, __line__, myfile)
# endif
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL load_adtotl (ng, tile, rold(ng), rold(ng), add)
          CALL tl_convolution (ng, tile, rold(ng), lweak, 2)
          CALL tl_variability (ng, tile, rold(ng), lweak)
# ifdef BALANCE_OPERATOR
          CALL tl_balance (ng, tile, rini, rold(ng))
# endif
        END DO
 
# ifdef POSTERIOR_ERROR_I
!
!  If computing the analysis error covariance matrix, copy TLM back
!  into ADM so that it can be written to Hessian NetCDF file.
!
        IF (lposterior) THEN
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
          END DO
        END IF
# endif
# ifdef PROFILE
        CALL wclock_off (ng, model, 82, __line__, myfile)
# endif
      END DO
 
# if defined ARRAY_MODES            || defined R4DVAR || \
     defined r4dvar_ana_sensitivity
!
!  Copy back to adjoint state arrays when done with the convolution.
!  Compute representer model initial conditions by adding convolved
!  adjoint solution to the reference nonlinear state (INI(ng)%name,
!  record Rbck).
!
      IF ((model.eq.irpm).and.                                          &
     &    (index(driver,'r4dvar').ne.0)) THEN
        bckrec=rbck
        DO ng=1,ngrids
          CALL get_state (ng, inlm, 9, ini(ng), bckrec, rnew(ng))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
        add=.false.
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
            CALL rp_ini_adjust (ng, tile, rnew(ng), rold(ng))
          END DO
        END DO
      END IF
 
# elif defined RBL4DVAR || defined RBL4DVAR_ANA_SENSITIVITY
#  ifndef RPCG
!
!  Copy back to adjoint state arrays when done with the convolution.
!  Compute nonlinear model initial conditions by adding convolved
!  adjoint solution to the reference nonlinear state (INI(ng)%name,
!  record Rbck).
!
      IF ((model.eq.inlm).and.                                          &
     &    (index(driver,'rbl4dvar').ne.0)) THEN
        bckrec=rbck
        DO ng=1,ngrids
          CALL get_state (ng, inlm, 9, ini(ng), bckrec, rnew(ng))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
        add=.false.
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
            CALL ini_adjust (ng, tile, rold(ng), rnew(ng))
          END DO
        END DO
      END IF
#  endif
# endif
!
!  Write out tangent linear model initial conditions and tangent
!  linear surface forcing adjustments for next inner loop into
!  ITL(ng)%name (record Rec1). The tangent model initial
!  conditions are set to the convolved adjoint solution.
!
      IF (model.eq.itlm) THEN
        DO ng=1,ngrids
# ifdef DISTRIBUTE
          CALL tl_wrt_ini (ng, myrank, rold(ng), rec1)
# else
          CALL tl_wrt_ini (ng, -1, rold(ng), rec1)
# endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
# if defined ARRAY_MODES            || defined R4DVAR || \
     defined r4dvar_ana_sensitivity
      ELSE IF (model.eq.irpm) THEN
        DO ng=1,ngrids
#  ifdef DISTRIBUTE
          CALL rp_wrt_ini (ng, myrank, rold(ng), rec2)
#  else
          CALL rp_wrt_ini (ng, -1, rold(ng), rec2)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
# elif (defined RBL4DVAR || defined RBL4DVAR_ANA_SENSITIVITY) && \
       !defined RPCG
      ELSE IF (model.eq.inlm) THEN
        DO ng=1,ngrids
#  ifdef DISTRIBUTE
          CALL wrt_ini (ng, myrank, rnew(ng))
#  else
          CALL wrt_ini (ng, -1, rnew(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, rold(ng), ini(ng)%Rindex)
#   else
          CALL wrt_frc_ad (ng, -1, rold(ng), ini(ng)%Rindex)
#   endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
        END DO
# endif
      END IF
 
# ifdef RPCG
!
!  Write out tangent linear model initial conditions and tangent
!  linear surface forcing and obc adjustments for next outer
!  loop into ITLname (record Rec2). The tangent model initial
!  conditions are set to the convolved adjoint solution.
!
      IF (model.eq.inlm) THEN
        DO ng=1,ngrids
#  ifdef DISTRIBUTE
          CALL tl_wrt_ini (ng, myrank, rold(ng), rec2)
#  else
          CALL tl_wrt_ini (ng, -1, rold(ng), rec2)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END IF
# endif
 
# ifdef POSTERIOR_ERROR_I
!
!  If weak constraint, convolve records 2:Nrec in ADM(ng)%name and
!  Write convolved adjoint solution into Hessian NetCDF file for use
!  later.
!
      IF (lposterior.and.(inner.ne.0)) THEN
        DO ng=1,ngrids
#  ifdef DISTRIBUTE
          CALL wrt_hessian (ng, myrank, rold(ng), rold(ng))
#  else
          CALL wrt_hessian (ng, -1, rold(ng), rold(ng))
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END IF
# endif
# ifdef RPCG
!
!  Now compute the new B^-1(x(k)-x(k-1)) term for the weak constraint
!  forcing terms. Records 6+nADrec/2 to nADrec+5 contain the sums so
!  far. This is ONLY done in the outer-loop not the inner-loops, and
!  is controlled by the flag "LaugWeak".
!
      IF (laugweak) THEN
        ltlm1=1
        ltlm2=2
        rec5=5
        DO ng=1,ngrids
          IF (nadj(ng).lt.ntimes(ng)) THEN
            nadrec=2*(1+ntimes(ng)/nadj(ng))
          ELSE
            nadrec=0
          END IF
          DO irec=1,nadrec/2
            jrec=rec5+nadrec/2+irec
!
!  First add the augmented piece which is computed from the previous
!  sum.
!
            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
!
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL sum_grad (ng, tile, ltlm1, ltlm2)
            END DO
!
!  Need to read adjoint netcdf file in reverse.
!
            adrec=(nrec(ng)-1)-(irec-1)
            CALL get_state (ng, itlm, 4, adm(ng), adrec, 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)
            END DO
!
!  Write the current sum of adjoint solutions into record jrec of the
!  ITL file.
!
#  ifdef DISTRIBUTE
            CALL tl_wrt_ini (ng, myrank, ltlm2, jrec)
#  else
            CALL tl_wrt_ini (ng, -1, ltlm2, jrec)
#  endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          END DO
        END DO
      END IF
# endif
# ifdef TIME_CONV
!
!  Apply the factorized space-time model error covariance
!  matrix of the form ATA' where A is the square-root of the
!  model error covariance matrix, and T is the time-correlation
!  matrix.
!
!  If weak constraint, convolve records 2-Nrec in ADM(ng)%name
!  and impose model error covariance. NOTE: We will not use the
!  convolved forcing increments generated here since these arrays
!  do not contain the complete solution and are redundant.
!  AMM: We might want to get rid of these unwanted records to
!  avoid any confusion in the future.
!
      DO ng=1,ngrids
        IF (nrec(ng).gt.3) THEN
#  ifdef PROFILE
          CALL wclock_on (ng, model, 82, __line__, myfile)
#  endif
          DO irec=1,nrec(ng)-1
!
!  Read adjoint solution. Since routine "get_state" loads data into the
!  ghost points, the adjoint solution is read in the tangent linear
!  state arrays by using iTLM instead of iADM in the calling arguments.
!
            adrec=irec
            CALL get_state (ng, itlm, 4, adm(ng), adrec, rold(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Load interior solution, read above, into adjoint state arrays.
!  Then, multiply adjoint solution by the background-error standard
!  deviations. Next, convolve resulting adjoint solution with the
!  squared-root adjoint diffusion operator which impose the model-error
!  spatial correlations. Notice that the spatial convolution is only
!  done for half of the diffusion steps (squared-root filter). Clear
!  tangent linear state arrays when done.
!
            lweak=.true.
            add=.false.
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
#  ifdef BALANCE_OPERATOR
              CALL ad_balance (ng, tile, rini, rold(ng))
#  endif
              CALL ad_variability (ng, tile, rold(ng), lweak)
              CALL ad_convolution (ng, tile, rold(ng), lweak, 2)
              CALL initialize_ocean (ng, tile, itlm)
              CALL initialize_forces (ng, tile, itlm)
#  ifdef ADJUST_BOUNDARY
              CALL initialize_boundary (ng, tile, itlm)
#  endif
            END DO
!
!  Overwrite ADM(ng)%name history NetCDF file with convolved adjoint
!  solution.
!
            kstp(ng)=rold(ng)
#  ifdef SOLVE3D
            nstp(ng)=rold(ng)
#  endif
#  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
          lwrtstate2d(ng)=.false.
          lwrttime(ng)=.true.
#  ifdef PROFILE
        CALL wclock_off (ng, model, 82, __line__, myfile)
#  endif
        END IF
      END DO
!
!  Apply the time correlation matrix to the newly convolved adjoint
!  fields.
!
      DO ng=1,ngrids
        tlf(ng)%Rindex=0
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL time_corr (ng, tile, iadm, adm(ng)%IOtype, adm(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END DO
!
      DO ng=1,ngrids
        IF (nrec(ng).gt.3) THEN
#  ifdef PROFILE
          CALL wclock_on (ng, model, 82, __line__, myfile)
#  endif
          adm(ng)%Rindex=0
          lwrtstate2d(ng)=.true.
          DO irec=nrec(ng)-1,1,-1
!
!  Read the latest TLF solution. We need to read the TLF file in reverse
!  order since the final solution is written into the adjoint netcdf
!  file which will be rewritten in ascending order of time in the TLF
!  file by wrt_impulse.
!  NOTE: model=6 here so as to only read the state variables.
!
            adrec=irec
            CALL get_state (ng, 6, 6, tlf(ng), adrec, rold(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
            lweak=.true.
            add=.false.
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL tl_convolution (ng, tile, rold(ng), lweak, 2)
              CALL tl_variability (ng, tile, rold(ng), lweak)
#  ifdef BALANCE_OPERATOR
              CALL tl_balance (ng, tile, rini, rold(ng))
#  endif
              CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
            END DO
!
!  Overwrite ADM(ng)%name history NetCDF file with convolved adjoint
!  solution.
!
            kstp(ng)=rold(ng)
#  ifdef SOLVE3D
            nstp(ng)=rold(ng)
#  endif
#  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
          lwrtstate2d(ng)=.false.
          lwrttime(ng)=.true.
#  ifdef PROFILE
          CALL wclock_off (ng, model, 82, __line__, myfile)
#  endif
        END IF
      END DO
 
# else
!
!  If weak constraint, convolve records 2-Nrec in ADM(ng)%name
!  and impose model error covariance. NOTE: We will not use the
!  convolved forcing increments generated here since these arrays
!  do not contain the complete solution and are redundant.
!  AMM: We might want to get rid of these unwanted records to
!  avoid any confusion in the future.
!
      DO ng=1,ngrids
        IF (nrec(ng).gt.3) THEN
#  ifdef PROFILE
          CALL wclock_on (ng, model, 82, __line__, myfile)
#  endif
          DO irec=1,nrec(ng)-1
!
!  Read adjoint solution. Since routine "get_state" loads data into the
!  ghost points, the adjoint solution is read in the tangent linear
!  state arrays by using iTLM instead of iADM in the calling arguments.
!
            adrec=irec
            CALL get_state (ng, itlm, 4, adm(ng), adrec, rold(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Load interior solution, read above, into adjoint state arrays.
!  Then, multiply adjoint solution by the background-error standard
!  deviations. Next, convolve resulting adjoint solution with the
!  squared-root adjoint diffusion operator which impose the model-error
!  spatial correlations. Notice that the spatial convolution is only
!  done for half of the diffusion steps (squared-root filter). Clear
!  tangent linear state arrays when done.
!
            lweak=.true.
            add=.false.
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
#  ifdef BALANCE_OPERATOR
              CALL ad_balance (ng, tile, rini, rold(ng))
#  endif
              CALL ad_variability (ng, tile, rold(ng), lweak)
              CALL ad_convolution (ng, tile, rold(ng), lweak, 2)
              CALL initialize_ocean (ng, tile, itlm)
              CALL initialize_forces (ng, tile, itlm)
#  ifdef ADJUST_BOUNDARY
              CALL initialize_boundary (ng, tile, itlm)
#  endif
            END DO
!
!  To insure symmetry, convolve resulting filtered adjoint solution
!  from above with the squared-root (half of steps) tangent linear
!  diffusion operator. Then, multiply result with its corresponding
!  background-error standard deviations.  Since the convolved solution
!  is in the adjoint state arrays, first copy to tangent linear state
!  arrays including the ghosts points. Copy back to adjoint state
!  arrays when done with the convolution for output purposes.
!
            lweak=.true.
            add=.false.
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL load_adtotl (ng, tile, rold(ng), rold(ng), add)
              CALL tl_convolution (ng, tile, rold(ng), lweak, 2)
              CALL tl_variability (ng, tile, rold(ng), lweak)
#  ifdef BALANCE_OPERATOR
              CALL tl_balance (ng, tile, rini, rold(ng))
#  endif
              CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
            END DO
!
!  Overwrite ADM(ng)%name history NetCDF file with final solution.
!
            kstp(ng)=rold(ng)
#  ifdef SOLVE3D
            nstp(ng)=rold(ng)
#  endif
#  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
          lwrtstate2d(ng)=.false.
          lwrttime(ng)=.true.
#  ifdef PROFILE
          CALL wclock_off (ng, model, 82, __line__, myfile)
#  endif
        END IF
      END DO
# endif
!
      RETURN
      END SUBROUTINE error_covariance
 
# ifdef SP4DVAR
!
!***********************************************************************
      SUBROUTINE saddlec (driver, Lselect, Rini, Rold, Rnew)
!***********************************************************************
!
!  Imported variable declarations.
!
      integer, intent(in) :: rini
      integer, intent(in) :: rold(ngrids)
      integer, intent(in) :: rnew(ngrids)
!
      logical, intent(in) :: lselect
!
      character (len=*), intent(in) :: driver
!
!  Local variable declarations.
!
      logical :: lweak, add
!
      integer :: inirec
      integer :: ng, tile
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", saddlec"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Specify spatial error covariance on the TLM control vector, at time
!  index "Rold", via convolutions of a pseudo-diffusion operator.  The
!  convolved control vector is loaded into time index "Rnew" and index
!  "Rold" is preserved.
!-----------------------------------------------------------------------
 
# ifdef BALANCE_OPERATOR
!
!  Read NLM initial condition in readiness for the balance operator.
!
      inirec=rini
      DO ng=1,ngrids
        CALL get_state (ng, inlm, 2, ini(ng), inirec, inirec)
        IF (exit_flag.ne.noerror) RETURN
        nrhs(ng)=rini
      END DO
# endif
!
!  Load "Rold" time index TLM control vector into ADM state arrays.
!  Apply balance operator, if activated. Then, scale control vector
!  with error covariance standard deviations. Next, convolve resulting
!  vector with the squared-root adjoint diffusion operator. Notice
!  that the spatial convolution is only done for half of the diffusion
!  steps (squared-root filter).
!
      lweak=lselect
      add=.false.
      DO ng=1,ngrids
# ifdef PROFILE
        CALL wclock_on (ng, iadm, 82, __line__, myfile)
# endif
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
# ifdef BALANCE_OPERATOR
          CALL ad_balance (ng, tile, rini, rold(ng))
# endif
          CALL ad_variability (ng, tile, rold(ng), lweak)
          CALL ad_convolution (ng, tile, rold(ng), lweak, 2)
! AMM Important.
          CALL initialize_ocean (ng, tile, itlm)
          CALL initialize_forces (ng, tile, itlm)
        END DO
# ifdef PROFILE
        CALL wclock_off (ng, iadm, 82, __line__, myfile)
# endif
      END DO
!
!  Since we wish to preserve what is in tl_var(Rold), load resulting
!  filtered solution from above, ad_var(Rold), into tl_var(Rnew).
!  Then, convolve with the squared-root (half of steps) tangent
!  linear diffusion operator. Next, scale results with error
!  covariance standard deviations. Apply balance operator, if
!  activated.
!
      add=.false.
      DO ng=1,ngrids
# ifdef PROFILE
        CALL wclock_on (ng, itlm, 82, __line__, myfile)
# endif
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL load_adtotl (ng, tile, rold(ng), rnew(ng), add)
          CALL tl_convolution (ng, tile, rnew(ng), lweak, 2)
          CALL tl_variability (ng, tile, rnew(ng), lweak)
# ifdef BALANCE_OPERATOR
          CALL tl_balance (ng, tile, rini, rnew(ng))
# endif
        END DO
# ifdef PROFILE
        CALL wclock_off (ng, itlm, 82, __line__, myfile)
# endif
      END DO
!
      RETURN
      END SUBROUTINE saddlec
# endif
#endif
      END MODULE convolve_mod