convolve_mod Module Reference

Functions/Subroutines
subroutine, public	convolve (driver, rini, rold, rnew)
subroutine, public	error_covariance (model, driver, outloop, innloop, rbck, rini, rold, rnew, rec1, rec2, lposterior)
subroutine, public	saddlec (driver, lselect, rini, rold, rnew)

Function/Subroutine Documentation

convolve()

subroutine, public convolve_mod::convolve	(	character (len=*), intent(in)	driver,
		integer, intent(in)	rini,
		integer, dimension(ngrids), intent(in)	rold,
		integer, dimension(ngrids), intent(in)	rnew )

Definition at line 115 of file convolve.F.

!***********************************************************************
!
!  Imported variable declarations.
!
      integer, intent(in) :: Rini
      integer, intent(in) :: Rold(Ngrids)
      integer, intent(in) :: Rnew(Ngrids)
!
      character (len=*), intent(in) :: driver
!
!  Local variable declarations.
!
      logical :: Lweak, add
!
      integer :: IniRec
      integer :: ng, tile
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Specify spatial error covariance on the TLM control vector, at time
!  index "Rold", via convolutions of a pseudo-diffusion operator.  The
!  convolved control vector is loaded into time index "Rnew" and index
!  "Rold" is preserved.
!-----------------------------------------------------------------------
 
# ifdef BALANCE_OPERATOR
!
!  Read NLM initial condition in readiness for the balance operator.
!
      inirec=rini
      DO ng=1,ngrids
        CALL get_state (ng, inlm, 2, ini(ng), inirec, inirec)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        nrhs(ng)=rini
      END DO
# endif
!
!  Load "Rold" time index TLM control vector into ADM state arrays.
!  Apply balance operator, if activated. Then, scale control vector
!  with error covariance standard deviations. Next, convolve resulting
!  vector with the squared-root adjoint diffusion operator. Notice
!  that the spatial convolution is only done for half of the diffusion
!  steps (squared-root filter).
!
      lweak=.false.
      add=.false.
      DO ng=1,ngrids
# ifdef PROFILE
        CALL wclock_on (ng, iadm, 82, __line__, myfile)
# endif
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
# ifdef BALANCE_OPERATOR
          CALL ad_balance (ng, tile, rini, rold(ng))
# endif
          CALL ad_variability (ng, tile, rold(ng), lweak)
          CALL ad_convolution (ng, tile, rold(ng), lweak, 2)
        END DO
# ifdef PROFILE
        CALL wclock_off (ng, iadm, 82, __line__, myfile)
# endif
      END DO
!
!  Since we wish to preserve what is in tl_var(Rold), load resulting
!  filtered solution from above, ad_var(Rold), into tl_var(Rnew).
!  Then, convolve with the squared-root (half of steps) tangent
!  linear diffusion operator. Next, scale results with error
!  covariance standard deviations. Apply balance operator, if
!  activated.
!
      add=.false.
      DO ng=1,ngrids
# ifdef PROFILE
        CALL wclock_on (ng, itlm, 82, __line__, myfile)
# endif
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL load_adtotl (ng, tile, rold(ng), rnew(ng), add)
          CALL tl_convolution (ng, tile, rnew(ng), lweak, 2)
          CALL tl_variability (ng, tile, rnew(ng), lweak)
# ifdef BALANCE_OPERATOR
          CALL tl_balance (ng, tile, rini, rnew(ng))
# endif
        END DO
# ifdef PROFILE
        CALL wclock_off (ng, itlm, 82, __line__, myfile)
# endif
      END DO
!
      RETURN

References ad_balance_mod::ad_balance(), ad_convolution_mod::ad_convolution(), ad_variability_mod::ad_variability(), mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), get_state_mod::get_state(), mod_param::iadm, mod_iounits::ini, mod_param::inlm, mod_param::itlm, mod_parallel::last_tile, ini_adjust_mod::load_adtotl(), ini_adjust_mod::load_tltoad(), mod_param::ngrids, mod_scalars::noerror, mod_stepping::nrhs, mod_iounits::sourcefile, tl_balance_mod::tl_balance(), tl_convolution_mod::tl_convolution(), tl_variability_mod::tl_variability(), wclock_off(), and wclock_on().

Referenced by r4dvar_mod::posterior_error(), and rbl4dvar_mod::posterior_error().

Here is the call graph for this function:

Here is the caller graph for this function:

