163
164
165
166
167 logical, intent(in) :: Lcgini
168
169 integer, intent(in) :: ng, model, outLoop, innLoop, NinnLoop
170
171
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
204
205
206# ifdef PROFILE
207
208
209
210 CALL wclock_on (ng, model, 85, __line__, myfile)
211# endif
212
213
214
215
216
217
218
219
220
221
222
223
224
225 ritzmaxerr=hevecerr
226 ltrans=.false.
227
228 master_thread : IF (master) THEN
229
230
231
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
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
255
256 DO iobs=1,ndatum(ng)
257 bckmodval(iobs)=nlmodval(iobs)
258 END DO
259# endif
260
261
262
263 IF ((outloop.eq.1).or.(.not.lhotstart)) THEN
264
265
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
273
274
275
276
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
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
299
300 IF (lprecond.and.(outloop.gt.1)) THEN
301 lscale=2
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
325
326 IF (lprecond.and.(outloop.gt.1)) THEN
327 lscale=2
328 ltrans=.false.
329 CALL rprecond (ng, lscale, ltrans, outloop, ninnloop, &
330 & admodval)
331 END IF
332
333
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
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
365
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
381
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
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
405
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
413
414 IF (lprecond.and.(outloop.gt.1)) THEN
415 lscale=2
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
440
441 IF (lprecond.and.(outloop.gt.1)) THEN
442 lscale=2
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
469
470
471 ELSE
472
473
474
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
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
496
497 IF (lprecond.and.(outloop.gt.1)) THEN
498 lscale=2
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
508
509 IF (lprecond.and.(outloop.gt.1)) THEN
510 lscale=2
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
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
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
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
580
581
582
583
584
585
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
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
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
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
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
654
655 DO i=1,innloop
656 zwork(i,3)=ze(i)
657 END DO
658# else
659
660
661
662
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
696
697
698 IF (lprecond.and.(outloop.gt.1)) THEN
699 lscale=2
700 ltrans=.false.
701 CALL rprecond (ng, lscale, ltrans, outloop, ninnloop, px)
702 END IF
703
704
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
713
714 IF (lprecond.and.(outloop.gt.1)) THEN
715 lscale=-2
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
727
728
729
730 eps=abs(cg_beta(innloop+1,outloop)*zwork(innloop,3))
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
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
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
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
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
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
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
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
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
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
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
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
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
986
987
988
989 IF (innloop.eq.ninnloop) THEN
990
991
992
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
1011
1012 IF (lprecond.and.(outloop.gt.1)) THEN
1013 lscale=2
1014 ltrans=.false.
1015 CALL rprecond (ng, lscale, ltrans, outloop, ninnloop, &
1016 & admodval)
1017 END IF
1018 END IF
1019
1020
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
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
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
1096
1097
1098
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
1110
1111
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
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
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
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
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
1173
1174 CALL rpevecs (ng, outloop, ninnloop)
1175 END IF
1176 END IF
1177
1178 END IF master_thread
1179
1180
1181
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
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
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
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
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
recursive subroutine wclock_off(ng, model, region, line, routine)
recursive subroutine wclock_on(ng, model, region, line, routine)