program titwrite.f 

        implicit double precision (o-z, a-h)
        include "parameters.dat"
        character fname*64, iline*64, iret*1

c       ---------------------------------------------------------
c            Arrays : Model atmosphere
c       ---------------------------------------------------------
        dimension zalt(nd3),pres(nd3),temp(nd3),dens(nd3,0:4)
     &   ,col(nd3,0:4),ifund(4),vsum(nd4,2)

c       ---------------------------------------------------------
c            Arrays : IR band parameters
c       ---------------------------------------------------------
     &  , arad(12,3),pcol(12,3),dgn(12,3),rdmol(3),rmcol(3)
     &  , rmass(3), qrot(3), wave(12,3)

c       ---------------------------------------------------------
c            Arrays for R-T calculation
c       ---------------------------------------------------------
     &  , tmp(nd3), prs(nd3), den(nd3), z(nd3), clm(nd3)
     &  , sr(nd3), source(nd3), wavs(1000), signal(1000)
     &  , cf(163,11), op_sum(nd3)

c       ---------------------------------------------------------
c                           Data for Titan 
c       ---------------------------------------------------------
        data rmn2 / 28.d0 /
        data rdn2 / 3.681d0 /
        data gm / 8.95d18 /
        data GMm / 4.16d-4 /
        data rplnt / 2575.d5 /
        data nlortz/3000/
        data wlortz/5.2d-2/
        data wdel/10.d0/
        data pmin / 1.d-5 /
        data pmax / 1.d+2 /
        data rtmax / 0.1d0 /
        data ifund/1,5,7,8/

c       --------------------------------------------------------
c              Read Yung and Allen model atmosphere file
c       --------------------------------------------------------
        open(unit=60,file='titan.atm',status='old')

           read(unit=60,fmt=*) ndat
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (zalt(n),n=1,ndat)
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (pres(n),n=1,ndat)
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (temp(n),n=1,ndat)
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (dens(n,0),n=1,ndat)
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (col(n,0),n=1,ndat)
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (dens(n,1),n=1,ndat)
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (col(n,1),n=1,ndat)
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (dens(n,2),n=1,ndat)
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (col(n,2),n=1,ndat)
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (dens(n,3),n=1,ndat)
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (col(n,3),n=1,ndat)
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (dens(n,4),n=1,ndat)
           read(unit=60,fmt=900) iline
           read(unit=60,fmt=901) (col(n,4),n=1,ndat)
        close(unit=60)

        natm = locate(pres,ndat,pmax) + 1

 900    format(a)
 901    format(1x,1p10e11.3)
 902    format(1x,i6)

        ntau = 51
        dlmin = dlog10(pmin)
        dlmax = dlog10(pmax)
        dlstep = (dlmax - dlmin) / float(ntau - 1)

        do nt = 1, ntau 

           dum1 = dlmin + (float(nt - 1) * dlstep)

           prs(nt) = 10.d0 ** dum1
           if (prs(nt) .le. pres(1)) then
              tmp(nt) = temp(1)
              den(nt) = dens(1,1)*prs(nt)/pres(1)
              clm(nt) = col(1,1)*prs(nt)/pres(1)
           else
              tmp(nt) = fintrp(pres,temp,natm,prs(nt))
              den(nt) = fintrp(pres,dens(1,1),natm,prs(nt))
              clm(nt) = fintrp(pres,col(1,1),natm,prs(nt))
              i1 = locate(pres,natm,prs(nt))
              i2 = i1 + 1
              pl1 = dlog(pres(i1))
              pl2 = dlog(pres(i2))
              pl = dlog(prs(nt))
              z(nt) = zalt(i1)+(pl-pl1)*(zalt(i2)-zalt(i1))
     &          /(pl2-pl1) 
           end if

        end do