error_covariance()

subroutine, public convolve_mod::error_covariance	(	integer, intent(in)	model,
		character (len=*), intent(in)	driver,
		integer, intent(in)	outloop,
		integer, intent(in)	innloop,
		integer, intent(in)	rbck,
		integer, intent(in)	rini,
		integer, dimension(ngrids), intent(in)	rold,
		integer, dimension(ngrids), intent(in)	rnew,
		integer, intent(in)	rec1,
		integer, intent(in)	rec2,
		logical, intent(in)	lposterior )

Definition at line 212 of file convolve.F.

!***********************************************************************
!
!  Imported variable declarations.
!
      logical, intent(in) :: Lposterior
!
      integer, intent(in) :: model, outLoop, innLoop
      integer, intent(in) :: Rbck, Rini, Rec1, Rec2
      integer, intent(in) :: Rold(Ngrids)
      integer, intent(in) :: Rnew(Ngrids)
!
      character (len=*), intent(in) :: driver
!
!  Local variable declarations.
!
      logical :: Lweak, add
!
      integer :: ADrec, BckRec, Fcount, IniRec
      integer :: i, irec, ng, tile
# ifdef RPCG
      integer :: LTLM1, LTLM2, Rec5, jrec, nADrec
# endif
      integer, dimension(Ngrids) :: Nrec
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", error_covariance"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Convolve adjoint trajectory.
!-----------------------------------------------------------------------
!
      IF (master) THEN
        IF ((outloop.lt.0).and.(innloop.lt.0)) THEN
          WRITE (stdout,10)
 10       FORMAT (/,' Convolving Adjoint Trajectory...')
        ELSE
          WRITE (stdout,20) outloop, innloop
 20       FORMAT (/,' Convolving Adjoint Trajectory: Outer = ',i3.3,    &
     &            ' Inner = ',i3.3)
        END IF
      END IF
!
!  Clear adjoint state arrays.
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL initialize_ocean (ng, tile, iadm)
        END DO
      END DO
!
!  Get number of records in adjoint NetCDF and initialize indices.
!
      DO ng=1,ngrids
        fcount=adm(ng)%Fcount
        nrec(ng)=adm(ng)%Nrec(fcount)
        adm(ng)%Nrec(fcount)=0
        adm(ng)%Rindex=0
        lwrtstate2d(ng)=.true.
        lwrttime(ng)=.false.
      END DO
!
!  Read in adjoint control vector to convolve from NetCDF file
!  ADM(ng)%name for time record Nrec, which correspond to the
!  initial conditions time. If adjusting open boundaries and/or
!  surface forcing, only record Nrec is read since it is the
!  only record for which adjoint forcing arrays are complete.
!
!  Notice that since routine "get_state" loads data into the
!  ghost points, the adjoint contol vector is read into the
!  tangent linear state arrays by using iTLM instead of iADM
!  in the calling arguments.
!
      DO ng=1,ngrids
        adrec=nrec(ng)
        frcrec(ng)=nrec(ng)
        CALL get_state (ng, itlm, 4, adm(ng), adrec, rold(ng))
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
 
# ifdef BALANCE_OPERATOR
!
!  Read NLM initial condition in readiness for the balance
!  operator.
!
      inirec=rini
      DO ng=1,ngrids
        CALL get_state (ng, inlm, 2, ini(ng), inirec, inirec)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        nrhs(ng)=rini
      END DO
# endif
!
!  Load ADM control read above from TLM into ADM state arrays.
!  Then, multiply control vector by the error covariance standard
!  deviations. Next, convolve resulting adjoint solution with the
!  squared-root adjoint diffusion operator to impose apriori error
!  hypothesis. Notice that the spatial convolution are only done
!  for half of the diffusion steps (squared-root filter). Clear
!  tangent linear state arrays when done.
!
      lweak=.false.
      add=.false.
      DO ng=1,ngrids
# ifdef PROFILE
        CALL wclock_on (ng, model, 82, __line__, myfile)
# endif
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
# ifdef BALANCE_OPERATOR
          CALL ad_balance (ng, tile, rini, rold(ng))
# endif
          CALL ad_variability (ng, tile, rold(ng), lweak)
          CALL ad_convolution (ng, tile, rold(ng), lweak, 2)
          CALL initialize_ocean (ng, tile, itlm)
          CALL initialize_forces (ng, tile, itlm)
