;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"


;;Diri  = "/B0/CMIP5/CMIP5/GFDL-ESM2M/3d/"
;;FLIST = systemfunc("ls -1 "+Diri+"ua_6hrLev_GFDL-ESM2M_historical_r1i1p1_1976*")
Diri    = "./"
FLIST   = systemfunc("ls "+Diri+"ua_6hrLev_GFDL-ESM2M_rcp85_r1i1p1.nc")

NFILES  = dimsizes(FLIST)
print(FLIST)

;;Diro  = "/R3/kelly/NCL/"
Diro    = "./"

j=0
do while (j.lt.NFILES)
   Filo = "ua.test_"+sprinti("%0.4i", j)+".p0.nc"
   system("/bin/rm -f "+Diro+Filo)                  ; rm any pre-exist files
   Fout = addfile(Diro+Filo, "c")
   data = addfile(FLIST,"r")

   copy_VarAtts(data, Fout) 

;**********************************************
; read needed variables from file
;************************************************
   ua = data->ua                                    ; select variable to save
   a  = data->a                                     ; get a coefficiants (hyam)
   b  = data->b                                     ; get b coefficiants (hybm)
   ps = data->ps                                    ; get sfc pressure (Pa)
   p0 = data->p0

   p0new = p0/100                                   ; vinth2p requires mb or hPa
   p0new@units = "hPa"
;************************************************
; For illustration, calculate the pressures at all hybrid levels and gridpoints
; Arbitrarily, choose a grid point. Here:  (0, 180)
;************************************************

   printVarSummary(ua)
   print(" ")
   print("min(ua)="+min(ua)+"   max(ua)="+max(ua))
   print(" ")
   print("min(ps)="+min(ps)+"   max(ps)="+max(ps))

   p = pres_hybrid_ccm(ps*0.01,p0new,a,b)              ; bottom-to-top; all hPa
   copy_VarCoords(ps, p(:,0,:,:))                      ; trick
   printVarSummary(p)

   print(" ")
   print("---> a,b,u:  original order is bottom-to-top ------")
   print( sprintf("%8.5f", a)+"  "+sprintf("%8.5f", b)+"  " \
         +sprintf("%8.2f", p(0,:,{0},{180}))          +"  " \
         +sprintf("%8.2f",ua(0,:,{0},{180})) )
   print("------------------------------------------")

;************************************************
; define other arguments required by vinth2p, change units
;************************************************
   interp = 3       ; type of interpolation: 1 = linear, 2 = log, 3 = loglog
   extrap = False   ; is extrapolation desired if data is outside the range of PS

   pnew   = (/ 850.0,700.0,500.0,300.0,200.0 /)       ; create an array of desired pressure levels:
   pnew   = pnew(::-1)
   pnew@long_name = "pressure levels"
   pnew@units     = "hPa"
   pnew!0         = "pnew"

   print(" ")
   print(pnew)

;************************************************
; As per documentation, the function requires data to be 
; ordered from top to bottom. Reverse the order using ::-1 syntax.
;************************************************
   ua = ua(:,::-1,:,:)
   a  = a(::-1)
   b  = b(::-1)

   pp = pres_hybrid_ccm(ps*0.01,p0new,a,b)              ; bottom-to-top
   copy_VarCoords(ps, pp(:,0,:,:))                      ; trick
   printVarSummary(pp)

print(" ")
print("---> a,b,p,u:  order is top-to-bottom ------")
print(sprintf("%8.5f", a)+"  "+sprintf("%8.5f", b)+"  " \
      +sprintf("%8.2f",pp(0,:,{0},{180}))          +"  " \
      +sprintf("%8.2f",ua(0,:,{0},{180})) )
print("------------------------------------------")
print(" ")

;************************************************
; calculate variable on pressure levels
; note, the 7th argument is not used, and so is set to 1.
;************************************************

   up = vinth2p(ua,a,b,pnew,ps,interp,p0new,1,extrap)
   copy_VarAtts(ua, up)
print(" ")
print("======================================")
printVarSummary(up)
print(" ")
print("pnew="+pnew+"   up="+up(0,:,{0},{180}))
print("======================================")

   filedimdef(Fout, "time", -1, True)
   Fout->U = up

print(" ")
print("======================================")
print("min(ua)="+min(ua)+"   max(ua)="+max(ua))
print("min(up)="+min(up)+"   max(up)="+max(up))
print("======================================")

j=j+1
end do