rpcg_lanczos.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE rpcg_lanczos_mod
 
#if defined WEAK_CONSTRAINT && defined RPCG
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group            Selime Gurol      !
!    Licensed under a MIT/X style license            Andrew M. Moore   !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  Weak Constraint 4-Dimensional Variational (4D-Var) Restricted       !
!        B-preconditioned Lanczos (RBLanczos) Algorithm                !
!                      (Gurol et al., 2014)                            !
!                                                                      !
!  4D-Var dual formulation (space spanned by observation vector)       !
!  minimization algorithm.  The RBLanczos is equivalent to the         !
!  Restricted B-preconditioned Conjugate Gradient (RBCG) used in       !
!  the primal formulation (Gurol et al., 2014).  It has similar        !
!  convergence properties to the primal formulation minimization       !
!  algorithms making weak-constraint 4D-Var very affordabl.            !
!                                                                      !
!  References:                                                         !
!                                                                      !
!    Gratton, S. and J. Tshimanga, 2009: An observation-space          !
!      formulation of variational assimilation using a restricted      !
!      preconditioned conjugate gradient algorithm, QJRMS, 135,        !
!      1573-1585.                                                      !
!                                                                      !
!    Gurol, S., A.T. Weaver, A.M. Moore, A. Piacentini, H.G. Arango,   !
!      S. Gratton, 2014: B-preconditioned minimization algorithms for  !
!      data assimilation with the dual formulation, QJRMS, 140,        !
!      539-556.                                                        !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_fourdvar
      Use mod_iounits
      USE mod_ncparam
      USE mod_scalars
!
# ifdef DISTRIBUTE
      USE distribute_mod, ONLY : mp_bcastf, mp_bcasti, mp_bcastl
# endif
      USE lapack_mod,     ONLY : dsteqr
      USE strings_mod,    ONLY : founderror
!
      implicit none
!
      PUBLIC  :: rpcg_lanczos
      PUBLIC  :: cg_read_rpcg
!
      PRIVATE :: cg_read_rpcg_nf90
# if defined PIO_LIB && defined DISTRIBUTE
      PRIVATE :: cg_read_rpcg_pio
# endif
      PRIVATE :: cg_write_rpcg
      PRIVATE :: cg_write_rpcg_nf90
# if defined PIO_LIB && defined DISTRIBUTE
      PRIVATE :: cg_write_rpcg_pio
# endif
      PRIVATE :: rpevecs
      PRIVATE :: rprecond
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE rpcg_lanczos (ng, model, outLoop, innLoop, NinnLoop,   &
     &                         Lcgini)
!***********************************************************************
!
!  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(dp) :: zbet, eps, preducv, preducy
      real(dp) :: jopt, jf, jmod, jdata, jb, jobs, jact, cff
!
      real(dp), dimension(NinnLoop) :: zu, zgam, zfact, ztemp3
      real(dp), dimension(NinnLoop) :: ztemp1, ztemp2, zu1, zu2
      real(dp), dimension(NinnLoop) :: ztemp4
      real(dp), dimension(Ndatum(ng)+1) :: px, pgrad, zw, zt
      real(dp), dimension(Ndatum(ng)+1) :: zhv, zht, zd
      real(dp), dimension(Ninner,3) :: zwork
      real(dp), dimension(2*(NinnLoop-1)) :: work
      real(dp), dimension(Ninner,Ninner) :: zgv
!
      character (len=13) :: string
 
      character (len=*), parameter :: myfile =                          &
     &  __FILE__
!
!=======================================================================
!  Weak constraint 4D-Var conjugate gradient, Lanczos-based algorithm.
!=======================================================================
 
# ifdef PROFILE
!
!  Turn on time clock.
!
      CALL wclock_on (ng, model, 85, __line__, myfile)
# endif
!
!  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 algorithm 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.
!
!  Initialize internal parameters.  This is needed here for output
!  reasons.
!
      jopt=0.0_dp
      jf=0.0_dp
      jdata=0.0_dp
      jmod=0.0_dp
      jb=jb0(outloop-1)
      jobs=0.0_dp
      jact=0.0_dp
      ritzmaxerr=hevecerr
      eps=0.0_dp
      preducv=0.0_dp
      preducy=0.0_dp
      zwork=0.0_dp
      zgv=0.0_dp
 
!
      master_thread : IF (master) THEN
!
!  Multiply TLmodVal by ObsScale to effectively remove any rejected
!  observations from the analysis.
!
        IF (innloop.gt.0) THEN
          DO iobs=1,ndatum(ng)
            tlmodval(iobs)=obsscale(iobs)*tlmodval(iobs)
          END DO
        END IF
!
!-----------------------------------------------------------------------
!  Initialization step, innloop = 0.
!-----------------------------------------------------------------------
!
        minimize : IF (innloop.eq.0) THEN
 
# if defined RBL4DVAR                 || \
     defined rbl4dvar_ana_sensitivity || \
     defined rbl4dvar_fct_sensitivity || \
     defined tl_rbl4dvar
!
          DO iobs=1,ndatum(ng)
            bckmodval(iobs)=nlmodval(iobs)
          END DO
# endif
!
!  For the first outer-loop, x0=0.
!
          IF (outloop.eq.1) THEN
            DO iobs=1,ndatum(ng)+1
              px(iobs)=0.0_dp
            END DO
            DO iobs=1,ndatum(ng)
              hbk(iobs,outloop)=0.0_dp
            END DO
          ELSE
!! NOTE: CODE WILL NOT WORK FOR R4D-Var AT THE MOMENT SINCE THIS IS
!! THE WRONG TLmodVal HERE!!!
            DO iobs=1,ndatum(ng)
              hbk(iobs,outloop)=-tlmodval(iobs)
            END DO
            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
          END IF
!
!  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)
!
!  Selime code:  r=b.
!
# if defined RBL4DVAR                 || \
     defined rbl4dvar_ana_sensitivity || \
     defined rbl4dvar_fct_sensitivity || \
     defined tl_rbl4dvar
              pgrad(iobs)=obsscale(iobs)*                               &
     &                    (obsval(iobs)-bckmodval(iobs))
# else
              pgrad(iobs)=obsscale(iobs)*                               &
     &                    (obsval(iobs)-tlmodval(iobs))
# endif
              IF (obserr(iobs).ne.0.0_r8) THEN
                 pgrad(iobs)=pgrad(iobs)/obserr(iobs)
              END IF
              vgrad0(iobs)=pgrad(iobs)
              vgrad0s(iobs)=vgrad0(iobs)
            END DO
!
!  AUGMENTED
!  Augment the residual vector.
!
            pgrad(ndatum(ng)+1)=1.0_dp
            vgrad0(ndatum(ng)+1)=pgrad(ndatum(ng)+1)
            vgrad0s(ndatum(ng)+1)=vgrad0(ndatum(ng)+1)
!
!  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
!
!  Selime code: zgrad0=z1.
!
!  Selime code: vgrad0=r.
!
            DO iobs=1,ndatum(ng)+1
              zgrad0(iobs,outloop)=pgrad(iobs)
            END DO
!
!  Selime code: ADmodVal=z1.
!
            DO iobs=1,ndatum(ng)
              admodval(iobs)=zgrad0(iobs,outloop)
            END DO
!
          preducv=1.0_dp
          preducy=1.0_dp
!
!-----------------------------------------------------------------------
!  Minimization step, innLoop > 0.
!-----------------------------------------------------------------------
!
        ELSE
!
!  COMPLETE THE CALCULATION OF PREVIOUS LANCZOS VECTOR.
!
!  Selime's code: t1=GDG'*z1=TLmodVal.
!
!  Complete the computation of the Lanczos vector from previous inner-loop
!
          IF (innloop.EQ.1) THEN
            cg_gnorm_v(outloop)=0.0_dp
            DO iobs=1,ndatum(ng)
              cg_gnorm_v(outloop)=cg_gnorm_v(outloop)+                  &
     &                            tlmodval(iobs)*vgrad0(iobs)
            END DO
!
!  AUGMENTED. Here Hbk=H(x0-xb). Hbk=0 and Jb0=0 when outer=1.
!  Jb0 is computed using sum_i (G'*w)_i, i.e. the sum of the outputs'
!  of the adjoint model at the end of AD_LOOP2.
!
            DO iobs=1,ndatum(ng)
              cg_gnorm_v(outloop)=cg_gnorm_v(outloop)+                  &
     &                            hbk(iobs,outloop)*                    &
     &                            zgrad0(ndatum(ng)+1,outloop)*         &
     &                            vgrad0(iobs)+                         &
     &                            hbk(iobs,outloop)*                    &
     &                            vgrad0(ndatum(ng)+1)*                 &
     &                            zgrad0(iobs,outloop)
            END DO
            cg_gnorm_v(outloop)=cg_gnorm_v(outloop)+                    &
     &                          jb0(outloop-1)*vgrad0(ndatum(ng)+1)*    &
     &                          zgrad0(ndatum(ng)+1,outloop)
!
            cg_gnorm_v(outloop)=sqrt(cg_gnorm_v(outloop))
!
!  Selime's code: beta_first=sqrt(r't1)=cg_Gnorm_v.
!                 cg_QG(1)=beta_first
!
            cg_qg(1,outloop)=cg_gnorm_v(outloop)
