phys/module_mp_radar.F

   1 !+---+-----------------------------------------------------------------+
   2 !..This set of routines facilitates computing radar reflectivity.
   3 !.. This module is more library code whereas the individual microphysics
   4 !.. schemes contains specific details needed for the final computation,
   5 !.. so refer to location within each schemes calling the routine named
   6 !.. rayleigh_soak_wetgraupel.
   7 !.. The bulk of this code originated from Ulrich Blahak (Germany) and
   8 !.. was adapted to WRF by G. Thompson.  This version of code is only
   9 !.. intended for use when Rayleigh scattering principles dominate and
  10 !.. is not intended for wavelengths in which Mie scattering is a
  11 !.. significant portion.  Therefore, it is well-suited to use with
  12 !.. 5 or 10 cm wavelength like USA NEXRAD radars.
  13 !.. This code makes some rather simple assumptions about water
  14 !.. coating on outside of frozen species (snow/graupel).  Fraction of
  15 !.. meltwater is simply the ratio of mixing ratio below melting level
  16 !.. divided by mixing ratio at level just above highest T>0C.  Also,
  17 !.. immediately 90% of the melted water exists on the ice's surface
  18 !.. and 10% is embedded within ice.  No water is "shed" at all in these
  19 !.. assumptions. The code is quite slow because it does the reflectivity
  20 !.. calculations based on 50 individual size bins of the distributions.
  21 !+---+-----------------------------------------------------------------+
  22
  23 MODULE module_mp_radar
  24
  25       USE module_wrf_error
  26
  27       PUBLIC :: rayleigh_soak_wetgraupel
  28       PUBLIC :: radar_init
  29       PRIVATE :: m_complex_water_ray
  30       PRIVATE :: m_complex_ice_maetzler
  31       PRIVATE :: m_complex_maxwellgarnett
  32       PRIVATE :: get_m_mix_nested
  33       PRIVATE :: get_m_mix
  34       PRIVATE :: WGAMMA
  35       PRIVATE :: GAMMLN
  36
  37
  38       INTEGER, PARAMETER, PUBLIC:: nrbins = 50
  39       DOUBLE PRECISION, DIMENSION(nrbins+1), PUBLIC:: xxDx
  40       DOUBLE PRECISION, DIMENSION(nrbins), PUBLIC:: xxDs,xdts,xxDg,xdtg
  41       DOUBLE PRECISION, PARAMETER, PUBLIC:: lamda_radar = 0.10           ! in meters
  42       DOUBLE PRECISION, PUBLIC:: K_w, PI5, lamda4
  43       COMPLEX*16, PUBLIC:: m_w_0, m_i_0
  44       DOUBLE PRECISION, DIMENSION(nrbins+1), PUBLIC:: simpson
  45       DOUBLE PRECISION, DIMENSION(3), PARAMETER, PUBLIC:: basis =       &
  46                            (/1.d0/3.d0, 4.d0/3.d0, 1.d0/3.d0/)
  47       REAL, DIMENSION(4), PUBLIC:: xcre, xcse, xcge, xcrg, xcsg, xcgg
  48       REAL, PUBLIC:: xam_r, xbm_r, xmu_r, xobmr
  49       REAL, PUBLIC:: xam_s, xbm_s, xmu_s, xoams, xobms, xocms
  50       REAL, PUBLIC:: xam_g, xbm_g, xmu_g, xoamg, xobmg, xocmg
  51       REAL, PUBLIC:: xorg2, xosg2, xogg2
  52
  53       INTEGER, PARAMETER, PUBLIC:: slen = 20
  54       CHARACTER(len=slen), PUBLIC::                                     &
  55               mixingrulestring_s, matrixstring_s, inclusionstring_s,    &
  56               hoststring_s, hostmatrixstring_s, hostinclusionstring_s,  &
  57               mixingrulestring_g, matrixstring_g, inclusionstring_g,    &
  58               hoststring_g, hostmatrixstring_g, hostinclusionstring_g
  59
  60 !..Single melting snow/graupel particle 90% meltwater on external sfc
  61       DOUBLE PRECISION, PARAMETER:: melt_outside_s = 0.9d0
  62       DOUBLE PRECISION, PARAMETER:: melt_outside_g = 0.9d0
  63
  64       CHARACTER*256:: radar_debug
  65
  66 CONTAINS
  67
  68 !+---+-----------------------------------------------------------------+
  69 !+---+-----------------------------------------------------------------+
  70 !+---+-----------------------------------------------------------------+
  71
  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
 191       COMPLEX*16 FUNCTION m_complex_water_ray(lambda,T)
 192
 193 !      Complex refractive Index of Water as function of Temperature T
 194 !      [deg C] and radar wavelength lambda [m]; valid for
 195 !      lambda in [0.001,1.0] m; T in [-10.0,30.0] deg C
 196 !      after Ray (1972)
 197
 198       IMPLICIT NONE
 199       DOUBLE PRECISION, INTENT(IN):: T,lambda
 200       DOUBLE PRECISION:: epsinf,epss,epsr,epsi
 201       DOUBLE PRECISION:: alpha,lambdas,sigma,nenner
 202       COMPLEX*16, PARAMETER:: i = (0d0,1d0)
 203       DOUBLE PRECISION, PARAMETER:: PIx=3.1415926535897932384626434d0
 204
 205       epsinf  = 5.27137d0 + 0.02164740d0 * T - 0.00131198d0 * T*T
 206       epss    = 78.54d+0 * (1.0 - 4.579d-3 * (T - 25.0)                 &
 207               + 1.190d-5 * (T - 25.0)*(T - 25.0)                        &
 208               - 2.800d-8 * (T - 25.0)*(T - 25.0)*(T - 25.0))
 209       alpha   = -16.8129d0/(T+273.16) + 0.0609265d0
 210       lambdas = 0.00033836d0 * exp(2513.98d0/(T+273.16)) * 1e-2
 211
 212       nenner = 1.d0+2.d0*(lambdas/lambda)**(1d0-alpha)*sin(alpha*PIx*0.5) &
 213              + (lambdas/lambda)**(2d0-2d0*alpha)
 214       epsr = epsinf + ((epss-epsinf) * ((lambdas/lambda)**(1d0-alpha)   &
 215            * sin(alpha*PIx*0.5)+1d0)) / nenner
 216       epsi = ((epss-epsinf) * ((lambdas/lambda)**(1d0-alpha)            &
 217            * cos(alpha*PIx*0.5)+0d0)) / nenner                           &
 218            + lambda*1.25664/1.88496
 219
 220       m_complex_water_ray = SQRT(CMPLX(epsr,-epsi))
 221
 222       END FUNCTION m_complex_water_ray
 223
 224 !+---+-----------------------------------------------------------------+
 225
 226       COMPLEX*16 FUNCTION m_complex_ice_maetzler(lambda,T)
 227
 228 !      complex refractive index of ice as function of Temperature T
 229 !      [deg C] and radar wavelength lambda [m]; valid for
 230 !      lambda in [0.0001,30] m; T in [-250.0,0.0] C
 231 !      Original comment from the Matlab-routine of Prof. Maetzler:
 232 !      Function for calculating the relative permittivity of pure ice in
 233 !      the microwave region, according to C. Maetzler, "Microwave
 234 !      properties of ice and snow", in B. Schmitt et al. (eds.) Solar
 235 !      System Ices, Astrophys. and Space Sci. Library, Vol. 227, Kluwer
 236 !      Academic Publishers, Dordrecht, pp. 241-257 (1998). Input:
 237 !      TK = temperature (K), range 20 to 273.15
 238 !      f = frequency in GHz, range 0.01 to 3000
 239
 240       IMPLICIT NONE
 241       DOUBLE PRECISION, INTENT(IN):: T,lambda
 242       DOUBLE PRECISION:: f,c,TK,B1,B2,b,deltabeta,betam,beta,theta,alfa
 243
 244       c = 2.99d8
 245       TK = T + 273.16
 246       f = c / lambda * 1d-9
 247
 248       B1 = 0.0207
 249       B2 = 1.16d-11
 250       b = 335.0d0
 251       deltabeta = EXP(-10.02 + 0.0364*(TK-273.16))
 252       betam = (B1/TK) * ( EXP(b/TK) / ((EXP(b/TK)-1)**2) ) + B2*f*f
 253       beta = betam + deltabeta
 254       theta = 300. / TK - 1.
 255       alfa = (0.00504d0 + 0.0062d0*theta) * EXP(-22.1d0*theta)
 256       m_complex_ice_maetzler = 3.1884 + 9.1e-4*(TK-273.16)
 257       m_complex_ice_maetzler = m_complex_ice_maetzler                   &
 258                              + CMPLX(0.0d0, (alfa/f + beta*f))
 259       m_complex_ice_maetzler = SQRT(CONJG(m_complex_ice_maetzler))
 260
 261       END FUNCTION m_complex_ice_maetzler
 262
 263 !+---+-----------------------------------------------------------------+
 264
 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 !+---+-----------------------------------------------------------------+
 361
 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(radar_debug,*) 'GET_M_MIX_NESTED: bad matrix: ', matrix
 388         CALL wrf_debug(150, radar_debug)
 389         cumulerror = cumulerror + 1
 390        else
 391         vol1 = volice / MAX(volice+volwater,1d-10)
 392         vol2 = 1.0d0 - vol1
 393         mtmp = get_m_mix (m_a, m_i, m_w, 0.0d0, vol1, vol2,             &
 394                          mixingrule, matrix, inclusion, error)
 395         cumulerror = cumulerror + error
 396
 397         if (hostmatrix .eq. 'air') then
 398          get_m_mix_nested = get_m_mix (m_a, mtmp, 2.0*m_a,              &
 399                          volair, (1.0d0-volair), 0.0d0, mixingrule,     &
 400                          hostmatrix, hostinclusion, error)
 401          cumulerror = cumulerror + error
 402         elseif (hostmatrix .eq. 'icewater') then
 403          get_m_mix_nested = get_m_mix (m_a, mtmp, 2.0*m_a,              &
 404                          volair, (1.0d0-volair), 0.0d0, mixingrule,     &
 405                          'ice', hostinclusion, error)
 406          cumulerror = cumulerror + error
 407         else
 408          write(radar_debug,*) 'GET_M_MIX_NESTED: bad hostmatrix: ',        &
 409                            hostmatrix
 410          CALL wrf_debug(150, radar_debug)
 411          cumulerror = cumulerror + 1
 412         endif
 413        endif
 414
 415       elseif (host .eq. 'ice') then
 416
 417        if (matrix .eq. 'ice') then
 418         write(radar_debug,*) 'GET_M_MIX_NESTED: bad matrix: ', matrix
 419         CALL wrf_debug(150, radar_debug)
 420         cumulerror = cumulerror + 1
 421        else
 422         vol1 = volair / MAX(volair+volwater,1d-10)
 423         vol2 = 1.0d0 - vol1
 424         mtmp = get_m_mix (m_a, m_i, m_w, vol1, 0.0d0, vol2,             &
 425                          mixingrule, matrix, inclusion, error)
 426         cumulerror = cumulerror + error
 427
 428         if (hostmatrix .eq. 'ice') then
 429          get_m_mix_nested = get_m_mix (mtmp, m_i, 2.0*m_a,              &
 430                          (1.0d0-volice), volice, 0.0d0, mixingrule,     &
 431                          hostmatrix, hostinclusion, error)
 432          cumulerror = cumulerror + error
 433         elseif (hostmatrix .eq. 'airwater') then
 434          get_m_mix_nested = get_m_mix (mtmp, m_i, 2.0*m_a,              &
 435                          (1.0d0-volice), volice, 0.0d0, mixingrule,     &
 436                          'air', hostinclusion, error)
 437          cumulerror = cumulerror + error
 438         else
 439          write(radar_debug,*) 'GET_M_MIX_NESTED: bad hostmatrix: ',        &
 440                            hostmatrix
 441          CALL wrf_debug(150, radar_debug)
 442          cumulerror = cumulerror + 1
 443         endif
 444        endif
 445
 446       elseif (host .eq. 'water') then
 447
 448        if (matrix .eq. 'water') then
 449         write(radar_debug,*) 'GET_M_MIX_NESTED: bad matrix: ', matrix
 450         CALL wrf_debug(150, radar_debug)
 451         cumulerror = cumulerror + 1
 452        else
 453         vol1 = volair / MAX(volice+volair,1d-10)
 454         vol2 = 1.0d0 - vol1
 455         mtmp = get_m_mix (m_a, m_i, m_w, vol1, vol2, 0.0d0,             &
 456                          mixingrule, matrix, inclusion, error)
 457         cumulerror = cumulerror + error
 458
 459         if (hostmatrix .eq. 'water') then
 460          get_m_mix_nested = get_m_mix (2*m_a, mtmp, m_w,                &
 461                          0.0d0, (1.0d0-volwater), volwater, mixingrule, &
 462                          hostmatrix, hostinclusion, error)
 463          cumulerror = cumulerror + error
 464         elseif (hostmatrix .eq. 'airice') then
 465          get_m_mix_nested = get_m_mix (2*m_a, mtmp, m_w,                &
 466                          0.0d0, (1.0d0-volwater), volwater, mixingrule, &
 467                          'ice', hostinclusion, error)
 468          cumulerror = cumulerror + error
 469         else
 470          write(radar_debug,*) 'GET_M_MIX_NESTED: bad hostmatrix: ',         &
 471                            hostmatrix
 472          CALL wrf_debug(150, radar_debug)
 473          cumulerror = cumulerror + 1
 474         endif
 475        endif
 476
 477       elseif (host .eq. 'none') then
 478
 479        get_m_mix_nested = get_m_mix (m_a, m_i, m_w,                     &
 480                        volair, volice, volwater, mixingrule,            &
 481                        matrix, inclusion, error)
 482        cumulerror = cumulerror + error
 483
 484       else
 485        write(radar_debug,*) 'GET_M_MIX_NESTED: unknown matrix: ', host
 486        CALL wrf_debug(150, radar_debug)
 487        cumulerror = cumulerror + 1
 488       endif
 489
 490       IF (cumulerror .ne. 0) THEN
 491        write(radar_debug,*) 'GET_M_MIX_NESTED: error encountered'
 492        CALL wrf_debug(150, radar_debug)
 493        get_m_mix_nested = CMPLX(1.0d0,0.0d0)
 494       endif
 495
 496       end function get_m_mix_nested
 497
 498 !+---+-----------------------------------------------------------------+
 499
 500       COMPLEX*16 FUNCTION get_m_mix (m_a, m_i, m_w, volair, volice,     &
 501                      volwater, mixingrule, matrix, inclusion, error)
 502
 503       IMPLICIT NONE
 504
 505       DOUBLE PRECISION, INTENT(in):: volice, volair, volwater
 506       COMPLEX*16, INTENT(in):: m_a, m_i, m_w
 507       CHARACTER(len=*), INTENT(in):: mixingrule, matrix, inclusion
 508       INTEGER, INTENT(out):: error
 509
 510       error = 0
 511       get_m_mix = CMPLX(1.0d0,0.0d0)
 512
 513       if (mixingrule .eq. 'maxwellgarnett') then
 514        if (matrix .eq. 'ice') then
 515         get_m_mix = m_complex_maxwellgarnett(volice, volair, volwater,  &
 516                            m_i, m_a, m_w, inclusion, error)
 517        elseif (matrix .eq. 'water') then
 518         get_m_mix = m_complex_maxwellgarnett(volwater, volair, volice,  &
 519                            m_w, m_a, m_i, inclusion, error)
 520        elseif (matrix .eq. 'air') then
 521         get_m_mix = m_complex_maxwellgarnett(volair, volwater, volice,  &
 522                            m_a, m_w, m_i, inclusion, error)
 523        else
 524         write(radar_debug,*) 'GET_M_MIX: unknown matrix: ', matrix
 525         CALL wrf_debug(150, radar_debug)
 526         error = 1
 527        endif
 528
 529       else
 530        write(radar_debug,*) 'GET_M_MIX: unknown mixingrule: ', mixingrule
 531        CALL wrf_debug(150, radar_debug)
 532        error = 2
 533       endif
 534
 535       if (error .ne. 0) then
 536        write(radar_debug,*) 'GET_M_MIX: error encountered'
 537        CALL wrf_debug(150, radar_debug)
 538       endif
 539
 540       END FUNCTION get_m_mix
 541
 542 !+---+-----------------------------------------------------------------+
 543
 544       COMPLEX*16 FUNCTION m_complex_maxwellgarnett(vol1, vol2, vol3,    &
 545                      m1, m2, m3, inclusion, error)
 546
 547       IMPLICIT NONE
 548
 549       COMPLEX*16 :: m1, m2, m3
 550       DOUBLE PRECISION :: vol1, vol2, vol3
 551       CHARACTER(len=*) :: inclusion
 552
 553       COMPLEX*16 :: beta2, beta3, m1t, m2t, m3t
 554       INTEGER, INTENT(out) :: error
 555
 556       error = 0
 557
 558       if (DABS(vol1+vol2+vol3-1.0d0) .gt. 1d-6) then
 559        write(radar_debug,*) 'M_COMPLEX_MAXWELLGARNETT: sum of the ',       &
 560               'partial volume fractions is not 1...ERROR'
 561        CALL wrf_debug(150, radar_debug)
 562        m_complex_maxwellgarnett=CMPLX(-999.99d0,-999.99d0)
 563        error = 1
 564        return
 565       endif
 566
 567       m1t = m1**2
 568       m2t = m2**2
 569       m3t = m3**2
 570
 571       if (inclusion .eq. 'spherical') then
 572        beta2 = 3.0d0*m1t/(m2t+2.0d0*m1t)
 573        beta3 = 3.0d0*m1t/(m3t+2.0d0*m1t)
 574       elseif (inclusion .eq. 'spheroidal') then
 575        beta2 = 2.0d0*m1t/(m2t-m1t) * (m2t/(m2t-m1t)*LOG(m2t/m1t)-1.0d0)
 576        beta3 = 2.0d0*m1t/(m3t-m1t) * (m3t/(m3t-m1t)*LOG(m3t/m1t)-1.0d0)
 577       else
 578        write(radar_debug,*) 'M_COMPLEX_MAXWELLGARNETT: ',                  &
 579                          'unknown inclusion: ', inclusion
 580        CALL wrf_debug(150, radar_debug)
 581        m_complex_maxwellgarnett=DCMPLX(-999.99d0,-999.99d0)
 582        error = 1
 583        return
 584       endif
 585
 586       m_complex_maxwellgarnett = &
 587        SQRT(((1.0d0-vol2-vol3)*m1t + vol2*beta2*m2t + vol3*beta3*m3t) / &
 588        (1.0d0-vol2-vol3+vol2*beta2+vol3*beta3))
 589
 590       END FUNCTION m_complex_maxwellgarnett
 591
 592 !+---+-----------------------------------------------------------------+
 593       REAL FUNCTION GAMMLN(XX)
 594 !     --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 0.
 595       IMPLICIT NONE
 596       REAL, INTENT(IN):: XX
 597       DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0
 598       DOUBLE PRECISION, DIMENSION(6), PARAMETER:: &
 599                COF = (/76.18009172947146D0, -86.50532032941677D0, &
 600                        24.01409824083091D0, -1.231739572450155D0, &
 601                       .1208650973866179D-2, -.5395239384953D-5/)
 602       DOUBLE PRECISION:: SER,TMP,X,Y
 603       INTEGER:: J
 604
 605       X=XX
 606       Y=X
 607       TMP=X+5.5D0
 608       TMP=(X+0.5D0)*LOG(TMP)-TMP
 609       SER=1.000000000190015D0
 610       DO 11 J=1,6
 611         Y=Y+1.D0
 612         SER=SER+COF(J)/Y
 613 11    CONTINUE
 614       GAMMLN=TMP+LOG(STP*SER/X)
 615       END FUNCTION GAMMLN
 616 !  (C) Copr. 1986-92 Numerical Recipes Software 2.02
 617 !+---+-----------------------------------------------------------------+
 618       REAL FUNCTION WGAMMA(y)
 619
 620       IMPLICIT NONE
 621       REAL, INTENT(IN):: y
 622
 623       WGAMMA = EXP(GAMMLN(y))
 624
 625       END FUNCTION WGAMMA
 626
 627 !+---+-----------------------------------------------------------------+
 628 END MODULE module_mp_radar
 629 !+---+-----------------------------------------------------------------+