'open landdat.pgrbb.ctl'
'define land=LANDsfc'
'close 1'

'set display color white'
'clear'
'set grads off'
*'open /gpfs/hps3/emc/ens/save/Hong.Guan/plot_coupexp/plume/fluxes/dswrf/CERESLa.grb_1p0.ctl'
*'open /gpfs/hps3/emc/ens/save/Hong.Guan/plot_coupexp/plume/fluxes/dswrf/CERES.grb_1p0.ctl'

'open /lfs/h2/emc/ens/save/hong.guan/plot_coupexp/plume/fluxes/dswrf/CERES.grb_1p0.ctl'
'q file'
'set t 32 300'
'define obsdsw=PRMSLsfc'
'define obsdlw=PVORTsfc'
'define obsulw=PTENDsfc'
'define obsusw=PRESsfc'
'close 1'

*'open /lfs/h2/emc/ens/noscrub/hong.guan/ERA5/fluxes/hourly/grib1_dat_2015/LHTFL_dat.ctl'
*'open /lfs/h2/emc/ens/noscrub/hong.guan/ERA5/fluxes/hourly/grib1_dat_2015/LHTFL_dat.ctl'
'open /lfs/h2/emc/ens/noscrub/hong.guan/lhsh.dat/SLHF_1p00_1516-daily.ctl'
'q file'
'set t 1 269'
'define obslhfl=lh/3600.'
'close 1'

*'open /lfs/h2/emc/ens/noscrub/hong.guan/ERA5/fluxes/hourly/grib1_dat_2015/SHTFL_dat.ctl'
'open /lfs/h2/emc/ens/noscrub/hong.guan/lhsh.dat/SSHF_1p00_1516-daily.ctl'
'q file'
'set t 1 269'
'define obsshfl=sh/3600.'
'define obs=obsdlw+obsdsw-obsulw-obsusw+obsshfl+obslhfl'
'close 1'


*'open /gpfs/hps3/emc/ens/save/Hong.Guan/plot_coupexp/plume/fluxes/dswrf/Bing_9months_LaNina/dswrf_sfc_GEFS_01jun2017_00z_d01_d269_mean.grb2.ctl'
*'open /gpfs/hps3/emc/ens/save/Hong.Guan/GEFSv12_daily_dell_wei/seasngefs.20150601/extract_dell/output_p2/dswrf_sfc_GEFS_01jun2015_00z_d01_d269_mean.grb2.ctl'
*'open /lfs/h2/emc/ens/noscrub/Hong.Guan/Bing_9months/gefs.20150601/00/dswrf_sfc_GEFS_01jun2015_00z_d01_d269_mean.grb2.ctl'
'open /lfs/h2/emc/ens/noscrub/hong.guan/cpls_avg/DSWRFsfc-daily.data.ctl'
'q file'
'set t 1 269'
'define meandsw=DSWRF'
'close 1'
*'open /lfs/h2/emc/ens/noscrub/Hong.Guan/Bing_9months/gefs.20150601/00/dlwrf_sfc_GEFS_01jun2015_00z_d01_d269_mean.grb2.ctl'
'open /lfs/h2/emc/ens/noscrub/hong.guan/cpls_avg/DLWRFsfc-daily.data.ctl'
'q file'
'set t 1 269'
'define meandlw=DLWRF'
'close 1'
'open /lfs/h2/emc/ens/noscrub/hong.guan/cpls_avg/USWRFsfc-daily.data.ctl'
*'open /lfs/h2/emc/ens/noscrub/Hong.Guan/Bing_9months/gefs.20150601/00/uswrf_sfc_GEFS_01jun2015_00z_d01_d269_mean.grb2.ctl'
'q file'
'set t 1 269'
'define meanusw=USWRF'
'close 1'
*'open /lfs/h2/emc/ens/noscrub/Hong.Guan/Bing_9months/gefs.20150601/00/ulwrf_sfc_GEFS_01jun2015_00z_d01_d269_mean.grb2.ctl'
'open /lfs/h2/emc/ens/noscrub/hong.guan/cpls_avg/ULWRFsfc-daily.data.ctl'
'q file'
'set t 1 269'
'define meanulw=ULWRF'
'close 1'
*'open /lfs/h2/emc/ens/noscrub/Hong.Guan/Bing_9months/gefs.20150601/00/lhfl_sfc_GEFS_01jun2015_00z_d01_d269_mean.grb2.ctl'
'open /lfs/h2/emc/ens/noscrub/hong.guan/cpls_avg/LHTFLsfc-daily.data.ctl'
'q file'
'set t 1 269'
'define meanlhfl=lhtfl'
'close 1'
*'open /lfs/h2/emc/ens/noscrub/Hong.Guan/Bing_9months/gefs.20150601/00/shfl_sfc_GEFS_01jun2015_00z_d01_d269_mean.grb2.ctl'
'open /lfs/h2/emc/ens/noscrub/hong.guan/cpls_avg/SHTFLsfc-daily.data.ctl'
'q file'
'set t 1 269'
'define meanshfl=shtfl'
'define fcst=meandlw+meandsw-meanusw-meanulw-meanlhfl-meanshfl'

