;=========================================
;loading the preamble for the ncl scripts
;=========================================
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/esmf/ESMF_regridding.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

;=============================================
; File name : extractTimeSeries.ncl

; purpose :

; Creation Date : 13-04-2016

; Last Modified :

; Created by : Nikhil
;=============================================
;begin the coding here
;=============================================
begin

dirnames = systemfunc("ls -d */")
dirDim = dimsizes(dirnames)

nDirs = dirDim(0)
if (nDirs .gt. 1) then
		do i=0,nDirs-1
			if (dirnames(i) .eq. "ts_output/") then
				system("rm -rf ./ts_output")
		  end if
		end do
else
	if (dirnames(0) .eq. "ts_output/") then
		system("rm -rf ./ts_output")
  end if
end if

system("mkdir ./ts_output")

;this loads all the files, in a list, this is a new feature in ncl
files = systemfunc("ls *.nc")
ncfile = addfiles(files, "r")

;get time from the list
times = ncfile[:]->Times

outTimes = wrf_times_c(times, 1)
outTimes2 = doubletoint(cd_convert(outTimes, "minutes since 1901-01-01 00:00:00"))

;print(outTimes2)

timeDim = dimsizes(times)

nsteps = timeDim(0)

;get the 2d coordinates
lat2d = ncfile[:]->XLAT(0, :, :)
lon2d = ncfile[:]->XLONG(0, :, :)

;not sure why i read the landmask
landmask = ncfile[:]->LANDMASK(0, :, :)

dim2D = dimsizes(lat2d)

;again not sure why i make this information available
nx = dim2D(1)
ny = dim2D(0)

;get the location from the tslist file

if (fileexists("tslist")) then
		tslist = readAsciiTable("tslist", 4, "float", 3)
end if
tsList = tslist(:,1:3)
;print(tsList)

;'UST','U10','V10','T2','Q2','SLP'

;get nearest i, j coordinate for a given lat, lon
location = getind_latlon2d(lat2d, lon2d, tsList(:,1), tsList(:,2))

locDims = dimsizes(location)
nLocs = locDims(0)

locNames = new(nLocs, "string")

if (fileexists("tslist")) then
		tsNames = readAsciiTable("tslist", 1, "string", 3)
end if

;make a array with name of output locations, which will be
;used to write the output of this script
do i=0,nLocs-1
	station_name = stringtochar(tsNames(i, :))
	locNames(i) = chartostring(station_name(:4))
	delete([/station_name/])
end do

tsData = new((/nsteps, nLocs, 6/), "float")

;get the data which is available in wrf output
do i=0,nsteps-1
		ust = ncfile[:]->UST(i, :, :)
		u10 = ncfile[:]->U10(i, :, :)
		v10 = ncfile[:]->V10(i, :, :)
		t2 = ncfile[:]->T2(i, :, :)
		do j=0, nLocs-1
			tsData(i, j, 0) = ust(location(j, 0), location(j, 1))
			tsData(i, j, 1) = u10(location(j, 0), location(j, 1))
			tsData(i, j, 2) = v10(location(j, 0), location(j, 1))
			tsData(i, j, 3) = t2(location(j, 0), location(j, 1))
		end do
end do

nfiles = ListCount(ncfile)

;calculate the parameters which are not in wrf output
do i=0, nfiles-1
		fileTimes = ncfile[i]->Times(:,:)
		fileTimeDim = dimsizes(fileTimes)
		nTimes = fileTimeDim(0)
		if (nTimes .gt. 1) then
				do j=0,nTimes-1
					;print((/i, j/))
					slp = wrf_user_getvar(ncfile[i], "slp", j)
				  rh2 = wrf_user_getvar(ncfile[i], "rh2", j)
					idx = j + i * nTimes
					do k=0, nLocs -1
						tsData(idx, k, 4) = slp(location(k, 0), location(k, 1))
						tsData(idx, k, 5) = rh2(location(k, 0), location(k, 1))
				  end do
				end do
		else
				  slp = wrf_user_getvar(ncfile[i], "slp", 0)
					rh2 = wrf_user_getvar(ncfile[i], "rh2", 0)
					do k=0, nLocs-1
						tsData(nsteps-1, k, 4) = slp(location(k, 0), location(k, 1))
						tsData(nsteps-1, k, 5) = rh2(location(k, 0), location(k, 1))
				  end do
		end if
		delete([/fileTimes, nTimes, fileTimeDim, slp, rh2/])
end do


;write the time series output to the text file
header = (/"---------------------------------------------------------------",\
						"    the time series data from wrfout using ncl and tslist     ",\
						"--------------------------------------------------------------",\
						"Time     U10        V10        Ust       RH2       T2      SLP"/)

hlist = [/header/]
do i=0, nLocs-1
		station_name = locNames(i) + ".txt"
		if (fileexists(station_name)) then
				system("rm ./"+station_name)
		end if
		temp = [/outTimes2, tsData(:, i, 1), tsData(:, i, 2), tsData(:, i, 0), \
				tsData(:, i, 5), tsData(:, i, 3), tsData(:, i, 4)/]
		write_table(station_name, "w", hlist, "%s")
		write_table(station_name, "a", temp, "%i,%3.6f,%3.6f,%3.6f,%3.6f,%3.6f,%4.6f")
		delete([/station_name/])
end do

;delete all the variable names
delete([/temp, hlist, header, tsData/])
delete([/ust, u10, v10, t2, tsList, nx, ny, dim2D, landmask, lon2d, lat2d/])

system("mv ./*.txt ./ts_output/.")

print("Script finished normally")

end; this ends begin