!
!  Normalize TLmodVal here and save in TLmodVal_S.
!
!  AUGMENTED.
!
            DO iobs=1,ndatum(ng)+1
!
!  Selime code: v1=r/beta_first=zcglwk(1).
!
              zcglwk(iobs,1,outloop)=vgrad0(iobs)/                      &
     &                               cg_gnorm_v(outloop)
!
!  Selime code: z1=z1/beta_first=vcglwk(1).
!
              vcglwk(iobs,1,outloop)=zgrad0(iobs,outloop)/              &
     &                               cg_gnorm_v(outloop)
            END DO
!
!  AUGMENTED
!
            DO iobs=1,ndatum(ng)
!
!  Selime code: t1=t1/beta_first=TLmodVal
!
              tlmodval_s(iobs,innloop,outloop)=tlmodval(iobs)
              tlmodval(iobs)=tlmodval(iobs)/cg_gnorm_v(outloop)
            END DO
!
            preducv=1.0_dp
            preducy=1.0_dp
!
!  Compute the actual penalty function terms.
!
            jobs=0.0_dp
!
!  NO AUGMENTED TERM HERE.
!
            jact=jb
            DO iobs=1,ndatum(ng)
              IF (obserr(iobs).ne.0.0_r8) THEN
                jact=jact+obserr(iobs)*vgrad0s(iobs)*vgrad0s(iobs)
              END IF
            END DO
            jobs=jact-jb
!
          ELSE
!
!  First grab the non-normalized Lanczos vector.
!
!  AUGMENTED.
!
            DO iobs=1,ndatum(ng)+1
!
!  Selime code: w=pgrad, t1=TLmodVal.
!
              pgrad(iobs)=zcglwk(iobs,innloop,outloop)
            END DO
!
!  Selime's code:  beta=pgrad*HBH'*pgrad=cg_beta
!
            cg_beta(innloop,outloop)=0.0_dp
            DO iobs=1,ndatum(ng)
              cg_beta(innloop,outloop)=cg_beta(innloop,outloop)+        &
     &                                 pgrad(iobs)*tlmodval(iobs)
            END DO
!
!  AUGMENTED.
!
            DO iobs=1,ndatum(ng)
              cg_beta(innloop,outloop)=cg_beta(innloop,outloop)+        &
     &                                 hbk(iobs,outloop)*               &
     &                                 vcglwk(ndatum(ng)+1,innloop,     &
     &                                        outloop)*                 &
     &                                 pgrad(iobs)+                     &
     &                                 hbk(iobs,outloop)*               &
     &                                 pgrad(ndatum(ng)+1)*             &
     &                                 vcglwk(iobs,innloop,outloop)
            END DO
            cg_beta(innloop,outloop)=cg_beta(innloop,outloop)+          &
     &                               jb0(outloop-1)*                    &
     &                               pgrad(ndatum(ng)+1)*               &
     &                               vcglwk(ndatum(ng)+1,innloop,       &
     &                                      outloop)
!
            cg_beta(innloop,outloop)=sqrt(cg_beta(innloop,outloop))
!
!  Normalize TLmodVal here and save in TLmodVal_S.
!
!  AUGMENTED
!
            DO iobs=1,ndatum(ng)+1
!
!  Selime code: v1=w/beta=zcglwk
!
              zcglwk(iobs,innloop,outloop)=pgrad(iobs)/                 &
     &                                     cg_beta(innloop,outloop)
!
!  Selime code: z1=w/beta=vcglwk
!
              vcglwk(iobs,innloop,outloop)=vcglwk(iobs,innloop,outloop)/&
     &                                     cg_beta(innloop,outloop)
            END DO
!
!  Selime code: t1=t1/beta=TLmodVal
!
!  AUGMENTED
!
            DO iobs=1,ndatum(ng)
              tlmodval_s(iobs,innloop,outloop)=tlmodval(iobs)
              tlmodval(iobs)=tlmodval(iobs)/cg_beta(innloop,outloop)
            END DO
!
!  Selime's code:  cg_QG=zgrad0*HBH'*zcglwk
!
            cg_qg(innloop,outloop)=0.0_dp
            DO iobs=1,ndatum(ng)
              cg_qg(innloop,outloop)=cg_qg(innloop,outloop)+            &
     &                               tlmodval(iobs)*                    &
     &                               vgrad0(iobs)
            END DO
!
!  AUGMENTED.
!
            DO iobs=1,ndatum(ng)
              cg_qg(innloop,outloop)=cg_qg(innloop,outloop)+            &
     &                               hbk(iobs,outloop)*                 &
     &                               vcglwk(ndatum(ng)+1,innloop,       &
     &                                      outloop)*                   &
     &                               vgrad0(iobs)+                      &
     &                               hbk(iobs,outloop)*                 &
     &                               vcglwk(iobs,innloop,outloop)*      &
     &                               vgrad0(ndatum(ng)+1)
            END DO
            cg_qg(innloop,outloop)=cg_qg(innloop,outloop)+              &
     &                             jb0(outloop-1)*                      &
     &                             vcglwk(ndatum(ng)+1,innloop,         &
     &                                    outloop)*                     &
     &                             vgrad0(ndatum(ng)+1)
!TEST
!TEST       IF (innLoop.gt.1) cg_QG(innLoop,outLoop)=0.0_dp
!TEST
!
!  Calculate the new solution based upon the new, orthonormalized
!  Lanczos vector. First, the tridiagonal system is solved by
!  decomposition and forward/back substitution.
!
            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
!
!  AUGMENTED.
!
            DO iobs=1,ndatum(ng)+1
              zw(iobs)=zgrad0(iobs,outloop)+                            &
     &                 cg_beta(innloop,outloop)*                        &
     &                 vcglwk(iobs,innloop,outloop)*                    &
     &                 zwork(innloop-1,3)
            END DO
 
!
!  AUGMENTED.
!
            DO iobs=1,ndatum(ng)+1
              px(iobs)=0.0_dp
              DO ivec=1,innloop-1
                px(iobs)=px(iobs)+                                      &
     &                   zcglwk(iobs,ivec,outloop)*zwork(ivec,3)
                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=1                 ! Spectral LMP
              laug=.true.
              CALL rprecond (ng, lscale, laug, outloop, ninnloop-1, px)
            END IF
!
!  Compute the reduction in the gradient Ax-b (in y-space).
!
!AMM        preducy=0.0_dp
!AMM        DO iobs=1,Ndatum(ng)
!AMM          preducy=preducy+zw(iobs)*zw(iobs)
!AMM        END DO
!AMM        preducy=SQRT(preducy)/cg_Gnorm_y(outLoop)
!
!  Compute the reduction in the gradient Ax-b (in v-space).
!
            IF (lprecond.and.(outloop.gt.1)) THEN
              lscale=1                 ! Spectral LMP
              laug=.true.
              CALL rprecond (ng, lscale, laug, outloop, ninnloop-1, zw)
            END IF
!
            preducv=0.0_dp
!
!  AUGMENTED.
!
            DO iobs=1,ndatum(ng)+1
              preducv=preducv+zw(iobs)*zw(iobs)
            END DO
            preducv=sqrt(preducv)/cg_gnorm_v(outloop)
!
!  Estimate the residual (in y-space) by:
!
!    ||Ax-b|| = cg_beta(k+1,outLoop)*zwork(innLoop,3)
!
          eps=abs(cg_beta(innloop,outloop)*zwork(innloop-1,3))
!
!  Estimate the various penalty function components based on
!  Bennett (2002), page 22 and pages 42-44:
!
!     Jopt  = value of the penalty function when true px
!             has been identified.
!     Jf    = initial data misfit for the first-guess
!     Jdata = data misfit for the state estimate
!     Jmod  = the model penalty function (the sum of the
!             dynamical, initial and boundary penalties)
!     Jb    = the ACTUAL model penalty function.
!     Jobs  = the ACTUAL data penalty function.
!     Jact  = the ACTUAL total penalty function.
!
!  NOTE:  Unlike I4D-Var where we compute the actual cost function
!         for each iteration that decreases with increasing
!         number of iterations, the various penalty terms
!         represent the optimal values (i.e. the values when
!         optimal state estimate has been identified). Therefore
!         the various penalty terms will only be correct when the
!         optimal state has been found. Therefore, as the
!         R4D-Var proceeds, Jopt, Jdata and Jmod will asymptote
!         to their final values and WILL NOT necessary decrease
!         with increasing iteration number.
!
            IF (outloop.eq.1) THEN
              DO iobs=1,ndatum(ng)
                IF (obserr(iobs).ne.0.0_r8) THEN
!                 Jopt=Jopt+px(iobs)*vgrad0(iobs)/SQRT(ObsErr(iobs))
                  jopt=jopt+px(iobs)*tlmodval_s(iobs,1,outloop)
!                 Jf=Jf+vgrad0(iobs)*vgrad0(iobs)/ObsErr(iobs)
                  jf=jf+vgrad0(iobs)*tlmodval_s(iobs,1,outloop)
                END IF
