libsim Versione 7.1.11
|
◆ inssor_ttr_mapper()
Sorts into increasing order (Insertion sort) Sorts XDONT into increasing order (Insertion sort) This subroutine uses insertion sort. It does not use any work array and is faster when XDONT is of very small size (< 20), or already almost sorted, so it is used in a final pass when the partial quicksorting has left a sequence of small subsets and that sorting is only necessary within each subset to complete the process. Michel Olagnon - Apr. 2000 Definizione alla linea 1626 del file stat_proc_engine.F90. 1627! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1628! authors:
1629! Davide Cesari <dcesari@arpa.emr.it>
1630! Paolo Patruno <ppatruno@arpa.emr.it>
1631
1632! This program is free software; you can redistribute it and/or
1633! modify it under the terms of the GNU General Public License as
1634! published by the Free Software Foundation; either version 2 of
1635! the License, or (at your option) any later version.
1636
1637! This program is distributed in the hope that it will be useful,
1638! but WITHOUT ANY WARRANTY; without even the implied warranty of
1639! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1640! GNU General Public License for more details.
1641
1642! You should have received a copy of the GNU General Public License
1643! along with this program. If not, see <http://www.gnu.org/licenses/>.
1644#include "config.h"
1645
1652IMPLICIT NONE
1653
1654! per ogni elemento i,j dell'output, definire n elementi di input ad
1655! esso contribuenti (arrayof_ttr_mapper) con le seguenti informazioni
1656TYPE ttr_mapper
1657 INTEGER :: it=imiss ! k
1658 INTEGER :: itr=imiss ! l
1659 INTEGER :: extra_info=imiss ! dtratio for intervals, extreme for point in time
1660 TYPE(datetime) :: time=datetime_miss ! time for point in time
1661END TYPE ttr_mapper
1662
1663INTERFACE OPERATOR (==)
1664 MODULE PROCEDURE ttr_mapper_eq
1665END INTERFACE
1666
1667INTERFACE OPERATOR (/=)
1668 MODULE PROCEDURE ttr_mapper_ne
1669END INTERFACE
1670
1671INTERFACE OPERATOR (>)
1672 MODULE PROCEDURE ttr_mapper_gt
1673END INTERFACE
1674
1675INTERFACE OPERATOR (<)
1676 MODULE PROCEDURE ttr_mapper_lt
1677END INTERFACE
1678
1679INTERFACE OPERATOR (>=)
1680 MODULE PROCEDURE ttr_mapper_ge
1681END INTERFACE
1682
1683INTERFACE OPERATOR (<=)
1684 MODULE PROCEDURE ttr_mapper_le
1685END INTERFACE
1686
1687#undef VOL7D_POLY_TYPE
1688#undef VOL7D_POLY_TYPES
1689#undef ENABLE_SORT
1690#define VOL7D_POLY_TYPE TYPE(ttr_mapper)
1691#define VOL7D_POLY_TYPES _ttr_mapper
1692#define ENABLE_SORT
1693#include "array_utilities_pre.F90"
1694
1695#define ARRAYOF_ORIGTYPE TYPE(ttr_mapper)
1696#define ARRAYOF_TYPE arrayof_ttr_mapper
1697#define ARRAYOF_ORIGEQ 1
1698#define ARRAYOF_ORIGGT 1
1699#include "arrayof_pre.F90"
1700! from arrayof
1701
1702
1703CONTAINS
1704
1705! simplified comparisons only on time
1706ELEMENTAL FUNCTION ttr_mapper_eq(this, that) RESULT(res)
1707TYPE(ttr_mapper),INTENT(IN) :: this, that
1708LOGICAL :: res
1709
1710res = this%time == that%time
1711
1712END FUNCTION ttr_mapper_eq
1713
1714ELEMENTAL FUNCTION ttr_mapper_ne(this, that) RESULT(res)
1715TYPE(ttr_mapper),INTENT(IN) :: this, that
1716LOGICAL :: res
1717
1718res = this%time /= that%time
1719
1720END FUNCTION ttr_mapper_ne
1721
1722ELEMENTAL FUNCTION ttr_mapper_gt(this, that) RESULT(res)
1723TYPE(ttr_mapper),INTENT(IN) :: this, that
1724LOGICAL :: res
1725
1726res = this%time > that%time
1727
1728END FUNCTION ttr_mapper_gt
1729
1730ELEMENTAL FUNCTION ttr_mapper_lt(this, that) RESULT(res)
1731TYPE(ttr_mapper),INTENT(IN) :: this, that
1732LOGICAL :: res
1733
1734res = this%time < that%time
1735
1736END FUNCTION ttr_mapper_lt
1737
1738ELEMENTAL FUNCTION ttr_mapper_ge(this, that) RESULT(res)
1739TYPE(ttr_mapper),INTENT(IN) :: this, that
1740LOGICAL :: res
1741
1742res = this%time >= that%time
1743
1744END FUNCTION ttr_mapper_ge
1745
1746ELEMENTAL FUNCTION ttr_mapper_le(this, that) RESULT(res)
1747TYPE(ttr_mapper),INTENT(IN) :: this, that
1748LOGICAL :: res
1749
1750res = this%time <= that%time
1751
1752END FUNCTION ttr_mapper_le
1753
1754#include "arrayof_post.F90"
1755#include "array_utilities_inc.F90"
1756
1757
1758! common operations for statistical processing by differences
1759SUBROUTINE recompute_stat_proc_diff_common(itime, itimerange, stat_proc, step, &
1760 otime, otimerange, map_tr, f, keep_tr, time_definition, full_steps, &
1761 start)
1762TYPE(datetime),INTENT(in) :: itime(:)
1763TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1764INTEGER,INTENT(in) :: stat_proc
1765TYPE(timedelta),INTENT(in) :: step
1766TYPE(datetime),POINTER :: otime(:)
1767TYPE(vol7d_timerange),POINTER :: otimerange(:)
1768INTEGER,ALLOCATABLE,INTENT(out) :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
1769INTEGER :: nitr
1770LOGICAL,ALLOCATABLE :: mask_timerange(:)
1771INTEGER,INTENT(in) :: time_definition
1772LOGICAL,INTENT(in),OPTIONAL :: full_steps
1773TYPE(datetime),INTENT(in),OPTIONAL :: start
1774
1775INTEGER :: i, j, k, l, dirtyrep
1776INTEGER :: steps
1777LOGICAL :: lfull_steps, useful
1778TYPE(datetime) :: lstart, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1779TYPE(vol7d_timerange) :: tmptimerange
1780TYPE(arrayof_datetime) :: a_otime
1781TYPE(arrayof_vol7d_timerange) :: a_otimerange
1782
1783! compute length of cumulation step in seconds
1785
1786lstart = datetime_miss
1787IF (PRESENT(start)) lstart = start
1788lfull_steps = optio_log(full_steps)
1789
1790! create a mask of suitable timeranges
1791ALLOCATE(mask_timerange(SIZE(itimerange)))
1792mask_timerange(:) = itimerange(:)%timerange == stat_proc &
1793 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
1794 .AND. itimerange(:)%p1 >= 0 &
1795 .AND. itimerange(:)%p2 > 0
1796
1797IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges
1798 mask_timerange(:) = mask_timerange(:) .AND. &
1799 (itimerange(:)%p1 == 0 .OR. & ! all analyses
1802ENDIF
1803! mask_timerange includes all candidate timeranges
1804
1805nitr = count(mask_timerange)
1806ALLOCATE(f(nitr))
1807j = 1
1808DO i = 1, nitr
1809 DO WHILE(.NOT.mask_timerange(j))
1810 j = j + 1
1811 ENDDO
1812 f(i) = j
1813 j = j + 1
1814ENDDO
1815
1816! now we have to evaluate time/timerage pairs which do not need processing
1817ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
1818CALL compute_keep_tr()
1819
1820ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
1821map_tr(:,:,:,:,:) = imiss
1822
1823! scan through all possible combinations of time and timerange
1824DO dirtyrep = 1, 2
1825 IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
1830 ENDIF
1831 DO l = 1, SIZE(itime)
1832 DO k = 1, nitr
1833 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1834 time_definition, pstart2, pend2, reftime2)
1835
1836 DO j = 1, SIZE(itime)
1837 DO i = 1, nitr
1838 useful = .false.
1839 CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
1840 time_definition, pstart1, pend1, reftime1)
1841 tmptimerange = vol7d_timerange_new(timerange=stat_proc)
1842
1843 IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
1844 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
1845 CALL time_timerange_set_period(tmptime, tmptimerange, &
1846 time_definition, pend1, pend2, reftime2)
1847 IF (lfull_steps) THEN
1849 useful = .true.
1850 ENDIF
1851 ELSE
1852 useful = .true.
1853 ENDIF
1854
1855 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
1856 CALL time_timerange_set_period(tmptime, tmptimerange, &
1857 time_definition, pstart2, pstart1, pstart1)
1858 IF (lfull_steps) THEN
1860 useful = .true.
1861 ENDIF
1862 ELSE
1863 useful = .true.
1864 ENDIF
1865 ENDIF
1866
1867 ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
1868 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
1869 CALL time_timerange_set_period(tmptime, tmptimerange, &
1870 time_definition, pend1, pend2, reftime2)
1871 IF (lfull_steps) THEN
1873 useful = .true.
1874 ENDIF
1875 ELSE
1876 useful = .true.
1877 ENDIF
1878
1879 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! |-=
1880 CALL time_timerange_set_period(tmptime, tmptimerange, &
1881 time_definition, pstart2, pstart1, reftime2)
1882 IF (lfull_steps) THEN
1884 useful = .true.
1885 ENDIF
1886 ELSE
1887 useful = .true.
1888 ENDIF
1889
1890 ENDIF
1891 ENDIF
1892 useful = useful .AND. tmptime /= datetime_miss .AND. &
1893 tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
1894
1895 IF (useful) THEN ! add a_otime, a_otimerange
1896 map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
1897 map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
1898 ENDIF
1899 ENDDO
1900 ENDDO
1901 ENDDO
1902 ENDDO
1903ENDDO
1904
1905! we have to repeat the computation with sorted arrays
1906CALL compute_keep_tr()
1907
1908otime => a_otime%array
1909otimerange => a_otimerange%array
1910! delete local objects keeping the contents
1913
1914#ifdef DEBUG
1915CALL l4f_log(l4f_debug, &
1920CALL l4f_log(l4f_debug, &
1923CALL l4f_log(l4f_debug, &
1925CALL l4f_log(l4f_debug, &
1927CALL l4f_log(l4f_debug, &
1929CALL l4f_log(l4f_debug, &
1931#endif
1932
1933CONTAINS
1934
1935SUBROUTINE compute_keep_tr()
1936
1937keep_tr(:,:,:) = imiss
1938DO l = 1, SIZE(itime)
1939 DO k = 1, nitr
1940 IF (itimerange(f(k))%p2 == steps) THEN
1941 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1942 time_definition, pstart2, pend2, reftime2)
1943 useful = .false.
1944 IF (reftime2 == pend2) THEN ! analysis
1947 useful = .true.
1948 ENDIF
1949 ELSE IF (lfull_steps) THEN
1951 useful = .true.
1952 ENDIF
1953 ELSE
1954 useful = .true.
1955 ENDIF
1956 ELSE ! forecast
1957 IF (lfull_steps) THEN
1959 useful = .true.
1960 ENDIF
1961 ELSE
1962 useful = .true.
1963 ENDIF
1964 ENDIF
1965 IF (useful) THEN
1966! CALL time_timerange_set_period(tmptime, tmptimerange, &
1967! time_definition, pstart2, pend2, reftime2)
1968 keep_tr(k,l,1) = append_unique(a_otime, itime(l))
1969 keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k)))
1970 ENDIF
1971 ENDIF
1972 ENDDO
1973ENDDO
1974
1975END SUBROUTINE compute_keep_tr
1976
1977END SUBROUTINE recompute_stat_proc_diff_common
1978
1979
1980! common operations for statistical processing by metamorphosis
1981SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
1982 otimerange, map_tr)
1983INTEGER,INTENT(in) :: istat_proc
1984TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1985INTEGER,INTENT(in) :: ostat_proc
1986TYPE(vol7d_timerange),POINTER :: otimerange(:)
1987INTEGER,POINTER :: map_tr(:)
1988
1989INTEGER :: i
1990LOGICAL :: tr_mask(SIZE(itimerange))
1991
1992IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1993 ALLOCATE(otimerange(0), map_tr(0))
1994 RETURN
1995ENDIF
1996
1997! useful timeranges
1998tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
1999 .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
2000ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
2001
2002otimerange = pack(itimerange, mask=tr_mask)
2003otimerange(:)%timerange = ostat_proc
2004map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
2005
2006END SUBROUTINE compute_stat_proc_metamorph_common
2007
2008
2009! common operations for statistical processing by aggregation
2010SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
2011 step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
2012TYPE(datetime),INTENT(in) :: itime(:)
2013TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
2014INTEGER,INTENT(in) :: stat_proc
2015INTEGER,INTENT(in) :: tri
2016TYPE(timedelta),INTENT(in) :: step
2017INTEGER,INTENT(in) :: time_definition
2018TYPE(datetime),POINTER :: otime(:)
2019TYPE(vol7d_timerange),POINTER :: otimerange(:)
2020TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
2021INTEGER,POINTER,OPTIONAL :: dtratio(:)
2022TYPE(datetime),INTENT(in),OPTIONAL :: start
2023LOGICAL,INTENT(in),OPTIONAL :: full_steps
2024
2025INTEGER :: i, j, k, l, na, nf, n
2026INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
2027INTEGER(kind=int_ll) :: stepms, mstepms
2028LOGICAL :: lforecast
2029TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
2030TYPE(arrayof_datetime) :: a_otime
2031TYPE(arrayof_vol7d_timerange) :: a_otimerange
2032TYPE(arrayof_integer) :: a_dtratio
2033LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
2034TYPE(ttr_mapper) :: lmapper
2035CHARACTER(len=8) :: env_var
2036LOGICAL :: climat_behavior
2037
2038
2039! enable bad behavior for climat database (only for instantaneous data)
2040env_var = ''
2041CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
2042climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
2043
2044! compute length of cumulation step in seconds
2046
2047! create a mask of suitable timeranges
2048ALLOCATE(mask_timerange(SIZE(itimerange)))
2049mask_timerange(:) = itimerange(:)%timerange == tri &
2050 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
2051
2052IF (PRESENT(dtratio)) THEN
2053 WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0)
2054 mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0
2055 ELSEWHERE
2056 mask_timerange(:) = .false.
2057 END WHERE
2058ELSE ! instantaneous
2059 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
2060ENDIF
2061
2062#ifdef DEBUG
2063CALL l4f_log(l4f_debug, &
2064 '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
2065 t2c(count(mask_timerange)))
2066#endif
2067
2068! euristically determine whether we are dealing with an
2069! analysis/observation or a forecast dataset
2070na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
2071nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
2072lforecast = nf >= na
2073#ifdef DEBUG
2074CALL l4f_log(l4f_debug, &
2076#endif
2077
2078! keep only timeranges of one type (really necessary?)
2079IF (lforecast) THEN
2080 CALL l4f_log(l4f_info, &
2081 'recompute_stat_proc_agg, processing in forecast mode')
2082ELSE
2083 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
2084 CALL l4f_log(l4f_info, &
2085 'recompute_stat_proc_agg, processing in analysis mode')
2086ENDIF
2087
2088#ifdef DEBUG
2089CALL l4f_log(l4f_debug, &
2090 '(re)compute_stat_proc_agg, number of useful timeranges: '// &
2091 t2c(count(mask_timerange)))
2092#endif
2093
2094IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
2095 ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
2096 IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
2097 RETURN
2098ENDIF
2099
2100! determine start and end of processing period, should work with p2==0
2101lstart = datetime_miss
2102IF (PRESENT(start)) lstart = start
2103lend = itime(SIZE(itime))
2104! compute some quantities
2105maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
2106maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
2107minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
2108IF (time_definition == 0) THEN ! reference time
2109 lend = lend + timedelta_new(sec=maxp1)
2110ENDIF
2111! extend interval at the end in order to include all the data in case
2112! frac_valid<1; must use < and not <= in the DO WHILE loops some lines
2113! below in order to exclude the last full interval which would be empty
2114lend = lend + step
2115IF (lstart == datetime_miss) THEN ! autodetect
2116 lstart = itime(1)
2117! if autodetected, adjust to obtain real absolute start of data
2118 IF (time_definition == 0) THEN ! reference time
2119 lstart = lstart + timedelta_new(sec=minp1mp2)
2120 ELSE ! verification time
2121! go back to start of longest processing interval
2122 lstart = lstart - timedelta_new(sec=maxp2)
2123 ENDIF
2124! full_steps is effective only in analysis mode and when start is not
2125! specified (start by itself in analysis mode implies full_steps with
2126! respect to start instead of absolute full steps)
2127 IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN
2129 ENDIF
2130ENDIF
2131
2132#ifdef DEBUG
2133CALL l4f_log(l4f_debug, &
2135#endif
2136
2137! create output time and timerange lists
2138
2139IF (lforecast) THEN ! forecast mode
2140 IF (time_definition == 0) THEN ! reference time
2142
2143! apply start shift to timerange, not to reftime
2144! why did we use itime(SIZE(itime)) (last ref time)?
2145! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
2147! allow starting before first reftime but restrict dtstart to -steps
2148! dstart = MAX(0, dstart)
2150 DO p1 = steps + dstart, maxp1, steps
2151 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
2152 ENDDO
2153
2154 ELSE ! verification time
2155
2156! verification time in forecast mode is the ugliest case, because we
2157! don't know a priori how many different (thus incompatible) reference
2158! times we have, so some assumption of regularity has to be made. For
2159! this reason msteps, the minimum step between two times, is
2160! computed. We choose to compute it as a difference between itime
2161! elements, it could be also computed as difference of itimerange%p1
2162! elements. But what if it is not modulo steps?
2163 mstepms = steps*1000_int_ll
2164 DO i = 2, SIZE(itime)
2166 IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN
2167 msteps = stepms/1000_int_ll
2169 ENDIF
2170 ENDDO
2171 msteps = mstepms/1000_int_ll
2172
2173 tmptime = lstart + step
2174 DO WHILE(tmptime < lend) ! < since lend has been extended
2175 CALL insert_unique(a_otime, tmptime)
2176 tmptime = tmptime + step
2177 ENDDO
2178
2179! in next loop, we used initially steps instead of msteps (see comment
2180! above), this gave much less combinations of time/timeranges with
2181! possible empty output; we start from msteps instead of steps in
2182! order to include partial initial intervals in case frac_valid<1;
2183! however, as a gemeral rule, for aggregation of forecasts it's better
2184! to use reference time
2185 DO p1 = msteps, maxp1, msteps
2186 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
2187 ENDDO
2188
2189 ENDIF
2190
2191ELSE ! analysis mode
2192 tmptime = lstart + step
2193 DO WHILE(tmptime < lend) ! < since lend has been extended
2194 CALL insert_unique(a_otime, tmptime)
2195 tmptime = tmptime + step
2196 ENDDO
2197 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
2198
2199ENDIF
2200
2203otime => a_otime%array
2204otimerange => a_otimerange%array
2207! delete local objects keeping the contents
2210
2211#ifdef DEBUG
2212CALL l4f_log(l4f_debug, &
2213 'recompute_stat_proc_agg, output time and timerange: '//&
2215#endif
2216
2217IF (PRESENT(dtratio)) THEN
2218! count the possible i/o interval ratios
2219 DO k = 1, SIZE(itimerange)
2220 IF (itimerange(k)%p2 /= 0) &
2221 CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
2222 ENDDO
2224 dtratio => a_dtratio%array
2226! delete local object keeping the contents
2228
2229#ifdef DEBUG
2230 CALL l4f_log(l4f_debug, &
2232 ' possible aggregation ratios, from '// &
2234#endif
2235
2236 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
2237 do_itimerange1: DO l = 1, SIZE(itimerange)
2238 IF (.NOT.mask_timerange(l)) cycle do_itimerange1
2239 do_itime1: DO k = 1, SIZE(itime)
2240 CALL time_timerange_get_period(itime(k), itimerange(l), &
2241 time_definition, pstart1, pend1, reftime1)
2242 do_otimerange1: DO j = 1, SIZE(otimerange)
2243 do_otime1: DO i = 1, SIZE(otime)
2244 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
2245 time_definition, pstart2, pend2, reftime2)
2246 IF (lforecast) THEN
2247 IF (reftime1 /= reftime2) cycle do_otime1
2248 ENDIF
2249
2250 IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
2252 lmapper%it = k
2253 lmapper%itr = l
2254 lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
2255 n = append(map_ttr(i,j), lmapper)
2256 cycle do_itime1 ! can contribute only to a single interval
2257 ENDIF
2258 ENDDO do_otime1
2259 ENDDO do_otimerange1
2260 ENDDO do_itime1
2261 ENDDO do_itimerange1
2262
2263ELSE
2264
2265 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
2266 do_itimerange2: DO l = 1, SIZE(itimerange)
2267 IF (.NOT.mask_timerange(l)) cycle do_itimerange2
2268 do_itime2: DO k = 1, SIZE(itime)
2269 CALL time_timerange_get_period(itime(k), itimerange(l), &
2270 time_definition, pstart1, pend1, reftime1)
2271 do_otimerange2: DO j = 1, SIZE(otimerange)
2272 do_otime2: DO i = 1, SIZE(otime)
2273 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
2274 time_definition, pstart2, pend2, reftime2)
2275 IF (lforecast) THEN
2276 IF (reftime1 /= reftime2) cycle do_otime2
2277 ENDIF
2278
2279 IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
2280 IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
2281 lmapper%it = k
2282 lmapper%itr = l
2283 IF (pstart1 == pstart2) THEN
2284 lmapper%extra_info = 1 ! start of interval
2285 ELSE IF (pend1 == pend2) THEN
2286 lmapper%extra_info = 2 ! end of interval
2287 ELSE
2288 lmapper%extra_info = imiss
2289 ENDIF
2290 lmapper%time = pstart1 ! = pend1, order by time?
2291 n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
2292! no CYCLE, a single input can contribute to multiple output intervals
2293 ENDIF
2294 ENDDO do_otime2
2295 ENDDO do_otimerange2
2296 ENDDO do_itime2
2297 ENDDO do_itimerange2
2298
2299ENDIF
2300
2301END SUBROUTINE recompute_stat_proc_agg_common
2302
2303
2304SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
2305 max_step, weights)
2306TYPE(datetime),INTENT(in) :: vertime(:)
2307TYPE(datetime),INTENT(in) :: pstart
2308TYPE(datetime),INTENT(in) :: pend
2309LOGICAL,INTENT(in) :: time_mask(:)
2310TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
2311DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
2312
2313INTEGER :: i, nt
2314TYPE(datetime),ALLOCATABLE :: lvertime(:)
2315TYPE(datetime) :: half, nexthalf
2316INTEGER(kind=int_ll) :: dt, tdt
2317
2318nt = count(time_mask)
2319ALLOCATE(lvertime(nt))
2320lvertime = pack(vertime, mask=time_mask)
2321
2322IF (PRESENT(max_step)) THEN
2323! new way
2324! max_step = timedelta_0
2325! DO i = 1, nt - 1
2326! IF (lvertime(i+1) - lvertime(i) > max_step) &
2327! max_step = lvertime(i+1) - lvertime(i)
2328! ENDDO
2329! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
2330! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
2331! old way
2332 IF (nt == 1) THEN
2333 max_step = pend - pstart
2334 ELSE
2335 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
2336 max_step = half - pstart
2337 DO i = 2, nt - 1
2338 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
2339 IF (nexthalf - half > max_step) max_step = nexthalf - half
2340 half = nexthalf
2341 ENDDO
2342 IF (pend - half > max_step) max_step = pend - half
2343 ENDIF
2344
2345ENDIF
2346
2347IF (PRESENT(weights)) THEN
2348 IF (nt == 1) THEN
2349 weights(1) = 1.0
2350 ELSE
2352 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
2354 weights(1) = dble(dt)/dble(tdt)
2355 DO i = 2, nt - 1
2356 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
2358 weights(i) = dble(dt)/dble(tdt)
2359 half = nexthalf
2360 ENDDO
2362 weights(nt) = dble(dt)/dble(tdt)
2363 ENDIF
2364ENDIF
2365
2366END SUBROUTINE compute_stat_proc_agg_sw
2367
2368! get start of period, end of period and reference time from time,
2369! timerange, according to time_definition.
2370SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
2371 pstart, pend, reftime)
2372TYPE(datetime),INTENT(in) :: time
2373TYPE(vol7d_timerange),INTENT(in) :: timerange
2374INTEGER,INTENT(in) :: time_definition
2375TYPE(datetime),INTENT(out) :: reftime
2376TYPE(datetime),INTENT(out) :: pstart
2377TYPE(datetime),INTENT(out) :: pend
2378
2379TYPE(timedelta) :: p1, p2
2380
2381
2382p1 = timedelta_new(sec=timerange%p1) ! end of period
2383p2 = timedelta_new(sec=timerange%p2) ! length of period
2384
2386! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
2387 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
2388 pstart = datetime_miss
2389 pend = datetime_miss
2390 reftime = datetime_miss
2391 RETURN
2392ENDIF
2393
2394IF (time_definition == 0) THEN ! time == reference time
2395 reftime = time
2396 pend = time + p1
2397 pstart = pend - p2
2398ELSE IF (time_definition == 1) THEN ! time == verification time
2399 pend = time
2400 pstart = time - p2
2401 reftime = time - p1
2402ELSE
2403 pstart = datetime_miss
2404 pend = datetime_miss
2405 reftime = datetime_miss
2406ENDIF
2407
2408END SUBROUTINE time_timerange_get_period
2409
2410
2411! get start of period, end of period and reference time from time,
2412! timerange, according to time_definition. step is taken separately
2413! from timerange and can be "popular"
2414SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
2415 pstart, pend, reftime)
2416TYPE(datetime),INTENT(in) :: time
2417TYPE(vol7d_timerange),INTENT(in) :: timerange
2418TYPE(timedelta),INTENT(in) :: step
2419INTEGER,INTENT(in) :: time_definition
2420TYPE(datetime),INTENT(out) :: reftime
2421TYPE(datetime),INTENT(out) :: pstart
2422TYPE(datetime),INTENT(out) :: pend
2423
2424TYPE(timedelta) :: p1
2425
2426
2427p1 = timedelta_new(sec=timerange%p1) ! end of period
2428
2430! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
2431 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
2432 pstart = datetime_miss
2433 pend = datetime_miss
2434 reftime = datetime_miss
2435 RETURN
2436ENDIF
2437
2438IF (time_definition == 0) THEN ! time == reference time
2439 reftime = time
2440 pend = time + p1
2441 pstart = pend - step
2442ELSE IF (time_definition == 1) THEN ! time == verification time
2443 pend = time
2444 pstart = time - step
2445 reftime = time - p1
2446ELSE
2447 pstart = datetime_miss
2448 pend = datetime_miss
2449 reftime = datetime_miss
2450ENDIF
2451
2452END SUBROUTINE time_timerange_get_period_pop
2453
2454
2455! set time, timerange%p1, timerange%p2 according to pstart, pend,
2456! reftime and time_definition.
2457SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
2458 pstart, pend, reftime)
2459TYPE(datetime),INTENT(out) :: time
2460TYPE(vol7d_timerange),INTENT(inout) :: timerange
2461INTEGER,INTENT(in) :: time_definition
2462TYPE(datetime),INTENT(in) :: reftime
2463TYPE(datetime),INTENT(in) :: pstart
2464TYPE(datetime),INTENT(in) :: pend
2465
2466TYPE(timedelta) :: p1, p2
2467INTEGER(kind=int_ll) :: dmsec
2468
2469
2470IF (time_definition == 0) THEN ! time == reference time
2471 time = reftime
2472 p1 = pend - reftime
2473 p2 = pend - pstart
2474ELSE IF (time_definition == 1) THEN ! time == verification time
2475 time = pend
2476 p1 = pend - reftime
2477 p2 = pend - pstart
2478ELSE
2479 time = datetime_miss
2480ENDIF
2481
2482IF (time /= datetime_miss) THEN
2484 timerange%p1 = int(dmsec/1000_int_ll)
2486 timerange%p2 = int(dmsec/1000_int_ll)
2487ELSE
2488 timerange%p1 = imiss
2489 timerange%p2 = imiss
2490ENDIF
2491
2492END SUBROUTINE time_timerange_set_period
2493
2494
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 |