*'open /lfs/h2/emc/ens/save/Hong.Guan/plot_coupexp/plume/fluxes/dswrf/dswrf_sfc_GEFS_01jun2015_00z_d01_d269_spr.grb2.ctl'
*'q file'
'set gxout shaded'

'set rgb  16    0    0  255'
'set rgb  17   55   55  255'
'set rgb  18  110  110  255'
'set rgb  19  165  165  255'
'set rgb  20  220  220  255'
* These are the RED shades
'set rgb  21  255  220  220'
'set rgb  22  255  165  165'
'set rgb  23  255  110  110'
'set rgb  24  255   55   55'
'set rgb  25  255    0    0'

*'set vpage 0.0 4.3 5 10.5'
'set t 1 269'
'set y 1'
'set lon 150 280'
*'set lon 150 290'
'set grads off'
'set csmooth on'
'set yflip off'
*'set clevs   -100 -75 -50 -25 0 25 50 75 100'
'set clevs   -80 -60 -40 -20 0 20 40 60 80'
*'set clevs   -40 -30 -20 -10 0 10 20 30 40'
'set ccols 16 17 18 19 20 21 22 23 24 25'

'set xlopts 1 5 0.18'
'set xlint 20'
'set ylopts 1 5 0.18'
*'d ave(maskout((mean-obs+meandlw-obsdlw-meanulw+obsulw-meanusw+obsusw-meanlhfl+obslhfl-meanshfl+obslhfl),0.999999-land),lat=-5,lat=5))'
'd ave(maskout((fcst-obs),0.999999-land),lat=-5,lat=5))'
*'draw xlab longitude'
*'draw ylab time'
'cbarn'
*'draw title DSWRF bias, lat=5S-5N, init. 2015060100'
'draw title net flux bias, lat=5S-5N, init. 2015060100'
*'printim DSWRF_20150601_bias.png x450 y500'
'printim NET_20150601_bias.png x450 y500'
'quit'

*'define point=aave(ave(maskout(mean,0.999999-land),t-3,t+3),lon=270,lon=290,lat=-20,lat=0)'
'define point=aave(ave(maskout(mean,0.999999-land),t-3,t+3),lon=220,lon=260,lat=-5,lat=5)'
'set gxout print'
'set prnopts %6.2f 1 1'
'd point'
*dummy=write('DLWRF_NINO3.4_mean.txt',result,append)
*dummy=write('DSWRF_SAcoast_mean.txt',result,append)
dummy=write('DSWRF_LON220260_mean.txt',result,append)

*'define point=aave(ave(maskout(mean-obs,0.999999-land),t-3,t+3),lon=270,lon=290,lat=-20,lat=0)'
'define point=aave(ave(maskout(mean-obs,0.999999-land),t-3,t+3),lon=220,lon=260,lat=-5,lat=5)'
'set gxout print'
'set prnopts %6.2f 1 1'
'd point'
*dummy=write('DLWRF_NINO3.4_mean.txt',result,append)
*dummy=write('DSWRF_SAcoast_meanbias.txt',result,append)
dummy=write('DSWRF_LON220260_meanbias.txt',result,append)

quit
*'set lon 220'
'define point=aave(ave(mean,t-3,t+3),lon=190,lon=240,lat=-5,lat=5)'
'set gxout print'
'set prnopts %6.2f 1 1'
'd point'
dummy=write('DSWRF_NINO3.4_mean.txt',result,append)

'quit'
'clear'
'set t 1 269'
'set y 1'
*'set lon 150 280'
'set lon 150 300'
'set grads off'
'set csmooth on'
'set yflip off'
'set clevs   -1.5 -1 -0.5 -0.25 0 0.25 0.5 1 1.5'
'set ccols 16 17 18 19 20 21 22 23 24 25'

'set xlopts 1 5 0.18'
'set xlint 20'
'set ylopts 1 5 0.18'
'd ave(maskout((mean-obs),0.999999-land),lat=-5,lat=5))'
'cbarn'
'draw title SST rmse, lat=5S-5N, init. 2015060100'
'print'
'printim SST_20150601_rmse.png x450 y500'

'clear'
'quit'

function digs(string,num)
  nc=0
  pt=""
  while(pt = "")
    nc=nc+1
    zzz=substr(string,nc,1)
    if(zzz = "."); break; endif
  endwhile
  end=nc+num
  str=substr(string,1,end)
return str