!AMM: For Jdata we need GDG'*px but we don't have this except at the end
!AMM            Jdata=Jdata+px(iobs)*px(iobs)
              END DO
!AMM          Jmod=Jopt-Jdata
            ELSE
              DO iobs=1,ndatum(ng)
!AMM            Jopt=Jopt-                                              &
!AMM &               (px(iobs)+cg_pxsave(iobs))*cg_innov(iobs)
                IF (obserr(iobs).ne.0.0_r8) THEN
                  jf=jf+cg_innov(iobs)*cg_innov(iobs)
                END IF
!AMM            Jdata=Jdata+                                            &
!AMM &                (px(iobs)+cg_pxsave(iobs))*                       &
!AMM &                (px(iobs)+cg_pxsave(iobs))
              END DO
              jmod=jopt-jdata
            END IF
!
!  Compute the actual penalty function terms.
!
            DO ivec=1,innloop-1
              IF (ivec.eq.1) THEN
                zfact(ivec)=1.0_dp/cg_gnorm_v(outloop)
              ELSE
                zfact(ivec)=1.0_dp/cg_beta(ivec,outloop)
              ENDIF
            END DO
!
!  AUGMENTED.
!
            DO iobs=1,ndatum(ng)+1
              zd(iobs)=vgrad0s(iobs)
            END DO
!
!  AUGMENTED.
!
            DO iobs=1,ndatum(ng)+1
              zhv(iobs)=zd(iobs)
            END DO
            DO ivec=1,innloop-1
              ztemp1(ivec)=0.0_dp
              ztemp3(ivec)=0.0_dp
              ztemp4(ivec)=0.0_dp
              DO iobs=1,ndatum(ng)
                ztemp1(ivec)=ztemp1(ivec)+                              &
     &                       tlmodval_s(iobs,ivec,outloop)*             &
     &                       zhv(iobs)*                                 &
     &                       zfact(ivec)
                ztemp3(ivec)=ztemp3(ivec)+                              &
     &                       vcglwk(iobs,ivec,outloop)*                 &
     &                       hbk(iobs,outloop)
                ztemp4(ivec)=ztemp4(ivec)+                              &
     &                       tlmodval_s(iobs,ivec,outloop)*             &
     &                       zgrad0(iobs,outloop)*                      &
     &                       zfact(ivec)
              END DO
!
!  AUGMENTED
!  Note the factor zfact is not needed in the next loop since it has
!  already been applied to vcglwk above.
!
              DO iobs=1,ndatum(ng)
                ztemp1(ivec)=ztemp1(ivec)+                              &
     &                       zhv(iobs)*                                 &
     &                       hbk(iobs,outloop)*                         &
     &                       vcglwk(ndatum(ng)+1,ivec,outloop)+         &
     &                       hbk(iobs,outloop)*                         &
     &                       vcglwk(iobs,ivec,outloop)*                 &
     &                       zhv(ndatum(ng)+1)
                ztemp4(ivec)=ztemp4(ivec)+                              &
     &                       zgrad0(iobs,outloop)*                      &
     &                       hbk(iobs,outloop)*                         &
     &                       vcglwk(ndatum(ng)+1,ivec,outloop)+         &
     &                       hbk(iobs,outloop)*                         &
     &                       vcglwk(iobs,ivec,outloop)*                 &
     &                       zgrad0(ndatum(ng)+1,outloop)
              END DO
              ztemp1(ivec)=ztemp1(ivec)+                                &
     &                     zhv(ndatum(ng)+1)*                           &
     &                     jb0(outloop-1)*                              &
     &                     vcglwk(ndatum(ng)+1,ivec,outloop)
              ztemp3(ivec)=ztemp3(ivec)+                                &
     &                     jb0(outloop-1)*                              &
     &                     vcglwk(ndatum(ng)+1,ivec,outloop)
              ztemp4(ivec)=ztemp4(ivec)+                                &
     &                     zgrad0(ndatum(ng)+1,outloop)*                &
     &                     jb0(outloop-1)*                              &
     &                     vcglwk(ndatum(ng)+1,ivec,outloop)
!TEST
              ztemp4(ivec)=0.0_dp
!TEST
            END DO
            zbet=cg_delta(1,outloop)
!TEST       zu1(1)=ztemp1(1)/zbet
!TEST       zu1(1)=ztemp4(1)/zbet
            zu1(1)=cg_qg(1,outloop)/zbet
!TEST
            DO ivec=2,innloop-1
              zgam(ivec)=cg_beta(ivec,outloop)/zbet
              zbet=cg_delta(ivec,outloop)-                              &
     &             cg_beta(ivec,outloop)*zgam(ivec)
!TEST         zu1(ivec)=(ztemp1(ivec)-                                  &
              zu1(ivec)=(ztemp4(ivec)-                                  &
     &                   cg_beta(ivec,outloop)*zu1(ivec-1))/zbet
            END DO
            ztemp2(innloop-1)=zu1(innloop-1)
 
            DO ivec=innloop-2,1,-1
              zu1(ivec)=zu1(ivec)-zgam(ivec+1)*zu1(ivec+1)
              ztemp2(ivec)=zu1(ivec)
            END DO
!
            jb=jb0(outloop-1)
            jact=jb
            DO iobs=1,ndatum(ng)
              IF (obserr(iobs).ne.0.0_r8) THEN
                jact=jact+obserr(iobs)*vgrad0s(iobs)*vgrad0s(iobs)
              END IF
            END DO
            DO ivec=1,innloop-1
              jb=jb+ztemp2(ivec)*ztemp2(ivec)
!TEST         Jact=Jact-ztemp1(ivec)*ztemp2(ivec)+                      &
!    &             2.0_dp*ztemp2(ivec)*ztemp3(ivec)
              jact=jact-ztemp1(ivec)*ztemp2(ivec)
            END DO
            DO iobs=1,ndatum(ng)
              jb=jb-2.0_dp*px(iobs)*hbk(iobs,outloop)
!TEST
              jact=jact+2.0_dp*px(iobs)*hbk(iobs,outloop)
!TEST
            END DO
            jb=jb-2.0_dp*px(ndatum(ng)+1)*jb0(outloop-1)
!TEST
            jact=jact+2.0_dp*px(ndatum(ng)+1)*jb0(outloop-1)
!TEST
            jobs=jact-jb
 
          END IF
!
!  NEXT CALCULATION
!
!  After the initialization, for every other inner loop, calculate a
!  new Lanczos vector, store it in the matrix, and update search.
!
!  AUGMENTED.
!
          DO iobs=1,ndatum(ng)
!
!  Selime code: w=t1.
!
            pgrad(iobs)=tlmodval(iobs)+                                 &
     &                  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
!
!  Selime code: w=t1/R. No augmented term here since Raug^-1=[R^-1 0].
!
              pgrad(iobs)=pgrad(iobs)/obserr(iobs)
            END IF
          END DO
          pgrad(ndatum(ng)+1)=0.0_dp
!
!  Selime code: w=t1/R+z1
!
          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
!
!  Selime code: w=w-v0*beta
!
!  AUGMENTED.
!
            DO iobs=1,ndatum(ng)+1
              pgrad(iobs)=pgrad(iobs)-                                  &
     &                    cg_beta(innloop,outloop)*                     &
     &                    zcglwk(iobs,innloop-1,outloop)
            END DO
          END IF
!
          cg_delta(innloop,outloop)=0.0_dp
!
!  Selime's code: alpha=w'*t1=cg_delta
!
!  AUGMENTED.
!
          DO iobs=1,ndatum(ng)
            cg_delta(innloop,outloop)=cg_delta(innloop,outloop)+        &
     &                                tlmodval(iobs)*                   &
     &                                pgrad(iobs)+                      &
     &                                hbk(iobs,outloop)*                &
     &                                vcglwk(ndatum(ng)+1,innloop,      &
     &                                       outloop)*                  &
     &                                pgrad(iobs)+                      &
     &                                hbk(iobs,outloop)*                &
     &                                vcglwk(iobs,innloop,outloop)*     &
     &                                pgrad(ndatum(ng)+1)
          END DO
          cg_delta(innloop,outloop)=cg_delta(innloop,outloop)+          &
     &                              jb0(outloop-1)*                     &
     &                              vcglwk(ndatum(ng)+1,innloop,        &
     &                                     outloop)*                    &
     &                              pgrad(ndatum(ng)+1)
!
!  Exit, if not a positive definite algorithm.
!
          IF (cg_delta(innloop,outloop).le.0.0_dp) THEN
            WRITE (stdout,20) cg_delta(innloop,outloop), outloop,       &
     &                        innloop
            exit_flag=8
            GO TO 10
          END IF
!
!  Selime code: w=w-alpha*v1
!
!  AUGMENTED.
!
          DO iobs=1,ndatum(ng)+1
            pgrad(iobs)=pgrad(iobs)-cg_delta(innloop,outloop)*          &
     &                  zcglwk(iobs,innloop,outloop)
          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_dp/cg_gnorm_v(outloop)
            ELSE
              zfact(ivec)=1.0_dp/cg_beta(ivec,outloop)
            ENDIF
          END DO
!
!  Selime's code: cg_dla=pgrad*HBH'*v1=pgrad*TLmodVal_S/zfact
!
!  AUGMENTED
!
          DO ivec=innloop,1,-1
            cg_dla(ivec,outloop)=0.0_dp
            DO iobs=1,ndatum(ng)
              cg_dla(ivec,outloop)=cg_dla(ivec,outloop)+                &
     &                             pgrad(iobs)*                         &
     &                             tlmodval_s(iobs,ivec,outloop)*       &
     &                             zfact(ivec)+pgrad(iobs)*             &
     &                             hbk(iobs,outloop)*                   &
     &                             vcglwk(ndatum(ng)+1,ivec,outloop)+   &
     &                             hbk(iobs,outloop)*                   &
     &                             vcglwk(iobs,ivec,outloop)*           &
     &                             pgrad(ndatum(ng)+1)
            END DO
            cg_dla(ivec,outloop)=cg_dla(ivec,outloop)+                  &
     &                           jb0(outloop-1)*                        &
     &                           vcglwk(ndatum(ng)+1,ivec,outloop)*     &
     &                           pgrad(ndatum(ng)+1)
            DO iobs=1,ndatum(ng)+1
              pgrad(iobs)=pgrad(iobs)-                                  &
     &                    cg_dla(ivec,outloop)*                         &
     &                    vcglwk(iobs,ivec,outloop)
            END DO
          END DO
!
!  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
            zcglwk(iobs,innloop+1,outloop)=pgrad(iobs)
          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
!
!  Selime code: z1=vcglwk.
!
!  AUGMENTED.
!
          DO iobs=1,ndatum(ng)+1
            vcglwk(iobs,innloop+1,outloop)=pgrad(iobs)
          END DO
!
!  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)
              admodval(iobs)=px(iobs)
            END DO
            DO iobs=1,ndatum(ng)+1
              fourdvar(ng)%cg_pxsum(iobs)=fourdvar(ng)%cg_pxsum(iobs)+  &
     &                                    fourdvar(ng)%cg_pxsave(iobs)
              fourdvar(ng)%cg_pxsave(iobs)=px(iobs)
