tl_congrad.F

Go to the documentation of this file.

#include "cppdefs.h"
 
#if (defined SENSITIVITY_4DVAR || \
     defined tl_rbl4dvar       || \
     defined tl_r4dvar)        && \
    !defined RPCG
 
      SUBROUTINE tl_congrad (ng, model, outLoop, innLoop, NinnLoop,     &
     &                       Lcgini)
!
!git $Id$
!=================================================== Andrew M. Moore ===
!  Copyright (c) 2002-2025 The ROMS Group           Hernan G. Arango   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  Weak Constraint 4-Dimensional Variational (4DVar) Pre-conditioned   !
!                     Conjugate Gradient Algorithm                     !
 
# if defined R4DVAR_ANA_SENSITIVITY || defined TL_R4DVAR
!                                                                      !
!  The indirect representer method solves the system:                  !
!                                                                      !
!              (R_n + Cobs) * Beta_n = h_n                             !
!                                                                      !
!              h_n = Xo - H * X_n                                      !
!                                                                      !
!  where R_n is the representer matrix, Cobs is the observation-error  !
!  covariance,  Beta_n  are the representer coefficients,  h_n is the  !
!  misfit between observations (Xo) and model (H*X_n),  and  H is the  !
!  linearized observation operator. Here, _n denotes iteration.        !
!                                                                      !
!  This system does not need to be solved explicitly by inverting the  !
!  symmetric stabilized representer matrix, P_n:                       !
!                                                                      !
!              P_n = R_n + Cobs                                        !
!                                                                      !
!  but by computing the action of P_n on any vector PSI, such that     !
!                                                                      !
!              P_n * PSI = R_n * PSI + Cobs * PSI                      !
!                                                                      !
!  The representer matrix is not explicitly computed but evaluated by  !
!  one integration backward of the adjoint model  and one integration  !
!  forward of the tangent linear model for any forcing vector PSI.     !
!                                                                      !
!  Notice that "ObsScale" vector is used for screenning observations.  !
!  This scale is one (zero) for good (bad) observations.               !
!                                                                      !
!  Currently, parallelization of this algorithm is not needed because  !
!  each parallel node has a full copy of the assimilation vectors.     !
!                                                                      !
!  This code solves Ax=b by minimizing the cost function               !
!  0.5*x*A*x-x*b, assuming an initial guess of x=x0. In this case the  !
!  gradient is Ax-b and the Hessian is A.                              !
!                                                                      !
!  Reference:                                                          !
!                                                                      !
!    Chua, B. S. and A. F. Bennett,  2001:  An inverse ocean modeling  !
!      sytem, Ocean Modelling, 3, 137-165.                             !
 
# elif defined RBL4DVAR_ANA_SENSITIVITY || defined TL_RBL4DVAR
!                                                                      !
!  Solve the system (following Courtier, 1997):                        !
!                                                                      !
!              (H M_n B (M_n)' H' + Cobs) * w_n = d_n                  !
!                                                                      !
!              d_n = yo - H * Xb_n                                     !
!                                                                      !
!  where  M_n is the tangent linear model matrix and  _n denotes a     !
!  sequence of outer-loop estimates, Cobs is the observation-error     !
!  covariance,  B is the background error  covariance,  d_n is the     !
!  misfit between observations (yo) and model (H * Xb_n), and H is     !
!  the linearized observation operator.                                !
!                                                                      !
!  The analysis increment is:                                          !
!                                                                      !
!             dx_n = B M' H' w_n                                       !
!                                                                      !
!  so that Xa = Xb + dx_n.                                             !
!                                                                      !
!  The system does not need to be  solved explicitly  by inverting     !
!  the symmetric matrix, P_n:                                          !
!                                                                      !
!              P_n = H M_n B (M_n)' H' + Cobs                          !
!                                                                      !
!  but by computing the action of P_n on any vector PSI, such that     !
!                                                                      !
!              P_n * PSI =  H M_n B (M_n)' H' * PSI + Cobs * PSI       !
!                                                                      !
!  The (H M_n B (M_n)' H') matrix is not  explicitly computed  but     !
!  evaluated by  one integration backward of the adjoint model and     !
!  one  integration  forward of the  tangent linear model  for any     !
!  forcing vector PSI.                                                 !
!                                                                      !
!  A preconditioned conjugate gradient algorithm is used to compute    !
!  an approximation PSI for w_n.                                       !
!                                                                      !
!  Reference:                                                          !
!                                                                      !
!    Courtier, P., 1997: Dual formulation of four-dimensional          !
!      variational assimilation, Quart. J. Roy. Meteor. Soc.,          !
!      123, 2449-2461.                                                 !
# endif
!                                                                      !
!  Lanczos Algorithm Reference:                                        !
!                                                                      !
!    Fisher, M., 1998: Minimization Algorithms for Variational Data    !
!      Assimilation. In Seminar on Recent Developments in Numerical    !
!      Methods for Atmospheric Modelling, 1998.                        !
!                                                                      !
!    Tchimanga, J., S. Gratton, A.T. Weaver, and A. Sartenaer, 2008:   !
!      Limited-memory preconditioners, with application to incremental !
!      four-dimensional variational ocean data assimilation, Q.J.R.    !
!      Meteorol. Soc., 134, 753-771.                                   !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_fourdvar
      USE mod_iounits
      USE mod_scalars
 
