phys/module_ra_rrtmg_aero_optical_util_cmaq.F

   1 ! Revision History:
   2 !  2016/02/23 David Wong extracted the complex number module and put it in a file
   3 !  2016/05/23 David Wong - replaced rrtmg_aero_optical_util_module with
   4 !                          cmaq_rrtmg_aero_optical_util_module to avoid duplication
   5 !                          of the same module name on WRF side of the two-way model
   6
   7 MODULE rrtmg_aero_optical_util_module
   8
   9   Integer      :: AERO_UTIL_LOG = 0
  10
  11   private
  12   public :: aero_optical, aero_optical2, aero_optical_CS, AERO_UTIL_LOG
  13
  14   interface ghintBH
  15     module procedure ghintBH_1, ghintBH_2, ghintBH_Odd
  16   end interface
  17
  18   interface ghintBH_CS
  19     module procedure ghintBH_CS_even, ghintBH_CS_odd
  20   end interface
  21
  22   Logical, Parameter :: Use_Odd_Quadrature = .True.
  23   Integer, Parameter :: Quadrature_Points = 3
  24
  25 !B.Hutzell One point quadature IGH = 1
  26
  27   real, parameter :: ghxi_1(1) = 0.00000000000
  28   real, parameter :: ghwi_1(1) = 1.77245385091
  29
  30 !B.Hutzell Three point quadature IGH = 3
  31   real, parameter :: ghxi_3(3) = (/ -1.22474487139,   &
  32                                      0.00000000000,   &
  33                                      1.22474487139 /)
  34
  35   real, parameter :: ghwi_3(3) = (/ 0.295408975151,   &
  36                                     1.181635900000,   &
  37                                     0.295408975151 /)
  38
  39 !B.Hutzell Five point quadature IGH = 5
  40   real(8), parameter :: ghxi_5(5) = (/ -2.02018287046d0,  &
  41                                        -0.958572464614d0, &
  42                                         0.00000000000d0,  &
  43                                         0.958572464614d0, &
  44                                         2.02018287046d0 /)
  45
  46   real(8), parameter :: ghwi_5(5) = (/ 0.019953242059d0,   &
  47                                        0.393619323152d0,   &
  48                                        0.945308720483d0,   &
  49                                        0.393619323152d0,   &
  50                                        0.019953242059d0 /)
  51
  52 !B.Hutzell Nine point quadature IGH = 9 points
  53 !No. Abscissas  Weight  Total Weight
  54   real, parameter :: ghxi_9(9) = (/ -3.19099320178,  &
  55                                     -2.26658058453,  &
  56                                     -1.46855328922,  &
  57                                     -0.72355101875,  &
  58                                      0.00000000000,  &
  59                                      0.72355101875,  &
  60                                      1.46855328922,  &
  61                                      2.26658058453,  &
  62                                      3.19099320178 /)
  63
  64   real, parameter :: ghwi_9(9) = (/ 3.96069772633E-5, &
  65                                     0.00494362428,    &
  66                                     0.08847452739,    &
  67                                     0.43265155900,    &
  68                                     0.72023521561,    &
  69                                     0.43265155900,    &
  70                                     0.08847452739,    &
  71                                     0.004943624275,   &
  72                                     3.96069772633E-5 /)
  73
  74   contains
  75
  76 ! ------------------------------------------------------------------
  77   subroutine getqext_BH (xx, crefin, qextalf, qscatalf, gscatalfg,SUCCESS)
  78
  79     implicit none
  80
  81     real, intent(in)     :: XX
  82     real, intent(out)    :: qextalf, qscatalf, gscatalfg
  83     complex, intent(in)  :: CREFIN
  84     logical, intent(out) :: success
  85
  86 ! local variables
  87     real( 8 ), parameter :: one_third = 1.0d0 / 3.0d0
  88
  89     integer              :: NXX
  90     integer              :: nstop, modulus
  91
  92     real :: QEXT, QSCA, QBACK, G_MIE, xx1
  93
  94     real( 8 )    :: x
  95     complex( 8 ) :: refractive_index
  96
  97     x = real( XX, 8 )
  98     refractive_index = dcmplx( real( CREFIN ), imag( CREFIN ) )
  99
 100     modulus = int( abs( x * refractive_index ) )
 101     nstop = int( x + 4.0d0 * x**one_third + 2.0d0 )
 102
 103     nxx = max( modulus, nstop ) + 15
 104
 105     xx1 = 1.0 / XX
 106
 107     CALL BHMIE_FLEXI (XX, NXX, NSTOP, CREFIN,QEXT,QSCA,QBACK,G_MIE, SUCCESS)
 108
 109     qextalf   = QEXT * xx1
 110     qscatalf  = QSCA * xx1
 111     gscatalfg = qscatalf * G_MIE
 112
 113   end subroutine getqext_bh
 114
 115 ! ------------------------------------------------------------------
 116   SUBROUTINE BHMIE (X, REFREL, QQEXT, QQSCA, QBACK, GSCA, SUCCESS)
 117
 118 ! FSB Changed the call vector to return only QEXT, QSCAT QBACK GSCA
 119 !     and ignore NANG, S1 and S2 and all calculations for them
 120
 121     implicit none
 122
 123 ! Arguments:
 124        real, intent(in)    :: X        ! X = pi*particle_diameter / Wavelength
 125        complex, intent(in) :: REFREL
 126
 127 !    REFREL = (complex refr. index of sphere)/(real index of medium)
 128 !    in the current use the index of refraction of the the medium
 129 !    i taken at 1.0 real.
 130 !
 131 !    Output
 132
 133     real,    intent(out) :: QQEXT, QQSCA, QBACK, GSCA
 134     logical, intent(out) :: SUCCESS
 135
 136 !     QQEXT   Efficiency factor for extinction
 137 !     QQSCA   Efficiency factor for scattering
 138 !     QQBACK  Efficiency factor for back scatter
 139 !     GSCA    asymmetry factor <cos>
 140 !     SUCCESS flag for successful calculation
 141 ! REFERENCE:
 142 !  Bohren, Craig F. and Donald R. Huffman, Absorption and
 143 !    Scattering of Light by Small Particles, Wiley-Interscience
 144 !    copyright 1983. Paperback Published 1998.
 145 ! FSB
 146 !    This code was originally listed in Appendix A. pp 477-482.
 147 !    As noted below, the original code was subsequently
 148 !    modified by Prof. Bruce T. Drain of Princetion University.
 149 !    The code was further modified for a specific application
 150 !    in a large three-dimensional code requiring as much
 151 !    computational efficiency as possible.
 152 !    Prof. Francis S. Binkowski of The University of North
 153 !    Carolina at Chapel Hill.
 154
 155 ! Declare parameters:
 156 ! Note: important that MXNANG be consistent with dimension of S1 and S2
 157 !       in calling routine!
 158
 159     integer, parameter :: MXNANG=10, NMXX=600000   ! FSB new limits
 160     real*8, parameter  :: PII = 3.1415916536D0
 161     real*8, parameter  :: ONE = 1.0D0, TWO = 2.0D0
 162
 163 ! Local variables:
 164     integer    :: NANG
 165     integer    :: N,NSTOP,NMX,NN
 166     real*8     :: QSCA, QEXT, DX1, DXX1
 167     real*8     :: CHI,CHI0,CHI1,DX,EN,P,PSI,PSI0,PSI1,XSTOP,YMOD
 168     real*8     :: TWO_N_M_ONE, TWO_N_P_ONE, EN1, FACTOR
 169     complex*16 :: AN,AN1,BN,BN1,DREFRL,XI,XI1,Y, Y1, DREFRL1
 170     complex*16 :: D(NMXX), FAC1, FAC2
 171     complex*16 :: XBACK
 172
 173 !***********************************************************************
 174 ! Subroutine BHMIE is the Bohren-Huffman Mie scattering subroutine
 175 !    to calculate scattering and absorption by a homogenous isotropic
 176 !    sphere.
 177 ! Given:
 178 !    X = 2*pi*a/lambda
 179 !    REFREL = (complex refr. index of sphere)/(real index of medium)
 180 !    real refractive index of medium taken as 1.0
 181 ! Returns:
 182 !    QEXT  = efficiency factor for extinction
 183 !    QSCA  = efficiency factor for scattering
 184 !    QBACK = efficiency factor for backscatter
 185 !            see Bohren & Huffman 1983 p. 122
 186 !    GSCA = <cos> asymmetry for scattering
 187 !
 188 ! Original program taken from Bohren and Huffman (1983), Appendix A
 189 ! Modified by Prof. Bruce T.Draine, Princeton Univ. Obs., 90/10/26
 190 ! in order to compute <cos(theta)>
 191 ! 91/05/07 (BTD): Modified to allow NANG=1
 192 ! 91/08/15 (BTD): Corrected error (failure to initialize P)
 193 ! 91/08/15 (BTD): Modified to enhance vectorizability.
 194 ! 91/08/15 (BTD): Modified to make NANG=2 if called with NANG=1
 195 ! 91/08/15 (BTD): Changed definition of QBACK.
 196 ! 92/01/08 (BTD): Converted to full double precision and double complex
 197 !                 eliminated 2 unneed lines of code
 198 !                 eliminated redundant variables (e.g. APSI,APSI0)
 199 !                 renamed RN -> EN = double precision N
 200 !                 Note that DOUBLE COMPLEX and DCMPLX are not part
 201 !                 of f77 standard, so this version may not be fully
 202 !                 portable.  In event that portable version is
 203 !                 needed, use src/bhmie_f77.f
 204 ! 93/06/01 (BTD): Changed AMAX1 to generic function MAX
 205 ! FSB April 09,2012 This code was modified by:
 206 ! Prof.  Francis S. Binkowski University of North Carolina at
 207 ! Chapel Hill, Institue for the Environment.
 208 !
 209 ! The modifications were made to enhance computation speed
 210 ! for use in a three-dimensional code. This was done by
 211 ! removing code that calculated angular scattering. The method
 212 ! of calculating QEXT, QBACK was also changed.
 213
 214 !***********************************************************************
 215 !*** Safety checks
 216
 217     SUCCESS = .TRUE.
 218     NANG = 2 ! FSB only this value
 219 ! IF(NANG.GT.MXNANG)STOP'***Error: NANG > MXNANG in bhmie'
 220 !      IF (NANG .LT. 2) NANG = 2
 221
 222     DX = REAL( X, 8 )
 223 ! FSB Define reciprocals so that divisions can be replaced by multiplications.
 224     DX1  = ONE / DX
 225     DXX1 = DX1 * DX1
 226     DREFRL = DCMPLX( REFREL )
 227     DREFRL1 = ONE / DREFRL
 228     Y = DX * DREFRL
 229     Y1 = ONE / Y
 230     YMOD = ABS(Y)
 231
 232 !*** Series expansion terminated after NSTOP terms
 233 !    Logarithmic derivatives calculated from NMX on down
 234     XSTOP = REAL( X + 4.0 * X**0.3333 + 2.0, 8)
 235     NMX  = INT( MAX(XSTOP,YMOD) ) + 15
 236
 237 ! BTD experiment 91/1/15: add one more term to series and compare results
 238 !      NMX=AMAX1(XSTOP,YMOD)+16
 239 ! test: compute 7001 wavelengths between .0001 and 1000 micron
 240 ! for a=1.0micron SiC grain.  When NMX increased by 1, only a single
 241 ! computed number changed (out of 4*7001) and it only changed by 1/8387
 242 ! conclusion: we are indeed retaining enough terms in series!
 243     NSTOP = INT( XSTOP )
 244     FACTOR = 1.0D0
 245
 246     IF (NMX .GT. NMXX) THEN
 247        WRITE(6,*)'Error: NMX > NMXX=',NMXX,' for |m|x=',YMOD
 248        SUCCESS = .FALSE.
 249        RETURN
 250     END IF
 251
 252 ! FSB all code relating to scattering angles is removed out for
 253 !     reasons of efficiency when running in a three-dimensional
 254 !     code. We only need QQSCA, QQEXT, GSCA AND QBACK
 255
 256
 257 !*** Logarithmic derivative D(J) calculated by downward recurrence
 258 !    beginning with initial value (0.,0.)
 259
 260     D(NMX) = DCMPLX(0.0D0,0.0D0)
 261     NN = NMX - 1
 262     DO N = 1,NN
 263        EN  = REAL(NMX - N + 1, 8 )
 264 ! FSB In the following division by Y has been replaced by
 265 !     multiplication by Y1, the reciprocal of Y.
 266        D(NMX-N) = ( EN * Y1 ) - (ONE / ( D(NMX-N+1) + EN * Y1))
 267     END DO
 268
 269 !*** Riccati-Bessel functions with real argument X
 270 !    calculated by upward recurrence
 271
 272     PSI0 =  COS(DX)
 273     PSI1 =  SIN(DX)
 274     CHI0 = -SIN(DX)
 275     CHI1 =  PSI0
 276     XI1  =  DCMPLX(PSI1,-CHI1)
 277     QSCA =  0.0D0
 278     GSCA =  0.0D0
 279     QEXT =  0.0D0
 280     P    = -ONE
 281     XBACK = (0.0d0,0.0d0)
 282
 283 ! FSB Start main loop
 284     DO N = 1,NSTOP
 285        EN        = REAL( N, 8)
 286        EN1       = ONE / EN
 287        TWO_N_M_ONE = TWO * EN - ONE
 288 ! for given N, PSI  = psi_n        CHI  = chi_n
 289 !              PSI1 = psi_{n-1}    CHI1 = chi_{n-1}
 290 !              PSI0 = psi_{n-2}    CHI0 = chi_{n-2}
 291 ! Calculate psi_n and chi_n
 292        PSI = TWO_N_M_ONE * PSI1 * DX1 - PSI0
 293        CHI = TWO_N_M_ONE * CHI1 * DX1 - CHI0
 294        XI  = DCMPLX(PSI,-CHI)
 295
 296 !*** Compute AN and BN:
 297 ! FSB Rearrange to get common terms
 298        FAC1 = D(N) * DREFRL1 + EN * DX1
 299        AN   = (FAC1) * PSI - PSI1
 300        AN   = AN / ( (FAC1 )* XI - XI1 )
 301        FAC2 = ( DREFRL * D(N) + EN * DX1)
 302        BN   = ( FAC2) * PSI -PSI1
 303        BN   = BN / ((FAC2) * XI - XI1 )
 304
 305 ! FSB calculate sum for QEXT as done by Wiscombe
 306 !     get common factor
 307        TWO_N_P_ONE = (TWO * EN + ONE)
 308        QEXT = QEXT + (TWO_N_P_ONE) * (REAL(AN) + REAL(BN) )
 309        QSCA = QSCA + (TWO_N_P_ONE) * ( ABS(AN)**2+ ABS(BN)**2 )
 310
 311 ! FSB calculate XBACK from B & H Page 122
 312        FACTOR = -1.0d0 * FACTOR  ! calculate (-1.0 ** N)
 313        XBACK = XBACK + (TWO_N_P_ONE) * factor * (AN - BN)
 314
 315 ! FSB calculate asymmetry factor
 316         GSCA = GSCA + REAL( ((TWO_N_P_ONE)/(EN * (EN + ONE))) *     &
 317               (REAL(AN)*REAL(BN)+IMAG(AN)*IMAG(BN)))
 318
 319        IF (N .GT. 1)THEN
 320           GSCA = GSCA + REAL( (EN - EN1) *                         &
 321                  (REAL(AN1)*REAL(AN) + IMAG(AN1)*IMAG(AN) +  &
 322                   REAL(BN1)*REAL(BN) + IMAG(BN1)*IMAG(BN)))
 323        ENDIF
 324
 325 !*** Store previous values of AN and BN for use in computation of g=<cos(theta)>
 326        AN1 = AN
 327        BN1 = BN
 328
 329 ! FSB set up for next iteration
 330        PSI0 = PSI1
 331        PSI1 = PSI
 332        CHI0 = CHI1
 333        CHI1 = CHI
 334        XI1  = DCMPLX(PSI1,-CHI1)
 335
 336     END DO   ! main  loop on n
 337
 338 !*** Have summed sufficient terms.
 339
 340 !    Now compute QQSCA,QQEXT,QBACK,and GSCA
 341     GSCA  = REAL( TWO / QSCA )  * GSCA
 342
 343 ! FSB in the following, divisions by DX * DX has been replaced by
 344 !      multiplication by DXX1 the reciprocal of 1.0 / (DX *DX)
 345     QQSCA = REAL( TWO * QSCA * DXX1 )
 346     QQEXT = REAL( TWO * QEXT * DXX1 )
 347     QBACK = REAL( REAL ( 0.5d0 * XBACK * CONJG(XBACK), 8 ) * DXX1 )  ! B&H Page 122
 348
 349   END subroutine BHMIE
 350
 351 ! ------------------------------------------------------------------
 352   subroutine aero_optical ( lamda_in, nmode, nr, ni, Vol,   &
 353                             dgn, sig, bext, bscat, g_bar,   &
 354                             success, modulus )
 355
 356 ! *** calculate the extinction and scattering coefficients and
 357 !     assymetry factors for each wavelength as a sum over the
 358 !     individual lognormal modes. Each mode may have a different
 359 !     set of refractive indices.
 360
 361    IMPLICIT NONE
 362 ! *** input variables
 363    real, intent(in)    :: lamda_in               ! wavelengths  [micro-m]
 364    INTEGER, intent(in) :: nmode                  ! number of lognormal modes
 365    real, intent(in)    :: nr( nmode), ni(nmode)  ! real and imaginary
 366                                                  ! refractive indices
 367    real, intent(in)    :: Vol(nmode)             ! modal aerosol volumes [m**3 /m**3]
 368    real, intent(in)    :: dgn(nmode)             ! geometric mean diameters
 369                                                  ! for number distribution [ m]
 370    real, intent(in)    :: sig(nmode)             ! geometric standard deviation
 371
 372    real, intent(in), optional :: modulus(nmode)  ! modulus of refracive index
 373
 374 ! *** output variables
 375    real, intent(out) :: bext    ! extinction coefficient [ 1 / m ]
 376    real, intent(out) :: bscat   ! scattering coefficient [ 1 / m ]
 377    real, intent(out) :: g_bar   ! assymetry factor for Mie and molecular scattering
 378    logical, intent(out) :: success ! flag for successful calculation
 379 ! *** internal variables
 380    INTEGER  :: j             ! loop index
 381 !     real     :: xlnsig(nmode) ! natural log of geometric standard deviations
 382    real     :: beta_Sc, bsc  !aerosol scattering coefficient
 383
 384    real     :: beta_Ex       ! aerosol extinction coefficients
 385    real     :: G             ! modal aerosol assymetry factors
 386    real     :: sum_g
 387    real     :: LSIGX
 388    real     :: lamdam1       ! 1/ lamda
 389    real     :: alphav        ! Mie size parameter
 390    real     :: vfac
 391    real     :: modalph
 392
 393    real, parameter :: pi = 3.14159265359
 394
 395    Logical, Save :: Initialize = .True.
 396
 397 ! *** coded 09/08/2004 by Dr. Francis S. Binkowski
 398 ! FSB Modified for RRTMG version December 2009.
 399 ! FSB modified 10/06/2004, 10/12/2004, 10/18/2005
 400 ! FSB 01/12/2006
 401 !     Formerly Carolina Environmental Program
 402 ! FSB now the Institute for the Environment
 403 !     University of North Carolina at Chapel Hill
 404 !     email: frank_binkowski@unc.edu
 405
 406
 407 ! *** initialize variables
 408     lamdam1 = 1.0e6 / lamda_in   ! lamda now in [ m ]
 409     bext    = 0.0
 410     bscat   = 0.0
 411     sum_g   = 0.0
 412
 413     DO j = 1, nmode
 414 !    calculate the extinction and scattering coefficients
 415 !    for each mode
 416        LSIGX = log(sig(j))
 417
 418 !     calculate Mie size parameter for volume distribution
 419 !     exp(3.0 * xlnsig*xlnsig)  converts dgn to dgv (volume diameter)
 420        alphav =  pi * dgn(j) * exp(3.0 * LSIGX * LSIGX) * lamdam1
 421
 422        if (present(modulus)) then
 423           modalph = alphav * modulus(j)
 424        end if
 425
 426        CALL ghintBH (nr(j), ni(j), alphav, LSIGX, beta_EX, beta_Sc, G, success)
 427
 428 ! *** ghintBH returns the normalized values
 429 !     Calculate the actual extinction and scattering coefficients
 430 !     by multplying by the modal volume and dividing by the wavelength
 431
 432        vfac  =  Vol(j) * lamdam1
 433
 434 ! *** sum to get total extinction and scattering
 435 !     and contribution to the overal assymetry factor
 436
 437        bext    = bext  + vfac * beta_Ex  ! [ 1 / m ]
 438        bsc     = vfac * beta_Sc
 439        bscat   = bscat + bsc
 440        sum_g   = sum_g + bsc * G
 441
 442     END DO  ! loop on modes
 443
 444 ! *** calculate combined assymetry factor for all modes
 445
 446     g_bar = sum_g / bscat ! changed to divide by bscat
 447
 448   END SUBROUTINE aero_optical
 449
 450 ! ------------------------------------------------------------------
 451   subroutine ghintBH_1 (nr, ni, alfv, xlnsig, Qext_GH, Qscat_GH, g_gh, success)
 452
 453 ! FSB *********** This is the newest (05_30_2012) version of GhintBH
 454 !      this version does the Mie method and calculates the optimum set of
 455 !      set of Gauss-Hermite abscissas and weights.
 456 ! FSB Calls Penndorf codes for alfv .le.  0.3
 457
 458 !      Dr. Francis S. Binkowski, The University of North Carolina
 459 !                                at Chapel Hill
 460 ! FSB this code file now contains all of the necessary subroutines that
 461 !     are called to perform an integral of the Bohren and Huffman
 462 !     Mie codes ( as updated by Prof. Bruce C. Drain of Princeton)
 463 !       calculates the extinction and scattering coefficients
 464 !       normalized by wavelength and total particle volume
 465 !       concentration for a log normal particle distribution
 466 !       with the logarithm of the geometric  standard deviation
 467 !       given by xlnsig. The integral of the
 468 !       asymmetry factor g is also calculated.
 469 ! FSB Change 12/20/2011 This code now has a choice of IGH based
 470 !     upon alfv and nr.
 471 !  *** Does Gauss-Hermite quadrature of Qext / alfa & Qscat / alfa
 472 !      and asymmetry factor <cos> over log normal distribution using
 473 !      symmetric  points.
 474
 475     implicit none
 476
 477     real, intent(in)     :: nr, ni     ! refractive indices
 478     real, intent(in)     :: alfv       ! Mie parameter for dgv
 479     real, intent(in)     :: xlnsig     ! log of geometric  standard deviation
 480     real, intent(out)    :: Qext_GH    ! normalized extinction efficiency
 481     real, intent(out)    :: Qscat_GH   ! normalized scattering efficiency
 482     real, intent(out)    :: g_GH       ! asymmetry factor <cos>
 483     logical, intent(out) :: success    ! flag for successful calculation
 484
 485     real    :: bext_P, bscat_P, babs_P, g_PCS, xlnsg2  ! see below for definition
 486
 487     real    :: aa1                ! see below for definition
 488     real    :: alfaip, alfaim     ! Mie parameters at abscissas
 489
 490 !  *** these are Qext/alfa and Qscat/alfv at the abscissas
 491     real    :: qalfip_e, qalfim_e ! extinction
 492     real    :: qalfip_s, qalfim_s ! scattering
 493     real    :: gsalfp, gsalfm     ! scattering times asymmetry factor
 494     integer :: IGH                ! index for GH quadrature
 495
 496 ! FSB define parameters
 497     real, parameter :: pi = 3.14159265
 498     real, parameter :: sqrtpi = 1.772454
 499     real, parameter :: sqrtpi1 = 1.0 / sqrtpi
 500     real, parameter :: sqrt2 = 1.414214
 501     real, parameter :: three_pi_two = 3.0 * pi / 2.0
 502     real, parameter :: const = three_pi_two * sqrtpi1
 503
 504     integer :: i
 505     complex :: crefin                  ! complex index of refraction
 506     real    :: sum_e,sum_s, xi,wxi,xf
 507     real    :: sum_sg
 508
 509 ! Gauss-Hermite abscissas and weights
 510 ! *** the following weights and abscissas are from Abramowitz
 511 !     Stegun, Table 25.10 page 924
 512 ! FSB full precision from Table 25.10
 513
 514 ! FSB ten-point  - IGH = 5
 515     real, parameter :: ghxi_10(5) = (/ 0.342901327223705,     &
 516                                        1.036610829789514,     &
 517                                        1.756683649299882,     &
 518                                        2.532731674232790,     &
 519                                        3.436159118837738 /)
 520
 521     real, parameter :: ghwi_10(5) = (/ 6.108626337353e-01,    &
 522                                        2.401386110823e-01,    &
 523                                        3.387439445548e-02,    &
 524                                        1.343645746781e-03,    &
 525                                        7.640432855233e-06 /)
 526
 527 ! FSB six-point - IGH = 3
 528     real, parameter :: ghxi_6(3) = (/ 0.436077411927617,      &
 529                                       1.335849074013597,      &
 530                                       2.350604973674492 /)
 531
 532     real, parameter :: ghwi_6(3) = (/ 7.246295952244e-01,     &
 533                                       1.570673203229e-01,     &
 534                                       4.530009905509e-03 /)
 535
 536 ! FSB two-point - IGH = 1
 537     real, parameter :: ghxi_2(1) = (/ 0.707106781186548 /)
 538
 539     real, parameter :: ghwi_2(1) = (/ 8.862269254528e-01 /)
 540
 541     real    :: GHXI(5), GHWI(5) ! weight and abscissas
 542     integer :: NMAX             ! number of weights and abscissa
 543
 544 ! FSB Check for valid range of Penndorf application.
 545     if ( alfv .le. 0.3) then
 546        xlnsg2 = xlnsig*xlnsig
 547        call pennfsb (nr,ni,alfv,xlnsg2,bext_P,bscat_P,babs_P,g_PCS)
 548        Qext_GH  = bext_P
 549        Qscat_GH = bscat_p
 550        g_GH     = g_PCS * exp(4.0 * xlnsg2) ! match GH integral
 551     else
 552
 553 ! FSB We need to do a full Mie calculation now
 554 !     Choose IGH. These choices are designed to improve
 555 !     the computational efficiency without sacrificing accuracy.
 556
 557        IGH=3 ! default value; six_point is sufficient generally
 558 ! six point
 559        NMAX = 3
 560
 561        if (nr .ge. 1.7) then
 562 ! 10 point
 563           IGH = 5 ! more points needed here
 564           NMAX = 5
 565        end if
 566
 567        if ( alfv .gt. 20.0 .or. alfv .lt. 0.5 ) then
 568           IGH  = 1 ! in  this range fewer points are needed
 569           NMAX = 1
 570        end if
 571
 572        if (IGH == 1) then
 573           GHXI(1)    = ghxi_2(1)
 574           GHWI(1)    = ghwi_2(1)
 575        else if (IGH == 3) then
 576           do i = 1, NMAX
 577              GHXI(i) = ghxi_6(i)
 578              GHWI(i) = ghwi_6(i)
 579           end do
 580        else
 581           do i = 1,NMAX
 582              GHXI(i) = ghxi_10(i)
 583              GHWI(i) = ghwi_10(i)
 584           end do
 585        end if ! set up number of abscissas and weights
 586
 587 ! FSB set  complex refractive index.
 588        crefin= cmplx(nr,ni)
 589
 590 ! FSB now start the integration code
 591        aa1 = sqrt2 * xlnsig   ! This 1.0 / Sqrt( A ) in derivation of the integral
 592                                  ! where A = 1.0 / ( 2.0 * xlnsg**2 )
 593
 594 ! Then alpha = alfv * exp[ u / sqrt(A) ]
 595 ! For Gauss-Hermite Quadrature u = xi
 596 ! Therefore, xf = exp( xi / sqrt(A) ),
 597 !  or xf = exp( xi * aa1 )
 598        sum_e  = 0.0
 599        sum_s  = 0.0
 600        sum_sg = 0.0
 601 ! FSB do NMAX calls to the MIE codes
 602        do i = 1,NMAX
 603           xi      = GHXI(i)
 604           wxi     = GHWI(i)
 605           xf      = exp( xi * aa1 )
 606           alfaip  = alfv * xf
 607           alfaim  = alfv / xf ! division cheaper than another exp()
 608 ! *** call subroutine to fetch the effficiencies
 609
 610           call getqext_BH (alfaip, crefin, qalfip_e, qalfip_s, gsalfp, success)
 611           call getqext_BH (alfaim, crefin, qalfim_e, qalfim_s, gsalfm, success)
 612
 613           sum_e  = sum_e + wxi  * ( qalfip_e + qalfim_e )
 614           sum_s  = sum_s + wxi  * ( qalfip_s + qalfim_s )
 615           sum_sg = sum_sg + wxi * ( gsalfp + gsalfm )
 616        end do
 617
 618        g_GH     = sum_sg / sum_s ! this is <cos>
 619        Qext_GH  = const * sum_e  !
 620        Qscat_GH = const * sum_s
 621     end if
 622
 623   end subroutine ghintBH_1
 624
 625 ! ------------------------------------------------------------------
 626   subroutine pennfsb (n, k, xx, lnsg2, bext, bscat, babs, g)
 627
 628 ! FSB a new version of Penndorf's equations. This version does
 629 !     analytical integration for Qext, Qscat, Qabs to generate
 630 !     bext, bscat, babs. Note that the expressions for Qext & Qscat
 631 !     hve been divide through by xx.
 632 !
 633 !     Reference:
 634 !     Caldas, M., V. Semiao, 2001, Radiative properties of small
 635 !                     particles: and extension of the Penndorff Model. Journal
 636 !                     of the Optical Society of America A, Vol. 18, No. 4,
 637 !                     pp 831-838.
 638
 639 !       Penndorf, R., 1962a,Scattering and extinction coefficients for small
 640 !                     absorbing and nonabsorbing aerosols,
 641 !                     J. Optical Society of America, 52, 896-904.
 642
 643 !       Penndorf, P., 1962b,Scattering and extinction coefficients for
 644 !                     small Spherical aerosols, J. Atmos. Sci., 19, p 193
 645
 646 ! FSB Coded by Dr. Francis S. Binkowski on October 25, 2011 by combining
 647 !     two previous versions to get a common code for the Penndorf and
 648 !     and Caldas & Semiao approaches. The Penndorf Qext, Qscat are much
 649 !     better than the versions from Caldas & Semiao despite claims to
 650 !     the contrary. The values of the asymmetry factor from Caldas & Semiao
 651 !     are better than can be obtained from Penndorf.
 652
 653 ! FSB This version does the analytical integral ove a lognormal
 654 !     size distribution.
 655
 656     implicit none
 657 !     input variables
 658     real, intent(in)  :: n, k     ! refractive index
 659     real, intent(in)  :: xx       ! pi * diameter / wavelength
 660     real, intent(in)  :: lnsg2    ! log(sigma_g)**2
 661     real, intent(out) :: bext     ! extinction coefficient
 662     real, intent(out) :: bscat    ! scattering coefficient
 663     real, intent(out) :: babs     ! absorption coefficient
 664     real, intent(out) :: g        ! asmmetry factor
 665
 666 !     internal variables
 667     complex*16  :: m, m2,m4,m6,m21,m22
 668     complex*16  :: P,Q,R,S,T,U,V,W
 669     complex*16  :: Qprime, Rprime,Sprime,Tprime
 670     complex*16  :: Uprime, Vprime, Wprime
 671     real*8      :: Qs, gQs, gpennCS
 672     real*8      :: P1,P2, Q1, Q2 , S2,V1, V2 ! see usage
 673     real*8      :: P1SQ, P2SQ  ! see usage
 674     real*8      :: y, y2, y3, y4, y6, y7,  y8, y9
 675     real*8      :: x, x2, x3, x4, x6, x7,  x8, x9
 676     real        :: mag, modalf
 677 ! FSB define useful numbers and fractions
 678     real, parameter :: pi = 3.14159265358979324d0
 679     real, parameter :: three_pi_two = 1.5d0 * pi
 680
 681     real*8, parameter :: one = 1.0d0
 682     real*8, parameter :: two = 2.0d0
 683     real*8, parameter :: three = 3.0d0
 684     real*8, parameter :: four = 4.0d0
 685     real*8, parameter :: five = 5.0d0
 686     real*8, parameter :: six = 6.0d0
 687     real*8, parameter :: eight = 8.0d0
 688     real*8, parameter :: nine = 9.0d0
 689     real*8, parameter :: fifteen = 15.0d0
 690     real*8, parameter :: fortyfive = 45.0d0
 691 !   real*8, parameter :: two5ths = two / five
 692     real*8, parameter :: twothrds = two / three
 693     real*8, parameter :: fourthirds = four / three
 694     real*8, parameter :: onefifteenth = one / fifteen
 695     real*8, parameter :: twofifteenths = two * onefifteenth
 696 !   real*8, parameter :: fourninths = four / nine
 697     real*8, parameter :: eightthirds = two * fourthirds
 698     real*8, parameter :: one_big = one / 31500.0d0
 699     real*8, parameter :: two_fortyfive = two / fortyfive
 700     real*8, parameter :: four_225 = four / 225.0d0
 701     real*8, parameter :: one_210 = one / 210.0d0
 702 !   real*8, parameter :: one_half = one / two
 703 !   real*8, parameter :: four_two = two
 704     real*8, parameter :: nine_two = 4.5d0
 705 !   real*8, parameter :: sixteen_two = eight
 706 !   real*8, parameter :: thirtysix_two = 36.0 / two
 707 !   real*8, parameter :: twentyfive_two = 25.0d0 / two
 708 !   real*8, parameter :: sixtyfour_two = 64.0d0 / two
 709 !   real*8, parameter :: fortynine_two = 49.0d0 / two
 710 !   real*8, parameter :: eightyone_two = 81.0d0 / two
 711     real*8            :: A,B,C,D,E, AA,BB,CC
 712
 713 ! FSB start code
 714     mag = sqrt( n * n + k * k )
 715     modalf = mag * xx
 716     y  = REAL( xx, 8 ) ! convert to real*8
 717 ! FSB get powers of y
 718     y2 = y * y
 719     y3 = y2 * y
 720     y4 = y3 * y
 721     y6 = y3 * y3
 722     y7 = y3 * y4
 723     y8 = y4 * y4
 724     y9 = y6 * y3
 725
 726 ! FSB Calculate integrals ove the lognormal distribution
 727 !     this is done term by term and the form is
 728 !     xn = yn * exp( (n**2) * lnsig2 /2.0d0)
 729
 730     x  = y
 731     x2 = y2 * exp( two              * lnsg2)
 732     x3 = y3 * exp( nine_two         * lnsg2)
 733     x4 = y4 ! * exp( eight            * lnsg2)
 734     x6 = y6 ! * exp( thirtysix_two    * lnsg2)
 735     x7 = y7 ! * exp( fortynine_two    * lnsg2)
 736     x8 = y8 ! * exp( fortynine_two    * lnsg2)
 737     x9 = y9 ! * exp( eightyone_two    * lnsg2)
 738
 739
 740 ! FSB explicitly calculate complex refrative index m
 741     m = dcmplx(n,-k)
 742 ! FSB get powers and functions of m
 743     m2 = m * m
 744     m4 = m2 * m2
 745     m6 = m2 * m4
 746     m21 = m2 - one
 747     m22 = m2 + two
 748
 749 ! FSB calculate Penndorf's definitions from Table II of Penndorf (1962a)
 750     P = m21 / m22
 751     Q = (m2 - two ) / m22
 752     S = m21 / ( two * m2 + three)
 753     V = m21
 754 ! FSB get real & imaginary parts following Penndorf's mptation
 755     P1 = real(P)
 756     P2 = -aimag(P)
 757     P1SQ = P1 * P1
 758     P2SQ = P2 * P2
 759
 760     Q1 = real(Q)
 761     Q2 = -aimag(Q)
 762     S2 = -aimag(S)
 763     V1 = real(V)
 764     v2 = -aimag(V)
 765
 766
 767 ! FSB Get bext from Penndorf (1962a) Equation (7) up to x4
 768 !     consistent with equation (8)
 769 !     We have then divided through by x and integrated analytically
 770     bext = REAL( four * P2 + ( 2.4d0 * (P1 * Q2 + P2 * Q1 ) +  twothrds * S2          &
 771          + twofifteenths * V2 ) * x2 + ( eightthirds * ( P1SQ - P2SQ ) ) * x3, 4 )
 772
 773 ! FSB get bscat from Penndorf Equation (9) up to x4
 774 !     we have divided through by x and integrated analytically
 775     bscat = REAL( eightthirds * ( P1SQ + P2SQ ) * x3 )
 776 ! FSB calculate babs
 777 !   babs = bext - bscat
 778
 779 ! FSB now get asymmetry factor from Caldas & Semiao (2001)
 780 !
 781 ! *** The following additional variables from Caldas & Semiao (2001)
 782 !     are defined in Equations 10a to 10h.
 783
 784     R = (m6 + 20.0d0*m4 -200.0d0*m2 + 200.0d0) / m22**2
 785     T = m21 / ( ( 2.0d0 * M2 + 3.0d0) **2 )
 786     U = m21 / (3.0d0 * M2 + 4.0d0 )
 787     W = m21 * ( 2.0d0 * m2 - 5.0d0)
 788
 789 ! *** further definitions from Caldas & Semiao (2001)
 790     Qprime = Q
 791     Rprime = 18.0d0 * R
 792     Sprime = 5.0d0 * S / P
 793     Tprime = 375.0d0 * T / P
 794 !   Uprime = 28.0d0 * U / P
 795     Vprime = V / P
 796     Wprime = 5.0d0 * W / P
 797
 798 ! FSB calculate gQs and Qs from Caldas & Semiao (2001)
 799 ! *** calculate Qs equation 13
 800 !      Qs = eightthirds * abs(P)**2                                        &
 801 !           * (x4 + onefifteenth * real(Qprime) * x6                       &
 802 !           + fourthirds * aimag(P) * x7                                   &
 803 !           + one_big * ( 35.0d0 * abs(Qprime)**2                          &
 804 !           + 20.0d0 * real(Rprime) + 35.0d0 * abs(Vprime)**2              &
 805 !           + 21.0d0 * abs(Sprime)**2 ) * x8                               &
 806 !           + two_fortyfive * aimag( Qprime * ( P - conjg(P) )) * x9 )
 807
 808 ! *** calculate gQs equation 15
 809
 810 !      gQs = four_225 * abs(P)**2 * (                                      &
 811 !              (5.0d0 * Real(Vprime) + 3.0d0 * real(Sprime) ) * x6         &
 812 !             + one_210 * ( 35.0d0 * real(Vprime*conjg(Qprime) )           &
 813 !             + 21.0d0 * real(Sprime * conjg(Qprime) )                     &
 814 !             + 10.0d0 * real(Wprime)- 6.0d0 * real(Tprime) ) * x8         &
 815 !             - twothrds * ( 5.0d0 * aimag(Vprime * conjg(P) )             &
 816 !             + 3.0d0 * aimag(Sprime * conjg(P) ) ) * x9    )
 817
 818 ! FSB recast into specific terms
 819     A = 1.0D0 * x4
 820     B = onefifteenth * real(Qprime) * x6
 821     C = fourthirds * aimag(P) * x7
 822     D = one_big * ( 35.0d0 * abs(Qprime)**2                             &
 823         + 20.0d0 * real(Rprime) + 35.0d0 * abs(Vprime)**2               &
 824         + 21.0d0 * abs(Sprime)**2 ) * x8
 825     E = two_fortyfive * aimag( Qprime * ( P - conjg(P) )) * x9
 826
 827     Qs = eightthirds * abs(P)**2 *( A + B + C + D + E )
 828
 829     AA = (5.0d0 * Real(Vprime) + 3.0d0 * real(Sprime) ) * x6
 830     BB = one_210 * ( 35.0d0 * real(Vprime*conjg(Qprime) )               &
 831          + 21.0d0 * real(Sprime * conjg(Qprime) )                       &
 832          + 10.0d0 * real(Wprime)- 6.0d0 * real(Tprime) ) * x8
 833     CC = twothrds * ( 5.0d0 * aimag(Vprime * conjg(P) )                 &
 834          + 3.0d0 * aimag(Sprime * conjg(P) ) ) * x9
 835
 836     gQs = four_225 * abs(P)**2 * ( AA + BB + CC )
 837
 838 ! FSB calculate asymmetry factor and adjust with empirical term.
 839     g = REAL(gQs / Qs)
 840 !  FSB now multiply by three_pi_two  get output  values
 841     bext  = three_pi_two * bext
 842     bscat = three_pi_two * bscat
 843 ! FSB calculate babs
 844     babs = bext - bscat
 845
 846   end subroutine pennfsb
 847
 848 ! ------------------------------------------------------------------
 849
 850 ! FSB a new version of Penndorf's equations. This version does
 851 !     analytical integration for Qext, Qscat, Qabs to generate
 852 !     bext, bscat, babs. Note that the expressions for Qext & Qscat
 853 !     hve been divide through by xx.
 854 !
 855 !     References:
 856 !       Penndorf, R., 1962a,Scattering and extinction coefficients for small
 857 !                     absorbing and nonabsorbing aerosols,
 858 !                     J. Optical Society of America, 52, 896-904.
 859
 860 !       Penndorf, P., 1962b,Scattering and extinction coefficients for
 861 !                     small Spherical aerosols, J. Atmos. Sci., 19, p 193
 862
 863 ! FSB Coded by Dr. Francis S. Binkowski on October 25, 2011 by combining
 864 !     two previous versions to get a common code for the Penndorf and
 865 !     and Caldas & Semiao approaches. The Penndorf Qext, Qscat are much
 866 !     better than the versions from Caldas & Semiao despite claims to
 867 !     the contrary. The values of the asymmetry factor from Caldas & Semiao
 868 !     are better than can be obtained from Penndorf.
 869 ! FSB Modified by FSB on 12/02/2013 to remove the Caldas & Semiao parts.
 870 !     This version is set up to calculate LW bext and bscat so that absorption
 871 !     in the LW can be calculated. This will work with RRTMG LW.
 872 ! FSB This version does the analytical integral ove a lognormal
 873 !     size distribution.
 874
 875   subroutine pennfsbLW (crefin, xx, lnsig, bext, bscat)
 876
 877     implicit none
 878 !   input variables
 879     complex, intent(in)  :: crefin  !complex refractive index
 880     real, intent(in)     :: xx ! pi * diameter / wavelength
 881     real, intent(in)     :: lnsig   ! log(sigma_g)
 882     real, intent(out)    :: bext  ! extinction coefficient
 883     real, intent(out)    :: bscat ! scattering coefficient
 884 !   internal variables
 885     real*8      :: Qext, Qscat
 886     real*8      :: lnsg2
 887 !   internal variables
 888     complex*16  :: m, m2,m21,m22
 889     complex*16  :: P,Q,S,V
 890     real*8      :: P1,P2, Q1, Q2 , S2,V1, V2 ! see usage
 891     real*8      :: P1SQ, P2SQ  ! see usage
 892     real*8      :: y, y2, y3
 893     real*8      ::    x2, x3
 894 ! FSB define useful numbers and fractions
 895     real*8, parameter :: one= 1.0d0
 896     real*8, parameter :: two = 2.0d0
 897     real*8, parameter :: three = 3.0d0
 898     real*8, parameter :: four = 4.0d0
 899     real*8, parameter :: eight = 8.0d0
 900     real*8, parameter :: nine = 9.0d0
 901     real*8, parameter :: fifteen = 15.0d0
 902     real*8, parameter :: twothrds = two/three
 903     real*8, parameter :: twofifteenths = two / fifteen
 904     real*8, parameter :: eightthirds = eight / three
 905     real*8, parameter :: nine_two = nine / two
 906
 907 ! FSB define useful numbers and fractions
 908     real, parameter :: pi = four*atan(one)
 909     real, parameter :: three_pi_two = 1.5d0 * pi
 910
 911 ! FSB start code
 912
 913 ! FSB Calculate integrals ove the lognormal distribution
 914 !     this is done term by term; the form is
 915 !     xn = yn * exp( (n**2) * lnsig2 /2.0d0)
 916     lnsg2 = lnsig * lnsig
 917     y  = xx
 918     y2 = y * y
 919     y3 = y * y2
 920     x2 = y2 * exp( two      * lnsg2)
 921     x3 = y3 * exp( nine_two * lnsg2)
 922
 923     m= conjg(crefin) ! Penndorf asuumes k is negative.
 924     m2 = m * m
 925     m21 = m2 - one
 926     m22 = m2 + two
 927 ! FSB calculate Penndorf's definitions from Table II of Penndorf (1962a)
 928       P = m21 / m22
 929       Q = (m2 - two ) / m22
 930       S = m21 / ( two * m2 + three)
 931       V = m21
 932 ! FSB get real & imaginary parts following Penndorf's mptation
 933       P1 = real(P)
 934       P2 = -aimag(P)
 935       P1SQ = P1 * P1
 936       P2SQ = P2 * P2
 937
 938       Q1 = real(Q)
 939       Q2 = -aimag(Q)
 940       S2 = -aimag(S)
 941       V1 = real(V)
 942       v2 = -aimag(V)
 943
 944 ! FSB Get bext from Penndorf (1962a) Equation (7) up to x4
 945 !     consistent with equation (8)
 946 !     We have then divided through by x and integrated analytically
 947       Qext = four * P2                                                  &
 948              + ( 2.4d0 * (P1 * Q2 + P2 * Q1 ) +  twothrds * S2          &
 949                        + twofifteenths * V2   ) * x2                    &
 950              + ( eightthirds * ( P1SQ - P2SQ ) ) * x3
 951
 952 ! FSB get bscat from Penndorf Equation (9) up to x4
 953 !     we have divided through by x and integrated analytically
 954       Qscat = eightthirds * ( P1SQ + P2SQ ) * x3
 955
 956 !  FSB now multiply by three_pi_two  get output  values
 957       bext  = three_pi_two * Qext
 958       bscat = three_pi_two * Qscat
 959
 960   end subroutine pennfsbLW
 961
 962 ! ------------------------------------------------------------------
 963   subroutine aero_optical2( lamda_in, crefin, Vol, dgn, &
 964                             sig, bext, bscat, gfac, success )
 965
 966 ! FSB NOTE: this subroutine calculates for single mode
 967
 968 ! *** calculate the extinction and scattering coefficients and
 969 !     assymetry factors for each wavelength as a sum over the
 970 !     individual lognormal modes. Each mode may have a different
 971 !     set of refractive indices.
 972
 973    IMPLICIT NONE
 974
 975 ! *** input variables
 976    real, intent(in)    :: lamda_in   ! wavelengths  [micro-m]
 977    complex, intent(in) :: crefin     ! Complex refractive index
 978    real, intent(in)    :: Vol        ! modal aerosol volumes [m**3 /m**3]
 979    real, intent(in)    :: dgn        ! geometric mean diameters
 980                                         ! for number distribution [ m]
 981    real, intent(in)    :: sig        ! geometric standard deviation
 982
 983 ! *** output variables
 984    real, intent(out)    :: bext         ! extinction coefficient [ 1 / m ]
 985    real, intent(out)    :: bscat        ! scattering coefficient [ 1 / m ]
 986    real, intent(out)    :: gfac         ! assymetry factor for Mie and molecular scattering
 987    logical, intent(out) :: success      ! flag for successful calculation
 988
 989 ! *** internal variables
 990 !  real :: xlnsig(nmode) ! natural log of geometric standard deviations
 991    real :: beta_Sc       ! aerosol scattering coefficient
 992
 993    real :: beta_Ex       ! aerosol extinction coefficients
 994    real :: G             ! modal aerosol assymetry factors
 995    real :: sum_g
 996    real :: LSIGX
 997    real :: lamdam1       ! 1/ lamda
 998    real :: alphav        ! Mie size parameter
 999    real :: vfac
1000    real, parameter :: pi = 3.14159265359
1001
1002     Logical, Save :: Initialize = .True.
1003
1004 ! FSB coded 04/15/2012 by Dr. Francis S. Binkowski
1005 !     modified from an earlier version
1006 !     Center for Environmental Modeling for PolicyDevelopment
1007 !     Institute for the Environment
1008 !     University of North Carolina at Chapel Hill
1009 !     email: frank_binkowski@unc.edu
1010
1011 ! *** initialize variables
1012     lamdam1 = 1.0e6 / lamda_in ! lamda now in [ m ]
1013     bext  = 0.0
1014     bscat = 0.0
1015     sum_g = 0.0
1016     LSIGX = log(sig)
1017
1018 !     calculate Mie size parameter for volume distribution
1019 !     exp(3.0 * xlnsig*xlnsig)  converts dgn to dgv (volume diameter)
1020     alphav =  pi * dgn * exp(3.0 * LSIGX * LSIGX) * lamdam1
1021
1022     If(Initialize .And. AERO_UTIL_LOG .GT. 0 )Then
1023        If( Use_Odd_Quadrature )then
1024            write(AERO_UTIL_LOG,99501)Quadrature_Points
1025        else
1026            write(AERO_UTIL_LOG,99504)
1027            Initialize = .False.
1028        End If
1029     End If
1030
1031     If( Use_Odd_Quadrature )then
1032         CALL ghintBH (Initialize, crefin, alphav, LSIGX, beta_EX, beta_Sc, G, success)
1033     Else
1034         CALL ghintBH (crefin, alphav, LSIGX, beta_EX, beta_Sc, G, success)
1035     End If
1036
1037 ! *** ghintBH returns the normalized values
1038 !     Calculate the actual extinction and scattering coefficients
1039 !     by multplying by the modal volume and dividing by the wavelength
1040
1041     vfac  =  Vol * lamdam1
1042     bext    = vfac * beta_Ex  ! [ 1 / m ]
1043     bscat   = vfac * beta_Sc  ! [ 1 / m ]
1044     gfac    = G
1045 99501  Format(I2,' Quadrature Points for Volume Averaged Aerosol Optics')
1046 99504  Format('Even Number Quadrature Points for Volume Averaged Aerosol Optics')
1047
1048   END SUBROUTINE aero_optical2
1049
1050 ! ------------------------------------------------------------------
1051   subroutine aero_optical_CS ( lamda_in, refcor,refshell, VOLCOR,   &
1052                                VOLSHELL, DGNCOR, DGNSHELL, SIG,     &
1053                                bext, bscat, gfac, succesS )
1054
1055 ! FSB NOTE: values for one mode are returend
1056 ! *** calculate the extinction and scattering coefficients and
1057 !     assymetry factors for each wavelength as a sum over the
1058 !     individual lognormal modes. Each mode may have a different
1059 !     set of refractive indices.
1060
1061     IMPLICIT NONE
1062 ! *** input variables
1063     real,intent(in)    :: lamda_in   ! wavelengths  [micro-m]
1064     complex,intent(in) :: refcor     ! Complex refractive index -core
1065     complex,intent(in) :: refshell   ! Complex refractive index -shell
1066     real,intent(in)    ::  VOLCOR    ! volume of core
1067     real,intent(in)    ::  VOLSHELL  ! volume of shell
1068     real,intent(in)    ::  DGNCOR    ! geometric mean diameters
1069                                      ! for number distribution [m]
1070     real,intent(in)    ::  DGNSHELL  ! geometric mean diameters
1071                                      ! for number distribution [m]
1072     real,intent(in)    ::  SIG       ! geometric standard deviation
1073
1074 ! *** output variables
1075     real,intent(out)     ::  bext      ! extinction coefficient [ 1 / m ]
1076     real,intent(out)     ::  bscat     ! scattering coefficient [ 1 / m ]
1077     real,intent(out)     ::  gfac      ! assymetry factor
1078     logical, intent(OUT) :: success    ! flag for successful calculation
1079
1080 ! *** internal variables
1081 !   real    :: xlnsig(nmode)    ! natural log of geometric standard deviations
1082     real    :: beta_Sc          ! aerosol scattering coefficient
1083
1084     real    :: beta_Ex          ! aerosol extinction coefficients
1085     real    :: G                ! modal aerosol assymetry factors
1086     real    :: LSIGX
1087     real    :: XX, YY           ! Mie size parameter
1088     real    :: expfac
1089     real    :: lamdam1          ! 1/ lamda
1090     real    :: vfac
1091
1092     Logical, Save :: Initialize = .True.
1093
1094     real, parameter :: pi = 3.14159265359
1095
1096 ! FSB coded 04/15/2012 by Dr. Francis S. Binkowski
1097 !     modified from an earlier version
1098 !     Center for Environmental Modeling for PolicyDevelopment
1099 !     Institute for the Environment
1100 !     University of North Carolina at Chapel Hill
1101 !     email: frank_binkowski@unc.edu
1102
1103
1104 ! *** initialize variables
1105     lamdam1 = 1.0e6 / lamda_in ! lamda now in [ m ]
1106
1107 !    calculate the extinction and scattering coefficients
1108     LSIGX  = log(SIG)
1109     expfac = pi * exp(3.0 * LSIGX * LSIGX) * lamdam1
1110
1111 !     calculate Mie size parameter for volume distribution
1112 !     exp(3.0 * xlnsig*xlnsig)  converts dgn to dgv (volume diameter)
1113     XX =  DGNCOR * expfac
1114     YY =  DGNSHELL * expfac
1115
1116     If(Initialize .And. AERO_UTIL_LOG .GT. 0 )Then
1117        If( Use_Odd_Quadrature )then
1118            write(AERO_UTIL_LOG,99500)Quadrature_Points
1119        else
1120            write(AERO_UTIL_LOG,99502)
1121            Initialize = .False.
1122        End If
1123     End If
1124
1125     If( Use_Odd_Quadrature )then
1126         CALL ghintBH_CS(Initialize,refcor,refshell,XX,YY,LSIGX,beta_EX,beta_Sc,G, success)
1127     Else
1128         CALL ghintBH_CS(refcor,refshell,XX,YY,LSIGX,beta_EX,beta_Sc,G, success)
1129     End If
1130
1131 ! FSB ghintBH_CS returns the normalized values
1132 !     Calculate the actual extinction and scattering coefficients
1133 !     by multplying by the modal volume and dividing by the wavelength.
1134 !     For the coated-sphere (core-shell) calculation use the combined
1135 !     volume
1136
1137     vfac   =  (VOLCOR + VOLSHELL) * lamdam1
1138     bext   = vfac * beta_Ex  ! [ 1 / m ]
1139     bscat  = vfac * beta_Sc  ! [ 1 / m ]
1140     gfac   = G
1141 99500  Format(I2,' Quadrature Points for Core-Shell Aerosol Optics')
1142 99502  Format('Even Number Quadrature Points for Core-Shell Aerosol Optics')
1143
1144   END SUBROUTINE aero_optical_CS
1145
1146 ! ------------------------------------------------------------------
1147
1148   subroutine aero_optical_LW (lamda_in, crefin, Vol,   &
1149                               dgn, sig, bext, bscat )
1150
1151 ! *** calculate the extinction and scattering coefficients and
1152 !     asymmetry factor for LW radiation.
1153 !     The calling code only uses
1154
1155    IMPLICIT NONE
1156 ! *** input variables
1157    real, intent(in)    :: lamda_in      ! wavelengths  [micro-m]
1158
1159    complex, intent(in) ::  crefin  ! complex refractive index
1160    real, intent(in)    :: Vol ! modal aerosol volumes [m**3 /m**3]
1161    real, intent(in)    :: dgn ! geometric mean diameter
1162                               ! for number distribution [ m]
1163    real, intent(in)    :: sig ! geometric standard deviation
1164
1165
1166 ! *** output variables
1167    real, intent(out) :: bext  ! extinction coefficient [ 1 / m ]
1168    real, intent(out) :: bscat ! scattering coefficient [ 1 / m ]
1169
1170
1171 ! *** internal variables
1172    real, parameter :: pi = 3.14159265359
1173
1174    real :: lamda  ! wavelength [ m]
1175    real :: beta_Sc       ! aerosol scattering coefficient
1176
1177    real :: beta_Ex       ! aerosol extinction coefficients
1178    real :: G             ! modal aerosol assymetry factors
1179    real :: VLX, DGX, SIGX, LSIGX
1180    real :: lamdam1 ! 1/ lamda
1181    real :: alphav ! Mie size parameter
1182    real vfac
1183
1184    logical :: success
1185
1186 ! *** coded 09/08/2004 by Dr. Francis S. Binkowski
1187 ! FSB Modified for RRTMG version December 2009.
1188 ! FSB modified 10/06/2004, 10/12/2004, 10/18/2005
1189 ! FSB 01/12/2006
1190 ! FSB 12/02/2013
1191 ! FSB Institute for the Environment
1192 !     University of North Carolina at Chapel Hill
1193 !     email: frank_binkowski@unc.edu
1194
1195 ! FSB Start Code:
1196 ! *** initialize variables
1197 !      lamda = 1.0e-6 * lamda_in ! lamda now in [ m ]
1198     bext  = 0.0
1199     bscat = 0.0
1200
1201 !   lamdam1 = 1.0 / lamda  ! 1 / [m]
1202     lamdam1 = 1.0e6 / lamda_in  ! 1 / [m]
1203
1204     VLX   = Vol
1205     DGX   = dgn
1206     SIGX  = sig
1207     LSIGX = log(SIGX)
1208
1209 !     calculate Mie size parameter for volume distribution
1210 !     exp(3.0 * xlnsig*xlnsig)  converts dgn to dgv (volume diameter)
1211     alphav = pi * DGX * exp(3.0 * LSIGX * LSIGX) * lamdam1
1212     vfac   = VLX * lamdam1
1213
1214   ! FSB Check for valid range of Penndorf application.
1215     if ( alphav .le. 0.3) then
1216        call pennfsbLW(crefin, alphav, LSIGX, beta_EX, beta_Sc)
1217        G = 0.0
1218     else
1219        CALL ghintBH(crefin, alphav, LSIGX, beta_EX, beta_Sc, G, success)
1220     end if
1221
1222     bext  = vfac * beta_Ex  ! [ 1 / m ]
1223     bscat = vfac * beta_Sc
1224
1225   END SUBROUTINE aero_optical_LW
1226
1227 ! ------------------------------------------------------------------
1228   subroutine ghintBH_2 (crefin,alfv,xlnsig,Qext_GH,Qscat_GH,g_gh, success)
1229
1230 ! *************** REVISED VERSION < NOTE
1231 ! FSB *********** This is the newest (04_14_2012) version of GhintBH
1232 !      this version does the Mie method and calculates the optimum set of
1233 !      set of Gauss-Hermite abscissas and weights.
1234 !      Dr. Francis S. Binkowski, The University of North Carolina
1235 !                                at Chapel Hill
1236 ! FSB this code file now contains all of the necessary subroutines that
1237 !     are called to perform an integral of the Bohren and Huffman
1238 !     Mie codes ( as updated by Prof. Bruce C. Drain of Princeton)
1239 !       calculates the extinction and scattering coefficients
1240 !       normalized by wavelength and total particle volume
1241 !       concentration for a log normal particle distribution
1242 !       with the logarithm of the geometric  standard deviation
1243 !       given by xlnsig. The integral of the
1244 !       asymmetry factor g is also calculated.
1245 ! FSB Change 12/20/2011 This code now has a choice of IGH based
1246 !     upon alfv and nr.
1247 ! FBB Changes Simplified code. Eliminated Penndorf code
1248 !  *** Does Gauss-Hermite quadrature of Qext / alfa & Qscat / alfa
1249 !      and asymmetry factor <cos> over log normal distribution using
1250 !      symmetric  points.
1251 !
1252     implicit none
1253
1254     complex, intent(in) :: crefin     ! complex index of refraction
1255     real, intent(in)    :: alfv       ! Mie parameter for dgv
1256     real, intent(in)    :: xlnsig     ! log of geometric  standard deviation
1257     real, intent(out)   :: Qext_GH    ! normalized extinction efficiency
1258     real, intent(out)   :: Qscat_GH   ! normalized scattering efficiency
1259     real, intent(out)   :: g_GH       ! asymmetry factor <cos>
1260     logical, intent(out) :: success   ! flag for successful calculation
1261
1262     real    :: nr                 ! real part of  refractive index
1263     real    :: aa1                ! see below for definition
1264     real    :: alfaip, alfaim     ! Mie parameters at abscissas
1265
1266 !  *** these are Qext/alfa and Qscat/alfv at the abscissas
1267     real    :: qalfip_e, qalfim_e ! extinction
1268     real    :: qalfip_s, qalfim_s ! scattering
1269     real    :: gsalfp, gsalfm     ! scattering times asymmetry factor
1270     integer :: IGH                ! index for GH quadrature
1271
1272 ! FSB define parameters
1273     real, parameter :: pi = 3.14159265
1274     real, parameter :: sqrtpi = 1.772454
1275     real, parameter :: sqrtpi1 = 1.0 / sqrtpi
1276     real, parameter :: sqrt2 = 1.414214
1277     real, parameter :: three_pi_two = 3.0 * pi / 2.0
1278     real, parameter :: const = three_pi_two * sqrtpi1
1279
1280     integer ::  i
1281     real    ::  sum_e,sum_s, xi,wxi,xf
1282     real    ::  sum_sg
1283
1284 ! Gauss-Hermite abscissas and weights
1285 ! *** the following weights and abscissas are from Abramowitz
1286 !     Stegun, Table 25.10 page 924
1287 ! FSB full precision from Table 25.10
1288
1289 ! FSB ten-point  - IGH = 5
1290     real, parameter :: ghxi_10(5) = (/ 0.342901327223705,     &
1291                                        1.036610829789514,     &
1292                                        1.756683649299882,     &
1293                                        2.532731674232790,     &
1294                                        3.436159118837738 /)
1295
1296     real, parameter :: ghwi_10(5) = (/ 6.108626337353e-01,    &
1297                                        2.401386110823e-01,    &
1298                                        3.387439445548e-02,    &
1299                                        1.343645746781e-03,    &
1300                                        7.640432855233e-06 /)
1301
1302 ! FSB six-point - IGH = 3
1303     real, parameter :: ghxi_6(3) = (/ 0.436077411927617,      &
1304                                       1.335849074013597,      &
1305                                       2.350604973674492 /)
1306
1307     real, parameter :: ghwi_6(3) = (/ 7.246295952244e-01,     &
1308                                       1.570673203229e-01,     &
1309                                       4.530009905509e-03 /)
1310
1311 ! FSB two-point - IGH = 1
1312     real, parameter :: ghxi_2(1) = (/ 0.707106781186548 /)
1313
1314     real, parameter :: ghwi_2(1) = (/ 8.862269254528e-01 /)
1315
1316     real    :: GHXI(5), GHWI(5) ! weight and abscissas
1317     integer :: NMAX             ! number of weights and abscissa
1318
1319
1320 ! start code
1321 ! FSB now choose IGH. These choices are designed to improve
1322 !     the computational efficiency without sacrificing accuracy.
1323
1324     nr = real(crefin)
1325
1326     IGH=3 ! default value; six_point is sufficient generally
1327 ! six point
1328     NMAX = 3
1329
1330     if (nr .ge. 1.7) then
1331 ! 10 point
1332        IGH = 5 ! more points needed here
1333        NMAX = 5
1334     end if
1335
1336     if( alfv .gt. 20.0 .or. alfv .lt. 0.5 ) then
1337        IGH  = 1 ! in  this range fewer points are needed
1338        NMAX = 1
1339     end if
1340
1341     if (IGH == 1) then
1342 ! two point
1343        GHXI(1)    = ghxi_2(1)
1344        GHWI(1)    = ghwi_2(1)
1345     else if (IGH == 3) then
1346        do i = 1, NMAX
1347           GHXI(i) = ghxi_6(i)
1348           GHWI(i) = ghwi_6(i)
1349        end do
1350     else
1351        do i = 1,NMAX
1352           GHXI(i) = ghxi_10(i)
1353           GHWI(i) = ghwi_10(i)
1354        end do
1355     end if ! set up number of abscissas and weights
1356
1357 ! FSB now start the integration code
1358     aa1 = sqrt2 * xlnsig    ! This 1.0 / Sqrt( A ) in derivation of the integral
1359                             ! where A = 1.0 / ( 2.0 * xlnsg**2 )
1360
1361 ! Then alpha = alfv * exp[ u / sqrt(A) ]
1362 ! For Gauss-Hermite Quadrature u = xi
1363 ! Therefore, xf = exp( xi / sqrt(A) ),
1364 !  or xf = exp( xi * aa1 )
1365
1366     sum_e  = 0.0
1367     sum_s  = 0.0
1368     sum_sg = 0.0
1369 ! FSB do NMAX calls to the MIE codes
1370     do i = 1,NMAX
1371        xi      = GHXI(i)
1372        wxi     = GHWI(i)
1373        xf      = exp( xi * aa1 )
1374        alfaip  = alfv * xf
1375        alfaim  = alfv / xf ! division cheaper than another exp()
1376 ! *** call subroutine to fetch the effficiencies
1377
1378        call getqext_BH(alfaip,crefin,qalfip_e,qalfip_s, gsalfp, success)
1379        call getqext_BH(alfaim,crefin,qalfim_e,qalfim_s, gsalfm, success)
1380
1381        sum_e  = sum_e + wxi  * ( qalfip_e + qalfim_e )
1382        sum_s  = sum_s + wxi  * ( qalfip_s + qalfim_s )
1383        sum_sg = sum_sg + wxi * ( gsalfp + gsalfm )
1384
1385     end do
1386
1387     g_GH     = sum_sg / sum_s ! this is <cos>
1388     Qext_GH  = const * sum_e  !
1389     Qscat_GH = const * sum_s
1390
1391   end subroutine ghintBH_2
1392
1393 ! ------------------------------------------------------------------
1394   subroutine ghintBH_CS_even (RCORE, RSHELL , XX, YY, xlnsig,  &
1395                               Qext_GH,Qscat_GH, g_gh, success)
1396
1397 ! FSB code for coated-sphere (core-shell) version
1398
1399 ! *************** REVISED VERSION < NOTE
1400 ! FSB *********** This is the newest (04_14_2012) version of ghintBH_CS
1401 !      for the coated-sphere (core-shell) method using BHCOAT
1402 !      this version does the Mie method and calculates the optimum set of
1403 !      set of Gauss-Hermite abscissas and weights.
1404 !      Dr. Francis S. Binkowski, The University of North Carolina
1405 !                                at Chapel Hill
1406
1407 ! FSB this code file now contains all of the necessary subroutines that
1408 !     are called to perform an integral of the Bohren and Huffman
1409 !     Mie codes ( as updated by Prof. Bruce C. Drain of Princeton)
1410 !       calculates the extinction and scattering coefficients
1411 !       normalized by wavelength and total particle volume
1412 !       concentration for a log normal particle distribution
1413 !       with the logarithm of the geometric  standard deviation
1414 !       given by xlnsig. The integral of the
1415 !       asymmetry factor g is also calculated.
1416 ! FSB Change 12/20/2011 This code now has a choice of IGH based
1417 !     upon alfv and nr.
1418 ! FBB Changes Simplified code. Eliminated Penndorf code
1419 !  *** Does Gauss-Hermite quadrature of Qext / alfa & Qscat / alfa
1420 !      and asymmetry factor <cos> over log normal distribution using
1421 !      symmetric  points.
1422 !
1423     implicit none
1424     complex, intent(in) :: RCORE      ! refractive index of core
1425     complex, intent(in) :: RSHELL     ! refractive index of shell
1426     real, intent(in)    :: XX         ! Mie parameter for core
1427     real, intent(in)    :: YY         ! Mie parameter for shell
1428     real, intent(in)    :: xlnsig     ! log of geometric  standard deviation
1429     real, intent(out)   :: Qext_GH    ! normalized extinction efficiency
1430     real, intent(out)   :: Qscat_GH   ! normalized scattering efficiency
1431     real, intent(out)   :: g_GH       ! asymmetry factor <cos>
1432     logical, intent(out) :: success   ! flag for successful calculation
1433
1434     real    :: nr                     ! real part of  refractive index
1435     real    :: aa1                    ! see below for definition
1436     real    :: XXP, XXM               ! Mie parameters at abscissas - CORE
1437     real    :: YYP, YYM               ! Mie parameters at abscissas - SHELL
1438
1439 ! FSB define parameters
1440    real, parameter :: pi = 3.14159265
1441    real, parameter :: sqrtpi = 1.772454
1442    real, parameter :: sqrtpi1 = 1.0 / sqrtpi
1443    real, parameter :: sqrt2 = 1.414214
1444    real, parameter :: three_pi_two = 3.0 * pi / 2.0
1445    real, parameter ::  const = three_pi_two * sqrtpi1
1446
1447 !  *** these are Qext/alfa and Qscat/alfv at the abscissas
1448     real    :: qalfip_e, qalfim_e     ! extinction
1449     real    :: qalfip_s, qalfim_s     ! scattering
1450     real    :: gsalfp, gsalfm         ! scattering times asymmetry factor
1451     integer :: IGH                    ! index for GH quadrature
1452     integer ::  i
1453     real    ::  sum_e,sum_s, xi,wxi,xf, temp
1454     real    ::  sum_sg
1455
1456 ! Gauss-Hermite abscissas and weights
1457 ! *** the following weights and abscissas are from Abramowitz
1458 !     Stegun, Table 25.10 page 924
1459 ! FSB full precision from Table 25.10
1460
1461 ! FSB ten-point  - IGH = 5
1462     real, parameter :: ghxi_10(5) = (/ 0.342901327223705,     &
1463                                        1.036610829789514,     &
1464                                        1.756683649299882,     &
1465                                        2.532731674232790,     &
1466                                        3.436159118837738 /)
1467
1468     real, parameter :: ghwi_10(5) = (/ 6.108626337353e-01,    &
1469                                        2.401386110823e-01,    &
1470                                        3.387439445548e-02,    &
1471                                        1.343645746781e-03,    &
1472                                        7.640432855233e-06 /)
1473
1474 ! FSB six-point - IGH = 3
1475     real, parameter :: ghxi_6(3) = (/ 0.436077411927617,      &
1476                                       1.335849074013597,      &
1477                                       2.350604973674492 /)
1478
1479     real, parameter :: ghwi_6(3) = (/ 7.246295952244e-01,     &
1480                                       1.570673203229e-01,     &
1481                                       4.530009905509e-03 /)
1482
1483 ! FSB two-point - IGH = 1
1484     real, parameter :: ghxi_2(1) = (/ 0.707106781186548 /)
1485
1486     real, parameter :: ghwi_2(1) = (/ 8.862269254528e-01 /)
1487
1488     real GHXI(5), GHWI(5) ! weight and abscissas
1489     integer NMAX ! number of weights and abscissa
1490
1491 ! start code
1492 ! FSB now choose IGH. These choices are designed to improve
1493 !     the computational efficiency without sacrificing accuracy.
1494
1495     nr = real(RSHELL)
1496
1497     IGH=3 ! default value; six_point is sufficient generally
1498 ! six point
1499     NMAX = 3
1500
1501     if (nr .ge. 1.7) then
1502 ! 10 point
1503        IGH = 5 ! more points needed here
1504        NMAX = 5
1505     end if
1506
1507     if ( XX .gt. 20.0 .or. XX .lt. 0.5 ) then
1508        IGH  = 1 ! in  this range fewer points are needed
1509        NMAX = 1
1510     end if
1511
1512     if (IGH == 1) then
1513 ! two point
1514        GHXI(1)    = ghxi_2(1)
1515        GHWI(1)    = ghwi_2(1)
1516     else if (IGH == 3) then
1517        do i = 1, NMAX
1518           GHXI(i) = ghxi_6(i)
1519           GHWI(i) = ghwi_6(i)
1520        end do
1521     else
1522        do i = 1,NMAX
1523           GHXI(i) = ghxi_10(i)
1524           GHWI(i) = ghwi_10(i)
1525        end do
1526     end if ! set up number of abscissas and weights
1527
1528 ! FSB now start the integration code
1529     aa1 = sqrt2 * xlnsig   ! This 1.0 / Sqrt( A ) in derivation of the integral
1530                            ! where A = 1.0 / ( 2.0 * xlnsg**2 )
1531
1532 ! Then alpha = alfv * exp[ u / sqrt(A) ]
1533 ! For Gauss-Hermite Quadrature u = xi
1534 ! Therefore, xf = exp( xi / sqrt(A) ),
1535 !  or xf = exp( xi * aa1 )
1536     sum_e  = 0.0
1537     sum_s  = 0.0
1538     sum_sg = 0.0
1539 ! FSB do NMAX calls to the MIE codes
1540     do i = 1,NMAX
1541        xi      = GHXI(i)
1542        wxi     = GHWI(i)
1543        xf      = exp( xi * aa1 )
1544        temp    = 1.0 / xf
1545        XXP     = XX * xf
1546        XXM     = XX * temp ! division cheaper than another exp()
1547        YYP     = YY * xf
1548        YYM     = YY * temp ! division cheaper than another exp()
1549 ! *** call subroutine to fetch the effficiencies
1550
1551        call getqsgBHCS(XXP,YYP,RCORE,RSHELL,qalfip_e,qalfip_s,gsalfp, success)
1552        call getqsgBHCS(XXM,YYM,RCORE,RSHELL,qalfim_e,qalfim_s,gsalfm, success)
1553
1554        sum_e  = sum_e  + wxi  * ( qalfip_e + qalfim_e )
1555        sum_s  = sum_s  + wxi  * ( qalfip_s + qalfim_s )
1556        sum_sg = sum_sg + wxi  * ( gsalfp   + gsalfm   )
1557     end do
1558
1559     g_GH     = sum_sg / sum_s ! this is <cos>
1560     Qext_GH  = const * sum_e  !
1561     Qscat_GH = const * sum_s
1562
1563   end subroutine ghintBH_CS_even
1564
1565 ! ------------------------------------------------------------------
1566   subroutine getqsgBHCS (XX,YY,RRFRL1,RRFRL2,qxtalf,qscalf,qsgalf, success)
1567     implicit none
1568
1569     real, intent(in)    :: XX, YY
1570     real, intent(out)   :: qxtalf, qscalf, qsgalf
1571     complex, intent(in) :: RRFRL1,RRFRL2            ! refractive indices Core , Shell
1572     logical, intent(out) :: success                 ! flag for successful calculation
1573
1574     real    :: QEXT, QSCA, QBACK, G_MIE
1575     real    :: xx1
1576     character (len = 20) :: mystr1, mystr2, mystr3, mystr4
1577
1578     xx1    = 1.0 / YY
1579
1580 !      if (     (xx *  real(RRFRL1) >= 30.0)   &
1581 !          .or. (xx * aimag(RRFRL1) >= 30.0)   &
1582 !          .or. (yy * aimag(RRFRL2) >= 30.0)) then
1583 !         print *, ' ==d== bhcoat error'
1584 !      end if
1585
1586     call BHCOAT (XX,YY,RRFRL1,RRFRL2,QEXT,QSCA,QBACK,G_MIE, SUCCESS)
1587
1588 !      if ((trim(mystr1) == ' NaN') .or.  &
1589 !          (trim(mystr2) == ' NaN') .or.  &
1590 !          (trim(mystr3) == ' NaN') .or.  &
1591 !          (trim(mystr4) == ' NaN')) then
1592 !          call BHCOAT (XX,YY,RRFRL1,RRFRL2,QEXT,QSCA,QBACK,G_MIE)
1593 !      end if
1594
1595     qxtalf = QEXT * xx1
1596     qscalf = QSCA * xx1
1597     qsgalf = qscalf * G_MIE
1598
1599   END subroutine getqsgBHCS
1600
1601 ! ------------------------------------------------------------------
1602   SUBROUTINE BHCOAT (XX, YY, RRFRL1, RRFRL2, QQEXT, QQSCA, QBACK, GGSCA, SUCCESS)
1603
1604     use complex_number_module
1605
1606     implicit none ! added by FSB
1607
1608 ! Arguments:
1609      real, intent(in)    :: XX,YY             ! Defined below
1610      complex, intent(in) :: RRFRL1,RRFRL2     ! Defined below
1611      real, intent(out)   :: QQEXT,QQSCA,QBACK ! Defined below
1612      real, intent(out)   :: GGSCA             ! asymmetry factor <cos> added by FSB
1613      logical,intent(out) :: success
1614
1615 ! Local variables:
1616
1617      real*8, parameter     :: DEL = 1.0D-08
1618      real*8, parameter     :: ONE = 1.0D0, TWO = 2.0D0
1619 !    complex*16, save :: II
1620 !    data II/(0.D0,1.D0)/
1621      type(complex_number) :: II
1622
1623      integer :: IFLAG,N,NSTOP
1624
1625      character (len = 128) :: mystr
1626
1627 !         -----------------------------------------------------------
1628 !              del is the inner sphere convergence criterion
1629 !         -----------------------------------------------------------
1630
1631      real*8 :: CHI0Y,CHI1Y,CHIY,PSI0Y,PSI1Y,PSIY,QEXT,RN,QSCA,X,Y,YSTOP,GSCA
1632      real*8 :: TWO_N_M_ONE, TWO_N_P_ONE
1633      real*8 :: RY, RYY, RNRY, RN1, factor
1634
1635 !    complex*16 :: AMESS1,AMESS2,AMESS3,AMESS4,AN,ANCAP,AN1, BN,BNCAP,BN1, BRACK,   &
1636      type(complex_number) :: AMESS1,AMESS2,AMESS3,AMESS4,AN,ANCAP,AN1, BN,BNCAP,BN1, BRACK,   &
1637                              CHI0X2,CHI0Y2,CHI1X2,CHI1Y2,CHIX2,CHIPX2,CHIPY2,CHIY2,CRACK,     &
1638                              D0X1,D0X2,D0Y2,D1X1,D1X2,D1Y2,DNBAR,GNBAR,                       &
1639                              REFREL,RFREL1,RFREL2, XBACK,XI0Y,XI1Y,XIY,                       &
1640                              X1,X2,Y2,RCX1, RCX2,RCY2, FAC1, FAC2
1641
1642 !***********************************************************************
1643 ! NOTES from Prof. Bruce T. Draine, Princeton University
1644 ! Subroutine BHCOAT calculates Q_ext, Q_sca, Q_back for coated sphere.
1645 ! All bessel functions computed by upward recurrence.
1646 ! Input:
1647 !        XX = 2*PI*RCORE*REFMED/WAVEL
1648 !        YY = 2*PI*RMANT*REFMED/WAVEL
1649 !        RFREL1 = REFCOR/REFMED
1650 !        RFREL2 = REFMAN/REFMED
1651 ! where  REFCOR = complex refr.index of core)
1652 !        REFMAN = complex refr.index of mantle)
1653 !        REFMED = real refr.index of medium)
1654 !        RCORE = radius of core
1655 !        RMANT = radius of mantle
1656 !        WAVEL = wavelength of light in ambient medium
1657 !
1658 ! Routine BHCOAT is taken from Bohren & Huffman (1983)
1659 ! Obtained from C.L.Joseph
1660 !
1661 ! History:
1662 ! 92/11/24 (BTD) Explicit declaration of all variables
1663 ! April 30,2012 (FSB) added additional code to optimize
1664 !  run time by finding common terms and replacing multiple
1665 !  divisions by multiplication by a reciprocal.
1666 ! April 09, 2012  code transferred from BTD's BMHMIE to
1667 ! calculate the asymmetry factor by  Prof. Francis S. Binkowski of
1668 ! The University of North Carolina at Chapel Hill.
1669 ! April 30,2012 (FSB) added additional code to optimize
1670 !  run time by finding common terms and replacing multiple
1671 !  divisions by multiplication by a reciprocal.
1672 ! July 16, 2010  more optimization by Dr. David Wong (DW) at US EPA
1673
1674 ! REFERENCE:
1675 !  Bohren, Craig F. and Donald R. Huffman, Absorption and
1676 !    Scattering of Light by Small Particles, Wiley-Interscience
1677 !    copyright 1983. Paperback Published 1998.
1678 !    This code was originally listed in Appendix B. pp 483-489.
1679 !    As noted above , the original code was subsequently
1680 !    modified by Prof. Bruce T. Draine of Princeton University.
1681 !
1682 ! FSB The background for this code is discussed in Borhen & Huffman (1983)
1683 ! on pages 181-183 ( Equations 8.2 ) and on pages 483-484.
1684 !***********************************************************************
1685 !
1686 ! Start Code
1687
1688      SUCCESS = .TRUE.
1689
1690 #if ( RWORDSIZE == 8 )
1691      II = c_set(0.0, 1.0)
1692 #else
1693      II = c_set(0.0D0, 1.0D0)
1694 #endif
1695
1696 ! this technique will make the second 4 byte in the 8 byte variable be 0
1697 ! rather than arbitrary digits to increase accuracy
1698      write (mystr, *) xx, yy, real(RRFRL1), aimag(RRFRL1), real(RRFRL2), aimag(RRFRL2)
1699      read  (mystr, *) x,  y,  RFREL1, RFREL2
1700
1701 !    X      = XX
1702 !    Y      = YY
1703      RY     = ONE / Y
1704      RYY    = RY * RY
1705 !    RFREL1%real_part = real(RRFRL1)
1706 !    RFREL1%imag_part = aimag(RRFRL1)
1707 !    RFREL2%real_part = real(RRFRL2)
1708 !    RFREL2%imag_part = aimag(RRFRL2)
1709      x1     = c_mul(x, rfrel1)
1710      x2     = c_mul(x, rfrel2)
1711      y2     = c_mul(y, rfrel2)
1712      RCX1   = c_div(ONE, X1)
1713      RCX2   = c_div(ONE, X2)
1714      RCY2   = c_div(ONE, Y2)
1715      refrel = c_div(rfrel2, rfrel1)
1716      ystop  = y + 4.0 * y**0.3333 + 2.0
1717      nstop  = INT( ystop )
1718
1719 !         -----------------------------------------------------------
1720 !              series terminated after nstop terms
1721 !         -----------------------------------------------------------
1722
1723 !   initialize variables
1724      d0x1   = c_div(c_cos(x1), c_sin(x1))
1725      d0x2   = c_div(c_cos(x2), c_sin(x2))
1726      d0y2   = c_div(c_cos(y2), c_sin(y2))
1727
1728      psi0y  = cos(y)
1729      psi1y  = sin(y)
1730      chi0y  = -sin(y)
1731      chi1y  = cos(y)
1732
1733      xi0y   = c_sub(psi0y, c_mul(chi0y, II))
1734      xi1y   = c_sub(psi1y, c_mul(chi1y, II))
1735
1736 #if ( RWORDSIZE == 8 )
1737      chi0y2 = c_mul(-1.0, c_SIN(y2))
1738 #else
1739      chi0y2 = c_mul(-1.0d0, c_SIN(y2))
1740 #endif
1741      chi1y2 = c_COS(y2)
1742 #if ( RWORDSIZE == 8 )
1743      chi0x2 = c_mul(-1.0, c_SIN(x2))
1744 #else
1745      chi0x2 = c_mul(-1.0d0, c_SIN(x2))
1746 #endif
1747      chi1x2 = c_COS(x2)
1748      qsca   = 0.0d0
1749      qext   = 0.0d0
1750      GSCA   = 0.0d0
1751 #if ( RWORDSIZE == 8 )
1752      xback  = c_set(0.0, 0.0)
1753 #else
1754      xback  = c_set(0.0d0, 0.0d0)
1755 #endif
1756      iflag  = 0
1757      factor = 1.0d0
1758
1759 ! FSB Start main loop
1760      DO n = 1, nstop
1761         rn = REAL( n, 8 )
1762         RN1 = ONE / RN
1763         TWO_N_M_ONE = TWO * RN - ONE
1764         TWO_N_P_ONE = TWO * RN + ONE
1765         psiy = (TWO_N_M_ONE)*psi1y*RY - psi0y
1766         chiy = (TWO_N_M_ONE)*chi1y*RY - chi0y
1767         xiy  = c_sub(psiy, c_mul(chiy, II))
1768         d1y2 = c_sub(c_div(ONE, c_sub(c_mul(rn, RCY2), d0y2)), c_mul(rn, RCY2))
1769
1770         IF (iflag .eq. 0) THEN
1771 ! *** Calculate inner sphere ancap, bncap
1772 !      and brack and crack
1773            d1x1   = c_sub(c_div(ONE, c_sub(c_mul(rn, RCX1), d0x1)), c_mul(rn, RCX1))
1774            d1x2   = c_sub(c_div(ONE, c_sub(c_mul(rn, RCX2), d0x2)), c_mul(rn, RCX2))
1775
1776            chix2  = c_sub(c_mul(c_mul(TWO*rn - ONE, chi1x2), RCX2), chi0x2)
1777            chiy2  = c_sub(c_mul(c_mul(TWO*rn - ONE, chi1y2), RCY2), chi0y2)
1778
1779            chipx2 = c_sub(chi1x2, c_mul(c_mul(rn, chix2), RCX2))
1780            chipy2 = c_sub(chi1y2, c_mul(c_mul(rn, chiy2), RCY2))
1781
1782            ANCAP = c_sub(c_mul(c_mul(REFREL, D1X1), CHIX2), CHIPX2)
1783            ANCAP = c_mul(ANCAP, c_sub(c_mul(CHIX2, D1X2), CHIPX2))
1784            ANCAP = c_div(c_sub(c_mul(REFREL, D1X1), D1X2), ANCAP)
1785
1786            brack  = c_mul(ancap, c_sub(c_mul(chiy2, d1y2), chipy2))
1787
1788            bncap  = c_sub(c_mul(refrel, d1x2), d1x1)
1789            bncap  = c_div(bncap, c_sub(c_mul(refrel, chipx2), c_mul(d1x1, chix2)))
1790            bncap  = c_div(bncap, c_sub(c_mul(chix2, d1x2), chipx2))
1791
1792            crack  = c_mul(bncap, c_sub(c_mul(chiy2, d1y2), chipy2))
1793 ! *** calculate convergence test expressions
1794 !     for inner sphere.
1795 ! *** see pages 483-485 of Bohren & Huffman for
1796 !     definitions.
1797            amess1 = c_mul(brack, chipy2)
1798            amess2 = c_mul(brack, chiy2)
1799            amess3 = c_mul(crack, chipy2)
1800            amess4 = c_mul(crack, chiy2)
1801
1802 ! Now test for convergence for inner sphere
1803 !  All four criteria must be satisfied. See page 484 of B & H
1804            IF (c_ABS(amess1) .LE. del*c_ABS(d1y2)  .AND.                          &
1805               (c_ABS(amess2) .LE. del)             .AND.                          &
1806               (c_ABS(amess3) .LE. del*c_ABS(d1y2)) .AND.                          &
1807               (c_ABS(amess4) .LE. del)                ) THEN
1808 !             convergence for inner sphere
1809 #if ( RWORDSIZE == 8 )
1810               brack = c_set(0.0,0.0)
1811               crack = c_set(0.0,0.0)
1812 #else
1813               brack = c_set(0.0D0,0.0D0)
1814               crack = c_set(0.0D0,0.0D0)
1815 #endif
1816               iflag = 1
1817 !         ELSE
1818 ! no convergence yet
1819 !            iflag = 0
1820           END IF
1821         END IF ! test on iflag .eq. 0
1822
1823 ! *** note usage of brack and crack See equations on
1824 !     Page 485  and discussion on pages 486 -487 of B & H
1825         dnbar = c_sub(d1y2, c_mul(brack, chipy2))
1826         dnbar = c_div(dnbar, c_sub(ONE, c_mul(brack, chiy2)))
1827         gnbar = c_sub(d1y2, c_mul(crack, chipy2))
1828         gnbar = c_div(gnbar, c_sub(ONE, c_mul(crack, chiy2)))
1829 !*** Store previous values of an and bn for use
1830 !    in computation of g=<cos(theta)>
1831         IF (N .GT. 1) THEN
1832            AN1 = an
1833            BN1 = bn
1834         END IF
1835 ! *** update an and bn
1836         RNRY = rn * RY
1837         FAC1 = c_add(c_div(dnbar, rfrel2), RNRY)
1838
1839         an = c_sub(c_mul(psiy, FAC1), psi1y)
1840         an = c_div(an, c_sub(c_mul(FAC1, xiy), xi1y))
1841         FAC2 = c_add(c_mul(rfrel2, gnbar), RNRY)
1842         bn = c_sub(c_mul(psiy, FAC2), psi1y)
1843         bn = c_div(bn, c_sub(c_mul(FAC2, xiy), xi1y))
1844
1845 ! *** Calculate sums for qsca, qext, xback
1846         qsca  = qsca + (TWO_N_P_ONE) * (c_ABS(an)**2 + c_ABS(bn)**2)
1847
1848         qext  = qext + TWO_N_P_ONE * (an%real_part + bn%real_part)
1849
1850         FACTOR = FACTOR * (-1.0D0)
1851         XBACK = c_add(XBACK, c_mul(TWO_N_P_ONE * FACTOR, c_sub(AN, BN)))
1852
1853 ! FSB calculate the sum for the asymmetry factor
1854
1855         GSCA = GSCA + ((TWO_N_P_ONE)/(RN* (RN + ONE)))*                       &
1856                (an%real_part*bn%real_part + an%imag_part*bn%imag_part)
1857
1858         IF (n .GT. 1) THEN
1859
1860            GSCA = GSCA + (RN - RN1) *                                         &
1861                  (AN1%real_part*AN%real_part + AN1%imag_part*AN%imag_part +   &
1862                   BN1%real_part*BN%real_part + BN1%imag_part*BN%imag_part)
1863
1864         END IF
1865 ! continue update for next interation
1866         psi0y  = psi1y
1867         psi1y  = psiy
1868         chi0y  = chi1y
1869         chi1y  = chiy
1870         xi1y   = c_sub(psi1y, c_mul(chi1y, II))
1871         chi0x2 = chi1x2
1872         chi1x2 = chix2
1873         chi0y2 = chi1y2
1874         chi1y2 = chiy2
1875         d0x1   = d1x1
1876         d0x2   = d1x2
1877         d0y2   = d1y2
1878      END DO  ! end of main loop
1879
1880 !*** Have summed sufficient terms.
1881 !    Now compute QQSCA,QQEXT,QBACK,and GSCA
1882      GGSCA = REAL( TWO * GSCA / qsca )
1883      QQSCA = REAL( TWO * qsca * RYY )
1884      QQEXT = REAL( TWO * qext * RYY )
1885
1886      QBACK = 0.5 * real((xback%real_part**2 + xback%imag_part**2) * RYY)
1887
1888   end subroutine BHCOAT
1889
1890 ! ------------------------------------------------------------------
1891   subroutine ghintBH_Odd (INIT, crefin,alfv,xlnsig,Qext_GH,Qscat_GH,g_gh, success )
1892
1893 ! *************** REVISED VERSION < NOTE
1894 ! FSB *********** This is the newest (04_14_2012) version of GhintBH
1895 !      this version does the Mie method and calculates the optimum set of
1896 !      set of Gauss-Hermite abscissas and weights.
1897 !      Dr. Francis S. Binkowski, The University of North Carolina
1898 !                                at Chapel Hill
1899 ! FSB this code file now contains all of the necessary subroutines that
1900 !     are called to perform an integral of the Bohren and Huffman
1901 !     Mie codes ( as updated by Prof. Bruce C. Drain of Princeton)
1902 !       calculates the extinction and scattering coefficients
1903 !       normalized by wavelength and total particle volume
1904 !       concentration for a log normal particle distribution
1905 !       with the logarithm of the geometric  standard deviation
1906 !       given by xlnsig. The integral of the
1907 !       asymmetry factor g is also calculated.
1908 ! FSB Change 12/20/2011 This code now has a choice of IGH based
1909 !     upon alfv and nr.
1910 ! FBB Changes Simplified code. Eliminated Penndorf code
1911 !  *** Does Gauss-Hermite quadrature of Qext / alfa & Qscat / alfa
1912 !      and asymmetry factor <cos> over log normal distribution using
1913 !      symmetric  points.
1914 !
1915     implicit none
1916
1917     logical, intent(INOUT)        :: INIT       ! initialize number of qudraure points
1918     complex, intent(in)           :: crefin     ! complex index of refraction
1919     real, intent(in)              :: alfv       ! Mie parameter for dgv
1920     real, intent(in)              :: xlnsig     ! log of geometric  standard deviation
1921     real, intent(out)             :: Qext_GH    ! normalized extinction efficiency
1922     real, intent(out)             :: Qscat_GH   ! normalized scattering efficiency
1923     real, intent(out)             :: g_GH       ! asymmetry factor <cos>
1924     logical, intent(out)          :: success    ! flag for successful calculation
1925
1926     real    :: nr                 ! real part of  refractive index
1927     real    :: aa1                ! see below for definition
1928     real    :: alfaip, alfaim     ! Mie parameters at abscissas
1929
1930 !  *** these are Qext/alfa and Qscat/alfv at the abscissas
1931     real    :: qalfip_e, qalfim_e ! extinction
1932     real    :: qalfip_s, qalfim_s ! scattering
1933     real    :: gsalfp, gsalfm     ! scattering times asymmetry factor
1934
1935 ! FSB define parameters
1936     real, parameter :: pi = 3.14159265
1937     real, parameter :: sqrtpi = 1.772454
1938     real, parameter :: sqrtpi1 = 1.0 / sqrtpi
1939     real, parameter :: sqrt2 = 1.414214
1940     real, parameter :: three_pi_two = 3.0 * pi / 2.0
1941     real, parameter :: const = three_pi_two * sqrtpi1
1942
1943     integer ::  i
1944     real    ::  sum_e,sum_s, xi,wxi,xf
1945     real    ::  sum_sg
1946
1947     real,    allocatable,  save  :: GHXI(:), GHWI(:) ! weight and abscissas
1948     integer, save  :: IGH                            ! number of weights and abscissa
1949     integer, save  :: NMAX                           ! optimumized number of weights and abscissa
1950
1951
1952 ! start code
1953 ! FSB now choose IGH. These choices are designed to improve
1954 !     the computational efficiency without sacrificing accuracy.
1955
1956     If( INIT )Then
1957
1958        Select Case( Quadrature_Points )
1959          Case( 1,3,9 )
1960            IGH = Quadrature_Points
1961          Case Default
1962            IGH = 3
1963        End Select
1964
1965        NMAX = Max( Int( IGH / 2 ), 0)
1966
1967        If( Allocated( GHXI ) .Or. Allocated( GHWI ) )Then
1968            Success = .False.
1969            Return
1970        End If
1971
1972        Allocate( GHXI( NMAX + 1 ), GHWI( NMAX + 1 ) )
1973
1974        Select Case ( IGH )
1975          Case ( 1 )
1976            GHXI(1)  = ghxi_1(1)
1977            GHWI(1)  = ghwi_1(1)
1978          Case ( 3 )
1979            do i = 1, NMAX + 1
1980              GHXI(i) = ghxi_3(i)
1981              GHWI(i) = ghwi_3(i)
1982            end do
1983          Case ( 9 )
1984            do i = 1, NMAX + 1
1985              GHXI(i) = ghxi_9(i)
1986              GHWI(i) = ghwi_9(i)
1987            end do
1988        end select
1989
1990        If( AERO_UTIL_LOG .GT. 0 )Then
1991            write(AERO_UTIL_LOG,*)'BHMIE: IGH,(NMAX + 1) = ',IGH,(NMAX + 1)
1992            do i = 1, NMAX + 1
1993              write(AERO_UTIL_LOG,*)'BHMIE: i, GHXI(i), GHWI(i) = ',i, GHXI(i), GHWI(i)
1994            end do
1995        End If
1996
1997        INIT = .False.
1998     Else
1999        If( .Not. Allocated( GHXI ) .Or. .Not. Allocated( GHWI ) )Then
2000            Success = .False.
2001            Return
2002        End If
2003     End If ! set up number of abscissas and weights
2004
2005     nr = real(crefin)
2006
2007 ! FSB now start the integration code
2008     aa1 = sqrt2 * xlnsig    ! This 1.0 / Sqrt( A ) in derivation of the integral
2009                             ! where A = 1.0 / ( 2.0 * xlnsg**2 )
2010
2011 ! Then alpha = alfv * exp[ u / sqrt(A) ]
2012 ! For Gauss-Hermite Quadrature u = xi
2013 ! Therefore, xf = exp( xi / sqrt(A) ),
2014 !  or xf = exp( xi * aa1 )
2015
2016 !start integration at zero point
2017     xi      = 0.0
2018     wxi     = GHWI(NMAX+1)
2019     xf      = 1.0
2020     alfaip  = alfv
2021 ! fetch the effficiencies at zero point
2022
2023     call getqext_BH(alfaip,crefin,qalfip_e,qalfip_s, gsalfp, success)
2024
2025     sum_e  = wxi * qalfip_e
2026     sum_s  = wxi * qalfip_s
2027     sum_sg = wxi * gsalfp
2028
2029 ! FSB do NMAX calls to the MIE codes
2030     do i = 1, NMAX
2031        xi      = GHXI(i)
2032        wxi     = GHWI(i)
2033        xf      = exp( xi * aa1 )
2034        alfaip  = alfv * xf
2035        alfaim  = alfv / xf ! division cheaper than another exp()
2036 ! *** call subroutine to fetch the effficiencies
2037
2038        call getqext_BH(alfaip,crefin,qalfip_e,qalfip_s, gsalfp, success)
2039        call getqext_BH(alfaim,crefin,qalfim_e,qalfim_s, gsalfm, success)
2040
2041        sum_e  = sum_e + wxi  * ( qalfip_e + qalfim_e )
2042        sum_s  = sum_s + wxi  * ( qalfip_s + qalfim_s )
2043        sum_sg = sum_sg + wxi * ( gsalfp + gsalfm )
2044
2045     end do
2046
2047     g_GH     = sum_sg / sum_s ! this is <cos>
2048     Qext_GH  = const * sum_e  !
2049     Qscat_GH = const * sum_s
2050
2051   end subroutine ghintBH_Odd
2052
2053 ! ------------------------------------------------------------------
2054   subroutine ghintBH_CS_Odd (INIT, RCORE, RSHELL , XX, YY, xlnsig,  &
2055                              Qext_GH,Qscat_GH, g_gh, success)
2056
2057 ! FSB code for coated-sphere (core-shell) version
2058
2059 ! *************** REVISED VERSION < NOTE
2060 ! FSB *********** This is the newest (04_14_2012) version of ghintBH_CS
2061 !      for the coated-sphere (core-shell) method using BHCOAT
2062 !      this version does the Mie method and calculates the optimum set of
2063 !      set of Gauss-Hermite abscissas and weights.
2064 !      Dr. Francis S. Binkowski, The University of North Carolina
2065 !                                at Chapel Hill
2066
2067 ! FSB this code file now contains all of the necessary subroutines that
2068 !     are called to perform an integral of the Bohren and Huffman
2069 !     Mie codes ( as updated by Prof. Bruce C. Drain of Princeton)
2070 !       calculates the extinction and scattering coefficients
2071 !       normalized by wavelength and total particle volume
2072 !       concentration for a log normal particle distribution
2073 !       with the logarithm of the geometric  standard deviation
2074 !       given by xlnsig. The integral of the
2075 !       asymmetry factor g is also calculated.
2076 ! FSB Change 12/20/2011 This code now has a choice of IGH based
2077 !     upon alfv and nr.
2078 ! FBB Changes Simplified code. Eliminated Penndorf code
2079 !  *** Does Gauss-Hermite quadrature of Qext / alfa & Qscat / alfa
2080 !      and asymmetry factor <cos> over log normal distribution using
2081 !      symmetric  points.
2082 !
2083     implicit none
2084
2085     logical, intent(inout) :: INIT       ! initialize number of qudraure points
2086     complex, intent(in)    :: RCORE      ! refractive index of core
2087     complex, intent(in)    :: RSHELL     ! refractive index of shell
2088     real, intent(in)       :: XX         ! Mie parameter for core
2089     real, intent(in)       :: YY         ! Mie parameter for shell
2090     real, intent(in)       :: xlnsig     ! log of geometric  standard deviation
2091     real, intent(out)      :: Qext_GH    ! normalized extinction efficiency
2092     real, intent(out)      :: Qscat_GH   ! normalized scattering efficiency
2093     real, intent(out)      :: g_GH       ! asymmetry factor <cos>
2094     logical, intent(out)   :: success   ! flag for successful calculation
2095
2096     real    :: nr                     ! real part of  refractive index
2097     real    :: aa1                    ! see below for definition
2098     real    :: XXP, XXM               ! Mie parameters at abscissas - CORE
2099     real    :: YYP, YYM               ! Mie parameters at abscissas - SHELL
2100
2101 ! FSB define parameters
2102    real, parameter :: pi = 3.14159265
2103    real, parameter :: sqrtpi = 1.772454
2104    real, parameter :: sqrtpi1 = 1.0 / sqrtpi
2105    real, parameter :: sqrt2 = 1.414214
2106    real, parameter :: three_pi_two = 3.0 * pi / 2.0
2107    real, parameter ::  const = three_pi_two * sqrtpi1
2108
2109 !  *** these are Qext/alfa and Qscat/alfv at the abscissas
2110     real    :: qalfip_e, qalfim_e     ! extinction
2111     real    :: qalfip_s, qalfim_s     ! scattering
2112     real    :: gsalfp, gsalfm         ! scattering times asymmetry factor
2113     integer ::  i
2114     real    ::  sum_e,sum_s, xi,wxi,xf, temp
2115     real    ::  sum_sg
2116
2117     real,    allocatable,  save  :: GHXI(:), GHWI(:) ! weight and abscissas
2118     integer,               save  :: IGH              ! number of weights and abscissa
2119     integer,               save  :: NMAX             ! optimized number of weights and abscissa
2120
2121 ! start code
2122 ! FSB now choose IGH. These choices are designed to improve
2123 !     the computational efficiency without sacrificing accuracy.
2124
2125     If( INIT )Then
2126
2127        Select Case( Quadrature_Points )
2128          Case( 1,3,9 )
2129            IGH = Quadrature_Points
2130          Case Default
2131            IGH = 3
2132        End Select
2133
2134        If( Allocated( GHXI ) .Or. Allocated( GHWI ) )Then
2135            Success = .False.
2136            Return
2137        End If
2138
2139        NMAX = Max( Int( IGH / 2 ), 0)
2140
2141        Allocate( GHXI( NMAX + 1 ), GHWI( NMAX + 1 ) )
2142
2143        Select Case ( IGH )
2144          Case ( 1 )
2145            GHXI(1)  = ghxi_1(1)
2146            GHWI(1)  = ghwi_1(1)
2147          Case ( 3 )
2148            do i = 1, NMAX + 1
2149              GHXI(i) = ghxi_3(i)
2150              GHWI(i) = ghwi_3(i)
2151            end do
2152          Case ( 9 )
2153            do i = 1, NMAX + 1
2154              GHXI(i) = ghxi_9(i)
2155              GHWI(i) = ghwi_9(i)
2156            end do
2157        end select
2158
2159        If( AERO_UTIL_LOG .GT. 0 )Then
2160            write(AERO_UTIL_LOG,*)'BHCoat: IGH,(NMAX + 1) = ',IGH,(NMAX + 1)
2161            do i = 1, NMAX + 1
2162              write(AERO_UTIL_LOG,*)'BHCoat: i, GHXI(i), GHWI(i) = ',i, GHXI(i), GHWI(i)
2163            end do
2164        End If
2165
2166        INIT = .False.
2167
2168     Else
2169        If( .Not. Allocated( GHXI ) .Or. .Not. Allocated( GHWI ) )Then
2170            Success = .False.
2171            Return
2172        End If
2173     End If ! set up number of abscissas and weights
2174
2175     nr = real(RSHELL)
2176
2177 ! FSB now start the integration code
2178     aa1 = sqrt2 * xlnsig   ! This 1.0 / Sqrt( A ) in derivation of the integral
2179                            ! where A = 1.0 / ( 2.0 * xlnsg**2 )
2180
2181 ! Then alpha = alfv * exp[ u / sqrt(A) ]
2182 ! For Gauss-Hermite Quadrature u = xi
2183 ! Therefore, xf = exp( xi / sqrt(A) ),
2184 !  or xf = exp( xi * aa1 )
2185
2186 !start integration at zero point
2187
2188     xi      = 0.0
2189     wxi     = GHWI(NMAX+1)
2190     xf      = 1.0
2191     XXP     = XX
2192     YYP     = YY
2193
2194 ! fetch the effficiencies at zero point
2195
2196     call getqsgBHCS(XXP,YYP,RCORE,RSHELL,qalfip_e,qalfip_s,gsalfp, success)
2197
2198     sum_e  = wxi  * qalfip_e
2199     sum_s  = wxi  * qalfip_s
2200     sum_sg = wxi  * gsalfp
2201
2202 ! FSB do NMAX calls to the MIE codes
2203     do i = 1, NMAX
2204        xi      = GHXI(i)
2205        wxi     = GHWI(i)
2206        xf      = exp( xi * aa1 )
2207        temp    = 1.0 / xf
2208        XXP     = XX * xf
2209        XXM     = XX * temp ! division cheaper than another exp()
2210        YYP     = YY * xf
2211        YYM     = YY * temp ! division cheaper than another exp()
2212 ! *** call subroutine to fetch the effficiencies
2213
2214        call getqsgBHCS(XXP,YYP,RCORE,RSHELL,qalfip_e,qalfip_s,gsalfp, success)
2215        call getqsgBHCS(XXM,YYM,RCORE,RSHELL,qalfim_e,qalfim_s,gsalfm, success)
2216
2217        sum_e  = sum_e  + wxi  * ( qalfip_e + qalfim_e )
2218        sum_s  = sum_s  + wxi  * ( qalfip_s + qalfim_s )
2219        sum_sg = sum_sg + wxi  * ( gsalfp   + gsalfm   )
2220     end do
2221
2222     g_GH     = sum_sg / sum_s ! this is <cos>
2223     Qext_GH  = const * sum_e  !
2224     Qscat_GH = const * sum_s
2225
2226   end subroutine ghintBH_CS_Odd
2227
2228 ! ------------------------------------------------------------------
2229   SUBROUTINE BHMIE_FLEXI (X, NMX, NSTOP, REFREL, QQEXT, QQSCA, QBACK, GSCA, SUCCESS)
2230
2231 ! FSB Changed the call vector to return only QEXT, QSCAT QBACK GSCA
2232 !     and ignore NANG, S1 and S2 and all calculations for them
2233
2234     implicit none
2235
2236 ! Arguments:
2237     real,    intent(in) :: X        ! X = pi*particle_diameter / Wavelength
2238     integer, intent(in) :: NMX      ! maximum number of terms in Mie series
2239     integer, intent(in) :: NSTOP    ! minumum number of terms in Mie series
2240     complex, intent(in) :: REFREL   ! refractive index
2241
2242 !    REFREL = (complex refr. index of sphere)/(real index of medium)
2243 !    in the current use the index of refraction of the the medium
2244 !    i taken at 1.0 real.
2245 !
2246 !    Output
2247
2248     real,    intent(out) :: QQEXT, QQSCA, QBACK, GSCA
2249     logical, intent(out) :: SUCCESS
2250
2251 !     QQEXT   Efficiency factor for extinction
2252 !     QQSCA   Efficiency factor for scattering
2253 !     QQBACK  Efficiency factor for back scatter
2254 !     GSCA    asymmetry factor <cos>
2255 !     SUCCESS flag for successful calculation
2256 ! REFERENCE:
2257 !  Bohren, Craig F. and Donald R. Huffman, Absorption and
2258 !    Scattering of Light by Small Particles, Wiley-Interscience
2259 !    copyright 1983. Paperback Published 1998.
2260 ! FSB
2261 !    This code was originally listed in Appendix A. pp 477-482.
2262 !    As noted below, the original code was subsequently
2263 !    modified by Prof. Bruce T. Drain of Princetion University.
2264 !    The code was further modified for a specific application
2265 !    in a large three-dimensional code requiring as much
2266 !    computational efficiency as possible.
2267 !    Prof. Francis S. Binkowski of The University of North
2268 !    Carolina at Chapel Hill.
2269
2270 ! Declare parameters:
2271 ! Note: important that MXNANG be consistent with dimension of S1 and S2
2272 !       in calling routine!
2273
2274     integer, parameter    :: MXNANG=10, NMXX=150000   ! FSB new limits
2275     integer, parameter    :: NANG  = 2
2276     real*8, parameter     :: PII = 3.1415916536D0
2277     real*8, parameter     :: ONE = 1.0D0, TWO = 2.0D0
2278     complex*16, parameter :: COMPLEX_DZERO = (0.0D0,0.0D0)
2279     complex,    parameter :: COMPLEX_ZERO  = (0.0,0.0)
2280
2281 ! Local variables:
2282     integer    :: N, NN
2283     real*8     :: QSCA, QEXT, DX1, DXX1
2284     real*8     :: CHI,CHI0,CHI1,DX,EN,P,PSI,PSI0,PSI1,XSTOP,YMOD
2285     real*8     :: TWO_N_M_ONE, TWO_N_P_ONE, EN1, FACTOR
2286     complex*16 :: AN,AN1,BN,BN1,DREFRL,XI,XI1,Y, Y1, DREFRL1
2287     complex*16 :: D(NMX)
2288     complex*16 :: FAC1, FAC2
2289     complex*16 :: XBACK
2290
2291 !***********************************************************************
2292 ! Subroutine BHMIE is the Bohren-Huffman Mie scattering subroutine
2293 !    to calculate scattering and absorption by a homogenous isotropic
2294 !    sphere.
2295 ! Given:
2296 !    X = 2*pi*a/lambda
2297 !    REFREL = (complex refr. index of sphere)/(real index of medium)
2298 !    real refractive index of medium taken as 1.0
2299 ! Returns:
2300 !    QEXT  = efficiency factor for extinction
2301 !    QSCA  = efficiency factor for scattering
2302 !    QBACK = efficiency factor for backscatter
2303 !            see Bohren & Huffman 1983 p. 122
2304 !    GSCA = <cos> asymmetry for scattering
2305 !
2306 ! Original program taken from Bohren and Huffman (1983), Appendix A
2307 ! Modified by Prof. Bruce T.Draine, Princeton Univ. Obs., 90/10/26
2308 ! in order to compute <cos(theta)>
2309 ! 91/05/07 (BTD): Modified to allow NANG=1
2310 ! 91/08/15 (BTD): Corrected error (failure to initialize P)
2311 ! 91/08/15 (BTD): Modified to enhance vectorizability.
2312 ! 91/08/15 (BTD): Modified to make NANG=2 if called with NANG=1
2313 ! 91/08/15 (BTD): Changed definition of QBACK.
2314 ! 92/01/08 (BTD): Converted to full double precision and double complex
2315 !                 eliminated 2 unneed lines of code
2316 !                 eliminated redundant variables (e.g. APSI,APSI0)
2317 !                 renamed RN -> EN = double precision N
2318 !                 Note that DOUBLE COMPLEX and DCMPLX are not part
2319 !                 of f77 standard, so this version may not be fully
2320 !                 portable.  In event that portable version is
2321 !                 needed, use src/bhmie_f77.f
2322 ! 93/06/01 (BTD): Changed AMAX1 to generic function MAX
2323 ! FSB April 09,2012 This code was modified by:
2324 ! Prof.  Francis S. Binkowski University of North Carolina at
2325 ! Chapel Hill, Institue for the Environment.
2326 !
2327 ! The modifications were made to enhance computation speed
2328 ! for use in a three-dimensional code. This was done by
2329 ! removing code that calculated angular scattering. The method
2330 ! of calculating QEXT, QBACK was also changed.
2331
2332 !***********************************************************************
2333 !*** Safety checks
2334
2335     SUCCESS = .TRUE.
2336 !    NANG = 2 ! FSB only this value
2337 ! IF(NANG.GT.MXNANG)STOP'***Error: NANG > MXNANG in bhmie'
2338 !   IF (NANG .LT. 2) NANG = 2
2339
2340     DX = REAL( X, 8 )
2341 ! FSB Define reciprocals so that divisions can be replaced by multiplications.
2342     DX1  = ONE / DX
2343     DXX1 = DX1 * DX1
2344     DREFRL = DCMPLX( REAL( REFREL ), IMAG( REFREL ) )
2345     DREFRL1 = ONE / DREFRL
2346     Y = DX * DREFRL
2347     Y1 = ONE / Y
2348 !    YMOD = ABS(Y)
2349
2350 !*** Series expansion terminated after NSTOP terms
2351 !    Logarithmic derivatives calculated from NMX on down
2352 !       XSTOP = X + 4.0 * X**0.3333 + 2.0
2353 !       NMX  = MAX(XSTOP,YMOD) + 15
2354
2355 ! BTD experiment 91/1/15: add one more term to series and compare results
2356 !      NMX=AMAX1(XSTOP,YMOD)+16
2357 ! test: compute 7001 wavelengths between .0001 and 1000 micron
2358 ! for a=1.0micron SiC grain.  When NMX increased by 1, only a single
2359 ! computed number changed (out of 4*7001) and it only changed by 1/8387
2360 ! conclusion: we are indeed retaining enough terms in series!
2361
2362     FACTOR = 1.0D0
2363
2364 !       IF (NMX .GT. NMXX) THEN
2365 !          WRITE(6,*)'Error: NMX > NMXX=',NMXX,' for |m|x=',YMOD
2366 !          SUCCESS = .FALSE.
2367 !          RETURN
2368 !       END IF
2369
2370 ! FSB all code relating to scattering angles is removed out for
2371 !     reasons of efficiency when running in a three-dimensional
2372 !     code. We only need QQSCA, QQEXT, GSCA AND QBACK
2373
2374 !*** Logarithmic derivative D(J) calculated by downward recurrence
2375 !    beginning with initial value (0.,0.)
2376
2377     D(NMX) = COMPLEX_DZERO
2378     NN = NMX - 1
2379     DO N = 1,NN
2380        EN  = REAL( NMX - N + 1, 8 )
2381 ! FSB In the following division by Y has been replaced by
2382 !     multiplication by Y1, the reciprocal of Y.
2383        D(NMX-N) = ( EN * Y1 ) - (ONE / ( D(NMX-N+1) + EN * Y1))
2384     END DO
2385
2386 !*** Riccati-Bessel functions with real argument X
2387 !    calculated by upward recurrence
2388
2389     PSI0 =  COS(DX)
2390     PSI1 =  SIN(DX)
2391     CHI0 = -SIN(DX)
2392     CHI1 =  PSI0
2393     XI1  =  DCMPLX(PSI1,-CHI1)
2394     QSCA =  0.0D0
2395     GSCA =  0.0D0
2396     QEXT =  0.0D0
2397     P    = -ONE
2398     XBACK = COMPLEX_DZERO
2399
2400 ! FSB Start main loop
2401     DO N = 1,NSTOP
2402        EN        = REAL( N, 8 )
2403        EN1       = ONE / EN
2404        TWO_N_M_ONE = TWO * EN - ONE
2405 ! for given N, PSI  = psi_n        CHI  = chi_n
2406 !              PSI1 = psi_{n-1}    CHI1 = chi_{n-1}
2407 !              PSI0 = psi_{n-2}    CHI0 = chi_{n-2}
2408 ! Calculate psi_n and chi_n
2409        PSI = TWO_N_M_ONE * PSI1 * DX1 - PSI0
2410        CHI = TWO_N_M_ONE * CHI1 * DX1 - CHI0
2411        XI  = DCMPLX(PSI,-CHI)
2412
2413 !*** Compute AN and BN:
2414 ! FSB Rearrange to get common terms
2415        FAC1 = D(N) * DREFRL1 + EN * DX1
2416        AN   = (FAC1) * PSI - PSI1
2417        AN   = AN / ( (FAC1 )* XI - XI1 )
2418        FAC2 = ( DREFRL * D(N) + EN * DX1)
2419        BN   = ( FAC2) * PSI -PSI1
2420        BN   = BN / ((FAC2) * XI - XI1 )
2421
2422 ! FSB calculate sum for QEXT as done by Wiscombe
2423 !     get common factor
2424        TWO_N_P_ONE = (TWO * EN + ONE)
2425        QEXT = QEXT + (TWO_N_P_ONE) * (REAL(AN) + REAL(BN) )
2426        QSCA = QSCA + (TWO_N_P_ONE) * ( ABS(AN)**2 + ABS(BN)**2 )
2427
2428 ! FSB calculate XBACK from B & H Page 122
2429        FACTOR = -1.0d0 * FACTOR  ! calculate (-1.0 ** N)
2430        XBACK = XBACK + (TWO_N_P_ONE) * factor * (AN - BN)
2431
2432 ! FSB calculate asymmetry factor
2433
2434        GSCA = GSCA + REAL((TWO_N_P_ONE)/(EN * (EN + ONE)) *     &
2435               (REAL(AN)*REAL(BN)+IMAG(AN)*IMAG(BN)))
2436
2437        IF (N .GT. 1)THEN
2438           GSCA = GSCA + REAL((EN - EN1) *                         &
2439                  (REAL(AN1)*REAL(AN) + IMAG(AN1)*IMAG(AN) +  &
2440                   REAL(BN1)*REAL(BN) + IMAG(BN1)*IMAG(BN)))
2441        ENDIF
2442
2443 !*** Store previous values of AN and BN for use in computation of g=<cos(theta)>
2444        AN1 = AN
2445        BN1 = BN
2446
2447 ! FSB set up for next iteration
2448        PSI0 = PSI1
2449        PSI1 = PSI
2450        CHI0 = CHI1
2451        CHI1 = CHI
2452        XI1  = DCMPLX(PSI1,-CHI1)
2453
2454     END DO   ! main  loop on n
2455
2456 !*** Have summed sufficient terms.
2457
2458 !    Now compute QQSCA,QQEXT,QBACK,and GSCA
2459     GSCA  = REAL( TWO / QSCA ) * GSCA
2460
2461 ! FSB in the following, divisions by DX * DX has been replaced by
2462 !      multiplication by DXX1 the reciprocal of 1.0 / (DX *DX)
2463     QQSCA = REAL( TWO * QSCA * DXX1 )
2464     QQEXT = REAL( TWO * QEXT * DXX1 )
2465     QBACK = REAL( REAL( 0.5D0 * XBACK * CONJG(XBACK), 8 ) * DXX1 ) ! B&H Page 122
2466
2467   END subroutine BHMIE_FLEXI
2468
2469 END MODULE rrtmg_aero_optical_util_module
2470
2471 ! ------------------------------------------------------------------
2472 ! FSB REvised  Mie calculations  02/09/2011
2473
2474 MODULE module_ra_rrtmg_sw_cmaq
2475
2476 contains
2477
2478 !------------------------------------------------------------------
2479    Subroutine get_aerosol_Optics_RRTMG_SW ( ns, nmode,delta_z, INMASS_ws,     &
2480                                             INMASS_in, INMASS_ec, INMASS_ss,  &
2481                                             INMASS_h2o,  INDGN, INSIG,        &
2482                                             tauaer, waer, gaer )
2483
2484 !FSB This version switches between BHCOAT to BHMIE depending upon whether
2485 !    EC is present or not. 04/15/2012.
2486
2487 !FSB this version does a core-shell calculation with BHCOAT 04/11/2012
2488 ! This version is set up to be used with RRTMG_SW <<<<<<<<
2489 !     wavelenght is calculated internally
2490 ! FSB This routine calculates the aerosol information ( tauaer, waer,
2491 !     gaer, needed to calculate the solar radiation)  The calling
2492 !     program specifies the location ( row, column, layer,
2493 !     layer thicknes, and wave length for the calculation.
2494 ! FSB 02/09/2011 Modifications made to subroutine ghintBH.
2495 ! FSB 04/14/2012 REmoved MODULUS, made changes to ghintBH.
2496 !     Put in option for core-shell (coated-sphere). 2
2497
2498 ! FSB Input variables:
2499
2500       use rrtmg_aero_optical_util_module
2501
2502       implicit none
2503
2504       integer,intent(in) ::  ns      ! index for wavelength should be
2505                                      ! between 1 and 14. <<< RRTMG_SW
2506       integer,intent(in) ::  nmode   ! should be 3 for WRF/CMAQ calculation
2507       real,intent(in)    ::  delta_z ! layer thickness [m]
2508 ! FSB mode types for WRF/CMAQ
2509 !     nmode = 1 Aitken
2510 !     nmode = 2 accumulation
2511 !     nmode = 3 coarse
2512 ! FSB modal mass concentration by species  [ ug / m**3]  NOTE:  MKS
2513       real, intent(in) ::  INMASS_ws(nmode)   ! water soluble
2514       real, intent(in) ::  INMASS_in(nmode)   ! insolugle
2515       real, intent(in) ::  INMASS_ec(nmode)   ! elemental carbon or soot like
2516       real, intent(in) ::  INMASS_ss(nmode)   ! sea salt
2517       real, intent(in) ::  INMASS_h2o(nmode)  ! water
2518 ! FSB particle size-distribution information
2519       real, intent(in) ::  INDGN( nmode)      ! geometric mean diameter [ m ] NOTE: MKS
2520       real, intent(in) ::  INSIG( nmode)      ! geometric standard deviation
2521
2522 !FSB output aerosol radiative properties  [dimensionless]
2523       real, intent(out) ::  tauaer   ! aerosol extinction optical depth
2524       real, intent(out) ::  waer     ! aerosol single scattering albedo
2525       real, intent(out) ::  gaer     ! aerosol assymetry parameter
2526
2527 ! FSB Internal variables
2528
2529       real    :: NR(nmode), NI(nmode)           ! refractive indices
2530       complex :: refcor(nmode), refshell(nmode) ! complex refracive indices
2531       complex :: crefin(nmode)                  ! complex refractive index
2532
2533 ! FSB special values for EC CORE-shell calculation
2534       real :: DGNSHELL(nmode)    ! modal geometric mean diameter [m]
2535       real :: DGNCORE (nmode)    ! modal geometric mean diameter [m]
2536
2537 ! FSB Modal volumes [ m**3 / m**3 ]
2538       real :: MVOL_ws(nmode)   ! water soluble
2539       real :: MVOL_in(nmode)   ! insolugle
2540       real :: MVOL_ec(nmode)   ! soot like
2541       real :: MVOL_ss(nmode)   !sea salt
2542       real :: MVOL_h2o(nmode)  ! water
2543 !     real :: VOL(nmode)       ! total modal volume [m** 3 / m**3]
2544 ! FSB special values for EC CORE-shell calculation
2545       real :: VOLCOR(nmode)    ! volume of EC core  [m** 3 / m**3]
2546       real :: VOLSHELL(nmode)  ! volume of shell    [m** 3 / m**3]
2547
2548       integer :: m             ! loop index
2549       real    :: bext          ! extinction coefficient [1 / m]
2550       real    :: bscat         ! scattering coefficient [1 / m]
2551       real    :: gfac          ! asymmetry factor
2552
2553       real    :: bextsum, bscatsum, bsgsum
2554
2555 ! FSB History variables by wavelength and mode
2556 !     real    :: bext_wm(ns,nmode)
2557 !     real    :: bscat_wm(ns,nmode)
2558 !     real    :: gfac_wm(ns,nmode)
2559
2560       real, parameter ::  one3rd = 1.0 / 3.0
2561       real  :: dfac      ! ratio of (volcor/vol) ** one3rd
2562                          ! used for calculating the diameter
2563                          ! of the EC core
2564
2565       logical :: succesS
2566
2567 !...component densities [ g/ cm**3 ] <<<<< cgs
2568
2569       real, parameter :: rhows = 1.8   !  bulk density of water soluble aerosol
2570
2571       real, parameter :: rhoin = 2.2   ! bulk density forinsoluble aerosol
2572
2573 !     real, parameter :: rhoec = 1.7   ! bulk density for soot aerosol
2574       real, parameter :: rhoec = 1.8   ! new value
2575
2576       real, parameter :: rhoh2o = 1.0  !  bulk density of aerosol water
2577
2578       real, parameter :: rhoss = 2.2   ! bulk density of seasalt
2579
2580 ! FSB scale factor for volume calculation
2581 !            1.0d-12 * [ cm**3 / g] -> [ m** 3 / ug ]
2582       real, parameter :: scalefactor = 1.0e-12
2583
2584 ! FSB scale factor for [1/g] to [1/ug]
2585       real, parameter ::  cug2g = 1.0e-06
2586
2587 ! FSB reciprocal component densities[ m ** 3 / ug ]
2588
2589       real, parameter :: rhows1 = scalefactor / rhows    !   water soluble aerosol
2590
2591       real, parameter :: rhoin1 = scalefactor /  rhoin   ! insoluble aerosol
2592
2593       real, parameter :: rhoec1 = scalefactor / rhoec    !  soot aerosol
2594
2595       real, parameter :: rhoh2o1 = scalefactor / rhoh2o  !  aerosol water
2596
2597       real, parameter :: rhoss1 = scalefactor / rhoss    !  seasalt
2598
2599       integer,parameter ::  nspint_sw = 14 ! number of spectral intervals for RRTMG_SW
2600
2601 ! FSB Band numbers and wavelengths for RRTMG_SW
2602       integer, parameter :: Band(nspint_sw) = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 /)
2603
2604       real, parameter :: LAMDA_SW(nspint_sw) = (/ 3.4615,  2.7885, 2.325,  2.046,  1.784,         &
2605                                                   1.4625,  1.2705, 1.0101, 0.7016, 0.53325,       &
2606                                                   0.38815, 0.299,  0.2316, 8.24 /) ! wavelength [ um ]
2607
2608 ! *** refractive indices
2609
2610 ! *** Except as otherwise noted reference values of refractive
2611 !     indices  for aerosol particles are from the OPAC Data base.
2612 !     Hess,  Koepke, and  Schult, Optical properties of
2613 !     aerosols and clouds: The software package OPAC, Bulletan of
2614 !     the American Meteorological Society, Vol 79, No 5,
2615 !     pp 831 - 844, May 1998.
2616 !     OPAC is a downloadable data set of optical properties of
2617 !     10 aerosol components, 6 water clouds and 3 cirrus clouds
2618 !     at UV, visible and IR wavelengths
2619 !     www.lrz-muenchen.de/~uh234an/www/radaer/opac.htm
2620
2621
2622 ! FSB water soluble
2623       real, parameter :: xnreal_ws(nspint_sw) = (/ 1.443, 1.420, 1.420, 1.420, 1.463, 1.510, 1.510,   &
2624                                                    1.520, 1.530, 1.530, 1.530, 1.530, 1.530, 1.710 /)
2625       real, parameter :: xnimag_ws(nspint_sw) = (/ 5.718E-3, 1.777E-2, 1.060E-2, 8.368E-3, 1.621E-2,  &
2626                                                    2.198E-2, 1.929E-2, 1.564E-2, 7.000E-3, 5.666E-3,  &
2627                                                    5.000E-3, 8.440E-3, 3.000E-2, 1.100E-1 /)
2628
2629 ! FSB sea salt
2630       real, parameter :: xnreal_ss(nspint_sw) = (/ 1.480, 1.534, 1.437, 1.448, 1.450, 1.462, 1.469,   &
2631                                                    1.470, 1.490, 1.500, 1.502, 1.510, 1.510, 1.510 /)
2632
2633       real, parameter :: xnimag_ss(nspint_sw) = (/ 1.758E-3, 7.462E-3, 2.950E-3, 1.276E-3, 7.944E-4,  &
2634                                                    5.382E-4, 3.754E-4, 1.498E-4, 2.050E-7, 1.184E-8,  &
2635                                                    9.938E-8, 2.060E-6, 5.000E-6, 1.000E-2 /)
2636
2637 ! FSB insoluble
2638       real, parameter :: xnreal_in(nspint_sw) = (/ 1.272, 1.168, 1.208, 1.253, 1.329, 1.418, 1.456,   &
2639                                                    1.518, 1.530, 1.530, 1.530, 1.530, 1.530, 1.470 /)
2640       real, parameter :: xnimag_in(nspint_sw) = (/ 1.165E-2, 1.073E-2, 8.650E-3, 8.092E-3, 8.000E-3,  &
2641                                                    8.000E-3, 8.000E-3, 8.000E-3, 8.000E-3, 8.000E-3,  &
2642                                                    8.000E-3, 8.440E-3, 3.000E-2,9.000E-2 /)
2643
2644 ! FSB 02/11/2012 These values are replaced.
2645 !      data xnreal_ec /1.877, 1.832, 1.813, 1.802, 1.791, 1.769, 1.761,   &
2646 !                      1.760, 1.750, 1.740, 1.750, 1.738, 1.620, 2.120/
2647 !      data xnimag_ec/ 5.563E-1, 5.273E-1, 5.030E-1, 4.918E-1, 4.814E-1,  &
2648 !                      4.585E-1, 4.508E-1, 4.404E-1, 4.300E-1, 4.400E-1,  &
2649 !                      4.600E-1, 4.696E-1, 4.500E-1, 5.700E-1/
2650
2651 !    New  Refractive indices for EC at RRTMG Wavelengths
2652 !       Source  lamda   xnreal_ec       xnimag_ec
2653 !       C&C Values
2654 !          3.4615       2.089   1.070
2655 !               2.7885  2.014   0.939
2656 !               2.325   1.962   0.843
2657 !               2.046   1.950   0.784
2658 !       Bond values
2659 !           1.784       1.940   0.760
2660 !               1.4625  1.930   0.749
2661 !                   1.2705      1.905   0.737
2662 !                   1.0101      1.870   0.726
2663 !  B&B  Values
2664 !          0.7016       1.85    0.71
2665 !          0.53325  1.85    0.71
2666 !                  0.38815  1.85        0.71
2667 !                  0.299        1.85    0.71
2668 !                  0.2316       1.85    0.71
2669 !        C & C values
2670 !          8.24     2.589       1.771
2671 !References:
2672 ! Bond  Personal Communication from Tami Bond
2673 ! B&B Bond, T.C. & R.W. Bergstrom (2006) Light absorption by
2674 ! Carbonaceous Particles: An investigative review,
2675 ! Aerosol Science and Technology. Vol. 40. pp 27-67
2676 !
2677 ! C&C Chang,H and T.T. Charalmpopoulos (1990) Determination of the
2678 ! wavelength dependence of refractive indices of flame soot,
2679 ! Proceeding of the Royal Society of London A, Vol. 430, pp 577-591.
2680 ! FSB new values
2681
2682 ! FSB elemental carbon - soot like
2683
2684       real, parameter :: xnreal_ec(nspint_sw) = (/ 2.089, 2.014, 1.962, 1.950, 1.940, 1.930, 1.905,   &
2685                                                    1.870, 1.85,  1.85,  1.85,  1.85,  1.85,  2.589 /)
2686       real, parameter :: xnimag_ec(nspint_sw) = (/ 1.070, 0.939, 0.843, 0.784, 0.760, 0.749, 0.737,   &
2687                                                    0.726, 0.71,  0.71,  0.71,  0.71,  0.71,  1.771 /)
2688
2689 ! FSB water
2690       real, save :: xnreal_h2o(nspint_sw) = (/ 1.408, 1.324, 1.277, 1.302, 1.312, 1.321, 1.323,       &
2691                                                1.327, 1.331, 1.334, 1.340, 1.349, 1.362, 1.260 /)
2692       real, save :: xnimag_h2o(nspint_sw) = (/ 1.420E-2, 1.577E-1, 1.516E-3, 1.159E-3, 2.360E-4,      &
2693                                                1.713E-4, 2.425E-5, 3.125E-6, 3.405E-8, 1.639E-9,      &
2694                                                2.955E-9, 1.635E-8, 3.350E-8, 6.220E-2 /)
2695
2696 ! FSB Begin code ======================================================
2697
2698       bextsum  = 0.0
2699       bscatsum = 0.0
2700       bsgsum   = 0.0
2701       do m = 1, nmode
2702 ! FSB calculate volumes [ m**3 / m**3 ]
2703 ! FSB the reciprocal densities have been scaled to [ m**3 / ug ]
2704
2705          MVOL_ws(m)  = rhows1  * INMASS_ws(m)
2706          MVOL_in(m)  = rhoin1  * INMASS_in(m)
2707          MVOL_ec(m)  = rhoec1  * INMASS_ec(m)
2708          MVOL_ss(m)  = rhoss1  * INMASS_ss(m)
2709          MVOL_h2o(m) = rhoh2o1 * INMASS_h2o(m)
2710
2711          VOLSHELL(m) = MVOL_ws(m) + MVOL_in(m) + MVOL_ss(m) + MVOL_h2o(m)
2712          VOLCOR(m)   = MVOL_ec(m)
2713 !        VOL(m)      = VOLSHELL(m) + VOLCOR(m) ! VOL is total volume
2714
2715          if ( VOLCOR(m) .gt. 0.0 ) then
2716 ! FSB EC is present
2717 !       calculate the ratio of core to shell volume
2718 !       take cube root for scaling the diameter of
2719 !       the core to that of the shell.
2720
2721 !           dfac        = ( VOLCOR(m) / VOL(m) ) ** one3rd
2722             dfac        = ( VOLCOR(m) / ( VOLSHELL(m) + VOLCOR(m) ) ) ** one3rd
2723 !           dfac        = ( VOLCOR(m) / ( VOLSHELL(m) + VOLCOR(m) ) )
2724 ! FSB Set shell and core diameters
2725             DGNSHELL(m) = INDGN(m)
2726             DGNCORE(M)  = dfac * INDGN(m)
2727 ! FSB note that VOLSHELL(m) is the total volume when EC is not present
2728          end if
2729
2730 ! internal mixture of non-EC species.
2731
2732 ! modal real refractive index No EC
2733          nr(m) = (MVOL_ws(m)  * xnreal_ws(ns) +              &
2734                   MVOL_in(m)  * xnreal_in(ns) +              &
2735                   MVOL_ss(m)  * xnreal_ss(ns) +              &
2736 !                 MVOL_h2o(m) * xnreal_h2o(ns)) / VOL(m)
2737                   MVOL_h2o(m) * xnreal_h2o(ns)) / VOLSHELL(m)
2738
2739 ! modal imaginary refractive index no EC
2740          ni(m) = (MVOL_ws(m)  * xnimag_ws(ns) +              &
2741                   MVOL_in(m)  * xnimag_in(ns) +              &
2742                   MVOL_ss(m)  * xnimag_ss(ns) +              &
2743 !                 MVOL_h2o(m) * xnimag_h2o(ns)) / VOL(m)
2744                   MVOL_h2o(m) * xnimag_h2o(ns)) / VOLSHELL(m)
2745
2746          if ( VOLCOR(m) .gt. 0.0) then
2747
2748 ! FSB calculate the complex refractive indices for the CORE and
2749 !     the SHELL for case when and EC core is assumed to exist.
2750
2751             refcor(m)   = cmplx( xnreal_ec(ns), xnimag_ec(ns) )
2752             refshell(m) = cmplx(nr(m), ni(m) )
2753 ! FSB do BHCOAT case
2754             CALL aero_optical_CS( LAMDA_SW(ns), refcor(m), refshell(m),   &
2755                                   VOLCOR(m),VOLSHELL(m), DGNCORE(m),      &
2756                                   DGNSHELL(m), INSIG(m),                  &
2757                                   bext, bscat, gfac, succesS )
2758 !        else if ( VOL(m) .gt. 0.0) then
2759          else if ( VOLSHELL(m) .gt. 0.0) then
2760 ! FSB do BHMIE case for the case when EC is not present.
2761             crefin(m) = cmplx(nr(m), ni(m) )
2762 !           CALL aero_optical2( LAMDA_SW(ns), crefin(m), VOL(m),          &
2763             CALL aero_optical2( LAMDA_SW(ns), crefin(m), VOLSHELL(m),     &
2764                                 INDGN(m), INSIG(m),                       &
2765                                 bext, bscat, gfac, success )
2766          else
2767             bext = 0.0
2768             bscat = 0.0
2769             gfac = 0.0
2770          end if
2771
2772 ! FSB sum for total values
2773          bextsum  = bextsum + bext
2774          bscatsum = bscatsum +bscat
2775          bsgsum   = bsgsum + bscat * gfac
2776 ! FSB get history
2777 !        bext_wm(ns,m)  = bext
2778 !        bscat_wm(ns,m) = bscat
2779 !        gfac_wm(ns,m)  = gfac
2780       end do ! loop over modes
2781
2782 ! FSB construct output variables
2783       tauaer = bextsum * delta_z
2784       waer   = bscatsum / bextsum
2785       gaer   = bsgsum / bscatsum
2786
2787    end subroutine get_aerosol_Optics_RRTMG_SW
2788
2789 END MODULE module_ra_rrtmg_sw_cmaq