1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/all/02.lai.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,879 @@
1.4 +;**************************************************************
1.5 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
1.6 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
1.7 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
1.8 +;**************************************************************
1.9 +procedure set_line(lines:string,nline:integer,newlines:string)
1.10 +begin
1.11 +; add line to ascci/html file
1.12 +
1.13 + nnewlines = dimsizes(newlines)
1.14 + if(nline+nnewlines-1.ge.dimsizes(lines))
1.15 + print("set_line: bad index, not setting anything.")
1.16 + return
1.17 + end if
1.18 + lines(nline:nline+nnewlines-1) = newlines
1.19 +; print ("lines = " + lines(nline:nline+nnewlines-1))
1.20 + nline = nline + nnewlines
1.21 + return
1.22 +end
1.23 +;**************************************************************
1.24 +; Main code.
1.25 +begin
1.26 +
1.27 + plot_type = "ps"
1.28 + plot_type_new = "png"
1.29 +
1.30 +;-----------------------------------------------------
1.31 +; edit table.html of current model for movel1_vs_model2
1.32 +
1.33 + if (isvar("compare")) then
1.34 + html_name2 = compare+"/table.html"
1.35 + html_new2 = html_name2 +".new"
1.36 + end if
1.37 +
1.38 +;------------------------------------------------------
1.39 +; edit table.html for current model
1.40 +
1.41 + html_name = model_name+"/table.html"
1.42 + html_new = html_name +".new"
1.43 +
1.44 +;------------------------------------------------------
1.45 +; read data: model
1.46 +
1.47 + fm = addfile(dirm+film10,"r")
1.48 + laimod = fm->TLAI
1.49 +
1.50 + delete (fm)
1.51 +
1.52 + dsizes = dimsizes(laimod)
1.53 + ntime = dsizes(0)
1.54 + nlat = dsizes(1)
1.55 + nlon = dsizes(2)
1.56 +
1.57 +;-----------------------------------
1.58 +; get landfrac data
1.59 +
1.60 + film_l = "lnd_"+ model_grid + ".nc"
1.61 + fm_l = addfile (dirs+film_l,"r")
1.62 + landfrac = fm_l->landfrac
1.63 +
1.64 + delete (fm_l)
1.65 +;----------------------------------
1.66 +; read biome data: model
1.67 +
1.68 + biome_name_mod = "Model PFT Class"
1.69 +
1.70 + film_c = "class_pft_"+model_grid+".nc"
1.71 + fm_c = addfile(dirs+film_c,"r")
1.72 + classmod = fm_c->CLASS_PFT
1.73 +
1.74 + delete (fm_c)
1.75 +
1.76 +; model data has 17 land-type classes
1.77 + nclass_mod = 17
1.78 +
1.79 +;----------------------------------------------------------
1.80 +; read data: ob
1.81 +
1.82 +;----------------------------------
1.83 +; read biome data: observed
1.84 +
1.85 + biome_name_ob = "MODIS LandCover"
1.86 +
1.87 + dir_c = diro + "lai/"
1.88 + filo_c = "land_class_"+model_grid+".nc"
1.89 + fo = addfile(dir_c+filo_c,"r")
1.90 + classob = tofloat(fo->LAND_CLASS)
1.91 +
1.92 + delete (fo)
1.93 +
1.94 +; input observed data has 20 land-type classes
1.95 + nclass_ob = 20
1.96 +
1.97 +;---------------------------------
1.98 +; read lai data: observed
1.99 +
1.100 + ;ob_name = "MODIS MOD 15A2 2000-2005"
1.101 + ob_name = "MODIS MOD 15A2 2000-2004"
1.102 +
1.103 + dir_l = diro + "lai/"
1.104 + ;filo = "LAI_2000-2005_MONS_"+model_grid+".nc"
1.105 + filo = "LAI_2000-2004_MONS_"+model_grid+".nc"
1.106 + fo = addfile(dir_l+filo,"r")
1.107 + laiob = fo->LAI
1.108 +
1.109 + delete (fo)
1.110 +
1.111 +;-------------------------------------------------
1.112 +; take into account landfrac
1.113 +
1.114 + laimod = laimod * conform(laimod,landfrac,(/1,2/))
1.115 + laiob = laiob * conform(laiob ,landfrac,(/1,2/))
1.116 +
1.117 + delete (landfrac)
1.118 +
1.119 +;************************************************
1.120 +; global res
1.121 +;************************************************
1.122 + resg = True ; Use plot options
1.123 + resg@cnFillOn = True ; Turn on color fill
1.124 + resg@gsnSpreadColors = True ; use full colormap
1.125 + resg@cnLinesOn = False ; Turn off contourn lines
1.126 + resg@mpFillOn = False ; Turn off map fill
1.127 + resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
1.128 +
1.129 +;************************************************
1.130 +; plot global biome class: (1) observed
1.131 +;************************************************
1.132 +
1.133 + resg@cnMinLevelValF = 1. ; Min level
1.134 + resg@cnMaxLevelValF = 19. ; Max level
1.135 + resg@cnLevelSpacingF = 1. ; interval
1.136 +
1.137 + classob@_FillValue = 1.e+36
1.138 + classob = where(classob.eq.0,classob@_FillValue,classob)
1.139 +
1.140 + plot_name = "global_class_ob"
1.141 + title = biome_name_ob
1.142 + resg@tiMainString = title
1.143 +
1.144 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.145 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.146 +
1.147 + plot = gsn_csm_contour_map_ce(wks,classob,resg)
1.148 +
1.149 + delete (wks)
1.150 +
1.151 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.152 + "rm "+plot_name+"."+plot_type)
1.153 +
1.154 +;************************************************
1.155 +; plot global biome class: (2) model
1.156 +;************************************************
1.157 +
1.158 + resg@cnMinLevelValF = 0. ; Min level
1.159 + resg@cnMaxLevelValF = 16. ; Max level
1.160 + resg@cnLevelSpacingF = 1. ; interval
1.161 +
1.162 + plot_name = "global_class_model"
1.163 + title = biome_name_mod
1.164 + resg@tiMainString = title
1.165 +
1.166 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.167 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.168 +
1.169 + plot = gsn_csm_contour_map_ce(wks,classmod,resg)
1.170 +
1.171 + delete (wks)
1.172 +
1.173 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.174 + "rm "+plot_name+"."+plot_type)
1.175 +
1.176 +;*******************************************************************
1.177 +; for html table : all 3 components (Mean, Max, Phase)
1.178 +;*******************************************************************
1.179 +
1.180 +; column (not including header column)
1.181 +
1.182 + component = (/"Mean","Max","Phase"/)
1.183 +
1.184 + col_head = (/"observed",model_name,"M_score" \
1.185 + ,"observed",model_name,"M_score" \
1.186 + ,"observed",model_name,"M_score" \
1.187 + /)
1.188 +
1.189 + n_comp = dimsizes(component)
1.190 + ncol = dimsizes(col_head )
1.191 +
1.192 +; row (not including header row)
1.193 +
1.194 +;----------------------------------------------------
1.195 +; using model biome class:
1.196 + row_head = (/"Not Vegetated" \
1.197 + ,"Needleleaf Evergreen Temperate Tree" \
1.198 + ,"Needleleaf Evergreen Boreal Tree" \
1.199 +; ,"Needleleaf Deciduous Boreal Tree" \
1.200 + ,"Broadleaf Evergreen Tropical Tree" \
1.201 + ,"Broadleaf Evergreen Temperate Tree" \
1.202 + ,"Broadleaf Deciduous Tropical Tree" \
1.203 + ,"Broadleaf Deciduous Temperate Tree" \
1.204 +; ,"Broadleaf Deciduous Boreal Tree" \
1.205 +; ,"Broadleaf Evergreen Shrub" \
1.206 + ,"Broadleaf Deciduous Temperate Shrub" \
1.207 + ,"Broadleaf Deciduous Boreal Shrub" \
1.208 + ,"C3 Arctic Grass" \
1.209 + ,"C3 Non-Arctic Grass" \
1.210 + ,"C4 Grass" \
1.211 + ,"Corn" \
1.212 +; ,"Wheat" \
1.213 + ,"All Biomes" \
1.214 + /)
1.215 +
1.216 + nrow = dimsizes(row_head)
1.217 +
1.218 +; arrays to be passed to table.
1.219 + text = new ((/nrow, ncol/),string )
1.220 +
1.221 +; total M_score
1.222 + M_total = 0.
1.223 +
1.224 +;********************************************************************
1.225 +; use land-type class to bin the data in equally spaced ranges
1.226 +;********************************************************************
1.227 +
1.228 +; using model biome class
1.229 + nclass = nclass_mod
1.230 +
1.231 + range = fspan(0,nclass,nclass+1)
1.232 +
1.233 +; Use this range information to grab all the values in a
1.234 +; particular range, and then take an average.
1.235 +
1.236 + nx = dimsizes(range) - 1
1.237 +
1.238 +; for 2 data: model and observed
1.239 + data_n = 2
1.240 +
1.241 +; using model biome class
1.242 +
1.243 + base = ndtooned(classmod)
1.244 +
1.245 +; output
1.246 +
1.247 + yvalues = new((/data_n,nx/),float)
1.248 + count = new((/data_n,nx/),float)
1.249 +
1.250 +;************************************************************************
1.251 +; go through all components
1.252 +;************************************************************************
1.253 +
1.254 + do m = 0,n_comp-1
1.255 +
1.256 +;===================
1.257 +; get data:
1.258 +;===================
1.259 +; (A) Mean
1.260 +
1.261 + if (m .eq. 0) then
1.262 + data_ob = dim_avg_Wrap(laiob (lat|:,lon|:,time|:))
1.263 + data_mod = dim_avg_Wrap(laimod(lat|:,lon|:,time|:))
1.264 + end if
1.265 +
1.266 +; (B) Max
1.267 +
1.268 + if (m .eq. 1) then
1.269 +
1.270 +; observed
1.271 + data_ob = laiob(0,:,:)
1.272 + data_ob@long_name = "Leaf Area Index Max"
1.273 +
1.274 + do j = 0,nlat-1
1.275 + do i = 0,nlon-1
1.276 + data_ob(j,i) = max(laiob(:,j,i))
1.277 + end do
1.278 + end do
1.279 +
1.280 +; model
1.281 + data_mod = laimod(0,:,:)
1.282 + data_mod@long_name = "Leaf Area Index Max"
1.283 +
1.284 + do j = 0,nlat-1
1.285 + do i = 0,nlon-1
1.286 + data_mod(j,i) = max(laimod(:,j,i))
1.287 + end do
1.288 + end do
1.289 +
1.290 + end if
1.291 +
1.292 +; (C) phase
1.293 +
1.294 + if (m .eq. 2) then
1.295 +
1.296 +; observed
1.297 + data_ob = laiob(0,:,:)
1.298 + data_ob@long_name = "Leaf Area Index Max Month"
1.299 +
1.300 + do j = 0,nlat-1
1.301 + do i = 0,nlon-1
1.302 + data_ob(j,i) = maxind(laiob(:,j,i)) + 1
1.303 + end do
1.304 + end do
1.305 +
1.306 +; model
1.307 + data_mod = laimod(0,:,:)
1.308 + data_mod@long_name = "Leaf Area Index Max Month"
1.309 +
1.310 + do j = 0,nlat-1
1.311 + do i = 0,nlon-1
1.312 + data_mod(j,i) = maxind(laimod(:,j,i)) + 1
1.313 + end do
1.314 + end do
1.315 +
1.316 + end if
1.317 +
1.318 +;==============================
1.319 +; put data into bins
1.320 +;==============================
1.321 +
1.322 +; Loop through each range, using base
1.323 +
1.324 + do i=0,nx-1
1.325 +
1.326 + if (i.ne.(nx-1)) then
1.327 + idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
1.328 + else
1.329 + idx = ind(base.ge.range(i))
1.330 + end if
1.331 +
1.332 +; loop through each dataset
1.333 +
1.334 + do n = 0,data_n-1
1.335 +
1.336 + if (n .eq. 0) then
1.337 + data = ndtooned(data_ob)
1.338 + end if
1.339 +
1.340 + if (n .eq. 1) then
1.341 + data = ndtooned(data_mod)
1.342 + end if
1.343 +
1.344 +; Calculate average
1.345 +
1.346 + if (.not.any(ismissing(idx))) then
1.347 + yvalues(n,i) = avg(data(idx))
1.348 + count(n,i) = dimsizes(idx)
1.349 + else
1.350 + yvalues(n,i) = yvalues@_FillValue
1.351 + count(n,i) = 0
1.352 + end if
1.353 +
1.354 +;#############################################################
1.355 +; using model biome class:
1.356 +;
1.357 +; set the following 4 classes to _FillValue:
1.358 +; (3)Needleleaf Deciduous Boreal Tree,
1.359 +; (8)Broadleaf Deciduous Boreal Tree,
1.360 +; (9)Broadleaf Evergreen Shrub,
1.361 +; (16)Wheat
1.362 +
1.363 + if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
1.364 + yvalues(n,i) = yvalues@_FillValue
1.365 + count(n,i) = 0
1.366 + end if
1.367 +;#############################################################
1.368 +
1.369 + delete(data)
1.370 + end do ; n-loop
1.371 +
1.372 + delete(idx)
1.373 + end do ; i-loop
1.374 +
1.375 +;=====================================
1.376 +; compute correlation coef and M score
1.377 +;=====================================
1.378 +
1.379 + score_max = 5.0
1.380 +
1.381 + u = yvalues(0,:)
1.382 + v = yvalues(1,:)
1.383 + u_count = count(0,:)
1.384 + v_count = count(1,:)
1.385 +
1.386 + good = ind(.not.ismissing(u) .and. .not.ismissing(v))
1.387 +
1.388 + uu = u(good)
1.389 + vv = v(good)
1.390 + uu_count = u_count(good)
1.391 + vv_count = v_count(good)
1.392 +
1.393 +; compute correlation coef
1.394 + cc = esccr(uu,vv,0)
1.395 +
1.396 + if (n .eq. 2) then
1.397 + bias = avg(abs(vv-uu))
1.398 + bias = where((bias.gt. 6.),12.-bias,bias)
1.399 + Mscore = ((6. - bias)/6.)*score_max
1.400 + else
1.401 + bias = sum(abs(vv-uu)/abs(vv+uu))
1.402 + Mscore = (1.- (bias/dimsizes(uu)))*score_max
1.403 + end if
1.404 +
1.405 + M_score = sprintf("%.2f", Mscore)
1.406 +
1.407 +; compute M_total
1.408 +
1.409 + M_total = M_total + Mscore
1.410 +
1.411 +;================================
1.412 +; output M_score to score sheet
1.413 +;===============================
1.414 +
1.415 + if (isvar("compare")) then
1.416 + system("sed -e '1,/M_lai_"+component(m)+"/s/M_lai_"+component(m)+"/"+M_score+"/' "+html_name2+" > "+html_new2+";"+ \
1.417 + "mv -f "+html_new2+" "+html_name2)
1.418 + end if
1.419 +
1.420 + system("sed s#M_lai_"+component(m)+"#"+M_score+"# "+html_name+" > "+html_new+";"+ \
1.421 + "mv -f "+html_new+" "+html_name)
1.422 +
1.423 +;==============================
1.424 +; output M_score to html table
1.425 +;==============================
1.426 +
1.427 + n = m*3
1.428 +
1.429 + do i=0,nrow-2
1.430 + text(i,n) = sprintf("%.1f",uu(i))
1.431 + text(i,n+1) = sprintf("%.1f",vv(i))
1.432 + text(i,n+2) = "-"
1.433 + end do
1.434 + text(nrow-1,n) = sprintf("%.1f",sum(uu*uu_count)/sum(uu_count))
1.435 + text(nrow-1,n+1) = sprintf("%.1f",sum(vv*vv_count)/sum(vv_count))
1.436 + text(nrow-1,n+2) = M_score
1.437 +
1.438 + delete (u)
1.439 + delete (v)
1.440 + delete (uu)
1.441 + delete (vv)
1.442 + delete (u_count)
1.443 + delete (v_count)
1.444 + delete (uu_count)
1.445 + delete (vv_count)
1.446 + delete (good)
1.447 +
1.448 +;========================================
1.449 +; global res changes for each component
1.450 +;========================================
1.451 + delta = 0.00001
1.452 +
1.453 + if (m .eq. 0) then
1.454 + resg@cnMinLevelValF = 0.
1.455 + resg@cnMaxLevelValF = 10.
1.456 + resg@cnLevelSpacingF = 1.
1.457 +
1.458 + data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)
1.459 + end if
1.460 +
1.461 + if (m .eq. 1) then
1.462 + resg@cnMinLevelValF = 0.
1.463 + resg@cnMaxLevelValF = 10.
1.464 + resg@cnLevelSpacingF = 1.
1.465 +
1.466 + data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)
1.467 + end if
1.468 +
1.469 + if (m .eq. 2) then
1.470 + resg@cnMinLevelValF = 1.
1.471 + resg@cnMaxLevelValF = 12.
1.472 + resg@cnLevelSpacingF = 1.
1.473 +
1.474 + data_mod = where(ismissing(data_ob).or.(data_mod.lt.delta),data_mod@_FillValue,data_mod)
1.475 + data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)
1.476 +
1.477 + end if
1.478 +
1.479 +;=========================
1.480 +; global contour : ob
1.481 +;=========================
1.482 +
1.483 + plot_name = "global_"+component(m)+"_ob"
1.484 + title = ob_name
1.485 + resg@tiMainString = title
1.486 +
1.487 + wks = gsn_open_wks (plot_type,plot_name) ; open workstation
1.488 + gsn_define_colormap(wks,"gui_default") ; choose colormap
1.489 +
1.490 + plot = gsn_csm_contour_map_ce(wks,data_ob,resg)
1.491 +
1.492 + delete (wks)
1.493 + delete (plot)
1.494 +
1.495 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.496 + "rm "+plot_name+"."+plot_type)
1.497 +
1.498 +;============================
1.499 +; global contour : model
1.500 +;============================
1.501 +
1.502 + plot_name = "global_"+component(m)+"_model"
1.503 + title = "Model " + model_name
1.504 + resg@tiMainString = title
1.505 +
1.506 + wks = gsn_open_wks (plot_type,plot_name)
1.507 + gsn_define_colormap(wks,"gui_default")
1.508 +
1.509 + plot = gsn_csm_contour_map_ce(wks,data_mod,resg)
1.510 +
1.511 + delete (wks)
1.512 + delete (plot)
1.513 +
1.514 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.515 + "rm "+plot_name+"."+plot_type)
1.516 +
1.517 +;================================
1.518 +; global contour: model vs ob
1.519 +;================================
1.520 +
1.521 + plot_name = "global_"+component(m)+"_model_vs_ob"
1.522 +
1.523 + wks = gsn_open_wks (plot_type,plot_name)
1.524 + gsn_define_colormap(wks,"gui_default")
1.525 +
1.526 + plot=new(3,graphic) ; create graphic array
1.527 +
1.528 + resg@gsnFrame = False ; Do not draw plot
1.529 + resg@gsnDraw = False ; Do not advance frame
1.530 +
1.531 +; plot correlation coef
1.532 +
1.533 + gRes = True
1.534 + gRes@txFontHeightF = 0.02
1.535 + gRes@txAngleF = 90
1.536 +
1.537 + correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")"
1.538 +
1.539 + gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
1.540 +
1.541 +; plot ob
1.542 +
1.543 + title = ob_name
1.544 + resg@tiMainString = title
1.545 +
1.546 + plot(0) = gsn_csm_contour_map_ce(wks,data_ob,resg)
1.547 +
1.548 +; plot model
1.549 +
1.550 + title = "Model "+ model_name
1.551 + resg@tiMainString = title
1.552 +
1.553 + plot(1) = gsn_csm_contour_map_ce(wks,data_mod,resg)
1.554 +
1.555 +; plot model-ob
1.556 +
1.557 + if (m .eq. 0) then
1.558 + resg@cnMinLevelValF = -2.
1.559 + resg@cnMaxLevelValF = 2.
1.560 + resg@cnLevelSpacingF = 0.4
1.561 + end if
1.562 +
1.563 + if (m .eq. 1) then
1.564 + resg@cnMinLevelValF = -6.
1.565 + resg@cnMaxLevelValF = 6.
1.566 + resg@cnLevelSpacingF = 1.
1.567 + end if
1.568 +
1.569 + if (m .eq. 2) then
1.570 + resg@cnMinLevelValF = -6.
1.571 + resg@cnMaxLevelValF = 6.
1.572 + resg@cnLevelSpacingF = 1.
1.573 + end if
1.574 +
1.575 + zz = data_mod
1.576 + zz = data_mod - data_ob
1.577 + title = "Model_"+model_name+" - Observed"
1.578 + resg@tiMainString = title
1.579 +
1.580 + plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
1.581 +
1.582 +; plot panel
1.583 +
1.584 + pres = True ; panel plot mods desired
1.585 + pres@gsnMaximize = True ; fill the page
1.586 +
1.587 + gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
1.588 +
1.589 + delete (wks)
1.590 + delete (plot)
1.591 +
1.592 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.593 + "rm "+plot_name+"."+plot_type)
1.594 +
1.595 + delete (data_ob)
1.596 + delete (data_mod)
1.597 +
1.598 + resg@gsnFrame = True ; Do advance frame
1.599 + resg@gsnDraw = True ; Do draw plot
1.600 +
1.601 + end do ; m-loop
1.602 +
1.603 +;**************************************************
1.604 +; html table
1.605 +;**************************************************
1.606 + output_html = "table_model_vs_ob.html"
1.607 +
1.608 + header_text = "<H1>LAI: Model "+model_name+" vs. MODIS observations</H1>"
1.609 +
1.610 + header = (/"<HTML>" \
1.611 + ,"<HEAD>" \
1.612 + ,"<TITLE>CLAMP metrics</TITLE>" \
1.613 + ,"</HEAD>" \
1.614 + ,header_text \
1.615 + /)
1.616 + footer = "</HTML>"
1.617 +
1.618 + table_header = (/ \
1.619 + "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
1.620 + ,"<tr>" \
1.621 + ," <th bgcolor=DDDDDD rowspan=2>Biome Class</th>" \
1.622 + ," <th bgcolor=DDDDDD colspan=3>"+component(0)+" (m2/m2)</th>" \
1.623 + ," <th bgcolor=DDDDDD colspan=3>"+component(1)+" (m2/m2)</th>" \
1.624 + ," <th bgcolor=DDDDDD colspan=3>"+component(2)+" (months)</th>" \
1.625 + ," <th bgcolor=DDDDDD rowspan=2>Annual Plot</th>" \
1.626 + ,"</tr>" \
1.627 + ,"<tr>" \
1.628 + ," <th bgcolor=DDDDDD >observed</th>" \
1.629 + ," <th bgcolor=DDDDDD >"+model_name+"</th>" \
1.630 + ," <th bgcolor=DDDDDD >M_score</th>" \
1.631 + ," <th bgcolor=DDDDDD >observed</th>" \
1.632 + ," <th bgcolor=DDDDDD >"+model_name+"</th>" \
1.633 + ," <th bgcolor=DDDDDD >M_score</th>" \
1.634 + ," <th bgcolor=DDDDDD >observed</th>" \
1.635 + ," <th bgcolor=DDDDDD >"+model_name+"</th>" \
1.636 + ," <th bgcolor=DDDDDD >M_score</th>" \
1.637 + ,"</tr>" \
1.638 + /)
1.639 + table_footer = "</table>"
1.640 + row_header = "<tr>"
1.641 + row_footer = "</tr>"
1.642 +
1.643 + lines = new(50000,string)
1.644 + nline = 0
1.645 +
1.646 + set_line(lines,nline,header)
1.647 + set_line(lines,nline,table_header)
1.648 +;-----------------------------------------------
1.649 +;row of table
1.650 +
1.651 + do n = 0,nrow-1
1.652 + set_line(lines,nline,row_header)
1.653 +
1.654 + txt1 = row_head(n)
1.655 + txt2 = text(n,0)
1.656 + txt3 = text(n,1)
1.657 + txt4 = text(n,2)
1.658 + txt5 = text(n,3)
1.659 + txt6 = text(n,4)
1.660 + txt7 = text(n,5)
1.661 + txt8 = text(n,6)
1.662 + txt9 = text(n,7)
1.663 + txt10 = text(n,8)
1.664 + txt11 = "<a href=./annual_biome_"+n+".png>model_vs_ob</a>"
1.665 +
1.666 + set_line(lines,nline,"<th>"+txt1+"</th>")
1.667 + set_line(lines,nline,"<th>"+txt2+"</th>")
1.668 + set_line(lines,nline,"<th>"+txt3+"</th>")
1.669 + set_line(lines,nline,"<th>"+txt4+"</th>")
1.670 + set_line(lines,nline,"<th>"+txt5+"</th>")
1.671 + set_line(lines,nline,"<th>"+txt6+"</th>")
1.672 + set_line(lines,nline,"<th>"+txt7+"</th>")
1.673 + set_line(lines,nline,"<th>"+txt8+"</th>")
1.674 + set_line(lines,nline,"<th>"+txt9+"</th>")
1.675 + set_line(lines,nline,"<th>"+txt10+"</th>")
1.676 + set_line(lines,nline,"<th>"+txt11+"</th>")
1.677 +
1.678 + set_line(lines,nline,row_footer)
1.679 + end do
1.680 +;-----------------------------------------------
1.681 + set_line(lines,nline,table_footer)
1.682 + set_line(lines,nline,footer)
1.683 +
1.684 +; Now write to an HTML file
1.685 +
1.686 + idx = ind(.not.ismissing(lines))
1.687 + if(.not.any(ismissing(idx))) then
1.688 + asciiwrite(output_html,lines(idx))
1.689 + else
1.690 + print ("error?")
1.691 + end if
1.692 +
1.693 + delete (yvalues)
1.694 + delete (count)
1.695 + delete (idx)
1.696 +
1.697 +;************************************************************************
1.698 +; go through all ntime
1.699 +;************************************************************************
1.700 +
1.701 +; for 2 data: model and observed
1.702 + data_n = 2
1.703 +
1.704 +; using model biome class
1.705 +
1.706 + base = ndtooned(classmod)
1.707 +
1.708 +; output
1.709 +
1.710 + yvalues = new((/data_n,nx,ntime/),float)
1.711 +
1.712 +;==============================
1.713 +; put data into bins
1.714 +;==============================
1.715 +
1.716 +; Loop through each range, using base
1.717 +
1.718 + do i=0,nx-1
1.719 +
1.720 + if (i.ne.(nx-1)) then
1.721 + idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
1.722 + else
1.723 + idx = ind(base.ge.range(i))
1.724 + end if
1.725 +
1.726 +; loop through each dataset
1.727 +
1.728 + do n = 0,data_n-1
1.729 +
1.730 + do m = 0,ntime-1
1.731 +
1.732 + if (n .eq. 0) then
1.733 + data = ndtooned(laiob (m,:,:))
1.734 + end if
1.735 +
1.736 + if (n .eq. 1) then
1.737 + data = ndtooned(laimod(m,:,:))
1.738 + end if
1.739 +
1.740 +; Calculate average
1.741 +
1.742 + if (.not.any(ismissing(idx))) then
1.743 + yvalues(n,i,m) = avg(data(idx))
1.744 + else
1.745 + yvalues(n,i,m) = yvalues@_FillValue
1.746 + end if
1.747 +
1.748 +;#############################################################
1.749 +; using model biome class:
1.750 +;
1.751 +; set the following 4 classes to _FillValue:
1.752 +; (3)Needleleaf Deciduous Boreal Tree,
1.753 +; (8)Broadleaf Deciduous Boreal Tree,
1.754 +; (9)Broadleaf Evergreen Shrub,
1.755 +; (16)Wheat
1.756 +
1.757 + if (i.eq.3 .or. i.eq.8 .or. i.eq.9 .or. i.eq.16) then
1.758 + yvalues(n,i,m) = yvalues@_FillValue
1.759 + end if
1.760 +;#############################################################
1.761 +
1.762 + end do ; m-loop
1.763 +
1.764 + delete(data)
1.765 + end do ; n-loop
1.766 +
1.767 + delete(idx)
1.768 + end do ; i-loop
1.769 +
1.770 + good = ind(.not.ismissing(yvalues(0,:,0)))
1.771 +
1.772 + n_biome = dimsizes(good)
1.773 +
1.774 +;----------------------------------------------------------------
1.775 +; data for tseries plot
1.776 +
1.777 + yvalues_g = new((/data_n,n_biome,ntime/),float)
1.778 +
1.779 + yvalues_g@units = ""
1.780 +
1.781 + yvalues_g = yvalues(:,good,:)
1.782 +
1.783 + delete (good)
1.784 +
1.785 +;*******************************************************************
1.786 +; res for line plot
1.787 +;*******************************************************************
1.788 +; for x-axis in xyplot
1.789 + mon = ispan(1,12,1)
1.790 + mon@long_name = "month"
1.791 +
1.792 + res = True ; plot mods desired
1.793 + res@xyLineThicknesses = (/2.0,2.0,2.0/) ; make 2nd lines thicker
1.794 + res@xyLineColors = (/"blue","red"/) ; change line color
1.795 +
1.796 +; res@tiMainFontHeightF = 0.025 ; size of title
1.797 +
1.798 +; Add a boxed legend using the more simple method
1.799 + res@pmLegendDisplayMode = "Always"
1.800 +; res@pmLegendWidthF = 0.1
1.801 + res@pmLegendWidthF = 0.08
1.802 + res@pmLegendHeightF = 0.06
1.803 +; res@pmLegendOrthogonalPosF = -1.17
1.804 +; res@pmLegendOrthogonalPosF = -1.00 ;(downward)
1.805 + res@pmLegendOrthogonalPosF = -0.30 ;(downward)
1.806 +
1.807 +; res@pmLegendParallelPosF = 0.18
1.808 + res@pmLegendParallelPosF = 0.23 ;(rightward)
1.809 +
1.810 +; res@lgPerimOn = False
1.811 + res@lgLabelFontHeightF = 0.015
1.812 + res@xyExplicitLegendLabels = (/"observed",model_name/)
1.813 +
1.814 +;************************************************************
1.815 +
1.816 + plot_data = new((/2,ntime/),float)
1.817 + plot_data@long_name = "Leaf Area Index"
1.818 +
1.819 +;----------------------------------------------
1.820 +; time series plot : per biome
1.821 +
1.822 + do m = 0, n_biome-1
1.823 +
1.824 + plot_name = "annual_biome_"+ m
1.825 +
1.826 + wks = gsn_open_wks (plot_type,plot_name)
1.827 +
1.828 + title = "LAI : "+ row_head(m)
1.829 + res@tiMainString = title
1.830 +
1.831 + plot_data(0,:) = yvalues_g(0,m,:)
1.832 + plot_data(1,:) = yvalues_g(1,m,:)
1.833 +
1.834 + plot = gsn_csm_xy(wks,mon,plot_data,res)
1.835 +
1.836 + delete (wks)
1.837 + delete (plot)
1.838 +
1.839 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.840 + "rm "+plot_name+"."+plot_type)
1.841 + end do
1.842 +
1.843 +;--------------------------------------------
1.844 +; time series plot: all biome
1.845 +
1.846 + plot_name = "annual_biome_"+ n_biome
1.847 +
1.848 + wks = gsn_open_wks (plot_type,plot_name)
1.849 +
1.850 + title = "LAI : "+ row_head(n_biome)
1.851 + res@tiMainString = title
1.852 +
1.853 + do k = 0,ntime-1
1.854 + plot_data(0,k) = avg(yvalues_g(0,:,k))
1.855 + plot_data(1,k) = avg(yvalues_g(1,:,k))
1.856 + end do
1.857 +
1.858 + plot = gsn_csm_xy(wks,mon,plot_data,res)
1.859 +
1.860 + delete (wks)
1.861 + delete (plot)
1.862 +
1.863 + system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
1.864 + "rm "+plot_name+"."+plot_type)
1.865 +
1.866 + delete (plot_data)
1.867 +
1.868 +;***************************************************************************
1.869 +; write total score to file
1.870 +;***************************************************************************
1.871 +
1.872 + asciiwrite("M_save.lai", M_total)
1.873 +
1.874 +;***************************************************************************
1.875 +; output plots
1.876 +;***************************************************************************
1.877 + output_dir = model_name+"/lai"
1.878 +
1.879 + system("mv *.png *.html " + output_dir)
1.880 +;***************************************************************************
1.881 +end
1.882 +