libsim Versione 7.2.1
|
◆ pack_distinct_sorted_ttr_mapper()
compatta gli elementi distinti di vect in un sorted array Definizione alla linea 932 del file stat_proc_engine.F90. 934! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
935! authors:
936! Davide Cesari <dcesari@arpa.emr.it>
937! Paolo Patruno <ppatruno@arpa.emr.it>
938
939! This program is free software; you can redistribute it and/or
940! modify it under the terms of the GNU General Public License as
941! published by the Free Software Foundation; either version 2 of
942! the License, or (at your option) any later version.
943
944! This program is distributed in the hope that it will be useful,
945! but WITHOUT ANY WARRANTY; without even the implied warranty of
946! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
947! GNU General Public License for more details.
948
949! You should have received a copy of the GNU General Public License
950! along with this program. If not, see <http://www.gnu.org/licenses/>.
951#include "config.h"
952
959IMPLICIT NONE
960
961! per ogni elemento i,j dell'output, definire n elementi di input ad
962! esso contribuenti (arrayof_ttr_mapper) con le seguenti informazioni
963TYPE ttr_mapper
964 INTEGER :: it=imiss ! k
965 INTEGER :: itr=imiss ! l
966 INTEGER :: extra_info=imiss ! dtratio for intervals, extreme for point in time
967 TYPE(datetime) :: time=datetime_miss ! time for point in time
968END TYPE ttr_mapper
969
970INTERFACE OPERATOR (==)
971 MODULE PROCEDURE ttr_mapper_eq
972END INTERFACE
973
974INTERFACE OPERATOR (/=)
975 MODULE PROCEDURE ttr_mapper_ne
976END INTERFACE
977
978INTERFACE OPERATOR (>)
979 MODULE PROCEDURE ttr_mapper_gt
980END INTERFACE
981
982INTERFACE OPERATOR (<)
983 MODULE PROCEDURE ttr_mapper_lt
984END INTERFACE
985
986INTERFACE OPERATOR (>=)
987 MODULE PROCEDURE ttr_mapper_ge
988END INTERFACE
989
990INTERFACE OPERATOR (<=)
991 MODULE PROCEDURE ttr_mapper_le
992END INTERFACE
993
994#undef VOL7D_POLY_TYPE
995#undef VOL7D_POLY_TYPES
996#undef ENABLE_SORT
997#define VOL7D_POLY_TYPE TYPE(ttr_mapper)
998#define VOL7D_POLY_TYPES _ttr_mapper
999#define ENABLE_SORT
1000#include "array_utilities_pre.F90"
1001
1002#define ARRAYOF_ORIGTYPE TYPE(ttr_mapper)
1003#define ARRAYOF_TYPE arrayof_ttr_mapper
1004#define ARRAYOF_ORIGEQ 1
1005#define ARRAYOF_ORIGGT 1
1006#include "arrayof_pre.F90"
1007! from arrayof
1008
1009
1010CONTAINS
1011
1012! simplified comparisons only on time
1013ELEMENTAL FUNCTION ttr_mapper_eq(this, that) RESULT(res)
1014TYPE(ttr_mapper),INTENT(IN) :: this, that
1015LOGICAL :: res
1016
1017res = this%time == that%time
1018
1019END FUNCTION ttr_mapper_eq
1020
1021ELEMENTAL FUNCTION ttr_mapper_ne(this, that) RESULT(res)
1022TYPE(ttr_mapper),INTENT(IN) :: this, that
1023LOGICAL :: res
1024
1025res = this%time /= that%time
1026
1027END FUNCTION ttr_mapper_ne
1028
1029ELEMENTAL FUNCTION ttr_mapper_gt(this, that) RESULT(res)
1030TYPE(ttr_mapper),INTENT(IN) :: this, that
1031LOGICAL :: res
1032
1033res = this%time > that%time
1034
1035END FUNCTION ttr_mapper_gt
1036
1037ELEMENTAL FUNCTION ttr_mapper_lt(this, that) RESULT(res)
1038TYPE(ttr_mapper),INTENT(IN) :: this, that
1039LOGICAL :: res
1040
1041res = this%time < that%time
1042
1043END FUNCTION ttr_mapper_lt
1044
1045ELEMENTAL FUNCTION ttr_mapper_ge(this, that) RESULT(res)
1046TYPE(ttr_mapper),INTENT(IN) :: this, that
1047LOGICAL :: res
1048
1049res = this%time >= that%time
1050
1051END FUNCTION ttr_mapper_ge
1052
1053ELEMENTAL FUNCTION ttr_mapper_le(this, that) RESULT(res)
1054TYPE(ttr_mapper),INTENT(IN) :: this, that
1055LOGICAL :: res
1056
1057res = this%time <= that%time
1058
1059END FUNCTION ttr_mapper_le
1060
1061#include "arrayof_post.F90"
1062#include "array_utilities_inc.F90"
1063
1064
1065! common operations for statistical processing by differences
1066SUBROUTINE recompute_stat_proc_diff_common(itime, itimerange, stat_proc, step, &
1067 otime, otimerange, map_tr, f, keep_tr, time_definition, full_steps, &
1068 start)
1069TYPE(datetime),INTENT(in) :: itime(:)
1070TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1071INTEGER,INTENT(in) :: stat_proc
1072TYPE(timedelta),INTENT(in) :: step
1073TYPE(datetime),POINTER :: otime(:)
1074TYPE(vol7d_timerange),POINTER :: otimerange(:)
1075INTEGER,ALLOCATABLE,INTENT(out) :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
1076INTEGER :: nitr
1077LOGICAL,ALLOCATABLE :: mask_timerange(:)
1078INTEGER,INTENT(in) :: time_definition
1079LOGICAL,INTENT(in),OPTIONAL :: full_steps
1080TYPE(datetime),INTENT(in),OPTIONAL :: start
1081
1082INTEGER :: i, j, k, l, dirtyrep
1083INTEGER :: steps
1084LOGICAL :: lfull_steps, useful
1085TYPE(datetime) :: lstart, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1086TYPE(vol7d_timerange) :: tmptimerange
1087TYPE(arrayof_datetime) :: a_otime
1088TYPE(arrayof_vol7d_timerange) :: a_otimerange
1089
1090! compute length of cumulation step in seconds
1092
1093lstart = datetime_miss
1094IF (PRESENT(start)) lstart = start
1095lfull_steps = optio_log(full_steps)
1096
1097! create a mask of suitable timeranges
1098ALLOCATE(mask_timerange(SIZE(itimerange)))
1099mask_timerange(:) = itimerange(:)%timerange == stat_proc &
1100 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
1101 .AND. itimerange(:)%p1 >= 0 &
1102 .AND. itimerange(:)%p2 > 0
1103
1104IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges
1105 mask_timerange(:) = mask_timerange(:) .AND. &
1106 (itimerange(:)%p1 == 0 .OR. & ! all analyses
1109ENDIF
1110! mask_timerange includes all candidate timeranges
1111
1112nitr = count(mask_timerange)
1113ALLOCATE(f(nitr))
1114j = 1
1115DO i = 1, nitr
1116 DO WHILE(.NOT.mask_timerange(j))
1117 j = j + 1
1118 ENDDO
1119 f(i) = j
1120 j = j + 1
1121ENDDO
1122
1123! now we have to evaluate time/timerage pairs which do not need processing
1124ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
1125CALL compute_keep_tr()
1126
1127ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
1128map_tr(:,:,:,:,:) = imiss
1129
1130! scan through all possible combinations of time and timerange
1131DO dirtyrep = 1, 2
1132 IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
1137 ENDIF
1138 DO l = 1, SIZE(itime)
1139 DO k = 1, nitr
1140 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1141 time_definition, pstart2, pend2, reftime2)
1142
1143 DO j = 1, SIZE(itime)
1144 DO i = 1, nitr
1145 useful = .false.
1146 CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
1147 time_definition, pstart1, pend1, reftime1)
1148 tmptimerange = vol7d_timerange_new(timerange=stat_proc)
1149
1150 IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
1151 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
1152 CALL time_timerange_set_period(tmptime, tmptimerange, &
1153 time_definition, pend1, pend2, reftime2)
1154 IF (lfull_steps) THEN
1156 useful = .true.
1157 ENDIF
1158 ELSE
1159 useful = .true.
1160 ENDIF
1161
1162 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
1163 CALL time_timerange_set_period(tmptime, tmptimerange, &
1164 time_definition, pstart2, pstart1, pstart1)
1165 IF (lfull_steps) THEN
1167 useful = .true.
1168 ENDIF
1169 ELSE
1170 useful = .true.
1171 ENDIF
1172 ENDIF
1173
1174 ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
1175 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
1176 CALL time_timerange_set_period(tmptime, tmptimerange, &
1177 time_definition, pend1, pend2, reftime2)
1178 IF (lfull_steps) THEN
1180 useful = .true.
1181 ENDIF
1182 ELSE
1183 useful = .true.
1184 ENDIF
1185
1186 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! |-=
1187 CALL time_timerange_set_period(tmptime, tmptimerange, &
1188 time_definition, pstart2, pstart1, reftime2)
1189 IF (lfull_steps) THEN
1191 useful = .true.
1192 ENDIF
1193 ELSE
1194 useful = .true.
1195 ENDIF
1196
1197 ENDIF
1198 ENDIF
1199 useful = useful .AND. tmptime /= datetime_miss .AND. &
1200 tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
1201
1202 IF (useful) THEN ! add a_otime, a_otimerange
1203 map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
1204 map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
1205 ENDIF
1206 ENDDO
1207 ENDDO
1208 ENDDO
1209 ENDDO
1210ENDDO
1211
1212! we have to repeat the computation with sorted arrays
1213CALL compute_keep_tr()
1214
1215otime => a_otime%array
1216otimerange => a_otimerange%array
1217! delete local objects keeping the contents
1220
1221#ifdef DEBUG
1222CALL l4f_log(l4f_debug, &
1227CALL l4f_log(l4f_debug, &
1230CALL l4f_log(l4f_debug, &
1232CALL l4f_log(l4f_debug, &
1234CALL l4f_log(l4f_debug, &
1236CALL l4f_log(l4f_debug, &
1238#endif
1239
1240CONTAINS
1241
1242SUBROUTINE compute_keep_tr()
1243
1244keep_tr(:,:,:) = imiss
1245DO l = 1, SIZE(itime)
1246 DO k = 1, nitr
1247 IF (itimerange(f(k))%p2 == steps) THEN
1248 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1249 time_definition, pstart2, pend2, reftime2)
1250 useful = .false.
1251 IF (reftime2 == pend2) THEN ! analysis
1254 useful = .true.
1255 ENDIF
1256 ELSE IF (lfull_steps) THEN
1258 useful = .true.
1259 ENDIF
1260 ELSE
1261 useful = .true.
1262 ENDIF
1263 ELSE ! forecast
1264 IF (lfull_steps) THEN
1266 useful = .true.
1267 ENDIF
1268 ELSE
1269 useful = .true.
1270 ENDIF
1271 ENDIF
1272 IF (useful) THEN
1273! CALL time_timerange_set_period(tmptime, tmptimerange, &
1274! time_definition, pstart2, pend2, reftime2)
1275 keep_tr(k,l,1) = append_unique(a_otime, itime(l))
1276 keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k)))
1277 ENDIF
1278 ENDIF
1279 ENDDO
1280ENDDO
1281
1282END SUBROUTINE compute_keep_tr
1283
1284END SUBROUTINE recompute_stat_proc_diff_common
1285
1286
1287! common operations for statistical processing by metamorphosis
1288SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
1289 otimerange, map_tr)
1290INTEGER,INTENT(in) :: istat_proc
1291TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1292INTEGER,INTENT(in) :: ostat_proc
1293TYPE(vol7d_timerange),POINTER :: otimerange(:)
1294INTEGER,POINTER :: map_tr(:)
1295
1296INTEGER :: i
1297LOGICAL :: tr_mask(SIZE(itimerange))
1298
1299IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1300 ALLOCATE(otimerange(0), map_tr(0))
1301 RETURN
1302ENDIF
1303
1304! useful timeranges
1305tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
1306 .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
1307ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
1308
1309otimerange = pack(itimerange, mask=tr_mask)
1310otimerange(:)%timerange = ostat_proc
1311map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
1312
1313END SUBROUTINE compute_stat_proc_metamorph_common
1314
1315
1316! common operations for statistical processing by aggregation
1317SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
1318 step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
1319TYPE(datetime),INTENT(in) :: itime(:)
1320TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1321INTEGER,INTENT(in) :: stat_proc
1322INTEGER,INTENT(in) :: tri
1323TYPE(timedelta),INTENT(in) :: step
1324INTEGER,INTENT(in) :: time_definition
1325TYPE(datetime),POINTER :: otime(:)
1326TYPE(vol7d_timerange),POINTER :: otimerange(:)
1327TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
1328INTEGER,POINTER,OPTIONAL :: dtratio(:)
1329TYPE(datetime),INTENT(in),OPTIONAL :: start
1330LOGICAL,INTENT(in),OPTIONAL :: full_steps
1331
1332INTEGER :: i, j, k, l, na, nf, n
1333INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
1334INTEGER(kind=int_ll) :: stepms, mstepms
1335LOGICAL :: lforecast
1336TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1337TYPE(arrayof_datetime) :: a_otime
1338TYPE(arrayof_vol7d_timerange) :: a_otimerange
1339TYPE(arrayof_integer) :: a_dtratio
1340LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
1341TYPE(ttr_mapper) :: lmapper
1342CHARACTER(len=8) :: env_var
1343LOGICAL :: climat_behavior
1344
1345
1346! enable bad behavior for climat database (only for instantaneous data)
1347env_var = ''
1348CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
1349climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
1350
1351! compute length of cumulation step in seconds
1353
1354! create a mask of suitable timeranges
1355ALLOCATE(mask_timerange(SIZE(itimerange)))
1356mask_timerange(:) = itimerange(:)%timerange == tri &
1357 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
1358
1359IF (PRESENT(dtratio)) THEN
1360 WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0)
1361 mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0
1362 ELSEWHERE
1363 mask_timerange(:) = .false.
1364 END WHERE
1365ELSE ! instantaneous
1366 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
1367ENDIF
1368
1369#ifdef DEBUG
1370CALL l4f_log(l4f_debug, &
1371 '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
1372 t2c(count(mask_timerange)))
1373#endif
1374
1375! euristically determine whether we are dealing with an
1376! analysis/observation or a forecast dataset
1377na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
1378nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
1379lforecast = nf >= na
1380#ifdef DEBUG
1381CALL l4f_log(l4f_debug, &
1383#endif
1384
1385! keep only timeranges of one type (really necessary?)
1386IF (lforecast) THEN
1387 CALL l4f_log(l4f_info, &
1388 'recompute_stat_proc_agg, processing in forecast mode')
1389ELSE
1390 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
1391 CALL l4f_log(l4f_info, &
1392 'recompute_stat_proc_agg, processing in analysis mode')
1393ENDIF
1394
1395#ifdef DEBUG
1396CALL l4f_log(l4f_debug, &
1397 '(re)compute_stat_proc_agg, number of useful timeranges: '// &
1398 t2c(count(mask_timerange)))
1399#endif
1400
1401IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1402 ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
1403 IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
1404 RETURN
1405ENDIF
1406
1407! determine start and end of processing period, should work with p2==0
1408lstart = datetime_miss
1409IF (PRESENT(start)) lstart = start
1410lend = itime(SIZE(itime))
1411! compute some quantities
1412maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
1413maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
1414minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
1415IF (time_definition == 0) THEN ! reference time
1416 lend = lend + timedelta_new(sec=maxp1)
1417ENDIF
1418! extend interval at the end in order to include all the data in case
1419! frac_valid<1; must use < and not <= in the DO WHILE loops some lines
1420! below in order to exclude the last full interval which would be empty
1421lend = lend + step
1422IF (lstart == datetime_miss) THEN ! autodetect
1423 lstart = itime(1)
1424! if autodetected, adjust to obtain real absolute start of data
1425 IF (time_definition == 0) THEN ! reference time
1426 lstart = lstart + timedelta_new(sec=minp1mp2)
1427 ELSE ! verification time
1428! go back to start of longest processing interval
1429 lstart = lstart - timedelta_new(sec=maxp2)
1430 ENDIF
1431! full_steps is effective only in analysis mode and when start is not
1432! specified (start by itself in analysis mode implies full_steps with
1433! respect to start instead of absolute full steps)
1434 IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN
1436 ENDIF
1437ENDIF
1438
1439#ifdef DEBUG
1440CALL l4f_log(l4f_debug, &
1442#endif
1443
1444! create output time and timerange lists
1445
1446IF (lforecast) THEN ! forecast mode
1447 IF (time_definition == 0) THEN ! reference time
1449
1450! apply start shift to timerange, not to reftime
1451! why did we use itime(SIZE(itime)) (last ref time)?
1452! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
1454! allow starting before first reftime but restrict dtstart to -steps
1455! dstart = MAX(0, dstart)
1457 DO p1 = steps + dstart, maxp1, steps
1458 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1459 ENDDO
1460
1461 ELSE ! verification time
1462
1463! verification time in forecast mode is the ugliest case, because we
1464! don't know a priori how many different (thus incompatible) reference
1465! times we have, so some assumption of regularity has to be made. For
1466! this reason msteps, the minimum step between two times, is
1467! computed. We choose to compute it as a difference between itime
1468! elements, it could be also computed as difference of itimerange%p1
1469! elements. But what if it is not modulo steps?
1470 mstepms = steps*1000_int_ll
1471 DO i = 2, SIZE(itime)
1473 IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN
1474 msteps = stepms/1000_int_ll
1476 ENDIF
1477 ENDDO
1478 msteps = mstepms/1000_int_ll
1479
1480 tmptime = lstart + step
1481 DO WHILE(tmptime < lend) ! < since lend has been extended
1482 CALL insert_unique(a_otime, tmptime)
1483 tmptime = tmptime + step
1484 ENDDO
1485
1486! in next loop, we used initially steps instead of msteps (see comment
1487! above), this gave much less combinations of time/timeranges with
1488! possible empty output; we start from msteps instead of steps in
1489! order to include partial initial intervals in case frac_valid<1;
1490! however, as a gemeral rule, for aggregation of forecasts it's better
1491! to use reference time
1492 DO p1 = msteps, maxp1, msteps
1493 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1494 ENDDO
1495
1496 ENDIF
1497
1498ELSE ! analysis mode
1499 tmptime = lstart + step
1500 DO WHILE(tmptime < lend) ! < since lend has been extended
1501 CALL insert_unique(a_otime, tmptime)
1502 tmptime = tmptime + step
1503 ENDDO
1504 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
1505
1506ENDIF
1507
1510otime => a_otime%array
1511otimerange => a_otimerange%array
1514! delete local objects keeping the contents
1517
1518#ifdef DEBUG
1519CALL l4f_log(l4f_debug, &
1520 'recompute_stat_proc_agg, output time and timerange: '//&
1522#endif
1523
1524IF (PRESENT(dtratio)) THEN
1525! count the possible i/o interval ratios
1526 DO k = 1, SIZE(itimerange)
1527 IF (itimerange(k)%p2 /= 0) &
1528 CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
1529 ENDDO
1531 dtratio => a_dtratio%array
1533! delete local object keeping the contents
1535
1536#ifdef DEBUG
1537 CALL l4f_log(l4f_debug, &
1539 ' possible aggregation ratios, from '// &
1541#endif
1542
1543 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1544 do_itimerange1: DO l = 1, SIZE(itimerange)
1545 IF (.NOT.mask_timerange(l)) cycle do_itimerange1
1546 do_itime1: DO k = 1, SIZE(itime)
1547 CALL time_timerange_get_period(itime(k), itimerange(l), &
1548 time_definition, pstart1, pend1, reftime1)
1549 do_otimerange1: DO j = 1, SIZE(otimerange)
1550 do_otime1: DO i = 1, SIZE(otime)
1551 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1552 time_definition, pstart2, pend2, reftime2)
1553 IF (lforecast) THEN
1554 IF (reftime1 /= reftime2) cycle do_otime1
1555 ENDIF
1556
1557 IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
1559 lmapper%it = k
1560 lmapper%itr = l
1561 lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
1562 n = append(map_ttr(i,j), lmapper)
1563 cycle do_itime1 ! can contribute only to a single interval
1564 ENDIF
1565 ENDDO do_otime1
1566 ENDDO do_otimerange1
1567 ENDDO do_itime1
1568 ENDDO do_itimerange1
1569
1570ELSE
1571
1572 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1573 do_itimerange2: DO l = 1, SIZE(itimerange)
1574 IF (.NOT.mask_timerange(l)) cycle do_itimerange2
1575 do_itime2: DO k = 1, SIZE(itime)
1576 CALL time_timerange_get_period(itime(k), itimerange(l), &
1577 time_definition, pstart1, pend1, reftime1)
1578 do_otimerange2: DO j = 1, SIZE(otimerange)
1579 do_otime2: DO i = 1, SIZE(otime)
1580 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1581 time_definition, pstart2, pend2, reftime2)
1582 IF (lforecast) THEN
1583 IF (reftime1 /= reftime2) cycle do_otime2
1584 ENDIF
1585
1586 IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
1587 IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
1588 lmapper%it = k
1589 lmapper%itr = l
1590 IF (pstart1 == pstart2) THEN
1591 lmapper%extra_info = 1 ! start of interval
1592 ELSE IF (pend1 == pend2) THEN
1593 lmapper%extra_info = 2 ! end of interval
1594 ELSE
1595 lmapper%extra_info = imiss
1596 ENDIF
1597 lmapper%time = pstart1 ! = pend1, order by time?
1598 n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
1599! no CYCLE, a single input can contribute to multiple output intervals
1600 ENDIF
1601 ENDDO do_otime2
1602 ENDDO do_otimerange2
1603 ENDDO do_itime2
1604 ENDDO do_itimerange2
1605
1606ENDIF
1607
1608END SUBROUTINE recompute_stat_proc_agg_common
1609
1610
1611SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
1612 max_step, weights)
1613TYPE(datetime),INTENT(in) :: vertime(:)
1614TYPE(datetime),INTENT(in) :: pstart
1615TYPE(datetime),INTENT(in) :: pend
1616LOGICAL,INTENT(in) :: time_mask(:)
1617TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
1618DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
1619
1620INTEGER :: i, nt
1621TYPE(datetime),ALLOCATABLE :: lvertime(:)
1622TYPE(datetime) :: half, nexthalf
1623INTEGER(kind=int_ll) :: dt, tdt
1624
1625nt = count(time_mask)
1626ALLOCATE(lvertime(nt))
1627lvertime = pack(vertime, mask=time_mask)
1628
1629IF (PRESENT(max_step)) THEN
1630! new way
1631! max_step = timedelta_0
1632! DO i = 1, nt - 1
1633! IF (lvertime(i+1) - lvertime(i) > max_step) &
1634! max_step = lvertime(i+1) - lvertime(i)
1635! ENDDO
1636! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
1637! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
1638! old way
1639 IF (nt == 1) THEN
1640 max_step = pend - pstart
1641 ELSE
1642 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
1643 max_step = half - pstart
1644 DO i = 2, nt - 1
1645 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
1646 IF (nexthalf - half > max_step) max_step = nexthalf - half
1647 half = nexthalf
1648 ENDDO
1649 IF (pend - half > max_step) max_step = pend - half
1650 ENDIF
1651
1652ENDIF
1653
1654IF (PRESENT(weights)) THEN
1655 IF (nt == 1) THEN
1656 weights(1) = 1.0
1657 ELSE
1659 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
1661 weights(1) = dble(dt)/dble(tdt)
1662 DO i = 2, nt - 1
1663 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
1665 weights(i) = dble(dt)/dble(tdt)
1666 half = nexthalf
1667 ENDDO
1669 weights(nt) = dble(dt)/dble(tdt)
1670 ENDIF
1671ENDIF
1672
1673END SUBROUTINE compute_stat_proc_agg_sw
1674
1675! get start of period, end of period and reference time from time,
1676! timerange, according to time_definition.
1677SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
1678 pstart, pend, reftime)
1679TYPE(datetime),INTENT(in) :: time
1680TYPE(vol7d_timerange),INTENT(in) :: timerange
1681INTEGER,INTENT(in) :: time_definition
1682TYPE(datetime),INTENT(out) :: reftime
1683TYPE(datetime),INTENT(out) :: pstart
1684TYPE(datetime),INTENT(out) :: pend
1685
1686TYPE(timedelta) :: p1, p2
1687
1688
1689p1 = timedelta_new(sec=timerange%p1) ! end of period
1690p2 = timedelta_new(sec=timerange%p2) ! length of period
1691
1693! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
1694 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
1695 pstart = datetime_miss
1696 pend = datetime_miss
1697 reftime = datetime_miss
1698 RETURN
1699ENDIF
1700
1701IF (time_definition == 0) THEN ! time == reference time
1702 reftime = time
1703 pend = time + p1
1704 pstart = pend - p2
1705ELSE IF (time_definition == 1) THEN ! time == verification time
1706 pend = time
1707 pstart = time - p2
1708 reftime = time - p1
1709ELSE
1710 pstart = datetime_miss
1711 pend = datetime_miss
1712 reftime = datetime_miss
1713ENDIF
1714
1715END SUBROUTINE time_timerange_get_period
1716
1717
1718! get start of period, end of period and reference time from time,
1719! timerange, according to time_definition. step is taken separately
1720! from timerange and can be "popular"
1721SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
1722 pstart, pend, reftime)
1723TYPE(datetime),INTENT(in) :: time
1724TYPE(vol7d_timerange),INTENT(in) :: timerange
1725TYPE(timedelta),INTENT(in) :: step
1726INTEGER,INTENT(in) :: time_definition
1727TYPE(datetime),INTENT(out) :: reftime
1728TYPE(datetime),INTENT(out) :: pstart
1729TYPE(datetime),INTENT(out) :: pend
1730
1731TYPE(timedelta) :: p1
1732
1733
1734p1 = timedelta_new(sec=timerange%p1) ! end of period
1735
1737! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
1738 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
1739 pstart = datetime_miss
1740 pend = datetime_miss
1741 reftime = datetime_miss
1742 RETURN
1743ENDIF
1744
1745IF (time_definition == 0) THEN ! time == reference time
1746 reftime = time
1747 pend = time + p1
1748 pstart = pend - step
1749ELSE IF (time_definition == 1) THEN ! time == verification time
1750 pend = time
1751 pstart = time - step
1752 reftime = time - p1
1753ELSE
1754 pstart = datetime_miss
1755 pend = datetime_miss
1756 reftime = datetime_miss
1757ENDIF
1758
1759END SUBROUTINE time_timerange_get_period_pop
1760
1761
1762! set time, timerange%p1, timerange%p2 according to pstart, pend,
1763! reftime and time_definition.
1764SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
1765 pstart, pend, reftime)
1766TYPE(datetime),INTENT(out) :: time
1767TYPE(vol7d_timerange),INTENT(inout) :: timerange
1768INTEGER,INTENT(in) :: time_definition
1769TYPE(datetime),INTENT(in) :: reftime
1770TYPE(datetime),INTENT(in) :: pstart
1771TYPE(datetime),INTENT(in) :: pend
1772
1773TYPE(timedelta) :: p1, p2
1774INTEGER(kind=int_ll) :: dmsec
1775
1776
1777IF (time_definition == 0) THEN ! time == reference time
1778 time = reftime
1779 p1 = pend - reftime
1780 p2 = pend - pstart
1781ELSE IF (time_definition == 1) THEN ! time == verification time
1782 time = pend
1783 p1 = pend - reftime
1784 p2 = pend - pstart
1785ELSE
1786 time = datetime_miss
1787ENDIF
1788
1789IF (time /= datetime_miss) THEN
1791 timerange%p1 = int(dmsec/1000_int_ll)
1793 timerange%p2 = int(dmsec/1000_int_ll)
1794ELSE
1795 timerange%p1 = imiss
1796 timerange%p2 = imiss
1797ENDIF
1798
1799END SUBROUTINE time_timerange_set_period
1800
1801
Restituiscono il valore dell'oggetto nella forma desiderata. Definition: datetime_class.F90:322 Functions that return a trimmed CHARACTER representation of the input variable. Definition: datetime_class.F90:349 Quick method to append an element to the array. Definition: stat_proc_engine.F90:365 Destructor for finalizing an array object. Definition: stat_proc_engine.F90:378 Method for inserting elements of the array at a desired position. Definition: stat_proc_engine.F90:356 Method for packing the array object reducing at a minimum the memory occupation, without destroying i... Definition: stat_proc_engine.F90:388 This module contains functions that are only for internal use of the library. Definition: stat_proc_engine.F90:211 Classe per la gestione di un volume completo di dati osservati. Definition: vol7d_class.F90:273 |