;;;;;;;;;;;;;;;;;;;;;;; EXAMPLE WITH FLUXES ;;;;;;;;;;;;;;;;;;;;;;;; ;;; This program is used with AN_RC_MOD.pro to compute a gray ;;; ;;; radiative-convective temperature profile for Titan using ;;; ;;; parameters for Titan, following those used in T.D. Robinson ;;; ;;; and D. C. Catling (2012) "An analytic radiative-convective ;;; ;;; model for planetary atmospheres", Astrophysical Journal, 757, ;;; ;;; 104. doi:10.1088/0004-637X/757/1/104. These IDL programs can ;;; ;;; be executed at the command line as follows: ;;; ;;; IDL> .r AN_RC_MOD.pro ;;; ;;; IDL> .r EXAMPLE_W_FLUXES.pro ;;; ;;; IDL> EXAMPLE_W_FLUXES ;;; ;;; NOTE: this program also uses a user-library function, ;;; ;;; INTEGRAL.PRO to compute downward fluxes, which must be ;;; ;;; accessible in the IDL path. It is available online here: ;;; ;;; http://www.astro.washington.edu/docs/idl/cgi-bin/getpro/library28.html?INTEGRAL ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; function used generate vertical grid with equal ;;; ;;; log-spaced grid points ;;; ;;; zmin = coordinate at top of grid ;;; ;;; zmax = coordinate at bottom of grid ;;; ;;; Nlvls = desired number of levels in grid ;;; ;;; z = output grid ;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FUNCTION GRIDGEN, zmin, zmax, Nlvls dlnz = (ALOG(zmax)-ALOG(zmin))/(Nlvls-1) ;divisions in log space z = FLTARR(Nlvls) ;array of Nlvls points for logarithm to base 10 z[*] = 0. ;initialize grid to zero z[Nlvls-1] = ALOG(zmax) ;assign the maximum value in log spacing FOR i=Nlvls-2, 0, -1 DO z[i] = z[i+1] - dlnz ;generate log spacing values z = EXP(z) ;convert from log(z) to actual values RETURN, z END PRO EXAMPLE_W_FLUXES p0 = 1.5 ;surface pressure [bar] T0 = 94. ;surface temperature [K] Ttoa = 175. ;temperature at stratopause [K] F1 = 1.5 ;flux absorbed in stratosphere [W/m**2] F2 = 1.1 ;flux absorbed in troposphere [W/m**2] F20 = 0.75 ;flux absorbed in troposphere (excluding surface) [W/m**2] Fi = 0. ;internal heat flux [W/m**2] gam = 1.4 ;ratio of specific heats, cp/cv alpha = 0.77 ;adjustment to dry adiabatic lapse rate n = 4./3 ;power for optical depth - pressure scaling (tau~p^n) ;call solver AN_RC_MOD, p0, T0, Ttoa, F1, F2, F20, Fi, gam, alpha, n, taurc, tau0, k1, k2, /LOUD ;print results PRINT, 'Gray Thermal Optical Depth of R-C Boundary: ', taurc PRINT, 'Gray Thermal Optical Depth of Atmosphere: ', tau0 PRINT, 'Attenuation Parameter for Channel 1 (k1): ', k1 PRINT, 'Attenuation Parameter for Channel 2 (k2): ', k2 ;generate p-T profile ;important constants and parameters sigma = 5.67e-8 ;Stefan-Boltzmann constant, J/s/m**2/K**4 D = 1.66 ;diffusivity factor (see, e.g., Rodgers & Walshaw 1966) beta = alpha*(gam - 1.)/gam ;beta defined in Eq. 11 of RC12 ;grid for gray thermal optical depth, separate log-spacing in ;the radiative region and in the convective region taumin = 1.e-4 ;minimum value of optical depth, tau, in radiative region taumax = taurc ;optical depth at bottom of radiative region Nlvlsr = 200 ;number of points in radiative grid taur = GRIDGEN(taumin,taumax,Nlvlsr) taumin = taur[Nlvlsr-2] ;minimum value of optical depth, tau, in convective region taumax = tau0 ;optical depth at reference pressure, p0 Nlvlsc = 25 ;number of points in convective region tauc = GRIDGEN(taumin,taumax,Nlvlsc) tau = [taur[0:Nlvlsr-3],tauc] ;combine grids Nlvls = N_ELEMENTS(tau) ;number of elements in combined grid ;corresponding pressure grid p = p0*(tau/tau0)^(1/n) ;temperature and flux profiles, initialized to zero T = FLTARR(Nlvls) ;temperature Fupth= FLTARR(Nlvls) ;upwelling thermal flux Fdnth= FLTARR(Nlvls) ;downwelling thermal flux Fsol = FLTARR(Nlvls) ;net solar flux Fconv= FLTARR(Nlvls) ;convective flux ;net solar flux in both radiative and convective regions, ;computed from RC12 Eq. 15 Fsol = F1*EXP(-DOUBLE(k1*tau)) + F2*EXP(-DOUBLE(k2*tau)) ;temperature and thermal fluxes in radiative region (RC12 Eq. 18, 19, 20), ;which only apply for optical depths smaller than taurc (i.e., above the ;r-c boundary) ; iRAD = WHERE(tau LE taurc) ;find the array subscripts for the radiative regime ; ;temperature (RC12 Eq. 18) sigmaT = F1/2*(1 + D/k1 + (k1/D - D/k1)*EXP(DOUBLE(-k1*tau))) + $ F2/2*(1 + D/k2 + (k2/D - D/k2)*EXP(DOUBLE(-k2*tau))) + $ Fi/2*(1 + D*tau) T[iRAD] = (sigmaT[iRad]/sigma)^0.25 ;set temperatures in the radiative regime ; ;upwelling thermal flux (RC12 Eq. 19) Fupth[iRAD] = F1/2*(1 + D/k1 + (1 - D/k1)*EXP(DOUBLE(-k1*tau[iRAD]))) + $ F2/2*(1 + D/k2 + (1 - D/k2)*EXP(DOUBLE(-k2*tau[iRAD]))) + $ Fi/2*(2 + D*tau[iRAD]) ; ;downwelling thermal flux (RC12 Eq. 20) Fdnth[iRAD] = F1/2*(1 + D/k1 - (1 + D/k1)*EXP(DOUBLE(-k1*tau[iRAD]))) + $ F2/2*(1 + D/k2 - (1 + D/k2)*EXP(DOUBLE(-k2*tau[iRAD]))) + $ Fi/2*D*tau[iRAD] ; ;temperature and thermal fluxes in convective region (RC12 Eq. 11, 13, 14), ;which only apply for optical depths larger than taurc (i.e., below the ;r-c boundary) ; iCON = WHERE(tau GE taurc) ;find the array subscripts for the convective regime ; ;temperature (RC12 Eq. 11) T[iCON] = T0*(tau[iCON]/tau0)^(beta/n) ;set temperatures in the convective regime ; ;upwelling thermal flux (RC12 Eq. 13) Fupth[iCON] = sigma*T0^4*EXP(D*tau[iCON])*( EXP(DOUBLE(-D*tau0)) + $ 1/(D*tau0)^(4*beta/n)*GAMMA(1+4*beta/n)*(IGAMMA(1+4*beta/n,D*tau0, /DOUBLE) - $ IGAMMA(1+4*beta/n,D*tau[iCON],/DOUBLE))) ;downwelling thermal flux (RC12 Eq. 14), which we break into two parts: (1) the ;integral portion of the solution, which computes the contributions from the atmospheric ;layers below the r-c boundary, and (2) the analytic portion of the solution, which computes ;the attenuated contribution of the thermal flux from the radiative region above ; ;begin by defining a vector of optical depth values in ;the convective region which we will be integrating over taumin = taurc ;minimum value of optical depth, tau taumax = tau0 ;optical depth at reference pressure p0 Nlvls = 200 ;number of points tauc = GRIDGEN(taumin,taumax,Nlvls) ; y = ((tauc/tau0)^(4*beta/n))*EXP(DOUBLE(D*tauc)) ;integrand in RC12 Eq. 14 FOR i=0L, N_ELEMENTS(iCON)-1 DO BEGIN ;compute the integral in RC12 Eq. 14 Fdnth[iCON[i]] = D*sigma*(T0^4)*EXP(DOUBLE(-D*tau[iCON[i]]))*INTEGRAL(tauc,y,XRANGE=[taurc,tau[iCON[i]]]) ENDFOR ;add analytic portion of RC12 Eq. 14 Fdnth[iCON] = Fdnth[iCON] + Fdnth[iCON[0]-1]*EXP(DOUBLE(-D*(tau[iCON]-taurc))) ;convective flux, computed from RC12 Eq. 22 Fconv = Fsol - (Fupth - Fdnth) ;plot p-T profile !P.MULTI=[0,2,2] ;four plots in one figure PLOT, T, p, XRANGE=[60,175], XSTYLE=1, YRANGE=[1.5,0.001], /YLOG, YSTYLE=1, XTITLE='Temperature [K]', $ YTITLE='Pressure [bar]', POSITION=[0.05,0.525,0.475,0.95] PLOT, Fsol, p, XRANGE=[0,3.], XSTYLE=1, YRANGE=[1.5,0.001], /YLOG, YSTYLE=1, XTITLE='Net Solar Flux [W/m!E2!N]', $ YTITLE='Pressure [bar]', POSITION=[0.525,0.525,0.95,0.95] PLOT, Fupth, p, XRANGE=[0,5.], XSTYLE=1, YRANGE=[1.5,0.001], /YLOG, YSTYLE=1, XTITLE='Up/Down Thermal Flux [W/m!E2!N]', $ YTITLE='Pressure [bar]', POSITION=[0.05,0.05,0.475,0.475] OPLOT, Fdnth, p PLOT, Fconv, p, XRANGE=[0,0.1], XSTYLE=1, YRANGE=[1.5,1.4], YSTYLE=1, XTITLE='Convective Flux [W/m!E2!N]', $ YTITLE='Pressure [bar]', POSITION=[0.525,0.05,0.95,0.475] !P.MULTI=[0,1,1] ;reset to one plot per figure END