#include "cppdefs.h"

#if defined RPCG && (defined SENSITIVITY_4DVAR || defined TL_W4DPSAS || \
    defined TL_W4DVAR)
      SUBROUTINE ad_rpcg_lanczos (ng, model, outLoop, innLoop,          &
     &                            NinnLoop, Lcgini)
!
!svn $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2020 The ROMS/TOMS Group       Selime Gurol      !
!    Licensed under a MIT/X style license            Andrew M. Moore   !
!    See License_ROMS.txt                                              !
!=======================================================================
!                                                                      !
!  Weak Constraint 4-Dimensional Variational (4DVar) Restricted        !
!      Pre-conditioned Conjugate Gradient (RPCG) Algorithm             !

# if defined TL_W4DVAR || defined W4DVAR_SENSITIVITY
!                                                                      !
!  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 stabilized representer matrix is A.        !
!                                                                      !
!  Reference:                                                          !
!                                                                      !
!    Chua, B. S. and A. F. Bennett,  2001:  An inverse ocean modeling  !
!      sytem, Ocean Modelling, 3, 137-165.                             !

# elif defined TL_W4DPSAS || defined W4DPSAS_SENSITIVITY || \
       defined W4DPSAS_FCT_SENSITIVITY
!                                                                      !
!  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.                                                 !
!                                                                      !
!  Reference:                                                          !
!                                                                      !
!    Courtier, P., 1997: Dual formulation of four-dimensional          !
!      variational assimilation, Quart. J. Roy. Meteor. Soc.,          !
!      123, 2449-2461.                                                 !
# endif
!                                                                      !
!  The restricted preconditioned conjugate gradient (RPCG) algorithm   !
!  is used to compute the solution. Specfically, the following system  !
!  is solved:                                                          !
!                                                                      !
!              (Cobs^-1 H M_n B (M_n)' H'+I)*w=Cobs^-1 d_n             !
!                                                                      !
!  All inner-products are defined with respect to a norm involving     !
!  H M_n B (M_n)' H', and the convergence of J in model space is       !
!  guaranteed to converge at the same rate as I4DVAR in primal space.  !
!                                                                      !
!  RPCG Reference:                                                     !
!                                                                      !
!    Gratton, S. and J. Tshimanga, 2009: An observation-space          !
!      formulation of variational assimilation using a restricted      !
!      preconditioned conjugate gradient algorithm. QJRMS, 135,        !
!      1573-1585.
!                                                                      !
!=======================================================================
!
      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, Laug
      integer :: i, ic, j, iobs, ivec, Lscale, info

      real(r8) :: zbet, eps, preducv, preducy
      real(r8) :: cff, fact, ad_dla, adfac

      real(r8), dimension(NinnLoop) :: zu, zgam, zfact, ad_zfact
      real(r8), dimension(NinnLoop) :: ad_zu, ad_zrhs
      real(r8), dimension(Ndatum(ng)+1) :: px, pgrad, pgrad_S
      real(r8), dimension(Ndatum(ng)+1) :: ad_px, ad_pgrad
      real(r8), dimension(Ninner,3) :: zwork, ad_zwork

      character (len=13) :: string
!
!=======================================================================
!  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
!
!  Set vgrad0=zgrad0 for the appropriat outer-loop.
!
        DO iobs=1,Ndatum(ng)+1
          vgrad0(iobs)=zgrad0(iobs,outLoop)
        END DO
!
! Clear arrays.
!
        ADmodVal=0.0_r8
        ad_pgrad=0.0_r8
        ad_px=0.0_r8
        ad_zu=0.0_r8
        ad_zrhs=0.0_r8
        ad_zfact=0.0_r8
        ad_dla=0.0_r8
        ad_zwork=0.0_r8
!!
!!      IF (innLoop.eq.NinnLoop) THEN
!!        DO iobs=1,Ndatum(ng)
!!          BCKmodVal(iobs)=NLmodVal(iobs)
!!        END DO
!!      END IF
!!
# if defined ANL_EOFS || defined TRACE_EOFS
!
        IF (innLoop.eq.NinnLoop) THEN
!
!  Clear adjoint arrays when innLoop=NinnLoop.
!
          DO j=1,NinnLoop+1
            DO iobs=1,Ndatum(ng)+1
              ad_zcglwk(iobs,j)=0.0_r8
              ad_vcglwk(iobs,j)=0.0_r8
            END DO
          END DO
          DO iobs=1,Ndatum(ng)+1
            ad_zgrad0(iobs)=0.0_r8
            ad_vgrad0(iobs)=0.0_r8
            ad_vgrad0s(iobs)=0.0_r8
          END DO
          DO iobs=1,Ndatum(ng)
#  ifdef IMPACT_INNER
            ad_ObsVal(iobs,innLoop)=0.0_r8
#  else
            ad_ObsVal(iobs)=0.0_r8
#  endif
            ad_ObsErr(iobs)=0.0_r8
            ad_ObsScale(iobs)=0.0_r8
            ad_Hbk(iobs)=0.0_r8
            ad_BCKmodVal(iobs)=0.0_r8
          END DO
          DO j=1,NinnLoop+1
            ad_cg_beta(j)=0.0_r8
            ad_cg_QG(j)=0.0_r8
          END DO
          DO j=1,NinnLoop
            ad_cg_delta(j)=0.0_r8
          END DO
          ad_cg_Gnorm(outLoop)=0.0_r8
          ad_cg_Gnorm_v(outLoop)=0.0_r8
          ad_Jb0(0)=0.0_r8
          ad_Jb0(outer)=0.0_r8
        END IF
# endif
!
!-----------------------------------------------------------------------
!  Initialization step, innloop = 0.
!-----------------------------------------------------------------------
!
        MINIMIZE : IF (innLoop.eq.0) THEN

# if defined W4DPSAS || defined TL_W4DPSAS || \
     defined W4DPSAS_SENSITIVITY || defined W4DPSAS_FCT_SENSITIVITY
!
          DO iobs=1,Ndatum(ng)
            BCKmodVal(iobs)=NLmodVal(iobs)
          END DO
# endif
!
!  Selime code: ADmodVal=z1.
!
            DO iobs=1,Ndatum(ng)
!>            tl_ADmodVal(iobs)=tl_zgrad0(iobs,outLoop)
!>
              ad_zgrad0(iobs)=ad_zgrad0(iobs)+TLmodVal(iobs)
              TLmodVal(iobs)=0.0_r8
            END DO
!
!  Selime code: zgrad0=z1.
!
!  Selime code: vgrad0=r.
!
            DO iobs=1,Ndatum(ng)+1
!>            tl_zgrad0(iobs)=tl_pgrad(iobs)
!>
              ad_pgrad(iobs)=ad_pgrad(iobs)+ad_zgrad0(iobs)
              ad_zgrad0(iobs)=0.0_r8
            END DO
!
!  If preconditioning, convert pgrad from v-space to y-space.
!
!  Selime code: z1=G*r.
!
!>          IF (Lprecond.and.(outLoop.gt.1)) THEN
!>            Lscale=1                 ! Spectral LMP
!>            Laug=.TRUE.
!>            CALL rprecond (ng, Lscale, Laug, outLoop, NinnLoop-1,     &
!>   &                       pgrad)
!>          END IF
!
!  AUGMENTED
!  Augment the residual vector.
!
!>          tl_vgrad0s(Ndatum(ng)+1)=tl_vgrad0(Ndatum(ng)+1)
!>
            ad_vgrad0(Ndatum(ng)+1)=ad_vgrad0 (Ndatum(ng)+1)+           &
     &                              ad_vgrad0s(Ndatum(ng)+1)
            ad_vgrad0s(Ndatum(ng)+1)=0.0_r8
!>          tl_vgrad0(Ndatum(ng)+1)=tl_pgrad(Ndatum(ng)+1)
!>
            ad_pgrad(Ndatum(ng)+1)=ad_pgrad (Ndatum(ng)+1)+             &
     &                             ad_vgrad0(Ndatum(ng)+1)
            ad_vgrad0(Ndatum(ng)+1)=0.0_r8
!>          tl_pgrad(Ndatum(ng)+1)=0.0_r8
!>
            ad_pgrad(Ndatum(ng)+1)=0.0_r8
!
!  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 residual is just: (obs-model).
!
            DO iobs=1,Ndatum(ng)
# if defined W4DPSAS             || defined TL_W4DPSAS || \
     defined W4DPSAS_SENSITIVITY || defined W4DPSAS_FCT_SENSITIVITY
!>            pgrad(iobs)=ObsScale(iobs)*                               &
!>   &                    (ObsVal(iobs)-BCKmodVal(iobs))
# else
!
!  NOTE: THIS CODE IS NOT CORRECT YET FOR W4DVAR
!
!>            pgrad(iobs)=ObsScale(iobs)*                               &
!>   &                    (ObsVal(iobs)-                                &
!>                         TLmodVal_S(iobs,innLoop,outLoop))
# endif
!>            tl_vgrad0s(iobs)=tl_vgrad0(iobs)
!>
              ad_vgrad0(iobs)=ad_vgrad0(iobs)+ad_vgrad0s(iobs)
              ad_vgrad0s(iobs)=0.0_r8
!>            tl_vgrad0(iobs)=tl_pgrad(iobs)
!>
              ad_pgrad(iobs)=ad_pgrad(iobs)+ad_vgrad0(iobs)
              ad_vgrad0(iobs)=0.0_r8