# ifdef DISTRIBUTE
!
      USE distribute_mod, ONLY : mp_bcastf, mp_bcasti, mp_bcastl
# endif
      implicit none
!
!  Imported variable declarations
!
      logical, intent(in) :: Lcgini
      integer, intent(in) :: ng, model, outLoop, innLoop, NinnLoop
!
!  Local variable declarations.
!
      logical :: Ltrans
 
      integer :: i, j, iobs, ivec, Lscale, info
 
      real(r8) :: dla, zbet
      real(r8) :: tl_dla
# ifdef MINRES
      real(r8) :: zsum, zck, zgk
      real(r8) :: tl_zsum, tl_zck, tl_zgk
# endif
 
      real(r8), dimension(NinnLoop) :: zu, zgam
      real(r8), dimension(NinnLoop) :: tl_zu, tl_zrhs
      real(r8), dimension(Ndatum(ng)) :: pgrad, zt
      real(r8), dimension(Ndatum(ng)) :: tl_px, tl_pgrad, tl_zt
# ifdef MINRES
      real(r8), dimension(innLoop,innLoop) :: ztriT, zLT, zLTt
      real(r8), dimension(innLoop,innLoop) :: tl_ztriT, tl_zLT
      real(r8), dimension(innLoop,innLoop) :: tl_zLTt
      real(r8), dimension(innLoop) :: tau, zwork1, ze, zeref
      real(r8), dimension(innLoop) :: tl_tau, tl_zwork1, tl_ze, tl_zeref
# endif
!
!=======================================================================
!  Weak constraint 4DVar conjugate gradient, Lanczos-based algorithm.
!=======================================================================
!
!  This conjugate gradient algorithm is not run in parallel since the
!  weak constraint is done in observation space. Mostly all the
!  variables are 1D arrays. Therefore, in parallel applications (only
!  distributed-memory is possible) the master node does the computations
!  and then broadcasts results to remaining nodes in the communicator.
!
!  This version of congrad solves A(x+x0)=b for x, by minimizing
!  J=0.5*x'Ax-x'(b+Ax0), where x0 is a first-guess estimate of the
!  representer coefficients from the previous outer-loop.
!  For the first outer-loop, x0=0. In the code x0=cg_pxsave and
!  x=px.
!
      ltrans=.false.
 
      master_thread : IF (master) THEN
!
!  Initialize cg_Gnorm. The TL of precond is not available.
!
        DO i=1,outloop
          cg_gnorm(i)=cg_gnorm_v(i)
        END DO
!
!  Initialize internal parameters.  This is needed here for output
!  reasons.
!
        IF (innloop.eq.0) THEN
 
# if defined RBL4DVAR || defined TL_RBL4DVAR
!
!  If this is the first inner-loop, save NLmodVal in BCKmodVal.
!
          DO iobs=1,ndatum(ng)
            bckmodval(iobs)=nlmodval(iobs)
          END DO
# endif
!
!  If this is the first outer-loop, clear the solution vector px.
!
          IF ((outloop.eq.1).or.(.not.lhotstart)) THEN
!
!  For the first outer-loop, x0=0.
!
            DO iobs=1,ndatum(ng)
              tl_px(iobs)=0.0_r8
              tl_cg_pxsave(iobs)=0.0_r8
            END DO
!
!  If this is the first pass of the inner loop, set up the vectors and
!  store the first guess. The initial starting guess is assumed to be
!  zero in which case the gradient is just: -(obs-model).
!  A first-level preconditioning is applied using the inverse of the
!  observation error standard deviations (i.e. sqrt(ObsErr)).
!
            DO iobs=1,ndatum(ng)
# if defined RBL4DVAR || defined TL_RBL4DVAR
!^            pgrad(iobs)=-ObsScale(iobs)*                              &
!^   &                    (ObsVal(iobs)-BCKmodVal(iobs))
!^
              tl_pgrad(iobs)=-obsscale(iobs)*tl_obsval(iobs)
# else
!^            pgrad(iobs)=-ObsScale(iobs)*                              &
!^   &                    (ObsVal(iobs)-TLmodVal(iobs))
!<>           tl_pgrad(iobs)=-ObsScale(iobs)*                           &
!<>  &                       (tl_ObsVal(iobs)-tl_TLmodVal(iobs))
!^
              tl_pgrad(iobs)=-obsscale(iobs)*                           &
     &                       (tl_obsval(iobs)-tlmodval(iobs))
# endif
!
! Convert pgrad from x-space to v-space.
!
              IF (obserr(iobs).NE.0.0_r8) THEN
!^              pgrad(iobs)=pgrad(iobs)/SQRT(ObsErr(iobs))
!^
                tl_pgrad(iobs)=tl_pgrad(iobs)/sqrt(obserr(iobs))
              END IF
!^            vgrad0(iobs)=pgrad(iobs)
!^
            END DO
!
!  If preconditioning, convert pgrad from v-space to y-space.
!
!^          IF (Lprecond.and.(outLoop.gt.1)) THEN
!^            Lscale=2                 ! SQRT spectral LMP
!^            Ltrans=.TRUE.
!^            CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop,     &
!^   &                       pgrad)
!^          END IF
!^
!^          cg_Gnorm(outLoop)=0.0_r8
!^
            tl_cg_gnorm(outloop)=0.0_r8
