libsim Versione 7.1.11
|
◆ index_sorted_ttr_mapper()
Cerca l'indice del primo o ultimo elemento di vect uguale a search. Definizione alla linea 1379 del file stat_proc_engine.F90. 1381! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1382! authors:
1383! Davide Cesari <dcesari@arpa.emr.it>
1384! Paolo Patruno <ppatruno@arpa.emr.it>
1385
1386! This program is free software; you can redistribute it and/or
1387! modify it under the terms of the GNU General Public License as
1388! published by the Free Software Foundation; either version 2 of
1389! the License, or (at your option) any later version.
1390
1391! This program is distributed in the hope that it will be useful,
1392! but WITHOUT ANY WARRANTY; without even the implied warranty of
1393! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1394! GNU General Public License for more details.
1395
1396! You should have received a copy of the GNU General Public License
1397! along with this program. If not, see <http://www.gnu.org/licenses/>.
1398#include "config.h"
1399
1406IMPLICIT NONE
1407
1408! per ogni elemento i,j dell'output, definire n elementi di input ad
1409! esso contribuenti (arrayof_ttr_mapper) con le seguenti informazioni
1410TYPE ttr_mapper
1411 INTEGER :: it=imiss ! k
1412 INTEGER :: itr=imiss ! l
1413 INTEGER :: extra_info=imiss ! dtratio for intervals, extreme for point in time
1414 TYPE(datetime) :: time=datetime_miss ! time for point in time
1415END TYPE ttr_mapper
1416
1417INTERFACE OPERATOR (==)
1418 MODULE PROCEDURE ttr_mapper_eq
1419END INTERFACE
1420
1421INTERFACE OPERATOR (/=)
1422 MODULE PROCEDURE ttr_mapper_ne
1423END INTERFACE
1424
1425INTERFACE OPERATOR (>)
1426 MODULE PROCEDURE ttr_mapper_gt
1427END INTERFACE
1428
1429INTERFACE OPERATOR (<)
1430 MODULE PROCEDURE ttr_mapper_lt
1431END INTERFACE
1432
1433INTERFACE OPERATOR (>=)
1434 MODULE PROCEDURE ttr_mapper_ge
1435END INTERFACE
1436
1437INTERFACE OPERATOR (<=)
1438 MODULE PROCEDURE ttr_mapper_le
1439END INTERFACE
1440
1441#undef VOL7D_POLY_TYPE
1442#undef VOL7D_POLY_TYPES
1443#undef ENABLE_SORT
1444#define VOL7D_POLY_TYPE TYPE(ttr_mapper)
1445#define VOL7D_POLY_TYPES _ttr_mapper
1446#define ENABLE_SORT
1447#include "array_utilities_pre.F90"
1448
1449#define ARRAYOF_ORIGTYPE TYPE(ttr_mapper)
1450#define ARRAYOF_TYPE arrayof_ttr_mapper
1451#define ARRAYOF_ORIGEQ 1
1452#define ARRAYOF_ORIGGT 1
1453#include "arrayof_pre.F90"
1454! from arrayof
1455
1456
1457CONTAINS
1458
1459! simplified comparisons only on time
1460ELEMENTAL FUNCTION ttr_mapper_eq(this, that) RESULT(res)
1461TYPE(ttr_mapper),INTENT(IN) :: this, that
1462LOGICAL :: res
1463
1464res = this%time == that%time
1465
1466END FUNCTION ttr_mapper_eq
1467
1468ELEMENTAL FUNCTION ttr_mapper_ne(this, that) RESULT(res)
1469TYPE(ttr_mapper),INTENT(IN) :: this, that
1470LOGICAL :: res
1471
1472res = this%time /= that%time
1473
1474END FUNCTION ttr_mapper_ne
1475
1476ELEMENTAL FUNCTION ttr_mapper_gt(this, that) RESULT(res)
1477TYPE(ttr_mapper),INTENT(IN) :: this, that
1478LOGICAL :: res
1479
1480res = this%time > that%time
1481
1482END FUNCTION ttr_mapper_gt
1483
1484ELEMENTAL FUNCTION ttr_mapper_lt(this, that) RESULT(res)
1485TYPE(ttr_mapper),INTENT(IN) :: this, that
1486LOGICAL :: res
1487
1488res = this%time < that%time
1489
1490END FUNCTION ttr_mapper_lt
1491
1492ELEMENTAL FUNCTION ttr_mapper_ge(this, that) RESULT(res)
1493TYPE(ttr_mapper),INTENT(IN) :: this, that
1494LOGICAL :: res
1495
1496res = this%time >= that%time
1497
1498END FUNCTION ttr_mapper_ge
1499
1500ELEMENTAL FUNCTION ttr_mapper_le(this, that) RESULT(res)
1501TYPE(ttr_mapper),INTENT(IN) :: this, that
1502LOGICAL :: res
1503
1504res = this%time <= that%time
1505
1506END FUNCTION ttr_mapper_le
1507
1508#include "arrayof_post.F90"
1509#include "array_utilities_inc.F90"
1510
1511
1512! common operations for statistical processing by differences
1513SUBROUTINE recompute_stat_proc_diff_common(itime, itimerange, stat_proc, step, &
1514 otime, otimerange, map_tr, f, keep_tr, time_definition, full_steps, &
1515 start)
1516TYPE(datetime),INTENT(in) :: itime(:)
1517TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1518INTEGER,INTENT(in) :: stat_proc
1519TYPE(timedelta),INTENT(in) :: step
1520TYPE(datetime),POINTER :: otime(:)
1521TYPE(vol7d_timerange),POINTER :: otimerange(:)
1522INTEGER,ALLOCATABLE,INTENT(out) :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
1523INTEGER :: nitr
1524LOGICAL,ALLOCATABLE :: mask_timerange(:)
1525INTEGER,INTENT(in) :: time_definition
1526LOGICAL,INTENT(in),OPTIONAL :: full_steps
1527TYPE(datetime),INTENT(in),OPTIONAL :: start
1528
1529INTEGER :: i, j, k, l, dirtyrep
1530INTEGER :: steps
1531LOGICAL :: lfull_steps, useful
1532TYPE(datetime) :: lstart, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1533TYPE(vol7d_timerange) :: tmptimerange
1534TYPE(arrayof_datetime) :: a_otime
1535TYPE(arrayof_vol7d_timerange) :: a_otimerange
1536
1537! compute length of cumulation step in seconds
1539
1540lstart = datetime_miss
1541IF (PRESENT(start)) lstart = start
1542lfull_steps = optio_log(full_steps)
1543
1544! create a mask of suitable timeranges
1545ALLOCATE(mask_timerange(SIZE(itimerange)))
1546mask_timerange(:) = itimerange(:)%timerange == stat_proc &
1547 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
1548 .AND. itimerange(:)%p1 >= 0 &
1549 .AND. itimerange(:)%p2 > 0
1550
1551IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges
1552 mask_timerange(:) = mask_timerange(:) .AND. &
1553 (itimerange(:)%p1 == 0 .OR. & ! all analyses
1556ENDIF
1557! mask_timerange includes all candidate timeranges
1558
1559nitr = count(mask_timerange)
1560ALLOCATE(f(nitr))
1561j = 1
1562DO i = 1, nitr
1563 DO WHILE(.NOT.mask_timerange(j))
1564 j = j + 1
1565 ENDDO
1566 f(i) = j
1567 j = j + 1
1568ENDDO
1569
1570! now we have to evaluate time/timerage pairs which do not need processing
1571ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
1572CALL compute_keep_tr()
1573
1574ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
1575map_tr(:,:,:,:,:) = imiss
1576
1577! scan through all possible combinations of time and timerange
1578DO dirtyrep = 1, 2
1579 IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
1584 ENDIF
1585 DO l = 1, SIZE(itime)
1586 DO k = 1, nitr
1587 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1588 time_definition, pstart2, pend2, reftime2)
1589
1590 DO j = 1, SIZE(itime)
1591 DO i = 1, nitr
1592 useful = .false.
1593 CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
1594 time_definition, pstart1, pend1, reftime1)
1595 tmptimerange = vol7d_timerange_new(timerange=stat_proc)
1596
1597 IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
1598 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
1599 CALL time_timerange_set_period(tmptime, tmptimerange, &
1600 time_definition, pend1, pend2, reftime2)
1601 IF (lfull_steps) THEN
1603 useful = .true.
1604 ENDIF
1605 ELSE
1606 useful = .true.
1607 ENDIF
1608
1609 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
1610 CALL time_timerange_set_period(tmptime, tmptimerange, &
1611 time_definition, pstart2, pstart1, pstart1)
1612 IF (lfull_steps) THEN
1614 useful = .true.
1615 ENDIF
1616 ELSE
1617 useful = .true.
1618 ENDIF
1619 ENDIF
1620
1621 ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
1622 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
1623 CALL time_timerange_set_period(tmptime, tmptimerange, &
1624 time_definition, pend1, pend2, reftime2)
1625 IF (lfull_steps) THEN
1627 useful = .true.
1628 ENDIF
1629 ELSE
1630 useful = .true.
1631 ENDIF
1632
1633 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! |-=
1634 CALL time_timerange_set_period(tmptime, tmptimerange, &
1635 time_definition, pstart2, pstart1, reftime2)
1636 IF (lfull_steps) THEN
1638 useful = .true.
1639 ENDIF
1640 ELSE
1641 useful = .true.
1642 ENDIF
1643
1644 ENDIF
1645 ENDIF
1646 useful = useful .AND. tmptime /= datetime_miss .AND. &
1647 tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
1648
1649 IF (useful) THEN ! add a_otime, a_otimerange
1650 map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
1651 map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
1652 ENDIF
1653 ENDDO
1654 ENDDO
1655 ENDDO
1656 ENDDO
1657ENDDO
1658
1659! we have to repeat the computation with sorted arrays
1660CALL compute_keep_tr()
1661
1662otime => a_otime%array
1663otimerange => a_otimerange%array
1664! delete local objects keeping the contents
1667
1668#ifdef DEBUG
1669CALL l4f_log(l4f_debug, &
1674CALL l4f_log(l4f_debug, &
1677CALL l4f_log(l4f_debug, &
1679CALL l4f_log(l4f_debug, &
1681CALL l4f_log(l4f_debug, &
1683CALL l4f_log(l4f_debug, &
1685#endif
1686
1687CONTAINS
1688
1689SUBROUTINE compute_keep_tr()
1690
1691keep_tr(:,:,:) = imiss
1692DO l = 1, SIZE(itime)
1693 DO k = 1, nitr
1694 IF (itimerange(f(k))%p2 == steps) THEN
1695 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1696 time_definition, pstart2, pend2, reftime2)
1697 useful = .false.
1698 IF (reftime2 == pend2) THEN ! analysis
1701 useful = .true.
1702 ENDIF
1703 ELSE IF (lfull_steps) THEN
1705 useful = .true.
1706 ENDIF
1707 ELSE
1708 useful = .true.
1709 ENDIF
1710 ELSE ! forecast
1711 IF (lfull_steps) THEN
1713 useful = .true.
1714 ENDIF
1715 ELSE
1716 useful = .true.
1717 ENDIF
1718 ENDIF
1719 IF (useful) THEN
1720! CALL time_timerange_set_period(tmptime, tmptimerange, &
1721! time_definition, pstart2, pend2, reftime2)
1722 keep_tr(k,l,1) = append_unique(a_otime, itime(l))
1723 keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k)))
1724 ENDIF
1725 ENDIF
1726 ENDDO
1727ENDDO
1728
1729END SUBROUTINE compute_keep_tr
1730
1731END SUBROUTINE recompute_stat_proc_diff_common
1732
1733
1734! common operations for statistical processing by metamorphosis
1735SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
1736 otimerange, map_tr)
1737INTEGER,INTENT(in) :: istat_proc
1738TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1739INTEGER,INTENT(in) :: ostat_proc
1740TYPE(vol7d_timerange),POINTER :: otimerange(:)
1741INTEGER,POINTER :: map_tr(:)
1742
1743INTEGER :: i
1744LOGICAL :: tr_mask(SIZE(itimerange))
1745
1746IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1747 ALLOCATE(otimerange(0), map_tr(0))
1748 RETURN
1749ENDIF
1750
1751! useful timeranges
1752tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
1753 .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
1754ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
1755
1756otimerange = pack(itimerange, mask=tr_mask)
1757otimerange(:)%timerange = ostat_proc
1758map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
1759
1760END SUBROUTINE compute_stat_proc_metamorph_common
1761
1762
1763! common operations for statistical processing by aggregation
1764SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
1765 step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
1766TYPE(datetime),INTENT(in) :: itime(:)
1767TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1768INTEGER,INTENT(in) :: stat_proc
1769INTEGER,INTENT(in) :: tri
1770TYPE(timedelta),INTENT(in) :: step
1771INTEGER,INTENT(in) :: time_definition
1772TYPE(datetime),POINTER :: otime(:)
1773TYPE(vol7d_timerange),POINTER :: otimerange(:)
1774TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
1775INTEGER,POINTER,OPTIONAL :: dtratio(:)
1776TYPE(datetime),INTENT(in),OPTIONAL :: start
1777LOGICAL,INTENT(in),OPTIONAL :: full_steps
1778
1779INTEGER :: i, j, k, l, na, nf, n
1780INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
1781INTEGER(kind=int_ll) :: stepms, mstepms
1782LOGICAL :: lforecast
1783TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1784TYPE(arrayof_datetime) :: a_otime
1785TYPE(arrayof_vol7d_timerange) :: a_otimerange
1786TYPE(arrayof_integer) :: a_dtratio
1787LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
1788TYPE(ttr_mapper) :: lmapper
1789CHARACTER(len=8) :: env_var
1790LOGICAL :: climat_behavior
1791
1792
1793! enable bad behavior for climat database (only for instantaneous data)
1794env_var = ''
1795CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
1796climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
1797
1798! compute length of cumulation step in seconds
1800
1801! create a mask of suitable timeranges
1802ALLOCATE(mask_timerange(SIZE(itimerange)))
1803mask_timerange(:) = itimerange(:)%timerange == tri &
1804 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
1805
1806IF (PRESENT(dtratio)) THEN
1807 WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0)
1808 mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0
1809 ELSEWHERE
1810 mask_timerange(:) = .false.
1811 END WHERE
1812ELSE ! instantaneous
1813 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
1814ENDIF
1815
1816#ifdef DEBUG
1817CALL l4f_log(l4f_debug, &
1818 '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
1819 t2c(count(mask_timerange)))
1820#endif
1821
1822! euristically determine whether we are dealing with an
1823! analysis/observation or a forecast dataset
1824na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
1825nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
1826lforecast = nf >= na
1827#ifdef DEBUG
1828CALL l4f_log(l4f_debug, &
1830#endif
1831
1832! keep only timeranges of one type (really necessary?)
1833IF (lforecast) THEN
1834 CALL l4f_log(l4f_info, &
1835 'recompute_stat_proc_agg, processing in forecast mode')
1836ELSE
1837 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
1838 CALL l4f_log(l4f_info, &
1839 'recompute_stat_proc_agg, processing in analysis mode')
1840ENDIF
1841
1842#ifdef DEBUG
1843CALL l4f_log(l4f_debug, &
1844 '(re)compute_stat_proc_agg, number of useful timeranges: '// &
1845 t2c(count(mask_timerange)))
1846#endif
1847
1848IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1849 ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
1850 IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
1851 RETURN
1852ENDIF
1853
1854! determine start and end of processing period, should work with p2==0
1855lstart = datetime_miss
1856IF (PRESENT(start)) lstart = start
1857lend = itime(SIZE(itime))
1858! compute some quantities
1859maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
1860maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
1861minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
1862IF (time_definition == 0) THEN ! reference time
1863 lend = lend + timedelta_new(sec=maxp1)
1864ENDIF
1865! extend interval at the end in order to include all the data in case
1866! frac_valid<1; must use < and not <= in the DO WHILE loops some lines
1867! below in order to exclude the last full interval which would be empty
1868lend = lend + step
1869IF (lstart == datetime_miss) THEN ! autodetect
1870 lstart = itime(1)
1871! if autodetected, adjust to obtain real absolute start of data
1872 IF (time_definition == 0) THEN ! reference time
1873 lstart = lstart + timedelta_new(sec=minp1mp2)
1874 ELSE ! verification time
1875! go back to start of longest processing interval
1876 lstart = lstart - timedelta_new(sec=maxp2)
1877 ENDIF
1878! full_steps is effective only in analysis mode and when start is not
1879! specified (start by itself in analysis mode implies full_steps with
1880! respect to start instead of absolute full steps)
1881 IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN
1883 ENDIF
1884ENDIF
1885
1886#ifdef DEBUG
1887CALL l4f_log(l4f_debug, &
1889#endif
1890
1891! create output time and timerange lists
1892
1893IF (lforecast) THEN ! forecast mode
1894 IF (time_definition == 0) THEN ! reference time
1896
1897! apply start shift to timerange, not to reftime
1898! why did we use itime(SIZE(itime)) (last ref time)?
1899! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
1901! allow starting before first reftime but restrict dtstart to -steps
1902! dstart = MAX(0, dstart)
1904 DO p1 = steps + dstart, maxp1, steps
1905 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1906 ENDDO
1907
1908 ELSE ! verification time
1909
1910! verification time in forecast mode is the ugliest case, because we
1911! don't know a priori how many different (thus incompatible) reference
1912! times we have, so some assumption of regularity has to be made. For
1913! this reason msteps, the minimum step between two times, is
1914! computed. We choose to compute it as a difference between itime
1915! elements, it could be also computed as difference of itimerange%p1
1916! elements. But what if it is not modulo steps?
1917 mstepms = steps*1000_int_ll
1918 DO i = 2, SIZE(itime)
1920 IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN
1921 msteps = stepms/1000_int_ll
1923 ENDIF
1924 ENDDO
1925 msteps = mstepms/1000_int_ll
1926
1927 tmptime = lstart + step
1928 DO WHILE(tmptime < lend) ! < since lend has been extended
1929 CALL insert_unique(a_otime, tmptime)
1930 tmptime = tmptime + step
1931 ENDDO
1932
1933! in next loop, we used initially steps instead of msteps (see comment
1934! above), this gave much less combinations of time/timeranges with
1935! possible empty output; we start from msteps instead of steps in
1936! order to include partial initial intervals in case frac_valid<1;
1937! however, as a gemeral rule, for aggregation of forecasts it's better
1938! to use reference time
1939 DO p1 = msteps, maxp1, msteps
1940 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1941 ENDDO
1942
1943 ENDIF
1944
1945ELSE ! analysis mode
1946 tmptime = lstart + step
1947 DO WHILE(tmptime < lend) ! < since lend has been extended
1948 CALL insert_unique(a_otime, tmptime)
1949 tmptime = tmptime + step
1950 ENDDO
1951 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
1952
1953ENDIF
1954
1957otime => a_otime%array
1958otimerange => a_otimerange%array
1961! delete local objects keeping the contents
1964
1965#ifdef DEBUG
1966CALL l4f_log(l4f_debug, &
1967 'recompute_stat_proc_agg, output time and timerange: '//&
1969#endif
1970
1971IF (PRESENT(dtratio)) THEN
1972! count the possible i/o interval ratios
1973 DO k = 1, SIZE(itimerange)
1974 IF (itimerange(k)%p2 /= 0) &
1975 CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
1976 ENDDO
1978 dtratio => a_dtratio%array
1980! delete local object keeping the contents
1982
1983#ifdef DEBUG
1984 CALL l4f_log(l4f_debug, &
1986 ' possible aggregation ratios, from '// &
1988#endif
1989
1990 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1991 do_itimerange1: DO l = 1, SIZE(itimerange)
1992 IF (.NOT.mask_timerange(l)) cycle do_itimerange1
1993 do_itime1: DO k = 1, SIZE(itime)
1994 CALL time_timerange_get_period(itime(k), itimerange(l), &
1995 time_definition, pstart1, pend1, reftime1)
1996 do_otimerange1: DO j = 1, SIZE(otimerange)
1997 do_otime1: DO i = 1, SIZE(otime)
1998 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1999 time_definition, pstart2, pend2, reftime2)
2000 IF (lforecast) THEN
2001 IF (reftime1 /= reftime2) cycle do_otime1
2002 ENDIF
2003
2004 IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
2006 lmapper%it = k
2007 lmapper%itr = l
2008 lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
2009 n = append(map_ttr(i,j), lmapper)
2010 cycle do_itime1 ! can contribute only to a single interval
2011 ENDIF
2012 ENDDO do_otime1
2013 ENDDO do_otimerange1
2014 ENDDO do_itime1
2015 ENDDO do_itimerange1
2016
2017ELSE
2018
2019 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
2020 do_itimerange2: DO l = 1, SIZE(itimerange)
2021 IF (.NOT.mask_timerange(l)) cycle do_itimerange2
2022 do_itime2: DO k = 1, SIZE(itime)
2023 CALL time_timerange_get_period(itime(k), itimerange(l), &
2024 time_definition, pstart1, pend1, reftime1)
2025 do_otimerange2: DO j = 1, SIZE(otimerange)
2026 do_otime2: DO i = 1, SIZE(otime)
2027 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
2028 time_definition, pstart2, pend2, reftime2)
2029 IF (lforecast) THEN
2030 IF (reftime1 /= reftime2) cycle do_otime2
2031 ENDIF
2032
2033 IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
2034 IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
2035 lmapper%it = k
2036 lmapper%itr = l
2037 IF (pstart1 == pstart2) THEN
2038 lmapper%extra_info = 1 ! start of interval
2039 ELSE IF (pend1 == pend2) THEN
2040 lmapper%extra_info = 2 ! end of interval
2041 ELSE
2042 lmapper%extra_info = imiss
2043 ENDIF
2044 lmapper%time = pstart1 ! = pend1, order by time?
2045 n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
2046! no CYCLE, a single input can contribute to multiple output intervals
2047 ENDIF
2048 ENDDO do_otime2
2049 ENDDO do_otimerange2
2050 ENDDO do_itime2
2051 ENDDO do_itimerange2
2052
2053ENDIF
2054
2055END SUBROUTINE recompute_stat_proc_agg_common
2056
2057
2058SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
2059 max_step, weights)
2060TYPE(datetime),INTENT(in) :: vertime(:)
2061TYPE(datetime),INTENT(in) :: pstart
2062TYPE(datetime),INTENT(in) :: pend
2063LOGICAL,INTENT(in) :: time_mask(:)
2064TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
2065DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
2066
2067INTEGER :: i, nt
2068TYPE(datetime),ALLOCATABLE :: lvertime(:)
2069TYPE(datetime) :: half, nexthalf
2070INTEGER(kind=int_ll) :: dt, tdt
2071
2072nt = count(time_mask)
2073ALLOCATE(lvertime(nt))
2074lvertime = pack(vertime, mask=time_mask)
2075
2076IF (PRESENT(max_step)) THEN
2077! new way
2078! max_step = timedelta_0
2079! DO i = 1, nt - 1
2080! IF (lvertime(i+1) - lvertime(i) > max_step) &
2081! max_step = lvertime(i+1) - lvertime(i)
2082! ENDDO
2083! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
2084! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
2085! old way
2086 IF (nt == 1) THEN
2087 max_step = pend - pstart
2088 ELSE
2089 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
2090 max_step = half - pstart
2091 DO i = 2, nt - 1
2092 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
2093 IF (nexthalf - half > max_step) max_step = nexthalf - half
2094 half = nexthalf
2095 ENDDO
2096 IF (pend - half > max_step) max_step = pend - half
2097 ENDIF
2098
2099ENDIF
2100
2101IF (PRESENT(weights)) THEN
2102 IF (nt == 1) THEN
2103 weights(1) = 1.0
2104 ELSE
2106 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
2108 weights(1) = dble(dt)/dble(tdt)
2109 DO i = 2, nt - 1
2110 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
2112 weights(i) = dble(dt)/dble(tdt)
2113 half = nexthalf
2114 ENDDO
2116 weights(nt) = dble(dt)/dble(tdt)
2117 ENDIF
2118ENDIF
2119
2120END SUBROUTINE compute_stat_proc_agg_sw
2121
2122! get start of period, end of period and reference time from time,
2123! timerange, according to time_definition.
2124SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
2125 pstart, pend, reftime)
2126TYPE(datetime),INTENT(in) :: time
2127TYPE(vol7d_timerange),INTENT(in) :: timerange
2128INTEGER,INTENT(in) :: time_definition
2129TYPE(datetime),INTENT(out) :: reftime
2130TYPE(datetime),INTENT(out) :: pstart
2131TYPE(datetime),INTENT(out) :: pend
2132
2133TYPE(timedelta) :: p1, p2
2134
2135
2136p1 = timedelta_new(sec=timerange%p1) ! end of period
2137p2 = timedelta_new(sec=timerange%p2) ! length of period
2138
2140! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
2141 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
2142 pstart = datetime_miss
2143 pend = datetime_miss
2144 reftime = datetime_miss
2145 RETURN
2146ENDIF
2147
2148IF (time_definition == 0) THEN ! time == reference time
2149 reftime = time
2150 pend = time + p1
2151 pstart = pend - p2
2152ELSE IF (time_definition == 1) THEN ! time == verification time
2153 pend = time
2154 pstart = time - p2
2155 reftime = time - p1
2156ELSE
2157 pstart = datetime_miss
2158 pend = datetime_miss
2159 reftime = datetime_miss
2160ENDIF
2161
2162END SUBROUTINE time_timerange_get_period
2163
2164
2165! get start of period, end of period and reference time from time,
2166! timerange, according to time_definition. step is taken separately
2167! from timerange and can be "popular"
2168SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
2169 pstart, pend, reftime)
2170TYPE(datetime),INTENT(in) :: time
2171TYPE(vol7d_timerange),INTENT(in) :: timerange
2172TYPE(timedelta),INTENT(in) :: step
2173INTEGER,INTENT(in) :: time_definition
2174TYPE(datetime),INTENT(out) :: reftime
2175TYPE(datetime),INTENT(out) :: pstart
2176TYPE(datetime),INTENT(out) :: pend
2177
2178TYPE(timedelta) :: p1
2179
2180
2181p1 = timedelta_new(sec=timerange%p1) ! end of period
2182
2184! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
2185 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
2186 pstart = datetime_miss
2187 pend = datetime_miss
2188 reftime = datetime_miss
2189 RETURN
2190ENDIF
2191
2192IF (time_definition == 0) THEN ! time == reference time
2193 reftime = time
2194 pend = time + p1
2195 pstart = pend - step
2196ELSE IF (time_definition == 1) THEN ! time == verification time
2197 pend = time
2198 pstart = time - step
2199 reftime = time - p1
2200ELSE
2201 pstart = datetime_miss
2202 pend = datetime_miss
2203 reftime = datetime_miss
2204ENDIF
2205
2206END SUBROUTINE time_timerange_get_period_pop
2207
2208
2209! set time, timerange%p1, timerange%p2 according to pstart, pend,
2210! reftime and time_definition.
2211SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
2212 pstart, pend, reftime)
2213TYPE(datetime),INTENT(out) :: time
2214TYPE(vol7d_timerange),INTENT(inout) :: timerange
2215INTEGER,INTENT(in) :: time_definition
2216TYPE(datetime),INTENT(in) :: reftime
2217TYPE(datetime),INTENT(in) :: pstart
2218TYPE(datetime),INTENT(in) :: pend
2219
2220TYPE(timedelta) :: p1, p2
2221INTEGER(kind=int_ll) :: dmsec
2222
2223
2224IF (time_definition == 0) THEN ! time == reference time
2225 time = reftime
2226 p1 = pend - reftime
2227 p2 = pend - pstart
2228ELSE IF (time_definition == 1) THEN ! time == verification time
2229 time = pend
2230 p1 = pend - reftime
2231 p2 = pend - pstart
2232ELSE
2233 time = datetime_miss
2234ENDIF
2235
2236IF (time /= datetime_miss) THEN
2238 timerange%p1 = int(dmsec/1000_int_ll)
2240 timerange%p2 = int(dmsec/1000_int_ll)
2241ELSE
2242 timerange%p1 = imiss
2243 timerange%p2 = imiss
2244ENDIF
2245
2246END SUBROUTINE time_timerange_set_period
2247
2248
Restituiscono il valore dell'oggetto nella forma desiderata. Definition: datetime_class.F90:328 Functions that return a trimmed CHARACTER representation of the input variable. Definition: datetime_class.F90:355 Quick method to append an element to the array. Definition: stat_proc_engine.F90:371 Destructor for finalizing an array object. Definition: stat_proc_engine.F90:384 Method for inserting elements of the array at a desired position. Definition: stat_proc_engine.F90:362 Method for packing the array object reducing at a minimum the memory occupation, without destroying i... Definition: stat_proc_engine.F90:394 This module contains functions that are only for internal use of the library. Definition: stat_proc_engine.F90:217 Classe per la gestione di un volume completo di dati osservati. Definition: vol7d_class.F90:279 |