program trans 

c **********************************************************************
c *                                                                    *
c *                                                                    *
c *                                                                    *
c *                                                    Roger Yelle     *
c *                                                    01 Nov 1990     *
c *                                                                    *
c **********************************************************************

        implicit double precision (o-z, a-h)
        include "parameters.dat"
        character fname*64
        real*8 nld

c       ---------------------------------------------------------
c            Arrays : Model atmosphere
c       ---------------------------------------------------------
        dimension zalt(nd3),pres(nd3),temp(nd3),dens(nd3,0:4)
     &   ,col(nd3,0:4),fmix(4)

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), wlortz(3)

c       ---------------------------------------------------------
c            Arrays : line listing
c       ---------------------------------------------------------
     &  , wav(nd4),str(nd4),wid(nd4),eng(nd4),nlines(4)
     &  , stot(4),opp0(4),wdopref(4),wmean(4),wstart(4),wend(4)
     &  , wstep(4),wmid(100)

c       ---------------------------------------------------------
c          Arrays: absorption coefficients
c       ---------------------------------------------------------
     &  ,opf(nd1,nd2),wpf(nd2),wp(20000),op(20000,nd1)
     &  ,tsmt(3,100),tlbl(3,100),rsmt(3,100),rlbl(3,100),osmt(3,100)

c       ---------------------------------------------------------
c            Arrays for R-T calculation
c       ---------------------------------------------------------
     &  , z(nd1),tmp(nd1),prs(nd1),den(nd1,0:4),clm(nd1,0:4),dtdc(3)

c       ---------------------------------------------------------
c                           Data for Titan 
c       ---------------------------------------------------------
        data rmn2 / 28.d0 /
        data rdn2 / 3.d0 /
        data gm / 8.95d18 /
        data GMm / 4.16d-4 /
        data rplnt / 2575.d5 /
        data rh / 9.5d0 /
        data nu /3/
        data nlortz/300/
        data wlortz/5.2d-2,1.d-1,1.d-1/
        data wdel/10.d0/
        data pmin / 1.d-5 /
        data pmax / 1.d+2 /
        data fmix / 1.d0, 1.d0, 1.d0, 1.d0/
        data nld/ 2.7d19/

c        write(unit=*, fmt=901) 
c        read (unit=*, fmt=900) fname
        fname='titan.atm'

c        write(unit=*,fmt=902) 
c        read (unit=*,fmt=*) cmin,pmax,ntau
        write(unit=*,fmt='(14h Enter ntau > ,$)')
        read (unit=*,fmt=*) ntau
        nmid = (ntau+1)/2

c        write(unit=*,fmt=903)
c        read (unit=*,fmt=*) tref
        tref = 150.d0

        write(unit=*,fmt=904)
        read (unit=*,fmt=*) otol

c       ----------------------------------------------------------
c              Read in IR vibration transition parameters
c       ----------------------------------------------------------
        call bandread(wave,arad,pcol,dgn,rmass,rdmol,qrot)

c       ----------------------------------------------------------
c          Calculate reduced mass and collision diameters
c           for hydrocarbon-N2 collisions
c       ----------------------------------------------------------
        do im = 1,3 
           rdmol(im) = 1.d-8*half*(rdn2+rdmol(im))
           rmcol(im) = amu*rmn2*rmass(im)/(rmn2+rmass(im))
        end do

c       ----------------------------------------------------------
c             Read in line listing for fundamentals
c       ----------------------------------------------------------
        write(unit=*,fmt='(14h Enter imol > ,$)')
        read (unit=*,fmt=*) imol
        if(imol .eq. 1) then
          jb = 1
        else if(imol .eq. 2) then
          jb = 5
        else if(imol .eq. 3) then
          jb = 7
        else if(imol .eq. 4) then
          jb = 8
        end if
        call sread(jb,wav,str,wid,eng,nlines(imol)
     &   ,stot(imol),wmean(imol),wstart(imol),wend(imol))
        wdopref(imol)=dsqrt(two*rkb*tref/rmass(imol)/amu)
     &   *wmean(imol)/cs
        opp0(imol) = stot(imol)/wdopref(imol) 
        wstep(imol) = wdopref(imol)
        write(unit=*,fmt=810) wstart(imol),wend(imol)
 810    format(' wstart, wend = ',2f9.3)

