diff -r 000000000000 -r 0c6405ab2ff4 fire/21.table+tseries.ncl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fire/21.table+tseries.ncl Mon Jan 26 22:08:20 2009 -0500 @@ -0,0 +1,258 @@ +;******************************************************** +;using model biome vlass +; +; required command line input parameters: +; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl +; +; histogram normalized by rain and compute correleration +;************************************************************** +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" +;************************************************************** +procedure set_line(lines:string,nline:integer,newlines:string) +begin +; add line to ascci/html file + + nnewlines = dimsizes(newlines) + if(nline+nnewlines-1.ge.dimsizes(lines)) + print("set_line: bad index, not setting anything.") + return + end if + lines(nline:nline+nnewlines-1) = newlines +; print ("lines = " + lines(nline:nline+nnewlines-1)) + nline = nline + nnewlines + return +end +;************************************************************** +; Main code. +begin + + plot_type = "ps" + plot_type_new = "png" + +;--------------------------------------------------- +; model name and grid + + model_grid = "T42" + + model_name = "cn" + model_name1 = "i01.06cn" + model_name2 = "i01.10cn" + +;--------------------------------------------------- +; get biome data: model + + biome_name_mod = "Model PFT Class" + + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" + film = "class_pft_"+model_grid+".nc" + fm = addfile(dirm+film,"r") + + classmod = fm->CLASS_PFT + + delete (fm) + +; model data has 17 land-type classes + + nclass_mod = 17 + +;-------------------------------------------------- +; get model data: landfrac and area + + dirm = "/fis/cgd/cseg/people/jeff/surface_data/" + film = "lnd_T42.nc" + fm = addfile (dirm+film,"r") + + landmask = fm->landmask + landfrac = fm->landfrac + area = fm->area + + delete (fm) + +; change area from km**2 to m**2 + area = area * 1.e6 + +;---------------------------------------------------- +; read data: time series, model + + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" + film = model_name2 + "_Fire_C_1979-2004_monthly.nc" + fm = addfile (dirm+film,"r") + + data_mod = fm->COL_FIRE_CLOSS(18:25,:,:,:) + + delete (fm) + +; Units for these variables are: +; g C/m^2/s + +; change unit to g C/m^2/month + + nsec_per_month = 60*60*24*30 + + data_mod = data_mod * nsec_per_month + + data_mod@unit = "gC/m2/month" +;---------------------------------------------------- +; read data: time series, observed + + dirm = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/" + film = "Fire_C_1997-2006_monthly_"+ model_grid+".nc" + fm = addfile (dirm+film,"r") + + data_ob = fm->FIRE_C(0:7,:,:,:) + + delete (fm) + + ob_name = "GFEDv2" + +; Units for these variables are: +; g C/m^2/month + +;--------------------------------------------------- +; take into account landfrac + + data_mod = data_mod * conform(data_mod, landfrac, (/2,3/)) + data_ob = data_ob * conform(data_ob, landfrac, (/2,3/)) + +;--------------------------------------------------- +; get time-mean + + x = dim_avg_Wrap(data_mod(lat|:,lon|:,month|:,year|:)) + data_mod_m = dim_avg_Wrap( x(lat|:,lon|:,month|:)) + delete (x) + + x = dim_avg_Wrap( data_ob(lat|:,lon|:,month|:,year|:)) + data_ob_m = dim_avg_Wrap( x(lat|:,lon|:,month|:)) + delete (x) + +; printVarSummary(y) + +;---------------------------------------------------- +; compute correlation coef + + landmask_1d = ndtooned(landmask) + data_mod_1d = ndtooned(data_mod_m) + data_ob_1d = ndtooned(data_ob_m) + + good = ind(landmask_1d .gt. 0.) +; print (dimsizes(good)) + + cc = esccr(data_mod_1d(good),data_ob_1d(good),0) +; print (cc) + + delete (landmask_1d) + delete (data_mod_1d) + delete (data_ob_id) + +;---------------------------------------------------- +; compute M_global + + score_max = 1. + + Mscore = cc * cc * score_max + + M_global = sprintf("%.2f", Mscore) + +;---------------------------------------------------- +; global res + + resg = True ; Use plot options + resg@cnFillOn = True ; Turn on color fill + resg@gsnSpreadColors = True ; use full colormap + resg@cnLinesOn = False ; Turn off contourn lines + resg@mpFillOn = False ; Turn off map fill + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals + +;---------------------------------------------------- +; global contour: model vs ob + + plot_name = "global_model_vs_ob" + + wks = gsn_open_wks (plot_type,plot_name) + gsn_define_colormap(wks,"gui_default") + + plot=new(3,graphic) ; create graphic array + + resg@gsnFrame = False ; Do not draw plot + resg@gsnDraw = False ; Do not advance frame + +;---------------------- +; plot correlation coef + + gRes = True + gRes@txFontHeightF = 0.02 + gRes@txAngleF = 90 + + correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")" + + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) + +;----------------------- +; plot ob + + data_ob_m = where(landmask .gt. 0., data_ob_m, data_ob_m@_FillValue) + + title = ob_name + resg@tiMainString = title + + resg@cnMinLevelValF = 1. + resg@cnMaxLevelValF = 10. + resg@cnLevelSpacingF = 1. + + plot(0) = gsn_csm_contour_map_ce(wks,data_ob_m,resg) + +;----------------------- +; plot model + + data_mod_m = where(landmask .gt. 0., data_mod_m, data_mod_m@_FillValue) + + title = "Model "+ model_name + resg@tiMainString = title + + resg@cnMinLevelValF = 1. + resg@cnMaxLevelValF = 10. + resg@cnLevelSpacingF = 1. + + plot(1) = gsn_csm_contour_map_ce(wks,data_mod_m,resg) + +;----------------------- +; plot model-ob + + resg@cnMinLevelValF = -8. + resg@cnMaxLevelValF = 2. + resg@cnLevelSpacingF = 1. + + zz = data_ob_m + zz = data_mod_m - data_ob_m + title = "Model_"+model_name+" - Observed" + resg@tiMainString = title + + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) + +; plot panel + + pres = True ; panel plot mods desired + pres@gsnMaximize = True ; fill the page + + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot + +exit + +; system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ +; "rm "+plot_name+"."+plot_type) + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) + + clear (wks) + delete (plot) + + delete (data_ob_m) + delete (data_mod_m) + delete (zz) + + resg@gsnFrame = True ; Do advance frame + resg@gsnDraw = True ; Do draw plot + +end +