libsim Versione 7.2.1

◆ map_distinct_ttr_mapper()

integer function, dimension(size(vect)) map_distinct_ttr_mapper ( type(ttr_mapper), dimension(:), intent(in)  vect,
logical, dimension(:), intent(in), optional  mask,
logical, intent(in), optional  back 
)

map distinct

Definizione alla linea 1114 del file stat_proc_engine.F90.

1115! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1116! authors:
1117! Davide Cesari <dcesari@arpa.emr.it>
1118! Paolo Patruno <ppatruno@arpa.emr.it>
1119
1120! This program is free software; you can redistribute it and/or
1121! modify it under the terms of the GNU General Public License as
1122! published by the Free Software Foundation; either version 2 of
1123! the License, or (at your option) any later version.
1124
1125! This program is distributed in the hope that it will be useful,
1126! but WITHOUT ANY WARRANTY; without even the implied warranty of
1127! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1128! GNU General Public License for more details.
1129
1130! You should have received a copy of the GNU General Public License
1131! along with this program. If not, see <http://www.gnu.org/licenses/>.
1132#include "config.h"
1133
1137MODULE stat_proc_engine
1139USE vol7d_class
1140IMPLICIT NONE
1141
1142! per ogni elemento i,j dell'output, definire n elementi di input ad
1143! esso contribuenti (arrayof_ttr_mapper) con le seguenti informazioni
1144TYPE ttr_mapper
1145 INTEGER :: it=imiss ! k
1146 INTEGER :: itr=imiss ! l
1147 INTEGER :: extra_info=imiss ! dtratio for intervals, extreme for point in time
1148 TYPE(datetime) :: time=datetime_miss ! time for point in time
1149END TYPE ttr_mapper
1150
1151INTERFACE OPERATOR (==)
1152 MODULE PROCEDURE ttr_mapper_eq
1153END INTERFACE
1154
1155INTERFACE OPERATOR (/=)
1156 MODULE PROCEDURE ttr_mapper_ne
1157END INTERFACE
1158
1159INTERFACE OPERATOR (>)
1160 MODULE PROCEDURE ttr_mapper_gt
1161END INTERFACE
1162
1163INTERFACE OPERATOR (<)
1164 MODULE PROCEDURE ttr_mapper_lt
1165END INTERFACE
1166
1167INTERFACE OPERATOR (>=)
1168 MODULE PROCEDURE ttr_mapper_ge
1169END INTERFACE
1170
1171INTERFACE OPERATOR (<=)
1172 MODULE PROCEDURE ttr_mapper_le
1173END INTERFACE
1174
1175#undef VOL7D_POLY_TYPE
1176#undef VOL7D_POLY_TYPES
1177#undef ENABLE_SORT
1178#define VOL7D_POLY_TYPE TYPE(ttr_mapper)
1179#define VOL7D_POLY_TYPES _ttr_mapper
1180#define ENABLE_SORT
1181#include "array_utilities_pre.F90"
1182
1183#define ARRAYOF_ORIGTYPE TYPE(ttr_mapper)
1184#define ARRAYOF_TYPE arrayof_ttr_mapper
1185#define ARRAYOF_ORIGEQ 1
1186#define ARRAYOF_ORIGGT 1
1187#include "arrayof_pre.F90"
1188! from arrayof
1189
1190
1191CONTAINS
1192
1193! simplified comparisons only on time
1194ELEMENTAL FUNCTION ttr_mapper_eq(this, that) RESULT(res)
1195TYPE(ttr_mapper),INTENT(IN) :: this, that
1196LOGICAL :: res
1197
1198res = this%time == that%time
1199
1200END FUNCTION ttr_mapper_eq
1201
1202ELEMENTAL FUNCTION ttr_mapper_ne(this, that) RESULT(res)
1203TYPE(ttr_mapper),INTENT(IN) :: this, that
1204LOGICAL :: res
1205
1206res = this%time /= that%time
1207
1208END FUNCTION ttr_mapper_ne
1209
1210ELEMENTAL FUNCTION ttr_mapper_gt(this, that) RESULT(res)
1211TYPE(ttr_mapper),INTENT(IN) :: this, that
1212LOGICAL :: res
1213
1214res = this%time > that%time
1215
1216END FUNCTION ttr_mapper_gt
1217
1218ELEMENTAL FUNCTION ttr_mapper_lt(this, that) RESULT(res)
1219TYPE(ttr_mapper),INTENT(IN) :: this, that
1220LOGICAL :: res
1221
1222res = this%time < that%time
1223
1224END FUNCTION ttr_mapper_lt
1225
1226ELEMENTAL FUNCTION ttr_mapper_ge(this, that) RESULT(res)
1227TYPE(ttr_mapper),INTENT(IN) :: this, that
1228LOGICAL :: res
1229
1230res = this%time >= that%time
1231
1232END FUNCTION ttr_mapper_ge
1233
1234ELEMENTAL FUNCTION ttr_mapper_le(this, that) RESULT(res)
1235TYPE(ttr_mapper),INTENT(IN) :: this, that
1236LOGICAL :: res
1237
1238res = this%time <= that%time
1239
1240END FUNCTION ttr_mapper_le
1241
1242#include "arrayof_post.F90"
1243#include "array_utilities_inc.F90"
1244
1245
1246! common operations for statistical processing by differences
1247SUBROUTINE recompute_stat_proc_diff_common(itime, itimerange, stat_proc, step, &
1248 otime, otimerange, map_tr, f, keep_tr, time_definition, full_steps, &
1249 start)
1250TYPE(datetime),INTENT(in) :: itime(:)
1251TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1252INTEGER,INTENT(in) :: stat_proc
1253TYPE(timedelta),INTENT(in) :: step
1254TYPE(datetime),POINTER :: otime(:)
1255TYPE(vol7d_timerange),POINTER :: otimerange(:)
1256INTEGER,ALLOCATABLE,INTENT(out) :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
1257INTEGER :: nitr
1258LOGICAL,ALLOCATABLE :: mask_timerange(:)
1259INTEGER,INTENT(in) :: time_definition
1260LOGICAL,INTENT(in),OPTIONAL :: full_steps
1261TYPE(datetime),INTENT(in),OPTIONAL :: start
1262
1263INTEGER :: i, j, k, l, dirtyrep
1264INTEGER :: steps
1265LOGICAL :: lfull_steps, useful
1266TYPE(datetime) :: lstart, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1267TYPE(vol7d_timerange) :: tmptimerange
1268TYPE(arrayof_datetime) :: a_otime
1269TYPE(arrayof_vol7d_timerange) :: a_otimerange
1270
1271! compute length of cumulation step in seconds
1272CALL getval(step, asec=steps)
1273
1274lstart = datetime_miss
1275IF (PRESENT(start)) lstart = start
1276lfull_steps = optio_log(full_steps)
1277
1278! create a mask of suitable timeranges
1279ALLOCATE(mask_timerange(SIZE(itimerange)))
1280mask_timerange(:) = itimerange(:)%timerange == stat_proc &
1281 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
1282 .AND. itimerange(:)%p1 >= 0 &
1283 .AND. itimerange(:)%p2 > 0
1284
1285IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges
1286 mask_timerange(:) = mask_timerange(:) .AND. &
1287 (itimerange(:)%p1 == 0 .OR. & ! all analyses
1288 mod(itimerange(:)%p1, steps) == 0 .OR. & ! end time is mod step
1289 mod(itimerange(:)%p1 - itimerange(:)%p2, steps) == 0) ! start time is mod step
1290ENDIF
1291! mask_timerange includes all candidate timeranges
1292
1293nitr = count(mask_timerange)
1294ALLOCATE(f(nitr))
1295j = 1
1296DO i = 1, nitr
1297 DO WHILE(.NOT.mask_timerange(j))
1298 j = j + 1
1299 ENDDO
1300 f(i) = j
1301 j = j + 1
1302ENDDO
1303
1304! now we have to evaluate time/timerage pairs which do not need processing
1305ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
1306CALL compute_keep_tr()
1307
1308ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
1309map_tr(:,:,:,:,:) = imiss
1310
1311! scan through all possible combinations of time and timerange
1312DO dirtyrep = 1, 2
1313 IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
1314 CALL packarray(a_otime)
1315 CALL packarray(a_otimerange)
1316 CALL sort(a_otime%array)
1317 CALL sort(a_otimerange%array)
1318 ENDIF
1319 DO l = 1, SIZE(itime)
1320 DO k = 1, nitr
1321 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1322 time_definition, pstart2, pend2, reftime2)
1323
1324 DO j = 1, SIZE(itime)
1325 DO i = 1, nitr
1326 useful = .false.
1327 CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
1328 time_definition, pstart1, pend1, reftime1)
1329 tmptimerange = vol7d_timerange_new(timerange=stat_proc)
1330
1331 IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
1332 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
1333 CALL time_timerange_set_period(tmptime, tmptimerange, &
1334 time_definition, pend1, pend2, reftime2)
1335 IF (lfull_steps) THEN
1336 IF (mod(reftime2, step) == timedelta_0) THEN
1337 useful = .true.
1338 ENDIF
1339 ELSE
1340 useful = .true.
1341 ENDIF
1342
1343 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
1344 CALL time_timerange_set_period(tmptime, tmptimerange, &
1345 time_definition, pstart2, pstart1, pstart1)
1346 IF (lfull_steps) THEN
1347 IF (mod(pstart1, step) == timedelta_0) THEN
1348 useful = .true.
1349 ENDIF
1350 ELSE
1351 useful = .true.
1352 ENDIF
1353 ENDIF
1354
1355 ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
1356 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
1357 CALL time_timerange_set_period(tmptime, tmptimerange, &
1358 time_definition, pend1, pend2, reftime2)
1359 IF (lfull_steps) THEN
1360 IF (mod(pend2-reftime2, step) == timedelta_0) THEN
1361 useful = .true.
1362 ENDIF
1363 ELSE
1364 useful = .true.
1365 ENDIF
1366
1367 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! |-=
1368 CALL time_timerange_set_period(tmptime, tmptimerange, &
1369 time_definition, pstart2, pstart1, reftime2)
1370 IF (lfull_steps) THEN
1371 IF (mod(pstart1-reftime2, step) == timedelta_0) THEN
1372 useful = .true.
1373 ENDIF
1374 ELSE
1375 useful = .true.
1376 ENDIF
1377
1378 ENDIF
1379 ENDIF
1380 useful = useful .AND. tmptime /= datetime_miss .AND. &
1381 tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
1382
1383 IF (useful) THEN ! add a_otime, a_otimerange
1384 map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
1385 map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
1386 ENDIF
1387 ENDDO
1388 ENDDO
1389 ENDDO
1390 ENDDO
1391ENDDO
1392
1393! we have to repeat the computation with sorted arrays
1394CALL compute_keep_tr()
1395
1396otime => a_otime%array
1397otimerange => a_otimerange%array
1398! delete local objects keeping the contents
1399CALL delete(a_otime, nodealloc=.true.)
1400CALL delete(a_otimerange, nodealloc=.true.)
1401
1402#ifdef DEBUG
1403CALL l4f_log(l4f_debug, &
1404 'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr,1)))//', '// &
1405 t2c((SIZE(map_tr,2)))//', '// &
1406 t2c((SIZE(map_tr,3)))//', '// &
1407 t2c((SIZE(map_tr,4))))
1408CALL l4f_log(l4f_debug, &
1409 'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr))/2)//', '// &
1410 t2c(count(c_e(map_tr))/2))
1411CALL l4f_log(l4f_debug, &
1412 'recompute_stat_proc_diff, nitr: '//t2c(nitr))
1413CALL l4f_log(l4f_debug, &
1414 'recompute_stat_proc_diff, good timeranges: '//t2c(count(c_e(keep_tr))/2))
1415CALL l4f_log(l4f_debug, &
1416 'recompute_stat_proc_diff, output times: '//t2c(SIZE(otime)))
1417CALL l4f_log(l4f_debug, &
1418 'recompute_stat_proc_diff, output timeranges: '//t2c(SIZE(otimerange)))
1419#endif
1420
1421CONTAINS
1422
1423SUBROUTINE compute_keep_tr()
1424
1425keep_tr(:,:,:) = imiss
1426DO l = 1, SIZE(itime)
1427 DO k = 1, nitr
1428 IF (itimerange(f(k))%p2 == steps) THEN
1429 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1430 time_definition, pstart2, pend2, reftime2)
1431 useful = .false.
1432 IF (reftime2 == pend2) THEN ! analysis
1433 IF (c_e(lstart)) THEN ! in analysis mode start wins over full_steps
1434 IF (mod(reftime2-lstart, step) == timedelta_0) THEN
1435 useful = .true.
1436 ENDIF
1437 ELSE IF (lfull_steps) THEN
1438 IF (mod(reftime2, step) == timedelta_0) THEN
1439 useful = .true.
1440 ENDIF
1441 ELSE
1442 useful = .true.
1443 ENDIF
1444 ELSE ! forecast
1445 IF (lfull_steps) THEN
1446 IF (mod(itimerange(f(k))%p1, steps) == 0) THEN
1447 useful = .true.
1448 ENDIF
1449 ELSE
1450 useful = .true.
1451 ENDIF
1452 ENDIF
1453 IF (useful) THEN
1454! CALL time_timerange_set_period(tmptime, tmptimerange, &
1455! time_definition, pstart2, pend2, reftime2)
1456 keep_tr(k,l,1) = append_unique(a_otime, itime(l))
1457 keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k)))
1458 ENDIF
1459 ENDIF
1460 ENDDO
1461ENDDO
1462
1463END SUBROUTINE compute_keep_tr
1464
1465END SUBROUTINE recompute_stat_proc_diff_common
1466
1467
1468! common operations for statistical processing by metamorphosis
1469SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
1470 otimerange, map_tr)
1471INTEGER,INTENT(in) :: istat_proc
1472TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1473INTEGER,INTENT(in) :: ostat_proc
1474TYPE(vol7d_timerange),POINTER :: otimerange(:)
1475INTEGER,POINTER :: map_tr(:)
1476
1477INTEGER :: i
1478LOGICAL :: tr_mask(SIZE(itimerange))
1479
1480IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1481 ALLOCATE(otimerange(0), map_tr(0))
1482 RETURN
1483ENDIF
1484
1485! useful timeranges
1486tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
1487 .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
1488ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
1489
1490otimerange = pack(itimerange, mask=tr_mask)
1491otimerange(:)%timerange = ostat_proc
1492map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
1493
1494END SUBROUTINE compute_stat_proc_metamorph_common
1495
1496
1497! common operations for statistical processing by aggregation
1498SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
1499 step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
1500TYPE(datetime),INTENT(in) :: itime(:)
1501TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1502INTEGER,INTENT(in) :: stat_proc
1503INTEGER,INTENT(in) :: tri
1504TYPE(timedelta),INTENT(in) :: step
1505INTEGER,INTENT(in) :: time_definition
1506TYPE(datetime),POINTER :: otime(:)
1507TYPE(vol7d_timerange),POINTER :: otimerange(:)
1508TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
1509INTEGER,POINTER,OPTIONAL :: dtratio(:)
1510TYPE(datetime),INTENT(in),OPTIONAL :: start
1511LOGICAL,INTENT(in),OPTIONAL :: full_steps
1512
1513INTEGER :: i, j, k, l, na, nf, n
1514INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
1515INTEGER(kind=int_ll) :: stepms, mstepms
1516LOGICAL :: lforecast
1517TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1518TYPE(arrayof_datetime) :: a_otime
1519TYPE(arrayof_vol7d_timerange) :: a_otimerange
1520TYPE(arrayof_integer) :: a_dtratio
1521LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
1522TYPE(ttr_mapper) :: lmapper
1523CHARACTER(len=8) :: env_var
1524LOGICAL :: climat_behavior
1525
1526
1527! enable bad behavior for climat database (only for instantaneous data)
1528env_var = ''
1529CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
1530climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
1531
1532! compute length of cumulation step in seconds
1533CALL getval(timedelta_depop(step), asec=steps)
1534
1535! create a mask of suitable timeranges
1536ALLOCATE(mask_timerange(SIZE(itimerange)))
1537mask_timerange(:) = itimerange(:)%timerange == tri &
1538 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
1539
1540IF (PRESENT(dtratio)) THEN
1541 WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0)
1542 mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0
1543 ELSEWHERE
1544 mask_timerange(:) = .false.
1545 END WHERE
1546ELSE ! instantaneous
1547 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
1548ENDIF
1549
1550#ifdef DEBUG
1551CALL l4f_log(l4f_debug, &
1552 '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
1553 t2c(count(mask_timerange)))
1554#endif
1555
1556! euristically determine whether we are dealing with an
1557! analysis/observation or a forecast dataset
1558na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
1559nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
1560lforecast = nf >= na
1561#ifdef DEBUG
1562CALL l4f_log(l4f_debug, &
1563 'recompute_stat_proc_agg, na: '//t2c(na)//', nf: '//t2c(nf))
1564#endif
1565
1566! keep only timeranges of one type (really necessary?)
1567IF (lforecast) THEN
1568 CALL l4f_log(l4f_info, &
1569 'recompute_stat_proc_agg, processing in forecast mode')
1570ELSE
1571 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
1572 CALL l4f_log(l4f_info, &
1573 'recompute_stat_proc_agg, processing in analysis mode')
1574ENDIF
1575
1576#ifdef DEBUG
1577CALL l4f_log(l4f_debug, &
1578 '(re)compute_stat_proc_agg, number of useful timeranges: '// &
1579 t2c(count(mask_timerange)))
1580#endif
1581
1582IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1583 ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
1584 IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
1585 RETURN
1586ENDIF
1587
1588! determine start and end of processing period, should work with p2==0
1589lstart = datetime_miss
1590IF (PRESENT(start)) lstart = start
1591lend = itime(SIZE(itime))
1592! compute some quantities
1593maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
1594maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
1595minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
1596IF (time_definition == 0) THEN ! reference time
1597 lend = lend + timedelta_new(sec=maxp1)
1598ENDIF
1599! extend interval at the end in order to include all the data in case
1600! frac_valid<1; must use < and not <= in the DO WHILE loops some lines
1601! below in order to exclude the last full interval which would be empty
1602lend = lend + step
1603IF (lstart == datetime_miss) THEN ! autodetect
1604 lstart = itime(1)
1605! if autodetected, adjust to obtain real absolute start of data
1606 IF (time_definition == 0) THEN ! reference time
1607 lstart = lstart + timedelta_new(sec=minp1mp2)
1608 ELSE ! verification time
1609! go back to start of longest processing interval
1610 lstart = lstart - timedelta_new(sec=maxp2)
1611 ENDIF
1612! full_steps is effective only in analysis mode and when start is not
1613! specified (start by itself in analysis mode implies full_steps with
1614! respect to start instead of absolute full steps)
1615 IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN
1616 lstart = lstart - (mod(lstart, step)) ! round to step, (should be MODULO, not MOD)
1617 ENDIF
1618ENDIF
1619
1620#ifdef DEBUG
1621CALL l4f_log(l4f_debug, &
1622 'recompute_stat_proc_agg, processing period: '//t2c(lstart)//' - '//t2c(lend))
1623#endif
1624
1625! create output time and timerange lists
1626
1627IF (lforecast) THEN ! forecast mode
1628 IF (time_definition == 0) THEN ! reference time
1629 CALL insert(a_otime, itime) ! should I limit to elements itime >= lstart?
1630
1631! apply start shift to timerange, not to reftime
1632! why did we use itime(SIZE(itime)) (last ref time)?
1633! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
1634 CALL getval(lstart-itime(1), asec=dstart)
1635! allow starting before first reftime but restrict dtstart to -steps
1636! dstart = MAX(0, dstart)
1637 IF (dstart < 0) dstart = mod(dstart, steps)
1638 DO p1 = steps + dstart, maxp1, steps
1639 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1640 ENDDO
1641
1642 ELSE ! verification time
1643
1644! verification time in forecast mode is the ugliest case, because we
1645! don't know a priori how many different (thus incompatible) reference
1646! times we have, so some assumption of regularity has to be made. For
1647! this reason msteps, the minimum step between two times, is
1648! computed. We choose to compute it as a difference between itime
1649! elements, it could be also computed as difference of itimerange%p1
1650! elements. But what if it is not modulo steps?
1651 mstepms = steps*1000_int_ll
1652 DO i = 2, SIZE(itime)
1653 CALL getval(itime(i)-itime(i-1), amsec=stepms)
1654 IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN
1655 msteps = stepms/1000_int_ll
1656 IF (mod(steps, msteps) == 0) mstepms = stepms
1657 ENDIF
1658 ENDDO
1659 msteps = mstepms/1000_int_ll
1660
1661 tmptime = lstart + step
1662 DO WHILE(tmptime < lend) ! < since lend has been extended
1663 CALL insert_unique(a_otime, tmptime)
1664 tmptime = tmptime + step
1665 ENDDO
1666
1667! in next loop, we used initially steps instead of msteps (see comment
1668! above), this gave much less combinations of time/timeranges with
1669! possible empty output; we start from msteps instead of steps in
1670! order to include partial initial intervals in case frac_valid<1;
1671! however, as a gemeral rule, for aggregation of forecasts it's better
1672! to use reference time
1673 DO p1 = msteps, maxp1, msteps
1674 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1675 ENDDO
1676
1677 ENDIF
1678
1679ELSE ! analysis mode
1680 tmptime = lstart + step
1681 DO WHILE(tmptime < lend) ! < since lend has been extended
1682 CALL insert_unique(a_otime, tmptime)
1683 tmptime = tmptime + step
1684 ENDDO
1685 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
1686
1687ENDIF
1688
1689CALL packarray(a_otime)
1690CALL packarray(a_otimerange)
1691otime => a_otime%array
1692otimerange => a_otimerange%array
1693CALL sort(otime)
1694CALL sort(otimerange)
1695! delete local objects keeping the contents
1696CALL delete(a_otime, nodealloc=.true.)
1697CALL delete(a_otimerange, nodealloc=.true.)
1698
1699#ifdef DEBUG
1700CALL l4f_log(l4f_debug, &
1701 'recompute_stat_proc_agg, output time and timerange: '//&
1702 t2c(SIZE(otime))//', '//t2c(size(otimerange)))
1703#endif
1704
1705IF (PRESENT(dtratio)) THEN
1706! count the possible i/o interval ratios
1707 DO k = 1, SIZE(itimerange)
1708 IF (itimerange(k)%p2 /= 0) &
1709 CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
1710 ENDDO
1711 CALL packarray(a_dtratio)
1712 dtratio => a_dtratio%array
1713 CALL sort(dtratio)
1714! delete local object keeping the contents
1715 CALL delete(a_dtratio, nodealloc=.true.)
1716
1717#ifdef DEBUG
1718 CALL l4f_log(l4f_debug, &
1719 'recompute_stat_proc_agg, found '//t2c(size(dtratio))// &
1720 ' possible aggregation ratios, from '// &
1721 t2c(dtratio(1))//' to '//t2c(dtratio(SIZE(dtratio))))
1722#endif
1723
1724 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1725 do_itimerange1: DO l = 1, SIZE(itimerange)
1726 IF (.NOT.mask_timerange(l)) cycle do_itimerange1
1727 do_itime1: DO k = 1, SIZE(itime)
1728 CALL time_timerange_get_period(itime(k), itimerange(l), &
1729 time_definition, pstart1, pend1, reftime1)
1730 do_otimerange1: DO j = 1, SIZE(otimerange)
1731 do_otime1: DO i = 1, SIZE(otime)
1732 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1733 time_definition, pstart2, pend2, reftime2)
1734 IF (lforecast) THEN
1735 IF (reftime1 /= reftime2) cycle do_otime1
1736 ENDIF
1737
1738 IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
1739 mod(pstart1-pstart2, pend1-pstart1) == timedelta_0) THEN ! useful
1740 lmapper%it = k
1741 lmapper%itr = l
1742 lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
1743 n = append(map_ttr(i,j), lmapper)
1744 cycle do_itime1 ! can contribute only to a single interval
1745 ENDIF
1746 ENDDO do_otime1
1747 ENDDO do_otimerange1
1748 ENDDO do_itime1
1749 ENDDO do_itimerange1
1750
1751ELSE
1752
1753 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1754 do_itimerange2: DO l = 1, SIZE(itimerange)
1755 IF (.NOT.mask_timerange(l)) cycle do_itimerange2
1756 do_itime2: DO k = 1, SIZE(itime)
1757 CALL time_timerange_get_period(itime(k), itimerange(l), &
1758 time_definition, pstart1, pend1, reftime1)
1759 do_otimerange2: DO j = 1, SIZE(otimerange)
1760 do_otime2: DO i = 1, SIZE(otime)
1761 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1762 time_definition, pstart2, pend2, reftime2)
1763 IF (lforecast) THEN
1764 IF (reftime1 /= reftime2) cycle do_otime2
1765 ENDIF
1766
1767 IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
1768 IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
1769 lmapper%it = k
1770 lmapper%itr = l
1771 IF (pstart1 == pstart2) THEN
1772 lmapper%extra_info = 1 ! start of interval
1773 ELSE IF (pend1 == pend2) THEN
1774 lmapper%extra_info = 2 ! end of interval
1775 ELSE
1776 lmapper%extra_info = imiss
1777 ENDIF
1778 lmapper%time = pstart1 ! = pend1, order by time?
1779 n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
1780! no CYCLE, a single input can contribute to multiple output intervals
1781 ENDIF
1782 ENDDO do_otime2
1783 ENDDO do_otimerange2
1784 ENDDO do_itime2
1785 ENDDO do_itimerange2
1786
1787ENDIF
1788
1789END SUBROUTINE recompute_stat_proc_agg_common
1790
1791
1792SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
1793 max_step, weights)
1794TYPE(datetime),INTENT(in) :: vertime(:)
1795TYPE(datetime),INTENT(in) :: pstart
1796TYPE(datetime),INTENT(in) :: pend
1797LOGICAL,INTENT(in) :: time_mask(:)
1798TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
1799DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
1800
1801INTEGER :: i, nt
1802TYPE(datetime),ALLOCATABLE :: lvertime(:)
1803TYPE(datetime) :: half, nexthalf
1804INTEGER(kind=int_ll) :: dt, tdt
1805
1806nt = count(time_mask)
1807ALLOCATE(lvertime(nt))
1808lvertime = pack(vertime, mask=time_mask)
1809
1810IF (PRESENT(max_step)) THEN
1811! new way
1812! max_step = timedelta_0
1813! DO i = 1, nt - 1
1814! IF (lvertime(i+1) - lvertime(i) > max_step) &
1815! max_step = lvertime(i+1) - lvertime(i)
1816! ENDDO
1817! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
1818! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
1819! old way
1820 IF (nt == 1) THEN
1821 max_step = pend - pstart
1822 ELSE
1823 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
1824 max_step = half - pstart
1825 DO i = 2, nt - 1
1826 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
1827 IF (nexthalf - half > max_step) max_step = nexthalf - half
1828 half = nexthalf
1829 ENDDO
1830 IF (pend - half > max_step) max_step = pend - half
1831 ENDIF
1832
1833ENDIF
1834
1835IF (PRESENT(weights)) THEN
1836 IF (nt == 1) THEN
1837 weights(1) = 1.0
1838 ELSE
1839 CALL getval(pend - pstart, amsec=tdt)
1840 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
1841 CALL getval(half - pstart, amsec=dt)
1842 weights(1) = dble(dt)/dble(tdt)
1843 DO i = 2, nt - 1
1844 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
1845 CALL getval(nexthalf - half, amsec=dt)
1846 weights(i) = dble(dt)/dble(tdt)
1847 half = nexthalf
1848 ENDDO
1849 CALL getval(pend - half, amsec=dt)
1850 weights(nt) = dble(dt)/dble(tdt)
1851 ENDIF
1852ENDIF
1853
1854END SUBROUTINE compute_stat_proc_agg_sw
1855
1856! get start of period, end of period and reference time from time,
1857! timerange, according to time_definition.
1858SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
1859 pstart, pend, reftime)
1860TYPE(datetime),INTENT(in) :: time
1861TYPE(vol7d_timerange),INTENT(in) :: timerange
1862INTEGER,INTENT(in) :: time_definition
1863TYPE(datetime),INTENT(out) :: reftime
1864TYPE(datetime),INTENT(out) :: pstart
1865TYPE(datetime),INTENT(out) :: pend
1866
1867TYPE(timedelta) :: p1, p2
1868
1869
1870p1 = timedelta_new(sec=timerange%p1) ! end of period
1871p2 = timedelta_new(sec=timerange%p2) ! length of period
1872
1873IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
1874! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
1875 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
1876 pstart = datetime_miss
1877 pend = datetime_miss
1878 reftime = datetime_miss
1879 RETURN
1880ENDIF
1881
1882IF (time_definition == 0) THEN ! time == reference time
1883 reftime = time
1884 pend = time + p1
1885 pstart = pend - p2
1886ELSE IF (time_definition == 1) THEN ! time == verification time
1887 pend = time
1888 pstart = time - p2
1889 reftime = time - p1
1890ELSE
1891 pstart = datetime_miss
1892 pend = datetime_miss
1893 reftime = datetime_miss
1894ENDIF
1895
1896END SUBROUTINE time_timerange_get_period
1897
1898
1899! get start of period, end of period and reference time from time,
1900! timerange, according to time_definition. step is taken separately
1901! from timerange and can be "popular"
1902SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
1903 pstart, pend, reftime)
1904TYPE(datetime),INTENT(in) :: time
1905TYPE(vol7d_timerange),INTENT(in) :: timerange
1906TYPE(timedelta),INTENT(in) :: step
1907INTEGER,INTENT(in) :: time_definition
1908TYPE(datetime),INTENT(out) :: reftime
1909TYPE(datetime),INTENT(out) :: pstart
1910TYPE(datetime),INTENT(out) :: pend
1911
1912TYPE(timedelta) :: p1
1913
1914
1915p1 = timedelta_new(sec=timerange%p1) ! end of period
1916
1917IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
1918! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
1919 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
1920 pstart = datetime_miss
1921 pend = datetime_miss
1922 reftime = datetime_miss
1923 RETURN
1924ENDIF
1925
1926IF (time_definition == 0) THEN ! time == reference time
1927 reftime = time
1928 pend = time + p1
1929 pstart = pend - step
1930ELSE IF (time_definition == 1) THEN ! time == verification time
1931 pend = time
1932 pstart = time - step
1933 reftime = time - p1
1934ELSE
1935 pstart = datetime_miss
1936 pend = datetime_miss
1937 reftime = datetime_miss
1938ENDIF
1939
1940END SUBROUTINE time_timerange_get_period_pop
1941
1942
1943! set time, timerange%p1, timerange%p2 according to pstart, pend,
1944! reftime and time_definition.
1945SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
1946 pstart, pend, reftime)
1947TYPE(datetime),INTENT(out) :: time
1948TYPE(vol7d_timerange),INTENT(inout) :: timerange
1949INTEGER,INTENT(in) :: time_definition
1950TYPE(datetime),INTENT(in) :: reftime
1951TYPE(datetime),INTENT(in) :: pstart
1952TYPE(datetime),INTENT(in) :: pend
1953
1954TYPE(timedelta) :: p1, p2
1955INTEGER(kind=int_ll) :: dmsec
1956
1957
1958IF (time_definition == 0) THEN ! time == reference time
1959 time = reftime
1960 p1 = pend - reftime
1961 p2 = pend - pstart
1962ELSE IF (time_definition == 1) THEN ! time == verification time
1963 time = pend
1964 p1 = pend - reftime
1965 p2 = pend - pstart
1966ELSE
1967 time = datetime_miss
1968ENDIF
1969
1970IF (time /= datetime_miss) THEN
1971 CALL getval(p1, amsec=dmsec) ! end of period
1972 timerange%p1 = int(dmsec/1000_int_ll)
1973 CALL getval(p2, amsec=dmsec) ! length of period
1974 timerange%p2 = int(dmsec/1000_int_ll)
1975ELSE
1976 timerange%p1 = imiss
1977 timerange%p2 = imiss
1978ENDIF
1979
1980END SUBROUTINE time_timerange_set_period
1981
1982
1983END MODULE stat_proc_engine
Restituiscono il valore dell'oggetto nella forma desiderata.
Operatore di resto della divisione.
Functions that return a trimmed CHARACTER representation of the input variable.
Quick method to append an element to the array.
Destructor for finalizing an array object.
Method for inserting elements of the array at a desired position.
Method for packing the array object reducing at a minimum the memory occupation, without destroying i...
Classi per la gestione delle coordinate temporali.
This module contains functions that are only for internal use of the library.
Classe per la gestione di un volume completo di dati osservati.

Generated with Doxygen.