ROMS
Loading...
Searching...
No Matches
tl_rpcg_lanczos.F
Go to the documentation of this file.
1#include "cppdefs.h"
2
3#if (defined SENSITIVITY_4DVAR || \
4 defined tl_rbl4dvar || \
5 defined tl_r4dvar) && defined rpcg
6
7 SUBROUTINE tl_rpcg_lanczos (ng, model, outLoop, innLoop, &
8 & NinnLoop, Lcgini)
9!
10!git $Id$
11!================================================== Hernan G. Arango ===
12! Copyright (c) 2002-2025 The ROMS Group Selime Gurol !
13! Licensed under a MIT/X style license Andrew M. Moore !
14! See License_ROMS.md !
15!=======================================================================
16! !
17! Weak Constraint 4-Dimensional Variational (4DVar) Restricted !
18! Pre-conditioned Conjugate Gradient (RPCG) Algorithm !
19
20# if defined R4DVAR_ANA_SENSITIVITY || defined TL_R4DVAR
21! !
22! The indirect representer method solves the system: !
23! !
24! (R_n + Cobs) * Beta_n = h_n !
25! !
26! h_n = Xo - H * X_n !
27! !
28! where R_n is the representer matrix, Cobs is the observation-error !
29! covariance, Beta_n are the representer coefficients, h_n is the !
30! misfit between observations (Xo) and model (H*X_n), and H is the !
31! linearized observation operator. Here, _n denotes iteration. !
32! !
33! This system does not need to be solved explicitly by inverting the !
34! symmetric stabilized representer matrix, P_n: !
35! !
36! P_n = R_n + Cobs !
37! !
38! but by computing the action of P_n on any vector PSI, such that !
39! !
40! P_n * PSI = R_n * PSI + Cobs * PSI !
41! !
42! The representer matrix is not explicitly computed but evaluated by !
43! one integration backward of the adjoint model and one integration !
44! forward of the tangent linear model for any forcing vector PSI. !
45! !
46! Notice that "ObsScale" vector is used for screenning observations. !
47! This scale is one (zero) for good (bad) observations. !
48! !
49! Currently, parallelization of this algorithm is not needed because !
50! each parallel node has a full copy of the assimilation vectors. !
51! !
52! This code solves Ax=b by minimizing the cost function !
53! 0.5*x*A*x-x*b, assuming an initial guess of x=x0. In this case the !
54! gradient is Ax-b and the stabilized representer matrix is A. !
55! !
56! Reference: !
57! !
58! Chua, B. S. and A. F. Bennett, 2001: An inverse ocean modeling !
59! sytem, Ocean Modelling, 3, 137-165. !
60
61# elif defined RBL4DVAR_ANA_SENSITIVITY || \
62 defined rbl4dvar_fct_sensitivity || \
63 defined tl_rbl4dvar
64! !
65! Solve the system (following Courtier, 1997): !
66! !
67! (H M_n B (M_n)' H' + Cobs) * w_n = d_n !
68! !
69! d_n = yo - H * Xb_n !
70! !
71! where M_n is the tangent linear model matrix and _n denotes a !
72! sequence of outer-loop estimates, Cobs is the observation-error !
73! covariance, B is the background error covariance, d_n is the !
74! misfit between observations (yo) and model (H * Xb_n), and H is !
75! the linearized observation operator. !
76! !
77! The analysis increment is: !
78! !
79! dx_n = B M' H' w_n !
80! !
81! so that Xa = Xb + dx_n. !
82! !
83! The system does not need to be solved explicitly by inverting !
84! the symmetric matrix, P_n: !
85! !
86! P_n = H M_n B (M_n)' H' + Cobs !
87! !
88! but by computing the action of P_n on any vector PSI, such that !
89! !
90! P_n * PSI = H M_n B (M_n)' H' * PSI + Cobs * PSI !
91! !
92! The (H M_n B (M_n)' H') matrix is not explicitly computed but !
93! evaluated by one integration backward of the adjoint model and !
94! one integration forward of the tangent linear model for any !
95! forcing vector PSI. !
96! !
97! Reference: !
98! !
99! Courtier, P., 1997: Dual formulation of four-dimensional !
100! variational assimilation, Quart. J. Roy. Meteor. Soc., !
101! 123, 2449-2461. !
102# endif
103! !
104! The restricted preconditioned conjugate gradient (RPCG) algorithm !
105! is used to compute the solution. Specfically, the following system !
106! is solved: !
107! !
108! (Cobs^-1 H M_n B (M_n)' H'+I)*w=Cobs^-1 d_n !
109! !
110! All inner-products are defined with respect to a norm involving !
111! H M_n B (M_n)' H', and the convergence of J in model space is !
112! guaranteed to converge at the same rate as I4D-Var in primal space. !
113! !
114! RPCG Reference: !
115! !
116! Gratton, S. and J. Tshimanga, 2009: An observation-space !
117! formulation of variational assimilation using a restricted !
118! preconditioned conjugate gradient algorithm. QJRMS, 135, !
119! 1573-1585.
120! !
121!=======================================================================
122!
123 USE mod_param
124 USE mod_parallel
125 USE mod_fourdvar
126 USE mod_iounits
127 USE mod_scalars
128
129# ifdef DISTRIBUTE
130!
131 USE distribute_mod, ONLY : mp_bcastf, mp_bcasti, mp_bcastl
132# endif
133 implicit none
134!
135! Imported variable declarations
136!
137 logical, intent(in) :: Lcgini
138 integer, intent(in) :: ng, model, outLoop, innLoop, NinnLoop
139!
140! Local variable declarations.
141!
142 logical :: Ltrans, Laug
143 integer :: i, ic, j, iobs, ivec, Lscale, info
144
145 real(r8) :: zbet, eps, preducv, preducy
146 real(r8) :: Jopt, Jf, Jmod, Jdata, Jb, Jobs, Jact, cff
147 real(r8) :: tl_dla, fact
148
149 real(r8), dimension(NinnLoop) :: zu, zgam, zfact, tl_zfact
150 real(r8), dimension(NinnLoop) :: tl_zu, tl_zrhs
151 real(r8), dimension(Ndatum(ng)+1) :: px, pgrad
152 real(r8), dimension(Ndatum(ng)+1) :: tl_px, tl_pgrad
153 real(r8), dimension(Ninner,3) :: zwork, tl_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 based on the current outer-loop.
178!
179 DO iobs=1,ndatum(ng)+1
180 vgrad0(iobs)=zgrad0(iobs,outloop)
181 END DO
182!
183! Multiply TLmodVal by ObsScale to effectively remove any rejected
184! observations from the analysis.
185!
186 IF (innloop.gt.0) THEN
187 DO iobs=1,ndatum(ng)
188!^ TLmodVal(iobs)=ObsScale(iobs)*TLmodVal(iobs)
189!^
190 tlmodval(iobs)=obsscale(iobs)* &
191 & tlmodval(iobs)+ &
192 & tl_obsscale(iobs)* &
193 & tlmodval_s(iobs,innloop,outloop)
194 END DO
195 END IF
196!
197!-----------------------------------------------------------------------
198! Initialization step, innloop = 0.
199!-----------------------------------------------------------------------
200!
201 minimize : IF (innloop.eq.0) THEN
202
203# if defined RBL4DVAR_ANA_SENSITIVITY || \
204 defined rbl4dvar_fct_sensitivity || \
205 defined tl_rbl4dvar || \
206 defined tl_r4dvar
207!
208 DO iobs=1,ndatum(ng)
209 bckmodval(iobs)=nlmodval(iobs)
210 tl_bckmodval(iobs)=0.0_r8
211!!
212!! THIS IS A TEMPORARY FIX.
213!!
214!!AMM tl_ObsVal(iobs)=ObsVal(iobs)
215!!AMM tl_ObsErr(iobs)=ObsErr(iobs)
216!!
217 END DO
218# endif
219!
220! For the first outer-loop, x0=0.
221!
222 IF (outloop.eq.1) THEN
223 DO iobs=1,ndatum(ng)+1
224!^ px(iobs)=0.0_r8
225!^
226 tl_px(iobs)=0.0_r8
227 END DO
228 DO iobs=1,ndatum(ng)
229!^ Hbk(iobs,outLoop)=0.0_r8
230!^
231 tl_hbk(iobs)=0.0_r8
232 END DO
233 ELSE
234!!
235!! NOTE: CODE WON'T WORK FOR R4D-Var AT THE MOMENT SINCE THIS IS
236!! THE WRONG TLmodVal HERE!!!
237!!
238 DO iobs=1,ndatum(ng)
239!^ Hbk(iobs,outLoop)=-TLmodVal(iobs)
240!^ tl_Hbk(iobs)=-tl_TLmodVal(iobs)
241!^
242 tl_hbk(iobs)=-tlmodval(iobs)
243 END DO
244!^ IF (Lprecond) THEN
245!^ DO iobs=1,Ndatum(ng)+1
246!^ FOURDVAR(ng)%cg_Dpxsum(iobs)=FOURDVAR(ng)%cg_pxsum(iobs)
247!^ END DO
248!^ Lscale=1 ! Spectral LMP
249!^ Laug=.FALSE.
250!^ CALL rprecond (ng, Lscale, Laug, outLoop, NinnLoop-1, &
251!^ & FOURDVAR(ng)%cg_Dpxsum)
252!^ DO iobs=1,Ndatum(ng)+1
253!^ FOURDVAR(ng)%cg_Dpxsum(iobs)= &
254!^ & -FOURDVAR(ng)%cg_Dpxsum(iobs)+
255!^ & FOURDVAR(ng)%cg_pxsum(iobs)
256!^ END DO
257!^ END IF
258 END IF
259!
260! If this is the first pass of the inner loop, set up the vectors and
261! store the first guess. The initial starting guess is assumed to be
262! zero in which case the residual is just: (obs-model).
263!
264 DO iobs=1,ndatum(ng)
265!
266! Selime code: r=b.
267!
268# if defined RBL4DVAR || \
269 defined rbl4dvar_ana_sensitivity || \
270 defined rbl4dvar_fct_sensitivity || \
271 defined tl_rbl4dvar
272!^ pgrad(iobs)=ObsScale(iobs)* &
273!^ & (ObsVal(iobs)-BCKmodVal(iobs))
274!^
275 tl_pgrad(iobs)=obsscale(iobs)* &
276 & (tl_obsval(iobs)- &
277 & tl_bckmodval(iobs))+ &
278 & tl_obsscale(iobs)* &
279 & (obsval(iobs)-bckmodval(iobs))
280# else
281!!
282!! NOTE: THIS CODE IS NOT CORRECT YET FOR R4D-Var
283!!
284!^ pgrad(iobs)=ObsScale(iobs)* &
285!^ & (ObsVal(iobs)-TLmodVal(iobs))
286!^ pgrad(iobs)=ObsScale(iobs)* &
287!^ & (ObsVal(iobs)-TLmodVal_S(iobs,innLoop,outLoop))
288!^ tl_pgrad(iobs)=ObsScale(iobs)* &
289!^ & (tl_ObsVal(iobs)-tl_TLmodVal(iobs))
290!^
291 tl_pgrad(iobs)=obsscale(iobs)* &
292 & (tl_obsval(iobs)- &
293 & tlmodval(iobs))+ &
294 & tl_obsscale(iobs)* &
295 & (obsval(iobs)- &
296 & tlmodval_s(iobs,innloop,outloop))
297# endif
298 IF (obserr(iobs).NE.0.0_r8) THEN
299!^ pgrad(iobs)=pgrad(iobs)/ObsErr(iobs)
300!^
301 pgrad(iobs)=vgrad0(iobs)
302 tl_pgrad(iobs)=tl_pgrad(iobs)/obserr(iobs)- &
303 & tl_obserr(iobs)*pgrad(iobs)/obserr(iobs)
304 END IF
305!^ vgrad0(iobs)=pgrad(iobs)
306!^
307 tl_vgrad0(iobs)=tl_pgrad(iobs)
308!^ vgrad0s(iobs)=vgrad0(iobs)
309!^
310 tl_vgrad0s(iobs)=tl_vgrad0(iobs)
311 END DO
312!
313! AUGMENTED.
314!
315! Augment the residual vector.
316!
317 pgrad(ndatum(ng)+1)=1.0_r8
318 tl_pgrad(ndatum(ng)+1)=0.0_r8
319!^ vgrad0(Ndatum(ng)+1)=pgrad(Ndatum(ng)+1)
320!^
321 tl_vgrad0(ndatum(ng)+1)=tl_pgrad(ndatum(ng)+1)
322!^ vgrad0s(Ndatum(ng)+1)=vgrad0(Ndatum(ng)+1)
323!^
324 tl_vgrad0s(ndatum(ng)+1)=tl_vgrad0(ndatum(ng)+1)
325!
326! If preconditioning, convert pgrad from v-space to y-space.
327!
328! Selime code: z1=G*r.
329!
330!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
331!^ Lscale=1 ! Spectral LMP
332!^ Laug=.TRUE.
333!^ CALL rprecond (ng, Lscale, Laug, outLoop, NinnLoop-1, &
334!^ & pgrad)
335!^ END IF
336!^
337!
338! Selime code: zgrad0=z1.
339!
340! Selime code: vgrad0=r.
341!
342 DO iobs=1,ndatum(ng)+1
343!^ zgrad0(iobs,outLoop)=pgrad(iobs)
344!^
345 tl_zgrad0(iobs)=tl_pgrad(iobs)
346 END DO
347!
348! Selime code: ADmodVal=z1.
349!
350 DO iobs=1,ndatum(ng)
351!^ ADmodVal(iobs)=zgrad0(iobs,outLoop)
352!^ tl_ADmodVal(iobs)=tl_zgrad0(iobs,outLoop)
353!^
354 admodval(iobs)=tl_zgrad0(iobs)
355 END DO
356!
357 preducv=1.0_r8
358 preducy=1.0_r8
359!
360!-----------------------------------------------------------------------
361! Minimization step, innLoop > 0.
362!-----------------------------------------------------------------------
363!
364 ELSE
365!
366! COMPLETE THE CALCULATION OF PREVIOUS LANCZOS VECTOR.
367!
368! Selime code: t1=GDG'*z1=TLmodVal.
369!
370!
371! Complete the computation of the Lanczos vector from previous
372! inner-loop.
373!
374 IF (innloop.EQ.1) THEN
375!^ cg_Gnorm_v(outLoop)=0.0_r8
376!^
377 tl_cg_gnorm_v(outloop)=0.0_r8
378 DO iobs=1,ndatum(ng)
379!^ cg_Gnorm_v(outLoop)=cg_Gnorm_v(outLoop)+ &
380!^ & TLmodVal(iobs)*vgrad0(iobs)
381!^
382 tl_cg_gnorm_v(outloop)=tl_cg_gnorm_v(outloop)+ &
383 & tlmodval(iobs)* &
384 & vgrad0(iobs)+ &
385 & tlmodval_s(iobs,innloop,outloop)* &
386 & tl_vgrad0(iobs)
387 END DO
388!
389! AUGMENTED. Here Hbk=H(x0-xb). Hbk=0 and Jb0=0 when outer=1.
390! Jb0 is computed using sum_i (G'*w)_i, i.e. the sum of the outputs
391! of the adjoint model at the end of AD_LOOP2.
392!
393 DO iobs=1,ndatum(ng)
394!^ cg_Gnorm_v(outLoop)=cg_Gnorm_v(outLoop)+ &
395!^ & Hbk(iobs,outLoop)* &
396!^ & zgrad0(Ndatum(ng)+1,outLoop)* &
397!^ & vgrad0(iobs)+ &
398!^ & Hbk(iobs,outLoop)* &
399!^ & vgrad0(Ndatum(ng)+1)* &
400!^ & zgrad0(iobs,outLoop)
401!^
402 tl_cg_gnorm_v(outloop)=tl_cg_gnorm_v(outloop)+ &
403 & tl_hbk(iobs)* &
404 & zgrad0(ndatum(ng)+1,outloop)* &
405 & vgrad0(iobs)+ &
406 & hbk(iobs,outloop)* &
407 & tl_zgrad0(ndatum(ng)+1)* &
408 & vgrad0(iobs)+ &
409 & hbk(iobs,outloop)* &
410 & zgrad0(ndatum(ng)+1,outloop)* &
411 & tl_vgrad0(iobs)+ &
412 & tl_hbk(iobs)* &
413 & vgrad0(ndatum(ng)+1)* &
414 & zgrad0(iobs,outloop)+ &
415 & hbk(iobs,outloop)* &
416 & tl_vgrad0(ndatum(ng)+1)* &
417 & zgrad0(iobs,outloop)+ &
418 & hbk(iobs,outloop)* &
419 & vgrad0(ndatum(ng)+1)* &
420 & tl_zgrad0(iobs)
421 END DO
422!^ cg_Gnorm_v(outLoop)=cg_Gnorm_v(outLoop)+ &
423!^ & Jb0(outLoop-1)* &
424!^ & vgrad0(Ndatum(ng)+1)* &
425!^ & zgrad0(Ndatum(ng)+1,outLoop)
426!^
427 tl_cg_gnorm_v(outloop)=tl_cg_gnorm_v(outloop)+ &
428 & tl_jb0(outloop-1)* &
429 & vgrad0(ndatum(ng)+1)* &
430 & zgrad0(ndatum(ng)+1,outloop)+ &
431 & jb0(outloop-1)* &
432 & tl_vgrad0(ndatum(ng)+1)* &
433 & zgrad0(ndatum(ng)+1,outloop)+ &
434 & jb0(outloop-1)* &
435 & vgrad0(ndatum(ng)+1)* &
436 & tl_zgrad0(ndatum(ng)+1)
437!
438!^ cg_Gnorm_v(outLoop)=SQRT(cg_Gnorm_v(outLoop))
439!^
440 tl_cg_gnorm_v(outloop)=0.5_r8*tl_cg_gnorm_v(outloop)/ &
441 & cg_gnorm_v(outloop)
442!
443! Selime code: beta_first=sqrt(r't1)=cg_Gnorm_v.
444! cg_QG(1)=beta_first
445!
446!^ cg_QG(1,outLoop)=cg_Gnorm_v(outLoop)
447!^
448 tl_cg_qg(1)=tl_cg_gnorm_v(outloop)
449!
450! Normalize TLmodVal here and save in TLmodVal_S.
451!
452! AUGMENTED.
453!
454 DO iobs=1,ndatum(ng)+1
455!
456! Selime code: v1=r/beta_first=zcglwk(1).
457!
458!^ zcglwk(iobs,1,outLoop)=vgrad0(iobs)/ &
459!^ & cg_Gnorm_v(outLoop)
460!^
461 tl_zcglwk(iobs,1)=tl_vgrad0(iobs)/ &
462 & cg_gnorm_v(outloop)- &
463 & tl_cg_gnorm_v(outloop)* &
464 & zcglwk(iobs,1,outloop)/ &
465 & cg_gnorm_v(outloop)
466!
467! Selime code: z1=z1/beta_first=vcglwk(1).
468!
469!^ vcglwk(iobs,1,outLoop)=zgrad0(iobs,outLoop)/ &
470!^ & cg_Gnorm_v(outLoop)
471!^
472 tl_vcglwk(iobs,1)=tl_zgrad0(iobs)/ &
473 & cg_gnorm_v(outloop)- &
474 & tl_cg_gnorm_v(outloop)* &
475 & vcglwk(iobs,1,outloop)/ &
476 & cg_gnorm_v(outloop)
477 END DO
478!
479! AUGMENTED.
480!
481 DO iobs=1,ndatum(ng)
482!
483! Selime code: t1=t1/beta_first=TLmodVal
484!
485!^ TLmodVal_S(iobs,innLoop,outLoop)=TLmodVal(iobs)
486!^
487 tl_tlmodval_s(iobs,innloop,outloop)=tlmodval(iobs)
488!^ TLmodVal(iobs)=TLmodVal(iobs)/ &
489!^ & cg_Gnorm_v(outLoop)
490!^
491 tlmodval(iobs)=tlmodval(iobs)/ &
492 & cg_gnorm_v(outloop)- &
493 & tl_cg_gnorm_v(outloop)* &
494 & tlmodval_s(iobs,innloop,outloop)/ &
495 & (cg_gnorm_v(outloop)* &
496 & cg_gnorm_v(outloop))
497 END DO
498!
499 ELSE
500!
501! First grab the non-normalized Lanczos vector.
502!
503! AUGMENTED.
504!
505 DO iobs=1,ndatum(ng)+1
506!
507! Selime code: w=pgrad, t1=TLmodVal.
508!
509!
510! Need to multiply zcglwk by cg_beta when recomputing pgrad
511! since zcglwk is the normalized Lanczos vector.
512!
513 pgrad(iobs)=zcglwk(iobs,innloop,outloop)* &
514 & cg_beta(innloop,outloop)
515 tl_pgrad(iobs)=tl_zcglwk(iobs,innloop)
516 END DO
517!
518! Selime's code: beta=pgrad*HBH'*pgrad=cg_beta
519!
520!^ cg_beta(innLoop,outLoop)=0.0_r8
521!^
522 tl_cg_beta(innloop)=0.0_r8
523 DO iobs=1,ndatum(ng)
524!^ cg_beta(innLoop,outLoop)=cg_beta(innLoop,outLoop)+ &
525!^ & pgrad(iobs)* &
526!^ & TLmodVal(iobs)
527!^
528 tl_cg_beta(innloop)=tl_cg_beta(innloop)+ &
529 & tl_pgrad(iobs)* &
530 & tlmodval_s(iobs,innloop,outloop)+ &
531 & pgrad(iobs)* &
532 & tlmodval(iobs)
533 END DO
534!
535! AUGMENTED.
536!
537 DO iobs=1,ndatum(ng)
538!^ cg_beta(innLoop,outLoop)=cg_beta(innLoop,outLoop)+ &
539!^ & Hbk(iobs,outLoop)* &
540!^ & vcglwk(Ndatum(ng)+1,innLoop, &
541!^ & outLoop)* &
542!^ & pgrad(iobs)+ &
543!^ & Hbk(iobs,outLoop)* &
544!^ & pgrad(Ndatum(ng)+1)* &
545!^ & vcglwk(iobs,innLoop,outLoop)
546!^
547 tl_cg_beta(innloop)=tl_cg_beta(innloop)+ &
548 & tl_hbk(iobs)* &
549 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
550 & pgrad(iobs)+ &
551 & hbk(iobs,outloop)* &
552 & tl_vcglwk(ndatum(ng)+1,innloop)* &
553 & pgrad(iobs)+ &
554 & hbk(iobs,outloop)* &
555 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
556 & tl_pgrad(iobs)+ &
557 & tl_hbk(iobs)* &
558 & pgrad(ndatum(ng)+1)* &
559 & vcglwk(iobs,innloop,outloop)+ &
560 & hbk(iobs,outloop)* &
561 & tl_pgrad(ndatum(ng)+1)* &
562 & vcglwk(iobs,innloop,outloop)+ &
563 & hbk(iobs,outloop)* &
564 & pgrad(ndatum(ng)+1)* &
565 & tl_vcglwk(iobs,innloop)
566 END DO
567!^ cg_beta(innLoop,outLoop)=cg_beta(innLoop,outLoop)+ &
568!^ & Jb0(outLoop-1)* &
569!^ & pgrad(Ndatum(ng)+1)* &
570!^ & vcglwk(Ndatum(ng)+1,innLoop, &
571!^ & outLoop)
572!^
573 tl_cg_beta(innloop)=tl_cg_beta(innloop)+ &
574 & tl_jb0(outloop-1)* &
575 & pgrad(ndatum(ng)+1)* &
576 & vcglwk(ndatum(ng)+1,innloop,outloop)+ &
577 & jb0(outloop-1)* &
578 & tl_pgrad(ndatum(ng)+1)* &
579 & vcglwk(ndatum(ng)+1,innloop,outloop)+ &
580 & jb0(outloop-1)* &
581 & pgrad(ndatum(ng)+1)* &
582 & tl_vcglwk(ndatum(ng)+1,innloop)
583!
584!^ cg_beta(innLoop,outLoop)=SQRT(cg_beta(innLoop,outLoop))
585!^
586 tl_cg_beta(innloop)=0.5_r8*tl_cg_beta(innloop)/ &
587 & cg_beta(innloop,outloop)
588!
589! Normalize TLmodVal here and save in TLmodVal_S.
590!
591! AUGMENTED.
592!
593 DO iobs=1,ndatum(ng)+1
594!
595! Selime code: v1=w/beta=zcglwk
596!
597!^ zcglwk(iobs,innLoop,outLoop)=pgrad(iobs)/ &
598!^ & cg_beta(innLoop,outLoop)
599!^
600 tl_zcglwk(iobs,innloop)=tl_pgrad(iobs)/ &
601 & cg_beta(innloop,outloop)- &
602 & tl_cg_beta(innloop)* &
603 & zcglwk(iobs,innloop,outloop)/ &
604 & cg_beta(innloop,outloop)
605!
606! Selime code: z1=w/beta=vcglwk
607!
608!^ vcglwk(iobs,innLoop,outLoop)=vcglwk(iobs,innLoop, &
609!^ & outLoop)/ &
610!^ & cg_beta(innLoop,outLoop)
611!^
612 tl_vcglwk(iobs,innloop)=tl_vcglwk(iobs,innloop)/ &
613 & cg_beta(innloop,outloop)- &
614 & tl_cg_beta(innloop)* &
615 & vcglwk(iobs,innloop,outloop)/ &
616 & cg_beta(innloop,outloop)
617 END DO
618!
619! Selime code: t1=t1/beta=TLmodVal
620!
621! AUGMENTED.
622!
623 DO iobs=1,ndatum(ng)
624!^ TLmodVal_S(iobs,innLoop,outLoop)=TLmodVal(iobs)
625!^
626 tl_tlmodval_s(iobs,innloop,outloop)=tlmodval(iobs)
627!^ TLmodVal(iobs)=TLmodVal(iobs)/ &
628!^ & cg_beta(innLoop,outLoop)
629 tlmodval(iobs)=tlmodval(iobs)/cg_beta(innloop,outloop)- &
630 & tl_cg_beta(innloop)* &
631 & tlmodval_s(iobs,innloop,outloop)/ &
632 & (cg_beta(innloop,outloop)* &
633 & cg_beta(innloop,outloop))
634 END DO
635!
636! Selime code: cg_QG=zgrad0*HBH'*zcglwk
637!
638!^ cg_QG(innLoop,outLoop)=0.0_r8
639!^
640 tl_cg_qg(innloop)=0.0_r8
641 DO iobs=1,ndatum(ng)
642!^ cg_QG(innLoop,outLoop)=cg_QG(innLoop,outLoop)+ &
643!^ & TLmodVal(iobs)* &
644!^ & vgrad0(iobs)
645!^
646 tl_cg_qg(innloop)=tl_cg_qg(innloop)+ &
647 & tlmodval(iobs)* &
648 & vgrad0(iobs)+ &
649 & tlmodval_s(iobs,innloop,outloop)* &
650 & tl_vgrad0(iobs)/ &
651 & cg_beta(innloop,outloop)
652 END DO
653!
654! AUGMENTED.
655!
656 DO iobs=1,ndatum(ng)
657!^ cg_QG(innLoop,outLoop)=cg_QG(innLoop,outLoop)+ &
658!^ & Hbk(iobs,outLoop)* &
659!^ & vcglwk(Ndatum(ng)+1,innLoop, &
660!^ & outLoop)* &
661!^ & vgrad0(iobs)+ &
662!^ & Hbk(iobs,outLoop)* &
663!^ & vcglwk(iobs,innLoop,outLoop)* &
664!^ & vgrad0(Ndatum(ng)+1)
665!^
666 tl_cg_qg(innloop)=tl_cg_qg(innloop)+ &
667 & tl_hbk(iobs)* &
668 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
669 & vgrad0(iobs)+ &
670 & hbk(iobs,outloop)* &
671 & tl_vcglwk(ndatum(ng)+1,innloop)* &
672 & vgrad0(iobs)+ &
673 & hbk(iobs,outloop)* &
674 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
675 & tl_vgrad0(iobs)+ &
676 & tl_hbk(iobs)* &
677 & vcglwk(iobs,innloop,outloop)* &
678 & vgrad0(ndatum(ng)+1)+ &
679 & hbk(iobs,outloop)* &
680 & tl_vcglwk(iobs,innloop)* &
681 & vgrad0(ndatum(ng)+1)+ &
682 & hbk(iobs,outloop)* &
683 & vcglwk(iobs,innloop,outloop)* &
684 & tl_vgrad0(ndatum(ng)+1)
685 END DO
686!^ cg_QG(innLoop,outLoop)=cg_QG(innLoop,outLoop)+ &
687!^ & Jb0(outLoop-1)* &
688!^ & vcglwk(Ndatum(ng)+1,innLoop, &
689!^ & outLoop)* &
690!^ & vgrad0(Ndatum(ng)+1)
691!^
692 tl_cg_qg(innloop)=tl_cg_qg(innloop)+ &
693 & tl_jb0(outloop-1)* &
694 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
695 & vgrad0(ndatum(ng)+1)+ &
696 & jb0(outloop-1)* &
697 & tl_vcglwk(ndatum(ng)+1,innloop)* &
698 & vgrad0(ndatum(ng)+1)+ &
699 & jb0(outloop-1)* &
700 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
701 & tl_vgrad0(ndatum(ng)+1)
702!
703!!
704!! TEST IF (innLoop.gt.1) cg_QG(innLoop,outLoop)=0.0_r8
705!!
706 IF (innloop.eq.ninnloop) THEN
707!
708! Calculate the new solution based upon the new, orthonormalized
709! Lanczos vector. First, the tridiagonal system is solved by
710! decomposition and forward/back substitution.
711!
712 IF (ninnloop.eq.1) THEN
713 print *,'Illegal configuration!'
714 print *,'Ninner must be ge 2'
715 stop
716 END IF
717 IF (ninnloop.eq.2) THEN
718 zu(1)=cg_qg(1,outloop)/cg_delta(1,outloop)
719 tl_zrhs(1)=tl_cg_qg(1)-tl_cg_delta(1)*zu(1)
720 tl_zu(1)=tl_zrhs(1)/cg_delta(1,outloop)
721 zwork(1,3)=zu(1)
722 tl_zwork(1,3)=tl_zu(1)
723 ELSE
724!
725! Compute zu first.
726!
727 zbet=cg_delta(1,outloop)
728 zu(1)=cg_qg(1,outloop)/zbet
729 DO ivec=2,innloop-1
730 zgam(ivec)=cg_beta(ivec,outloop)/zbet
731 zbet=cg_delta(ivec,outloop)- &
732 & cg_beta(ivec,outloop)*zgam(ivec)
733 zu(ivec)=(cg_qg(ivec,outloop)- &
734 & cg_beta(ivec,outloop)* &
735 & zu(ivec-1))/zbet
736 END DO
737 zwork(innloop-1,3)=zu(innloop-1)
738!
739 DO ivec=innloop-2,1,-1
740 zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
741 zwork(ivec,3)=zu(ivec)
742 END DO
743!
744! Now compute tl_zrhs.
745!
746 tl_zrhs(1)=tl_cg_qg(1)- &
747 & tl_cg_delta(1)*zu(1)- &
748 & tl_cg_beta(2)*zu(2)
749 DO ivec=2,innloop-2
750 tl_zrhs(ivec)=tl_cg_qg(ivec)- &
751 & tl_cg_beta(ivec)*zu(ivec-1)- &
752 & tl_cg_delta(ivec)*zu(ivec)- &
753 & tl_cg_beta(ivec+1)*zu(ivec+1)
754 END DO
755 tl_zrhs(innloop-1)=tl_cg_qg(innloop-1)- &
756 & tl_cg_beta(innloop-1)*zu(innloop-2)- &
757 & tl_cg_delta(innloop-1)*zu(innloop-1)
758!
759! Now solve the TL tridiagonal system A*dx=b-dA*x
760!
761 zbet=cg_delta(1,outloop)
762 tl_zu(1)=tl_zrhs(1)/zbet
763 DO ivec=2,innloop-1
764 zgam(ivec)=cg_beta(ivec,outloop)/zbet
765 zbet=cg_delta(ivec,outloop)- &
766 & cg_beta(ivec,outloop)*zgam(ivec)
767 tl_zu(ivec)=(tl_zrhs(ivec)- &
768 & cg_beta(ivec,outloop)* &
769 & tl_zu(ivec-1))/zbet
770 END DO
771 tl_zwork(innloop-1,3)=tl_zu(innloop-1)
772!
773 DO ivec=innloop-2,1,-1
774!^ zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
775!^
776 tl_zu(ivec)=tl_zu(ivec)-zgam(ivec+1)*tl_zu(ivec+1)
777 tl_zwork(ivec,3)=tl_zu(ivec)
778 END DO
779 END IF
780!
781! AUGMENTED.
782!
783 DO iobs=1,ndatum(ng)+1
784 px(iobs)=0.0_r8
785 tl_px(iobs)=0.0_r8
786 DO ivec=1,innloop-1
787 px(iobs)=px(iobs)+ &
788 & zcglwk(iobs,ivec,outloop)*zwork(ivec,3)
789 tl_px(iobs)=tl_px(iobs)+ &
790 & tl_zcglwk(iobs,ivec)* &
791 & zwork(ivec,3)+ &
792 & zcglwk(iobs,ivec,outloop)* &
793 & tl_zwork(ivec,3)
794 END DO
795 END DO
796!
797! If preconditioning, convert px from y-space to v-space.
798! We will always keep px in v-space.
799!
800!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
801!^ Lscale=1 ! Spectral LMP
802!^ Laug=.TRUE.
803!^ CALL rprecond (ng, Lscale, Laug, outLoop, NinnLoop-1, &
804!^ & px)
805!^ END IF
806
807 END IF
808
809 END IF
810!
811! NEXT CALCULATION
812!
813! After the initialization, for every other inner loop, calculate a
814! new Lanczos vector, store it in the matrix, and update search.
815!
816! AUGMENTED.
817!
818 IF (innloop.eq.1) THEN
819 fact=1.0_r8/cg_gnorm_v(outloop)
820 ELSE
821 fact=1.0_r8/cg_beta(innloop,outloop)
822 END IF
823!
824 DO iobs=1,ndatum(ng)
825!
826! Selime code: w=t1.
827!
828!^ pgrad(iobs)=TLmodVal(iobs)+ &
829!^ & Hbk(iobs,outLoop)* &
830!^ & vcglwk(Ndatum(ng)+1,innLoop,outLoop)
831!^
832 pgrad(iobs)=tlmodval_s(iobs,innloop,outloop)* &
833 & fact+ &
834 & hbk(iobs,outloop)* &
835 & vcglwk(ndatum(ng)+1,innloop,outloop)
836 tl_pgrad(iobs)=tlmodval(iobs)+ &
837 & tl_hbk(iobs)* &
838 & vcglwk(ndatum(ng)+1,innloop,outloop)+ &
839 & hbk(iobs,outloop)* &
840 & tl_vcglwk(ndatum(ng)+1,innloop)
841 END DO
842!
843 DO iobs=1,ndatum(ng)
844!
845! Convert gradient contribution from x-space to v-space.
846!
847 IF (obserr(iobs).ne.0.0_r8) THEN
848!
849! Selime code: w=t1/R. No augmented term here since Raug^-1=[R^-1 0].
850!
851 pgrad(iobs)=pgrad(iobs)/obserr(iobs)
852 tl_pgrad(iobs)=tl_pgrad(iobs)/obserr(iobs)- &
853 & tl_obserr(iobs)*pgrad(iobs)/obserr(iobs)
854 END IF
855 END DO
856 pgrad(ndatum(ng)+1)=0.0_r8
857 tl_pgrad(ndatum(ng)+1)=0.0_r8
858!
859! Selime code: w=t1/R+z1
860!
861 DO iobs=1,ndatum(ng)
862 pgrad(iobs)=pgrad(iobs)+ &
863 & vcglwk(iobs,innloop,outloop)
864 tl_pgrad(iobs)=tl_pgrad(iobs)+ &
865 & tl_vcglwk(iobs,innloop)
866 END DO
867 pgrad(ndatum(ng)+1)=pgrad(ndatum(ng)+1)+ &
868 & vcglwk(ndatum(ng)+1,innloop,outloop)
869 tl_pgrad(ndatum(ng)+1)=tl_pgrad(ndatum(ng)+1)+ &
870 & tl_vcglwk(ndatum(ng)+1,innloop)
871!
872 IF (innloop.gt.1) THEN
873!
874! Selime code: w=w-v0*beta
875!
876! AUGMENTED.
877!
878 DO iobs=1,ndatum(ng)+1
879 pgrad(iobs)=pgrad(iobs)- &
880 & cg_beta(innloop,outloop)* &
881 & zcglwk(iobs,innloop-1,outloop)
882 tl_pgrad(iobs)=tl_pgrad(iobs)- &
883 & tl_cg_beta(innloop)* &
884 & zcglwk(iobs,innloop-1,outloop)- &
885 & cg_beta(innloop,outloop)* &
886 & tl_zcglwk(iobs,innloop-1)
887 END DO
888 END IF
889!
890!^ cg_delta(innLoop,outLoop)=0.0_r8
891!^
892 tl_cg_delta(innloop)=0.0_r8
893!
894! Selime code: alpha=w*t1=cg_delta
895!
896! AUGMENTED.
897!
898 IF (innloop.eq.1) THEN
899 fact=1.0_r8/cg_gnorm_v(outloop)
900 ELSE
901 fact=1.0_r8/cg_beta(innloop,outloop)
902 END IF
903!
904 DO iobs=1,ndatum(ng)
905!^ cg_delta(innLoop,outLoop)=cg_delta(innLoop,outLoop)+ &
906!^ & TLmodVal(iobs)* &
907!^ & pgrad(iobs)+ &
908!^ & Hbk(iobs,outLoop)* &
909!^ & vcglwk(Ndatum(ng)+1,innLoop, &
910!^ & outLoop)* &
911!^ & pgrad(iobs)+ &
912!^ & Hbk(iobs,outLoop)* &
913!^ & vcglwk(iobs,innLoop,outLoop)* &
914!^ & pgrad(Ndatum(ng)+1)
915!^
916 tl_cg_delta(innloop)=tl_cg_delta(innloop)+ &
917 & tlmodval(iobs)* &
918 & pgrad(iobs)+ &
919 & tlmodval_s(iobs,innloop,outloop)* &
920 & tl_pgrad(iobs)* &
921 & fact+ &
922 & tl_hbk(iobs)* &
923 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
924 & pgrad(iobs)+ &
925 & hbk(iobs,outloop)* &
926 & tl_vcglwk(ndatum(ng)+1,innloop)* &
927 & pgrad(iobs)+ &
928 & hbk(iobs,outloop)* &
929 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
930 & tl_pgrad(iobs)+ &
931 & tl_hbk(iobs)* &
932 & vcglwk(iobs,innloop,outloop)* &
933 & pgrad(ndatum(ng)+1)+ &
934 & hbk(iobs,outloop)* &
935 & tl_vcglwk(iobs,innloop)* &
936 & pgrad(ndatum(ng)+1)+ &
937 & hbk(iobs,outloop)* &
938 & vcglwk(iobs,innloop,outloop)* &
939 & tl_pgrad(ndatum(ng)+1)
940 END DO
941!^ cg_delta(innLoop,outLoop)=cg_delta(innLoop,outLoop)+ &
942!^ & Jb0(outLoop-1)* &
943!^ & vcglwk(Ndatum(ng)+1,innLoop, &
944!^ & outLoop)* &
945!^ & pgrad(Ndatum(ng)+1)
946!^
947 tl_cg_delta(innloop)=tl_cg_delta(innloop)+ &
948 & tl_jb0(outloop-1)* &
949 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
950 & pgrad(ndatum(ng)+1)+ &
951 & jb0(outloop-1)* &
952 & tl_vcglwk(ndatum(ng)+1,innloop)* &
953 & pgrad(ndatum(ng)+1)+ &
954 & jb0(outloop-1)* &
955 & vcglwk(ndatum(ng)+1,innloop,outloop)* &
956 & tl_pgrad(ndatum(ng)+1)
957!
958! Selime code: w=w-alpha*v1
959!
960! AUGMENTED.
961!
962 DO iobs=1,ndatum(ng)+1
963 pgrad(iobs)=pgrad(iobs)- &
964 & cg_delta(innloop,outloop)* &
965 & zcglwk(iobs,innloop,outloop)
966 tl_pgrad(iobs)=tl_pgrad(iobs)- &
967 & tl_cg_delta(innloop)* &
968 & zcglwk(iobs,innloop,outloop)- &
969 & cg_delta(innloop,outloop)* &
970 & tl_zcglwk(iobs,innloop)
971 END DO
972!
973! Orthonormalize against previous directions.
974!
975! Identify the appropriate Lanczos vector normalization coefficient.
976!
977 DO ivec=1,innloop
978 IF (ivec.eq.1) THEN
979 zfact(ivec)=1.0_r8/cg_gnorm_v(outloop)
980 tl_zfact(ivec)=-tl_cg_gnorm_v(outloop)* &
981 & zfact(ivec)/ &
982 & cg_gnorm_v(outloop)
983 ELSE
984 zfact(ivec)=1.0_r8/cg_beta(ivec,outloop)
985 tl_zfact(ivec)=-tl_cg_beta(ivec)* &
986 & zfact(ivec)/ &
987 & cg_beta(ivec,outloop)
988 ENDIF
989 END DO
990!
991! Selime code: cg_dla=pgrad*HBH'*v1=pgrad*TLmodVal_S/zfact
992!
993! AUGMENTED.
994!
995 DO ivec=innloop,1,-1
996!^ cg_dla(ivec,outLoop)=0.0_r8
997!^
998 tl_dla=0.0_r8
999 DO iobs=1,ndatum(ng)
1000!^ cg_dla(ivec,outLoop)=cg_dla(ivec,outLoop)+ &
1001!^ & pgrad(iobs)* &
1002!^ & TLmodVal_S(iobs,ivec,outLoop)* &
1003!^ & zfact(ivec)+ &
1004!^ & pgrad(iobs)* &
1005!^ & Hbk(iobs,outLoop)* &
1006!^ & vcglwk(Ndatum(ng)+1,ivec,outLoop)+ &
1007!^ & Hbk(iobs,outLoop)* &
1008!^ & vcglwk(iobs,ivec,outLoop)* &
1009!^ & pgrad(Ndatum(ng)+1)
1010!^
1011 tl_dla=tl_dla+ &
1012 & tl_pgrad(iobs)* &
1013 & tlmodval_s(iobs,ivec,outloop)* &
1014 & zfact(ivec)+ &
1015 & pgrad(iobs)* &
1016 & tl_tlmodval_s(iobs,ivec,outloop)* &
1017 & zfact(ivec)+ &
1018 & pgrad(iobs)* &
1019 & tlmodval_s(iobs,ivec,outloop)* &
1020 & tl_zfact(ivec)+ &
1021 & tl_pgrad(iobs)* &
1022 & hbk(iobs,outloop)* &
1023 & vcglwk(ndatum(ng)+1,ivec,outloop)+ &
1024 & pgrad(iobs)* &
1025 & tl_hbk(iobs)* &
1026 & vcglwk(ndatum(ng)+1,ivec,outloop)+ &
1027 & pgrad(iobs)* &
1028 & hbk(iobs,outloop)* &
1029 & tl_vcglwk(ndatum(ng)+1,ivec)+ &
1030 & tl_hbk(iobs)* &
1031 & vcglwk(iobs,ivec,outloop)* &
1032 & pgrad(ndatum(ng)+1)+ &
1033 & hbk(iobs,outloop)* &
1034 & tl_vcglwk(iobs,ivec)* &
1035 & pgrad(ndatum(ng)+1)+ &
1036 & hbk(iobs,outloop)* &
1037 & vcglwk(iobs,ivec,outloop)* &
1038 & tl_pgrad(ndatum(ng)+1)
1039 END DO
1040!^ cg_dla(ivec,outLoop)=cg_dla(ivec,outLoop)+ &
1041!^ & Jb0(outLoop-1)* &
1042!^ & vcglwk(Ndatum(ng)+1,ivec,outLoop)* &
1043!^ & pgrad(Ndatum(ng)+1)
1044!^
1045 tl_dla=tl_dla+ &
1046 & tl_jb0(outloop-1)* &
1047 & vcglwk(ndatum(ng)+1,ivec,outloop)* &
1048 & pgrad(ndatum(ng)+1)+ &
1049 & jb0(outloop-1)* &
1050 & tl_vcglwk(ndatum(ng)+1,ivec)* &
1051 & pgrad(ndatum(ng)+1)+ &
1052 & jb0(outloop-1)* &
1053 & vcglwk(ndatum(ng)+1,ivec,outloop)* &
1054 & tl_pgrad(ndatum(ng)+1)
1055 DO iobs=1,ndatum(ng)+1
1056 pgrad(iobs)=pgrad(iobs)- &
1057 & cg_dla(ivec,outloop)* &
1058 & vcglwk(iobs,ivec,outloop)
1059 tl_pgrad(iobs)=tl_pgrad(iobs)-tl_dla* &
1060 & vcglwk(iobs,ivec,outloop)- &
1061 & cg_dla(ivec,outloop)* &
1062 & tl_vcglwk(iobs,ivec)
1063 END DO
1064 END DO
1065!
1066! Save the non-normalized Lanczos vector. zcglwk is used as temporary
1067! storage. The calculation will be completed at the start of the next
1068! inner-loop.
1069!
1070! Selime code: w=zcglwk.
1071!
1072! AUGMENTED.
1073!
1074 DO iobs=1,ndatum(ng)+1
1075!^ zcglwk(iobs,innLoop+1,outLoop)=pgrad(iobs)
1076!^
1077 tl_zcglwk(iobs,innloop+1)=tl_pgrad(iobs)
1078 END DO
1079!
1080! If preconditioning, convert pgrad from v-space to y-space.
1081!
1082!
1083! Selime code: z1=Gw.
1084!
1085!^ IF (Lprecond.and.(outLoop.gt.1)) THEN
1086!^ Lscale=1 ! Spectral LMP
1087!^ Laug=.TRUE.
1088!^ CALL rprecond (ng, Lscale, Laug, outLoop, NinnLoop-1, pgrad)
1089!^ END IF
1090!
1091! Selime code: z1=vcglwk.
1092!
1093! AUGMENTED.
1094!
1095 DO iobs=1,ndatum(ng)+1
1096!^ vcglwk(iobs,innLoop+1,outLoop)=pgrad(iobs)
1097!^
1098 tl_vcglwk(iobs,innloop+1)=tl_pgrad(iobs)
1099 END DO
1100!
1101! Put the new trial solution into the adjoint vector for the next
1102! loop or put the final solution into the adjoint vector on the
1103! final inner-loop.
1104!
1105 IF (innloop.eq.ninnloop) THEN
1106!
1107! Note: px is already in v-space so there is no need to convert
1108! if preconditioning; cg_pxsave is also in v-space.
1109!
1110 DO iobs=1,ndatum(ng)
1111!^ ADmodVal(iobs)=px(iobs)
1112!^ tl_ADmodVal(iobs)=tl_px(iobs)
1113!^
1114 admodval(iobs)=tl_px(iobs)
1115 END DO
1116 DO iobs=1,ndatum(ng)+1
1117!^ FOURDVAR(ng)%cg_pxsum(iobs)=FOURDVAR(ng)%cg_pxsum(iobs)+ &
1118!^ & FOURDVAR(ng)%cg_pxsave(iobs)
1119!^ FOURDVAR(ng)%cg_pxsave(iobs)=px(iobs)
1120!^
1121 fourdvar(ng)%tl_cg_pxsave(iobs)=tl_px(iobs)
1122!AMM FOURDVAR(ng)%cg_pxsum(iobs)=FOURDVAR(ng)%cg_pxsum(iobs)+ &
1123!AMM & px(iobs)
1124 END DO
1125 ELSE
1126!
1127! Selime code: z1=ADmodVal.
1128!
1129 DO iobs=1,ndatum(ng)
1130!^ ADmodVal(iobs)=vcglwk(iobs,innLoop+1,outLoop)
1131!^ tl_ADmodVal(iobs)=tl_vcglwk(iobs,innLoop+1)
1132!^
1133 admodval(iobs)=tl_vcglwk(iobs,innloop+1)
1134 END DO
1135!
1136 END IF
1137!
1138 END IF minimize
1139
1140 END IF master_thread
1141
1142# ifdef DISTRIBUTE
1143!
1144!-----------------------------------------------------------------------
1145! Broadcast new solution to other nodes.
1146!-----------------------------------------------------------------------
1147!
1148 CALL mp_bcastf (ng, model, admodval)
1149 CALL mp_bcastf (ng, model, vgrad0(:))
1150 CALL mp_bcastf (ng, model, tl_obsval)
1151 CALL mp_bcastf (ng, model, tl_obserr)
1152 CALL mp_bcastf (ng, model, tl_obsscale)
1153# if defined RBL4DVAR || defined R4DVAR || \
1154 defined sensitivity_4dvar || defined tl_rbl4dvar || \
1155 defined tl_r4dvar
1156 CALL mp_bcastf (ng, model, tl_hbk(:))
1157 CALL mp_bcastf (ng, model, tl_jb0(outloop-1))
1158# endif
1159 CALL mp_bcasti (ng, model, info)
1160 CALL mp_bcastf (ng, model, tl_cg_beta(:))
1161 CALL mp_bcastf (ng, model, tl_cg_qg(:))
1162 CALL mp_bcastf (ng, model, tl_cg_delta(:))
1163 CALL mp_bcastf (ng, model, tl_zgrad0(:))
1164 CALL mp_bcastf (ng, model, tl_zcglwk(:,:))
1165 CALL mp_bcastf (ng, model, tl_vcglwk(:,:))
1166 CALL mp_bcastf (ng, model, fourdvar(ng)%tl_cg_pxsave(:))
1167# endif
1168!
1169 RETURN
1170 END SUBROUTINE tl_rpcg_lanczos
1171
1172#else
1173 SUBROUTINE tl_rpcg_lanczos
1174 RETURN
1175 END SUBROUTINE tl_rpcg_lanczos
1176#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::tl_vcglwk
real(r8), dimension(:,:), allocatable tl_vcglwk
Definition mod_fourdvar.F:330
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::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::tl_obsval
real(r8), dimension(:), allocatable tl_obsval
Definition mod_fourdvar.F:340
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::tl_vgrad0s
real(r8), dimension(:), allocatable tl_vgrad0s
Definition mod_fourdvar.F:332
mod_fourdvar::obsval
real(r8), dimension(:), allocatable obsval
Definition mod_fourdvar.F:232
mod_fourdvar::tl_cg_qg
real(dp), dimension(:), allocatable tl_cg_qg
Definition mod_fourdvar.F:683
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::tl_zgrad0
real(r8), dimension(:), allocatable tl_zgrad0
Definition mod_fourdvar.F:328
mod_fourdvar::jb0
real(r8), dimension(:), allocatable jb0
Definition mod_fourdvar.F:291
mod_fourdvar::tl_cg_delta
real(dp), dimension(:), allocatable tl_cg_delta
Definition mod_fourdvar.F:682
mod_fourdvar::tl_zcglwk
real(r8), dimension(:,:), allocatable tl_zcglwk
Definition mod_fourdvar.F:327
mod_fourdvar::tl_cg_beta
real(dp), dimension(:), allocatable tl_cg_beta
Definition mod_fourdvar.F:651
mod_fourdvar::tl_obsscale
real(r8), dimension(:), allocatable tl_obsscale
Definition mod_fourdvar.F:335
mod_fourdvar::admodval
real(r8), dimension(:), allocatable admodval
Definition mod_fourdvar.F:239
mod_fourdvar::zcglwk
real(r8), dimension(:,:,:), allocatable zcglwk
Definition mod_fourdvar.F:288
mod_fourdvar::tl_jb0
real(r8), dimension(:), allocatable tl_jb0
Definition mod_fourdvar.F:333
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::tl_obserr
real(r8), dimension(:), allocatable tl_obserr
Definition mod_fourdvar.F:334
mod_fourdvar::tl_vgrad0
real(r8), dimension(:), allocatable tl_vgrad0
Definition mod_fourdvar.F:331
mod_fourdvar::tl_cg_gnorm_v
real(dp), dimension(:), allocatable tl_cg_gnorm_v
Definition mod_fourdvar.F:685
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
tl_rpcg_lanczos
subroutine tl_rpcg_lanczos(ng, model, outloop, innloop, ninnloop, lcgini)
Definition tl_rpcg_lanczos.F:9