CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
module_mp_radar.F90
1
3
24!+---+-----------------------------------------------------------------+
25
26 MODULE module_mp_radar
27
28 PUBLIC :: rayleigh_soak_wetgraupel
29 PUBLIC :: radar_init
30 PRIVATE :: m_complex_water_ray
31 PRIVATE :: m_complex_ice_maetzler
32 PRIVATE :: m_complex_maxwellgarnett
33 PRIVATE :: get_m_mix_nested
34 PRIVATE :: get_m_mix
35 PRIVATE :: wgamma
36 PRIVATE :: gammln
37
38
39 INTEGER, PARAMETER, PUBLIC:: nrbins = 50
40 DOUBLE PRECISION, DIMENSION(nrbins+1), PUBLIC:: xxdx
41 DOUBLE PRECISION, DIMENSION(nrbins), PUBLIC:: xxds,xdts,xxdg,xdtg
42 DOUBLE PRECISION, PARAMETER, PUBLIC:: lamda_radar = 0.10 ! in meters
43 DOUBLE PRECISION, PUBLIC:: k_w, pi5, lamda4
44 COMPLEX*16, PUBLIC:: m_w_0, m_i_0
45 DOUBLE PRECISION, DIMENSION(nrbins+1), PUBLIC:: simpson
46 DOUBLE PRECISION, DIMENSION(3), PARAMETER, PUBLIC:: basis = &
47 (/1.d0/3.d0, 4.d0/3.d0, 1.d0/3.d0/)
48 REAL, DIMENSION(4), PUBLIC:: xcre, xcse, xcge, xcrg, xcsg, xcgg
49 REAL, PUBLIC:: xam_r, xbm_r, xmu_r, xobmr
50 REAL, PUBLIC:: xam_s, xbm_s, xmu_s, xoams, xobms, xocms
51 REAL, PUBLIC:: xam_g, xbm_g, xmu_g, xoamg, xobmg, xocmg
52 REAL, PUBLIC:: xorg2, xosg2, xogg2
53
54 INTEGER, PARAMETER, PUBLIC:: slen = 20
55 CHARACTER(len=slen), PUBLIC:: &
56 mixingrulestring_s, matrixstring_s, inclusionstring_s, &
57 hoststring_s, hostmatrixstring_s, hostinclusionstring_s, &
58 mixingrulestring_g, matrixstring_g, inclusionstring_g, &
59 hoststring_g, hostmatrixstring_g, hostinclusionstring_g
60
62 DOUBLE PRECISION, PARAMETER:: melt_outside_s = 0.9d0
63 DOUBLE PRECISION, PARAMETER:: melt_outside_g = 0.9d0
64
65 CONTAINS
66
67!+---+-----------------------------------------------------------------+
68!+---+-----------------------------------------------------------------+
69!+---+-----------------------------------------------------------------+
70
72 subroutine radar_init
73
74 IMPLICIT NONE
75 INTEGER:: n
76 pi5 = 3.14159*3.14159*3.14159*3.14159*3.14159
77 lamda4 = lamda_radar*lamda_radar*lamda_radar*lamda_radar
78 m_w_0 = m_complex_water_ray(lamda_radar, 0.0d0)
79 m_i_0 = m_complex_ice_maetzler(lamda_radar, 0.0d0)
80 k_w = (abs( (m_w_0*m_w_0 - 1.0) /(m_w_0*m_w_0 + 2.0) ))**2
81
82 do n = 1, nrbins+1
83 simpson(n) = 0.0d0
84 enddo
85 do n = 1, nrbins-1, 2
86 simpson(n) = simpson(n) + basis(1)
87 simpson(n+1) = simpson(n+1) + basis(2)
88 simpson(n+2) = simpson(n+2) + basis(3)
89 enddo
90
91 do n = 1, slen
92 mixingrulestring_s(n:n) = char(0)
93 matrixstring_s(n:n) = char(0)
94 inclusionstring_s(n:n) = char(0)
95 hoststring_s(n:n) = char(0)
96 hostmatrixstring_s(n:n) = char(0)
97 hostinclusionstring_s(n:n) = char(0)
98 mixingrulestring_g(n:n) = char(0)
99 matrixstring_g(n:n) = char(0)
100 inclusionstring_g(n:n) = char(0)
101 hoststring_g(n:n) = char(0)
102 hostmatrixstring_g(n:n) = char(0)
103 hostinclusionstring_g(n:n) = char(0)
104 enddo
105
106 mixingrulestring_s = 'maxwellgarnett'
107 hoststring_s = 'air'
108 matrixstring_s = 'water'
109 inclusionstring_s = 'spheroidal'
110 hostmatrixstring_s = 'icewater'
111 hostinclusionstring_s = 'spheroidal'
112
113 mixingrulestring_g = 'maxwellgarnett'
114 hoststring_g = 'air'
115 matrixstring_g = 'water'
116 inclusionstring_g = 'spheroidal'
117 hostmatrixstring_g = 'icewater'
118 hostinclusionstring_g = 'spheroidal'
119
120!..Create bins of snow (from 100 microns up to 2 cm).
121 xxdx(1) = 100.d-6
122 xxdx(nrbins+1) = 0.02d0
123 do n = 2, nrbins
124 xxdx(n) = dexp(dfloat(n-1)/dfloat(nrbins) &
125 *dlog(xxdx(nrbins+1)/xxdx(1)) +dlog(xxdx(1)))
126 enddo
127 do n = 1, nrbins
128 xxds(n) = dsqrt(xxdx(n)*xxdx(n+1))
129 xdts(n) = xxdx(n+1) - xxdx(n)
130 enddo
131
132!..Create bins of graupel (from 100 microns up to 5 cm).
133 xxdx(1) = 100.d-6
134 xxdx(nrbins+1) = 0.05d0
135 do n = 2, nrbins
136 xxdx(n) = dexp(dfloat(n-1)/dfloat(nrbins) &
137 *dlog(xxdx(nrbins+1)/xxdx(1)) +dlog(xxdx(1)))
138 enddo
139 do n = 1, nrbins
140 xxdg(n) = dsqrt(xxdx(n)*xxdx(n+1))
141 xdtg(n) = xxdx(n+1) - xxdx(n)
142 enddo
143
144
145!..The calling program must set the m(D) relations and gamma shape
146!.. parameter mu for rain, snow, and graupel. Easily add other types
147!.. based on the template here. For majority of schemes with simpler
148!.. exponential number distribution, mu=0.
149
150 xcre(1) = 1. + xbm_r
151 xcre(2) = 1. + xmu_r
152 xcre(3) = 1. + xbm_r + xmu_r
153 xcre(4) = 1. + 2.*xbm_r + xmu_r
154 do n = 1, 4
155 xcrg(n) = wgamma(xcre(n))
156 enddo
157 xorg2 = 1./xcrg(2)
158
159 xcse(1) = 1. + xbm_s
160 xcse(2) = 1. + xmu_s
161 xcse(3) = 1. + xbm_s + xmu_s
162 xcse(4) = 1. + 2.*xbm_s + xmu_s
163 do n = 1, 4
164 xcsg(n) = wgamma(xcse(n))
165 enddo
166 xosg2 = 1./xcsg(2)
167
168 xcge(1) = 1. + xbm_g
169 xcge(2) = 1. + xmu_g
170 xcge(3) = 1. + xbm_g + xmu_g
171 xcge(4) = 1. + 2.*xbm_g + xmu_g
172 do n = 1, 4
173 xcgg(n) = wgamma(xcge(n))
174 enddo
175 xogg2 = 1./xcgg(2)
176
177 xobmr = 1./xbm_r
178 xoams = 1./xam_s
179 xobms = 1./xbm_s
180 xocms = xoams**xobms
181 xoamg = 1./xam_g
182 xobmg = 1./xbm_g
183 xocmg = xoamg**xobmg
184
185
186 end subroutine radar_init
187
188!+---+-----------------------------------------------------------------+
189!+---+-----------------------------------------------------------------+
190
195 COMPLEX*16 FUNCTION m_complex_water_ray(lambda,T)
196
197 IMPLICIT NONE
198 DOUBLE PRECISION, INTENT(IN):: t,lambda
199 DOUBLE PRECISION:: epsinf,epss,epsr,epsi
200 DOUBLE PRECISION:: alpha,lambdas,sigma,nenner
201 COMPLEX*16, PARAMETER:: i = (0d0,1d0)
202 DOUBLE PRECISION, PARAMETER:: pix=3.1415926535897932384626434d0
203
204 epsinf = 5.27137d0 + 0.02164740d0 * t - 0.00131198d0 * t*t
205 epss = 78.54d+0 * (1.0 - 4.579d-3 * (t - 25.0) &
206 + 1.190d-5 * (t - 25.0)*(t - 25.0) &
207 - 2.800d-8 * (t - 25.0)*(t - 25.0)*(t - 25.0))
208 alpha = -16.8129d0/(t+273.16) + 0.0609265d0
209 lambdas = 0.00033836d0 * exp(2513.98d0/(t+273.16)) * 1e-2
210
211 nenner = 1.d0+2.d0*(lambdas/lambda)**(1d0-alpha)*sin(alpha*pix*0.5) &
212 + (lambdas/lambda)**(2d0-2d0*alpha)
213 epsr = epsinf + ((epss-epsinf) * ((lambdas/lambda)**(1d0-alpha) &
214 * sin(alpha*pix*0.5)+1d0)) / nenner
215 epsi = ((epss-epsinf) * ((lambdas/lambda)**(1d0-alpha) &
216 * cos(alpha*pix*0.5)+0d0)) / nenner &
217 + lambda*1.25664/1.88496
218
219 m_complex_water_ray = sqrt(cmplx(epsr,-epsi))
220
221 END FUNCTION m_complex_water_ray
222
223!+---+-----------------------------------------------------------------+
224
225 COMPLEX*16 FUNCTION m_complex_ice_maetzler(lambda,T)
226
227! complex refractive index of ice as function of Temperature T
228! [deg C] and radar wavelength lambda [m]; valid for
229! lambda in [0.0001,30] m; T in [-250.0,0.0] C
230! Original comment from the Matlab-routine of Prof. Maetzler:
231! Function for calculating the relative permittivity of pure ice in
232! the microwave region, according to C. Maetzler, "Microwave
233! properties of ice and snow", in B. Schmitt et al. (eds.) Solar
234! System Ices, Astrophys. and Space Sci. Library, Vol. 227, Kluwer
235! Academic Publishers, Dordrecht, pp. 241-257 (1998). Input:
236! TK = temperature (K), range 20 to 273.15
237! f = frequency in GHz, range 0.01 to 3000
238
239 IMPLICIT NONE
240 DOUBLE PRECISION, INTENT(IN):: t,lambda
241 DOUBLE PRECISION:: f,c,tk,b1,b2,b,deltabeta,betam,beta,theta,alfa
242
243 c = 2.99d8
244 tk = t + 273.16
245 f = c / lambda * 1d-9
246
247 b1 = 0.0207
248 b2 = 1.16d-11
249 b = 335.0d0
250 deltabeta = exp(-10.02 + 0.0364*(tk-273.16))
251 betam = (b1/tk) * ( exp(b/tk) / ((exp(b/tk)-1)**2) ) + b2*f*f
252 beta = betam + deltabeta
253 theta = 300. / tk - 1.
254 alfa = (0.00504d0 + 0.0062d0*theta) * exp(-22.1d0*theta)
255 m_complex_ice_maetzler = 3.1884 + 9.1e-4*(tk-273.16)
256 m_complex_ice_maetzler = m_complex_ice_maetzler &
257 + cmplx(0.0d0, (alfa/f + beta*f))
258 m_complex_ice_maetzler = sqrt(conjg(m_complex_ice_maetzler))
259
260 END FUNCTION m_complex_ice_maetzler
261
262!+---+-----------------------------------------------------------------+
263
265 subroutine rayleigh_soak_wetgraupel (x_g, a_geo, b_geo, fmelt, &
266 meltratio_outside, m_w, m_i, lambda, C_back, &
267 mixingrule,matrix,inclusion, &
268 host,hostmatrix,hostinclusion)
269
270 IMPLICIT NONE
271
272 DOUBLE PRECISION, INTENT(in):: x_g, a_geo, b_geo, fmelt, lambda, &
273 meltratio_outside
274 DOUBLE PRECISION, INTENT(out):: c_back
275 COMPLEX*16, INTENT(in):: m_w, m_i
276 CHARACTER(len=*), INTENT(in):: mixingrule, matrix, inclusion, &
277 host, hostmatrix, hostinclusion
278
279 COMPLEX*16:: m_core, m_air
280 DOUBLE PRECISION:: d_large, d_g, rhog, x_w, xw_a, fm, fmgrenz, &
281 volg, vg, volair, volice, volwater, &
282 meltratio_outside_grenz, mra
283 INTEGER:: error
284 DOUBLE PRECISION, PARAMETER:: pix=3.1415926535897932384626434d0
285
286! refractive index of air:
287 m_air = (1.0d0,0.0d0)
288
289! Limiting the degree of melting --- for safety:
290 fm = dmax1(dmin1(fmelt, 1.0d0), 0.0d0)
291! Limiting the ratio of (melting on outside)/(melting on inside):
292 mra = dmax1(dmin1(meltratio_outside, 1.0d0), 0.0d0)
293
294! ! The relative portion of meltwater melting at outside should increase
295! ! from the given input value (between 0 and 1)
296! ! to 1 as the degree of melting approaches 1,
297! ! so that the melting particle "converges" to a water drop.
298! ! Simplest assumption is linear:
299 mra = mra + (1.0d0-mra)*fm
300
301 x_w = x_g * fm
302
303 d_g = a_geo * x_g**b_geo
304
305 if (d_g .ge. 1d-12) then
306
307 vg = pix/6. * d_g**3
308 rhog = dmax1(dmin1(x_g / vg, 900.0d0), 10.0d0)
309 vg = x_g / rhog
310
311 meltratio_outside_grenz = 1.0d0 - rhog / 1000.
312
313 if (mra .le. meltratio_outside_grenz) then
314 !..In this case, it cannot happen that, during melting, all the
315 !.. air inclusions within the ice particle get filled with
316 !.. meltwater. This only happens at the end of all melting.
317 volg = vg * (1.0d0 - mra * fm)
318
319 else
320 !..In this case, at some melting degree fm, all the air
321 !.. inclusions get filled with meltwater.
322 fmgrenz=(900.0-rhog)/(mra*900.0-rhog+900.0*rhog/1000.)
323
324 if (fm .le. fmgrenz) then
325 !.. not all air pockets are filled:
326 volg = (1.0 - mra * fm) * vg
327 else
328 !..all air pockets are filled with meltwater, now the
329 !.. entire ice sceleton melts homogeneously:
330 volg = (x_g - x_w) / 900.0 + x_w / 1000.
331 endif
332
333 endif
334
335 d_large = (6.0 / pix * volg) ** (1./3.)
336 volice = (x_g - x_w) / (volg * 900.0)
337 volwater = x_w / (1000. * volg)
338 volair = 1.0 - volice - volwater
339
340 !..complex index of refraction for the ice-air-water mixture
341 !.. of the particle:
342 m_core = get_m_mix_nested(m_air, m_i, m_w, volair, volice, &
343 volwater, mixingrule, host, matrix, inclusion, &
344 hostmatrix, hostinclusion, error)
345 if (error .ne. 0) then
346 c_back = 0.0d0
347 return
348 endif
349
350 !..Rayleigh-backscattering coefficient of melting particle:
351 c_back = (abs((m_core**2-1.0d0)/(m_core**2+2.0d0)))**2 &
352 * pi5 * d_large**6 / lamda4
353
354 else
355 c_back = 0.0d0
356 endif
357
358 end subroutine rayleigh_soak_wetgraupel
359
360!+---+-----------------------------------------------------------------+
362 complex*16 function get_m_mix_nested (m_a, m_i, m_w, volair, &
363 volice, volwater, mixingrule, host, matrix, &
364 inclusion, hostmatrix, hostinclusion, cumulerror)
365
366 IMPLICIT NONE
367
368 DOUBLE PRECISION, INTENT(in):: volice, volair, volwater
369 COMPLEX*16, INTENT(in):: m_a, m_i, m_w
370 CHARACTER(len=*), INTENT(in):: mixingrule, host, matrix, &
371 inclusion, hostmatrix, hostinclusion
372 INTEGER, INTENT(out):: cumulerror
373
374 DOUBLE PRECISION:: vol1, vol2
375 COMPLEX*16:: mtmp
376 INTEGER:: error
377
378 !..Folded: ( (m1 + m2) + m3), where m1,m2,m3 could each be
379 !.. air, ice, or water
380
381 cumulerror = 0
382 get_m_mix_nested = cmplx(1.0d0,0.0d0)
383
384 if (host .eq. 'air') then
385
386 if (matrix .eq. 'air') then
387 write(*,*) 'GET_M_MIX_NESTED: bad matrix: ', matrix
388 cumulerror = cumulerror + 1
389 else
390 vol1 = volice / max(volice+volwater,1d-10)
391 vol2 = 1.0d0 - vol1
392 mtmp = get_m_mix(m_a, m_i, m_w, 0.0d0, vol1, vol2, &
393 mixingrule, matrix, inclusion, error)
394 cumulerror = cumulerror + error
395
396 if (hostmatrix .eq. 'air') then
397 get_m_mix_nested = get_m_mix(m_a, mtmp, 2.0*m_a, &
398 volair, (1.0d0-volair), 0.0d0, mixingrule, &
399 hostmatrix, hostinclusion, error)
400 cumulerror = cumulerror + error
401 elseif (hostmatrix .eq. 'icewater') then
402 get_m_mix_nested = get_m_mix(m_a, mtmp, 2.0*m_a, &
403 volair, (1.0d0-volair), 0.0d0, mixingrule, &
404 'ice', hostinclusion, error)
405 cumulerror = cumulerror + error
406 else
407 write(*,*) 'GET_M_MIX_NESTED: bad hostmatrix: ', &
408 hostmatrix
409 cumulerror = cumulerror + 1
410 endif
411 endif
412
413 elseif (host .eq. 'ice') then
414
415 if (matrix .eq. 'ice') then
416 write(*,*) 'GET_M_MIX_NESTED: bad matrix: ', matrix
417 cumulerror = cumulerror + 1
418 else
419 vol1 = volair / max(volair+volwater,1d-10)
420 vol2 = 1.0d0 - vol1
421 mtmp = get_m_mix(m_a, m_i, m_w, vol1, 0.0d0, vol2, &
422 mixingrule, matrix, inclusion, error)
423 cumulerror = cumulerror + error
424
425 if (hostmatrix .eq. 'ice') then
426 get_m_mix_nested = get_m_mix(mtmp, m_i, 2.0*m_a, &
427 (1.0d0-volice), volice, 0.0d0, mixingrule, &
428 hostmatrix, hostinclusion, error)
429 cumulerror = cumulerror + error
430 elseif (hostmatrix .eq. 'airwater') then
431 get_m_mix_nested = get_m_mix(mtmp, m_i, 2.0*m_a, &
432 (1.0d0-volice), volice, 0.0d0, mixingrule, &
433 'air', hostinclusion, error)
434 cumulerror = cumulerror + error
435 else
436 write(*,*) 'GET_M_MIX_NESTED: bad hostmatrix: ', &
437 hostmatrix
438 cumulerror = cumulerror + 1
439 endif
440 endif
441
442 elseif (host .eq. 'water') then
443
444 if (matrix .eq. 'water') then
445 write(*,*) 'GET_M_MIX_NESTED: bad matrix: ', matrix
446 cumulerror = cumulerror + 1
447 else
448 vol1 = volair / max(volice+volair,1d-10)
449 vol2 = 1.0d0 - vol1
450 mtmp = get_m_mix(m_a, m_i, m_w, vol1, vol2, 0.0d0, &
451 mixingrule, matrix, inclusion, error)
452 cumulerror = cumulerror + error
453
454 if (hostmatrix .eq. 'water') then
455 get_m_mix_nested = get_m_mix(2*m_a, mtmp, m_w, &
456 0.0d0, (1.0d0-volwater), volwater, mixingrule, &
457 hostmatrix, hostinclusion, error)
458 cumulerror = cumulerror + error
459 elseif (hostmatrix .eq. 'airice') then
460 get_m_mix_nested = get_m_mix(2*m_a, mtmp, m_w, &
461 0.0d0, (1.0d0-volwater), volwater, mixingrule, &
462 'ice', hostinclusion, error)
463 cumulerror = cumulerror + error
464 else
465 write(*,*) 'GET_M_MIX_NESTED: bad hostmatrix: ', &
466 hostmatrix
467 cumulerror = cumulerror + 1
468 endif
469 endif
470
471 elseif (host .eq. 'none') then
472
473 get_m_mix_nested = get_m_mix(m_a, m_i, m_w, &
474 volair, volice, volwater, mixingrule, &
475 matrix, inclusion, error)
476 cumulerror = cumulerror + error
477
478 else
479 write(*,*) 'GET_M_MIX_NESTED: unknown matrix: ', host
480 cumulerror = cumulerror + 1
481 endif
482
483 IF (cumulerror .ne. 0) THEN
484 write(*,*) 'GET_M_MIX_NESTED: error encountered'
485 get_m_mix_nested = cmplx(1.0d0,0.0d0)
486 endif
487
488 end function get_m_mix_nested
489
490!+---+-----------------------------------------------------------------+
492 COMPLEX*16 FUNCTION get_m_mix (m_a, m_i, m_w, volair, volice, &
493 volwater, mixingrule, matrix, inclusion, error)
494
495 IMPLICIT NONE
496
497 DOUBLE PRECISION, INTENT(in):: volice, volair, volwater
498 COMPLEX*16, INTENT(in):: m_a, m_i, m_w
499 CHARACTER(len=*), INTENT(in):: mixingrule, matrix, inclusion
500 INTEGER, INTENT(out):: error
501
502 error = 0
503 get_m_mix = cmplx(1.0d0,0.0d0)
504
505 if (mixingrule .eq. 'maxwellgarnett') then
506 if (matrix .eq. 'ice') then
507 get_m_mix = m_complex_maxwellgarnett(volice, volair, volwater, &
508 m_i, m_a, m_w, inclusion, error)
509 elseif (matrix .eq. 'water') then
510 get_m_mix = m_complex_maxwellgarnett(volwater, volair, volice, &
511 m_w, m_a, m_i, inclusion, error)
512 elseif (matrix .eq. 'air') then
513 get_m_mix = m_complex_maxwellgarnett(volair, volwater, volice, &
514 m_a, m_w, m_i, inclusion, error)
515 else
516 write(*,*) 'GET_M_MIX: unknown matrix: ', matrix
517 error = 1
518 endif
519
520 else
521 write(*,*) 'GET_M_MIX: unknown mixingrule: ', mixingrule
522 error = 2
523 endif
524
525 if (error .ne. 0) then
526 write(*,*) 'GET_M_MIX: error encountered'
527 endif
528
529 END FUNCTION get_m_mix
530
531!+---+-----------------------------------------------------------------+
533 COMPLEX*16 FUNCTION m_complex_maxwellgarnett(vol1, vol2, vol3, &
534 m1, m2, m3, inclusion, error)
535
536 IMPLICIT NONE
537
538 COMPLEX*16 :: m1, m2, m3
539 DOUBLE PRECISION :: vol1, vol2, vol3
540 CHARACTER(len=*) :: inclusion
541
542 COMPLEX*16 :: beta2, beta3, m1t, m2t, m3t
543 INTEGER, INTENT(out) :: error
544
545 error = 0
546
547 if (dabs(vol1+vol2+vol3-1.0d0) .gt. 1d-6) then
548 write(*,*) 'M_COMPLEX_MAXWELLGARNETT: sum of the ', &
549 'partial volume fractions is not 1...ERROR'
550 m_complex_maxwellgarnett=cmplx(-999.99d0,-999.99d0)
551 error = 1
552 return
553 endif
554
555 m1t = m1**2
556 m2t = m2**2
557 m3t = m3**2
558
559 if (inclusion .eq. 'spherical') then
560 beta2 = 3.0d0*m1t/(m2t+2.0d0*m1t)
561 beta3 = 3.0d0*m1t/(m3t+2.0d0*m1t)
562 elseif (inclusion .eq. 'spheroidal') then
563 beta2 = 2.0d0*m1t/(m2t-m1t) * (m2t/(m2t-m1t)*log(m2t/m1t)-1.0d0)
564 beta3 = 2.0d0*m1t/(m3t-m1t) * (m3t/(m3t-m1t)*log(m3t/m1t)-1.0d0)
565 else
566 write(*,*) 'M_COMPLEX_MAXWELLGARNETT: ', &
567 'unknown inclusion: ', inclusion
568 m_complex_maxwellgarnett=dcmplx(-999.99d0,-999.99d0)
569 error = 1
570 return
571 endif
572
573 m_complex_maxwellgarnett = &
574 sqrt(((1.0d0-vol2-vol3)*m1t + vol2*beta2*m2t + vol3*beta3*m3t) / &
575 (1.0d0-vol2-vol3+vol2*beta2+vol3*beta3))
576
577 END FUNCTION m_complex_maxwellgarnett
578
579!+---+-----------------------------------------------------------------+
581 REAL function gammln(xx)
582! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0.
583 IMPLICIT NONE
584 REAL, INTENT(IN):: xx
585 DOUBLE PRECISION, PARAMETER:: stp = 2.5066282746310005d0
586 DOUBLE PRECISION, DIMENSION(6), PARAMETER:: &
587 cof = (/76.18009172947146d0, -86.50532032941677d0, &
588 24.01409824083091d0, -1.231739572450155d0, &
589 .1208650973866179d-2, -.5395239384953d-5/)
590 DOUBLE PRECISION:: ser,tmp,x,y
591 INTEGER:: j
592
593 x=xx
594 y=x
595 tmp=x+5.5d0
596 tmp=(x+0.5d0)*log(tmp)-tmp
597 ser=1.000000000190015d0
598 DO 11 j=1,6
599 y=y+1.d0
600 ser=ser+cof(j)/y
60111 CONTINUE
602 gammln=tmp+log(stp*ser/x)
603 END FUNCTION gammln
604! (C) Copr. 1986-92 Numerical Recipes Software 2.02
605!+---+-----------------------------------------------------------------+
607 REAL function wgamma(y)
608
609 IMPLICIT NONE
610 REAL, INTENT(IN):: y
611
612 wgamma = exp(gammln(y))
613
614 END FUNCTION wgamma
615
616!+---+-----------------------------------------------------------------+
617 END MODULE module_mp_radar
618!+---+-----------------------------------------------------------------+
error
subroutine error(parameters, swdown,fsa,fsr,fira,fsh,fcev, fgev,fctr,ssoil,beg_wb,canliq,canice, sneqv,wa,smc,dzsnso,prcp,ecan, etran,edir,runsrf,runsub,dt,nsoil, nsnow,ist,errwat, iloc,jloc,fveg, sav,sag,fsrv,fsrg,zwt,pah, ifdef ccpp
check surface energy balance and water balance.
Definition module_sf_noahmplsm.F90:1469
module_mp_radar
This module is more library code whereas the individual microphysics schemes contains specific detail...
Definition module_mp_radar.F90:26