OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
afrt_rt1.f
Go to the documentation of this file.
1 c
2 c program main1
3 c
4 c***********************************************************************
5 c
6 c function- this program is the first part of the radiative transfer
7 c program.
8 c IN THIS VERSION OF THE CODE, THE INPUT READ STATEMENTS
9 C HAVE BEEN MODIFIED FOR OBPG PRODUCTION RUNS
10 c***********************************************************************
11 
12  SUBROUTINE af_rt1_read(infile,phsdir,outdir,
13  1 ozone_lut,atm_lut,coeff_lut,nsd,iair,iprin,ipsudo,ioptn,idust,
14  2 iaprof,iblyr,icld,nch,nt55,nrh,kset,ilam1,ilam2,isd1,isd2,
15  3 itau1,itau2,kznum,krhum,deltau,psrfc,ht_dust,sigma_dust,
16  4 pcldtop,pcldbtm,taum55,xlamb,c0,c1,c2,beta,rho,ppo,x101,
17  5 temp101,htdd,ppodd,dxdd,znpdd,nortau,ioznab,nrmww,
18  6 oqst,oqtt,oqst5,oqtt5)
19 
20  include 'afrt_rt1.cmn'
21 c
22 c parameter (klam=100,ktau=10,klyr=200,kair=3,kmac=100,
23 c 1 kozn=100,koznp=101,kdlyr=160,kdlvl=161,kmacp=101,
24 c 2 nsc=1000,krh=10,ksd=100)
25  integer*4 iair(3),krhum(ksd),iprin,ipsudo,ioptn,
26  1 idust,iaprof,iblyr,icld,nch,nt55,nrh,kset,ilam1,ilam2,
27  2 itau1,itau2,kznum,nortau,ioznab,nrmww,nsd
28  real*8 wvlth2,deltau,psrfc,ht_dust,sigma_dust,pcldtop,pcldbtm,nmnum
29  real*8 taum55(ktau),beta(klam),xlamb(klam),c0(klam),c1(klam),c2(klam),
30  1 rho(klam),ppo(koznp),x101(koznp),temp101(koznp),
31  2 htdd(kmac),ppodd(kmac),dxdd(kmac),znpdd(kmac),watrdd(kmac)
32  real*8 oqst(ktau,ksd,klam),oqtt(ktau,ksd,klam),
33  1 oqst5(ktau,ksd,klam),oqtt5(ktau,ksd,klam)
34  character*255 infile,ozone_lut,atm_lut,coeff_lut,phsdir,outdir,tmat,
35  1 tmat55,txt
36  character*90 tmata,tmat55a
37  character*4 cndex
38  character*2 dndex,endex,wwcc,krhc,ksetc
39  character*255 root
40  data root /'$RTDATA'/
41 c
42 c**********************************************************************
43 c
44  n=koznp
45 c nm=kmac
46  nm=85
47  nd=kdlvl
48 c
49 c dir1='../rt1/'
50 c dir2='../rt1/output/'
51 c dir3='../phs/output/'
52 c
53 c open(5,file='rt1.dir',status='old')
54 c read(5,'(a)')infile
55 c read(5,'(a)')outdir
56 c read(5,'(a)')phsdir
57 c
58  len1=index(infile,' ')-1
59  len2=index(outdir,' ')-1
60  len3=index(phsdir,' ')-1
61  len4=index(ozone_lut,' ')-1
62  len5=index(atm_lut,' ')-1
63  len6=index(coeff_lut,' ')-1
64 c
65  call getenv('RTDATA',root)
66  lenr = lenstr(root)
67 c
68 c open(unit=6,file='/tmp/dump',status='unknown')
69  open(unit=4,file=infile(1:len1),status='old')
70  open(unit=7,file=ozone_lut(1:len4),status='old')
71  open(unit=8,file=atm_lut(1:len5),status='old')
72  open(unit=15,file=coeff_lut(1:len6),status='old')
73 c
74 c read flags and other processing parameters
75 c
76  read(4,1002)txt
77 c read(4,1000)(iair(i),i=1,3),iprin,ipsudo,ioptn,idust,
78 c 1 iaprof,iblyr,icld
79 c**********************************************************
80  read(4,1000)(iair(i),i=1,3)
81  iprin=0
82  ipsudo=0
83  ioptn=0
84  idust=0
85  iaprof=1
86  iblyr=0
87  icld=0
88 c**********************************************************
89 c if idust=2 then ht_dust > top of the uniform dust layer and
90 c sigma_dust > bottom of the uniform dust layer
91 c if icld=1 then the program will replace aerosols with cloud
92 c if icld=0 then check idust,iaprof,iblyr flags (default iaprof=1)
93  if(idust.eq.0 .and. iaprof.eq.0 .and. iblyr.eq.0)then
94  iaprof=1
95  idust=0
96  iblyr=0
97  endif
98  if(idust.ge.1 .and. iaprof.eq.1 .and. iblyr.eq.1)then
99  iaprof=1
100  idust=0
101  iblyr=0
102  endif
103  if(iaprof.eq.0)then
104  if(idust.ge.1 .and. iblyr.eq.1)then
105  iaprof=0
106  idust=0
107  iblyr=1
108  endif
109  endif
110 c
111  read(4,1002)txt
112 c read(4,1003)nch,nt55,nortau,ioznab,nrmww
113 c**********************************************************
114  read(4,1003)nch,nt55,nrh,kset
115  nortau=0
116  ioznab=1
117  nrmww=0
118  if(kset.eq.0)nrh=1
119 c**********************************************************
120  call convtc(nrmww,2,wwcc)
121 c read starting and ending parameters for wavelength (1-15),
122 c size dist(1-11) and optical thickness (1-6)
123 c
124  read(4,1065)ilam1,ilam2
125  read(4,1065)isd1,isd2
126  read(4,1065)itau1,itau2
127 c
128  read(4,1002)txt
129 c read(4,1015)deltau,psrfc,ht_dust,sigma_dust
130 c**********************************************************
131  read(4,1015)deltau,psrfc
132  ht_dust=3.0
133  sigma_dust=1.0
134 c**********************************************************
135  read(4,1002)txt
136  read(4,1005)(taum55(i),i=1,nt55)
137 c
138  if(nrh.gt.0)then
139  read(4,1002)txt
140  read(4,1018)(krhum(i),i=1,nrh)
141  endif
142  if(nrh.eq.0)then
143  nrh=1
144  krhum(1)=0
145  endif
146 c
147  if(icld.eq.1)then
148  read(4,1002)txt
149  read(4,1015)pcldtop,pcldbtm !pressure in units of atm.
150  endif
151  close(4)
152 
153 c write(6,1005)(taum55(i),i=1,nt55)
154 c
155 c read wavelenghts, rayleigh scattering coeff., and absorption
156 c coefficients for ozone, h2o, co2 and molecular depolarization
157 c factors.
158 c ioznab=0 input is ozone optical thickness
159 c ioznab=1 input is ozone absorption coeff at temp t
160 c ioznab=2 input is ozone absorption at t=0c and the temp.
161 c dependent coeff. the input is read from a dataset:
162 c coeff.dat.
163 c
164  read(15,1002)txt
165  do kz=1,nch
166 c read(15,*)kznum,xlamb(kz),c0(kz),c1(kz),c2(kz),
167 c 1 beta(kz),rho(kz),ttozn(kz)
168  read(15,*)kznum,xlamb(kz),c0(kz),c1(kz),c2(kz),
169  1 beta(kz),rho(kz)
170 c write(6,111)kznum,xlamb(kz),c0(kz),c1(kz),c2(kz),
171 c 1 beta(kz),rho(kz)
172 111 format(i5,1p6e11.3)
173 
174  enddo
175  close(15)
176 c
177  do i=1,koznp
178  read(7,*)ppo(i),x101(i),temp101(i)
179  enddo
180  close(7)
181 
182  ndvsrf=k
183 c
184 c read height, pressure, ozone, and aerosol number density from
185 c a datatset which is used currently in the vpd program
186 c the dataset name is atmos_prof_vpd
187 c
188  do i=1,nm
189 c read (8,1220) htdd(i),ppodd(i),dxdd(i),znpdd(i)
190  read (8,*) htdd(i),ppodd(i),dxdd(i),znpdd(i)
191 c read (8,*) nmnum,htdd(i),ppodd(i),dxdd(i),znpdd(i)
192  watrdd(i)=0.0
193  enddo
194  close(8)
195 c
196 c
197 c create output datasets for modis aerosol look-up tables.
198 c
199  do 380 il=ilam1,ilam2,1
200 c if(xlamb(il).lt.1.0)then
201  lambz=xlamb(il)*1.0d3+0.01
202  wvlth=lambz*1.0e-3
203 c else
204 c lambz=xlamb(il)*1.0e4+0.1
205 c endif
206  do 390 isd=isd1,isd2,1
207  ih = (isd-1)*nrh/nsd + 1
208  do 400 itau=itau1,itau2,1
209 c
210 cbaf call convtc(lambz,4,cndex)
211  call convtc(idnint(xlamb(il)*1.0d3+0.01),4,cndex)
212  call convtc(isd,2,dndex)
213  call convtc(itau,2,endex)
214  call convtc(krhum(isd),2,krhc)
215  call convtc(kset,2,ksetc)
216 
217  if( kset.eq.0)krhc='00' !for rayleigh atmosphere set rh=0
218 
219  tmata='phs_twl'//cndex//'sd'//dndex//'rh'//krhc
220  k04=index(tmata,' ')-1
221  tmat=tmata(1:k04)//'_set'//ksetc//'.dat'
222 c
223  tmat55a='phs_twl'//wwcc//'sd'//dndex//'rh'//krhc
224  k05=index(tmat55a,' ')-1
225  tmat55=tmat55a(1:k05)//'_set'//ksetc//'.dat'
226 c
227  ntm=index(tmat,' ')-1
228  n55=index(tmat55,' ')-1
229 
230  oqst(itau,isd,il) = 1.0e-15
231  oqtt(itau,isd,il) = 1.0e-15
232 c
233 90 continue
234 110 continue
235  rewind(4)
236  rewind(7)
237  rewind(8)
238  rewind(15)
239 400 continue
240 390 continue
241 380 continue
242 c
243  return
244 c
245 c..................format statements....................................
246 c
247  1000 format(10i5)
248  1002 format(a)
249  1003 format(10i5)
250  1005 format(1p10e11.4)
251  1015 format(1p4e15.5)
252  1018 format(10i5)
253  1065 format(2i5)
254  1090 format(a)
255  1095 format(i5,11x,i5)
256  1105 format(1x,'input wavelength and phase matrix wavelength'
257  1 ' do not match.'/1x,'input wavelength=',1pe12.4,
258  2 'phase matrix wavelength=',1pe12.4)
259  1100 format(1p4e11.3,1p2e12.4)
260  end
261 
262 c
263 c**********************************************************************
264 c
265  SUBROUTINE af_rt1_process(outdir,ailam,aisd,airh,aitau,
266  1 iair,iprin,ipsudo,ioptn,idust,iaprof,iblyr,icld,nch,nt55,nsd,
267  2 nrh,kset,kznum,krhum,adeltau,psrfc,ht_dust,sigma_dust,
268  3 pcldtop,pcldbtm,taum55,xlamb,c0,c1,c2,beta,rho,ppo,x101,
269  4 temp101,htdd,ppodd,dxdd,znpdd,nortau,ioznab,nrmww,aqst,
270  5 aqtt,aqst5,aqtt5,ohtlvl,opplvl,ofr,ofn,ofa,oftot,odtrr,
271  6 odtmm,odtaa,odtot,otrp,otmp,otap,osalb,oturbhl,odx,ox,
272  7 otozn101,ohtdv,opdv,otaur,otaum,otaua,otau,owvlth,opsrfc,
273  8 orho,oozamtp,otautot,odtau1,onolyr,otcar,otwat,
274  9 otoznp,oifc,onmodl )
275 
276  implicit real*8 (a-h,o-z)
277  include 'afrt_rt1.cmn'
278 c
279  integer*4 ailam,aitau,aisd,airh,iair(3),krhum(ksd),iprin,
280  1 ipsudo,ioptn,idust,iaprof,iblyr,icld,nch,nt55,nsd,nrh,kset,
281  2 kznum,nortau,ioznab,nrmww,nmsz,asd
282  real*8 adeltau,psrfc,ht_dust,sigma_dust,pcldtop,pcldbtm
283  real*8 taum55(ktau),beta(klam),xlamb(klam),c0(klam),c1(klam),c2(klam),
284  1 rho(klam),ppo(koznp),x101(koznp),temp101(koznp),
285  2 htdd(kmac),ppodd(kmac),dxdd(kmac),znpdd(kmac),watrdd(kmac)
286  real*8 aqst,aqtt,aqst5,aqtt5,ohtlvl(klyr),opplvl(klyr),ofr(klyr),
287  1 ofn(klyr),ofa(klyr),oftot(klyr)
288  real*8 odtrr(klyr),odtmm(klyr),odtaa(klyr),odtot(klyr),
289  1 osalb(klyr),oturbhl(klyr),otrp,otmp,otap,
290  2 odx(koznp),ox(koznp),otozn101(koznp)
291  real*8 ohtdv(kdlvl),opdv(kdlvl),otaur(klyr),otaum(klyr),
292  1 otaua(klyr),otau(klyr),owvlth,opsrfc
293  real*8 orho,oozamtp,otautot,odtau1,otcar,otwat,otoznp
294  integer*4 onolyr,oifc,onmodl
295 
296  real*8 ttozn(klam),th2o(klam),tco2(klam),ppod(kmac),htd(kmac),znp(kmac),
297  1 dxd(kmac),watr(kmac),dx(kozn),x(koznp),pdv(kdlvl),htdv(kdlvl),
298  2 dtrr(klyr),dtmm(klyr),abozn(klam), dtaa(klyr),dtot(klyr),
299  3 temp1(klyr),temp2(klyr),xxd(kmacp),tozn101(koznp),xd(kmac)
300  real*8 wvlth,tr,tmt,ta,twat,tcar,rrho,oznabs,qst,qtt,
301  1 xozn,tautot,dtau1,tozn,p2km,h2km,cc0,cc1,cc2
302  real*8 fr(klyr),fn(klyr),fa(klyr),ftot(klyr),taur(klyr),
303  1 taum(klyr),taua(klyr),tauw(klyr),tauc(klyr),
304  2 tauozn(klyr),tau(klyr),lgppod(kmac),lgpdv(kdlvl),
305  3 htlvl(klyr),pplvl(klyr)
306  real*8 ht(kdlvl),hho(kdlvl),zpp(kdlvl),zppt(kdlvl),dnm(kdlvl)
307  integer*4 n,nm,ifc,nolyr,kndx,ndvsrf
308 
309  character*255 outdir,out2,out3,out23,txt
310  character*90 out2a,out3a,out23a
311  character*4 cndex
312  character*2 dndex,endex,wwcc,krhc,ksetc
313  character*255 root
314  data root /'$RTDATA'/
315 c
316 c**********************************************************************
317 c
318  deltau = adeltau
319  n=koznp
320 c nm=kmac
321  nd=kdlvl
322  nm = 85
323 
324  ioptn = nt55
325  do i=1,nch
326  th2o(i)=0.0d0
327  tco2(i)=0.0d0
328  enddo
329 c
330  do i=1,nch
331  th2o(i)=0.0d0
332  tco2(i)=0.0d0
333  enddo
334 
335 c express dx in atm-cm
336  do i=1,kozn
337  dx(i)=(x101(i+1)-x101(i))*1.0d-3
338  enddo
339 c
340 c compute pressure corresponding to 161 level of dave's vpd program
341 c use the emperical formula given in vpd report
342 c
343  a1=-3.4219d-3
344  a2=6.7056d-2
345  do i=1,nd
346  pdv(i)=1.0d0*dexp(a1*(nd-i)*(1.0d0+a2*(nd-i)))
347  enddo
348 c
349 c if(psrfc .lt. 1.0d0) then move the nearest dave pressure level
350 c to psrfc and redefine nd
351 c
352  do i=1,kdlvl
353  k=i
354  if(pdv(i).ge.psrfc)then
355  goto 116
356  endif
357  enddo
358 116 continue
359 
360  ndvsrf=k
361 
362  if(icld.eq.1)then
363  call cloud_prof(pcldtop,pcldbtm,ppodd,htdd,znpdd)
364  else
365  if(idust.ge.1)then
366  call dust_prof(znpdd,htdd)
367  endif
368 c
369  if(iblyr.eq.1)then
370 c compute the pressure coressponding to 2km
371  h2km=2.0d0
372  do i=1,nm
373  temp1(i)=dlog(ppodd(i))
374  enddo
375  call linear(h2km,htdd,temp1,nm,xlp2km)
376  p2km=dexp(xlp2km)/ppodd(1)
377 c write(6,181)h2km,xlp2km,p2km
378 181 format('h2km,xlp2km,p2km',1p3e11.3)
379  diff=dabs(pdv(1)-p2km)
380  do i=1,nd
381  diffx=dabs(pdv(i)-p2km)
382  if(diffx.le.diff)then
383  kndx=i
384  diff=diffx
385  endif
386  enddo
387  pdv(kndx)=p2km
388  endif
389  endif
390 c
391 c store top of the atmosphere in buffer element 1 and bottom in nm
392 c
393  do 30 i=1,nm
394  k = nm - i + 1
395  ppod(i) = ppodd(k) / ppodd(1)
396  htd(i)=htdd(k)
397  znp(i)=znpdd(k)
398  dxd(i)=dxdd(k)
399  watr(i)=watrdd(k)
400 30 continue
401 c
402 c create output datasets for modis aerosol look-up tables.
403 c
404  il = ailam
405  isd = aisd
406  ih = airh
407  itau = aitau
408 c
409  nmsz = 10
410  asd = mod(isd-1,nmsz) + 1
411  lambz=xlamb(il)*1.0d3+0.01
412  call convtc(lambz,4,cndex)
413  call convtc(idnint(xlamb(il)*1.0d3+0.01),4,cndex)
414  call convtc(asd,2,dndex)
415  call convtc(itau,2,endex)
416  call convtc(krhum(isd),2,krhc)
417  call convtc(kset,2,ksetc)
418 c
419 c write(*,*) ksetc,' ',endex,' ',dndex,' ',krhc,' ',cndex
420 c
421  if( kset.eq.0)krhc='00' !for rayleigh atmosphere set rh=0
422 
423  out2a='rt1_wl'//cndex//'sd'//dndex//'rh'//krhc
424  k01=index(out2a,' ')-1
425  out2=out2a(1:k01)//'ta'//endex//'_set'//ksetc//'.out'
426 c
427  out3a='rt1_spwl'//cndex//'sd'//dndex//'rh'//krhc
428  k02=index(out3a,' ')-1
429  out3=out3a(1:k02)//'ta'//endex//'_set'//ksetc//'.dat'
430 c
431  out23a='rt1_wl'//cndex//'sd'//dndex//'rh'//krhc
432  k03=index(out23a,' ')-1
433  out23=out23a(1:k03)//'ta'//endex//'_set'//ksetc//'.dat'
434 c
435  n02=index(out2,' ')-1
436  n03=index(out3,' ')-1
437  n23=index(out23,' ')-1
438 c
439  len2=index(outdir,' ')-1
440  open(2,file=outdir(1:len2)//'/'//out2(1:n02),status='unknown')
441  if(ipsudo.eq.1)then
442  open(3,file=outdir(1:len2)//'/'//out3(1:n03),status='unknown')
443  endif
444 c
445  open(unit=23,file=outdir(1:len2)//'/'//out23(1:n23),
446  . status='unknown')
447 c
448 c convert wavelength in cm
449 c
450  wvlth=xlamb(il)/10000.0d0
451  tr=beta(il)
452  tozn=ttozn(il)
453  cc0=c0(il)
454  cc1=c1(il)
455  cc2=c2(il)
456  twat=th2o(il)
457  tcar=tco2(il)
458  rrho=rho(il)
459 c if(ipol.eq.0)rrho=0.0d0
460  ifc = 0
461  if (iair(2) .eq. 1) ifc=1
462  if(taum55(itau).lt. 1.0d-6)ifc=0
463 c
464 c....set scattering and extinction coefficients equal to 1.0e-15
465 c
466  qst5 = aqst5
467  qtt5 = aqtt5
468  qst = aqst
469  qtt = aqtt
470 
471 c compute the aerosol optical thickness for the
472 c wavelength corresponding to taum55(itau)
473 c
474  if (iair(2) .ne. 1) go to 90
475 c wvlth2=wvlth2*1.0e-4
476  wvlth2=wvlth
477 
478  dwv=dabs(wvlth2-wvlth)
479  if(dwv .ge.1.0e-8)then
480  write(2,1105)wvlth,wvlth2
481  write(*,1105)wvlth,wvlth2
482  stop
483  endif
484  if(nortau.eq.1)then
485  tmt=(qtt/qtt5)*taum55(itau)
486  write (2,1165) qtt5,qst5
487  write (2,1170) qtt,qst
488  write(2,1175)taum55(itau),tmt
489  write(2,1177)wvlth*1.0d4
490  else
491  tmt=taum55(itau)
492  write(2,1170)qtt,qst
493  write(2,1176)tmt
494  write(2,1177)wvlth*1.0d4
495  endif
496 c
497  90 continue
498 c
499  write (2,1040)
500 c
501 c write(*,*)'ready to call tauin'
502 c write(*,*)'nd,psrfc',nd,psrfc
503 c write(*,*)' Main...ioznab,cc0,cc1,cc2',ioznab,cc0,cc1,cc2
504  call tauint( itau,isd,ih,il,ifc,iair,iprin,ipsudo,ioptn,idust,iaprof,
505  1 iblyr,icld,nd,nm,ndvsrf,psrfc,deltau,qst,qtt,tmt,znp,zpp,
506  2 pdv,fr,fn,fa,ftot,tau,taum,taua,taur,tauw,tauc,tauozn,ppod,lgppod,lgpdv,
507  3 x,dx,xd,dxd,htlvl,pplvl,htd,dtrr,dtmm,dtaa,dtot,htdv,tr,tcar,wvlth,rrho,
508  4 ohtlvl,opplvl,ofr,ofn,ofa,oftot,odtrr,odtmm,odtaa,odtot,osalb,
509  5 oturbhl,otrp,otmp,otap,odx,ox,temp101,x101,ppo,otozn101,ohtdv,
510  6 opdv,otaur,otaum,otaua,otau,owvlth,opsrfc,orho,oozamtp,otautot,odtau1,
511  7 onolyr,otcar,otwat,otoznp,oifc,onmodl,aqst,aqtt,aqst5,aqtt5 )
512 c
513 110 continue
514  close(2)
515  close(3)
516  close(23)
517 c
518  return
519 c
520 c..................format statements....................................
521 c
522  1020 format(2f10.3)
523  1022 format(t5,'wvlth',t17,'psrfc',t29,'rho',t41,'xozn',t53,'tautot',
524  1 t65,'deltau',t74,'nolyr')
525  1023 format(1p6e12.4,i4)
526  1024 format(t5,'tr',t17,'tmt',t29,'ta',t41,'tcar',t53,'twat',
527  1 t65,'tozn',t75,'ifc')
528  1025 format(1p6e12.4,i4)
529  1026 format(t3,'no',t9,'dtrr',t21,'dtmm',t33,'dtaa',t45,'dtot')
530  1027 format(i4,1p4e12.4)
531  1030 format(5e15.5)
532  1035 format(3f10.3)
533  1040 format (/t5,'radiative transfer calculations based on',
534  1 ' the clear atmosphere model')
535  1050 format (/t5,'radiative transfer calculations based on',
536  1 ' the thin cloud model')
537  1060 format(4e15.5,2i5)
538  1070 format(4e15.5)
539  1090 format(a)
540  1095 format(i5,11x,i5)
541  1105 format(1x,'input wavelength and phase matrix wavelength'
542  1 ' do not match.'/1x,'input wavelength=',1pe12.4,
543  2 'phase matrix wavelength=',1pe12.4)
544  1100 format(1p4e11.3,1p2e12.4)
545  1120 format(1p2d18.12,1p2d19.12,0pf6.2)
546  1130 format(i4,25f4.1)
547  1140 format(i4,10f4.2)
548  1165 format('volume extinction coefficient (550 nm)',t40,'=',1pe15.5,
549  1 '(cm-1)'/'volume scattering coefficient (550 nm)', t40,'=',
550  2 1pe15.5, '(cm-1)' / )
551  1170 format('volume extinction coefficient',t40,'=',1pe15.5,
552  1 '(cm-1)' / 'volume scattering coefficient', t40,'=',
553  2 1pe15.5, '(cm-1)' / )
554  1175 format('aerosol opt. thickness (550 nm)=',1pe15.5/
555  1 'aerosol opt. thickness =',1pe15.5)
556  1176 format('aerosol opt. thickness =',1pe15.5)
557  1177 format('wavelength (um) =',1pe15.5)
558  1180 format(4d18.10,5x,i3)
559  1220 format(f5.0,4d12.4)
560  1230 format(5x,f5.1,3f10.3)
561  end
562 
563 c***********************************************************************
564  subroutine tauint(itau,isd,ih,il,ifc,iair,iprin,ipsudo,ioptn,idust,iaprof,
565  1 iblyr,icld,nd,nm,ndvsrf,psrfc,deltau,qst,qtt,tmt,znp,zpp,
566  2 pdv,fr,fn,fa,ftot,tau,taum,taua,taur,tauw,tauc,tauozn,ppod,
567  3 lgppod,lgpdv,x,dx,xd,dxd,htlvl,pplvl,htd,dtrr,dtmm,dtaa,
568  4 dtot,htdv,tr,tcar,wvlth,rrho, ohtlvl,opplvl,ofr,ofn,ofa,
569  5 oftot,odtrr,odtmm,odtaa,odtot,osalb,oturbhl,otrp,otmp,otap,odx,ox,
570  6 temp101,x101,ppo,otozn101,ohtdv,opdv,otaur,otaum,otaua,otau,owvlth,
571  7 opsrfc,orho,oozamtp,otautot,odtau1,onolyr,otcar,otwat,
572  8 otoznp,oifc,onmodl,aqst,aqtt,aqst5,aqtt5 )
573 c
574 c***********************************************************************
575 c
576 c name- tauint
577 c
578 c***********************************************************************
579 c
580  include 'afrt_rt1.cmn'
581 c
582  integer*4 itau,isd,ih,il,nd,ndvsrf,n,nm,ifc,nolyr,kndx,
583  1 iair(3),krhum(ksd),iprin,ipsudo,ioptn,idust,iaprof,iblyr,icld
584  real*8 temp101(koznp),x101(koznp),ppo(koznp),aqst,aqtt,aqst5,tmt,
585  1 aqtt5,ohtlvl(klyr),opplvl(klyr),ofr(klyr),ofn(klyr),ofa(klyr),
586  2 oftot(klyr),psrfc,deltau,qst,qtt,tr,wvlth,rrho
587  real*8 odtrr(klyr),odtmm(klyr),odtaa(klyr),odtot(klyr),znp(kmac),
588  1 zpp(kdlvl),osalb(klyr),oturbhl(klyr),otrp,otmp,otap,
589  2 odx(koznp),ox(koznp),otozn101(koznp)
590  real*8 ohtdv(kdlvl),opdv(kdlvl),otaur(klyr),otaum(klyr),
591  1 otaua(klyr),otau(klyr),owvlth,opsrfc
592  real*8 orho,oozamtp,otautot,odtau1,otcar,otwat,otoznp,pp,pplg,argi
593  integer*4 onolyr,oifc,onmodl
594  real*8 pdv(kdlvl),fr(klyr),fn(klyr),fa(klyr),ftot(klyr),tau(klyr),
595  1 taum(klyr),taua(klyr),taur(klyr),tauw(klyr),tauc(klyr),
596  2 tauozn(klyr),ppod(kmac),lgppod(kmac),lgpdv(kdlvl),x(koznp),
597  3 dx(kozn),xd(kmac),dxd(kmac),htlvl(klyr),pplvl(klyr),htd(kmac),
598  4 dtrr(klyr),dtmm(klyr),dtaa(klyr),dtot(klyr),htdv(kdlvl),tcar
599 
600  nd=ndvsrf
601  pdv(nd)=psrfc
602 
603 c write(*,*)'tauin nd,pdv(nd)',nd,pdv(nd)
604 c....set the elements to zero
605  do i=1,klyr
606  fa(i) = 0.0d0
607  fn(i) = 0.0d0
608  fr(i) = 0.0d0
609  taum(i) = 0.0d0
610  taur(i) = 0.0d0
611  taua(i) = 0.0d0
612  tauw(i) = 0.0d0
613  tauc(i) = 0.0d0
614  tauozn(i)=0.0d0
615  enddo
616 c
617  do i=1,nm
618  lgppod(i)=dlog(ppod(i))
619  enddo
620 c
621  const = qst/qtt
622  do i=1,nd
623  lgpdv(i)=dlog(pdv(i))
624  enddo
625 c
626  if (iair(3).eq.1)then
627  x(1) = 0.0d0
628  if(ioznab.eq.0)then
629 c
630 c determine total ozone in the profile(note that the
631 c ozone density is gm/m**3 and the height is in km.
632 c also, the input is not given at equal delta height)
633 c first determine the scale height and from it the ozone
634 c amount above level 1
635 c
636  denm=dlog(dxd(1)/dxd(2))
637  o3h=-(htd(1)-htd(2))/denm
638  xd(1)=o3h*(1.0d+3*dxd(1)/2.1429d-2)
639  nm1=kmac-1
640  do i=1,nm1
641  xd(i+1)=xd(i)+0.5d0*1.0d+3*(htd(i)-htd(i+1))*
642  1 (dxd(i)+dxd(i+1))/2.1429d-2
643  enddo
644  ozamt=xd(kmac)
645 c write(6,*)'ozone amount(in vpd profile)=',ozamt
646  oznabs=tozn/ozamt
647 c determine ozone optical thickness profile as a
648 c function of pdv
649  call ozn3
650  toznp=tozn
651  ozamtp=ozamt
652  if(toznp.gt.0d0)then
653  write (2,1070) ozamtp
654  endif
655  endif
656  if(ioznab.ge.1)then
657  nm1 = n - 1
658  x(1)=x101(1)
659  do i=1,nm1
660  x(i+1) = x(i) + dx(i)
661  enddo
662  ozamt = x(n)
663 c write(*,*)'ozone amount in the input profile',ozamt
664 c determine ozone optical thickness profile as a
665 c function of pdv
666  call ozn3_abs(temp101,ppo)
667  toznp=tauozn(nd)
668  ozamtp=x(n)
669  if(toznp.gt.0d0)then
670  write (2,1070) ozamtp
671  endif
672  endif
673 c
674  endif
675 c
676 c determine water vapor optical thickness profile
677 c
678 c write(6,*)'ready to call sub water'
679 c call water(twat,pdv,watr,ppod,htd,nm,nd,tauw)
680 c
681 c determine co2 optical thickness profile
682 c
683 c write(6,*)'ready to call co2'
684  call co2(tcar,pdv,tauc,nd)
685 c
686  do i=1,nd
687  taua(i)=tauozn(i)+tauw(i)+tauc(i)
688  enddo
689 c write(6,198)(tauc(i),i=1,nd)
690 198 format('tauc',1p6e11.3)
691 c write(6,199)(taua(i),i=1,nd)
692 199 format('taua',1p6e11.3)
693 c
694 c write(6,*)'ready to call part2'
695 c write(*,*)'nd,pdv(nd)',nd,pdv(nd)
696 c write(6,*)iair
697 c write(6,*)'iaprof',iaprof
698 c
699  if (iair(2).eq.1)then
700  if(idust.ge.1 .or. iaprof.eq.1)then
701  call part1(nm,pdv,ppod,htd,znp,zpp,qtt,tmt,taum)
702  endif
703  if(iblyr.eq.1)then
704  call part2(nm,pdv,ppod,htd,znp,zpp,qtt,tmt,taum)
705  endif
706  endif
707 c
708  if (iair(1).eq.1)then
709  do i=1,nd
710  taur(i) = tr * pdv(i)
711  enddo
712  endif
713 c
714  do i=1,nd
715  tau(i)=taum(i)+taur(i)+taua(i)
716  enddo
717  do i=1,nd
718  call linear (lgpdv(i),lgppod,htd,nm,htdv(i))
719 c write(6,115)i,htdv(i),pdv(i),taur(i),taum(i),taua(i),tau(i)
720  enddo
721 c
722  if(ipsudo.eq.1)then
723  write(3,121)
724  write(3,122)nd
725  write(3,128)
726  do i=1,nd
727  write(3,115)i,htdv(i),pdv(i),taur(i),taum(i),taua(i),tau(i)
728  ohtdv(i) = htdv(i)
729  opdv(i) = pdv(i)
730  otaur(i) = taur(i)
731  otaum(i) = taum(i)
732  otaua(i) = taua(i)
733  otau(i) = tau(i)
734  enddo
735  onmodl = nd
736  endif
737 c
738  if(tau(nd).lt.0.02d0)then
739  deltau=tau(nd)/2.0d0
740  endif
741 c
742  if(ioptn.le.0 .or. ioptn.gt.4)ioptn=0
743 c
744  if(ioptn.eq.0)then
745  ntau=tau(nd)/deltau
746  ntau=ntau+1
747  dtau1=tau(nd)/dfloat(ntau)
748  ntaup=ntau+1
749 c
750 c write(6,*)'ntau,dtau1,deltau',ntau,dtau1,deltau
751 c fr(1)=tr*ppod(1)
752 c fn(1)=taum(1)
753 c fa(1)=taua(1)
754 c
755 c assign a value of 1.0e-10 fr(1),fn(1) and fa(1) (top of the atm.)
756 c
757  if(taur(nd).gt.0.0d0)fr(1)=1.0d-10
758  if(taum(nd).gt.0.0d0)fn(1)=1.0d-10
759  if(taua(nd).gt.0.0d0)fa(1)=1.0d-10
760  if(ifc.eq.0)fn(1)=0.0d0
761 c
762  htlvl(1)=htd(1)
763  pplvl(1)=ppod(1)
764  ftot(1)=fr(1)+fn(1)+fa(1)
765  do i=1,ntau
766  fi = i
767  argi=fi*dtau1
768  call linear (argi,tau,pdv,nd,pp)
769  pplvl(i+1) = pp
770  pplg=dlog(pp)
771  call linear (pplg,lgpdv,htdv,nd,htlvl(i+1))
772  if(htlvl(i+1).lt.0.0d0)htlvl(i+1)=0.0d0
773  if (iair(1).eq.1) call linear (pp,pdv,taur,nd,fr(i+1))
774  if (iair(2).eq.1) call linear (pp,pdv,taum,nd,fn(i+1))
775  if (iair(3).eq.1) call linear (pp,pdv,taua,nd,fa(i+1))
776  ftot(i+1)=fr(i+1)+fn(i+1)+fa(i+1)
777  write (2,1110)i,htlvl(i),pplvl(i),fr(i),fn(i),fa(i),ftot(i)
778 
779  ohtlvl(i) = htlvl(i)
780  opplvl(i) = pplvl(i)
781  ofr(i) = fr(i)
782  ofn(i) = fn(i)
783  ofa(i) = fa(i)
784  oftot(i) = ftot(i)
785 
786  enddo
787  write (2,1110)ntaup,htlvl(ntaup),pplvl(ntaup),fr(ntaup),
788  1 fn(ntaup),fa(ntaup),ftot(ntaup)
789 
790  ohtlvl(ntaup) = htlvl(ntaup)
791  opplvl(ntaup) = pplvl(ntaup)
792  ofr(ntaup) = fr(ntaup)
793  ofn(ntaup) = fn(ntaup)
794  ofa(ntaup) = fa(ntaup)
795  oftot(ntaup) = ftot(ntaup)
796 
797  endif
798 c
799  if(ioptn.gt.0)then
800  jdl=1
801  if(ioptn.eq.1)jdl=1
802  if(ioptn.eq.2)jdl=2
803  if(ioptn.eq.3)jdl=4
804  if(ioptn.eq.4)jdl=8
805 c
806  k=0
807  a1=-3.4219d-3
808  a2=6.7056d-2
809  do i=1,kdlvl,jdl
810  k=k+1
811  pplvl(k)=1.0d0*dexp(a1*(kdlvl-i)*(1.0d0+a2*(kdlvl-i)))
812  if(pplvl(k).gt.psrfc)then
813  goto 120
814  endif
815  enddo
816 120 continue
817  pplvl(k)=psrfc
818  ntau=k-1
819  il=1
820 c fr(1)=tr*ppod(1)
821 c fn(1)=taum(1)
822 c fa(1)=taua(1)
823 c
824 c assign a value of 1.0e-10 fr(1),fn(1) and fa(1) (top of the atm.)
825 c
826  if(taur(nd).gt.0.0d0)fr(1)=1.0d-10
827  if(taum(nd).gt.0.0d0)fn(1)=1.0d-10
828  if(taua(nd).gt.0.0d0)fa(1)=1.0d-10
829  if(ifc.eq.0)fn(1)=0.0d0
830 c
831  htlvl(1)=htd(1)
832  ftot(1)=fr(1)+fn(1)+fa(1)
833  write (2,1110)il,htlvl(1),pplvl(1),fr(1),fn(1),fa(1),ftot(1)
834 
835  ohtlvl(1) = htlvl(1)
836  opplvl(1) = pplvl(1)
837  ofr(1) = fr(1)
838  ofn(1) = fn(1)
839  ofa(1) = fa(1)
840  oftot(1) = ftot(1)
841 
842  do i=2,k
843  pp=pplvl(i)
844  pplg=dlog(pp)
845  call linear (pplg,lgpdv,htdv,nd,htlvl(i))
846  if (iair(1).eq.1) call linear (pp,pdv,taur,nd,fr(i))
847  if (iair(2).eq.1) call linear (pp,pdv,taum,nd,fn(i))
848  if (iair(3).eq.1) call linear (pp,pdv,taua,nd,fa(i))
849  ftot(i)=fr(i)+fn(i)+fa(i)
850  write (2,1110)i,htlvl(i),pplvl(i),fr(i),fn(i),fa(i),ftot(i)
851 
852  ohtlvl(i) = htlvl(i)
853  opplvl(i) = pplvl(i)
854  ofr(i) = fr(i)
855  ofn(i) = fn(i)
856  ofa(i) = fa(i)
857  oftot(i) = ftot(i)
858 
859  enddo
860  endif
861 c
862 c
863  do i=1,ntau
864  dtrr(i) = fr(i+1) - fr(i)
865  dtmm(i) = fn(i+1) - fn(i)
866  dtaa(i) = fa(i+1) - fa(i)
867  if(dtrr(i).lt.0.0d0)dtrr(i)=0.0d0
868  if(dtmm(i).lt.0.0d0)dtmm(i)=0.0d0
869  if(dtaa(i).lt.0.0d0)dtaa(i)=0.0d0
870  enddo
871  write(2,1130)
872 c
873  write (2,1190)
874  write (2,1200)
875 c
876  trp=0.0
877  tmp=0.0
878  tap=0.0
879  nolyr=ntau
880  do i=1,ntau
881  trp=trp+dtrr(i)
882  tmp=tmp+dtmm(i)
883  tap=tap+dtaa(i)
884  tots = dtrr(i) + dtmm(i)*const
885  dtot(i) = dtrr(i) + dtmm(i) + dtaa(i)
886  salb = tots/dtot(i)
887  turbhl=(dtmm(i)*const)/tots
888  write (2,1210) i,pplvl(i+1),dtrr(i),dtmm(i),dtaa(i),dtot(i),
889  1 salb,turbhl
890 
891  opplvl(i) = pplvl(i+1)
892  odtrr(i) = dtrr(i)
893  odtmm(i) = dtmm(i)
894  odtaa(i) = dtaa(i)
895  odtot(i) = dtot(i)
896  osalb(i) = salb
897  oturbhl(i) = turbhl
898 
899  enddo
900 c
901  tautot = trp + tmp + tap
902 c
903  write(2,1230)trp,tmp,tap,tautot
904 
905  otrp = trp
906  otmp = tmp
907  otap = tap
908  otautot = tautot
909 
910  write(23,1022)
911  write(23,1023)wvlth,psrfc,rrho,ozamtp,tautot,dtau1,nolyr
912 
913  owvlth = wvlth
914  opsrfc = psrfc
915  orho = rrho
916  oozamtp = ozamtp
917  otautot = tautot
918  odtau1 = dtau1
919  onolyr = nolyr
920 
921  write(23,1024)
922  write(23,1025)trp,tmp,tap,tcar,twat,toznp,ifc
923 
924  otrp = trp
925  otmp = tmp
926  otap = tap
927  otcar = tcar
928  otwat = twat
929  otoznp = toznp
930  oifc = ifc
931 
932  write(23,1026)
933  do i=1,nolyr
934  write(23,1027)i,htlvl(i+1),pplvl(i+1),dtrr(i),dtmm(i),
935  1 dtaa(i),dtot(i)
936 
937  ohtlvl(i) = htlvl(i+1)
938  opplvl(i) = pplvl(i+1)
939  odtrr(i) = dtrr(i)
940  odtmm(i) = dtmm(i)
941  odtaa(i) = dtaa(i)
942  odtot(i) = dtot(i)
943 
944  enddo
945 c
946  return
947 c
948 c................format statements......................................
949 c
950  115 format(i4,1x,1p6e12.4)
951  121 format(t1,'tot_lvl')
952  122 format(i5)
953  128 format(t3,'lvl',t9,'ht(km)',t21,'p(atm)',t33,'taur',t45,'taum',
954  1 t57,'taua',t69,'tau_tot')
955  1000 format (15x,'total optical thickness at 3.6 km = ',1pe15.4/
956  1 15x,'new dtau=',1pe15.4/ 15x, 'new no. of layers=',i5//)
957  1010 format(/,10x,'number of layers above the cloud = ',i5,10x,
958  1 'layer thickness = ',f10.6)
959  1011 format(/,10x,'number of layers below the cloud = ',i5,10x,
960  1 'layer thickness = ',f10.6)
961  1012 format(/,10x,'number of layers in the cloud = ',i5,10x,
962  1 'layer thickness = ',f10.6)
963  1015 format (1x,t50,'clear atmosphere')
964  1016 format (1x,t50,'thin cloud')
965  1017 format (1x,t50,'thick cloud')
966  1020 format (/,t50,'cumulative optical thicknesses'///t11,
967  1 'press.',t27,'rayleigh', t47, 'mie',t67,'ozone',t87,'total'/
968  2 t11,'(in atm)'/)
969  1022 format(t5,'wvlth',t17,'psrfc',t29,'rho',t41,'xozn',t53,'tautot',
970  1 t65,'deltau',t74,'nolyr')
971  1023 format(1p6e12.4,i4)
972  1024 format(t5,'tr',t17,'tmt',t29,'ta',t41,'tcar',t53,'twat',
973  1 t65,'tozn',t75,'ifc')
974  1025 format(1p6e12.4,i4)
975  1026 format(t3,'no',t9,'htl_b',t21,'pl_b',t33,'dtrr',t45,'dtmm',
976  1 t57,'dtaa',t69,'dtot')
977  1027 format(i4,1p6e12.4)
978  1030 format (5x,'delta tau =',1pe15.5,5x,'no. of layers =',i3/)
979  1040 format ('interpolated values of press and',
980  1 ' tau(rayleigh) etc. corresponding to'/'delta tau=',1pe15.5//
981  2 t10,'height',t23,'press.',t31,'tau(rayleigh)',t46, 'tau(mie)',t56,
982  3 'tau(ozone)',t68,'tau(tot)'/t10,'(in km)',t22,'(in atm)'/)
983  1050 format ('values of tau(rayleigh) etc. in full and half layers'/
984  1'(the first and the last layer are half layers)'/)
985  1060 format (//t15,'division of atmosphere into layers of equal optical
986  1 thickness for pupose of calculating radiative transfer:'//)
987  1070 format (15x,'total ozone(top to psrfc) in atm-cm',1pe15.5/)
988  1080 format (i4,1x,1p6e12.4)
989  1090 format (i3,2x,1p4e11.3,0pf8.4)
990  1100 format (t2,'no',t8,'dtau(ray)',t19,'dtau(mie)',t30,
991  1 ' dtau(o3)',t41,'dtau(tot)',t53,'turb'/)
992  1110 format (i4,1x,1p6e12.4)
993  1120 format (3e15.5)
994  1130 format (' ')
995  1140 format (t10,'sum of the total dtau = ',f10.5,
996  1 /t10,'sum of the total dtaur = ', f10.5,
997  2 /t10,'sum of the total dtaum(scat.)= ', f10.5,
998  3 /t10,'sum of the total dtaua = ', f10.5)
999  1190 format(/t5,'values of tau(rayleigh ) etc. in each layer'/)
1000  1200 format(t2,'no',t6,'press(atm)'t17,'dtau(ray)',t28,
1001  1 'dtau(mie)',t39, 'dtau(o3)',t49,' dtau(tot)',
1002  2 t60,'alb_sscat',t71,'turb'/)
1003  1210 format(i3,1p5e11.3,0p2f8.4)
1004  1230 format(/t10,'tau_r =',e12.4/t10,'tau_m =',e12.4/
1005  1 t10,'tau_a =',e12.4/t10,'tau_t =',e12.4)
1006  end
1007 c
1008 c......subroutine part2..............................................
1009 c
1010  subroutine part2(nm,pdv,ppod,htd,znp,zpp,qtt,tmt,taum)
1011 c
1012  include 'afrt_rt1.cmn'
1013  real*8 sumxx(nsc+1),hxx(nsc+1),e(nsc+1),u(nsc+1)
1014  real*8 ttozn(klam),th2o(klam),tco2(klam),ppod(kmac),htd(kmac),znp(kmac),
1015  1 dxd(kmac),watr(kmac),dx(kozn),x(koznp),pdv(kdlvl),htdv(kdlvl),
1016  2 dtrr(klyr),dtmm(klyr),abozn(klam), dtaa(klyr),dtot(klyr),
1017  3 temp1(klyr),xxd(kmacp),tozn101(koznp),xd(kmac),temp2(klyr)
1018  real*8 wvlth,tr,tmt,ta,twat,tcar,rrho,oznabs,qst,qtt,
1019  1 xozn,tautot,dtau1,tozn,p2km,h2km,cc0,cc1,cc2
1020  real*8 fr(klyr),fn(klyr),fa(klyr),ftot(klyr),taur(klyr),
1021  1 taum(klyr),taua(klyr),tauw(klyr),tauc(klyr),
1022  2 tauozn(klyr),tau(klyr),lgppod(kmac),lgpdv(kdlvl),
1023  3 htlvl(klyr),pplvl(klyr)
1024  real*8 ht(kdlvl),hho(kdlvl),zpp(kdlvl),zppt(kdlvl),dnm(kdlvl)
1025  integer*4 n,nm,ifc,nolyr,kndx,ndvsrf,ndv
1026 c
1027 c determine the geo. height and num. density corresponding to pdv
1028 c
1029 c write(*,*)'welcome to sub part2'
1030 c write(*,*)'pdv(1)',pdv(1)
1031  ndv=kdlvl
1032  do i=1,ndv
1033  do j=1,nm-1
1034  if(pdv(i).ge.ppod(j) .and. pdv(i).lt.ppod(j+1))then
1035  xn1=dlog(pdv(i))-dlog(ppod(j))
1036  xd1=dlog(ppod(j+1))-dlog(ppod(j))
1037  hho(i)=htd(j)+(xn1/xd1)*(htd(j+1)-htd(j))
1038  prod1=(hho(i)-htd(j))/(htd(j+1)-htd(j))
1039  zpp(i)=znp(j)*(znp(j+1)/znp(j))**prod1
1040  endif
1041  enddo
1042  hho(ndv)=0.0d0
1043  zpp(ndv)=znp(nm)
1044 c write(6,226)i,pdv(i),hho(i),zpp(i)
1045 226 format('i,pdv,hho,zpp',i3,1x,1p3e11.3)
1046  enddo
1047 c
1048 c write(*,*)'pdv(1)',pdv(1)
1049 c integrate the number density
1050 c
1051  zppt(1)=0.0d0
1052  do i=1,ndv-1
1053  dnm(i)=0.5d0*(zpp(i)+zpp(i+1))*(hho(i)-hho(i+1))*1.0d+5
1054  zppt(i+1)=zppt(i)+dnm(i)
1055 c write(6,151)i+1,zppt(i+1)
1056 151 format('integrated # of part i+1,zppt(i+1)',i3,1pe11.3)
1057  enddo
1058 c
1059 c....convert the total number of particles in to optical depth
1060 c
1061 c tp=zppt(ndv)*qtt
1062 c normalize the total aerosol optical thickness at psrfc to the
1063 c input values
1064  tp=zppt(ndvsrf)*qtt
1065 c
1066  if(tmt.le.10.0d0)then
1067 c if(tmt.le.0.2d0)then
1068 c scale the optical depth 'tp' to the input value 'tmt'
1069 c
1070  c = tmt/tp
1071 c
1072  do i=1,ndv
1073  taum(i) = c * zppt(i) * qtt
1074  enddo
1075  return
1076  else
1077 c
1078 c scale the optical depth 'tp' to a value of 0.2 and determine
1079 c total aerosol optical thickness up to 2km
1080 c
1081  c=0.20d0/tp
1082  do i=1,ndv
1083  zpp(i)=zpp(i)*c
1084  zppt(i)=zppt(i)*c
1085  taum(i)=zppt(i)*qtt
1086 c write(6,119)i,zppt(i),taum(i)
1087 119 format('i,zppt(i),taum(i)',i3,1x,1p2e11.3)
1088  enddo
1089 c
1090 c
1091 c total aerosol opt and numbers from 2km to the ground
1092  t2b=tmt-taum(kndx)
1093  xn2b=t2b/qtt
1094  yn2b=zppt(ndv)-zppt(kndx)
1095 c
1096 c determine the scale height between 0 and 2km in the input profile
1097 c
1098  h0=2.0d0/(dlog(znp(nm)/znp(nm-2)))
1099 c
1100 c determine the low and high value of scale height that bracket
1101 c the total number of particles 'xn2b' from 0-2km
1102 c
1103  hx=h0*1.1d0
1104 100 continue
1105  sum2x=0.0d0
1106  do j=kndx+1,ndv
1107  zpp(j)=zpp(kndx)*dexp((hho(kndx)-hho(j))/hx)
1108  dnm(j-1)=0.5d0*(zpp(j-1)+zpp(j))*(hho(j-1)-hho(j))*1.0d+5
1109  sum2x=sum2x+dnm(j-1)
1110  enddo
1111  if(sum2x.gt.xn2b)then
1112  hg=hx
1113  goto 200
1114  endif
1115  hl=hx
1116  hx=hx*0.90d0
1117  if(hx.lt.0.10d0)then
1118 c write(6,600)hx
1119 600 format('the value of the aerosol scale height is less',
1120  1 ' than 0.10',1pe11.3)
1121  stop
1122  endif
1123  goto 100
1124 200 continue
1125 c write(6,*)'h0,hl,hx',h0,hl,hx
1126 c compute sumxx as a function of scale height hxx
1127 c
1128  nscp=nsc+1
1129  dhx=(hg-hl)/dfloat(nsc)
1130 c
1131  do i=1,nscp
1132  hx=hl+(i-1)*dhx
1133  sum2x=0.0d0
1134  do j=kndx+1,ndv
1135  zpp(j)=zpp(kndx)*dexp((hho(kndx)-hho(j))/hx)
1136  dnm(j-1)=0.5d0*(zpp(j-1)+zpp(j))*(hho(j-1)-hho(j))*
1137  1 1.0d+5
1138  sum2x=sum2x+dnm(j-1)
1139  enddo
1140  sumxx(i)=sum2x
1141  hxx(i)=hx
1142  enddo
1143 c
1144 c interpolate for scale height
1145 c
1146  in=1
1147  il=1
1148  iu=1
1149  call spline(xn2b,sumxx,hxx,nscp,in,hxz,il,iu,vl,vu,e,u)
1150 c
1151  do j=kndx+1,ndv
1152  zpp(j)=zpp(kndx)*dexp((hho(kndx)-hho(j))/hxz)
1153  dnm(j-1)=0.5d0*(zpp(j-1)+zpp(j))*(hho(j-1)-hho(j))*
1154  1 1.0d+5
1155  zppt(j)=zppt(j-1)+dnm(j-1)
1156  taum(j)=zppt(j)*qtt
1157  enddo
1158 c
1159  endif
1160 c
1161 c
1162 c
1163  return
1164  end
1165 c...............subroutine water...............................
1166 c
1167  subroutine water(twat,pdv,watr,ppod,ht,nm,ndv,tauw)
1168 c
1169  implicit real*8 (a-h,o-z)
1170  real*8 pdv(1000),watr(1000),ppod(1000),tauw(1000),u(1000)
1171  real*8 dts(1000),dt(1000),tauww(1000),e(1000),ht(1000)
1172 c
1173 c convert the water vapor density from gm/m**3 to gm/cm**2/km
1174 c (see mcclachy 72 p89)
1175  do i=1,nm
1176  watr(i)=watr(i)*1.0d-1
1177  enddo
1178 c
1179 c determine the total water amount
1180  nm1 = nm - 1
1181  tp = 0.0d0
1182  do 50 i=1,nm1
1183  tp = tp + (watr(i)+watr(i+1))/2.0d0
1184  50 continue
1185 c
1186 c scale the water vapor amount 'tp' to the input optical thickness
1187 c value 'twat'
1188 c
1189  c = twat/tp
1190 c
1191  do 60 i=1,nm
1192  dts(i) = c * watr(i)
1193  60 continue
1194  do 70 i=1,nm1
1195  dt(i) = (dts(i+1) + dts(i)) / 2.0d0
1196  70 continue
1197  tauww(1) = 0.
1198  do 80 i=1,nm
1199 c.....cumulative water vapor optical depth
1200  tauww(i+1) = tauww(i) + dt(i)
1201  80 continue
1202  in = 1
1203  il = 1
1204  iu = 1
1205 c
1206  do 85 i=1,ndv
1207  call spline (pdv(i),ppod,tauww,nm,in,tauw(i),il,iu,vl,vu,e,u)
1208  ttww=dabs(tauw(i))
1209  if(tauw(i).lt.0.0d0 .and. ttww.lt.1.0d-8)then
1210  tauw(i)=0.0d0
1211  endif
1212  85 continue
1213  return
1214  end
1215 c***********************************************************************
1216  subroutine ozn3
1217 c
1218  include 'afrt_rt1.cmn'
1219 c
1220  real*8 ee(klyr),uu(klyr)
1221  real*8 ppod(kmac),pdv(kdlvl),tauozn(klyr),oznabs
1222  integer*4 n,ndvsrf
1223 c**********************************************************************
1224  in=1
1225  il=1
1226  iu=1
1227  do i=1,ndvsrf
1228  call spline(pdv(i),ppod,xd,n,in,tempo3,il,iu,vl,vu,ee,uu)
1229  tauozn(i)=tempo3*oznabs
1230  enddo
1231 c
1232  return
1233  end
1234 c***********************************************************************
1235 c
1236  subroutine co2(tcar,pdv,tauc,ndv)
1237 c
1238  implicit real*8 (a-h,o-z)
1239  include 'afrt_rt1.cmn'
1240 c
1241  real*8 pdv(kdlvl),tauc(kdlvl)
1242 c
1243 c normalize the surface pressure to tcar
1244 c
1245  cc=tcar/pdv(ndv)
1246 c
1247  do i=1,ndv
1248  tauc(i)=pdv(i)*cc
1249  enddo
1250 c
1251  return
1252  end
1253 c***********************************************************************
1254  subroutine linear(xp,x,y,n,yp)
1255 c***********************************************************************
1256 c subroutine performs linear interpolation. It assumes that
1257 c x is a monotonically increasing parameter
1258 c***********************************************************************
1259  implicit real*8 (a-h,o-z)
1260  real*8 x(n),y(n)
1261 c
1262  il=1
1263  ih=n
1264 100 continue
1265  m=(il+ih)/2
1266  if(xp.gt.x(m))then
1267  il=m
1268  else
1269  ih=m
1270  endif
1271  id=ih-il
1272  if(id.eq.1)goto 200
1273  goto 100
1274 200 continue
1275 c
1276  slope=(y(ih)-y(il))/(x(ih)-x(il))
1277  yp=y(il)+slope*(xp-x(il))
1278 c
1279  return
1280  end
1281 c**********************************************************************
1282 c......subroutine part1..............................................
1283 c
1284  subroutine part1(nm,pdv,ppod,htd,znp,zpp,qtt,tmt,taum)
1285 c
1286  include 'afrt_rt1.cmn'
1287  integer*4 nm
1288  real*8 oznabs,qtt,tmt
1289  real*8 ppod(kmac),htd(kmac),znp(kmac),taum(klyr),pdv(kdlvl)
1290  real*8 ht(kdlvl),hho(kdlvl),zpp(kdlvl),zppt(kdlvl),dnm(kdlvl)
1291 c
1292 c determine the geo. height and num. density corresponding to pdv
1293 c
1294 c write(6,*)'welcome to part1'
1295  ndv=kdlvl
1296  do i=1,ndv
1297  do j=1,nm-1
1298  if(pdv(i).ge.ppod(j) .and. pdv(i).lt.ppod(j+1))then
1299  xn1=dlog(pdv(i))-dlog(ppod(j))
1300  xd1=dlog(ppod(j+1))-dlog(ppod(j))
1301  hho(i)=htd(j)+(xn1/xd1)*(htd(j+1)-htd(j))
1302  prod1=(hho(i)-htd(j))/(htd(j+1)-htd(j))
1303  if(znp(j).gt.0.0)then
1304  zpp(i)=znp(j)*(znp(j+1)/znp(j))**prod1
1305  else
1306  zpp(i)=0.0d0
1307  endif
1308  endif
1309  enddo
1310  hho(ndv)=0.0d0
1311  zpp(ndv)=znp(nm)
1312 c write(6,226)i,pdv(i),hho(i),zpp(i)
1313 226 format('i,pdv,hho,zpp',i3,1x,1p3e11.3)
1314  enddo
1315 c
1316 c integrate the number density
1317 c
1318  zppt(1)=0.0d0
1319  do i=1,ndv-1
1320  dnm(i)=0.5d0*(zpp(i)+zpp(i+1))*(hho(i)-hho(i+1))*1.0d+5
1321  zppt(i+1)=zppt(i)+dnm(i)
1322  enddo
1323 c
1324 c....convert the total number of particles in to optical depth
1325 c
1326  tp=zppt(ndv)*qtt
1327 c
1328 c scale the optical depth 'tp' to the input value 'tmt'
1329 c
1330  c = tmt/tp
1331 c
1332  do i=1,ndv
1333  taum(i) = c * zppt(i) * qtt
1334  enddo
1335  return
1336  end
1337 c*********************************************************************
1338 c
1339  subroutine dust_prof(znpdd,htdd)
1340 c
1341 c*********************************************************************
1342 c
1343  include 'afrt_rt1.cmn'
1344 c
1345  real*8 znpdd(kmac),htdd(kmac)
1346 
1347 c write(*,*)'idust',idust
1348  if(idust.eq.2)then
1349  top_cld=ht_dust
1350  btm_cld=sigma_dust
1351  do i=1,nm
1352  if(htdd(i).le.top_cld .and. htdd(i).gt.btm_cld)then
1353  znpdd(i)=1.0d0
1354  else
1355  znpdd(i)=0.0d0
1356  endif
1357  enddo
1358  else
1359  do i=1,nm
1360  znpdd(i)=exp(-0.5d0*((htdd(i)-ht_dust)/sigma_dust)**2)
1361  enddo
1362  endif
1363 c
1364  return
1365  end
1366 c***************************************************************************
1367  subroutine ozn3_abs(temp101,ppo)
1368 c
1369  include 'afrt_rt1.cmn'
1370 c
1371  real*8 ee(klyr),uu(klyr),temp101(koznp),ppo(koznp),tozn101(koznp),
1372  1 tauozn(klyr),dx(kozn),x(kozn),pdv(kdlvl)
1373 c determine the ozone optical thickness using temperature
1374 c dependent coefficient. if ioznab=1 then set temp. at all
1375 c levels to 227.16k (-46.0 c)
1376 c
1377  if(ioznab.eq.1)then
1378  do i=1,n
1379  temp101(i)=227.16
1380  enddo
1381  endif
1382 c
1383  alfatemp=cc0+cc1*(temp101(1)-273.16)+
1384  1 cc2*(temp101(1)-273.16)**2.0
1385  tozn101(1)=alfatemp*dx(1)
1386 c write(*,*)'alfatemp,tozn101(1)',alfatemp,tozn101(1)
1387  do i=2,n
1388  alfatemp=cc0+cc1*(0.5d0*(temp101(i)+temp101(i+1))-273.16)+
1389  1 cc2*(0.5d0*(temp101(i)+temp101(i+1))-273.16)**2.0
1390  tozn101(i)=tozn101(i-1)+alfatemp*dx(i)
1391  enddo
1392 c
1393  in=1
1394  il=1
1395  iu=1
1396  do j=1,n
1397  write(2,100)j,ppo(j),dx(j),x(j),tozn101(j)
1398 100 format('i,ppo,dx,x,tozn101',i4,1x,1p4e11.3)
1399  enddo
1400  do i=1,ndvsrf
1401  call spline(pdv(i),ppo,tozn101,n,in,tauo3z,il,iu,vl,vu,ee,uu)
1402  tauozn(i)=tauo3z
1403  enddo
1404 c
1405  return
1406  end
1407 c***********************************************************************
1408 
1409 c**********************************************************************
1410  subroutine cloud_prof(pcldtop,pcldbtm,ppodd,htdd,znpdd)
1411 c
1412  implicit real*8 (a-h,o-z)
1413  include 'afrt_rt1.cmn'
1414 c
1415  real*8 ppodd(kmac),znpdd(kmac),htdd(kmac),pcldtop,pcldbtm,hcldtop,
1416  1 hcldbtm,temp1(klyr),temp2(klyr),qxldtoplg,qcldbtmlg
1417  integer*4 i,k,nm
1418 c**********************************************************************
1419 c compute the height corresponding to the top and
1420 c bottom of the cloud layer
1421 c reverse the order of pressure and height. in linear, x should
1422 c be monotonically increasing in value.
1423 
1424  do i=1,nm
1425  k=nm-i+1
1426  temp1(i)=dlog(ppodd(k))
1427  temp2(i)=htdd(k)
1428  enddo
1429  qcldtop=(pcldtop*1.0d3)
1430  qcldbtm=(pcldbtm*1.0d3)
1431  qcldtoplg=dlog(qcldtop)
1432  qcldbtmlg=dlog(qcldbtm)
1433  call linear(qcldtoplg,temp1,temp2,nm,hcldtop)
1434  call linear(qcldbtmlg,temp1,temp2,nm,hcldbtm)
1435 c
1436 c write(*,*)'hcldtop,hcldbtm',hcldtop,hcldbtm
1437 c
1438  diff=dabs(ppodd(1)-qcldtop)
1439 c write(*,*)'ppodd(1),qcldtop,diff',ppodd(1),qcldtop,diff
1440  do i=1,nm
1441  diffx=dabs(ppodd(i)-qcldtop)
1442  if(diffx.le.diff)then
1443  kndx=i
1444  diff=diffx
1445  endif
1446  enddo
1447  ppodd(kndx)=qcldtop
1448  diff=dabs(ppodd(1)-qcldbtm)
1449  do i=1,nm
1450  diffx=dabs(ppodd(i)-qcldbtm)
1451  if(diffx.le.diff)then
1452  kmdx=i
1453  diff=diffx
1454  endif
1455  enddo
1456  ppodd(kmdx)=qcldbtm
1457 c
1458 c reset the number density
1459 c
1460  do i=1,nm
1461  znpdd(i)=0.0d0
1462  if(i.ge.kmdx .and. i.le.kndx)znpdd(i)=1.0d0
1463 c write(*,*)'i,znpdd(i)',i,znpdd(i)
1464  enddo
1465 c write(*,*)'kndx,kmdx',kndx,kmdx
1466 c write(*,*)'qcldtop,qcldbtm',qcldtop,qcldbtm
1467 c
1468  return
1469  end
1470 c***********************************************************************
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
af_rt1_read
subroutine af_rt1_read(infile, phsdir, outdir, ozone_lut, atm_lut, coeff_lut, nsd, iair, iprin, ipsudo, ioptn, idus iaprof, iblyr, icld, nch, nt55, nrh, kset, ilam1, ilam2, isd1, isd2, itau1, itau2, kznum, krhum, deltau, psrfc, ht_dust, sigma_dust, pcldtop, pcldbtm, taum55, xlamb, c0, c1, c2, beta, rho, ppo, x101, temp101, htdd, ppodd, dxdd, znpdd, nortau, ioznab, nrmww, oqst, oqtt, oqst5, oqtt5)
Definition: afrt_rt1.f:19
tauint
subroutine tauint(itau, isd, ih, il, ifc, iair, iprin, ipsudo, ioptn, idust iblyr, icld, nd, nm, ndvsrf, psrfc, deltau, qst, qtt, tmt, znp, zpp, pdv, fr, fn, fa, ftot, tau, taum, taua, taur, tauw, tauc, tauozn, ppod, lgppod, lgpdv, x, dx, xd, dxd, htlvl, pplvl, htd, dtrr, dtmm, dtaa, dtot, htdv, tr, tcar, wvlth, rrho, ohtlvl, opplvl, ofr, ofn, ofa, oftot, odtrr, odtmm, odtaa, odtot, osalb, oturbhl, otrp, otmp, otap, o temp101, x101, ppo, otozn101, ohtdv, opdv, otaur, otaum, otaua, otau, opsrfc, orho, oozamtp, otautot, odtau1, onolyr, otcar, otwat, otoznp, oifc, onmodl, aqst, aqtt, aqst5, aqtt5)
Definition: afrt_rt1.f:573
spline
subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)
Definition: phs.f:1348
dust_prof
subroutine dust_prof(znpdd, htdd)
Definition: afrt_rt1.f:1340
lenstr
integer *4 function lenstr(string)
Definition: lenstr.f:2
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
part2
subroutine part2(nm, pdv, ppod, htd, znp, zpp, qtt, tmt, taum)
Definition: afrt_rt1.f:1011
part1
subroutine part1(nm, pdv, ppod, htd, znp, zpp, qtt, tmt, taum)
Definition: afrt_rt1.f:1285
ozn3_abs
subroutine ozn3_abs(temp101, ppo)
Definition: afrt_rt1.f:1368
diff
subroutine diff(x, conec, n, dconecno, dn, dconecmk, units, u, inno, i, outno, o, input, deriv)
Definition: ffnet.f:205
af_rt1_process
subroutine af_rt1_process(outdir, ailam, aisd, airh, aitau, iair, iprin, ipsudo, ioptn, idust, iaprof, iblyr, icld, nch, nt55, nsd nrh, kset, kznum, krhum, adeltau, psrfc, ht_dust, sigma_dust, pcldtop, pcldbtm, taum55, xlamb, c0, c1, c2, beta, rho, ppo, x101, temp101, htdd, ppodd, dxdd, znpdd, nortau, ioznab, nrmww, aqst, aqtt, aqst5, aqtt5, ohtlvl, opplvl, ofr, ofn, ofa, oftot, odtrr, odtmm, odtaa, odtot, otrp, otmp, otap, osalb, oturbhl, odx, ox, otozn101, ohtdv, opdv, otaur, otaum, otaua, otau, owvlth, opsrfc, orho, oozamtp, otautot, odtau1, onolyr, otcar, otwat, otoznp, oifc, onmodl)
Definition: afrt_rt1.f:275
ozn3
subroutine ozn3
Definition: afrt_rt1.f:1217
cloud_prof
subroutine cloud_prof(pcldtop, pcldbtm, ppodd, htdd, znpdd)
Definition: afrt_rt1.f:1411
water
subroutine water(twat, pdv, watr, ppod, ht, nm, ndv, tauw)
Definition: afrt_rt1.f:1168
linear
subroutine linear(xp, x, y, n, yp)
Definition: afrt_rt1.f:1255
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

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