**This program calculates the evaporative source for different pentads for**
**different years, based on the information from files titled master_xxx**
**It then calculates an average of the 16-20 periods of active/break phases**
**and creates an image and also writes it to a file**

'reinit'

*Set color bar colors**
'set rgb 24 0 0 255'
'set rgb 23 20 20 255'
'set rgb 22 40 40 255'
'set rgb 21 60 60 255'
'set rgb 20 80 80 255'
'set rgb 19 100 100 255'
'set rgb 18 120 120 255'
'set rgb 17 140 140 255'
'set rgb 16 160 160 255'
'set rgb 26 180 180 225'
'set rgb 25 200 200 225'

*Set domain to India**
lon0 = 30
lon1 = 110
lat0 = -10
lat1 = 40

set = 1 
**neutral/dry/wet
active = 1 
**active/break

WHILE (active <=2)
  WHILE (set <=3)

****Conditional statements to get correct master file read**
    if (active = 1)
      if (set = 1)
        master = master_neu
      endif
      if (set = 2)
        master = master_dry
      endif
      if (set = 3)
        master = master_wet
      endif
    endif
    if (active = 2)
      if (set = 1)
        master = master_neu_break
      endif
      if (set = 2)
        master = master_dry_break
      endif
      if (set = 3)
        master = master_wet_break
      endif
    endif

    'c'

****Read the file 'master' for information about**
****interannual and intraseasonal variability**
    master1 = read(master)
    line = sublin(master1, 2)
    type = subwrd(line, 1)
    spell = subwrd(line,2)
    master2 = read(master)
    master3 = read(master)
    periods = subwrd(master3, 6)
    say 'This is for 'type' years with 'periods' 'spell' periods'

****This loops works through all periods of a given set of**
****wet/dry/neutral or active/break situations**
****The 'master' file is read to determine how many periods**
****there will be**
    counter = 1
    WHILE (counter <= periods)
      say '# 'counter
      periodfinder = read(master)
      t0 = subwrd(periodfinder, 2)
      tf = subwrd(periodfinder, 3)
      byear = subwrd(periodfinder, 4)
        
**********************************************************
******This is just the images_specific.gs program copied**
******into this program, because it would have been very**
******complicated to do it any other way, I believe*******
**********************************************************

      tmpoutfile = 'test.txt'

      gridarea1  = read(gridarea)
      filenumlin = sublin(gridarea1,2)
      filenum    = subwrd(filenumlin,3)
      gridcor1   = read(gridarea)
      cornumlin  = sublin(gridcor1,2)
      cornum     = subwrd(cornumlin,3)
      gridcor2   = read(gridarea)
      gridcor    = sublin(gridcor2, 2)
 
*    say 'Pentads used in tracing: 't0'-'tf
*    say 'Number of files/points = 'filenum
*    say 'Number of corners      = 'cornum

      xycount = 1
      WHILE (xycount <= filenum)
        xylin1 = read(gridarea)
        xylin  = sublin(xylin1, 2)
        x1     = subwrd(xylin, 1)  
        y1     = subwrd(xylin, 2)
        x      = substr(x1, 2, 3)
        y      = substr(y1, 2, 3)
        file   = 'source.'byear'.'x1'.'y1'.grd'
        line1  = 'dset /elninost2/ppantina/Traj/cfsr/output.'byear'/'file
        line2  = '***'
        '!/bin/rm source.ctl'
        err = write(tmpoutfile, line1)
        err = write(tmpoutfile, line2, append)
        err = close(tmpoutfile)
        '!cat 'tmpoutfile' /elninost2/ppantina/Traj/cfsr/output.'byear'/source.ctl.template > source.ctl'

*       say 'Create Control File: 'file
        'open source.ctl'
        'define ps = 't0
        'define pf = 'tf
        'define year = 'byear
        'set lat -60 60'
        'set lon   0 360'
        if (xycount = 1)
          'define saavg't0''byear' = sum(s, t='t0', t='tf')*5'
        else
          'define saavg't0''byear' = saavg't0''byear' + sum(s, t='t0', t='tf')*5'
        endif
 
       'close 1'
        xycount = xycount+1

      ENDWHILE **xycount
 
******Define a new variable equal to saavg that is consequecutive**
******so that it is easier to loop for the stddev calculation** 
      'open source.ctl'
      'define saavg'set''active''counter' = saavg't0''byear
      'close 1'
  
*      'c'
*      'set lat 'lat0' 'lat1
*      'set lon 'lon0' 'lon1
*      'd saavg't0''byear
*      'run images_corners.gs'
*      'draw title Evaporate Source for Pentads 't0' to 'tf' from 'byear

      err = close(gridarea)

*************************************
*************************************
**This is the end of the subroutine**
*************************************
*************************************

******Open the source.ctl file again to calculate an average, first by summing**
******up the individual 16-20 periods so that we can see how they look overall**
      'open source.ctl'
      IF (counter = 1)
        'define avg = saavg't0''byear  
      ELSE
        'define avg = avg + saavg't0''byear
      ENDIF
 
      'close 1'
      counter = counter +1 
    ENDWHILE **counter

    counter = 1
 
****Now calculate the average and standard deviation of the average variable**
    'open source.ctl'
    'define avg'set''active' = avg/'periods

****This loop sums up the saavg variables so that they may be subtracted**
****from the average and then summed for stddev**    
    WHILE (counter <= periods)
      IF (counter = 1)
        'define sum'set''active' = pow((saavg'set''active''counter' - avg'set''active'),2)'
      ELSE
        'define sum'set''active' = sum'set''active' + pow((saavg'set''active''counter' - avg'set''active'),2)'
      ENDIF
  
      counter = counter + 1
    ENDWHILE **counter(2)

****Define stddev**
    'define stddev'set''active' = sqrt((sum'set''active')/('periods'-1))'


   
****Display the average plots**
*    'set parea 1 10 1 7.5'
*    'set display color white' 
*    'c'
*    'set grid off'
*    'set grads off'
*    'set mpdset hires'
*    'set lon 'lon0' 'lon1
*    'set lat 'lat0' 'lat1
*    'set grads off'
*    'set xlint 10'
*    'set ylint 5'
*    'set gxout shaded'
*    'set clevs 0 5 10 15 20 25 30 35 40 45 50'
*    'set ccols  0  0 25 26 16 17  18  19  20  21  22  23 24'
    
*    'd avg'set''active

*    'run /usr/people/ppantina/scripts/common/cbarn.gs'
*    'run images_corners.gs'
*    'draw title Average Evaporate Source for 'spell' Periods in 'type' Years'
*    'printim images/'spell'_'type'_spells.png x1280 y1024'

**Write to a file to save them at the same time!**
    'set gxout fwrite'
   
    IF (set =1)
      IF (active = 1)
        'set fwrite /elninost2/ppantina/Traj/compare/cfsr_spells.grd'
      ELSE
        'set fwrite -ap /elninost2/ppantina/Traj/compare/cfsr_spells.grd'
      ENDIF
    ELSE
       'set fwrite -ap /elninost2/ppantina/Traj/compare/cfsr_spells.grd'
    ENDIF
    
    'set x 1 192'
    'set y 1 94'
    'd avg'set''active
    'd stddev'set''active
*    'set gxout shaded'

    'disable fwrite'
    'close 1'

    err = close(master)

    set = set +1

  ENDWHILE **set

  active = active+1

  set = 1

ENDWHILE **active


**Now do some caluclates to find the the average of **
**the active/break periods, as well as the stddev of the 6 cases**
'open source.ctl'
'set x 1 192'
'set y 1 94'
'define avgactive = (avg11+avg21+avg31)/3'
'define avgbreak = (avg12+avg22+avg32)/3'

**Here, counter 'm' is equivalent to 'set'**
**and counter 'n' is equivalent to 'active'**
**First calculate stddev of active periods**
m=1
WHILE (m<=3)
  IF(m=1)
    endn = 19
  ENDIF
  IF(m=2)
    endn = 19
  ENDIF
  IF(m=3)
    endn = 18
  ENDIF
  
  n=1
  WHILE (n<=endn)
    IF (n = 1)
      'define sum1 = pow((saavg'm'1'n' - avgactive),2)'
    ELSE
      'define sum1 = sum1 + pow((saavg'm'1'n' - avgactive),2)'
    ENDIF
    n = n+1
  ENDWHILE **n

m = m+1
ENDWHILE **m

