1 ;********************************************************
2 ;using model biome vlass
4 ; required command line input parameters:
5 ; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
7 ; histogram normalized by rain and compute correleration
8 ;**************************************************************
9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
11 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
12 ;**************************************************************
13 procedure set_line(lines:string,nline:integer,newlines:string)
15 ; add line to ascci/html file
17 nnewlines = dimsizes(newlines)
18 if(nline+nnewlines-1.ge.dimsizes(lines))
19 print("set_line: bad index, not setting anything.")
22 lines(nline:nline+nnewlines-1) = newlines
23 ; print ("lines = " + lines(nline:nline+nnewlines-1))
24 nline = nline + nnewlines
27 ;**************************************************************
34 ;---------------------------------------------------
40 model_name1 = "i01.06cn"
41 model_name2 = "i01.10cn"
43 ;---------------------------------------------------
44 ; get biome data: model
46 biome_name_mod = "Model PFT Class"
48 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
49 film = "class_pft_"+model_grid+".nc"
50 fm = addfile(dirm+film,"r")
52 classmod = fm->CLASS_PFT
56 ; model data has 17 land-type classes
60 ;--------------------------------------------------
61 ; get model data: landmask, landfrac and area
63 dirm = "/fis/cgd/cseg/people/jeff/surface_data/"
65 fm = addfile (dirm+film,"r")
67 landmask = fm->landmask
68 landfrac = fm->landfrac
73 ; change area from km**2 to m**2
76 ;----------------------------------------------------
77 ; read data: time series, model
79 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
80 film = model_name2 + "_Fire_C_1979-2004_monthly.nc"
81 fm = addfile (dirm+film,"r")
83 data_mod = fm->COL_FIRE_CLOSS(18:25,:,:,:)
87 ; Units for these variables are:
90 ; change unit to g C/m^2/month
92 nsec_per_month = 60*60*24*30
94 data_mod = data_mod * nsec_per_month
96 data_mod@unit = "gC/m2/month"
97 ;----------------------------------------------------
98 ; read data: time series, observed
100 dirm = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"
101 film = "Fire_C_1997-2006_monthly_"+ model_grid+".nc"
102 fm = addfile (dirm+film,"r")
104 data_ob = fm->FIRE_C(0:7,:,:,:)
110 ; Units for these variables are:
113 ;---------------------------------------------------
114 ; take into account landfrac
116 area = area * landfrac
117 data_mod = data_mod * conform(data_mod, landfrac, (/2,3/))
118 data_ob = data_ob * conform(data_ob, landfrac, (/2,3/))
122 ;-------------------------------------------------------------
125 ; column (not including header column)
127 col_head = (/"Model Fire_Flux (PgC/yr)" \
128 ,"Observed Fire_Flux (PgC/yr)" \
129 ,"Correlation Coefficient" \
130 ,"Ratio model/observed" \
135 ncol = dimsizes(col_head)
137 ; row (not including header row)
139 ; using model biome class:
140 row_head = (/"Not Vegetated" \
141 ,"Needleleaf Evergreen Temperate Tree" \
142 ,"Needleleaf Evergreen Boreal Tree" \
143 ; ,"Needleleaf Deciduous Boreal Tree" \
144 ,"Broadleaf Evergreen Tropical Tree" \
145 ,"Broadleaf Evergreen Temperate Tree" \
146 ,"Broadleaf Deciduous Tropical Tree" \
147 ,"Broadleaf Deciduous Temperate Tree" \
148 ; ,"Broadleaf Deciduous Boreal Tree" \
149 ; ,"Broadleaf Evergreen Shrub" \
150 ,"Broadleaf Deciduous Temperate Shrub" \
151 ,"Broadleaf Deciduous Boreal Shrub" \
153 ,"C3 Non-Arctic Grass" \
159 nrow = dimsizes(row_head)
161 ; arrays to be passed to table.
162 text = new ((/nrow, ncol/),string )
164 ;*****************************************************************
166 ;*****************************************************************
168 x = dim_avg_Wrap(data_mod(lat|:,lon|:,month|:,year|:))
169 data_mod_m = dim_avg_Wrap( x(lat|:,lon|:,month|:))
172 x = dim_avg_Wrap( data_ob(lat|:,lon|:,month|:,year|:))
173 data_ob_m = dim_avg_Wrap( x(lat|:,lon|:,month|:))
176 ;----------------------------------------------------
177 ; compute correlation coef
179 landmask_1d = ndtooned(landmask)
180 data_mod_1d = ndtooned(data_mod_m)
181 data_ob_1d = ndtooned(data_ob_m)
183 good = ind(landmask_1d .gt. 0.)
184 ; print (dimsizes(good))
186 cc = esccr(data_mod_1d(good),data_ob_1d(good),0)
194 ;----------------------------------------------------
199 Mscore = cc * cc * score_max
201 M_global = sprintf("%.2f", Mscore)
203 ;----------------------------------------------------
206 resg = True ; Use plot options
207 resg@cnFillOn = True ; Turn on color fill
208 resg@gsnSpreadColors = True ; use full colormap
209 resg@cnLinesOn = False ; Turn off contourn lines
210 resg@mpFillOn = False ; Turn off map fill
211 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
213 ;----------------------------------------------------
214 ; global contour: model vs ob
216 plot_name = "global_model_vs_ob"
218 wks = gsn_open_wks (plot_type,plot_name)
219 gsn_define_colormap(wks,"gui_default")
221 plot=new(3,graphic) ; create graphic array
223 resg@gsnFrame = False ; Do not draw plot
224 resg@gsnDraw = False ; Do not advance frame
226 ;----------------------
227 ; plot correlation coef
230 gRes@txFontHeightF = 0.02
233 correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")"
235 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
237 ;-----------------------
240 data_ob_m = where(landmask .gt. 0., data_ob_m, data_ob_m@_FillValue)
243 resg@tiMainString = title
245 resg@cnMinLevelValF = 1.
246 resg@cnMaxLevelValF = 10.
247 resg@cnLevelSpacingF = 1.
249 plot(0) = gsn_csm_contour_map_ce(wks,data_ob_m,resg)
251 ;-----------------------
254 data_mod_m = where(landmask .gt. 0., data_mod_m, data_mod_m@_FillValue)
256 title = "Model "+ model_name
257 resg@tiMainString = title
259 resg@cnMinLevelValF = 1.
260 resg@cnMaxLevelValF = 10.
261 resg@cnLevelSpacingF = 1.
263 plot(1) = gsn_csm_contour_map_ce(wks,data_mod_m,resg)
265 ;-----------------------
268 resg@cnMinLevelValF = -8.
269 resg@cnMaxLevelValF = 2.
270 resg@cnLevelSpacingF = 1.
273 zz = data_mod_m - data_ob_m
274 title = "Model_"+model_name+" - Observed"
275 resg@tiMainString = title
277 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
281 pres = True ; panel plot mods desired
282 pres@gsnMaximize = True ; fill the page
284 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
286 ; system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
287 ; "rm "+plot_name+"."+plot_type)
296 resg@gsnFrame = True ; Do advance frame
297 resg@gsnDraw = True ; Do draw plot
299 ;*******************************************************************
300 ; (B) Time series : per biome
301 ;*******************************************************************
305 dsizes = dimsizes(data_mod)
308 ntime = nyear * nmonth
313 ;-------------------------------------------
314 ; Calculate "nice" bins for binning the data
316 ; using model biome class
319 range = fspan(0,nclass,nclass+1)
322 ; Use this range information to grab all the values in a
323 ; particular range, and then take an average.
325 nx = dimsizes(range) - 1
327 ;-------------------------------------------
330 ; using observed biome class
331 ; base = ndtooned(classob)
332 ; using model biome class
333 base = ndtooned(classmod)
337 area_bin = new((/nx/),float)
338 yvalues = new((/ntime,data_n,nx/),float)
340 ; Loop through each range, using base.
344 if (i.ne.(nx-1)) then
345 idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
347 idx = ind(base.ge.range(i))
349 ;---------------------
352 data = ndtooned(area)
354 if (.not.any(ismissing(idx))) then
355 area_bin(i) = sum(data(idx))
357 area_bin(i) = area_bin@_FillValue
360 ;#############################################################
361 ; using model biome class:
362 ; set the following 4 classes to _FillValue:
363 ; (3)Needleleaf Deciduous Boreal Tree,
364 ; (8)Broadleaf Deciduous Boreal Tree,
365 ; (9)Broadleaf Evergreen Shrub,
368 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
369 area_bin(i) = area_bin@_FillValue
371 ;#############################################################
375 ;---------------------
376 ; for data_mod and data_ob
387 data = ndtooned(data_ob(m,k,:,:))
391 data = ndtooned(data_mod(m,k,:,:))
396 if (.not.any(ismissing(idx))) then
397 yvalues(t,n,i) = avg(data(idx))
399 yvalues(t,n,i) = yvalues@_FillValue
402 ;#############################################################
403 ; using model biome class:
404 ; set the following 4 classes to _FillValue:
405 ; (3)Needleleaf Deciduous Boreal Tree,
406 ; (8)Broadleaf Deciduous Boreal Tree,
407 ; (9)Broadleaf Evergreen Shrub,
410 if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
411 yvalues(t,n,i) = yvalues@_FillValue
413 ;#############################################################
428 ;----------------------------------------------------------------
431 good = ind(.not.ismissing(area_bin))
433 area_g = area_bin(good)
435 n_biome = dimsizes(good)
437 ;----------------------------------------------------------------
438 ; data for tseries plot
440 yvalues_g = new((/ntime,data_n,n_biome/),float)
442 yvalues_g@units = "TgC/month"
444 ; change unit to Tg C/month
445 ; change unit from g to Tg (Tera gram)
448 yvalues_g = yvalues(:,:,good) * conform(yvalues_g,area_g,2) * factor_unit
450 ;-------------------------------------------------------------------
451 ; general settings for line plot
454 res@xyDashPatterns = (/0,0/) ; make lines solid
455 res@xyLineThicknesses = (/2.0,2.0/) ; make lines thicker
456 res@xyLineColors = (/"blue","red"/) ; line color
458 res@trXMinF = year_start
459 res@trXMaxF = year_end + 1
461 res@vpHeightF = 0.4 ; change aspect ratio of plot
465 res@tiMainFontHeightF = 0.025 ; size of title
467 res@tmXBFormat = "f" ; not to add trailing zeros
469 ; res@gsnMaximize = True
471 ;----------------------------------------------
472 ; Add a boxed legend using the simple method
474 res@pmLegendDisplayMode = "Always"
475 ; res@pmLegendWidthF = 0.1
476 res@pmLegendWidthF = 0.08
477 res@pmLegendHeightF = 0.06
478 res@pmLegendOrthogonalPosF = -1.17
479 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
480 ; res@pmLegendOrthogonalPosF = -0.30 ;(downward)
482 ; res@pmLegendParallelPosF = 0.18
483 res@pmLegendParallelPosF = 0.23 ;(rightward)
484 res@pmLegendParallelPosF = 0.73 ;(rightward)
485 res@pmLegendParallelPosF = 0.83 ;(rightward)
487 ; res@lgPerimOn = False
488 res@lgLabelFontHeightF = 0.015
489 res@xyExplicitLegendLabels = (/"observed",model_name/)
491 ;*******************************************************************
492 ; (A) time series plot: monthly ( 2 lines per plot)
493 ;*******************************************************************
495 ; x-axis in time series plot
497 timeI = new((/ntime/),integer)
498 timeF = new((/ntime/),float)
499 timeI = ispan(1,ntime,1)
500 timeF = year_start + (timeI-1)/12.
501 timeF@long_name = "year"
503 plot_data = new((/2,ntime/),float)
504 plot_data@long_name = "TgC/month"
506 ;----------------------------------------------
507 ; time series : per biome
511 plot_name = "monthly_biome_"+ m
513 wks = gsn_open_wks (plot_type,plot_name)
515 title = "Fire : "+ row_head(m)
516 res@tiMainString = title
518 plot_data(0,:) = yvalues_g(:,0,m)
519 plot_data(1,:) = yvalues_g(:,1,m)
521 plot = gsn_csm_xy(wks,timeF,plot_data,res)
523 ; system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
524 ; "rm "+plot_name+"."+plot_type)
531 ;--------------------------------------------
532 ; time series: global
534 plot_name = "monthly_global"
536 wks = gsn_open_wks (plot_type,plot_name)
538 title = "Fire : "+ row_head(n_biome)
539 res@tiMainString = title
542 plot_data(0,k) = sum(yvalues_g(k,0,:))
543 plot_data(1,k) = sum(yvalues_g(k,1,:))
546 plot = gsn_csm_xy(wks,timeF,plot_data,res)
548 ; system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
549 ; "rm "+plot_name+"."+plot_type)