1 ; ***********************************************
2 ; combine all scatter and all histogram
3 ; ***********************************************
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
7 ;************************************************
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)
22 ;************************************************
24 ;************************************************
26 dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
28 f81 = addfile (dir81+fil81,"r")
32 rain81 = tofloat(f81->PREC_ANN)
36 id81@long_name = "SITE_ID"
38 ;change longitude from (-180,180) to (0,360)
39 ;for model data interpolation
41 do i= 0,dimsizes(x81)-1
42 if (x81(i) .lt. 0.) then
47 ;-------------------------------------
49 dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
50 fil933 = "data.933.nc"
51 f933 = addfile (dir933+fil933,"r")
59 id933@long_name = "SITE_ID"
61 ;change longitude from (-180,180) to (0,360)
62 ;for model data interpolation
64 do i= 0,dimsizes(x933)-1
65 if (x933(i) .lt. 0.) then
66 x933(i) = x933(i)+ 360.
70 ;----------------------------------------
72 dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
73 filglobe = "Npp_T31_mean.nc"
74 fglobe = addfile (dirglobe+filglobe,"r")
76 nppglobe = fglobe->NPP
78 ;************************************************
80 ;************************************************
81 model_name = "b30.061n"
83 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
84 film = "b30.061n_1995-2004_ANN_climo_lnd.nc"
85 fm = addfile (dirm+film,"r")
92 nppmod = nppmod0(0,:,:)
93 rainmod = rainmod0(0,:,:)
97 nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0)
99 nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)
101 rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)
103 rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)
105 ;************************************************
106 ; Units for these variables are:
110 ; npp81 : g C/m^2/year
111 ; nppmod81: g C/m^2/s
112 ; nppglobe: g C/m^2/year
114 ; We want to convert these to "m/year" and "g C/m^2/year".
116 nsec_per_year = 60*60*24*365
118 rain81 = rain81 / 1000.
119 rainmod81 = (rainmod81/ 1000.) * nsec_per_year
120 nppmod81 = nppmod81 * nsec_per_year
122 rain933 = rain933 / 1000.
123 rainmod933 = (rainmod933/ 1000.) * nsec_per_year
124 nppmod933 = nppmod933 * nsec_per_year
126 nppmod = nppmod * nsec_per_year
128 npp81@units = "gC/m^2/yr"
129 nppmod81@units = "gC/m^2/yr"
130 npp933@units = "gC/m^2/yr"
131 nppmod933@units = "gC/m^2/yr"
132 nppmod@units = "gC/m^2/yr"
133 nppglobe@units = "gC/m^2/yr"
134 rain81@units = "m/yr"
135 rainmod81@units = "m/yr"
136 rain933@units = "m/yr"
137 rainmod933@units = "m/yr"
139 npp81@long_name = "NPP (gC/m2/year)"
140 npp933@long_name = "NPP (gC/m2/year)"
141 nppmod81@long_name = "NPP (gC/m2/year)"
142 nppmod933@long_name = "NPP (gC/m2/year)"
143 nppmod@long_name = "NPP (gC/m2/year)"
144 nppglobe@long_name = "NPP (gC/m2/year)"
145 rain81@long_name = "PREC (m/year)"
146 rain933@long_name = "PREC (m/year)"
147 rainmod81@long_name = "PREC (m/year)"
148 rainmod933@long_name = "PREC (m/year)"
150 ;(A) plot scatter ob 81
152 plot_name = "scatter_ob_81"
153 title = "Observed 81 sites"
155 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
156 res = True ; plot mods desired
157 res@tiMainString = title ; add title
158 res@xyMarkLineModes = "Markers" ; choose which have markers
159 res@xyMarkers = 16 ; choose type of marker
160 res@xyMarkerColor = "red" ; Marker color
161 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
162 res@tmLabelAutoStride = True ; nice tick mark labels
164 plot = gsn_csm_xy (wks,id81,npp81,res) ; create plot
167 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
168 system("rm "+plot_name+"."+plot_type)
169 system("rm "+plot_name+"-1."+plot_type_new)
170 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
174 ;(B) plot scatter ob 933
176 plot_name = "scatter_ob_933"
177 title = "Observed 933 sites"
179 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
180 res = True ; plot mods desired
181 res@tiMainString = title ; add title
182 res@xyMarkLineModes = "Markers" ; choose which have markers
183 res@xyMarkers = 16 ; choose type of marker
184 res@xyMarkerColor = "red" ; Marker color
185 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
186 res@tmLabelAutoStride = True ; nice tick mark labels
188 plot = gsn_csm_xy (wks,id933,npp933,res) ; create plot
191 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
192 system("rm "+plot_name+"."+plot_type)
193 system("rm "+plot_name+"-1."+plot_type_new)
194 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
198 ;(C) plot scatter model 81
200 plot_name = "scatter_model_81"
201 title = "Model "+ model_name +" 81 sites"
203 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
204 res = True ; plot mods desired
205 res@tiMainString = title ; add title
206 res@xyMarkLineModes = "Markers" ; choose which have markers
207 res@xyMarkers = 16 ; choose type of marker
208 res@xyMarkerColor = "red" ; Marker color
209 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
210 res@tmLabelAutoStride = True ; nice tick mark labels
212 plot = gsn_csm_xy (wks,id81,nppmod81,res) ; create plot
215 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
216 system("rm "+plot_name+"."+plot_type)
217 system("rm "+plot_name+"-1."+plot_type_new)
218 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
222 ;(D) plot scatter model 933
224 plot_name = "scatter_model_933"
225 title = "Model "+ model_name +" 933 sites"
227 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
228 res = True ; plot mods desired
229 res@tiMainString = title ; add title
230 res@xyMarkLineModes = "Markers" ; choose which have markers
231 res@xyMarkers = 16 ; choose type of marker
232 res@xyMarkerColor = "red" ; Marker color
233 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
234 res@tmLabelAutoStride = True ; nice tick mark labels
236 plot = gsn_csm_xy (wks,id933,nppmod933,res) ; create plot
239 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
240 system("rm "+plot_name+"."+plot_type)
241 system("rm "+plot_name+"-1."+plot_type_new)
242 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
246 ;(E) scatter model vs ob 81
248 plot_name = "scatter_model_vs_ob_81"
249 title = model_name +" vs ob 81 sites"
251 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
252 res = True ; plot mods desired
253 res@tiMainString = title ; add title
254 res@xyMarkLineModes = "Markers" ; choose which have markers
255 res@xyMarkers = 16 ; choose type of marker
256 res@xyMarkerColor = "red" ; Marker color
257 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
258 res@tmLabelAutoStride = True ; nice tick mark labels
259 ;-------------------------------
260 ;compute correlation coef. and M
261 ccr81 = esccr(nppmod81,npp81,0)
265 bias = sum(abs(nppmod81-npp81)/(nppmod81+npp81))
266 M81 = (1. - (bias/dimsizes(y81)))*5.
270 tRes@txFontHeightF = 0.025
272 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
274 gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
275 ;--------------------------------
276 plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot
279 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
280 system("rm "+plot_name+"."+plot_type)
281 system("rm "+plot_name+"-1."+plot_type_new)
282 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
286 ;(F) scatter model vs ob 933
288 plot_name = "scatter_model_vs_ob_933"
289 title = model_name +" vs ob 933 sites"
291 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
292 res = True ; plot mods desired
293 res@tiMainString = title ; add title
294 res@xyMarkLineModes = "Markers" ; choose which have markers
295 res@xyMarkers = 16 ; choose type of marker
296 res@xyMarkerColor = "red" ; Marker color
297 res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
298 res@tmLabelAutoStride = True ; nice tick mark labels
299 ;-------------------------------
300 ;compute correlation coef. and M
301 ccr933 = esccr(nppmod933,npp933,0)
305 bias = sum(abs(nppmod933-npp933)/(nppmod933+npp933))
306 M933 = (1. - (bias/dimsizes(y933)))*5.
310 tRes@txFontHeightF = 0.025
312 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
314 gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
315 ;--------------------------------
316 plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot
319 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
320 system("rm "+plot_name+"."+plot_type)
321 system("rm "+plot_name+"-1."+plot_type_new)
322 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
325 ;**************************************************************************
328 ;--------------------------------------------------------------------
331 RAIN1_1D = ndtooned(rain81)
332 RAIN2_1D = ndtooned(rainmod81)
333 NPP1_1D = ndtooned(npp81)
334 NPP2_1D = ndtooned(nppmod81)
336 ; Calculate some "nice" bins for binning the data in equally spaced
339 nbins = 15 ; Number of bins to use.
341 nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
342 nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
343 range = fspan(nicevals(0),nicevals(1),nvals)
345 ; Use this range information to grab all the values in a
346 ; particular range, and then take an average.
350 xvalues = new((/2,nx/),typeof(RAIN1_1D))
351 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
352 dx = xvalues(0,1) - xvalues(0,0) ; range width
353 dx4 = dx/4 ; 1/4 of the range
354 xvalues(1,:) = xvalues(0,:) - dx/5.
355 yvalues = new((/2,nx/),typeof(RAIN1_1D))
356 mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
357 mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
361 ; See if we are doing model or observational data.
371 ; Loop through each range and check for values.
374 if (i.ne.(nr-2)) then
376 ; print("In range ["+range(i)+","+range(i+1)+")")
377 idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
380 ; print("In range ["+range(i)+",)")
381 idx = ind(range(i).le.data)
384 ; Calculate average, and get min and max.
386 if(.not.any(ismissing(idx))) then
387 yvalues(nd,i) = avg(npp_data(idx))
388 mn_yvalues(nd,i) = min(npp_data(idx))
389 mx_yvalues(nd,i) = max(npp_data(idx))
390 count = dimsizes(idx)
393 yvalues(nd,i) = yvalues@_FillValue
394 mn_yvalues(nd,i) = yvalues@_FillValue
395 mx_yvalues(nd,i) = yvalues@_FillValue
398 ; Print out information.
400 ; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
401 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
404 ; Clean up for next time in loop.
411 ;----------------------------------------
412 ;compute correlation coeff and M score
416 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
420 ccr81h = esccr(uu,vv,0)
424 bias = sum(abs(vv-uu)/(vv+uu))
425 M81h = (1.- (bias/dimsizes(uu)))*5.
432 ;----------------------------------------------------------------------
436 res@gsnMaximize = True
439 res@xyMarkLineMode = "Markers"
440 res@xyMarkerSizeF = 0.014
442 res@xyMarkerColors = (/"Brown","Blue"/)
443 res@trYMinF = min(mn_yvalues) - 10.
444 res@trYMaxF = max(mx_yvalues) + 10.
446 res@tiYAxisString = "NPP (g C/m2/year)"
447 res@tiXAxisString = "Precipitation (m/year)"
449 max_bar = new((/2,nx/),graphic)
450 min_bar = new((/2,nx/),graphic)
451 max_cap = new((/2,nx/),graphic)
452 min_cap = new((/2,nx/),graphic)
455 line_colors = (/"brown","blue"/)
456 ;=================================================================
457 ; histogram ob 81 site only
459 plot_name = "histogram_ob_81"
460 title = "Observed 81 site"
461 res@tiMainString = title
463 wks = gsn_open_wks (plot_type,plot_name)
465 xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),res)
467 ;-------------------------------
468 ;Attach the vertical bar and the horizontal cap line
471 lnres@gsLineColor = line_colors(nd)
474 if(.not.ismissing(mn_yvalues(nd,i)).and. \
475 .not.ismissing(mx_yvalues(nd,i))) then
477 ; Attach the vertical bar, both above and below the marker.
481 y2 = mn_yvalues(nd,i)
482 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
484 y2 = mx_yvalues(nd,i)
485 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
487 ; Attach the horizontal cap line, both above and below the marker.
489 x1 = xvalues(nd,i) - dx4
490 x2 = xvalues(nd,i) + dx4
491 y1 = mn_yvalues(nd,i)
492 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
494 y1 = mx_yvalues(nd,i)
495 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
503 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
504 ; system("rm "+plot_name+"."+plot_type)
505 ; system("rm "+plot_name+"-1."+plot_type_new)
506 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
509 ;===========================================================================
510 ; histogram model vs ob 81 site
512 plot_name = "histogram_mod-ob_81"
513 title = model_name+ " vs Observed 81 site"
514 res@tiMainString = title
516 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
518 ;-----------------------------
519 ; Add a boxed legend using the more simple method, which won't have
520 ; vertical lines going through the markers.
522 res@pmLegendDisplayMode = "Always"
523 ; res@pmLegendWidthF = 0.1
524 res@pmLegendWidthF = 0.08
525 res@pmLegendHeightF = 0.05
526 res@pmLegendOrthogonalPosF = -1.17
527 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
528 ; res@pmLegendParallelPosF = 0.18
529 res@pmLegendParallelPosF = 0.23 ;(rightward)
531 ; res@lgPerimOn = False
532 res@lgLabelFontHeightF = 0.015
533 res@xyExplicitLegendLabels = (/"observed",model_name/)
534 ;-----------------------------
536 tRes@txFontHeightF = 0.025
538 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
540 gsn_text_ndc(wks,correlation_text,0.5,0.8,tRes)
542 xy = gsn_csm_xy(wks,xvalues,yvalues,res)
543 ;-------------------------------
544 ;Attach the vertical bar and the horizontal cap line
547 lnres@gsLineColor = line_colors(nd)
550 if(.not.ismissing(mn_yvalues(nd,i)).and. \
551 .not.ismissing(mx_yvalues(nd,i))) then
553 ; Attach the vertical bar, both above and below the marker.
557 y2 = mn_yvalues(nd,i)
558 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
560 y2 = mx_yvalues(nd,i)
561 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
563 ; Attach the horizontal cap line, both above and below the marker.
565 x1 = xvalues(nd,i) - dx4
566 x2 = xvalues(nd,i) + dx4
567 y1 = mn_yvalues(nd,i)
568 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
570 y1 = mx_yvalues(nd,i)
571 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
578 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
579 ;system("rm "+plot_name+"."+plot_type)
580 ;system("rm "+plot_name+"-1."+plot_type_new)
581 ;system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
599 ;**************************************************************************
602 ;--------------------------------------------------------------------
605 RAIN1_1D = ndtooned(rain933)
606 RAIN2_1D = ndtooned(rainmod933)
607 NPP1_1D = ndtooned(npp933)
608 NPP2_1D = ndtooned(nppmod933)
610 ; Calculate some "nice" bins for binning the data in equally spaced
613 nbins = 15 ; Number of bins to use.
615 nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
616 nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
617 range = fspan(nicevals(0),nicevals(1),nvals)
619 ; Use this range information to grab all the values in a
620 ; particular range, and then take an average.
624 xvalues = new((/2,nx/),typeof(RAIN1_1D))
625 xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
626 dx = xvalues(0,1) - xvalues(0,0) ; range width
627 dx4 = dx/4 ; 1/4 of the range
628 xvalues(1,:) = xvalues(0,:) - dx/5.
629 yvalues = new((/2,nx/),typeof(RAIN1_1D))
630 mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
631 mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
635 ; See if we are doing model or observational data.
645 ; Loop through each range and check for values.
648 if (i.ne.(nr-2)) then
650 ; print("In range ["+range(i)+","+range(i+1)+")")
651 idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
654 ; print("In range ["+range(i)+",)")
655 idx = ind(range(i).le.data)
658 ; Calculate average, and get min and max.
660 if(.not.any(ismissing(idx))) then
661 yvalues(nd,i) = avg(npp_data(idx))
662 mn_yvalues(nd,i) = min(npp_data(idx))
663 mx_yvalues(nd,i) = max(npp_data(idx))
664 count = dimsizes(idx)
667 yvalues(nd,i) = yvalues@_FillValue
668 mn_yvalues(nd,i) = yvalues@_FillValue
669 mx_yvalues(nd,i) = yvalues@_FillValue
672 ; Print out information.
674 ; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
675 ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
678 ; Clean up for next time in loop.
685 ;----------------------------------------
686 ;compute correlation coeff and M score
690 good = ind(.not.ismissing(u) .and. .not.ismissing(v))
694 ccr933h = esccr(uu,vv,0)
698 bias = sum(abs(vv-uu)/(vv+uu))
699 M933h = (1.- (bias/dimsizes(uu)))*5.
706 ;----------------------------------------------------------------------
710 res@gsnMaximize = True
713 res@xyMarkLineMode = "Markers"
714 res@xyMarkerSizeF = 0.014
716 res@xyMarkerColors = (/"Brown","Blue"/)
717 res@trYMinF = min(mn_yvalues) - 10.
718 res@trYMaxF = max(mx_yvalues) + 10.
720 res@tiYAxisString = "NPP (g C/m2/year)"
721 res@tiXAxisString = "Precipitation (m/year)"
723 max_bar = new((/2,nx/),graphic)
724 min_bar = new((/2,nx/),graphic)
725 max_cap = new((/2,nx/),graphic)
726 min_cap = new((/2,nx/),graphic)
729 line_colors = (/"brown","blue"/)
730 ;=================================================================
731 ; histogram ob 933 site only
733 plot_name = "histogram_ob_933"
734 title = "Observed 933 site"
735 res@tiMainString = title
737 wks = gsn_open_wks (plot_type,plot_name)
739 xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),res)
741 ;-------------------------------
742 ;Attach the vertical bar and the horizontal cap line
745 lnres@gsLineColor = line_colors(nd)
748 if(.not.ismissing(mn_yvalues(nd,i)).and. \
749 .not.ismissing(mx_yvalues(nd,i))) then
751 ; Attach the vertical bar, both above and below the marker.
755 y2 = mn_yvalues(nd,i)
756 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
758 y2 = mx_yvalues(nd,i)
759 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
761 ; Attach the horizontal cap line, both above and below the marker.
763 x1 = xvalues(nd,i) - dx4
764 x2 = xvalues(nd,i) + dx4
765 y1 = mn_yvalues(nd,i)
766 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
768 y1 = mx_yvalues(nd,i)
769 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
777 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
778 ; system("rm "+plot_name+"."+plot_type)
779 ; system("rm "+plot_name+"-1."+plot_type_new)
780 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
783 ;===========================================================================
784 ; histogram model vs ob 933 site
786 plot_name = "histogram_mod-ob_933"
787 title = model_name+ " vs Observed 933 site"
788 res@tiMainString = title
790 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
792 ;-----------------------------
793 ; Add a boxed legend using the more simple method, which won't have
794 ; vertical lines going through the markers.
796 res@pmLegendDisplayMode = "Always"
797 ; res@pmLegendWidthF = 0.1
798 res@pmLegendWidthF = 0.08
799 res@pmLegendHeightF = 0.05
800 res@pmLegendOrthogonalPosF = -1.17
801 ; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
802 ; res@pmLegendParallelPosF = 0.18
803 res@pmLegendParallelPosF = 0.23 ;(rightward)
805 ; res@lgPerimOn = False
806 res@lgLabelFontHeightF = 0.015
807 res@xyExplicitLegendLabels = (/"observed",model_name/)
808 ;-----------------------------
810 tRes@txFontHeightF = 0.025
812 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
814 gsn_text_ndc(wks,correlation_text,0.5,0.8,tRes)
816 xy = gsn_csm_xy(wks,xvalues,yvalues,res)
817 ;-------------------------------
818 ;Attach the vertical bar and the horizontal cap line
821 lnres@gsLineColor = line_colors(nd)
824 if(.not.ismissing(mn_yvalues(nd,i)).and. \
825 .not.ismissing(mx_yvalues(nd,i))) then
827 ; Attach the vertical bar, both above and below the marker.
831 y2 = mn_yvalues(nd,i)
832 min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
834 y2 = mx_yvalues(nd,i)
835 max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
837 ; Attach the horizontal cap line, both above and below the marker.
839 x1 = xvalues(nd,i) - dx4
840 x2 = xvalues(nd,i) + dx4
841 y1 = mn_yvalues(nd,i)
842 min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
844 y1 = mx_yvalues(nd,i)
845 max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
852 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
853 ; system("rm "+plot_name+"."+plot_type)
854 ; system("rm "+plot_name+"-1."+plot_type_new)
855 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
858 ;------------------------------------------------------------------------
862 resw = True ; Use plot options
863 resw@cnFillOn = True ; Turn on color fill
864 resw@gsnSpreadColors = True ; use full colormap
865 ; resw@cnFillMode = "RasterFill" ; Turn on raster color
866 ; resw@lbLabelAutoStride = True
867 resw@cnLinesOn = False ; Turn off contourn lines
868 resw@mpFillOn = False ; Turn off map fill
870 resw@gsnSpreadColors = True ; use full colormap
871 resw@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
872 resw@cnMinLevelValF = 0. ; Min level
873 resw@cnMaxLevelValF = 2200. ; Max level
874 resw@cnLevelSpacingF = 200. ; interval
875 ;------------------------------------------------------------------------
878 delta = 0.00000000001
879 nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
881 plot_name = "global_ob"
882 title = "Observed MODIS MOD 17"
883 resw@tiMainString = title
885 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
886 gsn_define_colormap(wks,"gui_default") ; choose colormap
889 plot = gsn_csm_contour_map_ce(wks,nppglobe,resw)
892 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
893 system("rm "+plot_name+"."+plot_type)
894 system("rm "+plot_name+"-1."+plot_type_new)
895 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
898 ;------------------------------------------------------------------------
899 ;global contour model
901 plot_name = "global_model"
902 title = "Model "+ model_name
903 resw@tiMainString = title
905 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
906 gsn_define_colormap(wks,"gui_default") ; choose colormap
908 plot = gsn_csm_contour_map_ce(wks,nppmod,resw)
911 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
912 system("rm "+plot_name+"."+plot_type)
913 system("rm "+plot_name+"-1."+plot_type_new)
914 system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)
917 ;------------------------------------------------------------------------
918 ;global contour model vs ob
920 plot_name = "global_model_vs_ob"
922 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
923 gsn_define_colormap(wks,"gui_default") ; choose colormap
926 plot=new(3,graphic) ; create graphic array
928 resw@gsnFrame = False ; Do not draw plot
929 resw@gsnDraw = False ; Do not advance frame
933 title = "Observed MODIS MOD 17"
934 resw@tiMainString = title
936 plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resw) ; for observed
940 title = "Model "+ model_name
941 resw@tiMainString = title
943 plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resw) ; for model
948 zz = nppmod - nppglobe
949 title = "Model_"+model_name+" - Observed"
951 resw@cnMinLevelValF = -500 ; Min level
952 resw@cnMaxLevelValF = 500. ; Max level
953 resw@cnLevelSpacingF = 50. ; interval
954 resw@tiMainString = title
956 plot(2) = gsn_csm_contour_map_ce(wks,zz,resw)
958 pres = True ; panel plot mods desired
959 pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
960 ; indiv. plots in panel
961 pres@gsnMaximize = True ; fill the page
963 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
965 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
966 system("rm "+plot_name+"."+plot_type)
967 ; system("rm "+plot_name+"-1."+plot_type_new)
968 ; system("mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new)