# ifdef ADJUST_BOUNDARY
          CALL initialize_boundary (ng, tile, itlm)
# endif
        END DO
# ifdef PROFILE
        CALL wclock_off (ng, model, 82, __line__, myfile)
# endif
      END DO
!
!  To insure symmetry, convolve resulting filtered adjoint solution
!  from above with the squared-root (half of steps) tangent linear
!  diffusion operator. Then, multiply result with its corresponding
!  error covariance standard deviations. Since the convolved solution
!  is in the adjoint state arrays, first copy to tangent linear state
!  arrays including the ghosts points.
!
      lweak=.false.
      add=.false.
      DO ng=1,ngrids
# ifdef PROFILE
        CALL wclock_on (ng, model, 82, __line__, myfile)
# endif
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL load_adtotl (ng, tile, rold(ng), rold(ng), add)
          CALL tl_convolution (ng, tile, rold(ng), lweak, 2)
          CALL tl_variability (ng, tile, rold(ng), lweak)
# ifdef BALANCE_OPERATOR
          CALL tl_balance (ng, tile, rini, rold(ng))
# endif
        END DO
 
# ifdef POSTERIOR_ERROR_I
!
!  If computing the analysis error covariance matrix, copy TLM back
!  into ADM so that it can be written to Hessian NetCDF file.
!
        IF (lposterior) THEN
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
          END DO
        END IF
# endif
# ifdef PROFILE
        CALL wclock_off (ng, model, 82, __line__, myfile)
# endif
      END DO
 
# if defined ARRAY_MODES            || defined R4DVAR || \
     defined r4dvar_ana_sensitivity
!
!  Copy back to adjoint state arrays when done with the convolution.
!  Compute representer model initial conditions by adding convolved
!  adjoint solution to the reference nonlinear state (INI(ng)%name,
!  record Rbck).
!
      IF ((model.eq.irpm).and.                                          &
     &    (index(driver,'r4dvar').ne.0)) THEN
        bckrec=rbck
        DO ng=1,ngrids
          CALL get_state (ng, inlm, 9, ini(ng), bckrec, rnew(ng))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
        add=.false.
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
            CALL rp_ini_adjust (ng, tile, rnew(ng), rold(ng))
          END DO
        END DO
      END IF
 
# elif defined RBL4DVAR || defined RBL4DVAR_ANA_SENSITIVITY
#  ifndef RPCG
!
!  Copy back to adjoint state arrays when done with the convolution.
!  Compute nonlinear model initial conditions by adding convolved
!  adjoint solution to the reference nonlinear state (INI(ng)%name,
!  record Rbck).
!
      IF ((model.eq.inlm).and.                                          &
     &    (index(driver,'rbl4dvar').ne.0)) THEN
        bckrec=rbck
        DO ng=1,ngrids
          CALL get_state (ng, inlm, 9, ini(ng), bckrec, rnew(ng))
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
!
        add=.false.
        DO ng=1,ngrids
          DO tile=first_tile(ng),last_tile(ng),+1
            CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
            CALL ini_adjust (ng, tile, rold(ng), rnew(ng))
          END DO
        END DO
      END IF
#  endif
# endif
!
!  Write out tangent linear model initial conditions and tangent
!  linear surface forcing adjustments for next inner loop into
!  ITL(ng)%name (record Rec1). The tangent model initial
!  conditions are set to the convolved adjoint solution.
!
      IF (model.eq.itlm) THEN
        DO ng=1,ngrids
# ifdef DISTRIBUTE
          CALL tl_wrt_ini (ng, myrank, rold(ng), rec1)
# else
          CALL tl_wrt_ini (ng, -1, rold(ng), rec1)
# endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
# if defined ARRAY_MODES            || defined R4DVAR || \
     defined r4dvar_ana_sensitivity
      ELSE IF (model.eq.irpm) THEN
        DO ng=1,ngrids
#  ifdef DISTRIBUTE
          CALL rp_wrt_ini (ng, myrank, rold(ng), rec2)
#  else
          CALL rp_wrt_ini (ng, -1, rold(ng), rec2)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
# elif (defined RBL4DVAR || defined RBL4DVAR_ANA_SENSITIVITY) && \
       !defined RPCG
      ELSE IF (model.eq.inlm) THEN
        DO ng=1,ngrids
#  ifdef DISTRIBUTE
          CALL wrt_ini (ng, myrank, rnew(ng))
