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