OB.DAAC Logo
NASA Logo
Ocean Color Science Software

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