* To be substituted by a master script
domain = 'DOMAIN'
deltat = 24
say ' '
say 'Processing 'deltat'-h precip for domain 'domain

*East Africa domain
if (domain = 'eaf')
  lats = "-14.89 14.89"
  lons = "22.98 52.02"
  labela = "2.7 1.4"
  labelb = "2.7 1.0"
  llint = 2
endif
*Kenya domain
if (domain = 'ken')
  lats = "-10.61 8.37"
  lons = "28.13 47.20"
  labela = "2.7 1.4"
  labelb = "2.7 1.0"
  llint = 2
endif

* time is the grads valid time after time = 1 (TBD in script)
* for the period ending the 24-h accumulation period
* ndays = # of days past grads start time
startday = 1
endday = 3

day=startday
while (day < endday)

'environ_shaded.ch'
'c'
'open imerge.ctl'
'open imergl.ctl'
'set lat ' lats
'set lon ' lons

*******************
* IMERG plots loop
*******************
time = day*24*2 + 12*2
endtime = time + 1
tback = time - deltat*2

while (time < endtime)
  'set t ' tback
  cdateback = mydate()
  chrback = subwrd(cdateback,1)
  cdyback = subwrd(cdateback,2)
  cmoback = subwrd(cdateback,3)
  cyrback = subwrd(cdateback,4)
  remain = math_fmod(tback,2)
  if (remain = 0)
    hhback = substr(chrback,1,2)
    cmnback = ":00"
    minback = "00"
  else
    hhback = chrback
    cmnback = ""
    minback = "30"
  endif
  cmmback = cmo2cmm(cmoback)
  say 'IMERG tback = 'cdateback ' hhmm = ' hhback minback

  'set t ' time
  cdate = mydate()
  chr = subwrd(cdate,1)
  cdy = subwrd(cdate,2)
  cmo = subwrd(cdate,3)
  cyr = subwrd(cdate,4)
  remain = math_fmod(time,2)
  if (remain = 0)
    hh = substr(chr,1,2)
    cmn = ":00"
    min = "00"
  else
    hh = chr
    cmn = ""
    min = "30"
  endif
  cmm = cmo2cmm(cmo)
  say 'Doing IMERG 'deltat'-h plots for time = 'cdate ' hhmm = 'hh min

********************
*Early NRT IMERG plot
********************
  'c'
  say 'Early'
  'set grads off'
* J.Case (6/29/2015) -- trying to display multi-panel images in one graphics file.
*  if (domain = 'pr')
** Center-left
*    'set vpage 0 5.5 2.12 6.37'
*  else
** Upper-left
*    'set vpage 0 5.5 4.25 8.5'
*  endif
  'environ_shaded.ch'
  'set grads off'
  'set xlint ' llint
  'set ylint ' llint
  'set xlopts 1 6 0.16'
  'set ylopts 1 6 0.16'
* mm
  'set clevs  1  2  5 10 15 20 25 50 75 100 125 150 200 250 300'
  'set ccols 0 33 35 37 39 43 45 47 49 53 55  57  59  65  67  69'
  'd sum(prate.1,t='tback',t='time-1',1)/2'
* AHPS inches scale and approximate color RGBs
*  'set clevs  0.01 0.10 0.25 0.50 0.75 1.0 1.5 2.0 2.5 3.0 4.0 5.0 6.0 8.0 10.0'
*  'set ccols 0   43   46    4   33   36  39  22  23  24  25  27  29  64   55   0'
* J.Case (7/30/2015) -- bug fix. Summing over one too many times (24.5 hours)
* Since IMERG is valid from t to t + 30 min, sum from tback to time-1
**  'd sum(prate.1,t='tback',t='time',1)/50.8'
*  'd sum(prate.1,t='tback',t='time-1',1)/50.8'

* draw Kenya county shapefile
  'set line 1 1 3'
  'draw shp Kenya_admin_2014_WGS84 1 46'

  'run cbar.gs'
  'set string 1 c 10 0'
  'set strsiz 0.19'
* mm accumulation
  'draw string 5.5 10.75 Early IMERG 'deltat'-h Accumulated Precipitation (mm) ending: '
  'set strsiz 0.14'
  'draw string 5.5 10.45 Valid for period ending 'chr'z 'cdy' 'cmo' 'cyr
* in accumulation
*  'draw string 5.5 8.1 Early IMERG 'deltat'-h prec (in) ending: 'chr''cmn'z 'cdy' 'cmo' 'cyr' '

