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