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