1MODULE localizationFunction_mod
2 ! MODULE localizationFunction_mod (prefix='lfn' category='2. B and R matrices')
3 !
4 !:Purpose: To store various functions for horizontal and vertical covariance
5 ! localization. Length-scale estimation for function available in
6 ! this module is also possible through curve fitting using
7 ! pre-computed optimal separation-distance localization values.
8 !
9 use earthConstants_mod
10 use utilities_mod
11 implicit none
12 save
13 private
14
15 ! public procedures
16 public :: lfn_Setup, lfn_Response, lfn_lengthscale, lfn_createBiPerFunction
17
18 logical :: initialized = .false.
19
20 character(len=128) :: LocFunction
21
22CONTAINS
23
24 !--------------------------------------------------------------------------
25 ! lfn_Setup
26 !--------------------------------------------------------------------------
27 subroutine lfn_Setup(LocFunctionWanted)
28 implicit none
29
30 ! Arguments:
31 character(len=*), intent(in) :: LocFunctionWanted
32
33 select case(trim(LocFunctionWanted))
34 case ('FifthOrder')
35 write(*,*)
36 write(*,*) 'lfn_setup: Using Gaspari and Cohn 5th order piecewise rationale function'
37 LocFunction='FifthOrder'
38 case default
39 write(*,*)
40 write(*,*) 'Unsupported HORIZONTAL localization function : ', trim(LocFunctionWanted)
41 call utl_abort('lfn_setup')
42 end select
43
44 initialized = .true.
45
46 end subroutine lfn_Setup
47
48 !--------------------------------------------------------------------------
49 ! lfn_Response
50 !--------------------------------------------------------------------------
51 function lfn_Response(distance,lengthscale) result(correlation)
52 implicit none
53
54 ! Arguments:
55 real(8), intent(in) :: distance
56 real(8), intent(in) :: lengthscale
57 ! Result:
58 real(8) :: correlation
59
60 if ( .not. initialized ) then
61 write(*,*)
62 write(*,*) 'The localixation function module was NOT initialized'
63 call utl_abort('lfn_horizResponse')
64 endif
65
66 select case(trim(LocFunction))
67 case ('FifthOrder')
68 correlation=lfn_FifthOrderFunction(distance,lengthscale)
69 case default
70 write(*,*)
71 write(*,*) 'Unsupported localization function : ', trim(LocFunction)
72 call utl_abort('lfn_Response')
73 end select
74
75 end function lfn_Response
76
77 !--------------------------------------------------------------------------
78 ! lfn_gradient
79 !--------------------------------------------------------------------------
80 function lfn_gradient(distance,lengthscale) result(gradient)
81 implicit none
82
83 ! Arguments:
84 real(8), intent(in) :: distance
85 real(8), intent(in) :: lengthscale
86 ! Result:
87 real(8) :: gradient
88
89 if ( .not. initialized ) then
90 write(*,*)
91 write(*,*) 'The localixation function module was NOT initialized'
92 call utl_abort('lfn_gradient')
93 endif
94
95 select case(trim(LocFunction))
96 case ('FifthOrder')
97 gradient=lfn_FifthOrderFunctionGradient(distance,lengthscale)
98 case default
99 write(*,*)
100 write(*,*) 'Unsupported localization function : ', trim(LocFunction)
101 call utl_abort('lfn_gradient')
102 end select
103
104 end function lfn_gradient
105
106 !--------------------------------------------------------------------------
107 ! lfn_FifthOrderFunction
108 !--------------------------------------------------------------------------
109 function lfn_FifthOrderFunction(distance,lengthscale) result(correlation)
110 implicit none
111
112 ! Arguments:
113 real(8), intent(in) :: distance
114 real(8), intent(in) :: lengthscale
115 ! Result:
116 real(8) :: correlation
117
118 ! Locals:
119 real(8) :: halflength
120
121 halflength = lengthscale/2.d0
122
123 if ( distance <= halflength ) then
124 correlation = -0.250d0*(distance/halflength)**5 &
125 + 0.5d0*(distance/halflength)**4 &
126 + 0.625d0*(distance/halflength)**3 &
127 - (5.0d0/3.0d0)*(distance/halflength)**2 &
128 + 1.0d0
129 else if ( distance <= (2.0d0*halflength) ) then
130 correlation = (1.0d0/12.0d0)*(distance/halflength)**5 &
131 - 0.5d0*(distance/halflength)**4 &
132 + 0.625d0*(distance/halflength)**3 &
133 + (5.0d0/3.0d0)*(distance/halflength)**2 &
134 - 5.0d0*(distance/halflength) &
135 + 4.0d0 &
136 - (2.0d0/3.0d0)*(halflength/distance)
137 else
138 correlation = 0.d0
139 endif
140
141 end function lfn_FifthOrderFunction
142
143 !--------------------------------------------------------------------------
144 ! lfn_FifthOrderFunctionGradient
145 !--------------------------------------------------------------------------
146 function lfn_FifthOrderFunctionGradient(distance,lengthscale) result(gradient)
147 implicit none
148
149 ! Arguments:
150 real(8), intent(in) :: distance
151 real(8), intent(in) :: lengthscale
152 ! Result:
153 real(8) :: gradient
154
155 ! Locals:
156 real(8) :: halflength
157
158 halflength = lengthscale/2.d0
159
160 ! Compute df/dhalflength
161 if ( distance <= halflength ) then
162 gradient = 5.d0*0.250d0*(distance/halflength)**5 / halflength &
163 - 4.d0*0.5d0*(distance/halflength)**4 / halflength &
164 - 3.d0*0.625d0*(distance/halflength)**3 / halflength &
165 + (10.0d0/3.0d0)*(distance/halflength)**2 / halflength
166 else if ( distance <= (2.0d0*halflength) ) then
167 gradient = (-5.0d0/12.0d0)*(distance/halflength)**5 / halflength &
168 + 4.0d0*0.5d0*(distance/halflength)**4 / halflength &
169 - 3.0d0*0.625d0*(distance/halflength)**3 / halflength &
170 - (10.0d0/3.0d0)*(distance/halflength)**2 / halflength &
171 + 5.0d0*(distance/halflength) / halflength &
172 - (2.0d0/3.0d0)*(1.0d0/distance)
173 else
174 gradient = 0.d0
175 endif
176
177 ! Compute df/dlengthscale (= df/dhalflength * dhalflength/dlengthscale)
178 gradient = gradient / 2.d0
179
180 end function lfn_FifthOrderFunctionGradient
181
182!--------------------------------------------------------------------------
183! lfn_CreateBiPerFunction
184!--------------------------------------------------------------------------
185 SUBROUTINE lfn_CreateBiPerFunction(gridpoint,CorrelLength,dlon,ni,nj,nk)
186 implicit none
187
188 ! Arguments:
189 integer, intent(in) :: ni
190 integer, intent(in) :: nj
191 integer, intent(in) :: nk
192 real(8), intent(in) :: dlon ! in radian
193 real(8), intent(in) :: CorrelLength(nk) ! in km
194 real(8), intent(out) :: gridpoint(ni,nj,nk)
195
196 ! Locals:
197 integer :: i, j, k, iref, jref
198 real(8) :: distance, distance_ref
199
200 gridpoint(:,:,:) = 0.d0
201
202 distance_ref = dlon * ec_ra
203
204 !
205 !- Create a bi-periodic correlation function by centering the function in each 4 corners
206 !
207
208 !- Lower-Left Corner (reference corner)
209 iref = 1
210 jref = 1
211 do j = 1, nj
212 do i = 1, ni
213 distance = distance_ref * sqrt( real((i-iref)**2 + (j-jref)**2,8) )
214 do k = 1, nk
215 gridpoint(i,j,k) = gridpoint(i,j,k) + lfn_response(distance,1000.d0*CorrelLength(k))
216 enddo
217 enddo
218 enddo
219
220 !- Upper-Left Corner
221 iref = 1
222 jref = nj+1 ! the lam spectral transform is periodic at nj+1
223 do j = 1, nj
224 do i = 1, ni
225 distance = distance_ref * sqrt( real((i-iref)**2 + (j-jref)**2,8) )
226 do k = 1, nk
227 gridpoint(i,j,k) = gridpoint(i,j,k) + lfn_response(distance,1000.d0*CorrelLength(k))
228 enddo
229 enddo
230 enddo
231
232 !- Lower-Right Corner
233 iref = ni+1 ! the lam spectral transform is periodic at ni+1
234 jref = 1
235 do j = 1, nj
236 do i = 1, ni
237 distance = distance_ref * sqrt( real((i-iref)**2 + (j-jref)**2,8) )
238 do k = 1, nk
239 gridpoint(i,j,k) = gridpoint(i,j,k) + lfn_response(distance,1000.d0*CorrelLength(k))
240 enddo
241 enddo
242 enddo
243
244 !- Upper-Right Corner
245 iref = ni+1 ! the lam spectral transform is periodic at ni+1
246 jref = nj+1 ! the lam spectral transform is periodic at nj+1
247 do j = 1, nj
248 do i = 1, ni
249 distance = distance_ref * sqrt( real((i-iref)**2 + (j-jref)**2,8) )
250 do k = 1, nk
251 gridpoint(i,j,k) = gridpoint(i,j,k) + lfn_response(distance,1000.d0*CorrelLength(k))
252 enddo
253 enddo
254 enddo
255
256 END SUBROUTINE Lfn_CreateBiPerFunction
257
258 !--------------------------------------------------------------------------
259 ! LFN_LENGTHSCALE
260 !--------------------------------------------------------------------------
261 subroutine lfn_lengthscale(lengthscale, rmse, localizationValues, &
262 distance, weight, numbins)
263 implicit none
264
265 ! Arguments:
266 integer, intent(in) :: numBins
267 real(8), intent(in) :: localizationValues(numBins)
268 real(8), intent(in) :: distance(numBins)
269 real(8), intent(in) :: weight(numBins)
270 real(8), intent(inout) :: lengthscale
271 real(8), intent(out) :: rmse
272
273 ! Locals:
274 integer :: ierr, ilist
275 real(8) :: minv, fit
276
277 ! 1. Initialization
278 ilist = 1
279 ierr = 0
280 minv = 0.01
281
282 ! 2. Apply Curve fitting algorithm
283 call lfn_curveFit(numBins, distance, & ! IN
284 localizationValues, & ! IN
285 weight, numBins, & ! IN
286 lengthscale, & ! INOUT
287 minv, ilist, & ! IN
288 ierr, & ! INOUT
289 fit) ! OUT
290
291 ! 3. root mean squared deviation output as 'degree of fit'
292 rmse = sqrt(fit)
293
294 end subroutine lfn_lengthscale
295
296 !--------------------------------------------------------------------------
297 ! LFN_CURVEFIT
298 !--------------------------------------------------------------------------
299 subroutine lfn_curveFit(Nn, x, y, w, nmax, param, minv, ilist, ierr, ssq)
300 !
301 !:Purpose: This subroutine computes the lengthscale (sl, here: param) of a
302 ! FUNCTION that best fits (in a least-square perspective) a
303 ! data set. The method follows routine CURVEFIT from:
304 ! Heeswijk, M.V., and C.G. Fox, 1988: Iterative Method and Fortran
305 ! Code for Nonlinear Curve Fitting, Computers and Geosciences, 14,
306 ! 4, pp. 489-503.
307 !
308 implicit none
309
310 ! Arguments:
311 INTEGER, INTENT(IN) :: Nn
312 REAL(8), INTENT(IN) :: x(Nn)
313 REAL(8), INTENT(IN) :: y(Nn)
314 REAL(8), INTENT(IN) :: w(Nn)
315 INTEGER, INTENT(IN) :: nmax
316 REAL(8), INTENT(INOUT) :: param
317 REAL(8), INTENT(IN) :: minv
318 INTEGER, INTENT(IN) :: ilist
319 INTEGER, INTENT(INOUT) :: ierr
320 REAL(8), INTENT(OUT) :: ssq
321
322 ! Locals:
323 INTEGER, PARAMETER :: icon = 100 ! max iteration
324 INTEGER, PARAMETER :: lbis = 10 ! max bisection
325 REAL(8) :: d(Nn)
326 REAL(8) :: g(Nn)
327 REAL(8) :: m
328 REAL(8) :: gtd
329 REAL(8) :: gtg
330 REAL(8) :: gtginv
331 REAL(8) :: ssqold, rootmsq
332 REAL(8) :: w_t
333 INTEGER :: iter,nbis,i ! Loop counters.
334 CHARACTER(LEN=60) :: FMT1
335 CHARACTER(LEN=12), PARAMETER :: RoutineName = 'lfn_curveFit'
336
337 !-------------------------------------------------------------------------------
338
339 !-------------------------------------------------------------------------------
340 ! 1. Initialize variables
341 !-------------------------------------------------------------------------------
342 ierr = 0 ! Error flags
343 iter = 0 ! Iteration counter
344 nbis = 0 ! Bisection counter
345 m = 9.9E9 ! dummy value to initialize m
346 ssqold = 9.9E9 ! dummy value to initialize ssqold
347
348 !-------------------------------------------------------------------------------
349 ! 2. Perform curve fitting using a minimization (least-square) approach
350 !-------------------------------------------------------------------------------
351 !PRINT*, RoutineName, ' fitting ', trim(locFunction)
352
353 DO WHILE ( (ABS(m) > minv) .AND. (iter <= icon) .AND. (nbis <= lbis) )
354
355 ! 2.1 Form vector D; i.e. a data point minus the isolated Taylor Series
356 ! term for that data point. In addition calculate the sum of the
357 ! square of the residuals
358 ssq = 0.d0
359 DO i = 1, nmax
360 d(i) = y(i) - lfn_response(x(i), param)
361 ssq = ssq + d(i)**2 * w(i)
362 END DO
363 rootmsq = sqrt(ssq)
364
365 ! 2.2 If we are converging
366 IF (ssq < ssqold) THEN
367
368 ! The matrix formation in the following program segment follows
369 ! equation 3.12 on page 41 of Geophysical Data Analysis: Discrete
370 ! Inverse Theory by Menke (1984)
371
372 ! 2.2.1 Form matrix G
373 DO i = 1, nmax
374 g(i)=lfn_gradient(x(i), param) ! curve derivative wrt to param
375 END DO
376
377 ! 2.2.2 Form matrix GTD
378 gtd = 0.d0
379 DO i = 1, nmax
380 gtd = gtd + d(i) * g(i) * w(i)
381 END DO
382
383 ! 2.2.3 Form matrix GTG
384 gtg = 0.d0
385 DO i = 1, nmax
386 gtg = gtg + g(i) * g(i) * w(i)
387 END DO
388
389 ! 2.2.4 Find the inverse of matrix GTG
390 if (gtg /= 0.d0) then
391 gtginv = 1.d0 / gtg
392 else
393 write(*,*) 'Lfn_Curvefit: gtg = 0 ! Aborting the iteration process'
394 DO i = 1, nmax
395 write(*,*) g(i), w(i), x(i), param
396 END DO
397 gtginv = tiny(1.d0)
398 !call utl_abort('Lfn_Curvefit')
399 end if
400
401 ! 2.2.5 Find vector M (i.e. increment of PARAM)
402 m = gtginv * gtd
403
404 ssqold = ssq ! keep
405
406 iter = iter + 1
407
408 ! 2.3 If we are diverging
409 ELSE
410
411 m = m / 2.d0 ! Bisect the correction
412 nbis = nbis + 1
413
414 END IF
415
416 ! 2.4 Update the parameter
417 param = param + m
418
419 ! 2.5 Print some info
420 FMT1 = "(A12,2X,A2,I6,2X,A2,I3,2X,A2,G10.4,2X,A5,G10.4,2X,A2,G10.4)"
421
422 IF (ilist == 1) WRITE(*,FMT1) RoutineName,'i=',iter,'n=',nbis,'m=',m, &
423 'rmsd=',rootmsq,'s=',param
424
425 END DO
426
427 PRINT*
428 PRINT*," Estimated lengthscale = ",param
429
430 !-------------------------------------------------------------------------------
431 ! 3. Compute the final mean squared deviation
432 !-------------------------------------------------------------------------------
433 ssq = 0.d0
434 w_t = 0.d0
435 DO i = 1, nmax
436 d(i) = y(i) - lfn_response(x(i), param)
437 w_t = w_t + w(i)
438 ssq = ssq + w(i) * d(i)**2
439 END DO
440 ssq = ssq / w_t
441
442 end SUBROUTINE Lfn_Curvefit
443
444END MODULE localizationFunction_mod
