; code to generate NO2 plots ; Figure 1: monthly map averaged over satellite period ; Figure 2: annual mean map averaged over satellite period and ; Figure 3: time-latitude cross section plot ; Figure 4: time series plot of the monthly mean ; Figure 5: seasonal cycle plot ; Figure 6: boxplot ;;;TVPLOT add ; /NOZINTERP -> set to switch off the vertical ; interpolation. This can be used only if the ZARR is ; evenly spaced. ; ; modified by J. Liu 08/27/2020, convert to 1 page plot ; Modified by J. Liu 10/28/2020, ; 1: fomat title of each plot, ; 2: add choice for c180q grid, res eq 0.5 ; the MR2GMI res eq 0.52, c180, 0.625x0.5 ; 3: update OMI NO2 data into V4 @GetColor.pro @cgboxplot_jliu.pro @integpr.pro pro no2_comp path='/gpfsm/dnb04/projects/p76/pub/'; MERRA2-GMI run run_name='MERRA2_GMI'; run name collection='inst0_3d_ovp_Nv' res=0.52; c360, c180, c90, c48 ystart=2005 yend=2018 ;@exp_specifics.inc ; for example: ;expdir='/gpfsm/dnb04/projects/p80/mmanyin/benchmarks/bench_i329_gmi_free_c180' ;run_name='bench_i329_gmi_free_c180' ;collection='inst0_3d_ovp_Nv' ;ystart=2005 ;yend=2014 ;figdir='/discover/nobackup/jliu13/CHEM_EVA/FIG/1page' ;res=1 nlevin=72 ; new test for bench_10-14-1_gmi_free_c180_72lev run ;expdir='/discover/nobackup/projects/gmao/ccmdev/bvanaart/run/bench_10-14-1_gmi_free_c180_72lev' ;run_name='bench_10-14-1_gmi_free_c180_72lev' ;collection='inst0_3d_ovp_Nv' ;ystart=2005 ;yend=2013 figdir='/discover/nobackup/jliu13/CHEM_EVA/FIG/'+run_name+'/' ;res=0.5 ;nlevin=72 ; collection is expected to include these fields: ; OVP14_NO2 ; OVP14_PS ; OVP14_TROPP ;path= expdir + '/' + collection + '/' ;add species='NO2' sptype='TRC' ob_platform='OMI' sat_name='OMI_trno2' specmod="OVP14_NO2" factor=1e-15; for unit change unit_name='(1e15 molec/cm2)' yrname= STRING(INDGEN(yend-ystart+1)+ystart) mname=['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12'] nyrs=n_elements (yrname) nmons=n_elements(mname) avg_name=strtrim(yrname[0], 2)+'-'+strtrim(yrname[nyrs-1],2); time range ; /add min_map=0 max_map=5 min_diff=-3 max_diff=3 yrange1=[-1,-1,-1,-2, -1, -1.5,-0.5, -0.5, -2, -0.5, -2, -2, -0.5, -2, -2, -0.5];*0.45 yrange2=[1, 2, 2, 5, 3, 16, 4, 2, 8, 10, 5, 10, 1.5, 4, 4, 1.0];*0.45 rname=['90S60S', '60S20S', '20S20N','20N60N', '60N90N', 'E.Asia', 'India',$ 'MidEast', 'Europe', 'E.US', 'S.Africa', 'C.Africa', 'Indonesia', 'W.US', $ 'S.America', 'Australia']; jliu, 08/27/2020 if res eq 0.5 then begin ModelInfo = CTM_type('generic', res=[0.5, 0.5], halfpolar=1, CENTER180=1) GridInfo = ctm_grid(ModelInfo) IMX = GridInfo.IMX ; Maximum I (longitude) dimension JMX = GridInfo.JMX ; Maximum J (latitude) dimension) lonmod = GridInfo.XMid ; Array of longitude centers latmod = GridInfo.YMid ; Array of latitude centers resname='c180q' resgrid='0.5x0.5' nlonmod=n_elements(lonmod) nlatmod=n_elements(latmod) endif else if res eq 0.52 then begin ModelInfo = CTM_type('generic', res=[0.625, 0.5], halfpolar=1, CENTER180=1) GridInfo = ctm_grid(ModelInfo) IMX = GridInfo.IMX ; Maximum I (longitude) dimension JMX = GridInfo.JMX ; Maximum J (latitude) dimension) lonmod = GridInfo.XMid ; Array of longitude centers latmod = GridInfo.YMid ; Array of latitude centers resname='c180' resgrid='0.50x0.625' nlonmod=n_elements(lonmod) nlatmod=n_elements(latmod) endif else if res eq 1 then begin ;-----------------------monthly mean----------------- ModelInfo = CTM_type('generic', res=[1, 1],halfpolar=1, CENTER180=1 ) GridInfo = ctm_grid(ModelInfo) IMX = GridInfo.IMX ; Maximum I (longitude) dimension JMX = GridInfo.JMX ; Maximum J (latitude) dimension) lonmod = GridInfo.XMid ; Array of longitude centers latmod = GridInfo.YMid ; Array of latitude centers resname='c90' resgrid='1.00x1.00' nlonmod=n_elements(lonmod) nlatmod=n_elements(latmod) endif ;;;;;;;;;;;;;;;;;;;;;;;; ; readin the model path ;;;;;;;;;;;;;;;;;;;;;;;;; ; read in am, bm If nlevin eq 72 then filename='/discover/nobackup/jliu13/CHEM_EVA/data/parameters/MR2GMI_ambm.sav' if nlevin eq 132 then filename= '/discover/nobackup/jliu13/CHEM_EVA/data/parameters/ambm_132.sav' restore, filename datacol_ob_mon=fltarr(nlonmod, nlatmod, nyrs, nmons)-999d0;lon, lat, year, month datacol_md_mon=fltarr(nlonmod, nlatmod, nyrs, nmons)-999d0; for y=0, nyrs-1 do begin for m=0, nmons-1 do begin ; add filedata='/discover/nobackup/jliu13/CHEM_EVA/data/'+species+'/'+resname+'/'+sat_name+'_'+resgrid+'_'+strtrim(yrname[y],2)+mname[m]+'_V4.sav' ;save,data, nlonmid, nlatmid, filename=savefile ;print, filedata ; Junhua if file_exist(filedata) eq 1 then begin print, 'reading omi data' restore, filedata ;if res eq 0.52 then data=NO2 datacol_ob_mon[*,*, y, m]=data endif else begin datacol_ob_mon[*,*, y, m]=!Values.F_NaN print, 'no data' endelse ;;;;;;;;DONE WITH OMI DATA then READ IN MODEL DATA ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if res eq 0.52 then begin ; specific foor MERRA2-GMI run filemod=path+'Y'+strtrim(yrname[y],2)+'/M'+mname[m]+'/'+run_name+'.'+collection+'.monthly.'+strtrim(yrname[y],2)+mname[m]+'.nc4' endif else begin filemod=path+run_name+'.'+collection+'.monthly.'+strtrim(yrname[y],2)+mname[m]+'.nc4' endelse ;filemod=path+'Y'+strtrim(yrname[y],2)+'/M'+mname[m]+'/MERRA2_GMI.inst0_3d_ovp_Nv.monthly.'+strtrim(yrname[y],2)+mname[m]+'.nc4' idu=ncdf_open(filemod, /nowrite) ; read in model species: O3, CO, SO2, CH2O, NO2 loid=ncdf_varid(idu,specmod) ncdf_varget,idu,loid,moddata loid=ncdf_varid(idu,"OVP14_PS"); surface pressure ncdf_varget,idu,loid,PS loid=ncdf_varid(idu,"OVP14_TROPP"); tropopause ncdf_varget,idu,loid,tropp for i=0, nlonmod-1 do begin for j=0, nlatmod-1 do begin trop=tropp[i, j]/100. if trop le 700 then begin ; omi has data; troppopause seem has some missing value Pres=am+bm*ps[i, j]/100. modtmp=reform(moddata[i, j, *]); mol/mol e-12 modtmp=reverse(modtmp); reverse saved data indle=where (Pres ge trop); add this line if calculated is tropospheric column Pressure=Pres[indle] modt=modtmp[indle] pressureSurface=ps[i, j]/100. tt=integPr(modt, pressure, pressureSurface,constMR=0 ) datacol_md_mon[i, j, y, m] = tt;/2.6867*1e-16; convert from molec/cm2 to DU for O3 and SO2 endif endfor ; i endfor; j endfor endfor ; eliminate all the negative data and do the same for model datacol_md_mon(where (datacol_md_mon eq -999d0))=!Values.F_NAN datacol_ob_mon(where (datacol_ob_mon lt 0))=!Values.F_NAN; datacol_md_mon(where (finite(datacol_ob_mon) eq 0))=!Values.F_NAN ; There are some negative values jliu 10/28/2020 ;savefile='/discover/nobackup/jliu13/CHEM_EVA/sav/'+sat_name+'_'+run_name+'_'+resname+'_'+species+'_monthly_mean.sav' ;save, datacol_ob_mon, datacol_md_mon, yrname, lonmod, latmod, filename=savefile ;;; for plotting ;restore, savefile ; seasonal cycle averaged over satellite period ob_mon=mean(datacol_ob_mon, 3, /NaN); get the monthly mean md_mon=mean(datacol_md_mon, 3, /NaN); get the monthly mean ob_mon[where (ob_mon eq 0)]=!Values.F_NAN md_mon[where (md_mon eq 0)]=!Values.F_NAN ; annnual mean md_anu=mean(md_mon, 3, /NaN) ob_anu=mean(ob_mon, 3, /NaN) ob_anu[where (ob_anu eq 0)]=!Values.F_NAN md_anu[where (md_anu eq 0)]=!Values.F_NAN ; monthly zonal mean month-latitude ob_za=mean(datacol_ob_mon, 1, /NaN); get the monthly mean md_za=mean(datacol_md_mon, 1, /NaN); get the monthly mean ob_za[where (ob_za eq 0)]=!Values.F_NAN md_za[where (md_za eq 0)]=!Values.F_NAN ob_zats=fltarr(nyrs*nmons, nlatmod) md_zats=fltarr(nyrs*nmons, nlatmod) for ys=0, nyrs-1 do begin for m=0, nmons-1 do begin for la=0, nlatmod-1 do begin ob_zats[ys*12+m, la]=ob_za[la, ys, m] md_zats[ys*12+m, la]=md_za[la, ys, m] endfor endfor endfor ;md_zats=reform(transpose(md_za), nyrs*nmons, nlatmod); correct has been double checked Datename = ['January','February','March','April','May','June','July', 'August','September','October','November','December'] ;================================================ ; plot settings ;;;;map plotting ; FIGURE 1: monthly mean map limit=[-70, -180, 70, 180]; map boundary for m = 0, nmons-1 do begin ;fps = figdir + '/' + 'FIG1_'+sat_name+'_'+run_name+'_clim_monthlymean_3panels.ps' ;fps = figdir + '/' + 'FIG1_'+sat_name+'_'+run_name+'_clim_'+Datename[m]+'_3panels.ps' fps = figdir + '/' + run_name+'.'+species+'.'+sptype+'.'+ob_platform+'.clim_map.'+Datename[m]+'_pt.ps' open_device, olddevice, /ps, filename=fps, Bits=8, $ /color, /portrait ; jLiu 08/27/2020 !X.OMARGIN = [6, 6] !Y.OMARGIN = [4, 4] !X.MARGIN = [6, 6] !Y.MARGIN = [6, 6] !P.CHARTHICK = 1.5 !P.Charsize = 2 !X.THICK = 2 !Y.THICK = 2 !P.font = 1.5 nrow=3 ncol=1 !P.Multi = [1, ncol, nrow, 0, 0] multipanel,rows=nrow,cols=ncol mdplot=reform(md_mon[*, *, m])*factor obplot=reform(ob_mon[*,*, m])*factor t1=sat_name+' '+datename[m]+ ' ('+avg_name+' mean)' t2=run_name+' '+datename[m]+ ' ('+avg_name+' mean)' t3= 'DIFF(MOD-OMI)' ;================================================ ; myct, 33, Ncolors=18 myct, 33, Ncolor=12 tvmap, obplot,lonmod,latmod,/sample, title=t1,TCSFAC=1.5, $ limit=limit, CsFac=1, TitleCsFac = 2, Margin = [ 0.04, 0.04, 0.04, 0.04 ], $ /countries, /continent, /coasts,grid=1, cbar=1,CBVERTICAL=1,cbposition=[0.9, 0.01, 0.92,0.98],divisions=7, Thick=2, $ NoGXLabels=0, NoGYLabels=0, WindowPosition=WindowPosition, $ mindata=min_map, maxdata=max_map, charsize=1.5, Nan_color=!myct.gray; Map_Set, 0, 0, Position=WindowPosition, Limit=Limit, Grid=0, Color=1, /noerase label = unit_name xyouts, !x.window[0]+0.72, !y.window[1]-0.23, label, $ charsize =1,color=1, alignment=0,/normal myct, 33, Ncolor=12 ;================================================ tvmap, mdplot,lonmod,latmod,/sample, title=t2,TCSFAC=1.5, $ limit=limit, CsFac=1, TitleCsFac = 2, Margin =[ 0.04, 0.04, 0.04, 0.04 ] , $ /countries, /continent, /coasts,grid=1, cbar=1,CBVERTICAL=1,cbposition=[0.9, 0.01, 0.92,0.98], divisions=7, Thick=2, $ NoGXLabels=0, NoGYLabels=0, WindowPosition=WindowPosition, $ mindata=min_map, maxdata=max_map, charsize=1.5, Nan_color=!myct.gray Map_Set, 0, 0, Position=WindowPosition, Limit=Limit, Grid=0, Color=1, /noerase label = unit_name xyouts, !x.window[0]+0.72, !y.window[1]-0.23, label, $ charsize =1,color=1, alignment=0,/normal DIFF=mdplot-obplot ;================================================ myct2, -1;, Ncolor=12 tvmap, DIFF,lonmod,latmod,/ sample, title=t3,TCSFAC=1.5, $ limit=limit, CsFac=1, TitleCsFac = 2, Margin =[ 0.04, 0.04, 0.04, 0.04 ], $ /countries, /continent, /coasts,grid=1, cbar=1,CBVERTICAL=1,divisions=8,cbposition=[0.9, 0.01, 0.92,0.98], Thick=2, $ NoGXLabels=0, NoGYLabels=0, WindowPosition=WindowPosition, $ mindata=min_diff, maxdata=max_diff, charsize=1.5, Nan_color=!myct.gray; Map_Set, 0, 0, Position=WindowPosition, Limit=Limit, Grid=0, Color=1, /noerase label = unit_name xyouts, !x.window[0]+0.72, !y.window[1]-0.23, label, $ charsize =1,color=1, alignment=0,/normal close_device;, /TIMESTAMP set_plot,olddevice close, /all endfor ;;;;;;FIgure 2 annual map;;;;;;;;; ;fps = figdir + '/' + 'FIG2_'+sat_name+'_'+run_name+'_clim_annualmean_3panels.ps' fps = figdir + '/' + run_name+'.'+species+'.'+sptype+'.'+ob_platform+'.clim_map.annual_pt.ps' open_device, olddevice, /ps, filename=fps, Bits=8, $ /color, /portrait !X.OMARGIN = [6, 6] !Y.OMARGIN = [4, 4] !X.MARGIN = [6, 6] !Y.MARGIN = [6, 6] !P.CHARTHICK = 1.5 !P.Charsize = 2 !X.THICK = 2 !Y.THICK = 2 !P.font = 1.5 nrow=3 ncol=1 !P.Multi = [1, ncol, nrow, 0, 0] multipanel,rows=nrow,cols=ncol limit=[-70, -180, 70, 180] mdplot=reform(md_anu[*, *])*factor obplot=reform(ob_anu[*,*])*factor t1=sat_name+ ' annual mean ('+avg_name+')' t2=run_name+ ' annual mean ('+avg_name+')' t3= 'DIFF(MOD-OMI)' ;================================================ myct, 33, Ncolors=12 tvmap, obplot,lonmod,latmod,/sample, title=t1,TCSFAC=1.5, $ limit=limit, CsFac=1, TitleCsFac = 2, Margin = [ 0.04, 0.04, 0.04, 0.04 ], $ /countries, /continent, /coasts,grid=1, cbar=1,CBVERTICAL=1,cbposition=[0.9, 0.01, 0.92,0.98], divisions=7, Thick=2, $ NoGXLabels=0, NoGYLabels=0, WindowPosition=WindowPosition, $ mindata=min_map, maxdata=max_map, charsize=1.5, Nan_color=!myct.gray; Map_Set, 0, 0, Position=WindowPosition, Limit=Limit, Grid=0, Color=1, /noerase label = unit_name xyouts, !x.window[0]+0.72, !y.window[1]-0.23, label, $ charsize =1,color=1, alignment=0,/normal ;================================================ tvmap, mdplot,lonmod,latmod,/sample, title=t2,TCSFAC=1.5, $ limit=limit, CsFac=1, TitleCsFac = 2, Margin =[ 0.04, 0.04, 0.04, 0.04 ] , $ /countries, /continent, /coasts,grid=1, cbar=1,CBVERTICAL=1,cbposition=[0.9, 0.01, 0.92,0.98],divisions=7, Thick=2, $ NoGXLabels=0, NoGYLabels=0, WindowPosition=WindowPosition, $ mindata=min_map, maxdata=max_map, charsize=1.5, Nan_color=!myct.gray; Map_Set, 0, 0, Position=WindowPosition, Limit=Limit, Grid=0, Color=1, /noerase label = unit_name xyouts, !x.window[0]+0.72, !y.window[1]-0.23, label, $ charsize =1,color=1, alignment=0,/normal ;================================================ myct2, -1 DIFF=mdplot-obplot ;DIFF[where (obplot eq -999d0 and mdplot eq -999d0 )]=-999d0 ;min_valid=min(DIFF[where (DIFF ne -999)]) tvmap, DIFF,lonmod ,latmod,/sample, title=t3,TCSFAC=1.5, $ limit=limit, CsFac=1, TitleCsFac = 2, Margin =[ 0.04, 0.04, 0.04, 0.04 ], $ /countries, /continent, /coasts,grid=1, cbar=1,CBVERTICAL=1,cbposition=[0.9, 0.01, 0.92,0.98],divisions=8, Thick=2, $ NoGXLabels=0, NoGYLabels=0, WindowPosition=WindowPosition, $ mindata=min_diff, maxdata=max_diff, charsize=1.5, Nan_color=!myct.gray Map_Set, 0, 0, Position=WindowPosition, Limit=Limit, Grid=0, Color=1, /noerase label = unit_name xyouts, !x.window[0]+0.72, !y.window[1]-0.23, label, $ charsize =1,color=1, alignment=0,/normal close_device;, /TIMESTAMP set_plot,olddevice close, /all ;fps = figdir + '/' + 'FIG3_'+sat_name+'_'+run_name+'_time_lat_zonalmean_3panel.ps' fps = figdir + '/' + run_name+'.'+species+'.'+sptype+'.'+ob_platform+'.time_lat_zonalmean_pt.ps' open_device, olddevice, /ps, filename=fps, Bits=8, $ /color, /portrait !X.OMARGIN = [6, 6] !Y.OMARGIN = [4, 4] !X.MARGIN = [6, 6] !Y.MARGIN = [6, 6] !P.CHARTHICK = 1.5 !P.Charsize = 2 !X.THICK = 2 !Y.THICK = 2 !P.font = 1.5 nrow=3 ncol=1 !P.Multi = [1, ncol, nrow, 0, 0] multipanel,rows=nrow,cols=ncol limit=[-70, -180, 70, 180] mdplot=reform(md_anu[*, *])*factor obplot=reform(ob_anu[*,*])*factor t1=sat_name+ ' annual mean ('+avg_name+')' t2=run_name+ ' annual mean ('+avg_name+')' t3= 'DIFF(MOD-OMI)' ; plot the IAV of zonal mean : X month from 2005 jan to 2014, Y: latitude grid indlatp=where (latmod ge -60 and latmod le 60) mdplot=md_zats[*, indlatp]*factor obplot=ob_zats[*, indlatp]*factor latp=latmod[indlatp] t1=sat_name+ ' ('+avg_name+')' t2=run_name+ ' ('+avg_name+')' t3= 'DIFF(MOD-OMI)' X=findgen(nyrs*nmons)/12.+ystart yname=['60N', '40N', '20N', '0', '20S', '40S', '60S'] myct, 33, Ncolors=12 ;obplot[where (finite(obplot) eq 0)]=-999d0 ;mdplot[where (finite(mdplot) eq 0)]=-999d0 tvplot, obplot,X ,latp, title=t1,CSFAC=2.5,yrange=[-60, 60], $ Margin =[ 0.04, 0.04, 0.04, 0.04 ] , cbar=1,CBVERTICAL=1,cbposition=[0.9, 0.01, 0.92,0.98], divisions=7, Thick=2,ytitle='Latitude',xtitle='Year', Nan_color=!myct.gray, mindata=min_map/2, maxdata=max_map/2, /NOZINTERP ;, min_valid=-10 Map_Set, 0, 0, Position=WindowPosition, Grid=0, Color=1, /noerase ; ;label = unit_name ;xyouts, !x.window[0]+0.72, !y.window[1]-0.23, label, $ ; charsize =1,color=1, alignment=0,/normal tvplot, mdplot,X ,latp, title=t2,CSFAC=2.5,yrange=[-60, 60], /NOZINTERP , $ Margin =[ 0.04, 0.04, 0.04, 0.04 ] , $ grid=1, cbar=1,CBVERTICAL=1,cbposition=[0.9, 0.01, 0.92,0.98],divisions=7, Thick=2, mindata=min_map/2, maxdata=max_map/2, ytitle='Latitude',xtitle='Year', Nan_color=!myct.gray;, B ;Map_Set, 0, 0, Position=WindowPosition, Grid=0, Color=1, /noerase ;label = unit_name ; xyouts, !x.window[0]+0.72, !y.window[1]-0.23, label, $ ; charsize =1,color=1, alignment=0,/normal ;================================================ ;================================================ myct2, -1 DIFF=mdplot-obplot ;DIFF[where (obplot eq -999 and mdplot eq -999)]=-999d0 tvplot, diff,X ,latp,title=t3, csfac=2.5,yrange=[-60, 60], /NOZINTERP , $ TitleCsFac = 1.5, Margin =[ 0.04, 0.04, 0.04, 0.04 ], ytitle='Latitude',xtitle='Year', $ grid=1, cbar=1,CBVERTICAL=1,cbposition=[0.9, 0.01, 0.92,0.98],divisions=8, Thick=2, $ mindata=min_diff/2, maxdata=max_diff/2, nan_color=!myct.gray;, min_valid=-100; ;Map_Set, 0, 0, Position=WindowPosition, Limit=Limit, Grid=0, Color=1, /noerase label = unit_name xyouts, !x.window[0]+0.72, !y.window[1]-0.03, label, $ charsize =1,color=1, alignment=0,/normal ; close_device;, /TIMESTAMP set_plot,olddevice close, /all ;;; FIG4-6 time series plots of selected regions. ; region info is saved in mkfig_map_region.pro savefile='/discover/nobackup/jliu13/CHEM_EVA/sav/regions_info.sav' ;save, regions, limit1, limit2, limit3, limit4, filename=savefile restore, savefile nr=n_elements(limit1) ob_mts=fltarr(nr, nyrs, nmons, 2)-999d0; 12 regions, 12 years and 12 months, 0 is mean and 1 is standard deviation md_mts=fltarr(nr, nyrs, nmons, 2)-999d0; 12 regions, 12 years and 12 months, 0 is mean and 1 is standard deviation ob_ats=fltarr(nr, nyrs, 2)-999d0; 12 regions, 12 years and 12 months, 0 is mean and 1 is standard deviation md_ats=fltarr(nr, nyrs, 2)-999d0; 12 regions, 12 yeaer, 0 is mean and 1 is standard deviation ob_season=fltarr(nr, nmons, 2)-999d0; 12 regions, 12 years and 12 months, 0 is mean and 1 is standard deviation md_season=fltarr(nr, nmons, 2)-999d0; 12 regions, 12 months, 0 is mean and 1 is standard deviation for r=0, nr-1 do begin indlon= where (lonmod ge limit2[r] and lonmod le limit4[r], nlon) indlat=where (latmod ge limit1[r] and latmod le limit3[r], nlat) obtmp1=reform(datacol_ob_mon[indlon, *, *,*]) obtmp2=reform(obtmp1[*, indlat, *,*]) mdtmp1=reform(datacol_md_mon[indlon, *, *,*]) mdtmp2=reform(mdtmp1[*, indlat, *,*]) ;;;;;; 1111111calculate the monthly mean and its std obtmp=fltarr(nyrs, nmons, 2)-999d0 mdtmp=fltarr(nyrs, nmons, 2)-999d0 for y=0, nyrs-1 do begin for m=0, nmons-1 do begin obtmp3=obtmp2[*, *, y, m] indob=where (finite (obtmp3) eq 1, n1) if (indob[0] ne -1 and n1 ge 2) then begin obtmp[y, m, 0]=mean(obtmp3[indob]) obtmp[y, m, 1]=stdev(obtmp3[indob]) endif else begin obtmp[y, m, *]=!Values.F_NAN endelse mdtmp3=mdtmp2[*, *, y, m] indmd=where (finite (mdtmp3) eq 1, n2) if (indmd[0] ne -1 and n2 ge 2) then begin mdtmp[y, m, 0]=mean(mdtmp3[indmd]) mdtmp[y, m, 1]=stdev(mdtmp3[indmd]) endif else begin mdtmp[y, m, *]=!Values.F_NAN endelse endfor endfor ob_mts[r, *, *, *]=obtmp md_mts[r, *, *, *]=mdtmp undefine, obtmp undefine, mdtmp ;;;;;; 2222222 calculate the annual mean and its std obtmp=fltarr(nyrs, 2)-999d0; 12 year mdtmp=fltarr(nyrs, 2)-999d0 for y=0, nyrs-1 do begin obtmp3=obtmp2[*, *, y, *] indob=where (finite (obtmp3) eq 1, n1) if (indob[0] ne -1 and n1 ge 2) then begin obtmp[y, 0]=mean(obtmp3[indob]) obtmp[y, 1]=stdev(obtmp3[indob]) endif else begin obtmp[y, *]=!Values.F_NAN endelse mdtmp3=mdtmp2[*, *, y, *] indmd=where (finite (mdtmp3) eq 1, n2) if (indmd[0] ne -1 and n2 ge 2) then begin mdtmp[y, 0]=mean(mdtmp3[indmd]) mdtmp[y, 1]=stdev(mdtmp3[indmd]) endif else begin mdtmp[y, *]=!Values.F_NAN endelse endfor ob_ats[r, *, *]=obtmp md_ats[r, *, *]=mdtmp undefine, obtmp undefine, mdtmp ;;;;;; 33333333 calculate the seasonal cycle mean and its std obtmp=fltarr(nmons, 2)-999d0; 12 year mdtmp=fltarr(nmons, 2)-999d0 for m=0, nmons-1 do begin obtmp3=obtmp2[*, *, *, m] indob=where (finite (obtmp3) eq 1, n1) if (indob[0] ne -1 and n1 ge 2) then begin obtmp[m, 0]=mean(obtmp3[indob]) obtmp[m, 1]=stdev(obtmp3[indob]) endif else begin obtmp[m, *]=!Values.F_NAN endelse mdtmp3=mdtmp2[*, *, *, m] indmd=where (finite (mdtmp3) eq 1, n2) if (indmd[0] ne -1 and n2 ge 2) then begin mdtmp[m, 0]=mean(mdtmp3[indmd]) mdtmp[m, 1]=stdev(mdtmp3[indmd]) endif else begin mdtmp[m, *]=!Values.F_NAN endelse endfor ob_season[r, *, *]=obtmp md_season[r, *, *]=mdtmp undefine, obtmp undefine, mdtmp endfor obm_ts=fltarr(nr, nyrs*nmons, 2)-999d0; 12 regions, 12 years and 12 months, 0 is mean and 1 is standard deviation mdm_ts=fltarr(nr, nyrs*nmons, 2)-999d0; 12 regions, 12 years and 12 months, 0 is mean and 1 is standard deviation ;;;;; convert to monthly time series for r=0, nr-1 do begin for t=0, 1 do begin obtmp=reform (transpose(ob_mts[r, *, *, t]), nyrs*nmons) obm_ts[r, *, t]=obtmp mdtmp=reform (transpose (md_mts[r, *, *, t]), nyrs*nmons) mdm_ts[r, *, t]=mdtmp endfor endfor X1=findgen(nyrs*nmons)/12.+ystart X2=findgen(nyrs)+ystart X3=findgen(12) 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] rrr=[1, 2, 3, 5, 6, 7, 8, 9, 10, 11] for r=0, nr-1 do begin ;fps = figdir + '/' + 'FIG4_'+sat_name+'_'+run_name+'_mon_tsplot_12regions.ps' ; jliu 08/27/2020 ;fps = figdir + '/' + 'FIG4_'+sat_name+'_'+run_name+'_mon_tsplot_'+rname[r]+'.ps' fps = figdir + '/' + run_name+'.'+species+'.'+sptype+'.'+ob_platform+'.mon_tsplot.'+rname[r]+'_ls.ps' open_device, olddevice, /ps, filename=fps, $ Bits=8, /color, /landscape !X.OMARGIN=[4,4] !Y.OMARGIN=[3,3] !X.MARGIN=[3,3] !Y.MARGIN=[2,2] !P.CHARTHICK=1.2 !P.Charsize = 1.2 !X.THICK=4 !Y.THICK=4 !P.font = 1.1 ;!P.Multi = [0, 2, 1, 0 ,0] multipanel,rows=1,cols=1 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) 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=[-10, 10] ; Draw the line plot with no data cgPlot, time[indo], obsp[indo], Title=regions[r], XTitle=xtitle, YTitle=species, $ Position=position, /NoData,xrange=xrange,Charsize=1.5, YRange=[yrange1[r], yrange2[r]], 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', $ ; SymSize=0.21;, Thick=6 cgPlots, time[indo], obsp[indo], Color='black',Thick=6, linestyle=0;,EQN_SAMPLES=100;, PSym=-16, SymColor='black', $ ; SymSize=0.21;, Thick=6 cgText, 0.5, 0.94, 'Monthly mean time series of '+species+' (OMI vs '+run_name+')', size=2,ALIGNMENT=0.5, /NORMAL,COLOR='blue' Legend, line=[0], color=[1], nlines=6, lcolor=[1], textcolor=1, label = ['OMI'], thick=3, $ position= [0.50, 0.83, 0.51, 0.84] legend, line=[0], color=[2], textcolor=2, lcolor=[2],label=[run_name],/ADD, thick=3 ;cgPS_Close close_device ;, /TIMESTAMP set_plot,olddevice close,/all print, 'DONE' endfor ;;;;;;;;;seasonal cycle time=findgen(nmons) blue=GetColor('blue', 102) Datename2 = ['J','F','M','A','M', 'J','J', 'A','S', 'O','N', 'D'] xticknamev=time[[1, 2, 3,4,5, 6, 7, 8, 9,10,11,12]-1] for r=0, nr-1 do begin ;fps = figdir + '/' + 'FIG5_'+sat_name+'_'+run_name+'_seasonalcycle_12regions.ps' ;fps = figdir + '/' + 'FIG5_'+sat_name+'_'+run_name+'_seasonalcycle_'+rname[r]+'.ps' fps = figdir + '/' + run_name+'.'+species+'.'+sptype+'.'+ob_platform+'.seasonalcycle.'+rname[r]+'_ls.ps' open_device, olddevice, /ps, filename=fps, $ Bits=8, /color, /landscape !X.OMARGIN=[4,4] !Y.OMARGIN=[3,3] !X.MARGIN=[3,3] !Y.MARGIN=[2,2] !P.CHARTHICK=1.2 !P.Charsize = 2 !X.THICK=4 !Y.THICK=4 !P.font = 1.1 ;!P.Multi = [0, 2, 1, 0 ,0] multipanel,rows=1,cols=1 obsp=reform(ob_season[r,*, 0])*factor mdp=reform(md_season[r,*, 0])*factor errorobs=reform(ob_season[r,*, 1])*factor indo=where (Finite(obsp) EQ 1, nindo) indm=where (Finite(mdp) EQ 1) 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 yrange=[-2, 10] ; Draw the line plot with no data cgPlot, time[indo], obsp[indo], Charsize=1.5, Title=regions[r], XTitle=xtitle, YTitle=species, $ Position=position, /NoData, YRange=[yrange1[r], yrange2[r]], 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', $ ; SymSize=0.21;, Thick=6 cgPlots, time[indo], obsp[indo], Color='black',Thick=6, linestyle=0;,EQN_SAMPLES=100;, PSym=-16, SymColor='black', $ ; SymSize=0.21;, Thick=6 cgText, 0.5, 0.94, 'Seasonal cycle of '+species+' (05-16)', size=2,ALIGNMENT=0.5, /NORMAL,COLOR='blue' Legend, line=[0], color=[1], nlines=6, lcolor=[1], textcolor=1, label = ['OMI'], thick=3, $ position= [0.50, 0.83, 0.51, 0.84] legend, line=[0], color=[2], textcolor=2, lcolor=[2],label=[run_name],/ADD, thick=3 ;cgPS_Close close_device ;, /TIMESTAMP set_plot,olddevice close,/all print, 'DONE' endfor maxdata_mon = [0.5, 0.5, 0.5, 0.5, 2, 2] ;fps = '/discover/nobackup/jliu13/CHEM_EVA/FIG/boxplot_'+species+'.ps' maxdata=[1, 1, 2, 3, 4, 8,4,3, 6,8, 6, 8] mindata=[0,0,-1, -1, -1, -1, -1, -1, -1,-1,-1, -1] for r=0, nr-1 do begin ;fps = figdir + '/' + 'FIG6_'+sat_name+'_'+run_name+'_boxplot_12regions.ps' ;fps = figdir + '/' + 'FIG6_'+sat_name+'_'+run_name+'_boxplot_'+rname[r]+'.ps' fps = figdir + '/' + run_name+'.'+species+'.'+sptype+'.'+ob_platform+'.boxplot.'+rname[r]+'_ls.ps' open_device, olddevice, /ps, filename=fps, Bits=8, $ /color, /landscape !X.OMARGIN = [8, 8] !Y.OMARGIN = [4, 4] !X.MARGIN = [3.5, 3.5] !Y.MARGIN = [2, 2] !P.CHARTHICK = 1.5 !P.Charsize = 1.5 !X.THICK = 2 !Y.THICK = 2 !P.font = 1.1 multipanel,rows=1,cols=1 cgPlot, [0,1], /nodata, yrange=[yrange1[r]-0.1,yrange2[r]+0.1], Title=sat_name+' at '+regions[r], xrange=[0,25], $ xtickformat='(A1)', ytitle=species+' '+unit_name, YStyle=1 index = indgen(12)*2 width = ((!X.CRange[1] - !X.Crange[0]) / (20)) * 0.4 indlon= where (lonmod ge limit2[r] and lonmod le limit4[r], nlon) indlat=where (latmod ge limit1[r] and latmod le limit3[r], nlat) obtmp1=reform(datacol_ob_mon[indlon, *, *,*]) obtmp2=reform(obtmp1[*, indlat, *,*]) mdtmp1=reform(datacol_md_mon[indlon, *, *,*]) mdtmp2=reform(mdtmp1[*, indlat, *,*]) ;;;;;; 33333333 Get data over each region for each month for m=0, 11 do begin obtmp3=obtmp2[*, *, *, m] indob=where (finite (obtmp3) eq 1, n1) if (indob[0] ne -1 and n1 ge 2) then begin obtmp=obtmp3[indob]*factor cgBoxPlot_jliu, obtmp,/overplot, XLOCATION=index[m]+0.75, WIDTH=width, $ BOXCOLOR='black', /FILLBOX endif mdtmp3=mdtmp2[*, *, *, m] indmd=where (finite (mdtmp3) eq 1, n2) if (indmd[0] ne -1 and n2 ge 2) then begin mdtmp=mdtmp3[indmd]*factor;/1e15 cgBoxPlot_jliu, mdtmp,/overplot, XLOCATION=index[m]+1.25, WIDTH=width, $ BOXCOLOR='red', /FILLBOX endif endfor ; labels = ['J', 'F', 'M', 'A', 'M', 'J', 'J', 'A', 'S','O', 'N', 'D'] if yrange2[r] ge 10 then begin for j=0,11 do cgText, (index+1)[j], yrange1[r]-1, labels[j],size=2, Alignment=0.5, color='black' cgText,index[7],yrange2[r]-1, sat_name, Alignment=0, size=1.2, color='black' cgText, index[7], yrange2[r]-1.5,run_name, Alignment=0,size=1.2, color='red' endif else if yrange2[r] lt 10 and yrange2[r] gt 2 then begin for j=0,11 do cgText, (index+1)[j], yrange1[r]-0.5, labels[j],size=2, Alignment=0.5, color='black' cgText,index[7],yrange2[r]-0.2, sat_name, Alignment=0, size=1.2, color='black' cgText, index[7], yrange2[r]-0.5,run_name, Alignment=0,size=1.2, color='red' endif else begin for j=0,11 do cgText, (index+1)[j], yrange1[r]-0.3, labels[j],size=2, Alignment=0.5, color='black' cgText,index[7],yrange2[r]-0.1, sat_name, Alignment=0, size=1.2, color='black' cgText, index[7], yrange2[r]-0.2,run_name, Alignment=0,size=1.2, color='red' endelse close,/all close_device ;, /TIMESTAMP endfor stop stop end