!^          vgnorm(outLoop)=0.0_r8
!^
            DO iobs=1,ndatum(ng)
!^            zgrad0(iobs,outLoop)=pgrad(iobs)
!^
              tl_zgrad0(iobs)=tl_pgrad(iobs)
!^            cg_Gnorm(outLoop)=cg_Gnorm(outLoop)+                      &
!^   &                          zgrad0(iobs)*zgrad0(iobs)
!^
              tl_cg_gnorm(outloop)=tl_cg_gnorm(outloop)+                &
     &                             2.0_r8*tl_zgrad0(iobs)*              &
     &                             zgrad0(iobs,outloop)
!^            vgnorm(outLoop)=vgnorm(outLoop)+vgrad0(iobs)*vgrad0(iobs)
!^
            END DO
!^          cg_Gnorm(outLoop)=SQRT(cg_Gnorm(outLoop))
!^
            tl_cg_gnorm(outloop)=0.5_r8*tl_cg_gnorm(outloop)/           &
     &                           cg_gnorm(outloop)
!^          vgnorm(outLoop)=SQRT(vgnorm(outLoop))
!^
            DO iobs=1,ndatum(ng)
!^            zcglwk(iobs,1,outLoop)=pgrad(iobs)/cg_Gnorm(outLoop)
!^
              tl_zcglwk(iobs,1)=(tl_pgrad(iobs)-                        &
     &                           tl_cg_gnorm(outloop)*                  &
     &                           zcglwk(iobs,1,outloop))/               &
     &                          cg_gnorm(outloop)
!^            ADmodVal(iobs)=zcglwk(iobs,1,outLoop)
!<>           tl_ADmodVal(iobs)=tl_zcglwk(iobs,1)
!^
              admodval(iobs)=tl_zcglwk(iobs,1)
            END DO
!
!  If preconditioning, convert ADmodVal from y-space to v-space.
!
!^          IF (Lprecond.and.(outLoop.gt.1)) THEN
!^            Lscale=2                 ! SQRT spectral LMP
!^            Ltrans=.FALSE.
!^            CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop,     &
!^   &                       ADmodVal)
!^          END IF
!
! Convert ADmodVal from v-space to x-space.
!
            DO iobs=1,ndatum(ng)
              IF (obserr(iobs).NE.0.0_r8) THEN
!^              ADmodVal(iobs)=ADmodVal(iobs)/SQRT(ObsErr(iobs))
!<>             tl_ADmodVal(iobs)=tl_ADmodVal(iobs)/SQRT(ObsErr(iobs))
!^
                admodval(iobs)=admodval(iobs)/sqrt(obserr(iobs))
              END IF
            END DO
 
!^          cg_QG(1,outLoop)=0.0_r8
!^
            tl_cg_qg(1)=0.0_r8
            DO iobs=1,ndatum(ng)
!^            cg_QG(1,outLoop)=cg_QG(1,outLoop)+                        &
!^   &                         zcglwk(iobs,1,outLoop)*zgrad0(iobs)
!^
              tl_cg_qg(1)=tl_cg_qg(1)+                                  &
     &                    tl_zcglwk(iobs,1)*zgrad0(iobs,outloop)+       &
     &                    zcglwk(iobs,1,outloop)*tl_zgrad0(iobs)
            END DO
 
          ELSE
            IF (lcgini) THEN
!
!  When outer>1 we need to evaluate Ax0 so for inner=0 we use
!  cg_pxsave in v-space as the starting vector.
!
              DO iobs=1,ndatum(ng)
!^              ADmodVal(iobs)=cg_pxsave(iobs)
!<>             tl_ADmodVal(iobs)=tl_cg_pxsave(iobs)
!^
                admodval(iobs)=tl_cg_pxsave(iobs)
# if defined RBL4DVAR || defined TL_RBL4DVAR
!^              cg_innov(iobs)=-ObsScale(iobs)*                         &
!^   &                         (ObsVal(iobs)-BCKmodVal(iobs))
!<>             tl_cg_innov(iobs)=0.0_r8
!^
                tl_cg_innov(iobs)=-obsscale(iobs)*tl_obsval(iobs)
# else
!^              cg_innov(iobs)=-ObsScale(iobs)*                         &
!^   &                  (ObsVal(iobs)-TLmodVal(iobs))
!<>             tl_cg_innov(iobs)=-ObsScale(iobs)*                      &
!<>  &                            (tl_ObsVal(iobs)-tl_TLmodVal(iobs))
!^
                tl_cg_innov(iobs)=-obsscale(iobs)*                      &
     &                            (tl_obsval(iobs)-tlmodval(iobs))
# endif
              END DO
!
!  Convert ADmodVal from v-space to x-space and cg_innov (the
!  contribution to the starting gradient) from x-space to v-space.
!
              DO iobs=1,ndatum(ng)
                IF (obserr(iobs).NE.0.0_r8) THEN
!^                ADmodVal(iobs)=ADmodVal(iobs)/SQRT(ObsErr(iobs))
!<>               tl_ADmodVal(iobs)=tl_ADmodVal(iobs)/SQRT(ObsErr(iobs))
!^
                  admodval(iobs)=admodval(iobs)/sqrt(obserr(iobs))