!
              IF (ObsErr(iobs).NE.0.0_r8) THEN
!>               pgrad(iobs)=pgrad(iobs)/ObsErr(iobs)
!>
                 pgrad(iobs)=vgrad0(iobs)
!>               tl_pgrad(iobs)=tl_pgrad(iobs)/ObsErr(iobs)-            &
!>   &                          tl_ObsErr(iobs)*pgrad(iobs)/ObsErr(iobs)
!>
                 ad_ObsErr(iobs)=ad_ObsErr(iobs)-ad_pgrad(iobs)*        &
     &                           pgrad(iobs)/ObsErr(iobs)
                 ad_pgrad(iobs)=ad_pgrad(iobs)/ObsErr(iobs)
              END IF
# if defined W4DPSAS             || defined TL_W4DPSAS || \
     defined W4DPSAS_SENSITIVITY || defined W4DPSAS_FCT_SENSITIVITY
!>            tl_pgrad(iobs)=ObsScale(iobs)*                            &
!>   &                       (tl_ObsVal(iobs)-tl_BCKmodVal(iobs))+      &
!>   &                       tl_ObsScale(iobs)*                         &
!>   &                       (ObsVal(iobs)-BCKmodVal(iobs))
#  ifdef IMPACT_INNER
              ad_ObsVal(iobs,innLoop)=ad_ObsVal(iobs,innLoop)+          &
     &                                ObsScale(iobs)*ad_pgrad(iobs)
#  else
              ad_ObsVal(iobs)=ad_ObsVal(iobs)+                          &
     &                        ObsScale(iobs)*ad_pgrad(iobs)
#  endif
              ad_BCKmodVal(iobs)=ad_BCKmodVal(iobs)-                    &
     &                           ObsScale(iobs)*ad_pgrad(iobs)
              ad_ObsScale(iobs)=ad_ObsScale(iobs)+                      &
     &                          ad_pgrad(iobs)*                         &
     &                          (ObsVal(iobs)-BCKmodVal(iobs))
              ad_pgrad(iobs)=0.0_r8
# else
!  NOTE: THIS CODE IS NOT CORRECT YET FOR W4DVAR
!
!>            tl_pgrad(iobs)=ObsScale(iobs)*                            &
!>   &                       (tl_ObsVal(iobs)-tl_TLmodVal(iobs))
!>            tl_pgrad(iobs)=ObsScale(iobs)*                            &
!>   &                       (tl_ObsVal(iobs)-TLmodVal(iobs))+          &
!>   &                       tl_ObsScale(iobs)*                         &
!>   &                       (ObsVal(iobs)-                             &
!>   &                        TLmodVal_S(iobs,innLoop,outLoop))
!>
#  ifdef IMPACT_INNER
              ad_ObsVal(iobs,innLoop)=ad_ObsVal(iobs,innLoop)+          &
     &                                ObsScale(iobs)*ad_pgrad(iobs)
#  else
              ad_ObsVal(iobs)=ad_ObsVal(iobs)+                          &
     &                        ObsScale(iobs)*ad_pgrad(iobs)
#  endif
              ADmodVal(iobs)=ADmodVal(iobs)-                            &
     &                       ObsScale(iobs)*ad_pgrad(iobs)
              ad_ObsScale(iobs)=ad_ObsScale(iobs)+
     &                          ad_pgrad(iobs)*                         &
     &                          (ObsVal(iobs)-                          &
     &                           TLmodVal_S(iobs,innLoop,outLoop))
              ad_pgrad(iobs)=0.0_r8
# endif
            END DO
!
!  For the first outer-loop, x0=0.
!
          IF (outLoop.eq.1) THEN
            DO iobs=1,Ndatum(ng)+1
!>            tl_px(iobs)=0.0_r8
!>
              ad_px(iobs)=0.0_r8
            END DO
            DO iobs=1,Ndatum(ng)
!>           tl_Hbk(iobs)=0.0_r8
!>
             ad_Hbk(iobs)=0.0_r8
            END DO
          ELSE
!>          IF (Lprecond) THEN
!>            DO iobs=1,Ndatum(ng)+1
!>              FOURDVAR(ng)%cg_Dpxsum(iobs)=FOURDVAR(ng)%cg_pxsum(iobs)
!>            END DO
!>            Lscale=1                 ! Spectral LMP
!>            Laug=.FALSE.
!>            CALL rprecond (ng, Lscale, Laug, outLoop, NinnLoop-1,     &
!>     &                     FOURDVAR(ng)%cg_Dpxsum)
!>            DO iobs=1,Ndatum(ng)+1
!>              FOURDVAR(ng)%cg_Dpxsum(iobs)=                           &
!>     &                               -FOURDVAR(ng)%cg_Dpxsum(iobs)+   &
!>     &                                FOURDVAR(ng)%cg_pxsum(iobs)
!>            END DO
!>          END IF
!!
!!  NOTE: CODE WON'T WORK FOR W4DVAR AT THE MOMENT SINCE THIS IS
!!  THE WRONG TLmodVal HERE!!!
!!
            DO iobs=1,Ndatum(ng)
!>            tl_Hbk(iobs)=-tl_TLmodVal(iobs)
!>
              ADmodVal(iobs)=ADmodVal(iobs)-ad_Hbk(iobs)
              ad_Hbk(iobs)=0.0_r8
            END DO
          END IF
# if defined W4DPSAS_SENSITIVITY || defined W4DPSAS_FCT_SENSITIVITY
!
          DO iobs=1,Ndatum(ng)
!>          tl_BCKmodVal(iobs)=0.0_r8
!>
            ad_BCKmodVal(iobs)=0.0_r8
          END DO
# endif
!
!-----------------------------------------------------------------------
!  Minimization step, innLoop > 0.
!-----------------------------------------------------------------------
!
        ELSE
!
!  ADJOINT NEXT CALCULATION
!
!  Put the new trial solution into the adjoint vector for the next
!  loop or put the final solution into the adjoint vector 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)+1
!>            FOURDVAR(ng)%tl_cg_pxsave(iobs)=tl_px(iobs)
!>
              ad_px(iobs)=ad_px(iobs)+FOURDVAR(ng)%ad_cg_pxsave(iobs)
              FOURDVAR(ng)%ad_cg_pxsave(iobs)=0.0_r8
            END DO

            DO iobs=1,Ndatum(ng)
!>            tl_ADmodVal(iobs)=tl_px(iobs)
!>
              ad_px(iobs)=ad_px(iobs)+TLmodVal(iobs)
              TLmodVal(iobs)=0.0_r8
            END DO
          ELSE
!
!  Selime code: z1=ADmodVal.
!
            DO iobs=1,Ndatum(ng)
!>            tl_ADmodVal(iobs)=tl_vcglwk(iobs,innLoop+1)
              ad_vcglwk(iobs,innLoop+1)=ad_vcglwk(iobs,innLoop+1)+      &
     &                                  TLmodVal(iobs)
              TLmodVal(iobs)=0.0_r8
            END DO
!
          END IF
!
!  Selime code: z1=vcglwk.
!
!  AUGMENTED.
!
          DO iobs=1,Ndatum(ng)+1
!>          tl_vcglwk(iobs,innLoop+1)=tl_pgrad(iobs)
!>
            ad_pgrad(iobs)=ad_pgrad(iobs)+ad_vcglwk(iobs,innLoop+1)
            ad_vcglwk(iobs,innLoop+1)=0.0_r8
          END DO

!
!  If preconditioning, convert pgrad from v-space to y-space.
!
!
!  Selime code: z1=Gw.
!
!>        IF (Lprecond.and.(outLoop.gt.1)) THEN
!>          Lscale=1                 ! Spectral LMP
!>          Laug=.TRUE.
!>          CALL rprecond (ng, Lscale, Laug, outLoop, NinnLoop-1, pgrad)
!>        END IF
!
!  Save the non-normalized Lanczos vector. zcglwk is used as temporary
!  storage. The calculation will be completed at the start of the next
!  inner-loop.
!
!  Selime code: w=zcglwk.
!
!  AUGMENTED
!
          DO iobs=1,Ndatum(ng)+1
!>          tl_zcglwk(iobs,innLoop+1)=tl_pgrad(iobs)
!>
            ad_pgrad(iobs)=ad_pgrad(iobs)+ad_zcglwk(iobs,innLoop+1)
            ad_zcglwk(iobs,innLoop+1)=0.0_r8
          END DO
!
!  Selime code: cg_dla=pgrad*HBH'*v1=pgrad*TLmodVal_S/zfact
!
!  AUGMENTED
!
!  Compute basic state pgrad.
!
          IF (innLoop.eq.1) THEN
            fact=1.0_r8/cg_Gnorm_v(outLoop)
          ELSE
            fact=1.0_r8/cg_beta(innLoop,outLoop)
          END IF
          DO iobs=1,Ndatum(ng)
            pgrad(iobs)=TLmodVal_S(iobs,innLoop,outLoop)*fact+           &
     &         Hbk(iobs,outLoop)*vcglwk(Ndatum(ng)+1,innLoop,outLoop)
          END DO
          DO iobs=1,Ndatum(ng)
            IF (ObsErr(iobs).NE.0.0_r8) THEN
              pgrad(iobs)=pgrad(iobs)/ObsErr(iobs)
            END IF
          END DO
          pgrad(Ndatum(ng)+1)=0.0_r8
          DO iobs=1,Ndatum(ng)
            pgrad(iobs)=pgrad(iobs)+vcglwk(iobs,innLoop,outLoop)
          END DO
          pgrad(Ndatum(ng)+1)=pgrad(Ndatum(ng)+1)+                       &
     &                             vcglwk(Ndatum(ng)+1,innLoop,outLoop)
          IF (innLoop.gt.1) THEN
            DO iobs=1,Ndatum(ng)+1
              pgrad(iobs)=pgrad(iobs)-cg_beta(innLoop,outLoop)*         &
     &                                 zcglwk(iobs,innLoop-1,outLoop)
            END DO
          END IF
          DO iobs=1,Ndatum(ng)+1
            pgrad(iobs)=pgrad(iobs)-cg_delta(innLoop,outLoop)*          &
     &                  zcglwk(iobs,innLoop,outLoop)
          END DO
