ROMS
Loading...
Searching...
No Matches
congrad.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE congrad_mod
3
4#if defined WEAK_CONSTRAINT && \
5 !(defined RPCG || \
6 defined r_symmetry || \
7 defined sp4dvar)
8!
9!git $Id$
10!================================================== Hernan G. Arango ===
11! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
12! Licensed under a MIT/X style license !
13! See License_ROMS.md !
14!=======================================================================
15! !
16! Weak Constraint 4-Dimensional Variational (4D-Var) Pre-conditioned !
17! Conjugate Gradient Algorithm !
18
19# if defined R4DVAR || defined R4DVAR_ANA_SENSITIVITY || \
20 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 Hessian 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 || defined RBL4DVAR_ANA_SENSITIVITY || \
62 defined tl_rbl4dvar
63! !
64! Solve the system (following Courtier, 1997): !
65! !
66! (H M_n B (M_n)' H' + Cobs) * w_n = d_n !
67! !
68! d_n = yo - H * Xb_n !
69! !
70! where M_n is the tangent linear model matrix and _n denotes a !
71! sequence of outer-loop estimates, Cobs is the observation-error !
72! covariance, B is the background error covariance, d_n is the !
73! misfit between observations (yo) and model (H * Xb_n), and H is !
74! the linearized observation operator. !
75! !
76! The analysis increment is: !
77! !
78! dx_n = B M' H' w_n !
79! !
80! so that Xa = Xb + dx_n. !
81! !
82! The system does not need to be solved explicitly by inverting !
83! the symmetric matrix, P_n: !
84! !
85! P_n = H M_n B (M_n)' H' + Cobs !
86! !
87! but by computing the action of P_n on any vector PSI, such that !
88! !
89! P_n * PSI = H M_n B (M_n)' H' * PSI + Cobs * PSI !
90! !
91! The (H M_n B (M_n)' H') matrix is not explicitly computed but !
92! evaluated by one integration backward of the adjoint model and !
93! one integration forward of the tangent linear model for any !
94! forcing vector PSI. !
95! !
96! A preconditioned conjugate gradient algorithm is used to compute !
97! an approximation PSI for w_n. !
98! !
99! Reference: !
100! !
101! Courtier, P., 1997: Dual formulation of four-dimensional !
102! variational assimilation, Quart. J. Roy. Meteor. Soc., !
103! 123, 2449-2461. !
104# endif
105! !
106! Lanczos Algorithm Reference: !
107! !
108! Fisher, M., 1998: Minimization Algorithms for Variational Data !
109! Assimilation. In Seminar on Recent Developments in Numerical !
110! Methods for Atmospheric Modelling, 1998. !
111! !
112! Tchimanga, J., S. Gratton, A.T. Weaver, and A. Sartenaer, 2008: !
113! Limited-memory preconditioners, with application to incremental !
114! four-dimensional variational ocean data assimilation, Q.J.R. !
115! Meteorol. Soc., 134, 753-771. !
116! !
117!=======================================================================
118!
119 USE mod_param
120 USE mod_parallel
121 USE mod_fourdvar
122 Use mod_iounits
123 USE mod_ncparam
124 USE mod_scalars
125!
126# ifdef DISTRIBUTE
127 USE distribute_mod, ONLY : mp_bcastf, mp_bcasti, mp_bcastl
128# endif
129 USE lapack_mod, ONLY : dsteqr
130 USE strings_mod, ONLY : founderror
131!
132 implicit none
133!
134! Interface for mixed precision minimization.
135!
136 INTERFACE rprecond
137# ifdef SINGLE_PRECISION
138 MODULE PROCEDURE rprecond_dp
139# endif
140 MODULE PROCEDURE rprecond_r8
141 END INTERFACE rprecond
142!
143 PUBLIC :: congrad
144 PUBLIC :: cg_read_congrad
145!
146 PRIVATE :: cg_read_congrad_nf90
147# if defined PIO_LIBRARY && defined DISTRIBUTE
148 PRIVATE :: cg_read_congrad_pio
149# endif
150 PRIVATE :: cg_write_congrad
151 PRIVATE :: cg_write_congrad_nf90
152# if defined PIO_LIBRARY && defined DISTRIBUTE
153 PRIVATE :: cg_write_congrad_pio
154# endif
155!
156 PRIVATE
157!
158 CONTAINS
159!
160!***********************************************************************
161 SUBROUTINE congrad (ng, model, outLoop, innLoop, NinnLoop, &
162 & Lcgini)
163!***********************************************************************
164!
165! Imported variable declarations
166!
167 logical, intent(in) :: lcgini
168!
169 integer, intent(in) :: ng, model, outloop, innloop, ninnloop
170!
171! Local variable declarations.
172!
173 logical :: ltrans
174!
175 integer :: i, ic, j, iobs, ivec, lscale, info
176!
177 real(dp) :: zbet, eps, preducv, preducy
178# ifdef MINRES
179 real(dp) :: zsum, zck, zgk
180# endif
181 real(dp) :: jopt, jf, jmod, jdata, jb, jobs, jact, cff
182
183 real(dp), dimension(NinnLoop) :: zu, zgam
184 real(dp), dimension(NinnLoop) :: ztemp1, ztemp2, zu1, zu2
185 real(dp), dimension(Ndatum(ng)) :: px, pgrad, zw, zt
186 real(dp), dimension(Ndatum(ng)) :: zhv, zht, zd
187 real(dp), dimension(Ninner,3) :: zwork
188 real(dp), dimension(2*(NinnLoop-1)) :: work
189 real(dp), dimension(Ninner,Ninner) :: zgv
190# ifdef MINRES
191 real(dp), dimension(innLoop,innLoop) :: ztrit, zlt, zltt
192 real(dp), dimension(innLoop) :: tau, zwork1, ze, zeref
193# endif
194!
195 character (len=13) :: string
196
197 character (len=*), parameter :: myfile = &
198 & __FILE__
199!
200 sourcefile=myfile
201!
202!=======================================================================
203! Weak constraint 4D-Var conjugate gradient, Lanczos-based algorithm.
204!=======================================================================
205
206# ifdef PROFILE
207!
208! Turn on time clock.
209!
210 CALL wclock_on (ng, model, 85, __line__, myfile)
211# endif
212!
213! This conjugate gradient algorithm is not run in parallel since the
214! weak constraint is done in observation space. Mostly all the
215! variables are 1D arrays. Therefore, in parallel applications (only
216! distributed-memory is possible) the master node does the computations
217! and then broadcasts results to remaining nodes in the communicator.
218!
219! This version of congrad solves A(x+x0)=b for x, by minimizing
220! J=0.5*x'Ax-x'(b+Ax0), where x0 is a first-guess estimate of the
221! representer coefficients from the previous outer-loop.
222! For the first outer-loop, x0=0. In the code x0=cg_pxsave and
223! x=px.
224!
225 ritzmaxerr=hevecerr
226 ltrans=.false.
227!
228 master_thread : IF (master) THEN
229!
230! Initialize internal parameters. This is needed here for output
231! reasons.
232!
233 jopt=0.0_dp
234 jf=0.0_dp
235 jdata=0.0_dp
236 jmod=0.0_dp
237 jb=0.0_dp
238 jobs=0.0_dp
239 jact=0.0_dp
240 eps=0.0_dp
241 preducv=0.0_dp
242 preducy=0.0_dp
243!
244!-----------------------------------------------------------------------
245! Initialization step, innloop = 0.
246!-----------------------------------------------------------------------
247!
248 minimize : IF (innloop.eq.0) THEN
249
250# if defined RBL4DVAR || \
251 defined rbl4dvar_ana_sensitivity || \
252 defined tl_rbl4dvar
253!
254! If this is the first inner-loop, save NLmodVal in BCKmodVal.
255!
256 DO iobs=1,ndatum(ng)
257 bckmodval(iobs)=nlmodval(iobs)
258 END DO
259# endif
260!
261! If this is the first outer-loop, clear the solution vector px.
262!
263 IF ((outloop.eq.1).or.(.not.lhotstart)) THEN
264!
265! For the first outer-loop, x0=0.
266!
267 DO iobs=1,ndatum(ng)
268 px(iobs)=0.0_dp
269 cg_pxsave(iobs)=0.0_dp
270 END DO
271!
272! If this is the first pass of the inner loop, set up the vectors and
273! store the first guess. The initial starting guess is assumed to be
274! zero in which case the gradient is just: -(obs-model).
275! A first-level preconditioning is applied using the inverse of the
276! observation error standard deviations (i.e. sqrt(ObsErr)).
277!
278 DO iobs=1,ndatum(ng)
279# if defined RBL4DVAR || \
280 defined rbl4dvar_ana_sensitivity || \
281 defined tl_rbl4dvar
282 pgrad(iobs)=-obsscale(iobs)* &
283 & (obsval(iobs)-bckmodval(iobs))
284# else
285 pgrad(iobs)=-obsscale(iobs)* &
286 & (obsval(iobs)-tlmodval(iobs))
287# endif
288!
289! Convert pgrad from x-space to v-space.
290!
291 IF (obserr(iobs).ne.0.0_r8) THEN
292 pgrad(iobs)=pgrad(iobs)/sqrt(obserr(iobs))
293 END IF
294 vgrad0(iobs)=pgrad(iobs)
295 vgrad0s(iobs)=-vgrad0(iobs)
296 END DO
297!
298! If preconditioning, convert pgrad from v-space to y-space.
299!
300 IF (lprecond.and.(outloop.gt.1)) THEN
301 lscale=2 ! SQRT spectral LMP
302 ltrans=.true.
303 CALL rprecond (ng, lscale, ltrans, outloop, ninnloop, &
304 & pgrad)
305 END IF
306!
307 cg_gnorm_y(outloop)=0.0_dp
308 cg_gnorm_v(outloop)=0.0_dp
309 DO iobs=1,ndatum(ng)
310 zgrad0(iobs,outloop)=pgrad(iobs)
311 cg_gnorm_y(outloop)=cg_gnorm_y(outloop)+ &
312 & zgrad0(iobs,outloop)* &
313 & zgrad0(iobs,outloop)
314 cg_gnorm_v(outloop)=cg_gnorm_v(outloop)+ &
315 & vgrad0(iobs)*vgrad0(iobs)
316 END DO
317 cg_gnorm_y(outloop)=sqrt(cg_gnorm_y(outloop))
318 cg_gnorm_v(outloop)=sqrt(cg_gnorm_v(outloop))
319 DO iobs=1,ndatum(ng)
320 zcglwk(iobs,1,outloop)=pgrad(iobs)/cg_gnorm_y(outloop)
321 admodval(iobs)=zcglwk(iobs,1,outloop)
322 END DO
323!
324! If preconditioning, convert ADmodVal from y-space to v-space.
325!
326 IF (lprecond.and.(outloop.gt.1)) THEN
327 lscale=2 ! SQRT spectral LMP
328 ltrans=.false.
329 CALL rprecond (ng, lscale, ltrans, outloop, ninnloop, &
330 & admodval)
331 END IF
332!
333! Convert ADmodVal from v-space to x-space.
334!
335 DO iobs=1,ndatum(ng)
336 IF (obserr(iobs).ne.0.0_r8) THEN
337 admodval(iobs)=admodval(iobs)/sqrt(obserr(iobs))
338 END IF
339 END DO
340
341 cg_qg(1,outloop)=0.0_dp
342 DO iobs=1,ndatum(ng)
343 cg_qg(1,outloop)=cg_qg(1,outloop)+ &
344 & zcglwk(iobs,1,outloop)* &
345 & zgrad0(iobs,outloop)
346 END DO
347
348 preducv=1.0_dp
349 preducy=1.0_dp
350!
351! Compute the actual penalty function terms.
352!
353 jb=0.0_dp
354 jobs=0.0_dp
355 DO iobs=1,ndatum(ng)
356 jobs=jobs+vgrad0s(iobs)*vgrad0s(iobs)
357 END DO
358 jact=jb+jobs
359
360 ELSE
361
362 IF (lcgini) THEN
363!
364! When outer>1 we need to evaluate Ax0 so for inner=0 we use
365! cg_pxsave in v-space as the starting vector.
366!
367 DO iobs=1,ndatum(ng)
368 admodval(iobs)=cg_pxsave(iobs)
369# if defined RBL4DVAR || \
370 defined rbl4dvar_ana_sensitivity || \
371 defined tl_rbl4dvar
372 cg_innov(iobs)=-obsscale(iobs)* &
373 & (obsval(iobs)-bckmodval(iobs))
374# else
375 cg_innov(iobs)=-obsscale(iobs)* &
376 & (obsval(iobs)-tlmodval(iobs))
377# endif
378 END DO
379!
380! Convert ADmodVal from v-space to x-space and cg_innov (the
381! contribution to the starting gradient) from x-space to v-space.
382!
383 DO iobs=1,ndatum(ng)
384 IF (obserr(iobs).ne.0.0_r8) THEN
385 admodval(iobs)=admodval(iobs)/sqrt(obserr(iobs))
386 cg_innov(iobs)=cg_innov(iobs)/sqrt(obserr(iobs))
387 END IF
388 END DO
389
390 ELSE
391
392 DO iobs=1,ndatum(ng)
393!
394! Convert gradient contribution from x-space to v-space.
395!
396 pgrad(iobs)=obsscale(iobs)*tlmodval(iobs)
397 gdgw(iobs)=pgrad(iobs)
398 IF (obserr(iobs).ne.0.0_r8) THEN
399 pgrad(iobs)=pgrad(iobs)/sqrt(obserr(iobs))
400 vgrad0s(iobs)=pgrad(iobs)+cg_innov(iobs)+ &
401 & cg_pxsave(iobs)
402 END IF
403!
404! Add I*x0=cg_pxsave contribution to the gradient and the term
405! -b=cg_innov (everything is in v-space at this point).
406!
407 pgrad(iobs)=pgrad(iobs)+obsscale(iobs)* &
408 & (cg_pxsave(iobs)+cg_innov(iobs))
409 vgrad0(iobs)=pgrad(iobs)
410 END DO
411!
412! If preconditioning, convert pgrad from v-space to y-space.
413!
414 IF (lprecond.and.(outloop.gt.1)) THEN
415 lscale=2 ! SQRT spectral LMP
416 ltrans=.true.
417 CALL rprecond (ng, lscale, ltrans, outloop, ninnloop, &
418 & pgrad)
419 END IF
420
421 cg_gnorm_y(outloop)=0.0_dp
422 cg_gnorm_v(outloop)=0.0_dp
423 DO iobs=1,ndatum(ng)
424 zgrad0(iobs,outloop)=pgrad(iobs)
425 cg_gnorm_y(outloop)=cg_gnorm_y(outloop)+ &
426 & zgrad0(iobs,outloop)* &
427 & zgrad0(iobs,outloop)
428 cg_gnorm_v(outloop)=cg_gnorm_v(outloop)+ &
429 & vgrad0(iobs)*vgrad0(iobs)
430 END DO
431 cg_gnorm_y(outloop)=sqrt(cg_gnorm_y(outloop))
432 cg_gnorm_v(outloop)=sqrt(cg_gnorm_v(outloop))
433!
434 DO iobs=1,ndatum(ng)
435 zcglwk(iobs,1,outloop)=pgrad(iobs)/cg_gnorm_y(outloop)
436 admodval(iobs)=zcglwk(iobs,1,outloop)
437 END DO
438!
439! If preconditioning, convert ADmodVal from y-space to v-space.
440!
441 IF (lprecond.and.(outloop.gt.1)) THEN
442 lscale=2 ! SQRT spectral LMP
443 ltrans=.false.
444 CALL rprecond (ng, lscale, ltrans, outloop, ninnloop, &
445 & admodval)
446 END IF
447!
448 DO iobs=1,ndatum(ng)
449 IF (obserr(iobs).ne.0.0_r8) THEN
450 admodval(iobs)=admodval(iobs)/sqrt(obserr(iobs))
451 END IF
452 END DO
453 cg_qg(1,outloop)=0.0_dp
454 DO iobs=1,ndatum(ng)
455 cg_qg(1,outloop)=cg_qg(1,outloop)+ &
456 & zcglwk(iobs,1,outloop)* &
457 & zgrad0(iobs,outloop)
458 END DO
459
460 END IF
461
462 END IF
463
464 preducv=1.0_dp
465 preducy=1.0_dp
466!
467!-----------------------------------------------------------------------
468! Minimization step, innLoop > 0.
469!-----------------------------------------------------------------------
470!
471 ELSE
472!
473! After the initialization, for every other inner loop, calculate a
474! new Lanczos vector, store it in the matrix, and update search.
475!
476 DO iobs=1,ndatum(ng)
477 pgrad(iobs)=obsscale(iobs)*tlmodval(iobs)
478# if defined RBL4DVAR || defined R4DVAR || \
479 defined sensitivity_4dvar || defined tl_rbl4dvar || \
480 defined tl_r4dvar
481 tlmodval_s(iobs,innloop,outloop)=tlmodval(iobs)
482# endif
483!
484! Convert gradient contribution from x-space to v-space.
485!
486 IF (obserr(iobs).ne.0.0_r8) THEN
487 pgrad(iobs)=pgrad(iobs)/sqrt(obserr(iobs))
488 END IF
489 END DO
490
491 DO iobs=1,ndatum(ng)
492 zt(iobs)=zcglwk(iobs,innloop,outloop)
493 END DO
494!
495! If preconditioning, convert zt from y-space to v-space.
496!
497 IF (lprecond.and.(outloop.gt.1)) THEN
498 lscale=2 ! SQRT spectral LMP
499 ltrans=.false.
500 CALL rprecond (ng, lscale, ltrans, outloop, ninnloop, zt)
501 END IF
502!
503 DO iobs=1,ndatum(ng)
504 pgrad(iobs)=pgrad(iobs)+obsscale(iobs)*zt(iobs)
505 END DO
506!
507! If preconditioning, convert pgrad from v-space to y-space.
508!
509 IF (lprecond.and.(outloop.gt.1)) THEN
510 lscale=2 ! SQRT spectral LMP
511 ltrans=.true.
512 CALL rprecond (ng, lscale, ltrans, outloop, ninnloop, pgrad)
513 END IF
514!
515 cg_delta(innloop,outloop)=0.0_dp
516 DO iobs=1,ndatum(ng)
517 cg_delta(innloop,outloop)=cg_delta(innloop,outloop)+ &
518 & zcglwk(iobs,innloop,outloop)* &
519 & pgrad(iobs)
520 END DO
521!
522! Exit, if not a positive definite algorithm.
523!
524 IF (cg_delta(innloop,outloop).le.0.0_dp) THEN
525 WRITE (stdout,20) cg_delta(innloop,outloop), outloop, &
526 & innloop
527 exit_flag=8
528 GO TO 10
529 END IF
530!
531! Compute the new Lanczos vector.
532!
533 DO iobs=1,ndatum(ng)
534 pgrad(iobs)=pgrad(iobs)-cg_delta(innloop,outloop)* &
535 & zcglwk(iobs,innloop,outloop)
536 END DO
537 IF (innloop.gt.1) THEN
538 DO iobs=1,ndatum(ng)
539 pgrad(iobs)=pgrad(iobs)-cg_beta(innloop,outloop)* &
540 & zcglwk(iobs,innloop-1,outloop)
541 END DO
542 END IF
543!
544! Orthonormalize against previous Lanczos vectors.
545!
546 DO ivec=innloop,1,-1
547 cg_dla(ivec,outloop)=0.0_dp
548 DO iobs=1,ndatum(ng)
549 cg_dla(ivec,outloop)=cg_dla(ivec,outloop)+pgrad(iobs)* &
550 & zcglwk(iobs,ivec,outloop)
551 END DO
552 DO iobs=1,ndatum(ng)
553 pgrad(iobs)=pgrad(iobs)- &
554 & cg_dla(ivec,outloop)* &
555 & zcglwk(iobs,ivec,outloop)
556 END DO
557 END DO
558!
559 cg_beta(innloop+1,outloop)=0.0_dp
560 DO iobs=1,ndatum(ng)
561 cg_beta(innloop+1,outloop)=cg_beta(innloop+1,outloop)+ &
562 & pgrad(iobs)*pgrad(iobs)
563 END DO
564 cg_beta(innloop+1,outloop)=sqrt(cg_beta(innloop+1,outloop))
565!
566 DO iobs=1,ndatum(ng)
567 zcglwk(iobs,innloop+1,outloop)=pgrad(iobs)/ &
568 & cg_beta(innloop+1,outloop)
569 END DO
570!
571 cg_qg(innloop+1,outloop)=0.0_dp
572 DO iobs=1,ndatum(ng)
573 cg_qg(innloop+1,outloop)=cg_qg(innloop+1,outloop)+ &
574 & zcglwk(iobs,innloop+1,outloop)* &
575 & zgrad0(iobs,outloop)
576 END DO
577# ifdef MINRES
578!
579! Use the minimum residual method as recommended by El Akkraoui and
580! Gauthier (2010, QJRMS, 136, 107-115) and described by Paige and
581! Saunders ("Sparse Indefinite Systems of Linear Equations", 1975,
582! SIAM Journal on Numerical Analysis, 617-619). Specifically we refer
583! to equations 6.10 and 6.11 of this paper.
584!
585! Perform a LQ factorization of the tridiagonal matrix.
586!
587 ztrit=0.0_dp
588 DO i=1,innloop
589 ztrit(i,i)=cg_delta(i,outloop)
590 END DO
591 DO i=1,innloop-1
592 ztrit(i,i+1)=cg_beta(i+1,outloop)
593 END DO
594 DO i=2,innloop
595 ztrit(i,i-1)=cg_beta(i,outloop)
596 END DO
597 CALL sqlq(innloop, ztrit, tau, zwork1)
598!
599! Isolate L=zLT and its transpose.
600!
601 zlt=0.0_dp
602 zltt=0.0_dp
603 DO i=1,innloop
604 DO j=1,i
605 zlt(i,j)=ztrit(i,j)
606 END DO
607 END DO
608 DO j=1,innloop
609 DO i=1,innloop
610 zltt(i,j)=zlt(j,i)
611 END DO
612 END DO
613!
614! Now form ze=beta1*Q*e1.
615!
616 ze=0.0_dp
617 DO i=1,innloop
618 ze(i)=-cg_qg(i,outloop)
619 END DO
620 DO i=1,innloop
621 DO j=1,innloop
622 zeref(j)=0.0_dp
623 END DO
624 zeref(i)=1.0_dp
625 DO j=i+1,innloop
626 zeref(j)=ztrit(i,j)
627 END DO
628 zsum=0.0_dp
629 DO j=1,innloop
630 zsum=zsum+ze(j)*zeref(j)
631 END DO
632 DO j=1,innloop
633 ze(j)=ze(j)-tau(i)*zsum*zeref(j)
634 END DO
635 END DO
636!
637! Now form ze=D*ze (using equation 5.6 and 6.5 also).
638!
639 zgk=sqrt(zlt(innloop,innloop)*zlt(innloop,innloop)+ &
640 & cg_beta(innloop+1,outloop)*cg_beta(innloop+1,outloop))
641 zck=zlt(innloop,innloop)/zgk
642 ze(innloop)=zck*ze(innloop)
643!
644! Now compute ze=inv(L')*ze.
645!
646 DO j=innloop,1,-1
647 ze(j)=ze(j)/zltt(j,j)
648 DO i=1,j-1
649 ze(i)=ze(i)-ze(j)*zltt(i,j)
650 END DO
651 END DO
652!
653! Copy the solution ze into zwork.
654!
655 DO i=1,innloop
656 zwork(i,3)=ze(i)
657 END DO
658# else
659!
660! Calculate the new solution based upon the new, orthonormalized
661! Lanczos vector. First, the tridiagonal system is solved by
662! decomposition and forward/back substitution.
663!
664 zbet=cg_delta(1,outloop)
665 zu(1)=-cg_qg(1,outloop)/zbet
666 DO ivec=2,innloop
667 zgam(ivec)=cg_beta(ivec,outloop)/zbet
668 zbet=cg_delta(ivec,outloop)-cg_beta(ivec,outloop)*zgam(ivec)
669 zu(ivec)=(-cg_qg(ivec,outloop)-cg_beta(ivec,outloop)* &
670 & zu(ivec-1))/zbet
671 END DO
672 zwork(innloop,3)=zu(innloop)
673
674 DO ivec=innloop-1,1,-1
675 zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
676 zwork(ivec,3)=zu(ivec)
677 END DO
678# endif
679 DO iobs=1,ndatum(ng)
680 zw(iobs)=zgrad0(iobs,outloop)+ &
681 & cg_beta(innloop+1,outloop)* &
682 & zcglwk(iobs,innloop+1,outloop)* &
683 & zwork(innloop,3)
684 END DO
685
686 DO iobs=1,ndatum(ng)
687 px(iobs)=0.0_dp
688 DO ivec=1,innloop
689 px(iobs)=px(iobs)+zcglwk(iobs,ivec,outloop)*zwork(ivec,3)
690 zw(iobs)=zw(iobs)-zcglwk(iobs,ivec,outloop)* &
691 & cg_qg(ivec,outloop)
692 END DO
693 END DO
694!
695! If preconditioning, convert px from y-space to v-space.
696! We will always keep px in v-space.
697!
698 IF (lprecond.and.(outloop.gt.1)) THEN
699 lscale=2 ! SQRT spectral LMP
700 ltrans=.false.
701 CALL rprecond (ng, lscale, ltrans, outloop, ninnloop, px)
702 END IF
703!
704! Compute the reduction in the gradient Ax-b (in y-space).
705!
706 preducy=0.0_dp
707 DO iobs=1,ndatum(ng)
708 preducy=preducy+zw(iobs)*zw(iobs)
709 END DO
710 preducy=sqrt(preducy)/cg_gnorm_y(outloop)
711!
712! Compute the reduction in the gradient Ax-b (in v-space).
713!
714 IF (lprecond.and.(outloop.gt.1)) THEN
715 lscale=-2 ! SQRT spectral LMP
716 ltrans=.true.
717 CALL rprecond (ng, lscale, ltrans, outloop, ninnloop, zw)
718 END IF
719!
720 preducv=0.0_dp
721 DO iobs=1,ndatum(ng)
722 preducv=preducv+zw(iobs)*zw(iobs)
723 END DO
724 preducv=sqrt(preducv)/cg_gnorm_v(outloop)
725!
726! Estimate the residual (in y-space) by:
727!
728! ||Ax-b|| = cg_beta(k+1,outLoop)*zwork(innLoop,3)
729!
730 eps=abs(cg_beta(innloop+1,outloop)*zwork(innloop,3))
731!
732! Estimate the various penalty function components based on
733! Bennett (2002), page 22 and pages 42-44:
734!
735! Jopt = value of the penalty function when true px
736! has been identified.
737! Jf = initial data misfit for the first-guess
738! Jdata = data misfit for the state estimate
739! Jmod = the model penalty function (the sum of the
740! dynamical, initial and boundary penalties)
741! Jb = the ACTUAL model penalty function.
742! Jobs = the ACTUAL data penalty function.
743! Jact = the ACTUAL total penalty function.
744!
745! NOTE: Unlike I4D-Var where we compute the actual cost function
746! for each iteration that decreases with increasing
747! number of iterations, the various penalty terms
748! represent the optimal values (i.e. the values when
749! optimal state estimate has been identified). Therefore
750! the various penalty terms will only be correct when the
751! optimal state has been found. Therefore, as the
752! R4D-Var proceeds, Jopt, Jdata and Jmod will asymptote
753! to their final values and WILL NOT necessary decrease
754! with increasing iteration number.
755!
756 IF (outloop.eq.1.or.(.not.lhotstart)) THEN
757 DO iobs=1,ndatum(ng)
758 jopt=jopt-px(iobs)*vgrad0(iobs)
759 IF (obserr(iobs).ne.0.0_r8) THEN
760 jf=jf+vgrad0(iobs)*vgrad0(iobs)
761 END IF
762 jdata=jdata+px(iobs)*px(iobs)
763 END DO
764 jmod=jopt-jdata
765 ELSE
766 DO iobs=1,ndatum(ng)
767 jopt=jopt-(px(iobs)+cg_pxsave(iobs))*cg_innov(iobs)
768 IF (obserr(iobs).ne.0.0_r8) THEN
769 jf=jf+cg_innov(iobs)*cg_innov(iobs)
770 END IF
771 jdata=jdata+ &
772 & (px(iobs)+cg_pxsave(iobs))* &
773 & (px(iobs)+cg_pxsave(iobs))
774 END DO
775 jmod=jopt-jdata
776 END IF
777!
778! Compute the actual penalty function terms.
779!
780 DO iobs=1,ndatum(ng)
781 zd(iobs)=vgrad0s(iobs)*sqrt(obserr(iobs))
782 END DO
783 DO iobs=1,ndatum(ng)
784 IF (obserr(iobs).ne.0.0_r8) THEN
785 zhv(iobs)=zd(iobs)/sqrt(obserr(iobs))
786 END IF
787 END DO
788 DO ivec=1,innloop
789 ztemp1(ivec)=0.0_dp
790 DO iobs=1,ndatum(ng)
791 ztemp1(ivec)=ztemp1(ivec)+ &
792 & zcglwk(iobs,ivec,outloop)*zhv(iobs)
793 END DO
794 END DO
795# ifdef MINRES
796!
797! Now form ze=Q*ztemp1.
798!
799 ze=0.0_dp
800 DO i=1,innloop
801 ze(i)=ztemp1(i)
802 END DO
803 DO i=1,innloop
804 DO j=1,innloop
805 zeref(j)=0.0_dp
806 END DO
807 zeref(i)=1.0_dp
808 DO j=i+1,innloop
809 zeref(j)=ztrit(i,j)
810 END DO
811 zsum=0.0_dp
812 DO j=1,innloop
813 zsum=zsum+ze(j)*zeref(j)
814 END DO
815 DO j=1,innloop
816 ze(j)=ze(j)-tau(i)*zsum*zeref(j)
817 END DO
818 END DO
819!
820! Now form ze=D*ze (using equation 5.6 and 6.5 also).
821!
822 zgk=sqrt(zlt(innloop,innloop)*zlt(innloop,innloop)+ &
823 & cg_beta(innloop+1,outloop)*cg_beta(innloop+1,outloop))
824 zck=zlt(innloop,innloop)/zgk
825 ze(innloop)=zck*ze(innloop)
826!
827! Now compute ze=inv(L')*ze.
828!
829 DO j=innloop,1,-1
830 ze(j)=ze(j)/zltt(j,j)
831 DO i=1,j-1
832 ze(i)=ze(i)-ze(j)*zltt(i,j)
833 END DO
834 END DO
835!
836! Copy the solution ze into ztemp2.
837!
838 DO i=1,innloop
839 ztemp2(i)=ze(i)
840 END DO
841!
842# else
843 zbet=cg_delta(1,outloop)
844 zu1(1)=ztemp1(1)/zbet
845 DO ivec=2,innloop
846 zgam(ivec)=cg_beta(ivec,outloop)/zbet
847 zbet=cg_delta(ivec,outloop)- &
848 & cg_beta(ivec,outloop)*zgam(ivec)
849 zu1(ivec)=(ztemp1(ivec)- &
850 & cg_beta(ivec,outloop)*zu1(ivec-1))/zbet
851 END DO
852 ztemp2(innloop)=zu1(innloop)
853
854 DO ivec=innloop-1,1,-1
855 zu1(ivec)=zu1(ivec)-zgam(ivec+1)*zu1(ivec+1)
856 ztemp2(ivec)=zu1(ivec)
857 END DO
858# endif
859 DO iobs=1,ndatum(ng)
860 zhv(iobs)=0.0_dp
861 END DO
862 DO iobs=1,ndatum(ng)
863 DO ivec=1,innloop
864 zhv(iobs)=zhv(iobs)+ &
865 & ztemp2(ivec)*tlmodval_s(iobs,ivec,outloop)
866 END DO
867 END DO
868!
869! Compute Jb.
870!
871 DO ivec=1,innloop
872 ztemp1(ivec)=0.0_dp
873 DO iobs=1,ndatum(ng)
874 IF (obserr(iobs).ne.0.0_r8) THEN
875 ztemp1(ivec)=ztemp1(ivec)+ &
876 & zcglwk(iobs,ivec,outloop)*zhv(iobs)/ &
877 & sqrt(obserr(iobs))
878 END IF
879 END DO
880 END DO
881# ifdef MINRES
882!
883! Now form ze=Q*ztemp1.
884!
885 ze=0.0_dp
886 DO i=1,innloop
887 ze(i)=ztemp1(i)
888 END DO
889 DO i=1,innloop
890 DO j=1,innloop
891 zeref(j)=0.0_dp
892 END DO
893 zeref(i)=1.0_dp
894 DO j=i+1,innloop
895 zeref(j)=ztrit(i,j)
896 END DO
897 zsum=0.0_dp
898 DO j=1,innloop
899 zsum=zsum+ze(j)*zeref(j)
900 END DO
901 DO j=1,innloop
902 ze(j)=ze(j)-tau(i)*zsum*zeref(j)
903 END DO
904 END DO
905!
906! Now form ze=D*ze (using equation 5.6 and 6.5 also).
907!
908 zgk=sqrt(zlt(innloop,innloop)*zlt(innloop,innloop)+ &
909 & cg_beta(innloop+1,outloop)*cg_beta(innloop+1,outloop))
910 zck=zlt(innloop,innloop)/zgk
911 ze(innloop)=zck*ze(innloop)
912!
913! Now compute ze=inv(L')*ze.
914!
915 DO j=innloop,1,-1
916 ze(j)=ze(j)/zltt(j,j)
917 DO i=1,j-1
918 ze(i)=ze(i)-ze(j)*zltt(i,j)
919 END DO
920 END DO
921!
922! Copy the solution ze into ztemp2.
923!
924 DO i=1,innloop
925 ztemp2(i)=ze(i)
926 END DO
927!
928# else
929 zbet=cg_delta(1,outloop)
930 zu2(1)=ztemp1(1)/zbet
931 DO ivec=2,innloop
932 zgam(ivec)=cg_beta(ivec,outloop)/zbet
933 zbet=cg_delta(ivec,outloop)- &
934 & cg_beta(ivec,outloop)*zgam(ivec)
935 zu2(ivec)=(ztemp1(ivec)- &
936 & cg_beta(ivec,outloop)*zu2(ivec-1))/zbet
937 END DO
938 ztemp2(innloop)=zu2(innloop)
939 DO ivec=innloop-1,1,-1
940 zu2(ivec)=zu2(ivec)-zgam(ivec+1)*zu2(ivec+1)
941 ztemp2(ivec)=zu2(ivec)
942 END DO
943# endif
944 DO iobs=1,ndatum(ng)
945 zht(iobs)=0.0_dp
946 END DO
947 DO iobs=1,ndatum(ng)
948 DO ivec=1,innloop
949 IF (obserr(iobs).ne.0.0_r8) THEN
950 zht(iobs)=zht(iobs)+ &
951 & ztemp2(ivec)*zcglwk(iobs,ivec,outloop)/ &
952 & sqrt(obserr(iobs))
953 END IF
954 END DO
955 END DO
956 jb=0.0_dp
957 DO iobs=1,ndatum(ng)
958 IF (obserr(iobs).ne.0.0_r8) THEN
959 cff=cg_pxsave(iobs)/sqrt(obserr(iobs))
960 jb=jb+zd(iobs)*zht(iobs)-2.0_dp*cff*zhv(iobs)+ &
961 & cff*gdgw(iobs)
962 END IF
963 END DO
964!
965! Compute Jobs.
966!
967 jobs=0.0_dp
968 IF (outloop.eq.1.or.(.not.lhotstart)) THEN
969 DO iobs=1,ndatum(ng)
970 IF (obserr(iobs).ne.0.0_r8) THEN
971 cff=zhv(iobs)-zd(iobs)
972 jobs=jobs+cff*cff/obserr(iobs)
973 END IF
974 END DO
975 ELSE
976 DO iobs=1,ndatum(ng)
977 IF (obserr(iobs).ne.0.0_r8) THEN
978 cff=gdgw(iobs)-zhv(iobs)+cg_innov(iobs)*sqrt(obserr(iobs))
979 jobs=jobs+cff*cff/obserr(iobs)
980 END IF
981 END DO
982 END IF
983 jact=jb+jobs
984!
985! Put the new trial solution into the adjoint vector for the next
986! loop or put the final solution into the adjoint vector on the
987! final inner-loop.
988!
989 IF (innloop.eq.ninnloop) THEN
990!
991! Note: px is already in v-space so there is no need to convert
992! if preconditioning; cg_pxsave is also in v-space.
993!
994 DO iobs=1,ndatum(ng)
995 admodval(iobs)=px(iobs)
996 END DO
997 IF (lhotstart) THEN
998 DO iobs=1,ndatum(ng)
999 admodval(iobs)=admodval(iobs)+cg_pxsave(iobs)
1000 END DO
1001 DO iobs=1,ndatum(ng)
1002 cg_pxsave(iobs)=admodval(iobs)
1003 END DO
1004 END IF
1005 ELSE
1006 DO iobs=1,ndatum(ng)
1007 admodval(iobs)=zcglwk(iobs,innloop+1,outloop)
1008 END DO
1009!
1010! If preconditioning, convert ADmodVal from y-space to v-space.
1011!
1012 IF (lprecond.and.(outloop.gt.1)) THEN
1013 lscale=2 ! SQRT spectral LMP
1014 ltrans=.false.
1015 CALL rprecond (ng, lscale, ltrans, outloop, ninnloop, &
1016 & admodval)
1017 END IF
1018 END IF
1019!
1020! Convert ADmodVal from v-space to x-space.
1021!
1022 DO iobs=1,ndatum(ng)
1023 IF (obserr(iobs).ne.0.0_r8) THEN
1024 admodval(iobs)=admodval(iobs)/sqrt(obserr(iobs))
1025 END IF
1026 END DO
1027
1028 END IF minimize
1029!
1030!-----------------------------------------------------------------------
1031! Report minimization progress.
1032!-----------------------------------------------------------------------
1033!
1034 WRITE (stdout,30) ndatum(ng), &
1035 & outloop, innloop, eps, &
1036 & outloop, innloop, preducy, &
1037 & outloop, innloop, preducv, &
1038 & outloop, innloop, jf, &
1039 & outloop, innloop, jdata, &
1040 & outloop, innloop, jmod, &
1041 & outloop, innloop, jopt, &
1042 & outloop, innloop, jb, &
1043 & outloop, innloop, jobs, &
1044 & outloop, innloop, jact, &
1045 & outloop, innloop
1046 DO ivec=1,innloop
1047 WRITE (stdout,40) ivec, cg_delta(ivec,outloop), &
1048 & cg_beta(ivec,outloop), zwork(ivec,3)
1049 END DO
1050 WRITE (stdout,50) innloop+1, cg_beta(innloop+1,outloop)
1051!
1052! If preconditioning, report the Ritz eigenvalues used.
1053!
1054 IF ((lhessianev.or.lprecond).and.(outloop.gt.1)) THEN
1055 IF (nritzev.gt.0) THEN
1056 WRITE (stdout,60) outloop, innloop, nconvritz
1057 ic=0
1058 DO ivec=1,ninnloop
1059 IF (ivec.gt.(ninnloop-nconvritz)) THEN
1060 string=' '
1061 ic=ic+1
1062 WRITE (stdout,70) ivec, cg_ritz(ivec,outloop-1), &
1063 & cg_ritzerr(ivec,outloop-1), &
1064 & trim(adjustl(string)), &
1065 & 'Used=', ic
1066 ELSE
1067 string=' '
1068 WRITE (stdout,80) ivec, cg_ritz(ivec,outloop-1), &
1069 & cg_ritzerr(ivec,outloop-1), &
1070 & trim(adjustl(string))
1071 END IF
1072 END DO
1073 ELSE
1074 WRITE (stdout,90) outloop, innloop, ritzmaxerr
1075 ic=0
1076 DO ivec=1,ninnloop
1077 IF (cg_ritzerr(ivec,outloop-1).le.ritzmaxerr) THEN
1078 string='converged'
1079 ic=ic+1
1080 WRITE (stdout,70) ivec, cg_ritz(ivec,outloop-1), &
1081 & cg_ritzerr(ivec,outloop-1), &
1082 & trim(adjustl(string)), &
1083 & 'Good=', ic
1084 ELSE
1085 string='not converged'
1086 WRITE (stdout,80) ivec, cg_ritz(ivec,outloop-1), &
1087 & cg_ritzerr(ivec,outloop-1), &
1088 & trim(adjustl(string))
1089 END IF
1090 END DO
1091 END IF
1092 END IF
1093!
1094!-----------------------------------------------------------------------
1095! If last inner loop, innLoop = NinnLoop.
1096!-----------------------------------------------------------------------
1097!
1098! Compute the eigenvalues and eigenvectors of the tridiagonal matrix.
1099!
1100 IF (innloop.eq.ninnloop) THEN
1101 IF (lhessianev.or.lprecond) THEN
1102 DO ivec=1,innloop
1103 cg_ritz(ivec,outloop)=cg_delta(ivec,outloop)
1104 END DO
1105 DO ivec=1,innloop-1
1106 zwork(ivec,1)=cg_beta(ivec+1,outloop)
1107 END DO
1108!
1109! Use the LAPACK routine DSTEQR to compute the eigenvectors and
1110! eigenvalues of the tridiagonal matrix. If applicable, the
1111! eigenpairs are computed by master thread only.
1112!
1113 CALL dsteqr ('I', innloop, cg_ritz(1,outloop), zwork(1,1), &
1114 & zgv, ninner, work, info)
1115 IF (info.ne.0) THEN
1116 WRITE (stdout,*) ' CONGRAD - Error in DSTEQR: info = ', &
1117 & info
1118 exit_flag=8
1119 GO TO 10
1120 END IF
1121!
1122 DO j=1,ninner
1123 DO i=1,ninner
1124 cg_zv(i,j,outloop)=zgv(i,j)
1125 END DO
1126 END DO
1127!
1128! Estimate the Ritz value error bounds.
1129!
1130 DO ivec=1,innloop
1131 cg_ritzerr(ivec,outloop)=abs(cg_beta(innloop+1,outloop)* &
1132 & cg_zv(innloop,ivec,outloop))
1133 END DO
1134!
1135! Check for exploding or negative Ritz values.
1136!
1137 DO ivec=1,innloop
1138 IF (cg_ritzerr(ivec,outloop).lt.0.0_dp) THEN
1139 WRITE (stdout,*) ' CONGRAD - negative Ritz value found.'
1140 exit_flag=8
1141 GO TO 10
1142 END IF
1143 END DO
1144!
1145! Change the error bounds for the acceptable eigenvectors.
1146!
1147 DO ivec=1,ninnloop
1148 cg_ritzerr(ivec,outloop)=cg_ritzerr(ivec,outloop)/ &
1149 & cg_ritz(ninnloop,outloop)
1150 END DO
1151!
1152! Report Eigenvalues and their relative accuracy.
1153!
1154 WRITE (stdout,100) outloop, innloop, ritzmaxerr
1155 ic=0
1156 DO ivec=1,ninnloop
1157 IF (cg_ritzerr(ivec,outloop).le.ritzmaxerr) THEN
1158 string='converged'
1159 ic=ic+1
1160 WRITE (stdout,70) ivec, cg_ritz(ivec,outloop), &
1161 & cg_ritzerr(ivec,outloop), &
1162 & trim(adjustl(string)), &
1163 & 'Good=', ic
1164 ELSE
1165 string='not converged'
1166 WRITE (stdout,80) ivec, cg_ritz(ivec,outloop), &
1167 & cg_ritzerr(ivec,outloop), &
1168 & trim(adjustl(string))
1169 END IF
1170 END DO
1171!
1172! Calculate the converged eigenvectors (vcglev) of the [R_n + Cobs].
1173!
1174 CALL rpevecs (ng, outloop, ninnloop)
1175 END IF
1176 END IF
1177
1178 END IF master_thread
1179!
1180!-----------------------------------------------------------------------
1181! Come here if not a possite definite algorithm.
1182!-----------------------------------------------------------------------
1183!
1184 10 CONTINUE
1185# ifdef DISTRIBUTE
1186 CALL mp_bcasti (ng, model, exit_flag)
1187# endif
1188 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1189
1190# ifdef DISTRIBUTE
1191!
1192!-----------------------------------------------------------------------
1193! Broadcast new solution to other nodes.
1194!-----------------------------------------------------------------------
1195!
1196 CALL mp_bcastf (ng, model, admodval)
1197# if defined RBL4DVAR || defined R4DVAR || \
1198 defined sensitivity_4dvar || defined tl_rbl4dvar || \
1199 defined tl_r4dvar
1200!! CALL mp_bcastf (ng, model, TLmodVal_S(:,:,outLoop)) ! not needed
1201# endif
1202 CALL mp_bcasti (ng, model, info)
1203 CALL mp_bcastf (ng, model, cg_beta(:,outloop))
1204 CALL mp_bcastf (ng, model, cg_qg(:,outloop))
1205 CALL mp_bcastf (ng, model, cg_delta(:,outloop))
1206 CALL mp_bcastf (ng, model, cg_dla(:,outloop))
1207 CALL mp_bcastf (ng, model, cg_gnorm_y(:))
1208 CALL mp_bcastf (ng, model, cg_gnorm_v(:))
1209 CALL mp_bcastf (ng, model, cg_pxsave(:))
1210 CALL mp_bcastf (ng, model, cg_innov(:))
1211 CALL mp_bcastf (ng, model, zgrad0(:,outloop))
1212 CALL mp_bcastf (ng, model, zcglwk(:,:,outloop))
1213 IF ((lhessianev.or.lprecond).and.(innloop.eq.ninnloop)) THEN
1214 CALL mp_bcastf (ng, model, cg_ritz(:,outloop))
1215 CALL mp_bcastf (ng, model, cg_ritzerr(:,outloop))
1216 CALL mp_bcastf (ng, model, cg_zv(:,:,outloop))
1217 CALL mp_bcastf (ng, model, vcglev(:,:,outloop))
1218 END IF
1219# endif
1220!
1221!-----------------------------------------------------------------------
1222! Write out conjugate gradient vectors into 4D-Var NetCDF file.
1223!-----------------------------------------------------------------------
1224!
1225 CALL cg_write_congrad (ng, model, innloop, outloop, ninnloop, &
1226 & jf, jdata, jmod, jopt, jb, jobs, jact, &
1227 & preducv, preducy)
1228
1229# ifdef PROFILE
1230!
1231! Turn off time clock.
1232!
1233 CALL wclock_off (ng, model, 85, __line__, myfile)
1234# endif
1235!
1236 20 FORMAT (/,' CONGRAD - Fatal error, not possitive definite', &
1237 & ' algorithm:',/, &
1238 & /,11x,'cg_delta = ',1p,e15.8,0p,3x,'(',i3.3,', ',i3.3,')')
1239 30 FORMAT (/,' CONGRAD - Conjugate Gradient Information:',/, &
1240 & /,11x,'Ndatum = ',i7.7,/,/, &
1241 & 1x,'(',i3.3,',',i3.3,'): ', &
1242 & 'Residual estimate, eps = ', &
1243 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1244 & 'Reduction in gradient norm, Greduc y-space = ', &
1245 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1246 & 'Reduction in gradient norm, Greduc v-space = ', &
1247 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1248 & 'First guess initial data misfit, Jf = ', &
1249 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1250 & 'State estimate data misfit, Jdata = ', &
1251 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1252 & 'Model penalty function, Jmod = ', &
1253 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1254 & 'Optimal penalty function, Jopt = ', &
1255 & 1p,e14.7,/,/,1x,'(',i3.3,',',i3.3,'): ', &
1256 & 'Actual Model penalty function, Jb = ', &
1257 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1258 & 'Actual data penalty function, Jobs = ', &
1259 & 1p,e14.7,/,/,1x,'(',i3.3,',',i3.3,'): ', &
1260 & 'Actual total penalty function, Jact = ', &
1261 & 1p,e14.7,/,/,1x,'(',i3.3,',',i3.3,'): ', &
1262 & 'Lanczos vectors - cg_delta, cg_beta, zwork:',/)
1263 40 FORMAT (6x,i3.3,4x,3(1p,e15.8,5x))
1264 50 FORMAT (6x,i3.3,24x,1p,e15.8)
1265 60 FORMAT (/,1x,'(',i3.3,',',i3.3,'): ', &
1266 & 'Ritz eigenvalues used in preconditioning, ', &
1267 & 'nConvRitz = ',i3.3,/)
1268 70 FORMAT (6x,i3.3,2x,1p,e14.7,2x,1p,e14.7,2x,a,5x,'('a,i3.3,')')
1269 80 FORMAT (6x,i3.3,2x,1p,e14.7,2x,1p,e14.7,2x,a)
1270 90 FORMAT (/,1x,'(',i3.3,',',i3.3,'): ', &
1271 & 'Ritz eigenvalues used in preconditioning, ', &
1272 & 'RitzMaxErr = ',1p,e12.5,/)
1273100 FORMAT (/,1x,'(',i3.3,',',i3.3,'): ', &
1274 & 'New Ritz eigenvalues and their accuracy, ', &
1275 & 'RitzMaxErr = ',1p,e12.5,/)
1276!
1277 RETURN
1278 END SUBROUTINE congrad
1279!
1280!***********************************************************************
1281 SUBROUTINE rpevecs (ng, outLoop, NinnLoop)
1282!***********************************************************************
1283! !
1284! This routine computes the converged eigenvectors of (R_n+Cobs). !
1285! !
1286!***********************************************************************
1287!
1288! Imported variable declarations
1289!
1290 integer, intent(in) :: ng, outLoop, NinnLoop
1291!
1292! Local variable declarations.
1293!
1294 integer :: iobs, ivec, jn, jk, ij, ingood
1295!
1296 real(dp) :: dla
1297!
1298!-----------------------------------------------------------------------
1299! Compute and store converged eigenvectors of (R_n+Cobs).
1300!-----------------------------------------------------------------------
1301!
1302! Count and collect the converged eigenvalues.
1303!
1304 ingood=0
1305 DO ivec=ninnloop,1,-1
1306 ingood=ingood+1
1307 ritz(ingood)=cg_ritz(ivec,outloop)
1308 END DO
1309!
1310! Calculate the converged eigenvectors of (R_n+Cobs).
1311!
1312 ij=0
1313 DO jn=ninnloop,1,-1
1314 ij=ij+1
1315 DO iobs=1,ndatum(ng)
1316 vcglev(iobs,ij,outloop)=0.0_dp
1317 END DO
1318 DO jk=1,ninnloop
1319 DO iobs=1,ndatum(ng)
1320 vcglev(iobs,ij,outloop)=vcglev(iobs,ij,outloop)+ &
1321 & zcglwk(iobs,jk,outloop)* &
1322 & cg_zv(jk,jn,outloop)
1323 END DO
1324 END DO
1325 DO jk=1,ij-1
1326 dla=0.0_dp
1327 DO iobs=1,ndatum(ng)
1328 dla=dla+vcglev(iobs,jk,outloop)*vcglev(iobs,ij,outloop)
1329 END DO
1330 DO iobs=1,ndatum(ng)
1331 vcglev(iobs,ij,outloop)=vcglev(iobs,ij,outloop)- &
1332 & dla*vcglev(iobs,jk,outloop)
1333 END DO
1334 END DO
1335 dla=0.0_dp
1336 DO iobs=1,ndatum(ng)
1337 dla=dla+vcglev(iobs,ij,outloop)*vcglev(iobs,ij,outloop)
1338 END DO
1339 DO iobs=1,ndatum(ng)
1340 vcglev(iobs,ij,outloop)=vcglev(iobs,ij,outloop)/sqrt(dla)
1341 END DO
1342 END DO
1343!
1344 RETURN
1345 END SUBROUTINE rpevecs
1346
1347# ifdef SINGLE_PRECISION
1348!
1349!***********************************************************************
1350 SUBROUTINE rprecond_dp (ng, Lscale, Ltrans, outLoop, NinnLoop, py)
1351!***********************************************************************
1352! !
1353! This routine is the preconditioner in observation space. !
1354! !
1355!***********************************************************************
1356!
1357! Imported variable declarations
1358!
1359 logical, intent(in) :: Ltrans
1360!
1361 integer, intent(in) :: ng, Lscale, outLoop, NinnLoop
1362!
1363 real(dp), intent(inout) :: py(Ndatum(ng))
1364!
1365! Local variable declarations.
1366!
1367 integer :: is, ie, inc, iss, i
1368 integer :: nol, nols, nole, ninc
1369 integer :: ingood
1370 integer :: iobs, nvec
1371!
1372 real(dp) :: dla, fac, facritz
1373!
1374! Set the do-loop indices for the sequential preconditioner
1375! loop.
1376!
1377 IF (ltrans) THEN
1378 nols=outloop-1
1379 nole=1
1380 ninc=-1
1381 ELSE
1382 nols=1
1383 nole=outloop-1
1384 ninc=1
1385 END IF
1386!
1387! Apply the preconditioners derived from all previous outer-loops
1388! sequentially.
1389!
1390 DO nol=nols,nole,ninc
1391!
1392! Determine the number of Ritz vectors to use.
1393! For Lritz=.TRUE., choose HvecErr to be larger.
1394!
1395 ingood=0
1396 DO i=1,ninnloop
1397 IF (cg_ritzerr(i,nol).le.ritzmaxerr) THEN
1398 ingood=ingood+1
1399 END IF
1400 END DO
1401 IF (nritzev.gt.0) THEN
1402 nconvritz=nritzev
1403 ingood=nconvritz
1404 ELSE
1405 nconvritz=ingood
1406 END IF
1407!
1408 IF (lscale.gt.0) THEN
1409 is=1
1410 ie=ingood
1411 inc=1
1412 ELSE
1413 is=ingood
1414 ie=1
1415 inc=-1
1416 END IF
1417!
1418 IF (ltrans) THEN
1419 iss=is
1420 is=ie
1421 ie=iss
1422 inc=-inc
1423 END IF
1424!
1425! Lscale determines the form of the preconditioner:
1426!
1427! 1= LMP (NOT CODED YET)
1428! -1= Inverse LMP (NOT CODED YET)
1429! 2= Square root LMP
1430! -2= Inverse square root LMP
1431!
1432 DO nvec=is,ie,inc
1433!
1434 fac=0.0_dp
1435!
1436 IF (lritz) THEN
1437!
1438! Note that the primitive Ritz vectors cg_zv are stored in order of
1439! ascending Ritz values.
1440!
1441 facritz=cg_beta(ninnloop+1,nol)* &
1442 & cg_zv(ninnloop,ninnloop+1-nvec,nol)
1443!
1444! Compute dot-product of py with q_k+1 from previous outer-loop.
1445!
1446 IF (.not.ltrans) THEN
1447 dla=0.0_dp
1448 DO iobs=1,ndatum(ng)
1449 dla=dla+zcglwk(iobs,ninnloop+1,nol)*py(iobs)
1450 END DO
1451 facritz=facritz*dla
1452 END IF
1453
1454 END IF
1455!
1456! Compute the dot-product of py with the Ritz vectors from previous
1457! outer-loop.
1458!
1459! Note that vcglev are arranged in order of descending Ritz values,
1460! while the Ritz values themselves that are stored in cg_Ritz are in
1461! ascending order.
1462!
1463 dla=0.0_dp
1464 DO iobs=1,ndatum(ng)
1465 dla=dla+vcglev(iobs,nvec,nol)*py(iobs)
1466 END DO
1467!
1468! NOTE: Lscale=1 and Lscale=-1 cases are not complete yet.
1469!
1470 IF (lscale.eq.-1) THEN
1471 fac=(cg_ritz(ninnloop+1-nvec,nol)-1.0_dp)*dla
1472 ELSE IF (lscale.eq.1) THEN
1473 fac=(1.0_dp/cg_ritz(ninnloop+1-nvec,nol)-1.0_dp)*dla
1474 ELSE IF (lscale.eq.-2) THEN
1475 fac=(sqrt(cg_ritz(ninnloop+1-nvec,nol))-1.0_dp)*dla
1476 ELSE IF (lscale.eq.2) THEN
1477 fac=(1.0_dp/sqrt(cg_ritz(ninnloop+1-nvec,nol))-1.0_dp)*dla
1478 END IF
1479!
1480 IF (.not.ltrans) THEN
1481 IF (lritz.and.(lscale.eq.-2)) THEN
1482 fac=fac+facritz/sqrt(cg_ritz(ninnloop+1-nvec,nol))
1483 END IF
1484 IF (lritz.and.(lscale.eq.2)) THEN
1485 fac=fac-facritz/cg_ritz(ninnloop+1-nvec,nol)
1486 END IF
1487 END IF
1488!
1489 DO iobs=1,ndatum(ng)
1490 py(iobs)=py(iobs)+fac*vcglev(iobs,nvec,nol)
1491 END DO
1492!
1493 IF (lritz.and.ltrans) THEN
1494
1495 IF (lscale.eq.2) THEN
1496 fac=-facritz*dla/cg_ritz(ninnloop+1-nvec,nol)
1497 END IF
1498 IF (lscale.eq.-2) THEN
1499 fac=facritz*dla/sqrt(cg_ritz(ninnloop+1-nvec,nol))
1500 END IF
1501
1502 DO iobs=1,ndatum(ng)
1503 py(iobs)=py(iobs)+fac*zcglwk(iobs,ninnloop+1,nol)
1504 END DO
1505
1506 END IF
1507
1508 END DO
1509 END DO
1510!
1511 RETURN
1512 END SUBROUTINE rprecond_dp
1513# endif
1514!
1515!***********************************************************************
1516 SUBROUTINE rprecond_r8 (ng, Lscale, Ltrans, outLoop, NinnLoop, py)
1517!***********************************************************************
1518! !
1519! This routine is the preconditioner in observation space. !
1520! !
1521!***********************************************************************
1522!
1523! Imported variable declarations
1524!
1525 logical, intent(in) :: Ltrans
1526!
1527 integer, intent(in) :: ng, Lscale, outLoop, NinnLoop
1528!
1529 real(r8), intent(inout) :: py(Ndatum(ng))
1530!
1531! Local variable declarations.
1532!
1533 integer :: is, ie, inc, iss, i
1534 integer :: nol, nols, nole, ninc
1535 integer :: ingood
1536 integer :: iobs, nvec
1537!
1538 real(r8) :: dla, fac, facritz
1539!
1540! Set the do-loop indices for the sequential preconditioner
1541! loop.
1542!
1543 IF (ltrans) THEN
1544 nols=outloop-1
1545 nole=1
1546 ninc=-1
1547 ELSE
1548 nols=1
1549 nole=outloop-1
1550 ninc=1
1551 END IF
1552!
1553! Apply the preconditioners derived from all previous outer-loops
1554! sequentially.
1555!
1556 DO nol=nols,nole,ninc
1557!
1558! Determine the number of Ritz vectors to use.
1559! For Lritz=.TRUE., choose HvecErr to be larger.
1560!
1561 ingood=0
1562 DO i=1,ninnloop
1563 IF (cg_ritzerr(i,nol).le.ritzmaxerr) THEN
1564 ingood=ingood+1
1565 END IF
1566 END DO
1567 IF (nritzev.gt.0) THEN
1568 nconvritz=nritzev
1569 ingood=nconvritz
1570 ELSE
1571 nconvritz=ingood
1572 END IF
1573!
1574 IF (lscale.gt.0) THEN
1575 is=1
1576 ie=ingood
1577 inc=1
1578 ELSE
1579 is=ingood
1580 ie=1
1581 inc=-1
1582 END IF
1583!
1584 IF (ltrans) THEN
1585 iss=is
1586 is=ie
1587 ie=iss
1588 inc=-inc
1589 END IF
1590!
1591! Lscale determines the form of the preconditioner:
1592!
1593! 1= LMP (NOT CODED YET)
1594! -1= Inverse LMP (NOT CODED YET)
1595! 2= Square root LMP
1596! -2= Inverse square root LMP
1597!
1598 DO nvec=is,ie,inc
1599!
1600 fac=0.0_dp
1601!
1602 IF (lritz) THEN
1603!
1604! Note that the primitive Ritz vectors cg_zv are stored in order of
1605! ascending Ritz values.
1606!
1607 facritz=cg_beta(ninnloop+1,nol)* &
1608 & cg_zv(ninnloop,ninnloop+1-nvec,nol)
1609!
1610! Compute dot-product of py with q_k+1 from previous outer-loop.
1611!
1612 IF (.not.ltrans) THEN
1613 dla=0.0_dp
1614 DO iobs=1,ndatum(ng)
1615 dla=dla+zcglwk(iobs,ninnloop+1,nol)*py(iobs)
1616 END DO
1617 facritz=facritz*dla
1618 END IF
1619
1620 END IF
1621!
1622! Compute the dot-product of py with the Ritz vectors from previous
1623! outer-loop.
1624!
1625! Note that vcglev are arranged in order of descending Ritz values,
1626! while the Ritz values themselves that are stored in cg_Ritz are in
1627! ascending order.
1628!
1629 dla=0.0_dp
1630 DO iobs=1,ndatum(ng)
1631 dla=dla+vcglev(iobs,nvec,nol)*py(iobs)
1632 END DO
1633!
1634! NOTE: Lscale=1 and Lscale=-1 cases are not complete yet.
1635!
1636 IF (lscale.eq.-1) THEN
1637 fac=(cg_ritz(ninnloop+1-nvec,nol)-1.0_dp)*dla
1638 ELSE IF (lscale.eq.1) THEN
1639 fac=(1.0_dp/cg_ritz(ninnloop+1-nvec,nol)-1.0_dp)*dla
1640 ELSE IF (lscale.eq.-2) THEN
1641 fac=(sqrt(cg_ritz(ninnloop+1-nvec,nol))-1.0_dp)*dla
1642 ELSE IF (lscale.eq.2) THEN
1643 fac=(1.0_dp/sqrt(cg_ritz(ninnloop+1-nvec,nol))-1.0_dp)*dla
1644 END IF
1645!
1646 IF (.not.ltrans) THEN
1647 IF (lritz.and.(lscale.eq.-2)) THEN
1648 fac=fac+facritz/sqrt(cg_ritz(ninnloop+1-nvec,nol))
1649 END IF
1650 IF (lritz.and.(lscale.eq.2)) THEN
1651 fac=fac-facritz/cg_ritz(ninnloop+1-nvec,nol)
1652 END IF
1653 END IF
1654!
1655 DO iobs=1,ndatum(ng)
1656 py(iobs)=py(iobs)+fac*vcglev(iobs,nvec,nol)
1657 END DO
1658!
1659 IF (lritz.and.ltrans) THEN
1660
1661 IF (lscale.eq.2) THEN
1662 fac=-facritz*dla/cg_ritz(ninnloop+1-nvec,nol)
1663 END IF
1664 IF (lscale.eq.-2) THEN
1665 fac=facritz*dla/sqrt(cg_ritz(ninnloop+1-nvec,nol))
1666 END IF
1667
1668 DO iobs=1,ndatum(ng)
1669 py(iobs)=py(iobs)+fac*zcglwk(iobs,ninnloop+1,nol)
1670 END DO
1671
1672 END IF
1673
1674 END DO
1675 END DO
1676!
1677 RETURN
1678 END SUBROUTINE rprecond_r8
1679!
1680!***********************************************************************
1681 SUBROUTINE cg_write_congrad (ng, model, innLoop, outLoop, &
1682 & NinnLoop, Jf, Jdata, Jmod, Jopt, &
1683 & Jb, Jobs, Jact, &
1684 & preducv, preducy)
1685!***********************************************************************
1686! !
1687! This routine writes conjugate gradient vectors into DAV file !
1688! using either the standard NetCDF library or the Parallel-IO (PIO) !
1689! library. It saves all the variables needed for split or restart !
1690! schemes. !
1691! !
1692!***********************************************************************
1693!
1694! Imported variable declarations
1695!
1696 integer, intent(in) :: ng, model, innloop, outloop, ninnloop
1697!
1698 real(dp), intent(in) :: jf, jdata, jmod, jopt, preducv, preducy
1699 real(dp), intent(in) :: jb, jobs, jact
1700!
1701! Local variable declarations.
1702!
1703 character (len=*), parameter :: myfile = &
1704 & __FILE__//", cg_write_congrad"
1705!
1706 sourcefile=myfile
1707!
1708!-----------------------------------------------------------------------
1709! Write out conjugate gradient variables.
1710!-----------------------------------------------------------------------
1711!
1712 SELECT CASE (dav(ng)%IOtype)
1713 CASE (io_nf90)
1714 CALL cg_write_congrad_nf90 (ng, model, innloop, outloop, &
1715 & ninnloop, jf, jdata, jmod, jopt, &
1716 & jb, jobs, jact, &
1717 & preducv, preducy)
1718
1719# if defined PIO_LIB && defined DISTRIBUTE
1720 CASE (io_pio)
1721 CALL cg_write_congrad_pio (ng, model, innloop, outloop, &
1722 & ninnloop, jf, jdata, jmod, jopt, &
1723 & jb, jobs, jact, &
1724 & preducv, preducy)
1725# endif
1726 CASE DEFAULT
1727 IF (master) WRITE (stdout,10) dav(ng)%IOtype
1728 exit_flag=3
1729 END SELECT
1730 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1731!
1732 10 FORMAT (' CG_WRITE_CONGRAD - Illegal output file type,', &
1733 & ' io_type = ',i0, &
1734 & /,20x,'Check KeyWord ''OUT_LIB'' in ''roms.in''.')
1735!
1736 RETURN
1737 END SUBROUTINE cg_write_congrad
1738!
1739!***********************************************************************
1740 SUBROUTINE cg_write_congrad_nf90 (ng, model, innLoop, outLoop, &
1741 & NinnLoop, Jf, Jdata, Jmod, &
1742 & Jopt, Jb, Jobs, Jact, &
1743 & preducv, preducy)
1744!***********************************************************************
1745!
1746 USE mod_netcdf
1747!
1748! Imported variable declarations
1749!
1750 integer, intent(in) :: ng, model, innloop, outloop, ninnloop
1751!
1752 real(dp), intent(in) :: jf, jdata, jmod, jopt, preducv, preducy
1753 real(dp), intent(in) :: jb, jobs, jact
1754!
1755! Local variable declarations.
1756!
1757 integer :: nconv, status
1758!
1759 character (len=*), parameter :: myfile = &
1760 & __FILE__//", cg_write_congrad_nf90"
1761!
1762 sourcefile=myfile
1763!
1764!-----------------------------------------------------------------------
1765! Write out conjugate gradient variables.
1766!-----------------------------------------------------------------------
1767!
1768! Write out outer and inner iteration.
1769!
1770 CALL netcdf_put_ivar (ng, model, dav(ng)%name, &
1771 & 'outer', outer, &
1772 & (/0/), (/0/), &
1773 & ncid = dav(ng)%ncid)
1774 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1775
1776 CALL netcdf_put_ivar (ng, model, dav(ng)%name, &
1777 & 'inner', inner, &
1778 & (/0/), (/0/), &
1779 & ncid = dav(ng)%ncid)
1780 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1781!
1782! Write out number of converged Ritz eigenvalues.
1783!
1784 IF (innloop.eq.ninner) THEN
1785 CALL netcdf_put_ivar (ng, model, dav(ng)%name, &
1786 & 'nConvRitz', nconvritz, &
1787 & (/outloop/), (/1/), &
1788 & ncid = dav(ng)%ncid)
1789 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1790 END IF
1791!
1792! Write out converged Ritz eigenvalues.
1793!
1794 IF (innloop.eq.ninner) THEN
1795 nconv=max(1,nconvritz)
1796 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1797 & 'Ritz', ritz(1:nconv), &
1798 & (/1,outloop/), (/nconv,1/), &
1799 & ncid = dav(ng)%ncid)
1800 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1801 END IF
1802!
1803! Write out conjugate gradient norm.
1804!
1805 IF (innloop.gt.0) THEN
1806 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1807 & 'cg_beta', cg_beta, &
1808 & (/1,1/), (/ninner+1,nouter/), &
1809 & ncid = dav(ng)%ncid)
1810 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1811 END IF
1812!
1813! Write out Lanczos algorithm coefficients.
1814!
1815 IF (innloop.gt.0) THEN
1816 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1817 & 'cg_delta', cg_delta, &
1818 & (/1,1/), (/ninner,nouter/), &
1819 & ncid = dav(ng)%ncid)
1820 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1821 END IF
1822!
1823! Write normalization coefficients for Lanczos vectos.
1824!
1825 IF (innloop.gt.0) THEN
1826 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1827 & 'cg_dla', cg_dla, &
1828 & (/1,1/), (/ninner,nouter/), &
1829 & ncid = dav(ng)%ncid)
1830 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1831 END IF
1832!
1833! Initial gradient normalization factor.
1834!
1835 IF (innloop.eq.0) THEN
1836 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1837 & 'cg_Gnorm_y', cg_gnorm_y, &
1838 & (/1/), (/nouter/), &
1839 & ncid = dav(ng)%ncid)
1840 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1841
1842 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1843 & 'cg_Gnorm_v', cg_gnorm_v, &
1844 & (/1/), (/nouter/), &
1845 & ncid = dav(ng)%ncid)
1846 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1847 END IF
1848!
1849! Lanczos vector normalization factor.
1850!
1851 IF (innloop.gt.0) THEN
1852 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1853 & 'cg_QG', cg_qg, &
1854 & (/1,1/), (/ninner+1,nouter/), &
1855 & ncid = dav(ng)%ncid)
1856 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1857 END IF
1858!
1859! Reduction in the gradient norm.
1860!
1861 IF (innloop.gt.0) THEN
1862 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1863 & 'cg_Greduc_y', preducy, &
1864 & (/innloop,outloop/), (/1,1/), &
1865 & ncid = dav(ng)%ncid)
1866 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1867
1868 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1869 & 'cg_Greduc_v', preducv, &
1870 & (/innloop,outloop/), (/1,1/), &
1871 & ncid = dav(ng)%ncid)
1872 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1873 END IF
1874!
1875! Eigenvalues of Lanczos recurrence relationship.
1876!
1877 IF (innloop.gt.0) THEN
1878 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1879 & 'cg_Ritz', cg_ritz, &
1880 & (/1,1/), (/ninner,nouter/), &
1881 & ncid = dav(ng)%ncid)
1882 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1883 END IF
1884!
1885! Eigenvalues relative error.
1886!
1887 IF (innloop.gt.0) THEN
1888 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1889 & 'cg_RitzErr', cg_ritzerr, &
1890 & (/1,1/), (/ninner,nouter/), &
1891 & ncid = dav(ng)%ncid)
1892 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1893 END IF
1894!
1895! Eigenvectors of Lanczos recurrence relationship.
1896!
1897 IF (innloop.gt.0) THEN
1898 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1899 & 'cg_zv', cg_zv(:,:,outloop), &
1900 & (/1,1,outloop/), (/ninner,ninner,1/), &
1901 & ncid = dav(ng)%ncid)
1902 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1903 END IF
1904!
1905! First guess initial data misfit.
1906!
1907 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1908 & 'Jf', jf, &
1909 & (/innloop+1,outloop/), (/1,1/), &
1910 & ncid = dav(ng)%ncid)
1911 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1912!
1913! State estimate data misfit.
1914!
1915 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1916 & 'Jdata', jdata, &
1917 & (/innloop+1,outloop/), (/1,1/), &
1918 & ncid = dav(ng)%ncid)
1919 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1920!
1921! Model penalty function.
1922!
1923 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1924 & 'Jmod', jmod, &
1925 & (/innloop+1,outloop/), (/1,1/), &
1926 & ncid = dav(ng)%ncid)
1927 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1928!
1929! Optimal penalty function.
1930!
1931 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1932 & 'Jopt', jopt, &
1933 & (/innloop+1,outloop/), (/1,1/), &
1934 & ncid = dav(ng)%ncid)
1935 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1936!
1937! Actual model penalty function.
1938!
1939 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1940 & 'Jb', jb, &
1941 & (/innloop+1,outloop/), (/1,1/), &
1942 & ncid = dav(ng)%ncid)
1943 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1944!
1945! Actual data penalty function.
1946!
1947 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1948 & 'Jobs', jobs, &
1949 & (/innloop+1,outloop/), (/1,1/), &
1950 & ncid = dav(ng)%ncid)
1951 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1952!
1953! Actual total penalty function.
1954!
1955 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1956 & 'Jact', jact, &
1957 & (/innloop+1,outloop/), (/1,1/), &
1958 & ncid = dav(ng)%ncid)
1959 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1960!
1961! First-guess estimate of representer coefficients.
1962!
1963 IF (innloop.eq.ninnloop) THEN
1964 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1965 & 'cg_pxsave', cg_pxsave, &
1966 & (/1,outloop/), (/ndatum(ng),1/), &
1967 & ncid = dav(ng)%ncid)
1968 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1969 END IF
1970!
1971! Initial gradient for minimization.
1972!
1973 IF (innloop.eq.0) THEN
1974 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1975 & 'zgrad0', zgrad0(:,outloop), &
1976 & (/1,outloop/), (/ndatum(ng),1/), &
1977 & ncid = dav(ng)%ncid)
1978 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1979 END IF
1980!
1981! Lanczos vectors.
1982!
1983 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1984 & 'zcglwk', zcglwk(:,innloop+1,outloop), &
1985 & (/1,innloop+1,outloop/), &
1986 & (/ndatum(ng),1,1/), &
1987 & ncid = dav(ng)%ncid)
1988 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1989!
1990! Saved TLmodVal_S.
1991!
1992 IF (innloop.gt.0) THEN
1993 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1994 & 'TLmodVal_S', &
1995 & tlmodval_s(:,innloop,outloop), &
1996 & (/1,innloop,outloop/), &
1997 & (/ndatum(ng),1,1/), &
1998 & ncid = dav(ng)%ncid)
1999 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2000 END IF
2001!
2002!-----------------------------------------------------------------------
2003! Synchronize model/observation NetCDF file to disk.
2004!-----------------------------------------------------------------------
2005!
2006 CALL netcdf_sync (ng, model, dav(ng)%name, dav(ng)%ncid)
2007!
2008 RETURN
2009 END SUBROUTINE cg_write_congrad_nf90
2010
2011# if defined PIO_LIB && defined DISTRIBUTE
2012!
2013!***********************************************************************
2014 SUBROUTINE cg_write_congrad_pio (ng, model, innLoop, outLoop, &
2015 & NinnLoop, Jf, Jdata, Jmod, &
2016 & Jopt, Jb, Jobs, Jact, &
2017 & preducv, preducy)
2018!***********************************************************************
2019!
2020 USE mod_pio_netcdf
2021!
2022! Imported variable declarations
2023!
2024 integer, intent(in) :: ng, model, innLoop, outLoop, NinnLoop
2025!
2026 real(dp), intent(in) :: Jf, Jdata, Jmod, Jopt, preducv, preducy
2027 real(dp), intent(in) :: Jb, Jobs, Jact
2028!
2029! Local variable declarations.
2030!
2031 integer :: Nconv, status
2032!
2033 character (len=*), parameter :: MyFile = &
2034 & __FILE__//", cg_write_congrad_pio"
2035!
2036 sourcefile=myfile
2037!
2038!-----------------------------------------------------------------------
2039! Write out conjugate gradient variables.
2040!-----------------------------------------------------------------------
2041!
2042! Write out outer and inner iteration.
2043!
2044 CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name, &
2045 & 'outer', outer, &
2046 & (/0/), (/0/), &
2047 & piofile = dav(ng)%pioFile)
2048 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2049
2050 CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name, &
2051 & 'inner', inner, &
2052 & (/0/), (/0/), &
2053 & piofile = dav(ng)%pioFile)
2054 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2055!
2056! Write out number of converged Ritz eigenvalues.
2057!
2058 IF (innloop.eq.ninner) THEN
2059 CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name, &
2060 & 'nConvRitz', nconvritz, &
2061 & (/outloop/), (/1/), &
2062 & piofile = dav(ng)%pioFile)
2063 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2064 END IF
2065!
2066! Write out converged Ritz eigenvalues.
2067!
2068 IF (innloop.eq.ninner) THEN
2069 nconv=max(1,nconvritz)
2070 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2071 & 'Ritz', ritz(1:nconv), &
2072 & (/1,outloop/), (/nconv,1/), &
2073 & piofile = dav(ng)%pioFile)
2074 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2075 END IF
2076!
2077! Write out conjugate gradient norm.
2078!
2079 IF (innloop.gt.0) THEN
2080 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2081 & 'cg_beta', cg_beta, &
2082 & (/1,1/), (/ninner+1,nouter/), &
2083 & piofile = dav(ng)%pioFile)
2084 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2085 END IF
2086!
2087! Write out Lanczos algorithm coefficients.
2088!
2089 IF (innloop.gt.0) THEN
2090 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2091 & 'cg_delta', cg_delta, &
2092 & (/1,1/), (/ninner,nouter/), &
2093 & piofile = dav(ng)%pioFile)
2094 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2095 END IF
2096!
2097! Write normalization coefficients for Lanczos vectos.
2098!
2099 IF (innloop.gt.0) THEN
2100 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2101 & 'cg_dla', cg_dla, &
2102 & (/1,1/), (/ninner,nouter/), &
2103 & piofile = dav(ng)%pioFile)
2104 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2105 END IF
2106!
2107! Initial gradient normalization factor.
2108!
2109 IF (innloop.eq.0) THEN
2110 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2111 & 'cg_Gnorm_y', cg_gnorm_y, &
2112 & (/1/), (/nouter/), &
2113 & piofile = dav(ng)%pioFile)
2114 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2115
2116 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2117 & 'cg_Gnorm_v', cg_gnorm_v, &
2118 & (/1/), (/nouter/), &
2119 & piofile = dav(ng)%pioFile)
2120 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2121 END IF
2122!
2123! Lanczos vector normalization factor.
2124!
2125 IF (innloop.gt.0) THEN
2126 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2127 & 'cg_QG', cg_qg, &
2128 & (/1,1/), (/ninner+1,nouter/), &
2129 & piofile = dav(ng)%pioFile)
2130 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2131 END IF
2132!
2133! Reduction in the gradient norm.
2134!
2135 IF (innloop.gt.0) THEN
2136 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2137 & 'cg_Greduc_y', preducy, &
2138 & (/innloop,outloop/), (/1,1/), &
2139 & piofile = dav(ng)%pioFile)
2140 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2141
2142 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2143 & 'cg_Greduc_v', preducv, &
2144 & (/innloop,outloop/), (/1,1/), &
2145 & piofile = dav(ng)%pioFile)
2146 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2147 END IF
2148!
2149! Eigenvalues of Lanczos recurrence relationship.
2150!
2151 IF (innloop.gt.0) THEN
2152 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2153 & 'cg_Ritz', cg_ritz, &
2154 & (/1,1/), (/ninner,nouter/), &
2155 & piofile = dav(ng)%pioFile)
2156 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2157 END IF
2158!
2159! Eigenvalues relative error.
2160!
2161 IF (innloop.gt.0) THEN
2162 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2163 & 'cg_RitzErr', cg_ritzerr, &
2164 & (/1,1/), (/ninner,nouter/), &
2165 & piofile = dav(ng)%pioFile)
2166 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2167 END IF
2168!
2169! Eigenvectors of Lanczos recurrence relationship.
2170!
2171 IF (innloop.gt.0) THEN
2172 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2173 & 'cg_zv', cg_zv(:,:,outloop), &
2174 & (/1,1,outloop/), (/ninner,ninner,1/), &
2175 & piofile = dav(ng)%pioFile)
2176 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2177 END IF
2178!
2179! First guess initial data misfit.
2180!
2181 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2182 & 'Jf', jf, &
2183 & (/innloop+1,outloop/), (/1,1/), &
2184 & piofile = dav(ng)%pioFile)
2185 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2186!
2187! State estimate data misfit.
2188!
2189 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2190 & 'Jdata', jdata, &
2191 & (/innloop+1,outloop/), (/1,1/), &
2192 & piofile = dav(ng)%pioFile)
2193 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2194!
2195! Model penalty function.
2196!
2197 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2198 & 'Jmod', jmod, &
2199 & (/innloop+1,outloop/), (/1,1/), &
2200 & piofile = dav(ng)%pioFile)
2201 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2202!
2203! Optimal penalty function.
2204!
2205 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2206 & 'Jopt', jopt, &
2207 & (/innloop+1,outloop/), (/1,1/), &
2208 & piofile = dav(ng)%pioFile)
2209 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2210!
2211! Actual model penalty function.
2212!
2213 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2214 & 'Jb', jb, &
2215 & (/innloop+1,outloop/), (/1,1/), &
2216 & piofile = dav(ng)%pioFile)
2217 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2218!
2219! Actual data penalty function.
2220!
2221 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2222 & 'Jobs', jobs, &
2223 & (/innloop+1,outloop/), (/1,1/), &
2224 & piofile = dav(ng)%pioFile)
2225 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2226!
2227! Actual total penalty function.
2228!
2229 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2230 & 'Jact', jact, &
2231 & (/innloop+1,outloop/), (/1,1/), &
2232 & piofile = dav(ng)%pioFile)
2233 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2234!
2235! First-guess estimate of representer coefficients.
2236!
2237 IF (innloop.eq.ninnloop) THEN
2238 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2239 & 'cg_pxsave', cg_pxsave, &
2240 & (/1,outloop/), (/ndatum(ng),1/), &
2241 & piofile = dav(ng)%pioFile)
2242 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2243 END IF
2244!
2245! Initial gradient for minimization.
2246!
2247 IF (innloop.eq.0) THEN
2248 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2249 & 'zgrad0', zgrad0(:,outloop), &
2250 & (/1,outloop/), (/ndatum(ng),1/), &
2251 & piofile = dav(ng)%pioFile)
2252 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2253 END IF
2254!
2255! Lanczos vectors.
2256!
2257 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2258 & 'zcglwk', zcglwk(:,innloop+1,outloop), &
2259 & (/1,innloop+1,outloop/), &
2260 & (/ndatum(ng),1,1/), &
2261 & piofile = dav(ng)%pioFile)
2262 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2263!
2264! Saved TLmodVal_S.
2265!
2266 IF (innloop.gt.0) THEN
2267 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2268 & 'TLmodVal_S', &
2269 & tlmodval_s(:,innloop,outloop), &
2270 & (/1,innloop,outloop/), &
2271 & (/ndatum(ng),1,1/), &
2272 & piofile = dav(ng)%pioFile)
2273 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2274 END IF
2275!
2276!-----------------------------------------------------------------------
2277! Synchronize model/observation NetCDF file to disk.
2278!-----------------------------------------------------------------------
2279!
2280 CALL pio_netcdf_sync (ng, model, dav(ng)%name, dav(ng)%pioFile)
2281!
2282 RETURN
2283 END SUBROUTINE cg_write_congrad_pio
2284# endif
2285!
2286!***********************************************************************
2287 SUBROUTINE cg_read_congrad (ng, model, outLoop)
2288!***********************************************************************
2289! !
2290! This routine reads conjugate gradient variables for previous outer !
2291! loops using either the standard NetCDF library or the Parallel-IO !
2292! (PIO) library. !
2293! !
2294!***********************************************************************
2295!
2296! Imported variable declarations
2297!
2298 integer, intent(in) :: ng, model, outloop
2299!
2300! Local variable declarations.
2301!
2302 character (len=*), parameter :: myfile = &
2303 & __FILE__//", cg_read_congrad"
2304!
2305 sourcefile=myfile
2306!
2307!-----------------------------------------------------------------------
2308! Read in conjugate gradient variables.
2309!-----------------------------------------------------------------------
2310!
2311 SELECT CASE (dav(ng)%IOtype)
2312 CASE (io_nf90)
2313 CALL cg_read_congrad_nf90 (ng, model, outloop)
2314
2315# if defined PIO_LIB && defined DISTRIBUTE
2316 CASE (io_pio)
2317 CALL cg_read_congrad_pio (ng, model, outloop)
2318# endif
2319 CASE DEFAULT
2320 IF (master) WRITE (stdout,10) dav(ng)%IOtype
2321 exit_flag=3
2322 END SELECT
2323 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2324!
2325 10 FORMAT (' CG_READ_CONGRAD - Illegal output type, io_type = ',i0)
2326!
2327 RETURN
2328 END SUBROUTINE cg_read_congrad
2329!
2330!***********************************************************************
2331 SUBROUTINE cg_read_congrad_nf90 (ng, model, outLoop)
2332!***********************************************************************
2333!
2334 USE mod_netcdf
2335!
2336! Imported variable declarations
2337!
2338 integer, intent(in) :: ng, model, outloop
2339!
2340! Local variable declarations.
2341!
2342 integer :: status
2343!
2344 character (len=*), parameter :: myfile = &
2345 & __FILE__//", cg_read_congrad_nf90"
2346!
2347 sourcefile=myfile
2348!
2349!-----------------------------------------------------------------------
2350! If split 4D-Var and outer>1, Read in conjugate gradient variables
2351! four outerloop restart.
2352!-----------------------------------------------------------------------
2353!
2354! Open DAV NetCDF file for reading.
2355!
2356 IF (dav(ng)%ncid.eq.-1) THEN
2357 CALL netcdf_open (ng, model, dav(ng)%name, 1, dav(ng)%ncid)
2358 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2359 END IF
2360!
2361! Read in number of converget Ritz eigenvalues
2362!
2363 CALL netcdf_get_ivar (ng, model, dav(ng)%name, &
2364 & 'nConvRitz', nconvritz, &
2365 & ncid = dav(ng)%ncid, &
2366 & start = (/outloop-1/), &
2367 & total = (/1/))
2368 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2369!
2370! Read in Ritz eigenvalues from previous outer loop.
2371!
2372 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2373 & 'Ritz', ritz, &
2374 & ncid = dav(ng)%ncid, &
2375 & start = (/1,outloop-1/), &
2376 & total = (/ninner,1/))
2377 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2378!
2379! Read in conjugate gradient "beta" coefficients.
2380!
2381 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2382 & 'cg_beta', cg_beta, &
2383 & ncid = dav(ng)%ncid, &
2384 & start = (/1,1/), &
2385 & total = (/ninner+1,outloop-1/))
2386 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2387!
2388! Read in conjugate gradient "delta" coefficients.
2389!
2390 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2391 & 'cg_delta', cg_delta, &
2392 & ncid = dav(ng)%ncid, &
2393 & start = (/1,1/), &
2394 & total = (/ninner,outloop-1/))
2395 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2396!
2397! Read in normalization coefficients for Lanczos vectos..
2398!
2399 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2400 & 'cg_dla', cg_dla, &
2401 & ncid = dav(ng)%ncid, &
2402 & start = (/1,1/), &
2403 & total = (/ninner,outloop-1/))
2404 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2405!
2406! Read in initial gradient normalization factor.
2407!
2408 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2409 & 'cg_Gnorm_y', cg_gnorm_y, &
2410 & ncid = dav(ng)%ncid, &
2411 & start = (/1/), &
2412 & total = (/outloop-1/))
2413 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2414!
2415 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2416 & 'cg_Gnorm_v', cg_gnorm_v, &
2417 & ncid = dav(ng)%ncid, &
2418 & start = (/1/), &
2419 & total = (/outloop-1/))
2420 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2421!
2422! Read in Lanczos vector normalization factor.
2423!
2424 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2425 & 'cg_QG', cg_qg, &
2426 & ncid = dav(ng)%ncid, &
2427 & start = (/1,1/), &
2428 & total = (/ninner+1,outloop-1/))
2429 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2430!
2431! Read in eigenvalues of Lanczos recurrence relationship.
2432!
2433 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2434 & 'cg_Ritz', cg_ritz, &
2435 & ncid = dav(ng)%ncid, &
2436 & start = (/1,1/), &
2437 & total = (/ninner,outloop-1/))
2438 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2439!
2440! Read eigenvalues relative error.
2441!
2442 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2443 & 'cg_RitzErr', cg_ritzerr, &
2444 & ncid = dav(ng)%ncid, &
2445 & start = (/1,1/), &
2446 & total = (/ninner,outloop-1/))
2447 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2448!
2449! Read in eigenvectors of Lanczos recurrence relationship for previous
2450! outer loop.
2451!
2452 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2453 & 'cg_zv', cg_zv, &
2454 & ncid = dav(ng)%ncid, &
2455 & start = (/1,1,outloop-1/), &
2456 & total = (/ninner,ninner,1/))
2457 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2458!
2459! Read first-guess estimate of representer coefficients for previous
2460! outer loop.
2461!
2462 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2463 & 'cg_pxsave', cg_pxsave, &
2464 & ncid = dav(ng)%ncid, &
2465 & start = (/1,outloop-1/), &
2466 & total = (/ndatum(ng),1/))
2467 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2468!
2469! Read in initial gradient for minimization.
2470!
2471 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2472 & 'zgrad0', zgrad0, &
2473 & ncid = dav(ng)%ncid, &
2474 & start = (/1,1/), &
2475 & total = (/ndatum(ng),outloop-1/))
2476 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2477!
2478! Read in Lanczos vectors.
2479!
2480 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2481 & 'zcglwk', zcglwk, &
2482 & ncid = dav(ng)%ncid, &
2483 & start = (/1,1,1/), &
2484 & total = (/ndatum(ng),ninner+1,outloop-1/))
2485 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2486!
2487 RETURN
2488 END SUBROUTINE cg_read_congrad_nf90
2489
2490# if defined PIO_LIB && defined DISTRIBUTE
2491!
2492!***********************************************************************
2493 SUBROUTINE cg_read_congrad_pio (ng, model, outLoop)
2494!***********************************************************************
2495!
2496 USE mod_pio_netcdf
2497!
2498! Imported variable declarations
2499!
2500 integer, intent(in) :: ng, model, outLoop
2501!
2502! Local variable declarations.
2503!
2504 integer :: status
2505!
2506 character (len=*), parameter :: MyFile = &
2507 & __FILE__//", cg_read_congrad_pio"
2508!
2509 sourcefile=myfile
2510!
2511!-----------------------------------------------------------------------
2512! If split 4D-Var and outer>1, Read in conjugate gradient variables
2513! four outerloop restart.
2514!-----------------------------------------------------------------------
2515!
2516! Open DAV NetCDF file for reading.
2517!
2518 IF (dav(ng)%pioFile%fh.eq.-1) THEN
2519 CALL pio_netcdf_open (ng, model, dav(ng)%name, 1, &
2520 & dav(ng)%pioFile)
2521 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2522 END IF
2523!
2524! Read in number of converget Ritz eigenvalues
2525!
2526 CALL pio_netcdf_get_ivar (ng, model, dav(ng)%name, &
2527 & 'nConvRitz', nconvritz, &
2528 & piofile = dav(ng)%pioFile, &
2529 & start = (/outloop-1/), &
2530 & total = (/1/))
2531 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2532!
2533! Read in Ritz eigenvalues from previous outer loop.
2534!
2535 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2536 & 'Ritz', ritz, &
2537 & piofile = dav(ng)%pioFile, &
2538 & start = (/1,outloop-1/), &
2539 & total = (/ninner,1/))
2540 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2541!
2542! Read in conjugate gradient "beta" coefficients.
2543!
2544 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2545 & 'cg_beta', cg_beta, &
2546 & piofile = dav(ng)%pioFile, &
2547 & start = (/1,1/), &
2548 & total = (/ninner+1,outloop-1/))
2549 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2550!
2551! Read in conjugate gradient "delta" coefficients.
2552!
2553 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2554 & 'cg_delta', cg_delta, &
2555 & piofile = dav(ng)%pioFile, &
2556 & start = (/1,1/), &
2557 & total = (/ninner,outloop-1/))
2558 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2559!
2560! Read in normalization coefficients for Lanczos vectos..
2561!
2562 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2563 & 'cg_dla', cg_dla, &
2564 & piofile = dav(ng)%pioFile, &
2565 & start = (/1,1/), &
2566 & total = (/ninner,outloop-1/))
2567 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2568!
2569! Read in initial gradient normalization factor.
2570!
2571 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2572 & 'cg_Gnorm_y', cg_gnorm_y, &
2573 & piofile = dav(ng)%pioFile, &
2574 & start = (/1/), &
2575 & total = (/outloop-1/))
2576 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2577!
2578 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2579 & 'cg_Gnorm_v', cg_gnorm_v, &
2580 & piofile = dav(ng)%pioFile, &
2581 & start = (/1/), &
2582 & total = (/outloop-1/))
2583 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2584!
2585! Read in Lanczos vector normalization factor.
2586!
2587 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2588 & 'cg_QG', cg_qg, &
2589 & piofile = dav(ng)%pioFile, &
2590 & start = (/1,1/), &
2591 & total = (/ninner+1,outloop-1/))
2592 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2593!
2594! Read in eigenvalues of Lanczos recurrence relationship.
2595!
2596 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2597 & 'cg_Ritz', cg_ritz, &
2598 & piofile = dav(ng)%pioFile, &
2599 & start = (/1,1/), &
2600 & total = (/ninner,outloop-1/))
2601 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2602!
2603! Read eigenvalues relative error.
2604!
2605 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2606 & 'cg_RitzErr', cg_ritzerr, &
2607 & piofile = dav(ng)%pioFile, &
2608 & start = (/1,1/), &
2609 & total = (/ninner,outloop-1/))
2610 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2611!
2612! Read in eigenvectors of Lanczos recurrence relationship for previous
2613! outer loop.
2614!
2615 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2616 & 'cg_zv', cg_zv, &
2617 & piofile = dav(ng)%pioFile, &
2618 & start = (/1,1,outloop-1/), &
2619 & total = (/ninner,ninner,1/))
2620 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2621!
2622! Read first-guess estimate of representer coefficients for previous
2623! outer loop.
2624!
2625 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2626 & 'cg_pxsave', cg_pxsave, &
2627 & piofile = dav(ng)%pioFile, &
2628 & start = (/1,outloop-1/), &
2629 & total = (/ndatum(ng),1/))
2630 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2631!
2632! Read in initial gradient for minimization.
2633!
2634 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2635 & 'zgrad0', zgrad0, &
2636 & piofile = dav(ng)%pioFile, &
2637 & start = (/1,1/), &
2638 & total = (/ndatum(ng),outloop-1/))
2639 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2640!
2641! Read in Lanczos vectors.
2642!
2643 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2644 & 'zcglwk', zcglwk, &
2645 & piofile = dav(ng)%pioFile, &
2646 & start = (/1,1,1/), &
2647 & total = (/ndatum(ng),ninner+1, &
2648 & outloop-1/))
2649 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2650!
2651 RETURN
2652 END SUBROUTINE cg_read_congrad_pio
2653# endif
2654#endif
2655 END MODULE congrad_mod
congrad_mod::rprecond
Definition congrad.F:136
distribute_mod::mp_bcastf
Definition distribute.F:75
distribute_mod::mp_bcasti
Definition distribute.F:95
distribute_mod::mp_bcastl
Definition distribute.F:89
mod_netcdf::netcdf_get_fvar
Definition mod_netcdf.F:45
mod_netcdf::netcdf_get_ivar
Definition mod_netcdf.F:59
mod_netcdf::netcdf_put_fvar
Definition mod_netcdf.F:82
mod_netcdf::netcdf_put_ivar
Definition mod_netcdf.F:96
mod_pio_netcdf::pio_netcdf_get_fvar
Definition mod_pio_netcdf.F:52
mod_pio_netcdf::pio_netcdf_get_ivar
Definition mod_pio_netcdf.F:66
mod_pio_netcdf::pio_netcdf_put_fvar
Definition mod_pio_netcdf.F:94
mod_pio_netcdf::pio_netcdf_put_ivar
Definition mod_pio_netcdf.F:108
congrad_mod
Definition congrad.F:2
congrad_mod::rpevecs
subroutine rpevecs(ng, outloop, ninnloop)
Definition congrad.F:1282
congrad_mod::cg_write_congrad_pio
subroutine cg_write_congrad_pio(ng, model, innloop, outloop, ninnloop, jf, jdata, jmod, jopt, jb, jobs, jact, preducv, preducy)
Definition congrad.F:2018
congrad_mod::cg_write_congrad
subroutine, private cg_write_congrad(ng, model, innloop, outloop, ninnloop, jf, jdata, jmod, jopt, jb, jobs, jact, preducv, preducy)
Definition congrad.F:1685
congrad_mod::cg_write_congrad_nf90
subroutine, private cg_write_congrad_nf90(ng, model, innloop, outloop, ninnloop, jf, jdata, jmod, jopt, jb, jobs, jact, preducv, preducy)
Definition congrad.F:1744
congrad_mod::congrad
subroutine, public congrad(ng, model, outloop, innloop, ninnloop, lcgini)
Definition congrad.F:163
congrad_mod::rprecond_r8
subroutine rprecond_r8(ng, lscale, ltrans, outloop, ninnloop, py)
Definition congrad.F:1517
congrad_mod::cg_read_congrad_nf90
subroutine, private cg_read_congrad_nf90(ng, model, outloop)
Definition congrad.F:2332
congrad_mod::cg_read_congrad
subroutine, public cg_read_congrad(ng, model, outloop)
Definition congrad.F:2288
congrad_mod::cg_read_congrad_pio
subroutine cg_read_congrad_pio(ng, model, outloop)
Definition congrad.F:2494
congrad_mod::rprecond_dp
subroutine rprecond_dp(ng, lscale, ltrans, outloop, ninnloop, py)
Definition congrad.F:1351
distribute_mod
Definition distribute.F:3
lapack_mod
Definition lapack_mod.F:2
lapack_mod::dsteqr
subroutine, public dsteqr(compz, n, d, e, z, ldz, work, info)
Definition lapack_mod.F:66
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::gdgw
real(r8), dimension(:), allocatable gdgw
Definition mod_fourdvar.F:297
mod_fourdvar::cg_beta
real(dp), dimension(:,:), allocatable cg_beta
Definition mod_fourdvar.F:654
mod_fourdvar::cg_ritz
real(dp), dimension(:,:), allocatable cg_ritz
Definition mod_fourdvar.F:716
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::ritz
real(dp), dimension(:), allocatable ritz
Definition mod_fourdvar.F:664
mod_fourdvar::hevecerr
real(r8) hevecerr
Definition mod_fourdvar.F:858
mod_fourdvar::obsval
real(r8), dimension(:), allocatable obsval
Definition mod_fourdvar.F:232
mod_fourdvar::cg_innov
real(r8), dimension(:), allocatable cg_innov
Definition mod_fourdvar.F:299
mod_fourdvar::cg_ritzerr
real(dp), dimension(:,:), allocatable cg_ritzerr
Definition mod_fourdvar.F:717
mod_fourdvar::lprecond
logical lprecond
Definition mod_fourdvar.F:821
mod_fourdvar::lritz
logical lritz
Definition mod_fourdvar.F:827
mod_fourdvar::vcglev
real(r8), dimension(:,:,:), allocatable vcglev
Definition mod_fourdvar.F:293
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_zv
real(dp), dimension(:,:,:), allocatable cg_zv
Definition mod_fourdvar.F:669
mod_fourdvar::lhessianev
logical lhessianev
Definition mod_fourdvar.F:812
mod_fourdvar::nritzev
integer nritzev
Definition mod_fourdvar.F:841
mod_fourdvar::vgrad0s
real(r8), dimension(:), allocatable vgrad0s
Definition mod_fourdvar.F:296
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::ritzmaxerr
real(dp) ritzmaxerr
Definition mod_fourdvar.F:659
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::cg_gnorm_y
real(dp), dimension(:), allocatable cg_gnorm_y
Definition mod_fourdvar.F:698
mod_fourdvar::jobs
real(r8), dimension(:,:), allocatable jobs
Definition mod_fourdvar.F:255
mod_fourdvar::zgrad0
real(r8), dimension(:,:), allocatable zgrad0
Definition mod_fourdvar.F:294
mod_fourdvar::nconvritz
integer nconvritz
Definition mod_fourdvar.F:842
mod_iounits
Definition mod_iounits.F:2
mod_iounits::dav
type(t_io), dimension(:), allocatable dav
Definition mod_iounits.F:205
mod_iounits::stdout
integer stdout
Definition mod_iounits.F:195
mod_iounits::sourcefile
character(len=256) sourcefile
Definition mod_iounits.F:317
mod_ncparam
Definition mod_ncparam.F:2
mod_ncparam::io_nf90
integer, parameter io_nf90
Definition mod_ncparam.F:95
mod_ncparam::io_pio
integer, parameter io_pio
Definition mod_ncparam.F:96
mod_netcdf
Definition mod_netcdf.F:2
mod_netcdf::netcdf_open
subroutine, public netcdf_open(ng, model, ncname, omode, ncid)
Definition mod_netcdf.F:9293
mod_netcdf::netcdf_sync
subroutine, public netcdf_sync(ng, model, ncname, ncid)
Definition mod_netcdf.F:9505
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_pio_netcdf
Definition mod_pio_netcdf.F:2
mod_pio_netcdf::pio_netcdf_sync
subroutine, public pio_netcdf_sync(ng, model, ncname, piofile)
Definition mod_pio_netcdf.F:9275
mod_pio_netcdf::pio_netcdf_open
subroutine, public pio_netcdf_open(ng, model, ncname, omode, piofile)
Definition mod_pio_netcdf.F:9148
mod_scalars
Definition mod_scalars.F:2
mod_scalars::ninner
integer ninner
Definition mod_scalars.F:353
mod_scalars::nouter
integer nouter
Definition mod_scalars.F:354
mod_scalars::nconv
integer, dimension(:), allocatable nconv
Definition mod_scalars.F:798
mod_scalars::exit_flag
integer exit_flag
Definition mod_scalars.F:559
mod_scalars::inner
integer inner
Definition mod_scalars.F:360
mod_scalars::noerror
integer noerror
Definition mod_scalars.F:561
mod_scalars::outer
integer outer
Definition mod_scalars.F:361
strings_mod
Definition strings.F:2
strings_mod::founderror
logical function, public founderror(flag, noerr, line, routine)
Definition strings.F:52
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3