!AMM          FOURDVAR(ng)%cg_pxsum(iobs)=FOURDVAR(ng)%cg_pxsum(iobs)+  &
!AMM &                                    px(iobs)
            END DO
          ELSE
!
!  Selime code: z1=ADmodVal.
!
            DO iobs=1,ndatum(ng)
              admodval(iobs)=vcglwk(iobs,innloop+1,outloop)
            END DO
!
          END IF
!
        END IF minimize
!
!-----------------------------------------------------------------------
!  Report minimization progress.
!-----------------------------------------------------------------------
!
        IF (innloop.gt.0) THEN
          WRITE (stdout,30) ndatum(ng),                                 &
     &                      outloop, innloop-1, eps,                    &
     &                      outloop, innloop-1, preducy,                &
     &                      outloop, innloop-1, preducv,                &
     &                      outloop, innloop-1, jf,                     &
     &                      outloop, innloop-1, jdata,                  &
     &                      outloop, innloop-1, jmod,                   &
     &                      outloop, innloop-1, jopt,                   &
     &                      outloop, innloop-1, jb,                     &
     &                      outloop, innloop-1, jobs,                   &
     &                      outloop, innloop-1, jact,                   &
     &                      outloop, innloop-1
        END IF
        DO ivec=1,innloop
          WRITE (stdout,40) ivec, cg_delta(ivec,outloop),               &
     &                      cg_beta(ivec,outloop), zwork(ivec,3)
        END DO
        WRITE (stdout,50) innloop+1, cg_beta(innloop+1,outloop)
!
!  If preconditioning, report the Ritz eigenvalues used.
!
        IF ((lhessianev.or.lprecond).and.(outloop.gt.1)) THEN
          IF (nritzev.gt.0) THEN
            WRITE (stdout,60) outloop, innloop, nconvritz
            ic=0
            DO ivec=1,ninnloop-1
              IF (ivec.gt.(ninnloop-1-nconvritz)) THEN
                string='         '
                ic=ic+1
                WRITE (stdout,70) ivec, cg_ritz(ivec,outloop-1),        &
     &                            cg_ritzerr(ivec,outloop-1),           &
     &                            trim(adjustl(string)),                &
     &                            'Used=', ic
              ELSE
                string='             '
                WRITE (stdout,80) ivec, cg_ritz(ivec,outloop-1),        &
     &                            cg_ritzerr(ivec,outloop-1),           &
     &                            trim(adjustl(string))
              END IF
            END DO
          ELSE
            WRITE (stdout,90) outloop, innloop, ritzmaxerr
            ic=0
            DO ivec=1,ninnloop
              IF (cg_ritzerr(ivec,outloop-1).le.ritzmaxerr) THEN
                string='converged'
                ic=ic+1
                WRITE (stdout,70) ivec, cg_ritz(ivec,outloop-1),        &
     &                            cg_ritzerr(ivec,outloop-1),           &
     &                            trim(adjustl(string)),                &
     &                            'Good=', ic
              ELSE
                string='not converged'
                WRITE (stdout,80) ivec, cg_ritz(ivec,outloop-1),        &
     &                            cg_ritzerr(ivec,outloop-1),           &
     &                            trim(adjustl(string))
              END IF
            END DO
          END IF
        END IF
!
!-----------------------------------------------------------------------
!  If last inner loop, innLoop = NinnLoop.
!-----------------------------------------------------------------------
!
!  Compute the eigenvalues and eigenvectors of the tridiagonal matrix.
!
        IF (innloop.eq.ninnloop) THEN
          IF (lhessianev.or.lprecond) THEN
            DO ivec=1,innloop-1
              cg_ritz(ivec,outloop)=cg_delta(ivec,outloop)
            END DO
            DO ivec=1,innloop-2
              zwork(ivec,1)=cg_beta(ivec+1,outloop)
            END DO
!
!  Use the LAPACK routine DSTEQR to compute the eigenvectors and
!  eigenvalues of the tridiagonal matrix. If applicable, the
!  eigenpairs are computed by master thread only.
!
            CALL dsteqr ('I', innloop-1, cg_ritz(1,outloop), zwork(1,1),&
     &                   zgv, ninner, work, info)
            IF (info.ne.0) THEN
              WRITE (stdout,*) ' RPCG_LANCZOS - Error in DSTEQR:',      &
     &                         ' info = ',info
              exit_flag=8
              GO TO 10
            END IF
!
            DO j=1,ninner
              DO i=1,ninner
                cg_zv(i,j,outloop)=zgv(i,j)
              END DO
            END DO
!
!  Estimate the Ritz value error bounds.
!
            DO ivec=1,innloop-1
              cg_ritzerr(ivec,outloop)=abs(cg_beta(innloop,outloop)*    &
     &                                 cg_zv(innloop-1,ivec,outloop))
            END DO
!
!  Check for exploding or negative Ritz values.
!
            DO ivec=1,innloop-1
              IF (cg_ritzerr(ivec,outloop).lt.0.0_dp) THEN
                WRITE (stdout,*) ' RPCG_LANCZOS - negative Ritz value', &
     &                           ' found.'
                exit_flag=8
                GO TO 10
              END IF
            END DO
!
!  Change the error bounds for the acceptable eigenvectors.
!
            DO ivec=1,ninnloop-1
              cg_ritzerr(ivec,outloop)=cg_ritzerr(ivec,outloop)/        &
     &                                 cg_ritz(ninnloop-1,outloop)
            END DO
!
!  Report Eigenvalues and their relative accuracy.
!
            WRITE (stdout,100) outloop, innloop, ritzmaxerr
            ic=0
            DO ivec=1,ninnloop-1
              IF (cg_ritzerr(ivec,outloop).le.ritzmaxerr) THEN
                string='converged'
                ic=ic+1
                WRITE (stdout,70) ivec, cg_ritz(ivec,outloop),          &
     &                            cg_ritzerr(ivec,outloop),             &
     &                            trim(adjustl(string)),                &
     &                            'Good=', ic
              ELSE
                string='not converged'
                WRITE (stdout,80) ivec, cg_ritz(ivec,outloop),          &
     &                            cg_ritzerr(ivec,outloop),             &
     &                            trim(adjustl(string))
              END IF
            END DO
!
!  Calculate the converged eigenvectors (vcglev) of the [R_n + Cobs].
!
            CALL rpevecs (ng, outloop, ninnloop-1)
          END IF
        END IF
 
      END IF master_thread
!
!-----------------------------------------------------------------------
!  Come here if not a possite definite algorithm.
!-----------------------------------------------------------------------
!
 10   CONTINUE
# ifdef DISTRIBUTE
      CALL mp_bcasti (ng, model, exit_flag)
# endif
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
# ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Broadcast new solution to other nodes.
!-----------------------------------------------------------------------
!
      CALL mp_bcastf (ng, model, admodval)
#  if defined RBL4DVAR          || defined R4DVAR      || \
      defined sensitivity_4dvar || defined tl_rbl4dvar || \
      defined tl_r4dvar
!!    CALL mp_bcastf (ng, model, TLmodVal_S(:,:,outLoop))   ! not needed
      CALL mp_bcastf (ng, model, hbk(:,outloop))
      CALL mp_bcastf (ng, model, jb0(outloop-1))
#  endif
      CALL mp_bcasti (ng, model, info)
      CALL mp_bcastf (ng, model, cg_beta(:,outloop))
      CALL mp_bcastf (ng, model, cg_qg(:,outloop))
      CALL mp_bcastf (ng, model, cg_delta(:,outloop))
      CALL mp_bcastf (ng, model, cg_dla(:,outloop))
      CALL mp_bcastf (ng, model, cg_gnorm_y(:))
      CALL mp_bcastf (ng, model, cg_gnorm_v(:))
      CALL mp_bcastf (ng, model, fourdvar(ng)%cg_pxsave(:))
!AMM  CALL mp_bcastf (ng, model, cg_pxsave(:))
      CALL mp_bcastf (ng, model, cg_innov(:))
      CALL mp_bcastf (ng, model, zgrad0(:,outloop))
      CALL mp_bcastf (ng, model, zcglwk(:,:,outloop))
      CALL mp_bcastf (ng, model, vcglwk(:,:,outloop))
      IF ((lhessianev.or.lprecond).and.(innloop.eq.ninnloop)) THEN
        CALL mp_bcastf (ng, model, cg_ritz(:,outloop))
        CALL mp_bcastf (ng, model, cg_ritzerr(:,outloop))
        CALL mp_bcastf (ng, model, cg_zv(:,:,outloop))
        CALL mp_bcastf (ng, model, vcglev(:,:,outloop))
      END IF
      CALL mp_bcastf (ng, model, jf)
      CALL mp_bcastf (ng, model, jdata)
      CALL mp_bcastf (ng, model, jmod)
      CALL mp_bcastf (ng, model, jopt)
      CALL mp_bcastf (ng, model, jobs)
      CALL mp_bcastf (ng, model, jact)
      CALL mp_bcastf (ng, model, preducv)
      CALL mp_bcastf (ng, model, preducy)
# endif
!
!-----------------------------------------------------------------------
!  Write out conjugate gradient vectors into 4D-Var NetCDF file.
!-----------------------------------------------------------------------
!
# ifdef BGQC
!
!  Save the cost function values for background quality control decisions
!
      jact_s(innloop)=jact
!
# endif
      IF (innloop.gt.0) THEN
         CALL cg_write_rpcg (ng, model, innloop, outloop,               &
     &                       jf, jdata, jmod, jopt, jb, jobs, jact,     &
     &                       preducv, preducy)
      END IF
 
# ifdef PROFILE
!
!  Turn off time clock.
!
      CALL wclock_off (ng, model, 85, __line__, myfile)
# endif
!
 20   FORMAT (/,' RPCG_LANCZOS - Fatal error, not possitive definite',  &
     &        ' algorithm:',/,                                          &
     &        /,11x,'cg_delta = ',1p,e15.8,0p,3x,'(',i3.3,', ',i3.3,')')
 30   FORMAT (/,' RPCG_LANCZOS - Conjugate Gradient Information:',/,    &
     &        /,11x,'Ndatum     = ',i0,/,  /,                           &
     &        1x,'(',i3.3,',',i3.3,'): ',                               &
     &        'Residual estimate,                      eps = ',         &
     &        1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ',                    &
     &        'Reduction in gradient norm,  Greduc y-space = ',         &
     &        1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ',                    &
     &        'Reduction in gradient norm,  Greduc v-space = ',         &
     &        1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ',                    &
     &        'First guess initial data misfit,         Jf = ',         &
     &        1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ',                    &
     &        'State estimate data misfit,           Jdata = ',         &
     &        1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ',                    &
     &        'Model penalty function,                Jmod = ',         &
     &        1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ',                    &
     &        'Optimal penalty function,              Jopt = ',         &
     &        1p,e14.7,/,/,1x,'(',i3.3,',',i3.3,'): ',                  &
     &        'Actual Model penalty function,           Jb = ',         &
     &        1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ',                    &
     &        'Actual data penalty function,          Jobs = ',         &
     &        1p,e14.7,/,/,1x,'(',i3.3,',',i3.3,'): ',                  &
     &        'Actual total penalty function,         Jact = ',         &
     &        1p,e14.7,/,/,1x,'(',i3.3,',',i3.3,'): ',                  &
     &        'Lanczos vectors - cg_delta, cg_beta, zwork:',/)
 40   FORMAT (6x,i3.3,4x,3(1p,e15.8,5x))
 50   FORMAT (6x,i3.3,24x,1p,e15.8)
 60   FORMAT (/,1x,'(',i3.3,',',i3.3,'): ',                             &
     &        'Ritz eigenvalues used in preconditioning, ',             &
     &        'nConvRitz = ',i3.3,/)
 70   FORMAT (6x,i3.3,2x,1p,e14.7,2x,1p,e14.7,2x,a,5x,'('a,i3.3,')')
 80   FORMAT (6x,i3.3,2x,1p,e14.7,2x,1p,e14.7,2x,a)
 90   FORMAT (/,1x,'(',i3.3,',',i3.3,'): ',                             &
     &        'Ritz eigenvalues used in preconditioning, ',             &
     &        'RitzMaxErr = ',1p,e12.5,/)
100   FORMAT (/,1x,'(',i3.3,',',i3.3,'): ',                             &
     &        'New Ritz eigenvalues and their accuracy,  ',             &
     &        'RitzMaxErr = ',1p,e12.5,/)
!
      RETURN
      END SUBROUTINE rpcg_lanczos
!
!***********************************************************************
      SUBROUTINE rpevecs (ng, outLoop, NinnLoop)
!***********************************************************************
!                                                                      !
!  This routine computes the converged eigenvectors of (R_n+Cobs).     !
!                                                                      !
!***********************************************************************
!
!  Imported variable declarations
!
      integer, intent(in) :: ng, outloop, ninnloop
!
!  Local variable declarations.
!
      integer :: iobs, ivec, jn, jk, ij, jkk, ingood
!
      real(dp) :: dla, zfact
      real(dp), dimension(Ndatum(ng)+1) :: zt
!
!-----------------------------------------------------------------------
!  Compute and store converged eigenvectors of (R_n+Cobs).
!-----------------------------------------------------------------------
!
!  Count and collect the converged eigenvalues.
!
      ingood=0
      DO ivec=ninnloop,1,-1
        ingood=ingood+1
        ritz(ingood)=cg_ritz(ivec,outloop)
      END DO
!
!  Calculate the converged eigenvectors of (R_n+Cobs).
!
      ij=0
      DO jn=ninnloop,1,-1
        ij=ij+1
!                                                             AUGMENTED
        DO iobs=1,ndatum(ng)
          vcglev(iobs,ij,outloop)=0.0_dp
        END DO
        DO jk=1,ninnloop
          DO iobs=1,ndatum(ng)
            vcglev(iobs,ij,outloop)=vcglev(iobs,ij,outloop)+            &
     &                              (vcglwk(iobs,jk,outloop)-           &
     &                               fourdvar(ng)%cg_pxsum(iobs)*       &
     &                               vcglwk(ndatum(ng)+1,jk,outloop))*  &
     &                              cg_zv(jk,jn,outloop)
          END DO
        END DO
!
!  Orthonormalize.
!
        DO jk=1,ij-1
          DO iobs=1,ndatum(ng)
            zt(iobs)=0.0_dp
          END DO
          DO iobs=1,ndatum(ng)
            DO jkk=1,ninnloop
              IF (jkk.eq.1) THEN
                zfact=1.0_dp/cg_gnorm_v(outloop)
              ELSE
                zfact=1.0_dp/cg_beta(jkk,outloop)
              ENDIF
              zt(iobs)=zt(iobs)+                                        &
     &                 (tlmodval_s(iobs,jkk,outloop)*                   &
     &                  zfact+                                          &
     &                  hbk(iobs,outloop)*                              &
     &                  vcglwk(ndatum(ng)+1,jkk,outloop))*              &
     &                 cg_zv(jkk,jk,outloop)
            END DO
          END DO
          dla=0.0_dp
          DO iobs=1,ndatum(ng)
            dla=dla+zt(iobs)*vcglev(iobs,ij,outloop)
          END DO
          DO iobs=1,ndatum(ng)
            vcglev(iobs,ij,outloop)=vcglev(iobs,ij,outloop)-            &
     &                              dla*zt(iobs)
          END DO
        END DO
        dla=0.0_dp
!                                                            AUGMENTED
        DO iobs=1,ndatum(ng)
          dla=dla+vcglev(iobs,ij,outloop)*vcglev(iobs,ij,outloop)
        END DO
!                                                            AUGMENTED+1
        DO iobs=1,ndatum(ng)
          vcglev(iobs,ij,outloop)=vcglev(iobs,ij,outloop)/sqrt(dla)
        END DO
      END DO
!
      RETURN
      END SUBROUTINE rpevecs
!
!***********************************************************************
      SUBROUTINE rprecond (ng, Lscale, Laug, outLoop, NinnLoop, py)
!***********************************************************************
!                                                                      !
!  This routine is the preconditioner in observation space.            !
!                                                                      !
!***********************************************************************
!
!  Imported variable declarations
!
      logical, intent(in) :: laug
      integer, intent(in) :: ng, lscale, outloop, ninnloop
