**This script is the main script which is called by images_years_write.gs**
**It calculates the moisture source for 3 sets of neutral/dry/wet years, **
**as well as the 5 cases for each one of those sets... these variables **
**are then passed to the other script and written out, but this script**
**can be run solo and the variables displayed from here**

'reinit'

**The 'rain' counter is a counter for the 3 sets of**
**wet, dry, or neutral years... 'rain1' replaces 'rain'**
**with a the word... this is used for the average variables**

rain = 1
WHILE (rain <= 3)

**Conditional rain years**
  IF (rain = 1)
    rain1 = Neutral
    byear1 = 1984
    byear2 = 1985
    byear3 = 1989
    byear4 = 1996
    byear5 = 1997
  ENDIF
  IF (rain = 2)
    rain1 = Dry
    byear1 = 1979
    byear2 = 1982
    byear3 = 1986
    byear4 = 1987
    byear5 = 2002
  ENDIF
  IF (rain = 3)
    rain1 = Wet
    byear1 = 1983
    byear2 = 1988
    byear3 = 1990
    byear4 = 1994
    byear5 = 1998
  ENDIF

**The 'counter' keeps track of the 5 years in each given set**
**of wet/dry/neutral years**
  counter = 1
  WHILE (counter <= 5)

**Set up the correct years for the file sources**
    IF (counter = 1)
     byear = byear1
    ENDIF
    IF (counter = 2)
      byear = byear2
    ENDIF
    IF (counter = 3)
      byear = byear3
    ENDIF
    IF (counter = 4)
      byear = byear4
    ENDIF
    IF (counter = 5)
     byear = byear5
    ENDIF

*******************************************
**FAMILIAR SUBROUTINE**********************
*******************************************
**This file is constantly rewritten as we loop**
**thru all the different source files... it acts**
**as a template for the first few lines of source.ctl**
    tmpoutfile = 'test.txt'

**Get the information regarding pentad starting**
**and ending times... may need to change this in**
**the future, espeically if we change periods**
    periodfile = read(wcr.in)
    periodline = sublin(periodfile,2)
    period0    = subwrd(periodline,2)
    period1    = subwrd(periodline,3)
    err        = close(wcr.in)

**First read the number of x-y files that are needed**
**There should be 41 over India... also find the number**
**and dimensions of the corners to be plotted...**
**Remember that gridarea is user-defined and will need**
**to change if you change the plot area**
    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)
 
    t0 = period0
    tf = period1

    say 'Pentads used in tracing: 'period0'-'period1
    say 'Number of files/points = 'filenum
    say 'Number of corners      = 'cornum

**Now loop thru the different source files**
**by reading the gridarea file... this loop**
**reads each text line, and finds the correct**
**source by its x and y coordinates**
    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  = '***'

**Remove the original source file and replace it with**
**one calling the next source.grd file... in this way**
**we find s from each .grd file, store it, and sum it**
**up as saavg**  
      '!/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'

**Open the newly created source file, and add s from it**
**to previous total of s from other .grd files**
      say 'Create Control File: 'file
      'open source.ctl'
     
      IF (xycount = 1)
        'define saavg'rain''counter' = sum(s, t='t0', t='tf')*5'
      ELSE
        'define saavg'rain''counter' = saavg'rain''counter' + sum(s, t='t0', t='tf')*5' 
      ENDIF

      'close 1'
*******************************************
**END FAMILIAR SUBROUTINE******************
*******************************************

**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'

**Now display the individual years... these will appear to 'grow'**
**if grads is open, as each gridpoint is added to the previous**
**to create the sum total of moisture source**
     'c'
     'set gxout shaded'
     'set clevs 30 60 90 120 150 180 210 240 270'
     'set ccols  0 25 26 16 17  18  19  20  21  22  23 24'
     'd saavg'rain''counter

      xycount = xycount+1
    ENDWHILE **xycount

    err = close(gridarea)
    counter = counter + 1
  ENDWHILE **counter

**Here we calculate an average of a given set of years**
  'open source.ctl'
  'define avg'rain1' = (saavg'rain'1+saavg'rain'2+saavg'rain'3+saavg'rain'4+saavg'rain'5)/5'
  'close 1'

 rain = rain + 1

ENDWHILE **rain