libsim Versione 7.2.4
|
◆ map_inv_distinct_ttr_mapper()
map inv distinct Definizione alla linea 1210 del file stat_proc_engine.F90. 1212! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1213! authors:
1214! Davide Cesari <dcesari@arpa.emr.it>
1215! Paolo Patruno <ppatruno@arpa.emr.it>
1216
1217! This program is free software; you can redistribute it and/or
1218! modify it under the terms of the GNU General Public License as
1219! published by the Free Software Foundation; either version 2 of
1220! the License, or (at your option) any later version.
1221
1222! This program is distributed in the hope that it will be useful,
1223! but WITHOUT ANY WARRANTY; without even the implied warranty of
1224! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1225! GNU General Public License for more details.
1226
1227! You should have received a copy of the GNU General Public License
1228! along with this program. If not, see <http://www.gnu.org/licenses/>.
1229#include "config.h"
1230
1237IMPLICIT NONE
1238
1239! per ogni elemento i,j dell'output, definire n elementi di input ad
1240! esso contribuenti (arrayof_ttr_mapper) con le seguenti informazioni
1241TYPE ttr_mapper
1242 INTEGER :: it=imiss ! k
1243 INTEGER :: itr=imiss ! l
1244 INTEGER :: extra_info=imiss ! dtratio for intervals, extreme for point in time
1245 TYPE(datetime) :: time=datetime_miss ! time for point in time
1246END TYPE ttr_mapper
1247
1248INTERFACE OPERATOR (==)
1249 MODULE PROCEDURE ttr_mapper_eq
1250END INTERFACE
1251
1252INTERFACE OPERATOR (/=)
1253 MODULE PROCEDURE ttr_mapper_ne
1254END INTERFACE
1255
1256INTERFACE OPERATOR (>)
1257 MODULE PROCEDURE ttr_mapper_gt
1258END INTERFACE
1259
1260INTERFACE OPERATOR (<)
1261 MODULE PROCEDURE ttr_mapper_lt
1262END INTERFACE
1263
1264INTERFACE OPERATOR (>=)
1265 MODULE PROCEDURE ttr_mapper_ge
1266END INTERFACE
1267
1268INTERFACE OPERATOR (<=)
1269 MODULE PROCEDURE ttr_mapper_le
1270END INTERFACE
1271
1272#undef VOL7D_POLY_TYPE
1273#undef VOL7D_POLY_TYPES
1274#undef ENABLE_SORT
1275#define VOL7D_POLY_TYPE TYPE(ttr_mapper)
1276#define VOL7D_POLY_TYPES _ttr_mapper
1277#define ENABLE_SORT
1278#include "array_utilities_pre.F90"
1279
1280#define ARRAYOF_ORIGTYPE TYPE(ttr_mapper)
1281#define ARRAYOF_TYPE arrayof_ttr_mapper
1282#define ARRAYOF_ORIGEQ 1
1283#define ARRAYOF_ORIGGT 1
1284#include "arrayof_pre.F90"
1285! from arrayof
1286
1287
1288CONTAINS
1289
1290! simplified comparisons only on time
1291ELEMENTAL FUNCTION ttr_mapper_eq(this, that) RESULT(res)
1292TYPE(ttr_mapper),INTENT(IN) :: this, that
1293LOGICAL :: res
1294
1295res = this%time == that%time
1296
1297END FUNCTION ttr_mapper_eq
1298
1299ELEMENTAL FUNCTION ttr_mapper_ne(this, that) RESULT(res)
1300TYPE(ttr_mapper),INTENT(IN) :: this, that
1301LOGICAL :: res
1302
1303res = this%time /= that%time
1304
1305END FUNCTION ttr_mapper_ne
1306
1307ELEMENTAL FUNCTION ttr_mapper_gt(this, that) RESULT(res)
1308TYPE(ttr_mapper),INTENT(IN) :: this, that
1309LOGICAL :: res
1310
1311res = this%time > that%time
1312
1313END FUNCTION ttr_mapper_gt
1314
1315ELEMENTAL FUNCTION ttr_mapper_lt(this, that) RESULT(res)
1316TYPE(ttr_mapper),INTENT(IN) :: this, that
1317LOGICAL :: res
1318
1319res = this%time < that%time
1320
1321END FUNCTION ttr_mapper_lt
1322
1323ELEMENTAL FUNCTION ttr_mapper_ge(this, that) RESULT(res)
1324TYPE(ttr_mapper),INTENT(IN) :: this, that
1325LOGICAL :: res
1326
1327res = this%time >= that%time
1328
1329END FUNCTION ttr_mapper_ge
1330
1331ELEMENTAL FUNCTION ttr_mapper_le(this, that) RESULT(res)
1332TYPE(ttr_mapper),INTENT(IN) :: this, that
1333LOGICAL :: res
1334
1335res = this%time <= that%time
1336
1337END FUNCTION ttr_mapper_le
1338
1339#include "arrayof_post.F90"
1340#include "array_utilities_inc.F90"
1341
1342
1343! common operations for statistical processing by differences
1344SUBROUTINE recompute_stat_proc_diff_common(itime, itimerange, stat_proc, step, &
1345 otime, otimerange, map_tr, f, keep_tr, time_definition, full_steps, &
1346 start)
1347TYPE(datetime),INTENT(in) :: itime(:)
1348TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1349INTEGER,INTENT(in) :: stat_proc
1350TYPE(timedelta),INTENT(in) :: step
1351TYPE(datetime),POINTER :: otime(:)
1352TYPE(vol7d_timerange),POINTER :: otimerange(:)
1353INTEGER,ALLOCATABLE,INTENT(out) :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
1354INTEGER :: nitr
1355LOGICAL,ALLOCATABLE :: mask_timerange(:)
1356INTEGER,INTENT(in) :: time_definition
1357LOGICAL,INTENT(in),OPTIONAL :: full_steps
1358TYPE(datetime),INTENT(in),OPTIONAL :: start
1359
1360INTEGER :: i, j, k, l, dirtyrep
1361INTEGER :: steps
1362LOGICAL :: lfull_steps, useful
1363TYPE(datetime) :: lstart, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1364TYPE(vol7d_timerange) :: tmptimerange
1365TYPE(arrayof_datetime) :: a_otime
1366TYPE(arrayof_vol7d_timerange) :: a_otimerange
1367TYPE(timedelta) :: start_delta
1368
1369! compute length of cumulation step in seconds
1371
1372lstart = datetime_miss
1373IF (PRESENT(start)) lstart = start
1374lfull_steps = optio_log(full_steps)
1375
1376! create a mask of suitable timeranges
1377ALLOCATE(mask_timerange(SIZE(itimerange)))
1378mask_timerange(:) = itimerange(:)%timerange == stat_proc &
1379 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
1380 .AND. itimerange(:)%p1 >= 0 &
1381 .AND. itimerange(:)%p2 > 0
1382
1383IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges
1384 mask_timerange(:) = mask_timerange(:) .AND. &
1385 (itimerange(:)%p1 == 0 .OR. & ! all analyses
1388ENDIF
1389! mask_timerange includes all candidate timeranges
1390
1391nitr = count(mask_timerange)
1392ALLOCATE(f(nitr))
1393j = 1
1394DO i = 1, nitr
1395 DO WHILE(.NOT.mask_timerange(j))
1396 j = j + 1
1397 ENDDO
1398 f(i) = j
1399 j = j + 1
1400ENDDO
1401
1402! now we have to evaluate time/timerage pairs which do not need processing
1403ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
1404CALL compute_keep_tr()
1405
1406ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
1407map_tr(:,:,:,:,:) = imiss
1408
1409! scan through all possible combinations of time and timerange
1410DO dirtyrep = 1, 2
1411 IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
1416 ENDIF
1417 DO l = 1, SIZE(itime)
1418 DO k = 1, nitr
1419 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1420 time_definition, pstart2, pend2, reftime2)
1421
1422 DO j = 1, SIZE(itime)
1423 DO i = 1, nitr
1424 useful = .false.
1425 CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
1426 time_definition, pstart1, pend1, reftime1)
1427 tmptimerange = vol7d_timerange_new(timerange=stat_proc)
1428
1429 IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
1430 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
1431 CALL time_timerange_set_period(tmptime, tmptimerange, &
1432 time_definition, pend1, pend2, reftime2)
1433 IF (lfull_steps) THEN
1435 useful = .true.
1436 ENDIF
1437 ELSE
1438 useful = .true.
1439 ENDIF
1440
1441 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
1442 CALL time_timerange_set_period(tmptime, tmptimerange, &
1443 time_definition, pstart2, pstart1, pstart1)
1444 IF (lfull_steps) THEN
1446 useful = .true.
1447 ENDIF
1448 ELSE
1449 useful = .true.
1450 ENDIF
1451 ENDIF
1452
1453 ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
1454 IF (lfull_steps) THEN
1456! lstart shifts the interval for computing modulo step, this does not
1457! remove data before lstart but just shifts the phase
1458 start_delta = lstart-reftime2
1459 ELSE
1460 start_delta = timedelta_0
1461 ENDIF
1462 ENDIF
1463
1464 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
1465 CALL time_timerange_set_period(tmptime, tmptimerange, &
1466 time_definition, pend1, pend2, reftime2)
1467 IF (lfull_steps) THEN
1469 useful = .true.
1470 ENDIF
1471 ELSE
1472 useful = .true.
1473 ENDIF
1474! keep only data after lstart
1476 IF (lstart > pend1) useful = .false.
1477 ENDIF
1478
1479 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! |-=
1480 CALL time_timerange_set_period(tmptime, tmptimerange, &
1481 time_definition, pstart2, pstart1, reftime2)
1482 IF (lfull_steps) THEN
1484 useful = .true.
1485 ENDIF
1486 ELSE
1487 useful = .true.
1488 ENDIF
1489! keep only data after lstart
1491 IF (lstart > pstart2) useful = .false.
1492 ENDIF
1493
1494 ENDIF
1495 ENDIF
1496 useful = useful .AND. tmptime /= datetime_miss .AND. &
1497 tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
1498
1499 IF (useful) THEN ! add a_otime, a_otimerange
1500 map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
1501 map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
1502 ENDIF
1503 ENDDO
1504 ENDDO
1505 ENDDO
1506 ENDDO
1507ENDDO
1508
1509! we have to repeat the computation with sorted arrays
1510CALL compute_keep_tr()
1511
1512otime => a_otime%array
1513otimerange => a_otimerange%array
1514! delete local objects keeping the contents
1517
1518#ifdef DEBUG
1519CALL l4f_log(l4f_debug, &
1524CALL l4f_log(l4f_debug, &
1527CALL l4f_log(l4f_debug, &
1529CALL l4f_log(l4f_debug, &
1531CALL l4f_log(l4f_debug, &
1533CALL l4f_log(l4f_debug, &
1535#endif
1536
1537CONTAINS
1538
1539SUBROUTINE compute_keep_tr()
1540INTEGER :: start_deltas
1541
1542keep_tr(:,:,:) = imiss
1543DO l = 1, SIZE(itime)
1544 itrloop: DO k = 1, nitr
1545 IF (itimerange(f(k))%p2 == steps) THEN
1546 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1547 time_definition, pstart2, pend2, reftime2)
1548 useful = .false.
1549! keep only data after lstart
1551 IF (lstart > pstart2) cycle itrloop
1552 ENDIF
1553 IF (reftime2 == pend2) THEN ! analysis
1556 useful = .true.
1557 ENDIF
1558 ELSE IF (lfull_steps) THEN
1560 useful = .true.
1561 ENDIF
1562 ELSE
1563 useful = .true.
1564 ENDIF
1565 ELSE ! forecast
1566 IF (lfull_steps) THEN
1567! same as above for start_delta, but in seconds and not in timerange/timedelta units
1569 start_deltas = timedelta_getamsec(lstart-reftime2)/1000_int_ll
1570 ELSE
1571 start_deltas = 0
1572 ENDIF
1574 useful = .true.
1575 ENDIF
1576 ELSE
1577 useful = .true.
1578 ENDIF
1579 ENDIF
1580 IF (useful) THEN
1581! CALL time_timerange_set_period(tmptime, tmptimerange, &
1582! time_definition, pstart2, pend2, reftime2)
1583 keep_tr(k,l,1) = append_unique(a_otime, itime(l))
1584 keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k)))
1585 ENDIF
1586 ENDIF
1587 ENDDO itrloop
1588ENDDO
1589
1590END SUBROUTINE compute_keep_tr
1591
1592END SUBROUTINE recompute_stat_proc_diff_common
1593
1594
1595! common operations for statistical processing by metamorphosis
1596SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
1597 otimerange, map_tr)
1598INTEGER,INTENT(in) :: istat_proc
1599TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1600INTEGER,INTENT(in) :: ostat_proc
1601TYPE(vol7d_timerange),POINTER :: otimerange(:)
1602INTEGER,POINTER :: map_tr(:)
1603
1604INTEGER :: i
1605LOGICAL :: tr_mask(SIZE(itimerange))
1606
1607IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1608 ALLOCATE(otimerange(0), map_tr(0))
1609 RETURN
1610ENDIF
1611
1612! useful timeranges
1613tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
1614 .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
1615ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
1616
1617otimerange = pack(itimerange, mask=tr_mask)
1618otimerange(:)%timerange = ostat_proc
1619map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
1620
1621END SUBROUTINE compute_stat_proc_metamorph_common
1622
1623
1624! common operations for statistical processing by aggregation
1625SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
1626 step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
1627TYPE(datetime),INTENT(in) :: itime(:)
1628TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1629INTEGER,INTENT(in) :: stat_proc
1630INTEGER,INTENT(in) :: tri
1631TYPE(timedelta),INTENT(in) :: step
1632INTEGER,INTENT(in) :: time_definition
1633TYPE(datetime),POINTER :: otime(:)
1634TYPE(vol7d_timerange),POINTER :: otimerange(:)
1635TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
1636INTEGER,POINTER,OPTIONAL :: dtratio(:)
1637TYPE(datetime),INTENT(in),OPTIONAL :: start
1638LOGICAL,INTENT(in),OPTIONAL :: full_steps
1639
1640INTEGER :: i, j, k, l, na, nf, n
1641INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
1642INTEGER(kind=int_ll) :: stepms, mstepms
1643LOGICAL :: lforecast
1644TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1645TYPE(arrayof_datetime) :: a_otime
1646TYPE(arrayof_vol7d_timerange) :: a_otimerange
1647TYPE(arrayof_integer) :: a_dtratio
1648LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
1649TYPE(ttr_mapper) :: lmapper
1650CHARACTER(len=8) :: env_var
1651LOGICAL :: climat_behavior
1652
1653
1654! enable bad behavior for climat database (only for instantaneous data)
1655env_var = ''
1656CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
1657climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
1658
1659! compute length of cumulation step in seconds
1661
1662! create a mask of suitable timeranges
1663ALLOCATE(mask_timerange(SIZE(itimerange)))
1664mask_timerange(:) = itimerange(:)%timerange == tri &
1665 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
1666
1667IF (PRESENT(dtratio)) THEN
1668 WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0)
1669 mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0
1670 ELSEWHERE
1671 mask_timerange(:) = .false.
1672 END WHERE
1673ELSE ! instantaneous
1674 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
1675ENDIF
1676
1677#ifdef DEBUG
1678CALL l4f_log(l4f_debug, &
1679 '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
1680 t2c(count(mask_timerange)))
1681#endif
1682
1683! euristically determine whether we are dealing with an
1684! analysis/observation or a forecast dataset
1685na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
1686nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
1687lforecast = nf >= na
1688#ifdef DEBUG
1689CALL l4f_log(l4f_debug, &
1691#endif
1692
1693! keep only timeranges of one type (really necessary?)
1694IF (lforecast) THEN
1695 CALL l4f_log(l4f_info, &
1696 'recompute_stat_proc_agg, processing in forecast mode')
1697ELSE
1698 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
1699 CALL l4f_log(l4f_info, &
1700 'recompute_stat_proc_agg, processing in analysis mode')
1701ENDIF
1702
1703#ifdef DEBUG
1704CALL l4f_log(l4f_debug, &
1705 '(re)compute_stat_proc_agg, number of useful timeranges: '// &
1706 t2c(count(mask_timerange)))
1707#endif
1708
1709IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1710 ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
1711 IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
1712 RETURN
1713ENDIF
1714
1715! determine start and end of processing period, should work with p2==0
1716lstart = datetime_miss
1717IF (PRESENT(start)) lstart = start
1718lend = itime(SIZE(itime))
1719! compute some quantities
1720maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
1721maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
1722minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
1723IF (time_definition == 0) THEN ! reference time
1724 lend = lend + timedelta_new(sec=maxp1)
1725ENDIF
1726! extend interval at the end in order to include all the data in case
1727! frac_valid<1; must use < and not <= in the DO WHILE loops some lines
1728! below in order to exclude the last full interval which would be empty
1729lend = lend + step
1730IF (lstart == datetime_miss) THEN ! autodetect
1731 lstart = itime(1)
1732! if autodetected, adjust to obtain real absolute start of data
1733 IF (time_definition == 0) THEN ! reference time
1734 lstart = lstart + timedelta_new(sec=minp1mp2)
1735 ELSE ! verification time
1736! go back to start of longest processing interval
1737 lstart = lstart - timedelta_new(sec=maxp2)
1738 ENDIF
1739! apply full_steps in analysis mode and when start is not specified
1740! (start by itself in analysis mode implies full_steps with respect to
1741! start instead of absolute full steps), todo in forecast mode
1742 IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN
1744 ENDIF
1745ENDIF
1746
1747#ifdef DEBUG
1748CALL l4f_log(l4f_debug, &
1750#endif
1751
1752! create output time and timerange lists
1753
1754IF (lforecast) THEN ! forecast mode
1755 IF (time_definition == 0) THEN ! reference time
1757
1758! apply start shift to timerange, not to reftime
1759! why did we use itime(SIZE(itime)) (last ref time)?
1760! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
1762! allow starting before first reftime but restrict dtstart to -steps
1763! dstart = MAX(0, dstart)
1765 DO p1 = steps + dstart, maxp1, steps
1766 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1767 ENDDO
1768
1769 ELSE ! verification time
1770
1771! verification time in forecast mode is the ugliest case, because we
1772! don't know a priori how many different (thus incompatible) reference
1773! times we have, so some assumption of regularity has to be made. For
1774! this reason msteps, the minimum step between two times, is
1775! computed. We choose to compute it as a difference between itime
1776! elements, it could be also computed as difference of itimerange%p1
1777! elements. But what if it is not modulo steps?
1778 mstepms = steps*1000_int_ll
1779 DO i = 2, SIZE(itime)
1781 IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN
1782 msteps = stepms/1000_int_ll
1784 ENDIF
1785 ENDDO
1786 msteps = mstepms/1000_int_ll
1787
1788 tmptime = lstart + step
1789 DO WHILE(tmptime < lend) ! < since lend has been extended
1790 CALL insert_unique(a_otime, tmptime)
1791 tmptime = tmptime + step
1792 ENDDO
1793
1794! in next loop, we used initially steps instead of msteps (see comment
1795! above), this gave much less combinations of time/timeranges with
1796! possible empty output; we start from msteps instead of steps in
1797! order to include partial initial intervals in case frac_valid<1;
1798! however, as a gemeral rule, for aggregation of forecasts it's better
1799! to use reference time
1800 DO p1 = msteps, maxp1, msteps
1801 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1802 ENDDO
1803
1804 ENDIF
1805
1806ELSE ! analysis mode
1807 tmptime = lstart + step
1808 DO WHILE(tmptime < lend) ! < since lend has been extended
1809 CALL insert_unique(a_otime, tmptime)
1810 tmptime = tmptime + step
1811 ENDDO
1812 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
1813
1814ENDIF
1815
1818otime => a_otime%array
1819otimerange => a_otimerange%array
1822! delete local objects keeping the contents
1825
1826#ifdef DEBUG
1827CALL l4f_log(l4f_debug, &
1828 'recompute_stat_proc_agg, output time and timerange: '//&
1830#endif
1831
1832IF (PRESENT(dtratio)) THEN
1833! count the possible i/o interval ratios
1834 DO k = 1, SIZE(itimerange)
1835 IF (itimerange(k)%p2 /= 0) &
1836 CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
1837 ENDDO
1839 dtratio => a_dtratio%array
1841! delete local object keeping the contents
1843
1844#ifdef DEBUG
1845 CALL l4f_log(l4f_debug, &
1847 ' possible aggregation ratios, from '// &
1849#endif
1850
1851 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1852 do_itimerange1: DO l = 1, SIZE(itimerange)
1853 IF (.NOT.mask_timerange(l)) cycle do_itimerange1
1854 do_itime1: DO k = 1, SIZE(itime)
1855 CALL time_timerange_get_period(itime(k), itimerange(l), &
1856 time_definition, pstart1, pend1, reftime1)
1857 do_otimerange1: DO j = 1, SIZE(otimerange)
1858 do_otime1: DO i = 1, SIZE(otime)
1859 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1860 time_definition, pstart2, pend2, reftime2)
1861 IF (lforecast) THEN
1862 IF (reftime1 /= reftime2) cycle do_otime1
1863 ENDIF
1864
1865 IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
1867 lmapper%it = k
1868 lmapper%itr = l
1869 lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
1870 n = append(map_ttr(i,j), lmapper)
1871 cycle do_itime1 ! can contribute only to a single interval
1872 ENDIF
1873 ENDDO do_otime1
1874 ENDDO do_otimerange1
1875 ENDDO do_itime1
1876 ENDDO do_itimerange1
1877
1878ELSE
1879
1880 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1881 do_itimerange2: DO l = 1, SIZE(itimerange)
1882 IF (.NOT.mask_timerange(l)) cycle do_itimerange2
1883 do_itime2: DO k = 1, SIZE(itime)
1884 CALL time_timerange_get_period(itime(k), itimerange(l), &
1885 time_definition, pstart1, pend1, reftime1)
1886 do_otimerange2: DO j = 1, SIZE(otimerange)
1887 do_otime2: DO i = 1, SIZE(otime)
1888 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1889 time_definition, pstart2, pend2, reftime2)
1890 IF (lforecast) THEN
1891 IF (reftime1 /= reftime2) cycle do_otime2
1892 ENDIF
1893
1894 IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
1895 IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
1896 lmapper%it = k
1897 lmapper%itr = l
1898 IF (pstart1 == pstart2) THEN
1899 lmapper%extra_info = 1 ! start of interval
1900 ELSE IF (pend1 == pend2) THEN
1901 lmapper%extra_info = 2 ! end of interval
1902 ELSE
1903 lmapper%extra_info = imiss
1904 ENDIF
1905 lmapper%time = pstart1 ! = pend1, order by time?
1906 n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
1907! no CYCLE, a single input can contribute to multiple output intervals
1908 ENDIF
1909 ENDDO do_otime2
1910 ENDDO do_otimerange2
1911 ENDDO do_itime2
1912 ENDDO do_itimerange2
1913
1914ENDIF
1915
1916END SUBROUTINE recompute_stat_proc_agg_common
1917
1918
1919SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
1920 max_step, weights)
1921TYPE(datetime),INTENT(in) :: vertime(:)
1922TYPE(datetime),INTENT(in) :: pstart
1923TYPE(datetime),INTENT(in) :: pend
1924LOGICAL,INTENT(in) :: time_mask(:)
1925TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
1926DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
1927
1928INTEGER :: i, nt
1929TYPE(datetime),ALLOCATABLE :: lvertime(:)
1930TYPE(datetime) :: half, nexthalf
1931INTEGER(kind=int_ll) :: dt, tdt
1932
1933nt = count(time_mask)
1934ALLOCATE(lvertime(nt))
1935lvertime = pack(vertime, mask=time_mask)
1936
1937IF (PRESENT(max_step)) THEN
1938! new way
1939! max_step = timedelta_0
1940! DO i = 1, nt - 1
1941! IF (lvertime(i+1) - lvertime(i) > max_step) &
1942! max_step = lvertime(i+1) - lvertime(i)
1943! ENDDO
1944! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
1945! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
1946! old way
1947 IF (nt == 1) THEN
1948 max_step = pend - pstart
1949 ELSE
1950 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
1951 max_step = half - pstart
1952 DO i = 2, nt - 1
1953 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
1954 IF (nexthalf - half > max_step) max_step = nexthalf - half
1955 half = nexthalf
1956 ENDDO
1957 IF (pend - half > max_step) max_step = pend - half
1958 ENDIF
1959
1960ENDIF
1961
1962IF (PRESENT(weights)) THEN
1963 IF (nt == 1) THEN
1964 weights(1) = 1.0
1965 ELSE
1967 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
1969 weights(1) = dble(dt)/dble(tdt)
1970 DO i = 2, nt - 1
1971 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
1973 weights(i) = dble(dt)/dble(tdt)
1974 half = nexthalf
1975 ENDDO
1977 weights(nt) = dble(dt)/dble(tdt)
1978 ENDIF
1979ENDIF
1980
1981END SUBROUTINE compute_stat_proc_agg_sw
1982
1983! get start of period, end of period and reference time from time,
1984! timerange, according to time_definition.
1985SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
1986 pstart, pend, reftime)
1987TYPE(datetime),INTENT(in) :: time
1988TYPE(vol7d_timerange),INTENT(in) :: timerange
1989INTEGER,INTENT(in) :: time_definition
1990TYPE(datetime),INTENT(out) :: reftime
1991TYPE(datetime),INTENT(out) :: pstart
1992TYPE(datetime),INTENT(out) :: pend
1993
1994TYPE(timedelta) :: p1, p2
1995
1996
1997p1 = timedelta_new(sec=timerange%p1) ! end of period
1998p2 = timedelta_new(sec=timerange%p2) ! length of period
1999
2001! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
2002 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
2003 pstart = datetime_miss
2004 pend = datetime_miss
2005 reftime = datetime_miss
2006 RETURN
2007ENDIF
2008
2009IF (time_definition == 0) THEN ! time == reference time
2010 reftime = time
2011 pend = time + p1
2012 pstart = pend - p2
2013ELSE IF (time_definition == 1) THEN ! time == verification time
2014 pend = time
2015 pstart = time - p2
2016 reftime = time - p1
2017ELSE
2018 pstart = datetime_miss
2019 pend = datetime_miss
2020 reftime = datetime_miss
2021ENDIF
2022
2023END SUBROUTINE time_timerange_get_period
2024
2025
2026! get start of period, end of period and reference time from time,
2027! timerange, according to time_definition. step is taken separately
2028! from timerange and can be "popular"
2029SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
2030 pstart, pend, reftime)
2031TYPE(datetime),INTENT(in) :: time
2032TYPE(vol7d_timerange),INTENT(in) :: timerange
2033TYPE(timedelta),INTENT(in) :: step
2034INTEGER,INTENT(in) :: time_definition
2035TYPE(datetime),INTENT(out) :: reftime
2036TYPE(datetime),INTENT(out) :: pstart
2037TYPE(datetime),INTENT(out) :: pend
2038
2039TYPE(timedelta) :: p1
2040
2041
2042p1 = timedelta_new(sec=timerange%p1) ! end of period
2043
2045! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
2046 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
2047 pstart = datetime_miss
2048 pend = datetime_miss
2049 reftime = datetime_miss
2050 RETURN
2051ENDIF
2052
2053IF (time_definition == 0) THEN ! time == reference time
2054 reftime = time
2055 pend = time + p1
2056 pstart = pend - step
2057ELSE IF (time_definition == 1) THEN ! time == verification time
2058 pend = time
2059 pstart = time - step
2060 reftime = time - p1
2061ELSE
2062 pstart = datetime_miss
2063 pend = datetime_miss
2064 reftime = datetime_miss
2065ENDIF
2066
2067END SUBROUTINE time_timerange_get_period_pop
2068
2069
2070! set time, timerange%p1, timerange%p2 according to pstart, pend,
2071! reftime and time_definition.
2072SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
2073 pstart, pend, reftime)
2074TYPE(datetime),INTENT(out) :: time
2075TYPE(vol7d_timerange),INTENT(inout) :: timerange
2076INTEGER,INTENT(in) :: time_definition
2077TYPE(datetime),INTENT(in) :: reftime
2078TYPE(datetime),INTENT(in) :: pstart
2079TYPE(datetime),INTENT(in) :: pend
2080
2081TYPE(timedelta) :: p1, p2
2082INTEGER(kind=int_ll) :: dmsec
2083
2084
2085IF (time_definition == 0) THEN ! time == reference time
2086 time = reftime
2087 p1 = pend - reftime
2088 p2 = pend - pstart
2089ELSE IF (time_definition == 1) THEN ! time == verification time
2090 time = pend
2091 p1 = pend - reftime
2092 p2 = pend - pstart
2093ELSE
2094 time = datetime_miss
2095ENDIF
2096
2097IF (time /= datetime_miss) THEN
2099 timerange%p1 = int(dmsec/1000_int_ll)
2101 timerange%p2 = int(dmsec/1000_int_ll)
2102ELSE
2103 timerange%p1 = imiss
2104 timerange%p2 = imiss
2105ENDIF
2106
2107END SUBROUTINE time_timerange_set_period
2108
2109
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 |