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