!^                cg_innov(iobs)=cg_innov(iobs)/SQRT(ObsErr(iobs))
!^
                  tl_cg_innov(iobs)=tl_cg_innov(iobs)/sqrt(obserr(iobs))
                END IF
              END DO
 
            ELSE
 
              DO iobs=1,ndatum(ng)
!
!  Convert gradient contribution from x-space to v-space.
!
!^              pgrad(iobs)=ObsScale(iobs)*TLmodVal(iobs)
!^              tl_pgrad(iobs)=ObsScale(iobs)*tl_TLmodVal(iobs)
!^
                tl_pgrad(iobs)=obsscale(iobs)*tlmodval(iobs)
                IF (obserr(iobs).NE.0.0_r8) THEN
!^                pgrad(iobs)=pgrad(iobs)/SQRT(ObsErr(iobs))
!^
                  tl_pgrad(iobs)=tl_pgrad(iobs)/sqrt(obserr(iobs))
                END IF
!
!  Add I*x0=cg_pxsave contribution to the gradient and the term
!  -b=cg_innov (everything is in v-space at this point).
!
!^              pgrad(iobs)=pgrad(iobs)+ObsScale(iobs)*                 &
!^   &                      (cg_pxsave(iobs)+cg_innov(iobs))
!^
                tl_pgrad(iobs)=tl_pgrad(iobs)+obsscale(iobs)*           &
     &                         (tl_cg_pxsave(iobs)+tl_cg_innov(iobs))
!^              vgrad0(iobs)=pgrad(iobs)
!^
              END DO
!
!  If preconditioning, convert pgrad from v-space to y-space.
!
!^            IF (Lprecond.and.(outLoop.gt.1)) THEN
!^              Lscale=2                 ! SQRT spectral LMP
!^              Ltrans=.TRUE.
!^              CALL rprecond(ng, Lscale, Ltrans, outLoop, NinnLoop,    &
!^   &                        pgrad)
!^            END IF
!^
!^            cg_Gnorm(outLoop)=0.0_r8
!^
              tl_cg_gnorm(outloop)=0.0_r8
!^            vgnorm(outLoop)=0.0_r8
!^
              DO iobs=1,ndatum(ng)
!^              zgrad0(iobs,outLoop)=pgrad(iobs)
!^
                tl_zgrad0(iobs)=tl_pgrad(iobs)
!^              cg_Gnorm(outLoop)=cg_Gnorm(outLoop)+                    &
!^   &                            zgrad0(iobs,outLoop)*                 &
!^   &                            zgrad0(iobs,outLoop)
!^
                tl_cg_gnorm(outloop)=tl_cg_gnorm(outloop)+              &
     &                               2.0_r8*tl_zgrad0(iobs)*            &
     &                               zgrad0(iobs,outloop)
!^              vgnorm(outLoop)=vgnorm(outLoop)+                        &
!^   &                          vgrad0(iobs)*vgrad0(iobs)
!^
              END DO
!^            cg_Gnorm(outLoop)=SQRT(cg_Gnorm(outLoop))
!^
              tl_cg_gnorm(outloop)=0.5_r8*tl_cg_gnorm(outloop)/         &
     &                             cg_gnorm(outloop)
!^            vgnorm(outLoop)=SQRT(vgnorm(outLoop))
!^
              DO iobs=1,ndatum(ng)
!^              zcglwk(iobs,1,outLoop)=pgrad(iobs)/cg_Gnorm(outLoop)
!^
                tl_zcglwk(iobs,1)=(tl_pgrad(iobs)-                      &
     &                             tl_cg_gnorm(outloop)*                &
     &                             zcglwk(iobs,1,outloop))/             &
     &                            cg_gnorm(outloop)
!^              ADmodVal(iobs)=zcglwk(iobs,1,outLoop)
!<>             tl_ADmodVal(iobs)=tl_zcglwk(iobs,1)
!^
                admodval(iobs)=tl_zcglwk(iobs,1)
              END DO
!
!  If preconditioning, convert ADmodVal from y-space to v-space.
!
!^            IF (Lprecond.and.(outLoop.gt.1)) THEN
!^              Lscale=2                 ! SQRT spectral LMP
!^              Ltrans=.FALSE.
!^              CALL rprecond(ng, Lscale, Ltrans, outLoop, NinnLoop,    &
!^   &                        ADmodVal)
!^            END IF
!^
              DO iobs=1,ndatum(ng)
                IF (obserr(iobs).NE.0.0_r8) THEN
!^                ADmodVal(iobs)=ADmodVal(iobs)/SQRT(ObsErr(iobs))
!<>               tl_ADmodVal(iobs)=tl_ADmodVal(iobs)/SQRT(ObsErr(iobs))
!^
                  admodval(iobs)=admodval(iobs)/sqrt(obserr(iobs))
                END IF
              END DO
!^            cg_QG(1,outLoop)=0.0_r8
!^
              tl_cg_qg(1)=0.0_r8
              DO iobs=1,ndatum(ng)
