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