!
!  Save the intermediate value of pgrad.
!
          DO iobs=1,Ndatum(ng)+1
            pgrad_S(iobs)=pgrad(iobs)
          END DO
!
!  Orthonormalize against previous directions.
!
!  Identify the appropriate Lanczos vector normalization coefficient.
!
          DO ivec=1,innLoop
            IF (ivec.eq.1) THEN
              zfact(ivec)=1.0_r8/cg_Gnorm_v(outLoop)
            ELSE
              zfact(ivec)=1.0_r8/cg_beta(ivec,outLoop)
            ENDIF
          END DO
!
!  Selime code: cg_dla=pgrad*HBH'*v1=pgrad*TLmodVal_S/zfact
!
!  AUGMENTED

!>        DO ivec=innLoop,1,-1
!>
          DO ivec=1,innLoop
            DO iobs=1,Ndatum(ng)+1
              pgrad(iobs)=pgrad_S(iobs)
            END DO
            DO i=ivec+1,innLoop
              DO iobs=1,Ndatum(ng)+1
                pgrad(iobs)=pgrad(iobs)-                                &
     &                      cg_dla(i,outLoop)*                          &
     &                      vcglwk(iobs,i,outLoop)
              END DO
            END DO
            DO iobs=1,Ndatum(ng)+1
!>            tl_pgrad(iobs)=tl_pgrad(iobs)-tl_dla*                     &
!>   &                       vcglwk(iobs,ivec,outLoop)-                 &
!>   &                       cg_dla(ivec,outLoop)*                      &
!>   &                       tl_vcglwk(iobs,ivec)
!>
              ad_dla=ad_dla-ad_pgrad(iobs)*vcglwk(iobs,ivec,outLoop)
              ad_vcglwk(iobs,ivec)=ad_vcglwk(iobs,ivec)-                &
     &                             cg_dla(ivec,outLoop)*ad_pgrad(iobs)
            END DO
!>          tl_dla=tl_dla+                                              &
!>   &             tl_Jb0(outLoop-1)*                                   &
!>   &             vcglwk(Ndatum(ng)+1,ivec,outLoop)*                   &
!>   &             pgrad(Ndatum(ng)+1)+                                 &
!>   &             Jb0(outLoop-1)*tl_vcglwk(Ndatum(ng)+1,ivec)*         &
!>   &             pgrad(Ndatum(ng)+1)+                                 &
!>   &             Jb0(outLoop-1)*vcglwk(Ndatum(ng)+1,ivec,outLoop)*    &
!>   &             tl_pgrad(Ndatum(ng)+1)
!>
            ad_Jb0(outLoop-1)=ad_Jb0(outLoop-1)+                        &
     &                        vcglwk(Ndatum(ng)+1,ivec,outLoop)*        &
     &                        pgrad(Ndatum(ng)+1)*ad_dla
            ad_vcglwk(Ndatum(ng)+1,ivec)=ad_vcglwk(Ndatum(ng)+1,ivec)+  &
     &                                   Jb0(outLoop-1)*                &
     &                                   pgrad(Ndatum(ng)+1)*ad_dla
            ad_pgrad(Ndatum(ng)+1)=ad_pgrad(Ndatum(ng)+1)+              &
     &                             Jb0(outLoop-1)*                      &
     &                             vcglwk(Ndatum(ng)+1,ivec,outLoop)*   &
     &                             ad_dla
            DO iobs=1,Ndatum(ng)
!>            tl_dla=tl_dla+                                            &
!>   &               tl_pgrad(iobs)*                                    &
!>   &               TLmodVal_S(iobs,ivec,outLoop)*                     &
!>   &               zfact(ivec)+                                       &
!>   &               pgrad(iobs)*tl_TLmodVal_S(iobs,ivec,outLoop)*      &
!>   &               zfact(ivec)+                                       &
!>   &               pgrad(iobs)*TLmodVal_S(iobs,ivec,outLoop)*         &
!>   &               tl_zfact(ivec)+                                    &
!>   &               tl_pgrad(iobs)*                                    &
!>   &               Hbk(iobs,outLoop)*                                 &
!>   &               vcglwk(Ndatum(ng)+1,ivec,outLoop)+                 &
!>   &               pgrad(iobs)*tl_Hbk(iobs)*                          &
!>   &               vcglwk(Ndatum(ng)+1,ivec,outLoop)+                 &
!>   &               pgrad(iobs)*Hbk(iobs,outLoop)*                     &
!>   &               tl_vcglwk(Ndatum(ng)+1,ivec)+                      &
!>   &               tl_Hbk(iobs)*vcglwk(iobs,ivec,outLoop)*            &
!>   &               pgrad(Ndatum(ng)+1)+                               &
!>   &               Hbk(iobs,outLoop)*tl_vcglwk(iobs,ivec)*            &
!>   &               pgrad(Ndatum(ng)+1)+                               &
!>   &               Hbk(iobs,outLoop)*vcglwk(iobs,ivec,outLoop)*       &
!>   &               tl_pgrad(Ndatum(ng)+1)
!>
              ad_pgrad(iobs)=ad_pgrad(iobs)+                            &
     &                       TLmodVal_S(iobs,ivec,outLoop)*             &
     &                       zfact(ivec)*ad_dla
              ad_TLmodVal_S(iobs,ivec,outLoop)=                         &
     &                               ad_TLmodVal_S(iobs,ivec,outLoop)+  &
     &                               pgrad(iobs)*zfact(ivec)*ad_dla
              ad_zfact(ivec)=ad_zfact(ivec)+                            &
     &                       pgrad(iobs)*                               &
     &                       TLmodVal_S(iobs,ivec,outLoop)*ad_dla
              ad_pgrad(iobs)=ad_pgrad(iobs)+                            &
     &                       Hbk(iobs,outLoop)*                         &
     &                       vcglwk(Ndatum(ng)+1,ivec,outLoop)*ad_dla
              ad_Hbk(iobs)=ad_Hbk(iobs)+                                &
     &                     pgrad(iobs)*                                 &
     &                     vcglwk(Ndatum(ng)+1,ivec,outLoop)*ad_dla
              ad_vcglwk(Ndatum(ng)+1,ivec)=ad_vcglwk(Ndatum(ng)+1,ivec)+&
     &                                     pgrad(iobs)*                 &
     &                                     Hbk(iobs,outLoop)*ad_dla
              ad_Hbk(iobs)=ad_Hbk(iobs)+                                &
     &                     vcglwk(iobs,ivec,outLoop)*                   &
     &                     pgrad(Ndatum(ng)+1)*ad_dla
              ad_vcglwk(iobs,ivec)=ad_vcglwk(iobs,ivec)+                &
     &                             Hbk(iobs,outLoop)*                   &
     &                             pgrad(Ndatum(ng)+1)*ad_dla
              ad_pgrad(Ndatum(ng)+1)=ad_pgrad(Ndatum(ng)+1)+            &
     &                               Hbk(iobs,outLoop)*                 &
     &                               vcglwk(iobs,ivec,outLoop)*ad_dla
            END DO
!>          tl_dla=0.0_r8
!>
            ad_dla=0.0_r8
          END DO
!
          DO ivec=1,innLoop
            IF (ivec.eq.1) THEN
              zfact(ivec)=1.0_r8/cg_Gnorm_v(outLoop)
!>            tl_zfact(ivec)=-tl_Gnorm_v(outLoop)*zfact(ivec)/          &
!>   &                        cg_Gnorm_v(outLoop)
!>
              ad_cg_Gnorm_v(outLoop)=ad_cg_Gnorm_v(outLoop)-            &
     &                               ad_zfact(ivec)*zfact(ivec)/        &
     &                               cg_Gnorm_v(outLoop)
              ad_zfact(ivec)=0.0_r8
            ELSE
              zfact(ivec)=1.0_r8/cg_beta(ivec,outLoop)
!>            tl_zfact(ivec)=-tl_cg_beta(ivec)*zfact(ivec)/             &
!>   &                        cg_beta(ivec,outLoop)
              ad_cg_beta(ivec)=ad_cg_beta(ivec)-                        &
     &                         ad_zfact(ivec)*zfact(ivec)/              &
     &                         cg_beta(ivec,outLoop)
              ad_zfact(ivec)=0.0_r8
            ENDIF
          END DO
!
!  Selime code: w=w-alpha*v1
!
!  AUGMENTED.
!
          DO iobs=1,Ndatum(ng)+1
