*TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
*     Author : Martin V. Mathew, vmartinmathew@gmail.com
*     Prepared on: 18 July 2013
*LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL 

'reinit'


*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* 1) Generating input file name & opening the file
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
infile1 ='/users/martin/work/research/input_data/wind_fields/wind_stress/tropflux/taux_tropflux_1d_1979_2012.nc'
mnCnt=0;
say 'Opening input file '%infile1
'sdfopen 'infile1

*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* 2) Using 'q file' command to know details of time dimension in the file
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*ga-> q file
*File 1 :
*  Descriptor: tauy_tropflux_1d_1979_2012.nc
*  Binary: tauy_tropflux_1d_1979_2012.nc
*  Type = Gridded
*  Xsize = 350  Ysize = 60  Zsize = 1  Tsize = 12266  Esize = 1
*  Number of Variables = 1
*     tauy  0  t,y,x  meridional wind stress 

* Using intrinsic function to read information about maximum no of time from result of 'q file'
* Refer http://grads.iges.org/grads/gadoc/script.html#functions
 '  q file '
   anslin5=sublin(result,5)

   maxday=subwrd(anslin5,12) ;

say anslin5
say maxday

* Remove following line after the trial, this line is included to overcome
* some error in the input data
maxday=7000 

*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* 3) Identifying start and end points of each month.
* Note: I am assuming that your data is a daily data for a single station.
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
days=1
strtMnT=days
while(days<=maxday)
'set t 'days
'set lat 0'
'set lon 70'

* Grads command 'q dim' is used to identify start and end points of each month 

*ga-> q dim
*Default file number is: 1
*X is varying   Lon = 30.5 to 379.5   X = 1 to 350
*Y is varying   Lat = -29.5 to 29.5   Y = 1 to 60
*Z is fixed     Lev = 0  Z = 1
*T is fixed     Time = 11:59Z01JAN1979  T = 1
*E is fixed     Ens = 1  E = 1

if days=1
'q dim'
   anslin5=sublin(result,5)
   prvM=substr(anslin5,31,3)
   yr=substr(anslin5,34,4)
else ;*if days=1
prvM=curM
   yr=substr(anslin5,34,4)
endif;*if days=1

'q dim'
   anslin5=sublin(result,5)
   curM=substr(anslin5,31,3) ;* Pls make necessary modification in the last two numbers
                            ;* to make sure that they correspond to text showing month.
                            ;* Numbers given here is valid for text JAN in the above eg  
if (prvM != curM)
stopMnT=days-1
mnCnt=mnCnt+1
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* 4) Variables 'strtMnT' & 'stopMnT' give the time indices showing start & stop of 
*    each month. Now use grads functions 'max'& 'min' or 'amax' & 'amin' to identify
*    maximum value for data of each month. 
* http://www.iges.org/grads/gadoc/gadocindex.html
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Here I am using 'max' & 'min' assuming that there is only one station in the data.
'define minVal= min(taux, t= 'strtMnT', t= 'stopMnT',1)'
'define maxVal= max(taux, t= 'strtMnT', t= 'stopMnT',1)'
'd minVal'; MinVal1=subwrd(result,4);
'd maxVal'; MaxVal1=subwrd(result,4);
*say 'minVal MinVal1 'minVal' 'MinVal1
*say 'maxVal MaxVal1 'maxVal' 'MaxVal1
* Saving the monthly maximum and minimum of data into a compound variable
OutAry.mnCnt.1=yr ;* saving year
OutAry.mnCnt.2=prvM  ;* saving text showing months JAN, FEB, etc
OutAry.mnCnt.3=MinVal1 ;* minimum value of the month
OutAry.mnCnt.4=MaxVal1 ;* minimum value of the month

strtMnT=stopMnT+1
endif ;*if (prvM != curM)

days=days+1
endwhile ;*while(days<=maxday)

'close 1' ;* closing file
* Computations for all years are completed.
 
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* 5) Writing the data in compound variable 'OutAry' to output file. 
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* terminal display of the calculated data
time=1
while(time<=mnCnt)
say OutAry.time.1'  'OutAry.time.2'  'OutAry.time.3'  'OutAry.time.4
time=time+1
endwhile ;*while(days<=maxday)