*GrADS script to pull out a data from a station given the latitude and longitude
*by Lydia Rill 06/01/16
*format should be year,doy,srad,tmax,tmin,prcp,,,
function main(args) 
lat=subwrd(args,1)
lon=subwrd(args,2)
countdays=subwrd(args,3)
sID=subwrd(args,4)
sname=subwrd(args,5)
styear=subwrd(args,6)
stdoy=subwrd(args,7)
newtime=subwrd(args,8)

'reinit'
'open /Volumes/WebData/NLDAS2/DailyNLDAS/daily.ctl'
*say 'going to pull station data. newtime and countdays:' newtime' 'countdays
'set x 1'
'set y 1'
'set t 'newtime' 'countdays''
* missing data is marked as -9.99E8 in the grib files
'set undef -9999'
*float format, 1 value for each line, 1 space between values
'set gxout print'
'set prnopts %6.1f 6 1' 

t2=newtime
t1=0
'd gr2stn(TMAX,'lon','lat',-n)'
tmaxresult=result
'd gr2stn(TMIN,'lon','lat',-n)'
tminresult=result
'd gr2stn(PRCP,'lon','lat',-n)'
prcpresult=result
'd gr2stn(SRAD,'lon','lat',-n)'
sradresult=result
 the data for the lat lon' 

doy=stdoy
year=styear
prev=missing

while (t2<=countdays)

* this command creates station data by sampling a gridded dataset and interpolating to a given location 
*(interpolation done bi-linearly. uses 4 surrounding grid points)
* use gr2stn(b=data,lon,lat,-n) if you want the nearest  grid point instead
tmaxline=sublin(tmaxresult,4+2*t1)
tmax2=subwrd(tmaxline,1)
tminline=sublin(tminresult,4+2*t1)
tmin2=subwrd(tminline,1)
prcpline=sublin(prcpresult,4+2*t1)
prcp2=subwrd(prcpline,1)
sradline=sublin(sradresult,4+2*t1)
srad2=subwrd(sradline,1)

* Replace missing data with the average of the previous and next day. (except set precip to 0)
if (tmax2<-99)
prev=sublin(tmaxresult,4+(2*t1)-2)
p1=subwrd(prev,1)
next=sublin(tmaxresult,4+(2*t1)+2)
n1=subwrd(next,1)
if (t2<countdays)
  if (p1 > -99 & n1 > -99)
    tmax2 = (p1+n1)/2
  else
    tmax2=-9999
  endif
else
  tmax2=p1
endif
endif

if (tmin2<-99)
prev2=sublin(tminresult,4+(2*t1)-2)
p2=subwrd(prev2,1)
next2=sublin(tminresult,4+(2*t1)+2)
n2=subwrd(next2,1)
if (t2<countdays)
  if (p2 > -99 & n2 > -99)
    tmin2=(p2+n2)/2
  else
    tmin2=-9999
  endif
else
  tmin2=p2
endif
endif

if (prcp2<-99)
  prcp2=0
endif

if (srad2<-99)
prev4=sublin(sradresult,4+(2*t1)-2)
p4=subwrd(prev4,1)
next4=sublin(sradpresult,4+(2*t1)+2)
n4=subwrd(next4,1)
if (t2<countdays)
  if (p4 > -99 & n4 > -99)
    srad2=(p4+n4)/2
  else
    srad2=-9999
  endif
else
  srad2=p4
endif
endif

say year','doy','srad2','tmax2','tmin2','prcp2',,,'
t1=t1+1
t2=t2+1

*is this a leap year?
rc1=math_mod(year,400)
rc2=math_mod(year,4)
rc3=math_mod(year,100)
if ((rc1=0) | (rc2=0 & rc3!=0))
endday=366
else
endday=365
endif

if doy<endday
doy=doy+1
else
year=year+1
doy=1
endif

endwhile
'quit'