;;************************************************************
;;This script calculates spec.humidity from air.temp & RH
;;For NCEP-DOE Reanalysis
;;Will work on NCL version 4.0
;;************************************************************
  begin
   year = "1979"
   rh_f = addfile("../rh/rhum."+year+".nc", "r")
   temp_f = addfile("../temp/air."+year+".nc", "r")
   q_f = addfile("../shum/shum_"+year+".nc", "r")
   psfc_f= addfile("../pres/pres.sfc."+year+".nc","r")
   zsfc_f= addfile("hgt.sfc.nc","r")


   rhf = short2flt(rh_f->rhum(:,{1000:300},:,:))
   tempf = short2flt(temp_f->air(:,{1000:300},:,:))
   presf = short2flt(psfc_f->pres(:,:,:))
   qf = q_f->shum(:,{1000:300},:,:)
   zsfcf = short2flt(zsfc_f->hgt(:,:,:))
         
   p  = (/100000,92500,85000,70000,60000,50000,40000,30000/)   
   z = (/111,762,1457,3011,4205,5572,7182,9160/)


   printVarSummary(presf)
   printVarSummary(zsfcf)

   zsfcfconform  = conform(presf,zsfcf,(/0,1,2/))

  
;   printVarSummary(presfconform)
 
   zconform = conform(qf,z,1)
   printVarSummary(zconform)

   pconform = conform(qf,p,1)
    
;   printVarSummary(p)
;   printVarSummary(rh)
;   printVarSummary(temp)
 
;;**************************************************
;;Calculates CAPE(J/kg)
;;**************************************************
 cape = wrf_cape_3d(pconform,tempf,qf,zconform,zsfcfconform,presf,False)
 capeh = cape(0,:,:,:,:)

; printVarSummary(cape)


 
copy_VarCoords(qf,capeh)
;;**************************************************
;;Save as netcdf file
;;**************************************************
  cape@long_name = "CAPE"
  cape@units     = "J/kg"
  diro= "./"                  ; output dir
  filo= "cape_"+year+".nc"  ; output file  

  ptho = diro+filo
  system("/bin/rm -f "+ptho)
  ncdf = addfile(ptho,"c")
  
  fAtt = True
  fAtt@title         = "Convective Available Potential Energy"
  fAtt@source_file   = "NCEP-NCAR Ranalysis-R2"
  fAtt@Conventions   = "None"
  fAtt@creation_date = systemfunc("date")
  fileattdef(ncdf,fAtt)            ; copy file attributes
  filedimdef(ncdf,"time",-1,True)

ncdf->cape = capeh
end