program mesopause

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, iret*1

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

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

c       ---------------------------------------------------------
c            Arrays : line listing
c       ---------------------------------------------------------
     &  , wav(nd4),str(nd4),wid(nd4),eng(nd4),nlines(5),nl(5)
     &  , stot(5),opp0(5),wdopref(5),wmean(5),wstart(5),wend(5)
     &  , wava(nd4),stra(nd4),wida(nd4),enga(nd4),nlinesa(3),nla(3)
     &  , wstep(5)

c       ---------------------------------------------------------
c          Arrays: absorption coefficients
c       ---------------------------------------------------------
     &  , opf(nd1,nd2),wpf(nd2), bf1(nd2), bf2(nd2), bpf(nd1,nd2)
     &  , dbpf(nd1,nd2), dopf(nd1,nd2), oppsum(nd1,4)

c       ---------------------------------------------------------
c            Arrays : Solar UV Flux
c       ---------------------------------------------------------
     &  , soluvwave(nd5), soluvflux(nd5), crs(nd5,0:3)
     
c       ---------------------------------------------------------
c            Arrays for R-T calculation
c       ---------------------------------------------------------
     &  , ep1(nd1),ep2(nd1), bj(nd1,4),bpmean(nd1,4),dep(nd1)
     &  , tmp(nd1), prs(nd1), den(nd1,0:4), clm(nd1,0:4), dtdc(3)
     &  , s(nd1,0:4),ts(nd1,4), xu(nd6), wu(nd6), z(nd1), ep(nd1)

c       ---------------------------------------------------------
c           Arrays for heating calculation
c       ---------------------------------------------------------
     &  , h(nd1), huv(nd1), hcn(nd1), habs(nd1), hcsc(nd1)
     &  , hir(nd1,3), p(nd1,nd1), pcn(nd1,nd1), pir(nd1,nd1,4)
     &  , dtmp(nd1), indx(nd1), tmplst(nd1), cp(nd1), hlte(nd1)
     &  , plte(nd1,nd1), haer(nd1)

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 rh / 9.5d0 /
        data nu /3/
        data nlortz/300/
        data wlortz/5.2d-2,1.d-1,1.d-1,1.d-1/
        data wdel/10.d0/
        data pmin / 1.d-5 /
        data pmax / 1.d+2 /
        data rtmax / 0.1d0 /
        data ifund/1,5,7,8,9/

c       ----------------------------------------------------------
c                     Zero some arrays 
c       ----------------------------------------------------------
        do n = 1, nd1
           h(n) = zero
           hcn(n) = zero
           habs(n) = zero
           hcsc(n) = zero
           do im = 1, 4
              hir(n,im) = zero
              ts(n,im) = zero
           end do
           huv(n) = zero
           hlte(n) = zero
           do m = 1, nd1
              p(m,n) = zero
              pcn(m,n) = zero
              plte(m,n) = zero
              do im = 1, 4
                 pir(m,n,im) = zero
              end do
           end do
        end do

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

        write(unit=*,fmt='(12h Restart ? >,$)')
        read (unit=*,fmt='(a)') iret
        if((iret .eq. 'y').or.(iret .eq. 'Y')) then

           write(unit=*, fmt=901) 
           read (unit=*, fmt=900) fname

           call atmread1(fname,ntau,z,prs,tmp,den,clm)

        else 

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

           write(unit=*, fmt='(19h Enter file name > ,$)') 
           read (unit=*, fmt=900) fname
           call atmread2(fname,pmax,natm,zalt,pres,temp,dens,col)

           call cgrid(gm,rplnt,pmin,pmax,ntau,zalt,temp,pres,dens
     &      ,col,natm,z,tmp,prs,den,clm)

        end if
        zbot = z(ntau)
        write(unit=*,fmt='(15h Enter tbase > ,$)')
        read (unit=*,fmt=*) tbase

        write(unit=*,fmt=
     &   '(43h Enter CH4, C2H2, C2H6, HCN scale factors >,$)') 
        read (unit=*,fmt=*) fmix(1),fmix(2),fmix(3),fmix(4) 

        do im = 1, 4
           do n = 1, ntau
              den(n,im) = fmix(im)*den(n,im)
              clm(n,im) = fmix(im)*clm(n,im)
           end do
        end do

        write(unit=*,fmt='(37h Enter aerosol heating rate at pmax >,$)')
        read (unit=*,fmt=*) haer0
        do n = 1, ntau
           haer(n) = haer0*prs(n)/pmax
        end do

        write(unit=*,fmt='(27h Enter tolerance and tmin >,$)')
        read (unit=*,fmt=*) otol,tmin

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,4 
           rdmol(im) = 1.d-8*half*(rdn2+rdmol(im))
           rmcol(im) = amu*rmn2*rmass(im)/(rmn2+rmass(im))
        end do

c       ----------------------------------------------------------
c                      Read in solar uv heating rates
c       ----------------------------------------------------------
        call solread(soluvwave,soluvflux,crs,nsol)

c       ----------------------------------------------------------
c             Read in line listing for fundamentals
c       ----------------------------------------------------------
        write(unit=*,fmt=*) ' reading line listings'
        n = 1
        do imol = 1, 5
           nl(imol) = n
           jb = ifund(imol)
           call sread(jb,wav(n),str(n),wid(n),eng(n),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)
           n = n + nlines(imol)
           write(unit=*,fmt=810) imol,nl(imol),nlines(imol)
        end do
 810    format(' im = ',i3,' nl = ',i5,' nlines = ',i5)

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)
           write(unit=*,fmt=810) ib,nla(ib),nlinesa(ib)
        end do

        write(unit=*,fmt=
     &   '(46h Enter tolerance, iterations, and time step > ,$)')
        read (unit=*,fmt=*) htol, ntmax, tstep
        tinv = one/tstep

        write(unit=*,fmt=
     &   '(50h Enter tolerance and iterations for Newton loop > ,$)')
        read (unit=*,fmt=*) ttol, nnmax

c       ##########################################################
c                    Begin time step loop
c       ##########################################################

        do n = 1, ntau
           tmplst(n) = tmp(n)
        end do

        ntit = 0
 10     continue

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

c       **********************************************************
c                 Calculate Methane cascade heating 
c       **********************************************************
        imol = 1
        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(imol),qrot(imol),amu*rmass(imol)
     &   ,tref,rh,1306.d0,wdel,xu,wu,nu,z,tmp,prs,clm(1,imol)
     &   ,den(1,imol),ep1,ep2,ntau,wava,stra,enga,wida,nla,nlinesa
     &   ,hcsc,ts(1,imol),habs)

        write(unit=*,fmt=*) ' cascade done'

        do n = 1, ntau
           ts(n,imol) = ts(n,imol)/stot(imol)/four/pi
        end do

c       ##########################################################
c                    Begin Newton-Rhapson loop
c       ##########################################################

        nnit = 0
 20     continue

c       *************************************************************
c           Loop through species calculate heating in vib bands 
c       *************************************************************
        do imol = 1, 4

           do jt = 1, ntau
              ep(jt) = epsln(tmp(jt),den(jt,0),arad(1,imol)
     &        ,pcol(1,imol),rdmol(imol),rmcol(imol))
              dep(jt) = ep(jt)*depsln(tmp(jt),den(jt,0),arad(1,imol)
     &        ,pcol(1,imol),rdmol(imol),rmcol(imol))/tmp(jt)
           end do

           if(fmix(imol) .gt. zero) then

           call absorbk(imol,nlortz,wlortz(imol),qrot(imol)
     &     ,wav(nl(imol)),str(nl(imol)),eng(nl(imol)),wid(nl(imol))
     &     ,nlines(imol),otol,amu*rmass(imol),wstart(imol)
     &     ,wend(imol),wdel,wstep(imol),wdopref(imol),opp0(imol)
     &     ,tbase,tmp,prs,clm(1,imol),ntau,wpf,opf,dopf,nopf,bpf
     &     ,dbpf,bf1,bf2)

           call abscheck(wpf,opf,dopf,bpf,dbpf,bf1,bf2,nopf,ntau)   

           call irheat(tmin,stot(imol),opp0(imol),wpf,opf,dopf,bpf
     &      ,dbpf,bf1,bf2,nopf,tbase,dtdc(imol),xu,wu,nu,tmp,prs
     &      ,den(1,imol),clm(1,imol),ep,dep,ts(1,imol),ntau
     &      ,oppsum(1,imol),bpmean(1,imol),bj(1,imol),s(1,imol)
     &      ,hir(1,imol),pir(1,1,imol))

c           write(unit=*,fmt='(1x,e10.3,2x,f6.1,5e11.3)') 
c     &      (prs(n),tmp(n),oppsum(n,imol),bj(n,imol)/bpmean(n,imol)
c     &      ,s(n,imol),hir(n,imol),pir(n,n,imol),n=1,ntau)

           end if

        end do

c       **********************************************************
c                 Calculate HCN rotational cooling rate
c       **********************************************************
        if(fmix(4) .gt. zero ) then
        call radioheat(tmin,wav(nl(5)),str(nl(5)),eng(nl(5)),wid(nl(5))
     &   ,nlines(5),tbase,tmp,den(1,4),clm(1,4),prs,ntau,xu,wu,nu
     &   ,hlte,plte)
        end if

c       **********************************************************
c                    Calculate conduction heating rate
c       **********************************************************
        call conduct(rplnt,tbase,0.d0,z,tmp,ntau,hcn,pcn)

c       ----------------------------------------------------------
c                         Add them all up
c       ----------------------------------------------------------
        do n = 1, ntau
           h(n)=huv(n)+hcn(n)+hcsc(n)+hir(n,1)+hir(n,2)+hir(n,3)
     &         +hir(n,4)+hlte(n)+haer(n)
           do m = 1, ntau
              p(m,n) = pcn(m,n)+pir(m,n,1)+pir(m,n,2)+pir(m,n,3)
     &               + pir(m,n,4)+plte(m,n)
           end do
        end do

c       ----------------------------------------------------------
c                Calculate temperature correction
c       ----------------------------------------------------------
        do n = 1, ntau
           cp(n) = 2.5d0*rkb*den(n,0)
        end do
        dhmax = zero
        do n = 1, ntau
           dtmp(n) = h(n) - tinv*cp(n)*(tmp(n)-tmplst(n))
           p(n,n) = p(n,n) - cp(n)*tinv 
           dhmax = dmax1(dabs(dtmp(n)),dhmax)
        end do
        call ludcmp(p,ntau,nd1,indx,dinx,ierr)
        if(ierr .eq. 1) then
           write(unit=*,fmt=*) ' Error in mesopause '
           stop 100
        end if
        call lubksb(p,ntau,nd1,indx,dtmp)

        do n = 1, ntau
           if(dtmp(n) .gt. zero) then
              dtmp(n) = dmin1(rtmax*tmp(n),dtmp(n))
           else
              dtmp(n) = dmax1(-rtmax*tmp(n),dtmp(n))
           end if
           tmp(n) = tmp(n) - dtmp(n)
        end do

c       ----------------------------------------------------------
c                         Check convergence
c       ----------------------------------------------------------
        terr = zero
        do n = 1, ntau
           terr=dmax1(dabs(dtmp(n)),terr)
        end do
        nnit = nnit + 1

        write(unit=*,fmt='(8h nnit = ,i5,8h dtmp = ,e11.3
     &   ,9h dhmax = ,e11.3)') nnit,terr,dhmax

        call atmwrite(nnit,terr,ntit,herr,ntau,z,tmp,prs,den,clm)

        if((terr .gt. ttol).and.(nnit .lt. nnmax)) go to 20

        herr = zero
        do n = 1, ntau
           hmax = dmax1(dabs(huv(n)),dabs(hcn(n)),dabs(hcsc(n))
     &      ,dabs(hir(n,1)),dabs(hir(n,2)),dabs(hir(n,3))
     &      ,dabs(hir(n,4)),dabs(hlte(n)))
           herr=dmax1(dabs(h(n))/hmax,herr)
        end do
        ntit = ntit + 1

        do n = 1, ntau
           write(unit=*,fmt='(1x,f6.1,8e10.3)')
     &      tmp(n),h(n),huv(n),hcsc(n),(hir(n,im),im=1,4)
     &      ,hlte(n)
        end do
        write(unit=*,fmt='(8h ntit = ,i5,8h herr = ,e11.3)') ntit,herr
c        read (unit=*,fmt='(a)') icr

c       **********************************************************
c                Solve Hydrostatic equilibrium equation
c       **********************************************************
	call hydrost(clm,tmp,zbot,rplnt,GMm,prs,z,den,tmplst,ntau)

        do n = 1, ntau
           tmplst(n) = tmp(n)
        end do

        if((herr .gt. htol).and.(ntit .lt. ntmax)) go to 10

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),tmp(n),h(n)
     &         ,huv(n),hcn(n),hcsc(n)
     &         ,hir(n,1),hir(n,2),hir(n,3),hir(n,4),hlte(n)
           end do

        close(unit=30)

        write(unit=*, fmt=907) 
        read (unit=*, fmt=900) fname
        open (unit=30, file=fname, status='new') 

           do n = 1, ntau
              write(unit=30, fmt=908) 1.e-5*z(n),prs(n)
     &         ,(s(n,im),im=1,3),ts(n,1)/(one+ep1(n))
           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,1p4e11.3)
  909   format(' Enter filename for time constants > ',$)
  910   format(1x,0pf7.1,1x,0pf7.2,2x,1p5e10.3)

        end



cddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd
c           imol = 2
c           wlortz(imol) = zero 
c           wreach = five*wdopref(imol)
c           wstep(imol) = half*wdopref(imol)
c           wstart(imol) = wmean(imol) - wreach
c           wend(imol) = wmean(imol) + wreach 
c           nlines(imol) = 1
c           wav(nl(imol)) = wmean(imol)
c           opp0(imol) = 1.d-17
c           str(nl(imol)) = opp0(imol)*wdopref(imol)
c           eng(nl(imol)) = 0.d0
c           qrot(imol) = 0.d0
cddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd