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