undef("Q1Q2_yanai.ncl")
function Q1Q2_yanai(time[*]:numeric,p,u,v,H,T,omega,npr[1]:integer,opt)
;
; References:
;
; Yanai, M. (1961): 
; A detailed analysis of typhoon formation.
; J. Meteor. Soc. Japan , 39 , 187–214
; 
; Yanai et al (1973):
;  Determination of bulk properties of tropical cloud clusters 
;    from large-scale heat and moisture budgets.
;  J.  Atmos.  Sci.  , 30 , 611–627,
;  https://doi.org/10.1175/1520-0469(1973)030<0611:DOBPOT>2.0.CO;2
;
;********************************************
; Diabatic heating in the atmosphere is a combined consequence of 
; radiative fluxes, phase changes of water substance, and turbulent 
; flux of sensible heat from the earth's surface. 
; 
; In the tropics, it is the major driving force of the atmospheric circulation. 
; It responds to the vertical gradient of diabatic heating.
;********************************************
local dimp, dimu, dimv, dimH, dimT, dimo       \
    , rankp, ranku, rankv, rankH, rankT, ranko \
    , Cp, Lc, dTdt, ss, opt_adv, long_name, units, gridType, advT, Q1 \
    , dHdt, advH, loq, Q2
begin
;;             Use for error testing
;;dimp         = dimsizes(p)
;;dimu         = dimsizes(u)
;;dimv         = dimsizes(v)
;;dimH         = dimsizes(H)
;;dimT         = dimsizes(T)
;;dimo         = dimsizes(omega)

;;rankp        = dimsizes(dimp)
;;ranku        = dimsizes(dimu)
;;rankv        = dimsizes(dimv)
;;rankH        = dimsizes(dimH)
;;rankT        = dimsizes(dimT)
;;ranko        = dimsizes(dimo)

  Cp           = 1.00464e3 
  Cp@long_name = "specific heat of dry air at constant pressure"
  Cp@units     = "J/(K*kg)"  

  Lc           = 2.26e6    ; [J/kg]=[m2/s2]  Latent Heat of Condensation/Vaporization
  Lc@long_name = "Latent Heat of Condensation/Vaporization"
  Lc@units     = "J/kg"  ; ==> "m2/s2"
;*******************************************
;---Compute local dT/dt
;*******************************************

  dTdt = center_finite_diff_n (T,time,False,0,0)
  copy_VarCoords(T, dTdt)
  dTdt@longname = "Temperature: Local Time derivative"
  dTdt@units = "K/s"
  printVarSummary(dTdt)
  printMinMax(dTdt,0)
  print("-----")
   
;****************************************
;---Compute static stability
;****************************************
  ss  = static_stability (p , temp, 1, 0) 
  printVarSummary(ss)
  printMinMax(ss,0)
  print("-----")

;****************************************
;---Compute advection term: spherical harmonics
;---U*(dT/dlon) + V*(dT/dlat)
;****************************************
  opt_adv   = 0
  long_name = "temp advection"
  units     = "K/s"
  gridType  = 1  
  advT      = advect_variable(uwnd,vwnd,T,gridType,long_name,units,opt_adv) 
   
;****************************************
;---Apparent Heat Source
;****************************************

  Q1 = Cp*dTdt - Cp*(omega*ss-advT)
  Q1@long_name = "Total Diabatic Heating as the Apparent Heat Source"
  Q1@units     = "K/s"     
  copy_VarCoords(temp,Q1)
  printVarSummary(Q1) 
  printMinMax(Q1,0)
  print("-----")

;*******************************************
;---Compute local dH/dt
;*******************************************

  dHdt = center_finite_diff_n (H,time,False,0,0)
  copy_VarCoords(H, dHdt)
  dHdt@longname = "Specific Humidity: Local Time derivative"
  dHdt@units = "1/s"
  printVarSummary(dHdt)
  printMinMax(dHdt,0)
  print("-----")

;****************************************
;---Compute advection term: spherical harmonics
;---U*(dH/dlon) + V*(dH/dlat)
;****************************************
  long_name = "moisture advection"
  units     = "1/s"
  advH      = advect_variable(uwnd,vwnd,H,gridType,long_name,units,opt_adv) 

;****************************************
;---Compute vertical movement of H
;---U*(dH/dlon) + V*(dH/dlat)
;****************************************

  dHdp = center_finite_diff_n (H,p,False,1,0)
  copy_VarCoords(H, dHdt)
  dHdp@longname = "Specific Humidity: Local Vertical Derivative"
  dHdp@units = "[kg/kg]/Pa""
  printVarSummary(dHdp)
  printMinMax(dHdp,0)
  print("-----")

  dHdp  = omega*dHdp
  dHdp@longname = "Specific Humidity: Vertical Moisture Transport"
  dHdp@units = "???"
  printVarSummary(dHdp)
  printMinMax(dHdp,0)
  print("-----")
   
;****************************************
;---Apparent Heat Source
;---Turbulent flux of moisture
;****************************************

  dHdt  = -Lc*dHdt
  advH  = -Lc*advH
  loq   = -Lc*dHdp

  Q2 = dHdt + advH + loq
  Q2@long_name = "Apparent Moisture Sink"
  Q2@units     = "C/s"   ; "K/s"
  copy_VarCoords(temp,Q2)
  printVarSummary(Q2) 
  printMinMax(Q2,0)
  print("-----")

  return( [/Q1, Q2/] )
end