diff -r 000000000000 -r 0c6405ab2ff4 npp/99.all.ncl
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/npp/99.all.ncl Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,1358 @@
+;*****************************************************
+; combine scatter, histogram, global and zonal plots
+; compute all correlation coef and M score
+; add 1-to-1 line in scatter plots
+;*****************************************************
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
+;****************************************************************************
+procedure set_line(lines:string,nline:integer,newlines:string)
+begin
+; add line to ascci/html file
+
+ nnewlines = dimsizes(newlines)
+ if(nline+nnewlines-1.ge.dimsizes(lines))
+ print("set_line: bad index, not setting anything.")
+ return
+ end if
+ lines(nline:nline+nnewlines-1) = newlines
+; print ("lines = " + lines(nline:nline+nnewlines-1))
+ nline = nline + nnewlines
+ return
+end
+;****************************************************************************
+procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \
+ ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \
+ ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \
+ ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \
+ ,dx4[1]:numeric \
+ )
+begin
+; Calculaee "nice" bins for binning the data in equally spaced ranges.
+; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D
+; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4
+
+ nbins = 15 ; Number of bins to use.
+
+ nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
+ nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
+ range = fspan(nicevals(0),nicevals(1),nvals)
+
+; Use this range information to grab all the values in a
+; particular range, and then take an average.
+
+ nr = dimsizes(range)
+ nx = nr-1
+; print (nx)
+
+; xvalues = new((/2,nx/),typeof(RAIN1_1D))
+ xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
+ dx = xvalues(0,1) - xvalues(0,0) ; range width
+ dx4 = dx/4 ; 1/4 of the range
+ xvalues(1,:) = xvalues(0,:) - dx/5.
+; yvalues = new((/2,nx/),typeof(RAIN1_1D))
+; mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
+; mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
+
+ do nd=0,1
+
+; See if we are doing model or observational data.
+
+ if(nd.eq.0) then
+ data = RAIN1_1D
+ npp_data = NPP1_1D
+ else
+ data = RAIN2_1D
+ npp_data = NPP2_1D
+ end if
+
+; Loop through each range and check for values.
+
+ do i=0,nr-2
+ if (i.ne.(nr-2)) then
+; print("")
+; print("In range ["+range(i)+","+range(i+1)+")")
+ idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
+ else
+; print("")
+; print("In range ["+range(i)+",)")
+ idx = ind(range(i).le.data)
+ end if
+
+; Calculate average, and get min and max.
+
+ if(.not.any(ismissing(idx))) then
+ yvalues(nd,i) = avg(npp_data(idx))
+ mn_yvalues(nd,i) = min(npp_data(idx))
+ mx_yvalues(nd,i) = max(npp_data(idx))
+ count = dimsizes(idx)
+ else
+ count = 0
+ yvalues(nd,i) = yvalues@_FillValue
+ mn_yvalues(nd,i) = yvalues@_FillValue
+ mx_yvalues(nd,i) = yvalues@_FillValue
+ end if
+
+; Print out information.
+; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
+; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
+
+; Clean up for next time in loop.
+
+ delete(idx)
+ end do
+ delete(data)
+ delete(npp_data)
+ end do
+end
+;****************************************************************************
+; Main code.
+;****************************************************************************
+begin
+
+ plot_type = "ps"
+ plot_type_new = "png"
+
+;************************************************
+; read model data
+;************************************************
+
+;film = "i01.06cn_1798-2004_ANN_climo.nc"
+;model_name = "06cn"
+;model_grid = "T42"
+
+;film = "i01.06casa_1798-2004_ANN_climo.nc"
+;model_name = "06casa"
+;model_grid = "T42"
+
+;film = "i01.10casa_1948-2004_ANN_climo.nc"
+;model_name = "10casa"
+;model_grid = "T42"
+
+ film = "i01.10cn_1948-2004_ANN_climo.nc"
+ model_name = "10cn"
+ model_grid = "T42"
+
+ html_name = "table.html." + model_name
+ html_new = html_name +".new"
+ system("sed s#model_name#"+model_name+"# "+html_name+" > "+html_new+";"+ \
+ "mv -f "+html_new+" "+html_name)
+;--------------------------------------------------
+ dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ fm = addfile (dirm+film,"r")
+
+ nppmod0 = fm->NPP
+ rainmod0 = fm->RAIN
+ xm = fm->lon
+ ym = fm->lat
+
+;************************************************
+; read ob data
+;************************************************
+
+;(1) data at 81 sites
+ dir81 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
+ fil81 = "data.81.nc"
+ f81 = addfile (dir81+fil81,"r")
+
+ id81 = f81->SITE_ID
+ npp81 = f81->TNPP_C
+ rain81 = tofloat(f81->PREC_ANN)
+ x81 = f81->LONG_DD
+ y81 = f81->LAT_DD
+
+ id81@long_name = "SITE_ID"
+
+;change longitude from (-180,180) to (0,360)
+;for model data interpolation
+
+ do i= 0,dimsizes(x81)-1
+ if (x81(i) .lt. 0.) then
+ x81(i) = x81(i)+ 360.
+ end if
+ end do
+;print (x81)
+;-------------------------------------
+;(2) data at 933 sites
+ dir933 = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
+ fil933 = "data.933.nc"
+ f933 = addfile (dir933+fil933,"r")
+
+ id933 = f933->SITE_ID
+ npp933 = f933->TNPP_C
+ rain933 = f933->PREC
+ x933 = f933->LONG_DD
+ y933 = f933->LAT_DD
+
+ id933@long_name = "SITE_ID"
+
+;change longitude from (-180,180) to (0,360)
+;for model data interpolation
+
+ do i= 0,dimsizes(x933)-1
+ if (x933(i) .lt. 0.) then
+ x933(i) = x933(i)+ 360.
+ end if
+ end do
+;print (x933)
+;----------------------------------------
+;(3) global data, interpolated into model grid
+ dirglobe = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
+ filglobe = "Npp_"+model_grid+"_mean.nc"
+ fglobe = addfile (dirglobe+filglobe,"r")
+
+ nppglobe0 = fglobe->NPP
+ nppglobe = nppglobe0
+
+;***********************************************************************
+; interpolate model data into ob sites
+;***********************************************************************
+
+ nppmod = nppmod0(0,:,:)
+ rainmod = rainmod0(0,:,:)
+ delete (nppmod0)
+ delete (rainmod0)
+
+ nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0)
+
+ nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)
+
+ rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)
+
+ rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)
+
+;**********************************************************
+; unified units
+;**********************************************************
+; Units for these variables are:
+;
+; rain81 : mm/year
+; rainmod : mm/s
+; npp81 : g C/m^2/year
+; nppmod81: g C/m^2/s
+; nppglobe: g C/m^2/year
+;
+; We want to convert these to "m/year" and "g C/m^2/year".
+
+ nsec_per_year = 60*60*24*365
+
+ rain81 = rain81 / 1000.
+ rainmod81 = (rainmod81/ 1000.) * nsec_per_year
+ nppmod81 = nppmod81 * nsec_per_year
+
+ rain933 = rain933 / 1000.
+ rainmod933 = (rainmod933/ 1000.) * nsec_per_year
+ nppmod933 = nppmod933 * nsec_per_year
+
+ nppmod = nppmod * nsec_per_year
+
+ npp81@units = "gC/m^2/yr"
+ nppmod81@units = "gC/m^2/yr"
+ npp933@units = "gC/m^2/yr"
+ nppmod933@units = "gC/m^2/yr"
+ nppmod@units = "gC/m^2/yr"
+ nppglobe@units = "gC/m^2/yr"
+ rain81@units = "m/yr"
+ rainmod81@units = "m/yr"
+ rain933@units = "m/yr"
+ rainmod933@units = "m/yr"
+
+ npp81@long_name = "NPP (gC/m2/year)"
+ npp933@long_name = "NPP (gC/m2/year)"
+ nppmod81@long_name = "NPP (gC/m2/year)"
+ nppmod933@long_name = "NPP (gC/m2/year)"
+ nppmod@long_name = "NPP (gC/m2/year)"
+ nppglobe@long_name = "NPP (gC/m2/year)"
+ rain81@long_name = "PREC (m/year)"
+ rain933@long_name = "PREC (m/year)"
+ rainmod81@long_name = "PREC (m/year)"
+ rainmod933@long_name = "PREC (m/year)"
+
+;*******************************************************************
+;(A)-1 html table of site81 -- observed
+;*******************************************************************
+ output_html = "table_site81_ob.html"
+
+ header = (/"" \
+ ,"
" \
+ ,"CLAMP metrics" \
+ ,"" \
+ ,"Observation at 81 sites
" \
+ /)
+ footer = ""
+
+ table_header = (/ \
+ "" \
+ ,"" \
+ ," Site ID | " \
+ ," Latitude | " \
+ ," Longitude | " \
+ ," NPP(gC/m2.year) | " \
+ ," RAIN(m/year) | " \
+ ,"
" \
+ /)
+ table_footer = "
"
+ row_header = ""
+ row_footer = "
"
+
+ lines = new(50000,string)
+ nline = 0
+
+ set_line(lines,nline,header)
+ set_line(lines,nline,table_header)
+;-----------------------------------------------
+;row of table
+
+ nrow = dimsizes(id81)
+ do n = 0,nrow-1
+ set_line(lines,nline,row_header)
+
+ txt1 = sprintf("%5.0f", id81(n))
+ txt2 = sprintf("%5.2f", y81(n))
+ txt3 = sprintf("%5.2f", x81(n))
+ txt4 = sprintf("%5.2f", npp81(n))
+ txt5 = sprintf("%5.2f", rain81(n))
+
+ set_line(lines,nline,""+txt1+" | ")
+ set_line(lines,nline,""+txt2+" | ")
+ set_line(lines,nline,""+txt3+" | ")
+ set_line(lines,nline,""+txt4+" | ")
+ set_line(lines,nline,""+txt5+" | ")
+
+ set_line(lines,nline,row_footer)
+ end do
+;-----------------------------------------------
+ set_line(lines,nline,table_footer)
+ set_line(lines,nline,footer)
+
+; Now write to an HTML file.
+ idx = ind(.not.ismissing(lines))
+ if(.not.any(ismissing(idx))) then
+ asciiwrite(output_html,lines(idx))
+ else
+ print ("error?")
+ end if
+;*******************************************************************
+;(A)-2 html table of site933 -- observed
+;*******************************************************************
+ output_html = "table_site933_ob.html"
+
+ header = (/"" \
+ ,"" \
+ ,"CLAMP metrics" \
+ ,"" \
+ ,"Observation at 933 sites
" \
+ /)
+ footer = ""
+
+ table_header = (/ \
+ "" \
+ ,"" \
+ ," Site ID | " \
+ ," Latitude | " \
+ ," Longitude | " \
+ ," NPP(gC/m2.year) | " \
+ ," RAIN(m/year) | " \
+ ,"
" \
+ /)
+ table_footer = "
"
+ row_header = ""
+ row_footer = "
"
+
+ delete (lines)
+ lines = new(50000,string)
+ nline = 0
+
+ set_line(lines,nline,header)
+ set_line(lines,nline,table_header)
+;-----------------------------------------------
+;row of table
+
+ nrow = dimsizes(id933)
+ do n = 0,nrow-1
+ set_line(lines,nline,row_header)
+
+ txt1 = sprintf("%5.0f", id933(n))
+ txt2 = sprintf("%5.2f", y933(n))
+ txt3 = sprintf("%5.2f", x933(n))
+ txt4 = sprintf("%5.2f", npp933(n))
+ txt5 = sprintf("%5.2f", rain933(n))
+
+ set_line(lines,nline,""+txt1+" | ")
+ set_line(lines,nline,""+txt2+" | ")
+ set_line(lines,nline,""+txt3+" | ")
+ set_line(lines,nline,""+txt4+" | ")
+ set_line(lines,nline,""+txt5+" | ")
+
+ set_line(lines,nline,row_footer)
+ end do
+;-----------------------------------------------
+ set_line(lines,nline,table_footer)
+ set_line(lines,nline,footer)
+
+; Now write to an HTML file.
+ delete (idx)
+ idx = ind(.not.ismissing(lines))
+ if(.not.any(ismissing(idx))) then
+ asciiwrite(output_html,lines(idx))
+ else
+ print ("error?")
+ end if
+;*******************************************************************
+;(A)-3 html table of site81 -- model vs ob
+;*******************************************************************
+ output_html = "table_site81_model_vs_ob.html"
+
+ header_text = "Model "+model_name+" vs Observed at 81 sites
"
+
+ header = (/"" \
+ ,"" \
+ ,"CLAMP metrics" \
+ ,"" \
+ ,header_text \
+ /)
+ footer = ""
+
+ delete (table_header)
+ table_header_text = " "+model_name+" | "
+ table_header = (/ \
+ "" \
+ ,"" \
+ ," Site ID | " \
+ ," Latitude | " \
+ ," Longitude | " \
+ ," NPP(gC/m2.year) | " \
+ ," RAIN(m/year) | " \
+ ,"
" \
+ ,"" \
+ ," observed | " \
+ , table_header_text \
+ ," observed | " \
+ , table_header_text \
+ ,"
" \
+ /)
+ table_footer = "
"
+ row_header = ""
+ row_footer = "
"
+
+ delete (lines)
+ lines = new(50000,string)
+ nline = 0
+
+ set_line(lines,nline,header)
+ set_line(lines,nline,table_header)
+;-----------------------------------------------
+;row of table
+
+ nrow = dimsizes(id81)
+ do n = 0,nrow-1
+ set_line(lines,nline,row_header)
+
+ txt1 = sprintf("%5.0f", id81(n))
+ txt2 = sprintf("%5.2f", y81(n))
+ txt3 = sprintf("%5.2f", x81(n))
+ txt4 = sprintf("%5.2f", npp81(n))
+ txt5 = sprintf("%5.2f", nppmod81(n))
+ txt6 = sprintf("%5.2f", rain81(n))
+ txt7 = sprintf("%5.2f", rainmod81(n))
+
+ set_line(lines,nline,""+txt1+" | ")
+ set_line(lines,nline,""+txt2+" | ")
+ set_line(lines,nline,""+txt3+" | ")
+ set_line(lines,nline,""+txt4+" | ")
+ set_line(lines,nline,""+txt5+" | ")
+ set_line(lines,nline,""+txt6+" | ")
+ set_line(lines,nline,""+txt7+" | ")
+
+ set_line(lines,nline,row_footer)
+ end do
+;-----------------------------------------------
+ set_line(lines,nline,table_footer)
+ set_line(lines,nline,footer)
+
+; Now write to an HTML file.
+ delete (idx)
+ idx = ind(.not.ismissing(lines))
+ if(.not.any(ismissing(idx))) then
+ asciiwrite(output_html,lines(idx))
+ else
+ print ("error?")
+ end if
+
+;*******************************************************************
+;(A)-4 html table of site933 -- model vs ob
+;*******************************************************************
+ output_html = "table_site933_model_vs_ob.html"
+
+ header_text = "Model "+model_name+" vs Observed at 933 sites
"
+
+ header = (/"" \
+ ,"" \
+ ,"CLAMP metrics" \
+ ,"" \
+ ,header_text \
+ /)
+ footer = ""
+
+ delete (table_header)
+ table_header_text = " "+model_name+" | "
+ table_header = (/ \
+ "" \
+ ,"" \
+ ," Site ID | " \
+ ," Latitude | " \
+ ," Longitude | " \
+ ," NPP(gC/m2.year) | " \
+ ," RAIN(m/year) | " \
+ ,"
" \
+ ,"" \
+ ," observed | " \
+ , table_header_text \
+ ," observed | " \
+ , table_header_text \
+ ,"
" \
+ /)
+ table_footer = "
"
+ row_header = ""
+ row_footer = "
"
+
+ delete (lines)
+ lines = new(50000,string)
+ nline = 0
+
+ set_line(lines,nline,header)
+ set_line(lines,nline,table_header)
+;-----------------------------------------------
+;row of table
+
+ nrow = dimsizes(id933)
+ do n = 0,nrow-1
+ set_line(lines,nline,row_header)
+
+ txt1 = sprintf("%5.0f", id933(n))
+ txt2 = sprintf("%5.2f", y933(n))
+ txt3 = sprintf("%5.2f", x933(n))
+ txt4 = sprintf("%5.2f", npp933(n))
+ txt5 = sprintf("%5.2f", nppmod933(n))
+ txt6 = sprintf("%5.2f", rain933(n))
+ txt7 = sprintf("%5.2f", rainmod933(n))
+
+ set_line(lines,nline,""+txt1+" | ")
+ set_line(lines,nline,""+txt2+" | ")
+ set_line(lines,nline,""+txt3+" | ")
+ set_line(lines,nline,""+txt4+" | ")
+ set_line(lines,nline,""+txt5+" | ")
+ set_line(lines,nline,""+txt6+" | ")
+ set_line(lines,nline,""+txt7+" | ")
+
+ set_line(lines,nline,row_footer)
+ end do
+;-----------------------------------------------
+ set_line(lines,nline,table_footer)
+ set_line(lines,nline,footer)
+
+; Now write to an HTML file.
+ delete (idx)
+ idx = ind(.not.ismissing(lines))
+ if(.not.any(ismissing(idx))) then
+ asciiwrite(output_html,lines(idx))
+ else
+ print ("error?")
+ end if
+;***************************************************************************
+;(A)-5 scatter plot, model vs ob 81
+;***************************************************************************
+
+ plot_name = "scatter_model_vs_ob_81"
+ title = model_name +" vs ob 81 sites"
+
+ wks = gsn_open_wks (plot_type,plot_name) ; open workstation
+ res = True ; plot mods desired
+ res@tiMainString = title ; add title
+ res@xyMarkLineModes = "Markers" ; choose which have markers
+ res@xyMarkers = 16 ; choose type of marker
+ res@xyMarkerColor = "red" ; Marker color
+ res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
+ res@tmLabelAutoStride = True ; nice tick mark labels
+
+ res@gsnDraw = False
+ res@gsnFrame = False ; don't advance frame yet
+;-------------------------------
+;compute correlation coef. and M
+ ccr81 = esccr(nppmod81,npp81,0)
+;print (ccr81)
+
+;new eq
+ bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81)))
+ M81 = (1. - (bias/dimsizes(y81)))*5.
+
+ M_npp_S81 = sprintf("%.2f", M81)
+ system("sed s#M_npp_S81#"+M_npp_S81+"# "+html_name+" > "+html_new+";"+ \
+ "mv -f "+html_new+" "+html_name)
+ print (M_npp_S81)
+
+ tRes = True
+ tRes@txFontHeightF = 0.025
+
+ correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
+
+ gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
+;--------------------------------
+ plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot
+;-------------------------------
+; add polyline
+
+ dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)
+;-------------------------------
+ draw (plot)
+ frame(wks)
+ clear (wks)
+ delete (dum)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"-*."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+;***************************************************************************
+;(A)-6 scatter plot, model vs ob 933
+;***************************************************************************
+
+ plot_name = "scatter_model_vs_ob_933"
+ title = model_name +" vs ob 933 sites"
+
+ wks = gsn_open_wks (plot_type,plot_name) ; open workstation
+ res = True ; plot mods desired
+ res@tiMainString = title ; add title
+ res@xyMarkLineModes = "Markers" ; choose which have markers
+ res@xyMarkers = 16 ; choose type of marker
+ res@xyMarkerColor = "red" ; Marker color
+ res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
+ res@tmLabelAutoStride = True ; nice tick mark labels
+
+ res@gsnDraw = False
+ res@gsnFrame = False ; don't advance frame yet
+;-------------------------------
+;compute correlation coef. and M
+ ccr933 = esccr(nppmod933,npp933,0)
+;print (ccr933)
+
+;new eq
+ bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))
+ M933 = (1. - (bias/dimsizes(y933)))*5.
+
+ M_npp_S933 = sprintf("%.2f", M933)
+ system("sed s#M_npp_S933#"+M_npp_S933+"# "+html_name+" > "+html_new+";"+ \
+ "mv -f "+html_new+" "+html_name)
+ print (M_npp_S933)
+
+ tRes = True
+ tRes@txFontHeightF = 0.025
+
+ correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
+
+ gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
+;--------------------------------
+ plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot
+;-------------------------------
+; add polyline
+
+ dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)
+;-------------------------------
+ draw (plot)
+ frame(wks)
+ clear (wks)
+ delete (dum)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"-*."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+;**************************************************************************
+;(B) histogram 81
+;**************************************************************************
+;get data
+ RAIN1_1D = ndtooned(rain81)
+ RAIN2_1D = ndtooned(rainmod81)
+ NPP1_1D = ndtooned(npp81)
+ NPP2_1D = ndtooned(nppmod81)
+
+; number of bin
+ nx = 8
+
+ xvalues = new((/2,nx/),float)
+ yvalues = new((/2,nx/),float)
+ mn_yvalues = new((/2,nx/),float)
+ mx_yvalues = new((/2,nx/),float)
+ dx4 = new((/1/),float)
+
+ get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
+ ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
+
+;----------------------------------------
+;compute correlation coeff and M score
+
+ u = yvalues(0,:)
+ v = yvalues(1,:)
+
+ good = ind(.not.ismissing(u) .and. .not.ismissing(v))
+ uu = u(good)
+ vv = v(good)
+
+ ccr81h = esccr(uu,vv,0)
+;print (ccr81h)
+
+;new eq
+ bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
+ M81h = (1.- (bias/dimsizes(uu)))*5.
+
+ M_npp_H81 = sprintf("%.2f", M81h)
+ system("sed s#M_npp_H81#"+M_npp_H81+"# "+html_name+" > "+html_new+";"+ \
+ "mv -f "+html_new+" "+html_name)
+ print (M_npp_H81)
+
+ delete (u)
+ delete (v)
+ delete (uu)
+ delete (vv)
+;----------------------------------------------------------------------
+; histogram res
+ resh = True
+ resh@gsnMaximize = True
+ resh@gsnDraw = False
+ resh@gsnFrame = False
+ resh@xyMarkLineMode = "Markers"
+ resh@xyMarkerSizeF = 0.014
+ resh@xyMarker = 16
+ resh@xyMarkerColors = (/"Brown","Blue"/)
+ resh@trYMinF = min(mn_yvalues) - 10.
+ resh@trYMaxF = max(mx_yvalues) + 10.
+
+ resh@tiYAxisString = "NPP (g C/m2/year)"
+ resh@tiXAxisString = "Precipitation (m/year)"
+
+ max_bar = new((/2,nx/),graphic)
+ min_bar = new((/2,nx/),graphic)
+ max_cap = new((/2,nx/),graphic)
+ min_cap = new((/2,nx/),graphic)
+
+ lnres = True
+ line_colors = (/"brown","blue"/)
+
+;**************************************************************************
+;(B)-1 histogram plot, ob 81 site
+;**************************************************************************
+
+ plot_name = "histogram_ob_81"
+ title = "Observed 81 site"
+ resh@tiMainString = title
+
+ wks = gsn_open_wks (plot_type,plot_name)
+
+ xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
+
+;-------------------------------
+;Attach the vertical bar and the horizontal cap line
+
+ do nd=0,0
+ lnres@gsLineColor = line_colors(nd)
+ do i=0,nx-1
+
+ if(.not.ismissing(mn_yvalues(nd,i)).and. \
+ .not.ismissing(mx_yvalues(nd,i))) then
+;
+; Attach the vertical bar, both above and below the marker.
+;
+ x1 = xvalues(nd,i)
+ y1 = yvalues(nd,i)
+ y2 = mn_yvalues(nd,i)
+ plx = (/x1,x1/)
+ ply = (/y1,y2/)
+ min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
+
+ y2 = mx_yvalues(nd,i)
+ plx = (/x1,x1/)
+ ply = (/y1,y2/)
+ max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
+;
+; Attach the horizontal cap line, both above and below the marker.
+;
+ x1 = xvalues(nd,i) - dx4
+ x2 = xvalues(nd,i) + dx4
+ y1 = mn_yvalues(nd,i)
+ plx = (/x1,x2/)
+ ply = (/y1,y1/)
+ min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
+
+ y1 = mx_yvalues(nd,i)
+ plx = (/x1,x2/)
+ ply = (/y1,y1/)
+ max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
+ end if
+ end do
+ end do
+
+ draw(xy)
+ frame(wks)
+ clear (wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"-*."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+;****************************************************************************
+;(B)-2 histogram plot, model vs ob 81 site
+;****************************************************************************
+
+ plot_name = "histogram_model_vs_ob_81"
+ title = model_name+ " vs Observed 81 site"
+ resh@tiMainString = title
+
+ wks = gsn_open_wks (plot_type,plot_name) ; open workstation
+
+;-----------------------------
+; Add a boxed legend using the more simple method, which won't have
+; vertical lines going through the markers.
+
+ resh@pmLegendDisplayMode = "Always"
+; resh@pmLegendWidthF = 0.1
+ resh@pmLegendWidthF = 0.08
+ resh@pmLegendHeightF = 0.05
+ resh@pmLegendOrthogonalPosF = -1.17
+; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
+; resh@pmLegendParallelPosF = 0.18
+ resh@pmLegendParallelPosF = 0.88 ;(rightward)
+
+; resh@lgPerimOn = False
+ resh@lgLabelFontHeightF = 0.015
+ resh@xyExplicitLegendLabels = (/"observed",model_name/)
+;-----------------------------
+ tRes = True
+ tRes@txFontHeightF = 0.025
+
+ correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
+
+ gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
+
+ xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
+;-------------------------------
+;Attach the vertical bar and the horizontal cap line
+
+ do nd=0,1
+ lnres@gsLineColor = line_colors(nd)
+ do i=0,nx-1
+
+ if(.not.ismissing(mn_yvalues(nd,i)).and. \
+ .not.ismissing(mx_yvalues(nd,i))) then
+;
+; Attach the vertical bar, both above and below the marker.
+;
+ x1 = xvalues(nd,i)
+ y1 = yvalues(nd,i)
+ y2 = mn_yvalues(nd,i)
+ min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
+
+ y2 = mx_yvalues(nd,i)
+ max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
+;
+; Attach the horizontal cap line, both above and below the marker.
+;
+ x1 = xvalues(nd,i) - dx4
+ x2 = xvalues(nd,i) + dx4
+ y1 = mn_yvalues(nd,i)
+ min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
+
+ y1 = mx_yvalues(nd,i)
+ max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
+ end if
+ end do
+ end do
+
+ draw(xy)
+ frame(wks)
+ clear(wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"-*."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+
+ delete (RAIN1_1D)
+ delete (RAIN2_1D)
+ delete (NPP1_1D)
+ delete (NPP2_1D)
+;delete (range)
+ delete (xvalues)
+ delete (yvalues)
+ delete (mn_yvalues)
+ delete (mx_yvalues)
+ delete (good)
+ delete (max_bar)
+ delete (min_bar)
+ delete (max_cap)
+ delete (min_cap)
+
+;**************************************************************************
+;(B) histogram 933
+;**************************************************************************
+
+; get data
+ RAIN1_1D = ndtooned(rain933)
+ RAIN2_1D = ndtooned(rainmod933)
+ NPP1_1D = ndtooned(npp933)
+ NPP2_1D = ndtooned(nppmod933)
+
+; number of bin
+ nx = 10
+
+ xvalues = new((/2,nx/),float)
+ yvalues = new((/2,nx/),float)
+ mn_yvalues = new((/2,nx/),float)
+ mx_yvalues = new((/2,nx/),float)
+ dx4 = new((/1/),float)
+
+ get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
+ ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
+
+;----------------------------------------
+;compute correlation coeff and M score
+
+ u = yvalues(0,:)
+ v = yvalues(1,:)
+
+ good = ind(.not.ismissing(u) .and. .not.ismissing(v))
+ uu = u(good)
+ vv = v(good)
+
+ ccr933h = esccr(uu,vv,0)
+;print (ccr933h)
+
+;new eq
+ bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
+ M933h = (1.- (bias/dimsizes(uu)))*5.
+
+ M_npp_H933 = sprintf("%.2f", M933h)
+ system("sed s#M_npp_H933#"+M_npp_H933+"# "+html_name+" > "+html_new+";"+ \
+ "mv -f "+html_new+" "+html_name)
+ print (M_npp_H933)
+
+ delete (u)
+ delete (v)
+ delete (uu)
+ delete (vv)
+;----------------------------------------------------------------------
+; histogram res
+ delete (resh)
+ resh = True
+ resh@gsnMaximize = True
+ resh@gsnDraw = False
+ resh@gsnFrame = False
+ resh@xyMarkLineMode = "Markers"
+ resh@xyMarkerSizeF = 0.014
+ resh@xyMarker = 16
+ resh@xyMarkerColors = (/"Brown","Blue"/)
+ resh@trYMinF = min(mn_yvalues) - 10.
+ resh@trYMaxF = max(mx_yvalues) + 10.
+
+ resh@tiYAxisString = "NPP (g C/m2/year)"
+ resh@tiXAxisString = "Precipitation (m/year)"
+
+ max_bar = new((/2,nx/),graphic)
+ min_bar = new((/2,nx/),graphic)
+ max_cap = new((/2,nx/),graphic)
+ min_cap = new((/2,nx/),graphic)
+
+ lnres = True
+ line_colors = (/"brown","blue"/)
+
+;**************************************************************************
+;(B)-3 histogram plot, ob 933 site
+;**************************************************************************
+
+ plot_name = "histogram_ob_933"
+ title = "Observed 933 site"
+ resh@tiMainString = title
+
+ wks = gsn_open_wks (plot_type,plot_name)
+
+ xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
+
+;-------------------------------
+;Attach the vertical bar and the horizontal cap line
+
+ do nd=0,0
+ lnres@gsLineColor = line_colors(nd)
+ do i=0,nx-1
+
+ if(.not.ismissing(mn_yvalues(nd,i)).and. \
+ .not.ismissing(mx_yvalues(nd,i))) then
+
+; Attach the vertical bar, both above and below the marker.
+
+ x1 = xvalues(nd,i)
+ y1 = yvalues(nd,i)
+ y2 = mn_yvalues(nd,i)
+ min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
+
+ y2 = mx_yvalues(nd,i)
+ max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
+
+; Attach the horizontal cap line, both above and below the marker.
+
+ x1 = xvalues(nd,i) - dx4
+ x2 = xvalues(nd,i) + dx4
+ y1 = mn_yvalues(nd,i)
+ min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
+
+ y1 = mx_yvalues(nd,i)
+ max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
+ end if
+ end do
+ end do
+
+ draw(xy)
+ frame(wks)
+ delete (xy)
+ clear (wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"-*."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+;****************************************************************************
+;(B)-4 histogram plot, model vs ob 933 site
+;****************************************************************************
+
+ plot_name = "histogram_model_vs_ob_933"
+ title = model_name+ " vs Observed 933 site"
+ resh@tiMainString = title
+
+ wks = gsn_open_wks (plot_type,plot_name) ; open workstation
+
+;-----------------------------
+; Add a boxed legend using the more simple method, which won't have
+; vertical lines going through the markers.
+
+ resh@pmLegendDisplayMode = "Always"
+; resh@pmLegendWidthF = 0.1
+ resh@pmLegendWidthF = 0.08
+ resh@pmLegendHeightF = 0.05
+ resh@pmLegendOrthogonalPosF = -1.17
+; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
+; resh@pmLegendParallelPosF = 0.18
+ resh@pmLegendParallelPosF = 0.88 ;(rightward)
+
+; resh@lgPerimOn = False
+ resh@lgLabelFontHeightF = 0.015
+ resh@xyExplicitLegendLabels = (/"observed",model_name/)
+;-----------------------------
+ tRes = True
+ tRes@txFontHeightF = 0.025
+
+ correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
+
+ gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
+
+ xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
+;-------------------------------
+;Attach the vertical bar and the horizontal cap line
+
+ do nd=0,1
+ lnres@gsLineColor = line_colors(nd)
+ do i=0,nx-1
+
+ if(.not.ismissing(mn_yvalues(nd,i)).and. \
+ .not.ismissing(mx_yvalues(nd,i))) then
+;
+; Attach the vertical bar, both above and below the marker.
+;
+ x1 = xvalues(nd,i)
+ y1 = yvalues(nd,i)
+ y2 = mn_yvalues(nd,i)
+ min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
+
+ y2 = mx_yvalues(nd,i)
+ max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
+;
+; Attach the horizontal cap line, both above and below the marker.
+;
+ x1 = xvalues(nd,i) - dx4
+ x2 = xvalues(nd,i) + dx4
+ y1 = mn_yvalues(nd,i)
+ min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
+
+ y1 = mx_yvalues(nd,i)
+ max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
+ end if
+ end do
+ end do
+
+ draw(xy)
+ frame(wks)
+ delete(xy)
+ clear(wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"-*."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+;***************************************************************************
+;(C) global contour
+;***************************************************************************
+
+;res
+ resg = True ; Use plot options
+ resg@cnFillOn = True ; Turn on color fill
+ resg@gsnSpreadColors = True ; use full colormap
+; resg@cnFillMode = "RasterFill" ; Turn on raster color
+; resg@lbLabelAutoStride = True
+ resg@cnLinesOn = False ; Turn off contourn lines
+ resg@mpFillOn = False ; Turn off map fill
+
+ resg@gsnSpreadColors = True ; use full colormap
+ resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
+ resg@cnMinLevelValF = 0. ; Min level
+ resg@cnMaxLevelValF = 2200. ; Max level
+ resg@cnLevelSpacingF = 200. ; interval
+
+;****************************************************************************
+;(C)-1 global contour plot, ob
+;****************************************************************************
+
+ delta = 0.00000000001
+ nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
+
+ plot_name = "global_ob"
+ title = "Observed MODIS MOD 17"
+ resg@tiMainString = title
+
+ wks = gsn_open_wks (plot_type,plot_name) ; open workstation
+ gsn_define_colormap(wks,"gui_default") ; choose colormap
+
+ plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)
+
+ frame(wks)
+ clear (wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"-*."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+;****************************************************************************
+;(C)-2 global contour plot, model
+;****************************************************************************
+
+ plot_name = "global_model"
+ title = "Model "+ model_name
+ resg@tiMainString = title
+
+ wks = gsn_open_wks (plot_type,plot_name) ; open workstation
+ gsn_define_colormap(wks,"gui_default") ; choose colormap
+
+ plot = gsn_csm_contour_map_ce(wks,nppmod,resg)
+
+ frame(wks)
+ clear (wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"-*."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+;****************************************************************************
+;(C)-3 global contour plot, model vs ob
+;****************************************************************************
+
+ plot_name = "global_model_vs_ob"
+
+ wks = gsn_open_wks (plot_type,plot_name) ; open workstation
+ gsn_define_colormap(wks,"gui_default") ; choose colormap
+
+ delete (plot)
+ plot=new(3,graphic) ; create graphic array
+
+ resg@gsnFrame = False ; Do not draw plot
+ resg@gsnDraw = False ; Do not advance frame
+
+; compute correlation coef and M score
+
+ uu1 = ndtooned(nppmod)
+ vv1 = ndtooned(nppglobe)
+
+ delete (good)
+ good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))
+
+ ug = uu1(good)
+ vg = vv1(good)
+
+ ccrG = esccr(ug,vg,0)
+; print (ccrG)
+
+ MG = (ccrG*ccrG)* 5.0
+
+ M_npp_G = sprintf("%.2f", MG)
+ system("sed s#M_npp_G#"+M_npp_G+"# "+html_name+" > "+html_new+";"+ \
+ "mv -f "+html_new+" "+html_name)
+ print (M_npp_G)
+
+; plot correlation coef
+
+ gRes = True
+ gRes@txFontHeightF = 0.02
+ gRes@txAngleF = 90
+
+ correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
+
+ gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
+;--------------------------------------------------------------------
+;(a) ob
+
+ title = "Observed MODIS MOD 17"
+ resg@tiMainString = title
+
+ plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)
+
+;(b) model
+
+ title = "Model "+ model_name
+ resg@tiMainString = title
+
+ plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)
+
+;(c) model-ob
+
+ zz = nppmod
+ zz = nppmod - nppglobe
+ title = "Model_"+model_name+" - Observed"
+
+ resg@cnMinLevelValF = -500 ; Min level
+ resg@cnMaxLevelValF = 500. ; Max level
+ resg@cnLevelSpacingF = 50. ; interval
+ resg@tiMainString = title
+
+ plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
+
+ pres = True ; panel plot mods desired
+; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
+ ; indiv. plots in panel
+ pres@gsnMaximize = True ; fill the page
+
+ gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
+
+ frame (wks)
+ clear (wks)
+ delete (plot)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"-*."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+;***************************************************************************
+;(D)-1 zonal line plot, ob
+;***************************************************************************
+
+ vv = zonalAve(nppglobe)
+ vv@long_name = nppglobe@long_name
+
+ plot_name = "zonal_ob"
+ title = "Observed MODIS MOD 17"
+
+ resz = True
+ resz@tiMainString = title
+
+ wks = gsn_open_wks (plot_type,plot_name) ; open workstation
+
+ plot = gsn_csm_xy (wks,ym,vv,resz)
+
+ frame(wks)
+ clear (wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"-*."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+;****************************************************************************
+;(D)-2 zonal line plot, model vs ob
+;****************************************************************************
+
+ s = new ((/2,dimsizes(ym)/), float)
+
+ s(0,:) = zonalAve(nppglobe)
+ s(1,:) = zonalAve(nppmod)
+
+ s@long_name = nppglobe@long_name
+;-------------------------------------------
+; compute correlation coef and M score
+
+ ccrZ = esccr(s(1,:), s(0,:),0)
+; print (ccrZ)
+
+ MZ = (ccrZ*ccrZ)* 5.0
+
+ M_npp_Z = sprintf("%.2f", MZ)
+ system("sed s#M_npp_Z#"+M_npp_Z+"# "+html_name+" > "+html_new+";"+ \
+ "mv -f "+html_new+" "+html_name)
+ print (M_npp_Z)
+;-------------------------------------------
+ plot_name = "zonal_model_vs_ob"
+ title = "Zonal Average"
+ resz@tiMainString = title
+
+ wks = gsn_open_wks (plot_type,plot_name)
+
+; resz@vpHeightF = 0.4 ; change aspect ratio of plot
+; resz@vpWidthF = 0.7
+
+ resz@xyMonoLineColor = "False" ; want colored lines
+ resz@xyLineColors = (/"black","red"/) ; colors chosen
+; resz@xyLineThicknesses = (/3.,3./) ; line thicknesses
+ resz@xyLineThicknesses = (/2.,2./) ; line thicknesses
+ resz@xyDashPatterns = (/0.,0./) ; make all lines solid
+
+ resz@tiYAxisString = s@long_name ; add a axis title
+ resz@txFontHeightF = 0.0195 ; change title font heights
+
+; Legent
+ resz@pmLegendDisplayMode = "Always" ; turn on legend
+ resz@pmLegendSide = "Top" ; Change location of
+; resz@pmLegendParallelPosF = .45 ; move units right
+ resz@pmLegendParallelPosF = .82 ; move units right
+ resz@pmLegendOrthogonalPosF = -0.4 ; move units down
+
+ resz@pmLegendWidthF = 0.10 ; Change width and
+ resz@pmLegendHeightF = 0.10 ; height of legend.
+ resz@lgLabelFontHeightF = .02 ; change font height
+; resz@lgTitleOn = True ; turn on legend title
+; resz@lgTitleString = "Example" ; create legend title
+; resz@lgTitleFontHeightF = .025 ; font of legend title
+ resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels
+;--------------------------------------------------------------------
+ zRes = True
+ zRes@txFontHeightF = 0.025
+
+ correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"
+
+ gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)
+;--------------------------------------------------------------------
+
+ plot = gsn_csm_xy (wks,ym,s,resz)
+
+ frame(wks)
+ clear (wks)
+
+ system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
+ "mv "+plot_name+"-0."+plot_type_new+" "+plot_name+"."+plot_type_new+";"+ \
+ "rm "+plot_name+"-*."+plot_type_new+";"+ \
+ "rm "+plot_name+"."+plot_type)
+;***************************************************************************
+; tar up output plots
+;***************************************************************************
+
+ temp_name = "npp." + model_name
+ system("mkdir -p " + temp_name+";"+ \
+ "cp "+ html_name + " " +temp_name+"/table.html"+";"+ \
+ "mv table_*.html " + temp_name +";"+ \
+ "mv *.png " + temp_name +";"+ \
+ "tar cf "+ temp_name +".tar " + temp_name)
+;***************************************************************************
+end
+