program radiodriver 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,hta do n = 1, ntau haer(n) = haer0*dexp(-(z(n)-zbot)/hta) c 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 uvflux c ---------------------------------------------------------- call solread(soluvwave,soluvflux,crs,nsol) c ---------------------------------------------------------- c Read in line listing for fundamentals c ---------------------------------------------------------- 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) 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 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 herr = 10.d0*htol cccccccc do while((herr .gt. htol).and.(ntit .lt. ntmax)) c ********************************************************** c Calculate UV heating rate c ********************************************************** call uvheat(rh,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) 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 terr = 10.d0*ttol cccccccc do while((terr .gt. ttol).and.(nnit .lt. nnmax)) 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) end do 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)),dabs(haer(n))) herr=dmax1(dabs(h(n))/hmax,herr) ccccccc end do ntit = ntit + 1 cd do n = 1, ntau cd write(unit=*,fmt='(1x,f6.1,8e10.3)') cd & tmp(n),h(n),huv(n),hcsc(n),(hir(n,im),im=1,4) cd & ,hlte(n) cd 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 cccccccc end do c ---------------------------------------------------------- c Write output to Disc c ---------------------------------------------------------- write(unit=*, fmt=905) read (unit=*, fmt=900) fname open (unit=30, file=fname, status='unknown') 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) & ,haer(n) end do close(unit=30) write(unit=*, fmt=907) read (unit=*, fmt=900) fname open (unit=30, file=fname, status='unknown') 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=911) read (unit=*, fmt=900) fname open (unit=30, file=fname, status='unknown') write(unit=30, fmt=913) do n = 1, ntau write(unit=30, fmt=912) 1.e-5*z(n),tmp(n),prs(n) & ,(den(n,i),i=0,4) 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,1p10e11.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) 911 format(' Enter filename for model atmosphere > ',$) 912 format(1x,0pf7.1,1x,0pf7.2,2x,1p6e10.3) 913 format(' z T P [N2] [CH4]' & ,' [C2H2] [C2H6] [HCN]') 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