afrt_rt1.f

Go to the documentation of this file.

c                                                                       
c      program   main1                                                 
c                                                                       
c***********************************************************************
c                                                                       
c   function-  this program is the first part of the radiative transfer 
c              program.                                                 
c              IN THIS VERSION OF THE CODE, THE INPUT READ STATEMENTS
C              HAVE BEEN MODIFIED FOR OBPG PRODUCTION RUNS
c***********************************************************************
 
      SUBROUTINE af_rt1_read(infile,phsdir,outdir,
     1      ozone_lut,atm_lut,coeff_lut,nsd,iair,iprin,ipsudo,ioptn,idust,
     2      iaprof,iblyr,icld,nch,nt55,nrh,kset,ilam1,ilam2,isd1,isd2,
     3      itau1,itau2,kznum,krhum,deltau,psrfc,ht_dust,sigma_dust,
     4      pcldtop,pcldbtm,taum55,xlamb,c0,c1,c2,beta,rho,ppo,x101,
     5      temp101,htdd,ppodd,dxdd,znpdd,nortau,ioznab,nrmww,
     6      oqst,oqtt,oqst5,oqtt5)
 
      include 'afrt_rt1.cmn'
c
c      parameter (klam=100,ktau=10,klyr=200,kair=3,kmac=100,
c     1           kozn=100,koznp=101,kdlyr=160,kdlvl=161,kmacp=101,
c     2           nsc=1000,krh=10,ksd=100)
      integer*4 iair(3),krhum(ksd),iprin,ipsudo,ioptn,
     1      idust,iaprof,iblyr,icld,nch,nt55,nrh,kset,ilam1,ilam2,
     2      itau1,itau2,kznum,nortau,ioznab,nrmww,nsd
      real*8 wvlth2,deltau,psrfc,ht_dust,sigma_dust,pcldtop,pcldbtm,nmnum
      real*8 taum55(ktau),beta(klam),xlamb(klam),c0(klam),c1(klam),c2(klam),
     1       rho(klam),ppo(koznp),x101(koznp),temp101(koznp),
     2       htdd(kmac),ppodd(kmac),dxdd(kmac),znpdd(kmac),watrdd(kmac)
      real*8 oqst(ktau,ksd,klam),oqtt(ktau,ksd,klam),
     1       oqst5(ktau,ksd,klam),oqtt5(ktau,ksd,klam)
      character*255 infile,ozone_lut,atm_lut,coeff_lut,phsdir,outdir,tmat,
     1       tmat55,txt
      character*90 tmata,tmat55a
      character*4  cndex
      character*2  dndex,endex,wwcc,krhc,ksetc
      character*255 root
      data root /'$RTDATA'/
c
c**********************************************************************
c                                                                       
      n=koznp
c     nm=kmac
      nm=85
      nd=kdlvl
c
c     dir1='../rt1/'
c     dir2='../rt1/output/'
c     dir3='../phs/output/'
c
c      open(5,file='rt1.dir',status='old')
c      read(5,'(a)')infile
c      read(5,'(a)')outdir
c      read(5,'(a)')phsdir
c
      len1=index(infile,' ')-1
      len2=index(outdir,' ')-1
      len3=index(phsdir,' ')-1
      len4=index(ozone_lut,' ')-1
      len5=index(atm_lut,' ')-1
      len6=index(coeff_lut,' ')-1
c
      call getenv('RTDATA',root)
      lenr = lenstr(root)
c
c      open(unit=6,file='/tmp/dump',status='unknown')
      open(unit=4,file=infile(1:len1),status='old') 
      open(unit=7,file=ozone_lut(1:len4),status='old')
      open(unit=8,file=atm_lut(1:len5),status='old')
      open(unit=15,file=coeff_lut(1:len6),status='old')
c                                                                       
c     read flags and other processing parameters
c
      read(4,1002)txt
c     read(4,1000)(iair(i),i=1,3),iprin,ipsudo,ioptn,idust,
c    1            iaprof,iblyr,icld
c**********************************************************   
      read(4,1000)(iair(i),i=1,3)
      iprin=0
      ipsudo=0
      ioptn=0
      idust=0
      iaprof=1
      iblyr=0
      icld=0
c**********************************************************
c     if idust=2 then  ht_dust > top of the uniform dust layer and
c     sigma_dust > bottom of the uniform dust layer
c     if icld=1 then the program will replace aerosols with cloud
c     if icld=0 then check idust,iaprof,iblyr flags (default iaprof=1)
      if(idust.eq.0 .and. iaprof.eq.0 .and. iblyr.eq.0)then
         iaprof=1
         idust=0
         iblyr=0
      endif
      if(idust.ge.1 .and. iaprof.eq.1 .and. iblyr.eq.1)then
         iaprof=1
         idust=0
         iblyr=0
      endif
      if(iaprof.eq.0)then
         if(idust.ge.1 .and. iblyr.eq.1)then
            iaprof=0
            idust=0
            iblyr=1
         endif
      endif
c
      read(4,1002)txt
c     read(4,1003)nch,nt55,nortau,ioznab,nrmww
c**********************************************************
      read(4,1003)nch,nt55,nrh,kset
      nortau=0
      ioznab=1
      nrmww=0
      if(kset.eq.0)nrh=1
c**********************************************************
      call convtc(nrmww,2,wwcc)
c     read starting and ending parameters for wavelength (1-15),
c     size dist(1-11) and optical thickness (1-6)
c
      read(4,1065)ilam1,ilam2
      read(4,1065)isd1,isd2
      read(4,1065)itau1,itau2
c
      read(4,1002)txt
c     read(4,1015)deltau,psrfc,ht_dust,sigma_dust
c**********************************************************
      read(4,1015)deltau,psrfc
      ht_dust=3.0
      sigma_dust=1.0
c**********************************************************      
      read(4,1002)txt
      read(4,1005)(taum55(i),i=1,nt55)
c
      if(nrh.gt.0)then
         read(4,1002)txt
         read(4,1018)(krhum(i),i=1,nrh)
      endif
      if(nrh.eq.0)then
         nrh=1
         krhum(1)=0
      endif
c
      if(icld.eq.1)then
         read(4,1002)txt
         read(4,1015)pcldtop,pcldbtm   !pressure in units of atm.
      endif
      close(4)
 
c     write(6,1005)(taum55(i),i=1,nt55)
c                                                                       
c     read wavelenghts, rayleigh scattering coeff., and absorption
c     coefficients for ozone, h2o, co2 and molecular depolarization
c     factors.  
c     ioznab=0   input is ozone optical thickness
c     ioznab=1   input is ozone absorption coeff at temp t
c     ioznab=2   input is ozone absorption at t=0c and the temp.
c                dependent coeff.  the input is read from a dataset:
c                coeff.dat.
c
      read(15,1002)txt
      do kz=1,nch
c        read(15,*)kznum,xlamb(kz),c0(kz),c1(kz),c2(kz),
c    1             beta(kz),rho(kz),ttozn(kz)
         read(15,*)kznum,xlamb(kz),c0(kz),c1(kz),c2(kz),
     1             beta(kz),rho(kz)
c        write(6,111)kznum,xlamb(kz),c0(kz),c1(kz),c2(kz),
c    1             beta(kz),rho(kz)
111      format(i5,1p6e11.3)
 
      enddo
      close(15)
c
      do i=1,koznp
         read(7,*)ppo(i),x101(i),temp101(i)
      enddo
      close(7)
 
      ndvsrf=k
c
c     read height, pressure, ozone, and aerosol number density from
c     a datatset which is used currently in the vpd program
c     the dataset name is atmos_prof_vpd
c
      do i=1,nm
c        read (8,1220) htdd(i),ppodd(i),dxdd(i),znpdd(i)
         read (8,*) htdd(i),ppodd(i),dxdd(i),znpdd(i)
