program rates c ********************************************************************** c * * c * * c * Rates calculates the radiative relaxation rate of Titan's * c * atmosphere through eigenvalue analysis of the heating rate * c * jacobian. * c * * c * Roger Yelle * c * 29 MAY 1991 * c * * c ********************************************************************** implicit double precision (o-z, a-h) include "parameters.dat" character fname*64, iret*1 double complex evalue 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), haer(nd1) c --------------------------------------------------------- c Arrays for eigen analysis c --------------------------------------------------------- & , evalue(nd1), evector(nd1,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/ data wdel/10.d0/ data pmin / 1.d-5 / data pmax / 1.d+2 / data rtmax / 0.1d0 / data ifund/1,5,7,8/ 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, 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 tbase = 150.d0 write(unit=*, fmt=901) read (unit=*, fmt=900) fname call atmread1(fname,ntau,z,prs,tmp,den,clm) write(unit=*,fmt= & '(43h Enter CH4, C2H2, C2H6, HCN scale factors >,$)') read (unit=*,fmt=*) fmix(1),fmix(2),fmix(3),fmix(4) write(unit=*,fmt='(27h Enter tolerance and tmin >,$)') read (unit=*,fmt=*) otol,tmin 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 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 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, 4 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) cddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd c imol = 2 c nl(imol) = 1 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 c ************************************************************* c Loop through species calculate heating in fundamentals c ************************************************************* do imol = 1, 3 c imol = 2 do n = 1, ntau ts(n,imol) = 0.d0 end do 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 cooling rate c ********************************************************** if(fmix(4) .gt. zero ) then call radioheat(tmin,wav(nl(4)),str(nl(4)),eng(nl(4)),wid(nl(4)) & ,nlines(4),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) & +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) & + plte(m,n) end do end do do n = 1, ntau cp(n) = 2.5d0*rkb*den(n,0) end do do n = 1, ntau do m = 1, ntau p(m,n) = p(m,n)/cp(m) end do end do c ---------------------------------------------------------- c Calculate eigenvalues/eigenvectors of p matrix c ---------------------------------------------------------- write(unit=*,fmt=*) ' Calculating eignevalues ' call eigen(ntau,p,evalue,evector) 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=920) n,z(n) write(unit=30,fmt=907) z(n),(evector(n,j),j=1,ntau) end do close(unit=30) write(unit=*, fmt=921) read (unit=*, fmt=900) fname open (unit=30, file=fname, status='new') do n = 1, ntau write(unit=30,fmt=922) evalue(n) end do close(unit=30) c write(unit=30,fmt=907) (z(n),n=1,ntau) c j = 0 c 10 continue c j = j + 1 c if(dimag(evalue(j)) .eq. zero) then c write(unit=30,fmt=906) j,evalue(j) c write(unit=30,fmt=907) (evector(n,j),n=1,ntau) c else if(dimag(evalue(j)) .ne. zero) then c write(unit=30,fmt=906) j,evalue(j) c write(unit=30,fmt=907) (evector(n,j),n=1,ntau) c write(unit=30,fmt=907) (evector(n,j+1),n=1,ntau) c write(unit=30,fmt=906) j,evalue(j+1) c write(unit=30,fmt=907) (evector(n,j),n=1,ntau) c write(unit=30,fmt=907) (-evector(n,j+1),n=1,ntau) c j = j + 1 c end if c if(j .lt. ntau) goto 10 c 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 eigenvectors > ',$) 921 format(' Enter filename for eigenvalues > ',$) 906 format('#',i3,1p2e13.5) 920 format('#',i3,1pe13.5) 907 format(1x,1p10e11.3) 922 format(1x,1p2e15.7) end