1 ;********************************************************
2 ; histogram normalized by rain and compute correleration
3 ;********************************************************
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
8 procedure pminmax(data:numeric,name:string)
10 print ("min/max " + name + " = " + min(data) + "/" + max(data))
11 if(isatt(data,"units")) then
12 print (name + " units = " + data@units)
24 ;************************************************
26 ;************************************************
27 ;film = "b30.061n_1995-2004_MONS_climo_lnd.nc"
28 ;model_name = "b30.061n"
31 ;film = "newcn05_ncep_1i_MONS_climo_lnd.nc"
35 ;film = "i01.06cn_1798-2004_MONS_climo.nc"
39 ;film = "i01.06casa_1798-2004_MONS_climo.nc"
40 ;model_name = "06casa"
43 film = "i01.10cn_1948-2004_MONS_climo.nc"
47 ;film = "i01.10casa_1948-2004_MONS_climo.nc"
48 ;model_name = "10casa"
51 html_name = "table.html." + model_name
52 html_new = html_name +".new"
53 system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \
54 "mv -f "+html_new+" "+html_name)
55 ;------------------------------------------------
56 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
57 fm = addfile(dirm+film,"r")
61 ;************************************************
62 ; read in data: observed
63 ;************************************************
65 ob_name = "MODIS MOD 15A2 2000-2005"
67 diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
68 filo1 = "land_class_"+model_grid+".nc"
69 filo2 = "LAI_2000-2005_MONS_"+model_grid+".nc"
71 fo1 = addfile(diro+filo1,"r")
72 fo2 = addfile(diro+filo2,"r")
74 classob = tofloat(fo1->LAND_CLASS)
77 ;*******************************************************************
79 ;*******************************************************************
81 table_header_name = "LAI"
83 ; column (not including header column)
85 col_header1 = (/"Mean","Max","Phase","Growth"/)
86 col_header2 = (/"ob","model","M" \
92 ncol1 = dimsizes(col_header1)
93 ncol2 = dimsizes(col_header2)
96 ; row (not including header row)
97 row_header = (/"Water Bodies" \
98 ,"Evergreen Needleleaf Forests" \
99 ,"Evergreen Broadleaf Forests" \
100 ,"Deciduous Needleleaf Forest" \
101 ,"Deciduous Broadleaf Forests" \
103 ,"Closed Bushlands" \
105 ,"Woody Savannas (S. Hem.)" \
106 ,"Savannas (S. Hem.)" \
108 ,"Permanent Wetlands" \
110 ,"Urban and Built-Up" \
111 ,"Cropland/Natural Vegetation Mosaic" \
112 ,"Permanent Snow and Ice" \
113 ,"Barren or Sparsely Vegetated" \
115 ,"Woody Savannas (N. Hem.)" \
116 ,"Savannas (N. Hem.)" \
117 ,"All biome average" \
119 nrow = dimsizes(row_header)
121 ; arrays to be passed to table.
122 value = new ((/nrow, ncol/),string )
127 ncr1 = (/1,1/) ; 1 row, 1 column
128 xx1 = (/0.005,0.25/) ; Start and end X
129 yy1 = (/0.900,0.995/) ; Start and end Y
130 text1 = table_header_name
132 res1@txFontHeightF = 0.03
133 res1@gsFillColor = "CornFlowerBlue"
135 ; Column header (equally space in x2)
136 ncr21 = (/1,ncol1/) ; 1 rows, 4 columns
137 xx21 = (/xx1(1),0.995/) ; start from end of x1
138 yy21 = (/0.9475,0.995/) ; half of y1
141 res21@txFontHeightF = 0.015
142 res21@gsFillColor = "Gray"
144 ncr22 = (/1,ncol2/) ; 1 rows, 12 columns
145 xx22 = xx21 ; start from end of x1
146 yy22 = (/0.900,0.9475/) ; half of y1
149 res22@txFontHeightF = 0.015
150 res22@gsFillColor = "Gray"
152 ; Row header (equally space in y2)
153 ncr3 = (/nrow,1/) ; 20 rows, 1 columns
154 xx3 = xx1 ; same as x1
155 yy3 = (/1.0-table_length,0.900/) ; end at start of y1
158 res3@txFontHeightF = 0.01
159 res3@gsFillColor = "Gray"
162 ncr4 = (/nrow,ncol/) ; 5 rows, 5 columns
163 xx4 = xx21 ; Start and end x
164 yy4 = yy3 ; Start and end Y
165 text4 = new((/nrow,ncol/),string)
167 color_fill4 = new((/nrow,ncol/),string)
168 color_fill4 = "white"
169 color_fill4(nrow-1,:) = "CornFlowerBlue"
171 res4 = True ; Set up resource list
172 ; res4@gsnDebug = True ; Useful to print NDC row,col values used.
173 res4@txFontHeightF = 0.015
174 res4@gsFillColor = color_fill4
178 ;************************************************
179 ; plot global land class: observed
180 ;************************************************
183 resg = True ; Use plot options
184 resg@cnFillOn = True ; Turn on color fill
185 resg@gsnSpreadColors = True ; use full colormap
186 ; resg@cnFillMode = "RasterFill" ; Turn on raster color
187 ; resg@lbLabelAutoStride = True
188 resg@cnLinesOn = False ; Turn off contourn lines
189 resg@mpFillOn = False ; Turn off map fill
191 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
192 resg@cnMinLevelValF = 1. ; Min level
193 resg@cnMaxLevelValF = 19. ; Max level
194 resg@cnLevelSpacingF = 1. ; interval
197 classob@_FillValue = 1.e+36
198 classob = where(classob.eq.0,classob@_FillValue,classob)
200 plot_name = "global_class_ob"
202 resg@tiMainString = title
204 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
205 gsn_define_colormap(wks,"gui_default") ; choose colormap
207 plot = gsn_csm_contour_map_ce(wks,classob,resg)
210 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
211 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
212 "rm "+plot_name+"-*."+plot_type_new+";"+ \
213 "rm "+plot_name+"."+plot_type)
216 ;*******************************************************************
217 ; Calculate "nice" bins for binning the data in equally spaced ranges
218 ;********************************************************************
220 range = fspan(0,nclassn-1,nclassn)
223 ; Use this range information to grab all the values in a
224 ; particular range, and then take an average.
228 xvalues = new((/2,nx/),float)
229 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
230 dx = xvalues(0,1) - xvalues(0,0) ; range width
231 dx4 = dx/4 ; 1/4 of the range
232 xvalues(1,:) = xvalues(0,:) - dx/5.
233 ;-----------------------------------------------------------------
235 ;--------------------------------------------------------------------
238 laiob_mean = dim_avg_Wrap(laiob(lat|:,lon|:,time|:))
239 laimod_mean = dim_avg_Wrap(laimod(lat|:,lon|:,time|:))
241 DATA11_1D = ndtooned(classob)
242 DATA12_1D = ndtooned(laiob_mean)
243 DATA22_1D = ndtooned(laimod_mean)
245 yvalues = new((/2,nx/),float)
246 mn_yvalues = new((/2,nx/),float)
247 mx_yvalues = new((/2,nx/),float)
251 ; See if we are doing model or observational data.
261 ; Loop through each range and check for values.
264 if (i.ne.(nr-2)) then
266 ; print("In range ["+range(i)+","+range(i+1)+")")
267 idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1)))
270 ; print("In range ["+range(i)+",)")
271 idx = ind(data_ob.ge.range(i))
274 ; Calculate average, and get min and max.
276 if(.not.any(ismissing(idx))) then
277 yvalues(nd,i) = avg(data_mod(idx))
278 mn_yvalues(nd,i) = min(data_mod(idx))
279 mx_yvalues(nd,i) = max(data_mod(idx))
280 count = dimsizes(idx)
283 yvalues(nd,i) = yvalues@_FillValue
284 mn_yvalues(nd,i) = yvalues@_FillValue
285 mx_yvalues(nd,i) = yvalues@_FillValue
288 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
289 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
291 ; Clean up for next time in loop.
298 ;-----------------------------------------------------------------
299 ; compute correlation coef and M score
304 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
308 ccrMean = esccr(uu,vv,0)
312 bias = sum(abs(vv-uu)/(vv+uu))
313 Mmean = (1.- (bias/dimsizes(uu)))*5.
315 M_lai_mean = sprintf("%.2f", Mmean)
316 system("sed s#M_lai_mean#"+M_lai_mean+"# "+html_name+" > "+html_new+";"+ \
317 "mv -f "+html_new+" "+html_name)
321 text4(i,0) = sprintf("%5.2f",u(i))
322 text4(i,1) = sprintf("%5.2f",v(i))
325 text4(nrow-1,0) = sprintf("%5.2f",avg(u))
326 text4(nrow-1,1) = sprintf("%5.2f",avg(v))
327 text4(nrow-1,2) = sprintf("%5.2f",Mmean)
333 ;--------------------------------------------------------------------
337 resm@gsnMaximize = True
339 resm@gsnFrame = False
340 resm@xyMarkLineMode = "Markers"
341 resm@xyMarkerSizeF = 0.014
343 resm@xyMarkerColors = (/"Brown","Blue"/)
344 ; resm@trYMinF = min(mn_yvalues) - 10.
345 ; resm@trYMaxF = max(mx_yvalues) + 10.
346 resm@trYMinF = min(mn_yvalues) - 2
347 resm@trYMaxF = max(mx_yvalues) + 4
349 resm@tiYAxisString = "Mean LAI (Leaf Area Index)"
350 resm@tiXAxisString = "Land Cover Type"
352 max_bar = new((/2,nx/),graphic)
353 min_bar = new((/2,nx/),graphic)
354 max_cap = new((/2,nx/),graphic)
355 min_cap = new((/2,nx/),graphic)
358 line_colors = (/"brown","blue"/)
359 ;------------------------------------------------------------------
360 ; Start the graphics.
362 plot_name = "histogram_mean"
363 title = model_name + " vs Observed"
364 resm@tiMainString = title
366 wks = gsn_open_wks (plot_type,plot_name)
367 ;-----------------------------
368 ; Add a boxed legend using the more simple method
370 resm@pmLegendDisplayMode = "Always"
371 ; resm@pmLegendWidthF = 0.1
372 resm@pmLegendWidthF = 0.08
373 resm@pmLegendHeightF = 0.05
374 resm@pmLegendOrthogonalPosF = -1.17
375 ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
376 ; resm@pmLegendParallelPosF = 0.18
377 resm@pmLegendParallelPosF = 0.88 ;(rightward)
379 ; resm@lgPerimOn = False
380 resm@lgLabelFontHeightF = 0.015
381 resm@xyExplicitLegendLabels = (/"observed",model_name/)
382 ;-----------------------------
384 tRes@txFontHeightF = 0.025
386 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")"
388 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
390 xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
391 ;-------------------------------
392 ;Attach the vertical bar and the horizontal cap line
395 lnres@gsLineColor = line_colors(nd)
398 if(.not.ismissing(mn_yvalues(nd,i)).and. \
399 .not.ismissing(mx_yvalues(nd,i))) then
401 ; Attach the vertical bar, both above and below the marker.
405 y2 = mn_yvalues(nd,i)
406 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
408 y2 = mx_yvalues(nd,i)
409 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
411 ; Attach the horizontal cap line, both above and below the marker.
413 x1 = xvalues(nd,i) - dx4
414 x2 = xvalues(nd,i) + dx4
415 y1 = mn_yvalues(nd,i)
416 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
418 y1 = mx_yvalues(nd,i)
419 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
427 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
428 "rm "+plot_name+"."+plot_type)
445 ;-----------------------------------------------------------------
448 resg = True ; Use plot options
449 resg@cnFillOn = True ; Turn on color fill
450 resg@gsnSpreadColors = True ; use full colormap
451 ; resg@cnFillMode = "RasterFill" ; Turn on raster color
452 ; resg@lbLabelAutoStride = True
453 resg@cnLinesOn = False ; Turn off contourn lines
454 resg@mpFillOn = False ; Turn off map fill
456 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
457 resg@cnMinLevelValF = 0. ; Min level
458 resg@cnMaxLevelValF = 10. ; Max level
459 resg@cnLevelSpacingF = 1. ; interval
464 laiob_mean = where(ismissing(laiob_mean).and.(ismissing(laimod_mean).or.(laimod_mean.lt.delta)),0.,laiob_mean)
466 plot_name = "global_mean_ob"
468 resg@tiMainString = title
470 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
471 gsn_define_colormap(wks,"gui_default") ; choose colormap
473 plot = gsn_csm_contour_map_ce(wks,laiob_mean,resg)
476 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
477 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
478 "rm "+plot_name+"-*."+plot_type_new+";"+ \
479 "rm "+plot_name+"."+plot_type)
482 ;------------------------------------------------------------------------
483 ;global contour model
485 plot_name = "global_mean_model"
486 title = "Model " + model_name
487 resg@tiMainString = title
489 wks = gsn_open_wks (plot_type,plot_name)
490 gsn_define_colormap(wks,"gui_default") ; choose colormap
492 plot = gsn_csm_contour_map_ce(wks,laimod_mean,resg)
495 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
496 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
497 "rm "+plot_name+"-*."+plot_type_new+";"+ \
498 "rm "+plot_name+"."+plot_type)
501 ;-----------------------------------------------------------------
503 ;--------------------------------------------------------------------
507 laiob_max = laiob(0,:,:)
509 laiob_max@long_name = "Leaf Area Index Max"
511 dsizes_z = dimsizes(laiob)
518 laiob_max(j,i) = max(s)
522 ; print (min(y)+"/"+max(y))
525 ;-------------------------
527 laimod_max = laimod(0,:,:)
529 laimod_max@long_name = "Leaf Area Index Max"
531 dsizes_z = dimsizes(laimod)
538 laimod_max(j,i) = max(s)
542 ; print (min(laimod_max)+"/"+max(laimod_max))
545 ;------------------------
546 DATA11_1D = ndtooned(classob)
547 DATA12_1D = ndtooned(laiob_max)
548 DATA22_1D = ndtooned(laimod_max)
550 yvalues = new((/2,nx/),float)
551 mn_yvalues = new((/2,nx/),float)
552 mx_yvalues = new((/2,nx/),float)
556 ; See if we are doing model or observational data.
566 ; Loop through each range and check for values.
569 if (i.ne.(nr-2)) then
571 ; print("In range ["+range(i)+","+range(i+1)+")")
572 idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))
575 ; print("In range ["+range(i)+",)")
576 idx = ind(range(i).le.data_ob)
579 ; Calculate average, and get min and max.
581 if(.not.any(ismissing(idx))) then
582 yvalues(nd,i) = avg(data_mod(idx))
583 mn_yvalues(nd,i) = min(data_mod(idx))
584 mx_yvalues(nd,i) = max(data_mod(idx))
585 count = dimsizes(idx)
588 yvalues(nd,i) = yvalues@_FillValue
589 mn_yvalues(nd,i) = yvalues@_FillValue
590 mx_yvalues(nd,i) = yvalues@_FillValue
593 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
594 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
596 ; Clean up for next time in loop.
603 ;-----------------------------------------------------------------
604 ; compute correlation coef and M score
609 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
613 ccrMax = esccr(uu,vv,0)
617 bias = sum(abs(vv-uu)/(vv+uu))
618 Mmax = (1.- (bias/dimsizes(uu)))*5.
620 M_lai_max = sprintf("%.2f", Mmax)
621 system("sed s#M_lai_max#"+M_lai_max+"# "+html_name+" > "+html_new+";"+ \
622 "mv -f "+html_new+" "+html_name)
626 text4(i,3) = sprintf("%5.2f",u(i))
627 text4(i,4) = sprintf("%5.2f",v(i))
630 text4(nrow-1,3) = sprintf("%5.2f",avg(u))
631 text4(nrow-1,4) = sprintf("%5.2f",avg(v))
632 text4(nrow-1,5) = sprintf("%5.2f",Mmax)
638 ;--------------------------------------------------------------------
642 resm@gsnMaximize = True
644 resm@gsnFrame = False
645 resm@xyMarkLineMode = "Markers"
646 resm@xyMarkerSizeF = 0.014
648 resm@xyMarkerColors = (/"Brown","Blue"/)
649 ; resm@trYMinF = min(mn_yvalues) - 10.
650 ; resm@trYMaxF = max(mx_yvalues) + 10.
651 resm@trYMinF = min(mn_yvalues) - 2
652 resm@trYMaxF = max(mx_yvalues) + 4
654 resm@tiYAxisString = "Max LAI (Leaf Area Index)"
655 resm@tiXAxisString = "Land Cover Type"
657 max_bar = new((/2,nx/),graphic)
658 min_bar = new((/2,nx/),graphic)
659 max_cap = new((/2,nx/),graphic)
660 min_cap = new((/2,nx/),graphic)
663 line_colors = (/"brown","blue"/)
664 ;------------------------------------------------------------------
665 ; Start the graphics.
667 plot_name = "histogram_max"
668 title = model_name + " vs Observed"
669 resm@tiMainString = title
671 wks = gsn_open_wks (plot_type,plot_name)
672 ;-----------------------------
673 ; Add a boxed legend using the more simple method
675 resm@pmLegendDisplayMode = "Always"
676 ; resm@pmLegendWidthF = 0.1
677 resm@pmLegendWidthF = 0.08
678 resm@pmLegendHeightF = 0.05
679 resm@pmLegendOrthogonalPosF = -1.17
680 ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
681 ; resm@pmLegendParallelPosF = 0.18
682 resm@pmLegendParallelPosF = 0.88 ;(rightward)
684 ; resm@lgPerimOn = False
685 resm@lgLabelFontHeightF = 0.015
686 resm@xyExplicitLegendLabels = (/"observed",model_name/)
687 ;-----------------------------
689 tRes@txFontHeightF = 0.025
691 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")"
693 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
695 xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
696 ;-------------------------------
697 ;Attach the vertical bar and the horizontal cap line
700 lnres@gsLineColor = line_colors(nd)
703 if(.not.ismissing(mn_yvalues(nd,i)).and. \
704 .not.ismissing(mx_yvalues(nd,i))) then
706 ; Attach the vertical bar, both above and below the marker.
710 y2 = mn_yvalues(nd,i)
711 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
713 y2 = mx_yvalues(nd,i)
714 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
716 ; Attach the horizontal cap line, both above and below the marker.
718 x1 = xvalues(nd,i) - dx4
719 x2 = xvalues(nd,i) + dx4
720 y1 = mn_yvalues(nd,i)
721 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
723 y1 = mx_yvalues(nd,i)
724 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
732 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
733 "rm "+plot_name+"."+plot_type)
750 ;-----------------------------------------------------------------
753 resg = True ; Use plot options
754 resg@cnFillOn = True ; Turn on color fill
755 resg@gsnSpreadColors = True ; use full colormap
756 ; resg@cnFillMode = "RasterFill" ; Turn on raster color
757 ; resg@lbLabelAutoStride = True
758 resg@cnLinesOn = False ; Turn off contourn lines
759 resg@mpFillOn = False ; Turn off map fill
761 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
762 resg@cnMinLevelValF = 0. ; Min level
763 resg@cnMaxLevelValF = 10. ; Max level
764 resg@cnLevelSpacingF = 1. ; interval
769 laiob_max = where(ismissing(laiob_max).and.(ismissing(laimod_max).or.(laimod_max.lt.delta)),0.,laiob_max)
771 plot_name = "global_max_ob"
773 resg@tiMainString = title
775 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
776 gsn_define_colormap(wks,"gui_default") ; choose colormap
778 plot_max = gsn_csm_contour_map_ce(wks,laiob_max,resg)
781 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
782 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
783 "rm "+plot_name+"-*."+plot_type_new+";"+ \
784 "rm "+plot_name+"."+plot_type)
787 ;------------------------------------------------------------------------
788 ;global contour model
790 plot_name = "global_max_model"
791 title = "Model " + model_name
792 resg@tiMainString = title
794 wks = gsn_open_wks (plot_type,plot_name)
795 gsn_define_colormap(wks,"gui_default") ; choose colormap
797 plot_max = gsn_csm_contour_map_ce(wks,laimod_max,resg)
800 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
801 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
802 "rm "+plot_name+"-*."+plot_type_new+";"+ \
803 "rm "+plot_name+"."+plot_type)
806 ;------------------------------------------------------------------------
808 ;--------------------------------------------------------------------
812 laiob_phase = laiob(0,:,:)
814 laiob_phase@long_name = "Leaf Area Index Max Month"
816 dsizes_z = dimsizes(laiob)
823 laiob_phase(j,i) = maxind(s) + 1
827 ; print (min(laiob_phase)+"/"+max(laiob_phase))
830 ;-------------------------
832 laimod_phase = laimod(0,:,:)
834 laimod_phase@long_name = "Leaf Area Index Max Month"
836 dsizes_z = dimsizes(laimod)
843 laimod_phase(j,i) = maxind(s) + 1
847 ; print (min(laimod_phase)+"/"+max(laimod_phase))
850 ;------------------------
851 DATA11_1D = ndtooned(classob)
852 DATA12_1D = ndtooned(laiob_phase)
853 DATA22_1D = ndtooned(laimod_phase)
855 yvalues = new((/2,nx/),float)
856 mn_yvalues = new((/2,nx/),float)
857 mx_yvalues = new((/2,nx/),float)
861 ; See if we are doing model or observational data.
871 ; Loop through each range and check for values.
874 if (i.ne.(nr-2)) then
876 ; print("In range ["+range(i)+","+range(i+1)+")")
877 idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))
880 ; print("In range ["+range(i)+",)")
881 idx = ind(range(i).le.data_ob)
884 ; Calculate average, and get min and max.
886 if(.not.any(ismissing(idx))) then
887 yvalues(nd,i) = avg(data_mod(idx))
888 mn_yvalues(nd,i) = min(data_mod(idx))
889 mx_yvalues(nd,i) = max(data_mod(idx))
890 count = dimsizes(idx)
893 yvalues(nd,i) = yvalues@_FillValue
894 mn_yvalues(nd,i) = yvalues@_FillValue
895 mx_yvalues(nd,i) = yvalues@_FillValue
898 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
899 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
901 ; Clean up for next time in loop.
908 ;-----------------------------------------------------------------
909 ; compute correlation coef and M score
914 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
918 ccrPhase = esccr(uu,vv,0)
922 ; bias = abs(avg(vv)-avg(uu))
924 bias = avg(abs(vv-uu))
926 bias = where((bias.gt. 6.),12.-bias,bias)
927 Mphase = ((6. - bias)/6.)*5.
929 M_lai_phase = sprintf("%.2f", Mphase)
930 system("sed s#M_lai_phase#"+M_lai_phase+"# "+html_name+" > "+html_new+";"+ \
931 "mv -f "+html_new+" "+html_name)
935 text4(i,6) = sprintf("%5.2f",u(i))
936 text4(i,7) = sprintf("%5.2f",v(i))
939 text4(nrow-1,6) = sprintf("%5.2f",avg(u))
940 text4(nrow-1,7) = sprintf("%5.2f",avg(v))
941 text4(nrow-1,8) = sprintf("%5.2f",Mphase)
947 ;--------------------------------------------------------------------
951 resm@gsnMaximize = True
953 resm@gsnFrame = False
954 resm@xyMarkLineMode = "Markers"
955 resm@xyMarkerSizeF = 0.014
957 resm@xyMarkerColors = (/"Brown","Blue"/)
958 ; resm@trYMinF = min(mn_yvalues) - 10.
959 ; resm@trYMaxF = max(mx_yvalues) + 10.
960 resm@trYMinF = min(mn_yvalues) - 2
961 resm@trYMaxF = max(mx_yvalues) + 4
963 resm@tiYAxisString = "Max LAI (Leaf Area Index) Month"
964 resm@tiXAxisString = "Land Cover Type"
966 max_bar = new((/2,nx/),graphic)
967 min_bar = new((/2,nx/),graphic)
968 max_cap = new((/2,nx/),graphic)
969 min_cap = new((/2,nx/),graphic)
972 line_colors = (/"brown","blue"/)
973 ;------------------------------------------------------------------
974 ; Start the graphics.
976 plot_name = "histogram_phase"
977 title = model_name + " vs Observed"
978 resm@tiMainString = title
980 wks = gsn_open_wks (plot_type,plot_name)
981 ;-----------------------------
982 ; Add a boxed legend using the more simple method
984 resm@pmLegendDisplayMode = "Always"
985 ; resm@pmLegendWidthF = 0.1
986 resm@pmLegendWidthF = 0.08
987 resm@pmLegendHeightF = 0.05
988 resm@pmLegendOrthogonalPosF = -1.17
989 ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
990 ; resm@pmLegendParallelPosF = 0.18
991 resm@pmLegendParallelPosF = 0.88 ;(rightward)
993 ; resm@lgPerimOn = False
994 resm@lgLabelFontHeightF = 0.015
995 resm@xyExplicitLegendLabels = (/"observed",model_name/)
996 ;-----------------------------
998 tRes@txFontHeightF = 0.025
1000 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")"
1002 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1004 xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
1005 ;-------------------------------
1006 ;Attach the vertical bar and the horizontal cap line
1009 lnres@gsLineColor = line_colors(nd)
1012 if(.not.ismissing(mn_yvalues(nd,i)).and. \
1013 .not.ismissing(mx_yvalues(nd,i))) then
1015 ; Attach the vertical bar, both above and below the marker.
1019 y2 = mn_yvalues(nd,i)
1020 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1022 y2 = mx_yvalues(nd,i)
1023 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1025 ; Attach the horizontal cap line, both above and below the marker.
1027 x1 = xvalues(nd,i) - dx4
1028 x2 = xvalues(nd,i) + dx4
1029 y1 = mn_yvalues(nd,i)
1030 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1032 y1 = mx_yvalues(nd,i)
1033 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1041 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1042 "rm "+plot_name+"."+plot_type)
1059 ;-----------------------------------------------------------------
1062 resg = True ; Use plot options
1063 resg@cnFillOn = True ; Turn on color fill
1064 resg@gsnSpreadColors = True ; use full colormap
1065 ; resg@cnFillMode = "RasterFill" ; Turn on raster color
1066 ; resg@lbLabelAutoStride = True
1067 resg@cnLinesOn = False ; Turn off contourn lines
1068 resg@mpFillOn = False ; Turn off map fill
1070 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1071 resg@cnMinLevelValF = 1. ; Min level
1072 resg@cnMaxLevelValF = 12. ; Max level
1073 resg@cnLevelSpacingF = 1. ; interval
1078 laiob_phase = where(ismissing(laiob_phase).and.(ismissing(laimod_phase).or.(laimod_phase.lt.delta)),0.,laiob_phase)
1080 plot_name = "global_phase_ob"
1082 resg@tiMainString = title
1084 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1085 gsn_define_colormap(wks,"gui_default") ; choose colormap
1087 plot = gsn_csm_contour_map_ce(wks,laiob_phase,resg)
1090 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1091 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1092 "rm "+plot_name+"-*."+plot_type_new+";"+ \
1093 "rm "+plot_name+"."+plot_type)
1096 ;------------------------------------------------------------------------
1097 ;global contour model
1099 plot_name = "global_phase_model"
1100 title = "Model " + model_name
1101 resg@tiMainString = title
1103 wks = gsn_open_wks (plot_type,plot_name)
1104 gsn_define_colormap(wks,"gui_default") ; choose colormap
1107 plot = gsn_csm_contour_map_ce(wks,laimod_phase,resg)
1110 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1111 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1112 "rm "+plot_name+"-*."+plot_type_new+";"+ \
1113 "rm "+plot_name+"."+plot_type)
1116 ;-----------------------------------------------------------------
1118 ;--------------------------------------------------------------------
1121 day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/)
1124 laiob_grow = laiob(0,:,:)
1125 laiob_grow@long_name = "Days of Growing Season"
1127 dsizes_z = dimsizes(laiob)
1136 if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then
1137 nday = nday + day_of_data(k)
1141 laiob_grow(j,i) = nday
1145 ; print (min(laiob_grow)+"/"+max(laiob_grow))
1146 ;-------------------------
1148 laimod_grow = laimod(0,:,:)
1149 laimod_grow@long_name = "Days of Growing Season"
1151 dsizes_z = dimsizes(laimod)
1160 if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then
1161 nday = nday + day_of_data(k)
1165 laimod_grow(j,i) = nday
1169 ; print (min(laimod_grow)+"/"+max(laimod_grow))
1170 ;------------------------
1171 DATA11_1D = ndtooned(classob)
1172 DATA12_1D = ndtooned(laiob_grow)
1173 DATA22_1D = ndtooned(laimod_grow)
1175 yvalues = new((/2,nx/),float)
1176 mn_yvalues = new((/2,nx/),float)
1177 mx_yvalues = new((/2,nx/),float)
1181 ; See if we are doing model or observational data.
1185 data_mod = DATA12_1D
1188 data_mod = DATA22_1D
1191 ; Loop through each range and check for values.
1194 if (i.ne.(nr-2)) then
1196 ; print("In range ["+range(i)+","+range(i+1)+")")
1197 idx = ind((range(i).le.data_ob).and.(data_ob.lt.range(i+1)))
1200 ; print("In range ["+range(i)+",)")
1201 idx = ind(range(i).le.data_ob)
1204 ; Calculate average, and get min and max.
1206 if(.not.any(ismissing(idx))) then
1207 yvalues(nd,i) = avg(data_mod(idx))
1208 mn_yvalues(nd,i) = min(data_mod(idx))
1209 mx_yvalues(nd,i) = max(data_mod(idx))
1210 count = dimsizes(idx)
1213 yvalues(nd,i) = yvalues@_FillValue
1214 mn_yvalues(nd,i) = yvalues@_FillValue
1215 mx_yvalues(nd,i) = yvalues@_FillValue
1218 ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1219 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1221 ; Clean up for next time in loop.
1228 ;-----------------------------------------------------------------
1229 ; compute correlation coef and M score
1234 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1238 ccrGrow = esccr(uu,vv,0)
1242 bias = sum(abs(vv-uu)/(vv+uu))
1243 Mgrow = (1.- (bias/dimsizes(uu)))*5.
1245 M_lai_grow = sprintf("%.2f", Mgrow)
1246 system("sed s#M_lai_grow#"+M_lai_grow+"# "+html_name+" > "+html_new+";"+ \
1247 "mv -f "+html_new+" "+html_name)
1251 text4(i,9) = sprintf("%5.2f",u(i))
1252 text4(i,10) = sprintf("%5.2f",v(i))
1255 text4(nrow-1,9) = sprintf("%5.2f",avg(u))
1256 text4(nrow-1,10) = sprintf("%5.2f",avg(v))
1257 text4(nrow-1,11) = sprintf("%5.2f",Mgrow)
1263 ;--------------------------------------------------------------------
1267 resm@gsnMaximize = True
1268 resm@gsnDraw = False
1269 resm@gsnFrame = False
1270 resm@xyMarkLineMode = "Markers"
1271 resm@xyMarkerSizeF = 0.014
1273 resm@xyMarkerColors = (/"Brown","Blue"/)
1274 ; resm@trYMinF = min(mn_yvalues) - 10.
1275 ; resm@trYMaxF = max(mx_yvalues) + 10.
1276 resm@trYMinF = min(mn_yvalues) - 2
1277 resm@trYMaxF = max(mx_yvalues) + 4
1279 resm@tiYAxisString = "Days of Growing season"
1280 resm@tiXAxisString = "Land Cover Type"
1282 max_bar = new((/2,nx/),graphic)
1283 min_bar = new((/2,nx/),graphic)
1284 max_cap = new((/2,nx/),graphic)
1285 min_cap = new((/2,nx/),graphic)
1288 line_colors = (/"brown","blue"/)
1289 ;------------------------------------------------------------------
1290 ; Start the graphics.
1292 plot_name = "histogram_grow"
1293 title = model_name + " vs Observed"
1294 resm@tiMainString = title
1296 wks = gsn_open_wks (plot_type,plot_name)
1297 ;-----------------------------
1298 ; Add a boxed legend using the more simple method
1300 resm@pmLegendDisplayMode = "Always"
1301 ; resm@pmLegendWidthF = 0.1
1302 resm@pmLegendWidthF = 0.08
1303 resm@pmLegendHeightF = 0.05
1304 resm@pmLegendOrthogonalPosF = -1.17
1305 ; resm@pmLegendOrthogonalPosF = -1.00 ;(downward)
1306 ; resm@pmLegendParallelPosF = 0.18
1307 resm@pmLegendParallelPosF = 0.88 ;(rightward)
1309 ; resm@lgPerimOn = False
1310 resm@lgLabelFontHeightF = 0.015
1311 resm@xyExplicitLegendLabels = (/"observed",model_name/)
1312 ;-----------------------------
1314 tRes@txFontHeightF = 0.025
1316 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")"
1318 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1320 xy = gsn_csm_xy(wks,xvalues,yvalues,resm)
1321 ;-------------------------------
1322 ;Attach the vertical bar and the horizontal cap line
1325 lnres@gsLineColor = line_colors(nd)
1328 if(.not.ismissing(mn_yvalues(nd,i)).and. \
1329 .not.ismissing(mx_yvalues(nd,i))) then
1331 ; Attach the vertical bar, both above and below the marker.
1335 y2 = mn_yvalues(nd,i)
1336 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1338 y2 = mx_yvalues(nd,i)
1339 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1341 ; Attach the horizontal cap line, both above and below the marker.
1343 x1 = xvalues(nd,i) - dx4
1344 x2 = xvalues(nd,i) + dx4
1345 y1 = mn_yvalues(nd,i)
1346 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1348 y1 = mx_yvalues(nd,i)
1349 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1357 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1358 "rm "+plot_name+"."+plot_type)
1375 ;-----------------------------------------------------------------
1378 resg = True ; Use plot options
1379 resg@cnFillOn = True ; Turn on color fill
1380 resg@gsnSpreadColors = True ; use full colormap
1381 ; resg@cnFillMode = "RasterFill" ; Turn on raster color
1382 ; resg@lbLabelAutoStride = True
1383 resg@cnLinesOn = False ; Turn off contourn lines
1384 resg@mpFillOn = False ; Turn off map fill
1386 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1387 resg@cnMinLevelValF = 60. ; Min level
1388 resg@cnMaxLevelValF = 360. ; Max level
1389 resg@cnLevelSpacingF = 20. ; interval
1393 laiob_grow@_FillValue = 1.e+36
1394 laiob_grow = where(laiob_grow .lt. 10.,laiob_grow@_FillValue,laiob_grow)
1396 plot_name = "global_grow_ob"
1398 resg@tiMainString = title
1400 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1401 gsn_define_colormap(wks,"gui_default") ; choose colormap
1403 plot = gsn_csm_contour_map_ce(wks,laiob_grow,resg)
1406 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1407 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1408 "rm "+plot_name+"-*."+plot_type_new+";"+ \
1409 "rm "+plot_name+"."+plot_type)
1412 ;------------------------------------------------------------------------
1413 ;global contour model
1415 laimod_grow@_FillValue = 1.e+36
1416 laimod_grow = where(laimod_grow .lt. 10.,laimod_grow@_FillValue,laimod_grow)
1418 plot_name = "global_grow_model"
1419 title = "Model " + model_name
1420 resg@tiMainString = title
1422 wks = gsn_open_wks (plot_type,plot_name)
1423 gsn_define_colormap(wks,"gui_default") ; choose colormap
1426 plot = gsn_csm_contour_map_ce(wks,laimod_grow,resg)
1429 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1430 "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
1431 "rm "+plot_name+"-*."+plot_type_new+";"+ \
1432 "rm "+plot_name+"."+plot_type)
1435 ;------------------------------------------------------------------------
1436 ;global contour model vs ob
1438 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1439 resg@cnMinLevelValF = 0. ; Min level
1440 resg@cnMaxLevelValF = 10. ; Max level
1441 resg@cnLevelSpacingF = 1. ; interval
1443 plot_name = "global_mean_model_vs_ob"
1445 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1446 gsn_define_colormap(wks,"gui_default") ; choose colormap
1449 plot=new(3,graphic) ; create graphic array
1451 resg@gsnFrame = False ; Do not draw plot
1452 resg@gsnDraw = False ; Do not advance frame
1454 ; plot correlation coef
1457 gRes@txFontHeightF = 0.02
1460 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMean)+")"
1462 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1463 ;--------------------------------------------------------------------
1468 resg@tiMainString = title
1470 plot(0) = gsn_csm_contour_map_ce(wks,laiob_mean,resg)
1474 title = "Model "+ model_name
1475 resg@tiMainString = title
1477 plot(1) = gsn_csm_contour_map_ce(wks,laimod_mean,resg)
1482 zz = laimod_mean - laiob_mean
1483 title = "Model_"+model_name+" - Observed"
1484 resg@tiMainString = title
1486 resg@cnMinLevelValF = -2. ; Min level
1487 resg@cnMaxLevelValF = 2. ; Max level
1488 resg@cnLevelSpacingF = 0.4 ; interval
1491 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1493 pres = True ; panel plot mods desired
1494 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1495 ; indiv. plots in panel
1496 pres@gsnMaximize = True ; fill the page
1498 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1500 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1501 "rm "+plot_name+"."+plot_type)
1505 ;-----------------------------------------------------------------
1506 ;global contour model vs ob
1508 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1509 resg@cnMinLevelValF = 0. ; Min level
1510 resg@cnMaxLevelValF = 10. ; Max level
1511 resg@cnLevelSpacingF = 1. ; interval
1513 plot_name = "global_max_model_vs_ob"
1515 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1516 gsn_define_colormap(wks,"gui_default") ; choose colormap
1519 plot=new(3,graphic) ; create graphic array
1521 resg@gsnFrame = False ; Do not draw plot
1522 resg@gsnDraw = False ; Do not advance frame
1524 ; plot correlation coef
1527 gRes@txFontHeightF = 0.02
1530 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrMax)+")"
1532 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1533 ;--------------------------------------------------------------------
1537 resg@tiMainString = title
1539 plot(0) = gsn_csm_contour_map_ce(wks,laiob_max,resg)
1543 title = "Model "+ model_name
1544 resg@tiMainString = title
1546 plot(1) = gsn_csm_contour_map_ce(wks,laimod_max,resg)
1552 zz = laimod_max - laiob_max
1553 title = "Model_"+model_name+" - Observed"
1554 resg@tiMainString = title
1556 resg@cnMinLevelValF = -6. ; Min level
1557 resg@cnMaxLevelValF = 6. ; Max level
1558 resg@cnLevelSpacingF = 1. ; interval
1561 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1563 pres = True ; panel plot mods desired
1564 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1565 ; indiv. plots in panel
1566 pres@gsnMaximize = True ; fill the page
1568 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1570 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1571 "rm "+plot_name+"."+plot_type)
1575 ;-----------------------------------------------------------------
1576 ;global contour model vs ob
1578 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1579 resg@cnMinLevelValF = 1. ; Min level
1580 resg@cnMaxLevelValF = 12. ; Max level
1581 resg@cnLevelSpacingF = 1. ; interval
1583 plot_name = "global_phase_model_vs_ob"
1585 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1586 gsn_define_colormap(wks,"gui_default") ; choose colormap
1589 plot=new(3,graphic) ; create graphic array
1591 resg@gsnFrame = False ; Do not draw plot
1592 resg@gsnDraw = False ; Do not advance frame
1594 ; plot correlation coef
1597 gRes@txFontHeightF = 0.02
1600 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrPhase)+")"
1602 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1603 ;--------------------------------------------------------------------
1607 resg@tiMainString = title
1609 plot(0) = gsn_csm_contour_map_ce(wks,laiob_phase,resg)
1613 title = "Model "+ model_name
1614 resg@tiMainString = title
1616 plot(1) = gsn_csm_contour_map_ce(wks,laimod_phase,resg)
1622 zz = laimod_phase - laiob_phase
1623 title = "Model_"+model_name+" - Observed"
1624 resg@tiMainString = title
1626 resg@cnMinLevelValF = -6. ; Min level
1627 resg@cnMaxLevelValF = 6. ; Max level
1628 resg@cnLevelSpacingF = 1. ; interval
1631 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1633 pres = True ; panel plot mods desired
1634 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1635 ; indiv. plots in panel
1636 pres@gsnMaximize = True ; fill the page
1638 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1640 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1641 "rm "+plot_name+"."+plot_type)
1645 ;------------------------------------------------------------------
1646 ;global contour model vs ob
1648 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1649 resg@cnMinLevelValF = 60. ; Min level
1650 resg@cnMaxLevelValF = 360. ; Max level
1651 resg@cnLevelSpacingF = 20. ; interval
1653 plot_name = "global_grow_model_vs_ob"
1655 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1656 gsn_define_colormap(wks,"gui_default") ; choose colormap
1659 plot=new(3,graphic) ; create graphic array
1661 resg@gsnFrame = False ; Do not draw plot
1662 resg@gsnDraw = False ; Do not advance frame
1664 ; plot correlation coef
1667 gRes@txFontHeightF = 0.02
1670 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrGrow)+")"
1672 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1673 ;--------------------------------------------------------------------
1677 resg@tiMainString = title
1679 plot(0) = gsn_csm_contour_map_ce(wks,laiob_grow,resg)
1683 title = "Model "+ model_name
1684 resg@tiMainString = title
1686 plot(1) = gsn_csm_contour_map_ce(wks,laimod_grow,resg)
1692 zz = laimod_grow - laiob_grow
1693 title = "Model_"+model_name+" - Observed"
1694 resg@tiMainString = title
1696 resg@cnMinLevelValF = -100. ; Min level
1697 resg@cnMaxLevelValF = 100. ; Max level
1698 resg@cnLevelSpacingF = 20. ; interval
1700 plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1702 pres = True ; panel plot mods desired
1703 ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1704 ; indiv. plots in panel
1705 pres@gsnMaximize = True ; fill the page
1707 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1709 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1710 "rm "+plot_name+"."+plot_type)
1714 ;**************************************************
1716 ;**************************************************
1718 plot_name = "table_lai"
1719 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1721 gsn_table(wks,ncr1,xx1,yy1,text1,res1)
1722 gsn_table(wks,ncr21,xx21,yy21,text21,res21)
1723 gsn_table(wks,ncr22,xx22,yy22,text22,res22)
1724 gsn_table(wks,ncr3,xx3,yy3,text3,res3)
1725 gsn_table(wks,ncr4,xx4,yy4,text4,res4)
1729 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1730 "rm "+plot_name+"."+plot_type)
1731 ;-------------------------------------------------------------------
1732 temp_name = "lai." + model_name
1733 system("mkdir -p " + temp_name+";"+ \
1734 "cp "+ html_name + " " +temp_name+"/table.html"+";"+ \
1735 "mv *.png " + temp_name +";"+ \
1736 "tar cf "+ temp_name +".tar " + temp_name)
1737 ;-------------------------------------------------------------------