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