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