Line data Source code
1 : !> \file mo_edk.f90
2 : !> \brief \copybrief mo_edk
3 : !> \details \copydetails mo_edk
4 :
5 : !> \brief Main module for EDK.
6 : !> \details Executes the EDK setup.
7 : !> \author Luis Samaniego
8 : !> \date 22.03.2006
9 : !> \date 14.04.2006
10 : !> \date 14.07.2010
11 : !! - DEM ckeck
12 : !> \author Matthias Zink
13 : !> \date 05.05.2011
14 : !! - IWD if No stations < 2
15 : !> \date 07.02.2012
16 : !! - correct prec data < 0 --> = 0
17 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
18 : !! EDK is released under the LGPLv3+ license \license_note
19 : module mo_edk
20 :
21 : implicit none
22 :
23 : private
24 :
25 : public :: EDK, clean, WriteDataMeteo
26 :
27 : contains
28 :
29 : !> \brief External Drift Kriging
30 : !> \details Perform EDK (daily)
31 : !! Output: output gridsize = cellFactor * gridsize DEM
32 2664 : subroutine EDK(k, jStart, jEnd, edk_dist, MetSta, cell, doOK)
33 : use mo_kind, only : i4, dp, sp
34 : use mainVar, only : MeteoStation, noDataValue, thresholdDist
35 : use kriging, only : CellCoarser
36 : use runControl, only : correctNeg, distZero
37 : use mo_edk_types, only: dist_t
38 :
39 : implicit none
40 :
41 : ! input / output variables
42 : integer(i4), intent(in) :: k !< cell id
43 : integer(i4), intent(in) :: jStart, jEnd
44 : type(dist_t), intent(in) :: edk_dist !< distances matrix
45 : type(MeteoStation), intent(in) :: MetSta(:) !< MeteoStation input
46 : type(CellCoarser), intent(inout) :: cell !< cell specification
47 : logical, optional, intent(in) :: doOK !< switch do ordinary kriging
48 :
49 : ! local variables
50 : logical :: doOK_loc, calc_weights
51 : integer(i4) :: jd ! julian day
52 : integer(i4) :: i, j
53 : integer(i4) :: n_select, n_zero
54 : integer(i4) :: ii, jj
55 0 : logical :: selectNS(cell%nNS) ! selected neighborhood (no no-data-values) current time-step
56 0 : logical :: selectNS_old(cell%nNS) ! selected neighborhood (no no-data-values) previous time-step
57 2664 : integer(i4), allocatable :: selectNS_ids(:)
58 2664 : real(dp), allocatable :: lamda(:)
59 2664 : real(dp), allocatable :: weights(:)
60 : ! real(dp), allocatable :: weights_old(:)
61 2664 : real(dp) :: sumLamda
62 :
63 53280 : selectNS_old = .false.
64 :
65 : ! switch ordinary kriging off if not explicitly given
66 2664 : doOK_loc = .False.
67 2664 : if (present(doOK)) doOK_loc = doOK
68 : ! IF NK changed -> re-estimate weights
69 28944 : timeloop: do jd = jStart, jEnd
70 525600 : selectNS = .false. ! reset selected neighborhood
71 26280 : n_select = 0 ! counter for no-data-values in current Neighborhood
72 26280 : n_zero = 0 ! counter for zero data values in current Neighborhood
73 525600 : do i = 1, cell%nNS
74 499320 : j = cell%listNS(i)
75 525600 : if ( MetSta(j)%z(jd) /= noDataValue ) then
76 :
77 : !*** zero field ***********************************
78 420480 : if (MetSta(j)%z(jd) == 0.0_dp ) n_zero = n_zero + 1
79 : !**********************************************************
80 420480 : n_select = n_select + 1
81 420480 : selectNS(i) = .true.
82 : end if
83 : end do
84 : ! list of IDs of selected neighbors
85 26280 : if ( allocated(selectNS_ids) ) deallocate(selectNS_ids)
86 78840 : allocate(selectNS_ids(n_select))
87 26280 : selectNS_ids = pack(cell%listNS, mask=selectNS)
88 :
89 : ! check of selected stations are same, weights are allocated and edk was executed in the previous time step (edk_true = .True.)
90 360072 : if (all(selectNS .eqv. selectNS_old) .and. allocated(weights)) then
91 : calc_weights = .False. ! same Neighborhood
92 : else
93 8712 : calc_weights = .True. ! new Neighborhood (new stations or old station with missing data)
94 : end if
95 :
96 28944 : if (.not. ( n_zero == n_select .or. n_select == 1 .or. n_select == 2 ) ) then
97 : ! no value ! avoid indetermination
98 : ! avoid 0 value calculations n_zero == n_select
99 : ! avoid calculations where only 1 station is available
100 : ! avoid numerical instabilities n_select == 2 (may happen that the solver matrix becomes singular)
101 : ! if (jd > jStart) weights_old = weights
102 15624 : if (calc_weights) then
103 4248 : if ( allocated(weights) ) deallocate(weights)
104 4248 : call get_kriging_weights(weights, n_select, selectNS_ids, doOK_loc, edk_dist, k, cell%h, MetSta)
105 52128 : if (jd > jStart) selectNS_old = selectNS ! copy previous neighborhood to variable if weights are calculated
106 : end if
107 :
108 : !if ((jd > jStart) .and. (sum(weights) /= sum(weights_old)) .and. all(selectNS .eqv. selectNS_old) .and. (calc_weights .eqv. .False.)) then
109 : !print *, sum(weights), ' - ', sum(weights_old)
110 : !print *, (sum(weights) /= sum(weights_old))
111 : !print *, "weights not equal but neighborhood is same at"
112 : !print *, 'time:', jd, 'cell: ', k
113 : !print *, all(selectNS .eqv. selectNS_old), allocated(weights)
114 : !print *, selectNS
115 : !print *, selectNS_old
116 : !end if
117 : ! print *, 'sum weights: ', sum(weights), maxval(weights), maxloc(weights), calc_weights
118 : ! The BLUE of z is then:
119 15624 : cell%z(jd) = 0.
120 265608 : do i=1, n_select
121 249984 : ii = selectNS_ids(i)
122 265608 : cell%z(jd) = cell%z(jd) + real(weights(i) * MetSta(ii)%z(jd), sp)
123 : end do
124 :
125 10656 : else if (n_zero /= n_select .AND. n_select == 1) then
126 : ! only one station available /= 0
127 0 : ii = selectNS_ids(1)
128 0 : cell%z(jd) = real(MetSta(ii)%z(jd), sp)
129 : ! TODO: this is actually wrong (should be also calculated by distance)
130 :
131 : ! for precipitation, distant values are set to zero
132 0 : if (distZero .and. edk_dist%getCS(k,ii) .gt. thresholdDist) then
133 0 : cell%z(jd) = 0.0_sp
134 : end if
135 :
136 10656 : else if (n_zero /= n_select .AND. n_select == 2) then
137 : ! avoid numerical instabilities --> IWD insverse weighted squared distance
138 : ! matrix of solver may become instable if only two or three stations are available
139 0 : allocate (lamda(n_select))
140 0 : do i=1, n_select
141 0 : ii=selectNS_ids(i)
142 0 : lamda(i)=1/edk_dist%getCS(k,ii)/edk_dist%getCS(k,ii)
143 : end do
144 0 : sumLamda=sum(lamda)
145 0 : lamda=lamda/sum(lamda)
146 0 : cell%z = 0.0_dp
147 0 : do i=1, n_select
148 0 : ii=selectNS_ids(i)
149 0 : cell%z(jd) = cell%z(jd) + real(lamda(i) * MetSta(ii)%z(jd), sp)
150 : end do
151 0 : deallocate (lamda)
152 :
153 10656 : else if (n_zero == n_select) then ! also for empty neighborhood
154 : ! if all stations have the value 0
155 10656 : cell%z(jd) = 0.0_sp
156 :
157 : end if
158 : end do timeloop
159 :
160 : ! correct negative
161 28944 : if (correctNeg) cell%z = merge(0._sp, cell%z, (cell%z .gt. -9999._sp) .and. (cell%z .lt. 0.))
162 : !
163 2664 : end subroutine EDK
164 :
165 : !> \brief get kriging weights
166 8496 : subroutine get_kriging_weights(X, nNmax, Nk, doOK_loc, edk_dist, k, cell_h, MetSta)
167 :
168 2664 : use mo_kind, only : dp, i4
169 : use mainVar, only : MeteoStation
170 : use varfit, only : beta
171 : use mo_edk_setvario, only : tVar
172 : use mo_edk_types, only : dist_t
173 :
174 : implicit none
175 :
176 : real(dp), allocatable, intent(out) :: X(:)
177 : integer(i4), intent(in) :: nNmax
178 : integer(i4), intent(in) :: Nk(:)
179 : logical, intent(in) :: doOK_loc
180 : type(dist_t), intent(in) :: edk_dist !< distances
181 : integer(i4), intent(in) :: k !< cell index
182 : real(dp), intent(in) :: cell_h !< cell elevation
183 : type(MeteoStation), intent(in) :: MetSta(:) !< MeteoStation input
184 :
185 : ! local variables
186 : integer(i4) :: i, j, ii, jj, info
187 4248 : integer(i4), allocatable :: ipvt(:)
188 4248 : real(dp), allocatable :: A (:,:), B(:), C(:,:)
189 :
190 : !
191 : ! initialize matrices
192 4248 : if (doOK_loc) then
193 0 : allocate (A(nNmax+1,nNmax+1), B(nNmax+1), X(nNmax+1), ipvt(nNmax + 1))
194 : else
195 42480 : allocate (A(nNmax+2,nNmax+2), B(nNmax+2), X(nNmax+2), ipvt(nNmax + 2))
196 : end if
197 : !
198 : ! assembling the system of equations: OK
199 1457064 : A=0.0_dp
200 :
201 : ! loop over available stations nNmax in neighborhood
202 72216 : do i = 1, nNmax
203 67968 : ii = Nk(i)
204 :
205 577728 : do j = i + 1, nNmax
206 509760 : jj = Nk(j)
207 :
208 : ! available only the upper triangular matrix
209 577728 : if (jj > ii) then
210 509760 : A(i,j)=tVar(edk_dist%getSS(ii,jj),beta(1),beta(2),beta(3))
211 509760 : A(j,i)=A(i,j)
212 : else
213 0 : A(i,j)=tVar(edk_dist%getSS(jj,ii),beta(1),beta(2),beta(3))
214 0 : A(j,i)=A(i,j)
215 : end if
216 : end do
217 67968 : A(i,i) = tVar(0._dp, beta(1),beta(2),beta(3))
218 :
219 67968 : A(i,nNmax+1) = 1.0_dp
220 67968 : A(nNmax+1,i) = 1.0_dp
221 67968 : if (.not. doOK_loc) A(i,nNmax+2) = MetSta(ii)%h
222 67968 : if (.not. doOK_loc) A(nNmax+2,i) = MetSta(ii)%h
223 :
224 72216 : B(i)=tVar(edk_dist%getCS(k,ii), beta(1),beta(2),beta(3))
225 : end do
226 :
227 : !
228 4248 : B(nNmax+1) = 1.0_dp
229 4248 : if (.not. doOK_loc) B(nNmax+2) = cell_h
230 :
231 : ! NOTE: only the upper triangular matrix is needed!
232 : ! call D_LSASF (A, B, X)
233 : ! print *, '<<<<<<<<<<<<<<'
234 : ! print *, 'rhs = ', B(:10)
235 1469808 : C = A
236 84960 : X = B
237 : ! print *, 'A = ', C(:10, 1)
238 : ! print *, '<<<<<<<<<<<<<<'
239 4248 : call dgesv(size(A, 1), 1, A, size(A, 1), ipvt, B, size(A, 1), info)
240 : ! print *, '<<<<<<<<<<<<<<'
241 : ! print *, 'B = ', B(:10)
242 : ! print *, 'A = ', A(:10, 1)
243 : ! print *, '<<<<<<<<<<<<<<'
244 84960 : if (maxval(abs(matmul(C, B) - X)) .gt. 1e-10) print *, 'maximum error: ', maxval(abs(matmul(C, B) - X))
245 84960 : X = B
246 4248 : if (info .ne. 0_i4) then
247 0 : print *, '***WARNING: calculation of weights failed'
248 : end if
249 : ! print *, 'shape of A: ', shape(A)
250 : ! print *, 'dem of stations: ', A(:, 10)
251 : ! print *, 'maximum of dem at stations: ', maxval(MetSta(:)%h)
252 : ! print *, 'easting: ', cell%x
253 : ! print *, 'northing: ', cell%y
254 : ! print *, 'number of neighbors: ', nNmax
255 : ! print *, ''
256 : ! ! print *, 'ipvt: ', ipvt
257 : ! print *, 'info: ', info
258 72216 : if (abs(sum(X(:nNmax)) - 1._dp) .gt. 1.e-4) then
259 0 : print *, '***WARNING: sum of weights is not 1, calculation of weights failed'
260 0 : print *, 'sum of weights: ', sum(X(:nNmax))
261 : ! stop 'testing'
262 : end if
263 4248 : deallocate (A, B, C, ipvt)
264 4248 : end subroutine get_kriging_weights
265 :
266 : !> \brief free allocated space of all arrays
267 2 : subroutine clean
268 4248 : use mainVar
269 : use mo_kind, only: i4
270 : use kriging
271 : use runControl
272 : integer(i4) :: i
273 : !
274 : ! DEM will be reused...
275 : !
276 :
277 : ! Stations
278 40 : do i = 1, nSta
279 40 : if ( allocated( MetSta(i)%z ) ) deallocate( MetSta(i)%z )
280 : end do
281 40 : if ( allocated(MetSta) ) deallocate (MetSta)
282 :
283 122 : do i=1,nCell
284 122 : if ( allocated( cell(i)%listNS ) ) deallocate ( cell(i)%listNS )
285 : end do
286 122 : if ( allocated(cell)) deallocate (cell)
287 :
288 2 : call edk_dist%clean
289 :
290 2 : end subroutine clean
291 :
292 : !> \brief WRITE METEREOLOGIC VARIABLES
293 1 : subroutine WriteDataMeteo
294 2 : use mo_kind, only : i4, dp
295 : use mainVar
296 : use runControl
297 : use kriging
298 : use VarFit
299 : use mo_edk_setvario, only : tvar
300 : !
301 : implicit none
302 : integer(i4) :: i, j, k
303 : integer(i4) :: leap ! leap day either 0 or 1
304 : character(256) :: dummy
305 : character(256) :: fileName
306 : logical :: wasOpened
307 :
308 : ! print varfit cross-val depending on variogram estimation (true or false)
309 1 : if ( flagVario ) then !
310 1 : fileName=trim(dataPathOut)//'varFit.txt'
311 1 : inquire(201, OPENED = wasOpened)
312 1 : if (.not.wasOpened) then
313 1 : open (201, file=filename, status='unknown', action='write')
314 1 : write (201, 200) 'nugget', 'sill', 'range', 'Type', 'Easting', 'Northing','BIAS', 'RMSE', 'r'
315 : end if
316 4 : write (201, 201) (beta(i),i=1,nParam), vType, (xl+xr)*0.5_dp, (yd+yu)*0.5_dp, E(1), E(2), E(7)
317 1 : close(201)
318 :
319 1 : fileName = trim(dataPathOut)//trim('variogram.dat')
320 1 : open (21, file=trim(fileName), status='unknown', action='write')
321 1 : write(21,203)
322 151 : do k=1,nbins
323 209 : if (nh(k) >0) write(21,204) (gamma(k,j), j=1,2), tVar(gamma(k,1),beta(1),beta(2),beta(3)), nh(k)
324 : end do
325 1 : close(21)
326 : end if
327 :
328 : !
329 : ! formats
330 : 110 format (i4,'.bin')
331 : 200 format (3a12 , a6, 2a11 , 3a9)
332 : 201 format (3e12.4, i6, 2f11.1, 3f9.4)
333 : 203 format (19x,'h',12x,'gamma(h)',12x,'g_cal(h)', 6x, 'N(h)')
334 : 204 format (3es20.5,i10)
335 :
336 : !
337 1 : end subroutine WriteDataMeteo
338 :
339 : end module mo_edk
|
