*
*	This program is modified from scm_heating_comp_2d.gs
* 
*
*	Make 48h time series plots 
*

function main(args) 

'reinit'
say "scm_heating_comp_1d.gs " args 

*	Get arguments 

if(args='?') 
  say 'scm_heating_comp_1d.gs requires 3 arguments: experiment, site, and day'
  return 
else 
  experiment=subwrd(args,1); if(experiment=''); experiment="";endif 
  site=subwrd(args,2); if(site=''); site=""; endif 
  ssdate=subwrd(args,3); if(ssdate=''); ssdate=""; endif 
  say 'experiment :'experiment
  say 'site :'site
  say 'day :'ssdate
endif 

*'reinit'

compare=yes
*compare=no

* 	examples of experiment and sites 
*experiment=armgfs
*site='FL2'
*site='imet1'
*site='sgpC1'
*site='Manus'
*site='togacoare'
*site='WHOI2'
*site=EPIC6
*site=Peace

* 	open files 

'open /ptmp/wx22rs/scm_'experiment'_'site'_new_20091215/flx.ctl'
say 'open /ptmp/wx22rs/scm_'experiment'_'site'_new_20091215/flx.ctl'
if(compare=yes)
  'open /ptmp/wx22rs/ncepscm/'experiment'/'site'/flx.ctl'
   say 'open /ptmp/wx22rs/ncepscm/'experiment'/'site'/flx.ctl'
endif 

'set display color white'
'clear'

* 	experiment titles 
fhour=48 
finthr=3
nhour=fhour/finthr +1
titexp=GFS_SCM
titexp1=GFS_SCMP
titexp2=GFS_SCMN
expname='SCM_'experiment'_'site
ssdate=ssdate'_fcst'

'run /global/save/wx23ss/gscripts/rgbset.gs'
'set xlopts 1 4 0.15'
'set ylopts 1 4 0.15'
'set strsiz 0.4'

varsx=1
while(varsx<=53)

if (varsx=1)
varname=dusfc
tvarname='dusfc'
varnamegfs=ustress
tvarnamegfs='ustress'
sunit=''
aclevs='0.00001 0.00002 3.e-5 4.e-5 5.e-5 6.e-5 7.e-5 8.e-5 9.e-5 1.e-4'
endif 

if (varsx=2)
varname=dvsfc
tvarname='dvsfc'
varnamegfs=vstress
tvarnamegfs='vstress'
sunit=''
aclevs='-0.0002 -0.0001 -0.00007 -0.00005 -0.00003 -0.00001 0.00001 0.00003 0.00005'
endif

if (varsx=3)
varname=dtsfc
tvarname='dtsfc'
varnamegfs=shf
tvarnamegfs='shf'
sunit=''
aclevs='-5e-4 -4e-4 -3e-4 -2.0e-4 -1.0e-4 -5e-5 -2.5e-5 2.5e-5 5e-5 1.0e-4 2.0e-4 3e-4 4e-4 5e-4'
endif

if (varsx=4)
varname=dqsfc
tvarname='dqsfc'
varnamegfs=lhf
tvarnamegfs='lhf'
sunit=''
aclevs='-5e-4 -4e-4 -3e-4 -2.0e-4 -1.0e-4 -5e-5 -2.5e-5 2.5e-5 5e-5 1.0e-4 2.0e-4 3e-4 4e-4 5e-4'
endif

if (varsx=5)
varname=tsea
tvarname='Surface temperature'
sunit='(K)'
aclevs='-5e-7 -4e-7 -3e-7 -2.0e-7 -1.0e-7 -5e-8 -2.5e-8 2.5e-8 5e-8 1.0e-7 2.0e-7 3e-7 4e-7 5e-7'
endif

if (varsx=6)
varname=smc1
tvarname='Soil moisture content'
sunit='(volume fraction)'
aclevs='-5e-5 -4e-5 -3e-5 -2.0e-5 -1.0e-5 -5e-6 5e-6 1.0e-5 2.0e-5 3e-5 4e-5 5e-5'
endif

if (varsx=7)
varname=smcr2
tvarname='smcr2'
sunit='(?)'
aclevs='-5e-7 -4e-7 -3e-7 -2.0e-7 -1.0e-7 -5e-8 5e-8 1.0e-7 2.0e-7 3e-7 4e-7 5e-7'
endif

if (varsx=8)
varname=stc1
tvarname='soil temperature (K)'
sunit='(?)'
aclevs='-5e-7 -4e-7 -3e-7 -2.0e-7 -1.0e-7 -5e-8 5e-8 1.0e-7 2.0e-7 3e-7 4e-7 5e-7'
endif

if (varsx=9)
varname=stc2
tvarname='stc2'
sunit='(?)'
aclevs='-5e-7 -4e-7 -3e-7 -2.0e-7 -1.0e-7 -5e-8 -1e-8 1e-8 5e-8 1.0e-7 2.0e-7 3e-7 4e-7 5e-7'
endif

if (varsx=10)
varname=weasd
tvarname='water equivalent of accum snow depth'
sunit='(m)'
aclevs='0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0'
endif

if (varsx=11)
varname=dlwsfc
tvarname='Sfc downward long-wave heat flux'
varnamegfs= dlws
tvarnamegfs='dlws'
sunit='(W/m^2)'
aclevs='200 210 220 230 240 250 260 270 280 290 300 310'
endif

if (varsx=12)
varname=ulwsfc
tvarname='Sfc upward long-wave heat flux'
varnamegfs= ulws
tvarnamegfs='ulws'
sunit='(W/m^2)'
aclevs='-0.002 0.002 0.004 0.006 0.008 0.010 0.012 0.014 0.016 0.018 0.020'
endif

if (varsx=13)
varname=olr
tvarname='Outgoing long-wave radiation flux'
varnamegfs=ulwt
tvarnamegfs='ulwt'
sunit='(W/m^2)'
aclevs='-0.0002 -0.00015 -0.0001 -0.00005 0.0005 0.0001 0.00015 0.0002'
endif

if (varsx=14)
varname=uswrftoa
tvarname='Upward short-wave radiation flux at TOA'
varnamegfs=uswt
tvarnamegfs='uswt'
sunit='(W/m^2)'
aclevs='-0.0002 -0.00015 -0.0001 -0.00005 0.0005 0.0001 0.00015 0.0002'
endif

if (varsx=15)
varname=uswrfsfc
tvarname='Upward short wave radiation flux at sfc'
varnamegfs=usws
tvarnamegfs='usws'
sunit='(W/m^2)'
endif

if (varsx=16)
varname=dswrfsfc
tvarname='Downward short wave radiation flux at sfc'
varnamegfs=dsws
tvarnamegfs='dsws'
sunit='(W/m^2)'
endif

if (varsx=17)
varname=tcdchcl
tvarname='high cloud cover'
varnamegfs=cldh
tvarnamegfs='cldh'
sunit='(%)'
endif

if (varsx=18)
varname=preshcl
tvarname='Pressure at high cloud top'
sunit='(Pa)'
endif

if (varsx=19)
varname=preshbl
tvarname='Pressure at high cloud base'
sunit='(Pa)'
endif

if (varsx=20)
varname=tmphcl
tvarname='Temperature of high clouds'
sunit='(K)'
endif

if (varsx=21)
varname=tcdcmcl
tvarname='Middle cloud cover'
varnamegfs=cldm
tvarnamegfs='cldm'
sunit='(%)'
endif

if (varsx=22)
varname=presmcl
tvarname='Pressure at middle cloud top'
sunit='(Pa)'
endif

if (varsx=23)
varname=presmbl
tvarname='Pressure at middle cloud base'
sunit='(Pa)'
endif

if (varsx=24)
varname=tmpmcl
tvarname='Temperature of middle clouds'
sunit='(K)'
endif

if (varsx=25)
varname=tcdclcl
tvarname='Low cloud cover'
varnamegfs=cldl
tvarnamegfs='cldl'
sunit='(%)'
endif

if (varsx=26)
varname=preslcl
tvarname='Pressure at low cloud top'
sunit='(Pa)'
endif

if (varsx=27)
varname=preslbl
tvarname='Pressure at low cloud base'
sunit='(Pa)'
endif

if (varsx=28)
varname=tmplcl
tvarname='Temperature of low clouds'
sunit='(K)'
endif

if (varsx=29)
varname=prate
tvarname='Precipitation rate'
varnamegfs=precip
tvarnamegfs='precip'
sunit='(mm/day)'
endif

if (varsx=30)
varname=cprate
tvarname='Convective precipitation rate'
varnamegfs=cnvprc
tvarnamegfs='cnvprc'
sunit='(mm/day)'
endif

if (varsx=31)
varname=gflux
tvarname='Conductive heating flux'
varnamegfs=gflx
tvarnamegfs='gflx'
sunit='(W/m^2)'
endif

if (varsx=32)
varname=slmsk
tvarname='Landmask'
varnamegfs=lsmask
tvarnamegfs='lsmask'
sunit=''
endif

if (varsx=33)
varname=icec
tvarname='Ice fraction'
sunit='(%)'
endif

if (varsx=34)
varname=u10m
tvarname='u10m'
varnamegfs=u10
tvarnamegfs='u10'
sunit='(m/s)'
endif

if (varsx=35)
varname=v10m
tvarname='v10m'
varnamegfs=v10
tvarnamegfs='v10'
sunit='(m/s)'
endif

if (varsx=36)
varname=t2m
tvarname='Temperature at 2m AGL'
varnamegfs=t2
tvarnamegfs='t2'
sunit='(K)'
endif

if (varsx=37)
varname=q2m
tvarname='Specific humidity at 2m AGL'
varnamegfs=q2
tvarnamegfs='q2'
sunit='(Kg/Kg)'
endif

if (varsx=38)
varname=psfc
tvarname='Surface pressure'
varnamegfs=psfc1
tvarnamegfs='psfc1'
sunit='(Pa)'
endif

if (varsx=39)
varname=tmpmax
tvarname='Maximum temperature '
varnamegfs=tmax
tvarnamegfs='tmax'
sunit='(K)'
endif

if (varsx=40)
varname=tmpmin
tvarname='Minimum temperature '
varnamegfs=tmin
tvarnamegfs='tmin'
sunit='(K)'
endif

if (varsx=41)
varname=runoff
tvarname='Water runoff'
sunit='(m/s)'
endif

if (varsx=42)
varname=ep
tvarname='Potential eavporation'
sunit=''
endif

if (varsx=43)
varname=dugwd
tvarname='Zonal gravity wave stress'
varnamegfs=ugrw
tvarnamegfs='ugrw'
sunit='(N/m^2)'
endif

if (varsx=44)
varname=dvgwd
tvarname='Meridional gravity wave stress'
varnamegfs=vgrw
tvarnamegfs='vgrw'
sunit='(N/m^2)'
endif

if (varsx=45)
varname=hpbl
tvarname='pbl height'
varnamegfs=hpbl
tvarnamegfs='hpbl'
sunit='(m)'
endif

if (varsx=46)
varname=pwat
tvarname='Precipitable water'
sunit=''
endif

if (varsx=47)
varname=albedo
tvarname='Surface albedo'
sunit=''
endif

if (varsx=48)
varname=tcdcclm
tvarname='Total cloud cover'
sunit='(%)'
endif

if (varsx=49)
varname=tcdccnv
tvarname='Totoal convective cloud cover'
sunit='(%)'
endif

if (varsx=50)
varname=cvt
tvarname='pressure of convective cloud top'
sunit='(Pa)'
endif

if (varsx=51)
varname=cvb
tvarname='pressure of convective cloud base'
sunit='(Pa)'
endif

