Diffusivity angle script

This is a function used to compute diffusivity angle of infrared spectra for given atmospheric conditions.

Reference: Feng J, Huang Y. Diffusivity factor approximation for spectral outgoing longwave radiation, Journal of Atmospheric Sciences (submitted)

function [diff,alpha]=dfa(opts,Ht,H,gamma,wavenum,ts,mode_geo);

%%%% opts: total optical depth, nadir view

%%%% Ht: scale height (km) of optical depth

%%%% gamma: temperature lapse rate (K/km); poistively defined

%%%% ts: surface temperature (K)

%%%% mode_geo: 0 if plane-parral atmosphere; 1 if curved earth surface

R       =6371.23; %%% Earth radius (km)

hkB     =1.438833;%%% h/kB

%H      =100;     %%% Height of TOA (km)

r       =0.5772156;%% Euler’s constant

alpha   =(hkB*wavenum/ts^2)*gamma*Ht*(exp(hkB*wavenum/ts))/(exp(hkB*wavenum/ts)–1);

%alpha  =(hkB*wavenum/ts^2/2+1/ts)*gamma*Ht;

if opts<1.45

if mode_geo==0

        scale=1;

        tmpmin=acos(exp(–0.5))/pi*180;

elseif mode_geo==1

        scale=(1–H/R)^2;

else

        display(‘mode_geo must be 0 or 1!’)

end

        a=(1–2*alpha*opts–0.5*alpha*opts^2*(r+log(opts)–2)+alpha/12*opts^3)*scale–(1–alpha/3*opts^3);

        b=–0.5*alpha*opts^3+alpha*opts;

        c=–alpha/4*opts^2+alpha/4*opts^3;

        diff=acos((–b-sqrt(b^2–4*a*c))/2/a)/pi*180;

        if mode_geo==1

        tmp=(–b-sqrt(b^2–4*a*c))/2/a;

        diff=acos(sqrt(1–(1–tmp^2)*scale))/pi*180;

        end

end

if opts>=1.45

if mode_geo==0

        diff=acos(exp(–0.5))/pi*180;

elseif mode_geo==1

       scale=(1–H/R)^2;

       diff=acos(sqrt(1–(1-exp(–0.5)^2)*scale))/pi*180;

end

end