* Program to calculate water droplet radiative parameters for ATMS 749.
PROGRAM mieATMS749 IMPLICIT DOUBLE PRECISION (a-h,o-z) xlam=0.6328d0 ! Wavelength in um. xnr = 1.333d0 ! real part of refractive index. xni = 1.d-10 ! Imag part of refractive index.
OPEN(1,FILE='DropletsOptics.txt') WRITE(1,*) 'Wavelength_um nr ni D_um s_um^2 sext_um^2 sabs_um^2 SSA g' DO i=1,100 r=DFLOAT(i) D=2.0d0*r ! Droplet diameter in Microns CALL MIEV0(xlam,xnr,xni,D, ! Inputs * s,CEXT,CABS,SingScatAlbedo,AsymmetryParam) ! Outputs WRITE(1,*) xlam,xnr,xni,D,s,CEXT,CABS,SingScatAlbedo,AsymmetryParam ENDDO CLOSE(1) END
* THIS SUBROUTINE CONTAINS BOTH THE MIE AND ADT SOLUTIONS FOR A * SPHERE!
********************************************************************** * Subroutine to compute the sphere function extinction cross section. * Uses anomalous diffraction theory. ********************************************************************* SUBROUTINE ANOM(lam,xnr,xni,d, cext,cabs) IMPLICIT DOUBLE PRECISION (a-h,o-z) DOUBLE PRECISION lam,xnr,xni,d,cext,cabs DOUBLE PRECISION projarea DOUBLE COMPLEX i,rho,arg,func,ior pi = 4.d0*DATAN(1.d0) i = (0.d0,1.d0) projarea=pi*d*d/4.d0 xx = pi*d/lam ior = DCMPLX(xnr,-xni) * Extinction efficiency. rho = 2.d0 * xx * (ior - 1.d0) arg = i * rho CALL K(arg,func) dqext = 4.d0*DREAL(func) * Absorption efficiency. arg = -4.d0 * xx * DIMAG(ior) CALL K(arg,func) dqabs = 2.d0 * DREAL(func) cext=dqext*projarea cabs=dqabs*projarea RETURN END
********************************************************************** * Subroutine to compute the sphere function K(arg). ********************************************************************* SUBROUTINE K(arg,func) DOUBLE COMPLEX arg,func,argsq,arg4th IF (CDABS(arg).LT.0.1d0) THEN argsq = arg*arg arg4th = argsq*argsq func=arg/3.-argsq/8.+arg*argsq/30.-arg4th/144.+arg4th*arg/840. ELSE func = 0.5+CDEXP(-arg)/arg+ (CDEXP(-arg) - 1.)/( arg * arg ) END IF RETURN END
C **************************************************************** ! MIE SCATTERING PROGRAM ! FOR DOCUMENTATION SEE THE PAPER BY WISCOMBE. ! This version was modified to only put out cross sections. * ALL INPUTS AND OUTPUTS ARE DOUBLE PRECISION REAL*8. * The wavelength is DLam. It is in the same units as * SphereDiameter. * DMRR is real part of refractive index. * DMRI is imag part of refractive index > 0. * CEXT is extinction cross section. * CABS is absorption cross section. C **************************************************************** SUBROUTINE MIEV0(DLam,DMRR,DMRI,SphereDiameter, * PhysicalXsec,CEXT,CABS,SingScatAlbedo,AsymmetryParam) ! Outputs COMPLEX*16 IOR DOUBLE PRECISION XX REAL QEXT,QSCA,GQSC REAL*8 CEXT,CABS,CSCA,SingScatAlbedo,AsymmetryParam DOUBLE PRECISION PI,DLam,DMRR,DMRI,SphereDiameter DOUBLE PRECISION PhysicalXsec REAL N2CUT COMPLEX SBACK, S1(450), S2(450) REAL*8 NP1DN LOGICAL NOIMAG, NOANGS COMPLEX*16 FF,AK,DEN,NUM,NTN,DTD,ANM1,BNM1 COMPLEX*16 DZET,DZETN,DZTNP1,DBIGA,SP,SM,SPS,SMS COMPLEX*16 BIGA,ANPM,BNPM,ZINV COMPLEX*16 ANP,BNP,IORINV,ZETN,ZETNP1 COMPLEX*16 TT,AN,BN DIMENSION SP(450),SM(450),SPS(450),SMS(450), * PIN(450),PINM1(450),TAUN(450),TMP(450),BIGA(20150) DOUBLE PRECISION PIN,PINM1,TMP,TAUN,XMU(900) * ,EQBIGA,COEFF,FN,FNP1,RNP1,RN,TWONP1 * ,XINV,PSIN,PSINP1,CHIN,CHINP1,REIOR,RIORIV EQUIVALENCE (S1(1),PIN(1)),(S1(450),PINM1(1)), * (S2(1),TAUN(1)),(S2(450),TMP(1)) * ,(ZETN,DZETN),(ZETNP1,DZTNP1) DATA EPS1/1.E-2/, EPS2/1.E-8/ DATA MAXIT/10000/ C C PI=4.D0*DATAN(1.D0) XX=PI*SphereDiameter/DLAM IOR=DCMPLX(DMRR,-DMRI) PhysicalXsec=PI*(SphereDiameter**2)/4.D0 * WRITE(*,*) 'PhysicalXsec square Microns',PhysicalXsec N2CUT = 3.e-8 NUMANG=0 REIOR = IOR AIMIOR = -DIMAG(IOR) IF(XX.LT.0. .OR. REIOR.LE.0. .OR. AIMIOR.LT.0. .OR. NUMANG.LT.0 * .OR. NUMANG.GT.900) STOP 1000 NN2 = NUMANG+1 NN = NN2/2 NOANGS = NUMANG.EQ.0 NOIMAG = AIMIOR.LE.N2CUT C C CALCULATE NUMBER OF TERMS IN MIE SERIES (A LEAST UPPER BOUND) C USING EMPIRICAL FORMULAE WHICH WERE FITTED FOR SIZE PARAMETERS C UP TO 20,000 C 7 IF(XX.LE.8.0) NT = XX+4.*XX**(1./3.)+1. IF(XX.GT.8.0 ) NT = XX+6.00*XX**(1./3.)+2. NTP1 = NT+1 NTP2=NT+2 C MAKE SURE ARRAY BIGA WILL BE LARGE ENOUGH IF(NTP1.LE.20150) GO TO 10 WRITE(6,8000) NT, XX 8000 FORMAT(///' ESTIMATED LENGTH OF MIE SERIES NT=',I6, * ' FOR SIZE PARAM=',F12.2/' EXCEEDS BIGA DIMENSIONS') STOP 1001 C C C COMPUTE BIGA C 10 XINV = 1.0D0/XX ZINV = XINV/IOR C C PREPARE FOR DOWN-RECURRENCE--- C COMPUTE INITIAL HIGH-ORDER BIGA-N USING LENTZ METHOD C FF = NTP1*ZINV MM = -1 KK = 2*NT+3 AK = (MM*KK)*ZINV DEN = AK NUM = DEN + 1.0D0/FF KOUNT = 1 C 20 KOUNT = KOUNT+1 IF(KOUNT.GT.MAXIT) GO TO 40 IF(CDABS(NUM/AK).GT.EPS1 .AND. CDABS(DEN/AK).GT.EPS1) GO TO 30 C ILL-CONDITIONED CASE--STRIDE TWO TERMS INSTEAD OF ONE MM = -MM KK = KK+2 AK = (MM*KK)*ZINV NTN = AK*NUM + 1.0D0 DTD = AK*DEN + 1.0D0 FF = (NTN/DTD) * FF MM = -MM KK = KK+2 AK = (MM*KK)*ZINV NUM = AK + NUM/NTN DEN = AK + DEN/DTD KOUNT = KOUNT+1 GO TO 20 C 30 TT = NUM/DEN FF = TT*FF C CHECK FOR CONVERGENCE IF(ABS(REAL(TT)-1.0).LT.EPS2 .AND. ABS(DIMAG(TT)).LT.EPS2) GOTO 50 MM = -MM KK = KK+2 AK = (MM*KK)*ZINV NUM = AK + 1.0D0/NUM DEN = AK + 1.0D0/DEN GO TO 20 C 40 WRITE(6,8001) NT, XX, IOR, AK,NUM,DEN,TT,FF 8001 FORMAT(///' CONTINUED FRACTION FOR A-SUB-NT FAILED TO CONVERGE'/ * ' NT=',I6/' X=',E20.8/' REFR INDEX=',2E20.8/' AK=',2E20.8/ * ' NUM=',2E20.8/' DEN=',2E20.8/' TT=',2E20.8/' FF=',2E20.8) STOP 1002 C 50 BIGA(NT) = FF DBIGA = FF C C DOWNWARD RECURRENCE FOR BIGA-N C DO 70 N = 2,NT DBIGA = (NTP2-N)*ZINV - 1.0D0/(((NTP2-N)*ZINV)+DBIGA) 70 BIGA(NTP1-N) = DBIGA IORINV = 1.0D0/IOR RIORIV = 1.0D0/REIOR C INITIALIZE QUANTITIES USED FOR EFFICIENT CALCULATION OF C NUMERICAL COEFFICIENTS IN MIE SERIES FN = 1.0D0 RN = 1.0D0 MM = 1 C INITIALIZE RICATTI-BESSEL FUNCTION ZETA FOR UPWARD RECURRENCE PSIN = DSIN(XX) CHIN = DCOS(XX) PSINP1 = XINV*PSIN-CHIN CHINP1 = XINV*CHIN+PSIN ZETN = DCMPLX(PSIN,CHIN) ZETNP1 = DCMPLX(PSINP1,CHINP1) C INITIALIZE PREVIOUS COEFFICIENTS (A-SUB-N-1, B-SUB-N-1) C FOR USE IN ASYMMETRY FACTOR SERIES ANM1 = (0.0,0.0) BNM1 = (0.0,0.0) C INITIALIZE SUMS FOR EFFICIENCIES, ASYMMETRY FACTOR, AND C BACKSCATTERING AMPLITUDE QEXT = 0.0 QSCA = 0.0 GQSC = 0. SBACK = (0.,0.) C INITIALIZE ANGULAR FCN PIN AND SUMS FOR S+,S- AT ALL ANGLES C ***** VECTORIZABLE LOOP ***** DO 250 J = 1,NN SP(J) = (0.0D0,0.0D0) SM(J) = (0.0D0,0.0D0) SPS(J) = (0.0D0,0.0D0) SMS(J) = (0.0D0,0.0D0) PINM1(J) = 0.0D0 250 PIN(J) = 1.0D0 C DO 500 N = 1,NT C COMPUTE THE VARIOUS NUMERICAL COEFFICIENTS NEEDED FNP1 = FN+1.0D0 TWONP1 = FN+FNP1 RNP1 = 1.0D0/FNP1 COEFF = RN+RNP1 NP1DN = 1.0D0+RN C C CALCULATE THE MIE SERIES COEFFICIENTS LITTLE-A AND LITTLE-B C IF(NOIMAG) GO TO 300 C GENERAL CASE AN = ((IORINV*BIGA(N)+(FN*XINV))*PSINP1-PSIN)/ * ((IORINV*BIGA(N)+(FN*XINV))*ZETNP1-ZETN) BN = (( IOR*BIGA(N)+(FN*XINV))*PSINP1-PSIN)/ * (( IOR*BIGA(N)+(FN*XINV))*ZETNP1-ZETN) C INCREMENT SERIES FOR SCATTERING EFFICIENCY QSCA = QSCA + TWONP1*((REAL(AN))**2+(DIMAG(AN))**2 * +(REAL(BN))**2+(DIMAG(BN))**2) GO TO 350 C 300 CONTINUE C NO-ABSORPTION CASE EQBIGA=BIGA(N) AN = ((RIORIV*EQBIGA+(FN*XINV))*PSINP1-PSIN)/ * ((RIORIV*EQBIGA+(FN*XINV))*ZETNP1-ZETN) BN = (( REIOR*EQBIGA+(FN*XINV))*PSINP1-PSIN)/ * (( REIOR*EQBIGA+(FN*XINV))*ZETNP1-ZETN) C 350 CONTINUE C INCREMENT SUMS FOR ASYMMETRY FACTOR, EXTINCTION EFFICIENCY, AND C BACKSCATTERING AMPLITUDE GQSC = GQSC +(FN-RN)*REAL(ANM1*CONJG(AN)+BNM1*CONJG(BN)) * + COEFF*REAL(AN*CONJG(BN)) QEXT = QEXT + TWONP1*REAL(AN+BN) SBACK = SBACK+(MM*TWONP1)*(AN-BN) IF(NOANGS) GO TO 450 C C PUT MIE SERIES COEFFICIENTS IN FORM NEEDED FOR COMPUTING S+, S- C ANP = COEFF*(AN+BN) BNP = COEFF*(AN-BN) ANPM = MM*ANP BNPM = MM*BNP C C ***** VECTORIZABLE LOOP ***** DO 400 J = 1,NN C ADD UP SUMS WHILE UPWARD RECURSING ANGULAR FUNCTIONS LITTLE PI C AND LITTLE TAU TMP(J) = (XMU(J)*PIN(J)) - PINM1(J) TAUN(J) = FN*TMP(J) - PINM1(J) SP(J) = SP(J) + ANP*(PIN(J)+TAUN(J)) SMS(J) = SMS(J) + BNPM*(PIN(J)+TAUN(J)) SM(J) = SM(J) + BNP*(PIN(J)-TAUN(J)) SPS(J) = SPS(J) + ANPM*(PIN(J)-TAUN(J)) PINM1(J) = PIN(J) PIN(J) = (XMU(J)*PIN(J)) + NP1DN*TMP(J) 400 CONTINUE C C UPDATE RELEVANT QUANTITIES FOR NEXT PASS THROUGH LOOP 450 MM = - MM FN = FNP1 RN = RNP1 ANM1 = AN BNM1 = BN C CALCULATE RICATTI-BESSEL FUNCTIONS BY UPWARD RECURRENCE DZET = (TWONP1*XINV)*DZTNP1-DZETN DZETN = DZTNP1 DZTNP1 = DZET PSIN = PSINP1 PSINP1 = ZETNP1 500 CONTINUE C C C MULTIPLY SUMS BY APPROPRIATE FACTORS TO GET QEXT, QSCA, GQSC C QEXT = 2.*(XINV**2)*QEXT QSCA = 2.*(XINV**2)*QSCA IF(NOIMAG) QSCA = QEXT GQSC = 4.*(XINV**2)*GQSC SBACK = 0.5*SBACK
CEXT=DREAL(QEXT)*PhysicalXSEC CSCA=DREAL(QSCA)*PhysicalXSEC CABS=CEXT-CSCA AsymmetryParam=DREAL(GQSC/QSCA) SingScatAlbedo=CSCA/CEXT
IF(NOANGS) RETURN
DO 800 J=1,NN C S1(J)=0.5*(SP(J)+SM(J)) S1(J)=SM(J) C S2(J)=0.5*(SP(J)-SM(J)) S2(J)=SP(J) C S1(NN2-J)=0.5*(SPS(J)+SMS(J)) C800 S2(NN2-J)=0.5*(SPS(J)-SMS(J)) S1(NN2-J)=SMS(J) 800 S2(NN2-J)=SPS(J)
RETURN END