;+ ; NAME: ; rt_szss96_Ch4heat23 ; PURPOSE: (one line) ; calculate CH4 heating at 2.3 micron ; DESCRIPTION: ; calculate CH4 heating at 2.3 micron ; CATEGORY: ; RT ; CALLING SEQUENCE: ; h23 = rt_szss96_Ch4heat23(t, p, n_ch4, gamn4=gamn4) ; INPUTS: ; t - temperature (K) ; p - pressure (microbar) ; n_ch4 - number density of CH4 (cm^-3) ; OPTIONAL INPUT PARAMETERS: ; KEYWORD INPUT PARAMETERS: ; codetype ; 'FORTRAN' - call routines from plutot.f ; 'FORTEST' - native IDL routines that reproduce the fortran ; 'IDL' - idl code, without the restriction of needing to reproduce fortran ; KEYWORD OUTPUT PARAMETERS: ; none ; OUTPUTS: ; returns cooling rate in erg/cm^3/s ; gamn4 : cool-to-space Gamma ; COMMON BLOCKS: ; SIDE EFFECTS: ; RESTRICTIONS: ; PROCEDURE: ; ; EXAMPLE: ; ; atm = rt_szss96_atmread('pluto.restart', 3.d) ; t = atm.t ; p = atm.p ; n_ch4 = atm.n[*,1] ; l77 = rt_szss96_Ch4heat23(t, p, n_ch4) ; MODIFICATION HISTORY: ; 2008 Nov 9 - fixed documentation ;- function rt_szss96_Ch4heat23, t, p, n_ch4, ncol_ch4, gamn4=gamn4, codetype=codetype fcn = 'rt_szss96__Ch4heat23' if not keyword_set(codetype) then codetype = 'FORTRAN' case codetype of 'FORTRAN': begin szss96_get_ch4all,GG,WW,TREF if tref[0] ne 40.d then begin szss96_input6 endif szss96_get_ch4v34, sigx5 if sigx5[0,0] lt 1d-51 then begin szss96_input5 endif szss96_get_parair, R00,WTMOL if WTMOL lt 1. then begin R00 = 1150.0D3 WTMOL = 28.0D szss96_set_parair,R00, WTMOL endif ; convert from standard atm structure to ; fortran's expected format AN = double(n_ch4 * 1d6) BN = double(ncol_ch4 * 1d4) TEM = double(t) PRE = double(p / 10.d) KM = n_elements(PRE) ICONT = 2L IYELLE = 2L IZHU = 1L if not keyword_set(gamn4) then begin BC = 1.3804D-23 DN1 = n_ch4 * 1d6 DN = PRE/(BC * TEM) RHO = double( ( DN1/ DN ) * (16/28.) ) szss96_qtch4z,TEM,PRE,RHO,QT,KM,ICONT,GAMN4,IYELLE endif szss96_htch5,HT,AN,BN,TEM,PRE,KM,ICONT,GAMN4 ;convert from K/s to erg/cm^3/s RAIR=8314.0D0/WTMOL ; gas constant CP=RAIR*7.0D0/2.0D0 ; specific heat at the constant pressure RHO = PRE/(RAIR*TEM) ; kg m**(-3) FACTX1=CP*RHO*10.0D0 ; To convert K/sec to erg cm**(-3) s**(-1) HT2CGS = HT * FACTX1 ht2 = HT2CGS return, ht2 end 'FORTEST': begin szss96_get_parair, R00,WTMOL if WTMOL lt 1. then begin R00 = 1150.0D3 WTMOL = 28.0D szss96_set_parair,R00, WTMOL endif szss96_get_panda, ZKVT,XA2,XA3,XA4,XP2,XP3,XP4 if XA2 lt 1.d then begin szss96_init_panda endif szss96_get_ch4all,GG,WW,TREF if tref[0] ne 40.d then begin szss96_input6 endif szss96_get_ch4v34, sigx5 if sigx5[0,0] lt 1d-51 then begin szss96_input5 endif TEM = double(t) PRE = double(p / 10.d) AN = double(n_ch4 * 1d6) BN = double(ncol_ch4 * 1d4) BC = 1.3804D-23 DN1 = n_ch4 * 1d6 DN = PRE/(BC * TEM) RHO = double( ( DN1/ DN ) * (16/28.) ) KM = n_elements(TEM) MM = 4 LM = 4 JM = 30 HT = dblarr(KM) SIG1 = dblarr(jm, km) gamsl1 = dblarr(km) gamsl2 = dblarr(km) PI=3.14159265D0 XX4 = [0.109063D0, 0.518378D0, 1.052419D0, 1.461733D0] WW4 = [0.273205D0, 0.512194D0, 0.512194D0, 0.273205D0] wk1 = dblarr(km) for ii = 0, mm-1 do begin for jj = 0, mm-1 do begin FFAC1=COS(XX4[ii]) FFAC2=COS(XX4[jj]) XMU1=FFAC1*FFAC2 ;----------------------------- ; CALL HTCH5M(HT,AN,BN,TEM,PRE,KM,XMU1,GAMN4) XMU = XMU1 BK=1.38D-23 ; Boltzmann constant in J K^-1 RCPR=7.0D0/2.0D0 RAIR=8314.0D0/WTMOL ; gas constant CP=RAIR*RCPR ; specific heat at the constant pressure FF0=15.5335D0 ; solar flux [pi B] at [nu3 + nu4] band in [W m^-2 / m^-1] EPS1=2.40232D-8 ; geometric factor for the solar constant in steradian SIGDNU=7.58D5*2.677D-26 ; band strength [sigma * delta nu] in [m] at 40 K XKVT3=XP2*ZKVT ; collisional rate coefficient in [m**3 s**(-1)] XKVT4=XP4*ZKVT ; collisional rate coefficient at T=80 K in [m**3 s**(-1)] log_sigx = alog(SIGX5) for j = 0, jm-1 do begin ; DO 20 J=1,JM sig1[j,*] = $ interpol(log_sigx[j,*], tref, tref[0] > tem < tref[lm-1], /qua) endfor sig1 = exp(sig1) for k = 0, km-1 do begin fac1 = total( ww * sig1[*,k] ) fac2 = total( ww * sig1[*,k] * expint(2, (SIG1[*,K]*BN[K]) < 300.d ) ) fac3 = total( ww * sig1[*,k] * exp( - ((SIG1[*,K]*BN[K]/xmu) < 300.d) ) ) GAMSL1(K)=FAC2/FAC1 ; one-side heating/cooling 0.5D0*FAC2/FAC1 GAMSL2(K)=FAC3/FAC1 endfor RNU34=0.0903D0 ; = [4320 - 3*1310]/4320 AIRN = PRE / (BK*TEM ) ; air number density in m**(-3) PHIK3=XKVT3*AIRN/XA2 PHIK4=XKVT4*AIRN/XA4 PHIKS2=5.0D-2*ZKVT*AN/(2.0D0*XA4) ; testingxz PHIKS3=5.0D-2*ZKVT*AN/(3.0D0*XA4) ; testingxz XRNUX=(1.0D0-RNU34)/3.0D0 ; testingxz FAC11=(1.0D0/(1.0D0+1.0D0/PHIKS3))*$ (1.0D0+2.0D0/(1.0D0+1.0D0/PHIKS2))*XRNUX EPS2=PHIK3/(GAMSL1+PHIK3) ; non-LTE efficient factor FAC22=PHIK4/(GAMN4+PHIK4) EPS3=RNU34+XRNUX*FAC11*FAC22 HT=FF0*EPS1*EPS2*AN*SIGDNU*GAMSL2 ; solar heating rate in [W m^-3] HT=HT*EPS3 ; vibration-vibration cascade HT=HT/(RCPR*PRE/TEM) ; solar heating rate in [K s^-1] ;----------------------------- wk1 = wk1 + ht * FFAC2 * WW4[ii]*WW4[jj] endfor endfor FACT1=(2.0D0/PI)/2.0D0 HT=FACT1*WK1 ; globally averaged solar heating rate in [K s^-1] ;convert from K/s to erg/cm^3/s RAIR=8314.0D0/WTMOL ; gas constant CP=RAIR*7.0D0/2.0D0 ; specific heat at the constant pressure RHO = PRE/(RAIR*TEM) ; kg m**(-3) FACTX1=CP*RHO*10.0D0 ; To convert K/sec to erg cm**(-3) s**(-1) HT2CGS = HT * FACTX1 ht2 = HT2CGS return, ht2 end else: begin print, fcn, ': code type ', codetype, ' not implemented' return, -1 end endcase end