!>          tl_pgrad(iobs)=tl_pgrad(iobs)-                              &
!>   &                     tl_cg_delta(innLoop)*                        &
!>   &                     zcglwk(iobs,innLoop,outLoop)-                &
!>   &                     cg_delta(innLoop,outLoop)*                   &
!>   &                     tl_zcglwk(iobs,innLoop)
!>
            ad_cg_delta(innLoop)=ad_cg_delta(innLoop)-                  &
     &                           zcglwk(iobs,innLoop,outLoop)*          &
     &                           ad_pgrad(iobs)
            ad_zcglwk(iobs,innLoop)=ad_zcglwk(iobs,innLoop)-            &
     &                              cg_delta(innLoop,outLoop)*          &
     &                              ad_pgrad(iobs)
          END DO
!
!  Selime code: alpha=w'*t1=cg_delta
!
!  AUGMENTED.
!
!  Recompute basic state pgrad.
!
          IF (innLoop.eq.1) THEN
            fact=1.0_r8/cg_Gnorm_v(outLoop)
          ELSE
            fact=1.0_r8/cg_beta(innLoop,outLoop)
          END IF
          DO iobs=1,Ndatum(ng)
            pgrad(iobs)=TLmodVal_S(iobs,innLoop,outLoop)*fact+          &
     &                  Hbk(iobs,outLoop)*                              &
     &                  vcglwk(Ndatum(ng)+1,innLoop,outLoop)
          END DO
          DO iobs=1,Ndatum(ng)
            IF (ObsErr(iobs).NE.0.0_r8) THEN
              pgrad(iobs)=pgrad(iobs)/ObsErr(iobs)
            END IF
          END DO
          pgrad(Ndatum(ng)+1)=0.0_r8
          DO iobs=1,Ndatum(ng)
            pgrad(iobs)=pgrad(iobs)+vcglwk(iobs,innLoop,outLoop)
          END DO
          pgrad(Ndatum(ng)+1)=pgrad(Ndatum(ng)+1)+                      &
     &                        vcglwk(Ndatum(ng)+1,innLoop,outLoop)
          IF (innLoop.gt.1) THEN
            DO iobs=1,Ndatum(ng)+1
              pgrad(iobs)=pgrad(iobs)-                                  &
     &                    cg_beta(innLoop,outLoop)*                     &
     &                    zcglwk(iobs,innLoop-1,outLoop)
            END DO
          END IF
!
!>        tl_cg_delta(innLoop)=tl_cg_delta(innLoop)+                    &
!>   &                         tl_Jb0(outLoop-1)*                       &
!>   &                         vcglwk(Ndatum(ng)+1,innLoop,outLoop)*    &
!>   &                         pgrad(Ndatum(ng)+1)+                     &
!>   &                         Jb0(outLoop-1)*                          &
!>   &                         tl_vcglwk(Ndatum(ng)+1,innLoop)*         &
!>   &                         pgrad(Ndatum(ng)+1)+                     &
!>   &                         Jb0(outLoop-1)*                          &
!>   &                         vcglwk(Ndatum(ng)+1,innLoop,outLoop)*    &
!>   &                         tl_pgrad(Ndatum(ng)+1)
!>
          ad_Jb0(outLoop-1)=ad_Jb0(outLoop-1)+                          &
     &                      vcglwk(Ndatum(ng)+1,innLoop,outLoop)*       &
     &                      pgrad(Ndatum(ng)+1)*                        &
     &                      ad_cg_delta(innLoop)
          ad_vcglwk(Ndatum(ng)+1,innLoop)=                              &
     &                           ad_vcglwk(Ndatum(ng)+1,innLoop)+       &
     &                           Jb0(outLoop-1)*                        &
     &                           pgrad(Ndatum(ng)+1)*                   &
     &                           ad_cg_delta(innLoop)
          ad_pgrad(Ndatum(ng)+1)=ad_pgrad(Ndatum(ng)+1)+                &
     &                           Jb0(outLoop-1)*                        &
     &                           vcglwk(Ndatum(ng)+1,innLoop,outLoop)*  &
     &                           ad_cg_delta(innLoop)
!
          IF (innLoop.eq.1) THEN
            fact=1.0_r8/cg_Gnorm_v(outLoop)
          ELSE
            fact=1.0_r8/cg_beta(innLoop,outLoop)
          END IF
!
          DO iobs=1,Ndatum(ng)
!>          tl_cg_delta(innLoop)=tl_cg_delta(innLoop)+                  &
!>   &                           tl_TLmodVal(iobs)*                     &
!>   &                           pgrad(iobs)+                           &
!>   &                           TLmodVal_S(iobs,innLoop,outLoop)*      &
!>   &                           tl_pgrad(iobs)*                        &
!>   &                           fact+                                  &
!>   &                           tl_Hbk(iobs)*                          &
!>   &                           vcglwk(Ndatum(ng)+1,innLoop,outLoop)*  &
!>   &                           pgrad(iobs)+                           &
!>   &                           Hbk(iobs,outLoop)*                     &
!>   &                           tl_vcglwk(Ndatum(ng)+1,innLoop)*       &
!>   &                           pgrad(iobs)+                           &
!>   &                           Hbk(iobs,outLoop)*                     &
!>   &                           vcglwk(Ndatum(ng)+1,innLoop,outLoop)*  &
!>   &                           tl_pgrad(iobs)+                        &
!>   &                           tl_Hbk(iobs)*                          &
!>   &                           vcglwk(iobs,innLoop,outLoop)*          &
!>   &                           pgrad(Ndatum(ng)+1)+                   &
!>   &                           Hbk(iobs,outLoop)*                     &
!>   &                           tl_vcglwk(iobs,innLoop)*               &
!>   &                           pgrad(Ndatum(ng)+1)+                   &
!>   &                           Hbk(iobs,outLoop)*                     &
!>   &                           vcglwk(iobs,innLoop,outLoop)*          &
!>   &                           tl_pgrad(Ndatum(ng)+1)
!>
             ADmodVal(iobs)=ADmodVal(iobs)+                             &
     &                      pgrad(iobs)*ad_cg_delta(innLoop)
             ad_pgrad(iobs)=ad_pgrad(iobs)+                             &
     &                      TLmodVal_S(iobs,innLoop,outLoop)*           &
     &                      ad_cg_delta(innLoop)*                       &
     &                      fact
             ad_Hbk(iobs)=ad_Hbk(iobs)+                                 &
     &                    vcglwk(Ndatum(ng)+1,innLoop,outLoop)*         &
     &                    pgrad(iobs)*                                  &
     &                    ad_cg_delta(innLoop)
             ad_vcglwk(Ndatum(ng)+1,innLoop)=                           &
     &                              ad_vcglwk(Ndatum(ng)+1,innLoop)+    &
     &                              Hbk(iobs,outLoop)*                  &
     &                              pgrad(iobs)*                        &
     &                              ad_cg_delta(innLoop)
             ad_pgrad(iobs)=ad_pgrad(iobs)+                             &
     &                      Hbk(iobs,outLoop)*                          &
     &                      vcglwk(Ndatum(ng)+1,innLoop,outLoop)*       &
     &                      ad_cg_delta(innLoop)
             ad_Hbk(iobs)=ad_Hbk(iobs)+                                 &
     &                    vcglwk(iobs,innLoop,outLoop)*                 &
     &                    pgrad(Ndatum(ng)+1)*                          &
     &                    ad_cg_delta(innLoop)
             ad_vcglwk(iobs,innLoop)=ad_vcglwk(iobs,innLoop)+           &
     &                               Hbk(iobs,outLoop)*                 &
     &                               pgrad(Ndatum(ng)+1)*               &
     &                               ad_cg_delta(innLoop)
             ad_pgrad(Ndatum(ng)+1)=ad_pgrad(Ndatum(ng)+1)+             &
     &                              Hbk(iobs,outLoop)*                  &
     &                              vcglwk(iobs,innLoop,outLoop)*       &
     &                              ad_cg_delta(innLoop)
          END DO
!
!>        tl_cg_delta(innLoop)=0.0_r8
!>
          ad_cg_delta(innLoop)=0.0_r8

          IF (innLoop.gt.1) THEN
!
!  Selime code: w=w-v0*beta
!
!  AUGMENTED.
!
            DO iobs=1,Ndatum(ng)+1
!>            tl_pgrad(iobs)=tl_pgrad(iobs)-                            &
!>   &                       tl_cg_beta(innLoop)*                       &
!>   &                       zcglwk(iobs,innLoop-1,outLoop)-            &
!>   &                       cg_beta(innLoop,outLoop)*                  &
!>   &                       tl_zcglwk(iobs,innLoop-1)
!>
              ad_cg_beta(innLoop)=ad_cg_beta(innLoop)-                  &
     &                            zcglwk(iobs,innLoop-1,outLoop)*       &
     &                            ad_pgrad(iobs)
              ad_zcglwk(iobs,innLoop-1)=ad_zcglwk(iobs,innLoop-1)-      &
     &                                  cg_beta(innLoop,outLoop)*       &
     &                                  ad_pgrad(iobs)
            END DO
          END IF
!
!  Selime code: w=t1/R+z1
!
!>        tl_pgrad(Ndatum(ng)+1)=tl_pgrad(Ndatum(ng)+1)+                &
!>   &                           tl_vcglwk(Ndatum(ng)+1,innLoop)
!>
          ad_vcglwk(Ndatum(ng)+1,innLoop)=                              &
     &                                ad_vcglwk(Ndatum(ng)+1,innLoop)+  &
     &                                ad_pgrad(Ndatum(ng)+1)
          DO iobs=1,Ndatum(ng)
!>          tl_pgrad(iobs)=tl_pgrad(iobs)+tl_vcglwk(iobs,innLoop)
!>
            ad_vcglwk(iobs,innLoop)=ad_vcglwk(iobs,innLoop)+            &
     &                              ad_pgrad(iobs)
          END DO
!
!>        tl_pgrad(Ndatum(ng)+1)=0.0_r8
!>
          ad_pgrad(Ndatum(ng)+1)=0.0_r8
!
!  Recompute pgrad.
!
          IF (innLoop.eq.1) THEN
            fact=1.0_r8/cg_Gnorm_v(outLoop)
          ELSE
            fact=1.0_r8/cg_beta(innLoop,outLoop)
          END IF
          DO iobs=1,Ndatum(ng)
            pgrad(iobs)=TLmodVal_S(iobs,innLoop,outLoop)*fact+          &
     &                  Hbk(iobs,outLoop)*                              &
     &                  vcglwk(Ndatum(ng)+1,innLoop,outLoop)
          END DO
!
          DO iobs=1,Ndatum(ng)
!
!  Convert gradient contribution from x-space to v-space.
!
            IF (ObsErr(iobs).NE.0.0_r8) THEN
              pgrad(iobs)=pgrad(iobs)/ObsErr(iobs)
!
!  Selime code: w=t1/R. No augmented term here since Raug^-1=[R^-1 0].
!
!>            tl_pgrad(iobs)=tl_pgrad(iobs)/ObsErr(iobs)-               &
!>   &                       tl_ObsErr(iobs)*pgrad(iobs)/ObsErr(iobs)
!>
              ad_ObsErr(iobs)=ad_ObsErr(iobs)-                          &
     &                         pgrad(iobs)*ad_pgrad(iobs)/ObsErr(iobs)
              ad_pgrad(iobs)=ad_pgrad(iobs)/ObsErr(iobs)
            END IF
          END DO
!
          IF (innLoop.eq.1) THEN
            fact=1.0_r8/cg_Gnorm_v(outLoop)
          ELSE
            fact=1.0_r8/cg_beta(innLoop,outLoop)
          END IF
!
          DO iobs=1,Ndatum(ng)
!
!  Selime code: w=t1.
!
!>          tl_pgrad(iobs)=tl_TLmodVal(iobs)+                           &
!>   &                     tl_Hbk(iobs)*                                &
!>   &                     vcglwk(Ndatum(ng)+1,innLoop,outLoop)+        &
!>   &                     Hbk(iobs,outLoop)*                           &
!>   &                     tl_vcglwk(Ndatum(ng)+1,innLoop)
!>
            ADmodVal(iobs)=ADmodVal(iobs)+ad_pgrad(iobs)
            ad_Hbk(iobs)=ad_Hbk(iobs)+                                  &
     &                   vcglwk(Ndatum(ng)+1,innLoop,outLoop)*          &
     &                   ad_pgrad(iobs)
            ad_vcglwk(Ndatum(ng)+1,innLoop)=                            &
     &                                ad_vcglwk(Ndatum(ng)+1,innLoop)+  &
     &                                Hbk(iobs,outLoop)*                &
     &                                ad_pgrad(iobs)
            ad_pgrad(iobs)=0.0_r8
          END DO

!
!  ADJOINT OF COMPLETE THE CALCULATION OF PREVIOUS LANCZOS VECTOR.
!
!  Complete the computation of the Lanczos vector from previous
! inner-loop
!
          IF (innLoop.EQ.1) THEN
!
            DO iobs=1,Ndatum(ng)
!
!  Selime code: t1=t1/beta_first=TLmodVal
!
!>            tl_TLmodVal(iobs)=tl_TLmodVal(iobs)/cg_Gnorm_v(outLoop)-  &
!>   &                          tl_cg_Gnorm_v(outLoop)*                 &
!>   &                          TLmodVal_S(iobs,innLoop,outLoop)/       &
!>   &                          (cg_Gnorm_v(outLoop)*                   &
!>   &                           cg_Gnorm_v(outLoop))
!>
              ad_cg_Gnorm_v(outLoop)=ad_cg_Gnorm_v(outLoop)-            &
     &                               ADmodVal(iobs)*                    &
     &                               TLmodVal_S(iobs,innLoop,outLoop)/  &
     &                               (cg_Gnorm_v(outLoop)*              &
     &                                cg_Gnorm_v(outLoop))
              ADmodVal(iobs)=ADmodVal(iobs)/cg_Gnorm_v(outLoop)
!>            tl_TLmodVal_S(iobs,innLoop,outLoop)=tl_TLmodVal(iobs)
!>
              ADmodVal(iobs)=ADmodVal(iobs)+                            &
     &                       ad_TLmodVal_S(iobs,innLoop,outLoop)
              ad_TLmodVal_S(iobs,innLoop,outLoop)=0.0_r8
            END DO
!
!  Normalize TLmodVal here and save in TLmodVal_S.
!
!  AUGMENTED.
!
            DO iobs=1,Ndatum(ng)+1
!
!  Selime code: z1=z1/beta_first=vcglwk(1).
!
!>            tl_vcglwk(iobs,1)=tl_zgrad0(iobs)/cg_Gnorm_v(outLoop)-    &
!>   &                          tl_cg_Gnorm_v(outLoop)*                 &
!>   &                          vcglwk(iobs,1,outLoop)/                 &
!>   &                          cg_Gnorm_v(outLoop)
!>
              ad_zgrad0(iobs)=ad_zgrad0(iobs)+                          &
     &                        ad_vcglwk(iobs,1)/cg_Gnorm_v(outLoop)
              ad_cg_Gnorm_v(outLoop)=ad_cg_Gnorm_v(outLoop)-            &
     &                               vcglwk(iobs,1,outLoop)*            &
     &                               ad_vcglwk(iobs,1)/                 &
     &                               cg_Gnorm_v(outLoop)
              ad_vcglwk(iobs,1)=0.0_r8
!
!  Selime code: v1=r/beta_first=zcglwk(1).
!
!>            tl_zcglwk(iobs,1)=tl_vgrad0(iobs)/                        &
!>   &                          cg_Gnorm_v(outLoop)-                    &
!>   &                          tl_cg_Gnorm_v(outLoop)*                 &
!>   &                          zcglwk(iobs,1,outLoop)/                 &
!>   &                          cg_Gnorm_v(outLoop)
!>
              ad_vgrad0(iobs)=ad_vgrad0(iobs)+ad_zcglwk(iobs,1)/        &
     &                        cg_Gnorm_v(outLoop)
              ad_cg_Gnorm_v(outLoop)=ad_cg_Gnorm_v(outLoop)-            &
     &                               zcglwk(iobs,1,outLoop)*            &
     &                               ad_zcglwk(iobs,1)/                 &
     &                               cg_Gnorm_v(outLoop)
              ad_zcglwk(iobs,1)=0.0_r8
            END DO