!^              cg_QG(1,outLoop)=cg_QG(1,outLoop)+                      &
!^   &                           zcglwk(iobs,1,outLoop)*                &
!^   &                           zgrad0(iobs,outLoop)
!^
                tl_cg_qg(1)=tl_cg_qg(1)+                                &
     &                      tl_zcglwk(iobs,1)*zgrad0(iobs,outloop)+     &
     &                      zcglwk(iobs,1,outloop)*tl_zgrad0(iobs)
              END DO
 
            END IF
 
          END IF
 
        ELSE
!
!  After the initialization, for every other inner loop, calculate a
!  new Lanczos vector, store it in the matrix, and update search.
!
          DO iobs=1,ndatum(ng)
            pgrad(iobs)=obsscale(iobs)*tlmodval_s(iobs,innloop,outloop)
!<>         tl_pgrad(iobs)=ObsScale(iobs)*tl_TLmodVal(iobs)
            tl_pgrad(iobs)=obsscale(iobs)*tlmodval(iobs)
!
!  Convert gradient contribution from x-space to v-space.
!
            IF (obserr(iobs).NE.0.0_r8) THEN
              pgrad(iobs)=pgrad(iobs)/sqrt(obserr(iobs))
              tl_pgrad(iobs)=tl_pgrad(iobs)/sqrt(obserr(iobs))
            END IF
          END DO
 
          DO iobs=1,ndatum(ng)
            zt(iobs)=zcglwk(iobs,innloop,outloop)
            tl_zt(iobs)=tl_zcglwk(iobs,innloop)
          END DO
!
!  If preconditioning, convert zt from y-space to v-space.
!
!^        IF (Lprecond.and.(outLoop.gt.1)) THEN
!^          Lscale=2                 ! SQRT spectral LMP
!^          Ltrans=.FALSE.
!^          CALL rprecond(ng, Lscale, Ltrans, outLoop, NinnLoop, zt)
!^        END IF
!^
          DO iobs=1,ndatum(ng)
            pgrad(iobs)=pgrad(iobs)+obsscale(iobs)*zt(iobs)
            tl_pgrad(iobs)=tl_pgrad(iobs)+obsscale(iobs)*tl_zt(iobs)
          END DO
!
!  If preconditioning, convert pgrad from v-space to y-space.
!
!^        IF (Lprecond.and.(outLoop.gt.1)) THEN
!^          Lscale=2                 ! SQRT spectral LMP
!^          Ltrans=.TRUE.
!^          CALL rprecond(ng, Lscale, Ltrans, outLoop, NinnLoop, pgrad)
!^        END IF
!^
!^        cg_delta(innLoop,outLoop)=0.0_r8
!^
          tl_cg_delta(innloop)=0.0_r8
          DO iobs=1,ndatum(ng)
!^          cg_delta(innLoop,outLoop)=cg_delta(innLoop,outLoop)+        &
!^   &                                zcglwk(iobs,innLoop,outLoop)*     &
!^   &                                pgrad(iobs)
!^
            tl_cg_delta(innloop)=tl_cg_delta(innloop)+                  &
     &                           tl_zcglwk(iobs,innloop)*pgrad(iobs)+   &
     &                           zcglwk(iobs,innloop,outloop)*          &
     &                           tl_pgrad(iobs)
          END DO
!
!  Exit, if not a positive definite algorithm.
!
!^        IF (cg_delta(innLoop,outLoop).le.0.0_r8) THEN
!^          WRITE (stdout,20) cg_delta(innLoop,outLoop), outLoop,       &
!^   &                        innLoop
!^          exit_flag=8
!^        END IF
!^
!
!  Compute the new Lanczos vector.
!
          DO iobs=1,ndatum(ng)
            pgrad(iobs)=pgrad(iobs)-                                    &
                        cg_delta(innloop,outloop)*                      &
     &                  zcglwk(iobs,innloop,outloop)
            tl_pgrad(iobs)=tl_pgrad(iobs)-                              &
     &                     tl_cg_delta(innloop)*                        &
     &                     zcglwk(iobs,innloop,outloop)-                &
     &                     cg_delta(innloop,outloop)*                   &
     &                     tl_zcglwk(iobs,innloop)
          END DO
          IF (innloop.gt.1) THEN
            DO iobs=1,ndatum(ng)
              pgrad(iobs)=pgrad(iobs)-                                  &
     &                    cg_beta(innloop,outloop)*                     &
     &                    zcglwk(iobs,innloop-1,outloop)
              tl_pgrad(iobs)=tl_pgrad(iobs)-                            &
     &                       tl_cg_beta(innloop)*                       &
     &                       zcglwk(iobs,innloop-1,outloop)-            &
     &                       cg_beta(innloop,outloop)*                  &
     &                       tl_zcglwk(iobs,innloop-1)
            END DO
          END IF
!
!  Orthonormalize against previous Lanczos vectors.
!
          DO ivec=innloop,1,-1
!^          cg_dla(ivec,outLoop)=0.0_r8
!^
            tl_dla=0.0_r8
            DO iobs=1,ndatum(ng)
