libsim Versione 7.2.6
|
◆ 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!> This module contains functions that are only for internal use of
953!! the library. It should not be used by user procedures because it is
954!! subject to change
955!! \ingroup vol7d
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
1089TYPE(timedelta) :: start_delta
1090
1091! compute length of cumulation step in seconds
1093
1094lstart = datetime_miss
1095IF (PRESENT(start)) lstart = start
1096lfull_steps = optio_log(full_steps)
1097
1098! create a mask of suitable timeranges
1099ALLOCATE(mask_timerange(SIZE(itimerange)))
1100mask_timerange(:) = itimerange(:)%timerange == stat_proc &
1101 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
1102 .AND. itimerange(:)%p1 >= 0 &
1103 .AND. itimerange(:)%p2 > 0
1104
1105IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges
1106 mask_timerange(:) = mask_timerange(:) .AND. &
1107 (itimerange(:)%p1 == 0 .OR. & ! all analyses
1110ENDIF
1111! mask_timerange includes all candidate timeranges
1112
1113nitr = count(mask_timerange)
1114ALLOCATE(f(nitr))
1115j = 1
1116DO i = 1, nitr
1117 DO WHILE(.NOT.mask_timerange(j))
1118 j = j + 1
1119 ENDDO
1120 f(i) = j
1121 j = j + 1
1122ENDDO
1123
1124! now we have to evaluate time/timerage pairs which do not need processing
1125ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
1126CALL compute_keep_tr()
1127
1128ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
1129map_tr(:,:,:,:,:) = imiss
1130
1131! scan through all possible combinations of time and timerange
1132DO dirtyrep = 1, 2
1133 IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
1138 ENDIF
1139 DO l = 1, SIZE(itime)
1140 DO k = 1, nitr
1141 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1142 time_definition, pstart2, pend2, reftime2)
1143
1144 DO j = 1, SIZE(itime)
1145 DO i = 1, nitr
1146 useful = .false.
1147 CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
1148 time_definition, pstart1, pend1, reftime1)
1149 tmptimerange = vol7d_timerange_new(timerange=stat_proc)
1150
1151 IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
1152 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
1153 CALL time_timerange_set_period(tmptime, tmptimerange, &
1154 time_definition, pend1, pend2, reftime2)
1155 IF (lfull_steps) THEN
1157 useful = .true.
1158 ENDIF
1159 ELSE
1160 useful = .true.
1161 ENDIF
1162
1163 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
1164 CALL time_timerange_set_period(tmptime, tmptimerange, &
1165 time_definition, pstart2, pstart1, pstart1)
1166 IF (lfull_steps) THEN
1168 useful = .true.
1169 ENDIF
1170 ELSE
1171 useful = .true.
1172 ENDIF
1173 ENDIF
1174
1175 ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
1176 IF (lfull_steps) THEN
1178! lstart shifts the interval for computing modulo step, this does not
1179! remove data before lstart but just shifts the phase
1180 start_delta = lstart-reftime2
1181 ELSE
1182 start_delta = timedelta_0
1183 ENDIF
1184 ENDIF
1185
1186 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
1187 CALL time_timerange_set_period(tmptime, tmptimerange, &
1188 time_definition, pend1, pend2, reftime2)
1189 IF (lfull_steps) THEN
1191 useful = .true.
1192 ENDIF
1193 ELSE
1194 useful = .true.
1195 ENDIF
1196! keep only data after lstart
1198 IF (lstart > pend1) useful = .false.
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, reftime2)
1204 IF (lfull_steps) THEN
1206 useful = .true.
1207 ENDIF
1208 ELSE
1209 useful = .true.
1210 ENDIF
1211! keep only data after lstart
1213 IF (lstart > pstart2) useful = .false.
1214 ENDIF
1215
1216 ENDIF
1217 ENDIF
1218 useful = useful .AND. tmptime /= datetime_miss .AND. &
1219 tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
1220
1221 IF (useful) THEN ! add a_otime, a_otimerange
1222 map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
1223 map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
1224 ENDIF
1225 ENDDO
1226 ENDDO
1227 ENDDO
1228 ENDDO
1229ENDDO
1230
1231! we have to repeat the computation with sorted arrays
1232CALL compute_keep_tr()
1233
1234otime => a_otime%array
1235otimerange => a_otimerange%array
1236! delete local objects keeping the contents
1239
1240#ifdef DEBUG
1241CALL l4f_log(l4f_debug, &
1246CALL l4f_log(l4f_debug, &
1249CALL l4f_log(l4f_debug, &
1251CALL l4f_log(l4f_debug, &
1253CALL l4f_log(l4f_debug, &
1255CALL l4f_log(l4f_debug, &
1257#endif
1258
1259CONTAINS
1260
1261SUBROUTINE compute_keep_tr()
1262INTEGER :: start_deltas
1263
1264keep_tr(:,:,:) = imiss
1265DO l = 1, SIZE(itime)
1266 itrloop: DO k = 1, nitr
1267 IF (itimerange(f(k))%p2 == steps) THEN
1268 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1269 time_definition, pstart2, pend2, reftime2)
1270 useful = .false.
1271! keep only data after lstart
1273 IF (lstart > pstart2) cycle itrloop
1274 ENDIF
1275 IF (reftime2 == pend2) THEN ! analysis
1278 useful = .true.
1279 ENDIF
1280 ELSE IF (lfull_steps) THEN
1282 useful = .true.
1283 ENDIF
1284 ELSE
1285 useful = .true.
1286 ENDIF
1287 ELSE ! forecast
1288 IF (lfull_steps) THEN
1289! same as above for start_delta, but in seconds and not in timerange/timedelta units
1291 start_deltas = timedelta_getamsec(lstart-reftime2)/1000_int_ll
1292 ELSE
1293 start_deltas = 0
1294 ENDIF
1296 useful = .true.
1297 ENDIF
1298 ELSE
1299 useful = .true.
1300 ENDIF
1301 ENDIF
1302 IF (useful) THEN
1303! CALL time_timerange_set_period(tmptime, tmptimerange, &
1304! time_definition, pstart2, pend2, reftime2)
1305 keep_tr(k,l,1) = append_unique(a_otime, itime(l))
1306 keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k)))
1307 ENDIF
1308 ENDIF
1309 ENDDO itrloop
1310ENDDO
1311
1312END SUBROUTINE compute_keep_tr
1313
1314END SUBROUTINE recompute_stat_proc_diff_common
1315
1316
1317! common operations for statistical processing by metamorphosis
1318SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
1319 otimerange, map_tr)
1320INTEGER,INTENT(in) :: istat_proc
1321TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1322INTEGER,INTENT(in) :: ostat_proc
1323TYPE(vol7d_timerange),POINTER :: otimerange(:)
1324INTEGER,POINTER :: map_tr(:)
1325
1326INTEGER :: i
1327LOGICAL :: tr_mask(SIZE(itimerange))
1328
1329IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1330 ALLOCATE(otimerange(0), map_tr(0))
1331 RETURN
1332ENDIF
1333
1334! useful timeranges
1335tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
1336 .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
1337ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
1338
1339otimerange = pack(itimerange, mask=tr_mask)
1340otimerange(:)%timerange = ostat_proc
1341map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
1342
1343END SUBROUTINE compute_stat_proc_metamorph_common
1344
1345
1346! common operations for statistical processing by aggregation
1347SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
1348 step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
1349TYPE(datetime),INTENT(in) :: itime(:)
1350TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1351INTEGER,INTENT(in) :: stat_proc
1352INTEGER,INTENT(in) :: tri
1353TYPE(timedelta),INTENT(in) :: step
1354INTEGER,INTENT(in) :: time_definition
1355TYPE(datetime),POINTER :: otime(:)
1356TYPE(vol7d_timerange),POINTER :: otimerange(:)
1357TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
1358INTEGER,POINTER,OPTIONAL :: dtratio(:)
1359TYPE(datetime),INTENT(in),OPTIONAL :: start
1360LOGICAL,INTENT(in),OPTIONAL :: full_steps
1361
1362INTEGER :: i, j, k, l, na, nf, n
1363INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
1364INTEGER(kind=int_ll) :: stepms, mstepms
1365LOGICAL :: lforecast
1366TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1367TYPE(arrayof_datetime) :: a_otime
1368TYPE(arrayof_vol7d_timerange) :: a_otimerange
1369TYPE(arrayof_integer) :: a_dtratio
1370LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
1371TYPE(ttr_mapper) :: lmapper
1372CHARACTER(len=8) :: env_var
1373LOGICAL :: climat_behavior
1374
1375
1376! enable bad behavior for climat database (only for instantaneous data)
1377env_var = ''
1378CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
1379climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
1380
1381! compute length of cumulation step in seconds
1383
1384! create a mask of suitable timeranges
1385ALLOCATE(mask_timerange(SIZE(itimerange)))
1386mask_timerange(:) = itimerange(:)%timerange == tri &
1387 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
1388
1389IF (PRESENT(dtratio)) THEN
1390 WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0)
1391 mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0
1392 ELSEWHERE
1393 mask_timerange(:) = .false.
1394 END WHERE
1395ELSE ! instantaneous
1396 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
1397ENDIF
1398
1399#ifdef DEBUG
1400CALL l4f_log(l4f_debug, &
1401 '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
1402 t2c(count(mask_timerange)))
1403#endif
1404
1405! euristically determine whether we are dealing with an
1406! analysis/observation or a forecast dataset
1407na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
1408nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
1409lforecast = nf >= na
1410#ifdef DEBUG
1411CALL l4f_log(l4f_debug, &
1413#endif
1414
1415! keep only timeranges of one type (really necessary?)
1416IF (lforecast) THEN
1417 CALL l4f_log(l4f_info, &
1418 'recompute_stat_proc_agg, processing in forecast mode')
1419ELSE
1420 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
1421 CALL l4f_log(l4f_info, &
1422 'recompute_stat_proc_agg, processing in analysis mode')
1423ENDIF
1424
1425#ifdef DEBUG
1426CALL l4f_log(l4f_debug, &
1427 '(re)compute_stat_proc_agg, number of useful timeranges: '// &
1428 t2c(count(mask_timerange)))
1429#endif
1430
1431IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1432 ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
1433 IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
1434 RETURN
1435ENDIF
1436
1437! determine start and end of processing period, should work with p2==0
1438lstart = datetime_miss
1439IF (PRESENT(start)) lstart = start
1440lend = itime(SIZE(itime))
1441! compute some quantities
1442maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
1443maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
1444minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
1445IF (time_definition == 0) THEN ! reference time
1446 lend = lend + timedelta_new(sec=maxp1)
1447ENDIF
1448! extend interval at the end in order to include all the data in case
1449! frac_valid<1; must use < and not <= in the DO WHILE loops some lines
1450! below in order to exclude the last full interval which would be empty
1451lend = lend + step
1452IF (lstart == datetime_miss) THEN ! autodetect
1453 lstart = itime(1)
1454! if autodetected, adjust to obtain real absolute start of data
1455 IF (time_definition == 0) THEN ! reference time
1456 lstart = lstart + timedelta_new(sec=minp1mp2)
1457 ELSE ! verification time
1458! go back to start of longest processing interval
1459 lstart = lstart - timedelta_new(sec=maxp2)
1460 ENDIF
1461! apply full_steps in analysis mode and when start is not specified
1462! (start by itself in analysis mode implies full_steps with respect to
1463! start instead of absolute full steps), todo in forecast mode
1464 IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN
1466 ENDIF
1467ENDIF
1468
1469#ifdef DEBUG
1470CALL l4f_log(l4f_debug, &
1472#endif
1473
1474! create output time and timerange lists
1475
1476IF (lforecast) THEN ! forecast mode
1477 IF (time_definition == 0) THEN ! reference time
1479
1480! apply start shift to timerange, not to reftime
1481! why did we use itime(SIZE(itime)) (last ref time)?
1482! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
1484! allow starting before first reftime but restrict dtstart to -steps
1485! dstart = MAX(0, dstart)
1487 DO p1 = steps + dstart, maxp1, steps
1488 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1489 ENDDO
1490
1491 ELSE ! verification time
1492
1493! verification time in forecast mode is the ugliest case, because we
1494! don't know a priori how many different (thus incompatible) reference
1495! times we have, so some assumption of regularity has to be made. For
1496! this reason msteps, the minimum step between two times, is
1497! computed. We choose to compute it as a difference between itime
1498! elements, it could be also computed as difference of itimerange%p1
1499! elements. But what if it is not modulo steps?
1500 mstepms = steps*1000_int_ll
1501 DO i = 2, SIZE(itime)
1503 IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN
1504 msteps = stepms/1000_int_ll
1506 ENDIF
1507 ENDDO
1508 msteps = mstepms/1000_int_ll
1509
1510 tmptime = lstart + step
1511 DO WHILE(tmptime < lend) ! < since lend has been extended
1512 CALL insert_unique(a_otime, tmptime)
1513 tmptime = tmptime + step
1514 ENDDO
1515
1516! in next loop, we used initially steps instead of msteps (see comment
1517! above), this gave much less combinations of time/timeranges with
1518! possible empty output; we start from msteps instead of steps in
1519! order to include partial initial intervals in case frac_valid<1;
1520! however, as a gemeral rule, for aggregation of forecasts it's better
1521! to use reference time
1522 DO p1 = msteps, maxp1, msteps
1523 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1524 ENDDO
1525
1526 ENDIF
1527
1528ELSE ! analysis mode
1529 tmptime = lstart + step
1530 DO WHILE(tmptime < lend) ! < since lend has been extended
1531 CALL insert_unique(a_otime, tmptime)
1532 tmptime = tmptime + step
1533 ENDDO
1534 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
1535
1536ENDIF
1537
1540otime => a_otime%array
1541otimerange => a_otimerange%array
1544! delete local objects keeping the contents
1547
1548#ifdef DEBUG
1549CALL l4f_log(l4f_debug, &
1550 'recompute_stat_proc_agg, output time and timerange: '//&
1552#endif
1553
1554IF (PRESENT(dtratio)) THEN
1555! count the possible i/o interval ratios
1556 DO k = 1, SIZE(itimerange)
1557 IF (itimerange(k)%p2 /= 0) &
1558 CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
1559 ENDDO
1561 dtratio => a_dtratio%array
1563! delete local object keeping the contents
1565
1566#ifdef DEBUG
1567 CALL l4f_log(l4f_debug, &
1569 ' possible aggregation ratios, from '// &
1571#endif
1572
1573 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1574 do_itimerange1: DO l = 1, SIZE(itimerange)
1575 IF (.NOT.mask_timerange(l)) cycle do_itimerange1
1576 do_itime1: DO k = 1, SIZE(itime)
1577 CALL time_timerange_get_period(itime(k), itimerange(l), &
1578 time_definition, pstart1, pend1, reftime1)
1579 do_otimerange1: DO j = 1, SIZE(otimerange)
1580 do_otime1: DO i = 1, SIZE(otime)
1581 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1582 time_definition, pstart2, pend2, reftime2)
1583 IF (lforecast) THEN
1584 IF (reftime1 /= reftime2) cycle do_otime1
1585 ENDIF
1586
1587 IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
1589 lmapper%it = k
1590 lmapper%itr = l
1591 lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
1592 n = append(map_ttr(i,j), lmapper)
1593 cycle do_itime1 ! can contribute only to a single interval
1594 ENDIF
1595 ENDDO do_otime1
1596 ENDDO do_otimerange1
1597 ENDDO do_itime1
1598 ENDDO do_itimerange1
1599
1600ELSE
1601
1602 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1603 do_itimerange2: DO l = 1, SIZE(itimerange)
1604 IF (.NOT.mask_timerange(l)) cycle do_itimerange2
1605 do_itime2: DO k = 1, SIZE(itime)
1606 CALL time_timerange_get_period(itime(k), itimerange(l), &
1607 time_definition, pstart1, pend1, reftime1)
1608 do_otimerange2: DO j = 1, SIZE(otimerange)
1609 do_otime2: DO i = 1, SIZE(otime)
1610 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1611 time_definition, pstart2, pend2, reftime2)
1612 IF (lforecast) THEN
1613 IF (reftime1 /= reftime2) cycle do_otime2
1614 ENDIF
1615
1616 IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
1617 IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
1618 lmapper%it = k
1619 lmapper%itr = l
1620 IF (pstart1 == pstart2) THEN
1621 lmapper%extra_info = 1 ! start of interval
1622 ELSE IF (pend1 == pend2) THEN
1623 lmapper%extra_info = 2 ! end of interval
1624 ELSE
1625 lmapper%extra_info = imiss
1626 ENDIF
1627 lmapper%time = pstart1 ! = pend1, order by time?
1628 n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
1629! no CYCLE, a single input can contribute to multiple output intervals
1630 ENDIF
1631 ENDDO do_otime2
1632 ENDDO do_otimerange2
1633 ENDDO do_itime2
1634 ENDDO do_itimerange2
1635
1636ENDIF
1637
1638END SUBROUTINE recompute_stat_proc_agg_common
1639
1640
1641SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
1642 max_step, weights)
1643TYPE(datetime),INTENT(in) :: vertime(:)
1644TYPE(datetime),INTENT(in) :: pstart
1645TYPE(datetime),INTENT(in) :: pend
1646LOGICAL,INTENT(in) :: time_mask(:)
1647TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
1648DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
1649
1650INTEGER :: i, nt
1651TYPE(datetime),ALLOCATABLE :: lvertime(:)
1652TYPE(datetime) :: half, nexthalf
1653INTEGER(kind=int_ll) :: dt, tdt
1654
1655nt = count(time_mask)
1656ALLOCATE(lvertime(nt))
1657lvertime = pack(vertime, mask=time_mask)
1658
1659IF (PRESENT(max_step)) THEN
1660! new way
1661! max_step = timedelta_0
1662! DO i = 1, nt - 1
1663! IF (lvertime(i+1) - lvertime(i) > max_step) &
1664! max_step = lvertime(i+1) - lvertime(i)
1665! ENDDO
1666! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
1667! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
1668! old way
1669 IF (nt == 1) THEN
1670 max_step = pend - pstart
1671 ELSE
1672 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
1673 max_step = half - pstart
1674 DO i = 2, nt - 1
1675 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
1676 IF (nexthalf - half > max_step) max_step = nexthalf - half
1677 half = nexthalf
1678 ENDDO
1679 IF (pend - half > max_step) max_step = pend - half
1680 ENDIF
1681
1682ENDIF
1683
1684IF (PRESENT(weights)) THEN
1685 IF (nt == 1) THEN
1686 weights(1) = 1.0
1687 ELSE
1689 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
1691 weights(1) = dble(dt)/dble(tdt)
1692 DO i = 2, nt - 1
1693 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
1695 weights(i) = dble(dt)/dble(tdt)
1696 half = nexthalf
1697 ENDDO
1699 weights(nt) = dble(dt)/dble(tdt)
1700 ENDIF
1701ENDIF
1702
1703END SUBROUTINE compute_stat_proc_agg_sw
1704
1705! get start of period, end of period and reference time from time,
1706! timerange, according to time_definition.
1707SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
1708 pstart, pend, reftime)
1709TYPE(datetime),INTENT(in) :: time
1710TYPE(vol7d_timerange),INTENT(in) :: timerange
1711INTEGER,INTENT(in) :: time_definition
1712TYPE(datetime),INTENT(out) :: reftime
1713TYPE(datetime),INTENT(out) :: pstart
1714TYPE(datetime),INTENT(out) :: pend
1715
1716TYPE(timedelta) :: p1, p2
1717
1718
1719p1 = timedelta_new(sec=timerange%p1) ! end of period
1720p2 = timedelta_new(sec=timerange%p2) ! length of period
1721
1723! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
1724 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
1725 pstart = datetime_miss
1726 pend = datetime_miss
1727 reftime = datetime_miss
1728 RETURN
1729ENDIF
1730
1731IF (time_definition == 0) THEN ! time == reference time
1732 reftime = time
1733 pend = time + p1
1734 pstart = pend - p2
1735ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time
1736 pend = time
1737 pstart = time - p2
1738 reftime = time - p1
1739ELSE
1740 pstart = datetime_miss
1741 pend = datetime_miss
1742 reftime = datetime_miss
1743ENDIF
1744
1745END SUBROUTINE time_timerange_get_period
1746
1747
1748! get start of period, end of period and reference time from time,
1749! timerange, according to time_definition. step is taken separately
1750! from timerange and can be "popular"
1751SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
1752 pstart, pend, reftime)
1753TYPE(datetime),INTENT(in) :: time
1754TYPE(vol7d_timerange),INTENT(in) :: timerange
1755TYPE(timedelta),INTENT(in) :: step
1756INTEGER,INTENT(in) :: time_definition
1757TYPE(datetime),INTENT(out) :: reftime
1758TYPE(datetime),INTENT(out) :: pstart
1759TYPE(datetime),INTENT(out) :: pend
1760
1761TYPE(timedelta) :: p1
1762
1763
1764p1 = timedelta_new(sec=timerange%p1) ! end of period
1765
1767! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
1768 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
1769 pstart = datetime_miss
1770 pend = datetime_miss
1771 reftime = datetime_miss
1772 RETURN
1773ENDIF
1774
1775IF (time_definition == 0) THEN ! time == reference time
1776 reftime = time
1777 pend = time + p1
1778 pstart = pend - step
1779ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time
1780 pend = time
1781 pstart = time - step
1782 reftime = time - p1
1783ELSE
1784 pstart = datetime_miss
1785 pend = datetime_miss
1786 reftime = datetime_miss
1787ENDIF
1788
1789END SUBROUTINE time_timerange_get_period_pop
1790
1791
1792! set time, timerange%p1, timerange%p2 according to pstart, pend,
1793! reftime and time_definition.
1794SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
1795 pstart, pend, reftime)
1796TYPE(datetime),INTENT(out) :: time
1797TYPE(vol7d_timerange),INTENT(inout) :: timerange
1798INTEGER,INTENT(in) :: time_definition
1799TYPE(datetime),INTENT(in) :: reftime
1800TYPE(datetime),INTENT(in) :: pstart
1801TYPE(datetime),INTENT(in) :: pend
1802
1803TYPE(timedelta) :: p1, p2
1804INTEGER(kind=int_ll) :: dmsec
1805
1806
1807IF (time_definition == 0) THEN ! time == reference time
1808 time = reftime
1809 p1 = pend - reftime
1810 p2 = pend - pstart
1811ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time
1812 time = pend
1813 p1 = pend - reftime
1814 p2 = pend - pstart
1815ELSE
1816 time = datetime_miss
1817ENDIF
1818
1819IF (time /= datetime_miss) THEN
1821 timerange%p1 = int(dmsec/1000_int_ll)
1823 timerange%p2 = int(dmsec/1000_int_ll)
1824ELSE
1825 timerange%p1 = imiss
1826 timerange%p2 = imiss
1827ENDIF
1828
1829END SUBROUTINE time_timerange_set_period
1830
1831
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 |