c       ----------------------------------------------------------
c                      Read in  Model Atmosphere
c       ----------------------------------------------------------
        call atmread(fname,pmax,natm,zalt,pres,temp,dens,col
     &   ,tbot,dtdz)

c       ----------------------------------------------------------
c                            Set up grid 
c       ----------------------------------------------------------
        call cgrid(fmix,gm,rplnt,pmin,pmax,ntau,zalt,temp,pres
     &   ,dens,col,natm,z,tmp,prs,den,clm,dtdz,dtdc)
        clm1 = clm(1,imol)
        clm2 = clm(nmid,imol)
        clm3 = clm(ntau,imol)
c       ----------------------------------------------------------
c                   Calculate absorption coefficients
c       ----------------------------------------------------------
        sum01 = zero
        sum02 = zero
        sum03 = zero
        npsum = 0
        wmin = wstart(imol)
        nb = 1
        do while(wmin .lt. wend(imol))
           wmax = dmin1(wmin+10.d0,wend(imol))
           wmid(nb) = half*(wmin+wmax)
           nmax = nd2
           call smt(nlortz,wlortz(imol),qrot(imol),otol,amu*rmass(imol)
     &     ,wmin,wmax,wmid(nb),wstep(imol),wav,str,eng,wid
     &     ,nlines(imol),tmp,prs,clm(1,imol),ntau,nmax,wpf,opf,npf,wp
     &     ,op,np)

           nsum = nsum + npf
           write(unit=*,fmt=*) ' wmid = ',wmid(nb),' npf = ',npf
     &     ,' nsum = ',nsum

           sum0 = zero
           sum1 = zero
           sum2 = zero
           sum3 = zero
           do n = 1, npf
              sum0=sum0+wpf(n)
              sum1=sum1+wpf(n)*dexp(-opf(1,n)*clm1)
              sum2=sum2+wpf(n)*dexp(-opf(nmid,n)*clm2)
              sum3=sum3+wpf(n)*dexp(-opf(ntau,n)*clm3)
           end do
           tsmt(1,nb) = sum1/sum0
           tsmt(2,nb) = sum2/sum0
           tsmt(3,nb) = sum3/sum0

           sum0 = zero
           sum1 = zero
           sum2 = zero
           sum3 = zero
           do n = 1, np
              sum0=sum0+wp(n)
              sum1=sum1+wp(n)*dexp(-op(n,1)*clm1)
              sum2=sum2+wp(n)*dexp(-op(n,nmid)*clm2)
              sum3=sum3+wp(n)*dexp(-op(n,ntau)*clm3)
           end do
           tlbl(1,nb) = sum1/sum0
           tlbl(2,nb) = sum2/sum0
           tlbl(3,nb) = sum3/sum0

           sum1 = zero
           sum2 = zero
           sum3 = zero
           do n = 1, npf
              sum1=sum1+wpf(n)*opf(1,n)*dexp(-opf(1,n)*clm1)
              sum2=sum2+wpf(n)*opf(nmid,n)
     &         *dexp(-opf(nmid,n)*clm2)
              sum3=sum3+wpf(n)*opf(ntau,n)
     &         *dexp(-opf(ntau,n)*clm3)
           end do
           rsmt(1,nb) = sum1
           rsmt(2,nb) = sum2
           rsmt(3,nb) = sum3

           sum1 = zero
           sum2 = zero
           sum3 = zero
           do n = 1, np
              sum01=sum01+wp(n)*op(n,1)
              sum1=sum1+wp(n)*op(n,1)*dexp(-op(n,1)*clm1)
              sum02=sum02+wp(n)*op(n,nmid)
              sum2=sum2+wp(n)*op(n,nmid)*dexp(-op(n,nmid)
     &         *clm2)
              sum03=sum03+wp(n)*op(n,ntau)
              sum3=sum3+wp(n)*op(n,ntau)*dexp(-op(n,ntau)
     &         *clm3)
           end do
           rlbl(1,nb) = sum1
           rlbl(2,nb) = sum2
           rlbl(3,nb) = sum3

           sum1 = zero
           sum2 = zero
           sum3 = zero
           do n = 1, npf
              sum1=sum1+wpf(n)*opf(1,n)
              sum2=sum2+wpf(n)*opf(nmid,n)
              sum3=sum3+wpf(n)*opf(ntau,n)
           end do
           osmt(1,nb) = nld*sum1
           osmt(2,nb) = nld*sum2
           osmt(3,nb) = nld*sum3

           wmin = wmax
           nb = nb + 1
        end do
        nbmax = nb - 1

        do nb = 1, nbmax
           rlbl(1,nb) = rlbl(1,nb)/sum01
           rlbl(2,nb) = rlbl(2,nb)/sum02
           rlbl(3,nb) = rlbl(3,nb)/sum03
           rsmt(1,nb) = rsmt(1,nb)/sum01
           rsmt(2,nb) = rsmt(2,nb)/sum02
           rsmt(3,nb) = rsmt(3,nb)/sum03
        end do