#  else
          CALL wrt_ini (ng, -1, rnew(ng))
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
      defined adjust_boundary
#   ifdef DISTRIBUTE
          CALL wrt_frc_ad (ng, myrank, rold(ng), ini(ng)%Rindex)
#   else
          CALL wrt_frc_ad (ng, -1, rold(ng), ini(ng)%Rindex)
#   endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
#  endif
        END DO
# endif
      END IF
 
# ifdef RPCG
!
!  Write out tangent linear model initial conditions and tangent
!  linear surface forcing and obc adjustments for next outer
!  loop into ITLname (record Rec2). The tangent model initial
!  conditions are set to the convolved adjoint solution.
!
      IF (model.eq.inlm) THEN
        DO ng=1,ngrids
#  ifdef DISTRIBUTE
          CALL tl_wrt_ini (ng, myrank, rold(ng), rec2)
#  else
          CALL tl_wrt_ini (ng, -1, rold(ng), rec2)
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END IF
# endif
 
# ifdef POSTERIOR_ERROR_I
!
!  If weak constraint, convolve records 2:Nrec in ADM(ng)%name and
!  Write convolved adjoint solution into Hessian NetCDF file for use
!  later.
!
      IF (lposterior.and.(inner.ne.0)) THEN
        DO ng=1,ngrids
#  ifdef DISTRIBUTE
          CALL wrt_hessian (ng, myrank, rold(ng), rold(ng))
#  else
          CALL wrt_hessian (ng, -1, rold(ng), rold(ng))
#  endif
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END IF
# endif
# ifdef RPCG
!
!  Now compute the new B^-1(x(k)-x(k-1)) term for the weak constraint
!  forcing terms. Records 6+nADrec/2 to nADrec+5 contain the sums so
!  far. This is ONLY done in the outer-loop not the inner-loops, and
!  is controlled by the flag "LaugWeak".
!
      IF (laugweak) THEN
        ltlm1=1
        ltlm2=2
        rec5=5
        DO ng=1,ngrids
          IF (nadj(ng).lt.ntimes(ng)) THEN
            nadrec=2*(1+ntimes(ng)/nadj(ng))
          ELSE
            nadrec=0
          END IF
          DO irec=1,nadrec/2
            jrec=rec5+nadrec/2+irec
!
!  First add the augmented piece which is computed from the previous
!  sum.
!
            CALL get_state (ng, itlm, 4, itl(ng), jrec, ltlm1)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL aug_oper (ng, tile, ltlm1, ltlm2)
            END DO
!
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL sum_grad (ng, tile, ltlm1, ltlm2)
            END DO
!
!  Need to read adjoint netcdf file in reverse.
!
            adrec=(nrec(ng)-1)-(irec-1)
            CALL get_state (ng, itlm, 4, adm(ng), adrec, ltlm1)
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL sum_grad (ng, tile, ltlm1, ltlm2)
            END DO
!
!  Write the current sum of adjoint solutions into record jrec of the
!  ITL file.
!
#  ifdef DISTRIBUTE
            CALL tl_wrt_ini (ng, myrank, ltlm2, jrec)
#  else
            CALL tl_wrt_ini (ng, -1, ltlm2, jrec)
#  endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
          END DO
        END DO
      END IF
# endif
# ifdef TIME_CONV
!
!  Apply the factorized space-time model error covariance
!  matrix of the form ATA' where A is the square-root of the
!  model error covariance matrix, and T is the time-correlation
!  matrix.
!
!  If weak constraint, convolve records 2-Nrec in ADM(ng)%name
!  and impose model error covariance. NOTE: We will not use the
!  convolved forcing increments generated here since these arrays
!  do not contain the complete solution and are redundant.
!  AMM: We might want to get rid of these unwanted records to
!  avoid any confusion in the future.
!
      DO ng=1,ngrids
        IF (nrec(ng).gt.3) THEN
#  ifdef PROFILE
          CALL wclock_on (ng, model, 82, __line__, myfile)
#  endif
          DO irec=1,nrec(ng)-1