!
      real(dp), intent(inout) :: py(ndatum(ng)+1)            ! AUGMENTED
!
!  Local variable declarations.
!
      integer :: is, ie, inc, iss, i, jk
      integer :: nol, nols, nole, ninc
      integer :: ingood
      integer :: iobs, nvec, ivec
!
      real(dp) :: dla, dlar, fac, facritz, zfact
      real(dp), dimension(NinnLoop) :: zvect
      real(dp), dimension(Ndatum(ng)+1) :: zt
!
!  Set the do-loop indices for the sequential preconditioner
!  loop.
!
      nols=1
      nole=outloop-1
      ninc=1
!
!  Apply the preconditioners derived from all previous outer-loops
!  sequentially.
!
      DO nol=nols,nole,ninc
!
!  Determine the number of Ritz vectors to use.
!  For Lritz=.TRUE., choose HvecErr to be larger.
!
        ingood=0
        DO i=1,ninnloop
          IF (cg_ritzerr(i,nol).le.ritzmaxerr) THEN
            ingood=ingood+1
          END IF
        END DO
        IF (nritzev.gt.0) THEN
          nconvritz=nritzev
          ingood=nconvritz
        ELSE
          nconvritz=ingood
        END IF
!
        IF (lscale.gt.0) THEN
          is=1
          ie=ingood
          inc=1
        ELSE
          is=ingood
          ie=1
          inc=-1
        END IF
!
!  Lscale determines the form of the preconditioner:
!
!     1= LMP
!    -1= Inverse LMP
!
        DO nvec=is,ie,inc
          DO iobs=1,ndatum(ng)
            zt(iobs)=0.0_dp
          END DO
          DO iobs=1,ndatum(ng)
            DO jk=1,ninnloop
              IF (jk.eq.1) THEN
                zfact=1.0_dp/cg_gnorm_v(nol)
              ELSE
                zfact=1.0_dp/cg_beta(jk,nol)
              ENDIF
              zt(iobs)=zt(iobs)+                                        &
     &                 (tlmodval_s(iobs,jk,nol)*                        &
     &                  zfact+                                          &
     &                  hbk(iobs,nol)*                                  &
     &                  vcglwk(ndatum(ng)+1,jk,nol))*                   &
     &                 cg_zv(jk,nvec,nol)
            END DO
          END DO
!
          dla=0.0_dp
          DO iobs=1,ndatum(ng)
            dla=dla+zt(iobs)*py(iobs)
          END DO
          IF (lritz) THEN
            dlar=0.0_dp
            DO iobs=1,ndatum(ng)
              dlar=dlar+                                                &
     &             (tlmodval_s(iobs,ninnloop+1,nol)/                    &
     &              cg_beta(ninnloop+1,nol)+                            &
     &              hbk(iobs,nol)*                                      &
     &              vcglwk(ndatum(ng)+1,ninnloop+1,nol))*               &
     &             py(iobs)
            END DO
          END IF
!
!  NOTE: Lscale=1 and Lscale=-1 cases are not complete yet.
!
          IF (lscale.eq.-1) THEN
            fac=(cg_ritz(ninnloop+1-nvec,nol)-1.0_dp)*dla
          ELSE
            fac=(1.0_dp/cg_ritz(ninnloop+1-nvec,nol)-1.0_dp)*dla
          END IF
!
          DO iobs=1,ndatum(ng)
             py(iobs)=py(iobs)+fac*vcglev(iobs,nvec,nol)
          END DO
!
          IF (lritz) THEN
!
!  NOTE: We still need the code for Lscale=-1.
!
           IF (lscale.eq.1) THEN
             facritz=cg_beta(ninnloop+1,nol)*                           &
     &               cg_zv(ninnloop,ninnloop+1-nvec,nol)/               &
     &               cg_ritz(ninnloop+1-nvec,nol)
             DO iobs=1,ndatum(ng)
               py(iobs)=py(iobs)+                                       &
     &                  facritz*facritz*vcglev(iobs,nvec,nol)-          &
     &                  facritz*dlar*                                   &
     &                  (vcglwk(iobs,ninnloop+1,nol)-                   &
     &                   fourdvar(ng)%cg_pxsum(iobs)*                   &
     &                   vcglwk(ndatum(ng)+1,ninnloop+1,nol))
             END DO
             DO iobs=1,ndatum(ng)
               py(iobs)=py(iobs)-facritz*dla*vcglev(iobs,nvec,nol)
             END DO
            END IF
          END IF
!
        END DO
      END DO
!
      IF (laug) THEN
        DO iobs=1,ndatum(ng)
           py(iobs)=py(iobs)+                                           &
     &              fourdvar(ng)%cg_Dpxsum(iobs)*py(ndatum(ng)+1)
        END DO
      END IF
!
      RETURN
      END SUBROUTINE rprecond
!
!***********************************************************************
      SUBROUTINE cg_write_rpcg (ng, model, innLoop, outLoop,            &
     &                          Jf, Jdata, Jmod, Jopt, Jb, Jobs, Jact,  &
     &                          preducv, preducy)
!***********************************************************************
!                                                                      !
!  This routine writes conjugate gradient vectors into DAV file        !
!  using either the standard NetCDF library or the Parallel-IO (PIO)   !
!  library. It saves all the variables needed for split or restart     !
!  schemes.                                                            !
!                                                                      !
!***********************************************************************
!
!  Imported variable declarations
!
      integer, intent(in) :: ng, model, innloop, outloop
 
      real(dp), intent(in) :: jf, jdata, jmod, jopt, preducv, preducy
      real(dp), intent(in) :: jb, jobs, jact
!
!  Local variable declarations.
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", cg_write_rpcg"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Write out conjugate gradient variables.
!-----------------------------------------------------------------------
!
      SELECT CASE (dav(ng)%IOtype)
        CASE (io_nf90)
          CALL cg_write_rpcg_nf90 (ng, model, innloop, outloop,         &
     &                             jf, jdata, jmod, jopt, jb, jobs,     &
     &                             jact, preducv, preducy)
 
# if defined PIO_LIB && defined DISTRIBUTE
        CASE (io_pio)
          CALL cg_write_rpcg_pio (ng, model, innloop, outloop,          &
     &                            jf, jdata, jmod, jopt, jb, jobs,      &
     &                            jact, preducv, preducy)
# endif
        CASE DEFAULT
          IF (master) WRITE (stdout,10) dav(ng)%IOtype
          exit_flag=3
      END SELECT
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
  10  FORMAT (' CG_WRITE_RPCG - Illegal output file type, io_type = ',  &
     &        i0,/,17x,'Check KeyWord ''OUT_LIB'' in ''roms.in''.')
!
      RETURN
      END SUBROUTINE cg_write_rpcg
!
!***********************************************************************
      SUBROUTINE cg_write_rpcg_nf90 (ng, model, innLoop, outLoop,       &
     &                               Jf, Jdata, Jmod, Jopt, Jb, Jobs,   &
     &                               Jact, preducv, preducy)
!***********************************************************************
!
      USE mod_netcdf
!
!  Imported variable declarations
!
      integer, intent(in) :: ng, model, innloop, outloop
 
      real(dp), intent(in) :: jf, jdata, jmod, jopt, preducv, preducy
      real(dp), intent(in) :: jb, jobs, jact
!
!  Local variable declarations.
!
      integer :: nconv, status
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", cg_write_rpcg_nf90"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Write out conjugate gradient variables.
!-----------------------------------------------------------------------
!
!  Write out outer and inner iteration.
!
      CALL netcdf_put_ivar (ng, model, dav(ng)%name,                    &
     &                      'outer', outer,                             &
     &                      (/0/), (/0/),                               &
     &                      ncid = dav(ng)%ncid)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      CALL netcdf_put_ivar (ng, model, dav(ng)%name,                    &
     &                      'inner', inner,                             &
     &                      (/0/), (/0/),                               &
     &                      ncid = dav(ng)%ncid)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Write out number of converged Ritz eigenvalues.
