223 ;********************************************************** |
224 ;********************************************************** |
224 ; Units for these variables are: |
225 ; Units for these variables are: |
225 ; |
226 ; |
226 ; rain81 : mm/year |
227 ; rain81 : mm/year |
227 ; rainmod : mm/s |
228 ; rainmod : mm/s |
228 ; npp81 : g C/m^2/year |
229 ; npp81 : gC/m^2/year |
229 ; nppmod81: g C/m^2/s |
230 ; nppmod81: gC/m^2/s |
230 ; nppglobe: g C/m^2/year |
231 ; nppglobe: gC/m^2/year |
231 ; |
232 ; |
232 ; We want to convert these to "m/year" and "g C/m^2/year". |
233 ; We want to convert these to "m/year" and "gC/m^2/year". |
233 |
234 |
234 nsec_per_year = 60*60*24*365 |
235 nsec_per_year = 60*60*24*365 |
235 |
236 |
236 rain81 = rain81 / 1000. |
237 rain81 = rain81 / 1000. |
237 rainmod81 = (rainmod81/ 1000.) * nsec_per_year |
238 rainmod81 = (rainmod81/ 1000.) * nsec_per_year |
241 rainmod933 = (rainmod933/ 1000.) * nsec_per_year |
242 rainmod933 = (rainmod933/ 1000.) * nsec_per_year |
242 nppmod933 = nppmod933 * nsec_per_year |
243 nppmod933 = nppmod933 * nsec_per_year |
243 |
244 |
244 nppmod = nppmod * nsec_per_year |
245 nppmod = nppmod * nsec_per_year |
245 |
246 |
246 npp81@units = "gC/m^2/yr" |
247 npp81@units = "gC m~S~-2~N~ y~S~-1~N~" |
247 nppmod81@units = "gC/m^2/yr" |
248 nppmod81@units = "gC m~S~-2~N~ y~S~-1~N~" |
248 npp933@units = "gC/m^2/yr" |
249 npp933@units = "gC m~S~-2~N~ y~S~-1~N~" |
249 nppmod933@units = "gC/m^2/yr" |
250 nppmod933@units = "gC m~S~-2~N~ y~S~-1~N~" |
250 nppmod@units = "gC/m^2/yr" |
251 nppmod@units = "gC m~S~-2~N~ y~S~-1~N~" |
251 nppglobe@units = "gC/m^2/yr" |
252 nppglobe@units = "gC m~S~-2~N~ y~S~-1~N~" |
252 rain81@units = "m/yr" |
253 rain81@units = "m y~S~-1~N~" |
253 rainmod81@units = "m/yr" |
254 rainmod81@units = "m y~S~-1~N~" |
254 rain933@units = "m/yr" |
255 rain933@units = "m y~S~-1~N~" |
255 rainmod933@units = "m/yr" |
256 rainmod933@units = "m y~S~-1~N~" |
256 |
257 |
257 npp81@long_name = "Obs. NPP (gC/m2/year)" |
258 npp81@long_name = "Obs. NPP (gC m~S~-2~N~ y~S~-1~N~)" |
258 npp933@long_name = "Obs. NPP (gC/m2/year)" |
259 npp933@long_name = "Obs. NPP (gC m~S~-2~N~ y~S~-1~N~)" |
259 nppmod81@long_name = "Model NPP (gC/m2/year)" |
260 nppmod81@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)" |
260 nppmod933@long_name = "Model NPP (gC/m2/year)" |
261 nppmod933@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)" |
261 nppmod@long_name = "Model NPP (gC/m2/year)" |
262 nppmod@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)" |
262 nppglobe@long_name = "NPP (gC/m2/year)" |
263 nppglobe@long_name = "NPP (gC m~S~-2~N~ y~S~-1~N~)" |
263 rain81@long_name = "PREC (m/year)" |
264 rain81@long_name = "PREC (m y~S~-1~N~)" |
264 rain933@long_name = "PREC (m/year)" |
265 rain933@long_name = "PREC (m y~S~-1~N~)" |
265 rainmod81@long_name = "PREC (m/year)" |
266 rainmod81@long_name = "PREC (m y~S~-1~N~)" |
266 rainmod933@long_name = "PREC (m/year)" |
267 rainmod933@long_name = "PREC (m y~S~-1~N~)" |
267 |
268 |
268 ; change longitude from 0-360 to -180-180 |
269 ; change longitude from 0-360 to -180-180 |
269 x81 = where(x81 .gt. 180., x81 -360., x81 ) |
270 x81 = where(x81 .gt. 180., x81 -360., x81 ) |
270 x933 = where(x933 .gt. 180., x933-360., x933) |
271 x933 = where(x933 .gt. 180., x933-360., x933) |
271 |
272 |
287 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \ |
288 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \ |
288 ,"<tr>" \ |
289 ,"<tr>" \ |
289 ," <th bgcolor=DDDDDD >Site ID</th>" \ |
290 ," <th bgcolor=DDDDDD >Site ID</th>" \ |
290 ," <th bgcolor=DDDDDD >Latitude</th>" \ |
291 ," <th bgcolor=DDDDDD >Latitude</th>" \ |
291 ," <th bgcolor=DDDDDD >Longitude</th>" \ |
292 ," <th bgcolor=DDDDDD >Longitude</th>" \ |
292 ," <th bgcolor=DDDDDD >NPP(gC/m2/year)</th>" \ |
293 ," <th bgcolor=DDDDDD >NPP (gC m<small><sup>-2</sup></small> y<small><sup>-1</sup></small>)</th>" \ |
293 ," <th bgcolor=DDDDDD >PPT(m/year)</th>" \ |
294 ," <th bgcolor=DDDDDD >PPT (m y<small><sup>-1</sup></small>)</th>" \ |
294 ,"</tr>" \ |
295 ,"</tr>" \ |
295 /) |
296 /) |
296 table_footer = "</table>" |
297 table_footer = "</table>" |
297 row_header = "<tr>" |
298 row_header = "<tr>" |
298 row_footer = "</tr>" |
299 row_footer = "</tr>" |
352 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \ |
353 "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \ |
353 ,"<tr>" \ |
354 ,"<tr>" \ |
354 ," <th bgcolor=DDDDDD >Site ID</th>" \ |
355 ," <th bgcolor=DDDDDD >Site ID</th>" \ |
355 ," <th bgcolor=DDDDDD >Latitude</th>" \ |
356 ," <th bgcolor=DDDDDD >Latitude</th>" \ |
356 ," <th bgcolor=DDDDDD >Longitude</th>" \ |
357 ," <th bgcolor=DDDDDD >Longitude</th>" \ |
357 ," <th bgcolor=DDDDDD >NPP(gC/m2/year)</th>" \ |
358 ," <th bgcolor=DDDDDD >NPP (gC m<small><sup>-2</sup></small> y<small><sup>-1</sup></small>)</th>" \ |
358 ," <th bgcolor=DDDDDD >PPT(m/year)</th>" \ |
359 ," <th bgcolor=DDDDDD >PPT (m y<small><sup>-1</sup></small>)</th>" \ |
359 ,"</tr>" \ |
360 ,"</tr>" \ |
360 /) |
361 /) |
361 table_footer = "</table>" |
362 table_footer = "</table>" |
362 row_header = "<tr>" |
363 row_header = "<tr>" |
363 row_footer = "</tr>" |
364 row_footer = "</tr>" |
423 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \ |
424 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \ |
424 ,"<tr>" \ |
425 ,"<tr>" \ |
425 ," <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \ |
426 ," <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \ |
426 ," <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \ |
427 ," <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \ |
427 ," <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \ |
428 ," <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \ |
428 ," <th bgcolor=DDDDDD colspan=2>NPP(gC/m2.year)</th>" \ |
429 ," <th bgcolor=DDDDDD colspan=2>NPP (gC m<small><sup>-2</sup></small> y<small><sup>-1</sup></small>)</th>" \ |
429 ," <th bgcolor=DDDDDD colspan=2>RAIN(m/year)</th>" \ |
430 ," <th bgcolor=DDDDDD colspan=2>RAIN (m y<small><sup>-1</sup></small>)</th>" \ |
430 ,"</tr>" \ |
431 ,"</tr>" \ |
431 ,"<tr>" \ |
432 ,"<tr>" \ |
432 ," <th bgcolor=DDDDDD >observed</th>" \ |
433 ," <th bgcolor=DDDDDD >observed</th>" \ |
433 , table_header_text \ |
434 , table_header_text \ |
434 ," <th bgcolor=DDDDDD >observed</th>" \ |
435 ," <th bgcolor=DDDDDD >observed</th>" \ |
504 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \ |
505 "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \ |
505 ,"<tr>" \ |
506 ,"<tr>" \ |
506 ," <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \ |
507 ," <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \ |
507 ," <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \ |
508 ," <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \ |
508 ," <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \ |
509 ," <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \ |
509 ," <th bgcolor=DDDDDD colspan=2>NPP(gC/m2.year)</th>" \ |
510 ," <th bgcolor=DDDDDD colspan=2>NPP (gC m<small><sup>-2</sup></small> y<small><sup>-1</sup></small>)</th>" \ |
510 ," <th bgcolor=DDDDDD colspan=2>RAIN(m/year)</th>" \ |
511 ," <th bgcolor=DDDDDD colspan=2>RAIN (m y<small><sup>-1</sup></small>)</th>" \ |
511 ,"</tr>" \ |
512 ,"</tr>" \ |
512 ,"<tr>" \ |
513 ,"<tr>" \ |
513 ," <th bgcolor=DDDDDD >observed</th>" \ |
514 ," <th bgcolor=DDDDDD >observed</th>" \ |
514 , table_header_text \ |
515 , table_header_text \ |
515 ," <th bgcolor=DDDDDD >observed</th>" \ |
516 ," <th bgcolor=DDDDDD >observed</th>" \ |
567 ;*************************************************************************** |
568 ;*************************************************************************** |
568 ;(A)-5 scatter plot, model vs ob 81 |
569 ;(A)-5 scatter plot, model vs ob 81 |
569 ;*************************************************************************** |
570 ;*************************************************************************** |
570 |
571 |
571 plot_name = "scatter_model_vs_ob_81" |
572 plot_name = "scatter_model_vs_ob_81" |
572 title = model_name +" vs Class A observations (81 sites)" |
573 ;title = model_name +" vs Class A observations (81 sites)" |
|
574 title = model_name +" (1975-2000) vs Class A observations (81 sites)" |
573 |
575 |
574 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
576 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
575 res = True ; plot mods desired |
577 res = True ; plot mods desired |
576 res@tiMainString = title ; add title |
578 res@tiMainString = title ; add title |
577 res@xyMarkLineModes = "Markers" ; choose which have markers |
579 res@xyMarkLineModes = "Markers" ; choose which have markers |
585 ;------------------------------- |
587 ;------------------------------- |
586 ;compute correlation coef. and M |
588 ;compute correlation coef. and M |
587 ccr81 = esccr(nppmod81,npp81,0) |
589 ccr81 = esccr(nppmod81,npp81,0) |
588 ;print (ccr81) |
590 ;print (ccr81) |
589 |
591 |
590 score_max = 2.5 |
592 score_max = 1.0 |
591 |
593 |
592 bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81))) |
594 bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81))) |
593 M81s = (1. - (bias/dimsizes(y81)))*score_max |
595 M81s = (1. - (bias/dimsizes(y81)))*score_max |
594 M_npp_S81 = sprintf("%.2f", M81s) |
596 M_npp_S81 = sprintf("%.2f", M81s) |
595 |
597 |
616 ;------------------------------- |
618 ;------------------------------- |
617 draw (plot) |
619 draw (plot) |
618 delete (wks) |
620 delete (wks) |
619 delete (dum) |
621 delete (dum) |
620 |
622 |
621 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
623 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
622 "rm "+plot_name+"."+plot_type) |
624 "rm "+plot_name+"."+plot_type) |
|
625 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
623 |
626 |
624 ;*************************************************************************** |
627 ;*************************************************************************** |
625 ;(A)-6 scatter plot, model vs ob 933 |
628 ;(A)-6 scatter plot, model vs ob 933 |
626 ;*************************************************************************** |
629 ;*************************************************************************** |
627 |
630 |
628 plot_name = "scatter_model_vs_ob_933" |
631 plot_name = "scatter_model_vs_ob_933" |
629 title = model_name +" vs Class B Observations (933 sites)" |
632 title = model_name +" (1975-2000) vs Class B Observations (933 sites)" |
630 |
633 |
631 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
634 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
632 res = True ; plot mods desired |
635 res = True ; plot mods desired |
633 res@tiMainString = title ; add title |
636 res@tiMainString = title ; add title |
634 res@xyMarkLineModes = "Markers" ; choose which have markers |
637 res@xyMarkLineModes = "Markers" ; choose which have markers |
642 ;------------------------------- |
645 ;------------------------------- |
643 ;compute correlation coef. and M |
646 ;compute correlation coef. and M |
644 |
647 |
645 ccr933 = esccr(nppmod933,npp933,0) |
648 ccr933 = esccr(nppmod933,npp933,0) |
646 |
649 |
647 score_max = 2.5 |
650 score_max = 1.0 |
648 |
651 |
649 bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933))) |
652 bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933))) |
650 M933s = (1. - (bias/dimsizes(y933)))*score_max |
653 M933s = (1. - (bias/dimsizes(y933)))*score_max |
651 M_npp_S933 = sprintf("%.2f", M933s) |
654 M_npp_S933 = sprintf("%.2f", M933s) |
652 |
655 |
740 resh@xyMarker = 16 |
744 resh@xyMarker = 16 |
741 resh@xyMarkerColors = (/"Brown","Blue"/) |
745 resh@xyMarkerColors = (/"Brown","Blue"/) |
742 resh@trYMinF = min(mn_yvalues) - 10. |
746 resh@trYMinF = min(mn_yvalues) - 10. |
743 resh@trYMaxF = max(mx_yvalues) + 10. |
747 resh@trYMaxF = max(mx_yvalues) + 10. |
744 |
748 |
745 resh@tiYAxisString = "NPP (g C/m2/year)" |
749 resh@tiYAxisString = "NPP (gC m~S~-2~N~ y~S~-1~N~)" |
746 resh@tiXAxisString = "Precipitation (m/year)" |
750 resh@tiXAxisString = "Precipitation (m y~S~-1~N~)" |
747 |
751 |
748 max_bar = new((/2,nx/),graphic) |
752 max_bar = new((/2,nx/),graphic) |
749 min_bar = new((/2,nx/),graphic) |
753 min_bar = new((/2,nx/),graphic) |
750 max_cap = new((/2,nx/),graphic) |
754 max_cap = new((/2,nx/),graphic) |
751 min_cap = new((/2,nx/),graphic) |
755 min_cap = new((/2,nx/),graphic) |
762 resh@tiMainString = title |
766 resh@tiMainString = title |
763 |
767 |
764 wks = gsn_open_wks (plot_type,plot_name) |
768 wks = gsn_open_wks (plot_type,plot_name) |
765 |
769 |
766 xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh) |
770 xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh) |
|
771 |
|
772 ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009 |
|
773 print("Observations: " + xvalues(0,:) + " " + yvalues(0,:)) |
767 |
774 |
768 ;------------------------------- |
775 ;------------------------------- |
769 ;Attach the vertical bar and the horizontal cap line |
776 ;Attach the vertical bar and the horizontal cap line |
770 |
777 |
771 do nd=0,0 |
778 do nd=0,0 |
807 end do |
814 end do |
808 |
815 |
809 draw(xy) |
816 draw(xy) |
810 delete (wks) |
817 delete (wks) |
811 |
818 |
812 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
819 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
813 "rm "+plot_name+"."+plot_type) |
820 "rm "+plot_name+"."+plot_type) |
|
821 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
814 |
822 |
815 ;**************************************************************************** |
823 ;**************************************************************************** |
816 ;(B)-2 histogram plot, model vs ob 81 site |
824 ;(B)-2 histogram plot, model vs ob 81 site |
817 ;**************************************************************************** |
825 ;**************************************************************************** |
818 |
826 |
819 plot_name = "histogram_model_vs_ob_81" |
827 plot_name = "histogram_model_vs_ob_81" |
820 title = model_name+ " vs Class A Observations (81 sites)" |
828 title = model_name+ " (1975-2000) vs Class A Observations (81 sites)" |
821 resh@tiMainString = title |
829 resh@tiMainString = title |
822 |
830 |
823 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
831 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
824 |
832 |
825 ;----------------------------- |
833 ;----------------------------- |
845 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")" |
853 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")" |
846 |
854 |
847 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
855 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
848 |
856 |
849 xy = gsn_csm_xy(wks,xvalues,yvalues,resh) |
857 xy = gsn_csm_xy(wks,xvalues,yvalues,resh) |
|
858 |
|
859 ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009 |
|
860 print(model_name + ": " + xvalues(1,:) + " " + yvalues(1,:)) |
|
861 ;print("All: " + xvalues(:,:) + " " + yvalues(:,:)) |
850 ;------------------------------- |
862 ;------------------------------- |
851 ;Attach the vertical bar and the horizontal cap line |
863 ;Attach the vertical bar and the horizontal cap line |
852 |
864 |
853 do nd=0,1 |
865 do nd=0,1 |
854 lnres@gsLineColor = line_colors(nd) |
866 lnres@gsLineColor = line_colors(nd) |
965 resh@xyMarker = 16 |
978 resh@xyMarker = 16 |
966 resh@xyMarkerColors = (/"Brown","Blue"/) |
979 resh@xyMarkerColors = (/"Brown","Blue"/) |
967 resh@trYMinF = min(mn_yvalues) - 10. |
980 resh@trYMinF = min(mn_yvalues) - 10. |
968 resh@trYMaxF = max(mx_yvalues) + 10. |
981 resh@trYMaxF = max(mx_yvalues) + 10. |
969 |
982 |
970 resh@tiYAxisString = "NPP (g C/m2/year)" |
983 resh@tiYAxisString = "NPP (gC m~S~-2~N~ y~S~-1~N~)" |
971 resh@tiXAxisString = "Precipitation (m/year)" |
984 resh@tiXAxisString = "Precipitation (m y~S~-1~N~)" |
972 |
985 |
973 max_bar = new((/2,nx/),graphic) |
986 max_bar = new((/2,nx/),graphic) |
974 min_bar = new((/2,nx/),graphic) |
987 min_bar = new((/2,nx/),graphic) |
975 max_cap = new((/2,nx/),graphic) |
988 max_cap = new((/2,nx/),graphic) |
976 min_cap = new((/2,nx/),graphic) |
989 min_cap = new((/2,nx/),graphic) |
987 resh@tiMainString = title |
1000 resh@tiMainString = title |
988 |
1001 |
989 wks = gsn_open_wks (plot_type,plot_name) |
1002 wks = gsn_open_wks (plot_type,plot_name) |
990 |
1003 |
991 xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh) |
1004 xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh) |
|
1005 |
|
1006 ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009 |
|
1007 print("Observations: " + xvalues(0,:) + " " + yvalues(0,:)) |
992 |
1008 |
993 ;------------------------------- |
1009 ;------------------------------- |
994 ;Attach the vertical bar and the horizontal cap line |
1010 ;Attach the vertical bar and the horizontal cap line |
995 |
1011 |
996 do nd=0,0 |
1012 do nd=0,0 |
1025 |
1041 |
1026 draw(xy) |
1042 draw(xy) |
1027 delete (xy) |
1043 delete (xy) |
1028 delete (wks) |
1044 delete (wks) |
1029 |
1045 |
1030 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1046 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1031 "rm "+plot_name+"."+plot_type) |
1047 "rm "+plot_name+"."+plot_type) |
|
1048 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
1032 |
1049 |
1033 ;**************************************************************************** |
1050 ;**************************************************************************** |
1034 ;(B)-4 histogram plot, model vs ob 933 site |
1051 ;(B)-4 histogram plot, model vs ob 933 site |
1035 ;**************************************************************************** |
1052 ;**************************************************************************** |
1036 |
1053 |
1037 plot_name = "histogram_model_vs_ob_933" |
1054 plot_name = "histogram_model_vs_ob_933" |
1038 title = model_name+ " vs Class B Observations (933 sites)" |
1055 title = model_name+ " (1975-2000) vs Class B Observations (933 sites)" |
1039 resh@tiMainString = title |
1056 resh@tiMainString = title |
1040 |
1057 |
1041 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
1058 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
1042 |
1059 |
1043 ;----------------------------- |
1060 ;----------------------------- |
1063 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")" |
1080 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")" |
1064 |
1081 |
1065 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
1082 gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes) |
1066 |
1083 |
1067 xy = gsn_csm_xy(wks,xvalues,yvalues,resh) |
1084 xy = gsn_csm_xy(wks,xvalues,yvalues,resh) |
|
1085 |
|
1086 ; Added by Forrest Hoffman to print out values on Wed Feb 4 14:36:00 EST 2009 |
|
1087 print("Observations: " + xvalues(1,:) + " " + yvalues(1,:)) |
1068 ;------------------------------- |
1088 ;------------------------------- |
1069 ;Attach the vertical bar and the horizontal cap line |
1089 ;Attach the vertical bar and the horizontal cap line |
1070 |
1090 |
1071 do nd=0,1 |
1091 do nd=0,1 |
1072 lnres@gsLineColor = line_colors(nd) |
1092 lnres@gsLineColor = line_colors(nd) |
1100 |
1120 |
1101 draw(xy) |
1121 draw(xy) |
1102 delete(xy) |
1122 delete(xy) |
1103 delete(wks) |
1123 delete(wks) |
1104 |
1124 |
1105 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1125 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1106 "rm "+plot_name+"."+plot_type) |
1126 "rm "+plot_name+"."+plot_type) |
|
1127 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
|
1128 |
|
1129 ;*************************************************************************** |
|
1130 ; Read 2000-2004 dataset for MODIS comparison |
|
1131 ;*************************************************************************** |
|
1132 |
|
1133 ;------------------------------------------------------ |
|
1134 ; read model data |
|
1135 |
|
1136 fm = addfile (dirm+film9,"r") |
|
1137 |
|
1138 nppmod0 = fm->NPP |
|
1139 xm = fm->lon |
|
1140 ym = fm->lat |
|
1141 |
|
1142 delete (fm) |
|
1143 |
|
1144 nppmod = nppmod0(0,:,:) |
|
1145 delete (nppmod0) |
|
1146 nppmod = nppmod * nsec_per_year |
|
1147 nppmod@units = "gC m~S~-2~N~ y~S~-1~N~" |
|
1148 nppmod@long_name = "Model NPP (gC m~S~-2~N~ y~S~-1~N~)" |
1107 |
1149 |
1108 ;*************************************************************************** |
1150 ;*************************************************************************** |
1109 ;(C) global contour |
1151 ;(C) global contour |
1110 ;*************************************************************************** |
1152 ;*************************************************************************** |
1111 |
1153 |
1130 |
1172 |
1131 delta = 0.00000000001 |
1173 delta = 0.00000000001 |
1132 nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe) |
1174 nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe) |
1133 |
1175 |
1134 plot_name = "global_ob" |
1176 plot_name = "global_ob" |
1135 title = "Observed MODIS MOD 17" |
1177 title = "MODIS MOD17A3 (2000-2004)" |
1136 resg@tiMainString = title |
1178 resg@tiMainString = title |
1137 |
1179 |
1138 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
1180 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
1139 gsn_define_colormap(wks,"gui_default") ; choose colormap |
1181 gsn_define_colormap(wks,"gui_default") ; choose colormap |
1140 |
1182 |
1141 plot = gsn_csm_contour_map_ce(wks,nppglobe,resg) |
1183 plot = gsn_csm_contour_map_ce(wks,nppglobe,resg) |
1142 |
1184 |
1143 delete (wks) |
1185 delete (wks) |
1144 |
1186 |
1145 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1187 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1146 "rm "+plot_name+"."+plot_type) |
1188 "rm "+plot_name+"."+plot_type) |
|
1189 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
1147 ;**************************************************************************** |
1190 ;**************************************************************************** |
1148 ;(C)-2 global contour plot, model |
1191 ;(C)-2 global contour plot, model |
1149 ;**************************************************************************** |
1192 ;**************************************************************************** |
1150 |
1193 |
1151 plot_name = "global_model" |
1194 plot_name = "global_model" |
1152 title = "Model "+ model_name |
1195 title = "Model "+ model_name + " (2000-2004)" |
1153 resg@tiMainString = title |
1196 resg@tiMainString = title |
1154 |
1197 |
1155 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
1198 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
1156 gsn_define_colormap(wks,"gui_default") ; choose colormap |
1199 gsn_define_colormap(wks,"gui_default") ; choose colormap |
1157 |
1200 |
1158 plot = gsn_csm_contour_map_ce(wks,nppmod,resg) |
1201 plot = gsn_csm_contour_map_ce(wks,nppmod,resg) |
1159 |
1202 |
1160 delete (wks) |
1203 delete (wks) |
1161 |
1204 |
1162 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1205 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1163 "rm "+plot_name+"."+plot_type) |
1206 "rm "+plot_name+"."+plot_type) |
|
1207 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
1164 |
1208 |
1165 ;**************************************************************************** |
1209 ;**************************************************************************** |
1166 ;(C)-3 global contour plot, model vs ob |
1210 ;(C)-3 global contour plot, model vs ob |
1167 ;**************************************************************************** |
1211 ;**************************************************************************** |
1168 |
1212 |
1213 |
1257 |
1214 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
1258 gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) |
1215 ;-------------------------------------------------------------------- |
1259 ;-------------------------------------------------------------------- |
1216 ;(a) ob |
1260 ;(a) ob |
1217 |
1261 |
1218 title = "Observed MODIS MOD 17" |
1262 title = "MODIS MOD17A3 (2000-2004)" |
1219 resg@tiMainString = title |
1263 resg@tiMainString = title |
1220 |
1264 |
1221 plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg) |
1265 plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg) |
1222 |
1266 |
1223 ;(b) model |
1267 ;(b) model |
1224 |
1268 |
1225 title = "Model "+ model_name |
1269 title = "Model "+ model_name + " (2000-2004)" |
1226 resg@tiMainString = title |
1270 resg@tiMainString = title |
1227 |
1271 |
1228 plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg) |
1272 plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg) |
1229 |
1273 |
1230 ;(c) model-ob |
1274 ;(c) model-ob |
1231 |
1275 |
1232 zz = nppmod |
1276 zz = nppmod |
1233 zz = nppmod - nppglobe |
1277 zz = nppmod - nppglobe |
1234 title = "Model_"+model_name+" - Observed" |
1278 title = "Model "+model_name+" - MODIS MOD17A3" |
1235 |
1279 |
1236 resg@cnMinLevelValF = -500 ; Min level |
1280 resg@cnMinLevelValF = -500 ; Min level |
1237 resg@cnMaxLevelValF = 500. ; Max level |
1281 resg@cnMaxLevelValF = 500. ; Max level |
1238 resg@cnLevelSpacingF = 50. ; interval |
1282 resg@cnLevelSpacingF = 50. ; interval |
1239 resg@tiMainString = title |
1283 resg@tiMainString = title |
1248 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
1292 gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot |
1249 |
1293 |
1250 delete (wks) |
1294 delete (wks) |
1251 delete (plot) |
1295 delete (plot) |
1252 |
1296 |
1253 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1297 system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1254 "rm "+plot_name+"."+plot_type) |
1298 "rm "+plot_name+"."+plot_type) |
|
1299 ;system("convert -trim -density 150 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
1255 |
1300 |
1256 ;*************************************************************************** |
1301 ;*************************************************************************** |
1257 ;(D)-1 zonal line plot, ob |
1302 ;(D)-1 zonal line plot, ob |
1258 ;*************************************************************************** |
1303 ;*************************************************************************** |
1259 |
1304 |
1260 vv = zonalAve(nppglobe) |
1305 vv = zonalAve(nppglobe) |
1261 vv@long_name = nppglobe@long_name |
1306 vv@long_name = nppglobe@long_name |
1262 |
1307 |
1263 plot_name = "zonal_ob" |
1308 plot_name = "zonal_ob" |
1264 title = "Observed MODIS MOD 17" |
1309 title = "MODIS MOD17A3 (2000-2004)" |
1265 |
1310 |
1266 resz = True |
1311 resz = True |
1267 resz@tiMainString = title |
1312 resz@tiMainString = title |
1268 |
1313 |
1269 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
1314 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
1270 |
1315 |
1271 plot = gsn_csm_xy (wks,ym,vv,resz) |
1316 plot = gsn_csm_xy (wks,ym,vv,resz) |
1272 |
1317 |
1273 delete (wks) |
1318 delete (wks) |
1274 |
1319 |
1275 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1320 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1276 "rm "+plot_name+"."+plot_type) |
1321 "rm "+plot_name+"."+plot_type) |
|
1322 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
1277 |
1323 |
1278 ;**************************************************************************** |
1324 ;**************************************************************************** |
1279 ;(D)-2 zonal line plot, model vs ob |
1325 ;(D)-2 zonal line plot, model vs ob |
1280 ;**************************************************************************** |
1326 ;**************************************************************************** |
1281 |
1327 |
1286 |
1332 |
1287 s@long_name = nppglobe@long_name |
1333 s@long_name = nppglobe@long_name |
1288 ;------------------------------------------- |
1334 ;------------------------------------------- |
1289 ; compute correlation coef and M score |
1335 ; compute correlation coef and M score |
1290 |
1336 |
1291 score_max = 5.0 |
1337 score_max = 2.0 |
1292 |
1338 |
1293 ccrZ = esccr(s(1,:), s(0,:),0) |
1339 ccrZ = esccr(s(1,:), s(0,:),0) |
1294 MZ = (ccrZ*ccrZ)* score_max |
1340 MZ = (ccrZ*ccrZ)* score_max |
1295 M_npp_Z = sprintf("%.2f", MZ) |
1341 M_npp_Z = sprintf("%.2f", MZ) |
1296 |
1342 |
1301 |
1347 |
1302 system("sed s#M_npp_Z#"+M_npp_Z+"# "+html_name+" > "+html_new+";"+ \ |
1348 system("sed s#M_npp_Z#"+M_npp_Z+"# "+html_name+" > "+html_new+";"+ \ |
1303 "mv -f "+html_new+" "+html_name) |
1349 "mv -f "+html_new+" "+html_name) |
1304 ;------------------------------------------- |
1350 ;------------------------------------------- |
1305 plot_name = "zonal_model_vs_ob" |
1351 plot_name = "zonal_model_vs_ob" |
1306 title = "Zonal Average" |
1352 title = "Zonal Average (2000-2004)" |
1307 resz@tiMainString = title |
1353 resz@tiMainString = title |
1308 |
1354 |
1309 wks = gsn_open_wks (plot_type,plot_name) |
1355 wks = gsn_open_wks (plot_type,plot_name) |
1310 |
1356 |
1311 ; resz@vpHeightF = 0.4 ; change aspect ratio of plot |
1357 ; resz@vpHeightF = 0.4 ; change aspect ratio of plot |
1322 |
1368 |
1323 ; Legent |
1369 ; Legent |
1324 resz@pmLegendDisplayMode = "Always" ; turn on legend |
1370 resz@pmLegendDisplayMode = "Always" ; turn on legend |
1325 resz@pmLegendSide = "Top" ; Change location of |
1371 resz@pmLegendSide = "Top" ; Change location of |
1326 ; resz@pmLegendParallelPosF = .45 ; move units right |
1372 ; resz@pmLegendParallelPosF = .45 ; move units right |
1327 resz@pmLegendParallelPosF = .82 ; move units right |
1373 resz@pmLegendParallelPosF = .75 ; move units right |
1328 resz@pmLegendOrthogonalPosF = -0.4 ; move units down |
1374 resz@pmLegendOrthogonalPosF = -0.4 ; move units down |
1329 |
1375 |
1330 resz@pmLegendWidthF = 0.10 ; Change width and |
1376 resz@pmLegendWidthF = 0.10 ; Change width and |
1331 resz@pmLegendHeightF = 0.10 ; height of legend. |
1377 resz@pmLegendHeightF = 0.10 ; height of legend. |
1332 resz@lgLabelFontHeightF = .02 ; change font height |
1378 resz@lgLabelFontHeightF = .015 ; change font height |
1333 ; resz@lgTitleOn = True ; turn on legend title |
1379 ; resz@lgTitleOn = True ; turn on legend title |
1334 ; resz@lgTitleString = "Example" ; create legend title |
1380 ; resz@lgTitleString = "Example" ; create legend title |
1335 ; resz@lgTitleFontHeightF = .025 ; font of legend title |
1381 ; resz@lgTitleFontHeightF = .025 ; font of legend title |
1336 resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels |
1382 resz@xyExplicitLegendLabels = (/"MODIS MOD17A3",model_name/) ; explicit labels |
1337 ;-------------------------------------------------------------------- |
1383 ;-------------------------------------------------------------------- |
1338 zRes = True |
1384 zRes = True |
1339 zRes@txFontHeightF = 0.025 |
1385 zRes@txFontHeightF = 0.025 |
1340 |
1386 |
1341 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")" |
1387 correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")" |
1345 |
1391 |
1346 plot = gsn_csm_xy (wks,ym,s,resz) |
1392 plot = gsn_csm_xy (wks,ym,s,resz) |
1347 |
1393 |
1348 delete (wks) |
1394 delete (wks) |
1349 |
1395 |
1350 system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1396 system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ |
1351 "rm "+plot_name+"."+plot_type) |
1397 "rm "+plot_name+"."+plot_type) |
|
1398 ;system("convert -trim -density 100 "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new) |
1352 |
1399 |
1353 ;*************************************************************************** |
1400 ;*************************************************************************** |
1354 ; add total score and write to file |
1401 ; add total score and write to file |
1355 ;*************************************************************************** |
1402 ;*************************************************************************** |
1356 M_total = M81s + M81h + M933s + M933h + MG + MZ |
1403 M_total = M81s + M81h + M933s + M933h + MG + MZ |