!
!  Selime code: beta_first=sqrt(r't1)=cg_Gnorm_v.
!                  cg_QG(1)=beta_first
!
!>          tl_cg_QG(1)=tl_cg_Gnorm_v(outLoop)
!>
            ad_cg_Gnorm_v(outLoop)=ad_cg_Gnorm_v(outLoop)+ad_cg_QG(1)
            ad_cg_QG(1)=0.0_r8
!>          tl_cg_Gnorm_v(outLoop)=0.5_r8*tl_cg_Gnorm_v(outLoop)/       &
!>   &                             cg_Gnorm_v(outLoop)
!>
            ad_cg_Gnorm_v(outLoop)=0.5_r8*ad_cg_Gnorm_v(outLoop)/       &
     &                           cg_Gnorm_v(outLoop)
!
!>          tl_cg_Gnorm_v(outLoop)=tl_cg_Gnorm_v(outLoop)+              &
!>   &                             tl_Jb0(outLoop-1)*                   &
!>   &                             vgrad0(Ndatum(ng)+1)*                &
!>   &                             zgrad0(Ndatum(ng)+1,outLoop)+        &
!>   &                             Jb0(outLoop-1)*                      &
!>   &                             tl_vgrad0(Ndatum(ng)+1)*             &
!>   &                             zgrad0(Ndatum(ng)+1,outLoop)+        &
!>   &                             Jb0(outLoop-1)*                      &
!>   &                             vgrad0(Ndatum(ng)+1)*                &
!>   &                             tl_zgrad0(Ndatum(ng)+1)
!>
            ad_Jb0(outLoop-1)=ad_Jb0(outLoop-1)+                        &
     &                        vgrad0(Ndatum(ng)+1)*                     &
     &                        zgrad0(Ndatum(ng)+1,outLoop)*             &
     &                        ad_cg_Gnorm_v(outLoop)
            ad_vgrad0(Ndatum(ng)+1)=ad_vgrad0(Ndatum(ng)+1)+            &
     &                              Jb0(outLoop-1)*                     &
     &                              zgrad0(Ndatum(ng)+1,outLoop)*       &
     &                              ad_cg_Gnorm_v(outLoop)
            ad_zgrad0(Ndatum(ng)+1)=ad_zgrad0(Ndatum(ng)+1)+            &
     &                              Jb0(outLoop-1)*                     &
     &                              vgrad0(Ndatum(ng)+1)*               &
     &                              ad_cg_Gnorm_v(outLoop)
            DO iobs=1,Ndatum(ng)
!>              tl_cg_Gnorm_v(outLoop)=tl_cg_Gnorm_v(outLoop)+          &
!>   &                                 tl_Hbk(iobs)*                    &
!>   &                                 zgrad0(Ndatum(ng)+1,outLoop)*    &
!>   &                                 vgrad0(iobs)+                    &
!>   &                                 Hbk(iobs,outLoop)*               &
!>   &                                 tl_zgrad0(Ndatum(ng)+1)*         &
!>   &                                 vgrad0(iobs)+                    &
!>   &                                 Hbk(iobs,outLoop)*               &
!>   &                                 zgrad0(Ndatum(ng)+1,outLoop)*    &
!>   &                                 tl_vgrad0(iobs)+                 &
!>   &                                 tl_Hbk(iobs)*                    &
!>   &                                 vgrad0(Ndatum(ng)+1)*            &
!>   &                                 zgrad0(iobs,outLoop)+            &
!>   &                                 Hbk(iobs,outLoop)*               &
!>   &                                 tl_vgrad0(Ndatum(ng)+1)*         &
!>   &                                 zgrad0(iobs,outLoop)+            &
!>   &                                 Hbk(iobs,outLoop)*               &
!>   &                                 vgrad0(Ndatum(ng)+1)*            &
!>   &                                 tl_zgrad0(iobs)
!>
                ad_Hbk(iobs)=ad_Hbk(iobs)+                              &
     &                       zgrad0(Ndatum(ng)+1,outLoop)*              &
     &                       vgrad0(iobs)*                              &
     &                       ad_cg_Gnorm_v(outLoop)
                ad_zgrad0(Ndatum(ng)+1)=ad_zgrad0(Ndatum(ng)+1)+        &
     &                                  Hbk(iobs,outLoop)*              &
     &                                  vgrad0(iobs)*                   &
     &                                  ad_cg_Gnorm_v(outLoop)
                ad_vgrad0(iobs)=ad_vgrad0(iobs)+                        &
     &                          Hbk(iobs,outLoop)*                      &
     &                          zgrad0(Ndatum(ng)+1,outLoop)*           &
     &                          ad_cg_Gnorm_v(outLoop)
                ad_Hbk(iobs)=ad_Hbk(iobs)+                              &
     &                       vgrad0(Ndatum(ng)+1)*                      &
     &                       zgrad0(iobs,outLoop)*                      &
     &                       ad_cg_Gnorm_v(outLoop)
                ad_vgrad0(Ndatum(ng)+1)=ad_vgrad0(Ndatum(ng)+1)+        &
     &                                  Hbk(iobs,outLoop)*              &
     &                                  zgrad0(iobs,outLoop)*           &
     &                                  ad_cg_Gnorm_v(outLoop)
                ad_zgrad0(iobs)=ad_zgrad0(iobs)+                        &
     &                          Hbk(iobs,outLoop)*                      &
     &                          vgrad0(Ndatum(ng)+1)*                   &
     &                          ad_cg_Gnorm_v(outLoop)
            END DO
!
            DO iobs=1,Ndatum(ng)
!>            tl_cg_Gnorm_v(outLoop)=tl_cg_Gnorm_v(outLoop)+            &
!>   &                               tl_TLmodVal(iobs)*                 &
!>   &                               vgrad0(iobs)+                      &
!>   &                               TLmodVal_S(iobs,innLoop,outLoop)*  &
!>   &                               tl_vgrad0(iobs)
!>
              ADmodVal(iobs)=ADmodVal(iobs)+                            &
     &                       vgrad0(iobs)*                              &
     &                       ad_cg_Gnorm_v(outLoop)
              ad_vgrad0(iobs)=ad_vgrad0(iobs)+                          &
     &                        TLmodVal_S(iobs,innLoop,outLoop)*         &
     &                        ad_cg_Gnorm_v(outLoop)
            END DO
!>          tl_cg_Gnorm_v(outLoop)=0.0_r8
!>
            ad_cg_Gnorm_v(outLoop)=0.0_r8

          ELSE
!
            IF (innLoop.eq.NinnLoop) THEN
!
!  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=1                 ! Spectral LMP
!>              Laug=.TRUE.
!>              CALL rprecond(ng, Lscale, Laug, outLoop, NinnLoop-1, px)
!>            END IF
!
!  Recompute zu and zwork.
!
!
            zbet=cg_delta(1,outLoop)
            zu(1)=cg_QG(1,outLoop)/zbet
            DO ivec=2,innLoop-1
              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
            zwork(innLoop-1,3)=zu(innLoop-1)
!
            DO ivec=innLoop-2,1,-1
              zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
              zwork(ivec,3)=zu(ivec)
            END DO
!
            DO iobs=1,Ndatum(ng)+1
              DO ivec=1,innLoop-1
!>              tl_px(iobs)=tl_px(iobs)+                                &
!>   &                      tl_zcglwk(iobs,ivec)*                       &
!>   &                      zwork(ivec,3)+                              &
!>   &                      zcglwk(iobs,ivec,outLoop)*                  &
!>   &                      tl_zwork(ivec,3)
!>
                ad_zcglwk(iobs,ivec)=ad_zcglwk(iobs,ivec)+              &
     &                               zwork(ivec,3)*                     &
     &                               ad_px(iobs)
                ad_zwork(ivec,3)=ad_zwork(ivec,3)+                      &
     &                           zcglwk(iobs,ivec,outLoop)*             &
     &                           ad_px(iobs)
              END DO
!>            tl_px(iobs)=0.0_r8
!>
              ad_px(iobs)=0.0_r8
            END DO
!
            IF (NinnLoop.eq.1) THEN
              print *,'Illegal configuration!'
              print *,'Ninner must be ge 2'
              STOP
            END IF
            IF (NinnLoop.eq.2) THEN
              zu(1)=cg_QG(1,outLoop)/cg_delta(1,outLoop)
!>            tl_zwork(1,3)=tl_zu(1)
!>
              ad_zu(1)=ad_zu(1)+ad_zwork(1,3)
              ad_zwork(1,3)=0.0_r8
!>            tl_zu(1)=tl_zrhs(1)/cg_delta(1,outLoop)
!>
              ad_zrhs(1)=ad_zrhs(1)+ad_zu(1)/cg_delta(1,outLoop)
              ad_zu(1)=0.0_r8
!>            tl_zrhs(1)=tl_cg_QG(1)-tl_cg_delta(1)*zu(1)
!>
              ad_cg_QG(1)=ad_cg_QG(1)+ad_zrhs(1)
              ad_cg_delta(1)=ad_cg_delta(1)-zu(1)*ad_zrhs(1)
              ad_zrhs(1)=0.0_r8
            ELSE
!>            DO ivec=innLoop-2,1,-1
!>
              DO ivec=1,innLoop-2
!>              tl_zwork(ivec,3)=tl_zu(ivec)
!>
                ad_zu(ivec)=ad_zu(ivec)+ad_zwork(ivec,3)
                ad_zwork(ivec,3)=0.0_r8

!>              tl_zu(ivec)=tl_zu(ivec)-zgam(ivec+1)*tl_zu(ivec+1)
!>
                ad_zu(ivec+1)=ad_zu(ivec+1)-zgam(ivec+1)*ad_zu(ivec)
              END DO
!>            tl_zwork(innLoop-1,3)=tl_zu(innLoop-1)
!>
              ad_zu(innLoop-1)=ad_zu(innLoop-1)+ad_zwork(innLoop-1,3)
              ad_zwork(innLoop-1,3)=0.0_r8

!>            DO ivec=2,innLoop-1
!>
              DO ivec=innLoop-1,2,-1
                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
!>
                adfac=ad_zu(ivec)/zbet
                ad_zu(ivec-1)=ad_zu(ivec-1)-                            &
     &                        cg_beta(ivec,outLoop)*adfac
                ad_zrhs(ivec)=ad_zrhs(ivec)+adfac
                ad_zu(ivec)=0.0_r8
              END DO
              zbet=cg_delta(1,outLoop)
!>            tl_zu(1)=tl_zrhs(1)/zbet
!>
              ad_zrhs(1)=ad_zrhs(1)+ad_zu(1)/zbet
              ad_zu(1)=0.0_r8
!
!>            tl_zrhs(innLoop-1)=tl_cg_QG(innLoop-1)-                   &
!>   &                         tl_cg_beta(innLoop-1)*zu(innLoop-2)-     &
!>   &                         tl_cg_delta(innLoop-1)*zu(innLoop-1)
!>
              ad_cg_QG(innLoop-1)=ad_cg_QG(innLoop-1)+                  &
     &                            ad_zrhs(innLoop-1)
              ad_cg_beta(innLoop-1)=ad_cg_beta(innLoop-1)-              &
     &                              zu(innLoop-2)*                      &
     &                              ad_zrhs(innLoop-1)
              ad_cg_delta(innLoop-1)=ad_cg_delta(innLoop-1)-            &
     &                               zu(innLoop-1)*                     &
     &                               ad_zrhs(innLoop-1)
              ad_zrhs(innLoop-1)=0.0_r8
              DO ivec=2,innLoop-2
!>              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)
!>
                ad_cg_QG(ivec)=ad_cg_QG(ivec)+                          &
     &                         ad_zrhs(ivec)
                ad_cg_beta(ivec)=ad_cg_beta(ivec)-                      &
     &                           zu(ivec-1)*ad_zrhs(ivec)
                ad_cg_delta(ivec)=ad_cg_delta(ivec)-                    &
     &                            zu(ivec)*ad_zrhs(ivec)
                ad_cg_beta(ivec+1)=ad_cg_beta(ivec+1)-                  &
     &                             zu(ivec+1)*ad_zrhs(ivec)
                ad_zrhs(ivec)=0.0_r8
              END DO
