program atmosphere 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 c --------------------------------------------------------- c Arrays : Model atmosphere c --------------------------------------------------------- dimension zalt(nd3),pres(nd3),temp(nd3),dens(nd3,0:3) & ,col(nd3,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 --------------------------------------------------------- & , wav(nd4),str(nd4),wid(nd4),eng(nd4) 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,0:3), ep2(nd1,0:3), bj(nd1,0:3), bpmean(nd1,0:3) & , tmp(nd1), prs(nd1), den(nd1,0:3), clm(nd1,0:3), dtdc(3) & , s(nd1,0:3), ts(nd1,0:3), xu(nd6), wu(nd6), z(nd1) c --------------------------------------------------------- c Arrays for heating calculation c --------------------------------------------------------- & , h(nd1), huv(nd1), hcn(nd1), habs(nd1), hcsc(nd1) & , hir(nd1,0:3), pcn(nd1,nd1) c --------------------------------------------------------- c Data for Titan c --------------------------------------------------------- data rmn2 / 28.d0 / data rdn2 / 3.d0 / data gm / 8.95e18 / data rplnt / 2575.d5 / data rh / 9.5d0 / data nu /3/ data nlortz/300/ data wlortz/5.2d-8,1.d-7,1.d-7/ data wdel/10.d0/ c ---------------------------------------------------------- c Define optional run parameters c ---------------------------------------------------------- data pmin / 1.d-5 / data pmax / 1.d+3 / 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 c write(unit=*,fmt=904) c read (unit=*,fmt=*) iopt,otol otol = 0.5d0 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*(one/(one/rmass(im)+one/rmn2)) end do c ---------------------------------------------------------- c Read in solar uvflux c ---------------------------------------------------------- call solread(soluvwave,soluvflux,crs,nsol) c ---------------------------------------------------------- c Read in Model Atmosphere c ---------------------------------------------------------- call atmread(fname,pmax,natm,zalt,pres,temp,dens,col & ,tbot,dtdz) 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,dtdz,dtdc) do n = 1, ntau write(unit=*,fmt='(1x,f6.1,1x,e10.3,f6.1,2e11.3)') & 1.d-5*z(n),prs(n),tmp(n),den(n,0),clm(n,0) end do c read(unit=*,fmt='(a)') icr c ********************************************************** c Calculate UV heating rate c ********************************************************** call uvheat(rh,clm,den,ntau,soluvwave,soluvflux,crs,nsol & ,huv) do n = 1, ntau write(unit=*,fmt=801) 1.d-5*z(n), huv(n) end do c read(unit=*,fmt='(a)') icr 801 format(1x,0pf7.1,2x,1pe11.3) c ********************************************************** c Calculate conduction heating rate c ********************************************************** call conduct(rplnt,dtdz,0.d0,z,tmp,ntau,hcn,pcn) do n = 1, ntau write(unit=*,fmt=801) 1.d-5*z(n), hcn(n) end do c read(unit=*,fmt='(a)') icr c ********************************************************** c Calculate Methane cascade heating c ********************************************************** imol = 1 do jt = 1, ntau ep1(jt,imol) = epsln(tmp(jt),den(jt,0),arad(1,imol) & ,pcol(1,imol),rdmol(imol),rmcol(imol)) ep2(jt,imol) = 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),rmcol(imol),tref & ,rh,1306.d0,wdel,tmp,prs,clm(1,imol),den(1,imol),ep1(1,imol) & ,ep2(1,imol),ntau,hcsc,ts(1,imol),habs) do n = 1, ntau write(unit=*,fmt='(1x,f6.1,2x,2e11.3)') & 1.d-5*z(n),habs(n),hcsc(n) end do c ************************************************************* c Loop through species calculate heating in fundamentals c ************************************************************* do imol = 1, 3 if(imol .eq. 1) then jb = 1 else if(imol .eq. 2) then jb = 5 else if(imol .eq. 3) then jb = 7 end if c --------------------------------- c Read in line listing c --------------------------------- call sread(jb,wav,str,wid,eng,nlines,stot,wmean & ,wstart,wend) wdopref = dsqrt(two*rkb*tref/rmcol(imol))*wmean/cs opp0 = stot/wdopref c ---------------------------------------------------------- c Calculate cascade contribution c ---------------------------------------------------------- if(imol .eq. 1) then do n = 1, ntau ts(n,imol) = ts(n,imol)/opp0/four/pi end do else do n = 1, ntau ts(n,imol) = zero end do end if do jt = 1, ntau ep1(jt,imol) = epsln(tmp(jt),den(jt,0),arad(1,imol) & ,pcol(1,imol),rdmol(imol),rmcol(imol)) end do call irheat(otol,nlortz,wlortz(imol),wav & ,str,eng,nlines,wstart,wend,wdel,wave(1,imol),dgn(1,imol) & ,qrot(imol),rmcol(imol),tref,wdopref,opp0,dtdc(imol) & ,xu,wu,nu,tmp,prs,den(1,imol),clm(1,imol),ep1(1,imol) & ,ts(1,imol),ntau,bpmean(1,imol),bj(1,imol),s(1,imol) & ,hir(1,imol)) do n = 1, ntau write(unit=*,fmt='(1x,f6.1,2x,e11.3)') 1.d-5*z(n),hir(n,imol) end do end do 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) 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),habs(n),hcsc(n) & ,hir(n,1),hir(n,2),hir(n,3) end do close(unit=30) write(unit=*, fmt=907) read (unit=*, fmt=900) fname open (unit=30, file=fname, status='new') do n = 1, natm write(unit=30, fmt=908) 1.e-5*zalt(n),prs(n) & ,(ep1(n,im),s(n,im),im=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 iopt and otol > ',$) 905 format(' Enter filename for heating rates > ',$) 906 format(1x,0pf7.1,1x,0pf7.2,2x,1p8e11.3) 907 format(' Enter filename for source functions > ',$) 908 format(1x,0pf7.1,1x,1pe10.3,2x,3(1pe11.3,1x,1pe11.3,2x)) end