;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;                                                               ;
;  Takes netCDF data (of WRF model output, for example) and     ;
;  creates a new netCDF file with pressure level coordinates    ;
;  containing height and 2 m air temperature data.  This is     ;
;  helpful for examining pressure level data in WRF model       ;
;  output which instead relies on eta vertical coordinates.     ;
;  The use of height and 2 m temperature data is specifically   ;
;  for PV analysis.  Other variables can be added easily.       ;
;                                                               ;
;  Ryan Clare                                                   ;
;  University of Wisconsin-Madison                              ;
;  November, 2018                                               ;
;                                                               ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;===============================================================;
;  User inputs                                                  ;
;    The following variables should be adjusted according to    ;
;    which file you will be obtaining the data from, as well as ;
;    the domain size and the length of the date/time string.    ;
;===============================================================;

original_file = "Case12_T-0_no-alter_2005-03-06_06:00:00"
npres = 19              ; Number of desired pressure levels
;nlat = 143              ; Meridional domain size
;nlon = 209              ; Zonal domain size
nDSL = 19               ; Size of the DateStrLen
nlat = 30
nlon = 58

; Which directory the new file will be stored in
dirname = "/air/rclare/WRF_OUTs/ncl_netCDF_test/"

; Name for new file (if "False", it will auto-generate)
filename = False


;===============================================================;
;  Extract/create new data fields                               ;
;===============================================================;

; Extract data from original file
a = addfile(original_file, "r")
t = a->Times
lat = wrf_user_getvar(a, "lat", -1)
lon = wrf_user_getvar(a, "lon", -1)
z = wrf_user_getvar(a, "height", -1)
T = a->T2               ; 2 m air temperature

; Create pressure levels
Peas = ispan(1000,100,50)
Ps = new((/1,19/),"float")
Ps(0,0) = 1000.
Ps(0,1) = 950.
Ps(0,2) = 900.
Ps(0,3) = 850.
Ps(0,4) = 800.
Ps(0,5) = 750.
Ps(0,6) = 700.
Ps(0,7) = 650.
Ps(0,8) = 600.
Ps(0,9) = 550.
Ps(0,10) = 500.
Ps(0,11) = 450.
Ps(0,12) = 400.
Ps(0,13) = 350.
Ps(0,14) = 300.
Ps(0,15) = 250.
Ps(0,16) = 200.
Ps(0,17) = 150.
Ps(0,18) = 100.

; Add base pressure to pressure perturbation to get actual
pbs = a->PB
pp = a->P
P = pbs+pp

intopts = True
intopts@extrapolate = True
intopts@field_type = "z"

hgtFIELD = wrf_user_vert_interp(a,z,"pres",Peas,intopts)


;===============================================================;
;  Interpolate to rectilinear grid                              ;
;===============================================================;

lonD1 = fspan(-69.5,-126.5,nlon)
latD1 = fspan(24.5,53.5,nlat)

rekt_hgt = rcm2rgrid_Wrap(lat(0,:,:),lon(0,:,:),hgtFIELD,latD1,lonD1,1)
;rekt_lon = rcm2rgrid_Wrap(lat(0,:,:),lon(0,:,:),lon,latD1,lonD1,1)
;rekt_lat = rcm2rgrid_Wrap(lat(0,:,:),lon(0,:,:),lat,latD1,lonD1,1)
rekt_lon = new((/1,nlat,nlon/),"float")
rekt_lon(0,:,:) = 
T_rekt = rcm2rgrid_Wrap(lat(0,:,:),lon(0,:,:),T,latD1,lonD1,1)


;===============================================================;
;  Create new netCDF file                                       ;
;===============================================================;

if (filename .eq. False) then
 strs = str_split(original_file, "_")
 case_no = strs(0)
 case_date = str_join(strs(3:), "_")
 filo = case_no + "_" + case_date + ".nc"
else
 filo = filename
end if

diro = dirname
fout = addfile(diro + filo, "c")

setfileoption(fout,"DefineMode",True)


;===============================================================;
; Global attributes                                             ;
;===============================================================;

fAtt = True
fAtt@title = "netCDF test"
fAtt@source_file = original_file
fAtt@Conventions = "None"
fAtt@creation_date = systemfunc("date")
fileattdef(fout,fAtt)


;===============================================================;
; Variable attributes                                           ;
;===============================================================;

dimNames = (/"Time", "DateStrLen", "SOUTH_NORTH", "WEST_EAST", "pres"/)  
dimSizes = (/ -1   ,  nDSL, nlat,  nlon, npres /) 
dimUnlim = (/ True , False, False, False, False/)   
filedimdef(fout,dimNames,dimSizes,dimUnlim)

filevardef(fout,"Time",typeof(t),getvardims(t))
filevardef(fout,"pres",typeof(Ps),(/"Time", "pres"/))
filevardef(fout,"lat",typeof(lat),getvardims(rekt_lat))
filevardef(fout,"lon",typeof(lon),getvardims(rekt_lon))
filevardef(fout,"hgt",typeof(lon),(/"Time", "pres", "SOUTH_NORTH", "WEST_EAST"/))
;filevardef(fout,"hgt",typeof(hgtFIELD),getvardims(rekt_hgt))
filevardef(fout,"temp",typeof(T),getvardims(T_rekt))

PresAtt=0
PresAtt@units = "hPa"
PresAtt@long_name = "Pressure Levels"
PresAtt@MemoryOrder = "XY"
PresAtt@coordinates = "XLONG XLAT XTIME"
filevarattdef(fout, "pres", PresAtt )

;HgtAtt=0
;HgtAtt@units = "m"
;HgtAtt@long_name = "Geopotential height"
;HgtAtt@MemoryOrder = "XYZ"
;HgtAtt@coordinates = "XLONG XLAT XPRES XTIME"
;filevarattdef(fout, "hgt", HgtAtt)
filevarattdef(fout, "hgt", rekt_hgt)

filevarattdef(fout,"Time",t)
filevarattdef(fout,"lon",rekt_lon)
filevarattdef(fout,"lat",rekt_lat)
filevarattdef(fout,"temp",T)

setfileoption(fout,"DefineMode",False)


;===============================================================;
; Writing to the new file                                       ;
;===============================================================;

fout->Time   = (/t/)
fout->lon    = (/rekt_lon/)
fout->lat    = (/rekt_lat/)
fout->pres   = (/Ps/)
fout->hgt    = (/rekt_hgt/)
fout->temp   = (/T_rekt/)