!>            tl_zrhs(1)=tl_cg_QG(1)-                                   &
!>   &                   tl_cg_delta(1)*zu(1)-                          &
!>   &                   tl_cg_beta(2)*zu(2)
!>
              ad_cg_QG(1)=ad_cg_QG(1)+ad_zrhs(1)
              ad_cg_delta(1)=ad_cg_delta(1)-zu(1)*ad_zrhs(1)
              ad_cg_beta(2)=ad_cg_beta(2)-zu(2)*ad_zrhs(1)
              ad_zrhs(1)=0.0_r8
            END IF

           END IF

!>          tl_cg_QG(innLoop)=tl_cg_QG(innLoop)+                        &
!>   &                        tl_Jb0(outLoop-1)*                        &
!>   &                        vcglwk(Ndatum(ng)+1,innLoop,outLoop)*     &
!>   &                        vgrad0(Ndatum(ng)+1)+                     &
!>   &                        Jb0(outLoop-1)*                           &
!>   &                        tl_vcglwk(Ndatum(ng)+1,innLoop)*          &
!>   &                        vgrad0(Ndatum(ng)+1)+                     &
!>   &                        Jb0(outLoop-1)*                           &
!>   &                        vcglwk(Ndatum(ng)+1,innLoop,outLoop)*     &
!>   &                        tl_vgrad0(Ndatum(ng)+1)
!>
            ad_Jb0(outLoop-1)=ad_Jb0(outLoop-1)+                        &
     &                        vcglwk(Ndatum(ng)+1,innLoop,outLoop)*     &
     &                        vgrad0(Ndatum(ng)+1)*                     &
     &                        ad_cg_QG(innLoop)
            ad_vcglwk(Ndatum(ng)+1,innLoop)=                            &
     &                           ad_vcglwk(Ndatum(ng)+1,innLoop)+       &
     &                           Jb0(outLoop-1)*                        &
     &                           vgrad0(Ndatum(ng)+1)*ad_cg_QG(innLoop)
            ad_vgrad0(Ndatum(ng)+1)=ad_vgrad0(Ndatum(ng)+1)+            &
     &                              Jb0(outLoop-1)*                     &
     &                            vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
     &                              ad_cg_QG(innLoop)
!
            DO iobs=1,Ndatum(ng)
!>            tl_cg_QG(innLoop)=tl_cg_QG(innLoop)+                      &
!>   &                          tl_Hbk(iobs)*                           &
!>   &                          vcglwk(Ndatum(ng)+1,innLoop,outLoop)*   &
!>   &                          vgrad0(iobs)+                           &
!>   &                          Hbk(iobs,outLoop)*                      &
!>   &                          tl_vcglwk(Ndatum(ng)+1,innLoop)*        &
!>   &                          vgrad0(iobs)+                           &
!>   &                          Hbk(iobs,outLoop)*                      &
!>   &                          vcglwk(Ndatum(ng)+1,innLoop,outLoop)*   &
!>   &                          tl_vgrad0(iobs)+                        &
!>   &                          tl_Hbk(iobs)*                           &
!>   &                          vcglwk(iobs,innLoop,outLoop)*           &
!>   &                          vgrad0(Ndatum(ng)+1)+                   &
!>   &                          Hbk(iobs,outLoop)*                      &
!>   &                          tl_vcglwk(iobs,innLoop)*                &
!>   &                          vgrad0(Ndatum(ng)+1)+                   &
!>   &                          Hbk(iobs,outLoop)*                      &
!>   &                          vcglwk(iobs,innLoop,outLoop)*           &
!>   &                          tl_vgrad0(Ndatum(ng)+1)
!>
              ad_Hbk(iobs)=ad_Hbk(iobs)+                                &
     &                     vcglwk(Ndatum(ng)+1,innLoop,outLoop)*        &
     &                     vgrad0(iobs)*                                &
     &                     ad_cg_QG(innLoop)
              ad_vcglwk(Ndatum(ng)+1,innLoop)=                          &
     &                               ad_vcglwk(Ndatum(ng)+1,innLoop)+   &
     &                               Hbk(iobs,outLoop)*                 &
     &                               vgrad0(iobs)*                      &
     &                               ad_cg_QG(innLoop)
              ad_vgrad0(iobs)=ad_vgrad0(iobs)+                          &
     &                        Hbk(iobs,outLoop)*                        &
     &                        vcglwk(Ndatum(ng)+1,innLoop,outLoop)*     &
     &                        ad_cg_QG(innLoop)
              ad_Hbk(iobs)=ad_Hbk(iobs)+                                &
     &                     vcglwk(iobs,innLoop,outLoop)*                &
     &                     vgrad0(Ndatum(ng)+1)*                        &
     &                     ad_cg_QG(innLoop)
              ad_vcglwk(iobs,innLoop)=ad_vcglwk(iobs,innLoop)+          &
     &                                Hbk(iobs,outLoop)*                &
     &                                vgrad0(Ndatum(ng)+1)*             &
     &                                ad_cg_QG(innLoop)
              ad_vgrad0(Ndatum(ng)+1)=ad_vgrad0(Ndatum(ng)+1)+          &
     &                                Hbk(iobs,outLoop)*                &
     &                                vcglwk(iobs,innLoop,outLoop)*     &
     &                                ad_cg_QG(innLoop)
            END DO
!
!  Selime code:  cg_QG=zgrad0*HBH'*zcglwk
!
            DO iobs=1,Ndatum(ng)
!>            tl_cg_QG(innLoop)=tl_cg_QG(innLoop)+                      &
!>   &                          tl_TLmodVal(iobs)*                      &
!>   &                          vgrad0(iobs)+                           &
!>   &                          TLmodVal_S(iobs,innLoop,outLoop)*       &
!>   &                          tl_vgrad0(iobs)/                        &
!>   &                          cg_beta(innLoop,outLoop)
              ADmodVal(iobs)=ADmodVal(iobs)+                            &
     &                       vgrad0(iobs)*ad_cg_QG(innLoop)
              ad_vgrad0(iobs)=ad_vgrad0(iobs)+                          &
     &                        TLmodVal_S(iobs,innLoop,outLoop)*         &
     &                        ad_cg_QG(innLoop)/                        &
     &                        cg_beta(innLoop,outLoop)
            END DO
!>          tl_cg_QG(innLoop)=0.0_r8
!>
            ad_cg_QG(innLoop)=0.0_r8
!
!  Selime code: t1=t1/beta=TLmodVal
!
!  AUGMENTED
!
            DO iobs=1,Ndatum(ng)
!>            tl_TLmodVal(iobs)=tl_TLmodVal(iobs)/                      &
!>   &                          cg_beta(innLoop,outLoop)-               &
!>   &                          tl_cg_beta(innLoop)*                    &
!>   &                          TLmodVal_S(iobs,innLoop,outLoop)/       &
!>   &                          (cg_beta(innLoop,outLoop)*              &
!>   &                           cg_beta(innLoop,outLoop))
!>
              ad_cg_beta(innLoop)=ad_cg_beta(innLoop)-                  &
     &                            TLmodVal_S(iobs,innLoop,outLoop)*     &
     &                            ADmodVal(iobs)/                       &
     &                            (cg_beta(innLoop,outLoop)*            &
     &                             cg_beta(innLoop,outLoop))
              ADmodVal(iobs)=ADmodVal(iobs)/cg_beta(innLoop,outLoop)
!>            tl_TLmodVal_S(iobs,innLoop,outLoop)=tl_TLmodVal(iobs)
!>
              ADmodVal(iobs)=ADmodVal(iobs)+                            &
     &                       ad_TLmodVal_S(iobs,innLoop,outLoop)
              ad_TLmodVal_S(iobs,innLoop,outLoop)=0.0_r8
            END DO
!
!  Selime code: z1=w/beta=vcglwk
!
            DO iobs=1,Ndatum(ng)+1
!>            tl_vcglwk(iobs,innLoop)=tl_vcglwk(iobs,innLoop)/          &
!>   &                                cg_beta(innLoop,outLoop)-         &
!>   &                                tl_cg_beta(innLoop)*              &
!>   &                                vcglwk(iobs,innLoop,outLoop)/     &
!>   &                                cg_beta(innLoop,outLoop)
              ad_cg_beta(innLoop)=ad_cg_beta(innLoop)-                  &
     &                            vcglwk(iobs,innLoop,outLoop)*         &
     &                            ad_vcglwk(iobs,innLoop)/              &
     &                            cg_beta(innLoop,outLoop)
              ad_vcglwk(iobs,innLoop)=ad_vcglwk(iobs,innLoop)/          &
     &                                cg_beta(innLoop,outLoop)
!>            tl_zcglwk(iobs,innLoop)=tl_pgrad(iobs)/                   &
!>   &                                cg_beta(innLoop,outLoop)-         &
!>   &                                tl_cg_beta(innLoop)*              &
!>   &                                zcglwk(iobs,innLoop,outLoop)/     &
!>   &                                cg_beta(innLoop,outLoop)
              ad_cg_beta(innLoop)=ad_cg_beta(innLoop)-                  &
     &                            zcglwk(iobs,innLoop,outLoop)*         &
     &                            ad_zcglwk(iobs,innLoop)/              &
     &                            cg_beta(innLoop,outLoop)
              ad_pgrad(iobs)=ad_pgrad(iobs)+                            &
     &                       ad_zcglwk(iobs,innLoop)/                   &
     &                       cg_beta(innLoop,outLoop)
              ad_zcglwk(iobs,innLoop)=0.0_r8
            END DO
