;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; sub.calc_theta_e.ncl
; Carl Schreck (cjschrec@ncsu.edu)
; January 2017
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Description: Based on https://www.ncl.ucar.edu/Applications/Scripts/iso_3.ncl
; and https://www.ncl.ucar.edu/Support/talk_archives/2007/0572.html
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


undef ( "calc_theta_e" )
function calc_theta_e( \
  t : numeric, \ ; temperature
  q : numeric, \ ; specific humidity
  p : numeric, \ ; pressure (must be scalar OR same size as t & q)
  p0 : numeric \ ; must be same units as p (Pa or hPa)
)
begin ; calc_theta_e
  ; These are some parameters that could be useful to have up top
  Lv   = 2.51e6  ; latent heat of vaporization at the triple point [J/kg]; Bluestein, p203        
;;cpd  = 1005.7  ; specific heat at constant pressure for air [Bolton (1980): MWR]
  cpd  = 1004.64 ; specific heat at constant pressure for air 
  R    = 287.04  ; specific gas constant for air [J/(Kg-K)]
  kap  = R/cpd   ; 0.285 

  r = q/(1-q) ; convert spc humidity to mixing ratio

;************************************************
; Calculate Equivalent Potential Temperature 
;************************************************
; Technically, this is a nonlinear computation and it should have
;    been done at each time step and, then, averaged. But ...
;************************************************
; https://en.wikipedia.org/wiki/Equivalent_potential_temperature
;************************************************

                                 ; wikipedia approx
  t_e  = t + ( Lv / cpd ) * r    ;                 Lv     cpd         q
                                 ; Units check: ([J/kg]/[J/(kg-K)])[kg/kg]  => K
  
  retVal = t_e * ( p0 / p )^kap  ; p0 and pm must be same units
 
  retVal@long_name = "equivalent potential temperature"
  retVal@units     = "K"
  copy_VarCoords( t, retVal ) ; assign coordinates
  printVarSummary( retVal )
  printMinMax( retVal, True )

  return(retVal)
end; calc_theta_e