cheyenne5-/glade/scratch/najiby>cat Annualaverage_MOPITTCO_nyg_modified.ncl ; Annualaverage_MOPITTCO_nyg_modified.ncl ; ; This code loads MOPITT CO measurements ; and applies averaging kernels to model results ; an annual average in 2013 ; ; To use type at command line: ; ncl Annualaverage_MOPITTCO_rrb_modified.ncl ; ; rrb 20201210 ;******************************************************** begin ;------------------------------------------------------------------------- ; User Input ;------------------------------------------------------------------------- PLOT = True if (PLOT) then pltdir = "./" pltname = "CO_mop_map_coave2013" plttype = "x11" ; send graphics to X11 or png file ; Plotting definitions factor = 1.e18 units = "CO total column" ;max_contour = 360.0 ;min_contour = 0.0 ;spacing = 20.0 end if ;------------------------------------------------------------------------- ; End User Input ;------------------------------------------------------------------------- ;------------------------------------------------------------------------- ; Setup ;------------------------------------------------------------------------- ; mopitt input wkdir = "/glade/scratch/najiby/MOPITT_temp/" diri = "/glade/scratch/najiby/MOPITT_temp/" ; model input dir_mo = "/glade/work/najiby/f.e20.FCSD_DICE_WAFRICA/f.e20.FCSD_DICE_WAFRICA/atm/hist/" ;------------------------------------------------------------------------- ;-- read MOPITT data and define ;-------------------------------------------- ; find all MOP data: all_mop= systemfunc("ls -1 "+diri+"MOP03JM*") print(all_mop) nmop = dimsizes(all_mop) date_all = new(nmop,"string") ;-------------------------------------------- ; open all MOPITT files ;-------------------------------------------- fileIn = all_mop ;f = addfiles(fileIn, "r") ; ListSetType (f, "cat") ; concatenate (=default) ; Read files and collect variables f = addfiles(fileIn,"r") ListSetType (f, "join") ; concatenate (=default) ;Uncomment below to look at variable names ; cocol_retr_flip = lonFlip(cocol_retr) ; printVarSummary(cocol_retr_flip) ; cocol_retr_flip := cocol_retr_flip(:,YDim_MOP03|:,XDim_MOP03|:) ; printVarSummary( cocol_retr_flip ) ; cocol_retr_flip := cocol_retr_flip(:,::-1,:) ; printVarSummary( cocol_retr_flip ) vnames = getfilevarnames(f[0]) print(vnames) lat_mop = f[0]->YDim_MOP03 lat_mop@long_name = "Latitude" lon_mop = f[0]->XDim_MOP03 lon_mop@long_name = "Longitude" ; lat_mop = f[0]->Latitude_MOP03 ; lat_mop@long_name = "Latitude" ; lon_mop = f[0]->Longitude_MOP03 ; lon_mop@long_name = "Longitude" ; time_bound = f->time_bound ; prnit(time_bound) ; time = (/ time_bound(:,0) /) ; print(time) ; date = cd_calendar(time, 0) ; print(date) ; lat_mop = f[0]->YDim_MOP03 ; lon_mop = f[0]->XDim_MOP03 psurf = f[:]->SurfacePressureDay_MOP03 cocol_retr = f[:]->RetrievedCOTotalColumnDay_MOP03 cocol_a_priori = f[:]->APrioriCOTotalColumnDay_MOP03 cosurf_retr = f[:]->RetrievedCOSurfaceMixingRatioDay_MOP03 cosurf_a_priori = f[:]->APrioriCOSurfaceMixingRatioDay_MOP03 coprof_retr = f[:]->RetrievedCOMixingRatioProfileDay_MOP03 ; coprof_a_priori = f[:]->APrioriCOMixingRatioProfileDay_MOP03 ; avkernel = f[:]->TotalColumnAveragingKernelDay_MOP03 printVarSummary(cocol_retr) ; cocol_retr_flip = lonFlip(cocol_retr) ; printVarSummary(cocol_retr_flip) cocol_retr!0 = "time" cocol_retr_reorder := cocol_retr(time|:,YDim_MOP03|:,XDim_MOP03|:) printVarSummary(cocol_retr_reorder) cocol_retr_reorder_lonflip = lonFlip(cocol_retr_reorder) printVarSummary(cocol_retr_reorder_lonflip) ; cocol_retr_SN := cocol_retr_reorder(:,::-1,:) ; printVarSummary(cocol_retr_SN) ;------------------------------------------------------------------------- ; Load and collect equivalent Model data ;------------------------------------------------------------------------- dir_mo = "/glade/work/najiby/f.e20.FCSD_DICE_WAFRICA/f.e20.FCSD_DICE_WAFRICA/atm/hist/" file_mo = systemfunc("ls -1 "+dir_mo+"f.e20.FCSD_DICE_WAFRICA.cam.h0.2013*.nc") print(file_mo) filein = file_mo f = addfiles(filein, "r") ListSetType (f, "cat") ; concatenate (=default) ; Read model files and collect variables ;f = addfiles(filein,"r") ; ListSetType (f, "join") ; concatenate (=default) vnames1 = getfilevarnames(f[0]) print(vnames1) time =f[:]->time lev1 = f[:]->lev lat1 = f[:]->lat lon1 = f[:]->lon nlat1 = dimsizes(lat1) nlon1 = dimsizes(lon1) nlev1 = dimsizes(lev1) co_model = f[:]->CO printVarSummary(co_model) mod_yyyymmdd = cd_calendar(co_model&time, 2) model_lev = co_model&lev nlev1 = dimsizes(model_lev) ; load level information for interpolation ps = f[:]->PS hyam = f[0]->hyam hybm = f[0]->hybm p0 = f[0]->P0 mod_press = pres_hybrid_ccm(ps,p0,hyam,hybm) copy_VarCoords(co_model, mod_press) ; initialize array to collect model data co_model_select = new((/dimsizes(time), nlev1/),float) altp_model_select = new((/dimsizes(time), nlev1/),float) co_model_smooth = new(dimsizes(cocol_retr),float) ln_co_mo_ak = new(nlev1,float) ; natural log of a priori ln_cocol_a_priori = log(cocol_a_priori) ; MOPITT Prior MOPITT_ap_annavg = dim_avg_n(cocol_a_priori,0) ; Loop through Measurement data and find equivalent model data and Collect base model data fileIn = all_mop f = addfiles(fileIn,"r") ListSetType (f, "join") ; concatenate (=default) time_collect = f[:]@StartTime_MOP03 ;units value presumes use of TAI93 (International Atomic Time) format time_collect@units = "seconds since 1993-1-1 00:00:00" date_collect := cd_calendar(time_collect, 0) ; print(file_mo) ;-------------------------------------------- ;-- Collect equivalent profile ;-------------------------------------------- ; cocol_retr_flip = lonFlip(cocol_retr) ; printVarSummary(cocol_retr_flip) do n=0,nmop-1 time_ind = ind(mod_yyyymmdd.eq.time_collect(n)) co_model_select(n,:) = co_model(time_ind,:,{lat_mop(n)},{lon_mop(n)}) altp_model_select(n,:) = mod_press(time_ind,:,{lat_mop(n)},{lon_mop(n)}) end do ; Regrid to MOPITT grid co_mo_regrid = int2p_n_Wrap(altp_model_select(n,:)/100,co_model_select(n,:),psurf(n,:),0,0) ln_co_mo_regrid = log(co_mo_regrid) ; Read AK avkernel = f[:]->TotalColumnAveragingKernelDay_MOP03 coprof_a_priori = f[:]->APrioriCOMixingRatioProfileDay_MOP03 ;-- Apply averaging kernel ;-------------------------------------------- ; Need to use the usual multiplication because matrix multiplication can't handle missing values do m=0,nlev1-1 ln_co_mo_ak(m) = ln_cocol_a_priori(n,m) + sum((avkernel(n,m,:)) * (ln_co_mo_regrid - ln_cocol_a_priori(n,:))) end do co_mo_ak_dummy = exp(ln_co_mo_ak) co_mo_ak = where(ismissing(co_mo_regrid), co_mo_ak_dummy@_FillValue, co_mo_ak_dummy) ; Collect smoothed model data do n=0,nmop-1 co_model_smooth(n,:) = co_mo_ak end do ; end go through all MOPITT data profiles and create smoothed model value ;------------------------------------------------------------------------- ; Perform Annual Averages ;------------------------------------------------------------------------- ; MOPITT MOPITT_co_annavg = dim_avg_n(cocol_retr,0) MOPITT_co_annavg!0 = "lon" MOPITT_co_annavg&lon = lon_mop MOPITT_co_annavg!1 = "lat" MOPITT_co_annavg&lat = lat_mop MOPITT_co_annsd = dim_stddev_n(cocol_retr,0) MOPITT_alt_annavg = dim_avg_n(psurf,0) printVarSummary(MOPITT_co_annavg) printVarSummary(MOPITT_co_annsd) ; MODEL model_co_model_annavg = dim_avg_n(co_model,0) model_co_model_annsd = dim_stddev_n(co_model,0) model_altp_annavg = dim_avg_n(lev1,0)/100 printVarSummary(model_co_model_annavg) printVarSummary(model_co_model_annsd) ; convert pressure to meters convert_altp_annavg = stdatmus_p2tdz(model_altp_annavg) model_alt_annavg = convert_altp_annavg(2,:) ; Smoothed Model smooth_model_co_annavg = dim_avg_n(co_model_smooth,0) smooth_model_co_annsd = dim_stddev_n(co_model_smooth,0) if (PLOT) then ;-------------------------------------------- ;-- Plot ;-------------------------------------------- wks = gsn_open_wks(plttype, pltname) ;-------------- ;plot resources res = True res@gsnDraw = False res@gsnFrame = False res@gsnMaximize = False res@vpHeightF = 0.4 res@vpWidthF = 0.6 res@xyMarkLineMode = "Marklines" res@xyMarker = 16 ;-------------- ; limits res@trYMinF = 0. res@trYMaxF = 13000. res@trXMinF = 15. res@trXMaxF = 120. ;-------------- ; titles res@tiYAxisString = "Altitude (m)" res@tiXAxisString = "CO (ppb)" res@tiMainString = "2013 annual average of MOPITT CO compared with CAM-chem" res@tiMainFontHeightF = 0.018 res@gsnLeftString = "MOPITT retrievals over Lagos" ;-------------- ; legend res@pmLegendDisplayMode = "Always" ; turn on legend res@lgJustification = "BottomLeft" res@pmLegendParallelPosF = 0.62 ; move legend right res@pmLegendHeightF = 0.04 ; height of legend res@pmLegendWidthF = 0.08 ; width of legend res@lgBoxMinorExtentF = 0.4 ; Shorten the legend lines res@lgLabelFontHeightF = 0.012 ; label font height res@lgPerimOn = False ; turn off outline ;-------------- ;-- create plot ; MOPITT res@xyLineColor = "black" res@xyDashPattern = 2 res@xyLineThicknessF = 5.0 res@xyExplicitLegendLabels = "MOPITT" ; create explicit labels res@pmLegendOrthogonalPosF = -1.10 ; move legend down plot1 = gsn_csm_xy(wks,MOPITT_co_annavg/factor,model_altp_annavg,res) res@gsnLeftString = "" res@xyLineColor = "darkgreen" res@xyMarkerColor = "darkgreen" res@xyExplicitLegendLabels = "MOPITT a priori" ; create explicit labels res@pmLegendOrthogonalPosF = -1.05 ; move legend down plot1a = gsn_csm_xy(wks,MOPITT_ap_annavg/factor,MOPITT_alt_annavg,res) overlay(plot1, plot1a) ; Model res@xyLineColor = "blue" res@xyMarkerColor = "blue" res@xyDashPattern = 14 res@xyExplicitLegendLabels = "CAM-chem" ; create explicit labels res@pmLegendOrthogonalPosF = -1.00 ; move legend down plot2 = gsn_csm_xy(wks,model_co_model_annavg/factor,MOPITT_alt_annavg,res) ; Model smoothed res@xyLineColor = "red" res@xyMarkerColor = "red" res@xyDashPattern = 0 res@xyExplicitLegendLabels = "CAM-chem with TES AK" ; create explicit labels res@pmLegendOrthogonalPosF = -0.95 ; move legend down plot3 = gsn_csm_xy(wks,smooth_model_co_annavg/factor,MOPITT_alt_annavg,res) ;-------------- ;-- add standard ;-- deviation gsres = True ; poly res gsres@tfPolyDrawOrder = "Predraw" ; draw this first gsres@gsFillOpacityF = 0.2 ; MOPITT nalt = dimsizes(model_altp_annavg) xp = new( (/2*nalt/), float ) yp = new( (/2*nalt/), float ) do k=0,nalt-1 dx = MOPITT_co_annsd(k)/factor yp(k) = MOPITT_co_annavg(k)/factor + dx xp(k) = MOPITT_alt_annavg(k) xp(2*nalt-1-k) = MOPITT_alt_annavg(k) yp(2*nalt-1-k) = MOPITT_co_annavg(k)/factor - dx end do gsres@gsFillColor = "Grey" ; color chosen dummy = gsn_add_polygon (wks,plot1,yp,xp,gsres) ; model m_nalt = dimsizes(model_alt_annavg) m_xp = new( (/2*m_nalt/), float ) m_yp = new( (/2*m_nalt/), float ) do k=0,m_nalt-1 dx = model_co_model_annsd(k)/factor m_yp(k) = model_co_model_annavg(k)/factor + dx m_xp(k) = model_alt_annavg(k) m_xp(2*m_nalt-1-k)= model_alt_annavg(k) m_yp(2*m_nalt-1-k)= model_co_model_annavg(k)/factor - dx end do gsres@gsFillColor = "blue" ; color chosen dummy2 = gsn_add_polygon (wks,plot2,m_yp,m_xp,gsres) ; smoothed model sm_xp = new( (/2*nalt/), float ) sm_yp = new( (/2*nalt/), float ) do k=0,nalt-1 dx = smooth_model_co_annsd(k)/factor sm_yp(k) = smooth_model_co_annavg(k)/factor + dx sm_xp(k) = MOPITT_alt_annavg(k) sm_xp(2*nalt-1-k)= MOPITT_alt_annavg(k) sm_yp(2*nalt-1-k)= smooth_model_co_annavg(k)/factor - dx end do gsres@gsFillColor = "red" ; color chosen dummy3 = gsn_add_polygon (wks,plot3,sm_yp,sm_xp,gsres) overlay(plot1, plot2) overlay(plot1, plot3) draw(plot1) frame(wks) end if end