!
!  Read adjoint solution. Since routine "get_state" loads data into the
!  ghost points, the adjoint solution is read in the tangent linear
!  state arrays by using iTLM instead of iADM in the calling arguments.
!
            adrec=irec
            CALL get_state (ng, itlm, 4, adm(ng), adrec, rold(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Load interior solution, read above, into adjoint state arrays.
!  Then, multiply adjoint solution by the background-error standard
!  deviations. Next, convolve resulting adjoint solution with the
!  squared-root adjoint diffusion operator which impose the model-error
!  spatial correlations. Notice that the spatial convolution is only
!  done for half of the diffusion steps (squared-root filter). Clear
!  tangent linear state arrays when done.
!
            lweak=.true.
            add=.false.
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
#  ifdef BALANCE_OPERATOR
              CALL ad_balance (ng, tile, rini, rold(ng))
#  endif
              CALL ad_variability (ng, tile, rold(ng), lweak)
              CALL ad_convolution (ng, tile, rold(ng), lweak, 2)
              CALL initialize_ocean (ng, tile, itlm)
              CALL initialize_forces (ng, tile, itlm)
#  ifdef ADJUST_BOUNDARY
              CALL initialize_boundary (ng, tile, itlm)
#  endif
            END DO
!
!  Overwrite ADM(ng)%name history NetCDF file with convolved adjoint
!  solution.
!
            kstp(ng)=rold(ng)
#  ifdef SOLVE3D
            nstp(ng)=rold(ng)
#  endif
#  ifdef DISTRIBUTE
            CALL ad_wrt_his (ng, myrank)
#  else
            CALL ad_wrt_his (ng, -1)
#  endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
          lwrtstate2d(ng)=.false.
          lwrttime(ng)=.true.
#  ifdef PROFILE
        CALL wclock_off (ng, model, 82, __line__, myfile)
#  endif
        END IF
      END DO
!
!  Apply the time correlation matrix to the newly convolved adjoint
!  fields.
!
      DO ng=1,ngrids
        tlf(ng)%Rindex=0
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL time_corr (ng, tile, iadm, adm(ng)%IOtype, adm(ng)%name)
          IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
        END DO
      END DO
!
      DO ng=1,ngrids
        IF (nrec(ng).gt.3) THEN
#  ifdef PROFILE
          CALL wclock_on (ng, model, 82, __line__, myfile)
#  endif
          adm(ng)%Rindex=0
          lwrtstate2d(ng)=.true.
          DO irec=nrec(ng)-1,1,-1
!
!  Read the latest TLF solution. We need to read the TLF file in reverse
!  order since the final solution is written into the adjoint netcdf
!  file which will be rewritten in ascending order of time in the TLF
!  file by wrt_impulse.
!  NOTE: model=6 here so as to only read the state variables.
!
            adrec=irec
            CALL get_state (ng, 6, 6, tlf(ng), adrec, rold(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
            lweak=.true.
            add=.false.
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL tl_convolution (ng, tile, rold(ng), lweak, 2)
              CALL tl_variability (ng, tile, rold(ng), lweak)
#  ifdef BALANCE_OPERATOR
              CALL tl_balance (ng, tile, rini, rold(ng))
#  endif
              CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
            END DO
!
!  Overwrite ADM(ng)%name history NetCDF file with convolved adjoint
!  solution.
!
            kstp(ng)=rold(ng)
#  ifdef SOLVE3D
            nstp(ng)=rold(ng)
#  endif
#  ifdef DISTRIBUTE
            CALL ad_wrt_his (ng, myrank)
#  else
            CALL ad_wrt_his (ng, -1)
#  endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
          lwrtstate2d(ng)=.false.
          lwrttime(ng)=.true.
#  ifdef PROFILE
          CALL wclock_off (ng, model, 82, __line__, myfile)
#  endif
        END IF
      END DO
 
# else
!
!  If weak constraint, convolve records 2-Nrec in ADM(ng)%name
!  and impose model error covariance. NOTE: We will not use the
!  convolved forcing increments generated here since these arrays
!  do not contain the complete solution and are redundant.
!  AMM: We might want to get rid of these unwanted records to
!  avoid any confusion in the future.
!
      DO ng=1,ngrids
        IF (nrec(ng).gt.3) THEN
#  ifdef PROFILE
          CALL wclock_on (ng, model, 82, __line__, myfile)
#  endif
          DO irec=1,nrec(ng)-1
!
!  Read adjoint solution. Since routine "get_state" loads data into the
!  ghost points, the adjoint solution is read in the tangent linear
!  state arrays by using iTLM instead of iADM in the calling arguments.
!
            adrec=irec
            CALL get_state (ng, itlm, 4, adm(ng), adrec, rold(ng))
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Load interior solution, read above, into adjoint state arrays.
!  Then, multiply adjoint solution by the background-error standard
!  deviations. Next, convolve resulting adjoint solution with the
!  squared-root adjoint diffusion operator which impose the model-error
!  spatial correlations. Notice that the spatial convolution is only
!  done for half of the diffusion steps (squared-root filter). Clear
!  tangent linear state arrays when done.
!
            lweak=.true.
            add=.false.
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
#  ifdef BALANCE_OPERATOR
              CALL ad_balance (ng, tile, rini, rold(ng))
#  endif
              CALL ad_variability (ng, tile, rold(ng), lweak)
              CALL ad_convolution (ng, tile, rold(ng), lweak, 2)
              CALL initialize_ocean (ng, tile, itlm)
              CALL initialize_forces (ng, tile, itlm)
#  ifdef ADJUST_BOUNDARY
              CALL initialize_boundary (ng, tile, itlm)
#  endif
            END DO
!
!  To insure symmetry, convolve resulting filtered adjoint solution
!  from above with the squared-root (half of steps) tangent linear
!  diffusion operator. Then, multiply result with its corresponding
!  background-error standard deviations.  Since the convolved solution
!  is in the adjoint state arrays, first copy to tangent linear state
!  arrays including the ghosts points. Copy back to adjoint state
!  arrays when done with the convolution for output purposes.
!
            lweak=.true.
            add=.false.
            DO tile=first_tile(ng),last_tile(ng),+1
              CALL load_adtotl (ng, tile, rold(ng), rold(ng), add)
              CALL tl_convolution (ng, tile, rold(ng), lweak, 2)
              CALL tl_variability (ng, tile, rold(ng), lweak)
#  ifdef BALANCE_OPERATOR
              CALL tl_balance (ng, tile, rini, rold(ng))
#  endif
              CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
            END DO
!
!  Overwrite ADM(ng)%name history NetCDF file with final solution.
!
            kstp(ng)=rold(ng)
#  ifdef SOLVE3D
            nstp(ng)=rold(ng)
#  endif
#  ifdef DISTRIBUTE
            CALL ad_wrt_his (ng, myrank)
#  else
            CALL ad_wrt_his (ng, -1)
#  endif
            IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
          END DO
          lwrtstate2d(ng)=.false.
          lwrttime(ng)=.true.
#  ifdef PROFILE
          CALL wclock_off (ng, model, 82, __line__, myfile)
#  endif
        END IF
      END DO
# endif
!
      RETURN

References ad_balance_mod::ad_balance(), ad_convolution_mod::ad_convolution(), ad_variability_mod::ad_variability(), ad_wrt_his_mod::ad_wrt_his(), mod_iounits::adm, comp_jb0_mod::aug_oper(), mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), mod_scalars::frcrec, get_state_mod::get_state(), mod_param::iadm, mod_iounits::ini, ini_adjust_mod::ini_adjust(), mod_boundary::initialize_boundary(), mod_forces::initialize_forces(), mod_ocean::initialize_ocean(), mod_param::inlm, mod_scalars::inner, mod_param::irpm, mod_iounits::itl, mod_param::itlm, mod_stepping::kstp, mod_parallel::last_tile, mod_fourdvar::laugweak, ini_adjust_mod::load_adtotl(), ini_adjust_mod::load_tltoad(), mod_scalars::lwrtstate2d, mod_scalars::lwrttime, mod_parallel::master, mod_parallel::myrank, mod_scalars::nadj, mod_param::ngrids, mod_scalars::noerror, mod_stepping::nrhs, mod_stepping::nstp, mod_scalars::ntimes, ini_adjust_mod::rp_ini_adjust(), rp_wrt_ini_mod::rp_wrt_ini(), mod_iounits::sourcefile, mod_iounits::stdout, sum_grad_mod::sum_grad(), time_corr_mod::time_corr(), tl_balance_mod::tl_balance(), tl_convolution_mod::tl_convolution(), tl_variability_mod::tl_variability(), tl_wrt_ini_mod::tl_wrt_ini(), mod_iounits::tlf, wclock_off(), wclock_on(), wrt_ini_mod::wrt_frc_ad(), wrt_hessian_mod::wrt_hessian(), and wrt_ini_mod::wrt_ini().

Referenced by r4dvar_mod::increment(), rbl4dvar_mod::increment(), and roms_kernel_mod::roms_run().

Here is the call graph for this function:

Here is the caller graph for this function:

saddlec()

subroutine, public convolve_mod::saddlec	(	character (len=*), intent(in)	driver,
		logical, intent(in)	lselect,
		integer, intent(in)	rini,
		integer, dimension(ngrids), intent(in)	rold,
		integer, dimension(ngrids), intent(in)	rnew )

Definition at line 807 of file convolve.F.

!***********************************************************************
!
!  Imported variable declarations.
!
      integer, intent(in) :: Rini
      integer, intent(in) :: Rold(Ngrids)
      integer, intent(in) :: Rnew(Ngrids)
!
      logical, intent(in) :: Lselect
!
      character (len=*), intent(in) :: driver
!
!  Local variable declarations.
!
      logical :: Lweak, add
!
      integer :: IniRec
      integer :: ng, tile
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__//", saddlec"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Specify spatial error covariance on the TLM control vector, at time
!  index "Rold", via convolutions of a pseudo-diffusion operator.  The
!  convolved control vector is loaded into time index "Rnew" and index
!  "Rold" is preserved.
!-----------------------------------------------------------------------
 
# ifdef BALANCE_OPERATOR
!
!  Read NLM initial condition in readiness for the balance operator.
!
      inirec=rini
      DO ng=1,ngrids
        CALL get_state (ng, inlm, 2, ini(ng), inirec, inirec)
        IF (exit_flag.ne.noerror) RETURN
        nrhs(ng)=rini
      END DO
# endif
!
!  Load "Rold" time index TLM control vector into ADM state arrays.
!  Apply balance operator, if activated. Then, scale control vector
!  with error covariance standard deviations. Next, convolve resulting
!  vector with the squared-root adjoint diffusion operator. Notice
!  that the spatial convolution is only done for half of the diffusion
!  steps (squared-root filter).
!
      lweak=lselect
      add=.false.
      DO ng=1,ngrids
# ifdef PROFILE
        CALL wclock_on (ng, iadm, 82, __line__, myfile)
# endif
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL load_tltoad (ng, tile, rold(ng), rold(ng), add)
# ifdef BALANCE_OPERATOR
          CALL ad_balance (ng, tile, rini, rold(ng))
# endif
          CALL ad_variability (ng, tile, rold(ng), lweak)
          CALL ad_convolution (ng, tile, rold(ng), lweak, 2)
! AMM Important.
          CALL initialize_ocean (ng, tile, itlm)
          CALL initialize_forces (ng, tile, itlm)
        END DO
# ifdef PROFILE
        CALL wclock_off (ng, iadm, 82, __line__, myfile)
# endif
      END DO
!
!  Since we wish to preserve what is in tl_var(Rold), load resulting
!  filtered solution from above, ad_var(Rold), into tl_var(Rnew).
!  Then, convolve with the squared-root (half of steps) tangent
!  linear diffusion operator. Next, scale results with error
!  covariance standard deviations. Apply balance operator, if
!  activated.
!
      add=.false.
      DO ng=1,ngrids
# ifdef PROFILE
        CALL wclock_on (ng, itlm, 82, __line__, myfile)
# endif
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL load_adtotl (ng, tile, rold(ng), rnew(ng), add)
          CALL tl_convolution (ng, tile, rnew(ng), lweak, 2)
          CALL tl_variability (ng, tile, rnew(ng), lweak)
# ifdef BALANCE_OPERATOR
          CALL tl_balance (ng, tile, rini, rnew(ng))
# endif
        END DO
# ifdef PROFILE
        CALL wclock_off (ng, itlm, 82, __line__, myfile)
# endif
      END DO
!
      RETURN

References ad_balance_mod::ad_balance(), ad_convolution_mod::ad_convolution(), ad_variability_mod::ad_variability(), mod_scalars::exit_flag, mod_parallel::first_tile, get_state_mod::get_state(), mod_param::iadm, mod_iounits::ini, mod_forces::initialize_forces(), mod_ocean::initialize_ocean(), mod_param::inlm, mod_param::itlm, mod_parallel::last_tile, ini_adjust_mod::load_adtotl(), ini_adjust_mod::load_tltoad(), mod_param::ngrids, mod_scalars::noerror, mod_stepping::nrhs, mod_iounits::sourcefile, tl_balance_mod::tl_balance(), tl_convolution_mod::tl_convolution(), tl_variability_mod::tl_variability(), wclock_off(), and wclock_on().

Here is the call graph for this function: