;************************************************
;NCL tutorial script: sigma2pressure.ncl 
;************************************************
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/csm/shea_util.ncl"
;************************************************
begin
;************************************************

;************************************************
; file handling
;************************************************
   in  = addfile("ha.flm_ptl.OMEGA.14751525.jfm.nc","r")        ; open netcdf file

   system ("/usr/bin/rm ha.flm_ptl.OMEGA-p.14751525.jfm.nc")    ; remove if exist 
   out = addfile("ha.flm_ptl.OMEGA-p.14751525.jfm.nc","c")      ; open a new netcdfile

;************************************************
; read needed variables from file
;************************************************
   otime=in->time
   olat=in->lat
   olon=in->lon
   P0mb =1000.
   hyam = in->hyam                              ; get a coefficiants
   hybm = in->hybm                              ; get b coefficiants
   PS   = in->PS                                ; get pressure

;*********************************************
; chose intepolate type
; ptype = 1 specify desired pressure levels
; ptype = 2 increase pressure by a certain value
; vtype = 1 specify desired varibales
; vtype = 2 chose all variables in the file 
; note:
; if you chose ptype = 1 you should define lev1 in part 1
; if you chose ptype = 2 you should define PB,PT and Pspan in part 2
; if you chose vtype = 1 you should define vname1 in part 3
;***********************************************
  ptype = 1 
  vtype = 2 

;************************************
; 1. specify desired prssure levels
;*************************************
   lev1 = (/ 200.0,300.0,500.0,700.0,850.0 /)

;************************************
; 2. increase pressur by a certain value
;************************************
   PB = 900.0                  ; bottom pressure
   PT = 100.0                  ; top pressure
   Pspan = 50.0               ; the value that increased 
   Plevs = floattoint((PB-PT)/Pspan) + 1 
   lev2 = fspan(PT,PB,Plevs)

;**************************************
; 3. specify the desired variables
;**************************************
   vname1 = (/"OMEGA","PS"/)

;********************************************
; 4. select all the variables in the file
;********************************************
   vname2 = getfilevarnames(in)    

;*****************************************
; end chose intepolate type not change below
;*****************************************
   if (ptype .eq. 1) then
      lev = lev1
   else
      lev = lev2
   end if
   if (vtype .eq. 1) then
      vname = vname1
   else
      vname = vname2
   end if
;**********************************************
; define lev attributes
;**********************************************
   lev!0         = "lev"                    ; variable/dim name
   lev&lev     =  lev                       ; create coordinate variable
   lev@long_name = "pressure"               ; attach some attributes
   lev@units     = "hPa"
   lev@positive  = "down"

;***********************************************
; define other arguments required by vinth2p
;************************************************
; type of interpolation: 1 = linear, 2 = log, 3 = loglog
   interp = 2 

; is extrapolation desired if data is outside the range of PS
   extrap = False
;************************************************
; calculate variable on pressure levels
; note, the 7th argument is not used, and so is set to 1.
;************************************************
; write dimension variables to the output file
;************************************************
   filedimdef(out,"time",-1,True) ; define time unlimited
   out->time = otime
   out->lat = olat
   out->lon = olon
   out->lev = lev

;************************************************
; calculate variable on pressure levels and write it to
; the output file.
; note, the 7th argument is not used, and so is set to 1.
;************************************************
   do i=0,dimsizes(vname)-1
     dims  = getfilevardims(in,vname(i))
     nrank = dimsizes(dims)               ; determine the number of dimensions
     if (nrank .gt. 1) then               ; not do dimention variables
        if (nrank.eq.4) then              ; only interp 4D variables
          Vold = in->$vname(i)$           ; read each variable to memory
          Vnew  = vinth2p(Vold(:,:,:,:),hyam,hybm,lev,PS(:,:,:),interp,P0mb,1,extrap)
          copy_VarAtts (Vold,Vnew)
          Vnew@missing_value = Vnew@_FillValue
          delete(Vnew@_FillValue)
          delete(Vold)
        else
         Vnew  = in->$vname(i)$
         Vnew@missing_value = -9999.
        end if
        out->$vname(i)$ = Vnew 
        delete(Vnew)
      end if
      delete(dims)  
      delete(nrank)
   end do 
;*************************************
end
;*************************************