libsim Versione 7.2.6

◆ map_distinct_ttr_mapper()

integer function, dimension(size(vect)) map_distinct_ttr_mapper ( type(ttr_mapper), dimension(:), intent(in) vect,
logical, dimension(:), intent(in), optional mask,
logical, intent(in), optional back )

map distinct

Definizione alla linea 1114 del file stat_proc_engine.F90.

1115! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1116! authors:
1117! Davide Cesari <dcesari@arpa.emr.it>
1118! Paolo Patruno <ppatruno@arpa.emr.it>
1119
1120! This program is free software; you can redistribute it and/or
1121! modify it under the terms of the GNU General Public License as
1122! published by the Free Software Foundation; either version 2 of
1123! the License, or (at your option) any later version.
1124
1125! This program is distributed in the hope that it will be useful,
1126! but WITHOUT ANY WARRANTY; without even the implied warranty of
1127! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1128! GNU General Public License for more details.
1129
1130! You should have received a copy of the GNU General Public License
1131! along with this program. If not, see <http://www.gnu.org/licenses/>.
1132#include "config.h"
1133!> This module contains functions that are only for internal use of
1134!! the library. It should not be used by user procedures because it is
1135!! subject to change
1136!! \ingroup vol7d
1137MODULE stat_proc_engine
1139USE vol7d_class
1140IMPLICIT NONE
1141
1142! per ogni elemento i,j dell'output, definire n elementi di input ad
1143! esso contribuenti (arrayof_ttr_mapper) con le seguenti informazioni
1144TYPE ttr_mapper
1145 INTEGER :: it=imiss ! k
1146 INTEGER :: itr=imiss ! l
1147 INTEGER :: extra_info=imiss ! dtratio for intervals, extreme for point in time
1148 TYPE(datetime) :: time=datetime_miss ! time for point in time
1149END TYPE ttr_mapper
1150
1151INTERFACE OPERATOR (==)
1152 MODULE PROCEDURE ttr_mapper_eq
1153END INTERFACE
1154
1155INTERFACE OPERATOR (/=)
1156 MODULE PROCEDURE ttr_mapper_ne
1157END INTERFACE
1158
1159INTERFACE OPERATOR (>)
1160 MODULE PROCEDURE ttr_mapper_gt
1161END INTERFACE
1162
1163INTERFACE OPERATOR (<)
1164 MODULE PROCEDURE ttr_mapper_lt
1165END INTERFACE
1166
1167INTERFACE OPERATOR (>=)
1168 MODULE PROCEDURE ttr_mapper_ge
1169END INTERFACE
1170
1171INTERFACE OPERATOR (<=)
1172 MODULE PROCEDURE ttr_mapper_le
1173END INTERFACE
1174
1175#undef VOL7D_POLY_TYPE
1176#undef VOL7D_POLY_TYPES
1177#undef ENABLE_SORT
1178#define VOL7D_POLY_TYPE TYPE(ttr_mapper)
1179#define VOL7D_POLY_TYPES _ttr_mapper
1180#define ENABLE_SORT
1181#include "array_utilities_pre.F90"
1182
1183#define ARRAYOF_ORIGTYPE TYPE(ttr_mapper)
1184#define ARRAYOF_TYPE arrayof_ttr_mapper
1185#define ARRAYOF_ORIGEQ 1
1186#define ARRAYOF_ORIGGT 1
1187#include "arrayof_pre.F90"
1188! from arrayof
1189
1190
1191CONTAINS
1192
1193! simplified comparisons only on time
1194ELEMENTAL FUNCTION ttr_mapper_eq(this, that) RESULT(res)
1195TYPE(ttr_mapper),INTENT(IN) :: this, that
1196LOGICAL :: res
1197
1198res = this%time == that%time
1199
1200END FUNCTION ttr_mapper_eq
1201
1202ELEMENTAL FUNCTION ttr_mapper_ne(this, that) RESULT(res)
1203TYPE(ttr_mapper),INTENT(IN) :: this, that
1204LOGICAL :: res
1205
1206res = this%time /= that%time
1207
1208END FUNCTION ttr_mapper_ne
1209
1210ELEMENTAL FUNCTION ttr_mapper_gt(this, that) RESULT(res)
1211TYPE(ttr_mapper),INTENT(IN) :: this, that
1212LOGICAL :: res
1213
1214res = this%time > that%time
1215
1216END FUNCTION ttr_mapper_gt
1217
1218ELEMENTAL FUNCTION ttr_mapper_lt(this, that) RESULT(res)
1219TYPE(ttr_mapper),INTENT(IN) :: this, that
1220LOGICAL :: res
1221
1222res = this%time < that%time
1223
1224END FUNCTION ttr_mapper_lt
1225
1226ELEMENTAL FUNCTION ttr_mapper_ge(this, that) RESULT(res)
1227TYPE(ttr_mapper),INTENT(IN) :: this, that
1228LOGICAL :: res
1229
1230res = this%time >= that%time
1231
1232END FUNCTION ttr_mapper_ge
1233
1234ELEMENTAL FUNCTION ttr_mapper_le(this, that) RESULT(res)
1235TYPE(ttr_mapper),INTENT(IN) :: this, that
1236LOGICAL :: res
1237
1238res = this%time <= that%time
1239
1240END FUNCTION ttr_mapper_le
1241
1242#include "arrayof_post.F90"
1243#include "array_utilities_inc.F90"
1244
1245
1246! common operations for statistical processing by differences
1247SUBROUTINE recompute_stat_proc_diff_common(itime, itimerange, stat_proc, step, &
1248 otime, otimerange, map_tr, f, keep_tr, time_definition, full_steps, &
1249 start)
1250TYPE(datetime),INTENT(in) :: itime(:)
1251TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1252INTEGER,INTENT(in) :: stat_proc
1253TYPE(timedelta),INTENT(in) :: step
1254TYPE(datetime),POINTER :: otime(:)
1255TYPE(vol7d_timerange),POINTER :: otimerange(:)
1256INTEGER,ALLOCATABLE,INTENT(out) :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
1257INTEGER :: nitr
1258LOGICAL,ALLOCATABLE :: mask_timerange(:)
1259INTEGER,INTENT(in) :: time_definition
1260LOGICAL,INTENT(in),OPTIONAL :: full_steps
1261TYPE(datetime),INTENT(in),OPTIONAL :: start
1262
1263INTEGER :: i, j, k, l, dirtyrep
1264INTEGER :: steps
1265LOGICAL :: lfull_steps, useful
1266TYPE(datetime) :: lstart, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1267TYPE(vol7d_timerange) :: tmptimerange
1268TYPE(arrayof_datetime) :: a_otime
1269TYPE(arrayof_vol7d_timerange) :: a_otimerange
1270TYPE(timedelta) :: start_delta
1271
1272! compute length of cumulation step in seconds
1273CALL getval(step, asec=steps)
1274
1275lstart = datetime_miss
1276IF (PRESENT(start)) lstart = start
1277lfull_steps = optio_log(full_steps)
1278
1279! create a mask of suitable timeranges
1280ALLOCATE(mask_timerange(SIZE(itimerange)))
1281mask_timerange(:) = itimerange(:)%timerange == stat_proc &
1282 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
1283 .AND. itimerange(:)%p1 >= 0 &
1284 .AND. itimerange(:)%p2 > 0
1285
1286IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges
1287 mask_timerange(:) = mask_timerange(:) .AND. &
1288 (itimerange(:)%p1 == 0 .OR. & ! all analyses
1289 mod(itimerange(:)%p1, steps) == 0 .OR. & ! end time is mod step
1290 mod(itimerange(:)%p1 - itimerange(:)%p2, steps) == 0) ! start time is mod step
1291ENDIF
1292! mask_timerange includes all candidate timeranges
1293
1294nitr = count(mask_timerange)
1295ALLOCATE(f(nitr))
1296j = 1
1297DO i = 1, nitr
1298 DO WHILE(.NOT.mask_timerange(j))
1299 j = j + 1
1300 ENDDO
1301 f(i) = j
1302 j = j + 1
1303ENDDO
1304
1305! now we have to evaluate time/timerage pairs which do not need processing
1306ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
1307CALL compute_keep_tr()
1308
1309ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
1310map_tr(:,:,:,:,:) = imiss
1311
1312! scan through all possible combinations of time and timerange
1313DO dirtyrep = 1, 2
1314 IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
1315 CALL packarray(a_otime)
1316 CALL packarray(a_otimerange)
1317 CALL sort(a_otime%array)
1318 CALL sort(a_otimerange%array)
1319 ENDIF
1320 DO l = 1, SIZE(itime)
1321 DO k = 1, nitr
1322 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1323 time_definition, pstart2, pend2, reftime2)
1324
1325 DO j = 1, SIZE(itime)
1326 DO i = 1, nitr
1327 useful = .false.
1328 CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
1329 time_definition, pstart1, pend1, reftime1)
1330 tmptimerange = vol7d_timerange_new(timerange=stat_proc)
1331
1332 IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
1333 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
1334 CALL time_timerange_set_period(tmptime, tmptimerange, &
1335 time_definition, pend1, pend2, reftime2)
1336 IF (lfull_steps) THEN
1337 IF (mod(reftime2, step) == timedelta_0) THEN
1338 useful = .true.
1339 ENDIF
1340 ELSE
1341 useful = .true.
1342 ENDIF
1343
1344 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
1345 CALL time_timerange_set_period(tmptime, tmptimerange, &
1346 time_definition, pstart2, pstart1, pstart1)
1347 IF (lfull_steps) THEN
1348 IF (mod(pstart1, step) == timedelta_0) THEN
1349 useful = .true.
1350 ENDIF
1351 ELSE
1352 useful = .true.
1353 ENDIF
1354 ENDIF
1355
1356 ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
1357 IF (lfull_steps) THEN
1358 IF (c_e(lstart)) THEN
1359! lstart shifts the interval for computing modulo step, this does not
1360! remove data before lstart but just shifts the phase
1361 start_delta = lstart-reftime2
1362 ELSE
1363 start_delta = timedelta_0
1364 ENDIF
1365 ENDIF
1366
1367 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
1368 CALL time_timerange_set_period(tmptime, tmptimerange, &
1369 time_definition, pend1, pend2, reftime2)
1370 IF (lfull_steps) THEN
1371 IF (mod(pend2-reftime2-start_delta, step) == timedelta_0) THEN
1372 useful = .true.
1373 ENDIF
1374 ELSE
1375 useful = .true.
1376 ENDIF
1377! keep only data after lstart
1378 IF (c_e(lstart)) THEN
1379 IF (lstart > pend1) useful = .false.
1380 ENDIF
1381
1382 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! |-=
1383 CALL time_timerange_set_period(tmptime, tmptimerange, &
1384 time_definition, pstart2, pstart1, reftime2)
1385 IF (lfull_steps) THEN
1386 IF (mod(pstart1-reftime2-start_delta, step) == timedelta_0) THEN
1387 useful = .true.
1388 ENDIF
1389 ELSE
1390 useful = .true.
1391 ENDIF
1392! keep only data after lstart
1393 IF (c_e(lstart)) THEN
1394 IF (lstart > pstart2) useful = .false.
1395 ENDIF
1396
1397 ENDIF
1398 ENDIF
1399 useful = useful .AND. tmptime /= datetime_miss .AND. &
1400 tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
1401
1402 IF (useful) THEN ! add a_otime, a_otimerange
1403 map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
1404 map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
1405 ENDIF
1406 ENDDO
1407 ENDDO
1408 ENDDO
1409 ENDDO
1410ENDDO
1411
1412! we have to repeat the computation with sorted arrays
1413CALL compute_keep_tr()
1414
1415otime => a_otime%array
1416otimerange => a_otimerange%array
1417! delete local objects keeping the contents
1418CALL delete(a_otime, nodealloc=.true.)
1419CALL delete(a_otimerange, nodealloc=.true.)
1420
1421#ifdef DEBUG
1422CALL l4f_log(l4f_debug, &
1423 'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr,1)))//', '// &
1424 t2c((SIZE(map_tr,2)))//', '// &
1425 t2c((SIZE(map_tr,3)))//', '// &
1426 t2c((SIZE(map_tr,4))))
1427CALL l4f_log(l4f_debug, &
1428 'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr))/2)//', '// &
1429 t2c(count(c_e(map_tr))/2))
1430CALL l4f_log(l4f_debug, &
1431 'recompute_stat_proc_diff, nitr: '//t2c(nitr))
1432CALL l4f_log(l4f_debug, &
1433 'recompute_stat_proc_diff, good timeranges: '//t2c(count(c_e(keep_tr))/2))
1434CALL l4f_log(l4f_debug, &
1435 'recompute_stat_proc_diff, output times: '//t2c(SIZE(otime)))
1436CALL l4f_log(l4f_debug, &
1437 'recompute_stat_proc_diff, output timeranges: '//t2c(SIZE(otimerange)))
1438#endif
1439
1440CONTAINS
1441
1442SUBROUTINE compute_keep_tr()
1443INTEGER :: start_deltas
1444
1445keep_tr(:,:,:) = imiss
1446DO l = 1, SIZE(itime)
1447 itrloop: DO k = 1, nitr
1448 IF (itimerange(f(k))%p2 == steps) THEN
1449 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1450 time_definition, pstart2, pend2, reftime2)
1451 useful = .false.
1452! keep only data after lstart
1453 IF (c_e(lstart)) THEN
1454 IF (lstart > pstart2) cycle itrloop
1455 ENDIF
1456 IF (reftime2 == pend2) THEN ! analysis
1457 IF (c_e(lstart)) THEN ! in analysis mode start wins over full_steps
1458 IF (mod(reftime2-lstart, step) == timedelta_0) THEN
1459 useful = .true.
1460 ENDIF
1461 ELSE IF (lfull_steps) THEN
1462 IF (mod(reftime2, step) == timedelta_0) THEN
1463 useful = .true.
1464 ENDIF
1465 ELSE
1466 useful = .true.
1467 ENDIF
1468 ELSE ! forecast
1469 IF (lfull_steps) THEN
1470! same as above for start_delta, but in seconds and not in timerange/timedelta units
1471 IF (c_e(lstart)) THEN
1472 start_deltas = timedelta_getamsec(lstart-reftime2)/1000_int_ll
1473 ELSE
1474 start_deltas = 0
1475 ENDIF
1476 IF (mod(itimerange(f(k))%p1 - start_deltas, steps) == 0) THEN
1477 useful = .true.
1478 ENDIF
1479 ELSE
1480 useful = .true.
1481 ENDIF
1482 ENDIF
1483 IF (useful) THEN
1484! CALL time_timerange_set_period(tmptime, tmptimerange, &
1485! time_definition, pstart2, pend2, reftime2)
1486 keep_tr(k,l,1) = append_unique(a_otime, itime(l))
1487 keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k)))
1488 ENDIF
1489 ENDIF
1490 ENDDO itrloop
1491ENDDO
1492
1493END SUBROUTINE compute_keep_tr
1494
1495END SUBROUTINE recompute_stat_proc_diff_common
1496
1497
1498! common operations for statistical processing by metamorphosis
1499SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
1500 otimerange, map_tr)
1501INTEGER,INTENT(in) :: istat_proc
1502TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1503INTEGER,INTENT(in) :: ostat_proc
1504TYPE(vol7d_timerange),POINTER :: otimerange(:)
1505INTEGER,POINTER :: map_tr(:)
1506
1507INTEGER :: i
1508LOGICAL :: tr_mask(SIZE(itimerange))
1509
1510IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1511 ALLOCATE(otimerange(0), map_tr(0))
1512 RETURN
1513ENDIF
1514
1515! useful timeranges
1516tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
1517 .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
1518ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
1519
1520otimerange = pack(itimerange, mask=tr_mask)
1521otimerange(:)%timerange = ostat_proc
1522map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
1523
1524END SUBROUTINE compute_stat_proc_metamorph_common
1525
1526
1527! common operations for statistical processing by aggregation
1528SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
1529 step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
1530TYPE(datetime),INTENT(in) :: itime(:)
1531TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1532INTEGER,INTENT(in) :: stat_proc
1533INTEGER,INTENT(in) :: tri
1534TYPE(timedelta),INTENT(in) :: step
1535INTEGER,INTENT(in) :: time_definition
1536TYPE(datetime),POINTER :: otime(:)
1537TYPE(vol7d_timerange),POINTER :: otimerange(:)
1538TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
1539INTEGER,POINTER,OPTIONAL :: dtratio(:)
1540TYPE(datetime),INTENT(in),OPTIONAL :: start
1541LOGICAL,INTENT(in),OPTIONAL :: full_steps
1542
1543INTEGER :: i, j, k, l, na, nf, n
1544INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
1545INTEGER(kind=int_ll) :: stepms, mstepms
1546LOGICAL :: lforecast
1547TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1548TYPE(arrayof_datetime) :: a_otime
1549TYPE(arrayof_vol7d_timerange) :: a_otimerange
1550TYPE(arrayof_integer) :: a_dtratio
1551LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
1552TYPE(ttr_mapper) :: lmapper
1553CHARACTER(len=8) :: env_var
1554LOGICAL :: climat_behavior
1555
1556
1557! enable bad behavior for climat database (only for instantaneous data)
1558env_var = ''
1559CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
1560climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
1561
1562! compute length of cumulation step in seconds
1563CALL getval(timedelta_depop(step), asec=steps)
1564
1565! create a mask of suitable timeranges
1566ALLOCATE(mask_timerange(SIZE(itimerange)))
1567mask_timerange(:) = itimerange(:)%timerange == tri &
1568 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
1569
1570IF (PRESENT(dtratio)) THEN
1571 WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0)
1572 mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0
1573 ELSEWHERE
1574 mask_timerange(:) = .false.
1575 END WHERE
1576ELSE ! instantaneous
1577 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
1578ENDIF
1579
1580#ifdef DEBUG
1581CALL l4f_log(l4f_debug, &
1582 '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
1583 t2c(count(mask_timerange)))
1584#endif
1585
1586! euristically determine whether we are dealing with an
1587! analysis/observation or a forecast dataset
1588na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
1589nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
1590lforecast = nf >= na
1591#ifdef DEBUG
1592CALL l4f_log(l4f_debug, &
1593 'recompute_stat_proc_agg, na: '//t2c(na)//', nf: '//t2c(nf))
1594#endif
1595
1596! keep only timeranges of one type (really necessary?)
1597IF (lforecast) THEN
1598 CALL l4f_log(l4f_info, &
1599 'recompute_stat_proc_agg, processing in forecast mode')
1600ELSE
1601 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
1602 CALL l4f_log(l4f_info, &
1603 'recompute_stat_proc_agg, processing in analysis mode')
1604ENDIF
1605
1606#ifdef DEBUG
1607CALL l4f_log(l4f_debug, &
1608 '(re)compute_stat_proc_agg, number of useful timeranges: '// &
1609 t2c(count(mask_timerange)))
1610#endif
1611
1612IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1613 ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
1614 IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
1615 RETURN
1616ENDIF
1617
1618! determine start and end of processing period, should work with p2==0
1619lstart = datetime_miss
1620IF (PRESENT(start)) lstart = start
1621lend = itime(SIZE(itime))
1622! compute some quantities
1623maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
1624maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
1625minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
1626IF (time_definition == 0) THEN ! reference time
1627 lend = lend + timedelta_new(sec=maxp1)
1628ENDIF
1629! extend interval at the end in order to include all the data in case
1630! frac_valid<1; must use < and not <= in the DO WHILE loops some lines
1631! below in order to exclude the last full interval which would be empty
1632lend = lend + step
1633IF (lstart == datetime_miss) THEN ! autodetect
1634 lstart = itime(1)
1635! if autodetected, adjust to obtain real absolute start of data
1636 IF (time_definition == 0) THEN ! reference time
1637 lstart = lstart + timedelta_new(sec=minp1mp2)
1638 ELSE ! verification time
1639! go back to start of longest processing interval
1640 lstart = lstart - timedelta_new(sec=maxp2)
1641 ENDIF
1642! apply full_steps in analysis mode and when start is not specified
1643! (start by itself in analysis mode implies full_steps with respect to
1644! start instead of absolute full steps), todo in forecast mode
1645 IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN
1646 lstart = lstart - (mod(lstart, step)) ! round to step, (should be MODULO, not MOD)
1647 ENDIF
1648ENDIF
1649
1650#ifdef DEBUG
1651CALL l4f_log(l4f_debug, &
1652 'recompute_stat_proc_agg, processing period: '//t2c(lstart)//' - '//t2c(lend))
1653#endif
1654
1655! create output time and timerange lists
1656
1657IF (lforecast) THEN ! forecast mode
1658 IF (time_definition == 0) THEN ! reference time
1659 CALL insert(a_otime, itime) ! should I limit to elements itime >= lstart?
1660
1661! apply start shift to timerange, not to reftime
1662! why did we use itime(SIZE(itime)) (last ref time)?
1663! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
1664 CALL getval(lstart-itime(1), asec=dstart)
1665! allow starting before first reftime but restrict dtstart to -steps
1666! dstart = MAX(0, dstart)
1667 IF (dstart < 0) dstart = mod(dstart, steps)
1668 DO p1 = steps + dstart, maxp1, steps
1669 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1670 ENDDO
1671
1672 ELSE ! verification time
1673
1674! verification time in forecast mode is the ugliest case, because we
1675! don't know a priori how many different (thus incompatible) reference
1676! times we have, so some assumption of regularity has to be made. For
1677! this reason msteps, the minimum step between two times, is
1678! computed. We choose to compute it as a difference between itime
1679! elements, it could be also computed as difference of itimerange%p1
1680! elements. But what if it is not modulo steps?
1681 mstepms = steps*1000_int_ll
1682 DO i = 2, SIZE(itime)
1683 CALL getval(itime(i)-itime(i-1), amsec=stepms)
1684 IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN
1685 msteps = stepms/1000_int_ll
1686 IF (mod(steps, msteps) == 0) mstepms = stepms
1687 ENDIF
1688 ENDDO
1689 msteps = mstepms/1000_int_ll
1690
1691 tmptime = lstart + step
1692 DO WHILE(tmptime < lend) ! < since lend has been extended
1693 CALL insert_unique(a_otime, tmptime)
1694 tmptime = tmptime + step
1695 ENDDO
1696
1697! in next loop, we used initially steps instead of msteps (see comment
1698! above), this gave much less combinations of time/timeranges with
1699! possible empty output; we start from msteps instead of steps in
1700! order to include partial initial intervals in case frac_valid<1;
1701! however, as a gemeral rule, for aggregation of forecasts it's better
1702! to use reference time
1703 DO p1 = msteps, maxp1, msteps
1704 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1705 ENDDO
1706
1707 ENDIF
1708
1709ELSE ! analysis mode
1710 tmptime = lstart + step
1711 DO WHILE(tmptime < lend) ! < since lend has been extended
1712 CALL insert_unique(a_otime, tmptime)
1713 tmptime = tmptime + step
1714 ENDDO
1715 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
1716
1717ENDIF
1718
1719CALL packarray(a_otime)
1720CALL packarray(a_otimerange)
1721otime => a_otime%array
1722otimerange => a_otimerange%array
1723CALL sort(otime)
1724CALL sort(otimerange)
1725! delete local objects keeping the contents
1726CALL delete(a_otime, nodealloc=.true.)
1727CALL delete(a_otimerange, nodealloc=.true.)
1728
1729#ifdef DEBUG
1730CALL l4f_log(l4f_debug, &
1731 'recompute_stat_proc_agg, output time and timerange: '//&
1732 t2c(SIZE(otime))//', '//t2c(size(otimerange)))
1733#endif
1734
1735IF (PRESENT(dtratio)) THEN
1736! count the possible i/o interval ratios
1737 DO k = 1, SIZE(itimerange)
1738 IF (itimerange(k)%p2 /= 0) &
1739 CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
1740 ENDDO
1741 CALL packarray(a_dtratio)
1742 dtratio => a_dtratio%array
1743 CALL sort(dtratio)
1744! delete local object keeping the contents
1745 CALL delete(a_dtratio, nodealloc=.true.)
1746
1747#ifdef DEBUG
1748 CALL l4f_log(l4f_debug, &
1749 'recompute_stat_proc_agg, found '//t2c(size(dtratio))// &
1750 ' possible aggregation ratios, from '// &
1751 t2c(dtratio(1))//' to '//t2c(dtratio(SIZE(dtratio))))
1752#endif
1753
1754 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1755 do_itimerange1: DO l = 1, SIZE(itimerange)
1756 IF (.NOT.mask_timerange(l)) cycle do_itimerange1
1757 do_itime1: DO k = 1, SIZE(itime)
1758 CALL time_timerange_get_period(itime(k), itimerange(l), &
1759 time_definition, pstart1, pend1, reftime1)
1760 do_otimerange1: DO j = 1, SIZE(otimerange)
1761 do_otime1: DO i = 1, SIZE(otime)
1762 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1763 time_definition, pstart2, pend2, reftime2)
1764 IF (lforecast) THEN
1765 IF (reftime1 /= reftime2) cycle do_otime1
1766 ENDIF
1767
1768 IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
1769 mod(pstart1-pstart2, pend1-pstart1) == timedelta_0) THEN ! useful
1770 lmapper%it = k
1771 lmapper%itr = l
1772 lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
1773 n = append(map_ttr(i,j), lmapper)
1774 cycle do_itime1 ! can contribute only to a single interval
1775 ENDIF
1776 ENDDO do_otime1
1777 ENDDO do_otimerange1
1778 ENDDO do_itime1
1779 ENDDO do_itimerange1
1780
1781ELSE
1782
1783 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1784 do_itimerange2: DO l = 1, SIZE(itimerange)
1785 IF (.NOT.mask_timerange(l)) cycle do_itimerange2
1786 do_itime2: DO k = 1, SIZE(itime)
1787 CALL time_timerange_get_period(itime(k), itimerange(l), &
1788 time_definition, pstart1, pend1, reftime1)
1789 do_otimerange2: DO j = 1, SIZE(otimerange)
1790 do_otime2: DO i = 1, SIZE(otime)
1791 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1792 time_definition, pstart2, pend2, reftime2)
1793 IF (lforecast) THEN
1794 IF (reftime1 /= reftime2) cycle do_otime2
1795 ENDIF
1796
1797 IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
1798 IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
1799 lmapper%it = k
1800 lmapper%itr = l
1801 IF (pstart1 == pstart2) THEN
1802 lmapper%extra_info = 1 ! start of interval
1803 ELSE IF (pend1 == pend2) THEN
1804 lmapper%extra_info = 2 ! end of interval
1805 ELSE
1806 lmapper%extra_info = imiss
1807 ENDIF
1808 lmapper%time = pstart1 ! = pend1, order by time?
1809 n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
1810! no CYCLE, a single input can contribute to multiple output intervals
1811 ENDIF
1812 ENDDO do_otime2
1813 ENDDO do_otimerange2
1814 ENDDO do_itime2
1815 ENDDO do_itimerange2
1816
1817ENDIF
1818
1819END SUBROUTINE recompute_stat_proc_agg_common
1820
1821
1822SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
1823 max_step, weights)
1824TYPE(datetime),INTENT(in) :: vertime(:)
1825TYPE(datetime),INTENT(in) :: pstart
1826TYPE(datetime),INTENT(in) :: pend
1827LOGICAL,INTENT(in) :: time_mask(:)
1828TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
1829DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
1830
1831INTEGER :: i, nt
1832TYPE(datetime),ALLOCATABLE :: lvertime(:)
1833TYPE(datetime) :: half, nexthalf
1834INTEGER(kind=int_ll) :: dt, tdt
1835
1836nt = count(time_mask)
1837ALLOCATE(lvertime(nt))
1838lvertime = pack(vertime, mask=time_mask)
1839
1840IF (PRESENT(max_step)) THEN
1841! new way
1842! max_step = timedelta_0
1843! DO i = 1, nt - 1
1844! IF (lvertime(i+1) - lvertime(i) > max_step) &
1845! max_step = lvertime(i+1) - lvertime(i)
1846! ENDDO
1847! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
1848! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
1849! old way
1850 IF (nt == 1) THEN
1851 max_step = pend - pstart
1852 ELSE
1853 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
1854 max_step = half - pstart
1855 DO i = 2, nt - 1
1856 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
1857 IF (nexthalf - half > max_step) max_step = nexthalf - half
1858 half = nexthalf
1859 ENDDO
1860 IF (pend - half > max_step) max_step = pend - half
1861 ENDIF
1862
1863ENDIF
1864
1865IF (PRESENT(weights)) THEN
1866 IF (nt == 1) THEN
1867 weights(1) = 1.0
1868 ELSE
1869 CALL getval(pend - pstart, amsec=tdt)
1870 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
1871 CALL getval(half - pstart, amsec=dt)
1872 weights(1) = dble(dt)/dble(tdt)
1873 DO i = 2, nt - 1
1874 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
1875 CALL getval(nexthalf - half, amsec=dt)
1876 weights(i) = dble(dt)/dble(tdt)
1877 half = nexthalf
1878 ENDDO
1879 CALL getval(pend - half, amsec=dt)
1880 weights(nt) = dble(dt)/dble(tdt)
1881 ENDIF
1882ENDIF
1883
1884END SUBROUTINE compute_stat_proc_agg_sw
1885
1886! get start of period, end of period and reference time from time,
1887! timerange, according to time_definition.
1888SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
1889 pstart, pend, reftime)
1890TYPE(datetime),INTENT(in) :: time
1891TYPE(vol7d_timerange),INTENT(in) :: timerange
1892INTEGER,INTENT(in) :: time_definition
1893TYPE(datetime),INTENT(out) :: reftime
1894TYPE(datetime),INTENT(out) :: pstart
1895TYPE(datetime),INTENT(out) :: pend
1896
1897TYPE(timedelta) :: p1, p2
1898
1899
1900p1 = timedelta_new(sec=timerange%p1) ! end of period
1901p2 = timedelta_new(sec=timerange%p2) ! length of period
1902
1903IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
1904! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
1905 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
1906 pstart = datetime_miss
1907 pend = datetime_miss
1908 reftime = datetime_miss
1909 RETURN
1910ENDIF
1911
1912IF (time_definition == 0) THEN ! time == reference time
1913 reftime = time
1914 pend = time + p1
1915 pstart = pend - p2
1916ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time
1917 pend = time
1918 pstart = time - p2
1919 reftime = time - p1
1920ELSE
1921 pstart = datetime_miss
1922 pend = datetime_miss
1923 reftime = datetime_miss
1924ENDIF
1925
1926END SUBROUTINE time_timerange_get_period
1927
1928
1929! get start of period, end of period and reference time from time,
1930! timerange, according to time_definition. step is taken separately
1931! from timerange and can be "popular"
1932SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
1933 pstart, pend, reftime)
1934TYPE(datetime),INTENT(in) :: time
1935TYPE(vol7d_timerange),INTENT(in) :: timerange
1936TYPE(timedelta),INTENT(in) :: step
1937INTEGER,INTENT(in) :: time_definition
1938TYPE(datetime),INTENT(out) :: reftime
1939TYPE(datetime),INTENT(out) :: pstart
1940TYPE(datetime),INTENT(out) :: pend
1941
1942TYPE(timedelta) :: p1
1943
1944
1945p1 = timedelta_new(sec=timerange%p1) ! end of period
1946
1947IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
1948! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
1949 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
1950 pstart = datetime_miss
1951 pend = datetime_miss
1952 reftime = datetime_miss
1953 RETURN
1954ENDIF
1955
1956IF (time_definition == 0) THEN ! time == reference time
1957 reftime = time
1958 pend = time + p1
1959 pstart = pend - step
1960ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time
1961 pend = time
1962 pstart = time - step
1963 reftime = time - p1
1964ELSE
1965 pstart = datetime_miss
1966 pend = datetime_miss
1967 reftime = datetime_miss
1968ENDIF
1969
1970END SUBROUTINE time_timerange_get_period_pop
1971
1972
1973! set time, timerange%p1, timerange%p2 according to pstart, pend,
1974! reftime and time_definition.
1975SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
1976 pstart, pend, reftime)
1977TYPE(datetime),INTENT(out) :: time
1978TYPE(vol7d_timerange),INTENT(inout) :: timerange
1979INTEGER,INTENT(in) :: time_definition
1980TYPE(datetime),INTENT(in) :: reftime
1981TYPE(datetime),INTENT(in) :: pstart
1982TYPE(datetime),INTENT(in) :: pend
1983
1984TYPE(timedelta) :: p1, p2
1985INTEGER(kind=int_ll) :: dmsec
1986
1987
1988IF (time_definition == 0) THEN ! time == reference time
1989 time = reftime
1990 p1 = pend - reftime
1991 p2 = pend - pstart
1992ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time
1993 time = pend
1994 p1 = pend - reftime
1995 p2 = pend - pstart
1996ELSE
1997 time = datetime_miss
1998ENDIF
1999
2000IF (time /= datetime_miss) THEN
2001 CALL getval(p1, amsec=dmsec) ! end of period
2002 timerange%p1 = int(dmsec/1000_int_ll)
2003 CALL getval(p2, amsec=dmsec) ! length of period
2004 timerange%p2 = int(dmsec/1000_int_ll)
2005ELSE
2006 timerange%p1 = imiss
2007 timerange%p2 = imiss
2008ENDIF
2009
2010END SUBROUTINE time_timerange_set_period
2011
2012
2013END MODULE stat_proc_engine
Restituiscono il valore dell'oggetto nella forma desiderata.
Operatore di resto della divisione.
Functions that return a trimmed CHARACTER representation of the input variable.
Quick method to append an element to the array.
Destructor for finalizing an array object.
Method for inserting elements of the array at a desired position.
Method for packing the array object reducing at a minimum the memory occupation, without destroying i...
Classi per la gestione delle coordinate temporali.
This module contains functions that are only for internal use of the library.
Classe per la gestione di un volume completo di dati osservati.

Generated with Doxygen.