libsim Versione 7.2.6

◆ pack_distinct_ttr_mapper()

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

compatta gli elementi distinti di vect in un array

Definizione alla linea 965 del file stat_proc_engine.F90.

967! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
968! authors:
969! Davide Cesari <dcesari@arpa.emr.it>
970! Paolo Patruno <ppatruno@arpa.emr.it>
971
972! This program is free software; you can redistribute it and/or
973! modify it under the terms of the GNU General Public License as
974! published by the Free Software Foundation; either version 2 of
975! the License, or (at your option) any later version.
976
977! This program is distributed in the hope that it will be useful,
978! but WITHOUT ANY WARRANTY; without even the implied warranty of
979! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
980! GNU General Public License for more details.
981
982! You should have received a copy of the GNU General Public License
983! along with this program. If not, see <http://www.gnu.org/licenses/>.
984#include "config.h"
985!> This module contains functions that are only for internal use of
986!! the library. It should not be used by user procedures because it is
987!! subject to change
988!! \ingroup vol7d
989MODULE stat_proc_engine
991USE vol7d_class
992IMPLICIT NONE
993
994! per ogni elemento i,j dell'output, definire n elementi di input ad
995! esso contribuenti (arrayof_ttr_mapper) con le seguenti informazioni
996TYPE ttr_mapper
997 INTEGER :: it=imiss ! k
998 INTEGER :: itr=imiss ! l
999 INTEGER :: extra_info=imiss ! dtratio for intervals, extreme for point in time
1000 TYPE(datetime) :: time=datetime_miss ! time for point in time
1001END TYPE ttr_mapper
1002
1003INTERFACE OPERATOR (==)
1004 MODULE PROCEDURE ttr_mapper_eq
1005END INTERFACE
1006
1007INTERFACE OPERATOR (/=)
1008 MODULE PROCEDURE ttr_mapper_ne
1009END INTERFACE
1010
1011INTERFACE OPERATOR (>)
1012 MODULE PROCEDURE ttr_mapper_gt
1013END INTERFACE
1014
1015INTERFACE OPERATOR (<)
1016 MODULE PROCEDURE ttr_mapper_lt
1017END INTERFACE
1018
1019INTERFACE OPERATOR (>=)
1020 MODULE PROCEDURE ttr_mapper_ge
1021END INTERFACE
1022
1023INTERFACE OPERATOR (<=)
1024 MODULE PROCEDURE ttr_mapper_le
1025END INTERFACE
1026
1027#undef VOL7D_POLY_TYPE
1028#undef VOL7D_POLY_TYPES
1029#undef ENABLE_SORT
1030#define VOL7D_POLY_TYPE TYPE(ttr_mapper)
1031#define VOL7D_POLY_TYPES _ttr_mapper
1032#define ENABLE_SORT
1033#include "array_utilities_pre.F90"
1034
1035#define ARRAYOF_ORIGTYPE TYPE(ttr_mapper)
1036#define ARRAYOF_TYPE arrayof_ttr_mapper
1037#define ARRAYOF_ORIGEQ 1
1038#define ARRAYOF_ORIGGT 1
1039#include "arrayof_pre.F90"
1040! from arrayof
1041
1042
1043CONTAINS
1044
1045! simplified comparisons only on time
1046ELEMENTAL FUNCTION ttr_mapper_eq(this, that) RESULT(res)
1047TYPE(ttr_mapper),INTENT(IN) :: this, that
1048LOGICAL :: res
1049
1050res = this%time == that%time
1051
1052END FUNCTION ttr_mapper_eq
1053
1054ELEMENTAL FUNCTION ttr_mapper_ne(this, that) RESULT(res)
1055TYPE(ttr_mapper),INTENT(IN) :: this, that
1056LOGICAL :: res
1057
1058res = this%time /= that%time
1059
1060END FUNCTION ttr_mapper_ne
1061
1062ELEMENTAL FUNCTION ttr_mapper_gt(this, that) RESULT(res)
1063TYPE(ttr_mapper),INTENT(IN) :: this, that
1064LOGICAL :: res
1065
1066res = this%time > that%time
1067
1068END FUNCTION ttr_mapper_gt
1069
1070ELEMENTAL FUNCTION ttr_mapper_lt(this, that) RESULT(res)
1071TYPE(ttr_mapper),INTENT(IN) :: this, that
1072LOGICAL :: res
1073
1074res = this%time < that%time
1075
1076END FUNCTION ttr_mapper_lt
1077
1078ELEMENTAL FUNCTION ttr_mapper_ge(this, that) RESULT(res)
1079TYPE(ttr_mapper),INTENT(IN) :: this, that
1080LOGICAL :: res
1081
1082res = this%time >= that%time
1083
1084END FUNCTION ttr_mapper_ge
1085
1086ELEMENTAL FUNCTION ttr_mapper_le(this, that) RESULT(res)
1087TYPE(ttr_mapper),INTENT(IN) :: this, that
1088LOGICAL :: res
1089
1090res = this%time <= that%time
1091
1092END FUNCTION ttr_mapper_le
1093
1094#include "arrayof_post.F90"
1095#include "array_utilities_inc.F90"
1096
1097
1098! common operations for statistical processing by differences
1099SUBROUTINE recompute_stat_proc_diff_common(itime, itimerange, stat_proc, step, &
1100 otime, otimerange, map_tr, f, keep_tr, time_definition, full_steps, &
1101 start)
1102TYPE(datetime),INTENT(in) :: itime(:)
1103TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1104INTEGER,INTENT(in) :: stat_proc
1105TYPE(timedelta),INTENT(in) :: step
1106TYPE(datetime),POINTER :: otime(:)
1107TYPE(vol7d_timerange),POINTER :: otimerange(:)
1108INTEGER,ALLOCATABLE,INTENT(out) :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
1109INTEGER :: nitr
1110LOGICAL,ALLOCATABLE :: mask_timerange(:)
1111INTEGER,INTENT(in) :: time_definition
1112LOGICAL,INTENT(in),OPTIONAL :: full_steps
1113TYPE(datetime),INTENT(in),OPTIONAL :: start
1114
1115INTEGER :: i, j, k, l, dirtyrep
1116INTEGER :: steps
1117LOGICAL :: lfull_steps, useful
1118TYPE(datetime) :: lstart, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1119TYPE(vol7d_timerange) :: tmptimerange
1120TYPE(arrayof_datetime) :: a_otime
1121TYPE(arrayof_vol7d_timerange) :: a_otimerange
1122TYPE(timedelta) :: start_delta
1123
1124! compute length of cumulation step in seconds
1125CALL getval(step, asec=steps)
1126
1127lstart = datetime_miss
1128IF (PRESENT(start)) lstart = start
1129lfull_steps = optio_log(full_steps)
1130
1131! create a mask of suitable timeranges
1132ALLOCATE(mask_timerange(SIZE(itimerange)))
1133mask_timerange(:) = itimerange(:)%timerange == stat_proc &
1134 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p2 /= imiss &
1135 .AND. itimerange(:)%p1 >= 0 &
1136 .AND. itimerange(:)%p2 > 0
1137
1138IF (lfull_steps .AND. steps /= 0) THEN ! keep only timeranges defining intervals ending at integer forecast steps or analysis timeranges
1139 mask_timerange(:) = mask_timerange(:) .AND. &
1140 (itimerange(:)%p1 == 0 .OR. & ! all analyses
1141 mod(itimerange(:)%p1, steps) == 0 .OR. & ! end time is mod step
1142 mod(itimerange(:)%p1 - itimerange(:)%p2, steps) == 0) ! start time is mod step
1143ENDIF
1144! mask_timerange includes all candidate timeranges
1145
1146nitr = count(mask_timerange)
1147ALLOCATE(f(nitr))
1148j = 1
1149DO i = 1, nitr
1150 DO WHILE(.NOT.mask_timerange(j))
1151 j = j + 1
1152 ENDDO
1153 f(i) = j
1154 j = j + 1
1155ENDDO
1156
1157! now we have to evaluate time/timerage pairs which do not need processing
1158ALLOCATE(keep_tr(nitr, SIZE(itime), 2))
1159CALL compute_keep_tr()
1160
1161ALLOCATE(map_tr(nitr, SIZE(itime), nitr, SIZE(itime), 2))
1162map_tr(:,:,:,:,:) = imiss
1163
1164! scan through all possible combinations of time and timerange
1165DO dirtyrep = 1, 2
1166 IF (dirtyrep == 2) THEN ! dirty and expensive trick for sorting descriptors
1167 CALL packarray(a_otime)
1168 CALL packarray(a_otimerange)
1169 CALL sort(a_otime%array)
1170 CALL sort(a_otimerange%array)
1171 ENDIF
1172 DO l = 1, SIZE(itime)
1173 DO k = 1, nitr
1174 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1175 time_definition, pstart2, pend2, reftime2)
1176
1177 DO j = 1, SIZE(itime)
1178 DO i = 1, nitr
1179 useful = .false.
1180 CALL time_timerange_get_period(itime(j), itimerange(f(i)), &
1181 time_definition, pstart1, pend1, reftime1)
1182 tmptimerange = vol7d_timerange_new(timerange=stat_proc)
1183
1184 IF (reftime2 == pend2 .AND. reftime1 == pend1) THEN ! analysis
1185 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! =-|
1186 CALL time_timerange_set_period(tmptime, tmptimerange, &
1187 time_definition, pend1, pend2, reftime2)
1188 IF (lfull_steps) THEN
1189 IF (mod(reftime2, step) == timedelta_0) THEN
1190 useful = .true.
1191 ENDIF
1192 ELSE
1193 useful = .true.
1194 ENDIF
1195
1196 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! -=|
1197 CALL time_timerange_set_period(tmptime, tmptimerange, &
1198 time_definition, pstart2, pstart1, pstart1)
1199 IF (lfull_steps) THEN
1200 IF (mod(pstart1, step) == timedelta_0) THEN
1201 useful = .true.
1202 ENDIF
1203 ELSE
1204 useful = .true.
1205 ENDIF
1206 ENDIF
1207
1208 ELSE IF (reftime2 == reftime1) THEN ! forecast, same reftime
1209 IF (lfull_steps) THEN
1210 IF (c_e(lstart)) THEN
1211! lstart shifts the interval for computing modulo step, this does not
1212! remove data before lstart but just shifts the phase
1213 start_delta = lstart-reftime2
1214 ELSE
1215 start_delta = timedelta_0
1216 ENDIF
1217 ENDIF
1218
1219 IF (pstart2 == pstart1 .AND. pend2 > pend1) THEN ! |=-
1220 CALL time_timerange_set_period(tmptime, tmptimerange, &
1221 time_definition, pend1, pend2, reftime2)
1222 IF (lfull_steps) THEN
1223 IF (mod(pend2-reftime2-start_delta, step) == timedelta_0) THEN
1224 useful = .true.
1225 ENDIF
1226 ELSE
1227 useful = .true.
1228 ENDIF
1229! keep only data after lstart
1230 IF (c_e(lstart)) THEN
1231 IF (lstart > pend1) useful = .false.
1232 ENDIF
1233
1234 ELSE IF (pstart2 < pstart1 .AND. pend2 == pend1) THEN ! |-=
1235 CALL time_timerange_set_period(tmptime, tmptimerange, &
1236 time_definition, pstart2, pstart1, reftime2)
1237 IF (lfull_steps) THEN
1238 IF (mod(pstart1-reftime2-start_delta, step) == timedelta_0) THEN
1239 useful = .true.
1240 ENDIF
1241 ELSE
1242 useful = .true.
1243 ENDIF
1244! keep only data after lstart
1245 IF (c_e(lstart)) THEN
1246 IF (lstart > pstart2) useful = .false.
1247 ENDIF
1248
1249 ENDIF
1250 ENDIF
1251 useful = useful .AND. tmptime /= datetime_miss .AND. &
1252 tmptimerange /= vol7d_timerange_miss .AND. tmptimerange%p2 == steps
1253
1254 IF (useful) THEN ! add a_otime, a_otimerange
1255 map_tr(i,j,k,l,1) = append_unique(a_otime, tmptime)
1256 map_tr(i,j,k,l,2) = append_unique(a_otimerange, tmptimerange)
1257 ENDIF
1258 ENDDO
1259 ENDDO
1260 ENDDO
1261 ENDDO
1262ENDDO
1263
1264! we have to repeat the computation with sorted arrays
1265CALL compute_keep_tr()
1266
1267otime => a_otime%array
1268otimerange => a_otimerange%array
1269! delete local objects keeping the contents
1270CALL delete(a_otime, nodealloc=.true.)
1271CALL delete(a_otimerange, nodealloc=.true.)
1272
1273#ifdef DEBUG
1274CALL l4f_log(l4f_debug, &
1275 'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr,1)))//', '// &
1276 t2c((SIZE(map_tr,2)))//', '// &
1277 t2c((SIZE(map_tr,3)))//', '// &
1278 t2c((SIZE(map_tr,4))))
1279CALL l4f_log(l4f_debug, &
1280 'recompute_stat_proc_diff, map_tr: '//t2c((SIZE(map_tr))/2)//', '// &
1281 t2c(count(c_e(map_tr))/2))
1282CALL l4f_log(l4f_debug, &
1283 'recompute_stat_proc_diff, nitr: '//t2c(nitr))
1284CALL l4f_log(l4f_debug, &
1285 'recompute_stat_proc_diff, good timeranges: '//t2c(count(c_e(keep_tr))/2))
1286CALL l4f_log(l4f_debug, &
1287 'recompute_stat_proc_diff, output times: '//t2c(SIZE(otime)))
1288CALL l4f_log(l4f_debug, &
1289 'recompute_stat_proc_diff, output timeranges: '//t2c(SIZE(otimerange)))
1290#endif
1291
1292CONTAINS
1293
1294SUBROUTINE compute_keep_tr()
1295INTEGER :: start_deltas
1296
1297keep_tr(:,:,:) = imiss
1298DO l = 1, SIZE(itime)
1299 itrloop: DO k = 1, nitr
1300 IF (itimerange(f(k))%p2 == steps) THEN
1301 CALL time_timerange_get_period(itime(l), itimerange(f(k)), &
1302 time_definition, pstart2, pend2, reftime2)
1303 useful = .false.
1304! keep only data after lstart
1305 IF (c_e(lstart)) THEN
1306 IF (lstart > pstart2) cycle itrloop
1307 ENDIF
1308 IF (reftime2 == pend2) THEN ! analysis
1309 IF (c_e(lstart)) THEN ! in analysis mode start wins over full_steps
1310 IF (mod(reftime2-lstart, step) == timedelta_0) THEN
1311 useful = .true.
1312 ENDIF
1313 ELSE IF (lfull_steps) THEN
1314 IF (mod(reftime2, step) == timedelta_0) THEN
1315 useful = .true.
1316 ENDIF
1317 ELSE
1318 useful = .true.
1319 ENDIF
1320 ELSE ! forecast
1321 IF (lfull_steps) THEN
1322! same as above for start_delta, but in seconds and not in timerange/timedelta units
1323 IF (c_e(lstart)) THEN
1324 start_deltas = timedelta_getamsec(lstart-reftime2)/1000_int_ll
1325 ELSE
1326 start_deltas = 0
1327 ENDIF
1328 IF (mod(itimerange(f(k))%p1 - start_deltas, steps) == 0) THEN
1329 useful = .true.
1330 ENDIF
1331 ELSE
1332 useful = .true.
1333 ENDIF
1334 ENDIF
1335 IF (useful) THEN
1336! CALL time_timerange_set_period(tmptime, tmptimerange, &
1337! time_definition, pstart2, pend2, reftime2)
1338 keep_tr(k,l,1) = append_unique(a_otime, itime(l))
1339 keep_tr(k,l,2) = append_unique(a_otimerange, itimerange(f(k)))
1340 ENDIF
1341 ENDIF
1342 ENDDO itrloop
1343ENDDO
1344
1345END SUBROUTINE compute_keep_tr
1346
1347END SUBROUTINE recompute_stat_proc_diff_common
1348
1349
1350! common operations for statistical processing by metamorphosis
1351SUBROUTINE compute_stat_proc_metamorph_common(istat_proc, itimerange, ostat_proc, &
1352 otimerange, map_tr)
1353INTEGER,INTENT(in) :: istat_proc
1354TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1355INTEGER,INTENT(in) :: ostat_proc
1356TYPE(vol7d_timerange),POINTER :: otimerange(:)
1357INTEGER,POINTER :: map_tr(:)
1358
1359INTEGER :: i
1360LOGICAL :: tr_mask(SIZE(itimerange))
1361
1362IF (SIZE(itimerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1363 ALLOCATE(otimerange(0), map_tr(0))
1364 RETURN
1365ENDIF
1366
1367! useful timeranges
1368tr_mask(:) = itimerange(:)%timerange == istat_proc .AND. itimerange(:)%p2 /= imiss &
1369 .AND. itimerange(:)%p2 /= 0 ! .AND. itimerange(:)%p1 == 0
1370ALLOCATE(otimerange(count(tr_mask)), map_tr(count(tr_mask)))
1371
1372otimerange = pack(itimerange, mask=tr_mask)
1373otimerange(:)%timerange = ostat_proc
1374map_tr = pack((/(i,i=1,SIZE(itimerange))/), mask=tr_mask)
1375
1376END SUBROUTINE compute_stat_proc_metamorph_common
1377
1378
1379! common operations for statistical processing by aggregation
1380SUBROUTINE recompute_stat_proc_agg_common(itime, itimerange, stat_proc, tri, &
1381 step, time_definition, otime, otimerange, map_ttr, dtratio, start, full_steps)
1382TYPE(datetime),INTENT(in) :: itime(:)
1383TYPE(vol7d_timerange),INTENT(in) :: itimerange(:)
1384INTEGER,INTENT(in) :: stat_proc
1385INTEGER,INTENT(in) :: tri
1386TYPE(timedelta),INTENT(in) :: step
1387INTEGER,INTENT(in) :: time_definition
1388TYPE(datetime),POINTER :: otime(:)
1389TYPE(vol7d_timerange),POINTER :: otimerange(:)
1390TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
1391INTEGER,POINTER,OPTIONAL :: dtratio(:)
1392TYPE(datetime),INTENT(in),OPTIONAL :: start
1393LOGICAL,INTENT(in),OPTIONAL :: full_steps
1394
1395INTEGER :: i, j, k, l, na, nf, n
1396INTEGER :: steps, p1, maxp1, maxp2, minp1mp2, dstart, msteps
1397INTEGER(kind=int_ll) :: stepms, mstepms
1398LOGICAL :: lforecast
1399TYPE(datetime) :: lstart, lend, pstart1, pstart2, pend1, pend2, reftime1, reftime2, tmptime
1400TYPE(arrayof_datetime) :: a_otime
1401TYPE(arrayof_vol7d_timerange) :: a_otimerange
1402TYPE(arrayof_integer) :: a_dtratio
1403LOGICAL,ALLOCATABLE :: mask_timerange(:) ! improve !!!!
1404TYPE(ttr_mapper) :: lmapper
1405CHARACTER(len=8) :: env_var
1406LOGICAL :: climat_behavior
1407
1408
1409! enable bad behavior for climat database (only for instantaneous data)
1410env_var = ''
1411CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
1412climat_behavior = len_trim(env_var) > 0 .AND. .NOT.PRESENT(dtratio)
1413
1414! compute length of cumulation step in seconds
1415CALL getval(timedelta_depop(step), asec=steps)
1416
1417! create a mask of suitable timeranges
1418ALLOCATE(mask_timerange(SIZE(itimerange)))
1419mask_timerange(:) = itimerange(:)%timerange == tri &
1420 .AND. itimerange(:)%p1 /= imiss .AND. itimerange(:)%p1 >= 0
1421
1422IF (PRESENT(dtratio)) THEN
1423 WHERE(itimerange(:)%p2 > 0 .AND. itimerange(:)%p2 /= imiss) ! for avoiding FPE MOD(steps, 0)
1424 mask_timerange(:) = mask_timerange(:) .AND. mod(steps, itimerange(:)%p2) == 0
1425 ELSEWHERE
1426 mask_timerange(:) = .false.
1427 END WHERE
1428ELSE ! instantaneous
1429 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p2 == 0
1430ENDIF
1431
1432#ifdef DEBUG
1433CALL l4f_log(l4f_debug, &
1434 '(re)compute_stat_proc_agg, number of useful timeranges before choosing analysis/forecast: '// &
1435 t2c(count(mask_timerange)))
1436#endif
1437
1438! euristically determine whether we are dealing with an
1439! analysis/observation or a forecast dataset
1440na = count(mask_timerange(:) .AND. itimerange(:)%p1 == 0)
1441nf = count(mask_timerange(:) .AND. itimerange(:)%p1 > 0)
1442lforecast = nf >= na
1443#ifdef DEBUG
1444CALL l4f_log(l4f_debug, &
1445 'recompute_stat_proc_agg, na: '//t2c(na)//', nf: '//t2c(nf))
1446#endif
1447
1448! keep only timeranges of one type (really necessary?)
1449IF (lforecast) THEN
1450 CALL l4f_log(l4f_info, &
1451 'recompute_stat_proc_agg, processing in forecast mode')
1452ELSE
1453 mask_timerange(:) = mask_timerange(:) .AND. itimerange(:)%p1 == 0
1454 CALL l4f_log(l4f_info, &
1455 'recompute_stat_proc_agg, processing in analysis mode')
1456ENDIF
1457
1458#ifdef DEBUG
1459CALL l4f_log(l4f_debug, &
1460 '(re)compute_stat_proc_agg, number of useful timeranges: '// &
1461 t2c(count(mask_timerange)))
1462#endif
1463
1464IF (SIZE(itime) == 0 .OR. count(mask_timerange) == 0) THEN ! avoid segmentation fault in case of empty volume
1465 ALLOCATE(otime(0), otimerange(0), map_ttr(0,0))
1466 IF (PRESENT(dtratio)) ALLOCATE(dtratio(0))
1467 RETURN
1468ENDIF
1469
1470! determine start and end of processing period, should work with p2==0
1471lstart = datetime_miss
1472IF (PRESENT(start)) lstart = start
1473lend = itime(SIZE(itime))
1474! compute some quantities
1475maxp1 = maxval(itimerange(:)%p1, mask=mask_timerange)
1476maxp2 = maxval(itimerange(:)%p2, mask=mask_timerange)
1477minp1mp2 = minval(itimerange(:)%p1 - itimerange(:)%p2, mask=mask_timerange)
1478IF (time_definition == 0) THEN ! reference time
1479 lend = lend + timedelta_new(sec=maxp1)
1480ENDIF
1481! extend interval at the end in order to include all the data in case
1482! frac_valid<1; must use < and not <= in the DO WHILE loops some lines
1483! below in order to exclude the last full interval which would be empty
1484lend = lend + step
1485IF (lstart == datetime_miss) THEN ! autodetect
1486 lstart = itime(1)
1487! if autodetected, adjust to obtain real absolute start of data
1488 IF (time_definition == 0) THEN ! reference time
1489 lstart = lstart + timedelta_new(sec=minp1mp2)
1490 ELSE ! verification time
1491! go back to start of longest processing interval
1492 lstart = lstart - timedelta_new(sec=maxp2)
1493 ENDIF
1494! apply full_steps in analysis mode and when start is not specified
1495! (start by itself in analysis mode implies full_steps with respect to
1496! start instead of absolute full steps), todo in forecast mode
1497 IF (optio_log(full_steps) .AND. .NOT.lforecast) THEN
1498 lstart = lstart - (mod(lstart, step)) ! round to step, (should be MODULO, not MOD)
1499 ENDIF
1500ENDIF
1501
1502#ifdef DEBUG
1503CALL l4f_log(l4f_debug, &
1504 'recompute_stat_proc_agg, processing period: '//t2c(lstart)//' - '//t2c(lend))
1505#endif
1506
1507! create output time and timerange lists
1508
1509IF (lforecast) THEN ! forecast mode
1510 IF (time_definition == 0) THEN ! reference time
1511 CALL insert(a_otime, itime) ! should I limit to elements itime >= lstart?
1512
1513! apply start shift to timerange, not to reftime
1514! why did we use itime(SIZE(itime)) (last ref time)?
1515! CALL getval(lstart-itime(SIZE(itime)), asec=dstart)
1516 CALL getval(lstart-itime(1), asec=dstart)
1517! allow starting before first reftime but restrict dtstart to -steps
1518! dstart = MAX(0, dstart)
1519 IF (dstart < 0) dstart = mod(dstart, steps)
1520 DO p1 = steps + dstart, maxp1, steps
1521 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1522 ENDDO
1523
1524 ELSE ! verification time
1525
1526! verification time in forecast mode is the ugliest case, because we
1527! don't know a priori how many different (thus incompatible) reference
1528! times we have, so some assumption of regularity has to be made. For
1529! this reason msteps, the minimum step between two times, is
1530! computed. We choose to compute it as a difference between itime
1531! elements, it could be also computed as difference of itimerange%p1
1532! elements. But what if it is not modulo steps?
1533 mstepms = steps*1000_int_ll
1534 DO i = 2, SIZE(itime)
1535 CALL getval(itime(i)-itime(i-1), amsec=stepms)
1536 IF (stepms > 0_int_ll .AND. stepms < mstepms) THEN
1537 msteps = stepms/1000_int_ll
1538 IF (mod(steps, msteps) == 0) mstepms = stepms
1539 ENDIF
1540 ENDDO
1541 msteps = mstepms/1000_int_ll
1542
1543 tmptime = lstart + step
1544 DO WHILE(tmptime < lend) ! < since lend has been extended
1545 CALL insert_unique(a_otime, tmptime)
1546 tmptime = tmptime + step
1547 ENDDO
1548
1549! in next loop, we used initially steps instead of msteps (see comment
1550! above), this gave much less combinations of time/timeranges with
1551! possible empty output; we start from msteps instead of steps in
1552! order to include partial initial intervals in case frac_valid<1;
1553! however, as a gemeral rule, for aggregation of forecasts it's better
1554! to use reference time
1555 DO p1 = msteps, maxp1, msteps
1556 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, p1, steps))
1557 ENDDO
1558
1559 ENDIF
1560
1561ELSE ! analysis mode
1562 tmptime = lstart + step
1563 DO WHILE(tmptime < lend) ! < since lend has been extended
1564 CALL insert_unique(a_otime, tmptime)
1565 tmptime = tmptime + step
1566 ENDDO
1567 CALL insert_unique(a_otimerange, vol7d_timerange_new(stat_proc, 0, steps))
1568
1569ENDIF
1570
1571CALL packarray(a_otime)
1572CALL packarray(a_otimerange)
1573otime => a_otime%array
1574otimerange => a_otimerange%array
1575CALL sort(otime)
1576CALL sort(otimerange)
1577! delete local objects keeping the contents
1578CALL delete(a_otime, nodealloc=.true.)
1579CALL delete(a_otimerange, nodealloc=.true.)
1580
1581#ifdef DEBUG
1582CALL l4f_log(l4f_debug, &
1583 'recompute_stat_proc_agg, output time and timerange: '//&
1584 t2c(SIZE(otime))//', '//t2c(size(otimerange)))
1585#endif
1586
1587IF (PRESENT(dtratio)) THEN
1588! count the possible i/o interval ratios
1589 DO k = 1, SIZE(itimerange)
1590 IF (itimerange(k)%p2 /= 0) &
1591 CALL insert_unique(a_dtratio, steps/itimerange(k)%p2) ! guaranteed to be integer
1592 ENDDO
1593 CALL packarray(a_dtratio)
1594 dtratio => a_dtratio%array
1595 CALL sort(dtratio)
1596! delete local object keeping the contents
1597 CALL delete(a_dtratio, nodealloc=.true.)
1598
1599#ifdef DEBUG
1600 CALL l4f_log(l4f_debug, &
1601 'recompute_stat_proc_agg, found '//t2c(size(dtratio))// &
1602 ' possible aggregation ratios, from '// &
1603 t2c(dtratio(1))//' to '//t2c(dtratio(SIZE(dtratio))))
1604#endif
1605
1606 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1607 do_itimerange1: DO l = 1, SIZE(itimerange)
1608 IF (.NOT.mask_timerange(l)) cycle do_itimerange1
1609 do_itime1: DO k = 1, SIZE(itime)
1610 CALL time_timerange_get_period(itime(k), itimerange(l), &
1611 time_definition, pstart1, pend1, reftime1)
1612 do_otimerange1: DO j = 1, SIZE(otimerange)
1613 do_otime1: DO i = 1, SIZE(otime)
1614 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1615 time_definition, pstart2, pend2, reftime2)
1616 IF (lforecast) THEN
1617 IF (reftime1 /= reftime2) cycle do_otime1
1618 ENDIF
1619
1620 IF (pstart1 >= pstart2 .AND. pend1 <= pend2 .AND. &
1621 mod(pstart1-pstart2, pend1-pstart1) == timedelta_0) THEN ! useful
1622 lmapper%it = k
1623 lmapper%itr = l
1624 lmapper%extra_info = steps/itimerange(l)%p2 ! dtratio, guaranteed to be integer
1625 n = append(map_ttr(i,j), lmapper)
1626 cycle do_itime1 ! can contribute only to a single interval
1627 ENDIF
1628 ENDDO do_otime1
1629 ENDDO do_otimerange1
1630 ENDDO do_itime1
1631 ENDDO do_itimerange1
1632
1633ELSE
1634
1635 ALLOCATE(map_ttr(SIZE(otime),SIZE(otimerange)))
1636 do_itimerange2: DO l = 1, SIZE(itimerange)
1637 IF (.NOT.mask_timerange(l)) cycle do_itimerange2
1638 do_itime2: DO k = 1, SIZE(itime)
1639 CALL time_timerange_get_period(itime(k), itimerange(l), &
1640 time_definition, pstart1, pend1, reftime1)
1641 do_otimerange2: DO j = 1, SIZE(otimerange)
1642 do_otime2: DO i = 1, SIZE(otime)
1643 CALL time_timerange_get_period_pop(otime(i), otimerange(j), step, &
1644 time_definition, pstart2, pend2, reftime2)
1645 IF (lforecast) THEN
1646 IF (reftime1 /= reftime2) cycle do_otime2
1647 ENDIF
1648
1649 IF (climat_behavior .AND. pstart1 == pstart2) cycle do_otime2
1650 IF (pstart1 >= pstart2 .AND. pend1 <= pend2) THEN ! useful
1651 lmapper%it = k
1652 lmapper%itr = l
1653 IF (pstart1 == pstart2) THEN
1654 lmapper%extra_info = 1 ! start of interval
1655 ELSE IF (pend1 == pend2) THEN
1656 lmapper%extra_info = 2 ! end of interval
1657 ELSE
1658 lmapper%extra_info = imiss
1659 ENDIF
1660 lmapper%time = pstart1 ! = pend1, order by time?
1661 n = insert_sorted(map_ttr(i,j), lmapper, .true., .true.)
1662! no CYCLE, a single input can contribute to multiple output intervals
1663 ENDIF
1664 ENDDO do_otime2
1665 ENDDO do_otimerange2
1666 ENDDO do_itime2
1667 ENDDO do_itimerange2
1668
1669ENDIF
1670
1671END SUBROUTINE recompute_stat_proc_agg_common
1672
1673
1674SUBROUTINE compute_stat_proc_agg_sw(vertime, pstart, pend, time_mask, &
1675 max_step, weights)
1676TYPE(datetime),INTENT(in) :: vertime(:)
1677TYPE(datetime),INTENT(in) :: pstart
1678TYPE(datetime),INTENT(in) :: pend
1679LOGICAL,INTENT(in) :: time_mask(:)
1680TYPE(timedelta),OPTIONAL,INTENT(out) :: max_step
1681DOUBLE PRECISION,OPTIONAL,INTENT(out) :: weights(:)
1682
1683INTEGER :: i, nt
1684TYPE(datetime),ALLOCATABLE :: lvertime(:)
1685TYPE(datetime) :: half, nexthalf
1686INTEGER(kind=int_ll) :: dt, tdt
1687
1688nt = count(time_mask)
1689ALLOCATE(lvertime(nt))
1690lvertime = pack(vertime, mask=time_mask)
1691
1692IF (PRESENT(max_step)) THEN
1693! new way
1694! max_step = timedelta_0
1695! DO i = 1, nt - 1
1696! IF (lvertime(i+1) - lvertime(i) > max_step) &
1697! max_step = lvertime(i+1) - lvertime(i)
1698! ENDDO
1699! IF (lvertime(1) - pstart > max_step) max_step = lvertime(1) - pstart
1700! IF (pend - lvertime(nt) > max_step) max_step = pend - lvertime(nt)
1701! old way
1702 IF (nt == 1) THEN
1703 max_step = pend - pstart
1704 ELSE
1705 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
1706 max_step = half - pstart
1707 DO i = 2, nt - 1
1708 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
1709 IF (nexthalf - half > max_step) max_step = nexthalf - half
1710 half = nexthalf
1711 ENDDO
1712 IF (pend - half > max_step) max_step = pend - half
1713 ENDIF
1714
1715ENDIF
1716
1717IF (PRESENT(weights)) THEN
1718 IF (nt == 1) THEN
1719 weights(1) = 1.0
1720 ELSE
1721 CALL getval(pend - pstart, amsec=tdt)
1722 half = lvertime(1) + (lvertime(2) - lvertime(1))/2
1723 CALL getval(half - pstart, amsec=dt)
1724 weights(1) = dble(dt)/dble(tdt)
1725 DO i = 2, nt - 1
1726 nexthalf = lvertime(i) + (lvertime(i+1) - lvertime(i))/2
1727 CALL getval(nexthalf - half, amsec=dt)
1728 weights(i) = dble(dt)/dble(tdt)
1729 half = nexthalf
1730 ENDDO
1731 CALL getval(pend - half, amsec=dt)
1732 weights(nt) = dble(dt)/dble(tdt)
1733 ENDIF
1734ENDIF
1735
1736END SUBROUTINE compute_stat_proc_agg_sw
1737
1738! get start of period, end of period and reference time from time,
1739! timerange, according to time_definition.
1740SUBROUTINE time_timerange_get_period(time, timerange, time_definition, &
1741 pstart, pend, reftime)
1742TYPE(datetime),INTENT(in) :: time
1743TYPE(vol7d_timerange),INTENT(in) :: timerange
1744INTEGER,INTENT(in) :: time_definition
1745TYPE(datetime),INTENT(out) :: reftime
1746TYPE(datetime),INTENT(out) :: pstart
1747TYPE(datetime),INTENT(out) :: pend
1748
1749TYPE(timedelta) :: p1, p2
1750
1751
1752p1 = timedelta_new(sec=timerange%p1) ! end of period
1753p2 = timedelta_new(sec=timerange%p2) ! length of period
1754
1755IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
1756! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
1757 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
1758 pstart = datetime_miss
1759 pend = datetime_miss
1760 reftime = datetime_miss
1761 RETURN
1762ENDIF
1763
1764IF (time_definition == 0) THEN ! time == reference time
1765 reftime = time
1766 pend = time + p1
1767 pstart = pend - p2
1768ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time
1769 pend = time
1770 pstart = time - p2
1771 reftime = time - p1
1772ELSE
1773 pstart = datetime_miss
1774 pend = datetime_miss
1775 reftime = datetime_miss
1776ENDIF
1777
1778END SUBROUTINE time_timerange_get_period
1779
1780
1781! get start of period, end of period and reference time from time,
1782! timerange, according to time_definition. step is taken separately
1783! from timerange and can be "popular"
1784SUBROUTINE time_timerange_get_period_pop(time, timerange, step, time_definition, &
1785 pstart, pend, reftime)
1786TYPE(datetime),INTENT(in) :: time
1787TYPE(vol7d_timerange),INTENT(in) :: timerange
1788TYPE(timedelta),INTENT(in) :: step
1789INTEGER,INTENT(in) :: time_definition
1790TYPE(datetime),INTENT(out) :: reftime
1791TYPE(datetime),INTENT(out) :: pstart
1792TYPE(datetime),INTENT(out) :: pend
1793
1794TYPE(timedelta) :: p1
1795
1796
1797p1 = timedelta_new(sec=timerange%p1) ! end of period
1798
1799IF (time == datetime_miss .OR. .NOT.c_e(timerange%p1) .OR. .NOT.c_e(timerange%p2) .OR. &
1800! (timerange%p1 > 0 .AND. timerange%p1 < timerange%p2) .OR. &
1801 timerange%p1 < 0 .OR. timerange%p2 < 0) THEN ! is this too pedantic and slow?
1802 pstart = datetime_miss
1803 pend = datetime_miss
1804 reftime = datetime_miss
1805 RETURN
1806ENDIF
1807
1808IF (time_definition == 0) THEN ! time == reference time
1809 reftime = time
1810 pend = time + p1
1811 pstart = pend - step
1812ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time
1813 pend = time
1814 pstart = time - step
1815 reftime = time - p1
1816ELSE
1817 pstart = datetime_miss
1818 pend = datetime_miss
1819 reftime = datetime_miss
1820ENDIF
1821
1822END SUBROUTINE time_timerange_get_period_pop
1823
1824
1825! set time, timerange%p1, timerange%p2 according to pstart, pend,
1826! reftime and time_definition.
1827SUBROUTINE time_timerange_set_period(time, timerange, time_definition, &
1828 pstart, pend, reftime)
1829TYPE(datetime),INTENT(out) :: time
1830TYPE(vol7d_timerange),INTENT(inout) :: timerange
1831INTEGER,INTENT(in) :: time_definition
1832TYPE(datetime),INTENT(in) :: reftime
1833TYPE(datetime),INTENT(in) :: pstart
1834TYPE(datetime),INTENT(in) :: pend
1835
1836TYPE(timedelta) :: p1, p2
1837INTEGER(kind=int_ll) :: dmsec
1838
1839
1840IF (time_definition == 0) THEN ! time == reference time
1841 time = reftime
1842 p1 = pend - reftime
1843 p2 = pend - pstart
1844ELSE IF (time_definition == 1 .OR. time_definition == 2) THEN ! time == verification time
1845 time = pend
1846 p1 = pend - reftime
1847 p2 = pend - pstart
1848ELSE
1849 time = datetime_miss
1850ENDIF
1851
1852IF (time /= datetime_miss) THEN
1853 CALL getval(p1, amsec=dmsec) ! end of period
1854 timerange%p1 = int(dmsec/1000_int_ll)
1855 CALL getval(p2, amsec=dmsec) ! length of period
1856 timerange%p2 = int(dmsec/1000_int_ll)
1857ELSE
1858 timerange%p1 = imiss
1859 timerange%p2 = imiss
1860ENDIF
1861
1862END SUBROUTINE time_timerange_set_period
1863
1864
1865END 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.