c       ----------------------------------------------------------
c                        Write output to Disc
c       ----------------------------------------------------------
        open (unit=30, file='opp.out', status='new') 
           do nb = 1, nbmax
           write(unit=30,fmt='(1x,f10.4,2x,3e11.3)') 
     &      wmid(nb),(osmt(i,nb),i=1,3)
           end do
        close(unit=30)

        open (unit=30, file='kernel.out', status='new') 
           do nb = 1, nbmax
           write(unit=30,fmt='(1x,f10.4,2x,6e11.3)') 
     &      wmid(nb),(rlbl(i,nb),i=1,3)
     &      ,((rsmt(i,nb)-rlbl(i,nb))/rlbl(i,nb),i=1,3)
           end do
        close(unit=30)

        stop

c       ----------------------------------------------------------
c                              Formats
c       ----------------------------------------------------------
  900   format(a)
  901   format(' Enter model atmosphere file name > ',$)
  902   format(' Enter cmin, pmax, and ntau > ',$)
  903   format(' Enter reference temperature > ',$)
  904   format(' Enter otol > ',$)
  905   format(' Enter filename for heating rates > ',$)
  906   format(1x,0pf7.1,1x,0pf7.2,2x,1p9e11.3)
  907   format(' Enter filename for source functions > ',$)
  908   format(1x,0pf7.1,1x,1pe10.3,2x,1p3e11.3)
  909   format(' Enter filename for time constants > ',$)
  910   format(1x,0pf7.1,1x,0pf7.2,2x,1p5e10.3)

        end

        subroutine smt(nlortz,wlortz,efc,otol,rmass,wmin,wmax,wmean
     &  ,wdel,wav,str,eng,wid,nlines,tmp,prs,clm,ntau,nmax,wpp,opp
     &  ,ibinmax,wp,op,np)

        implicit double precision (o-z, a-h)
        include "parameters.dat"

        logical itest, istop

c       -------------------------------------------------------------
c                         Input Arrays
c       -------------------------------------------------------------
        dimension tmp(nd1),prs(nd1),clm(nd1),wav(1),str(1),eng(1),wid(1)

c       -------------------------------------------------------------
c                         Output Arrays
c       -------------------------------------------------------------
        dimension opp(nd1,nd2), wpp(nd2)

c      --------------------------------------------------------------
c                       Internal Arrays
c      --------------------------------------------------------------
       dimension wp(20000),op(20000,nd1),wdop(nd1),wlor(nd1)
     &  ,wreach(nd1),hckt(nd1),qw(nd1),lw1(nd1),lw2(nd1),tpt(nd1)
     &  ,tpb(nd1,nd2)

       tref = 296.d0

       do nt = 1, ntau

          wdop(nt) = wmean*dsqrt(two*rkb*tmp(nt)/rmass)/cs
          wlor(nt) = 1.d-6 * prs(nt) * dsqrt(tmp(nt)/tref)/wdop(nt)
          wreach(nt) = dmax1(nlortz*wlortz*wlor(nt)*wdop(nt)
     &     ,five*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) = min(locate(wav,nlines,w1)+1,nlines)
          w2 = wmin + wreach(nt)
          lw2(nt) = max(locate(wav,nlines,w2),1)
       end do

