ROMS
Loading...
Searching...
No Matches
ad_congrad.F
Go to the documentation of this file.
1#include "cppdefs.h"
2
3#if defined SENSITIVITY_4DVAR && !defined RPCG
4 SUBROUTINE ad_congrad (ng, model, outLoop, innLoop, NinnLoop, &
5 & Lcgini)
6!
7!git $Id$
8!=================================================== Andrew M. Moore ===
9! Copyright (c) 2002-2025 The ROMS Group Hernan G. Arango !
10! Licensed under a MIT/X style license !
11! See License_ROMS.md !
12!=======================================================================
13! !
14! Weak Constraint 4-Dimensional Variational (4DVar) Pre-conditioned !
15! Conjugate Gradient Algorithm !
16
17# if defined TL_R4DVAR || defined R4DVAR_ANA_SENSITIVITY
18! !
19! The indirect representer method solves the system: !
20! !
21! (R_n + Cobs) * Beta_n = h_n !
22! !
23! h_n = Xo - H * X_n !
24! !
25! where R_n is the representer matrix, Cobs is the observation-error !
26! covariance, Beta_n are the representer coefficients, h_n is the !
27! misfit between observations (Xo) and model (H*X_n), and H is the !
28! linearized observation operator. Here, _n denotes iteration. !
29! !
30! This system does not need to be solved explicitly by inverting the !
31! symmetric stabilized representer matrix, P_n: !
32! !
33! P_n = R_n + Cobs !
34! !
35! but by computing the action of P_n on any vector PSI, such that !
36! !
37! P_n * PSI = R_n * PSI + Cobs * PSI !
38! !
39! The representer matrix is not explicitly computed but evaluated by !
40! one integration backward of the adjoint model and one integration !
41! forward of the tangent linear model for any forcing vector PSI. !
42! !
43! Notice that "ObsScale" vector is used for screenning observations. !
44! This scale is one (zero) for good (bad) observations. !
45! !
46! Currently, parallelization of this algorithm is not needed because !
47! each parallel node has a full copy of the assimilation vectors. !
48! !
49! This code solves Ax=b by minimizing the cost function !
50! 0.5*x*A*x-x*b, assuming an initial guess of x=x0. In this case the !
51! gradient is Ax-b and the Hessian is A. !
52! !
53! Reference: !
54! !
55! Chua, B. S. and A. F. Bennett, 2001: An inverse ocean modeling !
56! sytem, Ocean Modelling, 3, 137-165. !
57
58# elif defined TL_RBL4DVAR || defined RBL4DVAR_ANA_SENSITIVITY
59! !
60! Solve the system (following Courtier, 1997): !
61! !
62! (H M_n B (M_n)' H' + Cobs) * w_n = d_n !
63! !
64! d_n = yo - H * Xb_n !
65! !
66! where M_n is the tangent linear model matrix and _n denotes a !
67! sequence of outer-loop estimates, Cobs is the observation-error !
68! covariance, B is the background error covariance, d_n is the !
69! misfit between observations (yo) and model (H * Xb_n), and H is !
70! the linearized observation operator. !
71! !
72! The analysis increment is: !
73! !
74! dx_n = B M' H' w_n !
75! !
76! so that Xa = Xb + dx_n. !
77! !
78! The system does not need to be solved explicitly by inverting !
79! the symmetric matrix, P_n: !
80! !
81! P_n = H M_n B (M_n)' H' + Cobs !
82! !
83! but by computing the action of P_n on any vector PSI, such that !
84! !
85! P_n * PSI = H M_n B (M_n)' H' * PSI + Cobs * PSI !
86! !
87! The (H M_n B (M_n)' H') matrix is not explicitly computed but !
88! evaluated by one integration backward of the adjoint model and !
89! one integration forward of the tangent linear model for any !
90! forcing vector PSI. !
91! !
92! A preconditioned conjugate gradient algorithm is used to compute !
93! an approximation PSI for w_n. !
94! !
95! Reference: !
96! !
97! Courtier, P., 1997: Dual formulation of four-dimensional !
98! variational assimilation, Quart. J. Roy. Meteor. Soc., !
99! 123, 2449-2461. !
100# endif
101! !
102! Lanczos Algorithm Reference: !
103! !
104! Fisher, M., 1998: Minimization Algorithms for Variational Data !
105! Assimilation. In Seminar on Recent Developments in Numerical !
106! Methods for Atmospheric Modelling, 1998. !
107! !
108! Tchimanga, J., S. Gratton, A.T. Weaver, and A. Sartenaer, 2008: !
109! Limited-memory preconditioners, with application to incremental !
110! four-dimensional variational ocean data assimilation, Q.J.R. !
111! Meteorol. Soc., 134, 753-771. !
112! !
113!=======================================================================
114!
115 USE mod_param
116 USE mod_parallel
117 USE mod_fourdvar
118 USE mod_iounits
119 USE mod_scalars
120
121# ifdef DISTRIBUTE
122!
123 USE distribute_mod, ONLY : mp_bcastf, mp_bcasti, mp_bcastl
124# endif
125 implicit none
126!
127! Imported variable declarations
128!
129 integer, intent(in) :: ng, model, outLoop, innLoop, NinnLoop
130 logical, intent(in) :: Lcgini
131!
132! Local variable declarations.
133!
134 logical :: Ltrans
135
136 integer :: i, j, iobs, ivec, Lscale, info
137
138 real(r8) :: dla, zbet
139 real(r8) :: adfac, ad_dla
140# ifdef MINRES
141 integer :: ii
142 real(r8) :: zsum, zck, zgk
143 real(r8) :: ad_zsum, ad_zck, ad_zgk
144# endif
145
146 real(r8), dimension(NinnLoop) :: zu, zgam
147 real(r8), dimension(NinnLoop) :: ad_zu, ad_zrhs
148 real(r8), dimension(Ndatum(ng)) :: pgrad, zt, pgrad_S
149 real(r8), dimension(Ndatum(ng)) :: ad_px, ad_pgrad, ad_zt
150# ifdef MINRES
151 real(r8), dimension(innLoop,innLoop) :: ztriT, zLT, zLTt
152 real(r8), dimension(innLoop,innLoop) :: ztriTs
153 real(r8), dimension(innLoop,innLoop) :: ad_ztriT, ad_zLT
154 real(r8), dimension(innLoop,innLoop) :: ad_zLTt
155 real(r8), dimension(innLoop) :: tau, zwork1, ze, zeref, zes
156 real(r8), dimension(innLoop) :: ad_tau, ad_zwork1, ad_ze, ad_zeref
157# endif
158!
159!-----------------------------------------------------------------------
160! Weak constraint 4DVar conjugate gradient, Lanczos-based algorithm.
161!-----------------------------------------------------------------------
162!
163! This conjugate gradient algorithm is not run in parallel since the
164! weak constraint is done in observation space. Mostly all the
165! variables are 1D arrays. Therefore, in parallel applications (only
166! distributed-memory is possible) the master node does the computations
167! and then broadcasts results to remaining nodes in the communicator.
168!
169 ltrans=.false.
170
171 master_thread : IF (master) THEN
172!
173! Initialize internal parameters. This is needed here for output
174! reasons.
175!
176!
177! Clear arrays.
178!
179 admodval=0.0_r8
180 ad_pgrad=0.0_r8
181 ad_px=0.0_r8
182 ad_zu=0.0_r8
183 ad_zrhs=0.0_r8
184 ad_zt=0.0_r8
185# ifdef MINRES
186 ad_ze=0.0_r8
187 ad_zsum=0.0_r8
188 ad_zck=0.0_r8
189 ad_zgk=0.0_r8
190 ad_ztrit=0.0_r8
191 ad_zlt=0.0_r8
192 ad_zltt=0.0_r8
193 ad_tau=0.0_r8
194 ad_zwork1=0.0_r8
195 ad_zeref=0.0_r8
196# endif
197!
198! Initialize cg_Gnorm. The adjoint of reprecond is not available.
199!
200 DO i=1,outloop
201 cg_gnorm(i)=cg_gnorm_v(i)
202 END DO
203!
204 IF (innloop.eq.0) THEN
205
206# ifdef RBL4DVAR
207!
208! If this is the first inner-loop, save NLmodVal in BCKmodVal.
209!
210 DO iobs=1,ndatum(ng)
211 bckmodval(iobs)=nlmodval(iobs)
212 END DO
213# endif
214!
215 IF ((outloop.eq.1).or.(.not.lhotstart)) THEN
216!
217 DO iobs=1,ndatum(ng)
218!^ tl_cg_QG(1)=tl_cg_QG(1)+ &
219!^ & tl_zcglwk(iobs,1)*zgrad0(iobs,outLoop)+ &
220!^ & zcglwk(iobs,1,outLoop)*tl_zgrad0(iobs)
221!^
222 ad_zcglwk(iobs,1)=ad_zcglwk(iobs,1)+ &
223 & zgrad0(iobs,outloop)*ad_cg_qg(1)
224 ad_zgrad0(iobs)=ad_zgrad0(iobs)+ &
225 & zcglwk(iobs,1,outloop)*ad_cg_qg(1)
226 END DO
227!^ tl_cg_QG(1)=0.0_r8
228!^
229 ad_cg_qg(1)=0.0_r8
230!
231! Convert ADmodVal from v-space to x-space.
232!
233 DO iobs=1,ndatum(ng)
234 IF (obserr(iobs).NE.0.0_r8) THEN
235!<> tl_ADmodVal(iobs)=tl_ADmodVal(iobs)/SQRT(ObsErr(iobs))
236!<> ad_ADmodVal(iobs)=ad_ADmodVal(iobs)/SQRT(ObsErr(iobs))
237!^
238 tlmodval(iobs)=tlmodval(iobs)/sqrt(obserr(iobs))
239 END IF
240 END DO
241!
242! If preconditioning, convert ADmodVal from y-space to v-space.
243!
244!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
245!^ Lscale=2 ! SQRT spectral LMP
246!^ Ltrans=.FALSE.
247!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, &
248!^ & ADmodVal)
249!^ END IF
250!^
251 DO iobs=1,ndatum(ng)
252!<> tl_ADmodVal(iobs)=tl_zcglwk(iobs,1)
253!<> ad_zcglwk(iobs,1)=ad_zcglwk(iobs,1)+ad_ADmodVal(iobs)
254!^
255 ad_zcglwk(iobs,1)=ad_zcglwk(iobs,1)+tlmodval(iobs)
256!<> ad_ADmodVal(iobs)=0.0_r8
257 tlmodval(iobs)=0.0_r8
258
259!^ tl_zcglwk(iobs,1)=(tl_pgrad(iobs)- &
260!^ & tl_Gnorm(outLoop)* &
261!^ & zcglwk(iobs,1,outLoop))/ &
262!^ & cg_Gnorm(outLoop)
263!^
264 adfac=ad_zcglwk(iobs,1)/cg_gnorm(outloop)
265 ad_cg_gnorm(outloop)=ad_cg_gnorm(outloop)- &
266 & zcglwk(iobs,1,outloop)*adfac
267 ad_pgrad(iobs)=adfac
268 ad_zcglwk(iobs,1)=0.0_r8
269 END DO
270!^ tl_cg_Gnorm(outLoop)=0.5_r8*tl_cg_Gnorm(outLoop)/ &
271!^ & cg_Gnorm(outLoop)
272!^
273 ad_cg_gnorm(outloop)=0.5_r8*ad_cg_gnorm(outloop)/ &
274 & cg_gnorm(outloop)
275
276 DO iobs=1,ndatum(ng)
277!^ tl_cg_Gnorm(outLoop)=tl_cg_Gnorm(outLoop)+ &
278!^ & 2.0_r8*tl_zgrad0(iobs)* &
279!^ & zgrad0(iobs,outLoop)
280!^
281 ad_zgrad0(iobs)=ad_zgrad0(iobs)+ &
282 & 2.0_r8*ad_cg_gnorm(outloop)* &
283 & zgrad0(iobs,outloop)
284!^ tl_zgrad0(iobs)=tl_pgrad(iobs)
285!^
286 ad_pgrad(iobs)=ad_pgrad(iobs)+ad_zgrad0(iobs)
287 ad_zgrad0(iobs)=0.0_r8
288 END DO
289!^ tl_cg_Gnorm(outLoop)=0.0_r8
290!^
291 ad_cg_gnorm(outloop)=0.0_r8
292!
293! If preconditioning, convert pgrad from v-space to y-space.
294!
295!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
296!^ Lscale=2 ! SQRT spectral LMP
297!^ Ltrans=.TRUE.
298!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, &
299!^ & pgrad)
300!^ END IF
301
302 DO iobs=1,ndatum(ng)
303!
304! Convert pgrad from x-space to v-space.
305!
306 IF (obserr(iobs).NE.0.0_r8) THEN
307!^ tl_pgrad(iobs)=tl_pgrad(iobs)/SQRT(ObsErr(iobs))
308!^
309 ad_pgrad(iobs)=ad_pgrad(iobs)/sqrt(obserr(iobs))
310 END IF
311# ifdef RBL4DVAR
312!^ tl_pgrad(iobs)=-ObsScale(iobs)*tl_ObsVal(iobs)
313!^
314 ad_obsval(iobs)=ad_obsval(iobs)-obsscale(iobs)* &
315 & ad_pgrad(iobs)
316 ad_pgrad(iobs)=0.0_r8
317# else
318!<> tl_pgrad(iobs)=-ObsScale(iobs)* &
319!<> & (tl_ObsVal(iobs)-tl_TLmodVal(iobs))
320!<> ad_TLmodVal(iobs)=ad_TLmodVal(iobs)+ &
321!<> & ObsScale(iobs)*ad_pgrad(iobs)
322!^
323 admodval(iobs)=admodval(iobs)+ &
324 & obsscale(iobs)*ad_pgrad(iobs)
325 ad_obsval(iobs)=ad_obsval(iobs)-obsscale(iobs)* &
326 & ad_pgrad(iobs)
327 ad_pgrad(iobs)=0.0_r8
328# endif
329 END DO
330
331# ifdef RBL4DVAR_ANA_SENSITIVITY
332!
333! For observation sensitivity, copy ADmodVal into TLmodVal.
334!
335! DO iobs=1,Ndatum(ng)
336! TLmodVal(iobs)=ADmodVal(iobs)
337! END DO
338# endif
339 ELSE
340
341 IF (lcgini) THEN
342!
343! Convert ADmodVal from v-space to x-space and cg_innov (the
344! contribution to the starting gradient) from x-space to v-space.
345!
346 DO iobs=1,ndatum(ng)
347 IF (obserr(iobs).NE.0.0_r8) THEN
348!^ tl_cg_innov(iobs)=tl_cg_innov(iobs)/SQRT(ObsErr(iobs))
349!^
350 ad_cg_innov(iobs)=ad_cg_innov(iobs)/sqrt(obserr(iobs))
351!<> tl_ADmodVal(iobs)=tl_ADmodVal(iobs)/SQRT(ObsErr(iobs))
352!<> ad_ADmodVal(iobs)=ad_ADmodVal(iobs)/SQRT(ObsErr(iobs))
353!^
354 tlmodval(iobs)=tlmodval(iobs)/sqrt(obserr(iobs))
355 END IF
356 END DO
357!
358! When outer>1 we need to evaluate Ax0 so for inner=0 we use
359! cg_pxsave in v-space as the starting vector.
360!
361 DO iobs=1,ndatum(ng)
362# ifdef RBL4DVAR
363!^ tl_cg_innov(iobs)=-ObsScale(iobs)*tl_ObsVal(iobs)
364!^
365 ad_obsval(iobs)=ad_obsval(iobs)-obsscale(iobs)* &
366 & ad_cg_innov(iobs)
367 ad_cg_innov(iobs)=0.0_r8
368# else
369!<> tl_cg_innov(iobs)=-ObsScale(iobs)* &
370!<> & (tl_ObsVal(iobs)-tl_TLmodVal(iobs))
371!<> ad_TLmodVal(iobs)=ad_TLmodVal(iobs)+ &
372!<> & ObsScale(iobs)*ad_cg_innov(iobs)
373!^
374 admodval(iobs)=admodval(iobs)+ &
375 & obsscale(iobs)*ad_cg_innov(iobs)
376 ad_obsval(iobs)=ad_obsval(iobs)-obsscale(iobs)* &
377 & ad_cg_innov(iobs)
378 ad_cg_innov(iobs)=0.0_r8
379# endif
380!<> tl_ADmodVal(iobs)=tl_cg_pxsave(iobs)
381!<> ad_cg_pxsave(iobs)=ad_cg_pxsave(iobs)+ad_ADmodVal(iobs)
382!^
383 ad_cg_pxsave(iobs)=ad_cg_pxsave(iobs)+tlmodval(iobs)
384!<> ad_ADmodVal(iobs)=0.0_r8
385 tlmodval(iobs)=0.0_r8
386 END DO
387
388 ELSE
389
390 DO iobs=1,ndatum(ng)
391!^ tl_cg_QG(1)=tl_cg_QG(1)+ &
392!^ & tl_zcglwk(iobs,1)*zgrad0(iobs,outLoop)+ &
393!^ & zcglwk(iobs,1,outLoop)*tl_zgrad0(iobs)
394!^
395 ad_zgrad0(iobs)=ad_zgrad0(iobs)+ &
396 & zcglwk(iobs,1,outloop)*ad_cg_qg(1)
397 ad_zcglwk(iobs,1)=ad_zcglwk(iobs,1)+ &
398 & zgrad0(iobs,outloop)*ad_cg_qg(1)
399 END DO
400!^ tl_cg_QG(1)=0.0_r8
401!^
402 ad_cg_qg(1)=0.0_r8
403
404 DO iobs=1,ndatum(ng)
405 IF (obserr(iobs).NE.0.0_r8) THEN
406!<> tl_ADmodVal(iobs)=tl_ADmodVal(iobs)/SQRT(ObsErr(iobs))
407!<> ad_ADmodVal(iobs)=ad_ADmodVal(iobs)/SQRT(ObsErr(iobs))
408!^
409 tlmodval(iobs)=tlmodval(iobs)/sqrt(obserr(iobs))
410 END IF
411 END DO
412!
413! If preconditioning, convert ADmodVal from y-space to v-space.
414!
415!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
416!^ Lscale=2 ! SQRT spectral LMP
417!^ Ltrans=.FALSE.
418!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, &
419!^ & ADmodVal)
420!^ END IF
421!^
422 DO iobs=1,ndatum(ng)
423!<> tl_ADmodVal(iobs)=tl_zcglwk(iobs,1)
424!<> ad_zcglwk(iobs,1)=ad_zcglwk(iobs,1)+ad_ADmodVal(iobs)
425!^
426 ad_zcglwk(iobs,1)=ad_zcglwk(iobs,1)+tlmodval(iobs)
427!<> ad_ADmodVal(iobs)=0.0_r8
428 tlmodval(iobs)=0.0_r8
429!^ tl_zcglwk(iobs,1)=(tl_pgrad(iobs)- &
430!^ & tl_cg_Gnorm(outLoop)* &
431!^ & zcglwk(iobs,1,outLoop))/ &
432!^ & cg_Gnorm(outLoop)
433!^
434 adfac=ad_zcglwk(iobs,1)/cg_gnorm(outloop)
435 ad_cg_gnorm(outloop)=ad_cg_gnorm(outloop)- &
436 & zcglwk(iobs,1,outloop)*adfac
437 ad_pgrad(iobs)=ad_pgrad(iobs)+adfac
438 ad_zcglwk(iobs,1)=0.0_r8
439 END DO
440!^ tl_cg_Gnorm(outLoop)=0.5_r8*tl_cg_Gnorm(outLoop)/ &
441!^ & cg_Gnorm(outLoop)
442!^
443 ad_cg_gnorm(outloop)=0.5_r8*ad_cg_gnorm(outloop)/ &
444 & cg_gnorm(outloop)
445
446 DO iobs=1,ndatum(ng)
447!^ tl_cg_Gnorm(outLoop)=tl_cg_Gnorm(outLoop)+ &
448!^ & 2.0_r8*tl_zgrad0(iobs)* &
449!^ & zgrad0(iobs,outLoop)
450!^
451 ad_zgrad0(iobs)=ad_zgrad0(iobs)+ &
452 & 2.0_r8*zgrad0(iobs,outloop)* &
453 & ad_cg_gnorm(outloop)
454!^ tl_zgrad0(iobs)=tl_pgrad(iobs)
455!^
456 ad_pgrad(iobs)=ad_pgrad(iobs)+ad_zgrad0(iobs)
457 ad_zgrad0(iobs)=0.0_r8
458 END DO
459!^ tl_cg_Gnorm(outLoop)=0.0_r8
460!^
461 ad_cg_gnorm(outloop)=0.0_r8
462!
463! If preconditioning, convert pgrad from v-space to y-space.
464!
465!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
466!^ Lscale=2 ! SQRT spectral LMP
467!^ Ltrans=.TRUE.
468!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, &
469!^ & pgrad)
470!^ END IF
471
472 DO iobs=1,ndatum(ng)
473!
474! Add I*x0=cg_pxsave contribution to the gradient and the term
475! -b=cg_innov (everything is in v-space at this point).
476!
477!^ tl_pgrad(iobs)=tl_pgrad(iobs)+ObsScale(iobs)* &
478!^ & (tl_cg_pxsave(iobs)+tl_cg_innov(iobs))
479!^
480 adfac=obsscale(iobs)*ad_pgrad(iobs)
481 ad_cg_innov(iobs)=ad_cg_innov(iobs)+adfac
482 ad_cg_pxsave(iobs)=ad_cg_pxsave(iobs)+adfac
483!
484! Convert gradient contribution from x-space to v-space.
485!
486 IF (obserr(iobs).NE.0.0_r8) THEN
487!^ tl_pgrad(iobs)=tl_pgrad(iobs)/SQRT(ObsErr(iobs))
488!^
489 ad_pgrad(iobs)=ad_pgrad(iobs)/sqrt(obserr(iobs))
490 END IF
491!<> tl_pgrad(iobs)=ObsScale(iobs)*tl_TLmodVal(iobs)
492!<> ad_TLmodVal(iobs)=ad_TLmodVal(iobs)+ &
493!<> & ObsScale(iobs)*ad_pgrad(iobs)
494!^
495 admodval(iobs)=admodval(iobs)+ &
496 & obsscale(iobs)*ad_pgrad(iobs)
497 ad_pgrad(iobs)=0.0_r8
498 END DO
499 END IF
500 END IF
501
502 ELSE
503!
504! Convert ADmodVal from v-space to x-space.
505!
506 DO iobs=1,ndatum(ng)
507 IF (obserr(iobs).NE.0.0_r8) THEN
508!<> tl_ADmodVal(iobs)=tl_ADmodVal(iobs)/SQRT(ObsErr(iobs))
509!<> ad_ADmodVal(iobs)=ad_ADmodVal(iobs)/SQRT(ObsErr(iobs))
510!^
511 tlmodval(iobs)=tlmodval(iobs)/sqrt(obserr(iobs))
512 END IF
513 END DO
514
515 IF (innloop.eq.ninnloop) THEN
516
517 IF (lhotstart) THEN
518 DO iobs=1,ndatum(ng)
519!<> tl_cg_pxsave(iobs)=tl_ADmodVal(iobs)
520!<> ad_ADmodVal(iobs)=ad_ADmodVal(iobs)+ad_cg_pxsave(iobs)
521!^
522 tlmodval(iobs)=tlmodval(iobs)+ad_cg_pxsave(iobs)
523 ad_cg_pxsave(iobs)=0.0_r8
524 END DO
525 DO iobs=1,ndatum(ng)
526!<> tl_ADmodVal(iobs)=tl_ADmodVal(iobs)+tl_cg_pxsave(iobs)
527!<> ad_cg_pxsave(iobs)=ad_ADmodVal(iobs)
528!^
529 ad_cg_pxsave(iobs)=tlmodval(iobs)
530 END DO
531 END IF
532!
533! Note: px is already in v-space so there is no need to convert
534! if preconditioning.
535!
536 DO iobs=1,ndatum(ng)
537!<> tl_ADmodVal(iobs)=tl_px(iobs)
538!<> ad_px(iobs)=ad_px(iobs)+ad_ADmodVal(iobs)
539!^
540 ad_px(iobs)=ad_px(iobs)+tlmodval(iobs)
541!<> ad_ADmodVal(iobs)=0.0_r8
542 tlmodval(iobs)=0.0_r8
543 END DO
544 ELSE
545!
546! If preconditioning, convert ADmodVal from y-space to v-space.
547!
548!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
549!^ Lscale=2 ! SQRT spectral LMP
550!^ Ltrans=.FALSE.
551!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, &
552!^ & ADmodVal)
553!^ END IF
554!^
555 DO iobs=1,ndatum(ng)
556!<> tl_ADmodVal(iobs)=tl_zcglwk(iobs,innLoop+1)
557!<> ad_zcglwk(iobs,innLoop+1)=ad_zcglwk(iobs,innLoop+1)+ &
558!<> & ad_ADmodVal(iobs)
559!^
560 ad_zcglwk(iobs,innloop+1)=ad_zcglwk(iobs,innloop+1)+ &
561 & tlmodval(iobs)
562!<> ad_ADmodVal(iobs)=0.0_r8
563 tlmodval(iobs)=0.0_r8
564 END DO
565 END IF
566
567!
568! Calculate the new solution based upon the new, orthonormalized
569! Lanczos vector. First, the tridiagonal system is solved by
570! decomposition and forward/back substitution.
571!
572 IF (innloop.eq.ninnloop) THEN
573
574# ifdef MINRES
575!
576! Compute zu first.
577!
578! Perform a LQ factorization of the tridiagonal matrix.
579!
580 ztrit=0.0_r8
581 DO i=1,innloop
582 ztrit(i,i)=cg_delta(i,outloop)
583 END DO
584 DO i=1,innloop-1
585 ztrit(i,i+1)=cg_beta(i+1,outloop)
586 END DO
587 DO i=2,innloop
588 ztrit(i,i-1)=cg_beta(i,outloop)
589 END DO
590!
591! Save a copy of ztriT since it is overwritten by sqlq.
592!
593 DO j=1,innloop
594 DO i=1,innloop
595 ztrits(i,j)=ztrit(i,j)
596 END DO
597 END DO
598!
599 CALL sqlq(innloop, ztrit, tau, zwork1)
600!
601! Isolate L=zLT and its transpose.
602!
603 zlt=0.0_r8
604 zltt=0.0_r8
605 DO i=1,innloop
606 DO j=1,i
607 zlt(i,j)=ztrit(i,j)
608 END DO
609 END DO
610 DO j=1,innloop
611 DO i=1,innloop
612 zltt(i,j)=zlt(j,i)
613 END DO
614 END DO
615!
616! Now form ze=beta1*Q*e1.
617!
618 ze=0.0_r8
619 DO i=1,innloop
620 ze(i)=-cg_qg(i,outloop)
621 END DO
622 DO i=1,innloop
623 DO j=1,innloop
624 zeref(j)=0.0_r8
625 END DO
626 zeref(i)=1.0_r8
627 DO j=i+1,innloop
628 zeref(j)=ztrit(i,j)
629 END DO
630 zsum=0.0_r8
631 DO j=1,innloop
632 zsum=zsum+ze(j)*zeref(j)
633 END DO
634 DO j=1,innloop
635 ze(j)=ze(j)-tau(i)*zsum*zeref(j)
636 END DO
637 END DO
638!
639! Save ze for use later.
640!
641 DO i=1,innloop
642 zes(i)=ze(i)
643 END DO
644!
645! Now form ze=D*ze (using equation 5.6 and 6.5 also).
646!
647 zgk=sqrt(zlt(innloop,innloop)*zlt(innloop,innloop)+ &
648 & cg_beta(innloop+1,outloop)*cg_beta(innloop+1,outloop))
649 zck=zlt(innloop,innloop)/zgk
650 ze(innloop)=zck*ze(innloop)
651!
652! Now compute ze=inv(L')*ze.
653!
654 DO j=innloop,1,-1
655 ze(j)=ze(j)/zltt(j,j)
656 DO i=1,j-1
657 ze(i)=ze(i)-ze(j)*zltt(i,j)
658 END DO
659 END DO
660!
661! Copy the solution ze into zu.
662!
663 DO i=1,innloop
664 zu(i)=ze(i)
665 END DO
666# else
667!
668! Compute zu and zgam first.
669!
670 zbet=cg_delta(1,outloop)
671 zu(1)=-cg_qg(1,outloop)/zbet
672 DO ivec=2,innloop
673 zgam(ivec)=cg_beta(ivec,outloop)/zbet
674 zbet=cg_delta(ivec,outloop)- &
675 & cg_beta(ivec,outloop)*zgam(ivec)
676 zu(ivec)=(-cg_qg(ivec,outloop)-cg_beta(ivec,outloop)* &
677 & zu(ivec-1))/zbet
678 END DO
679 DO ivec=innloop-1,1,-1
680 zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
681 END DO
682# endif
683 DO iobs=1,ndatum(ng)
684 DO ivec=1,innloop
685!^ tl_px(iobs)=tl_px(iobs)+ &
686!^ & tl_zcglwk(iobs,ivec)*zu(ivec)+ &
687!^ & zcglwk(iobs,ivec,outLoop)*tl_zu(ivec)
688!^
689 ad_zu(ivec)=ad_zu(ivec)+ &
690 & zcglwk(iobs,ivec,outloop)*ad_px(iobs)
691 ad_zcglwk(iobs,ivec)=ad_zcglwk(iobs,ivec)+ &
692 & zu(ivec)*ad_px(iobs)
693 END DO
694!^ tl_px(iobs)=0.0_r8
695!^
696 ad_px(iobs)=0.0_r8
697 END DO
698
699# ifdef MINRES
700 DO i=1,innloop
701!^ tl_zu(i)=tl_ze(i)
702!^
703 ad_ze(i)=ad_ze(i)+ad_zu(i)
704 ad_zu(i)=0.0_r8
705 END DO
706!
707! Adjoint solver the linear triangular system.
708!
709 DO j=1,innloop
710 DO i=j-1,1,-1
711!^ tl_ze(i)=tl_ze(i)-tl_ze(j)*zLTt(i,j)
712!^
713 ad_ze(j)=ad_ze(j)-ad_ze(i)*zltt(i,j)
714 END DO
715!^ tl_ze(j)=tl_ze(j)/zLTt(j,j)
716!^
717 ad_ze(j)=ad_ze(j)/zltt(j,j)
718 END DO
719!
720 DO i=1,innloop
721!^ tl_ze(i)=tl_ze(i)-tl_zrhs(i)
722!^
723 ad_zrhs(i)=ad_zrhs(i)-ad_ze(i)
724 DO j=1,innloop
725!^ tl_zrhs(i)=tl_zrhs(i)+tl_zLTt(i,j)*ze(j)
726!^
727 ad_zltt(i,j)=ad_zltt(i,j)+ze(j)*ad_zrhs(i)
728 END DO
729!^ tl_zrhs(i)=0.0_r8
730!^
731 ad_zrhs(i)=0.0_r8
732 END DO
733!
734! Use saved values zes of ze here.
735!
736!^ tl_ze(innLoop)=tl_zck*ze(innLoop)+zck*tl_ze(innLoop)
737!^
738 ad_zck=ad_zck+zes(innloop)*ad_ze(innloop)
739 ad_ze(innloop)=zck*ad_ze(innloop)
740!^ tl_zck=tl_zLT(innLoop,innLoop)/zgk-tl_zgk*zck/zgk
741!^
742 adfac=ad_zck/zgk
743 ad_zlt(innloop,innloop)=ad_zlt(innloop,innloop)+adfac
744 ad_zgk=ad_zgk-zck*adfac
745 ad_zck=0.0_r8
746 IF (zgk.GT.0.0_r8) THEN
747!^ tl_zgk=(tl_zLT(innLoop,innLoop)*zLT(innLoop,innLoop)+ &
748!^ & tl_cg_beta(innLoop+1)*cg_beta(innLoop+1,outLoop))/&
749!^ & zgk
750!^
751 adfac=ad_zgk/zgk
752 ad_zlt(innloop,innloop)=ad_zlt(innloop,innloop)+ &
753 & zlt(innloop,innloop)*adfac
754 ad_cg_beta(innloop+1)=ad_cg_beta(innloop+1)+ &
755 & cg_beta(innloop+1,outloop)*adfac
756 ad_zgk=0.0_r8
757 ELSE
758!^ tl_zgk=0.0_r8
759!^
760 ad_zgk=0.0_r8
761 ENDIF
762!
763! This is the adjoint multiplication by Q via the H matrices so must
764! be done in reverse.
765!
766 DO i=innloop,1,-1
767!
768! Recompute intermediate arrays.
769!
770 DO j=1,innloop
771 ze(j)=-cg_qg(j,outloop)
772 END DO
773 DO ii=1,i
774 DO j=1,innloop
775 zeref(j)=0.0_r8
776 END DO
777 zeref(ii)=1.0_r8
778 DO j=ii+1,innloop
779 zeref(j)=ztrit(ii,j)
780 END DO
781 zsum=0.0_r8
782 DO j=1,innloop
783 zsum=zsum+ze(j)*zeref(j)
784 END DO
785!
786! Save ze again.
787!
788 DO j=1,innloop
789 zes(j)=ze(j)
790 END DO
791 DO j=1,innloop
792 ze(j)=ze(j)-tau(ii)*zsum*zeref(j)
793 END DO
794 END DO
795!
796 DO j=1,innloop
797!^ tl_ze(j)=tl_ze(j)-tl_tau(i)*zsum*zeref(j)- &
798!^ & tau(i)*tl_zsum*zeref(j)- &
799!^ & tau(i)*zsum*tl_zeref(j)
800!^
801 ad_tau(i)=ad_tau(i)-zsum*zeref(j)*ad_ze(j)
802 ad_zsum=ad_zsum-tau(i)*zeref(j)*ad_ze(j)
803 ad_zeref(j)=ad_zeref(j)-tau(i)*zsum*ad_ze(j)
804 END DO
805 DO j=1,innloop
806!^ tl_zsum=tl_zsum+tl_ze(j)*zeref(j)+ze(j)*tl_zeref(j)
807!^
808 ad_ze(j)=ad_ze(j)+zeref(j)*ad_zsum
809 ad_zeref(j)=ad_zeref(j)+zes(j)*ad_zsum
810 END DO
811!^ tl_zsum=0.0_r8
812!^
813 ad_zsum=0.0_r8
814 DO j=i+1,innloop
815!^ tl_zeref(j)=tl_ztriT(i,j)
816 ad_ztrit(i,j)=ad_ztrit(i,j)+ad_zeref(j)
817 ad_zeref(j)=0.0_r8
818 END DO
819!^ tl_zeref(i)=0.0_r8
820!^
821 ad_zeref(i)=0.0_r8
822 DO j=1,innloop
823!^ tl_zeref(j)=0.0_r8
824!^
825 ad_zeref(j)=0.0_r8
826 END DO
827 END DO
828 DO i=1,innloop
829!^ tl_ze(i)=-tl_cg_QG(i,outLoop)
830!^
831 ad_cg_qg(i)=ad_cg_qg(i)-ad_ze(i)
832 ad_ze(i)=0.0_r8
833 END DO
834!
835 DO j=1,innloop
836 DO i=1,innloop
837!^ tl_zLTt(i,j)=tl_zLT(j,i)
838!^
839 ad_zlt(j,i)=ad_zlt(j,i)+ad_zltt(i,j)
840 ad_zltt(i,j)=0.0_r8
841 END DO
842 END DO
843 DO i=1,innloop
844 DO j=1,i
845!^ tl_zLT(i,j)=tl_ztriT(i,j)
846!^
847 ad_ztrit(i,j)=ad_ztrit(i,j)+ad_zlt(i,j)
848 ad_zlt(i,j)=0.0_r8
849 END DO
850 END DO
851!^ tl_zLT=0.0_r8
852!^
853 ad_zlt=0.0_r8
854!^ tl_zLTt=0.0_r8
855!^
856 ad_zltt=0.0_r8
857!
858! Grab a clean copy of ztriT.
859!
860 DO j=1,innloop
861 DO i=1,innloop
862 ztrit(i,j)=ztrits(i,j)
863 END DO
864 END DO
865!
866 CALL ad_sqlq (innloop, ztrit, ad_ztrit, tau, ad_tau, &
867 & zwork1, ad_zwork1)
868!
869 DO i=2,innloop
870!^ tl_ztriT(i,i-1)=tl_cg_beta(i)
871!^
872 ad_cg_beta(i)=ad_cg_beta(i)+ad_ztrit(i,i-1)
873 ad_ztrit(i,i-1)=0.0_r8
874 END DO
875 DO i=1,innloop-1
876!^ tl_ztriT(i,i+1)=tl_cg_beta(i+1)
877!^
878 ad_cg_beta(i+1)=ad_cg_beta(i+1)+ad_ztrit(i,i+1)
879 ad_ztrit(i,i+1)=0.0_r8
880 END DO
881 DO i=1,innloop
882!^ tl_ztriT(i,i)=tl_cg_delta(i)
883!^
884 ad_cg_delta(i)=ad_cg_delta(i)+ad_ztrit(i,i)
885 ad_ztrit(i,i)=0.0_r8
886 END DO
887!^ tl_ztriT=0.0_r8
888!^
889 ad_ztrit=0.0_r8
890# else
891 IF (ninnloop.eq.1) THEN
892 zu(1)=-cg_qg(1,outloop)/cg_delta(1,outloop)
893!^ tl_zu(1)=tl_zrhs(1)/cg_delta(1,outLoop)
894!^
895 ad_zrhs(1)=ad_zrhs(1)+ad_zu(1)/cg_delta(1,outloop)
896 ad_zu(1)=0.0_r8
897!^ tl_zrhs(1)=-tl_cg_QG(1)-tl_cg_delta(1)*zu(1)
898!^
899 ad_cg_qg(1)=ad_cg_qg(1)-ad_zrhs(1)
900 ad_cg_delta(1)=ad_cg_delta(1)-zu(1)*ad_zrhs(1)
901 ad_zrhs(1)=0.0_r8
902 ELSE
903!^ DO ivec=innLoop-1,1,-1
904!^
905 DO ivec=1,innloop-1
906!^ tl_zu(ivec)=tl_zu(ivec)-zgam(ivec+1)*tl_zu(ivec+1)
907!^
908 ad_zu(ivec+1)=ad_zu(ivec+1)-zgam(ivec+1)*ad_zu(ivec)
909 END DO
910!^ DO ivec=2,innLoop
911!^
912 DO ivec=innloop,2,-1
913 zbet=cg_delta(ivec,outloop)- &
914 & cg_beta(ivec,outloop)*zgam(ivec)
915!^ tl_zu(ivec)=(tl_zrhs(ivec)-cg_beta(ivec,outLoop)* &
916!^ & tl_zu(ivec-1))/zbet
917!^
918 adfac=ad_zu(ivec)/zbet
919 ad_zu(ivec-1)=ad_zu(ivec-1)- &
920 & cg_beta(ivec,outloop)*adfac
921 ad_zrhs(ivec)=ad_zrhs(ivec)+adfac
922 ad_zu(ivec)=0.0_r8
923 END DO
924 zbet=cg_delta(1,outloop)
925!^ tl_zu(1)=tl_zrhs(1)/zbet
926!^
927 ad_zrhs(1)=ad_zrhs(1)+ad_zu(1)/zbet
928 ad_zu(1)=0.0_r8
929!^ tl_zrhs(innLoop)=-tl_cg_QG(innLoop)- &
930!^ & tl_cg_beta(innLoop)*zu(innLoop-1)- &
931!^ & tl_cg_delta(innLoop)*zu(innLoop)
932!^
933 ad_cg_qg(innloop)=ad_cg_qg(innloop)-ad_zrhs(innloop)
934 ad_cg_beta(innloop)=ad_cg_beta(innloop)- &
935 & zu(innloop-1)*ad_zrhs(innloop)
936 ad_cg_delta(innloop)=ad_cg_delta(innloop)- &
937 & zu(innloop)*ad_zrhs(innloop)
938 ad_zrhs(innloop)=0.0_r8
939
940 DO ivec=2,innloop-1
941!^ tl_zrhs(ivec)=-tl_cg_QG(ivec)- &
942!^ & tl_cg_beta(ivec)*zu(ivec-1)- &
943!^ & tl_cg_delta(ivec)*zu(ivec)- &
944!^ & tl_cg_beta(ivec+1)*zu(ivec+1)
945!^
946 ad_cg_qg(ivec)=ad_cg_qg(ivec)-ad_zrhs(ivec)
947 ad_cg_beta(ivec)=ad_cg_beta(ivec)- &
948 & zu(ivec-1)*ad_zrhs(ivec)
949 ad_cg_delta(ivec)=ad_cg_delta(ivec)- &
950 & zu(ivec)*ad_zrhs(ivec)
951 ad_cg_beta(ivec+1)=ad_cg_beta(ivec+1)- &
952 & zu(ivec+1)*ad_zrhs(ivec)
953 ad_zrhs(ivec)=0.0_r8
954 END DO
955!^ tl_zrhs(1)=-tl_cg_QG(1)- &
956!^ & tl_cg_delta(1)*zu(1)- &
957!^ & tl_cg_beta(2)*zu(2)
958!^
959 ad_cg_qg(1)=ad_cg_qg(1)-ad_zrhs(1)
960 ad_cg_delta(1)=ad_cg_delta(1)-zu(1)*ad_zrhs(1)
961 ad_cg_beta(2)=ad_cg_beta(2)-zu(2)*ad_zrhs(1)
962 ad_zrhs(1)=0.0_r8
963 END IF
964# endif
965 END IF
966
967 DO iobs=1,ndatum(ng)
968!^ tl_cg_QG(innLoop+1)=tl_cg_QG(innLoop+1)+ &
969!^ & tl_zcglwk(iobs,innLoop+1)* &
970!^ & zgrad0(iobs,outLoop)+ &
971!^ & zcglwk(iobs,innLoop+1,outLoop)* &
972!^ & tl_zgrad0(iobs)
973!^
974 ad_zgrad0(iobs)=ad_zgrad0(iobs)+ &
975 & zcglwk(iobs,innloop+1,outloop)* &
976 & ad_cg_qg(innloop+1)
977 ad_zcglwk(iobs,innloop+1)=ad_zcglwk(iobs,innloop+1)+ &
978 & zgrad0(iobs,outloop)* &
979 & ad_cg_qg(innloop+1)
980 END DO
981!^ tl_cg_QG(innLoop+1)=0.0_r8
982!^
983 ad_cg_qg(innloop+1)=0.0_r8
984 DO iobs=1,ndatum(ng)
985!^ tl_zcglwk(iobs,innLoop+1)=(tl_pgrad(iobs)- &
986!^ & tl_cg_beta(innLoop+1)* &
987!^ & zcglwk(iobs,innLoop+1,outLoop))/ &
988!^ & cg_beta(innLoop+1,outLoop)
989!^
990 adfac=ad_zcglwk(iobs,innloop+1)/cg_beta(innloop+1,outloop)
991 ad_cg_beta(innloop+1)=ad_cg_beta(innloop+1)- &
992 & zcglwk(iobs,innloop+1,outloop)*adfac
993 ad_pgrad(iobs)=ad_pgrad(iobs)+adfac
994 ad_zcglwk(iobs,innloop+1)=0.0_r8
995 END DO
996!^ tl_cg_beta(innLoop+1)=0.5_r8*tl_cg_beta(innLoop+1)/ &
997!^ & cg_beta(innLoop+1,outLoop)
998!^
999 ad_cg_beta(innloop+1)=0.5_r8*ad_cg_beta(innloop+1)/ &
1000 & cg_beta(innloop+1,outloop)
1001!
1002! Compute appropriate value of pgrad.
1003!
1004 DO iobs=1,ndatum(ng)
1005 pgrad(iobs)=obsscale(iobs)*tlmodval_s(iobs,innloop,outloop)
1006!
1007! Convert gradient contribution from x-space to v-space.
1008!
1009 IF (obserr(iobs).NE.0.0_r8) THEN
1010 pgrad(iobs)=pgrad(iobs)/sqrt(obserr(iobs))
1011 END IF
1012 END DO
1013 DO iobs=1,ndatum(ng)
1014 zt(iobs)=zcglwk(iobs,innloop,outloop)
1015 END DO
1016!
1017! If preconditioning, convert zt from y-space to v-space.
1018!
1019!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
1020!^ Lscale=2 ! SQRT spectral LMP
1021!^ Ltrans=.FALSE.
1022!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, zt)
1023!^ END IF
1024!^
1025 DO iobs=1,ndatum(ng)
1026 pgrad(iobs)=pgrad(iobs)+obsscale(iobs)*zt(iobs)
1027 END DO
1028!
1029! If preconditioning, convert pgrad from v-space to y-space.
1030!
1031!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
1032!^ Lscale=2 ! SQRT spectral LMP
1033!^ Ltrans=.TRUE.
1034!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, pgrad)
1035!^ END IF
1036!^
1037!
1038! Compute the new Lanczos vector.
1039!
1040 DO iobs=1,ndatum(ng)
1041 pgrad(iobs)=pgrad(iobs)- &
1042 & cg_delta(innloop,outloop)* &
1043 & zcglwk(iobs,innloop,outloop)
1044 END DO
1045 IF (innloop.gt.1) THEN
1046 DO iobs=1,ndatum(ng)
1047 pgrad(iobs)=pgrad(iobs)- &
1048 & cg_beta(innloop,outloop)* &
1049 & zcglwk(iobs,innloop-1,outloop)
1050 END DO
1051 END IF
1052!
1053! Orthonormalize against previous Lanczos vectors.
1054!
1055 DO ivec=innloop,1,-1
1056 DO iobs=1,ndatum(ng)
1057 pgrad(iobs)=pgrad(iobs)- &
1058 & cg_dla(ivec,outloop)* &
1059 & zcglwk(iobs,ivec,outloop)
1060 END DO
1061 END DO
1062
1063 DO iobs=1,ndatum(ng)
1064!^ tl_cg_beta(innLoop+1)=tl_cg_beta(innLoop+1)+ &
1065!^ & 2.0_r8*tl_pgrad(iobs)*pgrad(iobs)
1066!^
1067 ad_pgrad(iobs)=ad_pgrad(iobs)+ &
1068 & 2.0_r8*pgrad(iobs)*ad_cg_beta(innloop+1)
1069 END DO
1070!^ tl_cg_beta(innLoop+1)=0.0_r8
1071!^
1072 ad_cg_beta(innloop+1)=0.0_r8
1073!
1074! Compute appropriate value of pgrad.
1075!
1076 DO iobs=1,ndatum(ng)
1077 pgrad(iobs)=obsscale(iobs)*tlmodval_s(iobs,innloop,outloop)
1078!
1079! Convert gradient contribution from x-space to v-space.
1080!
1081 IF (obserr(iobs).NE.0.0_r8) THEN
1082 pgrad(iobs)=pgrad(iobs)/sqrt(obserr(iobs))
1083 END IF
1084 END DO
1085 DO iobs=1,ndatum(ng)
1086 zt(iobs)=zcglwk(iobs,innloop,outloop)
1087 END DO
1088!
1089! If preconditioning, convert zt from y-space to v-space.
1090!
1091!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
1092!^ Lscale=2 ! SQRT spectral LMP
1093!^ Ltrans=.FALSE.
1094!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, zt)
1095!^ END IF
1096!^
1097 DO iobs=1,ndatum(ng)
1098 pgrad(iobs)=pgrad(iobs)+obsscale(iobs)*zt(iobs)
1099 END DO
1100!
1101! If preconditioning, convert pgrad from v-space to y-space.
1102!
1103!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
1104!^ Lscale=2 ! SQRT spectral LMP
1105!^ Ltrans=.TRUE.
1106!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, pgrad)
1107!^ END IF
1108!^
1109!
1110! Compute the new Lanczos vector.
1111!
1112 DO iobs=1,ndatum(ng)
1113 pgrad(iobs)=pgrad(iobs)- &
1114 & cg_delta(innloop,outloop)* &
1115 & zcglwk(iobs,innloop,outloop)
1116 END DO
1117 IF (innloop.gt.1) THEN
1118 DO iobs=1,ndatum(ng)
1119 pgrad(iobs)=pgrad(iobs)- &
1120 & cg_beta(innloop,outloop)* &
1121 & zcglwk(iobs,innloop-1,outloop)
1122 END DO
1123 END IF
1124!
1125! Save the intermediate value of pgrad.
1126!
1127 DO iobs=1,ndatum(ng)
1128 pgrad_s(iobs)=pgrad(iobs)
1129 END DO
1130!
1131! Orthonormalize against previous Lanczos vectors.
1132!
1133!^ DO ivec=innLoop,1,-1
1134!^
1135 DO ivec=1,innloop
1136
1137 DO iobs=1,ndatum(ng)
1138 pgrad(iobs)=pgrad_s(iobs)
1139 END DO
1140 DO i=ivec+1,innloop
1141 DO iobs=1,ndatum(ng)
1142 pgrad(iobs)=pgrad(iobs)- &
1143 & cg_dla(i,outloop)* &
1144 & zcglwk(iobs,i,outloop)
1145 END DO
1146 END DO
1147
1148 ad_dla=0.0_r8
1149 DO iobs=1,ndatum(ng)
1150!^ tl_pgrad(iobs)=tl_pgrad(iobs)- &
1151!^ & cg_dla(ivec,outLoop)*tl_zcglwk(iobs,ivec)- &
1152!^ & tl_dla*zcglwk(iobs,ivec,outLoop)
1153!^
1154 ad_zcglwk(iobs,ivec)=ad_zcglwk(iobs,ivec)- &
1155 & cg_dla(ivec,outloop)*ad_pgrad(iobs)
1156 ad_dla=ad_dla-zcglwk(iobs,ivec,outloop)*ad_pgrad(iobs)
1157 END DO
1158 DO iobs=1,ndatum(ng)
1159!^ tl_dla=tl_dla+ &
1160!^ & tl_pgrad(iobs)*zcglwk(iobs,ivec,outLoop)+ &
1161!^ & pgrad(iobs)*tl_zcglwk(iobs,ivec)
1162!^
1163 ad_zcglwk(iobs,ivec)=ad_zcglwk(iobs,ivec)+ &
1164 & pgrad(iobs)*ad_dla
1165 ad_pgrad(iobs)=ad_pgrad(iobs)+ &
1166 & zcglwk(iobs,ivec,outloop)*ad_dla
1167 END DO
1168!^ tl_dla=0.0_r8
1169!^
1170 ad_dla=0.0_r8
1171 END DO
1172
1173 IF (innloop.gt.1) THEN
1174 DO iobs=1,ndatum(ng)
1175!^ tl_pgrad(iobs)=tl_pgrad(iobs)- &
1176!^ & tl_cg_beta(innLoop)* &
1177!^ & zcglwk(iobs,innLoop-1,outLoop)- &
1178!^ & cg_beta(innLoop,outLoop)* &
1179!^ & tl_zcglwk(iobs,innLoop-1)
1180!^
1181 ad_zcglwk(iobs,innloop-1)=ad_zcglwk(iobs,innloop-1)- &
1182 & cg_beta(innloop,outloop)* &
1183 & ad_pgrad(iobs)
1184 ad_cg_beta(innloop)=ad_cg_beta(innloop)- &
1185 & zcglwk(iobs,innloop-1,outloop)* &
1186 & ad_pgrad(iobs)
1187 END DO
1188 END IF
1189 DO iobs=1,ndatum(ng)
1190!^ tl_pgrad(iobs)=tl_pgrad(iobs)- &
1191!^ & tl_cg_delta(innLoop)* &
1192!^ & zcglwk(iobs,innLoop,outLoop)- &
1193!^ & cg_delta(innLoop,outLoop)* &
1194!^ & tl_zcglwk(iobs,innLoop)
1195!^
1196 ad_zcglwk(iobs,innloop)=ad_zcglwk(iobs,innloop)- &
1197 & cg_delta(innloop,outloop)* &
1198 & ad_pgrad(iobs)
1199 ad_cg_delta(innloop)=ad_cg_delta(innloop)- &
1200 & zcglwk(iobs,innloop,outloop)* &
1201 & ad_pgrad(iobs)
1202 END DO
1203!
1204! Compute appropriate value of pgrad.
1205!
1206 DO iobs=1,ndatum(ng)
1207 pgrad(iobs)=obsscale(iobs)*tlmodval_s(iobs,innloop,outloop)
1208!
1209! Convert gradient contribution from x-space to v-space.
1210!
1211 IF (obserr(iobs).NE.0.0_r8) THEN
1212 pgrad(iobs)=pgrad(iobs)/sqrt(obserr(iobs))
1213 END IF
1214 END DO
1215 DO iobs=1,ndatum(ng)
1216 zt(iobs)=zcglwk(iobs,innloop,outloop)
1217 END DO
1218!
1219! If preconditioning, convert zt from y-space to v-space.
1220!
1221!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
1222!^ Lscale=2 ! SQRT spectral LMP
1223!^ Ltrans=.FALSE.
1224!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, zt)
1225!^ END IF
1226!^
1227 DO iobs=1,ndatum(ng)
1228 pgrad(iobs)=pgrad(iobs)+obsscale(iobs)*zt(iobs)
1229 END DO
1230!
1231! If preconditioning, convert pgrad from v-space to y-space.
1232!
1233!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
1234!^ Lscale=2 ! SQRT spectral LMP
1235!^ Ltrans=.TRUE.
1236!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, pgrad)
1237!^ END IF
1238!^
1239 DO iobs=1,ndatum(ng)
1240!^ tl_cg_delta(innLoop)=tl_cg_delta(innLoop)+ &
1241!^ & tl_zcglwk(iobs,innLoop)*pgrad(iobs)+ &
1242!^ & zcglwk(iobs,innLoop,outLoop)* &
1243!^ & tl_pgrad(iobs)
1244!^
1245 ad_pgrad(iobs)=ad_pgrad(iobs)+ &
1246 & zcglwk(iobs,innloop,outloop)* &
1247 & ad_cg_delta(innloop)
1248 ad_zcglwk(iobs,innloop)=ad_zcglwk(iobs,innloop)+ &
1249 & pgrad(iobs)*ad_cg_delta(innloop)
1250 END DO
1251!^ tl_cg_delta(innLoop)=0.0_r8
1252!^
1253 ad_cg_delta(innloop)=0.0_r8
1254!
1255! If preconditioning, convert pgrad from v-space to y-space.
1256!
1257!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
1258!^ Lscale=2 ! SQRT spectral LMP
1259!^ Ltrans=.TRUE.
1260!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, pgrad)
1261!^ END IF
1262!^
1263 DO iobs=1,ndatum(ng)
1264!^ tl_pgrad(iobs)=tl_pgrad(iobs)+ObsScale(iobs)*tl_zt(iobs)
1265!^
1266 ad_zt(iobs)=ad_zt(iobs)+obsscale(iobs)*ad_pgrad(iobs)
1267 END DO
1268!
1269! If preconditioning, convert zt from y-space to v-space.
1270!
1271!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
1272!^ Lscale=2 ! SQRT spectral LMP
1273!^ Ltrans=.FALSE.
1274!^ CALL rprecond (ng, Lscale, Ltrans, outLoop, NinnLoop, zt)
1275!^ END IF
1276!^
1277 DO iobs=1,ndatum(ng)
1278!^ tl_zt(iobs)=tl_zcglwk(iobs,innLoop)
1279!^
1280 ad_zcglwk(iobs,innloop)=ad_zcglwk(iobs,innloop)+ad_zt(iobs)
1281 ad_zt(iobs)=0.0_r8
1282 END DO
1283 DO iobs=1,ndatum(ng)
1284 IF (obserr(iobs).NE.0.0_r8) THEN
1285!^ tl_pgrad(iobs)=tl_pgrad(iobs)/SQRT(ObsErr(iobs))
1286!^
1287 ad_pgrad(iobs)=ad_pgrad(iobs)/sqrt(obserr(iobs))
1288 END IF
1289!^ tl_pgrad(iobs)=ObsScale(iobs)*tl_TLmodVal(iobs)
1290!<> ad_TLmodVal(iobs)=ad_TLmodVal(iobs)+ObsScale(iobs)* &
1291!<> & ad_pgrad(iobs)
1292!^
1293 admodval(iobs)=admodval(iobs)+obsscale(iobs)*ad_pgrad(iobs)
1294 ad_pgrad(iobs)=0.0_r8
1295 END DO
1296
1297 END IF
1298 END IF master_thread
1299
1300# ifdef DISTRIBUTE
1301!
1302! Broadcast new solution to other nodes.
1303!
1304 CALL mp_bcasti (ng, model, exit_flag)
1305 CALL mp_bcastf (ng, model, admodval)
1306 CALL mp_bcastf (ng, model, ad_obsval)
1307 CALL mp_bcastf (ng, model, tlmodval)
1308 CALL mp_bcasti (ng, model, info)
1309 CALL mp_bcastf (ng, model, ad_cg_beta(:))
1310 CALL mp_bcastf (ng, model, ad_cg_delta(:))
1311 CALL mp_bcastf (ng, model, ad_zcglwk(:,:))
1312 CALL mp_bcastf (ng, model, ad_cg_qg(:))
1313 CALL mp_bcastf (ng, model, ad_cg_gnorm(:))
1314 CALL mp_bcastf (ng, model, ad_cg_pxsave(:))
1315 CALL mp_bcastf (ng, model, ad_cg_innov(:))
1316 CALL mp_bcastf (ng, model, ad_zgrad0(:))
1317# endif
1318
1319 RETURN
1320 END SUBROUTINE ad_congrad
1321#else
1322 SUBROUTINE ad_congrad
1323 RETURN
1324 END SUBROUTINE ad_congrad
1325#endif
1326
ad_congrad
subroutine ad_congrad
Definition ad_congrad.F:1323
ad_sqlq
subroutine ad_sqlq(innloop, a, ad_a, tau, ad_tau, y, ad_y)
Definition ad_sqlq.F:7
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::cg_beta
real(dp), dimension(:,:), allocatable cg_beta
Definition mod_fourdvar.F:654
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_obsval
real(r8), dimension(:,:), allocatable ad_obsval
Definition mod_fourdvar.F:322
mod_fourdvar::ad_zgrad0
real(r8), dimension(:), allocatable ad_zgrad0
Definition mod_fourdvar.F:307
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::lhotstart
logical lhotstart
Definition mod_fourdvar.F:817
mod_fourdvar::cg_gnorm
real(dp), dimension(:), allocatable cg_gnorm
Definition mod_fourdvar.F:696
mod_fourdvar::ad_cg_gnorm
real(dp), dimension(:), allocatable ad_cg_gnorm
Definition mod_fourdvar.F:680
mod_fourdvar::ad_cg_pxsave
real(r8), dimension(:), allocatable ad_cg_pxsave
Definition mod_fourdvar.F:319
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::nlmodval
real(r8), dimension(:), allocatable nlmodval
Definition mod_fourdvar.F:240
mod_fourdvar::ad_cg_innov
real(r8), dimension(:), allocatable ad_cg_innov
Definition mod_fourdvar.F:318
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::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::exit_flag
integer exit_flag
Definition mod_scalars.F:559