*------------------------------------------------
* Here we are dealing with the point data only
'reinit'
'open var2_ens.ctl'
*
'!rm temp*.txt'
'!rm formatted*.txt'
*
'set e 1 30'
t=1
while(t<=40)
'set t 't
'set gxout print'
'set prnopts %0.2f 20 1'
'd tmp2m'
say result
dummy=write('temp.txt',result)
dd=close('temp.txt')
* read 2 times to set the read pointer to 2nd line
rstr=read('temp.txt')
rstr=read('temp.txt')
pp=sublin(rstr,2)
ensm=mean(pp)
p10=percentile(pp,0.1)
q1=percentile(pp,0.25)
q2=percentile(pp,0.5)
q3=percentile(pp,0.75)
p90=percentile(pp,0.9)
p95=percentile(pp,0.95)
mx=max(pp)
mn=min(pp)
vals=mx' ' mn' ' q1' ' q2' ' q3' ' p10' ' p90' ' p95' ' ensm
'set prnopts %0.2f 8 1'
dummy1=write('temp22.txt',vals,append)
* say mx' ' mn' ' q1' ' q2' ' q3
dd=close('temp.txt')
* '!rm temp.txt'
t=t+1
endwhile
*
'quit'
*
****************************************************************
*
*###############          STATPACK v1.0          ###############
*
****************************************************************
*  
* ----------------------------------------
* STATPACK v1.0: GrADS-aholic October 2013
* Report bugs to http://gradsaddict.blogspot.com/
* ----------------------------------------
*
* Description: A collection of functions designed to perform 
* basic statistical analysis on user input data arrays.
*
* User Instructions:
* ------------------
* 1. Copy paste the contents of this file into the bottom your GrADS script.
* 2. Arrays must take the format of a string, with array elements separated by spaces.
*    e.g., arr='1 2 3 4 5 6 7 8'
* 3. List of included functions below.
*
* ------------------
* List of Functions:
* ------------------
* mean(arr): Returns mean of input array "arr"
* stdev(arr): Returns Standard Deviation of input array "arr"
* sort(arr): Returns sorted array (lowest to highest) of array "arr"
* rank(arr): Returns an array size equal to array "arr" ranging from 1 to size(arr)
* percentile(arr,p): Returns the data percentile (p) of array "arr"
* size(arr): Returns the number of elements in array "arr"
* max(arr): Returns the maximum value in array "arr"
* min(arr): Returns the minimum value in array "arr"
* correlation(arr1,arr2): Returns the correlation coefficient "r" between "arr1" and "arr2"
* regression(arr1,arr2): Returns the regression coefficient between "arr1" and "arr2"
* ------------------
* ------------------
*
* Example: Calculate the median of the array '9 7 8 6 5 4 3 2 1' using the percentile function.
* ------------------
* data='9 7 8 6 5 4 3 2 1'
* median=percentile(data,0.5)
* say median
* 
* the value returned in variable "median" is equal to 5 in this example.
* ------------------


****************************************************************
*                                                              *
*                         FUNCTIONS                            *
*                                                              *
****************************************************************




function mean(arr)
  sum=0
  ck=1;c=1
  while(ck=1)
    num=subwrd(arr,c)
    chknum=valnum(num)
    if(chknum=0);ck=0;else;sum=sum+num;endif
    c=c+1
  endwhile
  mean=sum/(c-2)
  return mean

******************

function stdev(arr)

  mean=mean(arr)

  sum_err_sq=0
  ck=1;c=1
  while(ck=1)
    num=subwrd(arr,c)
    chknum=valnum(num)
    if(chknum=0);ck=0;else
      sum_err_sq=sum_err_sq+math_pow(num-mean,2)
    endif
    c=c+1
  endwhile

      stdev=sum_err_sq/(c-3)
      stdev=math_sqrt(stdev)

  return stdev

******************

