;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; AN_RC_MOD ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; ;;; ;;; Description: Solves for the variables that describe the 1-D thermal structure ;;; ;;; of a planet's atmosphere in radiative-convective equilibrium, ;;; ;;; according to the following article: 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 (hereafter RC12). Papers using this ;;; ;;; code (or a modification of it) should please cite the above paper. ;;; ;;; This solver program is distributed with an example application to ;;; ;;; Titan: EXAMPLE.PRO. These IDL programs can be executed at the command ;;; ;;; line as follows: ;;; ;;; IDL> .r AN_RC_MOD.pro ;;; ;;; IDL> .r EXAMPLE.pro ;;; ;;; IDL> EXAMPLE ;;; ;;; ;;; ;;; A numerical, quasi-Newtonian scheme is used to find the gray thermal ;;; ;;; optical depths of both the radiative-convective boundary and down to a ;;; ;;; reference pressure p0. Absorbed solar (or stellar) flux is divided ;;; ;;; into two channels, and the strength of the atmospheric attenuation in ;;; ;;; these two channels is parametrized in k1 and k2 (see RC12). In ;;; ;;; this solver, channel '1' is taken to be a "stratosphere" channel, ;;; ;;; wherein the value of k1 is set by adjusting the T-p profile to match ;;; ;;; a user-defined top-of-atmosphere (TOA) temperature (Ttoa). Channel '2' ;;; ;;; is taken as a "troposphere" channel, and the value of k2 is set by ;;; ;;; comparing the solar flux absorbed between the TOA and p0 (F20) to the ;;; ;;; total solar flux absorbed in this channel (F2). Other parameters are ;;; ;;; described below, and follow the notation of RC12. ;;; ;;; ;;; ;;; Inputs: ;;; ;;; p0 - reference pressure (e.g., surface pressure) [bar] ;;; ;;; T0 - temperature at reference pressure [K] ;;; ;;; Ttoa - temperature at top of atmosphere (assumed as stratopause) [K] ;;; ;;; F1 - flux absorbed in channel '1' (assumed to be "stratosphere") [W/m**2] ;;; ;;; F2 - flux absorbed in channel '2' (assumed to be "troposphere") [W/m**2] ;;; ;;; F20 - flux absorbed in channel '2' down to reference level p0 [W/m**2] ;;; ;;; Fi - internal heat flux [W/m**2] ;;; ;;; gamma - ratio of specific heats, cp/cv (see Eq. 8 of RC12) ;;; ;;; alpha - empirical adjustment to dry adiabatic lapse rate (see Eq. 10 of RC12) ;;; ;;; n - scaling parameter that relates tau ~ p^n, where tau is gray thermal ;;; ;;; optical depth (see Eq. 6 of RC12) ;;; ;;; LOUD - keyword to print diagnostic quantitates during iterative solving ;;; ;;; process ;;; ;;; Outputs: ;;; ;;; taurc - gray thermal optical depth of radiative-convective boundary ;;; ;;; tau0 - gray thermal optical depth down to reference pressure p0 ;;; ;;; k1 - attenuation parameter for channel '1' [dimensionless] ;;; ;;; k2 - attenuation parameter for channel '2' [dimensionless] ;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO AN_RC_MOD, p0, T0, Ttoa, F1, F2, F20, Fi, gam, alpha, n, taurc, tau0, k1, k2, LOUD=loud ;;; common values passed to numerical solver ;;; COMMON TO_SOLVER,sigma,D,F1C,F20C,F2C,FiC,T0C,p0C,TtoaC,betaC,nC,k1C,k2C,LOUDC ;;; important constants ;;; sigma = 5.67e-8 ;Stefan-Boltzmann constant, J/s/m**2/K**4 D = 1.66 ;diffusivity factor (see, e.g., Rodgers & Walshaw 1966) ;;; set common values using user inputs ;;; F1C = F1 ;flux absorbed in channel '1' (assumed to be "stratosphere") [W/m**2] F20C = F20 ;flux absorbed in channel '2' down to reference level p0 [W/m**2] F2C = F2 ;flux absorbed in channel '2' (assumed to be "troposphere") [W/m**2] FiC = Fi ;internal heat flux [W/m**2] T0C = T0 ;temperature at reference pressure [K] p0C = p0 ;reference pressure (e.g., surface pressure) [bar] TtoaC = Ttoa ;temperature at top of atmosphere [K] nC = n ;scaling power between gray optical depth and pressure betaC = alpha*(gam - 1.)/gam ;beta defined in Eq. 11 of RC12 ;;;; solve non-linear system ;;;; ;X[0] == optical depth of Radiative-Convective (R-C) boundary, taurc ;X[1] == total thermal optical depth down to reference pressure p0 (tau0) ;initial guess that R-C boundary is at D*taurc~1, and tau0~10*p0^n X = [1./D,10.*p0^n] ; initial guesses of taurc and tau0, respectively ;if LOUD, then print the optical depth of the R-C boundary, ;the upwelling thermal flux at this location from the radiative region ;and the convective region, and the temperature from the radiative region ;and the convective region IF KEYWORD_SET( LOUD ) THEN BEGIN LOUDC = LOUD PRINT, ' D*taurc Fup rad Fup conv T rad T conv' PRINT, ' [W/m**2] [W/m**2] [K] [K]' ENDIF ;use a quasi-Newtonian solver to perform non-linear fit RESULT = BROYDEN(X,'BROYFUNC_RCMOD') ;return values for R-C boundary optical depth and total thermal optical depth k1 = k1C ;attenuation parameter for channel '1' (assumed as "stratosphere") k2 = k2C ;attenuation parameter for channel '2' (assumed as "troposphere") taurc= RESULT[0] ; result for gray thermal optical depth of radiative-convective boundary tau0 = RESULT[1] ; result for gray thermal optical depth down to reference pressure p0 END ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; function used in call to IDL routine "BRODYDEN" ;;; ;;; NOTE: X[0] = r-c boundary optical depth ;;; ;;; X[1] = total thermal optical depth down to reference pressure ;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; FUNCTION BROYFUNC_RCMOD, X ;;; common values ;;; COMMON TO_SOLVER,sigma,D,F1C,F20C,F2C,FiC,T0C,p0C,TtoaC,betaC,nC,k1C,k2C,LOUDC F1 = F1C ;flux absorbed in channel '1' (assumed to be "stratosphere") [W/m**2] F20 = F20C ;flux absorbed in channel '2' down to reference level p0 [W/m**2] F2 = F2C ;flux absorbed in channel '2' (assumed to be "troposphere") [W/m**2] Fi = FiC ;internal heat flux [W/m**2] T0 = T0C ;temperature at reference pressure [K] p0 = p0C ;reference pressure (e.g., surface pressure) [bar] Ttoa = TtoaC ;temperature at top of atmosphere [K] beta = betaC ;beta defined in Eq. 11 of RC12 n = nC ;scaling power between gray optical depth and pressure ;;; variables ;;; taurc = X[0] ;current value of optical depth of r-c boundary tau0 = X[1] ;current value of total thermal optical depth down to reference pressure p0 ;;; compute fluxes and temperatures at r-c boundary ;;; ;first check if solver is attempting to use negative values for the ;optical depths (which is unphysical), and, if so, return large numbers ;to the solver to prevent it from trying such values as solutions IF (X[0] LE 0 OR X[1] LE 0) THEN BEGIN RETURN, [1.e3,1.e3] ; ENDIF ELSE BEGIN ;shortwave optical depth in channel "2" is given by ;LOG( F2/(F2-F20) ), and k2 is defined as the ratio of ;this to the thermal optical depth at p0 (i.e., tau0) IF (F2 GT 0) THEN BEGIN k2 = ALOG(F2/(F2-F20))/tau0 ENDIF ELSE BEGIN ;if the flux in channel "2" is zero, ;then no attenuation in this channel k2 = 0 ENDELSE ;k1 is computed by requiring the radiative temperature profile ;to equal the top-of-atmosphere temperature (Ttoa, set by user), ;which is accomplished by evaluating Eq. 18 in RC12 at tau=0 IF (F1 GT 0) THEN BEGIN k1 = D*( 2/F1*( sigma*Ttoa^4 - Fi/2 - F2/2*(1 + k2/D) ) - 1 ) ENDIF ELSE BEGIN ;if the flux in channel "1" is zero, ;then no attenuation in this channel k1 = 0 ENDELSE ;return values of k1, k2 to user k1C = k1 ;attenuation parameter for channel '1' (assumed as "stratosphere") k2C = k2 ;attenuation parameter for channel '2' (assumed as "troposphere") ;four cases: (1) both k1, k2 are greater than zero ; (2) k1 greater than zero, k2 equal to (or less than) zero ; (3) k2 greater than zero, k1 equal to (or less than) zero ; (4) both k1, k2, are equal to (or less than) zero ;;; case 1: both k1, k2 are greater than zero ;;; IF ( k1 GT 0 AND k2 GT 0 ) THEN BEGIN ;upwelling thermal flux at R-C boundary from radiative equilibrium (RC12 Eq. 19) Fup_r = F1/2*(1 + D/k1 + (1 - D/k1)*EXP(DOUBLE(-k1*taurc))) + $ F2/2*(1 + D/k2 + (1 - D/k2)*EXP(DOUBLE(-k2*taurc))) + $ Fi/2*(2 + D*taurc) ;temperature at R-C boundary from radiative equilibrium (RC12 Eq. 18) sigmaT_r = F1/2*(1 + D/k1 + (k1/D - D/k1)*EXP(DOUBLE(-k1*taurc))) + $ F2/2*(1 + D/k2 + (k2/D - D/k2)*EXP(DOUBLE(-k2*taurc))) + $ Fi/2*(1 + D*taurc) T_r = (sigmaT_r/sigma)^0.25 ENDIF ;;; case 2: k1 greater than zero, k2 equal to (or less than) zero ;;; IF ( k1 GT 0 AND k2 LE 0 ) THEN BEGIN k2C = 0. ;upwelling thermal flux at R-C boundary from radiative equilibrium (RC12 Eq. 19) Fup_r = F1/2*(1 + D/k1 + (1 - D/k1)*EXP(DOUBLE(-k1*taurc))) + $ F2/2*(2 + D*taurc) + $ Fi/2*(2 + D*taurc) ;temperature at R-C boundary from radiative equilibrium (RC12 Eq. 18) sigmaT_r = F1/2*(1 + D/k1 + (k1/D - D/k1)*EXP(DOUBLE(-k1*taurc))) + $ F2/2*(1 + D*taurc) + $ Fi/2*(1 + D*taurc) T_r = (sigmaT_r/sigma)^0.25 ENDIF ;;; case 3: k2 greater than zero, k1 equal to (or less than) zero ;;; IF ( k1 LE 0 AND k2 GT 0 ) THEN BEGIN k1C = 0. ;upwelling thermal flux at R-C boundary from radiative equilibrium (RC12 Eq. 19) Fup_r = F1/2*(2 + D*taurc) + $ F2/2*(1 + D/k2 + (1 - D/k2)*EXP(DOUBLE(-k2*taurc))) + $ Fi/2*(2 + D*taurc) ;temperature at R-C boundary from radiative equilibrium (RC12 Eq. 18) sigmaT_r = F1/2*(1 + D*taurc) + $ F2/2*(1 + D/k2 + (k2/D - D/k2)*EXP(DOUBLE(-k2*taurc))) + $ Fi/2*(1 + D*taurc) T_r = (sigmaT_r/sigma)^0.25 ENDIF ;;; case 4: both k1, k2, are equal to (or less than) zero ;;; IF ( k1 LE 0 AND k2 LE 0 ) THEN BEGIN k1C = 0. k2C = 0. ;upwelling thermal flux at R-C boundary from radiative equilibrium (RC12 Eq. 19) Fup_r = F1/2*(2 + D*taurc) + $ F2/2*(2 + D*taurc) + $ Fi/2*(2 + D*taurc) ;temperature at R-C boundary from radiative equilibrium (RC12 Eq. 18) sigmaT_r = F1/2*(1 + D*taurc) + $ F2/2*(1 + D*taurc) + $ Fi/2*(1 + D*taurc) T_r = (sigmaT_r/sigma)^0.25 ENDIF ;upwelling thermal flux at R-C boundary from convective region (RC12 Eq. 13) Fup_c = sigma*T0^4*EXP(D*taurc)*( 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*taurc,/DOUBLE))) ;temperature at R-C boundary from convective region (RC12 Eq. 11) T_c = T0*(taurc/tau0)^(beta/n) IF KEYWORD_SET( LOUDC ) THEN BEGIN PRINT, D*taurc, Fup_r, Fup_c, T_r, T_c ENDIF ;system is solved when the temperatures and upwelling fluxes are continuous ;at the R-C boundary, so the solver attempts to minimize the difference ;between these quantities as computed from the radiative expressions and the ;convective expressions RETURN, [T_r - T_c,Fup_r - Fup_c] ;essentially we solve simultaneous equations: ; T_r - T_c = 0 and Fup_r - Fup_c = 0 ENDELSE END ;;;