!$Id: RTS_mie_sourcecode.f90,v 1.1 2010/07/30 23:47:04 daven Exp $ SUBROUTINE Mie_main & ( max_Mie_angles, max_Mie_sizes, & ! D max_Mie_points, max_Mie_distpoints, & ! D do_external_angles, do_coeffct_angles, do_use_cutoff, & ! I idis, nr_parameters, startup, & ! I nblocks, nweights, cutoff, & ! I n_external_angles, external_angle_cosines, & ! I n_coeffct_angles, coeff_cosines, coeff_weights, & ! I m_complex, xparticle_limit, wavelength, rmax, rmin, & ! I mie_bulk, dist, fmat, & ! O message, trace, action, failmie ) ! O ! modules USE Mie_precision USE MIE_constants, ONLY : d_zero, d_half, d_one, d_two, d_three, c_i ! implicit none statement IMPLICIT NONE ! Dimensioning input INTEGER, INTENT (IN) :: max_Mie_angles, max_Mie_sizes INTEGER, INTENT (IN) :: max_Mie_points, max_Mie_distpoints ! input LOGICAL , INTENT (IN) :: do_external_angles LOGICAL , INTENT (IN) :: do_coeffct_angles LOGICAL , INTENT (IN) :: do_use_cutoff INTEGER , INTENT (IN) :: idis REAL (KIND=dp), INTENT (IN) :: nr_parameters(3) INTEGER , INTENT (IN) :: nblocks INTEGER , INTENT (IN) :: nweights REAL (KIND=dp), INTENT (IN) :: cutoff INTEGER , INTENT (IN) :: n_external_angles REAL (KIND=dp), INTENT (IN) :: external_angle_cosines(max_Mie_angles) LOGICAL , INTENT (INOUT) :: startup INTEGER , INTENT (INOUT) :: n_coeffct_angles REAL (KIND=dp), INTENT (INOUT) :: coeff_cosines(max_Mie_angles) REAL (KIND=dp), INTENT (INOUT) :: coeff_weights(max_Mie_angles) COMPLEX (KIND=dp), INTENT (IN) :: m_complex REAL (KIND=dp), INTENT (IN) :: xparticle_limit REAL (KIND=dp), INTENT (IN) :: wavelength ! output REAL (KIND=dp), INTENT (OUT) :: MIE_BULK(4) REAL (KIND=dp), INTENT (OUT) :: DIST(5) REAL (KIND=dp), INTENT (OUT) :: FMAT(4,max_Mie_angles) LOGICAL , INTENT (OUT) :: failmie CHARACTER*(*) , INTENT (OUT) :: message, trace, action REAL (KIND=dp), INTENT (INOUT) :: rmin, rmax ! Local Mie output REAL (KIND=dp), DIMENSION (max_Mie_sizes) :: q_ext, q_sca, asym COMPLEX (KIND=dp), DIMENSION (max_Mie_angles, max_Mie_sizes) :: splus,sminus ! local variables for Mie code CHARACTER*5 :: char5 LOGICAL :: do_angular_variation LOGICAL :: failmm, faild INTEGER :: i, angle, n_angles INTEGER :: iblock, n_sizes, kf REAL (KIND=dp) :: factor_0, factor_1, d_pi REAL (KIND=dp) :: rstart, rfinis, help, rblock REAL (KIND=dp) :: quad, quadr2, quadr3, quadr4 REAL (KIND=dp) :: ndens, gxsec, reff, volume, veff REAL (KIND=dp) :: Qext, Qsca, Qasy, ssalbedo REAL (KIND=dp) :: f(4) REAL (KIND=dp) :: angle_cosines(max_Mie_angles) ! redundant variables ! REAL (KIND=dp) :: xeff, xeff_d(3) ! LOGICAL :: fail COMPLEX (KIND=dp) :: sp, sm, csp, csm COMPLEX (KIND=dp) :: c_mi REAL (KIND=dp) :: xpart_root3, xparticle COMPLEX (KIND=dp) :: y_argument INTEGER :: limmax, limstop, limsize REAL (KIND=dp), DIMENSION (max_Mie_sizes) :: particle_sizes REAL (KIND=dp), DIMENSION (max_Mie_sizes) :: rquad, weights, nr ! start up ! -------- c_mi = - c_i n_sizes = nweights d_pi = 4.0_dp * ATAN(d_one) factor_0 = d_two * d_pi / wavelength factor_1 = wavelength / factor_0 ! Zeroing ! ------- trace = ' ' message = ' ' action = ' ' failmie = .FALSE. failmm = .false. Qext = d_zero Qsca = d_zero Qasy = d_zero ndens = d_zero gxsec = d_zero reff = d_zero veff = d_zero ! limiting radii calculation if ( do_use_cutoff ) then CALL rminmax ( idis, nr_parameters, cutoff, rmin, rmax, message, failmm ) IF ( failmm ) THEN failmie = .TRUE. trace = 'Trace : First Check in Mie Main. Failed to find radii extrema' action = 'Action : Consult with R. Spurr' RETURN END IF endif ! Check limiting radii if set externally if ( .not. do_use_cutoff ) then if ( rmin.lt.0.0d0 ) then failmm = .true. message = 'External Rmin value < 0, out of bounds' else if ( rmax .le. 0.0d0 ) then failmm = .true. message = 'External Rmax value =< 0, out of bounds' else if ( rmin .ge. rmax ) then failmm = .true. message = 'External Rmin >= Rmax, Cannot be possible!' endif if ( failmm ) then trace = 'Trace : First Check in Mie Main. User Rmin/Rmax wrong' action = 'Action : Change input values of Rmin and Rmax' RETURN END IF endif ! number of blocks rblock = ( rmax - rmin ) / DBLE(nblocks) ! limiting number of terms for coefficient computation xparticle = factor_0 * rmax y_argument = xparticle * m_complex limstop = 2 IF ( xparticle > 0.02) THEN xpart_root3 = xparticle ** ( d_one / d_three ) IF ( xparticle <= 8.0_dp ) THEN limstop = xparticle + 4.0_dp * xpart_root3 + d_two ELSE IF ( xparticle < 4200.0_dp ) THEN limstop = xparticle + 4.05_dp * xpart_root3 + d_two ELSE limstop = xparticle + 4.0_dp * xpart_root3 + d_two END IF END IF limmax = nint(max(DBLE(limstop),ABS(y_argument)) + 15.0_dp) ! Dimensioning and exception handling checks ! ------------------------------------------ ! return if size limit exceeded IF ( xparticle > xparticle_limit ) THEN failmie = .TRUE. limsize = int(xparticle) + 1 write(char5,'(i5)')limsize message = 'Message : error size parameter overflow' trace = 'Trace : Second check in Mie_main' action = 'Action : In configuration file, increase cutoff or '// & 'increase xparticle_limit to at least '//char5 RETURN END IF ! return if maximum number of terms too great IF ( limstop > max_Mie_points ) THEN failmie = .TRUE. write(char5,'(i5)')limstop message = 'Message : Insufficient dimensioning for maximum number of terms' trace = 'Trace : Third check in Mie_main' action = 'Action : Increase max_Mie_points in calling program to at least '//char5 RETURN END IF ! And again, Dave recurrence IF ( limmax > max_Mie_points ) THEN failmie = .TRUE. write(char5,'(i5)')limmax message = 'Message : Insufficient dimensioning for maximum number of terms (Dave recurrence)' trace = 'Trace : Fourth check in Mie_main' action = 'Action : Increase max_Mie_points in calling program to at least '//char5 RETURN END IF ! Compute the number of angles required for coefficient computation IF ( do_coeffct_angles .AND. startup ) THEN n_coeffct_angles = 2*limmax + 2 if ( n_coeffct_angles > max_Mie_angles ) then failmie = .true. write(char5,'(i5)')n_coeffct_angles message = 'Message : Dimensioning error for number of terms for coefficient computation' trace = 'Trace : Fifth check in Mie Main' action = 'Action : Increase value of max_Mie_angles in calling program to at least '//char5 return endif ENDIF ! Compute the angles required for coefficient computation ! ------------------------------------------------------- ! Quadrature (only need to do it once) IF ( do_coeffct_angles .AND. startup ) THEN n_coeffct_angles = 2*limmax + 2 CALL mie_gauleg ( max_Mie_angles, n_coeffct_angles, -1.0_dp, 1.0_dp, & ! Input coeff_cosines, coeff_weights ) ! Output startup = .FALSE. END IF ! Overall cosines IF ( do_coeffct_angles ) THEN n_angles = n_coeffct_angles DO angle = 1, n_angles angle_cosines(angle) = coeff_cosines(angle) END DO do_angular_variation = .TRUE. ELSE IF ( do_external_angles ) THEN n_angles = n_external_angles DO angle = 1, n_angles angle_cosines(angle) = external_angle_cosines(angle) END DO do_angular_variation = .TRUE. ELSE n_angles = 0 do_angular_variation = .FALSE. END IF ! zero the angular input, if flagged IF ( do_angular_variation ) THEN DO angle = 1, n_angles DO kf = 1, 4 fmat(kf,angle) = d_zero END DO END DO END IF ! start integration ! ----------------- DO iblock = 1, nblocks rstart = rmin + ( iblock-1) * rblock rfinis = rstart + rblock CALL mie_gauleg ( max_Mie_sizes, n_sizes, rstart, rfinis, rquad, weights ) ! prepare particle sizes DO i = 1, n_sizes particle_sizes(i) = factor_0 * rquad(i) ENDDO ! Call to coefficients ! WARNING - easy o get segmentation fault before this call ! Check dimensioning first, Use lots of memory !!!!!!! CALL mie_coeffs & ( max_Mie_angles, max_Mie_sizes, max_Mie_points, & ! Dimensioning do_angular_variation, & ! Input n_angles, n_sizes, m_complex, & ! Input particle_sizes, angle_cosines, & ! Input q_ext, q_sca, asym, splus, sminus ) ! Output ! size distribution and derivatives CALL sizedis & ( max_Mie_sizes, idis, nr_parameters, rquad, n_sizes, & nr, message, faild ) IF ( faild ) THEN failmie = faild write(char5,'(i5)')iblock trace = 'Trace : Sixth check in Mie_main. Subroutine sizedis failed for block number '//char5 action = 'Action : Consult with R. Spurr' RETURN END IF ! Integration over particle sizes within block ! -------------------------------------------- DO i = 1, n_sizes ! Number density, geometric cross-section, 3rd and 4th powers quad = nr(i) * weights(i) quadr2 = quad * rquad(i) * rquad(i) quadr3 = quadr2 * rquad(i) quadr4 = quadr3 * rquad(i) ndens = ndens + quad gxsec = gxsec + quadr2 reff = reff + quadr3 veff = veff + quadr4 ! Basic coefficients Qext = Qext + quad * q_ext(i) Qsca = Qsca + quad * q_sca(i) Qasy = Qasy + quad * asym(i) ! angular variation loop IF ( do_angular_variation ) THEN DO angle = 1, n_angles sp = splus(angle,i) sm = sminus(angle,i) csp = CONJG(sp) csm = CONJG(sm) f(1) = REAL ( sp * csp + sm * csm ) f(2) = - REAL ( sm * csp + sp * csm ) f(3) = REAL ( sp * csp - sm * csm ) f(4) = REAL ( ( sm * csp - sp * csm ) * c_mi ) DO kf = 1, 4 FMAT(kf,angle) = FMAT(kf,angle) + quad*f(kf) ENDDO END DO END IF ! Finish integration loops END DO END DO ! Final Assignations ! ------------------ ! F matrix stuff IF ( do_angular_variation ) THEN DO angle = 1, n_angles DO kf = 1, 4 FMAT(kf,angle) = d_half * FMAT(kf,angle) / Qsca END DO END DO END IF ! geometric cross-section gxsec = d_pi * gxsec ! asymmetry parameter Qasy = d_two * Qasy / Qsca ! basic coefficients Qsca = Qsca * factor_1 Qext = Qext * factor_1 Qsca = Qsca/gxsec Qext = Qext/gxsec ! single scattering albedo ssalbedo = Qsca/Qext ! geometrical quantities volume= (4.0_dp/3.0_dp) * d_pi * reff reff = d_pi * reff / gxsec ! Variance output help = d_pi / gxsec / reff / reff veff = help * veff veff = veff - d_one ! Particle size parameter output ! xeff = factor_0 * reff ! Final assignation MIE_BULK(1) = Qext MIE_BULK(2) = Qsca MIE_BULK(3) = Qasy MIE_BULK(4) = ssalbedo DIST(1) = ndens DIST(2) = gxsec DIST(3) = volume DIST(4) = reff ! DIST(5) = xeff DIST(5) = veff RETURN END SUBROUTINE Mie_main SUBROUTINE mie_coeffs & ! ( max_Mie_angles, max_Mie_sizes, max_Mie_points, & ! Dimensioning do_angular_variation, n_angles, n_sizes, m_complex, & ! Input particle_sizes, angle_cosines, & ! Input q_ext, q_sca, asym, splus, sminus ) ! Output ! name: ! mie_coeffs ! purpose: ! calculates the scattering parameters of a series of particles ! using the mie scattering theory. FOR USE WITH POLYDISPERSE COD ! inputs: ! particle_sizes: array of particle size parameters ! angle_cosines: array of angle cosines ! m_complex: the complex refractive index of the particles ! n_angles, n_sizes: number of scattering angles, number of particle sizes ! do_angular_variation: flag for S+/S- output ! outputs (1): ! q_ext: the extinction efficiency ! q_sca: the scattering efficiency ! asym: the asymmetry parameter ! splus: the first amplitude function ! sminus: the second amplitude function ! modification history ! g. thomas IDL Mie code (February 2004). Basic Monodisperse derivatives. ! r. spurr F90 Mie code ( October 2004). Extension all Polydisperse derivatives. ! r. spurr Exception handling (September 2008). Exception handling removed. ! modules USE Mie_precision USE MIE_constants, ONLY : d_zero, d_half, d_one, d_two, d_three, & c_zero, c_one, c_i ! implicit none statement IMPLICIT NONE ! Dimensioning input INTEGER, INTENT (IN) :: max_Mie_angles, max_Mie_sizes INTEGER, INTENT (IN) :: max_Mie_points ! input LOGICAL , INTENT (IN) :: do_angular_variation INTEGER , INTENT (IN) :: n_angles, n_sizes COMPLEX (KIND=dp), INTENT (IN) :: m_complex REAL (KIND=dp), DIMENSION (max_Mie_sizes), INTENT (IN) :: particle_sizes REAL (KIND=dp), DIMENSION (max_Mie_angles), INTENT (IN) :: angle_cosines ! output (1) REAL (KIND=dp), DIMENSION (max_Mie_sizes), INTENT (OUT) :: q_ext, q_sca, asym COMPLEX (KIND=dp), DIMENSION (max_Mie_angles, max_Mie_sizes), INTENT (OUT) :: splus, sminus ! local variables for Mie code INTEGER :: size, angle, n, nm1, nstop(max_Mie_sizes), nmax, maxstop REAL (KIND=dp) :: xparticle, xpart_root3 REAL (KIND=dp) :: xinv, xinvsq, two_d_xsq, xinv_dx REAL (KIND=dp) :: dn, dnp1, dnm1, dnsq, dnnp1, tnp1, tnm1, hnp1, hnm1 REAL (KIND=dp) :: cos_x, sin_x, psi0, psi1, chi0, chi1, psi, chi REAL (KIND=dp) :: s, t, tau_n, factor, forward, bckward COMPLEX (KIND=dp) :: inverse_m, y_argument, yinv, yinvsq, a1, zeta, zeta1 COMPLEX (KIND=dp) :: an, bn, an_star, bn_star, anm1, bnm1, bnm1_star COMPLEX (KIND=dp) :: biga_divs_m, biga_mult_m, noverx, aterm, bterm COMPLEX (KIND=dp) :: facplus, facminus, c_mi COMPLEX (KIND=dp) :: an_denom, bn_denom ! redundant variables ! REAL (KIND=dp), INTENT (IN) :: xparticle_limit ! subroutine argument ! REAL (KIND=dp) :: four_d_xsq ! CHARACTER*4 :: char4 ! COMPLEX (KIND=dp) :: common, an_denom_dm, s1, s2 ! INTEGER :: nmax_end ! local arrays REAL (KIND=dp), DIMENSION (max_Mie_angles) :: pi_n, pi_nm1 COMPLEX (KIND=dp), DIMENSION (max_Mie_points) :: biga REAL (KIND=dp), DIMENSION (max_Mie_angles,max_Mie_points) :: polyplus, polyminus ! Initial section ! --------------- maxstop = 0 c_mi = -c_i DO size = 1, n_sizes ! particle size xparticle = particle_sizes (size) ! assign number of terms and maximum IF ( xparticle < 0.02) THEN nstop(size) = 2 ELSE xpart_root3 = xparticle ** ( d_one / d_three ) IF ( xparticle <= 8.0_dp ) THEN nstop(size) = xparticle + 4.0_dp * xpart_root3 + d_two ELSE IF ( xparticle < 4200.0_dp ) THEN nstop(size) = xparticle + 4.05_dp * xpart_root3 + d_two ELSE nstop(size) = xparticle + 4.0_dp * xpart_root3 + d_two END IF END IF maxstop = max(nstop(size),maxstop) END DO ! phase function expansion polynomials ! ---> initialise phase function Legendre polynomials ! ---> Recurrence phase function Legendre polynomials IF ( do_angular_variation ) THEN DO angle = 1, n_angles pi_nm1(angle) = d_zero pi_n(angle) = d_one END DO DO n = 1, maxstop nm1 = n - 1 dn = dble(n) dnp1 = dn + d_one forward = dnp1 / dn DO angle = 1, n_angles s = angle_cosines(angle) * pi_n(angle) t = s - pi_nm1(angle) tau_n = dn*t - pi_nm1(angle) polyplus(angle,n) = pi_n(angle) + tau_n polyminus(angle,n) = pi_n(angle) - tau_n pi_nm1(angle) = pi_n(angle) pi_n(angle) = s + t*forward END DO END DO END IF ! start loop over particle sizes ! ------------------------------ DO size = 1, n_sizes ! initialize output asym(size) = d_zero q_ext(size) = d_zero q_sca(size) = d_zero ! some auxiliary quantities xparticle = particle_sizes (size) xinv = d_one / xparticle xinvsq = xinv * xinv two_d_xsq = d_two * xinvsq xinv_dx = - d_two * xinv inverse_m = c_one / m_complex y_argument = xparticle * m_complex yinv = d_one / y_argument yinvsq = yinv * yinv ! Biga = ratio derivative, recurrence due to J. Dave nmax = nint(max(dble(nstop(size)),abs(y_argument)) + 15.0_dp) biga(nmax) = c_zero DO n = nmax-1, 1,-1 a1 = dble(n+1) / y_argument biga(n) = a1 - c_one / (a1+biga(n+1)) END DO ! initialize Riccati-Bessel functions tnp1 = d_one cos_x = COS(xparticle) sin_x = SIN(xparticle) psi0 = cos_x psi1 = sin_x chi1 =-cos_x chi0 = sin_x zeta1 = CMPLX(psi1,chi1,kind=dp) ! initialise sp and sm IF ( do_angular_variation ) THEN DO angle = 1, n_angles splus(angle,size) = c_zero sminus(angle,size) = c_zero END DO END IF ! main loop DO n = 1, nstop(size) ! various factors dn = dble(n) dnp1 = dn + d_one dnm1 = dn - d_one tnp1 = tnp1 + d_two tnm1 = tnp1 - d_two dnsq = dn * dn dnnp1 = dnsq + dn factor = tnp1 / dnnp1 bckward = dnm1 / dn ! Ricatti - Bessel recurrence psi = tnm1 * psi1/xparticle - psi0 chi = tnm1 * chi1/xparticle - chi0 zeta = CMPLX(psi,chi,kind=dp) ! a(n) and b(n) biga_divs_m = biga(n) * inverse_m biga_mult_m = biga(n) * m_complex noverx = CMPLX(dn/xparticle,d_zero,kind=dp) aterm = biga_divs_m + noverx bterm = biga_mult_m + noverx an_denom = (aterm * zeta - zeta1) bn_denom = (bterm * zeta - zeta1) an = ( aterm*psi-psi1 ) / an_denom bn = ( bterm*psi-psi1 ) / bn_denom an_star = CONJG(an) bn_star = CONJG(bn) ! basic coefficients ! ------------------ ! Q coefficients q_ext(size) = q_ext(size) + tnp1 * REAL ( an + bn ) q_sca(size) = q_sca(size) + tnp1 * REAL ( an*CONJG(an) + bn*CONJG(bn) ) ! asymmetry parameter IF ( n > 1 ) THEN hnp1 = bckward * dnp1 hnm1 = tnm1 / (dnsq - dn) asym(size) = asym(size) & + hnp1 * REAL ( anm1*an_star + bnm1*bn_star) & + hnm1 * REAL ( anm1*bnm1_star) END IF ! Upgrades ! -------- ! upgrade an/bn recurrences (only for asymmetry parameter) anm1 = an bnm1 = bn bnm1_star = bn_star ! upgrade Ricatti-Bessel recurrences psi0 = psi1 psi1 = psi chi0 = chi1 chi1 = chi zeta1 = CMPLX(psi1,chi1,kind=dp) ! S+/S- function stuff ! -------------------- IF ( do_angular_variation ) THEN facplus = factor * ( an + bn ) facminus = factor * ( an - bn ) DO angle = 1, n_angles splus(angle,size) = splus(angle,size) + facplus * polyplus(angle,n) sminus(angle,size) = sminus(angle,size) + facminus * polyminus(angle,n) END DO END IF ! end sum loop END DO ! End loop and finish ! ------------------- ! end loop over particle sizes END DO ! finish RETURN END SUBROUTINE mie_coeffs ! Contains the following modules ! sizedist ! sizedist_nod ! gammafunction ! gauleg ! rminmax SUBROUTINE sizedis & ( max_Mie_distpoints, idis, par, radius, numradius, & nwithr, message, faild ) !************************************************************************ !* Calculate the size distribution n(r) for the numr radius values * !* contained in array r and return the results through the array nwithr* !* The size distributions are normalized such that the integral over * !* all r is equal to one. * !************************************************************************ ! modules USE Mie_precision USE MIE_constants, ONLY : d_zero, d_half, d_one, d_two, d_three IMPLICIT NONE !* subroutine arguments INTEGER , INTENT (IN) :: max_Mie_distpoints REAL (KIND=dp), INTENT (IN) :: par(3) INTEGER , INTENT (IN) :: idis, numradius CHARACTER*(*) , INTENT (OUT) :: message LOGICAL , INTENT (OUT) :: faild REAL (KIND=dp), DIMENSION (max_Mie_distpoints), INTENT (IN) :: radius REAL (KIND=dp), DIMENSION (max_Mie_distpoints), INTENT (OUT) :: nwithr !* local variables INTEGER :: i REAL (KIND=dp) :: pi,r,logr,root2p REAL (KIND=dp) :: alpha,alpha1,b,logb,arg1,arg2,arg,argsq,r3 REAL (KIND=dp) :: b1,b2,logb1,logb2,rc REAL (KIND=dp) :: logrg,logsi,logsi_inv,gamma,gamma1,rg REAL (KIND=dp) :: rmin,rmax,fac1,fac2,aperg, alpha2 REAL (KIND=dp) :: n1, n2, C, logC, logC1, logC2 REAL (KIND=dp) :: gammln, dummy CHARACTER*70 :: message_gamma LOGICAL :: fail character*1 :: cdis ! check faild = .FALSE. if (idis == 0 ) RETURN IF ( IDIS > 8 ) THEN faild = .TRUE. message = 'illegal index in sizedis' RETURN END IF ! setup pi = dacos(-1.d0) root2p = dsqrt(pi+pi) ! IDIS = 1 : TWO-PARAMETER GAMMA with alpha and b given IF ( idis == 1 ) THEN alpha = par(1) b = par(2) alpha1 = alpha + d_one logb = LOG(b) CALL gammafunction ( alpha1, .false., gammln, dummy, fail, message_gamma ) IF ( fail ) go to 240 logC = alpha1*logb - gammln DO i = 1, numradius r = radius(i) logr = LOG(r) arg1 = logC + alpha*logr nwithr(i) = EXP ( arg1 - b*r ) END DO ! IDIS = 2 : TWO-PARAMETER GAMMA with par(1)= reff and par(2)= veff given ELSE IF ( idis == 2 ) THEN alpha = d_one/par(2) - d_three b = d_one/(par(1)*par(2)) alpha1 = alpha + d_one logb = LOG(b) CALL gammafunction ( alpha1, .false., gammln, dummy, fail, message_gamma ) IF ( fail ) go to 240 logC = alpha1*logb - gammln DO i = 1, numradius r = radius(i) logr = LOG(r) arg1 = logC + alpha*logr nwithr(i) = EXP ( arg1 - b*r ) END DO ! IDIS = 3 : BIMODAL GAMMA with equal mode weights ELSE IF ( idis == 3 ) THEN alpha = d_one/par(3) - d_three b1 = d_one/(par(1)*par(3)) b2 = d_one/(par(2)*par(3)) alpha1 = alpha + d_one CALL gammafunction ( alpha1, .false., gammln, dummy, fail, message_gamma ) logb1 = LOG(b1) logb2 = LOG(b2) logC1 = alpha1*logb1 - gammln logC2 = alpha1*logb2 - gammln DO i = 1, numradius r = radius(i) logr = LOG(r) arg1 = logC1 + alpha*logr arg2 = logC2 + alpha*logr n1 = EXP(arg1 - b1*r) n2 = EXP(arg2 - b2*r) nwithr(i) = d_half * ( n1 + n2 ) END DO ! 4 LOG-NORMAL with rg and sigma given ELSE IF ( idis == 4 ) THEN logrg = dlog(par(1)) logsi = dabs(dlog(par(2))) logsi_inv = d_one / logsi C = logsi_inv / root2p DO i = 1, numradius r = radius(i) logr = LOG(r) arg = ( logr - logrg ) / logsi argsq = arg * arg nwithr(i) = C * dexp( - d_half * argsq ) / r END DO ! 5 LOG-NORMAL with reff and veff given * ELSE IF ( idis == 5 ) THEN alpha1 = d_one + par(2) alpha2 = dlog(alpha1) rg = par(1)/(d_one+par(2))**2.5_dp logrg = dlog(rg) logsi = dsqrt(alpha2) logsi_inv = d_one / logsi C = logsi_inv / root2p DO i = 1, numradius r = radius(i) logr = LOG(r) arg = ( logr - logrg ) / logsi argsq = arg * arg nwithr(i) = C * dexp( - d_half * argsq ) / r END DO ! 6 POWER LAW * ELSE IF ( idis == 6 ) THEN alpha = par(1) rmin = par(2) rmax = par(3) alpha1 = alpha - d_one fac1 = (rmax/rmin)**alpha1 fac2 = d_one / ( fac1 - d_one ) C = alpha1 * rmax**alpha1 * fac2 DO i = 1, numradius r = radius(i) if ( (r < rmax) .and. (r > rmin) ) then nwithr(i) = C*r**(-alpha) else nwithr(i) = d_zero endif END DO ! 7 MODIFIED GAMMA with alpha, rc and gamma given ELSE IF ( idis == 7 ) THEN alpha = par(1) rc = par(2) gamma = par(3) b = alpha / (gamma*rc**gamma) logb = LOG(b) alpha1 = alpha + d_one gamma1 = d_one / gamma aperg = alpha1/gamma CALL gammafunction ( aperg, .false., gammln, dummy, fail, message_gamma ) IF ( fail ) go to 240 logC = dlog(gamma) + aperg*logb - gammln DO i = 1, numradius r = radius(i) logr = LOG(r) arg1 = logC + alpha*logr r3 = b*r ** gamma nwithr(i) = EXP ( arg1 - r3 ) END DO ! 8 MODIFIED GAMMA with alpha, b and gamma given ELSE IF ( idis == 8 ) THEN alpha = par(1) b = par(2) gamma = par(3) alpha1 = alpha + d_one gamma1 = d_one / gamma logb = LOG(b) aperg = alpha1/gamma CALL gammafunction ( aperg, .false., gammln, dummy, fail, message_gamma ) IF ( fail ) go to 240 logC = dlog(gamma) + aperg*logb - gammln DO i = 1, numradius r = radius(i) logr = LOG(r) arg1 = logC + alpha*logr r3 = r ** gamma nwithr(i) = EXP ( arg1 - b*r3 ) END DO END IF ! normal return RETURN ! special return 240 CONTINUE faild = .TRUE. write(cdis,'(I1)')idis message = message_gamma(1:LEN(message_gamma))//', distribution : '//cdis RETURN END SUBROUTINE sizedis SUBROUTINE sizedis_nod & ( idis, par, numradius, radius, nwithr, message, failnod ) !************************************************************************ !* Calculate the size distribution n(r) for the numr radius values * !* contained in array r and return the results through the array nwithr* !* The size distributions are normalized such that the integral over * !* all r is equal to one. * !************************************************************************ ! modules USE Mie_precision USE MIE_constants, ONLY : d_zero, d_half, d_one, d_two, d_three IMPLICIT NONE !* subroutine arguments REAL (KIND=dp), INTENT (IN) :: par(3) INTEGER , INTENT (IN) :: idis, numradius CHARACTER*(*) , INTENT (OUT) :: message LOGICAL , INTENT (OUT) :: failnod REAL (KIND=dp), DIMENSION (numradius), INTENT (IN) :: radius REAL (KIND=dp), DIMENSION (numradius), INTENT (OUT) :: nwithr !* local variables LOGICAL :: deriv INTEGER :: i REAL (KIND=dp) :: pi,r,logr,root2p REAL (KIND=dp) :: alpha,alpha1,b,logb,arg1,arg2,arg,argsq,r3 REAL (KIND=dp) :: b1,b2,logb1,logb2,rc, gammln, dgammln REAL (KIND=dp) :: logrg,logsi,logsi_inv,gamma,gamma1,rg REAL (KIND=dp) :: rmin,rmax,fac1,fac2,aperg,alpha2 REAL (KIND=dp) :: n1, n2, C, logC, logC1, logC2 CHARACTER*70 :: message_gamma LOGICAL :: fail character*1 :: cdis ! check failnod = .FALSE. if (idis == 0 ) RETURN IF ( IDIS > 8 ) THEN failnod = .TRUE. message = 'illegal index in sizedis' RETURN END IF ! setup pi = dacos(-1.d0) root2p = dsqrt(pi+pi) deriv = .FALSE. ! IDIS = 1 : TWO-PARAMETER GAMMA with alpha and b given IF ( idis == 1 ) THEN alpha = par(1) b = par(2) alpha1 = alpha + d_one logb = LOG(b) CALL gammafunction ( alpha1, deriv, gammln, dgammln, fail, message_gamma ) IF ( fail ) go to 240 logC = alpha1*logb - gammln DO i = 1, numradius r = radius(i) logr = LOG(r) arg1 = logC + alpha*logr nwithr(i) = EXP ( arg1 - b*r ) END DO ! IDIS = 2 : TWO-PARAMETER GAMMA with par(1)= reff and par(2)= veff given ELSE IF ( idis == 2 ) THEN alpha = d_one/par(2) - d_three b = d_one/(par(1)*par(2)) alpha1 = alpha + d_one logb = LOG(b) CALL gammafunction ( alpha1, deriv, gammln, dgammln, fail, message_gamma ) IF ( fail ) go to 240 logC = alpha1*logb - gammln DO i = 1, numradius r = radius(i) logr = LOG(r) arg1 = logC + alpha*logr nwithr(i) = EXP ( arg1 - b*r ) END DO ! IDIS = 3 : BIMODAL GAMMA with equal mode weights ELSE IF ( idis == 3 ) THEN alpha = d_one/par(3) - d_three b1 = d_one/(par(1)*par(3)) b2 = d_one/(par(2)*par(3)) alpha1 = alpha + d_one CALL gammafunction ( alpha1, deriv, gammln, dgammln, fail, message_gamma ) logb1 = LOG(b1) logb2 = LOG(b2) logC1 = alpha1*logb1 - gammln logC2 = alpha1*logb2 - gammln DO i = 1, numradius r = radius(i) logr = LOG(r) arg1 = logC1 + alpha*logr arg2 = logC2 + alpha*logr n1 = EXP(arg1 - b1*r) n2 = EXP(arg2 - b2*r) nwithr(i) = d_half * ( n1 + n2 ) END DO ! 4 LOG-NORMAL with rg and sigma given ELSE IF ( idis == 4 ) THEN logrg = dlog(par(1)) logsi = dabs(dlog(par(2))) logsi_inv = d_one / logsi C = logsi_inv / root2p DO i = 1, numradius r = radius(i) logr = LOG(r) arg = ( logr - logrg ) / logsi argsq = arg * arg nwithr(i) = C * dexp( - d_half * argsq ) / r END DO ! 5 LOG-NORMAL with reff and veff given * ELSE IF ( idis == 5 ) THEN alpha1 = d_one + par(2) alpha2 = dlog(alpha1) rg = par(1)/(d_one+par(2))**2.5_dp logrg = dlog(rg) logsi = dsqrt(alpha2) logsi_inv = d_one / logsi C = logsi_inv / root2p DO i = 1, numradius r = radius(i) logr = LOG(r) arg = ( logr - logrg ) / logsi argsq = arg * arg nwithr(i) = C * dexp( - d_half * argsq ) / r END DO ! 6 POWER LAW * ELSE IF ( idis == 6 ) THEN alpha = par(1) rmin = par(2) rmax = par(3) alpha1 = alpha - d_one fac1 = (rmax/rmin)**alpha1 fac2 = d_one / ( fac1 - d_one ) C = alpha1 * rmax**alpha1 * fac2 DO i = 1, numradius r = radius(i) if ( (r < rmax) .and. (r > rmin) ) then nwithr(i) = C*r**(-alpha) else nwithr(i) = d_zero endif END DO ! 7 MODIFIED GAMMA with alpha, rc and gamma given ELSE IF ( idis == 7 ) THEN alpha = par(1) rc = par(2) gamma = par(3) b = alpha / (gamma*rc**gamma) logb = LOG(b) alpha1 = alpha + d_one gamma1 = d_one / gamma aperg = alpha1/gamma CALL gammafunction ( aperg, deriv, gammln, dgammln, fail, message_gamma ) IF ( fail ) go to 240 logC = dlog(gamma) + aperg*logb - gammln DO i = 1, numradius r = radius(i) logr = LOG(r) arg1 = logC + alpha*logr r3 = b*r ** gamma nwithr(i) = EXP ( arg1 - r3 ) END DO ! 8 MODIFIED GAMMA with alpha, b and gamma given ELSE IF ( idis == 8 ) THEN alpha = par(1) b = par(2) gamma = par(3) alpha1 = alpha + d_one gamma1 = d_one / gamma logb = LOG(b) aperg = alpha1/gamma CALL gammafunction ( aperg, deriv, gammln, dgammln, fail, message_gamma ) IF ( fail ) go to 240 logC = dlog(gamma) + aperg*logb - gammln DO i = 1, numradius r = radius(i) logr = LOG(r) arg1 = logC + alpha*logr r3 = r ** gamma nwithr(i) = EXP ( arg1 - b*r3 ) END DO END IF ! normal return RETURN ! special return 240 CONTINUE failnod = .TRUE. write(cdis,'(I1)')idis message = message_gamma(1:LEN(message_gamma))//', distribution : '//cdis RETURN END SUBROUTINE sizedis_nod SUBROUTINE rminmax( idis, par, cutoff, rmin, rmax, message, fail ) USE Mie_precision USE MIE_constants, ONLY : d_zero, d_half, d_one, d_two, d_three IMPLICIT NONE ! subroutine arguments REAL (KIND=dp), INTENT (IN) :: par(3), cutoff INTEGER , INTENT (IN) :: idis CHARACTER*(*) , INTENT (OUT) :: message REAL (KIND=dp), INTENT (OUT) :: rmin, rmax LOGICAL , INTENT (OUT) :: fail ! local variables REAL (KIND=dp) :: r(1), nwithr(1),ref,rref,sef,r0,r1,eps LOGICAL :: failnod !************************************************************************ !* Find the integration bounds rmin and rmax for the integration over * !* a size distribution. These bounds are chosen such that the size * !* distribution falls below the user specified cutoff. It is essential * !* that the size distribution is normalized such that the integral * !* over all r is equal to one ! * !* This is programmed rather clumsy and will in the future be changed * !************************************************************************ fail = .FALSE. eps = 1.0E-10_dp IF (idis == 0) THEN rmin= par(1) rmax= par(1) RETURN ELSE IF ( idis == 1) THEN sef = d_one/SQRT(par(2)+d_three) ref = d_one/(sef*sef*par(2)) rref = ref ELSE IF ( idis == 2) THEN ref = par(1) sef = SQRT(par(2)) rref= ref ELSE IF ( idis == 3) THEN sef = SQRT(par(3)) ref = MAX(par(1),par(2))+sef rref= MAX(par(1),par(2)) ELSE IF ( idis == 4) THEN sef = SQRT(EXP(LOG(par(2))**d_two)-d_one) ref = par(1)*(d_one+sef*sef)**0.4_dp rref= ref ELSE IF ( idis == 5) THEN ref = par(1) sef = SQRT(ref) rref= ref ELSE IF ( idis == 6) THEN rmin= par(2) rmax= par(3) RETURN ELSE IF ( idis == 7) THEN ref = par(2) sef = d_two*ref rref= d_half*ref ELSE IF ( idis == 8) THEN ref = (par(1)/(par(2)*par(3)))**par(3) sef = d_two*ref rref= d_half*ref END IF !************************************************************************ !* search for a value of r such that the size distribution !* is less than the cutoff. Start the search at ref+sef which * !* guarantees that such a value will be found on the TAIL of the * !* distribution. * !************************************************************************ r(1) = ref+sef r0 = ref 200 CONTINUE CALL sizedis_nod( idis, par, 1, r, nwithr, message, failnod ) IF ( failnod ) GO TO 899 IF ( nwithr(1) > cutoff) THEN r0 = r(1) r(1) = d_two*r(1) goto 200 END IF r1 = r(1) !************************************************************************ !* Now the size distribution assumes the cutoff value somewhere * !* between r0 and r1 Use bisection to find the corresponding r * !************************************************************************ 300 CONTINUE r(1) = d_half*(r0+r1) CALL sizedis_nod( idis, par, 1, r, nwithr, message, failnod ) IF ( failnod ) GO TO 899 IF ( nwithr(1) > cutoff) THEN r0 = r(1) ELSE r1 = r(1) END IF IF ((r1-r0) > eps) GOTO 300 rmax = d_half*(r0+r1) !************************************************************************ !* Search for a value of r on the low end of the size distribution * !* such that the distribution falls below the cutoff. There is no * !* guarantee that such a value exists, so use an extra test to see if * !* the search comes very near to r = 0 * !************************************************************************ r1 = rref r0 = d_zero 400 CONTINUE r(1) = d_half*r1 CALL sizedis_nod( idis, par, 1, r, nwithr, message, failnod ) IF ( failnod ) GO TO 899 IF ( nwithr(1) > cutoff) THEN r1 = r(1) IF (r1 > eps) GOTO 400 ELSE r0 = r(1) END IF !************************************************************************ !* Possibly the size distribution goes through cutoff between r0 * !* and r1 try to find the exact value of r where this happens by * !* bisection. * !* In case there is no solution, the algorithm will terminate soon. * !************************************************************************ 500 CONTINUE r(1) = d_half*(r0+r1) CALL sizedis_nod( idis, par, 1, r, nwithr, message, failnod ) IF ( failnod ) GO TO 899 IF ( nwithr(1) > cutoff) THEN r1 = r(1) ELSE r0 = r(1) END IF IF ( (r1-r0) > eps) GOTO 500 IF (r1 <= eps) THEN rmin = d_zero ELSE rmin = d_half*(r0+r1) END IF ! normal return RETURN ! error return 899 CONTINUE fail = .TRUE. RETURN END SUBROUTINE rminmax !************************************************************************ !************************************************************************ !************************************************************************ !************************************************************************ !************************************************************************ !************************************************************************ !************************************************************************ !************************************************************************ SUBROUTINE gammafunction & ( xarg, do_derivative, gammln, dgammln, gammafail, message ) !************************************************************************ !* Return the value of the natural logarithm of the gamma function. * !* and its derivative (the Digamma function) * !* The argument xarg must be real and positive. * !* This function is documented in : * !* * !* W.H. Press et al. 1986, 'Numerical Recipes' Cambridge Univ. Pr. * !* page 157 (ISBN 0-521-30811) * !* * !* When the argument xarg is between zero and one, the relation (6.1.4)* !* on page 156 of the book by Press is used. * !* V.L. Dolman April 18 1989 * !************************************************************************ ! modules USE Mie_precision USE MIE_constants, ONLY : d_zero, d_half, d_one, d_two IMPLICIT NONE !* subroutine arguments REAL (KIND=dp), INTENT (IN) :: xarg LOGICAL , INTENT (IN) :: do_derivative REAL (KIND=dp), INTENT (OUT) :: gammln, dgammln CHARACTER*(*) , INTENT (OUT) :: message LOGICAL , INTENT (OUT) :: gammafail !* local parameters and data REAL (KIND=dp), PARAMETER :: gammaf_eps = 1.d-10 REAL (KIND=dp), PARAMETER :: gammaf_fpf = 5.5D0 REAL (KIND=dp) :: cof(6),stp, c0 DATA cof / 76.18009172947146D0, & -86.50532032941677D0, & 24.01409824083091D0, & -1.231739572450155D0, & 0.1208650973866179D-2, & -0.5395239384953D-5 / DATA c0 / 1.000000000190015D0 / DATA stp / 2.5066282746310005D0 / !* local variables INTEGER :: j REAL (KIND=dp) :: x,xx,xxx,tmp,x1,x2,logtmp,pi,ser,dser,gtmp,dgtmp,pix !* initialize output message = ' ' gammln = d_zero dgammln = d_zero gammafail = .FALSE. !* check for bad input IF (xarg <= d_zero) THEN message = ' gammafunction: called with negative argument xarg' gammafail = .TRUE. RETURN END IF IF (ABS(xarg-d_one) < gammaf_eps) THEN message = ' gammafunction: argument too close to one ' gammafail = .TRUE. RETURN END IF !* set up pi = 4.0_dp * ATAN(d_one) IF (xarg .ge. d_one) THEN xxx = xarg ELSE xxx = xarg + d_two END IF !* Numerical Recipes stuff xx = xxx - d_one x1 = xx + gammaf_fpf x2 = xx + d_half logtmp = LOG(x1) tmp = x2*logtmp-x1 ser = c0 x = xx DO j =1, 6 x = x + d_one ser = ser+cof(j)/x END DO gtmp = tmp + LOG(stp*ser) !* derivative of gammln IF ( do_derivative ) THEN dser = d_zero x = xx DO j = 1, 6 x = x+d_one dser = dser+cof(j)/x/x END DO dgtmp = logtmp - (5.0_dp/x1) - (dser/ser) END IF !* assign output IF ( do_derivative ) THEN IF (xarg > d_one) THEN gammln = gtmp dgammln = dgtmp ELSE pix = pi*(d_one-xarg) gammln = LOG(pix/SIN(pix))-gtmp dgammln = - dgtmp - (pi*COS(pix)/SIN(pix)) - (d_one/(d_one-xarg)) END IF ELSE IF (xarg > d_one) THEN gammln = gtmp ELSE pix = pi*(d_one-xarg) gammln = LOG(pix/SIN(pix))-gtmp END IF END IF RETURN end SUBROUTINE gammafunction !************************************************************************ !************************************************************************ !************************************************************************ !************************************************************************ !************************************************************************ !************************************************************************ !************************************************************************ !************************************************************************ SUBROUTINE mie_gauleg(maxn,n,x1,x2,x,w) USE Mie_precision IMPLICIT NONE !* subroutine arguments INTEGER , INTENT (IN) :: maxn, n REAL (KIND=dp), INTENT (IN) :: x1,x2 REAL (KIND=dp), INTENT (OUT) :: x(maxn),w(maxn) INTEGER :: i, m, j REAL (KIND=dp) :: xm,xl,p1,p2,p3,pp,z,z1,eps eps=3.0e-14_dp m=(n+1)/2 xm=0.5_dp*(x2+x1) xl=0.5_dp*(x2-x1) DO i=1,m z=COS(3.1415926540_dp*(DBLE(i)-0.250_dp)/(n+0.50_dp)) 1 CONTINUE p1=1.0_dp p2=0.0_dp DO j=1,n p3=p2 p2=p1 p1=((2.0_dp*DBLE(j)-1.0_dp)*z*p2-(DBLE(j)-1.0_dp)*p3)/DBLE(j) END DO pp=n*(z*p1-p2)/(z*z-1.0_dp) z1=z z=z1-p1/pp if(dabs(z-z1) > eps) GO TO 1 x(i)=xm-xl*z x(n+1-i)=xm+xl*z w(i)=2.0_dp*xl/((1.0_dp-z*z)*pp*pp) w(n+1-i)=w(i) END DO RETURN END SUBROUTINE mie_gauleg SUBROUTINE develop ( max_Mie_angles, ncoeffs, nangles, & cosines, weights, FMAT, expcoeffs ) ! Based on the Meerhoff Mie code !************************************************************************ !* Calculate the expansion coefficients of the scattering matrix in * !* generalized spherical functions by numerical integration over the * !* scattering angle. * !************************************************************************ ! modules USE Mie_precision USE MIE_constants, ONLY : d_zero, d_half, d_one, d_two, d_three, d_four ! implicit none statement IMPLICIT NONE ! input INTEGER , INTENT (IN) :: max_Mie_angles INTEGER , INTENT (IN) :: ncoeffs, nangles REAL (KIND=dp), INTENT (IN) :: cosines(max_Mie_angles) REAL (KIND=dp), INTENT (IN) :: weights(max_Mie_angles) REAL (KIND=dp), INTENT (IN) :: FMAT(4,max_Mie_angles) ! output REAL (KIND=dp), INTENT (OUT) :: expcoeffs(6,0:max_Mie_angles) ! local variables REAL (KIND=dp) :: P00(max_Mie_angles,2) REAL (KIND=dp) :: P02(max_Mie_angles,2) REAL (KIND=dp) :: P22(max_Mie_angles,2) REAL (KIND=dp) :: P2m2(max_Mie_angles,2) REAL (KIND=dp) :: fmatw(4,max_Mie_angles) INTEGER :: i, j, l, lnew, lold, itmp INTEGER :: index_11, index_12, index_22, index_33, index_34, index_44 REAL (KIND=dp) :: dl, dl2, qroot6, fac1, fac2, fac3, fl,& sql4, sql41, twol1, tmp1, tmp2, denom, & alfap, alfam ! Initialization qroot6 = -0.25_dp*SQRT(6.0_dp) index_11 = 1 index_12 = 2 index_22 = 3 index_33 = 4 index_34 = 5 index_44 = 6 DO j = 0, ncoeffs DO i = 1, 6 expcoeffs(i,j) = d_zero END DO END DO ! Multiply the scattering matrix F with the weights w for all angles * ! We do this here because otherwise it should be done for each l * DO i = 1, 4 DO j = 1, nangles fmatw(i,j) = weights(j)*FMAT(i,j) END DO END DO ! Start loop over the coefficient index l * ! first update generalized spherical functions, then calculate coefs. * ! lold and lnew are pointer-like indices used in recurrence * lnew = 1 lold = 2 DO l = 0, ncoeffs IF (l == 0) THEN dl = d_zero DO i=1, nangles P00(i,lold) = d_one P00(i,lnew) = d_zero P02(i,lold) = d_zero P22(i,lold) = d_zero P2m2(i,lold)= d_zero P02(i,lnew) = d_zero P22(i,lnew) = d_zero P2m2(i,lnew)= d_zero END DO ELSE dl = DBLE(l) dl2 = dl * dl fac1 = (d_two*dl-d_one)/dl fac2 = (dl-d_one)/dl DO i=1, nangles P00(i,lold) = fac1*cosines(i)*P00(i,lnew) - fac2*P00(i,lold) END DO ENDIF IF (l == 2) THEN DO i=1, nangles P02(i,lold) = qroot6*(d_one-cosines(i)*cosines(i)) P22(i,lold) = 0.25_dp*(d_one+cosines(i))*(d_one+cosines(i)) P2m2(i,lold)= 0.25_dp*(d_one-cosines(i))*(d_one-cosines(i)) P02(i,lnew) = d_zero P22(i,lnew) = d_zero P2m2(i,lnew)= d_zero END DO sql41 = d_zero ELSE IF (l > 2) THEN sql4 = sql41 sql41 = dsqrt(dl2-d_four) twol1 = 2.D0*dl - d_one tmp1 = twol1/sql41 tmp2 = sql4/sql41 denom = (dl-d_one)*(dl2-d_four) fac1 = twol1*(dl-d_one)*dble(l)/denom fac2 = 4.D0*twol1/denom fac3 = dl*((dl-d_one)*(dl-d_one)-d_four)/denom DO i=1, nangles P02(i,lold) = tmp1*cosines(i)*P02(i,lnew) - tmp2*P02(i,lold) P22(i,lold) = (fac1*cosines(i)-fac2)*P22(i,lnew) - fac3*P22(i,lold) P2m2(i,lold)= (fac1*cosines(i)+fac2)*P2m2(i,lnew) - fac3*P2m2(i,lold) END DO END IF itmp = lnew lnew = lold lold = itmp alfap = d_zero alfam = d_zero fl = dl+d_half do i=1, nangles expcoeffs(index_11,l) = expcoeffs(index_11,l) + P00(i,lnew)*fmatw(1,i) alfap = alfap + P22(i,lnew) * (fmatw(1,i)+fmatw(3,i)) alfam = alfam + P2m2(i,lnew) * (fmatw(1,i)-fmatw(3,i)) expcoeffs(index_44,l) = expcoeffs(index_44,l) + P00(i,lnew)*fmatw(3,i) expcoeffs(index_12,l) = expcoeffs(index_12,l) + P02(i,lnew)*fmatw(2,i) expcoeffs(index_34,l) = expcoeffs(index_34,l) + P02(i,lnew)*fmatw(4,i) END DO expcoeffs(index_11,l) = fl*expcoeffs(index_11,l) expcoeffs(index_22,l) = fl*d_half*(alfap+alfam) expcoeffs(index_33,l) = fl*d_half*(alfap-alfam) expcoeffs(index_44,l) = fl*expcoeffs(index_44,l) expcoeffs(index_12,l) = fl*expcoeffs(index_12,l) expcoeffs(index_34,l) = fl*expcoeffs(index_34,l) END DO RETURN END SUBROUTINE develop SUBROUTINE expand ( max_Mie_angles, ncoeffs, nangles, cosines, expcoeffs, FMAT ) ! Based on the Meerhoff Mie code ! Use the expansion coefficients of the scattering matrix in ! generalized spherical functions to expand F matrix ! modules USE Mie_precision USE MIE_constants, ONLY : d_zero, d_one, d_two, d_four ! implicit none statement IMPLICIT NONE ! input INTEGER , INTENT (IN) :: max_Mie_angles INTEGER , INTENT (IN) :: ncoeffs, nangles REAL (KIND=dp), INTENT (IN) :: cosines(max_Mie_angles) REAL (KIND=dp), INTENT (IN) :: expcoeffs(6,0:max_Mie_angles) ! output REAL (KIND=dp), INTENT (OUT) :: FMAT(4,max_Mie_angles) ! local variables REAL (KIND=dp) :: P00(max_Mie_angles,2) REAL (KIND=dp) :: P02(max_Mie_angles,2) INTEGER :: i, j, l, lnew, lold, itmp INTEGER :: index_11, index_12, index_34, index_44 REAL (KIND=dp) :: dl, qroot6, fac1, fac2, sql4, sql41, tmp1, tmp2 ! Initialization qroot6 = -0.25_dp*SQRT(6.0_dp) index_11 = 1 index_12 = 2 index_34 = 5 index_44 = 6 ! Set scattering matrix F to zero DO j = 1, 4 DO i = 1, nangles FMAT(j,i) = d_zero END DO END DO ! Start loop over the coefficient index l ! first update generalized spherical functions, then calculate coefs. ! lold and lnew are pointer-like indices used in recurrence lnew = 1 lold = 2 DO l = 0, ncoeffs IF ( l == 0) THEN ! Adding paper Eqs. (76) and (77) with m=0 DO i=1, nangles P00(i,lold) = d_one P00(i,lnew) = d_zero P02(i,lold) = d_zero P02(i,lnew) = d_zero END DO ELSE dl = DBLE(l) fac1 = (d_two*dl-d_one)/dl fac2 = (dl-d_one)/dl ! Adding paper Eq. (81) with m=0 DO i=1, nangles P00(i,lold) = fac1*cosines(i)*P00(i,lnew) - fac2*P00(i,lold) END DO END IF IF ( l == 2) THEN ! Adding paper Eq. (78) ! sql4 contains the factor dsqrt((l+1)*(l+1)-4) needed in ! the recurrence Eqs. (81) and (82) DO i=1, nangles P02(i,lold) = qroot6*(d_one-cosines(i)*cosines(i)) P02(i,lnew) = d_zero END DO sql41 = d_zero ELSE IF ( l > 2) THEN ! Adding paper Eq. (82) with m=0 sql4 = sql41 sql41 = dsqrt(dl*dl-d_four) tmp1 = (d_two*dl-d_one)/sql41 tmp2 = sql4/sql41 DO i=1, nangles P02(i,lold) = tmp1*cosines(i)*P02(i,lnew) - tmp2*P02(i,lold) END DO END IF ! Switch indices so that lnew indicates the function with ! the present index value l, this mechanism prevents swapping ! of entire arrays. itmp = lnew lnew = lold lold = itmp ! Now add the l-th term to the scattering matrix. ! See de Haan et al. (1987) Eqs. (68)-(73). ! Remember for Mie scattering : F11 = F22 and F33 = F44 DO i=1, nangles FMAT(1,i) = FMAT(1,i) + expcoeffs(index_11,l)*P00(i,lnew) FMAT(2,i) = FMAT(2,i) + expcoeffs(index_12,l)*P02(i,lnew) FMAT(3,i) = FMAT(3,i) + expcoeffs(index_44,l)*P00(i,lnew) FMAT(4,i) = FMAT(4,i) + expcoeffs(index_34,l)*P02(i,lnew) END DO END DO RETURN END SUBROUTINE expand