;+ ; NAME: ; rt_mesopause ; PURPOSE: (one line) ; Reproduce the behavior of mesopause.f ; DESCRIPTION: ; Reproduce the behavior of mesopause.f ; CATEGORY: ; RT ; CALLING SEQUENCE: ; mesopause ; INPUTS: ; OPTIONAL INPUT PARAMETERS: ; none ; KEYWORD INPUT PARAMETERS: ; none ; KEYWORD OUTPUT PARAMETERS: ; none ; OUTPUTS: ; collision rate/deexcitation rate ; COMMON BLOCKS: ; None ; SIDE EFFECTS: ; RESTRICTIONS: ; None ; PROCEDURE: ; ; ; MODIFICATION HISTORY: ; Written 2004 May 8, Leslie Young SwRI ; based on mesopause.f, written by Roger Yelle, 1990 ;- pro rt_mesopause ; set prompt = 1 for interactive prompt = 0 iverbose = 2 ; replaces a call to parameters.h physconstants ; --------------------------------------------------------- ; Data for Titan ; --------------------------------------------------------- dat=rt_datatitan() tref = 150. ; initialize atm, zbot, and tbase iret='n' if prompt then $ read, iret, prompt = ' Restart ? >', format='(a)' if ( (iret eq 'y') or (iret eq 'Y') ) then begin fname = 'titan.atm' if prompt then $ read, fname, prompt = ' Enter model atmosphere file name > ' atm = rt_atmread1(fname) endif else begin nz = 51 if prompt then $ read, nz, prompt = ' Enter ntau > ' fname = 'titan.atm' if prompt then $ read, fname, prompt = ' Enter file name > ' atm = rt_atmread2(fname, dat.pmax) rt_cgrid, dat, nz, atm endelse zbot = atm.z[nz-1] tbase=185. if prompt then $ read, tbase, prompt = 'Enter tbase > ' fmix = [1.d,1.d,1.d,0.5d] if prompt then $ read, fmix, prompt=' Enter CH4, C2H2, C2H6, HCN scale factors >' for im = 1, 4 do atm.n[*,im] = fmix[im-1]*atm.n[*,im] for im = 1, 4 do atm.ncol[*,im] = fmix[im-1]*atm.ncol[*,im] if (iverbose ge 2) then begin print, 'z' print, atm.z[0:2], form='(3e11.3)' print, '..', atm.z[nz-2:nz-1], form='(A11,2e11.3)' print, 'density of N2' print, atm.n[0:2,0], form='(3e11.3)' print, '..', atm.n[nz-2:nz-1,0], form='(A11,2e11.3)' print, 'density of HCN' print, atm.n[0:2,4], form='(3e11.3)' print, '..', atm.n[nz-2:nz-1,4], form='(A11,2e11.3)' print, 'column density of N2' print, atm.ncol[0:2,0], form='(3e11.3)' print, '..', atm.ncol[nz-2:nz-1,0], form='(A11,2e11.3)' print, 'column density of HCN' print, atm.ncol[0:2,4], form='(3e11.3)' print, '..', atm.ncol[nz-2:nz-1,4], form='(A11,2e11.3)' endif ; Aerosol heating -- haer haer0 = 2.5d-5 & hta = 5.d6 if prompt then $ read, haer0,hta, prompt = 'Enter aerosol heating rate at pmax >' haer = haer0 * exp( -(atm.z-zbot)/hta) if (iverbose ge 2) then begin print, 'Calculate aerosol heating rate' print, haer[0:2], form='(3e11.3)' print, '..', haer[nz-2:nz-1], form='(A11,2e11.3)' endif otol=0.2 & tmin = 1.d-10 if prompt then $ read, otol, tmin, prompt = 'Enter tolerance and tmin >' ; ---------------------------------------------------------- ; Set up angular quadratures for R-T calculations ; ---------------------------------------------------------- gauleg, 0.d, 1.d, xu, wu, dat.nu ; ---------------------------------------------------------- ; Read in IR vibration transition parameters and ; Calculate reduced mass and collision diameters ; for hydrocarbon-N2 collisions ; ---------------------------------------------------------- band=rt_bandread() band.rdmol = 1.d-8 * 0.5d * (dat.rdn2 + band.rdmol) band.rmcol = !phys.m_u*dat.rmn2*band.rmass/(dat.rmn2+band.rmass) ; ---------------------------------------------------------- ; Read in solar uvflux ; ---------------------------------------------------------- sol = rt_solread() if (iverbose ge 2) then begin print, 'N2 cross section' print, sol.crs[0:2,0], form='(3e11.3)' print, '..', sol.crs[sol.nsol-2:sol.nsol-1,0], form='(A11,2e11.3)' print, 'C2H6 cross section' print, sol.crs[0:2,3], form='(3e11.3)' print, '..', sol.crs[sol.nsol-2:sol.nsol-1,3], form='(A11,2e11.3)' endif ; ---------------------------------------------------------- ; Read in line listing for fundamentals ; ---------------------------------------------------------- sfund = rt_sread(dat.ifund) wdopref = sqrt(2.*!phys.k*tref/(band.rmass*!phys.m_u)) * $ (sfund.wmean/!phys.c) sfund2 = {wdopref:wdopref, opp0:sfund.stot/wdopref, wstep:wdopref} ; ---------------------------------------------------------- ; Read in line listing for solar bands ; ---------------------------------------------------------- ssolar = rt_sread([1,2,3]) htol = 1.d-10 & ntmax=1 & tstep = 1.d20 if prompt then $ read, htol, ntmax, tstep, $ prompt = ' Enter tolerance, iterations, and time step > ' tinv = 1.d/tstep htol = 1.d-2 & nnmax = 1 if prompt then $ read, ttol, nnmax, $ prompt = ' Enter tolerance and iterations for Newton loop > ' ; ########################################################## ; Begin time step loop ; ########################################################## if (iverbose ge 2) then print, 'Begin time step loop' tlst = atm.t istop=0 ntit = 0 herr = 10.d0*htol while((herr gt htol) and (ntit lt ntmax) and (istop eq 0) ) do begin ; ********************************************************** ; Calculate UV heating rate ; ********************************************************** print, 'Calculate UV heating rate' huv = rt_uvheat(dat.rh, atm, sol) if (iverbose ge 2) then begin print, huv[0:2], form='(3e11.3)' print, '..', huv[nz-2:nz-1], form='(A11,2e11.3)' endif ; ********************************************************** ; Calculate Methane cascade heating ; ********************************************************** if (iverbose ge 2) then print, 'Methane cascade heating' im = 0 ib = 0 ep1 = rt_epsln(atm.t, atm.n[*,0],$ band.arad[ib,im], band.pcol[ib,im],$ band.rdmol[im], band.rmcol[im]) ib=3 ep2 = rt_epsln(atm.t, atm.n[*,0],$ band.arad[ib,im], band.pcol[ib,im],$ band.rdmol[im], band.rmcol[im]) if (iverbose ge 2) then begin print, 'ep1' print, ep1[0:2], form='(3e11.3)' print, '..', ep1[nz-2:nz-1], form='(A11,2e11.3)' print, 'ep2' print, ep2[0:2], form='(3e11.3)' print, '..', ep2[nz-2:nz-1], form='(A11,2e11.3)' ; FORTRAN OUTPUT ; ep1 ; 0.710E-07 0.981E-07 0.136E-06 ; .. 0.524E+00 0.723E+00 ; ep2 ; 0.867E-04 0.120E-03 0.167E-03 ; .. 0.640E+03 0.883E+03 endif ts = dblarr(atm.nz, 4) hcsc = rt_cascade(band.qrot[im],!phys.m_u*band.rmass[im], $ dat.rh,1306.d0,dat.wdel,xu,wu,dat.nu,$ atm,ep1,ep2,ssolar,ts) ts[*,im] = ts[*,im]/(sfund.stot[im]*4.*!pi) if (iverbose ge 2) then begin print, 'hcsc' print, hcsc[0:2], form='(3e11.3)' print, '..', hcsc[nz-2:nz-1], form='(A11,2e11.3)' endif ; FORTRAN OUTPUT ; 0.895E-15 0.146E-14 0.239E-14 ; .. 0.776E-07 0.992E-07 ; ########################################################## ; Begin Newton-Rhapson loop ; ########################################################## if (iverbose ge 2) then print, 'Begin Newton-Rhapson loop' nnit = 0 terr = 10.d0^tol while(terr gt ttol) and (nnit lt nnmax) and (istop eq 0)) do begin ; ************************************************************* ; Loop through species calculate heating in vib bands ; ************************************************************* if (iverbose ge 2) then print, 'calculate heating in vib band' for imol=0,3 do begin ib = 0 ep = rt_epsln(atm.t, atm.n[*,0],$ band.arad[ib,im], band.pcol[ib,im],$ band.rdmol[im], band.rmcol[im]) dep = drt_epsln(atm.t, atm.n[*,0],$ band.arad[ib,im], band.pcol[ib,im],$ band.rdmol[im], band.rmcol[im]) * (ep/atm.t) endfor ; loop over imol for vib heating ; ---------------------------------------------------------- ; make sure we only go through once, in micro mode ; ---------------------------------------------------------- istop=1 ; ---------------------------------------------------------- ; Check convergence ; ---------------------------------------------------------- nnit = nnit + 1 endwhile ; Newton-Rapson loop ntit = ntit+1 endwhile ; time step loop stop end