program mesonep

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
	logical first 
	real timeused,timeused1,timeused2,timeused3,timeused4, time(2)
	real dtime 

c       ---------------------------------------------------------
c            Arrays : Model atmosphere
c       ---------------------------------------------------------
        dimension zalt(nd3),pres(nd3),temp(nd3),dens(nd3,0:4)
     &   ,col(nd3,0:4), fmix(4), ifund(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),nl(4)
     &  , stot(4),opp0(4),wdopref(4),wmean(4),wstart(4),wend(4)
     &  , wava(nd4),stra(nd4),wida(nd4),enga(nd4),nlinesa(3),nla(3)
     &  , wstep(4)

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

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,3),bpmean(nd1,3),dep(nd1)
     &  , tmp(nd1), prs(nd1), den(nd1,0:4), clm(nd1,0:4), dtdc(3)
     &  , s(nd1,0:3),ts(nd1,3), 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,3)
     &  , dtmp(nd1), indx(nd1), tmplst(nd1), cp(nd1), hlte(nd1)
     &  , plte(nd1,nd1), hsou(nd1), psou(nd1,nd1),hast(nd1) 

c       ---------------------------------------------------------
c          Arrays for time constants
c       ---------------------------------------------------------
     &  , tcn(nd1),tir(nd1,3),tlte(nd1)

c       ---------------------------------------------------------
c                           Data for Neptune 
c       ---------------------------------------------------------
        data rmn2 / 2.d0 /
        data rdn2 / 2.968d0 /
        data GMm / 0.273d-1 /
        data rplnt / 24445.d5 /
        data rh / 30.06d0 /
        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+2 /
        data pmax / 1.d+5 /
        data rtmax / 0.1d0 /
        data ifund/1,5,7,8/

	open(40,file='time_nep.out') 
	write(40,*)'Time serise for the program runing' 
	timeused = dtime(time)
c       ----------------------------------------------------------
c                     Zero some arrays 
c       ----------------------------------------------------------
        do n = 1, nd1
           h(n) = zero
           hcn(n) = zero
           habs(n) = zero
           hcsc(n) = zero
	   hsou(n) = zero 
	   hast(n) = zero
           do im = 1, 3
              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, 3
                 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)

	write(*,*)'before cgrid',temp(51) 
           call cgrid(rplnt,pmin,pmax,ntau,zalt,temp,pres,dens
     &      ,col,natm,z,tmp,prs,den,clm)

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

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

	write(unit=*,911)
	read(*,*)zmax

        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=*) haer

        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-H2 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 solar uvflux
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, 3
           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)
        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
	
	timeused = dtime(time)
	write(40,1000)timeused 

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

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

        ntit = 0
 10     continue

	write(*,*)'ntit=',ntit 
	do n = 1,ntau
	   write(*,*)'tmplst=',tmplst(n),n
	end do
	timeused3 = dtime(time)
	write(40,1003)timeused3 
c       **********************************************************
c                    Calculate UV heating rate
c       **********************************************************
        call uvheat(rh,clm,den,ntau,soluvwave,soluvflux,crs,nsol
     &   ,huv)
	write(*,*)'Did uvheat'
	do n = 1,ntau
	  write(*,*)huv(n)
	end do

c	**********************************************************
c		     Calculate heating from wave < 1.5 um
c	**********************************************************

	call short(rh,ntau,xu,wu,den(1,1),clm(1,1),prs,hast)
	write(*,*)'Did short'
	do n= 1,ntau
	   write(*,*)'hast=',hast(n)
	end do

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),otol
     &   ,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)

c  	open(50,file='casheat.out')
  	write(*,*)'writing cascade heating rates' 
c  	do n = 1,ntau
c   	read (50,1008)hcsc(n),ts(n,imol),habs(n) 
c 	write(50,1008)hcsc(n),ts(n,imol),habs(n) 
c  	write(*,*)n,hcsc(n),ts(n,imol),habs(n)
c 	end do
c 	close(50)
        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
	timeused1 = dtime(time)
	write(40,1001)timeused1 
	write(*,*)'starting NR loop:',nnit 

c       *************************************************************
c           Loop through species calculate heating in fundamentals
c       *************************************************************
         do imol = 1, 3

           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)

	timeused = dtime(time)
	write(40,1005)timeused 


c  	open(60,file='opf_nep.out')
C 	do k = 1,nopf
C 	 write(60,1010)wpf(k),bf1(k),bf2(k)
C 	 do n = 1,ntau
c  	    write(60,1011)opf(n,k),dopf(n,k),bpf(n,k),dbpf(n,k) 
c  	 end do 
c  	end do
c  	close(60)

c  	open(60,file='opf_nep.out')
c 	k = 0
c  100	k = k+1
c  	 read(60,1010,end = 200)wpf(k),bf1(k),bf2(k)
c  	 do n = 1,ntau
c  	    read(60,1011,end=200)opf(n,k),dopf(n,k),bpf(n,k),dbpf(n,k) 
c 	if(k.eq.351 .and. n.eq.1)write(*,*)'in meso,dbpf=',dbpf(n,k) ,dbpf(n,k-2)
c  	 end do 
c  	goto 100
c  200	continue
c 	nopf = k-1