c         read (8,*) nmnum,htdd(i),ppodd(i),dxdd(i),znpdd(i)
         watrdd(i)=0.0
      enddo
      close(8)
c
c
c      create output datasets for modis aerosol look-up tables.
c
      do 380 il=ilam1,ilam2,1
c         if(xlamb(il).lt.1.0)then
            lambz=xlamb(il)*1.0d3+0.01
            wvlth=lambz*1.0e-3
c         else
c            lambz=xlamb(il)*1.0e4+0.1
c         endif
      do 390 isd=isd1,isd2,1
      ih = (isd-1)*nrh/nsd + 1
      do 400 itau=itau1,itau2,1
c
cbaf      call convtc(lambz,4,cndex)
      call convtc(idnint(xlamb(il)*1.0d3+0.01),4,cndex)
      call convtc(isd,2,dndex)
      call convtc(itau,2,endex)
      call convtc(krhum(isd),2,krhc)
      call convtc(kset,2,ksetc)
 
      if( kset.eq.0)krhc='00'  !for rayleigh atmosphere set rh=0
 
      tmata='phs_twl'//cndex//'sd'//dndex//'rh'//krhc
      k04=index(tmata,' ')-1
      tmat=tmata(1:k04)//'_set'//ksetc//'.dat'
c
      tmat55a='phs_twl'//wwcc//'sd'//dndex//'rh'//krhc
      k05=index(tmat55a,' ')-1
      tmat55=tmat55a(1:k05)//'_set'//ksetc//'.dat'
c
      ntm=index(tmat,' ')-1
      n55=index(tmat55,' ')-1
 
      oqst(itau,isd,il) = 1.0e-15
      oqtt(itau,isd,il) = 1.0e-15
c
90    continue
110   continue
      rewind(4)
      rewind(7)
      rewind(8)
      rewind(15)
400   continue
390   continue
380   continue
c
      return
c
c..................format statements....................................
c
 1000 format(10i5)
 1002 format(a)
 1003 format(10i5)
 1005 format(1p10e11.4)
 1015 format(1p4e15.5)
 1018 format(10i5)
 1065 format(2i5)
 1090 format(a)
 1095 format(i5,11x,i5)
 1105 format(1x,'input wavelength and phase matrix wavelength'
     1       ' do not match.'/1x,'input wavelength=',1pe12.4,
     2         'phase matrix wavelength=',1pe12.4)
 1100 format(1p4e11.3,1p2e12.4)
      end
 
c
c**********************************************************************
c
      SUBROUTINE af_rt1_process(outdir,ailam,aisd,airh,aitau,
     1      iair,iprin,ipsudo,ioptn,idust,iaprof,iblyr,icld,nch,nt55,nsd,
     2      nrh,kset,kznum,krhum,adeltau,psrfc,ht_dust,sigma_dust,
     3      pcldtop,pcldbtm,taum55,xlamb,c0,c1,c2,beta,rho,ppo,x101,
     4      temp101,htdd,ppodd,dxdd,znpdd,nortau,ioznab,nrmww,aqst,
     5      aqtt,aqst5,aqtt5,ohtlvl,opplvl,ofr,ofn,ofa,oftot,odtrr,
     6      odtmm,odtaa,odtot,otrp,otmp,otap,osalb,oturbhl,odx,ox,
     7      otozn101,ohtdv,opdv,otaur,otaum,otaua,otau,owvlth,opsrfc,
     8      orho,oozamtp,otautot,odtau1,onolyr,otcar,otwat,
     9      otoznp,oifc,onmodl )
 
      implicit real*8 (a-h,o-z)
      include 'afrt_rt1.cmn'
c
      integer*4 ailam,aitau,aisd,airh,iair(3),krhum(ksd),iprin,
     1      ipsudo,ioptn,idust,iaprof,iblyr,icld,nch,nt55,nsd,nrh,kset,
     2      kznum,nortau,ioznab,nrmww,nmsz,asd
      real*8 adeltau,psrfc,ht_dust,sigma_dust,pcldtop,pcldbtm
      real*8 taum55(ktau),beta(klam),xlamb(klam),c0(klam),c1(klam),c2(klam),
     1       rho(klam),ppo(koznp),x101(koznp),temp101(koznp),
     2       htdd(kmac),ppodd(kmac),dxdd(kmac),znpdd(kmac),watrdd(kmac)
      real*8 aqst,aqtt,aqst5,aqtt5,ohtlvl(klyr),opplvl(klyr),ofr(klyr),
     1       ofn(klyr),ofa(klyr),oftot(klyr)
      real*8 odtrr(klyr),odtmm(klyr),odtaa(klyr),odtot(klyr),
     1       osalb(klyr),oturbhl(klyr),otrp,otmp,otap,
     2       odx(koznp),ox(koznp),otozn101(koznp)
      real*8 ohtdv(kdlvl),opdv(kdlvl),otaur(klyr),otaum(klyr),
     1       otaua(klyr),otau(klyr),owvlth,opsrfc
      real*8 orho,oozamtp,otautot,odtau1,otcar,otwat,otoznp
      integer*4 onolyr,oifc,onmodl
 
      real*8 ttozn(klam),th2o(klam),tco2(klam),ppod(kmac),htd(kmac),znp(kmac),
     1       dxd(kmac),watr(kmac),dx(kozn),x(koznp),pdv(kdlvl),htdv(kdlvl),
     2       dtrr(klyr),dtmm(klyr),abozn(klam), dtaa(klyr),dtot(klyr),
     3       temp1(klyr),temp2(klyr),xxd(kmacp),tozn101(koznp),xd(kmac)
      real*8 wvlth,tr,tmt,ta,twat,tcar,rrho,oznabs,qst,qtt,
     1       xozn,tautot,dtau1,tozn,p2km,h2km,cc0,cc1,cc2
      real*8 fr(klyr),fn(klyr),fa(klyr),ftot(klyr),taur(klyr),
     1       taum(klyr),taua(klyr),tauw(klyr),tauc(klyr),
     2       tauozn(klyr),tau(klyr),lgppod(kmac),lgpdv(kdlvl),
     3       htlvl(klyr),pplvl(klyr)
      real*8 ht(kdlvl),hho(kdlvl),zpp(kdlvl),zppt(kdlvl),dnm(kdlvl)
      integer*4 n,nm,ifc,nolyr,kndx,ndvsrf
 
      character*255 outdir,out2,out3,out23,txt
      character*90 out2a,out3a,out23a
      character*4  cndex
      character*2  dndex,endex,wwcc,krhc,ksetc
      character*255 root
      data root /'$RTDATA'/
c
c**********************************************************************
c
      deltau = adeltau
      n=koznp
c     nm=kmac
      nd=kdlvl
      nm = 85
 
      ioptn = nt55
      do i=1,nch
         th2o(i)=0.0d0
         tco2(i)=0.0d0
      enddo
c
      do i=1,nch
         th2o(i)=0.0d0
         tco2(i)=0.0d0
      enddo   
 
c     express dx in atm-cm
      do i=1,kozn
         dx(i)=(x101(i+1)-x101(i))*1.0d-3
      enddo
c
c     compute pressure corresponding to 161 level of dave's vpd program
c     use the emperical formula given in vpd report
c
      a1=-3.4219d-3
      a2=6.7056d-2     
      do i=1,nd
         pdv(i)=1.0d0*dexp(a1*(nd-i)*(1.0d0+a2*(nd-i))) 
      enddo
c
c     if(psrfc .lt. 1.0d0) then move the nearest dave pressure level
c     to psrfc and redefine nd
c
      do i=1,kdlvl
         k=i
         if(pdv(i).ge.psrfc)then
            goto 116
         endif
      enddo
116   continue
 
      ndvsrf=k
 
      if(icld.eq.1)then
         call cloud_prof(pcldtop,pcldbtm,ppodd,htdd,znpdd)
      else
          if(idust.ge.1)then
            call dust_prof(znpdd,htdd)
          endif
c
          if(iblyr.eq.1)then
c            compute the pressure coressponding to 2km
             h2km=2.0d0
                do i=1,nm
                   temp1(i)=dlog(ppodd(i))
                enddo
             call linear(h2km,htdd,temp1,nm,xlp2km)
             p2km=dexp(xlp2km)/ppodd(1)
c            write(6,181)h2km,xlp2km,p2km
181          format('h2km,xlp2km,p2km',1p3e11.3)
             diff=dabs(pdv(1)-p2km)
             do i=1,nd
                diffx=dabs(pdv(i)-p2km)
                if(diffx.le.diff)then
                   kndx=i
                   diff=diffx
                endif
             enddo
             pdv(kndx)=p2km
          endif
      endif
c
c     store top of the atmosphere in buffer element 1 and bottom in nm
c
      do 30 i=1,nm                                                      
        k = nm - i + 1                                                  
        ppod(i) = ppodd(k) / ppodd(1)                                   
        htd(i)=htdd(k)
        znp(i)=znpdd(k)
        dxd(i)=dxdd(k)
        watr(i)=watrdd(k)
30    continue                                                          
c
c      create output datasets for modis aerosol look-up tables.
c
      il = ailam
      isd = aisd
      ih = airh
      itau = aitau
c
      nmsz = 10
      asd = mod(isd-1,nmsz) + 1
      lambz=xlamb(il)*1.0d3+0.01
      call convtc(lambz,4,cndex)
      call convtc(idnint(xlamb(il)*1.0d3+0.01),4,cndex) 
      call convtc(asd,2,dndex)
      call convtc(itau,2,endex)                                        
      call convtc(krhum(isd),2,krhc)
      call convtc(kset,2,ksetc)
c
c      write(*,*) ksetc,' ',endex,' ',dndex,' ',krhc,' ',cndex
c
      if( kset.eq.0)krhc='00'  !for rayleigh atmosphere set rh=0
 
      out2a='rt1_wl'//cndex//'sd'//dndex//'rh'//krhc
      k01=index(out2a,' ')-1
      out2=out2a(1:k01)//'ta'//endex//'_set'//ksetc//'.out'
c
      out3a='rt1_spwl'//cndex//'sd'//dndex//'rh'//krhc
      k02=index(out3a,' ')-1
      out3=out3a(1:k02)//'ta'//endex//'_set'//ksetc//'.dat'
c
      out23a='rt1_wl'//cndex//'sd'//dndex//'rh'//krhc
      k03=index(out23a,' ')-1
      out23=out23a(1:k03)//'ta'//endex//'_set'//ksetc//'.dat'
c
      n02=index(out2,' ')-1
      n03=index(out3,' ')-1
      n23=index(out23,' ')-1