c       ---------------------------------------------------------
c              Read mesosphere/thermosphere source functions
c       ---------------------------------------------------------
        open (unit=30, file='source_a.out', status='old') 

           do n = 1, ntau
              read(unit=30, fmt=*) zdum,pdum,sr(n)
           end do

        close(unit=30)

        open (unit=30, file='heat_a.out', status='old') 

           do n = 1, ntau
              read(unit=30, fmt=*) zdum,tmp(n)
           end do

        close(unit=30)

c       ---------------------------------------------------------
c           Splice them all together to get model atmosphere from
c           ground to exobase
c       ---------------------------------------------------------
        m = ntau
        nbreak = m
        do n = natm , ndat
           m = m + 1
	   z(m) = zalt(n)
           tmp(m) = temp(n)
           prs(m) = pres(n)
           sr(m) = one
           den(m) = dens(n,1)
        end do
        ntau = m

        clm(1) = 1.d10
        do n = 2, ntau
           clm(n) = clm(n-1)+half*(den(n-1)+den(n))*(z(n-1)-z(n))
        end do

        open(unit=50,file='dump.out',status='old')
           do n = 1, ntau
           write(unit=50,fmt='(1x,2f7.2,3e11.3)') 1.e-5*z(n),tmp(n)
     &      ,prs(n),clm(n),sr(n)
           end do
        close(unit=50)

c ********************************************************************
c *                 Calculate Methane nu-4 spectrum                  *
c ********************************************************************

c      -------------------------------------------------------------
c                       Read spectral line list
c      -------------------------------------------------------------
       call sread(1,wav,str,wid,eng,nlines,stot,wmean
     &  ,wstart,wend)

       write(unit=*,fmt=952) stot,nlines,wstart,wend
 952   format(' band strength = ',1pe12.4/
     &   ,' number of lines = ',i8/
     &   ,' wavenumber range = ',1p2e12.4)

c      ----------------------------------------------------
c                 Set up coarse wavenumber gird
c      ----------------------------------------------------
       write(unit=*, fmt=904) 
       read(unit=*, fmt=*) wstart, wend, wdel
       nwaves = ((wend - wstart) / wdel) + 0.4999d0
       wavs(1) = wstart + (half * wdel)
       do l = 2, nwaves
          wavs(l) = wavs(l - 1) + wdel
       end do

       write(unit=*, fmt='(33h Enter cosine of emission angle >,$)')
       read (unit=*, fmt=*) ue

c      -------------------------------------------------------
c                      Loop through coarse bins
c      -------------------------------------------------------
       do nt = 1, ntau
          op_sum(nt) = zero
       end do

c       do nt = 1, ntau
c          tmp(nt) = 175.d0
c       end do

       do l = 1, nwaves

          wmin = wavs(l) - half*wdel
          wmax = wmin + wdel

c         ----------------------------------------------------
c                   Calculate source functions
c         ----------------------------------------------------
          do n = 1, ntau
             source(n) = sr(n)*bplanck(wavs(l),tmp(n))
          end do

c         ----------------------------------------------------
c                    Calculate emergent intensity
c         ----------------------------------------------------
          tref = 296.d0
          wfc = wavs(l)*dsqrt(two*rkb/rmass(1)/amu)/cs
          wstep = wfc*dsqrt(tref)/two
          wlc = 1.d-6/dsqrt(tref)/wfc

          do nt = 1, ntau

             wdop(nt) = wfc*dsqrt(tmp(nt))
             wlor(nt) = wlc*prs(nt)
             wreach(nt) = nlortz*wlortz*wlor(nt)*wdop(nt)
     &        +10.d0*wdop(nt)
             hckt(nt) = hck*(one/tmp(nt) - one/tref)
             qw(nt) = (tref/tmp(nt))**efc/wdop(nt)

          end do

c         -----------------------------------------------------------
c           Find the range of lines needed to calculate absorption
c           coefficient at the first frequency
c         -----------------------------------------------------------

          do nt = 1, ntau
             w1 = wmin - wreach(nt)
             lw1(nt) = locate(wav,nlines,w1)
             w2 = wmin + wreach(nt)
             lw2(nt) = locate(wav,nlines,w2)+1
          end do

