program solar 

c **********************************************************************
c *                                                                    *
c *                                                                    *
c *                                                                    *
c *                                                    Roger Yelle     *
c *                                                    01 Nov 1990     *
c *                                                    Yiping Wang     *
c *                                                    11 Jan 1991     *
c **********************************************************************

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

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 : Solar UV Flux
c       ---------------------------------------------------------
     &  , soluvwave(nd5), soluvflux(nd5), crs(nd5,0:3)
     
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       ---------------------------------------------------------
     &  , wava(nd4),stra(nd4),wida(nd4),enga(nd4),nlinesa(3),nla(3)

c       ---------------------------------------------------------
c            Arrays for R-T calculation
c       ---------------------------------------------------------
     &  , ep1(nd1),ep2(nd1), bj(nd1,3),bpmean(nd1)
     &  , tmp(nd1), prs(nd1), den(nd1,0:4), clm(nd1,0:4), dtdc(3)
     &  , ts(nd1,3), xu(nd6), wu(nd6), z(nd1), ep(nd1)

c       ---------------------------------------------------------
c           Arrays for heating calculation
c       ---------------------------------------------------------
     &  ,  heat(nd1,3), habs(nd1,3), hcsc(nd1,3), huv(nd1,2)

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 /

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

        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
	tbase=175.
	imol=1

        do i = 1, 4
           fmix(i) = one
        end do
c       ----------------------------------------------------------
c           Set up angular quadratures for R-T calculations
c       ----------------------------------------------------------
        call gauleg(zero, one, xu, wu, nu)

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 solar bands 
c       ----------------------------------------------------------
        n = 1
        do ib = 1, 3
           nla(ib) = n
           call sread(ib+1,wava(n),stra(n),wida(n),enga(n)
     &      ,nlinesa(ib),stota,wmeana,wstarta,wenda)
           n = n + nlinesa(ib)
        end do

c       ----------------------------------------------------------
c                      Read in  Model Atmosphere
c       ----------------------------------------------------------
        call atmread2(fname,pmax,natm,zalt,pres,temp,dens,col)

c       ----------------------------------------------------------
c               Set up grid for heating calculation
c       ----------------------------------------------------------
        call cgrid(gm,rplnt,pmin,pmax,ntau,zalt,temp,pres,dens
     &   ,col,natm,z,tmp,prs,den,clm)

c       **********************************************************
c                    Calculate UV heating rate
c       **********************************************************
        call solread(soluvwave,soluvflux,crs,nsol)
        call uvheat(rh,clm,den,ntau,soluvwave,soluvflux,crs,nsol
     &   ,huv(1,1))

        call solread(soluvwave,soluvflux,crs,nsol)
        call uvheat(rh,clm,den,ntau,soluvwave,soluvflux,crs,nsol
     &   ,huv(1,2))

c       **********************************************************
c                 Calculate Methane cascade heating 
c       **********************************************************
        do jt = 1, ntau
           ep1(jt) = epsln(tmp(jt),den(jt,0),arad(1,imol)
     &     ,pcol(1,imol),rdmol(imol),rmcol(imol))
           ep2(jt) = epsln(tmp(jt),den(jt,0),arad(4,imol)
     &     ,pcol(4,imol),rdmol(imol),rmcol(imol))
        end do

        call cascade(nlortz,wlortz(1),qrot(imol),amu*rmass(1)
     &   ,tref,rh,1306.d0,wdel,xu,wu,nu,z,tmp,prs,clm(1,1)
     &   ,den(1,1),ep1,ep2,ntau,wava,stra,enga,wida,nla,nlinesa
     &   ,hcsc,ts,habs)

c       ----------------------------------------------------------
c                      Convert to degrees/day
c       ----------------------------------------------------------
        do n = 1, ntau
           do j = 1, 3
              hcsc(n,j) = (2.5d20/den(n,0))*hcsc(n,j)
           end do
           do j = 1, 2
              huv(n,j) = (2.5d20/den(n,0))*huv(n,j)
           end do
        end do

c       ----------------------------------------------------------
c                        Write output to Disc
c       ----------------------------------------------------------
        write(unit=*, fmt=905) 
        read (unit=*, fmt=900) fname
        open (unit=30, file=fname, status='new') 

           do n = 1, ntau
             write(unit=30,fmt=906) 1.e-5*z(n),prs(n),(hcsc(n,j),j=1,3)
     &        ,(huv(n,j),j=1,2)
           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,2x,0pf6.2,2x,1p5e11.3)
  907   format(1x,1p2e11.3)
        end