c
      len2=index(outdir,' ')-1
      open(2,file=outdir(1:len2)//'/'//out2(1:n02),status='unknown')
      if(ipsudo.eq.1)then
         open(3,file=outdir(1:len2)//'/'//out3(1:n03),status='unknown')
      endif
c
      open(unit=23,file=outdir(1:len2)//'/'//out23(1:n23),
     . status='unknown')
c                                                                       
c     convert wavelength in cm
c
      wvlth=xlamb(il)/10000.0d0                                          
      tr=beta(il)                                                       
      tozn=ttozn(il)
      cc0=c0(il)
      cc1=c1(il)
      cc2=c2(il)
      twat=th2o(il)
      tcar=tco2(il)
      rrho=rho(il)                                                      
c     if(ipol.eq.0)rrho=0.0d0
      ifc = 0                                                           
      if (iair(2) .eq. 1) ifc=1                                         
      if(taum55(itau).lt. 1.0d-6)ifc=0
c                                                                       
c....set  scattering and extinction coefficients equal to 1.0e-15          
c
      qst5 = aqst5
      qtt5 = aqtt5
      qst = aqst
      qtt = aqtt
 
c     compute the aerosol optical thickness for the 
c     wavelength corresponding to taum55(itau) 
c
      if (iair(2) .ne. 1) go to 90
c          wvlth2=wvlth2*1.0e-4
          wvlth2=wvlth
 
          dwv=dabs(wvlth2-wvlth)
          if(dwv .ge.1.0e-8)then
             write(2,1105)wvlth,wvlth2
             write(*,1105)wvlth,wvlth2
             stop
          endif
          if(nortau.eq.1)then
             tmt=(qtt/qtt5)*taum55(itau)
             write (2,1165) qtt5,qst5
             write (2,1170) qtt,qst
             write(2,1175)taum55(itau),tmt
             write(2,1177)wvlth*1.0d4
          else
             tmt=taum55(itau)
             write(2,1170)qtt,qst
             write(2,1176)tmt
             write(2,1177)wvlth*1.0d4
          endif
c
   90 continue
c                                                                       
       write (2,1040)                                   
c                                                                       
c     write(*,*)'ready to call tauin'
c     write(*,*)'nd,psrfc',nd,psrfc
c     write(*,*)' Main...ioznab,cc0,cc1,cc2',ioznab,cc0,cc1,cc2
      call tauint( itau,isd,ih,il,ifc,iair,iprin,ipsudo,ioptn,idust,iaprof,
     1  iblyr,icld,nd,nm,ndvsrf,psrfc,deltau,qst,qtt,tmt,znp,zpp,
     2  pdv,fr,fn,fa,ftot,tau,taum,taua,taur,tauw,tauc,tauozn,ppod,lgppod,lgpdv,
     3  x,dx,xd,dxd,htlvl,pplvl,htd,dtrr,dtmm,dtaa,dtot,htdv,tr,tcar,wvlth,rrho,
     4  ohtlvl,opplvl,ofr,ofn,ofa,oftot,odtrr,odtmm,odtaa,odtot,osalb,
     5  oturbhl,otrp,otmp,otap,odx,ox,temp101,x101,ppo,otozn101,ohtdv,
     6  opdv,otaur,otaum,otaua,otau,owvlth,opsrfc,orho,oozamtp,otautot,odtau1,
     7  onolyr,otcar,otwat,otoznp,oifc,onmodl,aqst,aqtt,aqst5,aqtt5 )
c                                                                       
110   continue                                                          
      close(2)                                                          
      close(3)
      close(23)
c
      return
c                                                                       
c..................format statements....................................
c                                                                       
 1020 format(2f10.3)
 1022 format(t5,'wvlth',t17,'psrfc',t29,'rho',t41,'xozn',t53,'tautot',
     1       t65,'deltau',t74,'nolyr')
 1023 format(1p6e12.4,i4)
 1024 format(t5,'tr',t17,'tmt',t29,'ta',t41,'tcar',t53,'twat',
     1       t65,'tozn',t75,'ifc')
 1025 format(1p6e12.4,i4)
 1026 format(t3,'no',t9,'dtrr',t21,'dtmm',t33,'dtaa',t45,'dtot')
 1027 format(i4,1p4e12.4)
 1030 format(5e15.5)                                                    
 1035 format(3f10.3)                                              
 1040 format (/t5,'radiative transfer calculations based on',         
     1 ' the clear atmosphere model')                                   
 1050 format (/t5,'radiative transfer calculations based on',         
     1 ' the thin cloud model')                                         
 1060 format(4e15.5,2i5)                                                
 1070 format(4e15.5)                                                    
 1090 format(a)                                                      
 1095 format(i5,11x,i5)
 1105 format(1x,'input wavelength and phase matrix wavelength'
     1       ' do not match.'/1x,'input wavelength=',1pe12.4,
     2         'phase matrix wavelength=',1pe12.4)
 1100 format(1p4e11.3,1p2e12.4)
 1120 format(1p2d18.12,1p2d19.12,0pf6.2)                                
 1130 format(i4,25f4.1)                                                 
 1140 format(i4,10f4.2)                                                 
 1165 format('volume extinction coefficient (550 nm)',t40,'=',1pe15.5,     
     1  '(cm-1)'/'volume scattering coefficient (550 nm)', t40,'=',    
     2 1pe15.5, '(cm-1)'     /    )                                     
 1170 format('volume extinction coefficient',t40,'=',1pe15.5,     
     1  '(cm-1)'   /   'volume scattering coefficient', t40,'=',    
     2 1pe15.5, '(cm-1)'     /    )                                     
 1175 format('aerosol opt. thickness (550 nm)=',1pe15.5/
     1       'aerosol opt. thickness         =',1pe15.5)
 1176 format('aerosol opt. thickness         =',1pe15.5)
 1177 format('wavelength (um)                =',1pe15.5)
 1180 format(4d18.10,5x,i3)                                             
 1220 format(f5.0,4d12.4)
 1230 format(5x,f5.1,3f10.3)
      end
 
c***********************************************************************
      subroutine tauint(itau,isd,ih,il,ifc,iair,iprin,ipsudo,ioptn,idust,iaprof,
     1      iblyr,icld,nd,nm,ndvsrf,psrfc,deltau,qst,qtt,tmt,znp,zpp,
     2      pdv,fr,fn,fa,ftot,tau,taum,taua,taur,tauw,tauc,tauozn,ppod,
     3      lgppod,lgpdv,x,dx,xd,dxd,htlvl,pplvl,htd,dtrr,dtmm,dtaa,
     4      dtot,htdv,tr,tcar,wvlth,rrho, ohtlvl,opplvl,ofr,ofn,ofa,
     5      oftot,odtrr,odtmm,odtaa,odtot,osalb,oturbhl,otrp,otmp,otap,odx,ox,
     6      temp101,x101,ppo,otozn101,ohtdv,opdv,otaur,otaum,otaua,otau,owvlth,
     7      opsrfc,orho,oozamtp,otautot,odtau1,onolyr,otcar,otwat,
     8      otoznp,oifc,onmodl,aqst,aqtt,aqst5,aqtt5 )
c                                                                       
c***********************************************************************
c                                                                       
c   name- tauint                                                        
c
c***********************************************************************
c                                                                       
      include 'afrt_rt1.cmn'
c                                                                       
      integer*4 itau,isd,ih,il,nd,ndvsrf,n,nm,ifc,nolyr,kndx,
     1        iair(3),krhum(ksd),iprin,ipsudo,ioptn,idust,iaprof,iblyr,icld
      real*8 temp101(koznp),x101(koznp),ppo(koznp),aqst,aqtt,aqst5,tmt,
     1       aqtt5,ohtlvl(klyr),opplvl(klyr),ofr(klyr),ofn(klyr),ofa(klyr),
     2       oftot(klyr),psrfc,deltau,qst,qtt,tr,wvlth,rrho
      real*8 odtrr(klyr),odtmm(klyr),odtaa(klyr),odtot(klyr),znp(kmac),
     1       zpp(kdlvl),osalb(klyr),oturbhl(klyr),otrp,otmp,otap,
     2       odx(koznp),ox(koznp),otozn101(koznp)
      real*8 ohtdv(kdlvl),opdv(kdlvl),otaur(klyr),otaum(klyr),
     1       otaua(klyr),otau(klyr),owvlth,opsrfc
      real*8 orho,oozamtp,otautot,odtau1,otcar,otwat,otoznp,pp,pplg,argi
      integer*4 onolyr,oifc,onmodl
      real*8 pdv(kdlvl),fr(klyr),fn(klyr),fa(klyr),ftot(klyr),tau(klyr),
     1       taum(klyr),taua(klyr),taur(klyr),tauw(klyr),tauc(klyr),
     2       tauozn(klyr),ppod(kmac),lgppod(kmac),lgpdv(kdlvl),x(koznp),
     3       dx(kozn),xd(kmac),dxd(kmac),htlvl(klyr),pplvl(klyr),htd(kmac),
     4       dtrr(klyr),dtmm(klyr),dtaa(klyr),dtot(klyr),htdv(kdlvl),tcar
 
      nd=ndvsrf
      pdv(nd)=psrfc
 
c     write(*,*)'tauin nd,pdv(nd)',nd,pdv(nd)
c....set the elements to zero                                           
      do i=1,klyr
         fa(i) = 0.0d0                                                      
         fn(i) = 0.0d0                                                      
         fr(i) = 0.0d0                                                      
         taum(i) = 0.0d0                                                    
         taur(i) = 0.0d0                                                    
         taua(i) = 0.0d0                                                    
         tauw(i) = 0.0d0
         tauc(i) = 0.0d0
         tauozn(i)=0.0d0
      enddo
c                                                                       
      do i=1,nm
         lgppod(i)=dlog(ppod(i))
      enddo
c
      const = qst/qtt                                                     
      do i=1,nd
         lgpdv(i)=dlog(pdv(i))
      enddo
c
      if (iair(3).eq.1)then
         x(1) = 0.0d0                                                        
         if(ioznab.eq.0)then
c
c           determine total ozone in the profile(note that the 
c           ozone density is gm/m**3 and the height is in km.
c           also, the  input is not given at equal delta height)
c           first determine the scale height and from it the ozone
c           amount above level 1
c
            denm=dlog(dxd(1)/dxd(2))
            o3h=-(htd(1)-htd(2))/denm
            xd(1)=o3h*(1.0d+3*dxd(1)/2.1429d-2)
            nm1=kmac-1
            do i=1,nm1
               xd(i+1)=xd(i)+0.5d0*1.0d+3*(htd(i)-htd(i+1))*
     1              (dxd(i)+dxd(i+1))/2.1429d-2
            enddo
            ozamt=xd(kmac)
c           write(6,*)'ozone amount(in vpd profile)=',ozamt
            oznabs=tozn/ozamt
c           determine ozone optical thickness profile as a 
c           function of pdv
            call ozn3
            toznp=tozn
            ozamtp=ozamt
            if(toznp.gt.0d0)then
               write (2,1070) ozamtp  
            endif
         endif
         if(ioznab.ge.1)then
            nm1 = n - 1
            x(1)=x101(1)
            do i=1,nm1
               x(i+1) = x(i) + dx(i)
            enddo
            ozamt = x(n)
c           write(*,*)'ozone amount in the input profile',ozamt
c           determine ozone optical thickness profile as a
c           function of pdv
            call ozn3_abs(temp101,ppo)
            toznp=tauozn(nd)
            ozamtp=x(n)
            if(toznp.gt.0d0)then
               write (2,1070) ozamtp
            endif
         endif
c          
      endif
c
c     determine water vapor optical thickness profile
c
c     write(6,*)'ready to call sub water'
c     call water(twat,pdv,watr,ppod,htd,nm,nd,tauw)
c
c     determine co2 optical thickness profile
c
c     write(6,*)'ready to call co2'
      call co2(tcar,pdv,tauc,nd)
c
      do i=1,nd                                                       
         taua(i)=tauozn(i)+tauw(i)+tauc(i)
      enddo
c        write(6,198)(tauc(i),i=1,nd)
198      format('tauc',1p6e11.3)
c        write(6,199)(taua(i),i=1,nd)
199      format('taua',1p6e11.3)
c
c     write(6,*)'ready to call part2'
c     write(*,*)'nd,pdv(nd)',nd,pdv(nd)
c     write(6,*)iair
c     write(6,*)'iaprof',iaprof
c
      if (iair(2).eq.1)then
         if(idust.ge.1 .or. iaprof.eq.1)then
            call part1(nm,pdv,ppod,htd,znp,zpp,qtt,tmt,taum)
         endif
         if(iblyr.eq.1)then
            call part2(nm,pdv,ppod,htd,znp,zpp,qtt,tmt,taum)
         endif
      endif
c
      if (iair(1).eq.1)then 
         do i=1,nd                                                      
            taur(i) = tr * pdv(i)                                           
         enddo
      endif
c
      do i=1,nd
         tau(i)=taum(i)+taur(i)+taua(i)
      enddo
      do i=1,nd                                                     
          call linear (lgpdv(i),lgppod,htd,nm,htdv(i))
c         write(6,115)i,htdv(i),pdv(i),taur(i),taum(i),taua(i),tau(i)
      enddo
c
      if(ipsudo.eq.1)then
        write(3,121)
        write(3,122)nd
        write(3,128) 
        do i=1,nd                                                         
           write(3,115)i,htdv(i),pdv(i),taur(i),taum(i),taua(i),tau(i)
           ohtdv(i) = htdv(i)
           opdv(i) = pdv(i)
           otaur(i) = taur(i)
           otaum(i) = taum(i)
           otaua(i) = taua(i)
           otau(i) = tau(i)
        enddo
        onmodl = nd
      endif                                                             
c
      if(tau(nd).lt.0.02d0)then
        deltau=tau(nd)/2.0d0
      endif
c
      if(ioptn.le.0 .or. ioptn.gt.4)ioptn=0
c
      if(ioptn.eq.0)then
        ntau=tau(nd)/deltau
        ntau=ntau+1
        dtau1=tau(nd)/dfloat(ntau)
        ntaup=ntau+1
c                                                                       
c       write(6,*)'ntau,dtau1,deltau',ntau,dtau1,deltau
c       fr(1)=tr*ppod(1)
c       fn(1)=taum(1)
c       fa(1)=taua(1)
c
c       assign a value of 1.0e-10 fr(1),fn(1) and fa(1) (top of the atm.)
c
        if(taur(nd).gt.0.0d0)fr(1)=1.0d-10
        if(taum(nd).gt.0.0d0)fn(1)=1.0d-10
        if(taua(nd).gt.0.0d0)fa(1)=1.0d-10
        if(ifc.eq.0)fn(1)=0.0d0
c
        htlvl(1)=htd(1)
        pplvl(1)=ppod(1)
        ftot(1)=fr(1)+fn(1)+fa(1)
        do i=1,ntau                                                   
           fi = i                                                         
           argi=fi*dtau1
           call linear (argi,tau,pdv,nd,pp)                                  
           pplvl(i+1) = pp                                                      
           pplg=dlog(pp)
           call linear (pplg,lgpdv,htdv,nd,htlvl(i+1))
           if(htlvl(i+1).lt.0.0d0)htlvl(i+1)=0.0d0
           if (iair(1).eq.1) call linear (pp,pdv,taur,nd,fr(i+1))              
           if (iair(2).eq.1) call linear (pp,pdv,taum,nd,fn(i+1))              
           if (iair(3).eq.1) call linear (pp,pdv,taua,nd,fa(i+1))              
           ftot(i+1)=fr(i+1)+fn(i+1)+fa(i+1) 
           write (2,1110)i,htlvl(i),pplvl(i),fr(i),fn(i),fa(i),ftot(i)
 
           ohtlvl(i) = htlvl(i)
           opplvl(i) = pplvl(i)
           ofr(i) = fr(i)
           ofn(i) = fn(i)
           ofa(i) = fa(i)
           oftot(i) = ftot(i)
 
        enddo
        write (2,1110)ntaup,htlvl(ntaup),pplvl(ntaup),fr(ntaup),
     1                fn(ntaup),fa(ntaup),ftot(ntaup)
 
        ohtlvl(ntaup) = htlvl(ntaup)
        opplvl(ntaup) = pplvl(ntaup)
        ofr(ntaup) = fr(ntaup)
        ofn(ntaup) = fn(ntaup)
        ofa(ntaup) = fa(ntaup)
        oftot(ntaup) = ftot(ntaup)
 
      endif
c
      if(ioptn.gt.0)then
        jdl=1
        if(ioptn.eq.1)jdl=1
        if(ioptn.eq.2)jdl=2
        if(ioptn.eq.3)jdl=4
        if(ioptn.eq.4)jdl=8
c
        k=0
        a1=-3.4219d-3
        a2=6.7056d-2
        do i=1,kdlvl,jdl
           k=k+1
           pplvl(k)=1.0d0*dexp(a1*(kdlvl-i)*(1.0d0+a2*(kdlvl-i)))
           if(pplvl(k).gt.psrfc)then
              goto 120
           endif
        enddo
120     continue
        pplvl(k)=psrfc
        ntau=k-1
        il=1
c       fr(1)=tr*ppod(1)
c       fn(1)=taum(1)
c       fa(1)=taua(1)
c
c       assign a value of 1.0e-10 fr(1),fn(1) and fa(1) (top of the atm.)
c
        if(taur(nd).gt.0.0d0)fr(1)=1.0d-10
        if(taum(nd).gt.0.0d0)fn(1)=1.0d-10
        if(taua(nd).gt.0.0d0)fa(1)=1.0d-10
        if(ifc.eq.0)fn(1)=0.0d0
c
        htlvl(1)=htd(1)
        ftot(1)=fr(1)+fn(1)+fa(1)
        write (2,1110)il,htlvl(1),pplvl(1),fr(1),fn(1),fa(1),ftot(1)
 
           ohtlvl(1) = htlvl(1)
           opplvl(1) = pplvl(1)
           ofr(1) = fr(1)
           ofn(1) = fn(1)
           ofa(1) = fa(1)
           oftot(1) = ftot(1)
 
        do i=2,k
           pp=pplvl(i)
           pplg=dlog(pp)
           call linear (pplg,lgpdv,htdv,nd,htlvl(i))
           if (iair(1).eq.1) call linear (pp,pdv,taur,nd,fr(i))
           if (iair(2).eq.1) call linear (pp,pdv,taum,nd,fn(i))
           if (iair(3).eq.1) call linear (pp,pdv,taua,nd,fa(i))
           ftot(i)=fr(i)+fn(i)+fa(i)
           write (2,1110)i,htlvl(i),pplvl(i),fr(i),fn(i),fa(i),ftot(i)
 
           ohtlvl(i) = htlvl(i)
           opplvl(i) = pplvl(i)
           ofr(i) = fr(i)
           ofn(i) = fn(i)
           ofa(i) = fa(i)
           oftot(i) = ftot(i)
 
        enddo
      endif
c
c                                                                       
      do i=1,ntau                                                   
         dtrr(i) = fr(i+1) - fr(i)                                         
         dtmm(i) = fn(i+1) - fn(i)
         dtaa(i) = fa(i+1) - fa(i)                                         
         if(dtrr(i).lt.0.0d0)dtrr(i)=0.0d0
         if(dtmm(i).lt.0.0d0)dtmm(i)=0.0d0
         if(dtaa(i).lt.0.0d0)dtaa(i)=0.0d0
      enddo
      write(2,1130)                                                     
c                                                                       
      write (2,1190)                                                    
      write (2,1200)                                                    
c                                                                       
      trp=0.0                                                           
      tmp=0.0                                                           
      tap=0.0                                                           
      nolyr=ntau
      do i=1,ntau                                                  
         trp=trp+dtrr(i)                                                   
         tmp=tmp+dtmm(i)                                                   
         tap=tap+dtaa(i)                                                   
         tots = dtrr(i) + dtmm(i)*const 
         dtot(i) = dtrr(i) + dtmm(i) + dtaa(i)                    
         salb = tots/dtot(i)                                               
         turbhl=(dtmm(i)*const)/tots
         write (2,1210) i,pplvl(i+1),dtrr(i),dtmm(i),dtaa(i),dtot(i),
     1                  salb,turbhl 
 
        opplvl(i) = pplvl(i+1)
        odtrr(i) = dtrr(i)
        odtmm(i) = dtmm(i)
        odtaa(i) = dtaa(i)
        odtot(i) = dtot(i)
        osalb(i) = salb
        oturbhl(i) = turbhl
 
      enddo
c                                                                       
      tautot = trp + tmp + tap                                          
c                                                                       
      write(2,1230)trp,tmp,tap,tautot                                   
 
      otrp = trp
      otmp = tmp
      otap = tap
      otautot = tautot
 
       write(23,1022)
       write(23,1023)wvlth,psrfc,rrho,ozamtp,tautot,dtau1,nolyr
 
       owvlth = wvlth
       opsrfc = psrfc
       orho = rrho
       oozamtp = ozamtp
       otautot = tautot
       odtau1 = dtau1
       onolyr = nolyr
 
       write(23,1024)
       write(23,1025)trp,tmp,tap,tcar,twat,toznp,ifc
 
       otrp = trp
       otmp = tmp
       otap = tap
       otcar = tcar
       otwat = twat
       otoznp = toznp
       oifc = ifc
 
       write(23,1026)
       do i=1,nolyr
          write(23,1027)i,htlvl(i+1),pplvl(i+1),dtrr(i),dtmm(i),
     1                  dtaa(i),dtot(i)
 
        ohtlvl(i) = htlvl(i+1)
        opplvl(i) = pplvl(i+1)
        odtrr(i) = dtrr(i)
        odtmm(i) = dtmm(i)
        odtaa(i) = dtaa(i)
        odtot(i) = dtot(i)
 
       enddo
c
      return                                                            
c                                                                       
c................format statements......................................
c                                                                       
 115  format(i4,1x,1p6e12.4)                                            
 121  format(t1,'tot_lvl')
 122  format(i5)
 128  format(t3,'lvl',t9,'ht(km)',t21,'p(atm)',t33,'taur',t45,'taum',
     1       t57,'taua',t69,'tau_tot')
 1000 format (15x,'total optical thickness at 3.6 km = ',1pe15.4/       
     1 15x,'new dtau=',1pe15.4/ 15x, 'new no. of layers=',i5//)         
 1010 format(/,10x,'number of layers above the cloud = ',i5,10x,        
     1 'layer thickness = ',f10.6)                                      
 1011 format(/,10x,'number of layers below the cloud = ',i5,10x,        
     1  'layer thickness = ',f10.6)                                     
 1012 format(/,10x,'number of layers in the cloud = ',i5,10x,           
     1  'layer thickness = ',f10.6)                                     
 1015 format (1x,t50,'clear atmosphere')                                
 1016 format (1x,t50,'thin cloud')                                      
 1017 format (1x,t50,'thick cloud')                                     
 1020 format (/,t50,'cumulative optical thicknesses'///t11,             
     1 'press.',t27,'rayleigh', t47, 'mie',t67,'ozone',t87,'total'/     
     2 t11,'(in atm)'/)                                                 
 1022 format(t5,'wvlth',t17,'psrfc',t29,'rho',t41,'xozn',t53,'tautot',
     1       t65,'deltau',t74,'nolyr')
 1023 format(1p6e12.4,i4)
 1024 format(t5,'tr',t17,'tmt',t29,'ta',t41,'tcar',t53,'twat',
     1       t65,'tozn',t75,'ifc')
 1025 format(1p6e12.4,i4)
 1026 format(t3,'no',t9,'htl_b',t21,'pl_b',t33,'dtrr',t45,'dtmm',
     1       t57,'dtaa',t69,'dtot')
 1027 format(i4,1p6e12.4)
 1030 format (5x,'delta tau =',1pe15.5,5x,'no. of layers =',i3/)   
 1040 format ('interpolated values of press  and',
     1 ' tau(rayleigh) etc. corresponding to'/'delta tau=',1pe15.5//     
     2 t10,'height',t23,'press.',t31,'tau(rayleigh)',t46, 'tau(mie)',t56,
     3 'tau(ozone)',t68,'tau(tot)'/t10,'(in km)',t22,'(in atm)'/)  
 1050 format ('values of tau(rayleigh) etc. in full and half layers'/
     1'(the first and the last layer are half layers)'/)
 1060 format (//t15,'division of atmosphere into layers of equal optical
     1 thickness for pupose of calculating radiative transfer:'//)      
 1070 format (15x,'total ozone(top to psrfc) in atm-cm',1pe15.5/)      
 1080 format (i4,1x,1p6e12.4)                                           
 1090 format (i3,2x,1p4e11.3,0pf8.4)                                          
 1100 format (t2,'no',t8,'dtau(ray)',t19,'dtau(mie)',t30,           
     1 ' dtau(o3)',t41,'dtau(tot)',t53,'turb'/)  
 1110 format (i4,1x,1p6e12.4)                                           
 1120 format (3e15.5)                                                   
 1130 format ('    ')                                                 
 1140 format (t10,'sum of the total dtau = ',f10.5,                 
     1 /t10,'sum of the total dtaur       = ', f10.5,                          
     2 /t10,'sum of the total dtaum(scat.)= ', f10.5,                          
     3 /t10,'sum of the total dtaua       = ', f10.5)                          
 1190 format(/t5,'values of tau(rayleigh ) etc. in each layer'/) 
 1200 format(t2,'no',t6,'press(atm)'t17,'dtau(ray)',t28,                
     1 'dtau(mie)',t39, 'dtau(o3)',t49,' dtau(tot)',              
     2  t60,'alb_sscat',t71,'turb'/)                                              
 1210 format(i3,1p5e11.3,0p2f8.4)                                            
 1230 format(/t10,'tau_r        =',e12.4/t10,'tau_m        =',e12.4/ 
     1       t10,'tau_a        =',e12.4/t10,'tau_t        =',e12.4) 
      end                                                               
c                                                                       
c......subroutine part2..............................................
c                                                                       
      subroutine part2(nm,pdv,ppod,htd,znp,zpp,qtt,tmt,taum)
c                                                                       
      include 'afrt_rt1.cmn'
      real*8 sumxx(nsc+1),hxx(nsc+1),e(nsc+1),u(nsc+1)
      real*8 ttozn(klam),th2o(klam),tco2(klam),ppod(kmac),htd(kmac),znp(kmac),
     1       dxd(kmac),watr(kmac),dx(kozn),x(koznp),pdv(kdlvl),htdv(kdlvl),
     2       dtrr(klyr),dtmm(klyr),abozn(klam), dtaa(klyr),dtot(klyr),
     3       temp1(klyr),xxd(kmacp),tozn101(koznp),xd(kmac),temp2(klyr)
      real*8 wvlth,tr,tmt,ta,twat,tcar,rrho,oznabs,qst,qtt,
     1       xozn,tautot,dtau1,tozn,p2km,h2km,cc0,cc1,cc2
      real*8 fr(klyr),fn(klyr),fa(klyr),ftot(klyr),taur(klyr),
     1       taum(klyr),taua(klyr),tauw(klyr),tauc(klyr),
     2       tauozn(klyr),tau(klyr),lgppod(kmac),lgpdv(kdlvl),
     3       htlvl(klyr),pplvl(klyr)
      real*8 ht(kdlvl),hho(kdlvl),zpp(kdlvl),zppt(kdlvl),dnm(kdlvl)
      integer*4 n,nm,ifc,nolyr,kndx,ndvsrf,ndv
c                                                                       
c     determine the geo. height and num. density corresponding to pdv
c                                                                       
c     write(*,*)'welcome to sub part2'
c     write(*,*)'pdv(1)',pdv(1)
      ndv=kdlvl
      do i=1,ndv
         do j=1,nm-1
            if(pdv(i).ge.ppod(j) .and. pdv(i).lt.ppod(j+1))then
              xn1=dlog(pdv(i))-dlog(ppod(j))
              xd1=dlog(ppod(j+1))-dlog(ppod(j))
              hho(i)=htd(j)+(xn1/xd1)*(htd(j+1)-htd(j))
              prod1=(hho(i)-htd(j))/(htd(j+1)-htd(j))
              zpp(i)=znp(j)*(znp(j+1)/znp(j))**prod1
            endif
         enddo
         hho(ndv)=0.0d0
         zpp(ndv)=znp(nm)
c     write(6,226)i,pdv(i),hho(i),zpp(i)
226   format('i,pdv,hho,zpp',i3,1x,1p3e11.3)
      enddo
c
c     write(*,*)'pdv(1)',pdv(1)
c     integrate the number density 
c                                                                       
      zppt(1)=0.0d0
      do i=1,ndv-1
         dnm(i)=0.5d0*(zpp(i)+zpp(i+1))*(hho(i)-hho(i+1))*1.0d+5
         zppt(i+1)=zppt(i)+dnm(i)
c        write(6,151)i+1,zppt(i+1)
151      format('integrated # of part i+1,zppt(i+1)',i3,1pe11.3)
      enddo
c
c....convert the total number of particles in to optical depth
c
c     tp=zppt(ndv)*qtt
c     normalize the total aerosol optical thickness at psrfc to the
c     input values 
      tp=zppt(ndvsrf)*qtt
c
      if(tmt.le.10.0d0)then
c     if(tmt.le.0.2d0)then
c       scale the optical depth 'tp' to the input value 'tmt'
c                                                                       
        c = tmt/tp
c                                                                       
        do  i=1,ndv                                                       
            taum(i) = c * zppt(i) * qtt                                         
        enddo
        return
      else
c
c     scale the optical depth 'tp' to a value of 0.2 and determine
c     total aerosol optical thickness up to 2km
c
        c=0.20d0/tp
        do i=1,ndv
           zpp(i)=zpp(i)*c
           zppt(i)=zppt(i)*c
           taum(i)=zppt(i)*qtt
c          write(6,119)i,zppt(i),taum(i)
119        format('i,zppt(i),taum(i)',i3,1x,1p2e11.3)
        enddo
c
c
c       total aerosol opt and numbers from 2km to the ground
        t2b=tmt-taum(kndx)
        xn2b=t2b/qtt
        yn2b=zppt(ndv)-zppt(kndx)
c
c       determine the scale height between 0 and 2km in the input profile
c
        h0=2.0d0/(dlog(znp(nm)/znp(nm-2)))
c
c       determine the low and high value of scale height that bracket
c       the total number of particles 'xn2b' from 0-2km
c
        hx=h0*1.1d0
100        continue
           sum2x=0.0d0
           do j=kndx+1,ndv
              zpp(j)=zpp(kndx)*dexp((hho(kndx)-hho(j))/hx)
              dnm(j-1)=0.5d0*(zpp(j-1)+zpp(j))*(hho(j-1)-hho(j))*1.0d+5
              sum2x=sum2x+dnm(j-1)
           enddo
           if(sum2x.gt.xn2b)then
              hg=hx
              goto 200
           endif
           hl=hx
           hx=hx*0.90d0
           if(hx.lt.0.10d0)then
c            write(6,600)hx
600          format('the value of the aerosol scale height is less',
     1              ' than 0.10',1pe11.3)
             stop
           endif
           goto 100
200        continue
c          write(6,*)'h0,hl,hx',h0,hl,hx
c          compute sumxx as a function of scale height hxx
c
           nscp=nsc+1
           dhx=(hg-hl)/dfloat(nsc)
c
           do i=1,nscp
              hx=hl+(i-1)*dhx
              sum2x=0.0d0
              do j=kndx+1,ndv
                 zpp(j)=zpp(kndx)*dexp((hho(kndx)-hho(j))/hx)
                 dnm(j-1)=0.5d0*(zpp(j-1)+zpp(j))*(hho(j-1)-hho(j))*
     1                    1.0d+5
                 sum2x=sum2x+dnm(j-1)
              enddo
              sumxx(i)=sum2x
              hxx(i)=hx
           enddo
c
c          interpolate for scale height
c
           in=1
           il=1
           iu=1
           call spline(xn2b,sumxx,hxx,nscp,in,hxz,il,iu,vl,vu,e,u)
c          
           do j=kndx+1,ndv
              zpp(j)=zpp(kndx)*dexp((hho(kndx)-hho(j))/hxz)
              dnm(j-1)=0.5d0*(zpp(j-1)+zpp(j))*(hho(j-1)-hho(j))*
     1                 1.0d+5
              zppt(j)=zppt(j-1)+dnm(j-1)
              taum(j)=zppt(j)*qtt
           enddo
c
      endif
c
c
c
      return
      end                                                               
c...............subroutine water...............................
c                                                                       
       subroutine water(twat,pdv,watr,ppod,ht,nm,ndv,tauw)                  
c                                                                       
      implicit real*8 (a-h,o-z)
      real*8 pdv(1000),watr(1000),ppod(1000),tauw(1000),u(1000)        
      real*8 dts(1000),dt(1000),tauww(1000),e(1000),ht(1000)
c                                                                       
c     convert the water vapor density from gm/m**3 to gm/cm**2/km
c     (see mcclachy 72 p89)
      do i=1,nm
         watr(i)=watr(i)*1.0d-1
      enddo
c                                                                       
c     determine the total water amount 
      nm1 = nm - 1                                                       
      tp = 0.0d0
      do 50 i=1,nm1                                                     
         tp = tp + (watr(i)+watr(i+1))/2.0d0 
   50 continue                                                          
c
c     scale the water vapor amount 'tp' to the input optical thickness
c     value 'twat'
c                                                                       
      c = twat/tp                                                          
c                                                                       
      do 60 i=1,nm                                                       
       dts(i) = c * watr(i)
   60 continue                                                          
      do 70 i=1,nm1                                                     
       dt(i) = (dts(i+1) + dts(i)) / 2.0d0
   70 continue                                                          
      tauww(1) = 0.                                                     
      do 80 i=1,nm                                                       
c.....cumulative water vapor optical depth 
       tauww(i+1) = tauww(i) + dt(i)                                    
   80 continue                                                          
      in = 1                                                            
      il = 1                                                            
      iu = 1                                                            
c                                                                       
      do 85 i=1,ndv                                                      
      call spline (pdv(i),ppod,tauww,nm,in,tauw(i),il,iu,vl,vu,e,u)      
      ttww=dabs(tauw(i))
      if(tauw(i).lt.0.0d0 .and. ttww.lt.1.0d-8)then
         tauw(i)=0.0d0
      endif
   85 continue                                                          
      return                                                            
      end                                                               
c***********************************************************************
      subroutine ozn3
c
      include 'afrt_rt1.cmn'
c
      real*8 ee(klyr),uu(klyr)
      real*8 ppod(kmac),pdv(kdlvl),tauozn(klyr),oznabs
      integer*4 n,ndvsrf
c**********************************************************************
      in=1
      il=1
      iu=1
      do i=1,ndvsrf
         call spline(pdv(i),ppod,xd,n,in,tempo3,il,iu,vl,vu,ee,uu)
         tauozn(i)=tempo3*oznabs
      enddo
c
      return
      end
c***********************************************************************
c
      subroutine co2(tcar,pdv,tauc,ndv)
c
      implicit real*8 (a-h,o-z)
      include 'afrt_rt1.cmn'
c
      real*8 pdv(kdlvl),tauc(kdlvl)
c
c     normalize the surface pressure to tcar
c
      cc=tcar/pdv(ndv)
c
      do i=1,ndv
         tauc(i)=pdv(i)*cc
      enddo
c
      return
      end
c***********************************************************************
      subroutine linear(xp,x,y,n,yp)                                    
c***********************************************************************
c     subroutine performs linear interpolation.  It assumes that        
c     x is a monotonically increasing parameter                         
c***********************************************************************
      implicit real*8 (a-h,o-z)
      real*8 x(n),y(n)                                                  
c                                                                       
      il=1                                                              
      ih=n                                                              
100   continue                                                          
      m=(il+ih)/2                                                       
      if(xp.gt.x(m))then                                                
         il=m                                                           
      else                                                              
         ih=m                                                           
      endif                                                             
      id=ih-il                                                          
      if(id.eq.1)goto 200                                               
      goto 100                                                          
200   continue                                                          
c                                                                       
      slope=(y(ih)-y(il))/(x(ih)-x(il))                                 
      yp=y(il)+slope*(xp-x(il))                                         
c                                                                       
      return                                                            
      end                                                               
c**********************************************************************
c......subroutine part1..............................................
c
      subroutine part1(nm,pdv,ppod,htd,znp,zpp,qtt,tmt,taum)
c
      include 'afrt_rt1.cmn'
      integer*4 nm
      real*8 oznabs,qtt,tmt
      real*8 ppod(kmac),htd(kmac),znp(kmac),taum(klyr),pdv(kdlvl)
      real*8 ht(kdlvl),hho(kdlvl),zpp(kdlvl),zppt(kdlvl),dnm(kdlvl)
c
c     determine the geo. height and num. density corresponding to pdv
c
c     write(6,*)'welcome to part1'
      ndv=kdlvl
      do i=1,ndv
         do j=1,nm-1
            if(pdv(i).ge.ppod(j) .and. pdv(i).lt.ppod(j+1))then
              xn1=dlog(pdv(i))-dlog(ppod(j))
              xd1=dlog(ppod(j+1))-dlog(ppod(j))
              hho(i)=htd(j)+(xn1/xd1)*(htd(j+1)-htd(j))
              prod1=(hho(i)-htd(j))/(htd(j+1)-htd(j))
              if(znp(j).gt.0.0)then
                 zpp(i)=znp(j)*(znp(j+1)/znp(j))**prod1
              else
                 zpp(i)=0.0d0
              endif
            endif
         enddo
         hho(ndv)=0.0d0
         zpp(ndv)=znp(nm)
c     write(6,226)i,pdv(i),hho(i),zpp(i)
226   format('i,pdv,hho,zpp',i3,1x,1p3e11.3)
      enddo
c
c     integrate the number density
c
      zppt(1)=0.0d0
      do i=1,ndv-1
         dnm(i)=0.5d0*(zpp(i)+zpp(i+1))*(hho(i)-hho(i+1))*1.0d+5
         zppt(i+1)=zppt(i)+dnm(i)
      enddo
c
c....convert the total number of particles in to optical depth
c
      tp=zppt(ndv)*qtt
c
c       scale the optical depth 'tp' to the input value 'tmt'
c
        c = tmt/tp
c
        do  i=1,ndv
            taum(i) = c * zppt(i) * qtt
        enddo
        return
      end
c*********************************************************************
c
      subroutine dust_prof(znpdd,htdd)
c
c*********************************************************************
c
      include 'afrt_rt1.cmn'
c
      real*8 znpdd(kmac),htdd(kmac)
 
c      write(*,*)'idust',idust
      if(idust.eq.2)then
        top_cld=ht_dust 
        btm_cld=sigma_dust
        do i=1,nm
           if(htdd(i).le.top_cld .and. htdd(i).gt.btm_cld)then
              znpdd(i)=1.0d0
           else
              znpdd(i)=0.0d0
           endif
        enddo
      else
        do i=1,nm
           znpdd(i)=exp(-0.5d0*((htdd(i)-ht_dust)/sigma_dust)**2)
        enddo
      endif
c
      return
      end
c***************************************************************************
      subroutine ozn3_abs(temp101,ppo)
c
      include 'afrt_rt1.cmn'
c
      real*8 ee(klyr),uu(klyr),temp101(koznp),ppo(koznp),tozn101(koznp),
     1       tauozn(klyr),dx(kozn),x(kozn),pdv(kdlvl)
c     determine the ozone optical thickness using temperature 
c     dependent coefficient. if ioznab=1 then set temp. at all
c     levels to 227.16k (-46.0 c)
c 
         if(ioznab.eq.1)then
            do i=1,n
               temp101(i)=227.16
            enddo
         endif
c
         alfatemp=cc0+cc1*(temp101(1)-273.16)+
     1            cc2*(temp101(1)-273.16)**2.0
         tozn101(1)=alfatemp*dx(1)
c        write(*,*)'alfatemp,tozn101(1)',alfatemp,tozn101(1)
      do i=2,n
         alfatemp=cc0+cc1*(0.5d0*(temp101(i)+temp101(i+1))-273.16)+
     1            cc2*(0.5d0*(temp101(i)+temp101(i+1))-273.16)**2.0
         tozn101(i)=tozn101(i-1)+alfatemp*dx(i) 
      enddo
c
      in=1
      il=1
      iu=1
      do j=1,n
         write(2,100)j,ppo(j),dx(j),x(j),tozn101(j)
100      format('i,ppo,dx,x,tozn101',i4,1x,1p4e11.3)
      enddo
      do i=1,ndvsrf
         call spline(pdv(i),ppo,tozn101,n,in,tauo3z,il,iu,vl,vu,ee,uu)
         tauozn(i)=tauo3z
      enddo
c
      return
      end
c***********************************************************************
                                                                        
c**********************************************************************
      subroutine cloud_prof(pcldtop,pcldbtm,ppodd,htdd,znpdd)
c
      implicit real*8 (a-h,o-z)
      include 'afrt_rt1.cmn'
c
      real*8 ppodd(kmac),znpdd(kmac),htdd(kmac),pcldtop,pcldbtm,hcldtop,
     1       hcldbtm,temp1(klyr),temp2(klyr),qxldtoplg,qcldbtmlg
      integer*4 i,k,nm
c**********************************************************************
c     compute the height corresponding to the top and
c     bottom of the cloud layer
c     reverse the order of pressure and height. in linear, x should
c     be monotonically increasing in value.
 
      do i=1,nm
         k=nm-i+1
         temp1(i)=dlog(ppodd(k))
         temp2(i)=htdd(k)
      enddo
      qcldtop=(pcldtop*1.0d3)
      qcldbtm=(pcldbtm*1.0d3)
      qcldtoplg=dlog(qcldtop)
      qcldbtmlg=dlog(qcldbtm)
      call linear(qcldtoplg,temp1,temp2,nm,hcldtop)
      call linear(qcldbtmlg,temp1,temp2,nm,hcldbtm)
c
c     write(*,*)'hcldtop,hcldbtm',hcldtop,hcldbtm
c
      diff=dabs(ppodd(1)-qcldtop)
c     write(*,*)'ppodd(1),qcldtop,diff',ppodd(1),qcldtop,diff
      do i=1,nm
         diffx=dabs(ppodd(i)-qcldtop)
         if(diffx.le.diff)then
            kndx=i
            diff=diffx
         endif
      enddo
      ppodd(kndx)=qcldtop
      diff=dabs(ppodd(1)-qcldbtm)
      do i=1,nm
         diffx=dabs(ppodd(i)-qcldbtm)
         if(diffx.le.diff)then
            kmdx=i
            diff=diffx
         endif
      enddo
      ppodd(kmdx)=qcldbtm
c
c     reset the number density
c
      do i=1,nm
      znpdd(i)=0.0d0
         if(i.ge.kmdx .and. i.le.kndx)znpdd(i)=1.0d0
c     write(*,*)'i,znpdd(i)',i,znpdd(i)
      enddo
c     write(*,*)'kndx,kmdx',kndx,kmdx
c     write(*,*)'qcldtop,qcldbtm',qcldtop,qcldbtm
c
      return
      end
c***********************************************************************