;   Script to write data to WPS intermediate file using in-built ncl function
;   Only works for NCL versions 6.3 or later
;   Returns several attributes attached to the return variable, each one a string array of length nfields 
;https://www.climatescience.org.au/sites/default/files/werner_post_processing.pdf

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
  
  DATE             = "2016-12-29"
  WPS_IM_root_name = "FILE3D"
  output_file_name = WPS_IM_root_name + ":" + DATE
  earth_radius     = 6367.470
  f = addfile("/home/student/Dev_WRF/DATA/asr15km.anl.3D.20161229.nc","r")
  U = f->UU
  V = f->VV
  T = f->TT
  R = f->RH
  Z = f->GHT
  xlat = f->XLAT
;  xlat = f->XLAT(0,:,:)
  xlon = f->XLONG
;  xlon = f->XLONG(0,:,:)
;  u!0 = "Time" 
;  u!1 = "num_metgrid_levels" 
;  u!2 = "south_north" 
;  u!3 = "east_west"
  printVarSummary(U)

;  system("/bin/rm " + output_file_name)
;  system("/home/student/Dev_WRF/DATA/asr15km.anl.3D.20161229.nc")

  FIELD_T          = "TT"
  UNITS_T          = "K"
  DESC_T           = "Temperature"

  FIELD_U          = "UU"
  UNITS_U          = "m s-1"
  DESC_U           = "Zonal Wind Speed"

  FIELD_V          = "VV"
  UNITS_V          = "m s-1"
  DESC_V           = "Meridional Wind Speed"

  FIELD_R          = "RH"
  UNITS_R          = "%K"
  DESC_R           = "Relative Humidity"

  FIELD_Z          = "GHT"
  UNITS_Z          = "m"
  DESC_Z           = "Geopotential Height"

 t = (/ 00, 03, 06, 09, 12, 15, 18, 21/)
 

 pnew = (/ 1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 725, 700, 675, \
            650, 625, 600, 575, 550, 525, 500, 450, 400, 350, 300, 250, 200 /)*100
 
  ntimes = dimsizes(t)
  NLEV = dimsizes(pnew)

  opt                   = True
  opt@projection        = 0                 ; "Equidistant_Lat_Lon"
  opt@date              = DATE
  opt@map_source        = "ASRV2"
  opt@startloc          = "SWCORNER"	    ; 8 chars exact
  opt@startlon          = xlon(0,0)
  opt@startlat          = xlat(0,0)
  opt@deltalon          = xlon(0,1)-xlon(0,0)
  opt@deltalat          = xlat(1,0)-xlat(0,0)
  opt@is_wind_earth_rel = False

;  do time=0,ntimes-1
;  opt@hour = t(time)
  	do ilev=0,NLEV-1
    	opt@level = pnew(ilev)

    	wrf_wps_write_int(WPS_IM_root_name,FIELD_U,UNITS_U,\
        	              DESC_U,U(ilev,:,:),opt)
    	wrf_wps_write_int(WPS_IM_root_name,FIELD_V,UNITS_V,\
                	      DESC_V,VonP(ilev,:,:),opt)
    	wrf_wps_write_int(WPS_IM_root_name,FIELD_T,UNITS_T,\
                      	      DESC_T,TonP(ilev,:,:),opt)
    	wrf_wps_write_int(WPS_IM_root_name,FIELD_R,UNITS_R,\
        	              DESC_R,RonP(ilev,:,:),opt)
    	wrf_wps_write_int(WPS_IM_root_name,FIELD_Z,UNITS_Z,\
        	              DESC_Z,ZonP(ilev,:,:),opt)
  	end do
;  end do