forrest@0: ;********************************************************
forrest@0: ; hardwired: co2_i = 283.1878
forrest@0: ; co2_f = 364.1252
forrest@0: ;
forrest@0: ; beta_4_ob = 0.6
forrest@0: ;**************************************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0: ;**************************************************************
forrest@0: procedure set_line(lines:string,nline:integer,newlines:string)
forrest@0: begin
forrest@0: ; add line to ascci/html file
forrest@0:
forrest@0: nnewlines = dimsizes(newlines)
forrest@0: if(nline+nnewlines-1.ge.dimsizes(lines))
forrest@0: print("set_line: bad index, not setting anything.")
forrest@0: return
forrest@0: end if
forrest@0: lines(nline:nline+nnewlines-1) = newlines
forrest@0: ; print ("lines = " + lines(nline:nline+nnewlines-1))
forrest@0: nline = nline + nnewlines
forrest@0: return
forrest@0: end
forrest@0: ;**************************************************************
forrest@0: ; Main code.
forrest@0: begin
forrest@0:
forrest@0: plot_type = "ps"
forrest@0: plot_type_new = "png"
forrest@0:
forrest@0: ;------------------------------------------------------
forrest@0: ; edit table.html of current model for movel1_vs_model2
forrest@0:
forrest@0: if (isvar("compare")) then
forrest@0: html_name2 = compare+"/table.html"
forrest@0: html_new2 = html_name2 +".new"
forrest@0: end if
forrest@0:
forrest@0: ;-------------------------------------------------------
forrest@0: ; edit table.html for current model
forrest@0:
forrest@0: html_name = model_name+"/table.html"
forrest@0: html_new = html_name +".new"
forrest@0:
forrest@0: ;-------------------------------------------------------
forrest@0: ; read model data
forrest@0:
forrest@0: ;###############################################################
forrest@0: ; hardwired for model data
forrest@0:
forrest@0: ; these values correspond to the start and the end of model data
forrest@0: co2_i = 283.1878
forrest@0: co2_f = 364.1252
forrest@0:
forrest@0: film_i = film6
forrest@0: film_f = film5
forrest@0: ;###############################################################
forrest@0:
forrest@0: fm_i = addfile (dirm+film_i,"r")
forrest@0: fm_f = addfile (dirm+film_f,"r")
forrest@0:
forrest@0: xm = fm_f->lon
forrest@0: ym = fm_f->lat
forrest@0:
forrest@0: npp_i = fm_i->NPP
forrest@0: npp_f = fm_f->NPP
forrest@0:
forrest@0: delete (fm_i)
forrest@0: delete (fm_f)
forrest@0:
forrest@0: ;Units for these variables are:
forrest@0: ;npp_i: g C/m^2/s
forrest@0:
forrest@0: nsec_per_year = 60*60*24*365
forrest@0:
forrest@0: npp_i = npp_i * nsec_per_year
forrest@0: npp_f = npp_f * nsec_per_year
forrest@0:
forrest@0: unit = "gC/m2/year"
forrest@0:
forrest@0: ;------------------------------
forrest@0: ; get landfrac data
forrest@0:
forrest@0: film_l = "lnd_" + model_grid + ".nc"
forrest@0: fm_l = addfile (dirs+film_l,"r")
forrest@0: landfrac = fm_l->landfrac
forrest@0:
forrest@0: npp_i(0,:,:) = npp_i(0,:,:) * landfrac(:,:)
forrest@0: npp_f(0,:,:) = npp_f(0,:,:) * landfrac(:,:)
forrest@0:
forrest@0: delete (fm_l)
forrest@0: delete (landfrac)
forrest@0:
forrest@0: ;-----------------------------
forrest@0: ; read biome data: model
forrest@0:
forrest@0: biome_name_mod = "Model PFT Class"
forrest@0:
forrest@0: film_c = "class_pft_"+model_grid+".nc"
forrest@0: fm_c = addfile (dirs+film_c,"r")
forrest@0: classmod = fm_c->CLASS_PFT
forrest@0:
forrest@0: delete (fm_c)
forrest@0:
forrest@0: ; model data has 17 land-type classes
forrest@0: nclass_mod = 17
forrest@0:
forrest@0: ;---------------------------------------------------
forrest@0: ; read data: observed at stations
forrest@0:
forrest@0: station = (/"DukeFACE" \
forrest@0: ,"AspenFACE" \
forrest@0: ,"ORNL-FACE" \
forrest@0: ,"POP-EUROFACE" \
forrest@0: /)
forrest@0:
forrest@0: lat_ob = (/ 35.58, 45.40, 35.54, 42.22/)
forrest@0: lon_ob = (/-79.05, -89.37, -84.20, 11.48/)
forrest@0: lon_obx = where(lon_ob.lt.0.,lon_ob+360.,lon_ob)
forrest@0:
forrest@0: n_sta = dimsizes(station)
forrest@0: beta_4_ob = new((/n_sta/),float)
forrest@0:
forrest@0: ;###################################################
forrest@0: ; this is a hardwired value
forrest@0: beta_4_ob = 0.60
forrest@0: ;###################################################
forrest@0: ;---------------------------------------------------
forrest@0: ; get model data at station
forrest@0:
forrest@0: npp_i_4 =linint2_points(xm,ym,npp_i,True,lon_obx,lat_ob,0)
forrest@0:
forrest@0: npp_f_4 =linint2_points(xm,ym,npp_f,True,lon_obx,lat_ob,0)
forrest@0:
forrest@0: ;---------------------------------------------------
forrest@0: ;compute beta_4
forrest@0:
forrest@0: score_max = 3.
forrest@0:
forrest@0: beta_4 = new((/n_sta/),float)
forrest@0:
forrest@0: beta_4 = ((npp_f_4/npp_i_4) - 1.)/log(co2_f/co2_i)
forrest@0:
forrest@0: beta_4_avg = avg(beta_4)
forrest@0:
forrest@0: bias = sum(abs(beta_4-beta_4_ob)/(abs(beta_4)+abs(beta_4_ob)))
forrest@0: Mbeta = (1. - (bias/n_sta))*score_max
forrest@0: M_beta = sprintf("%.2f", Mbeta)
forrest@0:
forrest@0: ;=========================
forrest@0: ; for html table - station
forrest@0: ;=========================
forrest@0:
forrest@0: output_html = "table_station.html"
forrest@0:
forrest@0: ; column (not including header column)
forrest@0:
forrest@0: col_head = (/"Latitude","Longitude","CO2_i","CO2_f","NPP_i","NPP_f","Beta_model","Beta_ob"/)
forrest@0:
forrest@0: ncol = dimsizes(col_head)
forrest@0:
forrest@0: ; row (not including header row)
forrest@0: row_head = (/"DukeFACE" \
forrest@0: ,"AspenFACE" \
forrest@0: ,"ORNL-FACE" \
forrest@0: ,"POP-EUROFACE" \
forrest@0: ,"All Station" \
forrest@0: /)
forrest@0: nrow = dimsizes(row_head)
forrest@0:
forrest@0: ; arrays to be passed to table.
forrest@0: text = new ((/nrow, ncol/),string )
forrest@0:
forrest@0: do i=0,nrow-2
forrest@0: text(i,0) = sprintf("%.1f",lat_ob(i))
forrest@0: text(i,1) = sprintf("%.1f",lon_ob(i))
forrest@0: text(i,2) = sprintf("%.1f",co2_i)
forrest@0: text(i,3) = sprintf("%.1f",co2_f)
forrest@0: text(i,4) = sprintf("%.1f",npp_i_4(0,i))
forrest@0: text(i,5) = sprintf("%.1f",npp_f_4(0,i))
forrest@0: text(i,6) = sprintf("%.2f",beta_4(i))
forrest@0: text(i,7) = "-"
forrest@0: end do
forrest@0: text(nrow-1,0) = "-"
forrest@0: text(nrow-1,1) = "-"
forrest@0: text(nrow-1,2) = "-"
forrest@0: text(nrow-1,3) = "-"
forrest@0: text(nrow-1,4) = "-"
forrest@0: text(nrow-1,5) = "-"
forrest@0: text(nrow-1,6) = sprintf("%.2f",beta_4_avg)
forrest@0: text(nrow-1,7) = sprintf("%.2f",avg(beta_4_ob))
forrest@0:
forrest@0: ;-----------
forrest@0: ; html table
forrest@0: ;-----------
forrest@0:
forrest@0: header_text = "Beta Factor: Model "+model_name+"
"
forrest@0:
forrest@0: header = (/"" \
forrest@0: ,"" \
forrest@0: ,"CLAMP metrics" \
forrest@0: ,"" \
forrest@0: ,header_text \
forrest@0: /)
forrest@0: footer = ""
forrest@0:
forrest@0: table_header = (/ \
forrest@0: "" \
forrest@0: ,"" \
forrest@0: ," | Station | " \
forrest@0: ," "+col_head(0)+" | " \
forrest@0: ," "+col_head(1)+" | " \
forrest@0: ," "+col_head(2)+" | " \
forrest@0: ," "+col_head(3)+" | " \
forrest@0: ," "+col_head(4)+"
("+unit+") | " \
forrest@0: ," "+col_head(5)+"
("+unit+") | " \
forrest@0: ," "+col_head(6)+" | " \
forrest@0: ," "+col_head(7)+" | " \
forrest@0: ,"
" \
forrest@0: /)
forrest@0: table_footer = "
"
forrest@0: row_header = ""
forrest@0: row_footer = "
"
forrest@0:
forrest@0: lines = new(50000,string)
forrest@0: nline = 0
forrest@0:
forrest@0: set_line(lines,nline,header)
forrest@0: set_line(lines,nline,table_header)
forrest@0: ;-----------------------------------------------
forrest@0: ;row of table
forrest@0:
forrest@0: do n = 0,nrow-1
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: txt1 = row_head(n)
forrest@0: txt2 = text(n,0)
forrest@0: txt3 = text(n,1)
forrest@0: txt4 = text(n,2)
forrest@0: txt5 = text(n,3)
forrest@0: txt6 = text(n,4)
forrest@0: txt7 = text(n,5)
forrest@0: txt8 = text(n,6)
forrest@0: txt9 = text(n,7)
forrest@0:
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0: set_line(lines,nline,""+txt4+" | ")
forrest@0: set_line(lines,nline,""+txt5+" | ")
forrest@0: set_line(lines,nline,""+txt6+" | ")
forrest@0: set_line(lines,nline,""+txt7+" | ")
forrest@0: set_line(lines,nline,""+txt8+" | ")
forrest@0: set_line(lines,nline,""+txt9+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: end do
forrest@0: ;-----------------------------------------------
forrest@0: set_line(lines,nline,table_footer)
forrest@0: set_line(lines,nline,footer)
forrest@0:
forrest@0: ; Now write to an HTML file.
forrest@0: idx = ind(.not.ismissing(lines))
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: asciiwrite(output_html,lines(idx))
forrest@0: else
forrest@0: print ("error?")
forrest@0: end if
forrest@0:
forrest@0: delete (col_head)
forrest@0: delete (row_head)
forrest@0: delete (text)
forrest@0: delete (table_header)
forrest@0: delete (idx)
forrest@0:
forrest@0: ;********************************************************************
forrest@0: ; use land-type class to bin the data in equally spaced ranges
forrest@0: ;********************************************************************
forrest@0:
forrest@0: ; using model biome class
forrest@0: nclass = nclass_mod
forrest@0:
forrest@0: range = fspan(0,nclass,nclass+1)
forrest@0:
forrest@0: ; Use this range information to grab all the values in a
forrest@0: ; particular range, and then take an average.
forrest@0:
forrest@0: nx = dimsizes(range) - 1
forrest@0:
forrest@0: ;==============================
forrest@0: ; put data into bins
forrest@0: ;==============================
forrest@0:
forrest@0: ; for model data and observed
forrest@0: data_n = 2
forrest@0:
forrest@0: ; using model biome class
forrest@0:
forrest@0: base = ndtooned(classmod)
forrest@0:
forrest@0: ; output
forrest@0:
forrest@0: yvalues = new((/data_n,nx/),float)
forrest@0: count = new((/data_n,nx/),float)
forrest@0:
forrest@0: ; Loop through each range, using base
forrest@0:
forrest@0: do i=0,nx-1
forrest@0:
forrest@0: if (i.ne.(nx-1)) then
forrest@0: idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
forrest@0: else
forrest@0: idx = ind(base.ge.range(i))
forrest@0: end if
forrest@0:
forrest@0: ; loop through each dataset
forrest@0:
forrest@0: do n = 0,data_n-1
forrest@0:
forrest@0: if (n .eq. 0) then
forrest@0: data = ndtooned(npp_i)
forrest@0: end if
forrest@0:
forrest@0: if (n .eq. 1) then
forrest@0: data = ndtooned(npp_f)
forrest@0: end if
forrest@0:
forrest@0: ; Calculate average
forrest@0:
forrest@0: if (.not.any(ismissing(idx))) then
forrest@0: yvalues(n,i) = avg(data(idx))
forrest@0: count(n,i) = dimsizes(idx)
forrest@0: else
forrest@0: yvalues(n,i) = yvalues@_FillValue
forrest@0: count(n,i) = 0
forrest@0: end if
forrest@0:
forrest@0: ;#############################################################
forrest@0: ; using model biome class:
forrest@0: ;
forrest@0: ; set the following 4 classes to _FillValue:
forrest@0: ; (3)Needleleaf Deciduous Boreal Tree,
forrest@0: ; (8)Broadleaf Deciduous Boreal Tree,
forrest@0: ; (9)Broadleaf Evergreen Shrub,
forrest@0: ; (16)Wheat
forrest@0:
forrest@0: if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
forrest@0: yvalues(n,i) = yvalues@_FillValue
forrest@0: count(n,i) = 0
forrest@0: end if
forrest@0: ;#############################################################
forrest@0:
forrest@0: delete(data)
forrest@0: end do ; n-loop
forrest@0:
forrest@0: delete(idx)
forrest@0: end do ; i-loop
forrest@0:
forrest@0: delete (base)
forrest@0: delete (npp_i)
forrest@0: delete (npp_f)
forrest@0:
forrest@0: ;============================
forrest@0: ;compute beta
forrest@0: ;============================
forrest@0:
forrest@0: u = yvalues(0,:)
forrest@0: v = yvalues(1,:)
forrest@0: u_count = count(0,:)
forrest@0: v_count = count(1,:)
forrest@0:
forrest@0: good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0:
forrest@0: uu = u(good)
forrest@0: vv = v(good)
forrest@0: uu_count = u_count(good)
forrest@0: vv_count = v_count(good)
forrest@0:
forrest@0: n_biome = dimsizes(uu)
forrest@0: beta_biome = new((/n_biome/),float)
forrest@0:
forrest@0: beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)
forrest@0:
forrest@0: beta_biome_avg = (sum(vv*vv_count)/sum(uu*uu_count) - 1.)/log(co2_f/co2_i)
forrest@0:
forrest@0: ;===========================
forrest@0: ; for html table - biome
forrest@0: ;===========================
forrest@0:
forrest@0: output_html = "table_biome.html"
forrest@0:
forrest@0: ; column (not including header column)
forrest@0:
forrest@0: col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta_model"/)
forrest@0:
forrest@0: ncol = dimsizes(col_head)
forrest@0:
forrest@0: ; row (not including header row)
forrest@0:
forrest@0: ;----------------------------------------------------
forrest@0: ; using model biome class:
forrest@0: ;
forrest@0: row_head = (/"Not Vegetated" \
forrest@0: ,"Needleleaf Evergreen Temperate Tree" \
forrest@0: ,"Needleleaf Evergreen Boreal Tree" \
forrest@0: ; ,"Needleleaf Deciduous Boreal Tree" \
forrest@0: ,"Broadleaf Evergreen Tropical Tree" \
forrest@0: ,"Broadleaf Evergreen Temperate Tree" \
forrest@0: ,"Broadleaf Deciduous Tropical Tree" \
forrest@0: ,"Broadleaf Deciduous Temperate Tree" \
forrest@0: ; ,"Broadleaf Deciduous Boreal Tree" \
forrest@0: ; ,"Broadleaf Evergreen Shrub" \
forrest@0: ,"Broadleaf Deciduous Temperate Shrub" \
forrest@0: ,"Broadleaf Deciduous Boreal Shrub" \
forrest@0: ,"C3 Arctic Grass" \
forrest@0: ,"C3 Non-Arctic Grass" \
forrest@0: ,"C4 Grass" \
forrest@0: ,"Corn" \
forrest@0: ; ,"Wheat" \
forrest@0: ,"All Biome" \
forrest@0: /)
forrest@0:
forrest@0: nrow = dimsizes(row_head)
forrest@0:
forrest@0: ; arrays to be passed to table.
forrest@0: text = new ((/nrow, ncol/),string )
forrest@0:
forrest@0: do i=0,nrow-2
forrest@0: text(i,0) = sprintf("%.1f",co2_i)
forrest@0: text(i,1) = sprintf("%.1f",co2_f)
forrest@0: text(i,2) = sprintf("%.1f",uu(i))
forrest@0: text(i,3) = sprintf("%.1f",vv(i))
forrest@0: text(i,4) = sprintf("%.2f",beta_biome(i))
forrest@0: end do
forrest@0: text(nrow-1,0) = "-"
forrest@0: text(nrow-1,1) = "-"
forrest@0: text(nrow-1,2) = "-"
forrest@0: text(nrow-1,3) = "-"
forrest@0: text(nrow-1,4) = sprintf("%.2f",beta_biome_avg)
forrest@0:
forrest@0: ;**************************************************
forrest@0: ; html table
forrest@0: ;**************************************************
forrest@0:
forrest@0: header_text = "Beta Factor: Model "+model_name+"
"
forrest@0:
forrest@0: header = (/"" \
forrest@0: ,"" \
forrest@0: ,"CLAMP metrics" \
forrest@0: ,"" \
forrest@0: ,header_text \
forrest@0: /)
forrest@0: footer = ""
forrest@0:
forrest@0: table_header = (/ \
forrest@0: "" \
forrest@0: ,"" \
forrest@0: ," | Biome Class | " \
forrest@0: ," "+col_head(0)+" | " \
forrest@0: ," "+col_head(1)+" | " \
forrest@0: ," "+col_head(2)+" | " \
forrest@0: ," "+col_head(3)+" | " \
forrest@0: ," "+col_head(4)+" | " \
forrest@0: ,"
" \
forrest@0: /)
forrest@0: table_footer = "
"
forrest@0: row_header = ""
forrest@0: row_footer = "
"
forrest@0:
forrest@0: lines = new(50000,string)
forrest@0: nline = 0
forrest@0:
forrest@0: set_line(lines,nline,header)
forrest@0: set_line(lines,nline,table_header)
forrest@0: ;-----------------------------------------------
forrest@0: ;row of table
forrest@0:
forrest@0: do n = 0,nrow-1
forrest@0: set_line(lines,nline,row_header)
forrest@0:
forrest@0: txt1 = row_head(n)
forrest@0: txt2 = text(n,0)
forrest@0: txt3 = text(n,1)
forrest@0: txt4 = text(n,2)
forrest@0: txt5 = text(n,3)
forrest@0: txt6 = text(n,4)
forrest@0:
forrest@0: set_line(lines,nline,""+txt1+" | ")
forrest@0: set_line(lines,nline,""+txt2+" | ")
forrest@0: set_line(lines,nline,""+txt3+" | ")
forrest@0: set_line(lines,nline,""+txt4+" | ")
forrest@0: set_line(lines,nline,""+txt5+" | ")
forrest@0: set_line(lines,nline,""+txt6+" | ")
forrest@0:
forrest@0: set_line(lines,nline,row_footer)
forrest@0: end do
forrest@0: ;-----------------------------------------------
forrest@0: set_line(lines,nline,table_footer)
forrest@0: set_line(lines,nline,footer)
forrest@0:
forrest@0: ; Now write to an HTML file
forrest@0:
forrest@0: idx = ind(.not.ismissing(lines))
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: asciiwrite(output_html,lines(idx))
forrest@0: else
forrest@0: print ("error?")
forrest@0: end if
forrest@0: delete (idx)
forrest@0:
forrest@0: ;**************************************************************************************
forrest@0: ; update score
forrest@0: ;**************************************************************************************
forrest@0: if (isvar("compare")) then
forrest@0: system("sed -e '1,/M_beta/s/M_beta/"+M_beta+"/' "+html_name2+" > "+html_new2+";"+ \
forrest@0: "mv -f "+html_new2+" "+html_name2)
forrest@0: end if
forrest@0:
forrest@0: system("sed s#M_beta#"+M_beta+"# "+html_name+" > "+html_new+";"+ \
forrest@0: "mv -f "+html_new+" "+html_name)
forrest@0:
forrest@0: ;***************************************************************************
forrest@0: ; get total score and write to file
forrest@0: ;***************************************************************************
forrest@0: M_total = Mbeta
forrest@0:
forrest@0: asciiwrite("M_save.beta", M_total)
forrest@0:
forrest@0: delete (M_total)
forrest@0:
forrest@0: ;***************************************************************************
forrest@0: ; output plot and html
forrest@0: ;***************************************************************************
forrest@0: output_dir = model_name+"/beta"
forrest@0:
forrest@0: system("mv *.html " + output_dir)
forrest@0: ;***************************************************************************
forrest@0:
forrest@0: end
forrest@0: