forrest@0: ;******************************************************** forrest@0: ;using model biome vlass forrest@0: ; forrest@0: ; required command line input parameters: forrest@0: ; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl forrest@0: ; forrest@0: ; histogram normalized by rain and compute correleration forrest@0: ;************************************************************** forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" forrest@0: ;************************************************************** forrest@0: procedure set_line(lines:string,nline:integer,newlines:string) forrest@0: begin forrest@0: ; add line to ascci/html file forrest@0: forrest@0: nnewlines = dimsizes(newlines) forrest@0: if(nline+nnewlines-1.ge.dimsizes(lines)) forrest@0: print("set_line: bad index, not setting anything.") forrest@0: return forrest@0: end if forrest@0: lines(nline:nline+nnewlines-1) = newlines forrest@0: ; print ("lines = " + lines(nline:nline+nnewlines-1)) forrest@0: nline = nline + nnewlines forrest@0: return forrest@0: end forrest@0: ;************************************************************** forrest@0: ; Main code. forrest@0: begin forrest@0: forrest@0: plot_type = "ps" forrest@0: plot_type_new = "png" forrest@0: forrest@0: ;--------------------------------------------------- forrest@0: ; model name and grid forrest@0: forrest@0: model_grid = "T42" forrest@0: forrest@0: model_name = "cn" forrest@0: model_name1 = "i01.06cn" forrest@0: model_name2 = "i01.10cn" forrest@0: forrest@0: ;--------------------------------------------------- forrest@0: ; get biome data: model forrest@0: forrest@0: biome_name_mod = "Model PFT Class" forrest@0: forrest@0: dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" forrest@0: film = "class_pft_"+model_grid+".nc" forrest@0: fm = addfile(dirm+film,"r") forrest@0: forrest@0: classmod = fm->CLASS_PFT forrest@0: forrest@0: delete (fm) forrest@0: forrest@0: ; model data has 17 land-type classes forrest@0: forrest@0: nclass_mod = 17 forrest@0: forrest@0: ;-------------------------------------------------- forrest@0: ; get model data: landmask, landfrac and area forrest@0: forrest@0: dirm = "/fis/cgd/cseg/people/jeff/surface_data/" forrest@0: film = "lnd_T42.nc" forrest@0: fm = addfile (dirm+film,"r") forrest@0: forrest@0: landmask = fm->landmask forrest@0: landfrac = fm->landfrac forrest@0: area = fm->area forrest@0: forrest@0: delete (fm) forrest@0: forrest@0: ; change area from km**2 to m**2 forrest@0: area = area * 1.e6 forrest@0: forrest@0: ;---------------------------------------------------- forrest@0: ; read data: time series, model forrest@0: forrest@0: dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" forrest@0: film = model_name2 + "_Fire_C_1979-2004_monthly.nc" forrest@0: fm = addfile (dirm+film,"r") forrest@0: forrest@0: data_mod = fm->COL_FIRE_CLOSS(18:25,:,:,:) forrest@0: forrest@0: delete (fm) forrest@0: forrest@0: ; Units for these variables are: forrest@0: ; g C/m^2/s forrest@0: forrest@0: ; change unit to g C/m^2/month forrest@0: forrest@0: nsec_per_month = 60*60*24*30 forrest@0: forrest@0: data_mod = data_mod * nsec_per_month forrest@0: forrest@0: data_mod@unit = "gC/m2/month" forrest@0: ;---------------------------------------------------- forrest@0: ; read data: time series, observed forrest@0: forrest@0: dirm = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/" forrest@0: film = "Fire_C_1997-2006_monthly_"+ model_grid+".nc" forrest@0: fm = addfile (dirm+film,"r") forrest@0: forrest@0: data_ob = fm->FIRE_C(0:7,:,:,:) forrest@0: forrest@0: delete (fm) forrest@0: forrest@0: ob_name = "GFEDv2" forrest@0: forrest@0: ; Units for these variables are: forrest@0: ; g C/m^2/month forrest@0: forrest@0: ;--------------------------------------------------- forrest@0: ; take into account landfrac forrest@0: forrest@0: area = area * landfrac forrest@0: data_mod = data_mod * conform(data_mod, landfrac, (/2,3/)) forrest@0: data_ob = data_ob * conform(data_ob, landfrac, (/2,3/)) forrest@0: forrest@0: delete (landfrac) forrest@0: forrest@0: ;------------------------------------------------------------- forrest@0: ; html table1 data forrest@0: forrest@0: ; column (not including header column) forrest@0: forrest@0: col_head = (/"Model Fire_Flux (PgC/yr)" \ forrest@0: ,"Observed Fire_Flux (PgC/yr)" \ forrest@0: ,"Correlation Coefficient" \ forrest@0: ,"Ratio model/observed" \ forrest@0: ,"M_score" \ forrest@0: ,"Timeseries plot" \ forrest@0: /) forrest@0: forrest@0: ncol = dimsizes(col_head) forrest@0: forrest@0: ; row (not including header row) forrest@0: forrest@0: ; using model biome class: forrest@0: row_head = (/"Not Vegetated" \ forrest@0: ,"Needleleaf Evergreen Temperate Tree" \ forrest@0: ,"Needleleaf Evergreen Boreal Tree" \ forrest@0: ; ,"Needleleaf Deciduous Boreal Tree" \ forrest@0: ,"Broadleaf Evergreen Tropical Tree" \ forrest@0: ,"Broadleaf Evergreen Temperate Tree" \ forrest@0: ,"Broadleaf Deciduous Tropical Tree" \ forrest@0: ,"Broadleaf Deciduous Temperate Tree" \ forrest@0: ; ,"Broadleaf Deciduous Boreal Tree" \ forrest@0: ; ,"Broadleaf Evergreen Shrub" \ forrest@0: ,"Broadleaf Deciduous Temperate Shrub" \ forrest@0: ,"Broadleaf Deciduous Boreal Shrub" \ forrest@0: ,"C3 Arctic Grass" \ forrest@0: ,"C3 Non-Arctic Grass" \ forrest@0: ,"C4 Grass" \ forrest@0: ,"Corn" \ forrest@0: ; ,"Wheat" \ forrest@0: ,"All Biome" \ forrest@0: /) forrest@0: nrow = dimsizes(row_head) forrest@0: forrest@0: ; arrays to be passed to table. forrest@0: text = new ((/nrow, ncol/),string ) forrest@0: forrest@0: ;***************************************************************** forrest@0: ; (A) get time-mean forrest@0: ;***************************************************************** forrest@0: forrest@0: x = dim_avg_Wrap(data_mod(lat|:,lon|:,month|:,year|:)) forrest@0: data_mod_m = dim_avg_Wrap( x(lat|:,lon|:,month|:)) forrest@0: delete (x) forrest@0: forrest@0: x = dim_avg_Wrap( data_ob(lat|:,lon|:,month|:,year|:)) forrest@0: data_ob_m = dim_avg_Wrap( x(lat|:,lon|:,month|:)) forrest@0: delete (x) forrest@0: forrest@0: ;---------------------------------------------------- forrest@0: ; compute correlation coef forrest@0: forrest@0: landmask_1d = ndtooned(landmask) forrest@0: data_mod_1d = ndtooned(data_mod_m) forrest@0: data_ob_1d = ndtooned(data_ob_m) forrest@0: forrest@0: good = ind(landmask_1d .gt. 0.) forrest@0: ; print (dimsizes(good)) forrest@0: forrest@0: cc = esccr(data_mod_1d(good),data_ob_1d(good),0) forrest@0: ; print (cc) forrest@0: forrest@0: delete (landmask_1d) forrest@0: delete (data_mod_1d) forrest@0: delete (data_ob_1d) forrest@0: delete (good) forrest@0: forrest@0: ;---------------------------------------------------- forrest@0: ; compute M_global forrest@0: forrest@0: score_max = 1. forrest@0: forrest@0: Mscore = cc * cc * score_max forrest@0: forrest@0: M_global = sprintf("%.2f", Mscore) forrest@0: forrest@0: ;---------------------------------------------------- forrest@0: ; global res forrest@0: forrest@0: resg = True ; Use plot options forrest@0: resg@cnFillOn = True ; Turn on color fill forrest@0: resg@gsnSpreadColors = True ; use full colormap forrest@0: resg@cnLinesOn = False ; Turn off contourn lines forrest@0: resg@mpFillOn = False ; Turn off map fill forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals forrest@0: forrest@0: ;---------------------------------------------------- forrest@0: ; global contour: model vs ob forrest@0: forrest@0: plot_name = "global_model_vs_ob" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: gsn_define_colormap(wks,"gui_default") forrest@0: forrest@0: plot=new(3,graphic) ; create graphic array forrest@0: forrest@0: resg@gsnFrame = False ; Do not draw plot forrest@0: resg@gsnDraw = False ; Do not advance frame forrest@0: forrest@0: ;---------------------- forrest@0: ; plot correlation coef forrest@0: forrest@0: gRes = True forrest@0: gRes@txFontHeightF = 0.02 forrest@0: gRes@txAngleF = 90 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) forrest@0: forrest@0: ;----------------------- forrest@0: ; plot ob forrest@0: forrest@0: data_ob_m = where(landmask .gt. 0., data_ob_m, data_ob_m@_FillValue) forrest@0: forrest@0: title = ob_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: resg@cnMinLevelValF = 1. forrest@0: resg@cnMaxLevelValF = 10. forrest@0: resg@cnLevelSpacingF = 1. forrest@0: forrest@0: plot(0) = gsn_csm_contour_map_ce(wks,data_ob_m,resg) forrest@0: forrest@0: ;----------------------- forrest@0: ; plot model forrest@0: forrest@0: data_mod_m = where(landmask .gt. 0., data_mod_m, data_mod_m@_FillValue) forrest@0: forrest@0: title = "Model "+ model_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: resg@cnMinLevelValF = 1. forrest@0: resg@cnMaxLevelValF = 10. forrest@0: resg@cnLevelSpacingF = 1. forrest@0: forrest@0: plot(1) = gsn_csm_contour_map_ce(wks,data_mod_m,resg) forrest@0: forrest@0: ;----------------------- forrest@0: ; plot model-ob forrest@0: forrest@0: resg@cnMinLevelValF = -8. forrest@0: resg@cnMaxLevelValF = 2. forrest@0: resg@cnLevelSpacingF = 1. forrest@0: forrest@0: zz = data_ob_m forrest@0: zz = data_mod_m - data_ob_m forrest@0: title = "Model_"+model_name+" - Observed" forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) forrest@0: forrest@0: ; plot panel forrest@0: forrest@0: pres = True ; panel plot mods desired forrest@0: pres@gsnMaximize = True ; fill the page forrest@0: forrest@0: gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot forrest@0: forrest@0: ; system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: ; "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: delete (plot) forrest@0: forrest@0: delete (data_ob_m) forrest@0: delete (data_mod_m) forrest@0: delete (zz) forrest@0: forrest@0: resg@gsnFrame = True ; Do advance frame forrest@0: resg@gsnDraw = True ; Do draw plot forrest@0: forrest@0: ;******************************************************************* forrest@0: ; (B) Time series : per biome forrest@0: ;******************************************************************* forrest@0: forrest@0: data_n = 2 forrest@0: forrest@0: dsizes = dimsizes(data_mod) forrest@0: nyear = dsizes(0) forrest@0: nmonth = dsizes(1) forrest@0: ntime = nyear * nmonth forrest@0: forrest@0: year_start = 1997 forrest@0: year_end = 2004 forrest@0: forrest@0: ;------------------------------------------- forrest@0: ; Calculate "nice" bins for binning the data forrest@0: forrest@0: ; using model biome class forrest@0: nclass = nclass_mod forrest@0: forrest@0: range = fspan(0,nclass,nclass+1) forrest@0: forrest@0: ; print (range) forrest@0: ; Use this range information to grab all the values in a forrest@0: ; particular range, and then take an average. forrest@0: forrest@0: nx = dimsizes(range) - 1 forrest@0: forrest@0: ;------------------------------------------- forrest@0: ; put data into bins forrest@0: forrest@0: ; using observed biome class forrest@0: ; base = ndtooned(classob) forrest@0: ; using model biome class forrest@0: base = ndtooned(classmod) forrest@0: forrest@0: ; output forrest@0: forrest@0: area_bin = new((/nx/),float) forrest@0: yvalues = new((/ntime,data_n,nx/),float) forrest@0: forrest@0: ; Loop through each range, using base. forrest@0: forrest@0: do i=0,nx-1 forrest@0: forrest@0: if (i.ne.(nx-1)) then forrest@0: idx = ind((base.ge.range(i)).and.(base.lt.range(i+1))) forrest@0: else forrest@0: idx = ind(base.ge.range(i)) forrest@0: end if forrest@0: ;--------------------- forrest@0: ; for area forrest@0: forrest@0: data = ndtooned(area) forrest@0: forrest@0: if (.not.any(ismissing(idx))) then forrest@0: area_bin(i) = sum(data(idx)) forrest@0: else forrest@0: area_bin(i) = area_bin@_FillValue forrest@0: end if forrest@0: forrest@0: ;############################################################# forrest@0: ; using model biome class: forrest@0: ; set the following 4 classes to _FillValue: forrest@0: ; (3)Needleleaf Deciduous Boreal Tree, forrest@0: ; (8)Broadleaf Deciduous Boreal Tree, forrest@0: ; (9)Broadleaf Evergreen Shrub, forrest@0: ; (16)Wheat forrest@0: forrest@0: if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then forrest@0: area_bin(i) = area_bin@_FillValue forrest@0: end if forrest@0: ;############################################################# forrest@0: forrest@0: delete (data) forrest@0: forrest@0: ;--------------------- forrest@0: ; for data_mod and data_ob forrest@0: forrest@0: do n = 0,data_n-1 forrest@0: forrest@0: t = -1 forrest@0: do m = 0,nyear-1 forrest@0: do k = 0,nmonth-1 forrest@0: forrest@0: t = t + 1 forrest@0: forrest@0: if (n.eq.0) then forrest@0: data = ndtooned(data_ob(m,k,:,:)) forrest@0: end if forrest@0: forrest@0: if (n.eq.1) then forrest@0: data = ndtooned(data_mod(m,k,:,:)) forrest@0: end if forrest@0: forrest@0: ; Calculate average forrest@0: forrest@0: if (.not.any(ismissing(idx))) then forrest@0: yvalues(t,n,i) = avg(data(idx)) forrest@0: else forrest@0: yvalues(t,n,i) = yvalues@_FillValue forrest@0: end if forrest@0: forrest@0: ;############################################################# forrest@0: ; using model biome class: forrest@0: ; set the following 4 classes to _FillValue: forrest@0: ; (3)Needleleaf Deciduous Boreal Tree, forrest@0: ; (8)Broadleaf Deciduous Boreal Tree, forrest@0: ; (9)Broadleaf Evergreen Shrub, forrest@0: ; (16)Wheat forrest@0: forrest@0: if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then forrest@0: yvalues(t,n,i) = yvalues@_FillValue forrest@0: end if forrest@0: ;############################################################# forrest@0: forrest@0: end do forrest@0: end do forrest@0: forrest@0: delete(data) forrest@0: end do forrest@0: forrest@0: delete(idx) forrest@0: end do forrest@0: forrest@0: delete (base) forrest@0: delete (data_mod) forrest@0: delete (data_ob) forrest@0: forrest@0: ;---------------------------------------------------------------- forrest@0: ; get area_good forrest@0: forrest@0: good = ind(.not.ismissing(area_bin)) forrest@0: forrest@0: area_g = area_bin(good) forrest@0: forrest@0: n_biome = dimsizes(good) forrest@0: forrest@0: ;---------------------------------------------------------------- forrest@0: ; data for tseries plot forrest@0: forrest@0: yvalues_g = new((/ntime,data_n,n_biome/),float) forrest@0: forrest@0: yvalues_g@units = "TgC/month" forrest@0: forrest@0: ; change unit to Tg C/month forrest@0: ; change unit from g to Tg (Tera gram) forrest@0: factor_unit = 1.e-12 forrest@0: forrest@0: yvalues_g = yvalues(:,:,good) * conform(yvalues_g,area_g,2) * factor_unit forrest@0: forrest@0: ;------------------------------------------------------------------- forrest@0: ; general settings for line plot forrest@0: forrest@0: res = True forrest@0: res@xyDashPatterns = (/0,0/) ; make lines solid forrest@0: res@xyLineThicknesses = (/2.0,2.0/) ; make lines thicker forrest@0: res@xyLineColors = (/"blue","red"/) ; line color forrest@0: forrest@0: res@trXMinF = year_start forrest@0: res@trXMaxF = year_end + 1 forrest@0: forrest@0: res@vpHeightF = 0.4 ; change aspect ratio of plot forrest@0: ; res@vpWidthF = 0.8 forrest@0: res@vpWidthF = 0.75 forrest@0: forrest@0: res@tiMainFontHeightF = 0.025 ; size of title forrest@0: forrest@0: res@tmXBFormat = "f" ; not to add trailing zeros forrest@0: forrest@0: ; res@gsnMaximize = True forrest@0: forrest@0: ;---------------------------------------------- forrest@0: ; Add a boxed legend using the simple method forrest@0: forrest@0: res@pmLegendDisplayMode = "Always" forrest@0: ; res@pmLegendWidthF = 0.1 forrest@0: res@pmLegendWidthF = 0.08 forrest@0: res@pmLegendHeightF = 0.06 forrest@0: res@pmLegendOrthogonalPosF = -1.17 forrest@0: ; res@pmLegendOrthogonalPosF = -1.00 ;(downward) forrest@0: ; res@pmLegendOrthogonalPosF = -0.30 ;(downward) forrest@0: forrest@0: ; res@pmLegendParallelPosF = 0.18 forrest@0: res@pmLegendParallelPosF = 0.23 ;(rightward) forrest@0: res@pmLegendParallelPosF = 0.73 ;(rightward) forrest@0: res@pmLegendParallelPosF = 0.83 ;(rightward) forrest@0: forrest@0: ; res@lgPerimOn = False forrest@0: res@lgLabelFontHeightF = 0.015 forrest@0: res@xyExplicitLegendLabels = (/"observed",model_name/) forrest@0: forrest@0: ;******************************************************************* forrest@0: ; (A) time series plot: monthly ( 2 lines per plot) forrest@0: ;******************************************************************* forrest@0: forrest@0: ; x-axis in time series plot forrest@0: forrest@0: timeI = new((/ntime/),integer) forrest@0: timeF = new((/ntime/),float) forrest@0: timeI = ispan(1,ntime,1) forrest@0: timeF = year_start + (timeI-1)/12. forrest@0: timeF@long_name = "year" forrest@0: forrest@0: plot_data = new((/2,ntime/),float) forrest@0: plot_data@long_name = "TgC/month" forrest@0: forrest@0: ;---------------------------------------------- forrest@0: ; time series : per biome forrest@0: forrest@0: do m = 0, n_biome-1 forrest@0: forrest@0: plot_name = "monthly_biome_"+ m forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: forrest@0: title = "Fire : "+ row_head(m) forrest@0: res@tiMainString = title forrest@0: forrest@0: plot_data(0,:) = yvalues_g(:,0,m) forrest@0: plot_data(1,:) = yvalues_g(:,1,m) forrest@0: forrest@0: plot = gsn_csm_xy(wks,timeF,plot_data,res) forrest@0: forrest@0: ; system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: ; "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: delete (plot) forrest@0: forrest@0: end do forrest@0: forrest@0: ;-------------------------------------------- forrest@0: ; time series: global forrest@0: forrest@0: plot_name = "monthly_global" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: forrest@0: title = "Fire : "+ row_head(n_biome) forrest@0: res@tiMainString = title forrest@0: forrest@0: do k = 0,ntime-1 forrest@0: plot_data(0,k) = sum(yvalues_g(k,0,:)) forrest@0: plot_data(1,k) = sum(yvalues_g(k,1,:)) forrest@0: end do forrest@0: forrest@0: plot = gsn_csm_xy(wks,timeF,plot_data,res) forrest@0: forrest@0: ; system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: ; "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: clear (wks) forrest@0: delete (plot) forrest@0: end forrest@0: