1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/all/01.npp.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,1370 @@
1.4 +;*****************************************************
1.5 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.6 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.7 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.8 +;****************************************************************************
1.9 +procedure set_line(lines:string,nline:integer,newlines:string)
1.10 +begin
1.11 +; add line to ascci/html file
1.12 +
1.13 + nnewlines = dimsizes(newlines)
1.14 + if(nline+nnewlines-1.ge.dimsizes(lines))
1.15 + print("set_line: bad index, not setting anything.")
1.16 + return
1.17 + end if
1.18 + lines(nline:nline+nnewlines-1) = newlines
1.19 +; print ("lines = " + lines(nline:nline+nnewlines-1))
1.20 + nline = nline + nnewlines
1.21 + return
1.22 +end
1.23 +;****************************************************************************
1.24 +procedure get_bin(RAIN1_1D[*]:numeric, NPP1_1D[*]:numeric \
1.25 + ,RAIN2_1D[*]:numeric, NPP2_1D[*]:numeric \
1.26 + ,xvalues[*][*]:numeric, yvalues[*][*]:numeric \
1.27 + ,mn_yvalues[*][*]:numeric, mx_yvalues[*][*]:numeric \
1.28 + ,dx4[1]:numeric \
1.29 + )
1.30 +begin
1.31 +; Calculaee "nice" bins for binning the data in equally spaced ranges.
1.32 +; input : RAIN1_1D, RAIN2_1D, NPP1_1D, NPP2_1D
1.33 +; output: xvalues, yvalues, mn_yvalues, mx_yvalues,dx4
1.34 +
1.35 + nbins = 15 ; Number of bins to use.
1.36 +
1.37 + nicevals = nice_mnmxintvl(min(RAIN1_1D),max(RAIN1_1D),nbins,True)
1.38 + nvals = floattoint((nicevals(1) - nicevals(0))/nicevals(2) + 1)
1.39 + range = fspan(nicevals(0),nicevals(1),nvals)
1.40 +
1.41 +; Use this range information to grab all the values in a
1.42 +; particular range, and then take an average.
1.43 +
1.44 + nr = dimsizes(range)
1.45 + nx = nr-1
1.46 +; print (nx)
1.47 +
1.48 +; xvalues = new((/2,nx/),typeof(RAIN1_1D))
1.49 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.50 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.51 + dx4 = dx/4 ; 1/4 of the range
1.52 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.53 +; yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.54 +; mn_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.55 +; mx_yvalues = new((/2,nx/),typeof(RAIN1_1D))
1.56 +
1.57 + do nd=0,1
1.58 +
1.59 +; See if we are doing model or observational data.
1.60 +
1.61 + if(nd.eq.0) then
1.62 + data = RAIN1_1D
1.63 + npp_data = NPP1_1D
1.64 + else
1.65 + data = RAIN2_1D
1.66 + npp_data = NPP2_1D
1.67 + end if
1.68 +
1.69 +; Loop through each range and check for values.
1.70 +
1.71 + do i=0,nr-2
1.72 + if (i.ne.(nr-2)) then
1.73 +; print("")
1.74 +; print("In range ["+range(i)+","+range(i+1)+")")
1.75 + idx = ind((range(i).le.data).and.(data.lt.range(i+1)))
1.76 + else
1.77 +; print("")
1.78 +; print("In range ["+range(i)+",)")
1.79 + idx = ind(range(i).le.data)
1.80 + end if
1.81 +
1.82 +; Calculate average, and get min and max.
1.83 +
1.84 + if(.not.any(ismissing(idx))) then
1.85 + yvalues(nd,i) = avg(npp_data(idx))
1.86 + mn_yvalues(nd,i) = min(npp_data(idx))
1.87 + mx_yvalues(nd,i) = max(npp_data(idx))
1.88 + count = dimsizes(idx)
1.89 + else
1.90 + count = 0
1.91 + yvalues(nd,i) = yvalues@_FillValue
1.92 + mn_yvalues(nd,i) = yvalues@_FillValue
1.93 + mx_yvalues(nd,i) = yvalues@_FillValue
1.94 + end if
1.95 +
1.96 +; Print out information.
1.97 +; print(data_types(nd) + ": " + count + " points, avg = " + yvalues(nd,i))
1.98 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.99 +
1.100 +; Clean up for next time in loop.
1.101 +
1.102 + delete(idx)
1.103 + end do
1.104 + delete(data)
1.105 + delete(npp_data)
1.106 + end do
1.107 +end
1.108 +;****************************************************************************
1.109 +; Main code.
1.110 +;****************************************************************************
1.111 +begin
1.112 +
1.113 + plot_type = "ps"
1.114 + plot_type_new = "png"
1.115 +
1.116 +;-----------------------------------------------------
1.117 +; edit table.html of current model for movel1_vs_model2
1.118 +
1.119 + if (isvar("compare")) then
1.120 + html_name2 = compare+"/table.html"
1.121 + html_new2 = html_name2 +".new"
1.122 + end if
1.123 +;------------------------------------------------------
1.124 +; edit table.html for current model
1.125 +
1.126 + html_name = model_name+"/table.html"
1.127 + html_new = html_name +".new"
1.128 +
1.129 +;------------------------------------------------------
1.130 +; read model data
1.131 +
1.132 + fm = addfile (dirm+film1,"r")
1.133 +
1.134 + nppmod0 = fm->NPP
1.135 + rainmod0 = fm->RAIN
1.136 + xm = fm->lon
1.137 + ym = fm->lat
1.138 +
1.139 + delete (fm)
1.140 +
1.141 +;----------------------------------------------------
1.142 +; read ob data
1.143 +
1.144 +;(1) data at 81 sites
1.145 +
1.146 + dir81 = diro + "npp/"
1.147 + fil81 = "data.81.nc"
1.148 + f81 = addfile (dir81+fil81,"r")
1.149 +
1.150 + id81 = f81->SITE_ID
1.151 + npp81 = f81->TNPP_C
1.152 + rain81 = tofloat(f81->PREC_ANN)
1.153 + x81 = f81->LONG_DD
1.154 + y81 = f81->LAT_DD
1.155 +
1.156 + delete (f81)
1.157 +
1.158 + id81@long_name = "SITE_ID"
1.159 +
1.160 +;change longitude from (-180,180) to (0,360)
1.161 +;for model data interpolation
1.162 +
1.163 + do i= 0,dimsizes(x81)-1
1.164 + if (x81(i) .lt. 0.) then
1.165 + x81(i) = x81(i)+ 360.
1.166 + end if
1.167 + end do
1.168 +
1.169 +;-------------------------------------
1.170 +;(2) data at 933 sites
1.171 +
1.172 + dir933 = diro + "npp/"
1.173 + fil933 = "data.933.nc"
1.174 + f933 = addfile (dir933+fil933,"r")
1.175 +
1.176 + id933 = f933->SITE_ID
1.177 + npp933 = f933->TNPP_C
1.178 + rain933 = f933->PREC
1.179 + x933 = f933->LONG_DD
1.180 + y933 = f933->LAT_DD
1.181 +
1.182 + delete (f933)
1.183 +
1.184 + id933@long_name = "SITE_ID"
1.185 +
1.186 +;change longitude from (-180,180) to (0,360)
1.187 +;for model data interpolation
1.188 +
1.189 + do i= 0,dimsizes(x933)-1
1.190 + if (x933(i) .lt. 0.) then
1.191 + x933(i) = x933(i)+ 360.
1.192 + end if
1.193 + end do
1.194 +
1.195 +;----------------------------------------
1.196 +;(3) global data, interpolated into model grid
1.197 +
1.198 + dirg = diro + "npp/"
1.199 + filg = "Npp_"+model_grid+"_mean.nc"
1.200 + fglobe = addfile (dirg+filg,"r")
1.201 +
1.202 + nppglobe0 = fglobe->NPP
1.203 + nppglobe = nppglobe0
1.204 +
1.205 + delete (fglobe)
1.206 +
1.207 +;***********************************************************************
1.208 +; interpolate model data into ob sites
1.209 +;***********************************************************************
1.210 +
1.211 + nppmod = nppmod0(0,:,:)
1.212 + rainmod = rainmod0(0,:,:)
1.213 + delete (nppmod0)
1.214 + delete (rainmod0)
1.215 +
1.216 + nppmod81 =linint2_points(xm,ym,nppmod,True,x81,y81,0)
1.217 +
1.218 + nppmod933 =linint2_points(xm,ym,nppmod,True,x933,y933,0)
1.219 +
1.220 + rainmod81 =linint2_points(xm,ym,rainmod,True,x81,y81,0)
1.221 +
1.222 + rainmod933=linint2_points(xm,ym,rainmod,True,x933,y933,0)
1.223 +
1.224 +;**********************************************************
1.225 +; unified units
1.226 +;**********************************************************
1.227 +; Units for these variables are:
1.228 +;
1.229 +; rain81 : mm/year
1.230 +; rainmod : mm/s
1.231 +; npp81 : g C/m^2/year
1.232 +; nppmod81: g C/m^2/s
1.233 +; nppglobe: g C/m^2/year
1.234 +;
1.235 +; We want to convert these to "m/year" and "g C/m^2/year".
1.236 +
1.237 + nsec_per_year = 60*60*24*365
1.238 +
1.239 + rain81 = rain81 / 1000.
1.240 + rainmod81 = (rainmod81/ 1000.) * nsec_per_year
1.241 + nppmod81 = nppmod81 * nsec_per_year
1.242 +
1.243 + rain933 = rain933 / 1000.
1.244 + rainmod933 = (rainmod933/ 1000.) * nsec_per_year
1.245 + nppmod933 = nppmod933 * nsec_per_year
1.246 +
1.247 + nppmod = nppmod * nsec_per_year
1.248 +
1.249 + npp81@units = "gC/m^2/yr"
1.250 + nppmod81@units = "gC/m^2/yr"
1.251 + npp933@units = "gC/m^2/yr"
1.252 + nppmod933@units = "gC/m^2/yr"
1.253 + nppmod@units = "gC/m^2/yr"
1.254 + nppglobe@units = "gC/m^2/yr"
1.255 + rain81@units = "m/yr"
1.256 + rainmod81@units = "m/yr"
1.257 + rain933@units = "m/yr"
1.258 + rainmod933@units = "m/yr"
1.259 +
1.260 + npp81@long_name = "Obs. NPP (gC/m2/year)"
1.261 + npp933@long_name = "Obs. NPP (gC/m2/year)"
1.262 + nppmod81@long_name = "Model NPP (gC/m2/year)"
1.263 + nppmod933@long_name = "Model NPP (gC/m2/year)"
1.264 + nppmod@long_name = "Model NPP (gC/m2/year)"
1.265 + nppglobe@long_name = "NPP (gC/m2/year)"
1.266 + rain81@long_name = "PREC (m/year)"
1.267 + rain933@long_name = "PREC (m/year)"
1.268 + rainmod81@long_name = "PREC (m/year)"
1.269 + rainmod933@long_name = "PREC (m/year)"
1.270 +
1.271 +; change longitude from 0-360 to -180-180
1.272 + x81 = where(x81 .gt. 180., x81 -360., x81 )
1.273 + x933 = where(x933 .gt. 180., x933-360., x933)
1.274 +
1.275 +;*******************************************************************
1.276 +;(A)-1 html table of site81 -- observed
1.277 +;*******************************************************************
1.278 +
1.279 + output_html = "table_site81_ob.html"
1.280 +
1.281 + header = (/"<HTML>" \
1.282 + ,"<HEAD>" \
1.283 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.284 + ,"</HEAD>" \
1.285 + ,"<H1>EMDI Observations Class A (81 sites)</H1>" \
1.286 + /)
1.287 + footer = "</HTML>"
1.288 +
1.289 + table_header = (/ \
1.290 + "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
1.291 + ,"<tr>" \
1.292 + ," <th bgcolor=DDDDDD >Site ID</th>" \
1.293 + ," <th bgcolor=DDDDDD >Latitude</th>" \
1.294 + ," <th bgcolor=DDDDDD >Longitude</th>" \
1.295 + ," <th bgcolor=DDDDDD >NPP(gC/m2/year)</th>" \
1.296 + ," <th bgcolor=DDDDDD >PPT(m/year)</th>" \
1.297 + ,"</tr>" \
1.298 + /)
1.299 + table_footer = "</table>"
1.300 + row_header = "<tr>"
1.301 + row_footer = "</tr>"
1.302 +
1.303 + lines = new(50000,string)
1.304 + nline = 0
1.305 +
1.306 + set_line(lines,nline,header)
1.307 + set_line(lines,nline,table_header)
1.308 +;-----------------------------------------------
1.309 +;row of table
1.310 +
1.311 + nrow = dimsizes(id81)
1.312 + do n = 0,nrow-1
1.313 + set_line(lines,nline,row_header)
1.314 +
1.315 + txt1 = sprintf("%5.0f", id81(n))
1.316 + txt2 = sprintf("%5.2f", y81(n))
1.317 + txt3 = sprintf("%5.2f", x81(n))
1.318 + txt4 = sprintf("%5.2f", npp81(n))
1.319 + txt5 = sprintf("%5.2f", rain81(n))
1.320 +
1.321 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.322 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.323 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.324 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.325 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.326 +
1.327 + set_line(lines,nline,row_footer)
1.328 + end do
1.329 +;-----------------------------------------------
1.330 + set_line(lines,nline,table_footer)
1.331 + set_line(lines,nline,footer)
1.332 +
1.333 +; Now write to an HTML file.
1.334 + idx = ind(.not.ismissing(lines))
1.335 + if(.not.any(ismissing(idx))) then
1.336 + asciiwrite(output_html,lines(idx))
1.337 + else
1.338 + print ("error?")
1.339 + end if
1.340 +
1.341 +;*******************************************************************
1.342 +;(A)-2 html table of site933 -- observed
1.343 +;*******************************************************************
1.344 + output_html = "table_site933_ob.html"
1.345 +
1.346 + header = (/"<HTML>" \
1.347 + ,"<HEAD>" \
1.348 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.349 + ,"</HEAD>" \
1.350 + ,"<H1>EMDI Observations Class B (933 sites)</H1>" \
1.351 + /)
1.352 + footer = "</HTML>"
1.353 +
1.354 + table_header = (/ \
1.355 + "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
1.356 + ,"<tr>" \
1.357 + ," <th bgcolor=DDDDDD >Site ID</th>" \
1.358 + ," <th bgcolor=DDDDDD >Latitude</th>" \
1.359 + ," <th bgcolor=DDDDDD >Longitude</th>" \
1.360 + ," <th bgcolor=DDDDDD >NPP(gC/m2/year)</th>" \
1.361 + ," <th bgcolor=DDDDDD >PPT(m/year)</th>" \
1.362 + ,"</tr>" \
1.363 + /)
1.364 + table_footer = "</table>"
1.365 + row_header = "<tr>"
1.366 + row_footer = "</tr>"
1.367 +
1.368 + delete (lines)
1.369 + lines = new(50000,string)
1.370 + nline = 0
1.371 +
1.372 + set_line(lines,nline,header)
1.373 + set_line(lines,nline,table_header)
1.374 +;-----------------------------------------------
1.375 +;row of table
1.376 +
1.377 + nrow = dimsizes(id933)
1.378 + do n = 0,nrow-1
1.379 + set_line(lines,nline,row_header)
1.380 +
1.381 + txt1 = sprintf("%5.0f", id933(n))
1.382 + txt2 = sprintf("%5.2f", y933(n))
1.383 + txt3 = sprintf("%5.2f", x933(n))
1.384 + txt4 = sprintf("%5.2f", npp933(n))
1.385 + txt5 = sprintf("%5.2f", rain933(n))
1.386 +
1.387 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.388 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.389 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.390 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.391 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.392 +
1.393 + set_line(lines,nline,row_footer)
1.394 + end do
1.395 +;-----------------------------------------------
1.396 + set_line(lines,nline,table_footer)
1.397 + set_line(lines,nline,footer)
1.398 +
1.399 +; Now write to an HTML file.
1.400 + delete (idx)
1.401 + idx = ind(.not.ismissing(lines))
1.402 + if(.not.any(ismissing(idx))) then
1.403 + asciiwrite(output_html,lines(idx))
1.404 + else
1.405 + print ("error?")
1.406 + end if
1.407 +
1.408 +;*******************************************************************
1.409 +;(A)-3 html table of site81 -- model vs ob
1.410 +;*******************************************************************
1.411 + output_html = "table_site81_model_vs_ob.html"
1.412 +
1.413 + header_text = "<H1>Model "+model_name+" vs Class A Observations (81 sites)</H1>"
1.414 +
1.415 + header = (/"<HTML>" \
1.416 + ,"<HEAD>" \
1.417 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.418 + ,"</HEAD>" \
1.419 + ,header_text \
1.420 + /)
1.421 + footer = "</HTML>"
1.422 +
1.423 + delete (table_header)
1.424 + table_header_text = " <th bgcolor=DDDDDD >"+model_name+"</th>"
1.425 + table_header = (/ \
1.426 + "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
1.427 + ,"<tr>" \
1.428 + ," <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \
1.429 + ," <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \
1.430 + ," <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \
1.431 + ," <th bgcolor=DDDDDD colspan=2>NPP(gC/m2.year)</th>" \
1.432 + ," <th bgcolor=DDDDDD colspan=2>RAIN(m/year)</th>" \
1.433 + ,"</tr>" \
1.434 + ,"<tr>" \
1.435 + ," <th bgcolor=DDDDDD >observed</th>" \
1.436 + , table_header_text \
1.437 + ," <th bgcolor=DDDDDD >observed</th>" \
1.438 + , table_header_text \
1.439 + ,"</tr>" \
1.440 + /)
1.441 + table_footer = "</table>"
1.442 + row_header = "<tr>"
1.443 + row_footer = "</tr>"
1.444 +
1.445 + delete (lines)
1.446 + lines = new(50000,string)
1.447 + nline = 0
1.448 +
1.449 + set_line(lines,nline,header)
1.450 + set_line(lines,nline,table_header)
1.451 +;-----------------------------------------------
1.452 +;row of table
1.453 +
1.454 + nrow = dimsizes(id81)
1.455 + do n = 0,nrow-1
1.456 + set_line(lines,nline,row_header)
1.457 +
1.458 + txt1 = sprintf("%5.0f", id81(n))
1.459 + txt2 = sprintf("%5.2f", y81(n))
1.460 + txt3 = sprintf("%5.2f", x81(n))
1.461 + txt4 = sprintf("%5.2f", npp81(n))
1.462 + txt5 = sprintf("%5.2f", nppmod81(n))
1.463 + txt6 = sprintf("%5.2f", rain81(n))
1.464 + txt7 = sprintf("%5.2f", rainmod81(n))
1.465 +
1.466 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.467 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.468 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.469 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.470 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.471 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.472 + set_line(lines,nline,"<th>"+txt7+"</th>")
1.473 +
1.474 + set_line(lines,nline,row_footer)
1.475 + end do
1.476 +;-----------------------------------------------
1.477 + set_line(lines,nline,table_footer)
1.478 + set_line(lines,nline,footer)
1.479 +
1.480 +; Now write to an HTML file.
1.481 + delete (idx)
1.482 + idx = ind(.not.ismissing(lines))
1.483 + if(.not.any(ismissing(idx))) then
1.484 + asciiwrite(output_html,lines(idx))
1.485 + else
1.486 + print ("error?")
1.487 + end if
1.488 +
1.489 +;*******************************************************************
1.490 +;(A)-4 html table of site933 -- model vs ob
1.491 +;*******************************************************************
1.492 + output_html = "table_site933_model_vs_ob.html"
1.493 +
1.494 + header_text = "<H1>Model "+model_name+" vs Class B Observations (933 sites)</H1>"
1.495 +
1.496 + header = (/"<HTML>" \
1.497 + ,"<HEAD>" \
1.498 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.499 + ,"</HEAD>" \
1.500 + ,header_text \
1.501 + /)
1.502 + footer = "</HTML>"
1.503 +
1.504 + delete (table_header)
1.505 + table_header_text = " <th bgcolor=DDDDDD >"+model_name+"</th>"
1.506 + table_header = (/ \
1.507 + "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
1.508 + ,"<tr>" \
1.509 + ," <th bgcolor=DDDDDD rowspan=2>Site ID</th>" \
1.510 + ," <th bgcolor=DDDDDD rowspan=2>Latitude</th>" \
1.511 + ," <th bgcolor=DDDDDD rowspan=2>Longitude</th>" \
1.512 + ," <th bgcolor=DDDDDD colspan=2>NPP(gC/m2.year)</th>" \
1.513 + ," <th bgcolor=DDDDDD colspan=2>RAIN(m/year)</th>" \
1.514 + ,"</tr>" \
1.515 + ,"<tr>" \
1.516 + ," <th bgcolor=DDDDDD >observed</th>" \
1.517 + , table_header_text \
1.518 + ," <th bgcolor=DDDDDD >observed</th>" \
1.519 + , table_header_text \
1.520 + ,"</tr>" \
1.521 + /)
1.522 + table_footer = "</table>"
1.523 + row_header = "<tr>"
1.524 + row_footer = "</tr>"
1.525 +
1.526 + delete (lines)
1.527 + lines = new(50000,string)
1.528 + nline = 0
1.529 +
1.530 + set_line(lines,nline,header)
1.531 + set_line(lines,nline,table_header)
1.532 +;-----------------------------------------------
1.533 +;row of table
1.534 +
1.535 + nrow = dimsizes(id933)
1.536 + do n = 0,nrow-1
1.537 + set_line(lines,nline,row_header)
1.538 +
1.539 + txt1 = sprintf("%5.0f", id933(n))
1.540 + txt2 = sprintf("%5.2f", y933(n))
1.541 + txt3 = sprintf("%5.2f", x933(n))
1.542 + txt4 = sprintf("%5.2f", npp933(n))
1.543 + txt5 = sprintf("%5.2f", nppmod933(n))
1.544 + txt6 = sprintf("%5.2f", rain933(n))
1.545 + txt7 = sprintf("%5.2f", rainmod933(n))
1.546 +
1.547 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.548 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.549 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.550 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.551 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.552 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.553 + set_line(lines,nline,"<th>"+txt7+"</th>")
1.554 +
1.555 + set_line(lines,nline,row_footer)
1.556 + end do
1.557 +;-----------------------------------------------
1.558 + set_line(lines,nline,table_footer)
1.559 + set_line(lines,nline,footer)
1.560 +
1.561 +; Now write to an HTML file.
1.562 + delete (idx)
1.563 + idx = ind(.not.ismissing(lines))
1.564 + if(.not.any(ismissing(idx))) then
1.565 + asciiwrite(output_html,lines(idx))
1.566 + else
1.567 + print ("error?")
1.568 + end if
1.569 +
1.570 +;***************************************************************************
1.571 +;(A)-5 scatter plot, model vs ob 81
1.572 +;***************************************************************************
1.573 +
1.574 + plot_name = "scatter_model_vs_ob_81"
1.575 + title = model_name +" vs Class A observations (81 sites)"
1.576 +
1.577 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.578 + res = True ; plot mods desired
1.579 + res@tiMainString = title ; add title
1.580 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.581 + res@xyMarkers = 16 ; choose type of marker
1.582 + res@xyMarkerColor = "red" ; Marker color
1.583 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.584 + res@tmLabelAutoStride = True ; nice tick mark labels
1.585 +
1.586 + res@gsnDraw = False
1.587 + res@gsnFrame = False ; don't advance frame yet
1.588 +;-------------------------------
1.589 +;compute correlation coef. and M
1.590 + ccr81 = esccr(nppmod81,npp81,0)
1.591 +;print (ccr81)
1.592 +
1.593 + score_max = 2.5
1.594 +
1.595 + bias = sum(abs(nppmod81-npp81)/(abs(nppmod81)+abs(npp81)))
1.596 + M81s = (1. - (bias/dimsizes(y81)))*score_max
1.597 + M_npp_S81 = sprintf("%.2f", M81s)
1.598 +
1.599 + if (isvar("compare")) then
1.600 + system("sed -e '1,/M_npp_S81/s/M_npp_S81/"+M_npp_S81+"/' "+html_name2+" > "+html_new2+";"+ \
1.601 + "mv -f "+html_new2+" "+html_name2)
1.602 + end if
1.603 +
1.604 + system("sed s#M_npp_S81#"+M_npp_S81+"# "+html_name+" > "+html_new+";"+ \
1.605 + "mv -f "+html_new+" "+html_name)
1.606 +
1.607 + tRes = True
1.608 + tRes@txFontHeightF = 0.025
1.609 +
1.610 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81)+")"
1.611 +
1.612 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.613 +;--------------------------------
1.614 + plot = gsn_csm_xy (wks,npp81,nppmod81,res) ; create plot
1.615 +;-------------------------------
1.616 +; add polyline
1.617 +
1.618 + dum=gsn_add_polyline(wks,plot,(/0,1200/),(/0,1200/),False)
1.619 +;-------------------------------
1.620 + draw (plot)
1.621 + delete (wks)
1.622 + delete (dum)
1.623 +
1.624 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.625 + "rm "+plot_name+"."+plot_type)
1.626 +
1.627 +;***************************************************************************
1.628 +;(A)-6 scatter plot, model vs ob 933
1.629 +;***************************************************************************
1.630 +
1.631 + plot_name = "scatter_model_vs_ob_933"
1.632 + title = model_name +" vs Class B Observations (933 sites)"
1.633 +
1.634 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.635 + res = True ; plot mods desired
1.636 + res@tiMainString = title ; add title
1.637 + res@xyMarkLineModes = "Markers" ; choose which have markers
1.638 + res@xyMarkers = 16 ; choose type of marker
1.639 + res@xyMarkerColor = "red" ; Marker color
1.640 + res@xyMarkerSizeF = 0.01 ; Marker size (default 0.01)
1.641 + res@tmLabelAutoStride = True ; nice tick mark labels
1.642 +
1.643 + res@gsnDraw = False
1.644 + res@gsnFrame = False ; don't advance frame yet
1.645 +;-------------------------------
1.646 +;compute correlation coef. and M
1.647 +
1.648 + ccr933 = esccr(nppmod933,npp933,0)
1.649 +
1.650 + score_max = 2.5
1.651 +
1.652 + bias = sum(abs(nppmod933-npp933)/(abs(nppmod933)+abs(npp933)))
1.653 + M933s = (1. - (bias/dimsizes(y933)))*score_max
1.654 + M_npp_S933 = sprintf("%.2f", M933s)
1.655 +
1.656 + if (isvar("compare")) then
1.657 + system("sed -e '1,/M_npp_S933/s/M_npp_S933/"+M_npp_S933+"/' "+html_name2+" > "+html_new2+";"+ \
1.658 + "mv -f "+html_new2+" "+html_name2)
1.659 + end if
1.660 +
1.661 + system("sed s#M_npp_S933#"+M_npp_S933+"# "+html_name+" > "+html_new+";"+ \
1.662 + "mv -f "+html_new+" "+html_name)
1.663 +
1.664 + tRes = True
1.665 + tRes@txFontHeightF = 0.025
1.666 +
1.667 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933)+")"
1.668 +
1.669 + gsn_text_ndc(wks,correlation_text,0.5,0.95,tRes)
1.670 +;--------------------------------
1.671 + plot = gsn_csm_xy (wks,npp933,nppmod933,res) ; create plot
1.672 +;-------------------------------
1.673 +; add polyline
1.674 +
1.675 + dum=gsn_add_polyline(wks,plot,(/0,1500/),(/0,1500/),False)
1.676 +;-------------------------------
1.677 + draw (plot)
1.678 + delete (wks)
1.679 + delete (dum)
1.680 +
1.681 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.682 + "rm "+plot_name+"."+plot_type)
1.683 +
1.684 +;**************************************************************************
1.685 +;(B) histogram 81
1.686 +;**************************************************************************
1.687 +;get data
1.688 + RAIN1_1D = ndtooned(rain81)
1.689 + RAIN2_1D = ndtooned(rainmod81)
1.690 + NPP1_1D = ndtooned(npp81)
1.691 + NPP2_1D = ndtooned(nppmod81)
1.692 +
1.693 +; number of bin
1.694 + nx = 8
1.695 +
1.696 + xvalues = new((/2,nx/),float)
1.697 + yvalues = new((/2,nx/),float)
1.698 + mn_yvalues = new((/2,nx/),float)
1.699 + mx_yvalues = new((/2,nx/),float)
1.700 + dx4 = new((/1/),float)
1.701 +
1.702 + get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
1.703 + ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
1.704 +
1.705 +;----------------------------------------
1.706 +;compute correlation coeff and M score
1.707 +
1.708 + u = yvalues(0,:)
1.709 + v = yvalues(1,:)
1.710 +
1.711 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.712 + uu = u(good)
1.713 + vv = v(good)
1.714 +
1.715 + ccr81h = esccr(uu,vv,0)
1.716 +
1.717 + score_max = 2.5
1.718 +
1.719 + bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
1.720 + M81h = (1.- (bias/dimsizes(uu)))*score_max
1.721 + M_npp_H81 = sprintf("%.2f", M81h)
1.722 +
1.723 + if (isvar("compare")) then
1.724 + system("sed -e '1,/M_npp_H81/s/M_npp_H81/"+M_npp_H81+"/' "+html_name2+" > "+html_new2+";"+ \
1.725 + "mv -f "+html_new2+" "+html_name2)
1.726 + end if
1.727 +
1.728 + system("sed s#M_npp_H81#"+M_npp_H81+"# "+html_name+" > "+html_new+";"+ \
1.729 + "mv -f "+html_new+" "+html_name)
1.730 +
1.731 + delete (u)
1.732 + delete (v)
1.733 + delete (uu)
1.734 + delete (vv)
1.735 +;----------------------------------------------------------------------
1.736 +; histogram res
1.737 + resh = True
1.738 + resh@gsnMaximize = True
1.739 + resh@gsnDraw = False
1.740 + resh@gsnFrame = False
1.741 + resh@xyMarkLineMode = "Markers"
1.742 + resh@xyMarkerSizeF = 0.014
1.743 + resh@xyMarker = 16
1.744 + resh@xyMarkerColors = (/"Brown","Blue"/)
1.745 + resh@trYMinF = min(mn_yvalues) - 10.
1.746 + resh@trYMaxF = max(mx_yvalues) + 10.
1.747 +
1.748 + resh@tiYAxisString = "NPP (g C/m2/year)"
1.749 + resh@tiXAxisString = "Precipitation (m/year)"
1.750 +
1.751 + max_bar = new((/2,nx/),graphic)
1.752 + min_bar = new((/2,nx/),graphic)
1.753 + max_cap = new((/2,nx/),graphic)
1.754 + min_cap = new((/2,nx/),graphic)
1.755 +
1.756 + lnres = True
1.757 + line_colors = (/"brown","blue"/)
1.758 +
1.759 +;**************************************************************************
1.760 +;(B)-1 histogram plot, ob 81 site
1.761 +;**************************************************************************
1.762 +
1.763 + plot_name = "histogram_ob_81"
1.764 + title = "Class A Observations (81 sites)"
1.765 + resh@tiMainString = title
1.766 +
1.767 + wks = gsn_open_wks (plot_type,plot_name)
1.768 +
1.769 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1.770 +
1.771 +;-------------------------------
1.772 +;Attach the vertical bar and the horizontal cap line
1.773 +
1.774 + do nd=0,0
1.775 + lnres@gsLineColor = line_colors(nd)
1.776 + do i=0,nx-1
1.777 +
1.778 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.779 + .not.ismissing(mx_yvalues(nd,i))) then
1.780 +;
1.781 +; Attach the vertical bar, both above and below the marker.
1.782 +;
1.783 + x1 = xvalues(nd,i)
1.784 + y1 = yvalues(nd,i)
1.785 + y2 = mn_yvalues(nd,i)
1.786 + plx = (/x1,x1/)
1.787 + ply = (/y1,y2/)
1.788 + min_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.789 +
1.790 + y2 = mx_yvalues(nd,i)
1.791 + plx = (/x1,x1/)
1.792 + ply = (/y1,y2/)
1.793 + max_bar(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.794 +;
1.795 +; Attach the horizontal cap line, both above and below the marker.
1.796 +;
1.797 + x1 = xvalues(nd,i) - dx4
1.798 + x2 = xvalues(nd,i) + dx4
1.799 + y1 = mn_yvalues(nd,i)
1.800 + plx = (/x1,x2/)
1.801 + ply = (/y1,y1/)
1.802 + min_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.803 +
1.804 + y1 = mx_yvalues(nd,i)
1.805 + plx = (/x1,x2/)
1.806 + ply = (/y1,y1/)
1.807 + max_cap(nd,i) = gsn_add_polyline(wks,xy,plx,ply,lnres)
1.808 + end if
1.809 + end do
1.810 + end do
1.811 +
1.812 + draw(xy)
1.813 + delete (wks)
1.814 +
1.815 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.816 + "rm "+plot_name+"."+plot_type)
1.817 +
1.818 +;****************************************************************************
1.819 +;(B)-2 histogram plot, model vs ob 81 site
1.820 +;****************************************************************************
1.821 +
1.822 + plot_name = "histogram_model_vs_ob_81"
1.823 + title = model_name+ " vs Class A Observations (81 sites)"
1.824 + resh@tiMainString = title
1.825 +
1.826 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.827 +
1.828 +;-----------------------------
1.829 +; Add a boxed legend using the more simple method, which won't have
1.830 +; vertical lines going through the markers.
1.831 +
1.832 + resh@pmLegendDisplayMode = "Always"
1.833 +; resh@pmLegendWidthF = 0.1
1.834 + resh@pmLegendWidthF = 0.08
1.835 + resh@pmLegendHeightF = 0.05
1.836 + resh@pmLegendOrthogonalPosF = -1.17
1.837 +; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.838 +; resh@pmLegendParallelPosF = 0.18
1.839 + resh@pmLegendParallelPosF = 0.88 ;(rightward)
1.840 +
1.841 +; resh@lgPerimOn = False
1.842 + resh@lgLabelFontHeightF = 0.015
1.843 + resh@xyExplicitLegendLabels = (/"observed",model_name/)
1.844 +;-----------------------------
1.845 + tRes = True
1.846 + tRes@txFontHeightF = 0.025
1.847 +
1.848 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr81h)+")"
1.849 +
1.850 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.851 +
1.852 + xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1.853 +;-------------------------------
1.854 +;Attach the vertical bar and the horizontal cap line
1.855 +
1.856 + do nd=0,1
1.857 + lnres@gsLineColor = line_colors(nd)
1.858 + do i=0,nx-1
1.859 +
1.860 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.861 + .not.ismissing(mx_yvalues(nd,i))) then
1.862 +;
1.863 +; Attach the vertical bar, both above and below the marker.
1.864 +;
1.865 + x1 = xvalues(nd,i)
1.866 + y1 = yvalues(nd,i)
1.867 + y2 = mn_yvalues(nd,i)
1.868 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.869 +
1.870 + y2 = mx_yvalues(nd,i)
1.871 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.872 +;
1.873 +; Attach the horizontal cap line, both above and below the marker.
1.874 +;
1.875 + x1 = xvalues(nd,i) - dx4
1.876 + x2 = xvalues(nd,i) + dx4
1.877 + y1 = mn_yvalues(nd,i)
1.878 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.879 +
1.880 + y1 = mx_yvalues(nd,i)
1.881 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.882 + end if
1.883 + end do
1.884 + end do
1.885 +
1.886 + draw(xy)
1.887 + delete(wks)
1.888 +
1.889 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.890 + "rm "+plot_name+"."+plot_type)
1.891 +
1.892 + delete (RAIN1_1D)
1.893 + delete (RAIN2_1D)
1.894 + delete (NPP1_1D)
1.895 + delete (NPP2_1D)
1.896 +;delete (range)
1.897 + delete (xvalues)
1.898 + delete (yvalues)
1.899 + delete (mn_yvalues)
1.900 + delete (mx_yvalues)
1.901 + delete (good)
1.902 + delete (max_bar)
1.903 + delete (min_bar)
1.904 + delete (max_cap)
1.905 + delete (min_cap)
1.906 +
1.907 +;**************************************************************************
1.908 +;(B) histogram 933
1.909 +;**************************************************************************
1.910 +
1.911 +; get data
1.912 + RAIN1_1D = ndtooned(rain933)
1.913 + RAIN2_1D = ndtooned(rainmod933)
1.914 + NPP1_1D = ndtooned(npp933)
1.915 + NPP2_1D = ndtooned(nppmod933)
1.916 +
1.917 +; number of bin
1.918 + nx = 10
1.919 +
1.920 + xvalues = new((/2,nx/),float)
1.921 + yvalues = new((/2,nx/),float)
1.922 + mn_yvalues = new((/2,nx/),float)
1.923 + mx_yvalues = new((/2,nx/),float)
1.924 + dx4 = new((/1/),float)
1.925 +
1.926 + get_bin(RAIN1_1D, NPP1_1D, RAIN2_1D, NPP2_1D \
1.927 + ,xvalues,yvalues,mn_yvalues,mx_yvalues,dx4)
1.928 +
1.929 +;----------------------------------------
1.930 +;compute correlation coeff and M score
1.931 +
1.932 + u = yvalues(0,:)
1.933 + v = yvalues(1,:)
1.934 +
1.935 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.936 + uu = u(good)
1.937 + vv = v(good)
1.938 +
1.939 + ccr933h = esccr(uu,vv,0)
1.940 +
1.941 + score_max = 2.5
1.942 +
1.943 + bias = sum(abs(vv-uu)/(abs(vv)+abs(uu)))
1.944 + M933h = (1.- (bias/dimsizes(uu)))*score_max
1.945 + M_npp_H933 = sprintf("%.2f", M933h)
1.946 +
1.947 + if (isvar("compare")) then
1.948 + system("sed -e '1,/M_npp_H933/s/M_npp_H933/"+M_npp_H933+"/' "+html_name2+" > "+html_new2+";"+ \
1.949 + "mv -f "+html_new2+" "+html_name2)
1.950 + end if
1.951 +
1.952 + system("sed s#M_npp_H933#"+M_npp_H933+"# "+html_name+" > "+html_new+";"+ \
1.953 + "mv -f "+html_new+" "+html_name)
1.954 +
1.955 + delete (u)
1.956 + delete (v)
1.957 + delete (uu)
1.958 + delete (vv)
1.959 +;----------------------------------------------------------------------
1.960 +; histogram res
1.961 + delete (resh)
1.962 + resh = True
1.963 + resh@gsnMaximize = True
1.964 + resh@gsnDraw = False
1.965 + resh@gsnFrame = False
1.966 + resh@xyMarkLineMode = "Markers"
1.967 + resh@xyMarkerSizeF = 0.014
1.968 + resh@xyMarker = 16
1.969 + resh@xyMarkerColors = (/"Brown","Blue"/)
1.970 + resh@trYMinF = min(mn_yvalues) - 10.
1.971 + resh@trYMaxF = max(mx_yvalues) + 10.
1.972 +
1.973 + resh@tiYAxisString = "NPP (g C/m2/year)"
1.974 + resh@tiXAxisString = "Precipitation (m/year)"
1.975 +
1.976 + max_bar = new((/2,nx/),graphic)
1.977 + min_bar = new((/2,nx/),graphic)
1.978 + max_cap = new((/2,nx/),graphic)
1.979 + min_cap = new((/2,nx/),graphic)
1.980 +
1.981 + lnres = True
1.982 + line_colors = (/"brown","blue"/)
1.983 +
1.984 +;**************************************************************************
1.985 +;(B)-3 histogram plot, ob 933 site
1.986 +;**************************************************************************
1.987 +
1.988 + plot_name = "histogram_ob_933"
1.989 + title = "Class B Observations (933 sites)"
1.990 + resh@tiMainString = title
1.991 +
1.992 + wks = gsn_open_wks (plot_type,plot_name)
1.993 +
1.994 + xy = gsn_csm_xy(wks,xvalues(0,:),yvalues(0,:),resh)
1.995 +
1.996 +;-------------------------------
1.997 +;Attach the vertical bar and the horizontal cap line
1.998 +
1.999 + do nd=0,0
1.1000 + lnres@gsLineColor = line_colors(nd)
1.1001 + do i=0,nx-1
1.1002 +
1.1003 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1004 + .not.ismissing(mx_yvalues(nd,i))) then
1.1005 +
1.1006 +; Attach the vertical bar, both above and below the marker.
1.1007 +
1.1008 + x1 = xvalues(nd,i)
1.1009 + y1 = yvalues(nd,i)
1.1010 + y2 = mn_yvalues(nd,i)
1.1011 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1012 +
1.1013 + y2 = mx_yvalues(nd,i)
1.1014 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1015 +
1.1016 +; Attach the horizontal cap line, both above and below the marker.
1.1017 +
1.1018 + x1 = xvalues(nd,i) - dx4
1.1019 + x2 = xvalues(nd,i) + dx4
1.1020 + y1 = mn_yvalues(nd,i)
1.1021 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1022 +
1.1023 + y1 = mx_yvalues(nd,i)
1.1024 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1025 + end if
1.1026 + end do
1.1027 + end do
1.1028 +
1.1029 + draw(xy)
1.1030 + delete (xy)
1.1031 + delete (wks)
1.1032 +
1.1033 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1034 + "rm "+plot_name+"."+plot_type)
1.1035 +
1.1036 +;****************************************************************************
1.1037 +;(B)-4 histogram plot, model vs ob 933 site
1.1038 +;****************************************************************************
1.1039 +
1.1040 + plot_name = "histogram_model_vs_ob_933"
1.1041 + title = model_name+ " vs Class B Observations (933 sites)"
1.1042 + resh@tiMainString = title
1.1043 +
1.1044 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1045 +
1.1046 +;-----------------------------
1.1047 +; Add a boxed legend using the more simple method, which won't have
1.1048 +; vertical lines going through the markers.
1.1049 +
1.1050 + resh@pmLegendDisplayMode = "Always"
1.1051 +; resh@pmLegendWidthF = 0.1
1.1052 + resh@pmLegendWidthF = 0.08
1.1053 + resh@pmLegendHeightF = 0.05
1.1054 + resh@pmLegendOrthogonalPosF = -1.17
1.1055 +; resh@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.1056 +; resh@pmLegendParallelPosF = 0.18
1.1057 + resh@pmLegendParallelPosF = 0.88 ;(rightward)
1.1058 +
1.1059 +; resh@lgPerimOn = False
1.1060 + resh@lgLabelFontHeightF = 0.015
1.1061 + resh@xyExplicitLegendLabels = (/"observed",model_name/)
1.1062 +;-----------------------------
1.1063 + tRes = True
1.1064 + tRes@txFontHeightF = 0.025
1.1065 +
1.1066 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccr933h)+")"
1.1067 +
1.1068 + gsn_text_ndc(wks,correlation_text,0.56,0.85,tRes)
1.1069 +
1.1070 + xy = gsn_csm_xy(wks,xvalues,yvalues,resh)
1.1071 +;-------------------------------
1.1072 +;Attach the vertical bar and the horizontal cap line
1.1073 +
1.1074 + do nd=0,1
1.1075 + lnres@gsLineColor = line_colors(nd)
1.1076 + do i=0,nx-1
1.1077 +
1.1078 + if(.not.ismissing(mn_yvalues(nd,i)).and. \
1.1079 + .not.ismissing(mx_yvalues(nd,i))) then
1.1080 +;
1.1081 +; Attach the vertical bar, both above and below the marker.
1.1082 +;
1.1083 + x1 = xvalues(nd,i)
1.1084 + y1 = yvalues(nd,i)
1.1085 + y2 = mn_yvalues(nd,i)
1.1086 + min_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1087 +
1.1088 + y2 = mx_yvalues(nd,i)
1.1089 + max_bar(nd,i) = gsn_add_polyline(wks,xy,(/x1,x1/),(/y1,y2/),lnres)
1.1090 +;
1.1091 +; Attach the horizontal cap line, both above and below the marker.
1.1092 +;
1.1093 + x1 = xvalues(nd,i) - dx4
1.1094 + x2 = xvalues(nd,i) + dx4
1.1095 + y1 = mn_yvalues(nd,i)
1.1096 + min_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1097 +
1.1098 + y1 = mx_yvalues(nd,i)
1.1099 + max_cap(nd,i) = gsn_add_polyline(wks,xy,(/x1,x2/),(/y1,y1/),lnres)
1.1100 + end if
1.1101 + end do
1.1102 + end do
1.1103 +
1.1104 + draw(xy)
1.1105 + delete(xy)
1.1106 + delete(wks)
1.1107 +
1.1108 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1109 + "rm "+plot_name+"."+plot_type)
1.1110 +
1.1111 +;***************************************************************************
1.1112 +;(C) global contour
1.1113 +;***************************************************************************
1.1114 +
1.1115 +;res
1.1116 + resg = True ; Use plot options
1.1117 + resg@cnFillOn = True ; Turn on color fill
1.1118 + resg@gsnSpreadColors = True ; use full colormap
1.1119 +; resg@cnFillMode = "RasterFill" ; Turn on raster color
1.1120 +; resg@lbLabelAutoStride = True
1.1121 + resg@cnLinesOn = False ; Turn off contourn lines
1.1122 + resg@mpFillOn = False ; Turn off map fill
1.1123 +
1.1124 + resg@gsnSpreadColors = True ; use full colormap
1.1125 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.1126 + resg@cnMinLevelValF = 0. ; Min level
1.1127 + resg@cnMaxLevelValF = 2200. ; Max level
1.1128 + resg@cnLevelSpacingF = 200. ; interval
1.1129 +
1.1130 +;****************************************************************************
1.1131 +;(C)-1 global contour plot, ob
1.1132 +;****************************************************************************
1.1133 +
1.1134 + delta = 0.00000000001
1.1135 + nppglobe = where(ismissing(nppglobe).and.(ismissing(nppmod).or.(nppmod.lt.delta)),0.,nppglobe)
1.1136 +
1.1137 + plot_name = "global_ob"
1.1138 + title = "Observed MODIS MOD 17"
1.1139 + resg@tiMainString = title
1.1140 +
1.1141 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1142 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1143 +
1.1144 + plot = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1.1145 +
1.1146 + delete (wks)
1.1147 +
1.1148 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1149 + "rm "+plot_name+"."+plot_type)
1.1150 +;****************************************************************************
1.1151 +;(C)-2 global contour plot, model
1.1152 +;****************************************************************************
1.1153 +
1.1154 + plot_name = "global_model"
1.1155 + title = "Model "+ model_name
1.1156 + resg@tiMainString = title
1.1157 +
1.1158 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1159 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1160 +
1.1161 + plot = gsn_csm_contour_map_ce(wks,nppmod,resg)
1.1162 +
1.1163 + delete (wks)
1.1164 +
1.1165 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1166 + "rm "+plot_name+"."+plot_type)
1.1167 +
1.1168 +;****************************************************************************
1.1169 +;(C)-3 global contour plot, model vs ob
1.1170 +;****************************************************************************
1.1171 +
1.1172 + plot_name = "global_model_vs_ob"
1.1173 +
1.1174 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1175 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.1176 +
1.1177 + delete (plot)
1.1178 + plot=new(3,graphic) ; create graphic array
1.1179 +
1.1180 + resg@gsnFrame = False ; Do not draw plot
1.1181 + resg@gsnDraw = False ; Do not advance frame
1.1182 +
1.1183 +; compute correlation coef and M score
1.1184 +
1.1185 + uu1 = ndtooned(nppmod)
1.1186 + vv1 = ndtooned(nppglobe)
1.1187 +
1.1188 + delete (good)
1.1189 + good = ind(.not.ismissing(uu1) .and. .not.ismissing(vv1))
1.1190 +
1.1191 + ug = uu1(good)
1.1192 + vg = vv1(good)
1.1193 +
1.1194 + ccrG = esccr(ug,vg,0)
1.1195 +
1.1196 + score_max = 5.0
1.1197 +
1.1198 + MG = (ccrG*ccrG)* score_max
1.1199 + M_npp_G = sprintf("%.2f", MG)
1.1200 +
1.1201 + if (isvar("compare")) then
1.1202 + system("sed -e '1,/M_npp_G/s/M_npp_G/"+M_npp_G+"/' "+html_name2+" > "+html_new2+";"+ \
1.1203 + "mv -f "+html_new2+" "+html_name2)
1.1204 + end if
1.1205 +
1.1206 + system("sed s#M_npp_G#"+M_npp_G+"# "+html_name+" > "+html_new+";"+ \
1.1207 + "mv -f "+html_new+" "+html_name)
1.1208 +
1.1209 +; plot correlation coef
1.1210 +
1.1211 + gRes = True
1.1212 + gRes@txFontHeightF = 0.02
1.1213 + gRes@txAngleF = 90
1.1214 +
1.1215 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrG)+")"
1.1216 +
1.1217 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.1218 +;--------------------------------------------------------------------
1.1219 +;(a) ob
1.1220 +
1.1221 + title = "Observed MODIS MOD 17"
1.1222 + resg@tiMainString = title
1.1223 +
1.1224 + plot(0) = gsn_csm_contour_map_ce(wks,nppglobe,resg)
1.1225 +
1.1226 +;(b) model
1.1227 +
1.1228 + title = "Model "+ model_name
1.1229 + resg@tiMainString = title
1.1230 +
1.1231 + plot(1) = gsn_csm_contour_map_ce(wks,nppmod,resg)
1.1232 +
1.1233 +;(c) model-ob
1.1234 +
1.1235 + zz = nppmod
1.1236 + zz = nppmod - nppglobe
1.1237 + title = "Model_"+model_name+" - Observed"
1.1238 +
1.1239 + resg@cnMinLevelValF = -500 ; Min level
1.1240 + resg@cnMaxLevelValF = 500. ; Max level
1.1241 + resg@cnLevelSpacingF = 50. ; interval
1.1242 + resg@tiMainString = title
1.1243 +
1.1244 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.1245 +
1.1246 + pres = True ; panel plot mods desired
1.1247 +; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around
1.1248 + ; indiv. plots in panel
1.1249 + pres@gsnMaximize = True ; fill the page
1.1250 +
1.1251 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.1252 +
1.1253 + delete (wks)
1.1254 + delete (plot)
1.1255 +
1.1256 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1257 + "rm "+plot_name+"."+plot_type)
1.1258 +
1.1259 +;***************************************************************************
1.1260 +;(D)-1 zonal line plot, ob
1.1261 +;***************************************************************************
1.1262 +
1.1263 + vv = zonalAve(nppglobe)
1.1264 + vv@long_name = nppglobe@long_name
1.1265 +
1.1266 + plot_name = "zonal_ob"
1.1267 + title = "Observed MODIS MOD 17"
1.1268 +
1.1269 + resz = True
1.1270 + resz@tiMainString = title
1.1271 +
1.1272 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.1273 +
1.1274 + plot = gsn_csm_xy (wks,ym,vv,resz)
1.1275 +
1.1276 + delete (wks)
1.1277 +
1.1278 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1279 + "rm "+plot_name+"."+plot_type)
1.1280 +
1.1281 +;****************************************************************************
1.1282 +;(D)-2 zonal line plot, model vs ob
1.1283 +;****************************************************************************
1.1284 +
1.1285 + s = new ((/2,dimsizes(ym)/), float)
1.1286 +
1.1287 + s(0,:) = zonalAve(nppglobe)
1.1288 + s(1,:) = zonalAve(nppmod)
1.1289 +
1.1290 + s@long_name = nppglobe@long_name
1.1291 +;-------------------------------------------
1.1292 +; compute correlation coef and M score
1.1293 +
1.1294 + score_max = 5.0
1.1295 +
1.1296 + ccrZ = esccr(s(1,:), s(0,:),0)
1.1297 + MZ = (ccrZ*ccrZ)* score_max
1.1298 + M_npp_Z = sprintf("%.2f", MZ)
1.1299 +
1.1300 + if (isvar("compare")) then
1.1301 + system("sed -e '1,/M_npp_Z/s/M_npp_Z/"+M_npp_Z+"/' "+html_name2+" > "+html_new2+";"+ \
1.1302 + "mv -f "+html_new2+" "+html_name2)
1.1303 + end if
1.1304 +
1.1305 + system("sed s#M_npp_Z#"+M_npp_Z+"# "+html_name+" > "+html_new+";"+ \
1.1306 + "mv -f "+html_new+" "+html_name)
1.1307 +;-------------------------------------------
1.1308 + plot_name = "zonal_model_vs_ob"
1.1309 + title = "Zonal Average"
1.1310 + resz@tiMainString = title
1.1311 +
1.1312 + wks = gsn_open_wks (plot_type,plot_name)
1.1313 +
1.1314 +; resz@vpHeightF = 0.4 ; change aspect ratio of plot
1.1315 +; resz@vpWidthF = 0.7
1.1316 +
1.1317 + resz@xyMonoLineColor = "False" ; want colored lines
1.1318 + resz@xyLineColors = (/"black","red"/) ; colors chosen
1.1319 +; resz@xyLineThicknesses = (/3.,3./) ; line thicknesses
1.1320 + resz@xyLineThicknesses = (/2.,2./) ; line thicknesses
1.1321 + resz@xyDashPatterns = (/0.,0./) ; make all lines solid
1.1322 +
1.1323 + resz@tiYAxisString = s@long_name ; add a axis title
1.1324 + resz@txFontHeightF = 0.0195 ; change title font heights
1.1325 +
1.1326 +; Legent
1.1327 + resz@pmLegendDisplayMode = "Always" ; turn on legend
1.1328 + resz@pmLegendSide = "Top" ; Change location of
1.1329 +; resz@pmLegendParallelPosF = .45 ; move units right
1.1330 + resz@pmLegendParallelPosF = .82 ; move units right
1.1331 + resz@pmLegendOrthogonalPosF = -0.4 ; move units down
1.1332 +
1.1333 + resz@pmLegendWidthF = 0.10 ; Change width and
1.1334 + resz@pmLegendHeightF = 0.10 ; height of legend.
1.1335 + resz@lgLabelFontHeightF = .02 ; change font height
1.1336 +; resz@lgTitleOn = True ; turn on legend title
1.1337 +; resz@lgTitleString = "Example" ; create legend title
1.1338 +; resz@lgTitleFontHeightF = .025 ; font of legend title
1.1339 + resz@xyExplicitLegendLabels = (/"Observed",model_name/) ; explicit labels
1.1340 +;--------------------------------------------------------------------
1.1341 + zRes = True
1.1342 + zRes@txFontHeightF = 0.025
1.1343 +
1.1344 + correlation_text = "(correlation coef = "+sprintf("%5.2f", ccrZ)+")"
1.1345 +
1.1346 + gsn_text_ndc(wks,correlation_text,0.50,0.77,zRes)
1.1347 +;--------------------------------------------------------------------
1.1348 +
1.1349 + plot = gsn_csm_xy (wks,ym,s,resz)
1.1350 +
1.1351 + delete (wks)
1.1352 +
1.1353 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.1354 + "rm "+plot_name+"."+plot_type)
1.1355 +
1.1356 +;***************************************************************************
1.1357 +; add total score and write to file
1.1358 +;***************************************************************************
1.1359 + M_total = M81s + M81h + M933s + M933h + MG + MZ
1.1360 +
1.1361 + asciiwrite("M_save.npp", M_total)
1.1362 +
1.1363 +;***************************************************************************
1.1364 +; output plots
1.1365 +;***************************************************************************
1.1366 + output_dir = model_name+"/npp"
1.1367 +
1.1368 + system("mv table_*.html " + output_dir +";"+ \
1.1369 + "mv *.png " + output_dir)
1.1370 +
1.1371 +;***************************************************************************
1.1372 +end
1.1373 +