;;;;; First read in the model and data information (monthly)

PRO test, md=md


sp=1
species='O3'
specmod='OVP14_O3'
factor=1

if md eq 2 then begin 

savefile='/gpfsm/dnb31/jliu13/CHEM_EVA/CCM_eva/sav/TOZsubv_bench_i329_19_gmi_free_c180_132lev_monthly.sav'

;save, md, data, lat1, lat2,filename=savefile
restore, savefile
mdlabel='bench132'
endif else if md eq 0 then begin 
savefile='/gpfsm/dnb31/jliu13/CHEM_EVA/CCM_eva/sav/TOZsubv_bench_i329_gmi_free_c180_monthly.sav'

;save, md, data, lat1, lat2,filename=savefile
restore, savefile
mdlabel='bench'
ystar=2005
yend=2013
nyrs=yend-ystar+1
nmons=12
endif else if md eq 1 then begin 
mdname='MERRA2_GMI'
mdlabel='MR2GMI'

savefile='/gpfsm/dnb31/jliu13/CHEM_EVA/CCM_eva/sav/TOZsubv_'+mdname+'_monthly.sav'

;save, md, data, lat1, lat2,filename=savefile
restore, savefile
;mdlabel='m2g'


ystar=2005
yend=2017
nyrs=yend-ystar+1
nmons=12
endif 





open_device, olddevice, /ps, filename='/discover/nobackup/jliu13/CHEM_EVA/CCM_eva/FIGS/FIG4_TOZ_'+mdlabel+'_tsplot_regions.ps', $
    Bits=8, /color, /portrait

!X.OMARGIN=[4,9] 
!Y.OMARGIN=[4,4]
!X.MARGIN=[3,3]
!Y.MARGIN=[3,3]
!P.CHARTHICK=1.7
!P.Charsize  = 1.7
!X.THICK=4
!Y.THICK=4
!P.font = 1.1
;!P.Multi = [0, 2, 1, 0 ,0]
multipanel,rows=3,cols=2



time=findgen(nyrs*nmons)/12.+ystar
blue=GetColor('blue', 102)

Datename2 = ['05','06','07','08','09', '10','11', '12','13']
xticknamev=time[[1, 13, 25, 37, 49, 61, 73,85,97,109]-1] 
y1=[500, 500, 400, 500, 600, 400, 320]
y0=[100, 200, 200, 200, 200, 200, 250]
nrs=n_elements(lat1)

for r=0, nrs-1 do begin 

obsp=reform(data[r,*]);*factor
 mdp=reform(md[r,*]);*factor
 
indo=where (Finite(obsp) EQ 1, nindo) 
  indm=where (Finite(mdp) EQ 1) 
resulto = REGRESS(reform(time[indo]), obsp[indo], SIGMA=sigma, CONST=const, $
   MEASURE_ERRORS=measure_errors)

consto=const
sigmao=sigma
slopeo=resulto
tyo=resulto[0]*time+consto[0]
resultm = REGRESS(reform(time[indm]), mdp[indm], SIGMA=sigma, CONST=const, $
   MEASURE_ERRORS=measure_errors)
constm=const
sigmam=sigma
slopem=resultm
tym=resultm[0]*time+constm[0]

corr=correlate(obsp[indo], mdp[indo])
meano=mean(obsp[indo])
meanm=mean(mdp[indm])

xrange=[ystar, yend+1]
yrange=[0, 500]
regions=['90S-60S', '60S-30S', '30S-30N', '30N-60N', '60N-90N', '90S-90N', '60S-60N']
; Draw the line plot with no data
   cgPlot, time[indo], obsp[indo], Title=regions[r], XTitle=xtitle, YTitle=specplotname, $
      Position=position, /NoData,xrange=xrange,Charsize=2, YRange=[y0[r], y1[r]], YStyle=1;,xticks=12,xtickv=xticknamev, xtickname=datename2
time2=time[indo]
 ;  ; Fill in the error estimates.
 
   cgPlots, time[indm], mdp[indm], Color='red', Thick=6, linestyle=0;,EQN_SAMPLES=100;, PSym=-16, SymColor='blue', $
cgPlots, time[indo], obsp[indo], Color='black',Thick=6, linestyle=0;,EQN_SAMPLES=100;, PSym=-16, SymColor='black', $
   cgPlots, time[indm], tym[indm], Color='red', Thick=3, linestyle=2;,EQN_SAMPLES=100;, PSym=-16, SymColor='blue', $
cgPlots, time[indo], tyo[indo], Color='black',Thick=3, linestyle=2;,EQN_SAMPLES=100;, PSym=-16, SymColor='black', $

label ='slope: '+ strtrim(string(slopeo,'(f8.1)'), 2)+cgsymbol('+-')+strtrim(string(sigmao,'(f7.1)'), 2)
    xyouts, !x.window[0]+0.02, !y.window[1]-0.02, label,   $
          charsize =1,color=1, alignment=0,/normal

label =strtrim(string(slopem,'(f8.1)'), 2)+cgsymbol('+-')+strtrim(string(sigmam,'(f7.1)'), 2)
    xyouts, !x.window[0]+0.18, !y.window[1]-0.02, label,   $
          charsize =1,color=2, alignment=0,/normal
bias=(meanm-meano)/meano*100

label ='r!S!E2!N = '+strtrim(string(corr*corr,'(f8.1)'), 2)+', mean ( '+strtrim(string(meano,'(f8.1)'),2)+', ';+strtrim(string(meanm,'(f8.3)'),2)+', '+string(bias,'(f8.1)')+'%)'
    xyouts, !x.window[0]+0.02, !y.window[1]-0.04, label,   $
          charsize =1,color=1, alignment=0,/normal

label = strtrim(string(meanm,'(f8.1)'),2)+', '
    xyouts, !x.window[0]+0.19, !y.window[1]-0.04, label,   $
          charsize =1,color=2, alignment=0,/normal

label = string(bias,'(f8.1)')+'%)'
    xyouts, !x.window[0]+0.23, !y.window[1]-0.04, label,   $
          charsize =1,color=1, alignment=0,/normal

     
Legend, line=[0], color=[1], nlines=6, lcolor=[1], textcolor=1, label = ['TOZ sbuv'], thick=3,  $
     position= [0.87, 0.90, 0.88, 0.91]
  legend, line=[0], color=[2], textcolor=2, lcolor=[2],label=['MOD'],/ADD, thick=3

endfor   
  

 
   ;cgPS_Close 
close_device   , /TIMESTAMP
set_plot,olddevice
close,/all
print, 'DONE'

stop

end