;***************************************************************
;This program is used to calculate velocity potential and streamfunction.
;Model variables are converted from hybrid sigma to pressure coordinates.
;Output to a NetCDF file is then done.
;***************************************************************
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
;***************************************************************
; variable and file handling
;***************************************************************
  diri  = "/chinook/zmangin/isugcm/isugcm_tri65_hourly/"
  diro  = "/home/zmangin/isu_gcm/ncl_plots/"
  f     = addfile(diri+"jja_1980-89_climo.nc","r")
  filo  = "isu65_JJA_sf_vp_u_v_q_climo.nc"
  system("/bin/rm -f " + diro + filo)    ; remove if exists
  fout  = addfile(diro+filo,"c")
;*************************************************
;define necessary variables from original netcdf file
;*************************************************
  u     = f->U     ;(t,lev,lat,lon)
  v     = f->V     ;(t,lev,lat,lon)
  q     = f->Q     ;(t,lev,lat,lon)

  time  = f->time
  lev   = f->lev     ;hybrid sigma levels 
  lat   = f->lat     ;64
  lon   = f->lon     ;128

  hyam   = f->hyam
  hybm   = f->hybm
  ps     = f->PS
  p0     = f->P0
  p0     = p0/100.

  ntime = dimsizes(time)
  nlev  = dimsizes(lev)
  nlat  = dimsizes(lat)
  nlon  = dimsizes(lon)
;***************************************************************  
; explicitly declare file definition mode. Improve efficiency.
;***************************************************************    
  setfileoption(fout,"DefineMode",True)
;***************************************************************  
; create global attributes of the file
;***************************************************************  
  fAtt               = True            ; assign file attributes
  fAtt@title         = "NCL Efficient Approach to netCDF Creation"  
  fAtt@source_file   = "/chinook/zmangin/isugcm/isugcm_tri65_hourly/jja_1980-89_climo.nc" 
  fAtt@Conventions   = "None"   
  fAtt@creation_date = systemfunc ("date")        
  fileattdef( fout, fAtt )            ; copy file attributes    
;***************************************************************
; predefine the coordinate variables and their dimensionality
; Note: to get an UNLIMITED record dimension, we set the dimensionality
; to -1 (or the actual size) and set the dimension name to True.
;*************************************************************** 
  dimNames = (/"time", "lat", "lon", "lev"/)  
  dimSizes = (/ -1   ,  nlat,  nlon, nlev /) 
  dimUnlim = (/ True , False, False, False/)    
  filedimdef(fout,dimNames,dimSizes,dimUnlim)
;*************************************************
;  calculates vars on new p surfaces
;************************************************
  interp = 2
  extrap = False
  
  lev_p  = new(16,typeof(lev),u@_FillValue)

  lev_p       = (/30.,50.,70.,100.,150.,200.,250.,300.,400.,500.,600.,700.,775.,850.,925.,1000./)
  lev_p@units = "hPa" 
  ;print(lev_p)
  
  nlev       = dimsizes(lev_p) 

  UonP       = vinth2p(u,hyam,hybm,lev_p,ps,interp,p0,1,extrap)
  VonP       = vinth2p(v,hyam,hybm,lev_p,ps,interp,p0,1,extrap)
  QonP       = vinth2p(q,hyam,hybm,lev_p,ps,interp,p0,1,extrap)

  ;printVarSummary(UonP)
  ;printVarSummary(VonP)
  ;printVarSummary(QonP)
  ;printMinMax(UonP,True)
;************************************************
;  calculates stream function and velocity potential
;************************************************
  uvmsg = 1e+36

  sf = new ( (/ntime,nlev,nlat,nlon /), float, uvmsg ) ; stream function
  vp = new ( (/ntime,nlev,nlat,nlon /), float, uvmsg ) ; velocity potential

  uv2sfvpg (UonP,VonP,sf,vp)
  ;printVarSummary(sf)
  ;printVarSummary(vp)

  ;copy_VarCoords(u, vp )
  ;copy_VarCoords(v, vp )
  sf@long_name  = "Streamfunction"
  sf@units      = "m2/s"
  sf!0          = "time"
  sf!1          = "lev_p"
  sf!2          = "lat"
  sf!3          = "lon"

  vp@long_name  = "Velocity Potential"
  vp@units      = "m2/s"
  copy_VarCoords(sf,vp)

  QonP@long_name  = "specific humidity"
  QonP@units      = "kg/kg"
  copy_VarCoords(sf,QonP)


  ;printVarSummary(lev_p)
  ;print(lev_p)
  ;printVarSummary(vp)
  ;print(vp)
;***************************************************************
;predefine the the dimensionality of the variables to be written out
;***************************************************************
; Here we are using NCL functions to facilitate defining 
; each variable's dimension name(s) and type. 
; The following could be replaced with explicit, user defined dimension 
; names different from those associated with the variable in memory. 
; Say, PS(time,lat,lon) in the NCL script. They could be redefined for the file via: 
; filevardef(fout, "PS"   ,typeof(PS) ,(/"TIME","latitude","longitude"/)) 
;*************************************************  
  filevardef(fout, "time"  ,typeof(time)   ,getvardims(time)) 
  ;filevardef(fout, "lev"  ,typeof(lev)    ,getvardims(lev))
  filevardef(fout, "lev_p" ,typeof(lev_p)  ,getvardims(lev_p))
  filevardef(fout, "lat"   ,typeof(lat)    ,getvardims(lat))                          
  filevardef(fout, "lon"   ,typeof(lon)    ,getvardims(lon))                          
  filevardef(fout, "sf"    ,typeof(sf)     ,getvardims(sf))    
  filevardef(fout, "vp"    ,typeof(vp)     ,getvardims(vp))                   
  filevardef(fout, "QonP"  ,typeof(QonP)   ,getvardims(QonP))  
  filevardef(fout, "UonP"  ,typeof(UonP)   ,getvardims(UonP))  ; different from name on script
  filevardef(fout, "VonP"  ,typeof(VonP)   ,getvardims(VonP))  ; different from name on script
;*************************************************
; Copy attributes associated with each variable to the file
; All attributes associated with each variable will be copied.
;*************************************************  
  filevarattdef(fout,"sf",sf)                         ; copy sf attributes
  filevarattdef(fout,"vp",vp)                         ; copy vp attributes
  filevarattdef(fout,"QonP",QonP)                     ; copy q attributes
  filevarattdef(fout,"UonP",UonP)                     ; copy u attributes
  filevarattdef(fout,"VonP",VonP)                     ; copy v attributes
  filevarattdef(fout,"time" ,time)                    ; copy time attributes
  ;filevarattdef(fout,"lev"  ,lev)                     ; copy lev attributes
  filevarattdef(fout,"lev_p",lev_p)                   ; copy lev_p attributes
  filevarattdef(fout,"lat"  ,lat)                     ; copy lat attributes
  filevarattdef(fout,"lon"  ,lon)                     ; copy lon attributes
;*************************************************
; explicitly exit file definition mode. **NOT REQUIRED**
;*************************************************  
 setfileoption(fout,"DefineMode",False)
;*************************************************
; output only the data values since the dimensionality and such have
; been predefined. The "(/", "/)" syntax tells NCL to only output the
; data values to the predefined locations on the file.
;*************************************************  
  fout->time    = (/time/)     
  ;fout->lev    = (/lev/)
  fout->lev_p   = (/lev_p/)
  fout->lat     = (/lat/)
  fout->lon     = (/lon/) 
  ;fout->sf     = (/sf/)
  ;fout->vp     = (/vp/)
  fout->QonP    = (/QonP/)
  fout->UonP    = (/UonP/)
  fout->VonP    = (/VonP/)
;*************************************************  
end