OB.DAAC Logo
NASA Logo
Ocean Color Science Software

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