* J.Case (7/29/2015) -- add label for # of valid times going into accumulation
  'set string 89 c 12 0'
  'set strsiz 0.22'
  if (day = 1)
    'draw string 'labela' NUMEA of 48'
  else
    'draw string 'labela' NUMEB of 48'
  endif
  'draw string 'labelb' times'
  'set string 1 c 9 0'
  'set strsiz 0.22'
  if (day = 1)
    'draw string 'labela' NUMEA of 48'
  else
    'draw string 'labela' NUMEB of 48'
  endif
  'draw string 'labelb' times'

  'gxprint imerge24h_'domain'_'cyr cmm cdy'_'hh min'Z.png png x1000 y800 white'

********************
*Late NRT IMERG plot
********************
  'c'
  say 'Late'
  'set grads off'
* J.Case (6/29/2015) -- trying to display multi-panel images in one graphics file.
*  if (domain = 'pr')
** Center-right
*    'set vpage 5.5 11.0 2.12 6.37'
*  else
** Upper-right
*    'set vpage 5.5 11.0 4.25 8.5'
*  endif
  'environ_shaded.ch'
  'set grads off'
  'set xlint ' llint
  'set ylint ' llint
  'set xlopts 1 6 0.16'
  'set ylopts 1 6 0.16'
* mm
  'set clevs  1  2  5 10 15 20 25 50 75 100 125 150 200 250 300'
  'set ccols 0 33 35 37 39 43 45 47 49 53 55  57  59  65  67  69'
  'd sum(prate.2,t='tback',t='time-1',1)/2'
* AHPS inches scale and approximate color RGBs
*  'set clevs  0.01 0.10 0.25 0.50 0.75 1.0 1.5 2.0 2.5 3.0 4.0 5.0 6.0 8.0 10.0'
*  'set ccols 0   43   46    4   33   36  39  22  23  24  25  27  29  64   55   0'
* J.Case (7/30/2015) -- bug fix. Summing over one too many times (24.5 hours)
* Since IMERG is valid from t to t + 30 min, sum from tback to time-1
**  'd sum(prate.2,t='tback',t='time',1)/50.8'
*  'd sum(prate.2,t='tback',t='time-1',1)/50.8'

* draw Kenya county shapefile
  'set line 1 1 3'
  'draw shp Kenya_admin_2014_WGS84 1 46'

  'run cbar.gs'
  'set string 1 c 10 0'
  'set strsiz 0.19'
* mm accumulation
  'draw string 5.5 10.75 Late IMERG 'deltat'-h Accumulated Precipitation (mm) ending: '
  'set strsiz 0.14'
  'draw string 5.5 10.45 Valid for period ending 'chr'z 'cdy' 'cmo' 'cyr
* in accumulation
*  'draw string 5.5 8.1 Late IMERG 'deltat'-h prec (in) ending: 'chr''cmn'z 'cdy' 'cmo' 'cyr' '

* J.Case (7/29/2015) -- add label for # of valid times going into accumulation
  'set string 89 c 12 0'
  'set strsiz 0.22'
  if (day = 1)
    'draw string 'labela' NUMLA of 48'
  else
    'draw string 'labela' NUMLB of 48'
  endif
  'draw string 'labelb' times'
  'set string 1 c 9 0'
  'set strsiz 0.22'
  if (day = 1)
    'draw string 'labela' NUMLA of 48'
  else
    'draw string 'labela' NUMLB of 48'
  endif
  'draw string 'labelb' times'

  'gxprint imergl24h_'domain'_'cyr cmm cdy'_'hh min'Z.png png x1000 y800 white'

  time = time + deltat*2
endwhile
'close 2'
'close 1'
*******************
* end IMERG loop
*******************

day = day + 1

endwhile
***************
* end day loop
***************

'quit'

************************************
function mydate
'query time'
sres = subwrd(result,3)
i=1
while (substr(sres,i,1)!='Z')
i=i+1
endwhile
hour = substr(sres,1,i-1)
isav = i
i = i + 1
while (substr(sres,i,1)>='0' & substr(sres,i,1)<='9')
i = i + 1
endwhile
day = substr(sres,isav+1,i-isav-1)
month = substr(sres,i,3)
year = substr(sres,i+3,4)
return (hour' 'day' 'month' 'year)
************************************

************************************
function cmo2cmm(cmo)
if (cmo = "JAN") ; cmm = "01" ; endif
if (cmo = "FEB") ; cmm = "02" ; endif
if (cmo = "MAR") ; cmm = "03" ; endif
if (cmo = "APR") ; cmm = "04" ; endif
if (cmo = "MAY") ; cmm = "05" ; endif
if (cmo = "JUN") ; cmm = "06" ; endif
if (cmo = "JUL") ; cmm = "07" ; endif
if (cmo = "AUG") ; cmm = "08" ; endif
if (cmo = "SEP") ; cmm = "09" ; endif
if (cmo = "OCT") ; cmm = "10" ; endif
if (cmo = "NOV") ; cmm = "11" ; endif
if (cmo = "DEC") ; cmm = "12" ; endif
return cmm
************************************

