1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/beta/04.biome.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,290 @@
1.4 +;********************************************************
1.5 +; required command line input parameters:
1.6 +; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
1.7 +;
1.8 +; histogram normalized by rain and compute correleration
1.9 +;**************************************************************
1.10 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
1.11 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
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 + nclass = 20
1.33 +
1.34 + plot_type = "ps"
1.35 + plot_type_new = "png"
1.36 +;************************************************
1.37 +; read data: model
1.38 +;************************************************
1.39 + co2_i = 283.1878
1.40 + co2_f = 364.1252
1.41 +
1.42 + model_grid = "T42"
1.43 +
1.44 +;model_name_i = "i01.07cn"
1.45 +;model_name_f = "i01.10cn"
1.46 +
1.47 + model_name_i = "i01.07casa"
1.48 + model_name_f = "i01.10casa"
1.49 +
1.50 + model_name = model_name_f
1.51 +
1.52 + dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.53 + film_i = model_name_i + "_1990-2004_ANN_climo.nc"
1.54 + film_f = model_name_f + "_1990-2004_ANN_climo.nc"
1.55 +
1.56 + fm_i = addfile (dirm+film_i,"r")
1.57 + fm_f = addfile (dirm+film_f,"r")
1.58 +
1.59 + npp_i = fm_i->NPP
1.60 + npp_f = fm_f->NPP
1.61 +
1.62 +;************************************************
1.63 +; read data: observed
1.64 +;************************************************
1.65 +
1.66 + ob_name = "MODIS MOD 15A2 2000-2005"
1.67 +
1.68 + diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
1.69 + filo = "land_class_"+model_grid+".nc"
1.70 +
1.71 + fo = addfile(diro+filo,"r")
1.72 +
1.73 + classob = tofloat(fo->LAND_CLASS)
1.74 +
1.75 +;*******************************************************************
1.76 +; Calculate "nice" bins for binning the data in equally spaced ranges
1.77 +;********************************************************************
1.78 + nclassn = nclass + 1
1.79 + range = fspan(0,nclassn-1,nclassn)
1.80 +; print (range)
1.81 +
1.82 +; Use this range information to grab all the values in a
1.83 +; particular range, and then take an average.
1.84 +
1.85 + nr = dimsizes(range)
1.86 + nx = nr-1
1.87 + xvalues = new((/2,nx/),float)
1.88 + xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
1.89 + dx = xvalues(0,1) - xvalues(0,0) ; range width
1.90 + dx4 = dx/4 ; 1/4 of the range
1.91 + xvalues(1,:) = xvalues(0,:) - dx/5.
1.92 +
1.93 +; get data
1.94 +
1.95 + DATA11_1D = ndtooned(classob)
1.96 + DATA12_1D = ndtooned(npp_i)
1.97 + DATA22_1D = ndtooned(npp_f)
1.98 +
1.99 + yvalues = new((/2,nx/),float)
1.100 + mn_yvalues = new((/2,nx/),float)
1.101 + mx_yvalues = new((/2,nx/),float)
1.102 +
1.103 + do nd=0,1
1.104 +
1.105 +; See if we are doing model or observational data.
1.106 +
1.107 + if(nd.eq.0) then
1.108 + data_ob = DATA11_1D
1.109 + data_mod = DATA12_1D
1.110 + else
1.111 + data_ob = DATA11_1D
1.112 + data_mod = DATA22_1D
1.113 + end if
1.114 +
1.115 +; Loop through each range and check for values.
1.116 +
1.117 + do i=0,nr-2
1.118 + if (i.ne.(nr-2)) then
1.119 +; print("")
1.120 +; print("In range ["+range(i)+","+range(i+1)+")")
1.121 + idx = ind((data_ob.ge.range(i)).and.(data_ob.lt.range(i+1)))
1.122 + else
1.123 +; print("")
1.124 +; print("In range ["+range(i)+",)")
1.125 + idx = ind(data_ob.ge.range(i))
1.126 + end if
1.127 +
1.128 +; Calculate average, and get min and max.
1.129 +
1.130 + if(.not.any(ismissing(idx))) then
1.131 + yvalues(nd,i) = avg(data_mod(idx))
1.132 + mn_yvalues(nd,i) = min(data_mod(idx))
1.133 + mx_yvalues(nd,i) = max(data_mod(idx))
1.134 + count = dimsizes(idx)
1.135 + else
1.136 + count = 0
1.137 + yvalues(nd,i) = yvalues@_FillValue
1.138 + mn_yvalues(nd,i) = yvalues@_FillValue
1.139 + mx_yvalues(nd,i) = yvalues@_FillValue
1.140 + end if
1.141 +
1.142 +; print(nd + ": " + count + " points, avg = " + yvalues(nd,i))
1.143 +; print("Min/Max: " + mn_yvalues(nd,i) + "/" + mx_yvalues(nd,i))
1.144 +
1.145 +; Clean up for next time in loop.
1.146 +
1.147 + delete(idx)
1.148 + end do
1.149 + delete(data_ob)
1.150 + delete(data_mod)
1.151 + end do
1.152 +;============================
1.153 +;compute beta
1.154 +;============================
1.155 +
1.156 + nsec_per_year = 60*60*24*365
1.157 +
1.158 + u = yvalues(0,:)
1.159 + v = yvalues(1,:)
1.160 +
1.161 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.162 + uu = u(good)* nsec_per_year
1.163 + vv = v(good)* nsec_per_year
1.164 +
1.165 + n_biome = dimsizes(uu)
1.166 +
1.167 + beta_biome = new((/n_biome/),float)
1.168 +
1.169 + beta_biome = ((vv/uu) - 1.)/log(co2_f/co2_i)
1.170 +
1.171 + beta_biome_avg = avg(beta_biome)
1.172 +
1.173 + print (beta_biome_avg)
1.174 +;*******************************************************************
1.175 +; for html table
1.176 +;*******************************************************************
1.177 +
1.178 +; column (not including header column)
1.179 +
1.180 + col_head = (/"CO2_i","CO2_f","NPP_i","NPP_f","Beta"/)
1.181 +
1.182 + ncol = dimsizes(col_head)
1.183 +
1.184 +; row (not including header row)
1.185 + row_head = (/"Water Bodies" \
1.186 + ,"Evergreen Needleleaf Forests" \
1.187 + ,"Evergreen Broadleaf Forests" \
1.188 + ,"Deciduous Needleleaf Forest" \
1.189 + ,"Deciduous Broadleaf Forests" \
1.190 + ,"Mixed Forests" \
1.191 + ,"Closed Bushlands" \
1.192 + ,"Open Bushlands" \
1.193 + ,"Woody Savannas (S. Hem.)" \
1.194 + ,"Savannas (S. Hem.)" \
1.195 + ,"Grasslands" \
1.196 + ,"Permanent Wetlands" \
1.197 + ,"Croplands" \
1.198 + ,"Cropland/Natural Vegetation Mosaic" \
1.199 + ,"Permanent Snow and Ice" \
1.200 + ,"Barren or Sparsely Vegetated" \
1.201 + ,"Woody Savannas (N. Hem.)" \
1.202 + ,"Savannas (N. Hem.)" \
1.203 + ,"All Biome" \
1.204 + /)
1.205 + nrow = dimsizes(row_head)
1.206 +
1.207 +; arrays to be passed to table.
1.208 + text4 = new ((/nrow, ncol/),string )
1.209 +
1.210 + do i=0,nrow-2
1.211 + text4(i,0) = sprintf("%.2f",co2_i)
1.212 + text4(i,1) = sprintf("%.2f",co2_f)
1.213 + text4(i,2) = sprintf("%.2f",uu(i))
1.214 + text4(i,3) = sprintf("%.2f",vv(i))
1.215 + text4(i,4) = sprintf("%.2f",beta_biome(i))
1.216 + end do
1.217 + text4(nrow-1,0) = "-"
1.218 + text4(nrow-1,1) = "-"
1.219 + text4(nrow-1,2) = "-"
1.220 + text4(nrow-1,3) = "-"
1.221 + text4(nrow-1,4) = sprintf("%.2f",beta_biome_avg)
1.222 +
1.223 +;**************************************************
1.224 +; html table
1.225 +;**************************************************
1.226 + output_html = "table_biome.html"
1.227 +
1.228 + header_text = "<H1>Beta Factor: Model "+model_name+"</H1>"
1.229 +
1.230 + header = (/"<HTML>" \
1.231 + ,"<HEAD>" \
1.232 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.233 + ,"</HEAD>" \
1.234 + ,header_text \
1.235 + /)
1.236 + footer = "</HTML>"
1.237 +
1.238 + table_header = (/ \
1.239 + "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
1.240 + ,"<tr>" \
1.241 + ," <th bgcolor=DDDDDD >Biome Class</th>" \
1.242 + ," <th bgcolor=DDDDDD >CO2_i</th>" \
1.243 + ," <th bgcolor=DDDDDD >CO2_f</th>" \
1.244 + ," <th bgcolor=DDDDDD >NPP_i</th>" \
1.245 + ," <th bgcolor=DDDDDD >NPP_f</th>" \
1.246 + ," <th bgcolor=DDDDDD >Beta</th>" \
1.247 + ,"</tr>" \
1.248 + /)
1.249 + table_footer = "</table>"
1.250 + row_header = "<tr>"
1.251 + row_footer = "</tr>"
1.252 +
1.253 + lines = new(50000,string)
1.254 + nline = 0
1.255 +
1.256 + set_line(lines,nline,header)
1.257 + set_line(lines,nline,table_header)
1.258 +;-----------------------------------------------
1.259 +;row of table
1.260 +
1.261 + do n = 0,nrow-1
1.262 + set_line(lines,nline,row_header)
1.263 +
1.264 + txt1 = row_head(n)
1.265 + txt2 = text4(n,0)
1.266 + txt3 = text4(n,1)
1.267 + txt4 = text4(n,2)
1.268 + txt5 = text4(n,3)
1.269 + txt6 = text4(n,4)
1.270 +
1.271 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.272 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.273 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.274 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.275 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.276 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.277 +
1.278 + set_line(lines,nline,row_footer)
1.279 + end do
1.280 +;-----------------------------------------------
1.281 + set_line(lines,nline,table_footer)
1.282 + set_line(lines,nline,footer)
1.283 +
1.284 +; Now write to an HTML file.
1.285 + idx = ind(.not.ismissing(lines))
1.286 + if(.not.any(ismissing(idx))) then
1.287 + asciiwrite(output_html,lines(idx))
1.288 + else
1.289 + print ("error?")
1.290 + end if
1.291 +
1.292 +end
1.293 +