!
!  Recompute pgrad.
!
            DO iobs=1,Ndatum(ng)+1
              pgrad(iobs)=zcglwk(iobs,innLoop,outLoop)*                 &
     &                    cg_beta(innLoop,outLoop)
            END DO
!
!>          tl_cg_beta(innLoop)=0.5_r8*tl_cg_beta(innLoop)/             &
!>   &                          cg_beta(innLoop,outLoop)
!>
            ad_cg_beta(innLoop)=0.5_r8*ad_cg_beta(innLoop)/             &
     &                         cg_beta(innLoop,outLoop)
!>          tl_cg_beta(innLoop)=tl_cg_beta(innLoop)+                    &
!>   &                          tl_Jb0(outLoop-1)*                      &
!>   &                          pgrad(Ndatum(ng)+1)*                    &
!>   &                          vcglwk(Ndatum(ng)+1,innLoop,outLoop)+   &
!>   &                          Jb0(outLoop-1)*                         &
!>   &                          tl_pgrad(Ndatum(ng)+1)*                 &
!>   &                          vcglwk(Ndatum(ng)+1,innLoop,outLoop)+   &
!>   &                          Jb0(outLoop-1)*                         &
!>   &                          pgrad(Ndatum(ng)+1)*                    &
!>   &                          tl_vcglwk(Ndatum(ng)+1,innLoop)
!>
            ad_Jb0(outLoop-1)=ad_Jb0(outLoop-1)+                        &
     &                        pgrad(Ndatum(ng)+1)*                      &
     &                        vcglwk(Ndatum(ng)+1,innLoop,outLoop)*     &
     &                        ad_cg_beta(innLoop)
            ad_pgrad(Ndatum(ng)+1)=ad_pgrad(Ndatum(ng)+1)+              &
     &                             Jb0(outLoop-1)*                      &
     &                             vcglwk(Ndatum(ng)+1,innLoop,outLoop)*&
     &                             ad_cg_beta(innLoop)
            ad_vcglwk(Ndatum(ng)+1,innLoop)=                            &
     &                            ad_vcglwk(Ndatum(ng)+1,innLoop)+      &
     &                            Jb0(outLoop-1)*                       &
     &                            pgrad(Ndatum(ng)+1)*                  &
     &                            ad_cg_beta(innLoop)
!
            DO iobs=1,Ndatum(ng)
!>            tl_cg_beta(innLoop)=tl_cg_beta(innLoop)+                  &
!>   &                            tl_Hbk(iobs)*                         &
!>   &                            vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
!>   &                            pgrad(iobs)+                          &
!>   &                            Hbk(iobs,outLoop)*                    &
!>   &                            tl_vcglwk(Ndatum(ng)+1,innLoop)*      &
!>   &                            pgrad(iobs)+                          &
!>   &                            Hbk(iobs,outLoop)*                    &
!>   &                            vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
!>   &                            tl_pgrad(iobs)+                       &
!>   &                            tl_Hbk(iobs)*                         &
!>   &                            pgrad(Ndatum(ng)+1)*                  &
!>   &                            vcglwk(iobs,innLoop,outLoop)+         &
!>   &                            Hbk(iobs,outLoop)*                    &
!>   &                            tl_pgrad(Ndatum(ng)+1)*               &
!>   &                            vcglwk(iobs,innLoop,outLoop)+         &
!>   &                            Hbk(iobs,outLoop)*                    &
!>   &                            pgrad(Ndatum(ng)+1)*                  &
!>   &                            tl_vcglwk(iobs,innLoop)
!>
              ad_Hbk(iobs)=ad_Hbk(iobs)+                                &
     &                     vcglwk(Ndatum(ng)+1,innLoop,outLoop)*        &
     &                     pgrad(iobs)*                                 &
     &                     ad_cg_beta(innLoop)
              ad_vcglwk(Ndatum(ng)+1,innLoop)=                          &
     &                              ad_vcglwk(Ndatum(ng)+1,innLoop)+    &
     &                              Hbk(iobs,outLoop)*                  &
     &                              pgrad(iobs)*                        &
     &                              ad_cg_beta(innLoop)
              ad_pgrad(iobs)=ad_pgrad(iobs)+                            &
     &                       Hbk(iobs,outLoop)*                         &
     &                       vcglwk(Ndatum(ng)+1,innLoop,outLoop)*      &
     &                       ad_cg_beta(innLoop)
              ad_Hbk(iobs)=ad_Hbk(iobs)+                                &
     &                      pgrad(Ndatum(ng)+1)*                        &
     &                      vcglwk(iobs,innLoop,outLoop)*               &
     &                      ad_cg_beta(innLoop)
              ad_pgrad(Ndatum(ng)+1)=ad_pgrad(Ndatum(ng)+1)+            &
     &                               Hbk(iobs,outLoop)*                 &
     &                               vcglwk(iobs,innLoop,outLoop)*      &
     &                               ad_cg_beta(innLoop)
              ad_vcglwk(iobs,innLoop)=ad_vcglwk(iobs,innLoop)+          &
     &                                Hbk(iobs,outLoop)*                &
     &                                pgrad(Ndatum(ng)+1)*              &
     &                                ad_cg_beta(innLoop)
            END DO
!
!  Selime code:  beta=pgrad*HBH'*pgrad=cg_beta
!
            DO iobs=1,Ndatum(ng)
!>            tl_cg_beta(innLoop)=tl_cg_beta(innLoop)+                  &
!>   &                            tl_pgrad(iobs)*                       &
!>   &                            TLmodVal_S(iobs,innLoop,outLoop)+     &
!>   &                            pgrad(iobs)*                          &
!>   &                            tl_TLmodVal(iobs)
              ad_pgrad(iobs)=ad_pgrad(iobs)+                            &
     &                       TLmodVal_S(iobs,innLoop,outLoop)*          &
     &                       ad_cg_beta(innLoop)
              ADmodVal(iobs)=ADmodVal(iobs)+                            &
     &                       pgrad(iobs)*ad_cg_beta(innLoop)
            END DO
!>          tl_cg_beta(innLoop)=0.0_r8
!>
            ad_cg_beta(innLoop)=0.0_r8
!
            DO iobs=1,Ndatum(ng)+1
!
!  Selime code: w=pgrad, t1=TLmodVal.
!
!>            tl_pgrad(iobs)=tl_zcglwk(iobs,innLoop)
!>
              ad_zcglwk(iobs,innLoop)=ad_zcglwk(iobs,innLoop)+          &
     &                                ad_pgrad(iobs)
              ad_pgrad(iobs)=0.0_r8
            END DO

          END IF

        END IF MINIMIZE
!
        IF (innLoop.gt.0) THEN
          DO iobs=1,Ndatum(ng)
!>          TLmodVal(iobs)=ObsScale(iobs)*                              &
!>   &                     TLmodVal(iobs)+                              &
!>   &                     tl_ObsScale(iobs)*                           &
!>   &                     TLmodVal_S(iobs,innLoop,outLoop)
!>
            ad_ObsScale(iobs)=ad_ObsScale(iobs)+                        &
     &                        TLmodVal_S(iobs,innLoop,outLoop)*         &
     &                        ADmodVal(iobs)
            ADmodVal(iobs)=ObsScale(iobs)*ADmodVal(iobs)
          END DO
        END IF

      END IF MASTER_THREAD

# ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Broadcast new solution to other nodes.
!-----------------------------------------------------------------------
!
      CALL mp_bcastf (ng, model, ADmodVal)
      CALL mp_bcastf (ng, model, vgrad0(:))
      CALL mp_bcastf (ng, model, TLmodVal)
#  ifdef IMPACT_INNER
      CALL mp_bcastf (ng, model, ad_ObsVal(:,:))
#  else
      CALL mp_bcastf (ng, model, ad_ObsVal)
#  endif
      CALL mp_bcastf (ng, model, ad_ObsErr)
      CALL mp_bcastf (ng, model, ad_ObsScale)
#  if defined SENSITIVITY_4DVAR || \
      defined W4DPSAS           || defined W4DVAR
      CALL mp_bcastf (ng, model, ad_Hbk(:))
      CALL mp_bcastf (ng, model, ad_Jb0(outLoop-1))
#  endif
      CALL mp_bcasti (ng, model, info)
      CALL mp_bcastf (ng, model, ad_cg_beta(:))
      CALL mp_bcastf (ng, model, ad_cg_QG(:))
      CALL mp_bcastf (ng, model, ad_cg_delta(:))
      CALL mp_bcastf (ng, model, ad_zgrad0(:))
      CALL mp_bcastf (ng, model, ad_zcglwk(:,:))
      CALL mp_bcastf (ng, model, ad_vcglwk(:,:))
      CALL mp_bcastf (ng, model, FOURDVAR(ng)%ad_cg_pxsave(:))
# endif
!
      RETURN
      END SUBROUTINE ad_rpcg_lanczos

#else
      SUBROUTINE ad_rpcg_lanczos
      RETURN
      END SUBROUTINE ad_rpcg_lanczos
#endif
