;;;;; 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 0 then begin mdname='bench_i329_gmi_free_c180' mdlabel='benchmark' ; omimls grid 1.25x1 -60-60 OModelInfo = CTM_type('generic', res=[1.25, 1], halfpolar=0, CENTER180=0) OGridInfo = ctm_grid(OModelInfo) OIMX = OGridInfo.IMX ; Maximum I (longitude) dimension OJMX = OGridInfo.JMX ; Maximum J (latitude) dimension) olonmid = OGridInfo.XMid ; Array of longitude centers olatmid = OGridInfo.YMid ; Array of latitude centers path='/gpfsm/dnb04/projects/p80/mmanyin/benchmarks/bench_i329_gmi_free_c180/inst0_3d_ovp_Nv/' ystart=2005 yend=2013 nyr=(yend-ystart+1) ; model grid NModelInfo = CTM_type('generic', res=[1, 1], halfpolar=1, CENTER180=1) NGridInfo = ctm_grid(NModelInfo) NIMX = NGridInfo.IMX ; Maximum I (longitude) dimension NJMX = NGridInfo.JMX ; Maximum J (latitude) dimension) nlonmid = NGridInfo.XMid ; Array of longitude centers nlatmid = NGridInfo.YMid ; Array of latitude centers nlon=n_elements(nlonmid) nlat=n_elements(nlatmid) OMIO3=fltarr(nlon, nlat, nyr, 12)-999d0 ; year month MODO3=fltarr(nlon, nlat, nyr, 12)-999d0 ; year month OMIO3raw=fltarr(288, 180, nyr, 12)-999d0 ; year month endif else if md eq 1 then begin mdname='MERRA2_GMI' mdlabel='MR2GMI' ; omimls grid 1.25x1 -60-60 OModelInfo = CTM_type('generic', res=[1.25, 1], halfpolar=0, CENTER180=0) OGridInfo = ctm_grid(OModelInfo) OIMX = OGridInfo.IMX ; Maximum I (longitude) dimension OJMX = OGridInfo.JMX ; Maximum J (latitude) dimension) olonmid = OGridInfo.XMid ; Array of longitude centers olatmid = OGridInfo.YMid ; Array of latitude centers path='/gpfsm/dnb04/projects/p76/pub/'; MERRA2-GMI run ystart=2005 yend=2017 NModelInfo = CTM_type('generic', res=[0.625, 0.5], halfpolar=1, CENTER180=1) NGridInfo = ctm_grid(NModelInfo) NIMX = NGridInfo.IMX ; Maximum I (longitude) dimension NJMX = NGridInfo.JMX ; Maximum J (latitude) dimension) nlonmid = NGridInfo.XMid ; Array of longitude centers nlatmid = NGridInfo.YMid ; Array of latitude centers nlon=n_elements(nlonmid) nlat=n_elements(nlatmid) nyr=yend-ystart+1 OMIO3=fltarr(nlon, nlat, nyr, 12)-999d0 ; year month MODO3=fltarr(nlon, nlat, nyr, 12)-999d0 ; year month OMIO3raw=fltarr(288, 180, nyr, 12)-999d0 ; year month endif else if md eq 2 then begin mdname='bench_i329_19_gmi_free_c180_132lev' mdlabel='benchmark132lev' ; omimls grid 1.25x1 -60-60 OModelInfo = CTM_type('generic', res=[1.25, 1], halfpolar=0, CENTER180=0) OGridInfo = ctm_grid(OModelInfo) OIMX = OGridInfo.IMX ; Maximum I (longitude) dimension OJMX = OGridInfo.JMX ; Maximum J (latitude) dimension) olonmid = OGridInfo.XMid ; Array of longitude centers olatmid = OGridInfo.YMid ; Array of latitude centers ;path='/gpfsm/dnb04/projects/p80/mmanyin/benchmarks/bench_i329_gmi_free_c180/inst0_3d_ovp_Nv/' path='/gpfsm/dnb02/projects/p128/mmanyin/bench_plots/inst0_3d_ovp_Nv/' ystart=2005 yend=2013 nyr=(yend-ystart+1) ; model grid NModelInfo = CTM_type('generic', res=[1, 1], halfpolar=1, CENTER180=1) NGridInfo = ctm_grid(NModelInfo) NIMX = NGridInfo.IMX ; Maximum I (longitude) dimension NJMX = NGridInfo.JMX ; Maximum J (latitude) dimension) nlonmid = NGridInfo.XMid ; Array of longitude centers nlatmid = NGridInfo.YMid ; Array of latitude centers nlon=n_elements(nlonmid) nlat=n_elements(nlatmid) OMIO3=fltarr(nlon, nlat, nyr, 12)-999d0 ; year month MODO3=fltarr(nlon, nlat, nyr, 12)-999d0 ; year month OMIO3raw=fltarr(288, 180, nyr, 12)-999d0 ; year month endif latm=nlatmid lonm=nlonmid ;;;;read in OMI O3 column data, which is derectly get from LOK mname=['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12'] yrname=['2005','2006', '2007', '2008', '2009', '2010', $ '2011', '2012', '2013','2014', '2015', '2016', '2017', '2018'] ; an extra copy of data was saved in ; /gpfsm/dnb31/jliu13/CHEM_EVA/data/O3/tco_omimls_60s60n_oct2004_to_apr2020.sav filedata='/gpfsm/dhome/jliu13/satellite/TOMSsubv_monthly_TOZ_9117.sav' restore, filedata lat=lats; 37 layser indye=where (year ge 2005 and year le yend) data=TO1Z[*, indye]; 36, 324 for yr=ystart, yend do begin for m=0, 11 do begin ambm=fltarr(2, 132) openr,lun,'ambm',/get_lun readf,lun,ambm close,/all ; close the data file am132=transpose(ambm[0, *]) bm132=transpose(ambm[1, *]) bm1=reverse(bm132) am1=reverse(am132) ;ad in am, bm filename='/discover/nobackup/jliu13/CHEM_EVA/data/parameters/MR2GMI_ambm.sav' restore, filename ;;;;;;;;;;;;;; ;;;;;;;;DONE WITH OMI DATA then READ IN MODEL DATA ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if md eq 0 or md eq 2 then begin filemod=path+mdname+'.inst0_3d_ovp_Nv.monthly.'+yrname[yr-ystart]+mname[m]+'.nc4' endif else if md eq 1 then begin filemod=path+'Y'+yrname[yr-ystart]+'/M'+mname[m]+'/'+mdname+'.inst0_3d_ovp_Nv.monthly.'+yrname[yr-ystart]+mname[m]+'.nc4' endif idu=ncdf_open(filemod, /nowrite) loid=ncdf_varid(idu,specmod) ncdf_varget,idu,loid,O3m oid=ncdf_varid(idu,"OVP14_PS"); surface pressure ncdf_varget,idu,oid,PSF for i=0, nlon-1 do begin for j=0, nlat-1 do begin Pres=am+bm*psf[i, j]/100. modtmp=reform(O3m[i, j, *]); mol/mol e-12 modtmp=reverse(modtmp); reverse saved data indle=where (modtmp gt 0); add this line if calculated is tropospheric column Pressure=Pres[indle] modt=modtmp[indle] pressureSurface=psf[i, j]/100. tt=integPr(modt, pressure, pressureSurface,constMR=0 ) modO3[i, j, yr-ystart, m] = tt/2.6867*1e-16; convert from molec/cm2 to DU for O3 and SO2 endfor ; i endfor; j endfor print, yr endfor modO3[where (modO3 le 0)]=!Values.F_NaN Otmp=mean(modO3, 1, /NaN) Otmp[where (Otmp eq 0)]=!Values.F_NaN lat1=[-90, -60, -30, 30, 60, -90, -60] lat2=[-60, -30, 30, 60, 90, 90, 60] n1=[0, 6, 12, 24, 30, 0, 6] n2=[5, 11, 23, 29, 35, 35, 29] nrs=n_elements(lat1) modO3lat=fltarr(nrs, nyr, 12)-999d0 dataO3lat=fltarr(nrs, nyr*12)-999d0 for rr=0, nrs-1 do begin indlat =where (latm ge lat1[rr] and latm le lat2[rr]) Otmp1=reform(Otmp[indlat, *,*]) modO3lat[rr, *, *]=mean(Otmp1, 1, /NaN) ns=n1[rr] nn=n2[rr] datatmp=data[ns:nn, *] dataO3lat[rr, *]=mean(datatmp, 1, /NaN) endfor dataO3lat[where (dataO3lat eq 0)]=!Values.F_NaN data=dataO3lat md=fltarr(nrs, nyr*12)-999d0 for r=0, nrs-1 do begin md[r, *]=reform(transpose(modO3lat[r, *, *]), nyr*12) endfor savefile='/gpfsm/dnb31/jliu13/CHEM_EVA/CCM_eva/sav/TOZsubv_'+mdname+'_monthly.sav' save, md, data, lat1, lat2,filename=savefile stop 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.+ystart blue=GetColor('blue', 102) Datename2 = ['05','06','07','08','09', '10','11', '12','13', '14','15'] xticknamev=time[[1, 13, 25, 37, 49, 61, 73,85,97,109,121,133]-1] for r=4, nr-1 do begin obsp=reform(obm_ts[r,*, 0]);*factor mdp=reform(mdm_ts[r,*, 0]);*factor errorobs=reform(obm_ts[r,*, 1]);*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]) high_error=obsp-obsp low_error=obsp-obsp high_error[indo] = obsp[indo]+errorobs[indo];*1.96 low_error[indo] = obsp[indo]-errorobs[indo];;*1.96 xrange=[ystart, yend+1] yrange=[-4, 4] ; 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=[10, 60], YStyle=1;,xticks=12,xtickv=xticknamev, xtickname=datename2 time2=time[indo] ; ; Fill in the error estimates. cgColorFill, [time2, Reverse(time2), time2[0]], $ [high_error[indo], Reverse(low_error[indo]), high_error[indo[0]]], $ Color='grey' 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.3)'), 2)+cgsymbol('+-')+strtrim(string(sigmao,'(f7.3)'), 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.3)'), 2)+cgsymbol('+-')+strtrim(string(sigmam,'(f7.3)'), 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.3)'), 2)+', mean ( '+strtrim(string(meano,'(f8.3)'),2)+', ';+strtrim(string(meanm,'(f8.3)'),2)+', '+string(bias,'(f8.3)')+'%)' xyouts, !x.window[0]+0.02, !y.window[1]-0.04, label, $ charsize =1,color=1, alignment=0,/normal label = strtrim(string(meanm,'(f8.3)'),2)+', ' xyouts, !x.window[0]+0.19, !y.window[1]-0.04, label, $ charsize =1,color=2, alignment=0,/normal label = string(bias,'(f8.3)')+'%)' xyouts, !x.window[0]+0.23, !y.window[1]-0.04, label, $ charsize =1,color=1, alignment=0,/normal cgText, 0.5, 0.94, 'Monthly mean time series of '+dname+', '+ mdlabel, size=2,ALIGNMENT=0.5, /NORMAL,COLOR='blue' Legend, line=[0], color=[1], nlines=6, lcolor=[1], textcolor=1, label = [dname], thick=3, $ position= [0.87, 0.90, 0.88, 0.91] legend, line=[0], color=[2], textcolor=2, lcolor=[2],label=[mdlabel],/ADD, thick=3 endfor ;cgPS_Close close_device , /TIMESTAMP set_plot,olddevice close,/all print, 'DONE' stop ;if md eq 0 then begin save, modO3, OMIO3raw, OMIO3, ystart, yend,lonr, latr, lonm, latm, filename=savefile ;endif else if md eq 1 then begin ;save, modO3,OMIO3y, ystart, yend,filename=savefile ;endif stop end