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