9
10
11
12
13
14
15
16
17
18
19
20# if defined R4DVAR_ANA_SENSITIVITY || defined TL_R4DVAR
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61# elif defined RBL4DVAR_ANA_SENSITIVITY || \
62 defined rbl4dvar_fct_sensitivity || \
63 defined tl_rbl4dvar
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102# endif
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
128
129# ifdef DISTRIBUTE
130
132# endif
133 implicit none
134
135
136
137 logical, intent(in) :: Lcgini
138 integer, intent(in) :: ng, model, outLoop, innLoop, NinnLoop
139
140
141
142 logical :: Ltrans, Laug
143 integer :: i, ic, j, iobs, ivec, Lscale, info
144
145 real(r8) :: zbet, eps, preducv, preducy
146 real(r8) :: Jopt, Jf, Jmod, Jdata, Jb, Jobs, Jact, cff
147 real(r8) :: tl_dla, fact
148
149 real(r8), dimension(NinnLoop) :: zu, zgam, zfact, tl_zfact
150 real(r8), dimension(NinnLoop) :: tl_zu, tl_zrhs
151 real(r8), dimension(Ndatum(ng)+1) :: px, pgrad
152 real(r8), dimension(Ndatum(ng)+1) :: tl_px, tl_pgrad
153 real(r8), dimension(Ninner,3) :: zwork, tl_zwork
154
155 character (len=13) :: string
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173 ltrans=.false.
174
175 master_thread :
IF (
master)
THEN
176
177
178
181 END DO
182
183
184
185
186 IF (innloop.gt.0) THEN
188
189
191 & tlmodval(iobs)+ &
193 & tlmodval_s(iobs,innloop,outloop)
194 END DO
195 END IF
196
197
198
199
200
201 minimize : IF (innloop.eq.0) THEN
202
203# if defined RBL4DVAR_ANA_SENSITIVITY || \
204 defined rbl4dvar_fct_sensitivity || \
205 defined tl_rbl4dvar || \
206 defined tl_r4dvar
207
210 tl_bckmodval(iobs)=0.0_r8
211
212
213
214
215
216
217 END DO
218# endif
219
220
221
222 IF (outloop.eq.1) THEN
224
225
226 tl_px(iobs)=0.0_r8
227 END DO
229
230
231 tl_hbk(iobs)=0.0_r8
232 END DO
233 ELSE
234
235
236
237
239
240
241
242 tl_hbk(iobs)=-tlmodval(iobs)
243 END DO
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258 END IF
259
260
261
262
263
265
266
267
268# if defined RBL4DVAR || \
269 defined rbl4dvar_ana_sensitivity || \
270 defined rbl4dvar_fct_sensitivity || \
271 defined tl_rbl4dvar
272
273
274
277 & tl_bckmodval(iobs))+ &
279 & (
obsval(iobs)-bckmodval(iobs))
280# else
281
282
283
284
285
286
287
288
289
290
293 & tlmodval(iobs))+ &
296 & tlmodval_s(iobs,innloop,outloop))
297# endif
298 IF (
obserr(iobs).NE.0.0_r8)
THEN
299
300
302 tl_pgrad(iobs)=tl_pgrad(iobs)/
obserr(iobs)- &
304 END IF
305
306
308
309
311 END DO
312
313
314
315
316
317 pgrad(
ndatum(ng)+1)=1.0_r8
318 tl_pgrad(
ndatum(ng)+1)=0.0_r8
319
320
322
323
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
343
344
346 END DO
347
348
349
351
352
353
355 END DO
356
357 preducv=1.0_r8
358 preducy=1.0_r8
359
360
361
362
363
364 ELSE
365
366
367
368
369
370
371
372
373
374 IF (innloop.EQ.1) THEN
375
376
379
380
381
383 & tlmodval(iobs)* &
385 & tlmodval_s(iobs,innloop,outloop)* &
387 END DO
388
389
390
391
392
394
395
396
397
398
399
400
401
403 & tl_hbk(iobs)* &
406 & hbk(iobs,outloop)* &
409 & hbk(iobs,outloop)* &
412 & tl_hbk(iobs)* &
415 & hbk(iobs,outloop)* &
418 & hbk(iobs,outloop)* &
421 END DO
422
423
424
425
426
437
438
439
442
443
444
445
446
447
449
450
451
452
453
455
456
457
458
459
460
464 &
zcglwk(iobs,1,outloop)/ &
466
467
468
469
470
471
475 &
vcglwk(iobs,1,outloop)/ &
477 END DO
478
479
480
482
483
484
485
486
487 tl_tlmodval_s(iobs,innloop,outloop)=tlmodval(iobs)
488
489
490
491 tlmodval(iobs)=tlmodval(iobs)/ &
494 & tlmodval_s(iobs,innloop,outloop)/ &
497 END DO
498
499 ELSE
500
501
502
503
504
506
507
508
509
510
511
512
513 pgrad(iobs)=
zcglwk(iobs,innloop,outloop)* &
516 END DO
517
518
519
520
521
524
525
526
527
529 & tl_pgrad(iobs)* &
530 & tlmodval_s(iobs,innloop,outloop)+ &
531 & pgrad(iobs)* &
532 & tlmodval(iobs)
533 END DO
534
535
536
538
539
540
541
542
543
544
545
546
548 & tl_hbk(iobs)* &
550 & pgrad(iobs)+ &
551 & hbk(iobs,outloop)* &
553 & pgrad(iobs)+ &
554 & hbk(iobs,outloop)* &
556 & tl_pgrad(iobs)+ &
557 & tl_hbk(iobs)* &
559 &
vcglwk(iobs,innloop,outloop)+ &
560 & hbk(iobs,outloop)* &
561 & tl_pgrad(
ndatum(ng)+1)* &
562 &
vcglwk(iobs,innloop,outloop)+ &
563 & hbk(iobs,outloop)* &
566 END DO
567
568
569
570
571
572
578 & tl_pgrad(
ndatum(ng)+1)* &
583
584
585
588
589
590
591
592
594
595
596
597
598
599
600 tl_zcglwk(iobs,innloop)=tl_pgrad(iobs)/ &
603 &
zcglwk(iobs,innloop,outloop)/ &
605
606
607
608
609
610
611
615 &
vcglwk(iobs,innloop,outloop)/ &
617 END DO
618
619
620
621
622
624
625
626 tl_tlmodval_s(iobs,innloop,outloop)=tlmodval(iobs)
627
628
629 tlmodval(iobs)=tlmodval(iobs)/
cg_beta(innloop,outloop)- &
631 & tlmodval_s(iobs,innloop,outloop)/ &
634 END DO
635
636
637
638
639
642
643
644
645
647 & tlmodval(iobs)* &
649 & tlmodval_s(iobs,innloop,outloop)* &
652 END DO
653
654
655
657
658
659
660
661
662
663
664
665
667 & tl_hbk(iobs)* &
670 & hbk(iobs,outloop)* &
673 & hbk(iobs,outloop)* &
676 & tl_hbk(iobs)* &
677 &
vcglwk(iobs,innloop,outloop)* &
679 & hbk(iobs,outloop)* &
682 & hbk(iobs,outloop)* &
683 &
vcglwk(iobs,innloop,outloop)* &
685 END DO
686
687
688
689
690
691
702
703
704
705
706 IF (innloop.eq.ninnloop) THEN
707
708
709
710
711
712 IF (ninnloop.eq.1) THEN
713 print *,'Illegal configuration!'
714 print *,'Ninner must be ge 2'
715 stop
716 END IF
717 IF (ninnloop.eq.2) THEN
720 tl_zu(1)=tl_zrhs(1)/
cg_delta(1,outloop)
721 zwork(1,3)=zu(1)
722 tl_zwork(1,3)=tl_zu(1)
723 ELSE
724
725
726
728 zu(1)=
cg_qg(1,outloop)/zbet
729 DO ivec=2,innloop-1
730 zgam(ivec)=
cg_beta(ivec,outloop)/zbet
732 &
cg_beta(ivec,outloop)*zgam(ivec)
733 zu(ivec)=(
cg_qg(ivec,outloop)- &
735 & zu(ivec-1))/zbet
736 END DO
737 zwork(innloop-1,3)=zu(innloop-1)
738
739 DO ivec=innloop-2,1,-1
740 zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
741 zwork(ivec,3)=zu(ivec)
742 END DO
743
744
745
749 DO ivec=2,innloop-2
754 END DO
755 tl_zrhs(innloop-1)=
tl_cg_qg(innloop-1)- &
758
759
760
762 tl_zu(1)=tl_zrhs(1)/zbet
763 DO ivec=2,innloop-1
764 zgam(ivec)=
cg_beta(ivec,outloop)/zbet
766 &
cg_beta(ivec,outloop)*zgam(ivec)
767 tl_zu(ivec)=(tl_zrhs(ivec)- &
769 & tl_zu(ivec-1))/zbet
770 END DO
771 tl_zwork(innloop-1,3)=tl_zu(innloop-1)
772
773 DO ivec=innloop-2,1,-1
774
775
776 tl_zu(ivec)=tl_zu(ivec)-zgam(ivec+1)*tl_zu(ivec+1)
777 tl_zwork(ivec,3)=tl_zu(ivec)
778 END DO
779 END IF
780
781
782
784 px(iobs)=0.0_r8
785 tl_px(iobs)=0.0_r8
786 DO ivec=1,innloop-1
787 px(iobs)=px(iobs)+ &
788 &
zcglwk(iobs,ivec,outloop)*zwork(ivec,3)
789 tl_px(iobs)=tl_px(iobs)+ &
791 & zwork(ivec,3)+ &
792 &
zcglwk(iobs,ivec,outloop)* &
793 & tl_zwork(ivec,3)
794 END DO
795 END DO
796
797
798
799
800
801
802
803
804
805
806
807 END IF
808
809 END IF
810
811
812
813
814
815
816
817
818 IF (innloop.eq.1) THEN
820 ELSE
821 fact=1.0_r8/
cg_beta(innloop,outloop)
822 END IF
823
825
826
827
828
829
830
831
832 pgrad(iobs)=tlmodval_s(iobs,innloop,outloop)* &
833 & fact+ &
834 & hbk(iobs,outloop)* &
836 tl_pgrad(iobs)=tlmodval(iobs)+ &
837 & tl_hbk(iobs)* &
839 & hbk(iobs,outloop)* &
841 END DO
842
844
845
846
847 IF (
obserr(iobs).ne.0.0_r8)
THEN
848
849
850
851 pgrad(iobs)=pgrad(iobs)/
obserr(iobs)
852 tl_pgrad(iobs)=tl_pgrad(iobs)/
obserr(iobs)- &
854 END IF
855 END DO
856 pgrad(
ndatum(ng)+1)=0.0_r8
857 tl_pgrad(
ndatum(ng)+1)=0.0_r8
858
859
860
862 pgrad(iobs)=pgrad(iobs)+ &
863 &
vcglwk(iobs,innloop,outloop)
864 tl_pgrad(iobs)=tl_pgrad(iobs)+ &
866 END DO
871
872 IF (innloop.gt.1) THEN
873
874
875
876
877
879 pgrad(iobs)=pgrad(iobs)- &
881 &
zcglwk(iobs,innloop-1,outloop)
882 tl_pgrad(iobs)=tl_pgrad(iobs)- &
884 &
zcglwk(iobs,innloop-1,outloop)- &
887 END DO
888 END IF
889
890
891
893
894
895
896
897
898 IF (innloop.eq.1) THEN
900 ELSE
901 fact=1.0_r8/
cg_beta(innloop,outloop)
902 END IF
903
905
906
907
908
909
910
911
912
913
914
915
917 & tlmodval(iobs)* &
918 & pgrad(iobs)+ &
919 & tlmodval_s(iobs,innloop,outloop)* &
920 & tl_pgrad(iobs)* &
921 & fact+ &
922 & tl_hbk(iobs)* &
924 & pgrad(iobs)+ &
925 & hbk(iobs,outloop)* &
927 & pgrad(iobs)+ &
928 & hbk(iobs,outloop)* &
930 & tl_pgrad(iobs)+ &
931 & tl_hbk(iobs)* &
932 &
vcglwk(iobs,innloop,outloop)* &
934 & hbk(iobs,outloop)* &
937 & hbk(iobs,outloop)* &
938 &
vcglwk(iobs,innloop,outloop)* &
940 END DO
941
942
943
944
945
946
957
958
959
960
961
963 pgrad(iobs)=pgrad(iobs)- &
965 &
zcglwk(iobs,innloop,outloop)
966 tl_pgrad(iobs)=tl_pgrad(iobs)- &
968 &
zcglwk(iobs,innloop,outloop)- &
971 END DO
972
973
974
975
976
977 DO ivec=1,innloop
978 IF (ivec.eq.1) THEN
981 & zfact(ivec)/ &
983 ELSE
984 zfact(ivec)=1.0_r8/
cg_beta(ivec,outloop)
986 & zfact(ivec)/ &
988 ENDIF
989 END DO
990
991
992
993
994
995 DO ivec=innloop,1,-1
996
997
998 tl_dla=0.0_r8
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011 tl_dla=tl_dla+ &
1012 & tl_pgrad(iobs)* &
1013 & tlmodval_s(iobs,ivec,outloop)* &
1014 & zfact(ivec)+ &
1015 & pgrad(iobs)* &
1016 & tl_tlmodval_s(iobs,ivec,outloop)* &
1017 & zfact(ivec)+ &
1018 & pgrad(iobs)* &
1019 & tlmodval_s(iobs,ivec,outloop)* &
1020 & tl_zfact(ivec)+ &
1021 & tl_pgrad(iobs)* &
1022 & hbk(iobs,outloop)* &
1024 & pgrad(iobs)* &
1025 & tl_hbk(iobs)* &
1027 & pgrad(iobs)* &
1028 & hbk(iobs,outloop)* &
1030 & tl_hbk(iobs)* &
1031 &
vcglwk(iobs,ivec,outloop)* &
1033 & hbk(iobs,outloop)* &
1036 & hbk(iobs,outloop)* &
1037 &
vcglwk(iobs,ivec,outloop)* &
1039 END DO
1040
1041
1042
1043
1044
1045 tl_dla=tl_dla+ &
1056 pgrad(iobs)=pgrad(iobs)- &
1057 &
cg_dla(ivec,outloop)* &
1058 &
vcglwk(iobs,ivec,outloop)
1059 tl_pgrad(iobs)=tl_pgrad(iobs)-tl_dla* &
1060 &
vcglwk(iobs,ivec,outloop)- &
1061 &
cg_dla(ivec,outloop)* &
1063 END DO
1064 END DO
1065
1066
1067
1068
1069
1070
1071
1072
1073
1075
1076
1077 tl_zcglwk(iobs,innloop+1)=tl_pgrad(iobs)
1078 END DO
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1096
1097
1098 tl_vcglwk(iobs,innloop+1)=tl_pgrad(iobs)
1099 END DO
1100
1101
1102
1103
1104
1105 IF (innloop.eq.ninnloop) THEN
1106
1107
1108
1109
1111
1112
1113
1115 END DO
1117
1118
1119
1120
1121 fourdvar(ng)%tl_cg_pxsave(iobs)=tl_px(iobs)
1122
1123
1124 END DO
1125 ELSE
1126
1127
1128
1130
1131
1132
1134 END DO
1135
1136 END IF
1137
1138 END IF minimize
1139
1140 END IF master_thread
1141
1142# ifdef DISTRIBUTE
1143
1144
1145
1146
1147
1153# if defined RBL4DVAR || defined R4DVAR || \
1154 defined sensitivity_4dvar || defined tl_rbl4dvar || \
1155 defined tl_r4dvar
1158# endif
1167# endif
1168
1169 RETURN
real(r8), dimension(:,:), allocatable tl_vcglwk
real(dp), dimension(:), allocatable cg_gnorm_v
real(r8), dimension(:,:,:), allocatable vcglwk
type(t_fourdvar), dimension(:), allocatable fourdvar
real(dp), dimension(:,:), allocatable cg_beta
real(r8), dimension(:), allocatable tl_obsval
integer, dimension(:), allocatable ndatum
real(r8), dimension(:,:), allocatable cg_dla
real(dp), dimension(:,:), allocatable cg_qg
real(r8), dimension(:), allocatable tl_vgrad0s
real(r8), dimension(:), allocatable obsval
real(dp), dimension(:), allocatable tl_cg_qg
real(r8), dimension(:), allocatable obsscale
real(r8), dimension(:), allocatable obserr
real(r8), dimension(:), allocatable tl_zgrad0
real(r8), dimension(:), allocatable jb0
real(dp), dimension(:), allocatable tl_cg_delta
real(r8), dimension(:,:), allocatable tl_zcglwk
real(dp), dimension(:), allocatable tl_cg_beta
real(r8), dimension(:), allocatable tl_obsscale
real(r8), dimension(:), allocatable admodval
real(r8), dimension(:,:,:), allocatable zcglwk
real(r8), dimension(:), allocatable tl_jb0
real(r8), dimension(:), allocatable vgrad0
real(r8), dimension(:), allocatable nlmodval
real(dp), dimension(:,:), allocatable cg_delta
real(r8), dimension(:), allocatable tl_obserr
real(r8), dimension(:), allocatable tl_vgrad0
real(dp), dimension(:), allocatable tl_cg_gnorm_v
real(r8), dimension(:,:), allocatable zgrad0