libsim Versione 7.2.4

◆ sort_ttr_mapper()

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

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.

Parametri
[in,out]xdontvector to sort inline

Definizione alla linea 1495 del file stat_proc_engine.F90.

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