!^            cg_dla(ivec,outLoop)=cg_dla(ivec,outLoop)+                &
!^   &                             pgrad(iobs)*                         &
!^   &                             zcglwk(iobs,ivec,outLoop)
!^
              tl_dla=tl_dla+                                            &
     &               tl_pgrad(iobs)*zcglwk(iobs,ivec,outloop)+          &
     &               pgrad(iobs)*tl_zcglwk(iobs,ivec)
            END DO
            DO iobs=1,ndatum(ng)
              pgrad(iobs)=pgrad(iobs)-                                  &
     &                    cg_dla(ivec,outloop)*                         &
     &                    zcglwk(iobs,ivec,outloop)
              tl_pgrad(iobs)=tl_pgrad(iobs)-                            &
     &                       cg_dla(ivec,outloop)*tl_zcglwk(iobs,ivec)- &
     &                       tl_dla*zcglwk(iobs,ivec,outloop)
            END DO
          END DO
!
!^        cg_beta(innLoop+1,outLoop)=0.0_r8
!^
          tl_cg_beta(innloop+1)=0.0_r8
          DO iobs=1,ndatum(ng)
!^          cg_beta(innLoop+1,outLoop)=cg_beta(innLoop+1,outLoop)+      &
!^   &                                 pgrad(iobs)*pgrad(iobs)
!^
            tl_cg_beta(innloop+1)=tl_cg_beta(innloop+1)+                &
     &                            2.0_r8*tl_pgrad(iobs)*pgrad(iobs)
          END DO
!^        cg_beta(innLoop+1,outLoop)=SQRT(cg_beta(innLoop+1,outLoop))
!^
          tl_cg_beta(innloop+1)=0.5_r8*tl_cg_beta(innloop+1)/           &
     &                          cg_beta(innloop+1,outloop)
!
          DO iobs=1,ndatum(ng)
!^          zcglwk(iobs,innLoop+1,outLoop)=pgrad(iobs)/                 &
!^   &                                     cg_beta(innLoop+1,outLoop)
!^
            tl_zcglwk(iobs,innloop+1)=(tl_pgrad(iobs)-                  &
     &                                 tl_cg_beta(innloop+1)*           &
     &                                 zcglwk(iobs,innloop+1,outloop))/ &
     &                                cg_beta(innloop+1,outloop)
          END DO
!
!^        cg_QG(innLoop+1,outLoop)=0.0_r8
!^
          tl_cg_qg(innloop+1)=0.0_r8
          DO iobs=1,ndatum(ng)
!^          cg_QG(innLoop+1,outLoop)=cg_QG(innLoop+1,outLoop)+          &
!^   &                               zcglwk(iobs,innLoop+1,outLoop)*
!^   &                               zgrad0(iobs)
!^
            tl_cg_qg(innloop+1)=tl_cg_qg(innloop+1)+                    &
     &                          tl_zcglwk(iobs,innloop+1)*              &
     &                          zgrad0(iobs,outloop)+                   &
     &                          zcglwk(iobs,innloop+1,outloop)*         &
     &                          tl_zgrad0(iobs)
          END DO
          IF (innloop.eq.ninnloop) THEN
# ifdef MINRES
!
!  Use the minimum residual method as described by Paige and Saunders
!  ("Sparse Indefinite Systems of Linear Equations", 1975, SIAM Journal
!  on Numerical Analysis, 617-619). Specifically we refer to equations
!  6.10 and 6.11 of this paper.
!
!  Perform a LQ factorization of the tridiagonal matrix.
!
          ztrit=0.0_r8
          tl_ztrit=0.0_r8
          DO i=1,innloop
            ztrit(i,i)=cg_delta(i,outloop)
            tl_ztrit(i,i)=tl_cg_delta(i)
          END DO
          DO i=1,innloop-1
            ztrit(i,i+1)=cg_beta(i+1,outloop)
            tl_ztrit(i,i+1)=tl_cg_beta(i+1)
          END DO
          DO i=2,innloop
            ztrit(i,i-1)=cg_beta(i,outloop)
            tl_ztrit(i,i-1)=tl_cg_beta(i)
          END DO
!
!  Note: tl_sqlq also computes the LQ factorization of ztriT.
!
          CALL tl_sqlq(innloop, ztrit, tl_ztrit, tau, tl_tau, zwork1,   &
     &                 tl_zwork1)
!
!   Isolate L=zLT and its transpose.
!
          zlt=0.0_r8
          tl_zlt=0.0_r8
          zltt=0.0_r8
          tl_zltt=0.0_r8
          DO i=1,innloop
            DO j=1,i
              zlt(i,j)=ztrit(i,j)
              tl_zlt(i,j)=tl_ztrit(i,j)
            END DO
          END DO
          DO j=1,innloop
            DO i=1,innloop
              zltt(i,j)=zlt(j,i)
              tl_zltt(i,j)=tl_zlt(j,i)
            END DO
          END DO
!
!   Now form ze=beta1*Q*e1.
!
          ze=0.0_r8
          tl_ze=0.0_r8
          DO i=1,innloop
            ze(i)=-cg_qg(i,outloop)
            tl_ze(i)=-tl_cg_qg(i)
          END DO
          DO i=1,innloop
            DO j=1,innloop
              zeref(j)=0.0_r8
              tl_zeref(j)=0.0_r8
            END DO
            zeref(i)=1.0_r8
            tl_zeref(i)=0.0_r8
            DO j=i+1,innloop
              zeref(j)=ztrit(i,j)
              tl_zeref(j)=tl_ztrit(i,j)
            END DO
            zsum=0.0_r8
            tl_zsum=0.0_r8
            DO j=1,innloop
              zsum=zsum+ze(j)*zeref(j)
              tl_zsum=tl_zsum+tl_ze(j)*zeref(j)+ze(j)*tl_zeref(j)
            END DO
            DO j=1,innloop
              ze(j)=ze(j)-tau(i)*zsum*zeref(j)
              tl_ze(j)=tl_ze(j)-tl_tau(i)*zsum*zeref(j)-                &
     &                          tau(i)*tl_zsum*zeref(j)-                &
     &                          tau(i)*zsum*tl_zeref(j)
            END DO
          END DO
