Line data Source code
1 : !> \file mo_edk_setvario.f90
2 : !> \brief \copybrief mo_edk_setvario
3 : !> \details \copydetails mo_edk_setvario
4 :
5 : !> \brief VARIOGRAM: Seting or estimating and fitting
6 : !> \details PURPOSE:
7 : !! 1) Set variagram from DB for a block, or
8 : !! 2.1) Estimate an empirical semi variogram for daily met. values
9 : !! 2.2) Fitting a teoretical variogram
10 : !! 2.3) Set variogram for a block
11 : !> \author Luis Samaniego
12 : !> \date 19.02.2004
13 : !! - main structure
14 : !> \date 12.04.2006
15 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
16 : !! EDK is released under the LGPLv3+ license \license_note
17 : module mo_edk_setvario
18 : use mo_kind, only: i4, dp
19 : implicit none
20 :
21 : private
22 :
23 : public :: setVario
24 : public :: dMatrix
25 : public :: tVar
26 :
27 : contains
28 :
29 : !> \brief find optimal variogram parameters, if wanted
30 2 : subroutine setVario(param)
31 : use runControl
32 : use VarFit
33 : use mainVar
34 : use mo_edk_empvar, only: EmpVar
35 : implicit none
36 : real(dp), intent(out) :: param(3) !< parameters (nugget, sill, range)
37 : integer(i4) :: jd, y, y0
38 : !
39 : ! Estimation
40 2 : if (flagVario) then
41 : ! Initialize
42 1 : nobs=0 ; m0=0.0_dp ; v0=0.0_dp
43 1 : y0 = 0
44 : !
45 366 : do jd= jStart, jEnd
46 365 : y = floor(float(jd-jStart)/365.)
47 365 : if (y > y0 ) then
48 0 : y0 = y
49 0 : print*, 'VarFit. Processing ', yStart+y0-1 , '...'
50 : end if
51 : ! update Variogram bins daily
52 366 : call EmpVar(jd, .true.)
53 : end do
54 : !
55 : ! variance
56 1 : v0 = v0 / real(nobs, dp)
57 : ! optimize parameters
58 : ! open(UNIT = 6,FILE = 'Report_OPT.sol', STATUS='REPLACE')
59 :
60 1 : print *, "ST: replace old GRG2 opti with something better"
61 1 : call opti(param)
62 : else
63 1 : print *, " ... no variogram estimation, parameters given."
64 : end if
65 2 : end subroutine setVario
66 :
67 : !> \brief Function to calulate variogram at given distance
68 : !> \return variogram value
69 651348 : real(dp) function tVar(h,c0,c,a)
70 2 : use VarFit, only : vType
71 : use mo_utils, only: eq
72 : real(dp), intent(in) :: h !< distance
73 : real(dp), intent(in) :: c0 !< nugget = beta(1) = XU(1)
74 : real(dp), intent(in) :: c !< sill = beta(2) = XU(2)
75 : real(dp), intent(in) :: a !< range = beta(3) = XU(3)
76 651348 : real(dp) :: r
77 : !
78 651348 : if ( a > 0.0_dp ) then
79 651348 : r = h/a
80 : else
81 : r = 0.0_dp
82 : end if
83 328500 : select case (vType)
84 : !
85 : case (1)
86 : ! composed: nugget + spherical + sill
87 328500 : if ( eq(h, 0.0_dp) ) then
88 33984 : tVar = c0 ! 0.0_dp
89 294516 : elseif ( h < a ) then
90 240642 : tVar = c0 + c * (1.5_dp * r - 0.5_dp * r**3)
91 : else
92 53874 : tVar = c0 + c
93 : end if
94 : !
95 : case (2)
96 : ! composed: nugget + exponential + sill
97 651348 : if ( eq(h, 0.0_dp) ) then
98 33984 : tVar = c0 ! 0.0_dp
99 288864 : elseif ( a > 0 ) then
100 288864 : tVar = c0 + c * (1.0_dp - exp(-r))
101 : else
102 0 : tVar = c0 + c
103 : end if
104 : end select
105 : !
106 651348 : end function tVar
107 :
108 : !> \brief DMATRIX:: To calculate distance between pairs.
109 2 : subroutine dMatrix
110 651348 : use mainVar
111 : use kriging
112 : use runControl
113 :
114 : implicit none
115 : integer(i4) :: i, j, k
116 : integer(i4) :: r, c, ii, jj
117 : integer(i4) :: delta, nTcell
118 2 : real(dp) :: xc, yc
119 2 : integer(i4), allocatable :: list(:)
120 : !
121 : ! Initialize variables
122 0 : if ( allocated(cell)) deallocate (cell)
123 : !
124 2 : print*, nCell, "cells", nSta, "stations"
125 2 : edk_dist%ncell = nCell
126 6 : allocate(edk_dist%cell_pos(nCell, 2))
127 126 : allocate ( cell(nCell) )
128 6 : allocate ( list(nSta) )
129 :
130 38 : do i=1,nSta-1
131 : ! distance matrix between stations: checked OK
132 380 : do j=i+1, nSta
133 342 : if (edk_dist%getSS(i,j) == 0.0_dp) then
134 : ! check if stations are closer than 5 meter
135 0 : print* , '--------------------------------------------------------------------------'
136 0 : print* , '!!! Warning: '
137 0 : print* , '!!! Stations: ', MetSta(i)%Id, ' and ', MetSta(j)%Id
138 0 : print* , '!!! have the same coordinates or are repeated. Check LUT. '
139 0 : print* , '!!! Rounded artefacts can be generated when stations have same coordinates'
140 0 : print* , '!!! and data at the time.'
141 0 : print* , '--------------------------------------------------------------------------'
142 : !stop
143 : end if
144 378 : if (edk_dist%getSS(i,j) > 0.0_dp .and. edk_dist%getSS(i,j) < 5.0_dp) then
145 0 : print* , '--------------------------------------------------------------------------'
146 0 : print* , '!!! Warning: '
147 0 : print* , '!!! Stations: ', MetSta(i)%Id, ' and ', MetSta(j)%Id
148 0 : print* , '!!! are closer than 5 meter distance. '
149 0 : print* , '--------------------------------------------------------------------------'
150 : end if
151 : end do
152 : end do
153 :
154 : ! cell coordinates and elevation : checked OK
155 : ! ***************************************
156 : ! cell numbering convention (1DIM first)
157 : ! c->1 2 ncol
158 : ! r
159 : ! ---+-----+------...+...-----+
160 : ! 1 nr+1
161 : ! 2
162 : ! ... k
163 : ! nr 2nr nCell
164 : ! ***************************************
165 : ! ii column in finer grid
166 : ! yy row in finer grid
167 : ! (xc,yc) coordinates of meteogrid
168 : ! ***************************************
169 : !
170 2 : r=1
171 2 : c=0
172 2 : if (DEMNcFlag /= 1) xc = gridMeteo%xllcorner + dble(gridMeteo%cellsize) * 0.5_dp
173 2 : delta = cellFactor / 2
174 2 : jj = delta
175 122 : do k=1,nCell
176 : ! advancing the counters
177 120 : if (r == 1) then
178 12 : c = c + 1
179 12 : if (c > 1) then
180 10 : jj = jj + cellFactor
181 : end if
182 : ii = delta
183 : else
184 108 : ii = ii + cellFactor
185 : end if
186 :
187 120 : if (DEMNcFlag == 1) then
188 0 : cell(k)%x = gridMeteo%easting(r,c)
189 0 : cell(k)%y = gridMeteo%northing(r,c)
190 : else
191 120 : if (r == 1) then
192 12 : if (c > 1) then
193 10 : xc = xc + dble(gridMeteo%cellsize)
194 : end if
195 12 : yc = gridMeteo%yllcorner + dble(gridMeteo%cellsize) * (dble(gridMeteo%nrows) - 0.5_dp)
196 : else
197 108 : yc = yc - dble(gridMeteo%cellsize)
198 : end if
199 120 : cell(k)%x = xc
200 120 : cell(k)%y = yc
201 : end if
202 360 : edk_dist%cell_pos(k,:) = [cell(k)%x, cell(k)%y]
203 :
204 : ! average of only four DEM cells around centre cell (from lower grid scale upto higher grid cell)
205 : !cell(k)%h = 0.25_dp*(G(ii,jj)%h + G(ii,jj+1)%h + G(ii+1,jj)%h + G(ii+1,jj+1)%h)
206 : !
207 : ! average of all DEM cells (from lower grid scale upto higher grid cell)
208 120 : if (cellFactor == 1) then
209 0 : cell(k)%h = G(ii+1,jj+1)%h
210 : else
211 196920 : nTcell = count(G( (ii-delta+1):(ii+delta) , (jj-delta+1):(jj+delta) )%h > grid%nodata_value )
212 120 : if (nTcell == 0) then
213 48 : cell(k)%h = gridMeteo%nodata_value
214 : else
215 0 : cell(k)%h = sum(G( (ii-delta+1):(ii+delta) , (jj-delta+1):(jj+delta) )%h, &
216 118152 : G( (ii-delta+1):(ii+delta) , (jj-delta+1):(jj+delta) )%h /= gridMeteo%nodata_value ) / dble(nTcell)
217 : end if
218 : end if
219 :
220 : ! advance the counters
221 120 : r=r+1
222 122 : if (r > gridMeteo%nrows) r = 1
223 : end do
224 :
225 2 : print*, " ... find neighborhoods"
226 :
227 : ! find the closest stations to cell i (any order): checked OK
228 122 : do i=1,nCell
229 120 : if ( modulo(i,100000) == 0 ) print*, " ... cells ready", i, "of", nCell
230 2400 : list = -9
231 2400 : do j=1,nSta
232 2400 : if (edk_dist%getCS(i,j) <= maxDist) list(j) = j
233 : end do
234 2400 : cell(i)%nNS = count(list > -9)
235 360 : allocate ( cell(i)%listNS( cell(i)%nNS ) )
236 4802 : cell(i)%listNS = pack(list, MASK = list >-9)
237 : end do
238 : !
239 2 : deallocate (list)
240 :
241 2 : end subroutine dMatrix
242 :
243 : !> \brief Optimization routine for variograms.
244 : !> \details Initialization of the Nonlinear Optimization Subroutine GRG2
245 : !! The function to be optimized is suplied in subroutine GCOMP
246 : !> \author Luis Samaniego
247 : !> \date 28.05.1999
248 1 : subroutine OPTI(pmin)
249 2 : use VarFit
250 : use mo_nelmin, only: nelminrange
251 :
252 : ! parameters for Nelder-Mead algorithm
253 : real(dp), intent(out) :: pmin(3) !< optimized parameter set
254 4 : real(dp) :: pstart(3) ! Starting point for the iteration.
255 9 : real(dp) :: prange(3, 2) ! Range of parameters (upper and lower bound).
256 1 : real(dp) :: varmin ! the terminating limit for the variance of the function values. varmin>0 is required
257 4 : real(dp) :: step(3) ! determines the size and shape of the initial simplex. The relative magnitudes of its elements should reflect the units of the variables. size(step)=size(start)
258 : integer(i4) :: konvge ! the convergence check is carried out every konvge iterations
259 : integer(i4) :: maxeval ! the maximum number of function evaluations. default: 1000
260 1 : real(dp) :: funcmin
261 : integer(i4) :: neval ! the number of function evaluations used.
262 : integer(i4) :: numrestart ! the number of restarts.
263 : integer(i4) :: ierror ! error indicator.
264 : ! 0: no errors detected.
265 : ! 1: varmin or konvge have an illegal value.
266 : ! 2: iteration terminated because maxeval was exceeded without convergence.
267 1 : real(dp), allocatable :: history(:)
268 :
269 :
270 : ! ! inputs for GRG2
271 : ! IMPLICIT DOUBLE PRECISION(A-H,O-Z), INTEGER(I,J,L,M,N)
272 : ! INTEGER*4 NCORE,NNVARS,NFUN,MAXBAS,MAXHES,LASTZ
273 : ! INTEGER*4 M,N,MP1,NPMP1,NBMAX,NNBMAX,NPNBMX,MPNBMX,NRTOT
274 : ! LOGICAL MAXIM,INPRNT,OTPRNT
275 : ! DIMENSION BLVAR(100),BUVAR(100),BLCON(10),BUCON(10)
276 : ! DIMENSION RAMCON(10),RAMVAR(10),INBIND(10),Z(5000)
277 : ! DIMENSION NONBAS(10),REDGR(10),DEFAUL(20),TITLE(19)
278 : ! DIMENSION XX(100),FCNS(10),RMULTS(10)
279 : ! DATA BIG/1.D31/
280 :
281 : ! Initialization of Nelder-Mead
282 1 : pstart = (/0.0, 1., 0.5/) ! Starting point for the iteration.
283 4 : prange(:, 1) = (/0., 0., 0./) ! Range of parameters (lower bound).
284 4 : prange(:, 2) = (/0.3, 5., 2./) ! Range of parameters (upper bound).
285 1 : varmin = 0.001 ! the terminating limit for the variance of the function values. varmin>0 is required
286 1 : step = (/0.15, 2.5, 1./) ! determines the size and shape of the initial simplex. The relative magnitudes of its elements should reflect the units of the variables. size(step)=size(start)
287 1 : konvge = 100 ! the convergence check is carried out every konvge iterations
288 1 : maxeval = 2000 ! the maximum number of function evaluations. default: 1000
289 :
290 : ! Call Nelder-Mead optimizer to reduce GCOMP
291 : ! pmin = nelmin(obj_func, pstart, varmin, step, konvge, maxeval, &
292 : ! funcmin, neval, numrestart, ierror, history)
293 : pmin = nelminrange(obj_func, pstart, prange, varmin, step, konvge, maxeval, &
294 1 : funcmin, neval, numrestart, ierror)
295 :
296 : ! scale up distance h
297 151 : where (gamma(:,1) > 0._dp) gamma(:,1) = gamma(:,1) * gmax(1)
298 :
299 : ! scale back parameters (only for range a)
300 1 : pmin(3)=pmin(3)*gmax(1)
301 5 : beta = pmin
302 :
303 1 : print *, '=============================='
304 1 : print *, ' Results of Nelder-Mead optimization '
305 1 : print *, neval, ' of ', maxeval
306 1 : print *, "funcmin: ", funcmin
307 1 : print *, "p_obj: ", pmin
308 1 : print *, 'error: ', ierror
309 1 : print *, 'varmin: ', varmin
310 : if (allocated(history)) print *, 'history: ', history(1), history(size(history))
311 1 : print *, 'gmax: ', gmax(1)
312 :
313 :
314 : ! estimate statistics
315 1 : call stats
316 : !
317 : ! scale up distance h
318 151 : where (gamma(:,1) > 0._dp) gamma(:,1) = gamma(:,1) * gmax(1)
319 :
320 1 : end subroutine OPTI
321 :
322 : !> \brief Function to be minimised for the nelder mead algorithm
323 : !> \return Target function value for given parameter set.
324 193 : function obj_func(p)
325 1 : use varfit, only : nbins, gamma, nh
326 : implicit none
327 : real(dp), dimension(:), intent(in) :: p
328 : real(dp) :: obj_func
329 193 : real(dp) :: gcal
330 : integer(i4) :: k
331 : !
332 193 : obj_func = 0.0_dp
333 : !
334 29143 : do k=1,nbins
335 29143 : if (gamma(k,1) > 0.0_dp) then
336 : ! gcal = tvar (gamma(k,1), xu(1), xu(2), xu(3))
337 5597 : gcal = tvar(gamma(k,1), p(1), p(2), p(3))
338 : !
339 : ! estimator l1 weighted
340 5597 : if (gamma(k,1) <= 1.0_dp ) obj_func = obj_func + dabs(gamma(k,2) - gcal) * dble(Nh(k))
341 : end if
342 : end do
343 193 : end function obj_func
344 :
345 : !> \brief statistics
346 1 : subroutine stats
347 193 : use varFit, only : E, beta, gamma
348 : implicit none
349 : integer(i4), parameter :: incx = 1
350 : integer(i4) :: k
351 : integer(i4) :: ne
352 1 : real(dp) :: SSE
353 1 : real(dp) :: zObsMean, zCalMean
354 1 : real(dp) :: zObsVar, zCalVar, sumP, NSE_denom
355 1 : real(dp), dimension(:), allocatable :: error, denom
356 1 : real(dp), dimension(:), allocatable :: zCal, zObs
357 : real(dp), parameter :: small = -9.999d3
358 :
359 : !
360 : !Initialize
361 151 : ne = count(gamma(:,1) > 0.0_dp)
362 6 : allocate (error(ne), denom(ne), zCal(ne), zObs(ne))
363 31 : zObs = gamma(1:ne,2)
364 30 : zCal = 0.0_dp
365 30 : do k=1,ne
366 30 : if (gamma(k,1) > 0._dp) zCal(k) = tvar(gamma(k,1),beta(1),beta(2),beta(3))
367 : end do
368 : !
369 31 : error = zObs-zCal
370 1 : if ( ne > 1 ) then
371 30 : zObsMean = sum(zObs)/real(ne, dp)
372 30 : zCalMean = sum(zCal)/real(ne, dp)
373 30 : sumP = dot_product(zObs,zCal)
374 30 : zObsVar = dot_product(zObs,zObs) - real(ne, dp) * zObsMean * zObsMean
375 30 : zCalVar = dot_product(zCal,zCal) - real(ne, dp) * zCalMean * zCalMean
376 30 : SSE = dot_product(error,error)
377 31 : denom = zObs - zObsMean
378 30 : NSE_denom = dot_product(denom,denom)
379 : else
380 : zObsMean = small
381 : zCalMean = small
382 : zObsVar = small
383 : zCalVar = small
384 : end if
385 : ! ****************
386 : ! Quality measures
387 : ! ****************
388 1 : if ( ne > 0 ) then
389 : !
390 : ! BIAS
391 1 : E(1) = zCalMean - zObsMean
392 1 : print *, 'BIAS: ', E(1)
393 : !
394 : ! MSE
395 1 : E(2) = SSE/real(ne, dp)
396 1 : print *, 'MSE: ', E(2)
397 : !
398 : ! RMSE
399 1 : if ( E(2) > 0.0_dp ) then
400 1 : E(3) = dsqrt(E(2))
401 : else
402 0 : E(3) = small
403 : end if
404 1 : print *, 'RMSE: ', E(3)
405 : !
406 : ! RRMSE
407 1 : if ( E(3) > 0.0_dp ) then
408 1 : E(4)=E(3)/zObsMean
409 : else
410 0 : E(4)= small
411 : end if
412 1 : print *, 'PRMSE:', E(4)
413 : !
414 : ! MAE
415 30 : E(5)= sum(abs(error))
416 1 : E(5)= E(5)/real(ne, dp)
417 1 : print *, 'MAE: ', E(5)
418 : !
419 : ! RMAE
420 1 : E(6)=E(5)/zObsMean
421 1 : print *, 'RMAE: ', E(6)
422 : !
423 : ! r
424 1 : E(7)= (sumP-real(ne, dp) * zCalMean * zObsMean) / dsqrt(zCalVar * zObsVar)
425 1 : print *, 'r: ', E(7)
426 : !
427 : ! NSE
428 1 : E(8)= 1.0_dp - (SSE/NSE_denom)
429 1 : print *, 'NSE: ', E(8)
430 : else
431 0 : E = small
432 : end if
433 1 : deallocate (error, denom, zCal, zObs)
434 1 : end subroutine stats
435 :
436 : end module mo_edk_setvario
|