c      --------------------------------------------------------------
c          ***** Everything is ready, start frequency loop *****
c      --------------------------------------------------------------

       wlast = wmin 
       ibinmax = 0
       iw = 0
       do while(wlast .lt. wmax)

          wnow = dmin1(wlast + wdel, wmax)
          wmid = half*(wnow+wlast)
          delw = wnow - wlast
          opmax = zero
          iw = iw + 1
          wp(iw) = delw

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)).and.
     &        (lw1(nt) .lt. nlines))
                lw1(nt) = lw1(nt) + 1
             end do

             do while((wav(lw2(nt))-wmid .le. wreach(nt)).and.
     &        (lw2(nt) .lt. nlines))
                lw2(nt) = lw2(nt) + 1
             end do

          end do

c         -----------------------------------------------------------
c            Loop backwards through altitudes, stop loop once the 
c            top of atmosphere  is reached (nt = 1) or if there are
c            no spectral lines within wreach of wmid
c         -----------------------------------------------------------

          nt = ntau  
          do while((nt .ge. 1).and.(lw1(nt).le.lw2(nt)))

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

             sumo = zero
             lw = lw1(nt) 
             do while(lw .le. lw2(nt))
                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)
                   opv = vg*str(lw)*dexp(-hckt(nt)*eng(lw))
                   sumo = sumo + opv
                end if
                lw = lw + 1
             end do
             op(iw,nt) = qw(nt)*sumo/sqpi

             nt = nt - 1

          end do

          tpt(1) = dexp(-two*op(iw,1)*(clm(2)-clm(1)))
          do nt = 2, ntau-1
             tpt(nt) = dexp(-two*op(iw,nt)*(clm(nt+1)-clm(nt-1)))
          end do
          tpt(ntau) = dexp(-two*op(iw,ntau)*(clm(ntau)-clm(ntau-1)))

c         -----------------------------------------------------------
c                      Find a home for this frequency
c          Loop backwards through bins and stop when a home is found
c         -----------------------------------------------------------

          istop = .true.
          ibin = ibinmax + 1 
          do while((ibin.gt.1).and. istop )

             ibin = ibin - 1

c            --------------------------------------------------------
c              does this frequency fit in any pre-existing bins ?
c            --------------------------------------------------------
              
             itest = .true.
             nt = 0
             do while(itest .and. (nt .lt. ntau))
                nt = nt + 1
                if(dabs(tpt(nt)-tpb(nt,ibin)).gt.otol) itest= .false. 
             end do

c            --------------------------------------------------------
c                 If a group is found, add into bin, redefine weight 
c                  and average
c            --------------------------------------------------------

             if(itest) then

                wsum = wpp(ibin) + delw
                do nt = 1, ntau
                  opp(nt,ibin)=(wpp(ibin)*opp(nt,ibin)+delw*op(iw,nt))
     &             /wsum
                end do
                wpp(ibin) = wsum

                tpb(1,ibin) = dexp(-two*opp(1,ibin)*(clm(2)-clm(1)))
                do nt = 2, ntau-1
                   tpb(nt,ibin) = dexp(-two*opp(nt,ibin)
     &                     *(clm(nt+1)-clm(nt-1)))
                end do
                tpb(ntau,ibin) = dexp(-two*opp(ntau,ibin)
     &                         *(clm(ntau)-clm(ntau-1)))

                istop = .false.

             end if

          end do

c         -----------------------------------------------------------
c                    If it didn't fit create a new bin
c         -----------------------------------------------------------
                
          if(istop) then

             ibinmax = ibinmax + 1
             do nt = 1, ntau
                opp(nt,ibinmax) = op(iw,nt)
             end do
             wpp(ibinmax) = delw

             do nt = 1, ntau
                tpb(nt,ibinmax) = tpt(nt)
             end do

          end if

c         -----------------------------------------------------------
c            Done with this frequency, update wlast and continue
c             with loop
c         -----------------------------------------------------------

          wlast = wnow

       end do
       np = iw

       return
       end