!
!   Now form ze=D*ze (using equation 5.6 and 6.5 also).
!
          zgk=sqrt(zlt(innloop,innloop)*zlt(innloop,innloop)+           &
     &        cg_beta(innloop+1,outloop)*cg_beta(innloop+1,outloop))
          IF (zgk.GT.0.0_r8) THEN
            tl_zgk=(tl_zlt(innloop,innloop)*zlt(innloop,innloop)+       &
     &              tl_cg_beta(innloop+1)*cg_beta(innloop+1,outloop))/  &
     &             zgk
          ELSE
            tl_zgk=0.0_r8
          ENDIF
          zck=zlt(innloop,innloop)/zgk
          tl_zck=tl_zlt(innloop,innloop)/zgk-tl_zgk*zck/zgk
          ze(innloop)=zck*ze(innloop)
          tl_ze(innloop)=tl_zck*ze(innloop)+zck*tl_ze(innloop)
!
!   Now compute tl_ze=inv(L')*(tl_ze-tl_L'*ze).
!
!   First solve for ze=inv(L')*ze.
!
          DO j=innloop,1,-1
            ze(j)=ze(j)/zltt(j,j)
            DO i=1,j-1
              ze(i)=ze(i)-ze(j)*zltt(i,j)
            END DO
          END DO
!
!   Next compute tl_rhs=tl_L'*ze then subtract from tl_ze.
!
          DO i=1,innloop
            tl_zrhs(i)=0.0_r8
            DO j=1,innloop
              tl_zrhs(i)=tl_zrhs(i)+tl_zltt(i,j)*ze(j)
            END DO
            tl_ze(i)=tl_ze(i)-tl_zrhs(i)
          END DO
!
!   Now solve the linear triangular system.
!
          DO j=innloop,1,-1
            tl_ze(j)=tl_ze(j)/zltt(j,j)
            DO i=1,j-1
              tl_ze(i)=tl_ze(i)-tl_ze(j)*zltt(i,j)
            END DO
          END DO
!
!   Copy the solution ze into zu.
!
          DO i=1,innloop
            zu(i)=ze(i)
            tl_zu(i)=tl_ze(i)
          END DO
# else
!
!  Calculate the new solution based upon the new, orthonormalized
!  Lanczos vector. First, the tridiagonal system is solved by
!  decomposition and forward/backward substitution.
!
            IF (ninnloop.eq.1) THEN
              zu(1)=-cg_qg(1,outloop)/cg_delta(1,outloop)
              tl_zrhs(1)=-tl_cg_qg(1)-tl_cg_delta(1)*zu(1)
              tl_zu(1)=tl_zrhs(1)/cg_delta(1,outloop)
            ELSE
!
!  Compute zu first.
!
              zbet=cg_delta(1,outloop)
              zu(1)=-cg_qg(1,outloop)/zbet
              DO ivec=2,innloop
                zgam(ivec)=cg_beta(ivec,outloop)/zbet
                zbet=cg_delta(ivec,outloop)-                            &
     &               cg_beta(ivec,outloop)*zgam(ivec)
                zu(ivec)=(-cg_qg(ivec,outloop)-cg_beta(ivec,outloop)*   &
     &                    zu(ivec-1))/zbet
              END DO
              DO ivec=innloop-1,1,-1
                zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
              END DO
!
!  Now compute tl_zrhs.
!
              tl_zrhs(1)=-tl_cg_qg(1)-                                  &
     &                   tl_cg_delta(1)*zu(1)-                          &
     &                   tl_cg_beta(2)*zu(2)
              DO ivec=2,innloop-1
                tl_zrhs(ivec)=-tl_cg_qg(ivec)-                          &
     &                        tl_cg_beta(ivec)*zu(ivec-1)-              &
     &                        tl_cg_delta(ivec)*zu(ivec)-               &
     &                        tl_cg_beta(ivec+1)*zu(ivec+1)
              END DO
              tl_zrhs(innloop)=-tl_cg_qg(innloop)-                      &
     &                         tl_cg_beta(innloop)*zu(innloop-1)-       &
     &                         tl_cg_delta(innloop)*zu(innloop)
!
!  Now solve the TL tridiagonal system A*dx=b-dA*x
!
              zbet=cg_delta(1,outloop)
              tl_zu(1)=tl_zrhs(1)/zbet
              DO ivec=2,innloop
                zgam(ivec)=cg_beta(ivec,outloop)/zbet
                zbet=cg_delta(ivec,outloop)-                            &
     &               cg_beta(ivec,outloop)*zgam(ivec)
                tl_zu(ivec)=(tl_zrhs(ivec)-cg_beta(ivec,outloop)*       &
     &                       tl_zu(ivec-1))/zbet
              END DO
 
              DO ivec=innloop-1,1,-1
!^              zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
!^
                tl_zu(ivec)=tl_zu(ivec)-zgam(ivec+1)*tl_zu(ivec+1)
              END DO
            END IF
 
!^          DO iobs=1,Ndatum(ng)
!^            zw(iobs)=zgrad0(iobs)+                                    &
!^   &                 cg_beta(innLoop+1,outLoop)*                      &
!^   &                 zcglwk(iobs,innLoop+1,outLoop)*zwork(innLoop,3)
!^          END DO
# endif
 
            DO iobs=1,ndatum(ng)
!^            px(iobs)=0.0_r8
!^
              tl_px(iobs)=0.0_r8
              DO ivec=1,innloop
!^              px(iobs)=px(iobs)+                                      &
!^   &                   zcglwk(iobs,ivec,outLoop)*zu(ivec)
!^
                tl_px(iobs)=tl_px(iobs)+                                &
     &                      tl_zcglwk(iobs,ivec)*zu(ivec)+              &
     &                      zcglwk(iobs,ivec,outloop)*tl_zu(ivec)
!^              zw(iobs)=zw(iobs)-                                      &
!^   &                   zcglwk(iobs,ivec,outLoop)*cg_QG(ivec,outLoop)
!^
              END DO
            END DO
!
!  If preconditioning, convert px from y-space to v-space.
!  We will always keep px in v-space.
!
!^          IF (Lprecond.and.(outLoop.gt.1)) THEN
!^            Lscale=2                 ! SQRT spectral LMP
!^            Ltrans=.FALSE.
!^            CALL rprecond(ng, Lscale, Ltrans, outLoop, NinnLoop, px)
!^          END IF
!^
          END IF
!
!  Put the new trial solution into the adjoint vector for the next loop
!  Put the final solution into the adjoint vector when converged
!  of on the final inner-loop.
!
          IF ((innloop.eq.ninnloop)) THEN
!
!  Note: px is already in v-space so there is no need to convert
!  if preconditioning. cg_pxsave is also in v-space.
!
            DO iobs=1,ndatum(ng)
!^            ADmodVal(iobs)=px(iobs)
!<>           tl_ADmodVal(iobs)=tl_px(iobs)
!^
              admodval(iobs)=tl_px(iobs)
            END DO
            IF (lhotstart) THEN
              DO iobs=1,ndatum(ng)
!^              ADmodVal(iobs)=ADmodVal(iobs)+cg_pxsave(iobs)
!<>             tl_ADmodVal(iobs)=tl_ADmodVal(iobs)+tl_cg_pxsave(iobs)
!^
                admodval(iobs)=admodval(iobs)+tl_cg_pxsave(iobs)
              END DO
              DO iobs=1,ndatum(ng)
!^              cg_pxsave(iobs)=ADmodVal(iobs)
!<>             tl_cg_pxsave(iobs)=tl_ADmodVal(iobs)
!^
                tl_cg_pxsave(iobs)=admodval(iobs)
              END DO
            END IF
          ELSE
            DO iobs=1,ndatum(ng)
!^            ADmodVal(iobs)=zcglwk(iobs,innLoop+1,outLoop)
!<>           tl_ADmodVal(iobs)=tl_zcglwk(iobs,innLoop+1)
!^
              admodval(iobs)=tl_zcglwk(iobs,innloop+1)
            END DO
!
!  If preconditioning, convert ADmodVal from y-space to v-space.
!
!^          IF (Lprecond.and.(outLoop.gt.1)) THEN
!^            Lscale=2                 ! SQRT spectral LMP
!^            Ltrans=.FALSE.
!^            CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop,     &
!^   &                       ADmodVal)
!^          END IF
!^
          END IF
