CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
surface_perturbation.F90
1
4
9
10 implicit none
11
12 private
13
14 public cdfnor, ppfbet
15
16 contains
17
18! mg, sfc-perts ***
19
20! the routines below are used in the percentile matching algorithm for the
21! albedo and vegetation fraction perturbations
22
26 subroutine cdfnor(z,cdfz)
27 use machine
28
29 implicit none
30 real(kind=kind_phys), intent(out) :: cdfz
31 real(kind=kind_phys),intent(in) :: z
32! local vars
33 integer iflag
34 real(kind=kind_phys) del,x,cdfx,eps
35
36 eps = 1.0e-5
37
38
39 ! definition of passed parameters !
40 ! z = value for which the normal CDF is to be computed
41 ! eps = the absolute accuracy requirment for the CDF
42 ! iflag = error indicator on output 0->no errors, 1->errorflag from
43 ! cdfgam, 2->errorflag from cdfgam
44 ! cdfz = the CDF of the standard normal distribution evaluated at z
45
46 del = 2.0*eps
47 if (z.eq.0.0) then
48 cdfz = 0.5
49 else
50 x = 0.5*z*z
51 call cdfgam(x,0.5,del,iflag, cdfx)
52 if (iflag.ne.0) return
53 if (z.gt.0.0) then
54 cdfz = 0.5+0.5*cdfx
55 else
56 cdfz = 0.5-0.5*cdfx
57 endif
58 endif
59
60 return
61 end
62
64 subroutine cdfgam(x,alpha,eps,iflag,cdfx)
65 use machine
66
67 implicit none
68 real(kind=kind_phys), intent(out) :: cdfx
69 real(kind=kind_phys),intent(in) :: x, alpha, eps
70! local vars
71 integer iflag,i,j,k, imax
72 logical LL
73 real(kind=kind_phys) dx, dgln, p,u,epsx,pdfl, eta, bl, uflo
74 data imax, uflo / 5000, 1.0e-37 /
75
76
77 ! definition of passed parameters !
78 ! x = value for which the CDF is to be computed
79 ! alpha = parameter of gamma function (>0)
80 ! eps = the absolute accuracy requirment for the CDF
81 ! iflag = error indicator on output 0->no errors, 1->either alpha or eps
82 ! is <= oflo, 2->number of terms evaluated in the infinite series exceeds
83 ! imax.
84 ! cdf = the CDF evaluated at x
85
86 cdfx = 0.0
87
88 if (alpha.le.uflo.or.eps.le.uflo) then
89 iflag=1
90 return
91 endif
92 iflag=0
93
94 ! check for special case of x
95 if (x.le.0) return
96
97 dx = x
98 call dgamln(alpha,dgln)
99 pdfl = (alpha-1.0)*log(dx)-dx-dgln
100 if (pdfl.lt.log(uflo)) then
101 if (x.ge.alpha) cdfx = 1.0
102 else
103 p = alpha
104 u = exp(pdfl)
105 ll = .true.
106 if (x.ge.p) then
107 k = int(p)
108 if (p.le.real(k)) k = k-1
109 eta = p - real(k)
110 call dgamln(eta,dgln)
111 bl = (eta-1)*log(dx)-dx-dgln
112 ll = bl.gt.log(eps)
113 endif
114 epsx = eps/x
115 if (ll) then
116 do i=0,imax
117 if (u.le.epsx*(p-x)) return
118 u = x*u/p
119 cdfx = cdfx+u
120 p = p+1.0
121 enddo
122 iflag = 2
123 else
124 do j=1,k
125 p=p-1.0
126 if (u.le.epsx*(x-p)) continue
127 cdfx = cdfx+u
128 u = p*u/x
129 enddo
130 cdfx = 1.0-cdfx
131 endif
132 endif
133 return
134 end subroutine cdfgam
135
137 subroutine dgamln(x,dgamlnout)
138
139 use machine
140 implicit none
141 real(kind=kind_phys), intent(in) :: x
142 real(kind=kind_phys), intent(out) :: dgamlnout
143! local vars
144 integer i, n
145 real(kind=kind_phys) absacc, b1, b2, b3, b4, b5, b6, b7, b8
146 real(kind=kind_phys) c, dx, q, r, xmin, xn
147 data xmin, absacc / 6.894d0, 1.0e-15 /
148 data c / 0.918938533204672741780329736d0 /
149 data b1 / 0.833333333333333333333333333d-1 /
150 data b2 / - 0.277777777777777777777777778d-2 /
151 data b3 / 0.793650793650793650793650794d-3 /
152 data b4 / - 0.595238095238095238095238095d-3 /
153 data b5 / 0.841750841750841750841750842d-3 /
154 data b6 / - 0.191752691752691752691752692d-2 /
155 data b7 / 0.641025641025641025641025641d-2 /
156 data b8 / - 0.295506535947712418300653595d-1 /
157
158 if (x.le.0.0) stop '*** x<=0.0 in function dgamln ***'
159 dx = x
160 n = max(0,int(xmin - dx + 1.0d0) )
161 xn = dx + n
162 r = 1.0d0/xn
163 q = r*r
164 dgamlnout = r*( b1+q*( b2+q*( b3+q*( b4+q*( b5+q*( b6+q*( b7+q*b8 ) ) ) ) ) ) ) +c + (xn-0.5d0)*log(xn)-xn
165
166 if (n.gt.0) then
167 q = 1.0d0
168 do i=0, n-1
169 q = q*(dx+i)
170 enddo
171 dgamlnout = dgamlnout-log(q)
172 endif
173
174 if (dgamlnout + absacc.eq.dgamlnout) then
175 print *,' ********* WARNING FROM FUNCTION DGAMLN *********'
176 print *,' REQUIRED ABSOLUTE ACCURACY NOT ATTAINED FOR X = ',x
177 endif
178 return
179 end subroutine dgamln
180
184 subroutine ppfbet(pr,p,q,iflag,x)
185 use machine
186 implicit none
187 real(kind=kind_phys), intent(in) :: pr, p, q
188 real(kind=kind_phys), intent(out) :: x
189 ! local variables
190 integer iflag, iter, itmax
191 real(kind=kind_phys) tol, a, b, fa, fb, fc, cdf, tol1
192 real(kind=kind_phys) c, d, e, xm, s, u, v, r, eps
193 data itmax, eps / 50, 1.0e-12 /
194
195 ! Compute beta distribution value corresponding to the
196 ! probability and distribution parameters a,b.
197 !
198 ! pr - a probability value in the interval [0,1]
199 ! p - the first parameter of the beta(p,q) distribution
200 ! q - the second parameter of the beta(p,q) distribution
201 ! iflag - erro indicator in output, 0-no errors, 1,2-error flags
202 ! from subroutine cdfbet, 3- pr<0 or pr>1, 4-p<=0 or
203 ! q<=0, 5-tol<1.E-8, 6-the cdfs at the endpoints have
204 ! the same sign and no value of x is defined, 7-maximum
205 ! iterations exceeded and current value of x returned
206
207 tol = 1.0e-5
208
209
210 iflag = 0
211 if (pr.lt.0.0.or.pr.gt.1.) then
212 iflag = 3
213 return
214 endif
215 if(min(p,q).le.0.) then
216 iflag =4
217 return
218 endif
219 if (tol.lt.1.0e-8) then
220 iflag = 5
221 return
222 endif
223 a = 0.
224 b = 1.
225 fa = -pr
226 fb = 1.-pr
227 if (fb*fa.gt.0.0) then
228 iflag = 6
229 return
230 endif
231
232 fc = fb
233 do iter =1,itmax
234 if (fb*fc.gt.0.) then
235 c=a
236 fc=fa
237 d = b-a
238 e=d
239 endif
240 if (abs(fc).lt.abs(fb)) then
241 a=b
242 b=c
243 c=a
244 fa=fb
245 fb=fc
246 fc=fa
247 endif
248
249 tol1 = 2.*eps*abs(b)+0.5*tol
250 xm = 0.5*(c-b)
251 if (abs(xm).le.tol1.or.fb.eq.0.0) then
252 x=b
253 return
254 endif
255 if (abs(e).ge.tol1.and.abs(fa).gt.abs(fb)) then
256 s = fb/fa
257 if (a.eq.c) then
258 u = 2.0*xm*s
259 v = 1.0-s
260 else
261 v = fa/fc
262 r = fb/fc
263 u = s*(2.0*xm*v*(v-r)-(b-a)*(r-1.0))
264 v = (v-1.0)*(r-1.0)*(s-1.0)
265 endif
266 if (u.gt.0.0) v = -v
267 u = abs(u)
268 if (2.0*u.lt.min(3.0*xm*v-abs(tol1*v),abs(e*v))) then
269 e = d
270 d = u/v
271 else
272 d = xm
273 e = d
274 endif
275
276 else
277
278 d=xm
279 e=d
280 endif
281
282 a = b
283 fa = fb
284 if (abs(d).gt.tol1) then
285 b = b+d
286 else
287 b = b+sign(tol1,xm)
288 endif
289 call cdfbet(b,p,q,eps,iflag,cdf)
290 if (iflag.ne.0) return
291 fb = cdf-pr
292 enddo
293 x = b
294
295 return
296 end subroutine ppfbet
297
301 subroutine cdfbet(x,p,q,eps,iflag,cdfx)
302 use machine
303
304 ! Computes the value of the cumulative beta distribution at a
305 ! single point x, given the distribution parameters p,q.
306 !
307 ! x - value at which the CDF is to be computed
308 ! p - first parameter of the beta function
309 ! q - second parameter of the beta function
310 ! eps - desired absolute accuracy
311
312 implicit none
313 real(kind=kind_phys), intent(in) :: x, p, q, eps
314 real(kind=kind_phys), intent(out) :: cdfx
315 ! local vars
316 integer iflag, jmax, j
317 logical LL
318 real(kind=kind_phys) dp, dq, gamln, yxeps, w, uflo
319 real(kind=kind_phys) xy, yx, pq, qp, pdfl, u, r, v
320 real(kind=kind_phys) tmp
321 data jmax, w, uflo / 5000, 20.0, 1.0e-30 /
322
323 cdfx = 0.0
324
325 if (p.le.uflo.or.q.le.uflo.or.eps.le.uflo) then
326 iflag = 1
327 endif
328 iflag = 0
329
330 if (x.le.0.0) return
331 if (x.ge.1.0) then
332 cdfx=1.0
333 else
334 ll = (p+w).ge.(p+q+2.0*w)*x
335 if (ll) then
336 xy = x
337 yx = 1.-xy
338 pq = p
339 qp = q
340 else
341 yx = x
342 xy = 1.-yx
343 qp = p
344 pq = q
345 endif
346
347 call gmln(pq,tmp)
348 dp = (pq-1.)*log(xy)-tmp
349 call gmln(qp,tmp)
350 dq = (qp-1.)*log(yx)-tmp
351 call gmln(pq+qp,tmp)
352 pdfl = tmp+dp+dq
353
354 if (pdfl.ge.log(uflo)) then
355 u = exp(pdfl)*xy/pq
356 r = xy/yx
357 do while (qp.gt.1.)
358 if (u.le.eps*(1.-(pq+qp)*xy/(pq+1.))) then
359 if (.not.ll) cdfx = 1.-cdfx
360 return
361 endif
362 cdfx = cdfx+u
363 pq = pq+1.
364 qp = qp-1.
365 u = qp*r*u/pq
366 enddo
367 v = yx*u
368 yxeps = yx*eps
369 do j = 0, jmax
370 if (v.le.yxeps) then
371 if (.not.ll) cdfx = 1.-cdfx
372 return
373 endif
374 cdfx = cdfx + v
375 pq = pq+1.
376 v = (pq+qp-1.)*xy*v/pq
377 enddo
378 iflag = 2
379 endif
380 if (.not.ll) cdfx = 1.-cdfx
381 endif
382
383 end subroutine cdfbet
384
388 subroutine gmln(x,y)
389 use machine
390 ! Computes the natural logarithm of the gamma distribution. Users
391 ! can set the absolute accuracy and corresponding xmin.
392
393 implicit none
394 real(kind=kind_phys), intent(in) :: x
395 real(kind=kind_phys), intent(out) :: y
396! local vars
397 integer i, n
398 real(kind=kind_phys) absacc, b1, b2, b3, b4, b5, b6, b7, b8
399 real(kind=kind_phys) c, dx, q, r, xmin, xn
400! data xmin, absacc / 6.894d0, 1.0E-15 /
401 data xmin, absacc / 1.357d0, 1.0e-3 /
402 data c / 0.918938533204672741780329736d0 /
403 data b1 / 0.833333333333333333333333333d-1 /
404 data b2 / - 0.277777777777777777777777778d-2 /
405 data b3 / 0.793650793650793650793650794d-3 /
406 data b4 / - 0.595238095238095238095238095d-3 /
407 data b5 / 0.841750841750841750841750842d-3 /
408 data b6 / - 0.191752691752691752691752692d-2 /
409 data b7 / 0.641025641025641025641025641d-2 /
410 data b8 / - 0.295506535947712418300653595d-1 /
411
412 if (x.le.0.0) stop '*** x<=0.0 in function gamln ***'
413 dx = x
414 n = max(0,int(xmin - dx + 1.0d0) )
415 xn = dx + n
416 r = 1.0d0/xn
417 q = r*r
418 y = r*( b1+q*( b2+q*( b3+q*( b4+q*( b5+q*( b6+q*( b7+q*b8 ) &
419 & )) ) ) ) ) +c + (xn-0.5d0)*log(xn)-xn
420
421 if (n.gt.0) then
422 q = 1.0d0
423 do i=0, n-1
424 q = q*(dx+i)
425 enddo
426 y = y-log(q)
427 endif
428
429 if (y + absacc.eq.y) then
430 print *,' ********* WARNING FROM FUNCTION GAMLN *********'
431 print *,' REQUIRED ABSOLUTE ACCURACY NOT ATTAINED FOR X = ',x
432 endif
433 return
434 end subroutine gmln
435
436! *** mg, sfc perts
437end module surface_perturbation
subroutine dgamln(x, dgamlnout)
subroutine gmln(x, y)
This subroutine computes the natural logarithm of the gamma distribution. Users can set the absolute ...
subroutine, public ppfbet(pr, p, q, iflag, x)
This subroutine computes the beta distribution value that matches the percentile from the random patt...
subroutine cdfgam(x, alpha, eps, iflag, cdfx)
subroutine, public cdfnor(z, cdfz)
This subrtouine calculates the CDF of the standard normal distribution evaluated at z.
subroutine cdfbet(x, p, q, eps, iflag, cdfx)
This subroutine computes the value of the cumulative beta distribution at a single point x,...