1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/all/07.beta.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,541 @@
1.4 +;********************************************************
1.5 +; hardwired: co2_i = 283.1878
1.6 +; co2_f = 364.1252
1.7 +;
1.8 +; beta_4_ob = 0.6
1.9 +;**************************************************************
1.10 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.11 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.12 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.13 +;**************************************************************
1.14 +procedure set_line(lines:string,nline:integer,newlines:string)
1.15 +begin
1.16 +; add line to ascci/html file
1.17 +
1.18 + nnewlines = dimsizes(newlines)
1.19 + if(nline+nnewlines-1.ge.dimsizes(lines))
1.20 + print("set_line: bad index, not setting anything.")
1.21 + return
1.22 + end if
1.23 + lines(nline:nline+nnewlines-1) = newlines
1.24 +; print ("lines = " + lines(nline:nline+nnewlines-1))
1.25 + nline = nline + nnewlines
1.26 + return
1.27 +end
1.28 +;**************************************************************
1.29 +; Main code.
1.30 +begin
1.31 +
1.32 + plot_type = "ps"
1.33 + plot_type_new = "png"
1.34 +
1.35 +;------------------------------------------------------
1.36 +; edit table.html of current model for movel1_vs_model2
1.37 +
1.38 + if (isvar("compare")) then
1.39 + html_name2 = compare+"/table.html"
1.40 + html_new2 = html_name2 +".new"
1.41 + end if
1.42 +
1.43 +;-------------------------------------------------------
1.44 +; edit table.html for current model
1.45 +
1.46 + html_name = model_name+"/table.html"
1.47 + html_new = html_name +".new"
1.48 +
1.49 +;-------------------------------------------------------
1.50 +; read model data
1.51 +
1.52 +;###############################################################
1.53 +; hardwired for model data
1.54 +
1.55 +; these values correspond to the start and the end of model data
1.56 + co2_i = 283.1878
1.57 + co2_f = 364.1252
1.58 +
1.59 + film_i = film6
1.60 + film_f = film5
1.61 +;###############################################################
1.62 +
1.63 + fm_i = addfile (dirm+film_i,"r")
1.64 + fm_f = addfile (dirm+film_f,"r")
1.65 +
1.66 + xm = fm_f->lon
1.67 + ym = fm_f->lat
1.68 +
1.69 + npp_i = fm_i->NPP
1.70 + npp_f = fm_f->NPP
1.71 +
1.72 + delete (fm_i)
1.73 + delete (fm_f)
1.74 +
1.75 +;Units for these variables are:
1.76 +;npp_i: g C/m^2/s
1.77 +
1.78 + nsec_per_year = 60*60*24*365
1.79 +
1.80 + npp_i = npp_i * nsec_per_year
1.81 + npp_f = npp_f * nsec_per_year
1.82 +
1.83 + unit = "gC/m2/year"
1.84 +
1.85 +;------------------------------
1.86 +; get landfrac data
1.87 +
1.88 + film_l = "lnd_" + model_grid + ".nc"
1.89 + fm_l = addfile (dirs+film_l,"r")
1.90 + landfrac = fm_l->landfrac
1.91 +
1.92 + npp_i(0,:,:) = npp_i(0,:,:) * landfrac(:,:)
1.93 + npp_f(0,:,:) = npp_f(0,:,:) * landfrac(:,:)
1.94 +
1.95 + delete (fm_l)
1.96 + delete (landfrac)
1.97 +
1.98 +;-----------------------------
1.99 +; read biome data: model
1.100 +
1.101 + biome_name_mod = "Model PFT Class"
1.102 +
1.103 + film_c = "class_pft_"+model_grid+".nc"
1.104 + fm_c = addfile (dirs+film_c,"r")
1.105 + classmod = fm_c->CLASS_PFT
1.106 +
1.107 + delete (fm_c)
1.108 +
1.109 +; model data has 17 land-type classes
1.110 + nclass_mod = 17
1.111 +
1.112 +;---------------------------------------------------
1.113 +; read data: observed at stations
1.114 +
1.115 + station = (/"DukeFACE" \
1.116 + ,"AspenFACE" \
1.117 + ,"ORNL-FACE" \
1.118 + ,"POP-EUROFACE" \
1.119 + /)
1.120 +
1.121 + lat_ob = (/ 35.58, 45.40, 35.54, 42.22/)
1.122 + lon_ob = (/-79.05, -89.37, -84.20, 11.48/)
1.123 + lon_obx = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)
1.124 +
1.125 + n_sta = dimsizes(station)
1.126 + beta_4_ob = new((/n_sta/),float)
1.127 +
1.128 +;###################################################
1.129 +; this is a hardwired value
1.130 + beta_4_ob = 0.60
1.131 +;###################################################
1.132 +;---------------------------------------------------
1.133 +; get model data at station
1.134 +
1.135 + npp_i_4 =linint2_points(xm,ym,npp_i,True,lon_obx,lat_ob,0)
1.136 +
1.137 + npp_f_4 =linint2_points(xm,ym,npp_f,True,lon_obx,lat_ob,0)
1.138 +
1.139 +;---------------------------------------------------
1.140 +;compute beta_4
1.141 +
1.142 + score_max = 3.
1.143 +
1.144 + beta_4 = new((/n_sta/),float)
1.145 +
1.146 + beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i)
1.147 +
1.148 + beta_4_avg = avg(beta_4)
1.149 +
1.150 + bias = sum(abs(beta_4-beta_4_ob)/(abs(beta_4)+abs(beta_4_ob)))
1.151 + Mbeta = (1. - (bias/n_sta))*score_max
1.152 + M_beta = sprintf("%.2f", Mbeta)
1.153 +
1.154 +;=========================
1.155 +; for html table - station
1.156 +;=========================
1.157 +
1.158 + output_html = "table_station.html"
1.159 +
1.160 +; column (not including header column)
1.161 +
1.162 + col_head = (/"Latitude","Longitude","CO2_i","CO2_f","NPP_i","NPP_f","Beta_model","Beta_ob"/)
1.163 +
1.164 + ncol = dimsizes(col_head)
1.165 +
1.166 +; row (not including header row)
1.167 + row_head = (/"DukeFACE" \
1.168 + ,"AspenFACE" \
1.169 + ,"ORNL-FACE" \
1.170 + ,"POP-EUROFACE" \
1.171 + ,"All Station" \
1.172 + /)
1.173 + nrow = dimsizes(row_head)
1.174 +
1.175 +; arrays to be passed to table.
1.176 + text = new ((/nrow, ncol/),string )
1.177 +
1.178 + do i=0,nrow-2
1.179 + text(i,0) = sprintf("%.1f",lat_ob(i))
1.180 + text(i,1) = sprintf("%.1f",lon_ob(i))
1.181 + text(i,2) = sprintf("%.1f",co2_i)
1.182 + text(i,3) = sprintf("%.1f",co2_f)
1.183 + text(i,4) = sprintf("%.1f",npp_i_4(0,i))
1.184 + text(i,5) = sprintf("%.1f",npp_f_4(0,i))
1.185 + text(i,6) = sprintf("%.2f",beta_4(i))
1.186 + text(i,7) = "-"
1.187 + end do
1.188 + text(nrow-1,0) = "-"
1.189 + text(nrow-1,1) = "-"
1.190 + text(nrow-1,2) = "-"
1.191 + text(nrow-1,3) = "-"
1.192 + text(nrow-1,4) = "-"
1.193 + text(nrow-1,5) = "-"
1.194 + text(nrow-1,6) = sprintf("%.2f",beta_4_avg)
1.195 + text(nrow-1,7) = sprintf("%.2f",avg(beta_4_ob))
1.196 +
1.197 +;-----------
1.198 +; html table
1.199 +;-----------
1.200 +
1.201 + header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
1.202 +
1.203 + header = (/"<HTML>" \
1.204 + ,"<HEAD>" \
1.205 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.206 + ,"</HEAD>" \
1.207 + ,header_text \
1.208 + /)
1.209 + footer = "</HTML>"
1.210 +
1.211 + table_header = (/ \
1.212 + "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
1.213 + ,"<tr>" \
1.214 + ," <th bgcolor=DDDDDD >Station</th>" \
1.215 + ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
1.216 + ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
1.217 + ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
1.218 + ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
1.219 + ," <th bgcolor=DDDDDD >"+col_head(4)+"<br>("+unit+")</th>" \
1.220 + ," <th bgcolor=DDDDDD >"+col_head(5)+"<br>("+unit+")</th>" \
1.221 + ," <th bgcolor=DDDDDD >"+col_head(6)+"</th>" \
1.222 + ," <th bgcolor=DDDDDD >"+col_head(7)+"</th>" \
1.223 + ,"</tr>" \
1.224 + /)
1.225 + table_footer = "</table>"
1.226 + row_header = "<tr>"
1.227 + row_footer = "</tr>"
1.228 +
1.229 + lines = new(50000,string)
1.230 + nline = 0
1.231 +
1.232 + set_line(lines,nline,header)
1.233 + set_line(lines,nline,table_header)
1.234 +;-----------------------------------------------
1.235 +;row of table
1.236 +
1.237 + do n = 0,nrow-1
1.238 + set_line(lines,nline,row_header)
1.239 +
1.240 + txt1 = row_head(n)
1.241 + txt2 = text(n,0)
1.242 + txt3 = text(n,1)
1.243 + txt4 = text(n,2)
1.244 + txt5 = text(n,3)
1.245 + txt6 = text(n,4)
1.246 + txt7 = text(n,5)
1.247 + txt8 = text(n,6)
1.248 + txt9 = text(n,7)
1.249 +
1.250 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.251 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.252 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.253 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.254 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.255 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.256 + set_line(lines,nline,"<th>"+txt7+"</th>")
1.257 + set_line(lines,nline,"<th>"+txt8+"</th>")
1.258 + set_line(lines,nline,"<th>"+txt9+"</th>")
1.259 +
1.260 + set_line(lines,nline,row_footer)
1.261 + end do
1.262 +;-----------------------------------------------
1.263 + set_line(lines,nline,table_footer)
1.264 + set_line(lines,nline,footer)
1.265 +
1.266 +; Now write to an HTML file.
1.267 + idx = ind(.not.ismissing(lines))
1.268 + if(.not.any(ismissing(idx))) then
1.269 + asciiwrite(output_html,lines(idx))
1.270 + else
1.271 + print ("error?")
1.272 + end if
1.273 +
1.274 + delete (col_head)
1.275 + delete (row_head)
1.276 + delete (text)
1.277 + delete (table_header)
1.278 + delete (idx)
1.279 +
1.280 +;********************************************************************
1.281 +; use land-type class to bin the data in equally spaced ranges
1.282 +;********************************************************************
1.283 +
1.284 +; using model biome class
1.285 + nclass = nclass_mod
1.286 +
1.287 + range = fspan(0,nclass,nclass+1)
1.288 +
1.289 +; Use this range information to grab all the values in a
1.290 +; particular range, and then take an average.
1.291 +
1.292 + nx = dimsizes(range) - 1
1.293 +
1.294 +;==============================
1.295 +; put data into bins
1.296 +;==============================
1.297 +
1.298 +; for model data and observed
1.299 + data_n = 2
1.300 +
1.301 +; using model biome class
1.302 +
1.303 + base = ndtooned(classmod)
1.304 +
1.305 +; output
1.306 +
1.307 + yvalues = new((/data_n,nx/),float)
1.308 + count = new((/data_n,nx/),float)
1.309 +
1.310 +; Loop through each range, using base
1.311 +
1.312 + do i=0,nx-1
1.313 +
1.314 + if (i.ne.(nx-1)) then
1.315 + idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
1.316 + else
1.317 + idx = ind(base.ge.range(i))
1.318 + end if
1.319 +
1.320 +; loop through each dataset
1.321 +
1.322 + do n = 0,data_n-1
1.323 +
1.324 + if (n .eq. 0) then
1.325 + data = ndtooned(npp_i)
1.326 + end if
1.327 +
1.328 + if (n .eq. 1) then
1.329 + data = ndtooned(npp_f)
1.330 + end if
1.331 +
1.332 +; Calculate average
1.333 +
1.334 + if (.not.any(ismissing(idx))) then
1.335 + yvalues(n,i) = avg(data(idx))
1.336 + count(n,i) = dimsizes(idx)
1.337 + else
1.338 + yvalues(n,i) = yvalues@_FillValue
1.339 + count(n,i) = 0
1.340 + end if
1.341 +
1.342 +;#############################################################
1.343 +; using model biome class:
1.344 +;
1.345 +; set the following 4 classes to _FillValue:
1.346 +; (3)Needleleaf Deciduous Boreal Tree,
1.347 +; (8)Broadleaf Deciduous Boreal Tree,
1.348 +; (9)Broadleaf Evergreen Shrub,
1.349 +; (16)Wheat
1.350 +
1.351 + if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
1.352 + yvalues(n,i) = yvalues@_FillValue
1.353 + count(n,i) = 0
1.354 + end if
1.355 +;#############################################################
1.356 +
1.357 + delete(data)
1.358 + end do ; n-loop
1.359 +
1.360 + delete(idx)
1.361 + end do ; i-loop
1.362 +
1.363 + delete (base)
1.364 + delete (npp_i)
1.365 + delete (npp_f)
1.366 +
1.367 +;============================
1.368 +;compute beta
1.369 +;============================
1.370 +
1.371 + u = yvalues(0,:)
1.372 + v = yvalues(1,:)
1.373 + u_count = count(0,:)
1.374 + v_count = count(1,:)
1.375 +
1.376 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.377 +
1.378 + uu = u(good)
1.379 + vv = v(good)
1.380 + uu_count = u_count(good)
1.381 + vv_count = v_count(good)
1.382 +
1.383 + n_biome = dimsizes(uu)
1.384 + beta_biome = new((/n_biome/),float)
1.385 +
1.386 + beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)
1.387 +
1.388 + beta_biome_avg = (sum(vv*vv_count)/sum(uu*uu_count) - 1.)/log(co2_f/co2_i)
1.389 +
1.390 +;===========================
1.391 +; for html table - biome
1.392 +;===========================
1.393 +
1.394 + output_html = "table_biome.html"
1.395 +
1.396 +; column (not including header column)
1.397 +
1.398 + col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta_model"/)
1.399 +
1.400 + ncol = dimsizes(col_head)
1.401 +
1.402 +; row (not including header row)
1.403 +
1.404 +;----------------------------------------------------
1.405 +; using model biome class:
1.406 +;
1.407 + row_head = (/"Not Vegetated" \
1.408 + ,"Needleleaf Evergreen Temperate Tree" \
1.409 + ,"Needleleaf Evergreen Boreal Tree" \
1.410 +; ,"Needleleaf Deciduous Boreal Tree" \
1.411 + ,"Broadleaf Evergreen Tropical Tree" \
1.412 + ,"Broadleaf Evergreen Temperate Tree" \
1.413 + ,"Broadleaf Deciduous Tropical Tree" \
1.414 + ,"Broadleaf Deciduous Temperate Tree" \
1.415 +; ,"Broadleaf Deciduous Boreal Tree" \
1.416 +; ,"Broadleaf Evergreen Shrub" \
1.417 + ,"Broadleaf Deciduous Temperate Shrub" \
1.418 + ,"Broadleaf Deciduous Boreal Shrub" \
1.419 + ,"C3 Arctic Grass" \
1.420 + ,"C3 Non-Arctic Grass" \
1.421 + ,"C4 Grass" \
1.422 + ,"Corn" \
1.423 +; ,"Wheat" \
1.424 + ,"All Biome" \
1.425 + /)
1.426 +
1.427 + nrow = dimsizes(row_head)
1.428 +
1.429 +; arrays to be passed to table.
1.430 + text = new ((/nrow, ncol/),string )
1.431 +
1.432 + do i=0,nrow-2
1.433 + text(i,0) = sprintf("%.1f",co2_i)
1.434 + text(i,1) = sprintf("%.1f",co2_f)
1.435 + text(i,2) = sprintf("%.1f",uu(i))
1.436 + text(i,3) = sprintf("%.1f",vv(i))
1.437 + text(i,4) = sprintf("%.2f",beta_biome(i))
1.438 + end do
1.439 + text(nrow-1,0) = "-"
1.440 + text(nrow-1,1) = "-"
1.441 + text(nrow-1,2) = "-"
1.442 + text(nrow-1,3) = "-"
1.443 + text(nrow-1,4) = sprintf("%.2f",beta_biome_avg)
1.444 +
1.445 +;**************************************************
1.446 +; html table
1.447 +;**************************************************
1.448 +
1.449 + header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
1.450 +
1.451 + header = (/"<HTML>" \
1.452 + ,"<HEAD>" \
1.453 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.454 + ,"</HEAD>" \
1.455 + ,header_text \
1.456 + /)
1.457 + footer = "</HTML>"
1.458 +
1.459 + table_header = (/ \
1.460 + "<table border=1 cellspacing=0 cellpadding=3 width=80%>" \
1.461 + ,"<tr>" \
1.462 + ," <th bgcolor=DDDDDD >Biome Class</th>" \
1.463 + ," <th bgcolor=DDDDDD >"+col_head(0)+"</th>" \
1.464 + ," <th bgcolor=DDDDDD >"+col_head(1)+"</th>" \
1.465 + ," <th bgcolor=DDDDDD >"+col_head(2)+"</th>" \
1.466 + ," <th bgcolor=DDDDDD >"+col_head(3)+"</th>" \
1.467 + ," <th bgcolor=DDDDDD >"+col_head(4)+"</th>" \
1.468 + ,"</tr>" \
1.469 + /)
1.470 + table_footer = "</table>"
1.471 + row_header = "<tr>"
1.472 + row_footer = "</tr>"
1.473 +
1.474 + lines = new(50000,string)
1.475 + nline = 0
1.476 +
1.477 + set_line(lines,nline,header)
1.478 + set_line(lines,nline,table_header)
1.479 +;-----------------------------------------------
1.480 +;row of table
1.481 +
1.482 + do n = 0,nrow-1
1.483 + set_line(lines,nline,row_header)
1.484 +
1.485 + txt1 = row_head(n)
1.486 + txt2 = text(n,0)
1.487 + txt3 = text(n,1)
1.488 + txt4 = text(n,2)
1.489 + txt5 = text(n,3)
1.490 + txt6 = text(n,4)
1.491 +
1.492 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.493 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.494 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.495 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.496 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.497 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.498 +
1.499 + set_line(lines,nline,row_footer)
1.500 + end do
1.501 +;-----------------------------------------------
1.502 + set_line(lines,nline,table_footer)
1.503 + set_line(lines,nline,footer)
1.504 +
1.505 +; Now write to an HTML file
1.506 +
1.507 + idx = ind(.not.ismissing(lines))
1.508 + if(.not.any(ismissing(idx))) then
1.509 + asciiwrite(output_html,lines(idx))
1.510 + else
1.511 + print ("error?")
1.512 + end if
1.513 + delete (idx)
1.514 +
1.515 +;**************************************************************************************
1.516 +; update score
1.517 +;**************************************************************************************
1.518 + if (isvar("compare")) then
1.519 + system("sed -e '1,/M_beta/s/M_beta/"+M_beta+"/' "+html_name2+" > "+html_new2+";"+ \
1.520 + "mv -f "+html_new2+" "+html_name2)
1.521 + end if
1.522 +
1.523 + system("sed s#M_beta#"+M_beta+"# "+html_name+" > "+html_new+";"+ \
1.524 + "mv -f "+html_new+" "+html_name)
1.525 +
1.526 +;***************************************************************************
1.527 +; get total score and write to file
1.528 +;***************************************************************************
1.529 + M_total = Mbeta
1.530 +
1.531 + asciiwrite("M_save.beta", M_total)
1.532 +
1.533 + delete (M_total)
1.534 +
1.535 +;***************************************************************************
1.536 +; output plot and html
1.537 +;***************************************************************************
1.538 + output_dir = model_name+"/beta"
1.539 +
1.540 + system("mv *.html " + output_dir)
1.541 +;***************************************************************************
1.542 +
1.543 +end
1.544 +