c          lwstart = min(locate(wav,nlines,wavs(l)-5.d0)+1,nlines)
c          lwend = max(locate(wav,nlines,wavs(l)+5.d0),1)
c         --------------------------------------------------------------
c                           Fine frequency loop 
c         --------------------------------------------------------------
          lw = 2200
          write(unit=*,fmt=803) wmin,wav(lw),wmax
 803      format(1x,3e15.7)
          read (unit=*,fmt=900) iret
          do nt = 1, ntau

             wlast = wmin
             vglast = zero
             sum = zero
             do while(wlast .lt. wmax)

                wnow = dmin1(wlast + wstep, wmax)
                wmid = half*(wnow+wlast)
                delw = wnow - wlast
                dw = dabs(wmid-wav(lw))
                if(dw .lt. wreach(nt)) then
                   xw = dw / wdop(nt)
                   avgt = wid(lw)*wlor(nt)
                   call voigt(xw,avgt,vgnow)
                   sum = sum + delw*half*(vgnow+vglast)
                   vglast = vgnow
                end if
                wlast = wnow

             end do
             sum = sum/sqpi/wdop(nt)
             write(unit=*,fmt=802) nt,sum

          end do
          read(unit=*,fmt=900) iret
 802      format(1x,i5,e13.5)

          lw1x = lw1(1)
          do lw = lw1x, nlines 
             vsum(lw,1) = zero
             vsum(lw,2) = zero
          end do

          sum_nu = zero
          wlast = wmin

          do while(wlast .lt. wmax)

             wnow = dmin1(wlast + wstep, wmax)
             wmid = half*(wnow+wlast)
             delw = wnow - wlast

c            --------------------------------------------------------
c              Find spectral lines which contribute to absorption 
c              coefficient at wmid
c            --------------------------------------------------------
             do nt = 1, ntau
                do while(wmid-wav(lw1(nt)) .gt. wreach(nt))
                   lw1(nt) = lw1(nt) + 1
                end do
                do while(wav(lw2(nt))-wmid .le. wreach(nt))
                   lw2(nt) = lw2(nt) + 1
                end do
             end do

c            --------------------------------------------------------
c                         Loop through altitudes
c            --------------------------------------------------------

             trans_last = one
             s_last = source(1)
             p_last = -20.d0 
             tau = zero
             sumz = zero
             nt = 1
             do while ((nt .le. ntau).and.(tau.lt.30.d0))

c               ------------------------------------------------------
c                 Calculate contribution of each line to absorption
c                  coefficient
c               ------------------------------------------------------

c                lw = 2184
c                dw = dabs(wmid-wav(lw))
c                xw = dw / wdop(nt)
c                avgt = wid(lw)*wlor(nt)
c                call voigt(xw,avgt,vg)
c                write(unit=*,fmt=801) nt,wmid,dw
c     &           ,wreach(nt),avgt,vg
c 801            format(i5,5e13.5)

                suml = zero
                do lw = lw1(nt), min(lw2(nt),nlines)
c                do lw = lwstart, lwend 
                   dw = dabs(wmid-wav(lw))
                   if(dw .le. wreach(nt)) then
                      xw = dw / wdop(nt)
                      avgt = wid(lw)*wlor(nt)
                      call voigt(xw,avgt,vg)
                      suml = suml + vg*str(lw)*dexp(-hckt(nt)*eng(lw))
                 
                      if(nt .eq. 1) then 
                      vsum(lw,1) = vsum(lw,1)+vg*delw/sqpi/wdop(1) 
c                        write(unit=*,fmt=*) lw,dw,vg
                      end if
                      if(nt .eq. 60) then 
                      vsum(lw,2) = vsum(lw,2)+vg*delw/sqpi/wdop(60) 