m=1
WHILE (m<=3)
  IF(m=1)
    endn = 16
  ENDIF
  IF(m=2)
    endn = 17
  ENDIF
  IF(m=3)
    endn = 16
  ENDIF
  
  n=1
  WHILE (n<=endn)
    IF (n = 1)
      'define sum2 = pow((saavg'm'2'n' - avgbreak),2)'
    ELSE
      'define sum2 = sum2 + pow((saavg'm'2'n' - avgbreak),2)'
    ENDIF
    n = n+1
  ENDWHILE **n

m = m+1
ENDWHILE **m

**
**
**
**
**
m=1
WHILE (m<=3)
  IF(m=1)
    endn = 19
  ENDIF
  IF(m=2)
    endn = 19
  ENDIF
  IF(m=3)
    endn = 18
  ENDIF
 
  n=1
  WHILE (n<=endn)
    IF (n = 1 & m=1)
      'define sum = saavg'm'1'n
    ELSE
      'define sum = sum + saavg'm'1'n
    ENDIF
    n = n+1
  ENDWHILE **n

m = m+1
ENDWHILE **m

m=1
WHILE (m<=3)
  IF(m=1)
    endn = 16
  ENDIF
  IF(m=2)
    endn = 17
  ENDIF
  IF(m=3)
    endn = 16
  ENDIF
   
  n=1
  WHILE (n<=endn)
    'define sum = sum + saavg'm'2'n
  n = n+1
  ENDWHILE **n
m = m+1
ENDWHILE **m

'avgtotal = (sum)/105'

m= 1
WHILE (m<=3)
  IF(m=1)
    endn = 19
  ENDIF
  IF(m=2)
    endn = 19
  ENDIF
  IF(m=3)
    endn = 18
  ENDIF
 
  n=1
  WHILE (n<=endn)
    IF (n = 1 & m=1)
      'define sum3 = pow((saavg'm'1'n' - avgtotal),2)'
    ELSE
      'define sum3 = sum3 + pow((saavg'm'1'n' - avgtotal),2)'
    ENDIF
    n = n+1
  ENDWHILE **n

m = m+1
ENDWHILE **m

m=1
WHILE (m<=3)
  IF(m=1)
    endn = 16
  ENDIF
  IF(m=2)
    endn = 17
  ENDIF
  IF(m=3)
    endn = 16
  ENDIF
   
  n=1
  WHILE (n<=endn)
    'define sum3 = sum3 + pow((saavg'm'2'n' - avgtotal),2)'
  n = n+1
  ENDWHILE **n

m = m+1
ENDWHILE **m

'define stddevtotal = sqrt(sum3/104)'

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

'stddevactive = sqrt(sum1/56)'
'stddevbreak = sqrt(sum2/49)' 



'avgwet = (avg31+avg32)/2'
counter = 1
'sum = 0'
WHILE (counter <=18)
  'define sum = sum + pow((saavg31'counter' - avgwet),2)'
  counter = counter +1
ENDWHILE
counter = 1
WHILE (counter <=16)
  'define sum = sum + pow((saavg32'counter' - avgwet),2)'
  counter = counter +1
ENDWHILE

'stddevwet = sqrt(sum/33)'

'avgdry = (avg21+avg22)/2'
counter = 1
'sum = 0'
WHILE (counter <=19)
  'define sum = sum + pow((saavg21'counter' - avgdry),2)'
  counter = counter +1
ENDWHILE
counter = 1
WHILE (counter <=17)
  'define sum = sum + pow((saavg22'counter' - avgdry),2)'
  counter = counter +1
ENDWHILE

'stddevdry = sqrt(sum/35)'

**Now display the final variables and append them to the file**
'set gxout fwrite'
'set fwrite -ap /elninost2/ppantina/Traj/compare/cfsr_spells.grd'
'set x 1 192'
'set y 1 94'
'd avgactive'
'd avgbreak'
'd stddevactive'
'd stddevbreak'
'd avgtotal'
'd stddevtotal'
'd avgwet'
'd avgdry'
'd stddevwet'
'd stddevdry'

'disable fwrite'
'close 1'

**This is here because if you display variables at the end of the run**
**to check for errors, you will 'write' to the file and get rogues**

'set gxout shaded'