;This script is a simple example of how to read data in and write it
;out to a new .nc file using NCL. The .nc file output from this script
;should be readable by GrADS or any other software.

;Questions? srausche@students.wisc.ediu


begin

; open file to write variable to (file should not already exist)
out = addfile("pme.nc","c")

; open data file to get data from
data = addfile("/grove1/people/behling/gen2lsx.GLOBS60sv.012.nc", "r")

; Define the variables and their dimensions.
; Variables have dimensions of (time,level,latitude,longitude) or (time,latitude,longitude)
; depending on whether it is a multi-level or sfc variable

dims = dimsizes(data->PRECIP)
print(dims)

pme = new((/dims(0),dims(1),dims(2)/),float)

;the following defines the coordinate vars in "pme"
;and makes time UNLIMITED
;taken from script at http://www.cgd.ucar.edu/cas/ncl-csm/IO_f2nc_2.shtml

dimNames = (/"time", "lat", "lon"/)   ; coordinate variables
dimSizes = (/  -1  ,  dims(1),  dims(2)     /)   ; coordinate dimensions
dimUnlim = (/ True , False, False /)   ; define unlimited d
filedimdef(out, dimNames  , dimSizes,  dimUnlim )

;put var names you want to define as "varNames2D"
varNames2D        = (/ "pme" /)
varTypes2D        = (/ "float" /)
filevardef(out, varNames2D, varTypes2D, (/"time", "lat" , "lon"/) )

;define variables
pme = data->PRECIP(:,:,:) - data->AVAP(:,:,:) 
time = data->time
lat = data->lat
lon = data->lon

; Assign metadata to variables

pme!0 = "time"
pme&time = time 
pme!1 = "lat"
pme&lat = lat 
pme!2 = "lon"
pme&lon = lon 

; write out data to file
; write variables
out->pme = pme 

; write metadata 
out->lat = lat 
out->lon = lon 
out->time = time

end