libsim Versione 7.2.1
|
◆ sort_ttr_mapper()
Sorts inline into ascending order - Quicksort Quicksort chooses a "pivot" in the set, and explores the array from both ends, looking for a value > pivot with the increasing index, for a value <= pivot with the decreasing index, and swapping them when it has found one of each. The array is then subdivided in 2 ([3]) subsets: { values <= pivot} {pivot} {values > pivot} One then call recursively the program to sort each subset. When the size of the subarray is small enough or the maximum level of recursion is gained, one uses an insertion sort that is faster for very small sets.
Definizione alla linea 1495 del file stat_proc_engine.F90. 1496! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1497! authors:
1498! Davide Cesari <dcesari@arpa.emr.it>
1499! Paolo Patruno <ppatruno@arpa.emr.it>
1500
1501! This program is free software; you can redistribute it and/or
1502! modify it under the terms of the GNU General Public License as
1503! published by the Free Software Foundation; either version 2 of
1504! the License, or (at your option) any later version.
1505
1506! This program is distributed in the hope that it will be useful,
1507! but WITHOUT ANY WARRANTY; without even the implied warranty of
1508! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1509! GNU General Public License for more details.
1510
1511! You should have received a copy of the GNU General Public License
1512! along with this program. If not, see <http://www.gnu.org/licenses/>.
1513#include "config.h"
1514
1521IMPLICIT NONE
1522
1523! per ogni elemento i,j dell'output, definire n elementi di input ad
1524! esso contribuenti (arrayof_ttr_mapper) con le seguenti informazioni
1525TYPE ttr_mapper
1526 INTEGER :: it=imiss ! k
1527 INTEGER :: itr=imiss ! l
1528 INTEGER :: extra_info=imiss ! dtratio for intervals, extreme for point in time
1529 TYPE(datetime) :: time=datetime_miss ! time for point in time
1530END TYPE ttr_mapper
1531
1532INTERFACE OPERATOR (==)
1533 MODULE PROCEDURE ttr_mapper_eq
1534END INTERFACE
1535
1536INTERFACE OPERATOR (/=)
1537 MODULE PROCEDURE ttr_mapper_ne
1538END INTERFACE
1539
1540INTERFACE OPERATOR (>)
1541 MODULE PROCEDURE ttr_mapper_gt
1542END INTERFACE
1543
1544INTERFACE OPERATOR (<)
1545 MODULE PROCEDURE ttr_mapper_lt
1546END INTERFACE
1547
1548INTERFACE OPERATOR (>=)
1549 MODULE PROCEDURE ttr_mapper_ge
1550END INTERFACE
1551
1552INTERFACE OPERATOR (<=)
1553 MODULE PROCEDURE ttr_mapper_le
1554END INTERFACE
1555
1556#undef VOL7D_POLY_TYPE
1557#undef VOL7D_POLY_TYPES
1558#undef ENABLE_SORT
1559#define VOL7D_POLY_TYPE TYPE(ttr_mapper)
1560#define VOL7D_POLY_TYPES _ttr_mapper
1561#define ENABLE_SORT
1562#include "array_utilities_pre.F90"
1563
1564#define ARRAYOF_ORIGTYPE TYPE(ttr_mapper)
1565#define ARRAYOF_TYPE arrayof_ttr_mapper
1566#define ARRAYOF_ORIGEQ 1
1567#define ARRAYOF_ORIGGT 1
1568#include "arrayof_pre.F90"
1569! from arrayof
1570
1571
1572CONTAINS
1573
1574! simplified comparisons only on time
1575ELEMENTAL FUNCTION ttr_mapper_eq(this, that) RESULT(res)
1576TYPE(ttr_mapper),INTENT(IN) :: this, that
1577LOGICAL :: res
1578
1579res = this%time == that%time
1580
1581END FUNCTION ttr_mapper_eq
1582
1583ELEMENTAL FUNCTION ttr_mapper_ne(this, that) RESULT(res)
1584TYPE(ttr_mapper),INTENT(IN) :: this, that
1585LOGICAL :: res
1586
1587res = this%time /= that%time
1588
1589END FUNCTION ttr_mapper_ne
1590
1591ELEMENTAL FUNCTION ttr_mapper_gt(this, that) RESULT(res)
1592TYPE(ttr_mapper),INTENT(IN) :: this, that
1593LOGICAL :: res
1594
1595res = this%time > that%time
1596
1597END FUNCTION ttr_mapper_gt
1598
1599ELEMENTAL FUNCTION ttr_mapper_lt(this, that) RESULT(res)
1600TYPE(ttr_mapper),INTENT(IN) :: this, that
1601LOGICAL :: res
1602
1603res = this%time < that%time
1604
1605END FUNCTION ttr_mapper_lt
1606
1607ELEMENTAL FUNCTION ttr_mapper_ge(this, that) RESULT(res)
1608TYPE(ttr_mapper),INTENT(IN) :: this, that
1609LOGICAL :: res
1610
1611res = this%time >= that%time
1612
1613END FUNCTION ttr_mapper_ge
1614
1615ELEMENTAL FUNCTION ttr_mapper_le(this, that) RESULT(res)
1616TYPE(ttr_mapper),INTENT(IN) :: this, that
1617LOGICAL :: res
1618
1619res = this%time <= that%time
1620
1621END FUNCTION ttr_mapper_le
1622
1623#include "arrayof_post.F90"
1624#include "array_utilities_inc.F90"
1625
1626
1627! common operations for statistical processing by differences
1628SUBROUTINE recompute_stat_proc_diff_common(itime, itimerange, stat_proc, step, &
1629 otime, otimerange, map_tr, f, keep_tr, time_definition, full_steps, &
1630 start)
1631TYPE(datetime),INTENT(in) :: itime(:)
1632TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1633INTEGER,INTENT(in) :: stat_proc
1634TYPE(timedelta),INTENT(in) :: step
1635TYPE(datetime),POINTER :: otime(:)
1636TYPE(vol7d_timerange),POINTER :: otimerange(:)
1637INTEGER,ALLOCATABLE,INTENT(out) :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
1638INTEGER :: nitr
1639LOGICAL,ALLOCATABLE :: mask_timerange(:)
1640INTEGER,INTENT(in) :: time_definition
1641LOGICAL,INTENT(in),OPTIONAL :: full_steps
1642TYPE(datetime),INTENT(in),OPTIONAL :: start
1643
1644INTEGER :: i, j, k, l, dirtyrep
1645INTEGER :: steps
1646LOGICAL :: lfull_steps, useful
1647TYPE(datetime) :: lstart, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1648TYPE(vol7d_timerange) :: tmptimerange
1649TYPE(arrayof_datetime) :: a_otime
1650TYPE(arrayof_vol7d_timerange) :: a_otimerange
1651
1652! compute length of cumulation step in seconds
1654
1655lstart = datetime_miss
1656IF (PRESENT(start)) lstart = start
1657lfull_steps = optio_log(full_steps)
1658
1659! create a mask of suitable timeranges
1660ALLOCATE(mask_timerange(SIZE(itimerange)))
1661mask_timerange(:) = itimerange(:)%timerange == stat_proc &
1662 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
1663 .AND. itimerange(:)%p1 >= 0 &
1664 .AND. itimerange(:)%p2 > 0
1665
1666IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges
1667 mask_timerange(:) = mask_timerange(:) .AND. &
1668 (itimerange(:)%p1 == 0 .OR. & ! all analyses
1671ENDIF
1672! mask_timerange includes all candidate timeranges
1673
1674nitr = count(mask_timerange)
1675ALLOCATE(f(nitr))
1676j = 1
1677DO i = 1, nitr
1678 DO WHILE(.NOT.mask_timerange(j))
1679 j = j + 1
1680 ENDDO
1681 f(i) = j
1682 j = j + 1
1683ENDDO
1684
1685! now we have to evaluate time/timerage pairs which do not need processing
1686ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
1687CALL compute_keep_tr()
1688
1689ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
1690map_tr(:,:,:,:,:) = imiss
1691
1692! scan through all possible combinations of time and timerange
1693DO dirtyrep = 1, 2
1694 IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
1699 ENDIF
1700 DO l = 1, SIZE(itime)
1701 DO k = 1, nitr
1702 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1703 time_definition, pstart2, pend2, reftime2)
1704
1705 DO j = 1, SIZE(itime)
1706 DO i = 1, nitr
1707 useful = .false.
1708 CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
1709 time_definition, pstart1, pend1, reftime1)
1710 tmptimerange = vol7d_timerange_new(timerange=stat_proc)
1711
1712 IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
1713 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
1714 CALL time_timerange_set_period(tmptime, tmptimerange, &
1715 time_definition, pend1, pend2, reftime2)
1716 IF (lfull_steps) THEN
1718 useful = .true.
1719 ENDIF
1720 ELSE
1721 useful = .true.
1722 ENDIF
1723
1724 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
1725 CALL time_timerange_set_period(tmptime, tmptimerange, &
1726 time_definition, pstart2, pstart1, pstart1)
1727 IF (lfull_steps) THEN
1729 useful = .true.
1730 ENDIF
1731 ELSE
1732 useful = .true.
1733 ENDIF
1734 ENDIF
1735
1736 ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
1737 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
1738 CALL time_timerange_set_period(tmptime, tmptimerange, &
1739 time_definition, pend1, pend2, reftime2)
1740 IF (lfull_steps) THEN
1742 useful = .true.
1743 ENDIF
1744 ELSE
1745 useful = .true.
1746 ENDIF
1747
1748 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! |-=
1749 CALL time_timerange_set_period(tmptime, tmptimerange, &
1750 time_definition, pstart2, pstart1, reftime2)
1751 IF (lfull_steps) THEN
1753 useful = .true.
1754 ENDIF
1755 ELSE
1756 useful = .true.
1757 ENDIF
1758
1759 ENDIF
1760 ENDIF
1761 useful = useful .AND. tmptime /= datetime_miss .AND. &
1762 tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
1763
1764 IF (useful) THEN ! add a_otime, a_otimerange
1765 map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
1766 map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
1767 ENDIF
1768 ENDDO
1769 ENDDO
1770 ENDDO
1771 ENDDO
1772ENDDO
1773
1774! we have to repeat the computation with sorted arrays
1775CALL compute_keep_tr()
1776
1777otime => a_otime%array
1778otimerange => a_otimerange%array
1779! delete local objects keeping the contents
1782
1783#ifdef DEBUG
1784CALL l4f_log(l4f_debug, &
1789CALL l4f_log(l4f_debug, &
1792CALL l4f_log(l4f_debug, &
1794CALL l4f_log(l4f_debug, &
1796CALL l4f_log(l4f_debug, &
1798CALL l4f_log(l4f_debug, &
1800#endif
1801
1802CONTAINS
1803
1804SUBROUTINE compute_keep_tr()
1805
1806keep_tr(:,:,:) = imiss
1807DO l = 1, SIZE(itime)
1808 DO k = 1, nitr
1809 IF (itimerange(f(k))%p2 == steps) THEN
1810 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1811 time_definition, pstart2, pend2, reftime2)
1812 useful = .false.
1813 IF (reftime2 == pend2) THEN ! analysis
1816 useful = .true.
1817 ENDIF
1818 ELSE IF (lfull_steps) THEN
1820 useful = .true.
1821 ENDIF
1822 ELSE
1823 useful = .true.
1824 ENDIF
1825 ELSE ! forecast
1826 IF (lfull_steps) THEN
1828 useful = .true.
1829 ENDIF
1830 ELSE
1831 useful = .true.
1832 ENDIF
1833 ENDIF
1834 IF (useful) THEN
1835! CALL time_timerange_set_period(tmptime, tmptimerange, &
1836! time_definition, pstart2, pend2, reftime2)
1837 keep_tr(k,l,1) = append_unique(a_otime, itime(l))
1838 keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k)))
1839 ENDIF
1840 ENDIF
1841 ENDDO
1842ENDDO
1843
1844END SUBROUTINE compute_keep_tr
1845
1846END SUBROUTINE recompute_stat_proc_diff_common
1847
1848
1849! common operations for statistical processing by metamorphosis
1850SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
1851 otimerange, map_tr)
1852INTEGER,INTENT(in) :: istat_proc
1853TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1854INTEGER,INTENT(in) :: ostat_proc
1855TYPE(vol7d_timerange),POINTER :: otimerange(:)
1856INTEGER,POINTER :: map_tr(:)
1857
1858INTEGER :: i
1859LOGICAL :: tr_mask(SIZE(itimerange))
1860
1861IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1862 ALLOCATE(otimerange(0), map_tr(0))
1863 RETURN
1864ENDIF
1865
1866! useful timeranges
1867tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
1868 .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
1869ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
1870
1871otimerange = pack(itimerange, mask=tr_mask)
1872otimerange(:)%timerange = ostat_proc
1873map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
1874
1875END SUBROUTINE compute_stat_proc_metamorph_common
1876
1877
1878! common operations for statistical processing by aggregation
1879SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
1880 step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
1881TYPE(datetime),INTENT(in) :: itime(:)
1882TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1883INTEGER,INTENT(in) :: stat_proc
1884INTEGER,INTENT(in) :: tri
1885TYPE(timedelta),INTENT(in) :: step
1886INTEGER,INTENT(in) :: time_definition
1887TYPE(datetime),POINTER :: otime(:)
1888TYPE(vol7d_timerange),POINTER :: otimerange(:)
1889TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
1890INTEGER,POINTER,OPTIONAL :: dtratio(:)
1891TYPE(datetime),INTENT(in),OPTIONAL :: start
1892LOGICAL,INTENT(in),OPTIONAL :: full_steps
1893
1894INTEGER :: i, j, k, l, na, nf, n
1895INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
1896INTEGER(kind=int_ll) :: stepms, mstepms
1897LOGICAL :: lforecast
1898TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1899TYPE(arrayof_datetime) :: a_otime
1900TYPE(arrayof_vol7d_timerange) :: a_otimerange
1901TYPE(arrayof_integer) :: a_dtratio
1902LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
1903TYPE(ttr_mapper) :: lmapper
1904CHARACTER(len=8) :: env_var
1905LOGICAL :: climat_behavior
1906
1907
1908! enable bad behavior for climat database (only for instantaneous data)
1909env_var = ''
1910CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
1911climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
1912
1913! compute length of cumulation step in seconds
1915
1916! create a mask of suitable timeranges
1917ALLOCATE(mask_timerange(SIZE(itimerange)))
1918mask_timerange(:) = itimerange(:)%timerange == tri &
1919 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
1920
1921IF (PRESENT(dtratio)) THEN
1922 WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0)
1923 mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0
1924 ELSEWHERE
1925 mask_timerange(:) = .false.
1926 END WHERE
1927ELSE ! instantaneous
1928 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
1929ENDIF
1930
1931#ifdef DEBUG
1932CALL l4f_log(l4f_debug, &
1933 '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
1934 t2c(count(mask_timerange)))
1935#endif
1936
1937! euristically determine whether we are dealing with an
1938! analysis/observation or a forecast dataset
1939na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
1940nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
1941lforecast = nf >= na
1942#ifdef DEBUG
1943CALL l4f_log(l4f_debug, &
1945#endif
1946
1947! keep only timeranges of one type (really necessary?)
1948IF (lforecast) THEN
1949 CALL l4f_log(l4f_info, &
1950 'recompute_stat_proc_agg, processing in forecast mode')
1951ELSE
1952 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
1953 CALL l4f_log(l4f_info, &
1954 'recompute_stat_proc_agg, processing in analysis mode')
1955ENDIF
1956
1957#ifdef DEBUG
1958CALL l4f_log(l4f_debug, &
1959 '(re)compute_stat_proc_agg, number of useful timeranges: '// &
1960 t2c(count(mask_timerange)))
1961#endif
1962
1963IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1964 ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
1965 IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
1966 RETURN
1967ENDIF
1968
1969! determine start and end of processing period, should work with p2==0
1970lstart = datetime_miss
1971IF (PRESENT(start)) lstart = start
1972lend = itime(SIZE(itime))
1973! compute some quantities
1974maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
1975maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
1976minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
1977IF (time_definition == 0) THEN ! reference time
1978 lend = lend + timedelta_new(sec=maxp1)
1979ENDIF
1980! extend interval at the end in order to include all the data in case
1981! frac_valid<1; must use < and not <= in the DO WHILE loops some lines
1982! below in order to exclude the last full interval which would be empty
1983lend = lend + step
1984IF (lstart == datetime_miss) THEN ! autodetect
1985 lstart = itime(1)
1986! if autodetected, adjust to obtain real absolute start of data
1987 IF (time_definition == 0) THEN ! reference time
1988 lstart = lstart + timedelta_new(sec=minp1mp2)
1989 ELSE ! verification time
1990! go back to start of longest processing interval
1991 lstart = lstart - timedelta_new(sec=maxp2)
1992 ENDIF
1993! full_steps is effective only in analysis mode and when start is not
1994! specified (start by itself in analysis mode implies full_steps with
1995! respect to start instead of absolute full steps)
1996 IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN
1998 ENDIF
1999ENDIF
2000
2001#ifdef DEBUG
2002CALL l4f_log(l4f_debug, &
2004#endif
2005
2006! create output time and timerange lists
2007
2008IF (lforecast) THEN ! forecast mode
2009 IF (time_definition == 0) THEN ! reference time
2011
2012! apply start shift to timerange, not to reftime
2013! why did we use itime(SIZE(itime)) (last ref time)?
2014! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
2016! allow starting before first reftime but restrict dtstart to -steps
2017! dstart = MAX(0, dstart)
2019 DO p1 = steps + dstart, maxp1, steps
2020 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
2021 ENDDO
2022
2023 ELSE ! verification time
2024
2025! verification time in forecast mode is the ugliest case, because we
2026! don't know a priori how many different (thus incompatible) reference
2027! times we have, so some assumption of regularity has to be made. For
2028! this reason msteps, the minimum step between two times, is
2029! computed. We choose to compute it as a difference between itime
2030! elements, it could be also computed as difference of itimerange%p1
2031! elements. But what if it is not modulo steps?
2032 mstepms = steps*1000_int_ll
2033 DO i = 2, SIZE(itime)
2035 IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN
2036 msteps = stepms/1000_int_ll
2038 ENDIF
2039 ENDDO
2040 msteps = mstepms/1000_int_ll
2041
2042 tmptime = lstart + step
2043 DO WHILE(tmptime < lend) ! < since lend has been extended
2044 CALL insert_unique(a_otime, tmptime)
2045 tmptime = tmptime + step
2046 ENDDO
2047
2048! in next loop, we used initially steps instead of msteps (see comment
2049! above), this gave much less combinations of time/timeranges with
2050! possible empty output; we start from msteps instead of steps in
2051! order to include partial initial intervals in case frac_valid<1;
2052! however, as a gemeral rule, for aggregation of forecasts it's better
2053! to use reference time
2054 DO p1 = msteps, maxp1, msteps
2055 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
2056 ENDDO
2057
2058 ENDIF
2059
2060ELSE ! analysis mode
2061 tmptime = lstart + step
2062 DO WHILE(tmptime < lend) ! < since lend has been extended
2063 CALL insert_unique(a_otime, tmptime)
2064 tmptime = tmptime + step
2065 ENDDO
2066 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
2067
2068ENDIF
2069
2072otime => a_otime%array
2073otimerange => a_otimerange%array
2076! delete local objects keeping the contents
2079
2080#ifdef DEBUG
2081CALL l4f_log(l4f_debug, &
2082 'recompute_stat_proc_agg, output time and timerange: '//&
2084#endif
2085
2086IF (PRESENT(dtratio)) THEN
2087! count the possible i/o interval ratios
2088 DO k = 1, SIZE(itimerange)
2089 IF (itimerange(k)%p2 /= 0) &
2090 CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
2091 ENDDO
2093 dtratio => a_dtratio%array
2095! delete local object keeping the contents
2097
2098#ifdef DEBUG
2099 CALL l4f_log(l4f_debug, &
2101 ' possible aggregation ratios, from '// &
2103#endif
2104
2105 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
2106 do_itimerange1: DO l = 1, SIZE(itimerange)
2107 IF (.NOT.mask_timerange(l)) cycle do_itimerange1
2108 do_itime1: DO k = 1, SIZE(itime)
2109 CALL time_timerange_get_period(itime(k), itimerange(l), &
2110 time_definition, pstart1, pend1, reftime1)
2111 do_otimerange1: DO j = 1, SIZE(otimerange)
2112 do_otime1: DO i = 1, SIZE(otime)
2113 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
2114 time_definition, pstart2, pend2, reftime2)
2115 IF (lforecast) THEN
2116 IF (reftime1 /= reftime2) cycle do_otime1
2117 ENDIF
2118
2119 IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
2121 lmapper%it = k
2122 lmapper%itr = l
2123 lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
2124 n = append(map_ttr(i,j), lmapper)
2125 cycle do_itime1 ! can contribute only to a single interval
2126 ENDIF
2127 ENDDO do_otime1
2128 ENDDO do_otimerange1
2129 ENDDO do_itime1
2130 ENDDO do_itimerange1
2131
2132ELSE
2133
2134 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
2135 do_itimerange2: DO l = 1, SIZE(itimerange)
2136 IF (.NOT.mask_timerange(l)) cycle do_itimerange2
2137 do_itime2: DO k = 1, SIZE(itime)
2138 CALL time_timerange_get_period(itime(k), itimerange(l), &
2139 time_definition, pstart1, pend1, reftime1)
2140 do_otimerange2: DO j = 1, SIZE(otimerange)
2141 do_otime2: DO i = 1, SIZE(otime)
2142 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
2143 time_definition, pstart2, pend2, reftime2)
2144 IF (lforecast) THEN
2145 IF (reftime1 /= reftime2) cycle do_otime2
2146 ENDIF
2147
2148 IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
2149 IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
2150 lmapper%it = k
2151 lmapper%itr = l
2152 IF (pstart1 == pstart2) THEN
2153 lmapper%extra_info = 1 ! start of interval
2154 ELSE IF (pend1 == pend2) THEN
2155 lmapper%extra_info = 2 ! end of interval
2156 ELSE
2157 lmapper%extra_info = imiss
2158 ENDIF
2159 lmapper%time = pstart1 ! = pend1, order by time?
2160 n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
2161! no CYCLE, a single input can contribute to multiple output intervals
2162 ENDIF
2163 ENDDO do_otime2
2164 ENDDO do_otimerange2
2165 ENDDO do_itime2
2166 ENDDO do_itimerange2
2167
2168ENDIF
2169
2170END SUBROUTINE recompute_stat_proc_agg_common
2171
2172
2173SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
2174 max_step, weights)
2175TYPE(datetime),INTENT(in) :: vertime(:)
2176TYPE(datetime),INTENT(in) :: pstart
2177TYPE(datetime),INTENT(in) :: pend
2178LOGICAL,INTENT(in) :: time_mask(:)
2179TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
2180DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
2181
2182INTEGER :: i, nt
2183TYPE(datetime),ALLOCATABLE :: lvertime(:)
2184TYPE(datetime) :: half, nexthalf
2185INTEGER(kind=int_ll) :: dt, tdt
2186
2187nt = count(time_mask)
2188ALLOCATE(lvertime(nt))
2189lvertime = pack(vertime, mask=time_mask)
2190
2191IF (PRESENT(max_step)) THEN
2192! new way
2193! max_step = timedelta_0
2194! DO i = 1, nt - 1
2195! IF (lvertime(i+1) - lvertime(i) > max_step) &
2196! max_step = lvertime(i+1) - lvertime(i)
2197! ENDDO
2198! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
2199! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
2200! old way
2201 IF (nt == 1) THEN
2202 max_step = pend - pstart
2203 ELSE
2204 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
2205 max_step = half - pstart
2206 DO i = 2, nt - 1
2207 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
2208 IF (nexthalf - half > max_step) max_step = nexthalf - half
2209 half = nexthalf
2210 ENDDO
2211 IF (pend - half > max_step) max_step = pend - half
2212 ENDIF
2213
2214ENDIF
2215
2216IF (PRESENT(weights)) THEN
2217 IF (nt == 1) THEN
2218 weights(1) = 1.0
2219 ELSE
2221 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
2223 weights(1) = dble(dt)/dble(tdt)
2224 DO i = 2, nt - 1
2225 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
2227 weights(i) = dble(dt)/dble(tdt)
2228 half = nexthalf
2229 ENDDO
2231 weights(nt) = dble(dt)/dble(tdt)
2232 ENDIF
2233ENDIF
2234
2235END SUBROUTINE compute_stat_proc_agg_sw
2236
2237! get start of period, end of period and reference time from time,
2238! timerange, according to time_definition.
2239SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
2240 pstart, pend, reftime)
2241TYPE(datetime),INTENT(in) :: time
2242TYPE(vol7d_timerange),INTENT(in) :: timerange
2243INTEGER,INTENT(in) :: time_definition
2244TYPE(datetime),INTENT(out) :: reftime
2245TYPE(datetime),INTENT(out) :: pstart
2246TYPE(datetime),INTENT(out) :: pend
2247
2248TYPE(timedelta) :: p1, p2
2249
2250
2251p1 = timedelta_new(sec=timerange%p1) ! end of period
2252p2 = timedelta_new(sec=timerange%p2) ! length of period
2253
2255! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
2256 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
2257 pstart = datetime_miss
2258 pend = datetime_miss
2259 reftime = datetime_miss
2260 RETURN
2261ENDIF
2262
2263IF (time_definition == 0) THEN ! time == reference time
2264 reftime = time
2265 pend = time + p1
2266 pstart = pend - p2
2267ELSE IF (time_definition == 1) THEN ! time == verification time
2268 pend = time
2269 pstart = time - p2
2270 reftime = time - p1
2271ELSE
2272 pstart = datetime_miss
2273 pend = datetime_miss
2274 reftime = datetime_miss
2275ENDIF
2276
2277END SUBROUTINE time_timerange_get_period
2278
2279
2280! get start of period, end of period and reference time from time,
2281! timerange, according to time_definition. step is taken separately
2282! from timerange and can be "popular"
2283SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
2284 pstart, pend, reftime)
2285TYPE(datetime),INTENT(in) :: time
2286TYPE(vol7d_timerange),INTENT(in) :: timerange
2287TYPE(timedelta),INTENT(in) :: step
2288INTEGER,INTENT(in) :: time_definition
2289TYPE(datetime),INTENT(out) :: reftime
2290TYPE(datetime),INTENT(out) :: pstart
2291TYPE(datetime),INTENT(out) :: pend
2292
2293TYPE(timedelta) :: p1
2294
2295
2296p1 = timedelta_new(sec=timerange%p1) ! end of period
2297
2299! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
2300 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
2301 pstart = datetime_miss
2302 pend = datetime_miss
2303 reftime = datetime_miss
2304 RETURN
2305ENDIF
2306
2307IF (time_definition == 0) THEN ! time == reference time
2308 reftime = time
2309 pend = time + p1
2310 pstart = pend - step
2311ELSE IF (time_definition == 1) THEN ! time == verification time
2312 pend = time
2313 pstart = time - step
2314 reftime = time - p1
2315ELSE
2316 pstart = datetime_miss
2317 pend = datetime_miss
2318 reftime = datetime_miss
2319ENDIF
2320
2321END SUBROUTINE time_timerange_get_period_pop
2322
2323
2324! set time, timerange%p1, timerange%p2 according to pstart, pend,
2325! reftime and time_definition.
2326SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
2327 pstart, pend, reftime)
2328TYPE(datetime),INTENT(out) :: time
2329TYPE(vol7d_timerange),INTENT(inout) :: timerange
2330INTEGER,INTENT(in) :: time_definition
2331TYPE(datetime),INTENT(in) :: reftime
2332TYPE(datetime),INTENT(in) :: pstart
2333TYPE(datetime),INTENT(in) :: pend
2334
2335TYPE(timedelta) :: p1, p2
2336INTEGER(kind=int_ll) :: dmsec
2337
2338
2339IF (time_definition == 0) THEN ! time == reference time
2340 time = reftime
2341 p1 = pend - reftime
2342 p2 = pend - pstart
2343ELSE IF (time_definition == 1) THEN ! time == verification time
2344 time = pend
2345 p1 = pend - reftime
2346 p2 = pend - pstart
2347ELSE
2348 time = datetime_miss
2349ENDIF
2350
2351IF (time /= datetime_miss) THEN
2353 timerange%p1 = int(dmsec/1000_int_ll)
2355 timerange%p2 = int(dmsec/1000_int_ll)
2356ELSE
2357 timerange%p1 = imiss
2358 timerange%p2 = imiss
2359ENDIF
2360
2361END SUBROUTINE time_timerange_set_period
2362
2363
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 |