!
!  Convert ADmodVal from v-space to x-space.
!
          DO iobs=1,ndatum(ng)
            IF (obserr(iobs).NE.0.0_r8) THEN
!^            ADmodVal(iobs)=ADmodVal(iobs)/SQRT(ObsErr(iobs))
!<>           tl_ADmodVal(iobs)=tl_ADmodVal(iobs)/SQRT(ObsErr(iobs))
!^
              admodval(iobs)=admodval(iobs)/sqrt(obserr(iobs))
            END IF
          END DO
        END IF
 
      END IF master_thread
 
# ifdef DISTRIBUTE
!
!  Broadcast new solution to other nodes.
!
      CALL mp_bcasti (ng, model, exit_flag)
      CALL mp_bcastf (ng, model, admodval)
      CALL mp_bcastf (ng, model, tl_cg_qg(:))
      CALL mp_bcastf (ng, model, tl_cg_gnorm(:))
      CALL mp_bcastf (ng, model, tl_cg_pxsave(:))
      CALL mp_bcastf (ng, model, tl_cg_innov(:))
      CALL mp_bcastf (ng, model, tl_cg_beta(:))
      CALL mp_bcastf (ng, model, tl_cg_delta(:))
      CALL mp_bcastf (ng, model, tl_zcglwk(:,:))
# endif
!
      RETURN
      END SUBROUTINE tl_congrad
#else
      SUBROUTINE tl_congrad
      RETURN
      END SUBROUTINE tl_congrad
#endif