1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/all/08.turnover.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,351 @@
1.4 +;********************************************************
1.5 +; hardwire: flux = flux/1200. (for casa only)
1.6 +;**************************************************************
1.7 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.8 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.10 +;**************************************************************
1.11 +procedure set_line(lines:string,nline:integer,newlines:string)
1.12 +begin
1.13 +; add line to ascci/html file
1.14 +
1.15 + nnewlines = dimsizes(newlines)
1.16 + if(nline+nnewlines-1.ge.dimsizes(lines))
1.17 + print("set_line: bad index, not setting anything.")
1.18 + return
1.19 + end if
1.20 + lines(nline:nline+nnewlines-1) = newlines
1.21 +; print ("lines = " + lines(nline:nline+nnewlines-1))
1.22 + nline = nline + nnewlines
1.23 + return
1.24 +end
1.25 +;**************************************************************
1.26 +; Main code.
1.27 +begin
1.28 +
1.29 + plot_type = "ps"
1.30 + plot_type_new = "png"
1.31 +
1.32 +;------------------------------------------------------
1.33 +; edit table.html of current model for movel1_vs_model2
1.34 +
1.35 + if (isvar("compare")) then
1.36 + html_name2 = compare+"/table.html"
1.37 + html_new2 = html_name2 +".new"
1.38 + end if
1.39 +
1.40 +;------------------------------------------------------
1.41 +; edit table.html for current model
1.42 +
1.43 + html_name = model_name+"/table.html"
1.44 + html_new = html_name +".new"
1.45 +
1.46 +;---------------------------------------------------------------
1.47 +;components
1.48 +
1.49 + component = (/"Leaf","Wood","Fine_Root","Litter","Coarse_Woody_Debris","Soil"/)
1.50 + n_comp = dimsizes(component)
1.51 +
1.52 + field_pool = (/"LEAFC","WOODC","FROOTC","LITTERC","CWDC","SOILC"/)
1.53 + field_flux = (/"LEAFC_ALLOC","WOODC_ALLOC","FROOTC_ALLOC","LITTERC_LOSS","CWDC_LOSS","SOILC_HR"/)
1.54 +
1.55 +;--------------------------------------------------
1.56 +; get landfrac data
1.57 +
1.58 + film_l = "lnd_"+ model_grid +".nc"
1.59 + fm_l = addfile (dirs+film_l,"r")
1.60 + landfrac = fm_l->landfrac
1.61 +
1.62 + delete (fm_l)
1.63 +;---------------------------------------------------
1.64 +; read biome data: model
1.65 +
1.66 + biome_name_mod = "Model PFT Class"
1.67 +
1.68 + film_c = "class_pft_"+ model_grid +".nc"
1.69 + fm_c = addfile (dirs+film_c,"r")
1.70 + classmod = fm_c->CLASS_PFT
1.71 +
1.72 + delete (fm_c)
1.73 +
1.74 +; model data has 17 land-type classes
1.75 + nclass_mod = 17
1.76 +
1.77 +;********************************************************************
1.78 +; use land-type class to bin the data in equally spaced ranges
1.79 +;********************************************************************
1.80 +
1.81 +; using model biome class
1.82 + nclass = nclass_mod
1.83 +
1.84 + range = fspan(0,nclass,nclass+1)
1.85 +
1.86 +; Use this range information to grab all the values in a
1.87 +; particular range, and then take an average.
1.88 +
1.89 + nx = dimsizes(range) - 1
1.90 +
1.91 +; for 2 data: pool and flux
1.92 + data_n = 2
1.93 +
1.94 +; using model biome class
1.95 +
1.96 + base = ndtooned(classmod)
1.97 +
1.98 +; output
1.99 +
1.100 + yvalues = new((/data_n,nx/),float)
1.101 + count = new((/data_n,nx/),float)
1.102 +
1.103 +;--------------------------------------------------
1.104 +; read model data, each component:
1.105 +
1.106 + fm = addfile (dirm+film4,"r")
1.107 +
1.108 + do k = 0,n_comp-1
1.109 +
1.110 + pool = fm->$field_pool(k)$
1.111 + flux = fm->$field_flux(k)$
1.112 +
1.113 +; Units for these variables are:
1.114 +; pool: g C/m^2
1.115 +; flux: g C/m^2/s
1.116 +
1.117 + nsec_per_year = 60*60*24*365
1.118 +
1.119 + flux = flux * nsec_per_year
1.120 +
1.121 + unit_p = "gC/m2"
1.122 + unit_f = "gC/m2/year"
1.123 + unit_t = "year"
1.124 +
1.125 +;#############################################################
1.126 +; casa only
1.127 +; all the plant pools (leaf, wood, and fine root) and
1.128 +; coarse woody debris (cwd) and litter pools for
1.129 +; CASA need to be divided by 1200. The soil flux
1.130 +; and turnover time are fine and do not need to be adjusted.
1.131 +
1.132 + if (BGC .eq. "casa") then
1.133 + if (k .ne. n_comp-1) then
1.134 + flux = flux/1200.
1.135 + end if
1.136 + end if
1.137 +;##############################################################
1.138 +
1.139 +; take into account landfrac
1.140 +
1.141 + pool = pool * conform(pool,landfrac,(/1,2/))
1.142 + flux = flux * conform(flux,landfrac,(/1,2/))
1.143 +
1.144 +; Loop through each range, using base
1.145 +
1.146 + do i=0,nx-1
1.147 +
1.148 + if (i.ne.(nx-1)) then
1.149 + idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
1.150 + else
1.151 + idx = ind(base.ge.range(i))
1.152 + end if
1.153 +
1.154 +; loop through each dataset
1.155 +
1.156 + do n = 0,data_n-1
1.157 +
1.158 + if (n .eq. 0) then
1.159 + data = ndtooned(pool)
1.160 + end if
1.161 +
1.162 + if (n .eq. 1) then
1.163 + data = ndtooned(flux)
1.164 + end if
1.165 +
1.166 +; Calculate average
1.167 +
1.168 + if (.not.any(ismissing(idx))) then
1.169 + yvalues(n,i) = avg(data(idx))
1.170 + count(n,i) = dimsizes(idx)
1.171 + else
1.172 + yvalues(n,i) = yvalues@_FillValue
1.173 + count(n,i) = 0
1.174 + end if
1.175 +
1.176 +;#############################################################
1.177 +; using model biome class:
1.178 +;
1.179 +; set the following 4 classes to _FillValue:
1.180 +; (3)Needleleaf Deciduous Boreal Tree,
1.181 +; (8)Broadleaf Deciduous Boreal Tree,
1.182 +; (9)Broadleaf Evergreen Shrub,
1.183 +; (16)Wheat
1.184 +
1.185 + if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
1.186 + yvalues(n,i) = yvalues@_FillValue
1.187 + count(n,i) = 0
1.188 + end if
1.189 +;#############################################################
1.190 +
1.191 + delete(data)
1.192 + end do ; n-loop
1.193 +
1.194 + delete(idx)
1.195 + end do ; i-loop
1.196 +
1.197 + delete (pool)
1.198 + delete (flux)
1.199 +
1.200 +;============================
1.201 +;compute turnover time
1.202 +;============================
1.203 +
1.204 + u = yvalues(0,:)
1.205 + v = yvalues(1,:)
1.206 + u_count = count(0,:)
1.207 + v_count = count(1,:)
1.208 +
1.209 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.210 +
1.211 + uu = u(good)
1.212 + vv = v(good)
1.213 + uu_count = u_count(good)
1.214 + vv_count = v_count(good)
1.215 +
1.216 + n_biome = dimsizes(uu)
1.217 + t_biome = new((/n_biome/),float)
1.218 +
1.219 + t_biome = uu/vv
1.220 +
1.221 + t_biome_avg = sum(uu*uu_count)/sum(vv*vv_count)
1.222 +
1.223 +;===========================
1.224 +; for html table - biome
1.225 +;===========================
1.226 +
1.227 + output_html = "table_"+component(k)+".html"
1.228 +
1.229 +; column (not including header column)
1.230 +
1.231 + col_head = (/component(k)+" Flux",component(k)+" Pool",component(k)+" Turnover Time"/)
1.232 +
1.233 + ncol = dimsizes(col_head)
1.234 +
1.235 +; row (not including header row)
1.236 +
1.237 +;----------------------------------------------------
1.238 +; using model biome class:
1.239 + row_head = (/"Not Vegetated" \
1.240 + ,"Needleleaf Evergreen Temperate Tree" \
1.241 + ,"Needleleaf Evergreen Boreal Tree" \
1.242 +; ,"Needleleaf Deciduous Boreal Tree" \
1.243 + ,"Broadleaf Evergreen Tropical Tree" \
1.244 + ,"Broadleaf Evergreen Temperate Tree" \
1.245 + ,"Broadleaf Deciduous Tropical Tree" \
1.246 + ,"Broadleaf Deciduous Temperate Tree" \
1.247 +; ,"Broadleaf Deciduous Boreal Tree" \
1.248 +; ,"Broadleaf Evergreen Shrub" \
1.249 + ,"Broadleaf Deciduous Temperate Shrub" \
1.250 + ,"Broadleaf Deciduous Boreal Shrub" \
1.251 + ,"C3 Arctic Grass" \
1.252 + ,"C3 Non-Arctic Grass" \
1.253 + ,"C4 Grass" \
1.254 + ,"Corn" \
1.255 +; ,"Wheat" \
1.256 + ,"All Biome" \
1.257 + /)
1.258 + nrow = dimsizes(row_head)
1.259 +
1.260 +; arrays to be passed to table.
1.261 + text = new ((/nrow, ncol/),string )
1.262 +
1.263 + do i=0,nrow-2
1.264 + text(i,0) = sprintf("%.1f",vv(i))
1.265 + text(i,1) = sprintf("%.1f",uu(i))
1.266 + text(i,2) = sprintf("%.2f",t_biome(i))
1.267 + end do
1.268 + text(nrow-1,0) = "-"
1.269 + text(nrow-1,1) = "-"
1.270 + text(nrow-1,2) = sprintf("%.2f",t_biome_avg)
1.271 +
1.272 +;**************************************************
1.273 +; html table
1.274 +;**************************************************
1.275 +
1.276 + header_text = "<H1>"+component(k)+" Turnover Time: Model "+model_name+"</H1>"
1.277 +
1.278 + header = (/"<HTML>" \
1.279 + ,"<HEAD>" \
1.280 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.281 + ,"</HEAD>" \
1.282 + ,header_text \
1.283 + /)
1.284 + footer = "</HTML>"
1.285 +
1.286 + table_header = (/ \
1.287 + "<table border=1 cellspacing=0 cellpadding=3 width=60%>" \
1.288 + ,"<tr>" \
1.289 + ," <th bgcolor=DDDDDD >Biome Class</th>" \
1.290 + ," <th bgcolor=DDDDDD >"+col_head(0)+"<br>("+unit_f+")</th>" \
1.291 + ," <th bgcolor=DDDDDD >"+col_head(1)+"<br>("+unit_p+")</th>" \
1.292 + ," <th bgcolor=DDDDDD >"+col_head(2)+"<br>("+unit_t+")</th>" \
1.293 + ,"</tr>" \
1.294 + /)
1.295 + table_footer = "</table>"
1.296 + row_header = "<tr>"
1.297 + row_footer = "</tr>"
1.298 +
1.299 + lines = new(50000,string)
1.300 + nline = 0
1.301 +
1.302 + set_line(lines,nline,header)
1.303 + set_line(lines,nline,table_header)
1.304 +;-----------------------------------------------
1.305 +; row of table
1.306 +
1.307 + do n = 0,nrow-1
1.308 + set_line(lines,nline,row_header)
1.309 +
1.310 + txt1 = row_head(n)
1.311 + txt2 = text(n,0)
1.312 + txt3 = text(n,1)
1.313 + txt4 = text(n,2)
1.314 +
1.315 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.316 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.317 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.318 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.319 +
1.320 + set_line(lines,nline,row_footer)
1.321 + end do
1.322 +;-----------------------------------------------
1.323 + set_line(lines,nline,table_footer)
1.324 + set_line(lines,nline,footer)
1.325 +
1.326 +; Now write to an HTML file
1.327 +
1.328 + idx = ind(.not.ismissing(lines))
1.329 + if(.not.any(ismissing(idx))) then
1.330 + asciiwrite(output_html,lines(idx))
1.331 + else
1.332 + print ("error?")
1.333 + end if
1.334 +
1.335 + delete (idx)
1.336 +
1.337 + delete (good)
1.338 + delete (t_biome)
1.339 + delete (text)
1.340 +
1.341 + end do ; k-loop
1.342 +
1.343 + delete (fm)
1.344 +
1.345 +;***************************************************************************
1.346 +; output plot and html
1.347 +;***************************************************************************
1.348 + output_dir = model_name+"/turnover"
1.349 +
1.350 + system("mv *.html " + output_dir)
1.351 +;******************************
1.352 +
1.353 +end
1.354 +