c  	close(60)

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

c  	open(60,file='opfp_nep.out')
c  	do k = 1,nopf
c  	 write(60,1010)wpf(k),bf1(k),bf2(k)
c  	 do n = 1,ntau
c  	    write(60,1011)opf(n,k),dopf(n,k),bpf(n,k),dbpf(n,k) 
c  	 end do 
c  	end do
c  	close(60)

	timeused = dtime(time)
	write(40,1006)timeused

	write(*,*)'I am getting irheat!',ntit,nnit,imol
	
           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))

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

           end if

	timeused = dtime(time)
	write(40,1007)timeused

         end do

	write(*,*)'I did irheat for all three molecules for 
     &   NR loop',nnit,'heating rates loop:',ntit 
c       **********************************************************
c          Calculate H2-h2 and h2-he collisional  cooling rate
c       **********************************************************
	first = .false. 
        if(fmix(4) .gt. zero ) then
	if(ntit.eq.0 .and. nnit.eq.0)first=.true. 
        call radioheat(first,tmin,tbase,tmp,den(1,0),clm(1,0),den(1,4)
     &  ,ntau,xu,wu,nu,hlte,plte)
        end if

c 	do n = 1,ntau
c 	  write(*,*)'hlte=',hlte(n)
c 	end do
	write(*,*)'I did radioheat for time=',nnit ,ntit
c	*********************************************************
c	     Adding in heating source with a Chapmann profile
c	*********************************************************
	if(haer .ne. zero)then
	  call hsource(ntau,z,tmp,rmn2,zmax,haer,hsou,psou) 
	end if 
	write(*,*) 'I did hsource',nnit,ntit 
	write(*,*)(hsou(n),n = 1,ntau)
c       **********************************************************
c                    Calculate conduction heating rate
c       **********************************************************
        call conduct(rplnt,tbase,0.d0,z,tmp,ntau,hcn,pcn)
	write(*,*)'I did conduct for time=',nnit ,ntit
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)
     &         +hlte(n)+hsou(n) + hast(n) 
           do m = 1, ntau
              p(m,n) = pcn(m,n)+pir(m,n,1)+pir(m,n,2)+pir(m,n,3)
     &               + plte(m,n) 
c 	write(*,*)'pmn=',p(m,n),'pcn=',pcn(m,n),'pir1=',pir(m,n,1)
           end do
        end do

c 	write(*,*)'psou=',(psou(n,n),n=1,ntau)
	write(*,*)'heating rates are OK!'
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)
	   write(*,*)'dtmp=',dtmp(n),'dhmax=',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 atm_nep(nnit,terr,ntit,herr,ntau,z,tmp,prs,den,clm)

	timeused2 = dtime(time)
	write(40,1002)timeused2 
        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(hlte(n)),dabs(hsou(n)),dabs(hast(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,3)
     &      ,hlte(n),hcn(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

	timeused4 = dtime(time)
	write(40,1004)timeused4 
        if((herr .gt. htol).and.(ntit .lt. ntmax)) go to 10

        do n = 1, ntau
           tcn(n) = dabs(cp(n)/pcn(n,n))
           tlte(n) = dabs(cp(n)/plte(n,n))
           do im = 1, 3
              tir(n,im) = dabs(cp(n)/(pir(n,n,im)+1.d-30))
           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),tmp(n),h(n)
     &         ,huv(n),hcn(n),hsou(n),hcsc(n)
     &         ,hir(n,1),hir(n,2),hir(n,3),hlte(n),hast(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)

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

           do n = 1, ntau
              write(unit=30, fmt=910) 1.e-5*z(n),tmp(n)
     &        ,tcn(n),(tir(n,im),im=1,3),tlte(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,0pf7.2,1p9e11.2,1pe10.2)
  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)
  911   format(' Enter the altitude for the Max. heating >',$)

 1000	format(/,' Time for reading everything in',
     &    / ' took', f10.3 ,'seconds.'/)
 1001	format(/,' Time used for starting NR loop again',
     &    / ' took', f10.3 ,'seconds.'/)
 1002	format(/,' Time used for finishing NR loop again',
     &    / ' took', f10.3 ,'seconds.'/)
 1003	format(/,' Time used for starting heating  loop again',
     &    / ' took', f10.3 ,'seconds.'/)
 1004	format(/,' Time used for finishing heating loop again',
     &    / ' took', f10.3 ,'seconds.'/)
 1005	format(/,' Time used from starting NR loop to finishing
     &    absorbk',
     &    / ' took', f10.3 ,'seconds.'/)

 1006	format(/,' Time used for doing abscheck once ',
     &    / ' took', f10.3 ,'seconds.'/)
 1007	format(//,' Time used for doing irheat once',
     &    / ' took', f10.3 ,'seconds.'/)
 1008	format(1x,1p3e12.3) 

 1010	format(1x,1p3e12.3)
 1011	format(1x,1p4e12.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