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