'reinit'
'run /usr/people/ppantina/scripts/common/tempcolors.gs'

**Calculate the stddev of wet/dry years since you were too lazy a LONG time**
**ago to calculate it and write it to the data file (oversight!!!)**
datatype = 1
WHILE (datatype <=3)
  IF (datatype =1)
    dataset = indian
    name = r2
  ENDIF
  IF (datatype =2)
    dataset = cfsr
    name = cfsr
  ENDIF
  IF (datatype =3)
    dataset = merra
    name = merra
  ENDIF
  
  'open 'dataset'.ctl'
  type = 2
  WHILE (type <=3)
    IF (type = 2)
      type1 = dry
    ENDIF
    IF (type = 3)
      type1 = wet
    ENDIF 

    counter = 1
    WHILE (counter <= 5)
      IF (counter = 1)
        'define sum = pow((saavg'type''counter' - avg'type1'),2)'
      ELSE
        'define sum = sum + pow((saavg'type''counter' - avg'type1'),2)'
      ENDIF
      counter = counter + 1
    ENDWHILE **counter

    'define stddev'type1''name' = sqrt((sum)/(5-1))'

    type = type +1
  ENDWHILE **type 

  datatype = datatype +1

  'close 1'
ENDWHILE **datatype



'open indian.ctl'
'open cfsr.ctl'
'open merra.ctl'
'open evapsource.ctl'
'set lat -10 40'
'set lon 30 110'
'set mpdset hires'
'set display color white'
'set grid off'
'c'

**************************
**Seasonal Average FIRST**
**************************
'define mask = const(maskout(average.4, average.4-5),0,-u)'

'set parea 1 4.25 6.75 9.75'
'set gxout shaded'
'set grads off'
'set clevs 25 50 75 100 125 150 175 200'
'set ccols  81 39 37 35 21 22 24 26 28'
'd avgrain.1'
'set gxout contour'
'set clab off'
'set ccolor 1'
'set clevs 25 50 75 100 125 150 175 200'
'd avgrain.1'
'set gxout shaded'
'run images_corners.gs' 

'set parea 1 4.25 3.75 6.75'
'set grads off'
'set clevs 25 50 75 100 125 150 175 200'
'set ccols  81 39 37 35 21 22 24 26 28'
'd avgrain.2'
'set gxout contour'
'set clab off'
'set ccolor 1'
'set clevs 25 50 75 100 125 150 175 200'
'd avgrain.2'
'set gxout shaded'
'run images_corners.gs'

'set parea 1 4.25 .75 3.75'
'set grads off'
'set clevs 25 50 75 100 125 150 175 200'
'set ccols  81 39 37 35 21 22 24 26 28'
'd avgrain.3'
'set gxout contour'
'set clab off'
'set ccolor 1'
'set clevs 25 50 75 100 125 150 175 200'
'd avgrain.3'
'set gxout shaded'
'run images_corners.gs'
'run /usr/people/ppantina/scripts/common/cbarn.gs .5 0 2.5 .5'
'draw string 4.25 0.5 mm'

m=15
n=45

'define s = (('m'-1)*pow(stddevrain.1,2))+(('n'-1)*pow(stddev.4,2)))/('m'+'n'-2)'
'define test = (avgrain.1-average.4-0)/sqrt((s*(1/'m'+1/'n')))'
'define difference = 100*(avgrain.1-average.4)/(mask)'
'define sigdiff = maskout(difference, abs(test)-1.67)'  
'set parea 4.75 8 6.75 9.75'
'set grads off'
'set clevs -100 -50 -25 -10 -5 5 10 25 50 100 '
'set ccols 49 47 45 43 41 81 21 23 25 27 29 '
'd difference'
'hatch sigdiff -300 300 -density 0.1 -int 0.05 -angle 45,135'
'run images_corners.gs'