!
      IF (innloop.eq.ninner) THEN
        CALL netcdf_put_ivar (ng, model, dav(ng)%name,                  &
     &                        'nConvRitz', nconvritz,                   &
     &                        (/outloop/), (/1/),                       &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Write out converged Ritz eigenvalues.
!
      IF (innloop.eq.ninner) THEN
        nconv=max(1,nconvritz)
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'Ritz', ritz(1:nconv),                    &
     &                        (/1,outloop/), (/nconv,1/),               &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Write out conjugate gradient norm.
!
      IF (innloop.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'cg_beta', cg_beta,                       &
     &                        (/1,1/), (/ninner+1,nouter/),             &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Write out Lanczos algorithm coefficients.
!
      IF (innloop.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'cg_delta', cg_delta,                     &
     &                        (/1,1/), (/ninner,nouter/),               &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Write normalization coefficients for Lanczos vectors.
!
      IF (innloop.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'cg_dla', cg_dla,                         &
     &                        (/1,1/), (/ninner,nouter/),               &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Initial gradient normalization factor.
!
      IF (innloop.eq.1) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'cg_Gnorm_y', cg_gnorm_y,                 &
     &                        (/1/), (/nouter/),                        &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'cg_Gnorm_v', cg_gnorm_v,                 &
     &                        (/1/), (/nouter/),                        &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Lanczos vector normalization factor.
!
      IF (innloop.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'cg_QG', cg_qg,                           &
     &                        (/1,1/), (/ninner+1,nouter/),             &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Reduction in the gradient norm.
!
      IF (innloop.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'cg_Greduc_y', preducy,                   &
     &                        (/innloop,outloop/), (/1,1/),             &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'cg_Greduc_v', preducv,                   &
     &                        (/innloop,outloop/), (/1,1/),             &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Eigenvalues of Lanczos recurrence relationship.
!
      IF (innloop.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'cg_Ritz', cg_ritz,                       &
     &                        (/1,1/), (/ninner,nouter/),               &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Eigenvalues relative error.
!
      IF (innloop.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'cg_RitzErr', cg_ritzerr,                 &
     &                        (/1,1/), (/ninner,nouter/),               &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Eigenvectors of Lanczos recurrence relationship.
!
      IF (innloop.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'cg_zv', cg_zv(:,:,outloop),              &
     &                        (/1,1,outloop/), (/ninner,ninner,1/),     &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  First guess initial data misfit.
!
      IF (innloop.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'Jf', jf,                                 &
     &                        (/innloop+1,outloop/), (/1,1/),           &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  State estimate data misfit.
!
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'Jdata', jdata,                           &
     &                        (/innloop+1,outloop/), (/1,1/),           &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Model penalty function.
!
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'Jmod', jmod,                             &
     &                        (/innloop+1,outloop/), (/1,1/),           &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Optimal penalty function.
!
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'Jopt', jopt,                             &
     &                        (/innloop+1,outloop/), (/1,1/),           &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Actual model penalty function.
!
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'Jb', jb,                                 &
     &                        (/innloop+1,outloop/), (/1,1/),           &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Actual data penalty function.
!
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'Jobs', jobs,                             &
     &                        (/innloop+1,outloop/), (/1,1/),           &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Actual total penalty function.
!
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'Jact', jact,                             &
     &                        (/innloop+1,outloop/), (/1,1/),           &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Initial gradient for minimization.
!
      IF (innloop.eq.1) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'zgrad0', zgrad0(:,outloop),              &
     &                        (/1,outloop/), (/ndatum(ng)+1,1/),        &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'vgrad0', vgrad0,                         &
     &                        (/1/), (/ndatum(ng)+1/),                  &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'Hbk', hbk(:,outloop),                    &
     &                        (/1,outloop/), (/ndatum(ng),1/),          &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Lanczos vectors.
!
      IF (innloop.gt.0) THEN
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'zcglwk', zcglwk(:,innloop,outloop),      &
     &                        (/1,innloop,outloop/),                    &
     &                        (/ndatum(ng)+1,1,1/),                     &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL netcdf_put_fvar (ng, model, dav(ng)%name,                  &
     &                        'vcglwk', vcglwk(:,innloop,outloop),      &
     &                        (/1,innloop,outloop/),                    &
     &                        (/ndatum(ng)+1,1,1/),                     &
     &                        ncid = dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Saved TLmodVal_S.
!
      CALL netcdf_put_fvar (ng, model, dav(ng)%name,                    &
     &                      'TLmodVal_S',                               &
     &                      tlmodval_s(:,innloop,outloop),              &
     &                      (/1,innloop,outloop/), (/ndatum(ng),1,1/),  &
     &                      ncid = dav(ng)%ncid)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!-----------------------------------------------------------------------
!  Synchronize model/observation NetCDF file to disk.
!-----------------------------------------------------------------------
!
      CALL netcdf_sync (ng, model, dav(ng)%name, dav(ng)%ncid)
!
      RETURN
      END SUBROUTINE cg_write_rpcg_nf90
 
# if defined PIO_LIB && defined DISTRIBUTE
!
!***********************************************************************
      SUBROUTINE cg_write_rpcg_pio (ng, model, innLoop, outLoop,        &
     &                              Jf, Jdata, Jmod, Jopt, Jb, Jobs,    &
     &                              Jact, preducv, preducy)
!***********************************************************************
!
      USE mod_pio_netcdf
!
!  Imported variable declarations
!
      integer, intent(in) :: ng, model, innloop, outloop
 
      real(dp), intent(in) :: jf, jdata, jmod, jopt, preducv, preducy
      real(dp), intent(in) :: jb, jobs, jact
!
!  Local variable declarations.
!
      integer :: nconv, status
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", cg_write_rpcg_pio"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Write out conjugate gradient variables.
!-----------------------------------------------------------------------
!
!  Write out outer and inner iteration.
!
      CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name,                &
     &                          'outer', outer,                         &
     &                          (/0/), (/0/),                           &
     &                          piofile = dav(ng)%pioFile)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
      CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name,                &
     &                          'inner', inner,                         &
     &                          (/0/), (/0/),                           &
     &                          piofile = dav(ng)%pioFile)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Write out number of converged Ritz eigenvalues.
!
      IF (innloop.eq.ninner) THEN
        CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name,              &
     &                            'nConvRitz', nconvritz,               &
     &                            (/outloop/), (/1/),                   &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Write out converged Ritz eigenvalues.
!
      IF (innloop.eq.ninner) THEN
        nconv=max(1,nconvritz)
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'Ritz', ritz(1:nconv),                &
     &                            (/1,outloop/), (/nconv,1/),           &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Write out conjugate gradient norm.
!
      IF (innloop.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'cg_beta', cg_beta,                   &
     &                            (/1,1/), (/ninner+1,nouter/),         &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Write out Lanczos algorithm coefficients.
!
      IF (innloop.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'cg_delta', cg_delta,                 &
     &                            (/1,1/), (/ninner,nouter/),           &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Write normalization coefficients for Lanczos vectors.
!
      IF (innloop.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'cg_dla', cg_dla,                     &
     &                            (/1,1/), (/ninner,nouter/),           &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Initial gradient normalization factor.
!
      IF (innloop.eq.1) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'cg_Gnorm_y', cg_gnorm_y,             &
     &                            (/1/), (/nouter/),                    &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'cg_Gnorm_v', cg_gnorm_v,             &
     &                            (/1/), (/nouter/),                    &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Lanczos vector normalization factor.
!
      IF (innloop.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'cg_QG', cg_qg,                       &
     &                            (/1,1/), (/ninner+1,nouter/),         &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Reduction in the gradient norm.
!
      IF (innloop.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'cg_Greduc_y', preducy,               &
     &                            (/innloop,outloop/), (/1,1/),         &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'cg_Greduc_v', preducv,               &
     &                            (/innloop,outloop/), (/1,1/),         &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Eigenvalues of Lanczos recurrence relationship.
!
      IF (innloop.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'cg_Ritz', cg_ritz,                   &
     &                            (/1,1/), (/ninner,nouter/),           &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Eigenvalues relative error.
!
      IF (innloop.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'cg_RitzErr', cg_ritzerr,             &
     &                            (/1,1/), (/ninner,nouter/),           &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Eigenvectors of Lanczos recurrence relationship.
!
      IF (innloop.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'cg_zv', cg_zv(:,:,outloop),          &
     &                            (/1,1,outloop/), (/ninner,ninner,1/), &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  First guess initial data misfit.
!
      IF (innloop.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'Jf', jf,                             &
     &                            (/innloop+1,outloop/), (/1,1/),       &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  State estimate data misfit.
!
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'Jdata', jdata,                       &
     &                            (/innloop+1,outloop/), (/1,1/),       &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Model penalty function.
!
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'Jmod', jmod,                         &
     &                            (/innloop+1,outloop/), (/1,1/),       &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Optimal penalty function.
!
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'Jopt', jopt,                         &
     &                            (/innloop+1,outloop/), (/1,1/),       &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Actual model penalty function.
!
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'Jb', jb,                             &
     &                            (/innloop+1,outloop/), (/1,1/),       &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Actual data penalty function.
!
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'Jobs', jobs,                         &
     &                            (/innloop+1,outloop/), (/1,1/),       &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Actual total penalty function.
!
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'Jact', jact,                         &
     &                            (/innloop+1,outloop/), (/1,1/),       &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Initial gradient for minimization.
!
      IF (innloop.eq.1) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'zgrad0', zgrad0(:,outloop),          &
     &                            (/1,outloop/), (/ndatum(ng)+1,1/),    &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'vgrad0', vgrad0,                     &
     &                            (/1/), (/ndatum(ng)+1/),              &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'Hbk', hbk(:,outloop),                &
     &                            (/1,outloop/), (/ndatum(ng),1/),      &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Lanczos vectors.
!
      IF (innloop.gt.0) THEN
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'zcglwk', zcglwk(:,innloop,outloop),  &
     &                            (/1,innloop,outloop/),                &
     &                            (/ndatum(ng)+1,1,1/),                 &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,              &
     &                            'vcglwk', vcglwk(:,innloop,outloop),  &
     &                            (/1,innloop,outloop/),                &
     &                            (/ndatum(ng)+1,1,1/),                 &
     &                            piofile = dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Saved TLmodVal_S.
!
      CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name,                &
     &                          'TLmodVal_S',                           &
     &                          tlmodval_s(:,innloop,outloop),          &
     &                          (/1,innloop,outloop/),                  &
     &                          (/ndatum(ng),1,1/),                     &
     &                          piofile = dav(ng)%pioFile)
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!-----------------------------------------------------------------------
!  Synchronize model/observation NetCDF file to disk.
!-----------------------------------------------------------------------
!
      CALL pio_netcdf_sync (ng, model, dav(ng)%name, dav(ng)%pioFile)
!
      RETURN
      END SUBROUTINE cg_write_rpcg_pio
# endif
!
!***********************************************************************
      SUBROUTINE cg_read_rpcg (ng, model, outLoop)
!***********************************************************************
!                                                                      !
!  This routine reads conjugate gradient variables for previous outer  !
!  loops using either the standard NetCDF library or the Parallel-IO   !
!  (PIO) library.                                                      !
!                                                                      !
!=======================================================================
!
!  Imported variable declarations
!
      integer, intent(in) :: ng, model, outloop
!
!  Local variable declarations.
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", cg_read_rpcg"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Read in conjugate gradient variables.
!-----------------------------------------------------------------------
!
      SELECT CASE (dav(ng)%IOtype)
        CASE (io_nf90)
          CALL cg_read_rpcg_nf90 (ng, model, outloop)
 
# if defined PIO_LIB && defined DISTRIBUTE
        CASE (io_pio)
          CALL cg_read_rpcg_pio (ng, model, outloop)
# endif
        CASE DEFAULT
          IF (master) WRITE (stdout,10) dav(ng)%IOtype
          exit_flag=3
      END SELECT
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
  10  FORMAT (' CG_READ_RPCG - Illegal output type, io_type = ',i0)
!
      RETURN
      END SUBROUTINE cg_read_rpcg
!
!***********************************************************************
      SUBROUTINE cg_read_rpcg_nf90 (ng, model, outLoop)
!***********************************************************************
!
      USE mod_netcdf
!
!  Imported variable declarations
!
      integer, intent(in) :: ng, model, outloop
!
!  Local variable declarations.
!
      integer :: status
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", cg_read_rpcg_nf90"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  If split 4D-Var and outer>1, Read in conjugate gradient variables
!  four outerloop restart.
!-----------------------------------------------------------------------
!
!  Open DAV NetCDF file for reading.
!
      IF (dav(ng)%ncid.eq.-1) THEN
        CALL netcdf_open (ng, model, dav(ng)%name, 1, dav(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Read in number of converget Ritz eigenvalues
!
      CALL netcdf_get_ivar (ng, model, dav(ng)%name,                    &
     &                      'nConvRitz', nconvritz,                     &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/outloop-1/),                      &
     &                      total = (/1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in Ritz eigenvalues.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'Ritz', ritz,                               &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,outloop-1/),                    &
     &                      total = (/ninner,1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in conjugate gradient "beta" coefficients.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'cg_beta', cg_beta,                         &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1/),                            &
     &                      total = (/ninner+1,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in conjugate gradient "delta" coefficients.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'cg_delta', cg_delta,                       &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1/),                            &
     &                      total = (/ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in normalization coefficients for Lanczos vectors.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'cg_dla', cg_dla,                           &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1/),                            &
     &                      total = (/ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in initial gradient normalization factor.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'cg_Gnorm_y', cg_gnorm_y,                   &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1/),                              &
     &                      total = (/outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'cg_Gnorm_v', cg_gnorm_v,                   &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1/),                              &
     &                      total = (/outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in Lanczos vector normalization factor.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'cg_QG', cg_qg,                             &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1/),                            &
     &                      total = (/ninner+1,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in eigenvalues of Lanczos recurrence relationship.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'cg_Ritz', cg_ritz,                         &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1/),                            &
     &                      total = (/ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read eigenvalues relative error.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'cg_RitzErr', cg_ritzerr,                   &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1/),                            &
     &                      total = (/ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in eigenvectors of Lanczos recurrence relationship.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'cg_zv', cg_zv,                             &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1,1/),                          &
     &                      total = (/ninner,ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in initial gradient for minimization.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'zgrad0', zgrad0,                           &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1/),                            &
     &                      total = (/ndatum(ng)+1,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'vgrad0', vgrad0,                           &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1/),                              &
     &                      total = (/ndatum(ng)+1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'Hbk', hbk,                                 &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1/),                            &
     &                      total = (/ndatum(ng),outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'Jb0', jb0(0:),                             &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1/),                              &
     &                      total = (/nouter+1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in Lanczos vectors.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'zcglwk', zcglwk,                           &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1,1/),                          &
     &                      total = (/ndatum(ng)+1,ninner+1,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'vcglwk', vcglwk,                           &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1,1/),                          &
     &                      total = (/ndatum(ng)+1,ninner+1,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in converged Lanczos eigenvectors.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'vcglev', vcglev,                           &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1,1/),                          &
     &                      total = (/ndatum(ng)+1,ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in TLmodVal_S.
!
      CALL netcdf_get_fvar (ng, model, dav(ng)%name,                    &
     &                      'TLmodVal_S', tlmodval_s,                   &
     &                      ncid = dav(ng)%ncid,                        &
     &                      start = (/1,1,1/),                          &
     &                      total = (/ndatum(ng),ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      RETURN
      END SUBROUTINE cg_read_rpcg_nf90
 
# if defined PIO_LIB && defined DISTRIBUTE
!
!***********************************************************************
      SUBROUTINE cg_read_rpcg_pio (ng, model, outLoop)
!***********************************************************************
!
      USE mod_pio_netcdf
!
!  Imported variable declarations
!
      integer, intent(in) :: ng, model, outloop
!
!  Local variable declarations.
!
      integer :: status
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//", cg_read_rpcg_pio"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  If split 4D-Var and outer>1, Read in conjugate gradient variables
!  four outerloop restart.
!-----------------------------------------------------------------------
!
!  Open DAV NetCDF file for reading.
!
      IF (dav(ng)%pioFile%fh.eq.-1) THEN
        CALL pio_netcdf_open (ng, model, dav(ng)%name, 1,               &
     &                        dav(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END IF
!
!  Read in number of converget Ritz eigenvalues
!
      CALL pio_netcdf_get_ivar (ng, model, dav(ng)%name,                &
     &                          'nConvRitz', nconvritz,                 &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/outloop-1/),                  &
     &                          total = (/1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in Ritz eigenvalues.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'Ritz', ritz,                           &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,outloop-1/),                &
     &                          total = (/ninner,1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in conjugate gradient "beta" coefficients.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'cg_beta', cg_beta,                     &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1/),                        &
     &                          total = (/ninner+1,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in conjugate gradient "delta" coefficients.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'cg_delta', cg_delta,                   &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1/),                        &
     &                          total = (/ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in normalization coefficients for Lanczos vectors.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'cg_dla', cg_dla,                       &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1/),                        &
     &                          total = (/ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in initial gradient normalization factor.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'cg_Gnorm_y', cg_gnorm_y,               &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1/),                          &
     &                          total = (/outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'cg_Gnorm_v', cg_gnorm_v,               &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1/),                          &
     &                          total = (/outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in Lanczos vector normalization factor.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'cg_QG', cg_qg,                         &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1/),                        &
     &                          total = (/ninner+1,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in eigenvalues of Lanczos recurrence relationship.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'cg_Ritz', cg_ritz,                     &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1/),                        &
     &                          total = (/ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read eigenvalues relative error.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'cg_RitzErr', cg_ritzerr,               &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1/),                        &
     &                          total = (/ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in eigenvectors of Lanczos recurrence relationship.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'cg_zv', cg_zv,                         &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1,1/),                      &
     &                          total = (/ninner,ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in initial gradient for minimization.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'zgrad0', zgrad0,                       &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1/),                        &
     &                          total = (/ndatum(ng)+1,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'vgrad0', vgrad0,                       &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1/),                          &
     &                          total = (/ndatum(ng)+1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'Hbk', hbk,                             &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1/),                        &
     &                          total = (/ndatum(ng),outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'Jb0', jb0(0:),                         &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1/),                          &
     &                          total = (/nouter+1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in Lanczos vectors.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'zcglwk', zcglwk,                       &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1,1/),                      &
     &                          total = (/ndatum(ng)+1,ninner+1,        &
     &                                    outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'vcglwk', vcglwk,                       &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1,1/),                      &
     &                          total = (/ndatum(ng)+1,ninner+1,        &
     &                                   outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in converged Lanczos eigenvectors.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'vcglev', vcglev,                       &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1,1/),                      &
     &                          total = (/ndatum(ng)+1,ninner,          &
     &                                    outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Read in TLmodVal_S.
!
      CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name,                &
     &                          'TLmodVal_S', tlmodval_s,               &
     &                          piofile = dav(ng)%pioFile,              &
     &                          start = (/1,1,1/),                      &
     &                          total = (/ndatum(ng),ninner,outloop-1/))
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
      RETURN
      END SUBROUTINE cg_read_rpcg_pio
# endif
#endif
      END MODULE rpcg_lanczos_mod