if (varsx=52)
varname=tcdcblc
tvarname='Cloud cover boundary layer'
sunit='(%)'
endif

if (varsx=53)
varname=dswrftoa
tvarname='Downward short wave radiation flux at TOA'
sunit='(W/m^2)'
endif

* 	define filenames 

giffile='./figures/'expname'.'varname'_'ssdate'_1d.gif'
grfile='./figures/'expname'.'varname'_'ssdate'_1d.gr'

*	plot results from the current version of SCM 

'set vpage 0 11 0 8.5'
'set dfile 1'
'set grads off'
'set t 1 'nhour
*'set gxout contour'
if(varname=prate | varname=cprate)
  'd 'varname'.1*86400'
  say 'prate | cprate'
else
  if(varname=tcdccnv|varname=tcdcclm|varname=tcdchcl|varname=tcdcmcl|varname=tcdclcl|varname=tcdcblc)
     'set vrange 0 105'
  endif
 'd 'varname'.1'
endif 
*'draw xlab fcst hour'
*'draw ylab sigma level'
'set strsiz 0.15 0.15'
'set string 1'
'draw string 3.5 8.25 'tvarname' 'sunit

*
*	need to compare with the previous version of SCM: RSUN, Jul 18, 2008
*
if(compare=yes)
compflag=0

if(varsx<5 ) 
compflag = 1 
endif  
if(varsx>10 & varsx <18)
compflag = 1
endif 
if(varsx=21 | varsx=25 )
compflag = 1
endif
if(varsx>28 & varsx <33)
compflag = 1
endif
if(varsx>33 & varsx <41)
compflag = 1
endif
if(varsx>42 & varsx <46)
compflag = 1
endif

if(compflag=1) 
*	make the plots using GFS output 

*'set dfile 2'
*'set grads off'
'set line 4'
if(varname=prate | varname=cprate)
  'd 'varnamegfs'.2*86400'
  say 'prate | cprate'
else
  if(varname=tcdccnv|varname=tcdcclm|varname=tcdchcl|varname=tcdcmcl|varname=tcdclcl|varname=tcdcblc)
     'set vrange 0 105'
  endif
 'd 'varnamegfs'.2'
endif 
*'draw xlab fcst hour'
*'draw ylab sigma level'
*'set strsiz 0.15 0.15'
*'set string 1'
*'draw string 3.5 8.25 'tvarname' 'sunit' (GFS20091215v)'

'printim 'giffile' gif x720 y500'
'enable print 'grfile
'print'
'disable print'
'clear'

*say 'type in c to continue or quit to exit'
*pull corquit
*corquit

* plot the difference in the results from the two versions of the SCM 

*       define filenames
clear 

giffile='./figures/'expname'.'varname'_'ssdate'_1d_diff.gif'
grfile='./figures/'expname'.'varname'_'ssdate'_1d_diff.gr'

'set vpage 0 11 0 8.5'
'set grads off'
'set t 1 'nhour
if(varname=prate | varname=cprate)
  'd ('varname'.1 - 'varnamegfs'.2)*86400'
  say 'prate | cprate'
else
  if(varname=tcdccnv|varname=tcdcclm|varname=tcdchcl|varname=tcdcmcl|varname=tcdclcl|varname=tcdcblc)
     'set vrange 0 105'
  endif
 'd 'varname'.1 - 'varnamegfs'.2'
endif
*'d 'varname'.1 - 'varnamegfs'.2'
'set strsiz 0.15 0.15'
'set string 1'
'draw string 3.5 8.25 diff 'tvarname' 'sunit

'printim 'giffile' gif x720 y500'
'enable print 'grfile
'print'
'disable print'
'clear'

* convert gr to ps 
*'gxps -c -i wgrdprs_0.trigger75.son.gr -o wgrdprs_0.trigger75.son.ps'

*say 'type in c to continue or quit to exit'
*pull corquit
*corquit

endif 
endif 

varsx=varsx+1
endwhile

*	close files 

if (compare=yes)
  'close 2' 
endif 
'close 1'