Line data Source code
1 : !> \file mo_smi.f90
2 : !> \brief \copybrief mo_smi
3 : !> \details \copydetails mo_smi
4 :
5 : !> \brief SOIL MOISTURE INDEX
6 : !> \details Estimate SMI based on monthly values of SM_mHM
7 : !> \author Luis Samaniego
8 : !> \date Feb 2011
9 : !> \author Stephan Thober
10 : !> \date 16.12.2014
11 : !! - added evalSMI capability
12 : !> \author Stephan Thober
13 : !> \date 07.11.2017
14 : !! - added inversion of CDFs
15 : !> \author Stephan Thober
16 : !> \date 02.08.2019
17 : !! - added non-full year
18 : !> \copyright Copyright 2005-\today, the CHS Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
19 : !! SMI is released under the LGPLv3+ license \license_note
20 : module mo_smi
21 : !
22 : use mo_smi_global_variables, only: period
23 : use mo_kind, only: dp
24 : !
25 : implicit none
26 : ! public routines
27 : public :: optimize_width
28 : public :: calSMI
29 : public :: invSMI
30 : !
31 :
32 : private
33 :
34 :
35 : contains
36 :
37 : !> \brief subroutine for estimating SMI for first array
38 1 : subroutine optimize_width( opt_h, silverman_h, SM, nCalendarStepsYear, per_kde)
39 : !$ use omp_lib
40 : use mo_kind, only : i4, sp
41 : use mo_kernel, only : kernel_density_h
42 : implicit none
43 :
44 : ! input variables
45 : logical, intent(in) :: silverman_h ! optimize kernel width
46 : real(sp), dimension(:,:), intent(in) :: SM
47 : integer(i4), intent(in) :: nCalendarStepsYear
48 : type(period), intent(in) :: per_kde
49 :
50 : ! output variables
51 : real(dp), dimension(:,:), intent(inout) :: opt_h
52 :
53 : ! local variables
54 : integer(i4) :: ii ! cell index
55 : integer(i4) :: mm ! month index
56 1 : real(dp), dimension(:), allocatable :: X
57 1 : logical, allocatable :: t_mask_kde(:)
58 1 : integer(i4), allocatable :: time_kde(:)
59 :
60 :
61 : call get_time_indizes(time_kde, per_kde, nCalendarStepsYear)
62 :
63 : !$OMP parallel default(shared) &
64 : !$OMP private(ii,mm,X, t_mask_kde)
65 : !$OMP do SCHEDULE(STATIC)
66 13 : do mm = 1, nCalendarStepsYear
67 :
68 4344 : t_mask_kde = (time_kde .eq. mm)
69 :
70 277 : do ii = 1, size( SM, 1 )
71 : #ifdef SMIDEBUG
72 : if ((mod(ii, 100) .eq. 0_i4) .and. (mm .eq. nCalendarStepsYear)) print *, ii, mm
73 : #endif
74 : !
75 95832 : allocate(X(count(t_mask_kde)))
76 8184 : X(:) = pack(SM( ii, :), t_mask_kde)
77 : !
78 : ! determine kernel width if these have not been read from file
79 264 : opt_h(ii,mm) = kernel_density_h( X(:), silverman = silverman_h )
80 276 : deallocate(X)
81 : end do
82 : end do
83 : !$OMP end do
84 : !$OMP end parallel
85 1 : end subroutine optimize_width
86 :
87 : !!======================================================
88 : !! ESTIMATE SMI
89 : !!======================================================
90 : !> \brief subroutine for calculating SMI for second array with pdf of first one
91 4 : subroutine calSMI( hh, sm_kde, sm_eval, nCalendarStepsYear, SMI, per_kde, per_eval)
92 :
93 : use mo_kind, only: i4, sp, dp
94 : use mo_kernel, only: kernel_cumdensity
95 : use mo_constants, only: nodata_dp, nodata_sp
96 :
97 : implicit none
98 :
99 : ! input variables
100 : real(dp), dimension(:,:), intent(in) :: hh
101 : real(sp), dimension(:,:), intent(in) :: sm_kde
102 : real(sp), dimension(:,:), intent(in) :: sm_eval
103 : integer(i4), intent(in) :: nCalendarStepsYear
104 : type(period), intent(in) :: per_kde
105 : type(period), intent(in) :: per_eval
106 :
107 : ! output variables
108 : real(sp), dimension(:,:), intent(out) :: SMI
109 :
110 : ! local variables
111 : integer(i4) :: mm ! loop index
112 : integer(i4) :: ii ! cell index
113 : real(dp), dimension(:), allocatable :: cdf
114 : real(dp), dimension(:), allocatable :: X_kde
115 : real(dp), dimension(:), allocatable :: X_eval
116 : integer(i4) :: n_time
117 4 : logical, allocatable :: t_mask_kde(:)
118 4 : integer(i4), allocatable :: time_kde(:)
119 4 : logical, allocatable :: t_mask_eval(:)
120 4 : integer(i4), allocatable :: time_eval(:)
121 : real(sp), dimension(:), allocatable :: dummy_1d_sp
122 :
123 : call get_time_indizes(time_kde, per_kde, nCalendarStepsYear)
124 : call get_time_indizes(time_eval, per_eval, nCalendarStepsYear)
125 : #ifdef SMIDEBUG
126 : print *, 'First 10 values of time_kde:'
127 : print *, time_kde(:10)
128 : print *, 'First 10 values of time_eval:'
129 : print *, time_eval(:10)
130 : #endif
131 :
132 405 : do mm = 1, nCalendarStepsYear ! time loop
133 :
134 12545307 : t_mask_kde = (time_kde .eq. mm)
135 119002 : t_mask_eval = (time_eval .eq. mm)
136 118601 : n_time = count(t_mask_eval)
137 : #ifdef SMIDEBUG
138 : print *, 'Processing: timestep, number of eval timesteps, number of kde timesteps'
139 : print *, mm, n_time, count(t_mask_kde)
140 : #endif
141 :
142 : ! check whether there is data for that timestep to be calculated
143 401 : if (n_time .eq. 0_i4) cycle
144 :
145 405 : call cellSMI(SM_kde, t_mask_kde, SM_eval, t_mask_eval, hh(:, mm), SMI)
146 : end do
147 :
148 : ! do leap days
149 4 : if (per_eval%n_leap_days .gt. 0) then
150 :
151 0 : mm = 60 ! take cdf of March first
152 0 : t_mask_kde = (time_kde .eq. mm)
153 0 : t_mask_eval = (time_eval .eq. -1_i4)
154 0 : n_time = count(t_mask_eval)
155 :
156 0 : call cellSMI(SM_kde, t_mask_kde, SM_eval, t_mask_eval, hh(:, mm), SMI)
157 : end if
158 :
159 4 : end subroutine calSMI
160 :
161 329 : subroutine cellSMI(SM_kde, t_mask_kde, SM_eval, t_mask_eval, hh, SMI)
162 :
163 : use mo_kind, only: i4, sp
164 : use mo_kernel, only: kernel_cumdensity
165 : use mo_constants, only: nodata_sp
166 :
167 : implicit none
168 :
169 : real(sp), intent(in) :: SM_kde(:, :), SM_eval(:, :)
170 : logical, intent(in) :: t_mask_kde(:), t_mask_eval(:)
171 : real(dp), intent(in) :: hh(:)
172 : real(sp), intent(inout) :: SMI(:, :)
173 :
174 : integer(i4) :: ii ! cell index
175 329 : real(dp), dimension(:), allocatable :: cdf
176 329 : real(dp), dimension(:), allocatable :: X_kde
177 : real(dp), dimension(:), allocatable :: X_eval
178 329 : real(sp), dimension(:), allocatable :: dummy_1d_sp
179 :
180 :
181 : !$OMP parallel default(shared) &
182 : !$OMP private(ii, X_kde, X_eval, cdf, dummy_1d_sp)
183 10311327 : allocate ( X_kde (count(t_mask_kde)))
184 99652 : allocate ( X_eval (count(t_mask_eval)))
185 99652 : allocate ( cdf (count(t_mask_eval)))
186 987 : allocate ( dummy_1d_sp (size(t_mask_eval, dim=1)))
187 :
188 : !$OMP do SCHEDULE(STATIC)
189 4567 : do ii = 1, size(SM_kde,1) ! cell loop
190 :
191 : #ifdef SMIDEBUG
192 : if ((mod(ii, 10000) .eq. 0_i4)) print *, ii
193 : #endif
194 :
195 123832718 : X_kde(:) = pack(real( SM_kde(ii,:), dp), t_mask_kde)
196 1274868 : X_eval(:) = pack(real(SM_eval(ii,:), dp), t_mask_eval)
197 :
198 23788 : cdf(:) = kernel_cumdensity(x_kde, hh(ii), xout=x_eval)
199 :
200 23788 : dummy_1d_sp = unpack(real(cdf(:), sp), t_mask_eval, nodata_sp)
201 1275197 : SMI(ii, :) = merge(dummy_1d_sp, SMI(ii, :), t_mask_eval)
202 :
203 : end do
204 : !$OMP end do
205 : !$OMP end parallel
206 :
207 : ! arrays are deallocated if openmp is active
208 329 : if (allocated(x_kde)) deallocate( x_kde)
209 329 : if (allocated(x_eval)) deallocate( X_eval)
210 329 : if (allocated(cdf)) deallocate( cdf)
211 329 : if (allocated(dummy_1d_sp)) deallocate( dummy_1d_sp )
212 :
213 329 : end subroutine cellSMI
214 :
215 : !> \brief create objective function of kernel_cumdensity and minimize it using
216 : !! nelmin because function is monotone
217 1 : subroutine invSMI(sm_kde, hh, SMI_invert, nCalendarStepsYear, per_kde, per_smi, SM_invert)
218 :
219 : use mo_kind, only: i4, sp, dp
220 : use mo_message, only: message
221 : use mo_constants, only: nodata_sp, nodata_dp
222 : use mo_kernel, only: kernel_cumdensity
223 :
224 : implicit none
225 :
226 : ! input variables
227 : real(dp), dimension(:,:), intent(in) :: hh
228 : real(sp), dimension(:,:), intent(in) :: sm_kde
229 : real(sp), dimension(:,:), intent(in) :: SMI_invert
230 : integer(i4), intent(in) :: nCalendarStepsYear
231 : type(period), intent(in) :: per_kde
232 : type(period), intent(in) :: per_smi
233 :
234 : ! output variables
235 : real(sp), dimension(:,:), allocatable, intent(out) :: SM_invert
236 :
237 : ! local variables
238 1 : logical, allocatable :: t_mask_kde(:), t_mask_smi(:)
239 1 : integer(i4), allocatable :: time_kde(:), time_smi(:)
240 : integer(i4) :: n_cells
241 : integer(i4) :: n_years_kde, n_years_invert
242 : integer(i4) :: ii, yy, mm ! loop variables
243 : integer(i4) :: xx_n_sample
244 : integer(i4), dimension(1) :: idx_invert
245 1 : real(sp), dimension(:), allocatable :: dummy_1d_sp
246 : real(dp) :: hh_kde
247 : real(dp) :: xx_min, xx_max, xx_h
248 1 : real(dp), dimension(:), allocatable :: y_inv, x_inv
249 1 : real(dp), dimension(:), allocatable :: xx_cdf, yy_cdf
250 1 : real(dp), dimension(:), allocatable :: xx_kde
251 :
252 : ! initialize extents
253 1 : n_cells = size(SMI_invert, 1)
254 1 : xx_n_sample = 2000_i4 ! gives precision of at least 0.0005 in SM
255 1 : allocate(xx_cdf(xx_n_sample))
256 1 : allocate(yy_cdf(xx_n_sample))
257 2001 : xx_cdf = nodata_dp
258 2001 : yy_cdf = nodata_dp
259 :
260 : ! initialize output array
261 4 : allocate(SM_invert(n_cells, size(SMI_invert, 2)))
262 116 : SM_invert = nodata_sp
263 :
264 : ! get time indizes
265 : call get_time_indizes(time_kde, per_kde, nCalendarStepsYear)
266 : call get_time_indizes(time_smi, per_smi, nCalendarStepsYear)
267 :
268 1 : call message('')
269 1 : call message(' start inversion of CDF')
270 13 : do mm = 1, nCalendarStepsYear
271 :
272 4344 : t_mask_kde = (time_kde .eq. mm)
273 84 : t_mask_smi = (time_smi .eq. mm)
274 72 : if (count(t_mask_smi) .eq. 0_i4) cycle
275 :
276 : !$OMP parallel default(shared) &
277 : !$OMP private(yy, xx_kde, hh_kde, y_inv, xx_min, xx_max, xx_h, xx_cdf, yy_cdf, idx_invert, x_inv, dummy_1d_sp)
278 1815 : allocate(xx_kde(count(t_mask_kde)))
279 40 : allocate(y_inv(count(t_mask_smi)))
280 40 : allocate(x_inv(count(t_mask_smi)))
281 40 : allocate(dummy_1d_sp(count(t_mask_smi)))
282 :
283 : !$OMP do SCHEDULE(STATIC)
284 115 : do ii = 1, n_cells
285 110 : if (modulo(ii, 1000) .eq. 0) print *, ii, n_cells
286 39710 : xx_kde(:) = pack(real(SM_kde( ii, : ), dp), t_mask_kde)
287 110 : hh_kde = hh(ii, mm)
288 660 : y_inv(:) = pack(real(SMI_invert(ii, :), dp), t_mask_smi)
289 :
290 : ! sample cdf
291 3520 : xx_min = max(0._dp, minval(xx_kde - 10._dp * hh_kde))
292 3520 : xx_max = min(1._dp, maxval(xx_kde + 10._dp * hh_kde))
293 110 : xx_h = (xx_max - xx_min) / real(xx_n_sample, dp)
294 220110 : do yy = 1, xx_n_sample
295 220110 : xx_cdf(yy) = xx_min + (yy - 1_i4) * xx_h
296 : end do
297 220110 : xx_cdf = merge(1._dp, xx_cdf, xx_cdf .gt. 1._dp)
298 220220 : yy_cdf = kernel_cumdensity(xx_kde, hh_kde, xout=xx_cdf)
299 :
300 220 : do yy = 1, size(y_inv, 1)
301 220220 : idx_invert = minloc(abs(y_inv(yy) - yy_cdf))
302 220 : x_inv(yy) = xx_cdf(idx_invert(1))
303 : end do
304 :
305 220 : dummy_1d_sp = unpack(real(x_inv(:), sp), t_mask_smi, nodata_sp)
306 665 : SM_invert(ii, :) = merge(dummy_1d_sp, SM_invert(ii, :), t_mask_smi)
307 :
308 : end do
309 : !$OMP end do
310 : !$OMP end parallel
311 :
312 5 : if (allocated(xx_kde)) deallocate(xx_kde)
313 5 : if (allocated(y_inv)) deallocate(y_inv)
314 5 : if (allocated(x_inv)) deallocate(x_inv)
315 6 : if (allocated(dummy_1d_sp)) deallocate(dummy_1d_sp)
316 :
317 : end do
318 1 : call message(' finish inversion of CDF... ok')
319 1 : call message('')
320 :
321 : ! free memory
322 1 : deallocate(xx_cdf, yy_cdf)
323 :
324 2 : end subroutine invSMI
325 :
326 :
327 11 : subroutine get_time_indizes(time, per, nCalendarStepsYear)
328 :
329 : use mo_kind, only: i4, dp
330 : use mo_smi_global_variables, only: period
331 : use mo_julian, only: dec2date, date2dec
332 :
333 : implicit none
334 :
335 : integer(i4), intent(in) :: nCalendarStepsYear
336 : type(period), intent(in) :: per
337 : integer(i4), allocatable, intent(out) :: time(:)
338 : integer(i4) :: ii, jj, start, dd, mm
339 :
340 : ! remove leap days
341 33 : allocate(time(size(per%time_points, dim=1)))
342 37174 : time(:) = 0_i4
343 :
344 11 : if (nCalendarStepsYear .eq. 12_i4) then
345 9 : start = per%m_start
346 : else
347 2 : start = int(per%j_start - date2dec(31, 12, per%y_start - 1), i4)
348 : ! account for removed leap days
349 2 : if ((int((date2dec(31, 12, per%y_start) - date2dec(1, 1, per%y_start) + 1), i4) .eq. 366) .and. per%m_start .gt. 2) then
350 0 : start = start - 1
351 : end if
352 : end if
353 :
354 :
355 11 : jj = start
356 37174 : do ii = 1, size(time)
357 37163 : time(ii) = jj
358 37163 : call dec2date(real(per%j_start + per%time_points(ii) - 1, dp), dd=dd, mm=mm)
359 37163 : if ((dd .eq. 29) .and. (mm .eq. 2) .and. (nCalendarStepsYear .ne. 12)) then
360 23 : time(ii) = -1
361 23 : cycle
362 : end if
363 37140 : jj = jj + 1_i4
364 37151 : if (jj .gt. nCalendarStepsYear) jj = 1_i4
365 : end do
366 :
367 11 : end subroutine get_time_indizes
368 :
369 : end module mo_smi
|
