libsim Versione 7.1.11
|
◆ vol7d_normalize_vcoord()
Metodo per normalizzare la coordinata verticale. Per ora la normalizzazione effettuata riporta i valori di pressione nella sezione dati alla coordinata verticale sostituendo quella eventualmente presente. Classicamente serve per i dati con coordinata verticale model layer (105) Essendo che la pressione varia nello spazio orizzontale e nel tempo questo metodo restituisce un solo profilo verticale.
Definizione alla linea 1748 del file vol7d_class_compute.F90. 1749! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
1750! authors:
1751! Davide Cesari <dcesari@arpa.emr.it>
1752! Paolo Patruno <ppatruno@arpa.emr.it>
1753
1754! This program is free software; you can redistribute it and/or
1755! modify it under the terms of the GNU General Public License as
1756! published by the Free Software Foundation; either version 2 of
1757! the License, or (at your option) any later version.
1758
1759! This program is distributed in the hope that it will be useful,
1760! but WITHOUT ANY WARRANTY; without even the implied warranty of
1761! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1762! GNU General Public License for more details.
1763
1764! You should have received a copy of the GNU General Public License
1765! along with this program. If not, see <http://www.gnu.org/licenses/>.
1766#include "config.h"
1767
1779IMPLICIT NONE
1780
1781CONTAINS
1782
1783
1854SUBROUTINE vol7d_compute_stat_proc(this, that, stat_proc_input, stat_proc, &
1855 step, start, full_steps, frac_valid, max_step, weighted, other)
1856TYPE(vol7d),INTENT(inout) :: this
1857TYPE(vol7d),INTENT(out) :: that
1858INTEGER,INTENT(in) :: stat_proc_input
1859INTEGER,INTENT(in) :: stat_proc
1860TYPE(timedelta),INTENT(in) :: step
1861TYPE(datetime),INTENT(in),OPTIONAL :: start
1862LOGICAL,INTENT(in),OPTIONAL :: full_steps
1863REAL,INTENT(in),OPTIONAL :: frac_valid
1864TYPE(timedelta),INTENT(in),OPTIONAL :: max_step
1865LOGICAL,INTENT(in),OPTIONAL :: weighted
1866TYPE(vol7d),INTENT(inout),OPTIONAL :: other
1867
1868TYPE(vol7d) :: that1, that2, other1
1869INTEGER :: steps
1870
1871IF (stat_proc_input == 254) THEN
1872 CALL l4f_log(l4f_info, 'computing statistical processing by aggregation '//&
1874
1875 CALL vol7d_compute_stat_proc_agg(this, that, stat_proc, &
1876 step, start, full_steps, max_step, weighted, other)
1877
1878ELSE IF (stat_proc == 254) THEN
1879 CALL l4f_log(l4f_info, &
1880 'computing instantaneous data from statistically processed '//&
1882
1883! compute length of cumulation step in seconds
1885
1886 IF (any(this%timerange(:)%p2 == steps)) THEN ! data is ready
1887 CALL vol7d_decompute_stat_proc(this, that, step, other, stat_proc_input)
1888 ELSE
1889 IF (any(this%timerange(:)%p2 == steps/2)) THEN ! need to average
1890! average twice on step interval, with a shift of step/2, check full_steps
1891 CALL vol7d_recompute_stat_proc_agg(this, that1, stat_proc_input, &
1892 step, full_steps=.false., frac_valid=1.0)
1893 CALL vol7d_recompute_stat_proc_agg(this, that2, stat_proc_input, &
1894 step, start=that1%time(1)+step/2, full_steps=.false., frac_valid=1.0)
1895! merge the result
1897! and process it
1898 CALL vol7d_decompute_stat_proc(that1, that, step, other, stat_proc_input)
1901 ELSE
1902! do nothing
1903 ENDIF
1904 ENDIF
1905
1906ELSE IF (stat_proc_input == stat_proc .OR. &
1907 (stat_proc == 0 .OR. stat_proc == 2 .OR. stat_proc == 3)) THEN
1908! avg, min and max can be computed from any input, with care
1909 CALL l4f_log(l4f_info, &
1910 'recomputing statistically processed data by aggregation and difference '//&
1912
1913 IF (PRESENT(other)) THEN
1914 CALL vol7d_recompute_stat_proc_agg(this, that1, stat_proc, &
1915 step, start, full_steps, frac_valid, &
1916 other=other, stat_proc_input=stat_proc_input)
1917 CALL vol7d_recompute_stat_proc_diff(this, that2, stat_proc, &
1918 step, full_steps, start, other=other1)
1920 ELSE
1921 CALL vol7d_recompute_stat_proc_agg(this, that1, stat_proc, &
1922 step, start, full_steps, frac_valid, stat_proc_input=stat_proc_input)
1923 CALL vol7d_recompute_stat_proc_diff(this, that2, stat_proc, step, full_steps, &
1924 start)
1925 ENDIF
1926
1929 that = that1
1930ELSE ! IF (stat_proc_input /= stat_proc) THEN
1931 IF ((stat_proc_input == 0 .AND. stat_proc == 1) .OR. &
1932 (stat_proc_input == 1 .AND. stat_proc == 0)) THEN
1933 CALL l4f_log(l4f_info, &
1934 'computing statistically processed data by integration/differentiation '// &
1936 CALL vol7d_compute_stat_proc_metamorph(this, that, stat_proc_input, &
1937 stat_proc)
1938 ELSE
1939 CALL l4f_log(l4f_error, &
1941 ' not implemented or does not make sense')
1942 CALL raise_error()
1943 ENDIF
1944ENDIF
1945
1946END SUBROUTINE vol7d_compute_stat_proc
1947
1948
1994SUBROUTINE vol7d_recompute_stat_proc_agg(this, that, stat_proc, &
1995 step, start, full_steps, frac_valid, other, stat_proc_input)
1996TYPE(vol7d),INTENT(inout) :: this
1997TYPE(vol7d),INTENT(out) :: that
1998INTEGER,INTENT(in) :: stat_proc
1999TYPE(timedelta),INTENT(in) :: step
2000TYPE(datetime),INTENT(in),OPTIONAL :: start
2001LOGICAL,INTENT(in),OPTIONAL :: full_steps
2002REAL,INTENT(in),OPTIONAL :: frac_valid
2003TYPE(vol7d),INTENT(inout),OPTIONAL :: other
2004INTEGER,INTENT(in),OPTIONAL :: stat_proc_input
2005
2006INTEGER :: tri
2007INTEGER :: i, j, n, n1, ndtr, i1, i3, i5, i6
2008INTEGER :: linshape(1)
2009REAL :: lfrac_valid, frac_c, frac_m
2010LOGICAL,ALLOCATABLE :: ttr_mask(:,:)
2011TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
2012INTEGER,POINTER :: dtratio(:)
2013
2014
2015IF (PRESENT(stat_proc_input)) THEN
2016 tri = stat_proc_input
2017ELSE
2018 tri = stat_proc
2019ENDIF
2020IF (PRESENT(frac_valid)) THEN
2021 lfrac_valid = frac_valid
2022ELSE
2023 lfrac_valid = 1.0
2024ENDIF
2025
2026! be safe
2027CALL vol7d_alloc_vol(this)
2028! initial check
2029
2030! cleanup the original volume
2031CALL vol7d_smart_sort(this, lsort_time=.true.) ! time-ordered volume needed
2033
2035CALL vol7d_alloc(that, nana=SIZE(this%ana), nlevel=SIZE(this%level), &
2036 nnetwork=SIZE(this%network))
2037IF (ASSOCIATED(this%dativar%r)) THEN
2038 CALL vol7d_alloc(that, ndativarr=SIZE(this%dativar%r))
2039 that%dativar%r = this%dativar%r
2040ENDIF
2041IF (ASSOCIATED(this%dativar%d)) THEN
2042 CALL vol7d_alloc(that, ndativard=SIZE(this%dativar%d))
2043 that%dativar%d = this%dativar%d
2044ENDIF
2045that%ana = this%ana
2046that%level = this%level
2047that%network = this%network
2048
2049! compute the output time and timerange and all the required mappings
2050CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, &
2051 step, this%time_definition, that%time, that%timerange, map_ttr, dtratio, &
2052 start, full_steps)
2053CALL vol7d_alloc_vol(that)
2054
2055ALLOCATE(ttr_mask(SIZE(this%time), SIZE(this%timerange)))
2056linshape = (/SIZE(ttr_mask)/)
2057! finally perform computations
2058IF (ASSOCIATED(this%voldatir)) THEN
2059 DO j = 1, SIZE(that%timerange)
2060 DO i = 1, SIZE(that%time)
2061
2062 DO i1 = 1, SIZE(this%ana)
2063 DO i3 = 1, SIZE(this%level)
2064 DO i6 = 1, SIZE(this%network)
2065 DO i5 = 1, SIZE(this%dativar%r)
2066
2067 frac_m = 0.
2068 DO n1 = SIZE(dtratio), 1, -1 ! precedence to longer periods
2069 IF (dtratio(n1) <= 0) cycle ! safety check
2070 ttr_mask = .false.
2071 DO n = 1, map_ttr(i,j)%arraysize
2072 IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN
2074 map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
2075 ttr_mask(map_ttr(i,j)%array(n)%it, &
2076 map_ttr(i,j)%array(n)%itr) = .true.
2077 ENDIF
2078 ENDIF
2079 ENDDO
2080
2081 ndtr = count(ttr_mask)
2082 frac_c = real(ndtr)/real(dtratio(n1))
2083
2084 IF (ndtr > 0 .AND. frac_c >= max(lfrac_valid, frac_m)) THEN
2085 frac_m = frac_c
2086 SELECT CASE(stat_proc)
2087 CASE (0, 200) ! average, vectorial mean
2088 that%voldatir(i1,i,i3,j,i5,i6) = &
2089 sum(this%voldatir(i1,:,i3,:,i5,i6), &
2090 mask=ttr_mask)/ndtr
2091 CASE (1, 4) ! accumulation, difference
2092 that%voldatir(i1,i,i3,j,i5,i6) = &
2093 sum(this%voldatir(i1,:,i3,:,i5,i6), &
2094 mask=ttr_mask)
2095 CASE (2) ! maximum
2096 that%voldatir(i1,i,i3,j,i5,i6) = &
2097 maxval(this%voldatir(i1,:,i3,:,i5,i6), &
2098 mask=ttr_mask)
2099 CASE (3) ! minimum
2100 that%voldatir(i1,i,i3,j,i5,i6) = &
2101 minval(this%voldatir(i1,:,i3,:,i5,i6), &
2102 mask=ttr_mask)
2103 CASE (6) ! stddev
2104 that%voldatir(i1,i,i3,j,i5,i6) = &
2105 stat_stddev( &
2106 reshape(this%voldatir(i1,:,i3,:,i5,i6), shape=linshape), &
2107 mask=reshape(ttr_mask, shape=linshape))
2108 END SELECT
2109 ENDIF
2110
2111 ENDDO ! dtratio
2112 ENDDO ! var
2113 ENDDO ! network
2114 ENDDO ! level
2115 ENDDO ! ana
2117 ENDDO ! otime
2118 ENDDO ! otimerange
2119ENDIF
2120
2121IF (ASSOCIATED(this%voldatid)) THEN
2122 DO j = 1, SIZE(that%timerange)
2123 DO i = 1, SIZE(that%time)
2124
2125 DO i1 = 1, SIZE(this%ana)
2126 DO i3 = 1, SIZE(this%level)
2127 DO i6 = 1, SIZE(this%network)
2128 DO i5 = 1, SIZE(this%dativar%d)
2129
2130 frac_m = 0.
2131 DO n1 = SIZE(dtratio), 1, -1 ! precedence to longer periods
2132 IF (dtratio(n1) <= 0) cycle ! safety check
2133 ttr_mask = .false.
2134 DO n = 1, map_ttr(i,j)%arraysize
2135 IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN
2137 map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
2138 ttr_mask(map_ttr(i,j)%array(n)%it, &
2139 map_ttr(i,j)%array(n)%itr) = .true.
2140 ENDIF
2141 ENDIF
2142 ENDDO
2143
2144 ndtr = count(ttr_mask)
2145 frac_c = real(ndtr)/real(dtratio(n1))
2146
2147 IF (ndtr > 0 .AND. frac_c >= max(lfrac_valid, frac_m)) THEN
2148 frac_m = frac_c
2149 SELECT CASE(stat_proc)
2150 CASE (0) ! average
2151 that%voldatid(i1,i,i3,j,i5,i6) = &
2152 sum(this%voldatid(i1,:,i3,:,i5,i6), &
2153 mask=ttr_mask)/ndtr
2154 CASE (1, 4) ! accumulation, difference
2155 that%voldatid(i1,i,i3,j,i5,i6) = &
2156 sum(this%voldatid(i1,:,i3,:,i5,i6), &
2157 mask=ttr_mask)
2158 CASE (2) ! maximum
2159 that%voldatid(i1,i,i3,j,i5,i6) = &
2160 maxval(this%voldatid(i1,:,i3,:,i5,i6), &
2161 mask=ttr_mask)
2162 CASE (3) ! minimum
2163 that%voldatid(i1,i,i3,j,i5,i6) = &
2164 minval(this%voldatid(i1,:,i3,:,i5,i6), &
2165 mask=ttr_mask)
2166 CASE (6) ! stddev
2167 that%voldatid(i1,i,i3,j,i5,i6) = &
2168 stat_stddev( &
2169 reshape(this%voldatid(i1,:,i3,:,i5,i6), shape=linshape), &
2170 mask=reshape(ttr_mask, shape=linshape))
2171 END SELECT
2172 ENDIF
2173
2174 ENDDO ! dtratio
2175 ENDDO ! var
2176 ENDDO ! network
2177 ENDDO ! level
2178 ENDDO ! ana
2180 ENDDO ! otime
2181 ENDDO ! otimerange
2182ENDIF
2183
2184DEALLOCATE(ttr_mask)
2185
2186CALL makeother()
2187
2188CONTAINS
2189
2190SUBROUTINE makeother()
2191IF (PRESENT(other)) THEN ! create volume with the remaining data for further processing
2193 ltimerange=(this%timerange(:)%timerange /= tri .OR. this%timerange(:)%p2 == imiss &
2194 .OR. this%timerange(:)%p2 == 0)) ! or MOD(steps, this%timerange(:)%p2) == 0
2195ENDIF
2196END SUBROUTINE makeother
2197
2198END SUBROUTINE vol7d_recompute_stat_proc_agg
2199
2200
2232SUBROUTINE vol7d_compute_stat_proc_agg(this, that, stat_proc, &
2233 step, start, full_steps, max_step, weighted, other)
2234TYPE(vol7d),INTENT(inout) :: this
2235TYPE(vol7d),INTENT(out) :: that
2236INTEGER,INTENT(in) :: stat_proc
2237TYPE(timedelta),INTENT(in) :: step
2238TYPE(datetime),INTENT(in),OPTIONAL :: start
2239LOGICAL,INTENT(in),OPTIONAL :: full_steps
2240TYPE(timedelta),INTENT(in),OPTIONAL :: max_step
2241LOGICAL,INTENT(in),OPTIONAL :: weighted
2242TYPE(vol7d),INTENT(inout),OPTIONAL :: other
2243!INTEGER,INTENT(in),OPTIONAL :: stat_proc_input !< to be used with care, type of statistical processing of data that has to be processed (from grib2 table), only data having timerange of this type will be recomputed, the actual statistical processing performed and which will appear in the output volume, is however determined by \a stat_proc argument
2244
2245TYPE(vol7d) :: v7dtmp
2246INTEGER :: tri
2247INTEGER :: i, j, n, ninp, ndtr, i1, i3, i5, i6, vartype, maxsize
2248TYPE(timedelta) :: lmax_step, act_max_step
2249TYPE(datetime) :: pstart, pend, reftime
2250TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
2251REAL,ALLOCATABLE :: tmpvolr(:)
2252DOUBLE PRECISION,ALLOCATABLE :: tmpvold(:), weights(:)
2253LOGICAL,ALLOCATABLE :: lin_mask(:)
2254LOGICAL :: lweighted
2255CHARACTER(len=8) :: env_var
2256
2257IF (PRESENT(max_step)) THEN
2258 lmax_step = max_step
2259ELSE
2260 lmax_step = timedelta_max
2261ENDIF
2262lweighted = optio_log(weighted)
2263tri = 254
2264! enable bad behavior for climat database
2265env_var = ''
2266CALL getenv('LIBSIM_CLIMAT_BEHAVIOR', env_var)
2267lweighted = lweighted .AND. len_trim(env_var) == 0
2268! only average can be weighted
2269lweighted = lweighted .AND. stat_proc == 0
2270
2271! be safe
2272CALL vol7d_alloc_vol(this)
2273! initial check
2274
2275! cleanup the original volume
2276CALL vol7d_smart_sort(this, lsort_time=.true.) ! time-ordered volume needed
2278! copy everything except time and timerange
2279CALL vol7d_copy(this, v7dtmp, ltime=(/.false./), ltimerange=(/.false./))
2280
2281! create new volume
2283! compute the output time and timerange and all the required mappings
2284CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, &
2285 step, this%time_definition, that%time, that%timerange, map_ttr, start=start, &
2286 full_steps=full_steps)
2287! merge with information from original volume
2288CALL vol7d_merge(that, v7dtmp)
2289
2290maxsize = maxval(map_ttr(:,:)%arraysize)
2291ALLOCATE(tmpvolr(maxsize), tmpvold(maxsize), lin_mask(maxsize), weights(maxsize))
2292do_otimerange: DO j = 1, SIZE(that%timerange)
2293 do_otime: DO i = 1, SIZE(that%time)
2294 ninp = map_ttr(i,j)%arraysize
2295 IF (ninp <= 0) cycle do_otime
2296! required for some computations
2297 CALL time_timerange_get_period(that%time(i), that%timerange(j), &
2298 that%time_definition, pstart, pend, reftime)
2299
2300 IF (ASSOCIATED(this%voldatir)) THEN
2301 DO i1 = 1, SIZE(this%ana)
2302 DO i3 = 1, SIZE(this%level)
2303 DO i6 = 1, SIZE(this%network)
2304 DO i5 = 1, SIZE(this%dativar%r)
2305! stat_proc difference treated separately here
2306 IF (stat_proc == 4) THEN
2307 IF (ninp >= 2) THEN
2308 IF (map_ttr(i,j)%array(1)%extra_info == 1 .AND. &
2309 map_ttr(i,j)%array(ninp)%extra_info == 2) THEN
2311 map_ttr(i,j)%array(1)%itr,i5,i6)) .AND. &
2312 c_e(this%voldatir(i1,map_ttr(i,j)%array(ninp)%it,i3, &
2313 map_ttr(i,j)%array(ninp)%itr,i5,i6))) THEN
2314 that%voldatir(i1,i,i3,j,i5,i6) = &
2315 this%voldatir(i1,map_ttr(i,j)%array(ninp)%it,i3, &
2316 map_ttr(i,j)%array(ninp)%itr,i5,i6) - &
2317 this%voldatir(i1,map_ttr(i,j)%array(1)%it,i3, &
2318 map_ttr(i,j)%array(1)%itr,i5,i6)
2319 ENDIF
2320 ENDIF
2321 ENDIF
2322 cycle
2323 ENDIF
2324! other stat_proc
2325 vartype = vol7d_vartype(this%dativar%r(i5))
2326 lin_mask = .false.
2327 ndtr = 0
2328 DO n = 1, ninp
2330 map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
2331 ndtr = ndtr + 1
2332 tmpvolr(ndtr) = this%voldatir(i1,map_ttr(i,j)%array(n)%it,i3, &
2333 map_ttr(i,j)%array(n)%itr,i5,i6)
2334 lin_mask(n) = .true.
2335 ENDIF
2336 ENDDO
2337 IF (ndtr == 0) cycle
2338 IF (lweighted) THEN
2339 CALL compute_stat_proc_agg_sw(map_ttr(i,j)%array(1:ninp)%time, &
2340 pstart, pend, lin_mask(1:ninp), act_max_step, weights)
2341 ELSE
2342 CALL compute_stat_proc_agg_sw(map_ttr(i,j)%array(1:ninp)%time, &
2343 pstart, pend, lin_mask(1:ninp), act_max_step)
2344 ENDIF
2345 IF (act_max_step > lmax_step) cycle
2346
2347 SELECT CASE(stat_proc)
2348 CASE (0) ! average
2349 IF (lweighted) THEN
2350 that%voldatir(i1,i,i3,j,i5,i6) = &
2351 sum(real(weights(1:ndtr))*tmpvolr(1:ndtr))
2352 ELSE
2353 that%voldatir(i1,i,i3,j,i5,i6) = &
2354 sum(tmpvolr(1:ndtr))/ndtr
2355 ENDIF
2356 CASE (2) ! maximum
2357 that%voldatir(i1,i,i3,j,i5,i6) = &
2358 maxval(tmpvolr(1:ndtr))
2359 CASE (3) ! minimum
2360 that%voldatir(i1,i,i3,j,i5,i6) = &
2361 minval(tmpvolr(1:ndtr))
2362 CASE (6) ! stddev
2363 that%voldatir(i1,i,i3,j,i5,i6) = &
2364 stat_stddev(tmpvolr(1:ndtr))
2365 CASE (201) ! mode
2366! mode only for angles at the moment, with predefined histogram
2367 IF (vartype == var_dir360) THEN
2368! remove undefined wind direction (=0), improve check?
2369! and reduce to interval [22.5,382.5[
2370 WHERE (tmpvolr(1:ndtr) == 0.0)
2371 tmpvolr(1:ndtr) = rmiss
2372 ELSE WHERE (tmpvolr(1:ndtr) < 22.5 .AND. tmpvolr(1:ndtr) > 0.0)
2373 tmpvolr(1:ndtr) = tmpvolr(1:ndtr) + 360.
2374 END WHERE
2375 that%voldatir(i1,i,i3,j,i5,i6) = &
2376 stat_mode_histogram(tmpvolr(1:ndtr), &
2377 8, 22.5, 382.5)
2378 ENDIF
2379 END SELECT
2380 ENDDO
2381 ENDDO
2382 ENDDO
2383 ENDDO
2384 ENDIF
2385
2386 IF (ASSOCIATED(this%voldatid)) THEN
2387 DO i1 = 1, SIZE(this%ana)
2388 DO i3 = 1, SIZE(this%level)
2389 DO i6 = 1, SIZE(this%network)
2390 DO i5 = 1, SIZE(this%dativar%d)
2391! stat_proc difference treated separately here
2392 IF (stat_proc == 4) THEN
2393 IF (n >= 2) THEN
2394 IF (map_ttr(i,j)%array(1)%extra_info == 1 .AND. &
2395 map_ttr(i,j)%array(ninp)%extra_info == 2) THEN
2397 map_ttr(i,j)%array(1)%itr,i5,i6)) .AND. &
2398 c_e(this%voldatid(i1,map_ttr(i,j)%array(ninp)%it,i3, &
2399 map_ttr(i,j)%array(ninp)%itr,i5,i6))) THEN
2400 that%voldatid(i1,i,i3,j,i5,i6) = &
2401 this%voldatid(i1,map_ttr(i,j)%array(ninp)%it,i3, &
2402 map_ttr(i,j)%array(ninp)%itr,i5,i6) - &
2403 this%voldatid(i1,map_ttr(i,j)%array(1)%it,i3, &
2404 map_ttr(i,j)%array(1)%itr,i5,i6)
2405 ENDIF
2406 ENDIF
2407 ENDIF
2408 cycle
2409 ENDIF
2410! other stat_proc
2411 vartype = vol7d_vartype(this%dativar%d(i5))
2412 lin_mask = .false.
2413 ndtr = 0
2414 DO n = 1, ninp
2416 map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
2417 ndtr = ndtr + 1
2418 tmpvold(ndtr) = this%voldatid(i1,map_ttr(i,j)%array(n)%it,i3, &
2419 map_ttr(i,j)%array(n)%itr,i5,i6)
2420 lin_mask(n) = .true.
2421 ENDIF
2422 ENDDO
2423 IF (ndtr == 0) cycle
2424 IF (lweighted) THEN
2425 CALL compute_stat_proc_agg_sw(map_ttr(i,j)%array(1:ninp)%time, &
2426 pstart, pend, lin_mask(1:ninp), act_max_step, weights)
2427 ELSE
2428 CALL compute_stat_proc_agg_sw(map_ttr(i,j)%array(1:ninp)%time, &
2429 pstart, pend, lin_mask(1:ninp), act_max_step)
2430 ENDIF
2431 IF (act_max_step > lmax_step) cycle
2432
2433 SELECT CASE(stat_proc)
2434 CASE (0) ! average
2435 IF (lweighted) THEN
2436 that%voldatid(i1,i,i3,j,i5,i6) = &
2437 sum(real(weights(1:ndtr))*tmpvold(1:ndtr))
2438 ELSE
2439 that%voldatid(i1,i,i3,j,i5,i6) = &
2440 sum(tmpvold(1:ndtr))/ndtr
2441 ENDIF
2442 CASE (2) ! maximum
2443 that%voldatid(i1,i,i3,j,i5,i6) = &
2444 maxval(tmpvold(1:ndtr))
2445 CASE (3) ! minimum
2446 that%voldatid(i1,i,i3,j,i5,i6) = &
2447 minval(tmpvold(1:ndtr))
2448 CASE (6) ! stddev
2449 that%voldatid(i1,i,i3,j,i5,i6) = &
2450 stat_stddev(tmpvold(1:ndtr))
2451 CASE (201) ! mode
2452! mode only for angles at the moment, with predefined histogram
2453 IF (vartype == var_dir360) THEN
2454! remove undefined wind direction (=0), improve check?
2455! and reduce to interval [22.5,382.5[
2456 WHERE (tmpvold(1:ndtr) == 0.0d0)
2457 tmpvold(1:ndtr) = dmiss
2458 ELSE WHERE (tmpvold(1:ndtr) < 22.5d0 .AND. tmpvold(1:ndtr) > 0.0d0)
2459 tmpvold(1:ndtr) = tmpvold(1:ndtr) + 360.0d0
2460 END WHERE
2461 that%voldatid(i1,i,i3,j,i5,i6) = &
2462 stat_mode_histogram(tmpvold(1:ndtr), &
2463 8, 22.5d0, 382.5d0)
2464 ENDIF
2465 END SELECT
2466 ENDDO
2467 ENDDO
2468 ENDDO
2469 ENDDO
2470 ENDIF
2471
2473 ENDDO do_otime
2474ENDDO do_otimerange
2475
2476DEALLOCATE(map_ttr)
2477DEALLOCATE(tmpvolr, tmpvold, lin_mask, weights)
2478
2479IF (PRESENT(other)) THEN ! create volume with the remaining data for further processing
2481 ltimerange=(this%timerange(:)%timerange /= tri))
2482ENDIF
2483
2484END SUBROUTINE vol7d_compute_stat_proc_agg
2485
2486
2502SUBROUTINE vol7d_decompute_stat_proc(this, that, step, other, stat_proc_input)
2503TYPE(vol7d),INTENT(inout) :: this
2504TYPE(vol7d),INTENT(out) :: that
2505TYPE(timedelta),INTENT(in) :: step
2506TYPE(vol7d),INTENT(inout),OPTIONAL :: other
2507INTEGER,INTENT(in),OPTIONAL :: stat_proc_input
2508
2509INTEGER :: i, tri, steps
2510
2511
2512IF (PRESENT(stat_proc_input)) THEN
2513 tri = stat_proc_input
2514ELSE
2515 tri = 0
2516ENDIF
2517! be safe
2518CALL vol7d_alloc_vol(this)
2519
2520! compute length of cumulation step in seconds
2522
2523! filter requested data
2525 ltimerange=(this%timerange(:)%timerange == tri .AND. &
2526 this%timerange(:)%p1 == 0 .AND. this%timerange(:)%p2 == steps))
2527
2528! convert timerange to instantaneous and go back half step in time
2529that%timerange(:)%timerange = 254
2530that%timerange(:)%p2 = 0
2531DO i = 1, SIZE(that%time(:))
2532 that%time(i) = that%time(i) - step/2
2533ENDDO
2534
2535IF (PRESENT(other)) THEN ! create volume with the remaining data for further processing
2537 ltimerange=(this%timerange(:)%timerange /= tri .OR. &
2538 this%timerange(:)%p1 /= 0 .OR. this%timerange(:)%p2 /= steps))
2539ENDIF
2540
2541END SUBROUTINE vol7d_decompute_stat_proc
2542
2543
2570SUBROUTINE vol7d_recompute_stat_proc_diff(this, that, stat_proc, step, full_steps, start, other)
2571TYPE(vol7d),INTENT(inout) :: this
2572TYPE(vol7d),INTENT(out) :: that
2573INTEGER,INTENT(in) :: stat_proc
2574TYPE(timedelta),INTENT(in) :: step
2575LOGICAL,INTENT(in),OPTIONAL :: full_steps
2576TYPE(datetime),INTENT(in),OPTIONAL :: start
2577TYPE(vol7d),INTENT(out),OPTIONAL :: other
2578
2579INTEGER :: i1, i3, i5, i6, i, j, k, l, nitr, steps
2580INTEGER,ALLOCATABLE :: map_tr(:,:,:,:,:), f(:), keep_tr(:,:,:)
2581LOGICAL,ALLOCATABLE :: mask_timerange(:)
2582LOGICAL,ALLOCATABLE :: mask_time(:)
2583TYPE(vol7d) :: v7dtmp
2584
2585
2586! be safe
2587CALL vol7d_alloc_vol(this)
2588! initialise the template of the output volume
2590
2591! compute length of cumulation step in seconds
2593
2594! compute the statistical processing relations, output time and
2595! timerange are defined here
2596CALL recompute_stat_proc_diff_common(this%time, this%timerange, stat_proc, step, &
2597 that%time, that%timerange, map_tr, f, keep_tr, &
2598 this%time_definition, full_steps, start)
2599nitr = SIZE(f)
2600
2601! complete the definition of the empty output template
2602CALL vol7d_alloc(that, nana=0, nlevel=0, nnetwork=0)
2603CALL vol7d_alloc_vol(that)
2604! shouldn't we exit here with an empty volume if stat_proc/=0,1 ?
2605ALLOCATE(mask_time(SIZE(this%time)), mask_timerange(SIZE(this%timerange)))
2606DO l = 1, SIZE(this%time)
2607 mask_time(l) = any(this%time(l) == that%time(:))
2608ENDDO
2609DO l = 1, SIZE(this%timerange)
2610 mask_timerange(l) = any(this%timerange(l) == that%timerange(:))
2611ENDDO
2612! create template for the output volume, keep all ana, level, network
2613! and variables; copy only the timeranges already satisfying the
2614! requested step, if any and only the times already existing in the
2615! output
2617 ltimerange=mask_timerange(:), ltime=mask_time(:))
2618! merge output created so far with template
2619CALL vol7d_merge(that, v7dtmp, lanasimple=.true., llevelsimple=.true.)
2620
2621! compute statistical processing
2622IF (ASSOCIATED(this%voldatir)) THEN
2623 DO l = 1, SIZE(this%time)
2624 DO k = 1, nitr
2625 DO j = 1, SIZE(this%time)
2626 DO i = 1, nitr
2628 DO i6 = 1, SIZE(this%network)
2629 DO i5 = 1, SIZE(this%dativar%r)
2630 DO i3 = 1, SIZE(this%level)
2631 DO i1 = 1, SIZE(this%ana)
2634
2635 IF (stat_proc == 0) THEN ! average
2636 that%voldatir( &
2637 i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6) = &
2638 (this%voldatir(i1,l,i3,f(k),i5,i6)*this%timerange(f(k))%p2 - &
2639 this%voldatir(i1,j,i3,f(i),i5,i6)*this%timerange(f(i))%p2)/ &
2640 steps ! optimize avoiding conversions
2641 ELSE IF (stat_proc == 1 .OR. stat_proc == 4) THEN ! acc, diff
2642 that%voldatir( &
2643 i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6) = &
2644 this%voldatir(i1,l,i3,f(k),i5,i6) - &
2645 this%voldatir(i1,j,i3,f(i),i5,i6)
2646 ENDIF
2647
2648 ENDIF
2649 ENDDO
2650 ENDDO
2651 ENDDO
2652 ENDDO
2653 ENDIF
2654 ENDDO
2655 ENDDO
2656 ENDDO
2657 ENDDO
2658ENDIF
2659
2660IF (ASSOCIATED(this%voldatid)) THEN
2661 DO l = 1, SIZE(this%time)
2662 DO k = 1, nitr
2663 DO j = 1, SIZE(this%time)
2664 DO i = 1, nitr
2666 DO i6 = 1, SIZE(this%network)
2667 DO i5 = 1, SIZE(this%dativar%d)
2668 DO i3 = 1, SIZE(this%level)
2669 DO i1 = 1, SIZE(this%ana)
2672! IF (.NOT.c_e(that%voldatid( &
2673! i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6))) THEN
2674
2675 IF (stat_proc == 0) THEN ! average
2676 that%voldatid( &
2677 i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6) = &
2678 (this%voldatid(i1,l,i3,f(k),i5,i6)*this%timerange(f(k))%p2 - &
2679 this%voldatid(i1,j,i3,f(i),i5,i6)*this%timerange(f(i))%p2)/ &
2680 steps ! optimize avoiding conversions
2681 ELSE IF (stat_proc == 1 .OR. stat_proc == 4) THEN ! acc, diff
2682 that%voldatid( &
2683 i1,map_tr(i,j,k,l,1),i3,map_tr(i,j,k,l,2),i5,i6) = &
2684 this%voldatid(i1,l,i3,f(k),i5,i6) - &
2685 this%voldatid(i1,j,i3,f(i),i5,i6)
2686 ENDIF
2687
2688! ENDIF
2689 ENDIF
2690 ENDDO
2691 ENDDO
2692 ENDDO
2693 ENDDO
2694 ENDIF
2695 ENDDO
2696 ENDDO
2697 ENDDO
2698 ENDDO
2699ENDIF
2700
2701! this should be avoided by sorting descriptors upstream
2702! descriptors now are sorted upstream with a dirty and expensive trick
2703! but the order may be scrambled in the call to vol7d_merge above
2704CALL vol7d_smart_sort(that, lsort_time=.true., lsort_timerange=.true.)
2705
2706CALL makeother(.true.)
2707
2708CONTAINS
2709
2710SUBROUTINE makeother(filter)
2711LOGICAL,INTENT(in) :: filter
2712IF (PRESENT(other)) THEN
2713 IF (filter) THEN ! create volume with the remaining data for further processing
2715 ltimerange=(this%timerange(:)%timerange /= stat_proc))
2716 ELSE
2718 ENDIF
2719ENDIF
2720END SUBROUTINE makeother
2721
2722END SUBROUTINE vol7d_recompute_stat_proc_diff
2723
2724
2752SUBROUTINE vol7d_compute_stat_proc_metamorph(this, that, stat_proc_input, stat_proc)
2753TYPE(vol7d),INTENT(inout) :: this
2754TYPE(vol7d),INTENT(out) :: that
2755INTEGER,INTENT(in) :: stat_proc_input
2756INTEGER,INTENT(in) :: stat_proc
2757
2758INTEGER :: j
2759LOGICAL,ALLOCATABLE :: tr_mask(:)
2760REAL,ALLOCATABLE :: int_ratio(:)
2761DOUBLE PRECISION,ALLOCATABLE :: int_ratiod(:)
2762
2763IF (.NOT.((stat_proc_input == 0 .AND. stat_proc == 1) .OR. &
2764 (stat_proc_input == 1 .AND. stat_proc == 0))) THEN
2765
2766 CALL l4f_log(l4f_warn, &
2767 'compute_stat_proc_metamorph, can only be applied to average->accumulated timerange and viceversa')
2768! return an empty volume, without signaling error
2770 CALL vol7d_alloc_vol(that)
2771 RETURN
2772ENDIF
2773
2774! be safe
2775CALL vol7d_alloc_vol(this)
2776
2777! useful timeranges
2778tr_mask = this%timerange(:)%timerange == stat_proc_input .AND. this%timerange(:)%p2 /= imiss &
2779 .AND. this%timerange(:)%p2 /= 0
2780
2781! initial check (necessary?)
2782IF (count(tr_mask) == 0) THEN
2783 CALL l4f_log(l4f_warn, &
2784 'vol7d_compute, no timeranges suitable for statistical processing by metamorphosis')
2786! CALL makeother()
2787 RETURN
2788ENDIF
2789
2790! copy required data and reset timerange
2791CALL vol7d_copy(this, that, ltimerange=tr_mask)
2792that%timerange(:)%timerange = stat_proc
2793! why next automatic f2003 allocation does not always work?
2794ALLOCATE(int_ratio(SIZE(that%timerange)), int_ratiod(SIZE(that%timerange)))
2795
2796IF (stat_proc == 0) THEN ! average -> integral
2797 int_ratio = 1./real(that%timerange(:)%p2)
2798 int_ratiod = 1./dble(that%timerange(:)%p2)
2799ELSE ! cumulation
2800 int_ratio = real(that%timerange(:)%p2)
2801 int_ratiod = dble(that%timerange(:)%p2)
2802ENDIF
2803
2804IF (ASSOCIATED(that%voldatir)) THEN
2805 DO j = 1, SIZE(that%timerange)
2807 that%voldatir(:,:,:,j,:,:) = that%voldatir(:,:,:,j,:,:)*int_ratio(j)
2808 ELSEWHERE
2809 that%voldatir(:,:,:,j,:,:) = rmiss
2810 END WHERE
2811 ENDDO
2812ENDIF
2813
2814IF (ASSOCIATED(that%voldatid)) THEN
2815 DO j = 1, SIZE(that%timerange)
2817 that%voldatid(:,:,:,j,:,:) = that%voldatid(:,:,:,j,:,:)*int_ratiod(j)
2818 ELSEWHERE
2819 that%voldatid(:,:,:,j,:,:) = rmiss
2820 END WHERE
2821 ENDDO
2822ENDIF
2823
2824
2825END SUBROUTINE vol7d_compute_stat_proc_metamorph
2826
2827
2828SUBROUTINE vol7d_recompute_stat_proc_agg_multiv(this, that, &
2829 step, start, frac_valid, multiv_proc)
2830TYPE(vol7d),INTENT(inout) :: this
2831TYPE(vol7d),INTENT(out) :: that
2832!INTEGER,INTENT(in) :: stat_proc !< type of statistical processing to be recomputed (from grib2 table), only data having timerange of this type will be recomputed and will appear in the output volume
2833TYPE(timedelta),INTENT(in) :: step
2834TYPE(datetime),INTENT(in),OPTIONAL :: start
2835REAL,INTENT(in),OPTIONAL :: frac_valid
2836!TYPE(vol7d),INTENT(inout),OPTIONAL :: other !< optional volume that, on exit, is going to contain the data that did not contribute to the statistical processing
2837!INTEGER,INTENT(in),OPTIONAL :: stat_proc_input !< to be used with care, type of statistical processing of data that has to be processed (from grib2 table), only data having timerange of this type will be recomputed, the actual statistical processing performed and which will appear in the output volume, is however determined by \a stat_proc argument
2838INTEGER,INTENT(in) :: multiv_proc
2839
2840INTEGER :: tri
2841INTEGER :: i, j, n, n1, ndtr, i1, i3, i5, i6
2842INTEGER :: linshape(1)
2843REAL :: lfrac_valid, frac_c, frac_m
2844LOGICAL,ALLOCATABLE :: ttr_mask(:,:)
2845TYPE(arrayof_ttr_mapper),POINTER :: map_ttr(:,:)
2846INTEGER,POINTER :: dtratio(:)
2847INTEGER :: stat_proc_input, stat_proc
2848
2849SELECT CASE(multiv_proc)
2850CASE (1) ! direction of maximum
2851 stat_proc_input = 205
2852 stat_proc=205
2853END SELECT
2854
2855tri = stat_proc_input
2856IF (PRESENT(frac_valid)) THEN
2857 lfrac_valid = frac_valid
2858ELSE
2859 lfrac_valid = 1.0
2860ENDIF
2861
2862! be safe
2863CALL vol7d_alloc_vol(this)
2864! initial check
2865
2866! cleanup the original volume
2867CALL vol7d_smart_sort(this, lsort_time=.true.) ! time-ordered volume needed
2869
2871CALL vol7d_alloc(that, nana=SIZE(this%ana), nlevel=SIZE(this%level), &
2872 nnetwork=SIZE(this%network))
2873IF (ASSOCIATED(this%dativar%r)) THEN
2874 CALL vol7d_alloc(that, ndativarr=SIZE(this%dativar%r))
2875 that%dativar%r = this%dativar%r
2876ENDIF
2877IF (ASSOCIATED(this%dativar%d)) THEN
2878 CALL vol7d_alloc(that, ndativard=SIZE(this%dativar%d))
2879 that%dativar%d = this%dativar%d
2880ENDIF
2881that%ana = this%ana
2882that%level = this%level
2883that%network = this%network
2884
2885! compute the output time and timerange and all the required mappings
2886CALL recompute_stat_proc_agg_common(this%time, this%timerange, stat_proc, tri, &
2887 step, this%time_definition, that%time, that%timerange, map_ttr, &
2888 dtratio=dtratio, start=start)
2889CALL vol7d_alloc_vol(that)
2890
2891ALLOCATE(ttr_mask(SIZE(this%time), SIZE(this%timerange)))
2892linshape = (/SIZE(ttr_mask)/)
2893! finally perform computations
2894IF (ASSOCIATED(this%voldatir)) THEN
2895 DO j = 1, SIZE(that%timerange)
2896 DO i = 1, SIZE(that%time)
2897
2898 DO i1 = 1, SIZE(this%ana)
2899 DO i3 = 1, SIZE(this%level)
2900 DO i6 = 1, SIZE(this%network)
2901 DO i5 = 1, SIZE(this%dativar%r)
2902
2903 frac_m = 0.
2904 DO n1 = SIZE(dtratio), 1, -1 ! precedence to longer periods
2905 IF (dtratio(n1) <= 0) cycle ! safety check
2906 ttr_mask = .false.
2907 DO n = 1, map_ttr(i,j)%arraysize
2908 IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN
2910 map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
2911 ttr_mask(map_ttr(i,j)%array(n)%it, &
2912 map_ttr(i,j)%array(n)%itr) = .true.
2913 ENDIF
2914 ENDIF
2915 ENDDO
2916
2917 ndtr = count(ttr_mask)
2918 frac_c = real(ndtr)/real(dtratio(n1))
2919
2920 IF (ndtr > 0 .AND. frac_c >= max(lfrac_valid, frac_m)) THEN
2921 frac_m = frac_c
2922 SELECT CASE(multiv_proc)
2923 CASE (1) ! average, vectorial mean
2924 that%voldatir(i1,i,i3,j,i5,i6) = &
2925 sum(this%voldatir(i1,:,i3,:,i5,i6), &
2926 mask=ttr_mask)/ndtr
2927 END SELECT
2928 ENDIF
2929
2930 ENDDO ! dtratio
2931 ENDDO ! var
2932 ENDDO ! network
2933 ENDDO ! level
2934 ENDDO ! ana
2936 ENDDO ! otime
2937 ENDDO ! otimerange
2938ENDIF
2939
2940IF (ASSOCIATED(this%voldatid)) THEN
2941 DO j = 1, SIZE(that%timerange)
2942 DO i = 1, SIZE(that%time)
2943
2944 DO i1 = 1, SIZE(this%ana)
2945 DO i3 = 1, SIZE(this%level)
2946 DO i6 = 1, SIZE(this%network)
2947 DO i5 = 1, SIZE(this%dativar%d)
2948
2949 frac_m = 0.
2950 DO n1 = SIZE(dtratio), 1, -1 ! precedence to longer periods
2951 IF (dtratio(n1) <= 0) cycle ! safety check
2952 ttr_mask = .false.
2953 DO n = 1, map_ttr(i,j)%arraysize
2954 IF (map_ttr(i,j)%array(n)%extra_info == dtratio(n1)) THEN
2956 map_ttr(i,j)%array(n)%itr,i5,i6))) THEN
2957 ttr_mask(map_ttr(i,j)%array(n)%it, &
2958 map_ttr(i,j)%array(n)%itr) = .true.
2959 ENDIF
2960 ENDIF
2961 ENDDO
2962
2963 ndtr = count(ttr_mask)
2964 frac_c = real(ndtr)/real(dtratio(n1))
2965
2966 IF (ndtr > 0 .AND. frac_c >= max(lfrac_valid, frac_m)) THEN
2967 frac_m = frac_c
2968 SELECT CASE(stat_proc)
2969 CASE (0) ! average
2970 that%voldatid(i1,i,i3,j,i5,i6) = &
2971 sum(this%voldatid(i1,:,i3,:,i5,i6), &
2972 mask=ttr_mask)/ndtr
2973 CASE (1, 4) ! accumulation, difference
2974 that%voldatid(i1,i,i3,j,i5,i6) = &
2975 sum(this%voldatid(i1,:,i3,:,i5,i6), &
2976 mask=ttr_mask)
2977 CASE (2) ! maximum
2978 that%voldatid(i1,i,i3,j,i5,i6) = &
2979 maxval(this%voldatid(i1,:,i3,:,i5,i6), &
2980 mask=ttr_mask)
2981 CASE (3) ! minimum
2982 that%voldatid(i1,i,i3,j,i5,i6) = &
2983 minval(this%voldatid(i1,:,i3,:,i5,i6), &
2984 mask=ttr_mask)
2985 CASE (6) ! stddev
2986 that%voldatid(i1,i,i3,j,i5,i6) = &
2987 stat_stddev( &
2988 reshape(this%voldatid(i1,:,i3,:,i5,i6), shape=linshape), &
2989 mask=reshape(ttr_mask, shape=linshape))
2990 END SELECT
2991 ENDIF
2992
2993 ENDDO ! dtratio
2994 ENDDO ! var
2995 ENDDO ! network
2996 ENDDO ! level
2997 ENDDO ! ana
2999 ENDDO ! otime
3000 ENDDO ! otimerange
3001ENDIF
3002
3003DEALLOCATE(ttr_mask)
3004
3005END SUBROUTINE vol7d_recompute_stat_proc_agg_multiv
3006
3023SUBROUTINE vol7d_fill_time(this, that, step, start, stopp, cyclicdt)
3024TYPE(vol7d),INTENT(inout) :: this
3025TYPE(vol7d),INTENT(inout) :: that
3026TYPE(timedelta),INTENT(in) :: step
3027TYPE(datetime),INTENT(in),OPTIONAL :: start
3028TYPE(datetime),INTENT(in),OPTIONAL :: stopp
3029TYPE(cyclicdatetime),INTENT(in),OPTIONAL :: cyclicdt
3030
3031TYPE(cyclicdatetime) :: lcyclicdt
3032TYPE(datetime) :: counter, lstart, lstop
3033INTEGER :: i, naddtime
3034
3035CALL safe_start_stop(this, lstart, lstop, start, stopp)
3037
3038lcyclicdt=cyclicdatetime_miss
3039if (present(cyclicdt)) then
3041end if
3042
3045
3046! Count the number of time levels required for completing the series
3047! valid also in the case (SIZE(this%time) == 0)
3048naddtime = 0
3049counter = lstart
3050i = 1
3051naddcount: DO WHILE(counter <= lstop)
3052 DO WHILE(i <= SIZE(this%time)) ! this%time(i) chases counter
3053 IF (counter < this%time(i)) THEN ! this%time(i) overtook counter
3054 i = max(i-1,1) ! go back if possible
3055 EXIT
3056 ELSE IF (counter == this%time(i) .OR. .NOT. counter == lcyclicdt) THEN ! found matching time
3057 counter = counter + step
3058 cycle naddcount
3059 ENDIF
3060 i = i + 1
3061 ENDDO
3062 naddtime = naddtime + 1
3063 counter = counter + step
3064ENDDO naddcount
3065
3066! old universal algorithm, not optimized, check that the new one is equivalent
3067!naddtime = 0
3068!counter = lstart
3069!DO WHILE(counter <= lstop)
3070! IF (.NOT.ANY(counter == this%time(:))) THEN
3071! naddtime = naddtime + 1
3072! ENDIF
3073! counter = counter + step
3074!ENDDO
3075
3076IF (naddtime > 0) THEN
3077
3079 CALL vol7d_alloc(that, ntime=naddtime)
3080 CALL vol7d_alloc_vol(that)
3081
3082 ! Repeat the count loop setting the time levels to be added
3083 naddtime = 0
3084 counter = lstart
3085 i = 1
3086 naddadd: DO WHILE(counter <= lstop)
3087 DO WHILE(i <= SIZE(this%time)) ! this%time(i) chases counter
3088 IF (counter < this%time(i)) THEN ! this%time(i) overtook counter
3089 i = max(i-1,1) ! go back if possible
3090 EXIT
3091 ELSE IF (counter == this%time(i) .OR. .NOT. counter == lcyclicdt) THEN ! found matching time
3092 counter = counter + step
3093 cycle naddadd
3094 ENDIF
3095 i = i + 1
3096 ENDDO
3097 naddtime = naddtime + 1
3098 that%time(naddtime) = counter ! only difference
3099 counter = counter + step
3100 ENDDO naddadd
3101
3103
3104ELSE
3105!! ? why sort all dimension ?
3106!! CALL vol7d_copy(this, that, lsort_time=.TRUE.)
3108ENDIF
3109
3110
3111END SUBROUTINE vol7d_fill_time
3112
3113
3125SUBROUTINE vol7d_filter_time(this, that, step, start, stopp, cyclicdt)
3126TYPE(vol7d),INTENT(inout) :: this
3127TYPE(vol7d),INTENT(inout) :: that
3128TYPE(timedelta),INTENT(in),optional :: step
3129TYPE(datetime),INTENT(in),OPTIONAL :: start
3130TYPE(datetime),INTENT(in),OPTIONAL :: stopp
3131TYPE(cyclicdatetime),INTENT(in),OPTIONAL :: cyclicdt
3132
3133TYPE(datetime) :: lstart, lstop
3134LOGICAL, ALLOCATABLE :: time_mask(:)
3135
3136CALL safe_start_stop(this, lstart, lstop, start, stopp)
3138
3141
3142ALLOCATE(time_mask(SIZE(this%time)))
3143
3144time_mask = this%time >= lstart .AND. this%time <= lstop
3145
3146IF (PRESENT(cyclicdt)) THEN
3148 time_mask = time_mask .AND. this%time == cyclicdt
3149 ENDIF
3150ENDIF
3151
3152IF (PRESENT(step)) THEN
3154 time_mask = time_mask .AND. mod(this%time - lstart, step) == timedelta_0
3155 ENDIF
3156ENDIF
3157
3158CALL vol7d_copy(this,that, ltime=time_mask)
3159
3160DEALLOCATE(time_mask)
3161
3162END SUBROUTINE vol7d_filter_time
3163
3164
3168SUBROUTINE vol7d_fill_data(this, step, start, stopp, tolerance)
3169TYPE(vol7d),INTENT(inout) :: this
3170TYPE(timedelta),INTENT(in) :: step
3171TYPE(datetime),INTENT(in),OPTIONAL :: start
3172TYPE(datetime),INTENT(in),OPTIONAL :: stopp
3173TYPE(timedelta),INTENT(in),optional :: tolerance
3174
3175TYPE(datetime) :: lstart, lstop
3176integer :: indana , indtime ,indlevel ,indtimerange ,inddativarr, indnetwork, iindtime
3177type(timedelta) :: deltato,deltat, ltolerance
3178
3179CALL safe_start_stop(this, lstart, lstop, start, stopp)
3181
3184
3185
3186ltolerance=step/2
3187
3188if (present(tolerance))then
3190end if
3191
3192
3193do indtime=1,size(this%time)
3194
3195 IF (this%time(indtime) < lstart .OR. this%time(indtime) > lstop .OR. &
3196 mod(this%time(indtime) - lstart, step) /= timedelta_0) cycle
3197 do indtimerange=1,size(this%timerange)
3198 if (this%timerange(indtimerange)%timerange /= 254) cycle
3199 do indnetwork=1,size(this%network)
3200 do inddativarr=1,size(this%dativar%r)
3201 do indlevel=1,size(this%level)
3202 do indana=1,size(this%ana)
3203
3204 !find the nearest data in time if data is missing
3205 if (.not. c_e(this%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork)))then
3206 deltato=timedelta_miss
3207
3208 !do iindtime=max(indtime-20,1),min(indtime+20,size(this%time)) !check on a chunk: 20 should be enought
3209
3210 do iindtime=indtime+1,size(this%time) !check forward
3211
3212 if (c_e(this%voldatir (indana, iindtime, indlevel, indtimerange, inddativarr, indnetwork )))then
3213 deltat=this%time(iindtime)-this%time(indtime)
3214
3215 if (deltat >= ltolerance) exit
3216
3217 if (deltat < deltato) then
3218 this%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork) = &
3219 this%voldatir(indana, iindtime, indlevel, indtimerange, inddativarr, indnetwork)
3220 deltato=deltat
3221 end if
3222 end if
3223 end do
3224
3225 do iindtime=indtime-1,1,-1 !check backward
3226
3227 if (c_e(this%voldatir (indana, iindtime, indlevel, indtimerange, inddativarr, indnetwork )))then
3228 if (iindtime < indtime) then
3229 deltat=this%time(indtime)-this%time(iindtime)
3230 else if (iindtime > indtime) then
3231 deltat=this%time(iindtime)-this%time(indtime)
3232 else
3233 cycle
3234 end if
3235
3236 if (deltat >= ltolerance) exit
3237
3238 if (deltat < deltato) then
3239 this%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork) = &
3240 this%voldatir(indana, iindtime, indlevel, indtimerange, inddativarr, indnetwork)
3241 deltato=deltat
3242 end if
3243 end if
3244 end do
3245
3246 end if
3247 end do
3248 end do
3249 end do
3250 end do
3251 end do
3252end do
3253
3254END SUBROUTINE vol7d_fill_data
3255
3256
3257! private utility routine for checking interval and start-stop times
3258! in input missing start-stop values are treated as not present
3259! in output missing start-stop values mean "do nothing"
3260SUBROUTINE safe_start_stop(this, lstart, lstop, start, stopp)
3261TYPE(vol7d),INTENT(inout) :: this
3262TYPE(datetime),INTENT(out) :: lstart
3263TYPE(datetime),INTENT(out) :: lstop
3264TYPE(datetime),INTENT(in),OPTIONAL :: start
3265TYPE(datetime),INTENT(in),OPTIONAL :: stopp
3266
3267lstart = datetime_miss
3268lstop = datetime_miss
3269! initial safety operation
3270CALL vol7d_alloc_vol(this)
3271IF (SIZE(this%time) == 0) RETURN ! avoid segmentation fault in case of empty volume
3272CALL vol7d_smart_sort(this, lsort_time=.true.)
3273
3274IF (PRESENT(start)) THEN
3276 lstart = start
3277 ELSE
3278 lstart = this%time(1)
3279 ENDIF
3280ELSE
3281 lstart = this%time(1)
3282ENDIF
3283IF (PRESENT(stopp)) THEN
3285 lstop = stopp
3286 ELSE
3287 lstop = this%time(SIZE(this%time))
3288 ENDIF
3289ELSE
3290 lstop = this%time(SIZE(this%time))
3291ENDIF
3292
3293END SUBROUTINE safe_start_stop
3294
3295
3302SUBROUTINE vol7d_normalize_vcoord(this,that,ana,time,timerange,network)
3303TYPE(vol7d),INTENT(INOUT) :: this
3304TYPE(vol7d),INTENT(OUT) :: that
3305integer,intent(in) :: time,ana,timerange,network
3306
3307character(len=1) :: type
3308integer :: ind
3309TYPE(vol7d_var) :: var
3310LOGICAL,allocatable :: ltime(:),ltimerange(:),lana(:),lnetwork(:)
3311logical,allocatable :: maschera(:)
3312
3313
3314allocate(ltime(size(this%time)))
3315allocate(ltimerange(size(this%timerange)))
3316allocate(lana(size(this%ana)))
3317allocate(lnetwork(size(this%network)))
3318
3319ltime=.false.
3320ltimerange=.false.
3321lana=.false.
3322lnetwork=.false.
3323
3324ltime(time)=.true.
3325ltimerange(timerange)=.true.
3326lana(ana)=.true.
3327lnetwork(network)=.true.
3328
3329call vol7d_copy(this, that,unique=.true.,&
3330 ltime=ltime,ltimerange=ltimerange,lana=lana,lnetwork=lnetwork )
3331
3333type=cmiss
3334!type="i"
3335ind = index(that%dativar, var, type=type)
3336
3337allocate(maschera(size(that%level)))
3338
3339maschera = (&
3340 (that%level%level1 == 105.and.that%level%level2 == 105) .or. &
3341 (that%level%level1 == 103 .and. that%level%level2 == imiss ) .or. &
3342 (that%level%level1 == 102 .and. that%level%level2 == imiss )) &
3343 .and. c_e(that%voldatic(1,1,:,1,ind,1))
3344
3345
3346select case (type)
3347
3348case("d")
3349
3350 where (maschera)
3351 that%level%level1 = 100
3352 that%level%l1 = int(realdat(that%voldatid(1,1,:,1,ind,1),that%dativar%d(ind)))
3353 that%level%l1 = int(that%voldatid(1,1,:,1,ind,1))
3354 that%level%level2 = imiss
3355 that%level%l2 = imiss
3356 end where
3357
3358case("r")
3359
3360 where (maschera)
3361 that%level%level1 = 100
3362 that%level%l1 = int(realdat(that%voldatir(1,1,:,1,ind,1),that%dativar%r(ind)))
3363 that%level%level2 = imiss
3364 that%level%l2 = imiss
3365 end where
3366
3367case("i")
3368
3369 where (maschera)
3370 that%level%level1 = 100
3371 that%level%l1 = int(realdat(that%voldatii(1,1,:,1,ind,1),that%dativar%i(ind)))
3372 that%level%level2 = imiss
3373 that%level%l2 = imiss
3374 end where
3375
3376case("b")
3377
3378 where (maschera)
3379 that%level%level1 = 100
3380 that%level%l1 = int(realdat(that%voldatib(1,1,:,1,ind,1),that%dativar%b(ind)))
3381 that%level%level2 = imiss
3382 that%level%l2 = imiss
3383 end where
3384
3385case("c")
3386
3387 where (maschera)
3388 that%level%level1 = 100
3389 that%level%l1 = int(realdat(that%voldatic(1,1,:,1,ind,1),that%dativar%c(ind)))
3390 that%level%level2 = imiss
3391 that%level%l2 = imiss
3392 end where
3393
3394end select
3395
3396deallocate(ltime)
3397deallocate(ltimerange)
3398deallocate(lana)
3399deallocate(lnetwork)
3400
3401END SUBROUTINE vol7d_normalize_vcoord
3402
3403
3404!!$!> Metodo per calcolare variabili derivate.
3405!!$!! TO DO !!
3406!!$SUBROUTINE vol7d_compute_var(this,that,var)
3407!!$TYPE(vol7d),INTENT(INOUT) :: this !< oggetto da normalizzare
3408!!$TYPE(vol7d),INTENT(OUT) :: that !< oggetto normalizzato
3409!!$
3410!!$character(len=1) :: type
3411!!$TYPE(vol7d_var),intent(in) :: var
3412
3413
3414!!$call init(var, btable="B10004") ! Pressure
3415!!$type=cmiss
3416!!$call vol7d_varvect_index(that%dativar,var , type=type,index_v=ind)
3417!!$
3418!!$select case (type)
3419!!$
3420!!$case("d")
3421!!$
3422!!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
3423!!$ that%level%level1 = 100
3424!!$ that%level%l1 = realdat(that%voldatid(1,1,:,1,ind,1),that%dativar%d(ind))
3425!!$ that%level%level2 = imiss
3426!!$ that%level%l2 = imiss
3427!!$ end where
3428!!$
3429!!$case("r")
3430!!$
3431!!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
3432!!$ that%level%level1 = 100
3433!!$ that%level%l1 = realdat(that%voldatir(1,1,:,1,ind,1),that%dativar%r(ind))
3434!!$ that%level%level2 = imiss
3435!!$ that%level%l2 = imiss
3436!!$ end where
3437!!$
3438!!$case("i")
3439!!$
3440!!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
3441!!$ that%level%level1 = 100
3442!!$ that%level%l1 = realdat(that%voldatii(1,1,:,1,ind,1),that%dativar%i(ind))
3443!!$ that%level%level2 = imiss
3444!!$ that%level%l2 = imiss
3445!!$ end where
3446!!$
3447!!$case("b")
3448!!$
3449!!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
3450!!$ that%level%level1 = 100
3451!!$ that%level%l1 = realdat(that%voldatib(1,1,:,1,ind,1),that%dativar%b(ind))
3452!!$ that%level%level2 = imiss
3453!!$ that%level%l2 = imiss
3454!!$ end where
3455!!$
3456!!$case("c")
3457!!$
3458!!$ where (that%level%level1 == 105.and.that%level%level2 == 105)
3459!!$ that%level%level1 = 100
3460!!$ that%level%l1 = realdat(that%voldatic(1,1,:,1,ind,1),that%dativar%c(ind))
3461!!$ that%level%level2 = imiss
3462!!$ that%level%l2 = imiss
3463!!$ end where
3464!!$
3465!!$end select
3466
3467!!$
3468!!$END SUBROUTINE vol7d_compute_var
3469!!$
3470
3471
3472
Restituiscono il valore dell'oggetto nella forma desiderata. Definition: datetime_class.F90:328 Costruttori per le classi datetime e timedelta. Definition: datetime_class.F90:317 Functions that return a trimmed CHARACTER representation of the input variable. Definition: datetime_class.F90:355 Restituiscono il valore dell'oggetto in forma di stringa stampabile. Definition: datetime_class.F90:333 Compute the mode of the random variable provided taking into account missing data. Definition: simple_stat.f90:160 Compute the standard deviation of the random variable provided, taking into account missing data. Definition: simple_stat.f90:69 Module for basic statistical computations taking into account missing data. Definition: simple_stat.f90:25 This module contains functions that are only for internal use of the library. Definition: stat_proc_engine.F90:217 Extension of vol7d_class with methods for performing simple statistical operations on entire volumes ... Definition: vol7d_class_compute.F90:220 Classe per la gestione di un volume completo di dati osservati. Definition: vol7d_class.F90:279 |