c                        write(unit=*,fmt=*) lw,dw,vg
                      end if

                   end if
                end do
                op_now = qw(nt)*suml/sqpi

                op_sum(nt) = op_sum(nt) + half*(op_now+op_last)*delw

                if(nt .eq. 1) then
                   tau = op_now*clm(1)
                else
                   tau = tau + op_now*(clm(nt)-clm(nt-1))
                end if

c                write(unit=*,fmt=800) lw1(nt),lw2(nt),nt,prs(nt)
c     &           ,op_now,tau,trans_now,s_now
c 800            format(1x,3i5,5e11.3)
                trans_now = dexp(-tau/ue)
                s_now = source(nt)
                p_now = dlog(prs(nt))                  

                sumz = sumz + half*(s_now+s_last)
     &            *(trans_last-trans_now)

                if(tau .gt. 1.d0) nnn = nt
c                cf(nt,l) = cf(nt,l) + half*delw*(s_now+s_last)
c     &            *(trans_last-trans_now)/(p_now-p_last)/(wmax-wmin)

                op_last = op_now
                trans_last = trans_now
                s_last = s_now
                p_last = p_now
                nt = nt + 1
c            --------------------------------------------------
c                          End altitude loop
c            --------------------------------------------------
             end do

c            ----------------------------------------------------
c                          Add in ground term
c            ----------------------------------------------------
             if(nt .eq. ntau) then
                sumz = sumz + source(ntau)*trans_now
             end if

             write(unit=*,fmt=806) wmid,sumz,prs(nnn),tmp(nnn)
 806     format(1x,4e13.5)
c            ----------------------------------------------------
c               Add altitude integral to frequency integration
c            ----------------------------------------------------
             sum_nu = sum_nu + sumz*delw
             wlast = wnow

c         -------------------------------------------------------
c                      End fine frequency loop
c         -------------------------------------------------------
          end do

c          read(unit=*,fmt=900) iret
          do lw = lw1x , lw2(60)
             write(unit=*,fmt=804) lw, vsum(lw,1),vsum(lw,2)
          end do
          read(unit=*,fmt=900) iret
 804      format(1x,i5,' vsum = ',2e13.5)

          signal(l) = sum_nu/(wmax-wmin)
          write(unit=*,fmt=920) wavs(l), signal(l)

c       ---------------------------------------------------------
c                     End coarse frequency loop
c       ---------------------------------------------------------
        end do

c       ---------------------------------------------------------
c                     Normalize contribution function
c       ---------------------------------------------------------
c        do l = 1, nwaves
c           do nt = 1, ntau
c              cf(nt,l) = cf(nt,l)/signal(l)
c           end do
c        end do

c       --------------------------------------------------------
c                    Write output file
c       --------------------------------------------------------
        write(unit=*,fmt=907)
        read (unit=*,fmt=900) fname
        open(unit=60,file=fname,status='new')
           do l = 1, nwaves
              write(unit=60,fmt=909) wavs(l), signal(l)
           end do
        close(unit=60)

c        open(unit=60,file='cf.out',status='old')
c           do n = 1, ntau
c              write(unit=60,fmt='(1x,1p11e11.3)') prs(n)
c     &         ,(cf(n,l),l=1,10)
c           end do
c        close(unit=60)

        open(unit=60,file='op.out',status='old')
           do n = 1, ntau
              write(unit=60,fmt='(1x,i5,0pf7.2,1p2e11.3)') 
     &         n,tmp(n),prs(n),op_sum(n)
           end do
        close(unit=60)

        stop

c       ----------------------------------------------------------
c                              Formats
c       ----------------------------------------------------------
  904   format(' Enter wmin, wmax, and delw > ',$)
  907   format(' Enter filename for output > ',$)
  908   format(1x,i6)
  909   format(1x,0pf9.3,2x,1pe13.5)
  910   format(1x,1p10e11.3)
  912   format(1x,0pf7.1,2x,1pe13.5)
  920   format(' wavenumber = ',0pf9.3,' intensity = ',1pe13.5)
	
	end