;;;;;;;;;;;;;;;;;;;;;;;;;;;;  EXAMPLE  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; 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.pro                              ;;;
;;;              IDL> EXAMPLE                                     ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
PRO EXAMPLE

  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
  
  ;grid for gray thermal optical depth
  taumin       = 1.e-4 ;minimum value of optical depth, tau
  taumax       = tau0  ;optical depth at reference pressure p0
  Nlvls        = 100   ;number of points
  dlntau       = (ALOG(taumax)-ALOG(taumin))/(Nlvls-1) ;divisions in log space of tau
  tau          = FLTARR(Nlvls) ;array of Nlvls points for logarithm to base 10 of optical depth
  tau[*]       = 0.            ;initialize array of log optical depth to zero
  tau[Nlvls-1] = ALOG(taumax)  ;assign the maximum value of log optical depth
  FOR i=Nlvls-2, 0, -1 DO tau[i] = tau[i+1] - dlntau ;generate the array
  tau          = EXP(tau) ;convert tau from log values to actual values
  
  ;corresponding pressure grid
  p = p0*(tau/tau0)^(1/n)
  
  ;temperature profile, initialized to zero
  T    = p  ;sets size and type of the temperature array equivalent to pressure array
  T[*] = 0. ;initialize all temperatures to zero 
  
  ;temperature in radiative region (RC12 Eq. 18), which 
  ;only applies for optical depths smaller than taurc (i.e., 
  ;above the R-C boundary)
  ;
  sigma   = 5.67e-8 ;Stefan-Boltzmann constant, J/s/m**2/K**4
  D       = 1.66    ;diffusivity factor (see, e.g., Rodgers & Walshaw 1966)
  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)
  iRAD    = WHERE(tau LE taurc) ;find the array subscripts for the radiative regime
  T[iRAD] = (sigmaT[iRad]/sigma)^0.25 ;set temperatures in the radiative regime
  ;
  
  ;temperature in convective region (RC12 Eq. 11), which 
  ;only applies for optical depths larger than taurc (i.e., 
  ;below the R-C boundary)
  ;
  beta    = alpha*(gam - 1.)/gam ;beta defined in Eq. 11 of RC12
  iCON    = WHERE(tau GT taurc)  ;find the array subscripts for the convective regime
  T[iCON] = T0*(tau[iCON]/tau0)^(beta/n) ;set temperatures in the convective regime
  ;
  
  ;plot p-T profile
  PLOT, T, P, XRANGE=[60,175], XSTYLE=1, YRANGE=[1.5,0.001], /YLOG, YSTYLE=1, XTITLE='Temperature [K]', $
              YTITLE='Pressure [bar]', POSITION=[0.12,0.12,0.95,0.95]
END