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