function sort(arr)
  ck=1;c=1
  while(ck=1)
    num=subwrd(arr,c)
    chknum=valnum(num)
    if(chknum=0);ck=0;endif
    c=c+1
  endwhile
  size=c-2

  ck=1
  cstart=1
  while(ck<=size)

    if(ck>1);arr=sorted_arr;endif
    sorted_arr=''
    c=cstart
    while(c<=size)
      n=c+1
      val1=subwrd(arr,c)
      val2=subwrd(arr,c+1)
     
      if(val1 > val2)
        val.c=val2;val.n=val1
      else
        val.c=val1;val.n=val2
      endif

      c=c+2
    endwhile

    c=size-1;n=size
    if(ck=1 & valnum(size/2)=1);val1=subwrd(arr,size-1);else;val1=val.c;endif
    val2=subwrd(arr,size)
    if(val1 > val2)
      val.c=val2;val.n=val1
    else
      val.c=val1;val.n=val2
    endif 

    c=1
    while(c<=size)
      sorted_arr=sorted_arr' 'val.c
      c=c+1
    endwhile
 
    if(cstart=1);cstart=2;else;cstart=1;endif
    ck=ck+1
  endwhile
 
  return sorted_arr

********************

function rank(arr)
  rank_arr=''
  c=1
  ck=1
  while(ck=1)
    num=subwrd(arr,c)
    chknum=valnum(num)
    if(chknum=0);ck=0;else;rank_arr=rank_arr' 'c;endif
    c=c+1
  endwhile
  
  return rank_arr

********************

function percentile(arr,p)
  
  if(p > 1 | p < 0)
    say "Not a valid number, please choose a number between 0 and 1."
  else

    sorted=sort(arr)
    size=size(arr)

    q1=p*(size-1)+1

    IR=math_int(q1)
    FR=q1-math_int(q1)

    if(q1=0);q1=1;endif

    if(q1=size);q1=subwrd(sorted,size);else
      q1=subwrd(sorted,math_int(q1))+(subwrd(sorted,math_int(q1)+1)-subwrd(sorted,math_int(q1)))*(q1-math_int(q1))
    endif

    percentile=q1

    return percentile
    
 endif

*******************


function size(arr)
  ck=1;c=1
  while(ck=1)
    num=subwrd(arr,c)
    chknum=valnum(num)
    if(chknum=0);ck=0;endif
    c=c+1
  endwhile
  size=c-2

return size

********************


function max(arr)
  size=size(arr)
  sorted=sort(arr)
  max=subwrd(sorted,size)
return max

********************

function min(arr)
  sorted=sort(arr)
  min=subwrd(sorted,1)
return min

********************

function correlation(arr1,arr2)


    size1=size(arr1)
    size2=size(arr2)
    if(size1 != size2)
      say "Correlation Coefficient must be determined from two equal sized vectors."
    else

      xmean=mean(arr1)
      ymean=mean(arr2)

      t=1
      denom=0
      numer=0
      ys=0
  
      while(t<=size1)
        xnum=subwrd(arr1,t)
        ynum=subwrd(arr2,t)
        denom=(xnum-xmean)*(xnum-xmean)+denom
        numer=(xnum-xmean)*(ynum-ymean)+numer
        ys=(ynum-ymean)*(ynum-ymean)+ys
        t=t+1
      endwhile

      r2=(numer*numer)/(denom*ys)
      r=math_sqrt(r2)

      return r
    endif

********************

function regression(arr1,arr2)


    size1=size(arr1)
    size2=size(arr2)

    if(size1 != size2)
      say "Regression Coefficient must be determined from two equal sized vectors."
    else

      xmean=mean(arr1)
      ymean=mean(arr2)
      t=1
      denom=0
      numer=0
      ys=0
  
      while(t<=size1)
        xnum=subwrd(arr1,t)
        ynum=subwrd(arr2,t)
        denom=(xnum-xmean)*(xnum-xmean)+denom
        numer=(xnum-xmean)*(ynum-ymean)+numer
        ys=(ynum-ymean)*(ynum-ymean)+ys
        t=t+1
      endwhile

      reg=numer/denom
      return reg
    endif

********************


****************************************************************
*                                                              *
*                        END FUNCTIONS                         *
*                                                              *
****************************************************************


****************************************************************
*
*
*
*###############                                 ###############
*
*
*
****************************************************************