forrest@0: ;********************************************************
forrest@0: ; required command line input parameters:
forrest@0: ; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
forrest@0: ;
forrest@0: ; histogram normalized by rain and compute correleration
forrest@0: ;**************************************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
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: nclass = 20
forrest@0:
forrest@0: plot_type = "ps"
forrest@0: plot_type_new = "png"
forrest@0: ;************************************************
forrest@0: ; read data: model
forrest@0: ;************************************************
forrest@0: co2_i = 283.1878
forrest@0: co2_f = 364.1252
forrest@0:
forrest@0: model_grid = "T42"
forrest@0:
forrest@0: ;model_name_i = "i01.07cn"
forrest@0: ;model_name_f = "i01.10cn"
forrest@0:
forrest@0: model_name_i = "i01.07casa"
forrest@0: model_name_f = "i01.10casa"
forrest@0:
forrest@0: model_name = model_name_f
forrest@0:
forrest@0: dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0: film_i = model_name_i + "_1990-2004_ANN_climo.nc"
forrest@0: film_f = model_name_f + "_1990-2004_ANN_climo.nc"
forrest@0:
forrest@0: fm_i = addfile (dirm+film_i,"r")
forrest@0: fm_f = addfile (dirm+film_f,"r")
forrest@0:
forrest@0: npp_i = fm_i->NPP
forrest@0: npp_f = fm_f->NPP
forrest@0:
forrest@0: ;************************************************
forrest@0: ; read data: observed
forrest@0: ;************************************************
forrest@0:
forrest@0: ob_name = "MODIS MOD 15A2 2000-2005"
forrest@0:
forrest@0: diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
forrest@0: filo = "land_class_"+model_grid+".nc"
forrest@0:
forrest@0: fo = addfile(diro+filo,"r")
forrest@0:
forrest@0: classob = tofloat(fo->LAND_CLASS)
forrest@0:
forrest@0: ;*******************************************************************
forrest@0: ; Calculate "nice" bins for binning the data in equally spaced ranges
forrest@0: ;********************************************************************
forrest@0: nclassn = nclass + 1
forrest@0: range = fspan(0,nclassn-1,nclassn)
forrest@0: ; print (range)
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: nr = dimsizes(range)
forrest@0: nx = nr-1
forrest@0: xvalues = new((/2,nx/),float)
forrest@0: xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
forrest@0: dx = xvalues(0,1) - xvalues(0,0) ; range width
forrest@0: dx4 = dx/4 ; 1/4 of the range
forrest@0: xvalues(1,:) = xvalues(0,:) - dx/5.
forrest@0:
forrest@0: ; get data
forrest@0:
forrest@0: DATA11_1D = ndtooned(classob)
forrest@0: DATA12_1D = ndtooned(npp_i)
forrest@0: DATA22_1D = ndtooned(npp_f)
forrest@0:
forrest@0: yvalues = new((/2,nx/),float)
forrest@0: mn_yvalues = new((/2,nx/),float)
forrest@0: mx_yvalues = new((/2,nx/),float)
forrest@0:
forrest@0: do nd=0,1
forrest@0:
forrest@0: ; See if we are doing model or observational data.
forrest@0:
forrest@0: if(nd.eq.0) then
forrest@0: data_ob = DATA11_1D
forrest@0: data_mod = DATA12_1D
forrest@0: else
forrest@0: data_ob = DATA11_1D
forrest@0: data_mod = DATA22_1D
forrest@0: end if
forrest@0:
forrest@0: ; Loop through each range and check for values.
forrest@0:
forrest@0: do i=0,nr-2
forrest@0: if (i.ne.(nr-2)) then
forrest@0: ; print("")
forrest@0: ; print("In range ["+range(i)+","+range(i+1)+")")
forrest@0: idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1)))
forrest@0: else
forrest@0: ; print("")
forrest@0: ; print("In range ["+range(i)+",)")
forrest@0: idx = ind(data_ob.ge.range(i))
forrest@0: end if
forrest@0:
forrest@0: ; Calculate average, and get min and max.
forrest@0:
forrest@0: if(.not.any(ismissing(idx))) then
forrest@0: yvalues(nd,i) = avg(data_mod(idx))
forrest@0: mn_yvalues(nd,i) = min(data_mod(idx))
forrest@0: mx_yvalues(nd,i) = max(data_mod(idx))
forrest@0: count = dimsizes(idx)
forrest@0: else
forrest@0: count = 0
forrest@0: yvalues(nd,i) = yvalues@_FillValue
forrest@0: mn_yvalues(nd,i) = yvalues@_FillValue
forrest@0: mx_yvalues(nd,i) = yvalues@_FillValue
forrest@0: end if
forrest@0:
forrest@0: ; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
forrest@0: ; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
forrest@0:
forrest@0: ; Clean up for next time in loop.
forrest@0:
forrest@0: delete(idx)
forrest@0: end do
forrest@0: delete(data_ob)
forrest@0: delete(data_mod)
forrest@0: end do
forrest@0: ;============================
forrest@0: ;compute beta
forrest@0: ;============================
forrest@0:
forrest@0: nsec_per_year = 60*60*24*365
forrest@0:
forrest@0: u = yvalues(0,:)
forrest@0: v = yvalues(1,:)
forrest@0:
forrest@0: good = ind(.not.ismissing(u) .and. .not.ismissing(v))
forrest@0: uu = u(good)* nsec_per_year
forrest@0: vv = v(good)* nsec_per_year
forrest@0:
forrest@0: n_biome = dimsizes(uu)
forrest@0:
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 = avg(beta_biome)
forrest@0:
forrest@0: print (beta_biome_avg)
forrest@0: ;*******************************************************************
forrest@0: ; for html table
forrest@0: ;*******************************************************************
forrest@0:
forrest@0: ; column (not including header column)
forrest@0:
forrest@0: col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta"/)
forrest@0:
forrest@0: ncol = dimsizes(col_head)
forrest@0:
forrest@0: ; row (not including header row)
forrest@0: row_head = (/"Water Bodies" \
forrest@0: ,"Evergreen Needleleaf Forests" \
forrest@0: ,"Evergreen Broadleaf Forests" \
forrest@0: ,"Deciduous Needleleaf Forest" \
forrest@0: ,"Deciduous Broadleaf Forests" \
forrest@0: ,"Mixed Forests" \
forrest@0: ,"Closed Bushlands" \
forrest@0: ,"Open Bushlands" \
forrest@0: ,"Woody Savannas (S. Hem.)" \
forrest@0: ,"Savannas (S. Hem.)" \
forrest@0: ,"Grasslands" \
forrest@0: ,"Permanent Wetlands" \
forrest@0: ,"Croplands" \
forrest@0: ,"Cropland/Natural Vegetation Mosaic" \
forrest@0: ,"Permanent Snow and Ice" \
forrest@0: ,"Barren or Sparsely Vegetated" \
forrest@0: ,"Woody Savannas (N. Hem.)" \
forrest@0: ,"Savannas (N. Hem.)" \
forrest@0: ,"All Biome" \
forrest@0: /)
forrest@0: nrow = dimsizes(row_head)
forrest@0:
forrest@0: ; arrays to be passed to table.
forrest@0: text4 = new ((/nrow, ncol/),string )
forrest@0:
forrest@0: do i=0,nrow-2
forrest@0: text4(i,0) = sprintf("%.2f",co2_i)
forrest@0: text4(i,1) = sprintf("%.2f",co2_f)
forrest@0: text4(i,2) = sprintf("%.2f",uu(i))
forrest@0: text4(i,3) = sprintf("%.2f",vv(i))
forrest@0: text4(i,4) = sprintf("%.2f",beta_biome(i))
forrest@0: end do
forrest@0: text4(nrow-1,0) = "-"
forrest@0: text4(nrow-1,1) = "-"
forrest@0: text4(nrow-1,2) = "-"
forrest@0: text4(nrow-1,3) = "-"
forrest@0: text4(nrow-1,4) = sprintf("%.2f",beta_biome_avg)
forrest@0:
forrest@0: ;**************************************************
forrest@0: ; html table
forrest@0: ;**************************************************
forrest@0: output_html = "table_biome.html"
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: ," CO2_i | " \
forrest@0: ," CO2_f | " \
forrest@0: ," NPP_i | " \
forrest@0: ," NPP_f | " \
forrest@0: ," Beta | " \
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 = text4(n,0)
forrest@0: txt3 = text4(n,1)
forrest@0: txt4 = text4(n,2)
forrest@0: txt5 = text4(n,3)
forrest@0: txt6 = text4(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: 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: end
forrest@0: