634
635
636
637
638
639
640
641
642
643
644
645
646
647# ifdef RPCG
648
653# endif
654
655
656
657 integer, intent(in) :: my_outer
658
659 real(dp), intent(in) :: RunInterval
660
661
662
663 logical :: Lcgini, Linner, Lposterior
664
665 integer :: i, ifile, lstr, my_inner, ng, tile
666 integer :: Fcount, InpRec
667# ifdef RPCG
668 integer :: ADrec, nLAST
669 integer :: irec, jrec, jrec1, jrec2
670 integer :: LiNL
671# endif
672# ifdef PROFILE
673 integer :: thread
674# endif
675
676 integer, dimension(Ngrids) :: Nrec
677# ifdef RPCG
678 integer, dimension(Ngrids) :: nADrec
679# endif
680# ifdef SPLIT_4DVAR
681
682 real(dp) :: stime
683# endif
684
685 character (len=8 ) :: driver = 'rbl4dvar'
686 character (len=10) :: suffix
687
688 character (len=*), parameter :: MyFile = &
689 & __FILE__//", increment"
690
691 sourcefile=myfile
692
693
694
695
696
697# ifdef PROFILE
698
699
700
701 DO ng=1,ngrids
702 DO thread=thread_range
703 CALL wclock_on (ng, itlm, 87, __line__, myfile)
704 END DO
705 END DO
706# endif
707
708# ifdef SPLIT_4DVAR
709
710
711
712
713
714
715
716
717 IF (my_outer.gt.1) THEN
718 nrun=1+(my_outer-1)*ninner
719 END IF
720
721
722
723
724 DO ng=1,ngrids
725 ldefini(ng)=.false.
726 CALL def_ini (ng)
727 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
728 END DO
729
730
731
732
733 IF (my_outer.gt.1) THEN
734 DO ng=1,ngrids
735 ldefitl(ng)=.false.
736 CALL tl_def_ini (ng)
737 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
738 END DO
739 END IF
740
741
742
743
744 IF (my_outer.gt.1) THEN
745 DO ng=1,ngrids
746 ldeftlf(ng)=.false.
747 CALL def_impulse (ng)
748 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
749 END DO
750 END IF
751
752
753
754
755 IF (my_outer.gt.1) THEN
756 DO ng=1,ngrids
757 ldefadj(ng)=.false.
758 CALL ad_def_his (ng, ldefadj(ng))
759 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
760 fcount=adm(ng)%Fcount
761 frcrec(ng)=adm(ng)%Nrec(fcount)
762 END DO
763 END IF
764
765
766
767
768
769
770 DO ng=1,ngrids
771 ldefmod(ng)=.false.
772 CALL def_mod (ng)
773 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
774 havenlmod(ng)=.false.
775 END DO
776
777
778
779
780
781 IF (my_outer.gt.1) THEN
782 DO ng=1,ngrids
783# ifdef RPCG
784 CALL cg_read_rpcg (ng, itlm, my_outer)
785# else
786 CALL cg_read_congrad (ng, itlm, my_outer)
787# endif
788 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
789 END DO
790 END IF
791
792
793
794
795
796
797
798
799 IF (my_outer.eq.1) THEN
800 ng=1
801 SELECT CASE (dav(ng)%IOtype)
802 CASE (io_nf90)
803 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
804 & vname(1,idnlmi), nlmodval, &
805 & ncid = dav(ng)%ncid)
806
807# if defined PIO_LIB && defined DISTRIBUTE
808 CASE (io_pio)
809 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
810 & vname(1,idnlmi), nlmodval, &
811 & piofile = dav(ng)%pioFile)
812
813# endif
814 END SELECT
815 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
816 END IF
817
818 ng=1
819 SELECT CASE (dav(ng)%IOtype)
820 CASE (io_nf90)
821 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
822 & vname(1,idobss), obsscale, &
823 & ncid = dav(ng)%ncid)
824
825# if defined PIO_LIB && defined DISTRIBUTE
826 CASE (io_pio)
827 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
828 & vname(1,idobss), obsscale, &
829 & piofile = dav(ng)%pioFile)
830# endif
831 END SELECT
832 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
833
834
835
836
837
838
839 ng=1
840 SELECT CASE (obs(ng)%IOtype)
841 CASE (io_nf90)
842 CALL netcdf_get_fvar (ng, itlm, obs(ng)%name, &
843 & vname(1,idoval), obsval)
844 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
845
846 CALL netcdf_get_fvar (ng, itlm, obs(ng)%name, &
847 & vname(1,idoerr), obserr)
848 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
849
850# if defined PIO_LIB && defined DISTRIBUTE
851 CASE (io_pio)
852 CALL pio_netcdf_get_fvar (ng, itlm, obs(ng)%name, &
853 & vname(1,idoval), obsval)
854 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
855
856 CALL pio_netcdf_get_fvar (ng, itlm, obs(ng)%name, &
857 & vname(1,idoerr), obserr)
858 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
859# endif
860 END SELECT
861
862
863
864
865
866 DO ng=1,ngrids
867 WRITE (his(ng)%name,10) trim(fwd(ng)%head), my_outer-1
868 lstr=len_trim(his(ng)%name)
869 his(ng)%base=his(ng)%name(1:lstr-3)
870 IF (his(ng)%Nfiles.gt.1) THEN
871 DO ifile=1,his(ng)%Nfiles
872 WRITE (suffix,"('_',i4.4,'.nc')") ifile
873 his(ng)%files(ifile)=trim(his(ng)%base)//trim(suffix)
874 END DO
875 his(ng)%name=trim(his(ng)%files(1))
876 ELSE
877 his(ng)%files(1)=trim(his(ng)%name)
878 END IF
879 END DO
880
881
882
883
884
885
886
887 DO ng=1,ngrids
888 WRITE (qck(ng)%name,10) trim(qck(ng)%head), my_outer-1
889 lstr=len_trim(qck(ng)%name)
890 qck(ng)%base=qck(ng)%name(1:lstr-3)
891 IF (qck(ng)%Nfiles.gt.1) THEN
892 DO ifile=1,qck(ng)%Nfiles
893 WRITE (suffix,"('_',i4.4,'.nc')") ifile
894 qck(ng)%files(ifile)=trim(qck(ng)%base)//trim(suffix)
895 END DO
896 qck(ng)%name=trim(qck(ng)%files(1))
897 ELSE
898 qck(ng)%files(1)=trim(qck(ng)%name)
899 END IF
900 END DO
901
902
903
904
905
906
907 inprec=1
908 DO ng=1,ngrids
909 SELECT CASE (his(ng)%IOtype)
910 CASE (io_nf90)
911 CALL netcdf_get_fvar (ng, itlm, his(ng)%name, &
912 & vname(1,idtime), stime, &
913 & start = (/inprec/), &
914 & total = (/1/))
915
916# if defined PIO_LIB && defined DISTRIBUTE
917 CASE (io_pio)
918 CALL pio_netcdf_get_fvar (ng, itlm, his(ng)%name, &
919 & vname(1,idtime), stime, &
920 & start = (/inprec/), &
921 & total = (/1/))
922# endif
923 END SELECT
924 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
925 initime(ng)=stime
926
927# ifdef ADJUST_BOUNDARY
928
929
930
931 obc_time(1,ng)=stime
932 DO i=2,nbrec(ng)
933 obc_time(i,ng)=obc_time(i-1,ng)+nobc(ng)*dt(ng)
934 END DO
935# endif
936
937# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
938
939
940
941 sf_time(1,ng)=stime
942 DO i=2,nfrec(ng)
943 sf_time(i,ng)=sf_time(i-1,ng)+nsff(ng)*dt(ng)
944 END DO
945# endif
946 END DO
947# endif
948
949
950
951
952
953
954
955
956
957
958
960 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
961 DO ng=1,ngrids
962 lreadfwd(ng)=.true.
963 END DO
964
965# ifdef FORWARD_FLUXES
966
967
968
969
970
971
972
973
974
975
976
977
978
980 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
981 DO ng=1,ngrids
982 lreadblk(ng)=.true.
983 lreadfrc(ng)=.false.
984 lreadqck(ng)=.false.
985 END DO
986# endif
987
988# ifdef RPCG
989
990
991
992
993 DO ng=1,ngrids
994 IF (nadj(ng).lt.ntimes(ng)) THEN
995 nadrec(ng)=2*(1+ntimes(ng)/nadj(ng))
996 ELSE
997 nadrec(ng)=0
998 END IF
999 END DO
1000# endif
1001
1002
1003
1004
1005
1006
1007 check_outer1 : IF (my_outer.eq.1) THEN
1008
1009
1010
1011 DO ng=1,ngrids
1012 ldefitl(ng)=.true.
1013 CALL tl_def_ini (ng)
1014 ldefitl(ng)=.false.
1015 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1016 END DO
1017
1018# ifdef RPCG
1019
1020
1021
1022
1023 DO ng=1,ngrids
1024 tdays(ng)=initime(ng)*sec2day
1025# ifdef DISTRIBUTE
1026 CALL tl_wrt_ini (ng, myrank, rec1, rec1)
1027# else
1028 CALL tl_wrt_ini (ng, -1, rec1, rec1)
1029# endif
1030 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1031# ifdef DISTRIBUTE
1032 CALL tl_wrt_ini (ng, myrank, rec1, rec2)
1033# else
1034 CALL tl_wrt_ini (ng, -1, rec1, rec2)
1035# endif
1036 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1037# ifdef DISTRIBUTE
1038 CALL tl_wrt_ini (ng, myrank, rec1, rec3)
1039# else
1040 CALL tl_wrt_ini (ng, -1, rec1, rec3)
1041# endif
1042 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1043# ifdef DISTRIBUTE
1044 CALL tl_wrt_ini (ng, myrank, rec1, rec4)
1045# else
1046 CALL tl_wrt_ini (ng, -1, rec1, rec4)
1047# endif
1048 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1049# ifdef DISTRIBUTE
1050 CALL tl_wrt_ini (ng, myrank, rec1, rec5)
1051# else
1052 CALL tl_wrt_ini (ng, -1, rec1, rec5)
1053# endif
1054 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1055
1056 IF (nadj(ng).lt.ntimes(ng)) THEN
1057 nlast=rec5
1058 DO irec=1,nadrec(ng)
1059# ifdef DISTRIBUTE
1060 CALL tl_wrt_ini (ng, myrank, rec1, nlast+irec)
1061# else
1062 CALL tl_wrt_ini (ng, -1, rec1, nlast+irec)
1063# endif
1064 IF (founderror(exit_flag, noerror, &
1065 & __line__, myfile)) RETURN
1066 END DO
1067 END IF
1068 END DO
1069# endif
1070
1071
1072
1073 DO ng=1,ngrids
1074 ldeftlf(ng)=.true.
1075 CALL def_impulse (ng)
1076 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1077 END DO
1078
1079# if defined POSTERIOR_EOFS || defined POSTERIOR_ERROR_I || \
1080 defined posterior_error_f
1081
1082
1083
1084
1085
1086 DO ng=1,ngrids
1087 ldefhss(ng)=.true.
1088 CALL def_hessian (ng)
1089 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1090 END DO
1091# endif
1092
1093# if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F
1094
1095
1096
1097
1098 DO ng=1,ngrids
1099 ldeferr(ng)=.true.
1100 CALL def_error (ng)
1101 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1102 END DO
1103# endif
1104
1105 END IF check_outer1
1106
1107# ifdef RPCG
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118 check_outer2 : IF (my_outer.gt.1) THEN
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129 DO ng=1,ngrids
1130 DO tile=first_tile(ng),last_tile(ng),+1
1131 CALL initialize_ocean (ng, tile, inlm)
1132 CALL initialize_forces (ng, tile, itlm)
1133# ifdef ADJUST_BOUNDARY
1134 CALL initialize_boundary (ng, tile, itlm)
1135# endif
1136 END DO
1137 END DO
1138
1139
1140
1141
1142 DO ng=1,ngrids
1143 wrtnlmod(ng)=.false.
1144 wrttlmod(ng)=.true.
1145 END DO
1146
1147
1148
1149
1150 DO ng=1,ngrids
1151 IF (frcrec(ng).gt.3) THEN
1152 frequentimpulse(ng)=.true.
1153 END IF
1154 END DO
1155
1156
1157
1158
1159 DO ng=1,ngrids
1160 itl(ng)%Rindex=rec3
1162 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1163 END DO
1164
1165
1166
1167 DO ng=1,ngrids
1168 IF (master) THEN
1169 WRITE (stdout,20) 'TL', ng, my_outer, 0, &
1170 & ntstart(ng), ntend(ng)
1171 END IF
1172 END DO
1173
1174# ifdef SOLVE3D
1176# else
1178# endif
1179 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1180
1181 DO ng=1,ngrids
1182 wrtnlmod(ng)=.false.
1183 wrttlmod(ng)=.false.
1184 END DO
1185 END IF check_outer2
1186# endif
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196 DO ng=1,ngrids
1197 DO tile=first_tile(ng),last_tile(ng),+1
1198 CALL initialize_ocean (ng, tile, inlm)
1199 CALL initialize_forces (ng, tile, itlm)
1200# ifdef ADJUST_BOUNDARY
1201 CALL initialize_boundary (ng, tile, itlm)
1202# endif
1203 END DO
1204 END DO
1205
1206# if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
1207
1208
1209
1210
1211 IF (balance(isfsur)) THEN
1212 DO ng=1,ngrids
1213 CALL get_state (ng, inlm, 2, ini(ng), lini, lini)
1214 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1215
1216 DO tile=first_tile(ng),last_tile(ng),+1
1217 CALL balance_ref (ng, tile, lini)
1218 CALL biconj (ng, tile, inlm, lini)
1219 END DO
1220 wrtzetaref(ng)=.true.
1221 END DO
1222 END IF
1223# endif
1224
1225
1226
1227 inner_loop : DO my_inner=0,ninner
1228 inner=my_inner
1229
1230# ifdef RPCG
1231 IF (inner.ne.ninner) THEN
1232 linner=.true.
1233 ELSE
1234 linner=.false.
1235 END IF
1236
1237
1238
1239
1240 IF ((inner.eq.0).and.(my_outer.gt.1)) THEN
1241 DO ng=1,ngrids
1242 SELECT CASE (dav(ng)%IOtype)
1243 CASE (io_nf90)
1244 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
1245 & 'TLmodel_value', tlmodval)
1246 IF (founderror(exit_flag, noerror, __line__, &
1247 & myfile)) RETURN
1248
1249 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
1250 & 'NLmodel_value', nlmodval, &
1251 & ncid = dav(ng)%ncid, &
1252 & start = (/1,my_outer-1/), &
1253 & total = (/ndatum(ng),1/))
1254 IF (founderror(exit_flag, noerror, __line__, &
1255 & myfile)) RETURN
1256
1257# if defined PIO_LIB && defined DISTRIBUTE
1258 CASE (io_pio)
1259 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
1260 & 'TLmodel_value', tlmodval)
1261 IF (founderror(exit_flag, noerror, __line__, &
1262 & myfile)) RETURN
1263
1264 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
1265 & 'NLmodel_value', nlmodval, &
1266 & piofile = dav(ng)%pioFile, &
1267 & start = (/1,my_outer-1/), &
1268 & total = (/ndatum(ng),1/))
1269 IF (founderror(exit_flag, noerror, __line__, &
1270 & myfile)) RETURN
1271# endif
1272 END SELECT
1273 END DO
1274 END IF
1275
1276 IF (inner.eq.0) THEN
1277 lcgini=.true.
1278 END IF
1279 DO ng=1,ngrids
1280 CALL rpcg_lanczos (ng, itlm, my_outer, inner, ninner, lcgini)
1281 END DO
1282# else
1283
1284
1285
1286
1287 IF (inner.eq.0) THEN
1288 lcgini=.true.
1289 DO ng=1,ngrids
1290 CALL congrad (ng, irpm, my_outer, inner, ninner, lcgini)
1291 END DO
1292 END IF
1293
1294
1295
1296 linner=.false.
1297 IF ((inner.ne.0).or.(nrun.ne.1)) THEN
1298 IF (((inner.eq.0).and.lhotstart).or.(inner.ne.0)) THEN
1299 linner=.true.
1300 END IF
1301 END IF
1302# endif
1303
1304
1305
1306 inner_compute : IF (linner) THEN
1307
1308
1309
1310
1311
1312
1313
1314
1315 DO ng=1,ngrids
1317 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1318 wrtmisfit(ng)=.false.
1319 END DO
1320
1321
1322
1323
1324 DO ng=1,ngrids
1325# ifdef WEAK_NOINTERP
1326 wrtforce(ng)=.false.
1327# else
1328 wrtforce(ng)=.true.
1329# endif
1330 IF (nrun.gt.1) ldefadj(ng)=.false.
1331 fcount=adm(ng)%load
1332 adm(ng)%Nrec(fcount)=0
1333 adm(ng)%Rindex=0
1334 END DO
1335
1336
1337
1338 DO ng=1,ngrids
1339 IF (master) THEN
1340 WRITE (stdout,20) 'AD', ng, my_outer, my_inner, &
1341 & ntstart(ng), ntend(ng)
1342 END IF
1343 END DO
1344
1345# ifdef SOLVE3D
1347# else
1349# endif
1350 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361 DO ng=1,ngrids
1362 lwrtstate2d(ng)=.true.
1363 END DO
1364
1365# ifndef WEAK_NOINTERP
1366
1367
1368
1369
1370
1371 DO ng=1,ngrids
1372# ifdef DISTRIBUTE
1373 CALL ad_wrt_his (ng, myrank)
1374# else
1375 CALL ad_wrt_his (ng, -1)
1376# endif
1377 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1378 END DO
1379# endif
1380
1381
1382
1383
1384 DO ng=1,ngrids
1385 wrtforce(ng)=.false.
1386# ifdef DISTRIBUTE
1387 CALL ad_wrt_his (ng, myrank)
1388# else
1389 CALL ad_wrt_his (ng, -1)
1390# endif
1391 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1392 END DO
1393
1394
1395
1396# ifdef POSTERIOR_ERROR_I
1397 lposterior=.true.
1398# else
1399 lposterior=.false.
1400# endif
1401# ifdef RPCG
1402
1403
1404
1405
1406
1407 laugweak=.false.
1408# endif
1409 CALL error_covariance (itlm, driver, my_outer, inner, &
1410 & lbck, lini, lold, lnew, &
1411 & rec1, rec2, lposterior)
1412 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1413
1414
1415
1416
1417
1418
1419
1420 IF (master) THEN
1421 WRITE (stdout,30) my_outer, inner
1422 END IF
1423 DO ng=1,ngrids
1424 tlf(ng)%Rindex=0
1425# ifdef DISTRIBUTE
1426 CALL wrt_impulse (ng, myrank, iadm, adm(ng)%name)
1427# else
1428 CALL wrt_impulse (ng, -1, iadm, adm(ng)%name)
1429# endif
1430 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1431 END DO
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442 DO ng=1,ngrids
1443 wrtnlmod(ng)=.false.
1444 wrttlmod(ng)=.true.
1445 END DO
1446
1447
1448
1449
1450 DO ng=1,ngrids
1451 IF (frcrec(ng).gt.3) THEN
1452 frequentimpulse(ng)=.true.
1453 END IF
1454 END DO
1455
1456
1457
1458 DO ng=1,ngrids
1459 itl(ng)%Rindex=rec1
1461 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1462 END DO
1463
1464
1465
1466
1467 IF ((my_outer.eq.1).and.(inner.eq.1)) THEN
1468 DO ng=1,ngrids
1469 wrtmisfit(ng)=.true.
1470 END DO
1471 END IF
1472
1473
1474
1475
1476
1477 DO ng=1,ngrids
1478 IF (master) THEN
1479 WRITE (stdout,20) 'TL', ng, my_outer, my_inner, &
1480 & ntstart(ng), ntend(ng)
1481 END IF
1482 END DO
1483
1484# ifdef SOLVE3D
1486# else
1488# endif
1489 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1490
1491 DO ng=1,ngrids
1492 wrtnlmod(ng)=.false.
1493 wrttlmod(ng)=.false.
1494 END DO
1495
1496# ifdef POSTERIOR_ERROR_F
1497
1498
1499
1500
1501 add=.false.
1502 DO ng=1,ngrids
1503 DO tile=first_tile(ng),last_tile(ng),+1
1504 CALL load_tltoad (ng, tile, lold(ng), lold(ng), add)
1505 END DO
1506 END DO
1507
1508
1509
1510
1511 IF (inner.ne.0) THEN
1512 DO ng=1,ngrids
1513# ifdef DISTRIBUTE
1514 CALL wrt_hessian (ng, tile, lold(ng), lold(ng))
1515# else
1516 CALL wrt_hessian (ng, -1, lold(ng), lold(ng))
1517# endif
1518 IF (founderror(exit_flag, noerror, &
1519 & __line__, myfile)) RETURN
1520 END DO
1521 END IF
1522# endif
1523
1524 nrun=nrun+1
1525
1526# ifndef RPCG
1527
1528
1529
1530
1531 DO ng=1,ngrids
1532 lcgini=.false.
1533 CALL congrad (ng, itlm, my_outer, inner, ninner, lcgini)
1534 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1535 END DO
1536# else
1537 lcgini=.false.
1538# endif
1539
1540 END IF inner_compute
1541
1542 END DO inner_loop
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553 DO ng=1,ngrids
1555 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1556 END DO
1557
1558
1559
1560
1561 DO ng=1,ngrids
1562# ifdef WEAK_NOINTERP
1563 wrtforce(ng)=.false.
1564# else
1565 wrtforce(ng)=.true.
1566# endif
1567 IF (nrun.gt.1) ldefadj(ng)=.false.
1568 fcount=adm(ng)%load
1569 adm(ng)%Nrec(fcount)=0
1570 adm(ng)%Rindex=0
1571 END DO
1572
1573
1574
1575
1576 DO ng=1,ngrids
1577 IF (master) THEN
1578 WRITE (stdout,20) 'AD', ng, my_outer, ninner, &
1579 & ntstart(ng), ntend(ng)
1580 END IF
1581 END DO
1582
1583# ifdef SOLVE3D
1585# else
1587# endif
1588 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1589
1590# ifndef WEAK_NOINTERP
1591
1592
1593
1594
1595
1596 DO ng=1,ngrids
1597# ifdef DISTRIBUTE
1598 CALL ad_wrt_his (ng, myrank)
1599# else
1600 CALL ad_wrt_his (ng, -1)
1601# endif
1602 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1603 END DO
1604# endif
1605
1606
1607
1608
1609 DO ng=1,ngrids
1610 wrtforce(ng)=.false.
1611# ifdef DISTRIBUTE
1612 CALL ad_wrt_his (ng, myrank)
1613# else
1614 CALL ad_wrt_his (ng, -1)
1615# endif
1616 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1617 END DO
1618
1619# ifdef RPCG
1620
1621
1622
1623
1624
1625 DO ng=1,ngrids
1626 fcount=adm(ng)%load
1627 nrec(ng)=adm(ng)%Nrec(fcount)
1628 END DO
1629# endif
1630
1631
1632
1633 lposterior=.false.
1634
1635# ifdef RPCG
1636
1637
1638
1639
1640
1641 laugweak=.true.
1642 linl=my_outer+1
1643 CALL error_covariance (inlm, driver, my_outer, inner, &
1644 & linl, lini, lold, lnew, &
1645 & rec1, rec2, lposterior)
1646 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1647# else
1648 CALL error_covariance (inlm, driver, my_outer, inner, &
1649 & lbck, lini, lold, lnew, &
1650 & rec1, rec2, lposterior)
1651 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1652# endif
1653
1654# ifdef RPCG
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668 laugweak=.false.
1669
1670# else
1671
1672
1673
1674
1675 DO ng=1,ngrids
1676# ifdef DISTRIBUTE
1677 CALL wrt_ini (ng, myrank, lnew(ng))
1678# else
1679 CALL wrt_ini (ng, -1, lnew(ng))
1680# endif
1681 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1682
1683# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
1684 defined adjust_boundary
1685# ifdef DISTRIBUTE
1686 CALL wrt_frc_ad (ng, myrank, lold(ng), ini(ng)%Rindex)
1687# else
1688 CALL wrt_frc_ad (ng, -1, lold(ng), ini(ng)%Rindex)
1689# endif
1690 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1691# endif
1692 END DO
1693# endif
1694
1695# ifdef RPCG
1696
1697
1698
1699
1700
1701 DO ng=1,ngrids
1702 CALL get_state (ng, itlm, 4, itl(ng), rec3, ltlm1)
1703 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1704 END DO
1705 DO ng=1,ngrids
1706 DO tile=first_tile(ng),last_tile(ng),+1
1707 CALL aug_oper (ng, tile, ltlm1, ltlm2)
1708 END DO
1709 END DO
1710
1711
1712
1713 DO ng=1,ngrids
1714# ifdef DISTRIBUTE
1715 CALL tl_wrt_ini (ng, myrank, ltlm2, rec5)
1716# else
1717 CALL tl_wrt_ini (ng, -1, ltlm2, rec5)
1718# endif
1719 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1720 END DO
1721
1722
1723
1724
1725 DO ng=1,ngrids
1726 CALL get_state (ng, itlm, 4, itl(ng), rec2, ltlm1)
1727 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1728 CALL get_state (ng, itlm, 4, itl(ng), rec5, ltlm2)
1729 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1730 END DO
1731
1732 DO ng=1,ngrids
1733 DO tile=first_tile(ng),last_tile(ng),+1
1734 CALL sum_grad (ng, tile, ltlm1, ltlm2)
1735 END DO
1736 END DO
1737
1738
1739
1740 DO ng=1,ngrids
1741# ifdef DISTRIBUTE
1742 CALL tl_wrt_ini (ng, myrank, ltlm2, rec2)
1743# else
1744 CALL tl_wrt_ini (ng, -1, ltlm2, rec2)
1745# endif
1746 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1747 END DO
1748
1749
1750
1751 linl=my_outer+1
1752 DO ng=1,ngrids
1753 CALL get_state (ng, inlm, 9, ini(ng), linl, lnew(ng))
1754 IF (timeiau(ng).eq.0.0_dp) THEN
1755 CALL get_state (ng, iadm, 4, itl(ng), rec2, ltlm1)
1756 END IF
1757 END DO
1758 DO ng=1,ngrids
1759 IF (timeiau(ng).eq.0.0_dp) THEN
1760 DO tile=first_tile(ng),last_tile(ng),+1
1762 END DO
1763 END IF
1764 END DO
1765
1766
1767
1768
1769 DO ng=1,ngrids
1770# ifdef DISTRIBUTE
1771 CALL wrt_ini (ng, myrank, lnew(ng))
1772# else
1773 CALL wrt_ini (ng, -1, lnew(ng))
1774# endif
1775 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1776
1777# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
1778 defined adjust_boundary
1779# ifdef DISTRIBUTE
1780 CALL wrt_frc_ad (ng, myrank, ltlm1, ini(ng)%Rindex)
1781# else
1782 CALL wrt_frc_ad (ng, -1, ltlm1, ini(ng)%Rindex)
1783# endif
1784 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1785# endif
1786 END DO
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796 DO ng=1,ngrids
1797 CALL get_state (ng, itlm, 4, itl(ng), rec4, ltlm1)
1798 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1799 END DO
1800
1801 DO ng=1,ngrids
1802 DO tile=first_tile(ng),last_tile(ng),+1
1803 CALL aug_oper (ng, tile, ltlm1, ltlm2)
1804 END DO
1805 END DO
1806
1807 DO ng=1,ngrids
1808 DO tile=first_tile(ng),last_tile(ng),+1
1809 CALL sum_grad (ng, tile, ltlm1, ltlm2)
1810 END DO
1811 END DO
1812
1813 DO ng=1,ngrids
1814 adrec=nrec(ng)
1815 CALL get_state (ng, itlm, 4, adm(ng), adrec, ltlm1)
1816 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1817 END DO
1818
1819 DO ng=1,ngrids
1820 DO tile=first_tile(ng),last_tile(ng),+1
1821 CALL sum_grad (ng, tile, ltlm1, ltlm2)
1822 END DO
1823 END DO
1824
1825
1826
1827
1828 DO ng=1,ngrids
1829# ifdef DISTRIBUTE
1830 CALL tl_wrt_ini (ng, myrank, ltlm2, rec4)
1831# else
1832 CALL tl_wrt_ini (ng, -1, ltlm2, rec4)
1833# endif
1834 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1835 END DO
1836# endif
1837
1838
1839
1840
1841
1842
1843
1844 IF (master) THEN
1845 WRITE (stdout,30) my_outer, inner
1846 END IF
1847 DO ng=1,ngrids
1848 tlf(ng)%Rindex=0
1849# ifdef DISTRIBUTE
1850 CALL wrt_impulse (ng, myrank, iadm, adm(ng)%name)
1851# else
1852 CALL wrt_impulse (ng, -1, iadm, adm(ng)%name)
1853# endif
1854 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1855 END DO
1856
1857# ifdef RPCG
1858
1859
1860
1861
1862
1863 DO ng=1,ngrids
1864 DO i=1,nadrec(ng)/2
1865 irec=i
1866 jrec=rec5+i
1867 CALL get_state (ng, itlm, 4, itl(ng), jrec, ltlm1)
1868 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1869
1870 DO tile=first_tile(ng),last_tile(ng),+1
1871 CALL aug_oper (ng, tile, ltlm1, ltlm2)
1872 END DO
1873
1874
1875
1876
1877
1878
1879
1880 CALL get_state (ng, 7, 4, tlf(ng), irec, ltlm1)
1881 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1882
1883 DO tile=first_tile(ng),last_tile(ng),+1
1884
1885 CALL sum_imp (ng, tile, ltlm2)
1886 END DO
1887
1888
1889
1890
1891 tlf(ng)%Rindex=0
1892# ifdef DISTRIBUTE
1893 CALL wrt_aug_imp (ng, myrank, itlm, ltlm2, i)
1894# else
1895 CALL wrt_aug_imp (ng, -1, itlm, ltlm2, i)
1896# endif
1897
1898 END DO
1899 END DO
1900
1901
1902
1903
1904
1905 DO ng=1,ngrids
1906 CALL get_state (ng, itlm, 4, itl(ng), rec2, ltlm1)
1907 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1908 CALL get_state (ng, itlm, 4, itl(ng), rec3, ltlm2)
1909 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1910 END DO
1911
1912 DO ng=1,ngrids
1913 DO tile=first_tile(ng),last_tile(ng),+1
1914 CALL sum_grad (ng, tile, ltlm1, ltlm2)
1915 END DO
1916 END DO
1917
1918
1919
1920 DO ng=1,ngrids
1921# ifdef DISTRIBUTE
1922 CALL tl_wrt_ini (ng, myrank, ltlm2, rec3)
1923# else
1924 CALL tl_wrt_ini (ng, -1, ltlm2, rec3)
1925# endif
1926 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1927 END DO
1928
1929# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
1930 defined adjust_boundary
1931
1932
1933
1934
1935
1936
1937 DO ng=1,ngrids
1938 CALL get_state (ng, iadm, 4, itl(ng), rec3, ltlm1)
1939 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1940 END DO
1941 DO ng=1,ngrids
1942# ifdef DISTRIBUTE
1943 CALL wrt_frc_ad (ng, myrank, ltlm1, ini(ng)%Rindex)
1944# else
1945 CALL wrt_frc_ad (ng, -1, ltlm1, ini(ng)%Rindex)
1946# endif
1947 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1948 END DO
1949# endif
1950
1951
1952
1953
1954 DO ng=1,ngrids
1955 DO i=1,nadrec(ng)/2
1956 irec=i
1957 jrec=rec5+i
1958 CALL get_state (ng, itlm, 4, itl(ng), jrec, ltlm1)
1959 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1960
1961 CALL get_state (ng, 7, 4, tlf(ng), irec, ltlm2)
1962 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1963
1964 DO tile=first_tile(ng),last_tile(ng),+1
1965
1966 CALL sum_imp (ng, tile, ltlm1)
1967 END DO
1968
1969
1970
1971# ifdef DISTRIBUTE
1972
1973 CALL tl_wrt_ini (ng, myrank, ltlm1, jrec)
1974# else
1975
1976 CALL tl_wrt_ini (ng, -1, ltlm1, jrec)
1977# endif
1978 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1979 END DO
1980 END DO
1981
1982
1983
1984
1985
1986 jb0(my_outer)=0.0_r8
1987 DO ng=1,ngrids
1988 CALL get_state (ng, iadm, 8, itl(ng), rec4, ltlm1)
1989 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1990 CALL get_state (ng, itlm, 8, itl(ng), rec3, ltlm1)
1991 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1992 END DO
1993
1994 DO ng=1,ngrids
1995 DO tile=first_tile(ng),last_tile(ng),+1
1996 CALL comp_jb0 (ng, tile, itlm, my_outer, ltlm1, ltlm1)
1997 END DO
1998 END DO
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011 DO ng=1,ngrids
2012 DO irec=1,nadrec(ng)/2
2013 jrec1=rec5+irec
2014 jrec2=rec5+nadrec(ng)/2+irec
2015 CALL get_state (ng, iadm, 8, itl(ng), jrec1, ltlm1)
2016 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2017 CALL get_state (ng, itlm, 8, itl(ng), jrec2, ltlm1)
2018 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2019
2020 DO tile=first_tile(ng),last_tile(ng),+1
2021 CALL comp_jb0 (ng, tile, itlm, my_outer, ltlm1, ltlm1)
2022 END DO
2023 END DO
2024 END DO
2025
2026
2027
2028
2029
2030 DO ng=1,ngrids
2031 DO irec=1,nadrec(ng)/2
2032 jrec=rec5+irec
2033 CALL get_state (ng, itlm, 8, itl(ng), jrec, ltlm1)
2034
2035 tlf(ng)%Rindex=0
2036# ifdef DISTRIBUTE
2037 CALL wrt_aug_imp (ng, myrank, itlm, ltlm1, irec)
2038# else
2039 CALL wrt_aug_imp (ng, -1, itlm, ltlm1, irec)
2040# endif
2041 END DO
2042 END DO
2043# endif
2044
2045# ifdef PROFILE
2046
2047
2048
2049 DO ng=1,ngrids
2050 DO thread=thread_range
2051 CALL wclock_off (ng, itlm, 87, __line__, myfile)
2052 END DO
2053 END DO
2054# endif
2055
2056 10 FORMAT (a,'_outer',i0,'.nc')
2057 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
2058 & ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3, &
2059 ', TimeSteps: ',i0,' - ',i0,')',/)
2060 30 FORMAT (/,' Converting Convolved Adjoint Trajectory to', &
2061 & ' Impulses: Outer = ',i3.3,' Inner = ',i3.3,/)
2062
2063 RETURN
subroutine ad_initial(ng)
subroutine ad_main3d(runinterval)
subroutine, public comp_jb0(ng, tile, model, outloop, ltlm, ladj)
subroutine, public aug_oper(ng, tile, linp, lout)
subroutine, public load_tltoad(ng, tile, linp, lout, add)
subroutine, public ini_adjust(ng, tile, linp, lout)
subroutine, public sum_grad(ng, tile, linp, lout)
subroutine, public sum_imp(ng, tile, lout)
subroutine tl_initial(ng)
subroutine tl_main3d(runinterval)