c*********************************************************************

        subroutine cascade(nlortz,wlortz,efc,rmass,tref,rh
     &   ,wav0,wdel,xu,wu,nu,z,tmp,prs,clm,den,ep1,ep2,ntau
     &   ,wav,str,eng,wid,nl,nlines,ht,th,ha)

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

c       ----------------------------------------------------
c                         Input Arrays
c       ----------------------------------------------------
        dimension z(1),tmp(1),prs(1),clm(1),den(1),ep1(1),ep2(1)
     &  ,xu(1),wu(1),nl(1),nlines(1),wav(1),str(1),eng(1),wid(1)

c       ----------------------------------------------------
c                         Output Arrays
c       ----------------------------------------------------
        dimension ht(nd1,3), th(nd1,3), ha(nd1,3)

c       ----------------------------------------------------
c                       Internal Arrays
c       ----------------------------------------------------
        dimension h(nd1),temp(nd1),pres(nd1),col(nd1),a(nd1,nd6)
     &  , x1(nd1),x2(nd1)

        data tmp0 / 296.d0 /
        data tmax / 7.d0 /

        do n = 1, ntau-1
           col(n) = dsqrt(clm(n)*clm(n+1))
           temp(n) = half*(tmp(n)+tmp(n+1))
           pres(n) = dsqrt(prs(n)*prs(n+1))
        end do
        col(ntau) = col(ntau-1)**2/col(ntau-2)
        temp(ntau) = temp(ntau-1)
        pres(ntau) = pres(ntau-1)**2/pres(ntau-2)

c       -------------------------------------------------------
c                      Loop through near IR bands
c       -------------------------------------------------------
        do jb = 1, 3

        do n = 1, ntau
           ha(n,jb) = zero
           ht(n,jb) = zero
           th(n,jb) = zero
        end do

        do n = 1, ntau
           x2(n) = ep2(n)/(one+ep2(n))
           x1(n) = ep1(n)/(one+ep1(n))
        end do

        wend = wav(nl(jb)+nlines(jb)-1)
        wlast = wav(nl(jb))
        do while (wlast .lt. wend)

           wnext = dmin1(wlast+wdel,wend)
           wmid = half*(wlast+wnext)
           delw = wnext - wlast
           solflux = solirflx(wmid)/rh**2

           lw1 = min(locate(wav(nl(jb)),nlines(jb),wlast)+1,nlines(jb))
           lw2 = max(locate(wav(nl(jb)),nlines(jb),wnext),1)
           nlw = lw2 - lw1 + 1

           do i = 1, nu
              call random(xu(i),pres,temp,col,ntau,str(lw1)
     &         ,wid(lw1),eng(lw1),nlw,efc,rmass,wmid,delw,a(1,i))
           end do
           
           do n = 1, ntau
              sum = zero
              do i = 1, nu
                if(n .gt. 1) then
                   sum=sum+xu(i)*wu(i)*(a(n,i)-a(n-1,i))
     &              /(col(n)-col(n-1))
                else
                   sum=sum+xu(i)*wu(i)*a(1,i)/col(1)
                endif
              end do
              h(n) = half*solflux*delw*sum
           end do
c           write(unit=*,fmt=*) ' IR heating: ',wmid,h(1)*den(1)

           do n = 1, ntau
              ha(n,jb) = ha(n,jb) + h(n)*den(n)
           end do

c          -----------------------------------------
c              Cascade heating rate (ergs/cm3/s)
c          -----------------------------------------
           fb = float(jb)
           dwr = (wmid-(fb+one)*wav0)/wmid
           rw = wav0/wmid
           fbr = fb*rw
           do n = 1, ntau
              ht(n,jb)=ht(n,jb)+x2(n)*(dwr+fbr*x1(n))*h(n)*den(n)
           end do

c          -----------------------------------------
c              nu-4 source function 
c          -----------------------------------------
           do n = 1, ntau
              th(n,jb) = th(n,jb) + rw*x2(n)*h(n)
           end do

           wlast = wnext

	end do

        end do


        return
        end


        function solirflx(ww)
        parameter (one = 1.d0, c1 = 8.37d-10, c2 = 2.4945d-4)
        real*8 ww,solirflx
        solirflx = c1*(ww**3)/(dexp(c2*ww) - one)
        return
        end