OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
afrt_phs.f
Go to the documentation of this file.
1 ***********************************************************************
2 c
3 c........................phase matrix main program.....................
4 c
5 c modis look-up table version v1.1
6 c
7 c the program can accept a junge dist.,dermendjian's modified gamma
8 c dist and lognormal distributions (upto three modes).
9 c ifunc=1 (junge dist.); ifunc=2 (dermendjian's modified gamma dist)
10 c ifunc=3 (lognormal uni-modal); ifunc=4 (lognormal bi_modal)
11 c ifunc=5 (lognormal tri_modal)
12 c maximum cut off radius for tables is 100 microns stored in
13 c ddd(15).
14 c
15 c
16  SUBROUTINE af_phs_read(infile,nlamb,nsd,ilam1,ilam2,isd1,isd2,
17  1 jfunc,jnorm,jrgm,kfunc,irh,iset,xtitle,
18  2 xww,xn1,xn2,xrg,xsig,xnpar,xdx,epar)
19 c***********************************************************************
20  implicit real *8 (a-h,o-z)
21  include 'afrt_phs.cmn'
22 c***********************************************************************
23  real*8 xww(nsdt,nw),xn1(nmd,nsdt,nw),xn2(nmd,nsdt,nw),
24  1 xrg(nmd,nsdt,nw),xsig(nmd,nsdt,nw),xnpar(nmd,nsdt,nw),
25  2 xdx(nmd,nsdt,nw),epar(nsdt,npar)
26  integer*4 jfunc(nsdt),jnorm(nsdt),jrgm(nsdt),kfunc(nsdt),
27  1 irh(nsdt),iset(nsdt)
28  integer*4 nlamb,nsd,ilam1,ilam2,isd1,isd2
29  character*80 xtitle(nsdt)
30 c
31  character*255 infile
32  character*80 stitle
33 c
34  data ddd/15*0.0d0/
35 c
36 c open(21,file='phs.dir',status='old')
37 c read(*,'(a)')infile
38 c read(*,'(a)')outdir
39 c
40  len1=index(infile,' ')-1
41 c
42  open(5, file=infile(1:len1),status='old')
43 c
44  nlamb = 0
45  ddd(15)=100.0d0
46  read(5,25)nlamb,nsd
47 c write(*,*) nlamb,nsd
48  read(5,25)ilam1,ilam2
49 c write(*,*) ilam1,ilam2
50  read(5,25)isd1,isd2
51 c write(*,*) isd1,isd2
52 c
53  do i=1,nsd,1
54  read(5,27)xtitle(i)
55  read(5,28)jfunc(i),kfunc(i)
56 c write(6,*)'jfunc,kfunc',jfunc(i),kfunc(i)
57  if(jfunc(i).gt.5)then
58  write(6,50)jfunc(i)
59  stop
60  endif
61  if(jfunc(i).le.2)then
62  read(5,29)stitle
63  read(5,30)(epar(i,ik),ik=1,6)
64  read(5,29)stitle
65  do j=1,nlamb,1
66  read(5,30)xww(i,j),xn1(1,i,j),xn2(1,i,j),xdx(1,i,j)
67 c write(*,*) xww(i,j),xn1(1,i,j),xn2(1,i,j),xdx(1,i,j)
68  enddo
69  else
70  read(5,28)jrgm(i),irh(i),iset(i)
71  read(5,29)stitle
72  jcard=1
73  jcard=jfunc(i)-2
74  do j=1,nlamb,1
75  do k=1,jcard,1
76  read(5,30)xww(i,j),xn1(k,i,j),xn2(k,i,j),xrg(k,i,j),
77  1 xsig(k,i,j),xnpar(k,i,j),xdx(k,i,j)
78 c write(*,*) xww(i,j),xn1(k,i,j),xn2(k,i,j),xrg(k,i,j),
79 c 1 xsig(k,i,j),xnpar(k,i,j),xdx(k,i,j)
80  enddo
81  enddo
82  endif
83  enddo
84  close (5)
85 c.....format statements........
86 25 format(3i5)
87 26 format(4i5)
88 27 format( a)
89 28 format(3i5)
90 29 format(a)
91 30 format(1p7e11.4)
92 35 format(1p8e10.3)
93 36 format(i5,1pe10.3)
94 40 format(10x,1p7e10.3)
95 50 format('jfunc is=',i2,' max. value allowed is 5')
96 c**********************************************************************
97  return
98  end
99 c***********************************************************************
100  SUBROUTINE af_phs_process(outdir,ailam,aisd,aisd1,aisd2,
101  1 jfunc,jnorm,jrgm,kfunc,irh,iset,
102  2 xtitle,xww,xn1,xn2,xrg, xsig,xnpar,xdx,epar,
103  3 ormin,ormax,odelr,odelx,
104  4 or11,or22,or33,or44,oreff,oveff,occnsml,obsr,
105  5 osalb,oasf,oqscat,oqext,oangl,ophfu,opol,ot,othd,
106  6 orbar,odnzp,odnp,odvp,osumnp,owt)
107 c***********************************************************************
108  implicit real *8 (a-h,o-z)
109  include 'afrt_phs.cmn'
110 c***********************************************************************
111  integer*4 ailam,aisd,aisd1,aisd2,isd0,jsd,nmsz,asd
112  integer*4 jfunc(nsdt),jnorm(nsdt),jrgm(nsdt),kfunc(nsdt),
113  1 irh(nsdt),iset(nsdt)
114  real*8 xww(nsdt,nw),xn1(nmd,nsdt,nw),xn2(nmd,nsdt,nw),
115  1 xrg(nmd,nsdt,nw),xsig(nmd,nsdt,nw),xnpar(nmd,nsdt,nw),
116  2 xdx(nmd,nsdt,nw),epar(nsdt,npar)
117  real*8 ormin(nmd),ormax(nmd),odelr(nmd),odelx(nmd),
118  1 or11,or22,or33,or44,
119  2 oreff,oveff,occnsml,obsr,
120  3 osalb,oasf,oqscat,oqext,
121  4 oangl(ntf),ophfu(ntf),opol(ntf),
122  5 ot(ntf,nstk),othd(ntf)
123  real*8 orbar(ndr),odnzp(ndr),odnp(ndr),
124  1 odvp(ndr),owt(ndr),osumnp(ndr)
125  character*255 xtitle(nsdt)
126  character*255 outdir
127  character*255 xname,outz,tmat,pf,dnr
128  character*4 cxwwz
129  character*2 ciband,crh,cisd,cset
130  integer*4 ixwwz
131 c
132 c
133  il=ailam
134  isd=aisd
135  isd1=aisd1
136  isd2=aisd2
137  khum=irh(isd)
138  nmsz=10
139  asd = mod(isd-1,nmsz) + 1
140 c
141  ixwwz=dnint(xww(isd,il)*1.0e3+0.01)
142 c if(xww(isd,il) .lt. 1.0)ixwwz=ixwwz/10
143  call convtc(il,2,ciband)
144  call convtc(asd,2,cisd)
145  call convtc(irh(isd),2,crh)
146  call convtc(iset(isd),2,cset)
147  call convtc(ixwwz,4,cxwwz)
148 
149 c write(*,*) ixwwz, xww(isd,il), cxwwz,' ',isd,' of ',isd2
150 c
151  xname='wl'//cxwwz//'sd'//cisd//'rh'//crh//'_set'//cset//'.dat'
152 c
153  len2=index(outdir,' ')-1
154  kx1=index(xname,' ')-1
155  outz=outdir(1:len2)//'/phs_o'//xname(1:kx1)
156  tmat=outdir(1:len2)//'/phs_t'//xname(1:kx1)
157  pf=outdir(1:len2)//'/phs_p'//xname(1:kx1)
158  dnr=outdir(1:len2)//'/phs_dn'//xname(1:kx1)
159 c
160  ko1=index(outz,' ')-1
161  kt1=index(tmat,' ')-1
162  kp1=index(pf,' ')-1
163  kn1=index(dnr,' ')-1
164 c
165 cbaf open(unit=2,file=outz(1:ko1),status='unknown')
166  open(unit=7,file=tmat(1:kt1),form='formatted',status='unknown')
167  open(unit=8,file=pf(1:kp1),form='formatted',status='unknown')
168  open(unit=4,file=dnr(1:kn1),form='formatted',status='unknown')
169 c
170  call pfunc(il,isd,jfunc,kfunc,jrgm,xtitle,xww,xn1,xn2,
171  1 xrg, xsig,xnpar,xdx,epar,ormin,ormax,odelr,odelx,
172  2 or11,or22,or33,or44, oreff,oveff,occnsml,obsr,
173  3 osalb,oasf,oqscat,oqext, oangl,ophfu,opol,ot,othd,
174  4 orbar,odnzp,odnp,odvp,osumnp,owt)
175 c
176 cbaf close(2)
177  close(7)
178  close(8)
179  close(4)
180 c
181 c**********************************************************************
182 c.....format statements........
183 25 format(3i5)
184 26 format(4i5)
185 27 format( a)
186 28 format(3i5)
187 29 format(a)
188 30 format(1p7e11.4)
189 35 format(1p8e10.3)
190 36 format(i5,1pe10.3)
191 40 format(10x,1p7e10.3)
192 50 format('jfunc is=',i2,' max. value allowed is 5')
193 c**********************************************************************
194  return
195  end
196 c***********************************************************************
197  subroutine flbfr(il,isd,ifunc,xrg,xsig,xnpar,epar)
198 c
199 c***********************************************************************
200 
201  implicit real *8 (a-h,o-z)
202 c
203  include 'afrt_phs.cmn'
204 c
205 c***********************************************************************
206 c
207  real*8 xrg(nmd,nsdt,nw),xsig(nmd,nsdt,nw),xnpar(nmd,nsdt,nw),
208  1 epar(nsdt,npar)
209 
210  if(ifunc.le.2)then
211  do i=1,6
212  ddd(i)=epar(isd,i)
213  enddo
214  endif
215 c
216  if(ifunc.ge.3 .and. ifunc.le.5)then
217  nmode=ifunc-2
218  do i=1,nmode
219  mm=(i-1)*3
220  ddd(mm+1)=xrg(i,isd,il)
221  ddd(mm+2)=xsig(i,isd,il)
222  ddd(mm+3)=xnpar(i,isd,il)
223  enddo
224 c ddd(2)=dlog(10.0d0**0.35d0)
225 c ddd(5)=dlog(10.0d0**0.40d0)
226  endif
227 c
228  return
229  end
230 c***********************************************************************
231  subroutine pfunc(il,isd,jfunc,kfunc,jrgm,xtitle,xww,xn1,xn2,
232  1 xrg,xsig,xnpar,xdx,epar,ormin,ormax,odelr,odelx,
233  2 or11,or22,or33,or44,oreff,oveff,occnsml,obsr, osalb,oasf,oqscat,
234  3 oqext,oangl,ophfu,opol,ot,othd,orbar,odnzp,odnp,odvp,osumnp,owt)
235 c
236 c***********************************************************************
237 
238  implicit real *8 (a-h,o-z)
239 c
240  include 'afrt_phs.cmn'
241 c
242 c***********************************************************************
243  integer*4 jfunc(nsdt),jnorm(nsdt),jrgm(nsdt),kfunc(nsdt),
244  1 irh(nsdt),iset(nsdt)
245  real*8 xlamb,xww(nsdt,nw),xn1(nmd,nsdt,nw),xn2(nmd,nsdt,nw)
246  real*8 xrg(nmd,nsdt,nw),xsig(nmd,nsdt,nw),xnpar(nmd,nsdt,nw),
247  1 xdx(nmd,nsdt,nw),epar(nsdt,npar)
248  real*8 ormin(nmd),ormax(nmd),odelr(nmd),odelx(nmd),
249  1 or11,or22,or33,or44,
250  2 oreff,oveff,occnsml,obsr,
251  3 osalb,oasf,oqscat,oqext,
252  4 oangl(ntf),ophfu(ntf),opol(ntf),
253  5 ot(ntf,nstk),othd(ntf)
254  real*8 orbar(ndr),odnzp(ndr),odnp(ndr),
255  1 odvp(ndr),owt(ndr),osumnp(ndr)
256  character*80 xtitle(nsdt)
257  character*80 title
258  data nd/5,6,3,6,9/
259 c
260 c
261  pi=dacos(-1.0d0)
262  ifunc=jfunc(isd)
263  mfunc=kfunc(isd)
264  irgm=jrgm(isd)
265  nterms=nd(ifunc)
266  xlamb=xww(isd,il)
267  title=xtitle(isd)
268 
269  call flbfr(il,isd,ifunc,xrg,xsig,xnpar,epar)
270 c
271  call norz(ifunc,nn,il,isd,irgm,title,xlamb,xn1,xn2,xnpar,xdx,
272  1 orbar,odnzp,odnp,odvp,osumnp,owt)
273 c
274 c if ifunc=3 then compute ccn
275 c
276  ccnsml=-777.0d0
277  if(ifunc.eq.3)then
278  call ccn(ifunc,irgm)
279  endif
280 c
281 c compute moments numerically
282 c
283  call moment(ifunc,nn)
284  izia=0
285  if(izia.eq.1)stop
286 c
287  jx=901
288  do i=1,1801
289  do j=1,4
290  t(i,j) = 0.0d0
291  enddo
292  enddo
293  zwt=0.0d0
294  thetd(1)=0.0d0
295  dthe=1.0d-1
296  do i=1,900
297  thetd(i+1)=thetd(i)+dthe
298  enddo
299 c
300  rfr1=xn1(1,isd,il)
301  rfi1=xn2(1,isd,il)
302  rfr2=xn1(2,isd,il)
303  rfi2=xn2(2,isd,il)
304  rfr3=xn1(3,isd,il)
305  rfi3=xn2(3,isd,il)
306  yyz1=xdx(1,isd,il)
307  yyz2=xdx(2,isd,il)
308  yyz3=xdx(3,isd,il)
309  xxz1=(yyz1*xlamb)/(2.0d0*pi)
310  xxz2=(yyz2*xlamb)/(2.0d0*pi)
311  xxz3=(yyz3*xlamb)/(2.0d0*pi)
312 c
313 c yyz=yyz1
314 c xxz=xxz1
315 c rfr=rfr1
316 c rfi=rfi1
317 c
318  qext=0.0d0
319  qscat=0.0d0
320  do 2 i=1,nn
321  x=(2.0d0*pi*rbar(i))/xlamb
322 c
323  if(ifunc.le.3)then
324  wt1=wt(i)
325  call dbmie(x,rfr1,rfi1,thetd,jx,qext1,qscat1,ctbrqs,
326  1 eltrmx,wt1,t)
327  factor=wt1*pi*rbar(i)**2*1.0e-8
328  qext=qext+qext1*factor
329  qscat=qscat+qscat1*factor
330  zwt = zwt + wt1
331  endif
332 c
333  if(ifunc.eq.4)then
334  wt1=dn1(i)*deltar(i)
335  wt2=dn2(i)*deltar(i)
336  if(rbar(i).ge.rmin1 .and. rbar(i).le.rmax1)then
337  call dbmie(x,rfr1,rfi1,thetd,jx,qext1,qscat1,ctbrqs,
338  1 eltrmx,wt1,t)
339  factor=wt1*pi*rbar(i)**2*1.0e-8
340  qext=qext+qext1*factor
341  qscat=qscat+qscat1*factor
342  zwt = zwt + wt1
343  endif
344  if(rbar(i).ge.rmin2 .and. rbar(i).le.rmax2)then
345  call dbmie(x,rfr2,rfi2,thetd,jx,qext1,qscat1,ctbrqs,
346  1 eltrmx,wt2,t)
347  factor=wt2*pi*rbar(i)**2*1.0e-8
348  qext=qext+qext1*factor
349  qscat=qscat+qscat1*factor
350  zwt = zwt + wt2
351  endif
352  endif
353 c
354  if(ifunc.eq.5)then
355  wt1=dn1(i)*deltar(i)
356  wt2=dn2(i)*deltar(i)
357  wt3=dn3(i)*deltar(i)
358  if(rbar(i).ge.rmin1 .and. rbar(i).le.rmax1)then
359  call dbmie(x,rfr1,rfi1,thetd,jx,qext1,qscat1,ctbrqs,
360  1 eltrmx,wt1,t)
361  factor=wt1*pi*rbar(i)**2*1.0e-8
362  qext=qext+qext1*factor
363  qscat=qscat+qscat1*factor
364  zwt = zwt + wt1
365  endif
366  if(rbar(i).ge.rmin2 .and. rbar(i).le.rmax2)then
367  call dbmie(x,rfr2,rfi2,thetd,jx,qext1,qscat1,ctbrqs,
368  1 eltrmx,wt2,t)
369  factor=wt2*pi*rbar(i)**2*1.0e-8
370  qext=qext+qext1*factor
371  qscat=qscat+qscat1*factor
372  zwt = zwt + wt2
373  endif
374  if(rbar(i).ge.rmin3 .and. rbar(i).le.rmax3)then
375  call dbmie(x,rfr3,rfi3,thetd,jx,qext1,qscat1,ctbrqs,
376  1 eltrmx,wt3,t)
377  factor=wt3*pi*rbar(i)**2*1.0e-8
378  qext=qext+qext1*factor
379  qscat=qscat+qscat1*factor
380  zwt = zwt + wt3
381  endif
382  endif
383 2 continue
384 c
385  alambd=xlamb*1.0d-4
386  const=((alambd**2)/(4.0*pi**2*qscat)) * 4.0*pi
387  do i=1,1801
388  rr(i,1)=t(i,1) * const
389  rr(i,2)=t(i,2) * const
390  rr(i,3)=t(i,3) * const
391  rr(i,4)=t(i,4) * const
392  if(i .le. jx) xx=thetd(i)
393  if(i .gt. jx) xx=180.0d0-thetd(1801-i+1)
394  enddo
395 c
396  do i=1,1801
397  angl(i) = dfloat(i-1)/10.0d0
398  if(i .le. jx) xx=(thetd(i)*pi)/180.0d0
399  if(i .gt. jx) xx=((180.0d0-thetd(1801-i+1))*pi)/180.0d0
400  y =dcos(xx)
401  rr(i,1)=const*(t(i,1)+t(i,2)*y**2+2.0d0*t(i,3)*y)
402  rr(i,2)=const*(t(i,2)+t(i,1)*y**2+2.0d0*t(i,3)*y)
403  rr(i,3)=const*((t(i,1)+t(i,2))*y+t(i,3)*(1.0d0+y**2))
404  rr(i,4)=const*(1.d0-y**2)*t(i,4)
405  phfu(i)=(rr(i,1)+rr(i,2))/2.0
406  pol(i)=(rr(i,2)-rr(i,1))/(rr(i,1)+rr(i,2))
407  enddo
408 c
409 c compute asymmetry factor and backscattering ratio
410 c
411  call asym(angl,phfu,asf,bsr)
412 c
413  salb=qscat/qext
414 c write onto the output data sets
415 c
416  goto(170,172,174,176,178),ifunc
417 170 continue
418  write(8,268)title
419  write(8,249)ifunc,mfunc
420  write(8,251)xlamb,rfr1,rfi1
421  write(8,252)ddd(10),ddd(11)
422  write(8,253)xxz1,yyz1
423  write(8,270)ddd(1),ddd(2),ddd(3)
424  write(7,268)title
425  write(7,249)ifunc,mfunc
426  write(7,251)xlamb,rfr1,rfi1
427  write(7,252)ddd(10),ddd(11)
428  write(7,253)xxz1,yyz1
429  write(7,270)ddd(1),ddd(2),ddd(3)
430  goto 180
431 172 continue
432  write(8,268)title
433  write(8,249)ifunc,mfunc
434  write(8,251)xlamb,rfr1,rfi1
435  write(8,252)ddd(10),ddd(11)
436  write(8,253)xxz1,yyz1
437  write(8,272)ddd(1),ddd(2),ddd(3),ddd(4)
438  write(7,268)title
439  write(7,249)ifunc,mfunc
440  write(7,251)xlamb,rfr1,rfi1
441  write(7,252)ddd(10),ddd(11)
442  write(7,253)xxz1,yyz1
443  write(7,272)ddd(1),ddd(2),ddd(3),ddd(4)
444  goto 180
445 174 continue
446  write(8,268)title
447  write(8,249)ifunc,mfunc
448  write(8,251)xlamb,rfr1,rfi1
449  write(8,252)ddd(10),ddd(11)
450  write(8,253)xxz1,yyz1
451  if(irgm.eq.1)then
452  write(8,275)ddd(1),ddd(2),ddd(3)
453  else
454  write(8,274)ddd(1),ddd(2),ddd(3)
455  endif
456  write(7,268)title
457  write(7,249)ifunc,mfunc
458  write(7,251)xlamb,rfr1,rfi1
459  write(7,252)ddd(10),ddd(11)
460  write(7,253)xxz1,yyz1
461  if(irgm.eq.1)then
462  write(7,275)ddd(1),ddd(2),ddd(3)
463  else
464  write(7,274)ddd(1),ddd(2),ddd(3)
465  endif
466  goto 180
467 176 continue
468  write(8,268)title
469  write(8,249)ifunc,mfunc
470  write(8,254)xlamb,rfr1,rfi1,rfr2,rfi2
471  write(8,255)rmin1,rmax1,rmin2,rmax2
472  write(8,256)xxz1,yyz1,xxz2,yyz2
473  if(irgm.eq.1)then
474  write(8,277)ddd(1),ddd(2),ddd(3),ddd(4),ddd(5),ddd(6)
475  else
476  write(8,276)ddd(1),ddd(2),ddd(3),ddd(4),ddd(5),ddd(6)
477  endif
478  write(7,268)title
479  write(7,249)ifunc,mfunc
480  write(7,254)xlamb,rfr1,rfi1,rfr2,rfi2
481  write(7,255)rmin1,rmax1,rmin2,rmax2
482  write(7,256)xxz1,yyz1,xxz2,yyz2
483  if(irgm.eq.1)then
484  write(7,277)ddd(1),ddd(2),ddd(3),ddd(4),ddd(5),ddd(6)
485  else
486  write(7,276)ddd(1),ddd(2),ddd(3),ddd(4),ddd(5),ddd(6)
487  endif
488  goto 180
489 178 continue
490  write(8,268)title
491  write(8,249)ifunc,mfunc
492  write(8,257)xlamb,rfr1,rfi1,rfr2,rfi2,rfr3,rfi3
493  write(8,258)rmin1,rmax1,rmin2,rmax2,rmin3,rmax3
494  write(8,259)xxz1,yyz1,xxz2,yyz2,xxz3,yyz3
495  if(irgm.eq.1)then
496  write(8,279)ddd(1),ddd(2),ddd(3),ddd(4),ddd(5),ddd(6),
497  1 ddd(7),ddd(8),ddd(9)
498  else
499  write(8,278)ddd(1),ddd(2),ddd(3),ddd(4),ddd(5),ddd(6),
500  1 ddd(7),ddd(8),ddd(9)
501  endif
502  write(7,268)title
503  write(7,249)ifunc,mfunc
504  write(7,257)xlamb,rfr1,rfi1,rfr2,rfi2,rfr3,rfi3
505  write(7,258)rmin1,rmax1,rmin2,rmax2,rmin3,rmax3
506  write(7,259)xxz1,yyz1,xxz2,yyz2,xxz3,yyz3
507  if(irgm.eq.1)then
508  write(7,279)ddd(1),ddd(2),ddd(3),ddd(4),ddd(5),ddd(6),
509  1 ddd(7),ddd(8),ddd(9)
510  else
511  write(7,278)ddd(1),ddd(2),ddd(3),ddd(4),ddd(5),ddd(6),
512  1 ddd(7),ddd(8),ddd(9)
513  endif
514 180 continue
515  write(8,281)r11,r22,r33,r44,reff,veff
516  write(8,280)ccnsml,bsr,salb,asf,qscat,qext
517  write(7,281)r11,r22,r33,r44,reff,veff
518  write(7,280)ccnsml,bsr,salb,asf,qscat,qext
519  ormin(1) = rmin1
520  ormin(2) = rmin2
521  ormin(3) = rmin3
522  ormax(1) = rmax1
523  ormax(2) = rmax2
524  ormax(3) = rmax3
525  odelr(1) = xxz1
526  odelr(2) = xxz2
527  odelr(3) = xxz3
528  odelx(1) = yyz1
529  odelx(2) = yyz2
530  odelx(3) = yyz3
531  or11 = r11
532  or22 = r22
533  or33 = r33
534  or44 = r44
535  oreff = reff
536  oveff = veff
537  occnsml = ccnsml
538  obsr = bsr
539  osalb = salb
540  oasf = asf
541  oqscat = qscat
542  oqext = qext
543  write(8,282)
544  do 50 i=1,1801
545  write(8,49)angl(i),phfu(i),pol(i)
546  write(7,47) (t(i,k),k=1,4),angl(i)
547  oangl(i) = angl(i)
548  ophfu(i) = phfu(i)
549  opol(i) = pol(i)
550  othd(i) = angl(i)
551  do 51 k=1,4
552  ot(i,k) = t(i,k)
553 51 continue
554 50 continue
555 c***********************************************************************
556 c.....format statements.......
557 c
558  47 format(1p,2d18.12,1p,2d19.12,0p,f6.2)
559  49 format(f5.1,1p,2e15.5)
560 249 format(t2,'ifunc',t16,'ifunc_sub'/i5,11x,i5)
561 251 format(t2,'wavelength',t16,'refr',t27,'refi'/1p3e11.3)
562 252 format(t5,'rmin',t16,'rmax',t24/1p2e11.3)
563 253 format(t4,'deltar',t15,'deltax'/1p2e11.3)
564 cbaf254 format(t2,'wavelength',t16,'refr1',t27,'refi1',t38,'refr2',
565 cbaf 1 t49,'refi2'/1p5e11.3)
566 254 format(t2,'wavelength',t16,'refr1',t27,'refi1',t38,'refr2',
567  1 t49,'refi2'/1p5e11.4)
568 255 format(t5,'rmin1',t16,'rmax1',t27,'rmin2',t38,'rmax2'/
569  1 1p4e11.3)
570 256 format(t4,'deltar1',t15,'deltax1',t26,'deltar2',
571  1 t37,'deltax2'/1p4e11.3)
572 257 format(t2,'wavelength',t15,'refr1',t26,'refi1',t27,'refr2',
573  1 t38,'refi2',t49,'refr3',t60,'refi3'/1p7e11.3)
574 258 format(t5,'rmin1',t16,'rmax1',t27,'rmin2',t38,'rmax2't49,'rmin3',
575  1 t60,'rmax3',t68/1p6e11.3)
576 259 format(t4,'deltar1',t15,'deltax1',t26,'deltar2',t37,'deltax2',
577  1 t48,'deltar3',t59,'deltax3'/1p6e11.3)
578 268 format(a)
579 271 format(t2,'wavelength',t16,'refr',t27,'refi't38,'rmin',
580  1 t49,'rmax',t51,'deltar',t62,'deltax'/1p7e11.3)
581 270 format(t5,'rc',t16,'c',t27,'nu'/1p3e11.3)
582 272 format(t5,'aa',t16,'b',t27,'alfa',t38,'gama'/1p4e11.3)
583 274 format(t5,'rg1',t16,'sig1',t27,'num1'/1p3e11.3)
584 275 format(t5,'rg1',t15,'sig1(ln)',t27,'num1'/1p3e11.3)
585 276 format(t5,'rg1't16,'sig1',t27,'num1',t38,'rg2',
586  1 t49,'sig2',t61,'num2'/1p6e11.3)
587 277 format(t5,'rg1't15,'sig1(ln)',t27,'num1',t38,'rg2',
588  1 t48,'sig2(ln)',t61,'num2'/1p6e11.3)
589 278 format(t3,'rg1',t11,'sig1',t19,'num1',t27,'rg2',t35,'sig2',
590  1 t43,'num2',t51,'rg3',t59,'sig3',t67,'num3'/9f8.5)
591 279 format(t3,'rg1',t10,'sig1(ln)',t19,'num1',t27,'rg2',
592  1 t34,'sig2(ln)',t43,'num2',t51,'rg3',t58,'sig3(ln)',
593  2 t67,'num3'/9f8.5)
594 280 format(t4,'ccnsml',t16,'bsr',t27,'salb',t38,'asf',
595  1 t49,'qscat',t61,'qext'/1p4e11.3,1p2e12.4)
596 281 format(t5,'r11',t16,'r22',t27,'r33',t38,'r44',
597  1 t49,'reff't60,'veff'/1p6e11.3)
598 282 format(t3,'angle',t11,'phs_func',t26,'deg_pol')
599 c***********************************************************************
600  return
601  end
602 c***********************************************************************
603 c
604  subroutine norz(ifunc,nn,il,isd,irgm,title,xlamb,xn1,xn2,xnpar,xdx,
605  1 orbar,odnzp,odnp,odvp,osumnp,owt)
606 c
607 c this program is designed for five diffrent size distributions
608 c other size distributions can be added easily by modifying
609 c appropriate portion of the code
610 c
611 c ifunc=1 (junge), =2(deirmendjian haze: l,m,c1 etc.),=3(lognormal)
612 c =4 (bimodal- in log(r) ), =5(tri-modal- in log(r)
613 c***********************************************************************
614  implicit real *8 (a-h,o-z)
615 c
616  include 'afrt_phs.cmn'
617 c
618 c***********************************************************************
619 c
620  real*8 xlamb,xn1(nmd,nsdt,nw),xn2(nmd,nsdt,nw),
621  1 xdx(nmd,nsdt,nw),xnpar(nmd,nsdt,nw),
622  2 orbar(ndr),odnzp(ndr),odnp(ndr),
623  3 odvp(ndr),owt(ndr),osumnp(ndr)
624 
625  character*80 title
626 c
627 c note buffer ddd(15) contains the cut-off value for rmax. ddd(15)
628 c is set in the main program
629 c title=xtitle(isd)
630 c xlamb=xww(isd,il)
631  kzia=1
632  goto(25,30,35,40,45),ifunc
633 25 continue
634 c refr1=xn1(1,isd,il)
635 c refi1=xn2(1,isd,il)
636  rmin=ddd(4)
637  rmax=ddd(5)
638  rc=ddd(1)
639  c=ddd(2)
640  xnu=ddd(3)
641  if(rmax.gt.ddd(15))rmax=ddd(15)
642  ddd(10)=rmin
643  ddd(11)=rmax
644  rmin1=rmin
645  rmax1=rmax
646  write(4,700)title
647  write(4,710)xlamb,ddd(10),ddd(11)
648  write(4,712)refr1,refi1
649  write(4,720)ddd(1),ddd(2),ddd(3)
650  goto 50
651 30 continue
652  refr1=xn1(1,isd,il)
653  refi1=xn2(1,isd,il)
654  a=ddd(1)
655  b=ddd(2)
656  alfa=ddd(3)
657  gama=ddd(4)
658  rmin=ddd(5)
659  rmax=ddd(6)
660  if(rmax.gt.ddd(15))rmax=ddd(15)
661  ddd(10)=rmin
662  ddd(11)=rmax
663  rmin1=rmin
664  rmax1=rmax
665  write(4,700)title
666  write(4,710)xlamb,ddd(10),ddd(11)
667  write(4,712)refr1,refi1
668  write(4,730)ddd(1),ddd(2),ddd(3),ddd(4)
669  goto 50
670 35 continue
671  refr1=xn1(1,isd,il)
672  refi1=xn2(1,isd,il)
673  rm1=ddd(1)
674  if(irgm.eq.1)then
675  sig11=ddd(2)
676  else
677  sig11=dlog(ddd(2))
678  endif
679  par1=ddd(3)
680  rmin1=dexp(dlog(rm1)-4.0d0*sig11)
681  if(rmin1 .lt. 0.00d0)rmin1=0.00d0
682  rmax1=dexp(dlog(rm1)+4.0d0*sig11)
683  if(rmax1.gt.ddd(15))rmax1=ddd(15)
684  ddd(10)=rmin1
685  ddd(11)=rmax1
686  rmin=rmin1
687  rmax=rmax1
688  write(4,700)title
689  write(4,710)xlamb,ddd(10),ddd(11)
690  write(4,712)refr1,refi1
691  if(irgm.eq.1)then
692  write(4,741)ddd(1),ddd(2),ddd(3)
693  else
694  write(4,740)ddd(1),ddd(2),ddd(3)
695  endif
696  goto 50
697 40 continue
698  refr1=xn1(1,isd,il)
699  refi1=xn2(1,isd,il)
700  refr2=xn1(2,isd,il)
701  refi2=xn2(2,isd,il)
702  rm1=ddd(1)
703  par1=ddd(3)
704  rm2=ddd(4)
705  par2=ddd(6)
706  if(irgm.eq.1)then
707  sig11=ddd(2)
708  sig21=ddd(5)
709  else
710  sig11=dlog(ddd(2))
711  sig21=dlog(ddd(5))
712  endif
713  rmin1=dexp(dlog(rm1)-4.0d0*sig11)
714  rmax1=dexp(dlog(rm1)+4.0d0*sig11)
715  if(rmin1 .lt. 0.0d0)rmin=0.0d0
716  rmin2=dexp(dlog(rm2)-4.0d0*sig21)
717  rmax2=dexp(dlog(rm2)+4.0d0*sig21)
718  if(rmax1.gt.ddd(15))rmax1=ddd(15)
719  if(rmax2.gt.ddd(15))rmax2=ddd(15)
720  rmin=rmin1
721  rmax=rmax2
722  ddd(10)=rmin
723  ddd(11)=rmax
724  write(4,700)title
725  write(4,710)xlamb,ddd(10),ddd(11)
726  write(4,714)refr1,refi1,refr2,refi2
727  if(irgm.eq.1)then
728  write(4,751)ddd(1),ddd(2),ddd(3),ddd(4),ddd(5),ddd(6)
729  else
730  write(4,750)ddd(1),ddd(2),ddd(3),ddd(4),ddd(5),ddd(6)
731  endif
732  goto 50
733  45 continue
734  refr1=xn1(1,isd,il)
735  refi1=xn2(1,isd,il)
736  refr2=xn1(2,isd,il)
737  refi2=xn2(2,isd,il)
738  refr3=xn1(3,isd,il)
739  refi3=xn2(3,isd,il)
740  rm1=ddd(1)
741  par1=ddd(3)
742  rm2=ddd(4)
743  par2=ddd(6)
744  rm3=ddd(7)
745  par3=ddd(9)
746  if(irgm.eq.1)then
747  sig11=ddd(2)
748  sig21=ddd(5)
749  sig31=ddd(8)
750  else
751  sig11=dlog(ddd(2))
752  sig21=dlog(ddd(5))
753  sig31=dlog(ddd(8))
754  endif
755  rmin1=dexp(dlog(rm1)-4.0d0*sig11)
756  rmax1=dexp(dlog(rm1)+4.0d0*sig11)
757  if(rmin1 .lt. 0.0d0)rmin=0.0d0
758  rmin2=dexp(dlog(rm2)-4.0d0*sig21)
759  rmax2=dexp(dlog(rm2)+4.0d0*sig21)
760  rmin3=dexp(dlog(rm3)-4.0d0*sig31)
761  rmax3=dexp(dlog(rm3)+4.0d0*sig31)
762  if(rmax1.gt.ddd(15))rmax1=ddd(15)
763  if(rmax2.gt.ddd(15))rmax2=ddd(15)
764  if(rmax3.gt.ddd(15))rmax3=ddd(15)
765  rmin=rmin1
766  rmax=rmax3
767  ddd(10)=rmin
768  ddd(11)=rmax
769  write(4,700)title
770  write(4,710)xlamb,ddd(10),ddd(11)
771  write(4,716)refr1,refi1,refr2,refi2,refr3,refi3
772  if(irgm.eq.1)then
773  write(4,760)ddd(1),ddd(2),ddd(3),ddd(4),ddd(5),ddd(6),
774  1 ddd(7),ddd(8),ddd(9)
775  else
776  write(4,760)ddd(1),ddd(2),ddd(3),ddd(4),ddd(5),ddd(6),
777  1 ddd(7),ddd(8),ddd(9)
778  endif
779  goto 50
780 50 continue
781  r(1)=rmin
782  if(ifunc .gt. 1)goto 53
783  delr=(xlamb*xdx(1,isd,il))/(2.0d0*pi)
784  n1=(rc-rmin)/delr+1.0d0
785  n2=(rmax-rc)/delr+1.0d0
786  xx1=(rc-rmin)/dfloat(n1)
787  xx2=(rmax-rc)/dfloat(n2)
788  do 51 i=1,n1
789  r(i+1)=r(i)+xx1
790  deltar(i)=xx1
791 51 continue
792  do 52 i=1,n2
793  r(i+1+n1)=r(i+n1)+xx2
794  deltar(n1+i)=xx2
795 52 continue
796  r(n1+n2+1)=rmax
797  nn=n1+n2
798  goto 57
799 53 continue
800  if(ifunc.le.3)then
801  delr=(xlamb*xdx(1,isd,il))/(2.0d0*pi)
802  nn=(rmax-rmin)/delr+1.0d0
803  xxx=(rmax-rmin)/dfloat(nn)
804  endif
805 c
806  if(ifunc.eq.4)then
807  delr1=(xlamb*xdx(1,isd,il))/(2.0d0*pi)
808  delr2=(xlamb*xdx(2,isd,il))/(2.0d0*pi)
809  nn1=(rmax1-rmin1)/delr1+1.0d0
810  xxx1=(rmax1-rmin1)/dfloat(nn1)
811  nn2=(rmax2-rmax1)/delr2+1.0d0
812  xxx2=(rmax2-rmax1)/dfloat(nn2)
813  nn=nn1+nn2
814  if (nn .gt. ndr) then
815  write(*,*) 'Increase ndr dimension in phs common, or decrease delx',nn,ndr
816  call exit(1)
817  endif
818  endif
819  if(ifunc.eq.5)then
820  delr1=(xlamb*xdx(1,isd,il))/(2.0d0*pi)
821  delr2=(xlamb*xdx(2,isd,il))/(2.0d0*pi)
822  delr3=(xlamb*xdx(3,isd,il))/(2.0d0*pi)
823  nn1=(rmax1-rmin1)/delr1+1.0d0
824  xxx1=(rmax1-rmin1)/dfloat(nn1)
825  nn2=(rmax2-rmax1)/delr2+1.0d0
826  xxx2=(rmax2-rmax1)/dfloat(nn2)
827  nn3=(rmax3-rmax2)/delr3+1.0d0
828  xxx3=(rmax3-rmax2)/dfloat(nn3)
829  nn=nn1+nn2+nn3
830  endif
831  do 55 i=1,nn
832  if(ifunc.le.3)then
833  kk=i
834  deltar(i)=xxx
835  r(i+1)=r(i)+xxx
836  if(r(i+1).gt.ddd(15))goto 56
837  endif
838  if(ifunc.eq.4)then
839  kk=i
840  xxx=xxx1
841  if(i.gt.nn1)xxx=xxx2
842  deltar(i)=xxx
843  r(i+1)=r(i)+xxx
844  if(r(i+1).gt.ddd(15))goto 56
845  endif
846  if(ifunc.eq.5)then
847  kk=i
848  xxx=xxx1
849  if(i.gt.nn1)xxx=xxx2
850  if(i.gt.nn2)xxx=xxx3
851  deltar(i)=xxx
852  r(i+1)=r(i)+xxx
853  if(r(i+1).gt.ddd(15))goto 56
854  endif
855 55 continue
856 56 continue
857  nn=kk
858  r(nn+1)=rmax
859 57 continue
860 c
861  goto(60,65,70,75,77),ifunc
862 60 continue
863  write(4,800)
864  sumn1(1)=0.0d0
865  do 64 i=1,nn
866  rbar(i)=(r(i+1)+r(i) )/2.0d0
867  xxx=r(i+1)-r(i)
868  if(r(i+1) .gt. rc)goto 62
869  dn1(i)=c
870  dnz1(i)=c*rbar(i)
871  wt(i)=dn1(i)*xxx
872  sumn1(i+1)=sumn1(i)+wt(i)
873  goto 63
874 62 continue
875  dn1(i)=c*(rc/rbar(i))**(xnu+1)
876  dnz1(i)=c*rbar(i)*(rc/rbar(i))**(xnu+1)
877  wt(i)=dn1(i)*xxx
878  sumn1(i+1)=sumn1(i)+wt(i)
879 63 continue
880  write(4,118)i,rbar(i),dnz1(i),dn1(i),wt(i),sumn1(i)
881  orbar(i) = rbar(i)
882  odnzp(i) = dnz1(i)
883  odnp(i) = dn1(i)
884  owt(i) = wt(i)
885  osumnp(i) = sumn1(i)
886 64 continue
887  goto 80
888 65 continue
889  write(4,810)
890  sumn1(1)=0.0d0
891  do 66 i=1,nn
892  rbar(i)=(r(i+1)+r(i) )/2.0d0
893  xxx=r(i+1)-r(i)
894  dn1(i)=a*rbar(i)**alfa*dexp(-b*rbar(i)**gama )
895  dnz1(i)=rbar(i)*dn1(i)
896  wt(i)=dn1(i)*xxx
897  sumn1(i+1)=sumn1(i)+wt(i)
898  write(4,118)i,rbar(i),dnz1(i),dn1(i),wt(i),sumn1(i)
899  orbar(i) = rbar(i)
900  odnzp(i) = dnz1(i)
901  odnp(i) = dn1(i)
902  owt(i) = wt(i)
903  osumnp(i) = sumn1(i)
904 66 continue
905  goto 80
906 70 continue
907  write(4,820)
908  sumn1(1)=0.0d0
909  do 71 i=1,nn
910  rbar(i)=(r(i+1)+r(i) )/2.0d0
911  xxx=r(i+1)-r(i)
912  xlgr=dlog(rbar(i))
913  xlgrm1=dlog(rm1)
914  if(irgm.eq.1)then
915  xsig11=ddd(2)
916  else
917  xsig11=dlog(ddd(2))
918  endif
919  par1=xnpar(1,isd,il)
920  consr1=par1*1.0d0/(dsqrt(2.0d0*pi)*xsig11)
921  dnz1(i)=consr1*dexp(-0.5d0*((xlgr-xlgrm1)/xsig11)**2)
922  dv1(i)=(4.0d0/3.0d0)*pi*rbar(i)**3*dnz1(i)
923  dn1(i)=(dnz1(i)/rbar(i))
924  wt(i)=dn1(i)*xxx
925  sumn1(i+1)=sumn1(i)+wt(i)
926 71 continue
927  do i=1,nn
928  write(4,118)i,rbar(i),dnz1(i),dn1(i),dv1(i),wt(i),sumn1(i)
929  orbar(i) = rbar(i)
930  odnzp(i) = dnz1(i)
931  odnp(i) = dn1(i)
932  odvp(i) = dv1(i)
933  owt(i) = wt(i)
934  osumnp(i) = sumn1(i)
935  enddo
936  goto 80
937 75 continue
938  write(4,830)
939  sumn1(1)=0.0d0
940  sumn2(1)=0.0d0
941  sumnp(1)=0.0d0
942  par1=xnpar(1,isd,il)
943  par2=xnpar(2,isd,il)
944  do 76 i=1,nn
945  rbar(i)=(r(i+1)+r(i) )/2.0d0
946  xxx=r(i+1)-r(i)
947  xlgr=dlog(rbar(i))
948  xlgrm1=dlog(rm1)
949  xlgrm2=dlog(rm2)
950  if(irgm.eq.1)then
951  xsig11=ddd(2)
952  xsig21=ddd(5)
953  else
954  xsig11=dlog(ddd(2))
955  xsig21=dlog(ddd(5))
956  endif
957  consr1=1.0d0/(dsqrt(2.0d0*pi)*xsig11)
958  consr2=1.0d0/(dsqrt(2.0d0*pi)*xsig21)
959  dnz1(i)=consr1*dexp(-0.5d0*((xlgr-xlgrm1)/xsig11)**2)
960  dnz2(i)=consr2*dexp(-0.5d0*((xlgr-xlgrm2)/xsig21)**2)
961  dn1(i)=dnz1(i)/rbar(i)
962  dn2(i)=dnz2(i)/rbar(i)
963  dv1(i)=(4.0d0/3.0d0)*pi*rbar(i)**3*dnz1(i)
964  dv2(i)=(4.0d0/3.0d0)*pi*rbar(i)**3*dnz2(i)
965  dn1(i)=par1*dn1(i)
966  dn2(i)=par2*dn2(i)
967  dnp(i)=dn1(i)+dn2(i)
968  dnz1(i)=par1*dnz1(i)
969  dnz2(i)=par2*dnz2(i)
970  dnzp(i)=dnz1(i)+dnz2(i)
971  dvp(i)=par1*dv1(i)+par2*dv2(i)
972  wt(i)=dn1(i)*xxx+dn2(i)*xxx
973  sumn1(i+1)=sumn1(i)+dn1(i)*xxx
974  sumn2(i+1)=sumn2(i)+dn2(i)*xxx
975  sumnp(i+1)=sumn1(i+1)+sumn2(i+1)
976  write(4,121)i,rbar(i),dnzp(i),dnp(i),dvp(i),wt(i),sumnp(i)
977  orbar(i) = rbar(i)
978  odnzp(i) = dnzp(i)
979  odnp(i) = dnp(i)
980  odvp(i) = dvp(i)
981  owt(i) = wt(i)
982  osumnp(i) = sumnp(i)
983 76 continue
984  goto 80
985 77 continue
986  write(4,840)
987  sumn1(1)=0.0d0
988  sumn2(1)=0.0d0
989  sumn3(1)=0.0d0
990  sumnp(1)=0.0d0
991  par1=xnpar(1,isd,il)
992  par2=xnpar(2,isd,il)
993  par3=xnpar(3,isd,il)
994  do 78 i=1,nn
995  rbar(i)=(r(i+1)+r(i) )/2.0d0
996  xxx=r(i+1)-r(i)
997  xlgr=dlog(rbar(i))
998  xlgrm1=dlog(rm1)
999  xlgrm2=dlog(rm2)
1000  xlgrm3=dlog(rm3)
1001  if(irgm.eq.1)then
1002  xsig11=ddd(2)
1003  xsig21=ddd(5)
1004  xsig31=ddd(8)
1005  else
1006  xsig11=dlog(ddd(2))
1007  xsig21=dlog(ddd(5))
1008  xsig31=dlog(ddd(8))
1009  endif
1010  consr1=par1/(dsqrt(2.0d0*pi)*xsig11)
1011  consr2=par2/(dsqrt(2.0d0*pi)*xsig21)
1012  consr3=par3/(dsqrt(2.0d0*pi)*xsig31)
1013  dnz1(i)=consr1*dexp(-0.5d0*((xlgr-xlgrm1)/xsig11)**2)
1014  dnz2(i)=consr2*dexp(-0.5d0*((xlgr-xlgrm2)/xsig21)**2)
1015  dnz3(i)=consr3*dexp(-0.5d0*((xlgr-xlgrm3)/xsig31)**2)
1016  dnzp(i)=dnz1(i)+dnz2(i)+dnz3(i)
1017  dv1(i)=par1*(4.0d0/3.0d0)*pi*rbar(i)**3*dnz1(i)
1018  dv2(i)=par2*(4.0d0/3.0d0)*pi*rbar(i)**3*dnz2(i)
1019  dv3(i)=par3*(4.0d0/3.0d0)*pi*rbar(i)**3*dnz3(i)
1020  dvp(i)=dv1(i)+dv2(i)+dv3(i)
1021  dn1(i)=(dnz1(i)/rbar(i))
1022  dn2(i)=(dnz2(i)/rbar(i))
1023  dn3(i)=(dnz3(i)/rbar(i))
1024  dnp(i)=dn1(i)+dn2(i)+dn3(i)
1025  wt(i)=dn1(i)*xxx+dn2(i)*xxx+dn3(i)*xxx
1026  sumn1(i+1)=sumn1(i)+dn1(i)*xxx
1027  sumn2(i+1)=sumn2(i)+dn2(i)*xxx
1028  sumn3(i+1)=sumn3(i)+dn3(i)*xxx
1029  sumnp(i+1)=sumn1(i+1)+sumn2(i+1)+sumn3(i+1)
1030  write(4,123)i,rbar(i),dnzp(i),dnp(i),dvp(i),wt(i),sumnp(i)
1031  orbar(i) = rbar(i)
1032  odnzp(i) = dnzp(i)
1033  odnp(i) = dnp(i)
1034  odvp(i) = dvp(i)
1035  owt(i) = wt(i)
1036  osumnp(i) = sumnp(i)
1037 78 continue
1038 80 continue
1039 c*****format statements*************************************************
1040 c
1041 118 format(i4,1x,1p6e11.3)
1042 121 format(i4,1x,1p6e11.3)
1043 123 format(i4,1x,1p6e11.3)
1044 700 format(a)
1045 710 format(t4,'lambda',t16,'rmin',t27,'rmax'/1p3e11.3)
1046 712 format(t5,'refr1',t16,'refi1'/1p2e11.3)
1047 714 format(t5,'refr1',t16,'refi1',t27,'refr2',t38,'refi2'/
1048  1 1p2e11.3,1x,1p2e11.3)
1049 716 format(t5,'refr1',t16,'refi1',t27,'refr2',t38,'refi2',t49,'refr3',
1050  1 t60,'refi3'/1p2e11.3,1x,1p2e11.3,1x,1p2e11.3)
1051 720 format(t5,'rc',t16,'c',t27,'nu_star'/1p3e11.3)
1052 730 format(t5,'aa',t16,'b',t27,'alfa',t38,'gama'/1p4e11.3)
1053 740 format(t5,'rg1',t16,'sig1',t27,'num1'/1p3e11.3)
1054 741 format(t5,'rg1',t15,'sig1(ln)',t27,'num1'/1p3e11.3)
1055 750 format(t5,'rg1',t16,'sig1',t27,'num1',t38,'rg2',t49,'sig2',t60,
1056  1 'num2'/1p6e11.3)
1057 751 format(t5,'rg1',t15,'sig1(ln)',t27,'num1',t38,'rg2',t48,
1058  1 'sig2(ln)',t60,'num2'/1p6e11.3)
1059 760 format(t4,'rg1',t12,'sig1',t20,'num1',t28,'rg2',t36,'sig2',
1060  1 t44,'num2',t52,'rg3',t60,'sig3',t68,'num3'/9f8.5)
1061 761 format(t4,'rg1',t11,'sig1(ln)',t20,'num1',t28,'rg2',t35,
1062  1 'sig2(ln)'t44,'num2',t52,'rg3',t59,'sig3(ln)',
1063  2 t68,'num3'/9f8.5)
1064 800 format(t2,'num',t10,'rbar',t21,'dnz',t32,'dn',t43,'wt',
1065  1 t53,'tot_num')
1066 810 format(t3,'num',t10,'rbar',t19,'dn/dlogr',t32,'dn/dr',t43,'wt',
1067  1 t53,'tot_num')
1068 820 format(t3,'num',t10,'rbar',t19,'dn/dlogr',t32,'dn/dr',
1069  1 t41,'dv/dlogr',t54,'wt',t64,'tot_num')
1070 830 format(t3,'num',t10,'rbar',t19,'dn/dlogr',t32,'dn/dr',
1071  1 t41,'dv/dlogr',t54,'wt',t64,'tot_num')
1072 840 format(t3,'num',t10,'rbar',t19,'dn/dlogr',t32,'dn/dr',
1073  1 t41,'dv/dlogr',t54,'wt',t64,'tot_num')
1074 c***********************************************************************
1075  return
1076  end
1077 c**********************************************************************
1078  subroutine ccn(ifunc,irgm)
1079 c
1080 c subroutine to compute ccn when rg is small (le 0.1)
1081 c
1082 c**********************************************************************
1083 c
1084  implicit real*8 (a-h,o-z)
1085  include 'afrt_phs.cmn'
1086 c***********************************************************************
1087 c write(*,*)'welcome to subroutine ccn'
1088  ccnsml=-777.0d0
1089  r0=0.03d0
1090  rg=ddd(1)
1091  if(irgm.eq.1)then
1092  sigln=ddd(2)
1093  else
1094  sigln=dlog(ddd(2))
1095  endif
1096  if(rg.le.0.1d0)then
1097  aa=(dlog(r0/rg))/(sigln*dsqrt(2.0d0))
1098  erfa1=erfz(aa)
1099  ccnsml=0.5d0*(1.0d0-erfa1)
1100  endif
1101 c
1102  return
1103  end
1104 c***********************************************************************
1105 c
1106  subroutine dbmie (x,rfr,rfi,thetd,jx,qext,qscat,ctbrqs,eltrmx,
1107  1 wat,tt )
1108 c subroutine for computing the parameters of the electromagnetic
1109 c radiation scattered by a sphere. this subroutine carries out all
1110 c computations in single precision arithmetic.
1111 c this subroutine computes the capital a function by making use of
1112 c downward recorrence relationship.
1113 c x 0 size parameter of the sphere,( 2 * pi * radius of the sphere)/
1114 c wavelength of the incident radiation).
1115 c rf0 refractive index of the material of the sphere. complex
1116 c quantity..form0 (rfr - i * rfi )
1117 c thetd(j)0 angle in degrees between the directions of the incident
1118 c and the scattered radiation. thetd(j) is ª or = 90.0.
1119 c if thetd(j) should happen to be greater than 90.0, enter with
1120 c supplementary value0 see comments below on eltrmx..
1121 c jx0 total number of thetd for which the computation arerequirde.
1122 c jx should not exceed 200 unless the dimensions statements
1123 c are appropriately modified.
1124 c main program should also have real thetd(200),eltrmx(4,200,2).
1125 c definitions for the following symbols can be found in ' light
1126 c scattering by small particles, h. c. van de hulst, john wiley +
1127 c sons, inc., new york, 1957 '.
1128 c qext82 effieciency factor for extinction, van de hulst, p.14 + 127
1129 c qscat82 effieciency factor for scattering,van de hulst,p.14 + 127.
1130 c ctbrqs0 average(cosine theta) * qscat,van de hulst, p. 128.
1131 c eltrmx(k,i,j)0 elements of the transformation matrix f,van de hul
1132 c st,p.34,45 + 125. i = 10 element m sub 2..i = 20element m sub 1..
1133 c i = 30 element s sub 21.. i = 40 element d sub 21...
1134 c eltrmx(1,i,j) represents the ith element of the matrix for
1135 c the angle thetd(j).. eltrmx(i,j,2) represents the ith element
1136 c of the matrix for the angle 180.0 - thetd(j) ..
1137 c*************************************************************************
1138  implicit real * 8 (a-h,o-z)
1139  real*8 x,rx,rfr,rfi,qext,qscat,t(5),ta(4),tb(2),tc(2)
1140  real*8 td(2),te(2),ctbrqs
1141  real*8 pip(3,1000)
1142  real*8 eltrmx(4,1000,2),pi(3,1000),tau(3,1000),cstht(1000),
1143  1 si2tht(1000),thetd(1000)
1144  complex*16 rf,rrf,rrfx,wm1,fna,fnb,tc1,tc2,wfn(2),acap(7000)
1145  complex*16 fnap,fnbp
1146 c ta(1)0 real part of wfn(1).. ta(2)0 imaginary part of wfn(1)..
1147 c ta(3)0 real part of wfn(2).. ta(4)0 imaginary part of wfn(2)..
1148 c tb(1)0 real part of fna...tb(2)0 imaginary part of fna...
1149 c tc(1)0 real part of fnb...tc(2)0 imaginary part of fnb...
1150 c td(1)0 real part of fnap.. td(2) imaginary part of fnap...
1151 c te(1)0 real part of fnbp... te(2)0 imaginary part of fnbp...
1152 c fnap + fnbp are the preceding values of fna + fnb respectively.
1153  real * 8 tt(1801,4),den,az,bz,rt1,rt2,rt3,rt4
1154  equivalence(wfn(1), ta(1)), (fna, tb(1)), (fnb, tc(1))
1155  equivalence(fnap, td(1)), (fnbp, te(1))
1156 c*************************************************************************
1157  if ( jx .ge. 1001 ) then
1158  write (6, 7)
1159  write (6, 6)
1160  stop 1
1161  endif
1162  rf=dcmplx(rfr,-rfi)
1163  rrf = 1.0d0/rf
1164  rx = 1.0d0/x
1165  rrfx = rrf * rx
1166  t(1)=(x**2)*(rfr**2+rfi**2)
1167  t(1)=dsqrt(t(1))
1168  nmx1 = 1.10d0 * t(1)
1169  if (nmx1 .ge. 7000) then
1170  write(6, 8)
1171  stop 2
1172  endif
1173  nmx2 = t(1)
1174  if (nmx1 .le. 150) then
1175  nmx1 = 150
1176  nmx2 = 135
1177  endif
1178  acap(nmx1 + 1 ) = ( 0.0d0, 0.0d0 )
1179  do n = 1, nmx1
1180  nn = nmx1-n+1
1181  acap(nn) = (nn+1)*rrfx-1.0d0/((nn+1)*rrfx+acap(nn+1))
1182  enddo
1183  do j = 1, jx
1184  if ( thetd(j) .lt. 0.0d0 ) thetd(j) = dabs(thetd(j))
1185  if ( thetd(j) .le. 0.0d0 ) then
1186  cstht(j) = 1.0d0
1187  si2tht(j) = 0.0d0
1188  else if (thetd(j) .lt. 90.0d0 ) then
1189  t(1) = (3.1415926535897932d0*thetd(j))/180.0d0
1190  cstht(j) = dcos(t(1))
1191  si2tht(j) = 1.0d0 - cstht(j)**2
1192  else if ( thetd(j) .eq. 90.0d0 ) then
1193  cstht(j) = 0.0d0
1194  si2tht(j) = 1.0d0
1195  else
1196  write (6, 5) thetd(j)
1197  write (6, 6)
1198  stop 3
1199  endif
1200  enddo
1201  do j = 1, jx
1202  pi(1,j) = 0.0d0
1203  pi(2,j) = 1.0d0
1204  tau(1,j) = 0.0d0
1205  pip(1,j)=0.
1206  pip(2,j)=0
1207  tau(2, j) = cstht(j)
1208  enddo
1209  t(1) = dcos(x)
1210  t(2) = dsin(x)
1211  wm1 = dcmplx(t(1),-t(2))
1212  wfn(1) = dcmplx(t(2),t(1))
1213  wfn(2) = rx * wfn(1)-wm1
1214  tc1 = acap(1)*rrf+rx
1215  tc2 = acap(1)*rf+rx
1216  fna = (tc1*ta(3)-ta(1))/(tc1*wfn(2)-wfn(1))
1217  fnb = (tc2*ta(3)-ta(1))/(tc2*wfn(2)-wfn(1))
1218  fnap = fna
1219  fnbp = fnb
1220  t(1) = 1.50d0
1221  tb(1) = t(1)*tb(1)
1222  tb(2) = t(1)*tb(2)
1223  tc(1) = t(1)*tc(1)
1224  tc(2) = t(1)*tc(2)
1225  j = 1
1226  do while (j.le.jx)
1227  if(cstht(j).ne.1.) then
1228  eltrmx(1,j,1)=tc(1)*pi(2,j)+tb(1)*tau(2,j)
1229  eltrmx(2,j,1)=tc(2)*pi(2,j)+tb(2)*tau(2,j)
1230  eltrmx(3,j,1)=tb(1)*pi(2,j)+tc(1)*tau(2,j)
1231  eltrmx(4,j,1)=tb(2)*pi(2,j)+tc(2)*tau(2,j)
1232  eltrmx(1,j,2)=tc(1)*pi(2,j)-tb(1)*tau(2,j)
1233  eltrmx(2,j,2)=tc(2)*pi(2,j)-tb(2)*tau(2,j)
1234  eltrmx(3,j,2)=tb(1)*pi(2,j)-tc(1)*tau(2,j)
1235  eltrmx(4,j,2)=tb(2)*pi(2,j)-tc(2)*tau(2,j)
1236  else
1237  eltrmx(1,j,1)=pip(2,j)*(tc(1)-tb(1))-pi(2,j)*tb(1)
1238  eltrmx(2,j,1)=pip(2,j)*(tc(2)-tb(2))-pi(2,j)*tb(2)
1239  eltrmx(3,j,1)=pip(2,j)*(tb(1)-tc(1))-pi(2,j)*tc(1)
1240  eltrmx(4,j,1)=pip(2,j)*(tb(2)-tc(2))-pi(2,j)*tc(2)
1241  eltrmx(1,j,2)=pip(2,j)*(tc(1)+tb(1))-pi(2,j)*tb(1)
1242  eltrmx(2,j,2)=pip(2,j)*(tc(2)+tb(2))-pi(2,j)*tb(2)
1243  eltrmx(3,j,2)=pip(2,j)*(tc(1)+tb(1))-pi(2,j)*tc(1)
1244  eltrmx(4,j,2)=pip(2,j)*(tc(2)+tb(2))-pi(2,j)*tc(2)
1245  endif
1246  j = j + 1
1247  enddo
1248  qext = 2.0d0 * ( tb(1) + tc(1))
1249  qscat =(tb(1)**2 + tb(2)**2 + tc(1)**2 + tc(2)**2)/0.75d0
1250  ctbrqs = 0.0d0
1251  n = 2
1252 65 t(1) = 2 * n - 1
1253  t(2) = n - 1
1254  t(3) = 2 * n + 1
1255  do j=1,jx
1256  pi(3,j)=(t(1)*pi(2,j)*cstht(j)-n*pi(1,j))/t(2)
1257  tau(3,j)=cstht(j)*(pi(3,j)-pi(1,j))-t(1)*si2tht(j)*pi(2,j)+tau(1,j)
1258  pip(3,j)=t(1)*pi(2,j)+pip(1,j)
1259  enddo
1260  wm1 = wfn(1)
1261  wfn(1) = wfn(2)
1262  wfn(2) = t(1)*rx*wfn(1) - wm1
1263  tc1 = acap(n)*rrf+n*rx
1264  tc2 = acap(n)*rf+n*rx
1265  fna = (tc1*ta(3)-ta(1))/(tc1*wfn(2)-wfn(1))
1266  fnb = (tc2*ta(3)-ta(1))/(tc2*wfn(2)-wfn(1))
1267  t(5) = n
1268  t(4) = t(1)/(t(5)*t(2))
1269  t(2) = (t(2)*(t(5) + 1.0d0))/t(5)
1270  ctbrqs = ctbrqs+t(2)*(td(1)*tb(1)+td(2)*tb(2)+te(1)*
1271  1 tc(1)+te(2)*tc(2))+t(4)*(td(1)*te(1)+td(2)*te(2))
1272  qext = qext+t(3)*(tb(1)+tc(1))
1273  t(4) = tb(1)**2 + tb(2)**2 + tc(1)**2 + tc(2)**2
1274  qscat = qscat+t(3)*t(4)
1275  t(2) = n*(n+1)
1276  t(1) = t(3)/t(2)
1277  k = (n/2)*2
1278  do j = 1, jx
1279  if(cstht(j).ne.1.) then
1280  eltrmx(1,j,1)=eltrmx(1,j,1)+t(1)*(tc(1)*pi(3,j)+tb(1)*tau(3,j))
1281  eltrmx(2,j,1)=eltrmx(2,j,1)+t(1)*(tc(2)*pi(3,j)+tb(2)*tau(3,j))
1282  eltrmx(3,j,1)=eltrmx(3,j,1)+t(1)*(tb(1)*pi(3,j)+tc(1)*tau(3,j))
1283  eltrmx(4,j,1)=eltrmx(4,j,1)+t(1)*(tb(2)*pi(3,j)+tc(2)*tau(3,j))
1284  az=1.0d0
1285  if(k.eq.n)az=-az
1286  eltrmx(1,j,2)=eltrmx(1,j,2)+t(1)*az*(tc(1)*pi(3,j)-tb(1)*tau(3,j))
1287  eltrmx(2,j,2)=eltrmx(2,j,2)+t(1)*az*(tc(2)*pi(3,j)-tb(2)*tau(3,j))
1288  eltrmx(3,j,2)=eltrmx(3,j,2)+t(1)*az*(tb(1)*pi(3,j)-tc(1)*tau(3,j))
1289  eltrmx(4,j,2)=eltrmx(4,j,2)+t(1)*az*(tb(2)*pi(3,j)-tc(2)*tau(3,j))
1290  else
1291  eltrmx(1,j,1)=eltrmx(1,j,1)+t(1)*(pip(3,j)*(tc(1)-tb(1))-pi(3,j)*
1292  1 tb(1))
1293  eltrmx(2,j,1)=eltrmx(2,j,1)+t(1)*(pip(3,j)*(tc(2)-tb(2))-pi(3,j)*
1294  1 tb(2))
1295  eltrmx(3,j,1)=eltrmx(3,j,1)+t(1)*(pip(3,j)*(tb(1)-tc(1))-pi(3,j)*
1296  1 tc(1))
1297  eltrmx(4,j,1)=eltrmx(4,j,1)+t(1)*(pip(3,j)*(tb(2)-tc(2))-pi(3,j)*
1298  1 tc(2))
1299  az=1.0d0
1300  bz=-1.0d0
1301  if(k.ne.n) then
1302  az=-az
1303  bz=-bz
1304  endif
1305  eltrmx(1,j,2)=eltrmx(1,j,2)+t(1)*(az*pip(3,j)*(tc(1)+tb(1))-bz*
1306  1 pi(3,j)*tb(1))
1307  eltrmx(2,j,2)=eltrmx(2,j,2)+t(1)*(az*pip(3,j)*(tc(2)+tb(2))-bz*
1308  1 pi(3,j)*tb(2))
1309  eltrmx(3,j,2)=eltrmx(3,j,2)+t(1)*(az*pip(3,j)*(tc(1)+tb(1))-bz*
1310  1 pi(3,j)*tc(1))
1311  eltrmx(4,j,2)=eltrmx(4,j,2)+t(1)*(az*pip(3,j)*(tc(2)+tb(2))-bz*
1312  1 pi(3,j)*tc(2))
1313  endif
1314  enddo
1315  if( t(4) .ge. 1.0d-14 ) then
1316  n = n + 1
1317  do j = 1, jx
1318  pi(1, j) = pi(2, j)
1319  pi(2, j) = pi(3, j)
1320  tau(1, j) = tau(2, j)
1321  tau(2, j) = tau(3, j)
1322  pip(1,j)=pip(2,j)
1323  pip(2,j)=pip(3,j)
1324  enddo
1325  fnap = fna
1326  fnbp = fnb
1327  if (n .le. nmx2) go to 65
1328  write(6, 8)
1329  stop 4
1330  endif
1331  do j = 1, jx
1332  m=j
1333  den=si2tht(j)
1334  do k=1,2
1335  do i=1,4
1336  t(i)=eltrmx(i,j,k)
1337  enddo
1338  if(cstht(j).ne.1.) then
1339  az=1.
1340  if(k.eq.2) az = -az
1341  rt3=(t(3)-az*cstht(j)*t(1))/den
1342  rt4=(t(4)-az*cstht(j)*t(2))/den
1343  rt1=(t(1)-az*cstht(j)*t(3))/den
1344  rt2=(t(2)-az*cstht(j)*t(4))/den
1345  else
1346  rt3=t(1)
1347  rt4=t(2)
1348  rt1=t(3)
1349  rt2=t(4)
1350  endif
1351  if(k.eq.2) m = 1802-m
1352  if(k.ne.2.or.m.ne.901 ) then
1353  tt(m,1) = tt(m,1)+(rt1**2+rt2**2)*wat
1354  tt(m,2) = tt(m,2)+(rt3**2+rt4**2)*wat
1355  tt(m,3) = tt(m,3)+(rt1*rt3+rt2*rt4)*wat
1356  tt(m,4) = tt(m,4)+(rt1*rt4-rt2*rt3)*wat
1357  endif
1358  enddo
1359  enddo
1360  t(1) = 2.0d0 * rx**2
1361  qext = qext * t(1)
1362  qscat = qscat * t(1)
1363  ctbrqs = 2.0d0 * ctbrqs * t(1)
1364 c
1365 c*****format statements**************************************************
1366 5 format(10x' the value of the scattering angle is greater than 90.0
1367  $ degrees. it is ',e15.4)
1368 6 format(//10x,' please read the comments'//)
1369 7 format(//10x,' the value of the argument jx is greater the 1000')
1370 8 format(//10x,'the upper limit for acap is not enough. suggest get
1371  1detailed output and modify subroutine'//)
1372  return
1373  end
1374 c*************************************************************************
1375  subroutine asym(x,y,f,bsr)
1376 c compute the asym. factor for a given phase function
1377 c x=angle
1378 c y=phase function
1379 c f=asym. factor
1380 c n=number of x-y pairs
1381 c*********************************************************************
1382 c
1383  parameter(m=1801)
1384  implicit real*8 (a-h,o-z)
1385  real*8 x(m),y(m)
1386  real*8 xmu(m),dxmu(m),xy(m),xb(m)
1387 c
1388  pi=dacos(-1.0d0)
1389  conv=pi/180.0d0
1390  mm1=m-1
1391  do i=1,m
1392  xmu(i)=cos(x(i)*conv)
1393  enddo
1394 c
1395  do i=1,m
1396  xy(i)=xmu(i)*y(i)
1397  xb(i)=dacos(xmu(i))*y(i)
1398  if(i.le.10)then
1399  endif
1400  enddo
1401  do i=1,mm1
1402  dxmu(i)=dabs(xmu(i)-xmu(i+1))
1403  enddo
1404 c
1405  sumn=0.0d0
1406  sumd=0.0d0
1407  sumb=0.0d0
1408  do i=1,mm1
1409  sumn=sumn+0.5d0*(xy(i)+xy(i+1))*dxmu(i)
1410  sumd=sumd+0.5d0*(y(i)+y(i+1))*dxmu(i)
1411  sumb=sumb+0.5d0*(xb(i)+xb(i+1))*dxmu(i)
1412  enddo
1413 c
1414  xnum=2.0d0*sumn*pi
1415  xden=2.0d0*sumd*pi
1416  f=xnum/xden
1417  bsr=(2.0d0*sumb)/xden
1418 c
1419 c write(6,100)sumn,sumd,sumb,xnum,xden,f,bsr
1420 100 format('asym sumn,sumd,sumb,xnum,xden,f,bsr',1p7e11.3)
1421 c
1422  return
1423  end
1424 c***********************************************************************
1425 c**************************************************************************
1426  function erfz(x)
1427  implicit real*8 (a-h,o-z)
1428  if(x.lt.0.d0)then
1429  erfz=-gammp(.5d0,x**2)
1430  else
1431  erfz=gammp(.5d0,x**2)
1432  endif
1433  return
1434  end
1435 c
1436  function gammp(a,x)
1437  implicit real*8 (a-h,o-z)
1438  if(x.lt.0.d0.or.a.le.0.d0) stop
1439  if(x.lt.a+1.d0)then
1440  call gser(gamser,a,x,gln)
1441  gammp=gamser
1442  else
1443  call gcf(gammcf,a,x,gln)
1444  gammp=1.d0-gammcf
1445  endif
1446  return
1447  end
1448 c
1449  subroutine gser(gamser,a,x,gln)
1450  implicit real*8 (a-h,o-z)
1451  parameter(itmax=100,eps=3.0d-7)
1452  gln=gammln(a)
1453  if(x.le.0.d0)then
1454  if(x.lt.0.d0) stop
1455  gamser=0.d0
1456  return
1457  endif
1458  ap=a
1459  sum=1.d0/a
1460  del=sum
1461  do 11 n=1,itmax
1462  ap=ap+1.d0
1463  del=del*x/ap
1464  sum=sum+del
1465  if(dabs(del).lt.dabs(sum)*eps)go to 1
1466 11 continue
1467  stop 'a too large, itmax too small'
1468 1 gamser=sum*dexp(-x+a*dlog(x)-gln)
1469  return
1470  end
1471 c
1472  subroutine gcf(gammcf,a,x,gln)
1473  implicit real*8 (a-h,o-z)
1474  parameter(itmax=100,eps=3.0d-7)
1475  gln=gammln(a)
1476  gold=0.d0
1477  a0=1.d0
1478  a1=x
1479  b0=0.d0
1480  b1=1.d0
1481  fac=1.d0
1482  do 11 n=1,itmax
1483  an=dfloat(n)
1484  ana=an-a
1485  a0=(a1+a0*ana)*fac
1486  b0=(b1+b0*ana)*fac
1487  anf=an*fac
1488  a1=x*a0+anf*a1
1489  b1=x*b0+anf*b1
1490  if(a1.ne.0.d0)then
1491  fac=1.d0/a1
1492  g=b1*fac
1493  if(dabs((g-gold)/g).lt.eps)go to 1
1494  gold=g
1495  endif
1496 11 continue
1497  stop 'a too large, itmax too small'
1498 1 continue
1499  gammcf=dexp(-x+a*dlog(x)-gln)*g
1500  return
1501  end
1502 c
1503  function gammln(xx)
1504  implicit real*8 (a-h,o-z)
1505  real*8 cof(6),stp,half,one,fpf,x,tmp,ser
1506  data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
1507  * -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
1508  data half,one,fpf/0.5d0,1.0d0,5.5d0/
1509  x=xx-one
1510  tmp=x+fpf
1511  tmp=(x+half)*log(tmp)-tmp
1512  ser=one
1513  do 11 j=1,6
1514  x=x+one
1515  ser=ser+cof(j)/x
1516 11 continue
1517  gammln=tmp+dlog(stp*ser)
1518  return
1519  end
1520 c***********************************************************************
1521 c
1522  subroutine moment(ifunc,nn)
1523 c
1524 c***********************************************************************
1525 
1526  implicit real *8 (a-h,o-z)
1527 c
1528  include 'afrt_phs.cmn'
1529 c
1530 c***********************************************************************
1531 c
1532 c write(*,*)'welcome to subroutine moment'
1533 c compute effective radius
1534 c
1535  prod1=0.0d0
1536  prod2=0.0d0
1537  prod3=0.0d0
1538  prod4=0.0d0
1539 c write(*,*)'ifunc,nn',ifunc,nn
1540 c write(*,*)'rmin1,rmax1',rmin1,rmax1
1541 c do i=1,10
1542 c write(*,120)rbar(i),dn1(i),deltar(i)
1543 c120 format('rbar,dn1,deltar',1p3e12.4)
1544 c enddo
1545  do i=1,nn
1546  if(rbar(i).gt.rmin1 .and. rbar(i).le.rmax1)then
1547  prod1=rbar(i)*dn1(i)*deltar(i)+prod1
1548  prod2=rbar(i)**2*dn1(i)*deltar(i)+prod2
1549  prod3=rbar(i)**3*dn1(i)*deltar(i)+prod3
1550  prod4=rbar(i)**4*dn1(i)*deltar(i)+prod4
1551  endif
1552  enddo
1553  if(ifunc.eq.4 .or. ifunc.eq.5)then
1554  do i=1,nn
1555  if(rbar(i).gt.rmin2 .and. rbar(i).le.rmax2)then
1556  prod1=rbar(i)*dn2(i)*deltar(i)+prod1
1557  prod2=rbar(i)**2*dn2(i)*deltar(i)+prod2
1558  prod3=rbar(i)**3*dn2(i)*deltar(i)+prod3
1559  prod4=rbar(i)**4*dn2(i)*deltar(i)+prod4
1560  endif
1561  enddo
1562  endif
1563  if(ifunc.le.5)then
1564  do i=1,nn
1565  if(rbar(i).gt.rmin3 .and. rbar(i).le.rmax3)then
1566  prod1=rbar(i)*dn3(i)*deltar(i)+prod1
1567  prod2=rbar(i)**2*dn3(i)*deltar(i)+prod2
1568  prod3=rbar(i)**3*dn3(i)*deltar(i)+prod3
1569  prod4=rbar(i)**4*dn3(i)*deltar(i)+prod4
1570  endif
1571  enddo
1572  endif
1573  reffz=prod3/prod2
1574  veffz=prod4/(reffz**2*prod2)
1575  r11=prod1
1576  r22=prod2
1577  r33=prod3
1578  r44=prod4
1579  reff=reffz
1580  veff=veffz
1581 c write(*,100)prod1,prod2,prod3,prod4,reffz,veffz
1582 100 format('prod1,prod2,prod3,prod4,reffz,veffz'/1p6e11.3)
1583  return
1584  end
1585 c***********************************************************************
1586 c
1587  subroutine xmntst(ifunc,irgm)
1588 c
1589 c subroutine to compute the moments and effective radius for log-
1590 c normal distributions
1591 c
1592 c**********************************************************************
1593 c
1594  implicit real*8 (a-h,o-z)
1595  include 'afrt_phs.cmn'
1596  real*8 q11(3),q22(3),q33(3),q44(3),qnn(3)
1597  real*8 amin(3),amax(3),agm(3),asigln(3)
1598 c***********************************************************************
1599 c compute ccnsml for lognormal unimodal distributions
1600 c define r0=0.03d0 (see tanre's write-up)
1601 c write(*,*)'welcome to subroutine xmntst'
1602  ccnsml=-777.0d0
1603  r0=0.03
1604  if(ifunc.eq.3)then
1605  rg=ddd(1)
1606  if(irgm.eq.1)then
1607  sigln=ddd(2)
1608  else
1609  sigln=dlog(ddd(2))
1610  endif
1611  if(rg.le.0.1d0)then
1612  aa=(dlog(r0/rg))/(sigln*dsqrt(2.0d0))
1613  erfa1=erf(aa)
1614  ccnsml=0.5d0*(1.0d0-erfa1)
1615  endif
1616  endif
1617 c
1618 c compute moments
1619  nmode=ifunc-2
1620  amin(1)=rmin1
1621  amin(2)=rmin2
1622  amin(3)=rmin3
1623  amax(1)=rmax1
1624  amax(2)=rmax2
1625  amax(3)=rmax3
1626  agm(1)=ddd(1)
1627  agm(2)=ddd(4)
1628  agm(3)=ddd(7)
1629  qnn(1)=ddd(3)
1630  qnn(2)=ddd(6)
1631  qnn(3)=ddd(9)
1632  if(irgm.eq.1)then
1633  asigln(1)=ddd(2)
1634  asigln(2)=ddd(5)
1635  asigln(3)=ddd(8)
1636  else
1637  asigln(1)=dlog(ddd(2))
1638  asigln(2)=dlog(ddd(5))
1639  asigln(3)=dlog(ddd(8))
1640  endif
1641  if(nmode.gt.3)return
1642 c
1643  do i=1,nmode,1
1644  q11(i)=agm(i)*dexp(0.5d0*asigln(i)**2)
1645  q22(i)=agm(i)**2*dexp(0.5d0*(2.0d0*asigln(i))**2)
1646  q33(i)=agm(i)**3*dexp(0.5d0*(3.0d0*asigln(i))**2)
1647  q44(i)=agm(i)**4*dexp(0.5d0*(4.0d0*asigln(i))**2)
1648  enddo
1649 c
1650  do i=1,nmode,1
1651  rgp1=dlog(agm(i))+asigln(i)**2
1652  rgp2=dlog(agm(i))+2.0d0*asigln(i)**2
1653  rgp3=dlog(agm(i))+3.0d0*asigln(i)**2
1654  rgp4=dlog(agm(i))+4.0d0*asigln(i)**2
1655 c
1656  erf1=erf((dlog(amax(i))-rgp1)/(asigln(i)*dsqrt(2.0d0)))
1657  erf2=erf((dlog(amax(i))-rgp2)/(asigln(i)*dsqrt(2.0d0)))
1658  erf3=erf((dlog(amax(i))-rgp3)/(asigln(i)*dsqrt(2.0d0)))
1659  erf4=erf((dlog(amax(i))-rgp4)/(asigln(i)*dsqrt(2.0d0)))
1660 c
1661  grf1=erf((dlog(amin(i))-rgp1)/(asigln(i)*dsqrt(2.0d0)))
1662  grf2=erf((dlog(amin(i))-rgp2)/(asigln(i)*dsqrt(2.0d0)))
1663  grf3=erf((dlog(amin(i))-rgp3)/(asigln(i)*dsqrt(2.0d0)))
1664  grf4=erf((dlog(amin(i))-rgp4)/(asigln(i)*dsqrt(2.0d0)))
1665 c
1666  q11(i)=0.5d0*q11(i)*(erf1-grf1)
1667  q22(i)=0.5d0*q22(i)*(erf2-grf2)
1668  q33(i)=0.5d0*q33(i)*(erf3-grf3)
1669  q44(i)=0.5d0*q44(i)*(erf4-grf4)
1670  enddo
1671 c
1672  r11=0.0d0
1673  r22=0.0d0
1674  r33=0.0d0
1675  r44=0.0d0
1676  do i=1,nmode,1
1677  r11=r11+qnn(i)*q11(i)
1678  r22=r22+qnn(i)*q22(i)
1679  r33=r33+qnn(i)*q33(i)
1680  r44=r44+qnn(i)*q44(i)
1681  enddo
1682 
1683 c compute effective radius
1684  reff=r33/r22
1685  veff=r44/(reff**2*r22)
1686 c
1687  write(*,100)r11,r22,r33,r44,reff,veff
1688 100 format('final: r11,r22,r33,r44,reff,veff'/1p6e11.3)
1689 c
1690  return
1691  end
1692 c***********************************************************************
gammln
function gammln(xx)
Definition: ocn.f:1270
index
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame index(0-based) was previously being used the frame number(1-based) should have been used. There were only a few minor changes to the uncertainty index(maximum of 1 digit). 3. Some inefficient arrays(Sigma_RVS_norm_sq) were eliminated and some code lines in Preprocess_L1A_Data were moved into Process_OBCEng_Emiss. There were no changes to the product. Required RAM was reduced by 20 MB. Now
erfz
function erfz(x)
Definition: ocn.f:1181
gser
subroutine gser(gamser, a, x, gln)
Definition: ocn.f:1210
ccn
subroutine ccn(ifunc, irgm)
Definition: phs.f:940
angl
subroutine angl
Definition: angl.f:2
gammp
function gammp(a, x)
Definition: ocn.f:1194
real
#define real
Definition: DbAlgOcean.cpp:26
const
subroutine const(NGAUSS, NMAX, MMAX, P, X, W, AN, ANN, S, SS, NP, EPS)
Definition: ampld.lp.f:924
gcf
subroutine gcf(gammcf, a, x, gln)
Definition: ocn.f:1236
fac
#define fac
Definition: atmocor1_land.c:111
xmntst
subroutine xmntst(ifunc, irgm)
Definition: phs.f:1599
parameter
README for MOD_PR03(V6.1.0) 2. POINTS OF CONTACT it can be either SDP Toolkit or MODIS Packet for Terra input files The orbit validation configuration parameter(LUN 600281) must be either "TRUE" or "FALSE". It needs to be "FALSE" when running in Near Real Time mode
flbfr
subroutine flbfr(il, isd, ifunc)
Definition: phs.f:148
moment
subroutine moment(ifunc, nn)
Definition: phs.f:1534
pfunc
subroutine pfunc(il, isd)
Definition: phs.f:179
rc
float rc(float x, float y)
dbmie
subroutine dbmie(x, rfr, rfi, thetd, jx, qext, qscat, ctbrqs, eltrmx, wat, tt)
Definition: phs.f:969
asym
subroutine asym(x, y, f, bsr)
Definition: phs.f:1228
pi
#define pi
Definition: vincenty.c:23
af_phs_read
subroutine af_phs_read(infile, nlamb, nsd, ilam1, ilam2, isd1, isd2, jfunc, jnorm, jrgm, kfunc, irh, iset, xtitle, xww, xn1, xn2, xrg, xsig, xnpar, xdx, epar)
Definition: afrt_phs.f:19
af_phs_process
subroutine af_phs_process(outdir, ailam, aisd, aisd1, aisd2, jfunc, jnorm, jrgm, kfunc, irh, iset, xtitle, xww, xn1, xn2, xrg, xsig, xnpar, xdx, epar, ormin, ormax, odelr, odelx, or11, or22, or33, or44, oreff, oveff, occnsml, obsr, osalb, oasf, oqscat, oqext, oangl, ophfu, opol, ot, othd, orbar, odnzp, odnp, odvp, osumnp, owt)
Definition: afrt_phs.f:107
f
#define f
Definition: l1_czcs_hdf.c:702
norz
subroutine norz(ifunc, nn, il, isd, irgm)
Definition: phs.f:504
convtc
subroutine convtc(num, nchar, loc)
Definition: ocn.f:156


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Level 3 Time Series Plotter

Request a SeaHawk Image

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster