;+ ; NAME: ; Stability_indices2 ; ; PURPOSE: ; To calculate the stability indices for a thermodynamic profile ; including lifted index, total totals, and showalter index ; CATEGORY: ; ; CALLING SEQUENCE: ; ; INPUTS: ; ; OPTIONAL INPUTS: ; ; KEYWORD PARAMETERS: ; ; OUTPUTS: ; ; OPTIONAL OUTPUTS: ; ; COMMON BLOCKS: ; ; SIDE EFFECTS: ; ; RESTRICTIONS: ; ; PROCEDURE: ; ; EXAMPLE: ; ; MODIFICATION HISTORY: ; Written by: Wayne F. Feltz, CIMSS/SSEC ; ; $Date:$ May 10, 1998 ; $Id:$ ; $Log:$ ; ; ;- pro stability_indices2, nl, p, alt, t, td, mixr, li, totalt, kindex, lcl, showwal,pres_lfc,$ pres_eq, cape, cin li=0 kindex=0 totalt=0 lcl=0 showwal=0 ; Check to make sure T and TD are in Kelvin if t(nl-1) le 150 then begin t=t+273.13 td=td+273.13 endif ;############################################################################# ; Compute Total Totals (t850 + td850 - 2*t500) ;############################################################################# result=abs(p - 500.0) t500_nl=where(min(result) eq result) result=abs(p - 700.0) t700_nl=where(min(result) eq result) result=abs(p - 850.0) t850_nl=where(min(result) eq result) totalt = (t( t850_nl )-273.16) + (td( t850_nl )-273.16) - $ 2.0 * (t( t500_nl )-273.16) ; END Total Totals (t850 + td850 - 2*t500) ;############################################################################# ;############################################################################# ; Compute k-index (t850 + td850) - t500 - (t700 - td700) ;############################################################################# result=abs(p - 700.0) t700_nl=where(min(result) eq result) kindex = (t(t850_nl)+td(t850_nl)) - t(t500_nl)-(t(t700_nl)-td(t700_nl))-273.15 ; End k-index (t850 + td850) - t500 - (t700 - td700) ;############################################################################# ;############################################################################# ; Compute Showwalter Index ;############################################################################# thetae_850=parcel_thetae(850.0,t(t850_nl),td(t850_nl),p(nl-1)) showwal=t(t500_nl)-temp_along_satadiabat(thetae_850, 500.0,p(nl-1)) ;############################################################################# ; End Showwalter Index ;############################################################################# ;############################################################################# ; Compute Lifted Index ;############################################################################# ; Find Parcel Surface Values parcel_boundary=p(nl-1)-100.0 p_parcel=p(nl-1)-50.0 t_sum=0.0 td_sum=0.0 weight_sum=0.0 p_flip=reverse(p) t_flip=reverse(t) td_flip=reverse(td) alt_flip=reverse(alt) mr_flip=reverse(mixr) for i=1,nl-1 do $ if p_flip(i) gt parcel_boundary then begin p_weight=(p_flip(i-1)-p_flip(i)) t_sum=t_sum+(t_flip(i-1)+t_flip(i))*p_weight/2.0 td_sum=td_sum+(td_flip(i-1)+td_flip(i))*p_weight/2.0 weight_sum=weight_sum+p_weight level_top=i+1 endif ;Compute Temperature and Dewpoint Temperature at parcel top t_top=t_flip(level_top)+(t_flip(level_top-1)-t_flip(level_top))*$ (parcel_boundary-p_flip(level_top))/(p_flip(level_top-1)-P_flip(level_top)) td_top=td_flip(level_top)+(td_flip(level_top-1)-td_flip(level_top))*$ (parcel_boundary-p_flip(level_top))/(p_flip(level_top-1)-P_flip(level_top)) ; Average these values in with other parcel attributes p_weight=(p_flip(level_top-1)-parcel_boundary) t_sum=t_sum+(t_flip(level_top-1)+t_top)*p_weight/2.0 td_sum=td_sum+(td_flip(level_top-1)+td_top)*p_weight/2.0 weight_sum=weight_sum+p_weight ; Now define parcel attributes including lifted condensation level t_parcel=t_sum/weight_sum td_parcel=td_sum/weight_sum lcl=pressure_at_lcl(p_parcel, t_parcel, td_parcel, p(nl-1)) ; Calculate the parcel theta E thetae_parcel=parcel_thetae(p_parcel,t_parcel,td_parcel,p(nl-1)) ; Calculate Lifted Index li = t(t500_nl)-temp_along_satadiabat(thetae_parcel, 500.0,p(nl-1)) ; End Lifted Index ;############################################################################# ;############################################################################# ; Compute Parcel Path ;############################################################################# ; First compute level of free convection thetae_sounding=fltarr(nl) tempdiff=fltarr(nl) index_pres_lfc=0 index_pres_eq=0 maxthetae_level=0 for i=0, nl-18 do begin thetae_sounding(i)=parcel_thetae(p_flip(i),t_flip(i),td_flip(i),p(nl-1)) endfor result=where(p_flip GT 500.) maxthetae_level=where(max(thetae_sounding(result)) EQ thetae_sounding(result)) maxthetae_level=fix(fix(maxthetae_level),0) maxthetae_temp=thetae_sounding(maxthetae_level) for i=0, nl-16 do begin tempdiff(i)=temp_along_satadiabat(maxthetae_temp,p_flip(i),p(nl-1))-t_flip(i) endfor for i=0, nl-16 do begin if (tempdiff(i) GE 0.8) AND (index_pres_lfc EQ 0) AND ( i GT maxthetae_level)$ then index_pres_lfc=i if (tempdiff(i) LT 0.8) AND (index_pres_lfc NE 0) AND (index_pres_eq EQ 0) $ then index_pres_eq=i endfor if (index_pres_lfc gt 0) AND (index_pres_lfc lt 32) then begin pres_lfc=p_flip(index_pres_lfc-1)+((p_flip(index_pres_lfc)-p_flip(index_pres_lfc-1))* $ (0.-tempdiff(index_pres_lfc-1))/(tempdiff(index_pres_lfc)-tempdiff(index_pres_lfc-1))) endif else begin pres_lfc=-999 endelse if index_pres_eq ne 0 then begin pres_eq=p_flip(index_pres_eq-1)+((p_flip(index_pres_eq)-p_flip(index_pres_eq-1))* $ (0.-tempdiff(index_pres_eq-1))/(tempdiff(index_pres_eq)-tempdiff(index_pres_eq-1))) endif else begin pres_eq=-999 endelse ;############################################################################# ; Compute CAPE ;############################################################################# cape=0.0 cin=0.0 partcape=fltarr(32) partcin=fltarr(32) if index_pres_lfc gt 0 then begin for i=index_pres_lfc, index_pres_eq-1 do begin alttop=alt_flip(i+1) altbot=alt_flip(i) altdiff=alttop-altbot tparceltop=temp_along_satadiabat(maxthetae_temp,p_flip(i+1),p(nl-1)) tparcelbot=temp_along_satadiabat(maxthetae_temp,p_flip(i),p(nl-1)) mrparceltop=wvmr(p_flip(i+1),tparceltop,tparceltop) mrparcelbot=wvmr(p_flip(i),tparcelbot,tparcelbot) vtemptop=virtemp(tparceltop,mrparceltop) vtempbot=virtemp(tparcelbot,mrparcelbot) parcel_tv_average=0.5*(vtemptop+vtempbot) vtempenvtop=virtemp(t_flip(i+1), mr_flip(i+1)) vtempenvbot=virtemp(t_flip(i), mr_flip(i)) env_tv_average=0.5*(vtempenvtop+vtempenvbot) cape=cape+((9.8/env_tv_average)*(parcel_tv_average-env_tv_average)*altdiff) partcape(i)=cape endfor for i=1, index_pres_lfc-1 do begin alttop=alt_flip(i+1) altbot=alt_flip(i) altdiff=alttop-altbot tparceltop=temp_along_satadiabat(thetae_sounding(i+1),p_flip(i+1),p(nl-1)) tparcelbot=temp_along_satadiabat(thetae_sounding(i),p_flip(i),p(nl-1)) mrparceltop=wvmr(p_flip(i+1),tparceltop,tparceltop) mrparcelbot=wvmr(p_flip(i),tparcelbot,tparcelbot) vtemptop=virtemp(tparceltop,mrparceltop) vtempbot=virtemp(tparcelbot,mrparcelbot) parcel_tv_average=0.5*(vtemptop+vtempbot) vtempenvtop=virtemp(t_flip(i+1), mr_flip(i+1)) vtempenvbot=virtemp(t_flip(i), mr_flip(i)) env_tv_average=0.5*(vtempenvtop+vtempenvbot) cin=cin+((9.8/env_tv_average)*(parcel_tv_average-env_tv_average)*altdiff) partcin(i)=cin endfor cape=max(partcape) cin=min(partcin) endif else begin cape=-9999 cin=-9999 endelse print, cape, cin,li ;############################################################################# end