libsim Versione 7.2.1

◆ inssor_ttr_mapper()

subroutine inssor_ttr_mapper ( type(ttr_mapper), dimension (:), intent(inout) xdont)

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
1643MODULE stat_proc_engine
1645USE vol7d_class
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
1776
1777! compute length of cumulation step in seconds
1778CALL getval(step, asec=steps)
1779
1780lstart = datetime_miss
1781IF (PRESENT(start)) lstart = start
1782lfull_steps = optio_log(full_steps)
1783
1784! create a mask of suitable timeranges
1785ALLOCATE(mask_timerange(SIZE(itimerange)))
1786mask_timerange(:) = itimerange(:)%timerange == stat_proc &
1787 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
1788 .AND. itimerange(:)%p1 >= 0 &
1789 .AND. itimerange(:)%p2 > 0
1790
1791IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges
1792 mask_timerange(:) = mask_timerange(:) .AND. &
1793 (itimerange(:)%p1 == 0 .OR. & ! all analyses
1794 mod(itimerange(:)%p1, steps) == 0 .OR. & ! end time is mod step
1795 mod(itimerange(:)%p1 - itimerange(:)%p2, steps) == 0) ! start time is mod step
1796ENDIF
1797! mask_timerange includes all candidate timeranges
1798
1799nitr = count(mask_timerange)
1800ALLOCATE(f(nitr))
1801j = 1
1802DO i = 1, nitr
1803 DO WHILE(.NOT.mask_timerange(j))
1804 j = j + 1
1805 ENDDO
1806 f(i) = j
1807 j = j + 1
1808ENDDO
1809
1810! now we have to evaluate time/timerage pairs which do not need processing
1811ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
1812CALL compute_keep_tr()
1813
1814ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
1815map_tr(:,:,:,:,:) = imiss
1816
1817! scan through all possible combinations of time and timerange
1818DO dirtyrep = 1, 2
1819 IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
1820 CALL packarray(a_otime)
1821 CALL packarray(a_otimerange)
1822 CALL sort(a_otime%array)
1823 CALL sort(a_otimerange%array)
1824 ENDIF
1825 DO l = 1, SIZE(itime)
1826 DO k = 1, nitr
1827 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1828 time_definition, pstart2, pend2, reftime2)
1829
1830 DO j = 1, SIZE(itime)
1831 DO i = 1, nitr
1832 useful = .false.
1833 CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
1834 time_definition, pstart1, pend1, reftime1)
1835 tmptimerange = vol7d_timerange_new(timerange=stat_proc)
1836
1837 IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
1838 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
1839 CALL time_timerange_set_period(tmptime, tmptimerange, &
1840 time_definition, pend1, pend2, reftime2)
1841 IF (lfull_steps) THEN
1842 IF (mod(reftime2, step) == timedelta_0) THEN
1843 useful = .true.
1844 ENDIF
1845 ELSE
1846 useful = .true.
1847 ENDIF
1848
1849 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
1850 CALL time_timerange_set_period(tmptime, tmptimerange, &
1851 time_definition, pstart2, pstart1, pstart1)
1852 IF (lfull_steps) THEN
1853 IF (mod(pstart1, step) == timedelta_0) THEN
1854 useful = .true.
1855 ENDIF
1856 ELSE
1857 useful = .true.
1858 ENDIF
1859 ENDIF
1860
1861 ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
1862 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
1863 CALL time_timerange_set_period(tmptime, tmptimerange, &
1864 time_definition, pend1, pend2, reftime2)
1865 IF (lfull_steps) THEN
1866 IF (mod(pend2-reftime2, step) == timedelta_0) THEN
1867 useful = .true.
1868 ENDIF
1869 ELSE
1870 useful = .true.
1871 ENDIF
1872
1873 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! |-=
1874 CALL time_timerange_set_period(tmptime, tmptimerange, &
1875 time_definition, pstart2, pstart1, reftime2)
1876 IF (lfull_steps) THEN
1877 IF (mod(pstart1-reftime2, step) == timedelta_0) THEN
1878 useful = .true.
1879 ENDIF
1880 ELSE
1881 useful = .true.
1882 ENDIF
1883
1884 ENDIF
1885 ENDIF
1886 useful = useful .AND. tmptime /= datetime_miss .AND. &
1887 tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
1888
1889 IF (useful) THEN ! add a_otime, a_otimerange
1890 map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
1891 map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
1892 ENDIF
1893 ENDDO
1894 ENDDO
1895 ENDDO
1896 ENDDO
1897ENDDO
1898
1899! we have to repeat the computation with sorted arrays
1900CALL compute_keep_tr()
1901
1902otime => a_otime%array
1903otimerange => a_otimerange%array
1904! delete local objects keeping the contents
1905CALL delete(a_otime, nodealloc=.true.)
1906CALL delete(a_otimerange, nodealloc=.true.)
1907
1908#ifdef DEBUG
1909CALL l4f_log(l4f_debug, &
1910 'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr,1)))//', '// &
1911 t2c((SIZE(map_tr,2)))//', '// &
1912 t2c((SIZE(map_tr,3)))//', '// &
1913 t2c((SIZE(map_tr,4))))
1914CALL l4f_log(l4f_debug, &
1915 'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr))/2)//', '// &
1916 t2c(count(c_e(map_tr))/2))
1917CALL l4f_log(l4f_debug, &
1918 'recompute_stat_proc_diff, nitr: '//t2c(nitr))
1919CALL l4f_log(l4f_debug, &
1920 'recompute_stat_proc_diff, good timeranges: '//t2c(count(c_e(keep_tr))/2))
1921CALL l4f_log(l4f_debug, &
1922 'recompute_stat_proc_diff, output times: '//t2c(SIZE(otime)))
1923CALL l4f_log(l4f_debug, &
1924 'recompute_stat_proc_diff, output timeranges: '//t2c(SIZE(otimerange)))
1925#endif
1926
1927CONTAINS
1928
1929SUBROUTINE compute_keep_tr()
1930
1931keep_tr(:,:,:) = imiss
1932DO l = 1, SIZE(itime)
1933 DO k = 1, nitr
1934 IF (itimerange(f(k))%p2 == steps) THEN
1935 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1936 time_definition, pstart2, pend2, reftime2)
1937 useful = .false.
1938 IF (reftime2 == pend2) THEN ! analysis
1939 IF (c_e(lstart)) THEN ! in analysis mode start wins over full_steps
1940 IF (mod(reftime2-lstart, step) == timedelta_0) THEN
1941 useful = .true.
1942 ENDIF
1943 ELSE IF (lfull_steps) THEN
1944 IF (mod(reftime2, step) == timedelta_0) THEN
1945 useful = .true.
1946 ENDIF
1947 ELSE
1948 useful = .true.
1949 ENDIF
1950 ELSE ! forecast
1951 IF (lfull_steps) THEN
1952 IF (mod(itimerange(f(k))%p1, steps) == 0) THEN
1953 useful = .true.
1954 ENDIF
1955 ELSE
1956 useful = .true.
1957 ENDIF
1958 ENDIF
1959 IF (useful) THEN
1960! CALL time_timerange_set_period(tmptime, tmptimerange, &
1961! time_definition, pstart2, pend2, reftime2)
1962 keep_tr(k,l,1) = append_unique(a_otime, itime(l))
1963 keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k)))
1964 ENDIF
1965 ENDIF
1966 ENDDO
1967ENDDO
1968
1969END SUBROUTINE compute_keep_tr
1970
1971END SUBROUTINE recompute_stat_proc_diff_common
1972
1973
1974! common operations for statistical processing by metamorphosis
1975SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
1976 otimerange, map_tr)
1977INTEGER,INTENT(in) :: istat_proc
1978TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1979INTEGER,INTENT(in) :: ostat_proc
1980TYPE(vol7d_timerange),POINTER :: otimerange(:)
1981INTEGER,POINTER :: map_tr(:)
1982
1983INTEGER :: i
1984LOGICAL :: tr_mask(SIZE(itimerange))
1985
1986IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1987 ALLOCATE(otimerange(0), map_tr(0))
1988 RETURN
1989ENDIF
1990
1991! useful timeranges
1992tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
1993 .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
1994ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
1995
1996otimerange = pack(itimerange, mask=tr_mask)
1997otimerange(:)%timerange = ostat_proc
1998map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
1999
2000END SUBROUTINE compute_stat_proc_metamorph_common
2001
2002
2003! common operations for statistical processing by aggregation
2004SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
2005 step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
2006TYPE(datetime),INTENT(in) :: itime(:)
2007TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
2008INTEGER,INTENT(in) :: stat_proc
2009INTEGER,INTENT(in) :: tri
2010TYPE(timedelta),INTENT(in) :: step
2011INTEGER,INTENT(in) :: time_definition
2012TYPE(datetime),POINTER :: otime(:)
2013TYPE(vol7d_timerange),POINTER :: otimerange(:)
2014TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
2015INTEGER,POINTER,OPTIONAL :: dtratio(:)
2016TYPE(datetime),INTENT(in),OPTIONAL :: start
2017LOGICAL,INTENT(in),OPTIONAL :: full_steps
2018
2019INTEGER :: i, j, k, l, na, nf, n
2020INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
2021INTEGER(kind=int_ll) :: stepms, mstepms
2022LOGICAL :: lforecast
2023TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
2024TYPE(arrayof_datetime) :: a_otime
2025TYPE(arrayof_vol7d_timerange) :: a_otimerange
2026TYPE(arrayof_integer) :: a_dtratio
2027LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
2028TYPE(ttr_mapper) :: lmapper
2029CHARACTER(len=8) :: env_var
2030LOGICAL :: climat_behavior
2031
2032
2033! enable bad behavior for climat database (only for instantaneous data)
2034env_var = ''
2035CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
2036climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
2037
2038! compute length of cumulation step in seconds
2039CALL getval(timedelta_depop(step), asec=steps)
2040
2041! create a mask of suitable timeranges
2042ALLOCATE(mask_timerange(SIZE(itimerange)))
2043mask_timerange(:) = itimerange(:)%timerange == tri &
2044 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
2045
2046IF (PRESENT(dtratio)) THEN
2047 WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0)
2048 mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0
2049 ELSEWHERE
2050 mask_timerange(:) = .false.
2051 END WHERE
2052ELSE ! instantaneous
2053 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
2054ENDIF
2055
2056#ifdef DEBUG
2057CALL l4f_log(l4f_debug, &
2058 '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
2059 t2c(count(mask_timerange)))
2060#endif
2061
2062! euristically determine whether we are dealing with an
2063! analysis/observation or a forecast dataset
2064na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
2065nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
2066lforecast = nf >= na
2067#ifdef DEBUG
2068CALL l4f_log(l4f_debug, &
2069 'recompute_stat_proc_agg, na: '//t2c(na)//', nf: '//t2c(nf))
2070#endif
2071
2072! keep only timeranges of one type (really necessary?)
2073IF (lforecast) THEN
2074 CALL l4f_log(l4f_info, &
2075 'recompute_stat_proc_agg, processing in forecast mode')
2076ELSE
2077 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
2078 CALL l4f_log(l4f_info, &
2079 'recompute_stat_proc_agg, processing in analysis mode')
2080ENDIF
2081
2082#ifdef DEBUG
2083CALL l4f_log(l4f_debug, &
2084 '(re)compute_stat_proc_agg, number of useful timeranges: '// &
2085 t2c(count(mask_timerange)))
2086#endif
2087
2088IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
2089 ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
2090 IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
2091 RETURN
2092ENDIF
2093
2094! determine start and end of processing period, should work with p2==0
2095lstart = datetime_miss
2096IF (PRESENT(start)) lstart = start
2097lend = itime(SIZE(itime))
2098! compute some quantities
2099maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
2100maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
2101minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
2102IF (time_definition == 0) THEN ! reference time
2103 lend = lend + timedelta_new(sec=maxp1)
2104ENDIF
2105! extend interval at the end in order to include all the data in case
2106! frac_valid<1; must use < and not <= in the DO WHILE loops some lines
2107! below in order to exclude the last full interval which would be empty
2108lend = lend + step
2109IF (lstart == datetime_miss) THEN ! autodetect
2110 lstart = itime(1)
2111! if autodetected, adjust to obtain real absolute start of data
2112 IF (time_definition == 0) THEN ! reference time
2113 lstart = lstart + timedelta_new(sec=minp1mp2)
2114 ELSE ! verification time
2115! go back to start of longest processing interval
2116 lstart = lstart - timedelta_new(sec=maxp2)
2117 ENDIF
2118! full_steps is effective only in analysis mode and when start is not
2119! specified (start by itself in analysis mode implies full_steps with
2120! respect to start instead of absolute full steps)
2121 IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN
2122 lstart = lstart - (mod(lstart, step)) ! round to step, (should be MODULO, not MOD)
2123 ENDIF
2124ENDIF
2125
2126#ifdef DEBUG
2127CALL l4f_log(l4f_debug, &
2128 'recompute_stat_proc_agg, processing period: '//t2c(lstart)//' - '//t2c(lend))
2129#endif
2130
2131! create output time and timerange lists
2132
2133IF (lforecast) THEN ! forecast mode
2134 IF (time_definition == 0) THEN ! reference time
2135 CALL insert(a_otime, itime) ! should I limit to elements itime >= lstart?
2136
2137! apply start shift to timerange, not to reftime
2138! why did we use itime(SIZE(itime)) (last ref time)?
2139! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
2140 CALL getval(lstart-itime(1), asec=dstart)
2141! allow starting before first reftime but restrict dtstart to -steps
2142! dstart = MAX(0, dstart)
2143 IF (dstart < 0) dstart = mod(dstart, steps)
2144 DO p1 = steps + dstart, maxp1, steps
2145 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
2146 ENDDO
2147
2148 ELSE ! verification time
2149
2150! verification time in forecast mode is the ugliest case, because we
2151! don't know a priori how many different (thus incompatible) reference
2152! times we have, so some assumption of regularity has to be made. For
2153! this reason msteps, the minimum step between two times, is
2154! computed. We choose to compute it as a difference between itime
2155! elements, it could be also computed as difference of itimerange%p1
2156! elements. But what if it is not modulo steps?
2157 mstepms = steps*1000_int_ll
2158 DO i = 2, SIZE(itime)
2159 CALL getval(itime(i)-itime(i-1), amsec=stepms)
2160 IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN
2161 msteps = stepms/1000_int_ll
2162 IF (mod(steps, msteps) == 0) mstepms = stepms
2163 ENDIF
2164 ENDDO
2165 msteps = mstepms/1000_int_ll
2166
2167 tmptime = lstart + step
2168 DO WHILE(tmptime < lend) ! < since lend has been extended
2169 CALL insert_unique(a_otime, tmptime)
2170 tmptime = tmptime + step
2171 ENDDO
2172
2173! in next loop, we used initially steps instead of msteps (see comment
2174! above), this gave much less combinations of time/timeranges with
2175! possible empty output; we start from msteps instead of steps in
2176! order to include partial initial intervals in case frac_valid<1;
2177! however, as a gemeral rule, for aggregation of forecasts it's better
2178! to use reference time
2179 DO p1 = msteps, maxp1, msteps
2180 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
2181 ENDDO
2182
2183 ENDIF
2184
2185ELSE ! analysis mode
2186 tmptime = lstart + step
2187 DO WHILE(tmptime < lend) ! < since lend has been extended
2188 CALL insert_unique(a_otime, tmptime)
2189 tmptime = tmptime + step
2190 ENDDO
2191 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
2192
2193ENDIF
2194
2195CALL packarray(a_otime)
2196CALL packarray(a_otimerange)
2197otime => a_otime%array
2198otimerange => a_otimerange%array
2199CALL sort(otime)
2200CALL sort(otimerange)
2201! delete local objects keeping the contents
2202CALL delete(a_otime, nodealloc=.true.)
2203CALL delete(a_otimerange, nodealloc=.true.)
2204
2205#ifdef DEBUG
2206CALL l4f_log(l4f_debug, &
2207 'recompute_stat_proc_agg, output time and timerange: '//&
2208 t2c(SIZE(otime))//', '//t2c(size(otimerange)))
2209#endif
2210
2211IF (PRESENT(dtratio)) THEN
2212! count the possible i/o interval ratios
2213 DO k = 1, SIZE(itimerange)
2214 IF (itimerange(k)%p2 /= 0) &
2215 CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
2216 ENDDO
2217 CALL packarray(a_dtratio)
2218 dtratio => a_dtratio%array
2219 CALL sort(dtratio)
2220! delete local object keeping the contents
2221 CALL delete(a_dtratio, nodealloc=.true.)
2222
2223#ifdef DEBUG
2224 CALL l4f_log(l4f_debug, &
2225 'recompute_stat_proc_agg, found '//t2c(size(dtratio))// &
2226 ' possible aggregation ratios, from '// &
2227 t2c(dtratio(1))//' to '//t2c(dtratio(SIZE(dtratio))))
2228#endif
2229
2230 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
2231 do_itimerange1: DO l = 1, SIZE(itimerange)
2232 IF (.NOT.mask_timerange(l)) cycle do_itimerange1
2233 do_itime1: DO k = 1, SIZE(itime)
2234 CALL time_timerange_get_period(itime(k), itimerange(l), &
2235 time_definition, pstart1, pend1, reftime1)
2236 do_otimerange1: DO j = 1, SIZE(otimerange)
2237 do_otime1: DO i = 1, SIZE(otime)
2238 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
2239 time_definition, pstart2, pend2, reftime2)
2240 IF (lforecast) THEN
2241 IF (reftime1 /= reftime2) cycle do_otime1
2242 ENDIF
2243
2244 IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
2245 mod(pstart1-pstart2, pend1-pstart1) == timedelta_0) THEN ! useful
2246 lmapper%it = k
2247 lmapper%itr = l
2248 lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
2249 n = append(map_ttr(i,j), lmapper)
2250 cycle do_itime1 ! can contribute only to a single interval
2251 ENDIF
2252 ENDDO do_otime1
2253 ENDDO do_otimerange1
2254 ENDDO do_itime1
2255 ENDDO do_itimerange1
2256
2257ELSE
2258
2259 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
2260 do_itimerange2: DO l = 1, SIZE(itimerange)
2261 IF (.NOT.mask_timerange(l)) cycle do_itimerange2
2262 do_itime2: DO k = 1, SIZE(itime)
2263 CALL time_timerange_get_period(itime(k), itimerange(l), &
2264 time_definition, pstart1, pend1, reftime1)
2265 do_otimerange2: DO j = 1, SIZE(otimerange)
2266 do_otime2: DO i = 1, SIZE(otime)
2267 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
2268 time_definition, pstart2, pend2, reftime2)
2269 IF (lforecast) THEN
2270 IF (reftime1 /= reftime2) cycle do_otime2
2271 ENDIF
2272
2273 IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
2274 IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
2275 lmapper%it = k
2276 lmapper%itr = l
2277 IF (pstart1 == pstart2) THEN
2278 lmapper%extra_info = 1 ! start of interval
2279 ELSE IF (pend1 == pend2) THEN
2280 lmapper%extra_info = 2 ! end of interval
2281 ELSE
2282 lmapper%extra_info = imiss
2283 ENDIF
2284 lmapper%time = pstart1 ! = pend1, order by time?
2285 n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
2286! no CYCLE, a single input can contribute to multiple output intervals
2287 ENDIF
2288 ENDDO do_otime2
2289 ENDDO do_otimerange2
2290 ENDDO do_itime2
2291 ENDDO do_itimerange2
2292
2293ENDIF
2294
2295END SUBROUTINE recompute_stat_proc_agg_common
2296
2297
2298SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
2299 max_step, weights)
2300TYPE(datetime),INTENT(in) :: vertime(:)
2301TYPE(datetime),INTENT(in) :: pstart
2302TYPE(datetime),INTENT(in) :: pend
2303LOGICAL,INTENT(in) :: time_mask(:)
2304TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
2305DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
2306
2307INTEGER :: i, nt
2308TYPE(datetime),ALLOCATABLE :: lvertime(:)
2309TYPE(datetime) :: half, nexthalf
2310INTEGER(kind=int_ll) :: dt, tdt
2311
2312nt = count(time_mask)
2313ALLOCATE(lvertime(nt))
2314lvertime = pack(vertime, mask=time_mask)
2315
2316IF (PRESENT(max_step)) THEN
2317! new way
2318! max_step = timedelta_0
2319! DO i = 1, nt - 1
2320! IF (lvertime(i+1) - lvertime(i) > max_step) &
2321! max_step = lvertime(i+1) - lvertime(i)
2322! ENDDO
2323! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
2324! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
2325! old way
2326 IF (nt == 1) THEN
2327 max_step = pend - pstart
2328 ELSE
2329 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
2330 max_step = half - pstart
2331 DO i = 2, nt - 1
2332 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
2333 IF (nexthalf - half > max_step) max_step = nexthalf - half
2334 half = nexthalf
2335 ENDDO
2336 IF (pend - half > max_step) max_step = pend - half
2337 ENDIF
2338
2339ENDIF
2340
2341IF (PRESENT(weights)) THEN
2342 IF (nt == 1) THEN
2343 weights(1) = 1.0
2344 ELSE
2345 CALL getval(pend - pstart, amsec=tdt)
2346 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
2347 CALL getval(half - pstart, amsec=dt)
2348 weights(1) = dble(dt)/dble(tdt)
2349 DO i = 2, nt - 1
2350 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
2351 CALL getval(nexthalf - half, amsec=dt)
2352 weights(i) = dble(dt)/dble(tdt)
2353 half = nexthalf
2354 ENDDO
2355 CALL getval(pend - half, amsec=dt)
2356 weights(nt) = dble(dt)/dble(tdt)
2357 ENDIF
2358ENDIF
2359
2360END SUBROUTINE compute_stat_proc_agg_sw
2361
2362! get start of period, end of period and reference time from time,
2363! timerange, according to time_definition.
2364SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
2365 pstart, pend, reftime)
2366TYPE(datetime),INTENT(in) :: time
2367TYPE(vol7d_timerange),INTENT(in) :: timerange
2368INTEGER,INTENT(in) :: time_definition
2369TYPE(datetime),INTENT(out) :: reftime
2370TYPE(datetime),INTENT(out) :: pstart
2371TYPE(datetime),INTENT(out) :: pend
2372
2373TYPE(timedelta) :: p1, p2
2374
2375
2376p1 = timedelta_new(sec=timerange%p1) ! end of period
2377p2 = timedelta_new(sec=timerange%p2) ! length of period
2378
2379IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
2380! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
2381 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
2382 pstart = datetime_miss
2383 pend = datetime_miss
2384 reftime = datetime_miss
2385 RETURN
2386ENDIF
2387
2388IF (time_definition == 0) THEN ! time == reference time
2389 reftime = time
2390 pend = time + p1
2391 pstart = pend - p2
2392ELSE IF (time_definition == 1) THEN ! time == verification time
2393 pend = time
2394 pstart = time - p2
2395 reftime = time - p1
2396ELSE
2397 pstart = datetime_miss
2398 pend = datetime_miss
2399 reftime = datetime_miss
2400ENDIF
2401
2402END SUBROUTINE time_timerange_get_period
2403
2404
2405! get start of period, end of period and reference time from time,
2406! timerange, according to time_definition. step is taken separately
2407! from timerange and can be "popular"
2408SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
2409 pstart, pend, reftime)
2410TYPE(datetime),INTENT(in) :: time
2411TYPE(vol7d_timerange),INTENT(in) :: timerange
2412TYPE(timedelta),INTENT(in) :: step
2413INTEGER,INTENT(in) :: time_definition
2414TYPE(datetime),INTENT(out) :: reftime
2415TYPE(datetime),INTENT(out) :: pstart
2416TYPE(datetime),INTENT(out) :: pend
2417
2418TYPE(timedelta) :: p1
2419
2420
2421p1 = timedelta_new(sec=timerange%p1) ! end of period
2422
2423IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
2424! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
2425 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
2426 pstart = datetime_miss
2427 pend = datetime_miss
2428 reftime = datetime_miss
2429 RETURN
2430ENDIF
2431
2432IF (time_definition == 0) THEN ! time == reference time
2433 reftime = time
2434 pend = time + p1
2435 pstart = pend - step
2436ELSE IF (time_definition == 1) THEN ! time == verification time
2437 pend = time
2438 pstart = time - step
2439 reftime = time - p1
2440ELSE
2441 pstart = datetime_miss
2442 pend = datetime_miss
2443 reftime = datetime_miss
2444ENDIF
2445
2446END SUBROUTINE time_timerange_get_period_pop
2447
2448
2449! set time, timerange%p1, timerange%p2 according to pstart, pend,
2450! reftime and time_definition.
2451SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
2452 pstart, pend, reftime)
2453TYPE(datetime),INTENT(out) :: time
2454TYPE(vol7d_timerange),INTENT(inout) :: timerange
2455INTEGER,INTENT(in) :: time_definition
2456TYPE(datetime),INTENT(in) :: reftime
2457TYPE(datetime),INTENT(in) :: pstart
2458TYPE(datetime),INTENT(in) :: pend
2459
2460TYPE(timedelta) :: p1, p2
2461INTEGER(kind=int_ll) :: dmsec
2462
2463
2464IF (time_definition == 0) THEN ! time == reference time
2465 time = reftime
2466 p1 = pend - reftime
2467 p2 = pend - pstart
2468ELSE IF (time_definition == 1) THEN ! time == verification time
2469 time = pend
2470 p1 = pend - reftime
2471 p2 = pend - pstart
2472ELSE
2473 time = datetime_miss
2474ENDIF
2475
2476IF (time /= datetime_miss) THEN
2477 CALL getval(p1, amsec=dmsec) ! end of period
2478 timerange%p1 = int(dmsec/1000_int_ll)
2479 CALL getval(p2, amsec=dmsec) ! length of period
2480 timerange%p2 = int(dmsec/1000_int_ll)
2481ELSE
2482 timerange%p1 = imiss
2483 timerange%p2 = imiss
2484ENDIF
2485
2486END SUBROUTINE time_timerange_set_period
2487
2488
2489END MODULE stat_proc_engine
Restituiscono il valore dell'oggetto nella forma desiderata.
Operatore di resto della divisione.
Functions that return a trimmed CHARACTER representation of the input variable.
Quick method to append an element to the array.
Destructor for finalizing an array object.
Method for inserting elements of the array at a desired position.
Method for packing the array object reducing at a minimum the memory occupation, without destroying i...
Classi per la gestione delle coordinate temporali.
This module contains functions that are only for internal use of the library.
Classe per la gestione di un volume completo di dati osservati.

Generated with Doxygen.