;;load "./static_stability.ncl"

;************************************************
; variable and file handling
;************************************************
  diri    = "/Users/shea/Data/netCDF/"
  f1      = addfile(diri+"temp_lag4tolead2.nc","r")    ; open netcdf file
  f2      = addfile(diri+"omega_lag4tolead2.nc","r")     
  f3      = addfile(diri+"uwnd_lag4tolead2.nc","r")
  f4      = addfile(diri+"vwnd_lag4tolead2.nc","r")
                                                       ; convenience: make all:  S->N
  temp    = f1->air(:,:,::-1,:)                        ; degK
  omega   = f2->omega(:,:,::-1,:)                      ; Pascal/s
  uwnd    = f3->uwnd(:,:,::-1,:)                       ; m/s
  vwnd    = f4->vwnd(:,:,::-1,:)                       ; m/s

  time    = f1->time                                   ; "hours since 1800-1-1"
  time    = time*3600                                          
  time@units  = "seconds since 1800-1-1 00:00:0.0"
  printVarSummary(time)
  print("---")
  t       = time                 ; Lyndz' name 

  p       = f1->level                                  ; hPa [*]
  p       = p*100                                      ; Pa  [100000,...,10000]
  p@units = "Pa"
  p!0     = "p"
  p&p     =  p                   ; not necessary
  printVarSummary(p)
  print("---")
  
  Cp      = 1.00464e3 ; specific heat of dry air at constant pressure 
  Cp@units= "J/(K*kg)"  
 
;*******************************************
;---Compute local dT/dt
;*******************************************

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

;****************************************
;---Compute advection term: spherical harmonics
;---U*(dX/dlon) + V*(dX/dlat)
;****************************************
  opt_adv   = 0
  long_name = "temp advection"
  units     = "K/s"
  gridType  = 1  
  advT      = advect_variable(uwnd,vwnd,temp,gridType,long_name,units,opt_adv)
  printVarSummary(advT)
  printMinMax(advT,0)
  print("-----")

;********************************************
; 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.
;********************************************

  Q1 = Cp*dTdt - Cp*(omega*s1+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("-----")
   
;***********************************************
; save to a netcdf file
;***********************************************
  diro = "./"
  filo = "heatsource_lag4tolead2.nc"
  ptho = diro+filo
  system("/bin/rm -f "+ptho)
  ncdf = addfile(ptho,"c")

  fAtt = True
  fAtt@title         = "Apparent Heat Source based on Yanai et al. 1973"
  fAtt@source_file   = "NCEP-NCAR Reanaysis"
  fAtt@Conventions   = "None"
  fAtt@creation_date = systemfunc("date")
  fileattdef(ncdf,fAtt)            ; copy file attributes
  filedimdef(ncdf,"time",-1,True)
  ncdf->Q1 = Q1