libsim Versione 7.2.6

◆ index_ttr_mapper()

integer function index_ttr_mapper ( type(ttr_mapper), dimension(:), intent(in) vect,
type(ttr_mapper), intent(in) search,
logical, dimension(:), intent(in), optional mask,
logical, intent(in), optional back,
integer, intent(in), optional cache )

Cerca l'indice del primo o ultimo elemento di vect uguale a search.

Definizione alla linea 1296 del file stat_proc_engine.F90.

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