OB.DAAC Logo
NASA Logo
Ocean Color Science Software

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