'define s = (('m'-1)*pow(stddevrain.2,2))+(('n'-1)*pow(stddev.4,2)))/('m'+'n'-2)'
'define test = (avgrain.2-average.4-0)/sqrt((s*(1/'m'+1/'n')))'
'define difference = 100*(avgrain.2-average.4)/(mask)'
'define sigdiff = maskout(difference, abs(test)-1.67)' 
'set parea 4.75 8 3.75 6.75'
'set grads off'
'set clevs -100 -50 -25 -10 -5 5 10 25 50 100 '
'set ccols 49 47 45 43 41 81 21 23 25 27 29 '
'd difference'
'set line 1 1 1'
'hatch sigdiff -300 300 -density 0.1 -int 0.05 -angle 45,135'
'run images_corners.gs'

'define s = (('m'-1)*pow(stddevrain.3,2))+(('n'-1)*pow(stddev.4,2)))/('m'+'n'-2)'
'define test = (avgrain.1-average.4-0)/sqrt((s*(1/'m'+1/'n')))'
'define difference = 100*(avgrain.3-average.4)/(mask)'
'define sigdiff = maskout(difference, abs(test)-1.67)'   
'set parea 4.75 8 0.75 3.75'
'set grads off'
'set clevs -100 -50 -25 -10 -5 5 10 25 50 100 '
'set ccols 49 47 45 43 41 81 21 23 25 27 29 '
'd difference'
'set line 1 1 1'
'hatch sigdiff -300 300 -density 0.1 -int 0.05 -angle 45,135'
'run images_corners.gs'
'run /usr/people/ppantina/scripts/common/cbarn.gs .5 0 6.5 .5'
'draw string 8.25 0.5 %'

**Titles and things**
'set string 1 tc 10 0'
'set strsiz .2'
'draw string 4.25 10.7 Seasonal'
'draw string 4.25 10.35 Evaporative Source'
'set strsiz .15'
'set string 1 tc 10 0'
'draw string 2.625 10 Seasonal Average'
'draw string 6.375 10 Difference'
'draw string 1 9.75 a)'
'draw string 1 6.75 b)'
'draw string 1 3.75 c)'
'draw string 4.75 9.75 d)'
'draw string 4.75 6.75 e)'
'draw string 4.75 3.75 f)'
'set string 1 tc 10 90'
'draw string 0.25 8.25 R2'
'draw string 0.25 5.25 CFSR'
'draw string 0.25 2.25 MERRA'
'set string 1 tc 10 0'
'printim images/paper_evap_1.png x1024 y1280'



pull pause

'c'

m=5
n=45

'set parea 1 4.25 6.75 9.75'
'define s = (('m'-1)*pow(stddevwetr2,2))+(('n'-1)*pow(stddev.4,2)))/('m'+'n'-2)'
'define test = (avgwet.1-average.4-0)/sqrt((s*(1/'m'+1/'n')))'
'define difference = 100*(avgwet.1-average.4)/(mask)'
'define sigdiff = maskout(difference, abs(test)-1.684)'  
'set grads off'
'set clevs -100 -50 -25 -10 -5 5 10 25 50 100 '
'set ccols 49 47 45 43 41 81 21 23 25 27 29 '
'd difference'
'hatch sigdiff -300 300 -density 0.1 -int 0.05 -angle 45,135'
'run images_corners.gs'
'set line 1 1 1'

'set parea 1 4.25 3.75 6.75'
'define s = (('m'-1)*pow(stddevwetcfsr,2))+(('n'-1)*pow(stddev.4,2)))/('m'+'n'-2)'
'define test = (avgwet.2-average.4-0)/sqrt((s*(1/'m'+1/'n')))'
'define difference = 100*(avgwet.2-average.4)/(mask)'
'define sigdiff = maskout(difference, abs(test)-1.684)'  
'set grads off'
'set clevs -100 -50 -25 -10 -5 5 10 25 50 100 '
'set ccols 49 47 45 43 41 81 21 23 25 27 29 '
'd difference'
'hatch sigdiff -300 300 -density 0.1 -int 0.05 -angle 45,135'
'run images_corners.gs'
'set line 1 1 1'

'set parea 1 4.25 .75 3.75'
'define s = (('m'-1)*pow(stddevwetmerra,2))+(('n'-1)*pow(stddev.4,2)))/('m'+'n'-2)'
'define test = (avgwet.3-average.4-0)/sqrt((s*(1/'m'+1/'n')))'
'define difference = 100*(avgwet.3-average.4)/(mask)'
'define sigdiff = maskout(difference, abs(test)-1.684)'  
'set grads off'
'set clevs -100 -50 -25 -10 -5 5 10 25 50 100 '
'set ccols 49 47 45 43 41 81 21 23 25 27 29 '
'd difference'
'hatch sigdiff -300 300 -density 0.1 -int 0.05 -angle 45,135'
'run images_corners.gs'
'set line 1 1 1'

'set parea 4.75 8 6.75 9.75'
'define s = (('m'-1)*pow(stddevdryr2,2))+(('n'-1)*pow(stddev.4,2)))/('m'+'n'-2)'
'define test = (avgdry.1-average.4-0)/sqrt((s*(1/'m'+1/'n')))'
'define difference = 100*(avgdry.1-average.4)/(mask)'
'define sigdiff = maskout(difference, abs(test)-1.684)'  
'set grads off'
'set clevs -100 -50 -25 -10 -5 5 10 25 50 100 '
'set ccols 49 47 45 43 41 81 21 23 25 27 29 '
'd difference'
'hatch sigdiff -300 300 -density 0.1 -int 0.05 -angle 45,135'
'run images_corners.gs'
'set line 1 1 1'

'set parea 4.75 8 3.75 6.75'
'define s = (('m'-1)*pow(stddevdrycfsr,2))+(('n'-1)*pow(stddev.4,2)))/('m'+'n'-2)'
'define test = (avgdry.2-average.4-0)/sqrt((s*(1/'m'+1/'n')))'
'define difference = 100*(avgdry.2-average.4)/(mask)'
'define sigdiff = maskout(difference, abs(test)-1.684)' 
'set grads off'
'set clevs -100 -50 -25 -10 -5 5 10 25 50 100 '
'set ccols 49 47 45 43 41 81 21 23 25 27 29 '
'd difference'
'set line 1 1 1'
'hatch sigdiff -300 300 -density 0.1 -int 0.05 -angle 45,135'
'run images_corners.gs'
'set line 1 1 1'

'set parea 4.75 8 0.75 3.75'
'define s = (('m'-1)*pow(stddevdrymerra,2))+(('n'-1)*pow(stddev.4,2)))/('m'+'n'-2)'
'define test = (avgdry.3-average.4-0)/sqrt((s*(1/'m'+1/'n')))'
'define difference = 100*(avgdry.3-average.4)/(mask)'
'define sigdiff = maskout(difference, abs(test)-1.684)'   
'set grads off'
'set clevs -100 -50 -25 -10 -5 5 10 25 50 100 '
'set ccols 49 47 45 43 41 81 21 23 25 27 29 '
'd difference'
'set line 1 1 1'
'hatch sigdiff -300 300 -density 0.1 -int 0.05 -angle 45,135'
'run images_corners.gs'
'set line 1 1 1'

**Titles and things**
'set string 1 tc 10 0'
'set strsiz .2'
'draw string 4.25 10.7 Interannual'
'draw string 4.25 10.35 Evaporative Source'
'run /usr/people/ppantina/scripts/common/cbarn.gs 1 0 4.25 .5'
'draw string 8.25 0.5 %'
'set strsiz .15'
'set string 1 tc 10 0'
'draw string 2.625 10 Wet Years'
'draw string 6.375 10 Dry Years'
'draw string 1 9.75 a)'
'draw string 1 6.75 b)'
'draw string 1 3.75 c)'
'draw string 4.75 9.75 d)'
'draw string 4.75 6.75 e)'
'draw string 4.75 3.75 f)'
'set string 1 tc 10 90'
'draw string 0.25 8.25 R2'
'draw string 0.25 5.25 CFSR'
'draw string 0.25 2.25 MERRA'
'set string 1 tc 10 0'
'printim images/paper_evap_2.png x1024 y1280'

pull pause