****************************************************************************************************
*
*        Program to extract 10 metres winds for selected stations (geographic points) from the Canadian 
*        Global Model of Environment Canada and write them toan ASCII file for use in loading data to 
*        a database.
*
*        Prepared by: M.A. MacLeod - Scotia Weather Services Inc - Halifax, Nova Scotia, Canada
*
*        Date: Dec 21. 2004
*
*****************************************************************************************************

*   Set initial counters

i = 1
t = 1
loopend = 1
stptm = 1
*   Read in the stations to be extracted from the GRIB fields
'open CMCHG.ctl'
stncount = read("awtsglbinput.txt")
rc = sublin(stncount,1) 
If (rc > 0)  
text1 = "Error in Input File - Code Number"%" "%rc
rc = write(cmcerror.txt,text1)
ENDIF

*   Read in the number of Stations to process

maxstn = sublin(stncount,2)
loopend = subwrd(maxstn,1)
* Read in the first Station to Process and preoduce the results to the output file
While (i <= loopend)

t = 1

stninfo = read("awtsglbinput.txt")

*  Read the input to check for errors in the reading
Filecode = sublin(stninfo,1)
 
If (Filecode > 0)  
text1 = "Error in Input File - Code Number"%" "%Filecode 
rc = write(cmcerror.txt,text1)
ENDIF

*  Extract the station information includidng name, latitude and longitude


stndata = sublin(stninfo,2)

stnname = subwrd(stndata,1)
lat1 = subwrd(stndata,2)
lon1 = subwrd(stndata,3)
strtm = subwrd(stndata,4)
stptm = subwrd(stndata,5)


val1 = "Stninfo = "%stndata%" Max Number of Stns = "%maxstn" Stn Id = "%stnname%" Latitude = "%lat1%" Longitude = "%lon1%" Start Time ="%strtm%" Stop Time ="%stptm

rc = write(cmcstnout.txt,val1)

fileset = stnname
'set lat  'lat1
'set lon 'lon1
'set mpdset hires '
 t = strtm 


say t

val3 = t" , "stptm
rc= write(cmcoutwrite.txt,val3)

* Loop through the time steps for the station 

While t <= stptm

'set t 't
say t
time='query time'
time
time1 = substr(result, 8, 2)
time2 = substr(result, 11, 2)
time3 = substr(result, 13, 3)
time4 = substr(result, 16, 4)

'set lev 1000 '
'define t10 = t - 273.14'
'd t10 '
Tmp1k = subwrd(result,4)

'define spd10 = mag(u,v) '
'd spd10 '
WndSpd1k = subwrd(result,4)

'define d1 = -atan2(u,-v)*180/3.1416'
'define d2 = const(maskout(d1,d1),0,-u)'
'define d3 = const(360+maskout(d1,-d2),0,-u)'
'define dir1k = d2+d3'
'display dir1k'

Wnddir1k = subwrd(result,4)


val = stnname%","%time3%","%time2%","%time4%","%time1%":00,"%Wnddir1k%","%WndSpd1k%","%Tmp1k
rc = write(cmcglb1fcst.txt, val)
rc = write(cmcglb1arch.ach,val,append)

t = t + 1
endwhile

* Increment counter for Input of Stations
i = i + 1
endwhile


close(cmcglb1fcst.txt)

*'!ftp nimbus < /home/wx_developer/GrADS_Scripts/cmc_ftp'