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 data: model forrest@0: forrest@0: fm = addfile(dirm+film10,"r") forrest@0: laimod = fm->TLAI forrest@0: forrest@0: delete (fm) forrest@0: forrest@0: dsizes = dimsizes(laimod) forrest@0: ntime = dsizes(0) forrest@0: nlat = dsizes(1) forrest@0: nlon = dsizes(2) 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: delete (fm_l) 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: ob forrest@0: forrest@0: ;---------------------------------- forrest@0: ; read biome data: observed forrest@0: forrest@0: biome_name_ob = "MODIS LandCover" forrest@0: forrest@0: dir_c = diro + "lai/" forrest@0: filo_c = "land_class_"+model_grid+".nc" forrest@0: fo = addfile(dir_c+filo_c,"r") forrest@0: classob = tofloat(fo->LAND_CLASS) forrest@0: forrest@0: delete (fo) forrest@0: forrest@0: ; input observed data has 20 land-type classes forrest@0: nclass_ob = 20 forrest@0: forrest@0: ;--------------------------------- forrest@0: ; read lai data: observed forrest@0: forrest@0: ;ob_name = "MODIS MOD 15A2 2000-2005" forrest@0: ob_name = "MODIS MOD 15A2 2000-2004" forrest@0: forrest@0: dir_l = diro + "lai/" forrest@0: ;filo = "LAI_2000-2005_MONS_"+model_grid+".nc" forrest@0: filo = "LAI_2000-2004_MONS_"+model_grid+".nc" forrest@0: fo = addfile(dir_l+filo,"r") forrest@0: laiob = fo->LAI forrest@0: forrest@0: delete (fo) forrest@0: forrest@0: ;------------------------------------------------- forrest@0: ; take into account landfrac forrest@0: forrest@0: laimod = laimod * conform(laimod,landfrac,(/1,2/)) forrest@0: laiob = laiob * conform(laiob ,landfrac,(/1,2/)) forrest@0: forrest@0: delete (landfrac) forrest@0: forrest@0: ;************************************************ forrest@0: ; global res forrest@0: ;************************************************ forrest@0: resg = True ; Use plot options forrest@0: resg@cnFillOn = True ; Turn on color fill forrest@0: resg@gsnSpreadColors = True ; use full colormap forrest@0: resg@cnLinesOn = False ; Turn off contourn lines forrest@0: resg@mpFillOn = False ; Turn off map fill forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals forrest@0: forrest@0: ;************************************************ forrest@0: ; plot global biome class: (1) observed forrest@0: ;************************************************ forrest@0: forrest@0: resg@cnMinLevelValF = 1. ; Min level forrest@0: resg@cnMaxLevelValF = 19. ; Max level forrest@0: resg@cnLevelSpacingF = 1. ; interval forrest@0: forrest@0: classob@_FillValue = 1.e+36 forrest@0: classob = where(classob.eq.0,classob@_FillValue,classob) forrest@0: forrest@0: plot_name = "global_class_ob" forrest@0: title = biome_name_ob forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: plot = gsn_csm_contour_map_ce(wks,classob,resg) forrest@0: forrest@0: delete (wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: ;************************************************ forrest@0: ; plot global biome class: (2) model forrest@0: ;************************************************ forrest@0: forrest@0: resg@cnMinLevelValF = 0. ; Min level forrest@0: resg@cnMaxLevelValF = 16. ; Max level forrest@0: resg@cnLevelSpacingF = 1. ; interval forrest@0: forrest@0: plot_name = "global_class_model" forrest@0: title = biome_name_mod forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: plot = gsn_csm_contour_map_ce(wks,classmod,resg) forrest@0: forrest@0: delete (wks) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: ;******************************************************************* forrest@0: ; for html table : all 3 components (Mean, Max, Phase) forrest@0: ;******************************************************************* forrest@0: forrest@0: ; column (not including header column) forrest@0: forrest@0: component = (/"Mean","Max","Phase"/) forrest@0: forrest@0: col_head = (/"observed",model_name,"M_score" \ forrest@0: ,"observed",model_name,"M_score" \ forrest@0: ,"observed",model_name,"M_score" \ forrest@0: /) forrest@0: forrest@0: n_comp = dimsizes(component) 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: 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 Biomes" \ 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: ; total M_score forrest@0: M_total = 0. 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: ; for 2 data: model 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: ;************************************************************************ forrest@0: ; go through all components forrest@0: ;************************************************************************ forrest@0: forrest@0: do m = 0,n_comp-1 forrest@0: forrest@0: ;=================== forrest@0: ; get data: forrest@0: ;=================== forrest@0: ; (A) Mean forrest@0: forrest@0: if (m .eq. 0) then forrest@0: data_ob = dim_avg_Wrap(laiob (lat|:,lon|:,time|:)) forrest@0: data_mod = dim_avg_Wrap(laimod(lat|:,lon|:,time|:)) forrest@0: end if forrest@0: forrest@0: ; (B) Max forrest@0: forrest@0: if (m .eq. 1) then forrest@0: forrest@0: ; observed forrest@0: data_ob = laiob(0,:,:) forrest@0: data_ob@long_name = "Leaf Area Index Max" forrest@0: forrest@0: do j = 0,nlat-1 forrest@0: do i = 0,nlon-1 forrest@0: data_ob(j,i) = max(laiob(:,j,i)) forrest@0: end do forrest@0: end do forrest@0: forrest@0: ; model forrest@0: data_mod = laimod(0,:,:) forrest@0: data_mod@long_name = "Leaf Area Index Max" forrest@0: forrest@0: do j = 0,nlat-1 forrest@0: do i = 0,nlon-1 forrest@0: data_mod(j,i) = max(laimod(:,j,i)) forrest@0: end do forrest@0: end do forrest@0: forrest@0: end if forrest@0: forrest@0: ; (C) phase forrest@0: forrest@0: if (m .eq. 2) then forrest@0: forrest@0: ; observed forrest@0: data_ob = laiob(0,:,:) forrest@0: data_ob@long_name = "Leaf Area Index Max Month" forrest@0: forrest@0: do j = 0,nlat-1 forrest@0: do i = 0,nlon-1 forrest@0: data_ob(j,i) = maxind(laiob(:,j,i)) + 1 forrest@0: end do forrest@0: end do forrest@0: forrest@0: ; model forrest@0: data_mod = laimod(0,:,:) forrest@0: data_mod@long_name = "Leaf Area Index Max Month" forrest@0: forrest@0: do j = 0,nlat-1 forrest@0: do i = 0,nlon-1 forrest@0: data_mod(j,i) = maxind(laimod(:,j,i)) + 1 forrest@0: end do forrest@0: end do forrest@0: forrest@0: end if forrest@0: forrest@0: ;============================== forrest@0: ; put data into bins forrest@0: ;============================== 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(data_ob) forrest@0: end if forrest@0: forrest@0: if (n .eq. 1) then forrest@0: data = ndtooned(data_mod) 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: ;===================================== forrest@0: ; compute correlation coef and M score forrest@0: ;===================================== forrest@0: forrest@0: score_max = 5.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: ; compute correlation coef forrest@0: cc = esccr(uu,vv,0) forrest@0: forrest@0: if (n .eq. 2) then forrest@0: bias = avg(abs(vv-uu)) forrest@0: bias = where((bias.gt. 6.),12.-bias,bias) forrest@0: Mscore = ((6. - bias)/6.)*score_max forrest@0: else forrest@0: bias = sum(abs(vv-uu)/abs(vv+uu)) forrest@0: Mscore = (1.- (bias/dimsizes(uu)))*score_max forrest@0: end if forrest@0: forrest@0: M_score = sprintf("%.2f", Mscore) forrest@0: forrest@0: ; compute M_total forrest@0: forrest@0: M_total = M_total + Mscore forrest@0: forrest@0: ;================================ forrest@0: ; output M_score to score sheet forrest@0: ;=============================== forrest@0: forrest@0: if (isvar("compare")) then forrest@0: system("sed -e '1,/M_lai_"+component(m)+"/s/M_lai_"+component(m)+"/"+M_score+"/' "+html_name2+" > "+html_new2+";"+ \ forrest@0: "mv -f "+html_new2+" "+html_name2) forrest@0: end if forrest@0: forrest@0: system("sed s#M_lai_"+component(m)+"#"+M_score+"# "+html_name+" > "+html_new+";"+ \ forrest@0: "mv -f "+html_new+" "+html_name) forrest@0: forrest@0: ;============================== forrest@0: ; output M_score to html table forrest@0: ;============================== forrest@0: forrest@0: n = m*3 forrest@0: forrest@0: do i=0,nrow-2 forrest@0: text(i,n) = sprintf("%.1f",uu(i)) forrest@0: text(i,n+1) = sprintf("%.1f",vv(i)) forrest@0: text(i,n+2) = "-" forrest@0: end do forrest@0: text(nrow-1,n) = sprintf("%.1f",sum(uu*uu_count)/sum(uu_count)) forrest@0: text(nrow-1,n+1) = sprintf("%.1f",sum(vv*vv_count)/sum(vv_count)) forrest@0: text(nrow-1,n+2) = M_score forrest@0: forrest@0: delete (u) forrest@0: delete (v) forrest@0: delete (uu) forrest@0: delete (vv) forrest@0: delete (u_count) forrest@0: delete (v_count) forrest@0: delete (uu_count) forrest@0: delete (vv_count) forrest@0: delete (good) forrest@0: forrest@0: ;======================================== forrest@0: ; global res changes for each component forrest@0: ;======================================== forrest@0: delta = 0.00001 forrest@0: forrest@0: if (m .eq. 0) then forrest@0: resg@cnMinLevelValF = 0. forrest@0: resg@cnMaxLevelValF = 10. forrest@0: resg@cnLevelSpacingF = 1. forrest@0: forrest@0: data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob) forrest@0: end if forrest@0: forrest@0: if (m .eq. 1) then forrest@0: resg@cnMinLevelValF = 0. forrest@0: resg@cnMaxLevelValF = 10. forrest@0: resg@cnLevelSpacingF = 1. forrest@0: forrest@0: data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob) forrest@0: end if forrest@0: forrest@0: if (m .eq. 2) then forrest@0: resg@cnMinLevelValF = 1. forrest@0: resg@cnMaxLevelValF = 12. forrest@0: resg@cnLevelSpacingF = 1. forrest@0: forrest@0: data_mod = where(ismissing(data_ob).or.(data_mod.lt.delta),data_mod@_FillValue,data_mod) forrest@0: data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob) forrest@0: forrest@0: end if forrest@0: forrest@0: ;========================= forrest@0: ; global contour : ob forrest@0: ;========================= forrest@0: forrest@0: plot_name = "global_"+component(m)+"_ob" forrest@0: title = ob_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) ; open workstation forrest@0: gsn_define_colormap(wks,"gui_default") ; choose colormap forrest@0: forrest@0: plot = gsn_csm_contour_map_ce(wks,data_ob,resg) forrest@0: forrest@0: delete (wks) forrest@0: delete (plot) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: ;============================ forrest@0: ; global contour : model forrest@0: ;============================ forrest@0: forrest@0: plot_name = "global_"+component(m)+"_model" forrest@0: title = "Model " + model_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: gsn_define_colormap(wks,"gui_default") forrest@0: forrest@0: plot = gsn_csm_contour_map_ce(wks,data_mod,resg) forrest@0: forrest@0: delete (wks) forrest@0: delete (plot) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: ;================================ forrest@0: ; global contour: model vs ob forrest@0: ;================================ forrest@0: forrest@0: plot_name = "global_"+component(m)+"_model_vs_ob" forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: gsn_define_colormap(wks,"gui_default") forrest@0: forrest@0: plot=new(3,graphic) ; create graphic array forrest@0: forrest@0: resg@gsnFrame = False ; Do not draw plot forrest@0: resg@gsnDraw = False ; Do not advance frame forrest@0: forrest@0: ; plot correlation coef forrest@0: forrest@0: gRes = True forrest@0: gRes@txFontHeightF = 0.02 forrest@0: gRes@txAngleF = 90 forrest@0: forrest@0: correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) forrest@0: forrest@0: ; plot ob forrest@0: forrest@0: title = ob_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(0) = gsn_csm_contour_map_ce(wks,data_ob,resg) forrest@0: forrest@0: ; plot model forrest@0: forrest@0: title = "Model "+ model_name forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(1) = gsn_csm_contour_map_ce(wks,data_mod,resg) forrest@0: forrest@0: ; plot model-ob forrest@0: forrest@0: if (m .eq. 0) then forrest@0: resg@cnMinLevelValF = -2. forrest@0: resg@cnMaxLevelValF = 2. forrest@0: resg@cnLevelSpacingF = 0.4 forrest@0: end if forrest@0: forrest@0: if (m .eq. 1) then forrest@0: resg@cnMinLevelValF = -6. forrest@0: resg@cnMaxLevelValF = 6. forrest@0: resg@cnLevelSpacingF = 1. forrest@0: end if forrest@0: forrest@0: if (m .eq. 2) then forrest@0: resg@cnMinLevelValF = -6. forrest@0: resg@cnMaxLevelValF = 6. forrest@0: resg@cnLevelSpacingF = 1. forrest@0: end if forrest@0: forrest@0: zz = data_mod forrest@0: zz = data_mod - data_ob forrest@0: title = "Model_"+model_name+" - Observed" forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) forrest@0: forrest@0: ; plot panel forrest@0: forrest@0: pres = True ; panel plot mods desired forrest@0: pres@gsnMaximize = True ; fill the page forrest@0: forrest@0: gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot forrest@0: forrest@0: delete (wks) forrest@0: delete (plot) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: delete (data_ob) forrest@0: delete (data_mod) forrest@0: forrest@0: resg@gsnFrame = True ; Do advance frame forrest@0: resg@gsnDraw = True ; Do draw plot forrest@0: forrest@0: end do ; m-loop forrest@0: forrest@0: ;************************************************** forrest@0: ; html table forrest@0: ;************************************************** forrest@0: output_html = "table_model_vs_ob.html" forrest@0: forrest@0: header_text = "

LAI: Model "+model_name+" vs. MODIS observations

" 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: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ,"" \ forrest@0: ,"" \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ," " \ forrest@0: ,"" \ forrest@0: /) forrest@0: table_footer = "
Biome Class"+component(0)+" (m2/m2)"+component(1)+" (m2/m2)"+component(2)+" (months)Annual Plot
observed"+model_name+"M_scoreobserved"+model_name+"M_scoreobserved"+model_name+"M_score
" 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: txt10 = text(n,8) forrest@0: txt11 = "model_vs_ob" 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: set_line(lines,nline,""+txt10+"") forrest@0: set_line(lines,nline,""+txt11+"") 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: forrest@0: delete (yvalues) forrest@0: delete (count) forrest@0: delete (idx) forrest@0: forrest@0: ;************************************************************************ forrest@0: ; go through all ntime forrest@0: ;************************************************************************ forrest@0: forrest@0: ; for 2 data: model 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,ntime/),float) forrest@0: forrest@0: ;============================== forrest@0: ; put data into bins forrest@0: ;============================== 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: do m = 0,ntime-1 forrest@0: forrest@0: if (n .eq. 0) then forrest@0: data = ndtooned(laiob (m,:,:)) forrest@0: end if forrest@0: forrest@0: if (n .eq. 1) then forrest@0: data = ndtooned(laimod(m,:,:)) 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,m) = avg(data(idx)) forrest@0: else forrest@0: yvalues(n,i,m) = yvalues@_FillValue 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,m) = yvalues@_FillValue forrest@0: end if forrest@0: ;############################################################# forrest@0: forrest@0: end do ; m-loop 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: good = ind(.not.ismissing(yvalues(0,:,0))) forrest@0: forrest@0: n_biome = dimsizes(good) forrest@0: forrest@0: ;---------------------------------------------------------------- forrest@0: ; data for tseries plot forrest@0: forrest@0: yvalues_g = new((/data_n,n_biome,ntime/),float) forrest@0: forrest@0: yvalues_g@units = "" forrest@0: forrest@0: yvalues_g = yvalues(:,good,:) forrest@0: forrest@0: delete (good) forrest@0: forrest@0: ;******************************************************************* forrest@0: ; res for line plot forrest@0: ;******************************************************************* forrest@0: ; for x-axis in xyplot forrest@0: mon = ispan(1,12,1) forrest@0: mon@long_name = "month" forrest@0: forrest@0: res = True ; plot mods desired forrest@0: res@xyLineThicknesses = (/2.0,2.0,2.0/) ; make 2nd lines thicker forrest@0: res@xyLineColors = (/"blue","red"/) ; change line color forrest@0: forrest@0: ; res@tiMainFontHeightF = 0.025 ; size of title forrest@0: forrest@0: ; Add a boxed legend using the more simple method forrest@0: res@pmLegendDisplayMode = "Always" forrest@0: ; res@pmLegendWidthF = 0.1 forrest@0: res@pmLegendWidthF = 0.08 forrest@0: res@pmLegendHeightF = 0.06 forrest@0: ; res@pmLegendOrthogonalPosF = -1.17 forrest@0: ; res@pmLegendOrthogonalPosF = -1.00 ;(downward) forrest@0: res@pmLegendOrthogonalPosF = -0.30 ;(downward) forrest@0: forrest@0: ; res@pmLegendParallelPosF = 0.18 forrest@0: res@pmLegendParallelPosF = 0.23 ;(rightward) forrest@0: forrest@0: ; res@lgPerimOn = False forrest@0: res@lgLabelFontHeightF = 0.015 forrest@0: res@xyExplicitLegendLabels = (/"observed",model_name/) forrest@0: forrest@0: ;************************************************************ forrest@0: forrest@0: plot_data = new((/2,ntime/),float) forrest@0: plot_data@long_name = "Leaf Area Index" forrest@0: forrest@0: ;---------------------------------------------- forrest@0: ; time series plot : per biome forrest@0: forrest@0: do m = 0, n_biome-1 forrest@0: forrest@0: plot_name = "annual_biome_"+ m forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: forrest@0: title = "LAI : "+ row_head(m) forrest@0: res@tiMainString = title forrest@0: forrest@0: plot_data(0,:) = yvalues_g(0,m,:) forrest@0: plot_data(1,:) = yvalues_g(1,m,:) forrest@0: forrest@0: plot = gsn_csm_xy(wks,mon,plot_data,res) forrest@0: forrest@0: delete (wks) forrest@0: delete (plot) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: end do forrest@0: forrest@0: ;-------------------------------------------- forrest@0: ; time series plot: all biome forrest@0: forrest@0: plot_name = "annual_biome_"+ n_biome forrest@0: forrest@0: wks = gsn_open_wks (plot_type,plot_name) forrest@0: forrest@0: title = "LAI : "+ row_head(n_biome) forrest@0: res@tiMainString = title forrest@0: forrest@0: do k = 0,ntime-1 forrest@0: plot_data(0,k) = avg(yvalues_g(0,:,k)) forrest@0: plot_data(1,k) = avg(yvalues_g(1,:,k)) forrest@0: end do forrest@0: forrest@0: plot = gsn_csm_xy(wks,mon,plot_data,res) forrest@0: forrest@0: delete (wks) forrest@0: delete (plot) forrest@0: forrest@0: system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \ forrest@0: "rm "+plot_name+"."+plot_type) forrest@0: forrest@0: delete (plot_data) forrest@0: forrest@0: ;*************************************************************************** forrest@0: ; write total score to file forrest@0: ;*************************************************************************** forrest@0: forrest@0: asciiwrite("M_save.lai", M_total) forrest@0: forrest@0: ;*************************************************************************** forrest@0: ; output plots forrest@0: ;*************************************************************************** forrest@0: output_dir = model_name+"/lai" forrest@0: forrest@0: system("mv *.png *.html " + output_dir) forrest@0: ;*************************************************************************** forrest@0: end forrest@0: