program cmie c c..program in order to determine the single scattering properties c for size distributions of spherical particles c ( includes backscattering ratio calculation ) c c three data-files are required: c sizset.d - amounts of particles for selected radii c refrac.d - refractive indices for all wavelengths c sizdis.d - limits and constants for size distributions c implicit real*8 (a-h,o-z) cmc __________________________________________________________________s parameter(nwl = 62, nrm = 36) ! suso, waso, soot, ssam, sscm c parameter(nwl = 62, nrm = 1) ! inso c parameter(nwl = 62, nrm = 8) ! dust dimension RH(nrm),rm(nrm),re(nrm),dd(nrm) dimension wl(nwl),rnw(nwl),riw(nwl),rns(nwl),ris(nwl) dimension QQex(nwl,nrm),QQsc(nwl,nrm), & QQbk(nwl,nrm),QQab(nwl,nrm), & SSal(nwl,nrm),SAal(nwl,nrm), & GG( nwl,nrm),Bex( nwl,nrm),xx(nwl,nrm) character*4 type character*39 header(8) cmc ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~e common /mie00/ rmin,rmax,psda(4),psdb(4),psdc(4),psdd(4),md,mn,nd common /mie01/ envo(75),erad(75),maxin common /mie02/ pi,x,rfr,rfi,crfr,crfi,crat,core,thetd(75),accu, & small common /mie03/ alambd,qext,qscat,ctbrqs,qbs,eltrmx(4,75,2),jx common /mie04/ alr,w0,verh,alwc,gasym,beta55,phafu(4,75,2) common /miemc/ qex,qsc,qbk,qab,reff c real*8 ain(75,2),polr(75,2) cmc _________________ c CHOOSE A TYPE type = 'waso' if (type.ne.'dust') rmax1= 0.5 !max radius at 0.9 RH cmc ~~~~~~~~~~~~~~~~~ data header /'Wl(um) Extinction efficiency ', & 'Wl(um) Scattering efficiency ', & 'Wl(um) Back-scatt efficiency ', & 'Wl(um) Absorption efficiency ', & 'Wl(um) Single-Scatt Albedo ', & 'Wl(um) 1.-Single-Scatt Albedo ', & 'Wl(um) Asymmetry factor ', & 'Wl(um) Mass extinction coefficient m2/g'/ c cm open ( 7,file='wavlen.dat',status='unknown') cm open ( 8,file='phasef.dat',status='unknown') cm open ( 9,file='result.dat',status='unknown') open (10,file='sizset.d',status='old') open (11,file='refrac.d',status='old') c lam = 0 pi = 4.0d+0 *datan(1.0d+0) c c..read 'sizset.d' data read (10,260) maxin,nphase core = 0.0d+0 if(nphase.ge.0) read (10,270) accu if(nphase.lt.0) read (10,280) accu,core if(nphase.lt.0) nphase = -nphase c if(maxin.le.0) goto 110 do 100 i = 1,maxin 100 read (10,290) envo(i),erad(i) rmin = erad(1) rmax = erad(maxin) goto 130 c c..read 'sizdis.d' data 110 open (12,file='sizdis.'//type//'.dat',status='old') write(6,*) 'sizdis.'//type//'.dat' read (12,300) md,mn,nd,rmin0,rmax0 do 120 i=1,4 psda(i) = 0.0d+0 psdb(i) = 0.0d+0 psdc(i) = 0.0d+0 psdd(i) = 0.0d+0 if(mn.lt.i) goto 120 if(md.eq.1) read (12,*) psda(i),psdb(i),psdc(i),psdd(i) if(md.eq.2) read (12,*) psda(i),psdb(i),psdc(i) 120 continue if(md.ne.(-maxin)) goto 240 if(maxin.lt.0. and. rmin0.le.0.0d+0) goto 240 close (12) cmc___________________________________________________s r0 = psdb(1) c== Read size at different RH: open(12,file='sizrad.'//type//'.dat',status='old') read (12,*) do n = 1,nrm read (12,*) RH(n),rm(n),dd(n) end do close (12) open(12,file='refrac.water.dat',status='old') read (12,'(/)') do l = 1,nwl read (12,'(f6.3,f7.3,E10.2)') & wl(l),rnw(l),riw(l) end do close (12) c= Read wavelengths, refractive indices for water and solute c= open(12,file='refrac.'//type//'.dat',status='old') read (12,'(/)') do l = 1,nwl read (12,'(f6.3,f7.3,E10.2)') & wl(l),rns(l),ris(l) end do close (12) cmc~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~e c 130 continue c c..read refractive indices for the particular wavelength creal = 0.0d+0 cimag = 0.0d+0 cmc if(core.eq.0.0d+0) read (11,*) wavel,rreal,rimag if(core.ne.0.0d+0) read (11,*) wavel,rreal,rimag, * creal,cimag cmc _________________________________________________________s if(core.eq.0.0d+0) then do n = 1,nrm !======================> Rm loop fs = r0**3/rm(n)**3 psdb(1) = rm(n) rmin = rmin0 *(rm(n)/r0) rmax = rmax0 *(rm(n)/r0) c------------------------------ cmc let rmax0 change with RH if (type.ne.'dust') then rmax00 = rmax0+(rmax1-rmax0)/(0.9-0.)*RH(n) c rmax=rmax00*(rm(n)/r0) rmax=rmax00 if (n.eq.1.or.n.eq.27) write(6,*) RH(n),rmax00,rmax end if c------------------------------ c do l = 1,nwl !======================> WL loop c if (nphase.le.1) then c..simple (3-point phasefunction) jx = 2 jxx = 3 thetd(1) = 0.0d+0 thetd(2) = 9.0d+1 else if (nphase.eq.2) then c..special forward/backward resolution (135-points) jx = 68 jxx = 2*jx -1 ajx = 5.0d+0 thetd(1) = 0.0d+0 do j = 2,jx aj = ajx if(j.le.68) aj = 5.0d+0 if(j.le.60) aj = 2.0d+0 if(j.le.45) aj = 1.0d+0 if(j.le.35) aj = 5.0d-1 if(j.le.21) aj = 2.0d-1 if(j.le.11) aj = 1.0d-1 thetd(j) = thetd(j-1) +aj end do if(thetd(jx).gt.9.0d+1) thetd(jx) = 9.0d+1 else c..read 'angles.d' data open (12,file='angles.d',status='old') read (12,*) jx read (12,*) (thetd(j),j=1,jx) close (12) jxx = 2*jx -1 end if wavel = wl(l) if (type.eq.'dust'.or.type.eq.'inso') then rreal = rns(l) rimag = ris(l) else rreal = fs*rns(l) + (1.-fs)*rnw(l) rimag = fs*ris(l) + (1.-fs)*riw(l) end if if (n.eq. 1.and.l.eq. 5) write(6,*) 'RH=',RH(n),'nr400=',rreal if (n.eq. 1.and.l.eq.13) write(6,*) 'RH=',RH(n),'nr800=',rreal if (n.eq.27.and.l.eq. 5) write(6,*) 'RH=',RH(n),'nr400=',rreal if (n.eq.27.and.l.eq.13) write(6,*) 'RH=',RH(n),'nr800=',rreal lam = lam+1 alambd = wavel rfr = rreal rfi = rimag crfr = creal crfi = cimag if(lam.eq.1) alr = alambd c c..call spectral subroutine .......................................... c call dsm1 c c..(1-w0) for wavelengths smaller than 3.5um cm ww = w0 cm if(alambd.lt.3.5d+0) ww = 1.0d+0 -w0 ww = 1.0d+0 -w0 c cm if(lam.gt.1) goto 180 cm beta66 = beta55 /1.0d+3 cm if(core.eq.0.0d+0) write (9,310) alwc,beta66 cm if(core.ne.0.0d+0) write (9,320) alwc,beta66,core cm goto 230 c cm180 if(core.eq.0.0d+0) write (9,330) ww,gasym,verh,alambd, cm * rfr,rfi cm if(core.ne.0.0d+0) write (9,340) ww,gasym,verh,alambd, cm * rfr,rfi,crfr,crfi cm230 continue QQex(l,n) = qex QQsc(l,n) = qsc QQbk(l,n) = qbk QQab(l,n) = qab SSal(l,n) = w0 SAal(l,n) = ww GG( l,n) = gasym c c..phasefunction (0.0-90.0) cm do 190 k = 1,2 cm do 190 j = 1,jx cm if(k.eq.1) ja = j cm if(k.eq.2) ja = jx -j+1 cm ain(ja,k) = (phafu(1,ja,k)+phafu(2,ja,k)) cm polr(ja,k) = (phafu(2,ja,k)-phafu(1,ja,k)) /ain(ja,k) cm ain(ja,k) = 0.5d+0 *ain(ja,k) cm190 continue c cm if(core.eq.0.0d+0) write (8,350) jxx,alambd,rfr,rfi cm if(core.ne.0.0d+0) write (8,360) jxx,alambd,rfr,rfi,core, cm * crfr,crfi cm write (8,370) cm do 200 j = 1,jx cm write (8,380) thetd(j),ain(j,1),(phafu(i,j,1),i=1,4),polr(j,1) cm200 continue c c..phasefunction (90.0-180.0) cm do 210 j = 1,jx cm thetd(j) = 1.8d+2 -thetd(j) cm210 continue c cm do 220 jj = 2,jx cm j = jx -jj +1 cm write (8,380) thetd(j),ain(j,2),(phafu(i,j,2),i=1,4),polr(j,2) cm220 continue end do !<====================== WL loop re(n) = reff fm = (re(n)**3-re(1)**3)*1./(re(1)**3*dd(1)) + 1. if (type.eq.'dust') fm = 1. Bex (:,n) = 3./4.*QQex(:,n)/dd(n)/re(n)*fm write(6,'(f4.2,6f8.3)') RH(n),rm(n),re(n),rmin,rmax, & QQex(1,n),Bex(1,n) end do !<====================== Rm loop end if write(6,'(a6,f5.3,1x,a6,f4.2,1x,a3,2f6.3)') & 'alpha=', -(log(Bex(5, 1)/Bex(13, 1)))/(log(wl(5)/wl(13))), & 'at RH=',RH( 1),'wl=',wl(5),wl(13) write(6,'(a6,f5.3,1x,a6,f4.2,1x,a3,2f6.3)') & 'alpha=', -(log(Bex(5,27)/Bex(13,27)))/(log(wl(5)/wl(13))), & 'at RH=',RH(27),'wl=',wl(5),wl(13) cmc ________________________________________________________s write(6,'(a7,f10.4)') 'r0 =',r0 write(6,'(a7,f10.3)') 'max wl=',wl(nwl) open(60,file='/misc/mcc03/chin/radiation/mie.'//type//'.out', & status='unknown') write(60,'(a58)') 'Mie parameters for '//type// & ' as a function of RH and wavelength' write(60,'(a6,36f10.2)') 'RH ',RH write(60,'(a6,36f10.4)') 'rm(um)',rm write(60,'(a6,36f10.4)') 're(um)',re open(61,file='/misc/mcc03/chin/radiation/mie.'//type//'.550', & status='unknown') write(61,'(a58)') 'Mie parameters for '//type// & ' as a function of RH and wavelength' write(61,'(a6,36f10.2)') 'RH ',RH write(61,'(a6,36f10.4)') 'rm(um)',rm write(61,'(a6,36f10.4)') 're(um)',re do nt = 1,8 if (nt.eq.1) xx(:,:) = QQex(:,:) if (nt.eq.2) xx(:,:) = QQsc(:,:) if (nt.eq.3) xx(:,:) = QQbk(:,:) if (nt.eq.4) xx(:,:) = QQab(:,:) if (nt.eq.5) xx(:,:) = SSal(:,:) if (nt.eq.6) xx(:,:) = SAal(:,:) if (nt.eq.7) xx(:,:) = GG( :,:) if (nt.eq.8) xx(:,:) = Bex( :,:) write(60,'(a39)') header(nt) write(61,'(a39)') header(nt) do l = 1,nwl if (l.eq.1) write(61,'(f6.3,1p36e10.3)') wl(l),xx(l,:) if (l.gt.1) write(60,'(f6.3,1p36e10.3)') wl(l),xx(l,:) end do end do close (60) close (61) cmc ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~e c cm230 if(wavel.gt.0.0d+0) goto 130 cm if(wavel.eq.0.0d+0) goto 250 if(wavel.eq.wl(nwl)) goto 250 240 write (6,390) 250 continue c c..formats c 260 format (//,i4,/,i4,41x,i4) 270 format ( /,e10.1,/) 280 format ( /,e10.1,23x,e11.3,/) 290 format (2x,e11.4,42x,e11.4) 300 format (/////,3(i6,/),2(/,e11.3),///) 310 format (/,3x,'liquid-water-content ',e10.3,'[g/m3]', * /,3x,'extinction-coefficient',e10.3,'[ 1/m]', * /,4x,'(1-) w0',2x,'aniso-fac',2x,'extin-rel', * 3x,'w [um]',2x,'refractive index') 320 format (/,3x,'liquid-water-content ',e10.3,'[g/m3]', * /,3x,'extinction-coefficient',e10.3,'[ 1/m]', * /,4x,'(1-) w0',2x,'aniso-fac',2x,'extin-rel', * 3x,'w [um]',2x,'refractive index', * 4x,'core: ',e9.3) 330 format (1x,e10.4,2e11.4, * 2x,f7.3,2x,'(',f5.3,',',e9.3,')') 340 format (1x,e10.4,2e11.4, * 2x,f7.3,2x,'(',f5.3,',',e9.3,')', * 2x,'(',f5.3,',',e9.3,')') 350 format (//,3x,i3,2x,f8.3,1x,'um',2x,'(',f5.3,',',e9.3,')') 360 format (//,3x,i3,2x,f8.3,1x,'um',2x,'(',f5.3,',',e9.3,')', * 2x,'core: ',e9.3,1x,'(',f5.3,',',e9.3,')') 370 format (/,8x,'angle intensity', * 1x,'p1/(4*pi) p2/(4*pi) p3/(4*pi) p4/(4*pi)', * 2x,'polaris.') 380 format ( 3x,f10.3,6e10.3) 390 format (3x,'wrong particle size distributions or lower limit') c close (7) close (8) close (9) close (10) close (11) c stop end c-------------------------------------------------------------------- subroutine dsm1 c c..determination of the spectral values c implicit real*8 (a-h,o-z) c common /mie00/ rmin,rmax,psda(4),psdb(4),psdc(4),psdd(4),md,mn,nd common /mie01/ envo(75),erad(75),maxin common /mie02/ pi,x,rfr,rfi,crfr,crfi,crat,core,thetd(75),accu, & small common /mie03/ alambd,qext,qscat,ctbrqs,qbs,eltrmx(4,75,2),jx common /mie04/ alr,w0,verh,alwc,gasym,beta55,phafu(4,75,2) common /miemc/ qex,qsc,qbk,qab,reff c real*8 pha(4,75,2) c welz = 1.0d+4 /alambd rhoh2o = 1.0d+0 xmax = 2.0d+0 *pi*rmax /alambd c c..set to zero an = 0.0d+0 ar = 0.0d+0 fc = 0.0d+0 ga = 0.0d+0 qe = 0.0d+0 qs = 0.0d+0 qb = 0.0d+0 ro = 0.0d+0 vo = 0.0d+0 do 400 k = 1,2 do 400 i = 1,4 do 400 j = 1,jx pha(i,j,k) = 0.0d+0 400 continue c c..number of intervals for 'maxin' data points inter = maxin if(maxin.eq.0) return if(maxin.gt.0) goto 410 c c..number of intervals for formulas (see below) inter = nd if(nd.eq.0) inter = 30 if(nd.lt.0) inter = -nd 410 continue c c..loop for intervals c do 680 n = 1,inter c c..determine mean radius (r) and size parameter (x) if(maxin.gt.0) goto 430 if(nd.lt.0) goto 420 c c..linear intervals rdel = rmax -rmin radd = rdel *dble(2*n-1) /dble(2*inter) ralo = rdel *dble(2*n-2) /dble(2*inter) rahi = rdel *dble(2*n ) /dble(2*inter) deltar = rahi -ralo c deltar = rdel /dble(inter-1) c radd = rdel *dble(n-1) /dble(inter-1) r = rmin +radd goto 440 c c..logarithmic intervals 420 rdel = dlog(rmax) -dlog(rmin) ralo = rdel *dble(2*n-2) /dble(2*inter) radd = rdel *dble(2*n-1) /dble(2*inter) rahi = rdel *dble(2*n ) /dble(2*inter) r = dexp(dlog(rmin) +radd) deltar = dexp(dlog(rmin) +rahi) -dexp(dlog(rmin) +ralo) goto 440 c c..for data column distributions 430 r = erad(n) deltar = 0.0 c 440 x = 2.0d+0 *pi*r /alambd c c..set maximum size-parameter to 500 c if(x.gt.500.0) x = 500.0 c c..determine the core-radius ratio if(core.ge.0.0) crat = core/r if(core.lt.0.0) crat = -core if(crat.gt.1.0d+0) crat = 1.0d+0 care = dabs(core) c c..call mie-parameter subroutine (dsm2=normal, dsm2c=normal/core) c c call dsm2 call dsm2c c c..determine particle density (#/m3) c if(maxin.lt.0) goto 450 anzahl = envo(n) goto 660 c c..modified gamma distributions [n(r) = a*r**b *exp(-c*r**d)] c 450 if(maxin.lt.-1) goto 470 anzahl = 0.0 do 460 i=1,mn envonr = 0.0 epsi = psdc(i) *r**psdd(i) if(epsi.lt.small) epsi = 0.0 if(epsi.gt.1.0d+2) goto 460 envonr = psda(i)* r**psdb(i) *dexp(-epsi) anzahl = anzahl +envonr *deltar 460 continue if(anzahl.eq.0.0) goto 680 goto 660 c c..ln-type distr. [n(r) = a/r *exp(-((ln(r)-ln(b))**2)/(2*ln(c)**2))] c with a = a'/(sqrt(2*pi)*ln(c)) c 470 if(maxin.lt.-2) goto 490 anzahl = 0.0 do 480 i=1,mn envonr = 0.0 epsi = ((dlog(r)-dlog(psdb(i)))**2) /(2.0*(dlog(psdc(i)))**2) if(epsi.lt.small) epsi = 0.0 if(epsi.gt.1.0d+2) goto 480 c envonr = psda(i)/r *dexp(-epsi) envonr = psda(i)/(r*dsqrt(2.0d+0*pi)*dlog(psdc(i))) *dexp(-epsi) anzahl = anzahl +envonr *deltar 480 continue if(anzahl.eq.0.0) goto 680 goto 660 490 continue c c c..space for any other particle size distribution c : : : : : c : : : : : c 660 an = anzahl +an ro = anzahl *r +ro c r2 = r *r ar = anzahl *r2 +ar ga = anzahl *r2 *ctbrqs +ga qe = anzahl *r2 *qext +qe qs = anzahl *r2 *qscat +qs qb = anzahl *r2 *qbs +qb c r3 = r *r *r fc = anzahl *r3 *r3 +fc vo = anzahl *r3 +vo c do 670 k = 1,2 do 670 i = 1,4 do 670 j = 1,jx pha(i,j,k) = pha(i,j,k) +alambd**2 * anzahl/pi *eltrmx(i,j,k) 670 continue c 680 continue c c..determine average values do 690 k = 1,2 do 690 i = 1,4 do 690 j = 1,jx phafu(i,j,k) = pha(i,j,k) /(4.0d+0 *pi**2 *qs) 690 continue c c..radii ran = ro/an c rar = dsqrt (ar/an) c rvo = (vo/an)**(1.0d+0 /3.0d+0) reff = vo/ar c qex = qe /ar qsc = qs /ar qbk = qb /ar qab = qex -qsc c xsect = pi *ar volume = (4.0d+0 /3.0d+0) *pi *vo c c..[um**3/m**3] to [g**3/m**3] alwc = 1.0d-12 *volume *rhoh2o c c..[um**2/m**3] to [1/km] optic = 2.0d-9 *pi *ar betaex = 1.0d-9 *pi *qe betasc = 1.0d-9 *pi *qs betabk = 1.0d-9 *pi *qb c betaab = betaex -betasc if(alambd.eq.alr) beta55 = betaex verh = betaex /beta55 w0 = betasc /betaex w0min = 1.0d+0 - w0 gasym = ga /qs c radar = 6.4d-5 *fc sicht = 3.91d+0 /betaex c c..rayleigh (back-)scattering coair = 2.68684d+25 alair = alambd /1.0d+6 c c..refractive index of air rfair = 2.8760305d+3 +1.62876d+1/alambd**2 +1.3638d-1/alambd**4 rfair = (1.0d+0 +rfair/1.0d+7)**2 -1.0d+0 c c..extinction and backscattering (1/m) exair = coair*(8.0d+0 *pi**3 *rfair**2 * /(3.0d+0 *coair**2 *alair**4)) *1.061 bkair = exair*(3.0d+0/(8.0d+0*pi)) bkrat = betabk /(1.0d+3 *bkair) c cm if(maxin.ge. 0) write (7,700) maxin,rmin,rmax,xmax cm if(maxin.eq.-1) write (7,710) mn, cm * psda(1),psdb(1),psdc(1),psdd(1), cm * psda(2),psdb(2),psdc(2),psdd(2), cm * psda(3),psdb(3),psdc(3),psdd(3), cm * psda(4),psdb(4),psdc(4),psdd(4) cm if(maxin.eq.-2) write (7,720) mn, cm * psda(1),psdb(1),psdc(1), cm * psda(2),psdb(2),psdc(2), cm * psda(3),psdb(3),psdc(3), cm * psda(4),psdb(4),psdc(4) cm if(maxin.lt.0 .and. nd.lt.0) write (7,730) inter,rmin,rmax,xmax cm if(maxin.lt.0 .and. nd.ge.0) write (7,740) inter,rmin,rmax,xmax cm if(core.eq.0.0) write (7,750) alambd,welz,rfr,rfi cm if(core.lt.0.0) write (7,760) alambd,welz, cm * rfr,rfi,crfr,crfi,care cm if(core.gt.0.0) write (7,770) alambd,welz, cm * rfr,rfi,crfr,crfi,care cm write (7,780) an,ran,reff,volume,alwc,qex,qsc,qbk,qab cm write (7,790) xsect,optic,w0min,w0,gasym, cm * betaex,betasc,betabk,betaab,bkrat,radar,sicht c c..formats c 700 format (//,3x,'mie -calculation - detailed results', * ///,3x,'for', * i3,' discrete concentration/radius data-pairs', * //,6x,'smallest particle size [um]',f13.4, * /,6x,'largest particle size [um]',f13.4, * /,6x,'maximum size-parameter (',f11.3,')') 710 format (//,3x,'mie -calculation - detailed results', * ///,3x,'for a (', * i2,'-modal) modified gamma size-distribution', * //,6x,'[n(r)/dr = a*r**b*exp(-c*r**d)', * 3x,'r = (b/(c*d))**(1/d)]', * //,5x,' a b c d ', * 4(/,5x,4(1x,e8.3))) 720 format (//,3x,'mie -calculation - detailed results', * ///,3x,'for a ( ', * i2,'-modal) log-normal size-distribution', * //,6x,'[n(r)/dr = a/r*exp(-((ln(r)-ln(b))**2)', * '/(2*ln(c)**2))]', * /,17x,'a = #/(sqrt(2*pi)*ln(c))', * //,5x,' a(#/m3) b(rm) c(std.d)', * 4(/,5x,3(1x,e8.3))) 730 format (/,4x,i4,' intervals - linear spacing', * /,6x,'lower particle size limit [um]',f13.4, * /,6x,'upper particle size limit [um]',f13.4, * /,6x,'maximum size-parameter (',f11.3,')') 740 format (/,4x,i4,' intervals - logarithmic spacing', * /,6x,'lower particle size limit [um]',f13.4, * /,6x,'upper particle size limit [um]',f13.4, * /,6x,'maximum size-parameter (',f11.3,')') 750 format (/,3x,'wavelength [um]',f13.4, * /,3x,'wavenumber [1/cm]',f13.4, * //,3x,'refractive index - real part ',f13.4, * /,3x,'refractive index - imag part ',e17.4) 760 format (/,3x,'wavelength [um]',f13.4, * /,3x,'wavenumber [1/cm]',f13.4, * //,3x,'shell refractive index - real part ',f13.4, * /,3x,'shell refractive index - imag part ',e17.4, * /,3x,'core refractive index - real part ',f13.4, * /,3x,'core refractive index - imag part ',e17.4, * /,3x,' (core-fraction:',f6.3,')') 770 format (/,3x,'wavelength [um]',f13.4, * /,3x,'wavenumber [1/cm]',f13.4, * //,3x,'shell refractive index - real part ',f13.4, * /,3x,'shell refractive index - imag part ',e17.4, * /,3x,'core refractive index - real part ',f13.4, * /,3x,'core refractive index - imag part ',e17.4, * /,3x,' (core-size:',f7.3,' um)') 780 format (/,3x,'number of droplets [#/m**3] ',e17.4, * //,3x,'calculated parameters :', * / 8x,'mean (rel to number) radius [um]',f13.4, * /,8x,'effective (vol/sur) radius [um]',f13.4, * //,8x,'drop-volume [um**3/m**3]',e17.4, * /,8x,'liquid water (1g/cm3) [g/m**3]',e17.4, * //,3x,'(average) extinction efficiency (q-ex)',e17.4, * /,3x,'(average) scattering efficiency (q-sc)',e17.4, * /,3x,'(average) back-scatt efficiency (q-bk)',e17.4, * /,3x,'(average) absorption efficiency (q-ab)',e17.4) 790 format (/,3x,'cross-section area (x-sect) [um2/m3]',e17.4, * /,3x,'geom.optics limit (2*pi*x-sect) [/km]',e17.4, * //,3x,'co-single-scattering-albedo (1-w0)',e17.4, * //,3x,' single-scattering-albedo (w0)',e17.4, * /,3x,' asymmetry-factor (g)',e17.4, * //,3x,'extinction-coefficient (beta-e) [/km]',e17.4, * /,3x,'scattering-coefficient (beta-s) [/km]',e17.4, * /,3x,'back-scatt.coefficient (beta-b) [/km]',e17.4, * /,3x,'absorption-coefficient (beta-a) [/km]',e17.4, * //,3x,'back-scatt.ratio -1 ,if multipl.by p0/p',e17.4, * //,3x,'radar reflectivity-factor (z) [mm6/m3]',e17.4, * //,3x,'visibility (koschmieder) [km]',e17.4,/) c return end c-------------------------------------------------------------------- subroutine dsm2 c c..subroutine for computing the parameters of the electromagtic c radiation scattered of spheres by determining the c capital a function, using the downward recurrance c relationship (original fortran-code by dave,1968) c c qext - efficiency factor for extinction c qscat - efficiency factor for scattering c ctbrqs - average of cos(angle) * qscat c eltrmx (1..) - real part of 1.complex amplitude c eltrmx (2..) - imag part of 1.complex amplitude c eltrmx (3..) - real part of 2.complex amplitude c eltrmx (4..) - imag part of 2.complex amplitude c implicit real*8 (a-h,o-z) c common /mie02/ pi,x,rfr,rfi,crfr,crfi,crat,core,thetd(75),accu, & small common /mie03/ alambd,qext,qscat,ctbrqs,qbs,eltrmx(4,75,2),jx c real*8 t(5),ta(4),tb(2),tc(2),td(2),te(2), * ph(3,75),tau(3,75),cstht(75),siztht(75) c complex*16 rf,rrf,rrfx,wm1,fna,fnb,tc1,tc2,wfn(2), * acap(50000),fnap,fnbp c equivalence (wfn(1),ta(1)),(fna,tb(1)),(fnb,tc(1)), * (fnap,td(1)),(fnbp,te(1)) c l = 1 n = 1 c if(jx.gt.75) goto 990 small = dsqrt(accu) small = dmin1(small,1.0d-15) c rf = dcmplx (rfr,-rfi) rrf = 1.0d+0/rf rx = 1.0d+0/x rrfx = rrf*rx t(1) = x*x *(rfr*rfr+rfi*rfi) t(1) = dsqrt (t(1)) n1 = idint (1.1d+0 *t(1)) n2 = idint (t(1)) if(n1.ge.50000) goto 950 if(n1.gt. 150) goto 800 n1 = 165 n2 = 150 800 n3 = n1 +1 acap(n3) = dcmplx (0.0d+0,0.0d+0) c do 810 nn = 1,n1 n3 = n1 -nn +1 n4 = n1 -nn +2 tc1 = dble (n4) *rrfx tc2 = 1.0d+0 / (tc1+acap(n4)) acap(n3) = tc1-tc2 810 continue c do 840 j = 1,jx angle = dabs (thetd(j)) if(angle.gt.0.0) goto 820 cstht(j) = 1.0d+0 siztht(j) = 0.0d+0 goto 840 820 if(angle.ge.9.0d+1) goto 830 cstht(j) = dcos(pi*angle/1.8d+2) siztht(j) = 1.0d+0 -cstht(j)*cstht(j) goto 840 830 if(angle.gt.9.0d+1) goto 970 cstht(j) = 0.0d+0 siztht(j) = 1.0d+0 840 continue c do 850 j = 1,jx ph(1,j) = 0.0d+0 ph(2,j) = 1.0d+0 tau(1,j) = 0.0d+0 tau(2,j) = cstht(j) 850 continue c t(1) = dcos (x) t(2) = dsin (x) c wm1 = dcmplx (t(1),-t(2)) wfn(1) = dcmplx (t(2), t(1)) wfn(2) = rx*wfn(1) -wm1 if(dabs(ta(3)).lt.accu) ta(3) = 0.0d+0 if(dabs(ta(4)).lt.accu) ta(4) = 0.0d+0 tc1 = acap(1)*rrf +rx tc2 = acap(1)*rf +rx fna = (tc1*ta(3)-ta(1)) /(tc1*wfn(2)-wfn(1)) fnb = (tc2*ta(3)-ta(1)) /(tc2*wfn(2)-wfn(1)) if(dabs(tb(1)).lt.accu) tb(1) = 0.0d+0 if(dabs(tb(2)).lt.accu) tb(2) = 0.0d+0 if(dabs(tc(1)).lt.accu) tc(1) = 0.0d+0 if(dabs(tc(2)).lt.accu) tc(2) = 0.0d+0 c fnap = fna fnbp = fnb t(1) = 1.5d+0 tb(1) = t(1)*tb(1) tb(2) = t(1)*tb(2) tc(1) = t(1)*tc(1) tc(2) = t(1)*tc(2) c do 860 j = 1,jx eltrmx(1,j,1) = tb(1)*ph(2,j) +tc(1)*tau(2,j) eltrmx(2,j,1) = tb(2)*ph(2,j) +tc(2)*tau(2,j) eltrmx(3,j,1) = tc(1)*ph(2,j) +tb(1)*tau(2,j) eltrmx(4,j,1) = tc(2)*ph(2,j) +tb(2)*tau(2,j) eltrmx(1,j,2) = tb(1)*ph(2,j) -tc(1)*tau(2,j) eltrmx(2,j,2) = tb(2)*ph(2,j) -tc(2)*tau(2,j) eltrmx(3,j,2) = tc(1)*ph(2,j) -tb(1)*tau(2,j) eltrmx(4,j,2) = tc(2)*ph(2,j) -tb(2)*tau(2,j) 860 continue c qext = 2.0d+0 *(tb(1)+tc(1)) qscat = (tb(1)**2 +tb(2)**2 +tc(1)**2 +tc(2)**2) /7.5d-1 ctbrqs = 0.0d+0 rmm = -1.0 qbsr = -2.0 * (tc(1) -tb(1)) qbsi = -2.0 * (tc(2) -tb(2)) c 870 n = n+1 c t(1) = dble (2*n -1) t(2) = dble ( n -1) t(3) = dble (2*n +1) c do 880 j = 1,jx ph(3,j) = (t(1)*ph(2,j)*cstht(j) -dble(n)*ph(1,j)) /t(2) tau(3,j) = cstht(j)*(ph(3,j)-ph(1,j)) -t(1)*siztht(j)*ph(2,j) * +tau(1,j) 880 continue c wm1 = wfn(1) wfn(1) = wfn(2) wfn(2) = t(1)*rx*wfn(1) -wm1 if(dabs(ta(3)).lt.accu) ta(3) = 0.0d+0 if(dabs(ta(4)).lt.accu) ta(4) = 0.0d+0 tc1 = acap(n)*rrf + dble(n)*rx tc2 = acap(n)*rf + dble(n)*rx fna = (tc1*ta(3) -ta(1)) /(tc1*wfn(2) -wfn(1)) fnb = (tc2*ta(3) -ta(1)) /(tc2*wfn(2) -wfn(1)) if(dabs(tb(1)).lt.accu) tb(1) = 0.0d+0 if(dabs(tb(2)).lt.accu) tb(2) = 0.0d+0 if(dabs(tc(1)).lt.accu) tc(1) = 0.0d+0 if(dabs(tc(2)).lt.accu) tc(2) = 0.0d+0 t(5) = dble (n) t(4) = t(1) /(t(5)*t(2)) t(2) = (t(2)*(t(5)+1.0d+0)) /t(5) ctbrqs = ctbrqs +t(2)*(tb(1)*td(1) +tb(2)*td(2) * +tc(1)*te(1) +tc(2)*te(2)) * +t(4)*(td(1)*te(1) +td(2)*te(2)) qext = qext +t(3)*(tb(1)+tc(1)) t(4) = tb(1)**2 +tb(2)**2 +tc(1)**2 +tc(2)**2 qscat = qscat +t(3)*t(4) rmm = -rmm qbsr = qbsr + t(3)*rmm*(tc(1) - tb(1)) qbsi = qbsi + t(3)*rmm*(tc(2) - tb(2)) t(2) = dble (n* (n+1)) t(1) = t(3) /t(2) k = (n/2) *2 c do 900 j = 1,jx eltrmx(1,j,1)= eltrmx(1,j,1) +t(1)*( tb(1)*ph(3,j)+tc(1)*tau(3,j)) eltrmx(2,j,1)= eltrmx(2,j,1) +t(1)*( tb(2)*ph(3,j)+tc(2)*tau(3,j)) eltrmx(3,j,1)= eltrmx(3,j,1) +t(1)*( tc(1)*ph(3,j)+tb(1)*tau(3,j)) eltrmx(4,j,1)= eltrmx(4,j,1) +t(1)*( tc(2)*ph(3,j)+tb(2)*tau(3,j)) if(k.eq.n) goto 890 eltrmx(1,j,2)= eltrmx(1,j,2) +t(1)*( tb(1)*ph(3,j)-tc(1)*tau(3,j)) eltrmx(2,j,2)= eltrmx(2,j,2) +t(1)*( tb(2)*ph(3,j)-tc(2)*tau(3,j)) eltrmx(3,j,2)= eltrmx(3,j,2) +t(1)*( tc(1)*ph(3,j)-tb(1)*tau(3,j)) eltrmx(4,j,2)= eltrmx(4,j,2) +t(1)*( tc(2)*ph(3,j)-tb(2)*tau(3,j)) goto 900 890 eltrmx(1,j,2)= eltrmx(1,j,2) +t(1)*(-tb(1)*ph(3,j)+tc(1)*tau(3,j)) eltrmx(2,j,2)= eltrmx(2,j,2) +t(1)*(-tb(2)*ph(3,j)+tc(2)*tau(3,j)) eltrmx(3,j,2)= eltrmx(3,j,2) +t(1)*(-tc(1)*ph(3,j)+tb(1)*tau(3,j)) eltrmx(4,j,2)= eltrmx(4,j,2) +t(1)*(-tc(2)*ph(3,j)+tb(2)*tau(3,j)) 900 continue c if(t(4).lt.small) goto 920 c do 910 j=1,jx ph(1,j) = ph(2,j) ph(2,j) = ph(3,j) tau(1,j) = tau(2,j) tau(2,j) = tau(3,j) 910 continue c fnap = fna fnbp = fnb if(n.lt.n2) goto 870 if(l.eq.10) goto 950 c n = 1 l = l +1 t(1) = x*x *(rfr*rfr+rfi*rfi) t(1) = dsqrt (t(1)) n1 = idint ((1.1d+0)**l *t(1)) n2 = idint ((1.1d+0)**(l-1) *t(1)) if(n1.ge.50000) goto 950 goto 800 c 920 do 940 j = 1,jx do 940 k = 1,2 do 930 i = 1,4 t(i) = eltrmx(i,j,k) 930 continue eltrmx(1,j,k) = t(3)*t(3) +t(4)*t(4) eltrmx(2,j,k) = t(1)*t(1) +t(2)*t(2) eltrmx(3,j,k) = t(1)*t(3) +t(2)*t(4) eltrmx(4,j,k) = t(2)*t(3) -t(1)*t(4) 940 continue c t(1) = 2.0d+0 *rx*rx qext = qext *t(1) qscat = qscat *t(1) ctbrqs = 2.0d+0 *ctbrqs *t(1) qbs = (qbsr*qbsr + qbsi*qbsi) *rx*rx /(4.0 *pi) goto 990 c 950 write (6,960) n2,l 960 format (/,3x,'upper limit (acap/n2) is not enough: n2=', * i7,' in loop:',i3) goto 990 c 970 write (6,980) angle 980 format (3x,'scattering angle is > 90.0! thetd:',e17.10) c 990 return end c-------------------------------------------------------------------- subroutine dsm2c c c..subroutine for computing the parameters of the electromagtic c radiation scattered of spheres by determining the c capital a function, using the downward recurrance c relationship (original fortran-code by dave,1968) c c [this routine allows a cores within the spheres!] c c qext - efficiency factor for extinction c qscat - efficiency factor for scattering c ctbrqs - average of cos(angle) * qscat c qbs - backscatter cross section c eltrmx (1..) - real part of 1.complex amplitude c eltrmx (2..) - imag part of 1.complex amplitude c eltrmx (3..) - real part of 2.complex amplitude c eltrmx (4..) - imag part of 2.complex amplitude c implicit real*8 (a-h,o-z) c common /mie02/ pi,x,rfr,rfi,crfr,crfi,crat,core,thetd(75),accu, & small common /mie03/ alambd,qext,qscat,ctbrqs,qbs,eltrmx(4,75,2),jx c real*8 t(5),ta(4),tb(2),tc(2),td(2),te(2),tx(2),ty(2), * ph(3,75),tau(3,75),cstht(75),siztht(75) c complex*16 rf,rrf,crf,rrfx,wm1,z1,z2,z3,tc1,tc2,wfn(2), * w(3,50000),acap(50000),fna,fnb,fnap,fnbp, * al1,al2,al3,u(8),dh1,dh2,dhx,ph24,ph21, * store,dummy,dumsq,fma,fmb c equivalence (wfn(1),ta(1)),(wfn(2),ta(3)), * (fna,tb(1)),(fnb,tc(1)), * (fnap,td(1)),(fnbp,te(1)),(z1,tx(1)),(z2,ty(1)) c l = 1 n = 1 c if(jx.gt.75) goto 990 small = dsqrt(accu) small = dmin1(small,1.0d-15) iflag = 1 if(crat.lt.small) iflag = 2 c rf = dcmplx (rfr,-rfi) crf = dcmplx (crfr,-crfi) rrf = 1.0d+0/rf rx = 1.0d+0/x rrfx = rrf*rx t(1) = x*x *(rfr*rfr+rfi*rfi) t(1) = dsqrt (t(1)) n1 = idint (1.1d+0 *t(1)) n2 = idint (t(1)) if(n1.ge.50000) goto 950 if(n1.gt. 150) goto 800 n1 = 165 n2 = 150 800 n3 = n1 +1 acap(n3) = dcmplx (0.0d+0,0.0d+0) c z1 = rf *x z2 = rf *x *crat z3 = crf *x *crat al0 = 2.0d+4 *pi /alambd al1 = al0 *crf al2 = al0 *rf al3 = dcmplx (al0,0.0d+0) w(1,n3) = dcmplx (0.0d+0,0.0d+0) w(2,n3) = dcmplx (0.0d+0,0.0d+0) w(3,n3) = dcmplx (0.0d+0,0.0d+0) c do 810 nn = 1,n1 n3 = n1 -nn +1 n4 = n1 -nn +2 tc1 = dble (n4) *rrfx tc2 = 1.0d+0 / (tc1+acap(n4)) acap(n3) = tc1- tc2 c if(iflag.eq.2) goto 810 tc1 = dble (n4) /x tc2 = 1.0d+0 / (tc1+w(1,n4)) w(1,n3) = tc1 - tc2 tc1 = dble (n4) /z3 tc2 = 1.0d+0 / (tc1+w(2,n4)) w(2,n3) = tc1 - tc2 tc1 = dble (n4) /z2 tc2 = 1.0d+0 / (tc1+w(3,n4)) w(3,n3) = tc1 - tc2 810 continue c do 840 j = 1,jx angle = dabs (thetd(j)) if(angle.gt.0.0) goto 820 cstht(j) = 1.0d+0 siztht(j) = 0.0d+0 goto 840 820 if(angle.ge.9.0d+1) goto 830 cstht(j) = dcos(pi*angle/1.8d+2) siztht(j) = 1.0d+0 -cstht(j)*cstht(j) goto 840 830 if(angle.gt.9.0d+1) goto 970 cstht(j) = 0.0d+0 siztht(j) = 1.0d+0 840 continue c do 850 j = 1,jx ph(1,j) = 0.0d+0 ph(2,j) = 1.0d+0 tau(1,j) = 0.0d+0 tau(2,j) = cstht(j) 850 continue c t(1) = dcos (x) t(2) = dsin (x) c wm1 = dcmplx (t(1),-t(2)) wfn(1) = dcmplx (t(2), t(1)) wfn(2) = rx*wfn(1) -wm1 c if(dabs(ta(3)).lt.accu) ta(3) = 0.0d+0 c if(dabs(ta(4)).lt.accu) ta(4) = 0.0d+0 tc1 = acap(n)*rrf +rx tc2 = acap(n)*rf +rx fna = (tc1*ta(3)-ta(1)) /(tc1*wfn(2)-wfn(1)) fnb = (tc2*ta(3)-ta(1)) /(tc2*wfn(2)-wfn(1)) c if(dabs(tb(1)).lt.accu) tb(1) = 0.0d+0 c if(dabs(tb(2)).lt.accu) tb(2) = 0.0d+0 c if(dabs(tc(1)).lt.accu) tc(1) = 0.0d+0 c if(dabs(tc(2)).lt.accu) tc(2) = 0.0d+0 c if(iflag.eq.2) goto 855 sinx = dsin (tx(1)) siny = dsin (ty(1)) cosx = dcos (tx(1)) cosy = dcos (ty(1)) etx = dexp (tx(2)) ety = dexp (ty(2)) etm = dexp (tx(2) -ty(2)) etp = dexp (tx(2) +ty(2)) c aa = siny * (etp +etm) bb = cosy * (etp -etm) cc = sinx * (etx*etx +1.0d+0) dd = cosx * (etx*etx -1.0d+0) denom = 1.0d+0 +etx*etx * (4.0d+0*sinx*sinx -2.0d+0 +etx*etx) realp = (aa*cc +bb*dd) /denom amagp = (bb*cc -aa*dd) /denom dummy = dcmplx (realp,amagp) c aa = siny *siny -5.0d-1 bb = cosy *siny ph24 = 5.0d-1 + dcmplx (aa,bb) *ety*ety c aa = sinx *siny -cosx *cosy bb = sinx *cosy +cosx *siny cc = sinx *siny +cosx *cosy dd = -sinx*cosy +cosx *siny ph21 = (dcmplx (aa,bb) *etx*ety +dcmplx (cc,dd) *etm) *5.0d-1 c dh2 = z2 / (1.0d+0 +(0.0d+0,1.0d+0) *z2) -1.0d+0 /z2 dh1 = z1 / (1.0d+0 +(0.0d+0,1.0d+0) *z1) -1.0d+0 /z1 dhx = x / (1.0d+0 +(0.0d+0,1.0d+0) * x) -1.0d+0 / x store = (dh2 +n /z2) * (w(3,n) +n /z2) ph24 = ph24 /store store = (dh1 +n /z1) * (w(3,n) +n /z2) ph21 = ph21 /store store = (acap(n) +n /z1) / (w(3,n) +n /z2) dummy = dummy * store dumsq = dummy * dummy c c note: the definitions of u(i) in this program are not the same as c the usubi defined in the article by toon and ackerman. the c corresponding terms are: c usub1 = u(1) usub2 = u(5) c usub3 = u(7) usub4 = dumsq c usub5 = u(2) usub6 = u(3) c usub7 = u(6) usub8 = u(4) c ratio of spherical bessel ftn to spherical henkel ftn = u(8) c u(1) = al3 *acap(n) -al2 *w(1,n) u(2) = al3 *acap(n) -al2 *dhx u(3) = al2 *acap(n) -al3 *w(1,n) u(4) = al2 *acap(n) -al3 *dhx u(5) = al1 * w(3,n) -al2 *w(2,n) u(6) = al2 * w(3,n) -al1 *w(2,n) u(7) = (0.0d+0,-1.0d+0) *(dummy *ph21 -ph24) u(8) = ta(3) /wfn(2) c fna = u(8) * (u(1)*u(5)*u(7) +al1*u(1) -dumsq*al3*u(5))/ * (u(2)*u(5)*u(7) +al1*u(2) -dumsq*al3*u(5)) fnb = u(8) * (u(3)*u(6)*u(7) +al2*u(3) -dumsq*al2*u(6))/ * (u(4)*u(6)*u(7) +al2*u(4) -dumsq*al2*u(6)) c 855 fnap = fna fnbp = fnb t(1) = 1.5d+0 tb(1) = t(1)*tb(1) tb(2) = t(1)*tb(2) tc(1) = t(1)*tc(1) tc(2) = t(1)*tc(2) c do 860 j = 1,jx eltrmx(1,j,1) = tb(1)*ph(2,j) +tc(1)*tau(2,j) eltrmx(2,j,1) = tb(2)*ph(2,j) +tc(2)*tau(2,j) eltrmx(3,j,1) = tc(1)*ph(2,j) +tb(1)*tau(2,j) eltrmx(4,j,1) = tc(2)*ph(2,j) +tb(2)*tau(2,j) eltrmx(1,j,2) = tb(1)*ph(2,j) -tc(1)*tau(2,j) eltrmx(2,j,2) = tb(2)*ph(2,j) -tc(2)*tau(2,j) eltrmx(3,j,2) = tc(1)*ph(2,j) -tb(1)*tau(2,j) eltrmx(4,j,2) = tc(2)*ph(2,j) -tb(2)*tau(2,j) 860 continue c qext = 2.0d+0 *(tb(1)+tc(1)) qscat = (tb(1)**2 +tb(2)**2 +tc(1)**2 +tc(2)**2) /7.5d-1 ctbrqs = 0.0d+0 rmm = -1.0 qbsr = -2.0 * (tc(1) -tb(1)) qbsi = -2.0 * (tc(2) -tb(2)) c 870 n = n+1 c t(1) = dble (2*n -1) t(2) = dble ( n -1) t(3) = dble (2*n +1) c do 880 j = 1,jx ph(3,j) = (t(1)*ph(2,j)*cstht(j) -dble(n)*ph(1,j)) /t(2) tau(3,j) = cstht(j)*(ph(3,j)-ph(1,j)) -t(1)*siztht(j)*ph(2,j) * +tau(1,j) 880 continue c wm1 = wfn(1) wfn(1) = wfn(2) wfn(2) = t(1)*rx*wfn(1) -wm1 c if(dabs(ta(3)).lt.accu) ta(3) = 0.0d+0 c if(dabs(ta(4)).lt.accu) ta(4) = 0.0d+0 tc1 = acap(n)*rrf + dble(n)*rx tc2 = acap(n)*rf + dble(n)*rx fna = (tc1*ta(3) -ta(1)) /(tc1*wfn(2) -wfn(1)) fnb = (tc2*ta(3) -ta(1)) /(tc2*wfn(2) -wfn(1)) c if(dabs(tb(1)).lt.accu) tb(1) = 0.0d+0 c if(dabs(tb(2)).lt.accu) tb(2) = 0.0d+0 c if(dabs(tc(1)).lt.accu) tc(1) = 0.0d+0 c if(dabs(tc(2)).lt.accu) tc(2) = 0.0d+0 c if(iflag.eq.2) goto 885 dhx = -n / x +1.0d+0 / (n / x -dhx) dh2 = -n /z2 +1.0d+0 / (n /z2 -dh2) dh1 = -n /z1 +1.0d+0 / (n /z1 -dh1) store = (dh2 + n /z2) * (w(3,n) +n /z2) ph24 = ph24 /store store = (dh1 + n /z1) * (w(3,n) +n /z2) ph21 = ph21 /store store = (acap(n) + n /z1) / (w(3,n) +n /z2) dummy = dummy * store dumsq = dummy * dummy c u(1) = al3 *acap(n) - al2 *w(1,n) u(2) = al3 *acap(n) - al2 *dhx u(3) = al2 *acap(n) - al3 *w(1,n) u(4) = al2 *acap(n) - al3 *dhx u(5) = al1 * w(3,n) - al2 *w(2,n) u(6) = al2 * w(3,n) - al1 *w(2,n) u(7) = (0.0d+0,-1.0d+0) * (dummy *ph21 -ph24) u(8) = ta(3) / wfn(2) c fma = u(8) * (u(1)*u(5)*u(7) +al1*u(1) -dumsq*al3*u(5)) / * (u(2)*u(5)*u(7) +al1*u(2) -dumsq*al3*u(5)) fmb = u(8) * (u(3)*u(6)*u(7) +al2*u(3) -dumsq*al2*u(6)) / * (u(4)*u(6)*u(7) +al2*u(4) -dumsq*al2*u(6)) c if(cdabs ((fna-fma)/fna) .lt.accu .and. * cdabs ((fnb-fmb)/fnb) .lt.accu) iflag = 2 if(iflag.eq.2) goto 885 fna = fma fnb = fmb c 885 t(5) = dble (n) t(4) = t(1) /(t(5)*t(2)) t(2) = (t(2)*(t(5)+1.0d+0)) /t(5) ctbrqs = ctbrqs +t(2)*(tb(1)*td(1) +tb(2)*td(2) * +tc(1)*te(1) +tc(2)*te(2)) * +t(4)*(td(1)*te(1) +td(2)*te(2)) qext = qext +t(3)*(tb(1)+tc(1)) t(4) = tb(1)**2 +tb(2)**2 +tc(1)**2 +tc(2)**2 qscat = qscat +t(3)*t(4) rmm = -rmm qbsr = qbsr + t(3)*rmm*(tc(1) - tb(1)) qbsi = qbsi + t(3)*rmm*(tc(2) - tb(2)) t(2) = dble (n* (n+1)) t(1) = t(3) /t(2) k = (n/2) *2 c do 900 j = 1,jx eltrmx(1,j,1)= eltrmx(1,j,1) +t(1)*( tb(1)*ph(3,j)+tc(1)*tau(3,j)) eltrmx(2,j,1)= eltrmx(2,j,1) +t(1)*( tb(2)*ph(3,j)+tc(2)*tau(3,j)) eltrmx(3,j,1)= eltrmx(3,j,1) +t(1)*( tc(1)*ph(3,j)+tb(1)*tau(3,j)) eltrmx(4,j,1)= eltrmx(4,j,1) +t(1)*( tc(2)*ph(3,j)+tb(2)*tau(3,j)) if(k.eq.n) goto 890 eltrmx(1,j,2)= eltrmx(1,j,2) +t(1)*( tb(1)*ph(3,j)-tc(1)*tau(3,j)) eltrmx(2,j,2)= eltrmx(2,j,2) +t(1)*( tb(2)*ph(3,j)-tc(2)*tau(3,j)) eltrmx(3,j,2)= eltrmx(3,j,2) +t(1)*( tc(1)*ph(3,j)-tb(1)*tau(3,j)) eltrmx(4,j,2)= eltrmx(4,j,2) +t(1)*( tc(2)*ph(3,j)-tb(2)*tau(3,j)) goto 900 890 eltrmx(1,j,2)= eltrmx(1,j,2) +t(1)*(-tb(1)*ph(3,j)+tc(1)*tau(3,j)) eltrmx(2,j,2)= eltrmx(2,j,2) +t(1)*(-tb(2)*ph(3,j)+tc(2)*tau(3,j)) eltrmx(3,j,2)= eltrmx(3,j,2) +t(1)*(-tc(1)*ph(3,j)+tb(1)*tau(3,j)) eltrmx(4,j,2)= eltrmx(4,j,2) +t(1)*(-tc(2)*ph(3,j)+tb(2)*tau(3,j)) 900 continue c if(t(4).lt.small) goto 920 c do 910 j=1,jx ph(1,j) = ph(2,j) ph(2,j) = ph(3,j) tau(1,j) = tau(2,j) tau(2,j) = tau(3,j) 910 continue c fnap = fna fnbp = fnb if(n.lt.n2) goto 870 if(l.eq.10) goto 950 c n = 1 l = l +1 t(1) = x*x *(rfr*rfr+rfi*rfi) t(1) = dsqrt (t(1)) n1 = idint ((1.1d+0)**l *t(1)) n2 = idint ((1.1d+0)**(l-1) *t(1)) if(n1.ge.50000) goto 950 goto 800 c 920 do 940 j = 1,jx do 940 k = 1,2 do 930 i = 1,4 t(i) = eltrmx(i,j,k) 930 continue eltrmx(1,j,k) = t(3)*t(3) +t(4)*t(4) eltrmx(2,j,k) = t(1)*t(1) +t(2)*t(2) eltrmx(3,j,k) = t(1)*t(3) +t(2)*t(4) eltrmx(4,j,k) = t(2)*t(3) -t(1)*t(4) 940 continue c t(1) = 2.0d+0 *rx*rx qext = qext *t(1) qscat = qscat *t(1) ctbrqs = 2.0d+0 *ctbrqs *t(1) qbs = (qbsr*qbsr + qbsi*qbsi) *rx*rx /(4.0 *pi) goto 990 c 950 write (6,960) n2,l 960 format (/,3x,'upper limit (acap/n2) is not enough: n2=', * i7,' in loop:',i3) goto 990 c 970 write (6,980) angle 980 format (3x,'scattering angle is > 90.0! thetd:',e17.10) c 990 return end c c ********************************************************************** c this subroutine computes mie scattering by a stratified sphere, c i.e. a particle consisting of a spherical core surrounded by a c spherical shell. the basic code used was that described in the c report: " subroutines for computing the parameters of the c electromagnetic radiation scattered by a sphere " j.v. dave, c i b m scientific center, palo alto , california. c report no. 320 - 3236 .. may 1968 . c c the modifications for stratified spheres are described in c toon and ackerman, appl. optics, in press, 1981 c c the parameters in the common statement are defined as follows : c alamda wavelength c x size parameter (2*pi*r/lamda) c crat ratio of core radius and outside radius c (rfr, rfi) refractive index of shell c (crfr, crfi) refractive index of core c thetd(j) angle (degree) between the directions of incident c and scattered radiation. thetd(75max) is < or= 90 c c the definitions for the following symbols can be found in"light c 'light scattering by small particles' c h.c.van de hulst, john wiley & sons, ny, 1957. c qext extinction efficiency p.14,127. c qscat scattering efficiency p.14,127. c ctbrqs average(cosine theta) * qscat, p.128 c eltrmx(i,j,k) elements of the transformation matrix f,p.34,45,125. c i=1: element m sub 2.. i=2: element m sub 1.. c i=3: element s sub 21.. i=4: element d sub 21.. c eltrmx(i,j,1) ith-element of matrix for angle thetd(j).. c eltrmx(i,j,2) ith-element of matrix for angle 180.0 - thetd(j) .. c qbs is the back scatter cross section. c c p.s. - in the original program the dimension of acap was 7000. c for conserving space this should be not much higher than c the value, n=1.1*(nreal**2 + nimag**2)**.5 * x + 1 c - the subroutine computes the capital a function c by making use of the downward recurrence relationship. c - be aware of the following 'equivalence' statements: c wfn(1) = (ta(1),ta(2)) c wfn(2) = (ta(3),ta(4)) c z1 = (tx(1),tx(2)) c z2 = (ty(1),ty(2)) c fna = (tb(1),tb(2)) c fnb = (tc(1),tc(2)) c fnap = (td(1),td(2)) c fnbp = (te(1),te(2)) c fnap, fnbp are the preceding values of fna, fnb. c ********************************************************************** c c note: the definitions of u(i) in this program are not the same as c the usubi defined in the article by toon and ackerman. the c corresponding terms are: c usub1 = u(1) usub2 = u(5) c usub3 = u(7) usub4 = dumsq c usub5 = u(2) usub6 = u(3) c usub7 = u(6) usub8 = u(4) c ratio of spherical bessel ftn to spherical henkal ftn = u(8) c c ********************************************************************** c c up to the statment number 910, eltrmx(i,j,k) means this: c - eltrmx(1,j,k): real part of the first complex amplitude. c - eltrmx(2,j,k): imaginary part of the first complex amplitude. c - eltrmx(3,j,k): real part of the second complex amplitude. c - eltrmx(4,j,k): imaginary part of the second complex amplitude. c k = 1 : for thetd(j) and k = 2 : for 180.0 - thetd(j) c definition of the complex amplitude: van de hulst,p.125. c c **********************************************************************