ROMS
Loading...
Searching...
No Matches
ad_rpcg_lanczos.F
Go to the documentation of this file.
1#include "cppdefs.h"
2
3#if defined RPCG && \
4 (defined sensitivity_4dvar || \
5 defined tl_rbl4dvar || \
6 defined tl_r4dvar)
7
8 SUBROUTINE ad_rpcg_lanczos (ng, model, outLoop, innLoop, &
9 & NinnLoop, Lcgini)
10!
11!git $Id$
12!================================================== Hernan G. Arango ===
13! Copyright (c) 2002-2025 The ROMS Group Selime Gurol !
14! Licensed under a MIT/X style license Andrew M. Moore !
15! See License_ROMS.md !
16!=======================================================================
17! !
18! Weak Constraint 4-Dimensional Variational (4DVar) Restricted !
19! Pre-conditioned Conjugate Gradient (RPCG) Algorithm !
20
21# if defined R4DVAR_ANA_SENSITIVITY || defined TL_R4DVAR
22! !
23! The indirect representer method solves the system: !
24! !
25! (R_n + Cobs) * Beta_n = h_n !
26! !
27! h_n = Xo - H * X_n !
28! !
29! where R_n is the representer matrix, Cobs is the observation-error !
30! covariance, Beta_n are the representer coefficients, h_n is the !
31! misfit between observations (Xo) and model (H*X_n), and H is the !
32! linearized observation operator. Here, _n denotes iteration. !
33! !
34! This system does not need to be solved explicitly by inverting the !
35! symmetric stabilized representer matrix, P_n: !
36! !
37! P_n = R_n + Cobs !
38! !
39! but by computing the action of P_n on any vector PSI, such that !
40! !
41! P_n * PSI = R_n * PSI + Cobs * PSI !
42! !
43! The representer matrix is not explicitly computed but evaluated by !
44! one integration backward of the adjoint model and one integration !
45! forward of the tangent linear model for any forcing vector PSI. !
46! !
47! Notice that "ObsScale" vector is used for screenning observations. !
48! This scale is one (zero) for good (bad) observations. !
49! !
50! Currently, parallelization of this algorithm is not needed because !
51! each parallel node has a full copy of the assimilation vectors. !
52! !
53! This code solves Ax=b by minimizing the cost function !
54! 0.5*x*A*x-x*b, assuming an initial guess of x=x0. In this case the !
55! gradient is Ax-b and the stabilized representer matrix is A. !
56! !
57! Reference: !
58! !
59! Chua, B. S. and A. F. Bennett, 2001: An inverse ocean modeling !
60! sytem, Ocean Modelling, 3, 137-165. !
61
62# elif defined RBL4DVAR_ANA_SENSITIVITY || \
63 defined rbl4dvar_fct_sensitivity || \
64 defined tl_rbl4dvar
65! !
66! Solve the system (following Courtier, 1997): !
67! !
68! (H M_n B (M_n)' H' + Cobs) * w_n = d_n !
69! !
70! d_n = yo - H * Xb_n !
71! !
72! where M_n is the tangent linear model matrix and _n denotes a !
73! sequence of outer-loop estimates, Cobs is the observation-error !
74! covariance, B is the background error covariance, d_n is the !
75! misfit between observations (yo) and model (H * Xb_n), and H is !
76! the linearized observation operator. !
77! !
78! The analysis increment is: !
79! !
80! dx_n = B M' H' w_n !
81! !
82! so that Xa = Xb + dx_n. !
83! !
84! The system does not need to be solved explicitly by inverting !
85! the symmetric matrix, P_n: !
86! !
87! P_n = H M_n B (M_n)' H' + Cobs !
88! !
89! but by computing the action of P_n on any vector PSI, such that !
90! !
91! P_n * PSI = H M_n B (M_n)' H' * PSI + Cobs * PSI !
92! !
93! The (H M_n B (M_n)' H') matrix is not explicitly computed but !
94! evaluated by one integration backward of the adjoint model and !
95! one integration forward of the tangent linear model for any !
96! forcing vector PSI. !
97! !
98! Reference: !
99! !
100! Courtier, P., 1997: Dual formulation of four-dimensional !
101! variational assimilation, Quart. J. Roy. Meteor. Soc., !
102! 123, 2449-2461. !
103# endif
104! !
105! The restricted preconditioned conjugate gradient (RPCG) algorithm !
106! is used to compute the solution. Specfically, the following system !
107! is solved: !
108! !
109! (Cobs^-1 H M_n B (M_n)' H'+I)*w=Cobs^-1 d_n !
110! !
111! All inner-products are defined with respect to a norm involving !
112! H M_n B (M_n)' H', and the convergence of J in model space is !
113! guaranteed to converge at the same rate as I4DVAR in primal space. !
114! !
115! RPCG Reference: !
116! !
117! Gratton, S. and J. Tshimanga, 2009: An observation-space !
118! formulation of variational assimilation using a restricted !
119! preconditioned conjugate gradient algorithm. QJRMS, 135, !
120! 1573-1585.
121! !
122!=======================================================================
123!
124 USE mod_param
125 USE mod_parallel
126 USE mod_fourdvar
127 USE mod_iounits
128 USE mod_scalars
129
130# ifdef DISTRIBUTE
131!
132 USE distribute_mod, ONLY : mp_bcastf, mp_bcasti, mp_bcastl
133# endif
134 implicit none
135!
136! Imported variable declarations
137!
138 logical, intent(in) :: Lcgini
139 integer, intent(in) :: ng, model, outLoop, innLoop, NinnLoop
140!
141! Local variable declarations.
142!
143 logical :: Ltrans, Laug
144 integer :: i, ic, j, iobs, ivec, Lscale, info
145
146 real(r8) :: zbet, eps, preducv, preducy
147 real(r8) :: cff, fact, ad_dla, adfac
148
149 real(r8), dimension(NinnLoop) :: zu, zgam, zfact, ad_zfact
150 real(r8), dimension(NinnLoop) :: ad_zu, ad_zrhs
151 real(r8), dimension(Ndatum(ng)+1) :: px, pgrad, pgrad_S
152 real(r8), dimension(Ndatum(ng)+1) :: ad_px, ad_pgrad
153 real(r8), dimension(Ninner,3) :: zwork, ad_zwork
154
155 character (len=13) :: string
156!
157!=======================================================================
158! Weak constraint 4DVar conjugate gradient, Lanczos-based algorithm.
159!=======================================================================
160!
161! This conjugate gradient algorithm is not run in parallel since the
162! weak constraint is done in observation space. Mostly all the
163! variables are 1D arrays. Therefore, in parallel applications (only
164! distributed-memory is possible) the master node does the computations
165! and then broadcasts results to remaining nodes in the communicator.
166!
167! This version of congrad solves A(x+x0)=b for x, by minimizing
168! J=0.5*x'Ax-x'(b+Ax0), where x0 is a first-guess estimate of the
169! representer coefficients from the previous outer-loop.
170! For the first outer-loop, x0=0. In the code x0=cg_pxsave and
171! x=px.
172!
173 ltrans=.false.
174!
175 master_thread : IF (master) THEN
176!
177! Set vgrad0=zgrad0 for the appropriat outer-loop.
178!
179 DO iobs=1,ndatum(ng)+1
180 vgrad0(iobs)=zgrad0(iobs,outloop)
181 END DO
182!
183! Clear arrays.
184!
185 admodval=0.0_r8
186 ad_pgrad=0.0_r8
187 ad_px=0.0_r8
188 ad_zu=0.0_r8
189 ad_zrhs=0.0_r8
190 ad_zfact=0.0_r8
191 ad_dla=0.0_r8
192 ad_zwork=0.0_r8
193!!
194!! IF (innLoop.eq.NinnLoop) THEN
195!! DO iobs=1,Ndatum(ng)
196!! BCKmodVal(iobs)=NLmodVal(iobs)
197!! END DO
198!! END IF
199!!
200# if defined ANL_EOFS || defined TRACE_EOFS
201!
202 IF (innloop.eq.ninnloop) THEN
203!
204! Clear adjoint arrays when innLoop=NinnLoop.
205!
206 DO j=1,ninnloop+1
207 DO iobs=1,ndatum(ng)+1
208 ad_zcglwk(iobs,j)=0.0_r8
209 ad_vcglwk(iobs,j)=0.0_r8
210 END DO
211 END DO
212 DO iobs=1,ndatum(ng)+1
213 ad_zgrad0(iobs)=0.0_r8
214 ad_vgrad0(iobs)=0.0_r8
215 ad_vgrad0s(iobs)=0.0_r8
216 END DO
217 DO iobs=1,ndatum(ng)
218# ifdef IMPACT_INNER
219 ad_obsval(iobs,innloop)=0.0_r8
220# else
221 ad_obsval(iobs)=0.0_r8
222# endif
223 ad_obserr(iobs)=0.0_r8
224 ad_obsscale(iobs)=0.0_r8
225 ad_hbk(iobs)=0.0_r8
226 ad_bckmodval(iobs)=0.0_r8
227 END DO
228 DO j=1,ninnloop+1
229 ad_cg_beta(j)=0.0_r8
230 ad_cg_qg(j)=0.0_r8
231 END DO
232 DO j=1,ninnloop
233 ad_cg_delta(j)=0.0_r8
234 END DO
235 ad_cg_gnorm(outloop)=0.0_r8
236 ad_cg_gnorm_v(outloop)=0.0_r8
237 ad_jb0(0)=0.0_r8
238 ad_jb0(outer)=0.0_r8
239 END IF
240# endif
241!
242!-----------------------------------------------------------------------
243! Initialization step, innloop = 0.
244!-----------------------------------------------------------------------
245!
246 minimize : IF (innloop.eq.0) THEN
247
248# if defined RBL4DVAR || \
249 defined rbl4dvar_ana_sensitivity || \
250 defined rbl4dvar_fct_sensitivity || \
251 defined tl_rbl4dvar
252!
253 DO iobs=1,ndatum(ng)
254 bckmodval(iobs)=nlmodval(iobs)
255 END DO
256# endif
257!
258! Selime code: ADmodVal=z1.
259!
260 DO iobs=1,ndatum(ng)
261!^ tl_ADmodVal(iobs)=tl_zgrad0(iobs,outLoop)
262!^
263 ad_zgrad0(iobs)=ad_zgrad0(iobs)+tlmodval(iobs)
264 tlmodval(iobs)=0.0_r8
265 END DO
266!
267! Selime code: zgrad0=z1.
268!
269! Selime code: vgrad0=r.
270!
271 DO iobs=1,ndatum(ng)+1
272!^ tl_zgrad0(iobs)=tl_pgrad(iobs)
273!^
274 ad_pgrad(iobs)=ad_pgrad(iobs)+ad_zgrad0(iobs)
275 ad_zgrad0(iobs)=0.0_r8
276 END DO
277!
278! If preconditioning, convert pgrad from v-space to y-space.
279!
280! Selime code: z1=G*r.
281!
282!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
283!^ Lscale=1 ! Spectral LMP
284!^ Laug=.TRUE.
285!^ CALL rprecond (ng, Lscale, Laug, outLoop, NinnLoop-1, &
286!^ & pgrad)
287!^ END IF
288!
289! AUGMENTED
290! Augment the residual vector.
291!
292!^ tl_vgrad0s(Ndatum(ng)+1)=tl_vgrad0(Ndatum(ng)+1)
293!^
294 ad_vgrad0(ndatum(ng)+1)=ad_vgrad0(ndatum(ng)+1)+ &
295 & ad_vgrad0s(ndatum(ng)+1)
296 ad_vgrad0s(ndatum(ng)+1)=0.0_r8
297!^ tl_vgrad0(Ndatum(ng)+1)=tl_pgrad(Ndatum(ng)+1)
298!^
299 ad_pgrad(ndatum(ng)+1)=ad_pgrad(ndatum(ng)+1)+ &
300 & ad_vgrad0(ndatum(ng)+1)
301 ad_vgrad0(ndatum(ng)+1)=0.0_r8
302!^ tl_pgrad(Ndatum(ng)+1)=0.0_r8
303!^
304 ad_pgrad(ndatum(ng)+1)=0.0_r8
305!
306! If this is the first pass of the inner loop, set up the vectors and
307! store the first guess. The initial starting guess is assumed to be
308! zero in which case the residual is just: (obs-model).
309!
310 DO iobs=1,ndatum(ng)
311# if defined RBL4DVAR || \
312 defined rbl4dvar_ana_sensitivity || \
313 defined rbl4dvar_fct_sensitivity || \
314 defined tl_rbl4dvar
315!^ pgrad(iobs)=ObsScale(iobs)* &
316!^ & (ObsVal(iobs)-BCKmodVal(iobs))
317# else
318!
319! NOTE: THIS CODE IS NOT CORRECT YET FOR R4D-Var
320!
321!^ pgrad(iobs)=ObsScale(iobs)* &
322!^ & (ObsVal(iobs)- &
323!^ TLmodVal_S(iobs,innLoop,outLoop))
324# endif
325!^ tl_vgrad0s(iobs)=tl_vgrad0(iobs)
326!^
327 ad_vgrad0(iobs)=ad_vgrad0(iobs)+ad_vgrad0s(iobs)
328 ad_vgrad0s(iobs)=0.0_r8
329!^ tl_vgrad0(iobs)=tl_pgrad(iobs)
330!^
331 ad_pgrad(iobs)=ad_pgrad(iobs)+ad_vgrad0(iobs)
332 ad_vgrad0(iobs)=0.0_r8
333!
334 IF (obserr(iobs).NE.0.0_r8) THEN
335!^ pgrad(iobs)=pgrad(iobs)/ObsErr(iobs)
336!^
337 pgrad(iobs)=vgrad0(iobs)
338!^ tl_pgrad(iobs)=tl_pgrad(iobs)/ObsErr(iobs)- &
339!^ & tl_ObsErr(iobs)*pgrad(iobs)/ObsErr(iobs)
340!^
341 ad_obserr(iobs)=ad_obserr(iobs)-ad_pgrad(iobs)* &
342 & pgrad(iobs)/obserr(iobs)
343 ad_pgrad(iobs)=ad_pgrad(iobs)/obserr(iobs)
344 END IF
345# if defined RBL4DVAR || \
346 defined rbl4dvar_ana_sensitivity || \
347 defined rbl4dvar_fct_sensitivity || \
348 defined tl_rbl4dvar
349!^ tl_pgrad(iobs)=ObsScale(iobs)* &
350!^ & (tl_ObsVal(iobs)-tl_BCKmodVal(iobs))+ &
351!^ & tl_ObsScale(iobs)* &
352!^ & (ObsVal(iobs)-BCKmodVal(iobs))
353# ifdef IMPACT_INNER
354 ad_obsval(iobs,innloop)=ad_obsval(iobs,innloop)+ &
355 & obsscale(iobs)*ad_pgrad(iobs)
356# else
357 ad_obsval(iobs)=ad_obsval(iobs)+ &
358 & obsscale(iobs)*ad_pgrad(iobs)
359# endif
360 ad_bckmodval(iobs)=ad_bckmodval(iobs)- &
361 & obsscale(iobs)*ad_pgrad(iobs)
362 ad_obsscale(iobs)=ad_obsscale(iobs)+ &
363 & ad_pgrad(iobs)* &
364 & (obsval(iobs)-bckmodval(iobs))
365 ad_pgrad(iobs)=0.0_r8
366# else
367! NOTE: THIS CODE IS NOT CORRECT YET FOR R4D-Var
368!
369!^ tl_pgrad(iobs)=ObsScale(iobs)* &
370!^ & (tl_ObsVal(iobs)-tl_TLmodVal(iobs))
371!^ tl_pgrad(iobs)=ObsScale(iobs)* &
372!^ & (tl_ObsVal(iobs)-TLmodVal(iobs))+ &
373!^ & tl_ObsScale(iobs)* &
374!^ & (ObsVal(iobs)- &
375!^ & TLmodVal_S(iobs,innLoop,outLoop))
376!^
377# ifdef IMPACT_INNER
378 ad_obsval(iobs,innloop)=ad_obsval(iobs,innloop)+ &
379 & obsscale(iobs)*ad_pgrad(iobs)
380# else
381 ad_obsval(iobs)=ad_obsval(iobs)+ &
382 & obsscale(iobs)*ad_pgrad(iobs)
383# endif
384 admodval(iobs)=admodval(iobs)- &
385 & obsscale(iobs)*ad_pgrad(iobs)
386 ad_obsscale(iobs)=ad_obsscale(iobs)+
387 & ad_pgrad(iobs)* &
388 & (obsval(iobs)- &
389 & tlmodval_s(iobs,innloop,outloop))
390 ad_pgrad(iobs)=0.0_r8
391# endif
392 END DO
393!
394! For the first outer-loop, x0=0.
395!
396 IF (outloop.eq.1) THEN
397 DO iobs=1,ndatum(ng)+1
398!^ tl_px(iobs)=0.0_r8
399!^
400 ad_px(iobs)=0.0_r8
401 END DO
402 DO iobs=1,ndatum(ng)
403!^ tl_Hbk(iobs)=0.0_r8
404!^
405 ad_hbk(iobs)=0.0_r8
406 END DO
407 ELSE
408!^ IF (Lprecond) THEN
409!^ DO iobs=1,Ndatum(ng)+1
410!^ FOURDVAR(ng)%cg_Dpxsum(iobs)=FOURDVAR(ng)%cg_pxsum(iobs)
411!^ END DO
412!^ Lscale=1 ! Spectral LMP
413!^ Laug=.FALSE.
414!^ CALL rprecond (ng, Lscale, Laug, outLoop, NinnLoop-1, &
415!^ & FOURDVAR(ng)%cg_Dpxsum)
416!^ DO iobs=1,Ndatum(ng)+1
417!^ FOURDVAR(ng)%cg_Dpxsum(iobs)= &
418!^ & -FOURDVAR(ng)%cg_Dpxsum(iobs)+ &
419!^ & FOURDVAR(ng)%cg_pxsum(iobs)
420!^ END DO
421!^ END IF
422!!
423!! NOTE: CODE WON'T WORK FOR R4D-Var AT THE MOMENT SINCE THIS IS
424!! THE WRONG TLmodVal HERE!!!
425!!
426 DO iobs=1,ndatum(ng)
427!^ tl_Hbk(iobs)=-tl_TLmodVal(iobs)
428!^
429 admodval(iobs)=admodval(iobs)-ad_hbk(iobs)
430 ad_hbk(iobs)=0.0_r8
431 END DO
432 END IF
433# if defined RBL4DVAR_ANA_SENSITIVITY || \
434 defined rbl4dvar_fct_sensitivity
435!
436 DO iobs=1,ndatum(ng)
437!^ tl_BCKmodVal(iobs)=0.0_r8
438!^
439 ad_bckmodval(iobs)=0.0_r8
440 END DO
441# endif
442!
443!-----------------------------------------------------------------------
444! Minimization step, innLoop > 0.
445!-----------------------------------------------------------------------
446!
447 ELSE
448!
449! ADJOINT NEXT CALCULATION
450!
451! Put the new trial solution into the adjoint vector for the next
452! loop or put the final solution into the adjoint vector on the
453! final inner-loop.
454!
455 IF (innloop.eq.ninnloop) THEN
456!
457! Note: px is already in v-space so there is no need to convert
458! if preconditioning; cg_pxsave is also in v-space.
459!
460 DO iobs=1,ndatum(ng)+1
461!^ FOURDVAR(ng)%tl_cg_pxsave(iobs)=tl_px(iobs)
462!^
463 ad_px(iobs)=ad_px(iobs)+fourdvar(ng)%ad_cg_pxsave(iobs)
464 fourdvar(ng)%ad_cg_pxsave(iobs)=0.0_r8
465 END DO
466
467 DO iobs=1,ndatum(ng)
468!^ tl_ADmodVal(iobs)=tl_px(iobs)
469!^
470 ad_px(iobs)=ad_px(iobs)+tlmodval(iobs)
471 tlmodval(iobs)=0.0_r8
472 END DO
473 ELSE
474!
475! Selime code: z1=ADmodVal.
476!
477 DO iobs=1,ndatum(ng)
478!^ tl_ADmodVal(iobs)=tl_vcglwk(iobs,innLoop+1)
479 ad_vcglwk(iobs,innloop+1)=ad_vcglwk(iobs,innloop+1)+ &
480 & tlmodval(iobs)
481 tlmodval(iobs)=0.0_r8
482 END DO
483!
484 END IF
485!
486! Selime code: z1=vcglwk.
487!
488! AUGMENTED.
489!
490 DO iobs=1,ndatum(ng)+1
491!^ tl_vcglwk(iobs,innLoop+1)=tl_pgrad(iobs)
492!^
493 ad_pgrad(iobs)=ad_pgrad(iobs)+ad_vcglwk(iobs,innloop+1)
494 ad_vcglwk(iobs,innloop+1)=0.0_r8
495 END DO
496
497!
498! If preconditioning, convert pgrad from v-space to y-space.
499!
500!
501! Selime code: z1=Gw.
502!
503!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
504!^ Lscale=1 ! Spectral LMP
505!^ Laug=.TRUE.
506!^ CALL rprecond (ng, Lscale, Laug, outLoop, NinnLoop-1, pgrad)
507!^ END IF
508!
509! Save the non-normalized Lanczos vector. zcglwk is used as temporary
510! storage. The calculation will be completed at the start of the next
511! inner-loop.
512!
513! Selime code: w=zcglwk.
514!
515! AUGMENTED
516!
517 DO iobs=1,ndatum(ng)+1
518!^ tl_zcglwk(iobs,innLoop+1)=tl_pgrad(iobs)
519!^
520 ad_pgrad(iobs)=ad_pgrad(iobs)+ad_zcglwk(iobs,innloop+1)
521 ad_zcglwk(iobs,innloop+1)=0.0_r8
522 END DO
523!
524! Selime code: cg_dla=pgrad*HBH'*v1=pgrad*TLmodVal_S/zfact
525!
526! AUGMENTED
527!
528! Compute basic state pgrad.
529!
530 IF (innloop.eq.1) THEN
531 fact=1.0_r8/cg_gnorm_v(outloop)
532 ELSE
533 fact=1.0_r8/cg_beta(innloop,outloop)
534 END IF
535 DO iobs=1,ndatum(ng)
536 pgrad(iobs)=tlmodval_s(iobs,innloop,outloop)*fact+ &
537 & hbk(iobs,outloop)*vcglwk(ndatum(ng)+1,innloop,outloop)
538 END DO
539 DO iobs=1,ndatum(ng)
540 IF (obserr(iobs).NE.0.0_r8) THEN
541 pgrad(iobs)=pgrad(iobs)/obserr(iobs)
542 END IF
543 END DO
544 pgrad(ndatum(ng)+1)=0.0_r8
545 DO iobs=1,ndatum(ng)
546 pgrad(iobs)=pgrad(iobs)+vcglwk(iobs,innloop,outloop)
547 END DO
548 pgrad(ndatum(ng)+1)=pgrad(ndatum(ng)+1)+ &
549 & vcglwk(ndatum(ng)+1,innloop,outloop)
550 IF (innloop.gt.1) THEN
551 DO iobs=1,ndatum(ng)+1
552 pgrad(iobs)=pgrad(iobs)-cg_beta(innloop,outloop)* &
553 & zcglwk(iobs,innloop-1,outloop)
554 END DO
555 END IF
556 DO iobs=1,ndatum(ng)+1
557 pgrad(iobs)=pgrad(iobs)-cg_delta(innloop,outloop)* &
558 & zcglwk(iobs,innloop,outloop)
559 END DO
560!
561! Save the intermediate value of pgrad.
562!
563 DO iobs=1,ndatum(ng)+1
564 pgrad_s(iobs)=pgrad(iobs)
565 END DO
566!
567! Orthonormalize against previous directions.
568!
569! Identify the appropriate Lanczos vector normalization coefficient.
570!
571 DO ivec=1,innloop
572 IF (ivec.eq.1) THEN
573 zfact(ivec)=1.0_r8/cg_gnorm_v(outloop)
574 ELSE
575 zfact(ivec)=1.0_r8/cg_beta(ivec,outloop)
576 ENDIF
577 END DO
578!
579! Selime code: cg_dla=pgrad*HBH'*v1=pgrad*TLmodVal_S/zfact
580!
581! AUGMENTED
582
583!^ DO ivec=innLoop,1,-1
584!^
585 DO ivec=1,innloop
586 DO iobs=1,ndatum(ng)+1
587 pgrad(iobs)=pgrad_s(iobs)
588 END DO
589 DO i=ivec+1,innloop
590 DO iobs=1,ndatum(ng)+1
591 pgrad(iobs)=pgrad(iobs)- &
592 & cg_dla(i,outloop)* &
593 & vcglwk(iobs,i,outloop)
594 END DO
595 END DO
596 DO iobs=1,ndatum(ng)+1
597!^ tl_pgrad(iobs)=tl_pgrad(iobs)-tl_dla* &
598!^ & vcglwk(iobs,ivec,outLoop)- &
599!^ & cg_dla(ivec,outLoop)* &
600!^ & tl_vcglwk(iobs,ivec)
601!^
602 ad_dla=ad_dla-ad_pgrad(iobs)*vcglwk(iobs,ivec,outloop)
603 ad_vcglwk(iobs,ivec)=ad_vcglwk(iobs,ivec)- &
604 & cg_dla(ivec,outloop)*ad_pgrad(iobs)
605 END DO
606!^ tl_dla=tl_dla+ &
607!^ & tl_Jb0(outLoop-1)* &
608!^ & vcglwk(Ndatum(ng)+1,ivec,outLoop)* &
609!^ & pgrad(Ndatum(ng)+1)+ &
610!^ & Jb0(outLoop-1)*tl_vcglwk(Ndatum(ng)+1,ivec)* &
611!^ & pgrad(Ndatum(ng)+1)+ &
612!^ & Jb0(outLoop-1)*vcglwk(Ndatum(ng)+1,ivec,outLoop)* &
613!^ & tl_pgrad(Ndatum(ng)+1)
614!^
615 ad_jb0(outloop-1)=ad_jb0(outloop-1)+ &
616 & vcglwk(ndatum(ng)+1,ivec,outloop)* &
617 & pgrad(ndatum(ng)+1)*ad_dla
618 ad_vcglwk(ndatum(ng)+1,ivec)=ad_vcglwk(ndatum(ng)+1,ivec)+ &
619 & jb0(outloop-1)* &
620 & pgrad(ndatum(ng)+1)*ad_dla
621 ad_pgrad(ndatum(ng)+1)=ad_pgrad(ndatum(ng)+1)+ &
622 & jb0(outloop-1)* &
623 & vcglwk(ndatum(ng)+1,ivec,outloop)* &
624 & ad_dla
625 DO iobs=1,ndatum(ng)
626!^ tl_dla=tl_dla+ &
627!^ & tl_pgrad(iobs)* &
628!^ & TLmodVal_S(iobs,ivec,outLoop)* &
629!^ & zfact(ivec)+ &
630!^ & pgrad(iobs)*tl_TLmodVal_S(iobs,ivec,outLoop)* &
631!^ & zfact(ivec)+ &
632!^ & pgrad(iobs)*TLmodVal_S(iobs,ivec,outLoop)* &
633!^ & tl_zfact(ivec)+ &
634!^ & tl_pgrad(iobs)* &
635!^ & Hbk(iobs,outLoop)* &
636!^ & vcglwk(Ndatum(ng)+1,ivec,outLoop)+ &
637!^ & pgrad(iobs)*tl_Hbk(iobs)* &
638!^ & vcglwk(Ndatum(ng)+1,ivec,outLoop)+ &
639!^ & pgrad(iobs)*Hbk(iobs,outLoop)* &
640!^ & tl_vcglwk(Ndatum(ng)+1,ivec)+ &
641!^ & tl_Hbk(iobs)*vcglwk(iobs,ivec,outLoop)* &
642!^ & pgrad(Ndatum(ng)+1)+ &
643!^ & Hbk(iobs,outLoop)*tl_vcglwk(iobs,ivec)* &
644!^ & pgrad(Ndatum(ng)+1)+ &
645!^ & Hbk(iobs,outLoop)*vcglwk(iobs,ivec,outLoop)* &
646!^ & tl_pgrad(Ndatum(ng)+1)
647!^
648 ad_pgrad(iobs)=ad_pgrad(iobs)+ &
649 & tlmodval_s(iobs,ivec,outloop)* &
650 & zfact(ivec)*ad_dla
651 ad_tlmodval_s(iobs,ivec,outloop)= &
652 & ad_tlmodval_s(iobs,ivec,outloop)+ &
653 & pgrad(iobs)*zfact(ivec)*ad_dla
654 ad_zfact(ivec)=ad_zfact(ivec)+ &
655 & pgrad(iobs)* &
656 & tlmodval_s(iobs,ivec,outloop)*ad_dla
657 ad_pgrad(iobs)=ad_pgrad(iobs)+ &
658 & hbk(iobs,outloop)* &
659 & vcglwk(ndatum(ng)+1,ivec,outloop)*ad_dla
660 ad_hbk(iobs)=ad_hbk(iobs)+ &
661 & pgrad(iobs)* &
662 & vcglwk(ndatum(ng)+1,ivec,outloop)*ad_dla
663 ad_vcglwk(ndatum(ng)+1,ivec)=ad_vcglwk(ndatum(ng)+1,ivec)+&
664 & pgrad(iobs)* &
665 & hbk(iobs,outloop)*ad_dla
666 ad_hbk(iobs)=ad_hbk(iobs)+ &
667 & vcglwk(iobs,ivec,outloop)* &
668 & pgrad(ndatum(ng)+1)*ad_dla
669 ad_vcglwk(iobs,ivec)=ad_vcglwk(iobs,ivec)+ &
670 & hbk(iobs,outloop)* &
671 & pgrad(ndatum(ng)+1)*ad_dla
672 ad_pgrad(ndatum(ng)+1)=ad_pgrad(ndatum(ng)+1)+ &
673 & hbk(iobs,outloop)* &
674 & vcglwk(iobs,ivec,outloop)*ad_dla
675 END DO
676!^ tl_dla=0.0_r8
677!^
678 ad_dla=0.0_r8
679 END DO
680!
681 DO ivec=1,innloop
682 IF (ivec.eq.1) THEN
683 zfact(ivec)=1.0_r8/cg_gnorm_v(outloop)
684!^ tl_zfact(ivec)=-tl_Gnorm_v(outLoop)*zfact(ivec)/ &
685!^ & cg_Gnorm_v(outLoop)
686!^
687 ad_cg_gnorm_v(outloop)=ad_cg_gnorm_v(outloop)- &
688 & ad_zfact(ivec)*zfact(ivec)/ &
689 & cg_gnorm_v(outloop)
690 ad_zfact(ivec)=0.0_r8
691 ELSE
692 zfact(ivec)=1.0_r8/cg_beta(ivec,outloop)
693!^ tl_zfact(ivec)=-tl_cg_beta(ivec)*zfact(ivec)/ &
694!^ & cg_beta(ivec,outLoop)
695 ad_cg_beta(ivec)=ad_cg_beta(ivec)- &
696 & ad_zfact(ivec)*zfact(ivec)/ &
697 & cg_beta(ivec,outloop)
698 ad_zfact(ivec)=0.0_r8
699 ENDIF
700 END DO
701!
702! Selime code: w=w-alpha*v1
703!
704! AUGMENTED.
705!
706 DO iobs=1,ndatum(ng)+1
707!^ tl_pgrad(iobs)=tl_pgrad(iobs)- &
708!^ & tl_cg_delta(innLoop)* &
709!^ & zcglwk(iobs,innLoop,outLoop)- &
710!^ & cg_delta(innLoop,outLoop)* &
711!^ & tl_zcglwk(iobs,innLoop)
712!^
713 ad_cg_delta(innloop)=ad_cg_delta(innloop)- &
714 & zcglwk(iobs,innloop,outloop)* &
715 & ad_pgrad(iobs)
716 ad_zcglwk(iobs,innloop)=ad_zcglwk(iobs,innloop)- &
717 & cg_delta(innloop,outloop)* &
718 & ad_pgrad(iobs)
719 END DO
720!
721! Selime code: alpha=w'*t1=cg_delta
722!
723! AUGMENTED.
724!
725! Recompute basic state pgrad.
726!
727 IF (innloop.eq.1) THEN
728 fact=1.0_r8/cg_gnorm_v(outloop)
729 ELSE
730 fact=1.0_r8/cg_beta(innloop,outloop)
731 END IF
732 DO iobs=1,ndatum(ng)
733 pgrad(iobs)=tlmodval_s(iobs,innloop,outloop)*fact+ &
734 & hbk(iobs,outloop)* &
735 & vcglwk(ndatum(ng)+1,innloop,outloop)
736 END DO
737 DO iobs=1,ndatum(ng)
738 IF (obserr(iobs).NE.0.0_r8) THEN
739 pgrad(iobs)=pgrad(iobs)/obserr(iobs)
740 END IF
741 END DO
742 pgrad(ndatum(ng)+1)=0.0_r8
743 DO iobs=1,ndatum(ng)
744 pgrad(iobs)=pgrad(iobs)+vcglwk(iobs,innloop,outloop)
745 END DO
746 pgrad(ndatum(ng)+1)=pgrad(ndatum(ng)+1)+ &
747 & vcglwk(ndatum(ng)+1,innloop,outloop)
748 IF (innloop.gt.1) THEN
749 DO iobs=1,ndatum(ng)+1
750 pgrad(iobs)=pgrad(iobs)- &
751 & cg_beta(innloop,outloop)* &
752 & zcglwk(iobs,innloop-1,outloop)
753 END DO
754 END IF
755!
756!^ tl_cg_delta(innLoop)=tl_cg_delta(innLoop)+ &
757!^ & tl_Jb0(outLoop-1)* &
758!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
759!^ & pgrad(Ndatum(ng)+1)+ &
760!^ & Jb0(outLoop-1)* &
761!^ & tl_vcglwk(Ndatum(ng)+1,innLoop)* &
762!^ & pgrad(Ndatum(ng)+1)+ &
763!^ & Jb0(outLoop-1)* &
764!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
765!^ & tl_pgrad(Ndatum(ng)+1)
766!^
767 ad_jb0(outloop-1)=ad_jb0(outloop-1)+ &
768 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
769 & pgrad(ndatum(ng)+1)* &
770 & ad_cg_delta(innloop)
771 ad_vcglwk(ndatum(ng)+1,innloop)= &
772 & ad_vcglwk(ndatum(ng)+1,innloop)+ &
773 & jb0(outloop-1)* &
774 & pgrad(ndatum(ng)+1)* &
775 & ad_cg_delta(innloop)
776 ad_pgrad(ndatum(ng)+1)=ad_pgrad(ndatum(ng)+1)+ &
777 & jb0(outloop-1)* &
778 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
779 & ad_cg_delta(innloop)
780!
781 IF (innloop.eq.1) THEN
782 fact=1.0_r8/cg_gnorm_v(outloop)
783 ELSE
784 fact=1.0_r8/cg_beta(innloop,outloop)
785 END IF
786!
787 DO iobs=1,ndatum(ng)
788!^ tl_cg_delta(innLoop)=tl_cg_delta(innLoop)+ &
789!^ & tl_TLmodVal(iobs)* &
790!^ & pgrad(iobs)+ &
791!^ & TLmodVal_S(iobs,innLoop,outLoop)* &
792!^ & tl_pgrad(iobs)* &
793!^ & fact+ &
794!^ & tl_Hbk(iobs)* &
795!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
796!^ & pgrad(iobs)+ &
797!^ & Hbk(iobs,outLoop)* &
798!^ & tl_vcglwk(Ndatum(ng)+1,innLoop)* &
799!^ & pgrad(iobs)+ &
800!^ & Hbk(iobs,outLoop)* &
801!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
802!^ & tl_pgrad(iobs)+ &
803!^ & tl_Hbk(iobs)* &
804!^ & vcglwk(iobs,innLoop,outLoop)* &
805!^ & pgrad(Ndatum(ng)+1)+ &
806!^ & Hbk(iobs,outLoop)* &
807!^ & tl_vcglwk(iobs,innLoop)* &
808!^ & pgrad(Ndatum(ng)+1)+ &
809!^ & Hbk(iobs,outLoop)* &
810!^ & vcglwk(iobs,innLoop,outLoop)* &
811!^ & tl_pgrad(Ndatum(ng)+1)
812!^
813 admodval(iobs)=admodval(iobs)+ &
814 & pgrad(iobs)*ad_cg_delta(innloop)
815 ad_pgrad(iobs)=ad_pgrad(iobs)+ &
816 & tlmodval_s(iobs,innloop,outloop)* &
817 & ad_cg_delta(innloop)* &
818 & fact
819 ad_hbk(iobs)=ad_hbk(iobs)+ &
820 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
821 & pgrad(iobs)* &
822 & ad_cg_delta(innloop)
823 ad_vcglwk(ndatum(ng)+1,innloop)= &
824 & ad_vcglwk(ndatum(ng)+1,innloop)+ &
825 & hbk(iobs,outloop)* &
826 & pgrad(iobs)* &
827 & ad_cg_delta(innloop)
828 ad_pgrad(iobs)=ad_pgrad(iobs)+ &
829 & hbk(iobs,outloop)* &
830 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
831 & ad_cg_delta(innloop)
832 ad_hbk(iobs)=ad_hbk(iobs)+ &
833 & vcglwk(iobs,innloop,outloop)* &
834 & pgrad(ndatum(ng)+1)* &
835 & ad_cg_delta(innloop)
836 ad_vcglwk(iobs,innloop)=ad_vcglwk(iobs,innloop)+ &
837 & hbk(iobs,outloop)* &
838 & pgrad(ndatum(ng)+1)* &
839 & ad_cg_delta(innloop)
840 ad_pgrad(ndatum(ng)+1)=ad_pgrad(ndatum(ng)+1)+ &
841 & hbk(iobs,outloop)* &
842 & vcglwk(iobs,innloop,outloop)* &
843 & ad_cg_delta(innloop)
844 END DO
845!
846!^ tl_cg_delta(innLoop)=0.0_r8
847!^
848 ad_cg_delta(innloop)=0.0_r8
849
850 IF (innloop.gt.1) THEN
851!
852! Selime code: w=w-v0*beta
853!
854! AUGMENTED.
855!
856 DO iobs=1,ndatum(ng)+1
857!^ tl_pgrad(iobs)=tl_pgrad(iobs)- &
858!^ & tl_cg_beta(innLoop)* &
859!^ & zcglwk(iobs,innLoop-1,outLoop)- &
860!^ & cg_beta(innLoop,outLoop)* &
861!^ & tl_zcglwk(iobs,innLoop-1)
862!^
863 ad_cg_beta(innloop)=ad_cg_beta(innloop)- &
864 & zcglwk(iobs,innloop-1,outloop)* &
865 & ad_pgrad(iobs)
866 ad_zcglwk(iobs,innloop-1)=ad_zcglwk(iobs,innloop-1)- &
867 & cg_beta(innloop,outloop)* &
868 & ad_pgrad(iobs)
869 END DO
870 END IF
871!
872! Selime code: w=t1/R+z1
873!
874!^ tl_pgrad(Ndatum(ng)+1)=tl_pgrad(Ndatum(ng)+1)+ &
875!^ & tl_vcglwk(Ndatum(ng)+1,innLoop)
876!^
877 ad_vcglwk(ndatum(ng)+1,innloop)= &
878 & ad_vcglwk(ndatum(ng)+1,innloop)+ &
879 & ad_pgrad(ndatum(ng)+1)
880 DO iobs=1,ndatum(ng)
881!^ tl_pgrad(iobs)=tl_pgrad(iobs)+tl_vcglwk(iobs,innLoop)
882!^
883 ad_vcglwk(iobs,innloop)=ad_vcglwk(iobs,innloop)+ &
884 & ad_pgrad(iobs)
885 END DO
886!
887!^ tl_pgrad(Ndatum(ng)+1)=0.0_r8
888!^
889 ad_pgrad(ndatum(ng)+1)=0.0_r8
890!
891! Recompute pgrad.
892!
893 IF (innloop.eq.1) THEN
894 fact=1.0_r8/cg_gnorm_v(outloop)
895 ELSE
896 fact=1.0_r8/cg_beta(innloop,outloop)
897 END IF
898 DO iobs=1,ndatum(ng)
899 pgrad(iobs)=tlmodval_s(iobs,innloop,outloop)*fact+ &
900 & hbk(iobs,outloop)* &
901 & vcglwk(ndatum(ng)+1,innloop,outloop)
902 END DO
903!
904 DO iobs=1,ndatum(ng)
905!
906! Convert gradient contribution from x-space to v-space.
907!
908 IF (obserr(iobs).NE.0.0_r8) THEN
909 pgrad(iobs)=pgrad(iobs)/obserr(iobs)
910!
911! Selime code: w=t1/R. No augmented term here since Raug^-1=[R^-1 0].
912!
913!^ tl_pgrad(iobs)=tl_pgrad(iobs)/ObsErr(iobs)- &
914!^ & tl_ObsErr(iobs)*pgrad(iobs)/ObsErr(iobs)
915!^
916 ad_obserr(iobs)=ad_obserr(iobs)- &
917 & pgrad(iobs)*ad_pgrad(iobs)/obserr(iobs)
918 ad_pgrad(iobs)=ad_pgrad(iobs)/obserr(iobs)
919 END IF
920 END DO
921!
922 IF (innloop.eq.1) THEN
923 fact=1.0_r8/cg_gnorm_v(outloop)
924 ELSE
925 fact=1.0_r8/cg_beta(innloop,outloop)
926 END IF
927!
928 DO iobs=1,ndatum(ng)
929!
930! Selime code: w=t1.
931!
932!^ tl_pgrad(iobs)=tl_TLmodVal(iobs)+ &
933!^ & tl_Hbk(iobs)* &
934!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)+ &
935!^ & Hbk(iobs,outLoop)* &
936!^ & tl_vcglwk(Ndatum(ng)+1,innLoop)
937!^
938 admodval(iobs)=admodval(iobs)+ad_pgrad(iobs)
939 ad_hbk(iobs)=ad_hbk(iobs)+ &
940 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
941 & ad_pgrad(iobs)
942 ad_vcglwk(ndatum(ng)+1,innloop)= &
943 & ad_vcglwk(ndatum(ng)+1,innloop)+ &
944 & hbk(iobs,outloop)* &
945 & ad_pgrad(iobs)
946 ad_pgrad(iobs)=0.0_r8
947 END DO
948
949!
950! ADJOINT OF COMPLETE THE CALCULATION OF PREVIOUS LANCZOS VECTOR.
951!
952! Complete the computation of the Lanczos vector from previous
953! inner-loop
954!
955 IF (innloop.EQ.1) THEN
956!
957 DO iobs=1,ndatum(ng)
958!
959! Selime code: t1=t1/beta_first=TLmodVal
960!
961!^ tl_TLmodVal(iobs)=tl_TLmodVal(iobs)/cg_Gnorm_v(outLoop)- &
962!^ & tl_cg_Gnorm_v(outLoop)* &
963!^ & TLmodVal_S(iobs,innLoop,outLoop)/ &
964!^ & (cg_Gnorm_v(outLoop)* &
965!^ & cg_Gnorm_v(outLoop))
966!^
967 ad_cg_gnorm_v(outloop)=ad_cg_gnorm_v(outloop)- &
968 & admodval(iobs)* &
969 & tlmodval_s(iobs,innloop,outloop)/ &
970 & (cg_gnorm_v(outloop)* &
971 & cg_gnorm_v(outloop))
972 admodval(iobs)=admodval(iobs)/cg_gnorm_v(outloop)
973!^ tl_TLmodVal_S(iobs,innLoop,outLoop)=tl_TLmodVal(iobs)
974!^
975 admodval(iobs)=admodval(iobs)+ &
976 & ad_tlmodval_s(iobs,innloop,outloop)
977 ad_tlmodval_s(iobs,innloop,outloop)=0.0_r8
978 END DO
979!
980! Normalize TLmodVal here and save in TLmodVal_S.
981!
982! AUGMENTED.
983!
984 DO iobs=1,ndatum(ng)+1
985!
986! Selime code: z1=z1/beta_first=vcglwk(1).
987!
988!^ tl_vcglwk(iobs,1)=tl_zgrad0(iobs)/cg_Gnorm_v(outLoop)- &
989!^ & tl_cg_Gnorm_v(outLoop)* &
990!^ & vcglwk(iobs,1,outLoop)/ &
991!^ & cg_Gnorm_v(outLoop)
992!^
993 ad_zgrad0(iobs)=ad_zgrad0(iobs)+ &
994 & ad_vcglwk(iobs,1)/cg_gnorm_v(outloop)
995 ad_cg_gnorm_v(outloop)=ad_cg_gnorm_v(outloop)- &
996 & vcglwk(iobs,1,outloop)* &
997 & ad_vcglwk(iobs,1)/ &
998 & cg_gnorm_v(outloop)
999 ad_vcglwk(iobs,1)=0.0_r8
1000!
1001! Selime code: v1=r/beta_first=zcglwk(1).
1002!
1003!^ tl_zcglwk(iobs,1)=tl_vgrad0(iobs)/ &
1004!^ & cg_Gnorm_v(outLoop)- &
1005!^ & tl_cg_Gnorm_v(outLoop)* &
1006!^ & zcglwk(iobs,1,outLoop)/ &
1007!^ & cg_Gnorm_v(outLoop)
1008!^
1009 ad_vgrad0(iobs)=ad_vgrad0(iobs)+ad_zcglwk(iobs,1)/ &
1010 & cg_gnorm_v(outloop)
1011 ad_cg_gnorm_v(outloop)=ad_cg_gnorm_v(outloop)- &
1012 & zcglwk(iobs,1,outloop)* &
1013 & ad_zcglwk(iobs,1)/ &
1014 & cg_gnorm_v(outloop)
1015 ad_zcglwk(iobs,1)=0.0_r8
1016 END DO
1017!
1018! Selime code: beta_first=sqrt(r't1)=cg_Gnorm_v.
1019! cg_QG(1)=beta_first
1020!
1021!^ tl_cg_QG(1)=tl_cg_Gnorm_v(outLoop)
1022!^
1023 ad_cg_gnorm_v(outloop)=ad_cg_gnorm_v(outloop)+ad_cg_qg(1)
1024 ad_cg_qg(1)=0.0_r8
1025!^ tl_cg_Gnorm_v(outLoop)=0.5_r8*tl_cg_Gnorm_v(outLoop)/ &
1026!^ & cg_Gnorm_v(outLoop)
1027!^
1028 ad_cg_gnorm_v(outloop)=0.5_r8*ad_cg_gnorm_v(outloop)/ &
1029 & cg_gnorm_v(outloop)
1030!
1031!^ tl_cg_Gnorm_v(outLoop)=tl_cg_Gnorm_v(outLoop)+ &
1032!^ & tl_Jb0(outLoop-1)* &
1033!^ & vgrad0(Ndatum(ng)+1)* &
1034!^ & zgrad0(Ndatum(ng)+1,outLoop)+ &
1035!^ & Jb0(outLoop-1)* &
1036!^ & tl_vgrad0(Ndatum(ng)+1)* &
1037!^ & zgrad0(Ndatum(ng)+1,outLoop)+ &
1038!^ & Jb0(outLoop-1)* &
1039!^ & vgrad0(Ndatum(ng)+1)* &
1040!^ & tl_zgrad0(Ndatum(ng)+1)
1041!^
1042 ad_jb0(outloop-1)=ad_jb0(outloop-1)+ &
1043 & vgrad0(ndatum(ng)+1)* &
1044 & zgrad0(ndatum(ng)+1,outloop)* &
1045 & ad_cg_gnorm_v(outloop)
1046 ad_vgrad0(ndatum(ng)+1)=ad_vgrad0(ndatum(ng)+1)+ &
1047 & jb0(outloop-1)* &
1048 & zgrad0(ndatum(ng)+1,outloop)* &
1049 & ad_cg_gnorm_v(outloop)
1050 ad_zgrad0(ndatum(ng)+1)=ad_zgrad0(ndatum(ng)+1)+ &
1051 & jb0(outloop-1)* &
1052 & vgrad0(ndatum(ng)+1)* &
1053 & ad_cg_gnorm_v(outloop)
1054 DO iobs=1,ndatum(ng)
1055!^ tl_cg_Gnorm_v(outLoop)=tl_cg_Gnorm_v(outLoop)+ &
1056!^ & tl_Hbk(iobs)* &
1057!^ & zgrad0(Ndatum(ng)+1,outLoop)* &
1058!^ & vgrad0(iobs)+ &
1059!^ & Hbk(iobs,outLoop)* &
1060!^ & tl_zgrad0(Ndatum(ng)+1)* &
1061!^ & vgrad0(iobs)+ &
1062!^ & Hbk(iobs,outLoop)* &
1063!^ & zgrad0(Ndatum(ng)+1,outLoop)* &
1064!^ & tl_vgrad0(iobs)+ &
1065!^ & tl_Hbk(iobs)* &
1066!^ & vgrad0(Ndatum(ng)+1)* &
1067!^ & zgrad0(iobs,outLoop)+ &
1068!^ & Hbk(iobs,outLoop)* &
1069!^ & tl_vgrad0(Ndatum(ng)+1)* &
1070!^ & zgrad0(iobs,outLoop)+ &
1071!^ & Hbk(iobs,outLoop)* &
1072!^ & vgrad0(Ndatum(ng)+1)* &
1073!^ & tl_zgrad0(iobs)
1074!^
1075 ad_hbk(iobs)=ad_hbk(iobs)+ &
1076 & zgrad0(ndatum(ng)+1,outloop)* &
1077 & vgrad0(iobs)* &
1078 & ad_cg_gnorm_v(outloop)
1079 ad_zgrad0(ndatum(ng)+1)=ad_zgrad0(ndatum(ng)+1)+ &
1080 & hbk(iobs,outloop)* &
1081 & vgrad0(iobs)* &
1082 & ad_cg_gnorm_v(outloop)
1083 ad_vgrad0(iobs)=ad_vgrad0(iobs)+ &
1084 & hbk(iobs,outloop)* &
1085 & zgrad0(ndatum(ng)+1,outloop)* &
1086 & ad_cg_gnorm_v(outloop)
1087 ad_hbk(iobs)=ad_hbk(iobs)+ &
1088 & vgrad0(ndatum(ng)+1)* &
1089 & zgrad0(iobs,outloop)* &
1090 & ad_cg_gnorm_v(outloop)
1091 ad_vgrad0(ndatum(ng)+1)=ad_vgrad0(ndatum(ng)+1)+ &
1092 & hbk(iobs,outloop)* &
1093 & zgrad0(iobs,outloop)* &
1094 & ad_cg_gnorm_v(outloop)
1095 ad_zgrad0(iobs)=ad_zgrad0(iobs)+ &
1096 & hbk(iobs,outloop)* &
1097 & vgrad0(ndatum(ng)+1)* &
1098 & ad_cg_gnorm_v(outloop)
1099 END DO
1100!
1101 DO iobs=1,ndatum(ng)
1102!^ tl_cg_Gnorm_v(outLoop)=tl_cg_Gnorm_v(outLoop)+ &
1103!^ & tl_TLmodVal(iobs)* &
1104!^ & vgrad0(iobs)+ &
1105!^ & TLmodVal_S(iobs,innLoop,outLoop)* &
1106!^ & tl_vgrad0(iobs)
1107!^
1108 admodval(iobs)=admodval(iobs)+ &
1109 & vgrad0(iobs)* &
1110 & ad_cg_gnorm_v(outloop)
1111 ad_vgrad0(iobs)=ad_vgrad0(iobs)+ &
1112 & tlmodval_s(iobs,innloop,outloop)* &
1113 & ad_cg_gnorm_v(outloop)
1114 END DO
1115!^ tl_cg_Gnorm_v(outLoop)=0.0_r8
1116!^
1117 ad_cg_gnorm_v(outloop)=0.0_r8
1118
1119 ELSE
1120!
1121 IF (innloop.eq.ninnloop) THEN
1122!
1123! If preconditioning, convert px from y-space to v-space.
1124! We will always keep px in v-space.
1125!
1126!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
1127!^ Lscale=1 ! Spectral LMP
1128!^ Laug=.TRUE.
1129!^ CALL rprecond(ng, Lscale, Laug, outLoop, NinnLoop-1, px)
1130!^ END IF
1131!
1132! Recompute zu and zwork.
1133!
1134!
1135 zbet=cg_delta(1,outloop)
1136 zu(1)=cg_qg(1,outloop)/zbet
1137 DO ivec=2,innloop-1
1138 zgam(ivec)=cg_beta(ivec,outloop)/zbet
1139 zbet=cg_delta(ivec,outloop)- &
1140 & cg_beta(ivec,outloop)*zgam(ivec)
1141 zu(ivec)=(cg_qg(ivec,outloop)-cg_beta(ivec,outloop)* &
1142 & zu(ivec-1))/zbet
1143 END DO
1144 zwork(innloop-1,3)=zu(innloop-1)
1145!
1146 DO ivec=innloop-2,1,-1
1147 zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
1148 zwork(ivec,3)=zu(ivec)
1149 END DO
1150!
1151 DO iobs=1,ndatum(ng)+1
1152 DO ivec=1,innloop-1
1153!^ tl_px(iobs)=tl_px(iobs)+ &
1154!^ & tl_zcglwk(iobs,ivec)* &
1155!^ & zwork(ivec,3)+ &
1156!^ & zcglwk(iobs,ivec,outLoop)* &
1157!^ & tl_zwork(ivec,3)
1158!^
1159 ad_zcglwk(iobs,ivec)=ad_zcglwk(iobs,ivec)+ &
1160 & zwork(ivec,3)* &
1161 & ad_px(iobs)
1162 ad_zwork(ivec,3)=ad_zwork(ivec,3)+ &
1163 & zcglwk(iobs,ivec,outloop)* &
1164 & ad_px(iobs)
1165 END DO
1166!^ tl_px(iobs)=0.0_r8
1167!^
1168 ad_px(iobs)=0.0_r8
1169 END DO
1170!
1171 IF (ninnloop.eq.1) THEN
1172 print *,'Illegal configuration!'
1173 print *,'Ninner must be ge 2'
1174 stop
1175 END IF
1176 IF (ninnloop.eq.2) THEN
1177 zu(1)=cg_qg(1,outloop)/cg_delta(1,outloop)
1178!^ tl_zwork(1,3)=tl_zu(1)
1179!^
1180 ad_zu(1)=ad_zu(1)+ad_zwork(1,3)
1181 ad_zwork(1,3)=0.0_r8
1182!^ tl_zu(1)=tl_zrhs(1)/cg_delta(1,outLoop)
1183!^
1184 ad_zrhs(1)=ad_zrhs(1)+ad_zu(1)/cg_delta(1,outloop)
1185 ad_zu(1)=0.0_r8
1186!^ tl_zrhs(1)=tl_cg_QG(1)-tl_cg_delta(1)*zu(1)
1187!^
1188 ad_cg_qg(1)=ad_cg_qg(1)+ad_zrhs(1)
1189 ad_cg_delta(1)=ad_cg_delta(1)-zu(1)*ad_zrhs(1)
1190 ad_zrhs(1)=0.0_r8
1191 ELSE
1192!^ DO ivec=innLoop-2,1,-1
1193!^
1194 DO ivec=1,innloop-2
1195!^ tl_zwork(ivec,3)=tl_zu(ivec)
1196!^
1197 ad_zu(ivec)=ad_zu(ivec)+ad_zwork(ivec,3)
1198 ad_zwork(ivec,3)=0.0_r8
1199
1200!^ tl_zu(ivec)=tl_zu(ivec)-zgam(ivec+1)*tl_zu(ivec+1)
1201!^
1202 ad_zu(ivec+1)=ad_zu(ivec+1)-zgam(ivec+1)*ad_zu(ivec)
1203 END DO
1204!^ tl_zwork(innLoop-1,3)=tl_zu(innLoop-1)
1205!^
1206 ad_zu(innloop-1)=ad_zu(innloop-1)+ad_zwork(innloop-1,3)
1207 ad_zwork(innloop-1,3)=0.0_r8
1208
1209!^ DO ivec=2,innLoop-1
1210!^
1211 DO ivec=innloop-1,2,-1
1212 zbet=cg_delta(ivec,outloop)- &
1213 & cg_beta(ivec,outloop)*zgam(ivec)
1214!^ tl_zu(ivec)=(tl_zrhs(ivec)-cg_beta(ivec,outLoop)* &
1215!^ & tl_zu(ivec-1))/zbet
1216!^
1217 adfac=ad_zu(ivec)/zbet
1218 ad_zu(ivec-1)=ad_zu(ivec-1)- &
1219 & cg_beta(ivec,outloop)*adfac
1220 ad_zrhs(ivec)=ad_zrhs(ivec)+adfac
1221 ad_zu(ivec)=0.0_r8
1222 END DO
1223 zbet=cg_delta(1,outloop)
1224!^ tl_zu(1)=tl_zrhs(1)/zbet
1225!^
1226 ad_zrhs(1)=ad_zrhs(1)+ad_zu(1)/zbet
1227 ad_zu(1)=0.0_r8
1228!
1229!^ tl_zrhs(innLoop-1)=tl_cg_QG(innLoop-1)- &
1230!^ & tl_cg_beta(innLoop-1)*zu(innLoop-2)- &
1231!^ & tl_cg_delta(innLoop-1)*zu(innLoop-1)
1232!^
1233 ad_cg_qg(innloop-1)=ad_cg_qg(innloop-1)+ &
1234 & ad_zrhs(innloop-1)
1235 ad_cg_beta(innloop-1)=ad_cg_beta(innloop-1)- &
1236 & zu(innloop-2)* &
1237 & ad_zrhs(innloop-1)
1238 ad_cg_delta(innloop-1)=ad_cg_delta(innloop-1)- &
1239 & zu(innloop-1)* &
1240 & ad_zrhs(innloop-1)
1241 ad_zrhs(innloop-1)=0.0_r8
1242 DO ivec=2,innloop-2
1243!^ tl_zrhs(ivec)=tl_cg_QG(ivec)- &
1244!^ & tl_cg_beta(ivec)*zu(ivec-1)- &
1245!^ & tl_cg_delta(ivec)*zu(ivec)- &
1246!^ & tl_cg_beta(ivec+1)*zu(ivec+1)
1247!^
1248 ad_cg_qg(ivec)=ad_cg_qg(ivec)+ &
1249 & ad_zrhs(ivec)
1250 ad_cg_beta(ivec)=ad_cg_beta(ivec)- &
1251 & zu(ivec-1)*ad_zrhs(ivec)
1252 ad_cg_delta(ivec)=ad_cg_delta(ivec)- &
1253 & zu(ivec)*ad_zrhs(ivec)
1254 ad_cg_beta(ivec+1)=ad_cg_beta(ivec+1)- &
1255 & zu(ivec+1)*ad_zrhs(ivec)
1256 ad_zrhs(ivec)=0.0_r8
1257 END DO
1258!^ tl_zrhs(1)=tl_cg_QG(1)- &
1259!^ & tl_cg_delta(1)*zu(1)- &
1260!^ & tl_cg_beta(2)*zu(2)
1261!^
1262 ad_cg_qg(1)=ad_cg_qg(1)+ad_zrhs(1)
1263 ad_cg_delta(1)=ad_cg_delta(1)-zu(1)*ad_zrhs(1)
1264 ad_cg_beta(2)=ad_cg_beta(2)-zu(2)*ad_zrhs(1)
1265 ad_zrhs(1)=0.0_r8
1266 END IF
1267
1268 END IF
1269
1270!^ tl_cg_QG(innLoop)=tl_cg_QG(innLoop)+ &
1271!^ & tl_Jb0(outLoop-1)* &
1272!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
1273!^ & vgrad0(Ndatum(ng)+1)+ &
1274!^ & Jb0(outLoop-1)* &
1275!^ & tl_vcglwk(Ndatum(ng)+1,innLoop)* &
1276!^ & vgrad0(Ndatum(ng)+1)+ &
1277!^ & Jb0(outLoop-1)* &
1278!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
1279!^ & tl_vgrad0(Ndatum(ng)+1)
1280!^
1281 ad_jb0(outloop-1)=ad_jb0(outloop-1)+ &
1282 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
1283 & vgrad0(ndatum(ng)+1)* &
1284 & ad_cg_qg(innloop)
1285 ad_vcglwk(ndatum(ng)+1,innloop)= &
1286 & ad_vcglwk(ndatum(ng)+1,innloop)+ &
1287 & jb0(outloop-1)* &
1288 & vgrad0(ndatum(ng)+1)*ad_cg_qg(innloop)
1289 ad_vgrad0(ndatum(ng)+1)=ad_vgrad0(ndatum(ng)+1)+ &
1290 & jb0(outloop-1)* &
1291 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
1292 & ad_cg_qg(innloop)
1293!
1294 DO iobs=1,ndatum(ng)
1295!^ tl_cg_QG(innLoop)=tl_cg_QG(innLoop)+ &
1296!^ & tl_Hbk(iobs)* &
1297!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
1298!^ & vgrad0(iobs)+ &
1299!^ & Hbk(iobs,outLoop)* &
1300!^ & tl_vcglwk(Ndatum(ng)+1,innLoop)* &
1301!^ & vgrad0(iobs)+ &
1302!^ & Hbk(iobs,outLoop)* &
1303!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
1304!^ & tl_vgrad0(iobs)+ &
1305!^ & tl_Hbk(iobs)* &
1306!^ & vcglwk(iobs,innLoop,outLoop)* &
1307!^ & vgrad0(Ndatum(ng)+1)+ &
1308!^ & Hbk(iobs,outLoop)* &
1309!^ & tl_vcglwk(iobs,innLoop)* &
1310!^ & vgrad0(Ndatum(ng)+1)+ &
1311!^ & Hbk(iobs,outLoop)* &
1312!^ & vcglwk(iobs,innLoop,outLoop)* &
1313!^ & tl_vgrad0(Ndatum(ng)+1)
1314!^
1315 ad_hbk(iobs)=ad_hbk(iobs)+ &
1316 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
1317 & vgrad0(iobs)* &
1318 & ad_cg_qg(innloop)
1319 ad_vcglwk(ndatum(ng)+1,innloop)= &
1320 & ad_vcglwk(ndatum(ng)+1,innloop)+ &
1321 & hbk(iobs,outloop)* &
1322 & vgrad0(iobs)* &
1323 & ad_cg_qg(innloop)
1324 ad_vgrad0(iobs)=ad_vgrad0(iobs)+ &
1325 & hbk(iobs,outloop)* &
1326 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
1327 & ad_cg_qg(innloop)
1328 ad_hbk(iobs)=ad_hbk(iobs)+ &
1329 & vcglwk(iobs,innloop,outloop)* &
1330 & vgrad0(ndatum(ng)+1)* &
1331 & ad_cg_qg(innloop)
1332 ad_vcglwk(iobs,innloop)=ad_vcglwk(iobs,innloop)+ &
1333 & hbk(iobs,outloop)* &
1334 & vgrad0(ndatum(ng)+1)* &
1335 & ad_cg_qg(innloop)
1336 ad_vgrad0(ndatum(ng)+1)=ad_vgrad0(ndatum(ng)+1)+ &
1337 & hbk(iobs,outloop)* &
1338 & vcglwk(iobs,innloop,outloop)* &
1339 & ad_cg_qg(innloop)
1340 END DO
1341!
1342! Selime code: cg_QG=zgrad0*HBH'*zcglwk
1343!
1344 DO iobs=1,ndatum(ng)
1345!^ tl_cg_QG(innLoop)=tl_cg_QG(innLoop)+ &
1346!^ & tl_TLmodVal(iobs)* &
1347!^ & vgrad0(iobs)+ &
1348!^ & TLmodVal_S(iobs,innLoop,outLoop)* &
1349!^ & tl_vgrad0(iobs)/ &
1350!^ & cg_beta(innLoop,outLoop)
1351 admodval(iobs)=admodval(iobs)+ &
1352 & vgrad0(iobs)*ad_cg_qg(innloop)
1353 ad_vgrad0(iobs)=ad_vgrad0(iobs)+ &
1354 & tlmodval_s(iobs,innloop,outloop)* &
1355 & ad_cg_qg(innloop)/ &
1356 & cg_beta(innloop,outloop)
1357 END DO
1358!^ tl_cg_QG(innLoop)=0.0_r8
1359!^
1360 ad_cg_qg(innloop)=0.0_r8
1361!
1362! Selime code: t1=t1/beta=TLmodVal
1363!
1364! AUGMENTED
1365!
1366 DO iobs=1,ndatum(ng)
1367!^ tl_TLmodVal(iobs)=tl_TLmodVal(iobs)/ &
1368!^ & cg_beta(innLoop,outLoop)- &
1369!^ & tl_cg_beta(innLoop)* &
1370!^ & TLmodVal_S(iobs,innLoop,outLoop)/ &
1371!^ & (cg_beta(innLoop,outLoop)* &
1372!^ & cg_beta(innLoop,outLoop))
1373!^
1374 ad_cg_beta(innloop)=ad_cg_beta(innloop)- &
1375 & tlmodval_s(iobs,innloop,outloop)* &
1376 & admodval(iobs)/ &
1377 & (cg_beta(innloop,outloop)* &
1378 & cg_beta(innloop,outloop))
1379 admodval(iobs)=admodval(iobs)/cg_beta(innloop,outloop)
1380!^ tl_TLmodVal_S(iobs,innLoop,outLoop)=tl_TLmodVal(iobs)
1381!^
1382 admodval(iobs)=admodval(iobs)+ &
1383 & ad_tlmodval_s(iobs,innloop,outloop)
1384 ad_tlmodval_s(iobs,innloop,outloop)=0.0_r8
1385 END DO
1386!
1387! Selime code: z1=w/beta=vcglwk
1388!
1389 DO iobs=1,ndatum(ng)+1
1390!^ tl_vcglwk(iobs,innLoop)=tl_vcglwk(iobs,innLoop)/ &
1391!^ & cg_beta(innLoop,outLoop)- &
1392!^ & tl_cg_beta(innLoop)* &
1393!^ & vcglwk(iobs,innLoop,outLoop)/ &
1394!^ & cg_beta(innLoop,outLoop)
1395 ad_cg_beta(innloop)=ad_cg_beta(innloop)- &
1396 & vcglwk(iobs,innloop,outloop)* &
1397 & ad_vcglwk(iobs,innloop)/ &
1398 & cg_beta(innloop,outloop)
1399 ad_vcglwk(iobs,innloop)=ad_vcglwk(iobs,innloop)/ &
1400 & cg_beta(innloop,outloop)
1401!^ tl_zcglwk(iobs,innLoop)=tl_pgrad(iobs)/ &
1402!^ & cg_beta(innLoop,outLoop)- &
1403!^ & tl_cg_beta(innLoop)* &
1404!^ & zcglwk(iobs,innLoop,outLoop)/ &
1405!^ & cg_beta(innLoop,outLoop)
1406 ad_cg_beta(innloop)=ad_cg_beta(innloop)- &
1407 & zcglwk(iobs,innloop,outloop)* &
1408 & ad_zcglwk(iobs,innloop)/ &
1409 & cg_beta(innloop,outloop)
1410 ad_pgrad(iobs)=ad_pgrad(iobs)+ &
1411 & ad_zcglwk(iobs,innloop)/ &
1412 & cg_beta(innloop,outloop)
1413 ad_zcglwk(iobs,innloop)=0.0_r8
1414 END DO
1415!
1416! Recompute pgrad.
1417!
1418 DO iobs=1,ndatum(ng)+1
1419 pgrad(iobs)=zcglwk(iobs,innloop,outloop)* &
1420 & cg_beta(innloop,outloop)
1421 END DO
1422!
1423!^ tl_cg_beta(innLoop)=0.5_r8*tl_cg_beta(innLoop)/ &
1424!^ & cg_beta(innLoop,outLoop)
1425!^
1426 ad_cg_beta(innloop)=0.5_r8*ad_cg_beta(innloop)/ &
1427 & cg_beta(innloop,outloop)
1428!^ tl_cg_beta(innLoop)=tl_cg_beta(innLoop)+ &
1429!^ & tl_Jb0(outLoop-1)* &
1430!^ & pgrad(Ndatum(ng)+1)* &
1431!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)+ &
1432!^ & Jb0(outLoop-1)* &
1433!^ & tl_pgrad(Ndatum(ng)+1)* &
1434!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)+ &
1435!^ & Jb0(outLoop-1)* &
1436!^ & pgrad(Ndatum(ng)+1)* &
1437!^ & tl_vcglwk(Ndatum(ng)+1,innLoop)
1438!^
1439 ad_jb0(outloop-1)=ad_jb0(outloop-1)+ &
1440 & pgrad(ndatum(ng)+1)* &
1441 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
1442 & ad_cg_beta(innloop)
1443 ad_pgrad(ndatum(ng)+1)=ad_pgrad(ndatum(ng)+1)+ &
1444 & jb0(outloop-1)* &
1445 & vcglwk(ndatum(ng)+1,innloop,outloop)*&
1446 & ad_cg_beta(innloop)
1447 ad_vcglwk(ndatum(ng)+1,innloop)= &
1448 & ad_vcglwk(ndatum(ng)+1,innloop)+ &
1449 & jb0(outloop-1)* &
1450 & pgrad(ndatum(ng)+1)* &
1451 & ad_cg_beta(innloop)
1452!
1453 DO iobs=1,ndatum(ng)
1454!^ tl_cg_beta(innLoop)=tl_cg_beta(innLoop)+ &
1455!^ & tl_Hbk(iobs)* &
1456!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
1457!^ & pgrad(iobs)+ &
1458!^ & Hbk(iobs,outLoop)* &
1459!^ & tl_vcglwk(Ndatum(ng)+1,innLoop)* &
1460!^ & pgrad(iobs)+ &
1461!^ & Hbk(iobs,outLoop)* &
1462!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)* &
1463!^ & tl_pgrad(iobs)+ &
1464!^ & tl_Hbk(iobs)* &
1465!^ & pgrad(Ndatum(ng)+1)* &
1466!^ & vcglwk(iobs,innLoop,outLoop)+ &
1467!^ & Hbk(iobs,outLoop)* &
1468!^ & tl_pgrad(Ndatum(ng)+1)* &
1469!^ & vcglwk(iobs,innLoop,outLoop)+ &
1470!^ & Hbk(iobs,outLoop)* &
1471!^ & pgrad(Ndatum(ng)+1)* &
1472!^ & tl_vcglwk(iobs,innLoop)
1473!^
1474 ad_hbk(iobs)=ad_hbk(iobs)+ &
1475 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
1476 & pgrad(iobs)* &
1477 & ad_cg_beta(innloop)
1478 ad_vcglwk(ndatum(ng)+1,innloop)= &
1479 & ad_vcglwk(ndatum(ng)+1,innloop)+ &
1480 & hbk(iobs,outloop)* &
1481 & pgrad(iobs)* &
1482 & ad_cg_beta(innloop)
1483 ad_pgrad(iobs)=ad_pgrad(iobs)+ &
1484 & hbk(iobs,outloop)* &
1485 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
1486 & ad_cg_beta(innloop)
1487 ad_hbk(iobs)=ad_hbk(iobs)+ &
1488 & pgrad(ndatum(ng)+1)* &
1489 & vcglwk(iobs,innloop,outloop)* &
1490 & ad_cg_beta(innloop)
1491 ad_pgrad(ndatum(ng)+1)=ad_pgrad(ndatum(ng)+1)+ &
1492 & hbk(iobs,outloop)* &
1493 & vcglwk(iobs,innloop,outloop)* &
1494 & ad_cg_beta(innloop)
1495 ad_vcglwk(iobs,innloop)=ad_vcglwk(iobs,innloop)+ &
1496 & hbk(iobs,outloop)* &
1497 & pgrad(ndatum(ng)+1)* &
1498 & ad_cg_beta(innloop)
1499 END DO
1500!
1501! Selime code: beta=pgrad*HBH'*pgrad=cg_beta
1502!
1503 DO iobs=1,ndatum(ng)
1504!^ tl_cg_beta(innLoop)=tl_cg_beta(innLoop)+ &
1505!^ & tl_pgrad(iobs)* &
1506!^ & TLmodVal_S(iobs,innLoop,outLoop)+ &
1507!^ & pgrad(iobs)* &
1508!^ & tl_TLmodVal(iobs)
1509 ad_pgrad(iobs)=ad_pgrad(iobs)+ &
1510 & tlmodval_s(iobs,innloop,outloop)* &
1511 & ad_cg_beta(innloop)
1512 admodval(iobs)=admodval(iobs)+ &
1513 & pgrad(iobs)*ad_cg_beta(innloop)
1514 END DO
1515!^ tl_cg_beta(innLoop)=0.0_r8
1516!^
1517 ad_cg_beta(innloop)=0.0_r8
1518!
1519 DO iobs=1,ndatum(ng)+1
1520!
1521! Selime code: w=pgrad, t1=TLmodVal.
1522!
1523!^ tl_pgrad(iobs)=tl_zcglwk(iobs,innLoop)
1524!^
1525 ad_zcglwk(iobs,innloop)=ad_zcglwk(iobs,innloop)+ &
1526 & ad_pgrad(iobs)
1527 ad_pgrad(iobs)=0.0_r8
1528 END DO
1529
1530 END IF
1531
1532 END IF minimize
1533!
1534 IF (innloop.gt.0) THEN
1535 DO iobs=1,ndatum(ng)
1536!^ TLmodVal(iobs)=ObsScale(iobs)* &
1537!^ & TLmodVal(iobs)+ &
1538!^ & tl_ObsScale(iobs)* &
1539!^ & TLmodVal_S(iobs,innLoop,outLoop)
1540!^
1541 ad_obsscale(iobs)=ad_obsscale(iobs)+ &
1542 & tlmodval_s(iobs,innloop,outloop)* &
1543 & admodval(iobs)
1544 admodval(iobs)=obsscale(iobs)*admodval(iobs)
1545 END DO
1546 END IF
1547
1548 END IF master_thread
1549
1550# ifdef DISTRIBUTE
1551!
1552!-----------------------------------------------------------------------
1553! Broadcast new solution to other nodes.
1554!-----------------------------------------------------------------------
1555!
1556 CALL mp_bcastf (ng, model, admodval)
1557 CALL mp_bcastf (ng, model, vgrad0(:))
1558 CALL mp_bcastf (ng, model, tlmodval)
1559# ifdef IMPACT_INNER
1560 CALL mp_bcastf (ng, model, ad_obsval(:,:))
1561# else
1562 CALL mp_bcastf (ng, model, ad_obsval)
1563# endif
1564 CALL mp_bcastf (ng, model, ad_obserr)
1565 CALL mp_bcastf (ng, model, ad_obsscale)
1566# if defined RBL4DVAR || \
1567 defined r4dvar || \
1568 defined sensitivity_4dvar
1569 CALL mp_bcastf (ng, model, ad_hbk(:))
1570 CALL mp_bcastf (ng, model, ad_jb0(outloop-1))
1571# endif
1572 CALL mp_bcasti (ng, model, info)
1573 CALL mp_bcastf (ng, model, ad_cg_beta(:))
1574 CALL mp_bcastf (ng, model, ad_cg_qg(:))
1575 CALL mp_bcastf (ng, model, ad_cg_delta(:))
1576 CALL mp_bcastf (ng, model, ad_zgrad0(:))
1577 CALL mp_bcastf (ng, model, ad_zcglwk(:,:))
1578 CALL mp_bcastf (ng, model, ad_vcglwk(:,:))
1579 CALL mp_bcastf (ng, model, fourdvar(ng)%ad_cg_pxsave(:))
1580# endif
1581!
1582 RETURN
1583 END SUBROUTINE ad_rpcg_lanczos
1584
1585#else
1586 SUBROUTINE ad_rpcg_lanczos
1587 RETURN
1588 END SUBROUTINE ad_rpcg_lanczos
1589#endif
distribute_mod::mp_bcastf
Definition distribute.F:75
distribute_mod::mp_bcasti
Definition distribute.F:95
distribute_mod::mp_bcastl
Definition distribute.F:89
distribute_mod
Definition distribute.F:3
mod_fourdvar
Definition mod_fourdvar.F:2
mod_fourdvar::cg_gnorm_v
real(dp), dimension(:), allocatable cg_gnorm_v
Definition mod_fourdvar.F:697
mod_fourdvar::vcglwk
real(r8), dimension(:,:,:), allocatable vcglwk
Definition mod_fourdvar.F:290
mod_fourdvar::ad_obserr
real(r8), dimension(:), allocatable ad_obserr
Definition mod_fourdvar.F:315
mod_fourdvar::fourdvar
type(t_fourdvar), dimension(:), allocatable fourdvar
Definition mod_fourdvar.F:211
mod_fourdvar::cg_beta
real(dp), dimension(:,:), allocatable cg_beta
Definition mod_fourdvar.F:654
mod_fourdvar::ad_obsscale
real(r8), dimension(:), allocatable ad_obsscale
Definition mod_fourdvar.F:316
mod_fourdvar::ad_zcglwk
real(r8), dimension(:,:), allocatable ad_zcglwk
Definition mod_fourdvar.F:306
mod_fourdvar::ndatum
integer, dimension(:), allocatable ndatum
Definition mod_fourdvar.F:409
mod_fourdvar::cg_dla
real(r8), dimension(:,:), allocatable cg_dla
Definition mod_fourdvar.F:300
mod_fourdvar::cg_qg
real(dp), dimension(:,:), allocatable cg_qg
Definition mod_fourdvar.F:706
mod_fourdvar::ad_vgrad0
real(r8), dimension(:), allocatable ad_vgrad0
Definition mod_fourdvar.F:312
mod_fourdvar::obsval
real(r8), dimension(:), allocatable obsval
Definition mod_fourdvar.F:232
mod_fourdvar::ad_obsval
real(r8), dimension(:,:), allocatable ad_obsval
Definition mod_fourdvar.F:322
mod_fourdvar::ad_cg_gnorm_v
real(dp), dimension(:), allocatable ad_cg_gnorm_v
Definition mod_fourdvar.F:686
mod_fourdvar::ad_zgrad0
real(r8), dimension(:), allocatable ad_zgrad0
Definition mod_fourdvar.F:307
mod_fourdvar::ad_jb0
real(r8), dimension(:), allocatable ad_jb0
Definition mod_fourdvar.F:314
mod_fourdvar::obsscale
real(r8), dimension(:), allocatable obsscale
Definition mod_fourdvar.F:227
mod_fourdvar::obserr
real(r8), dimension(:), allocatable obserr
Definition mod_fourdvar.F:225
mod_fourdvar::jb0
real(r8), dimension(:), allocatable jb0
Definition mod_fourdvar.F:291
mod_fourdvar::ad_cg_gnorm
real(dp), dimension(:), allocatable ad_cg_gnorm
Definition mod_fourdvar.F:680
mod_fourdvar::ad_vgrad0s
real(r8), dimension(:), allocatable ad_vgrad0s
Definition mod_fourdvar.F:313
mod_fourdvar::admodval
real(r8), dimension(:), allocatable admodval
Definition mod_fourdvar.F:239
mod_fourdvar::ad_cg_qg
real(dp), dimension(:), allocatable ad_cg_qg
Definition mod_fourdvar.F:679
mod_fourdvar::zcglwk
real(r8), dimension(:,:,:), allocatable zcglwk
Definition mod_fourdvar.F:288
mod_fourdvar::vgrad0
real(r8), dimension(:), allocatable vgrad0
Definition mod_fourdvar.F:295
mod_fourdvar::nlmodval
real(r8), dimension(:), allocatable nlmodval
Definition mod_fourdvar.F:240
mod_fourdvar::cg_delta
real(dp), dimension(:,:), allocatable cg_delta
Definition mod_fourdvar.F:691
mod_fourdvar::ad_cg_beta
real(dp), dimension(:), allocatable ad_cg_beta
Definition mod_fourdvar.F:649
mod_fourdvar::ad_cg_delta
real(dp), dimension(:), allocatable ad_cg_delta
Definition mod_fourdvar.F:678
mod_fourdvar::ad_vcglwk
real(r8), dimension(:,:), allocatable ad_vcglwk
Definition mod_fourdvar.F:311
mod_fourdvar::zgrad0
real(r8), dimension(:,:), allocatable zgrad0
Definition mod_fourdvar.F:294
mod_iounits
Definition mod_iounits.F:2
mod_parallel
Definition mod_parallel.F:2
mod_parallel::master
logical master
Definition mod_parallel.F:40
mod_param
Definition mod_param.F:2
mod_scalars
Definition mod_scalars.F:2
mod_scalars::outer
integer outer
Definition mod_scalars.F:361