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: 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: ; get model data: landfrac and area forrest@0: forrest@0: film_l = "lnd_"+ model_grid + ".nc" forrest@0: fm_l = addfile (dirs+film_l,"r") forrest@0: forrest@0: landfrac = fm_l->landfrac forrest@0: area0 = fm_l->area forrest@0: forrest@0: delete (fm_l) forrest@0: forrest@0: ; change area from km**2 to m**2 forrest@0: area0 = area0 * 1.e6 forrest@0: forrest@0: ;----------------------------- forrest@0: ; take into account landfrac forrest@0: forrest@0: area0 = area0 * landfrac forrest@0: forrest@0: ;-------------------------------------------- forrest@0: ; read model data forrest@0: forrest@0: fm = addfile (dirm+film1,"r") forrest@0: forrest@0: if (BGC .eq. "cn") then forrest@0: data1 = fm->LIVESTEMC forrest@0: data2 = fm->DEADSTEMC forrest@0: data3 = fm->LEAFC forrest@0: datamod0 = data1(0,:,:) forrest@0: datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:) forrest@0: end if forrest@0: forrest@0: if (BGC .eq. "casa") then forrest@0: factor_WOODC = 0.7 forrest@0: data1 = fm->WOODC forrest@0: data2 = fm->LEAFC forrest@0: datamod0 = data1(0,:,:) forrest@0: datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:) forrest@0: end if forrest@0: forrest@0: ; unit: gC/m2 forrest@0: forrest@0: xm = fm->lon forrest@0: ym = fm->lat forrest@0: forrest@0: delete (fm) forrest@0: forrest@0: ;------------------------------------------------ forrest@0: ; read amazon mask data forrest@0: forrest@0: dir_m = diro + "biomass/" forrest@0: fil_m = "amazon_mask_"+ model_grid + ".nc" forrest@0: fm = addfile (dir_m+fil_m,"r") forrest@0: forrest@0: mask_amazon0 = fm->mask_amazon forrest@0: forrest@0: delete (fm) forrest@0: forrest@0: ;------------------------------------------------ forrest@0: ; read ob data forrest@0: forrest@0: ob_name = "LC15_Amazon_Biomass" forrest@0: forrest@0: dir_b = diro + "biomass/" forrest@0: fil_b = "amazon_biomass_"+model_grid+".nc" forrest@0: fo = addfile (dir_b+fil_b,"r") forrest@0: forrest@0: dataob = fo->BIOMASS forrest@0: xo = fo->lon forrest@0: yo = fo->lat forrest@0: forrest@0: delete (fo) forrest@0: forrest@0: ;************************************************ forrest@0: ; Units for these variables are: forrest@0: ; dataob : MgC/ha forrest@0: ; datamod0 : gC/m2 forrest@0: ; We want to convert these to KgC/m2 forrest@0: ; ha = 100m*100m = 10,000 m2 forrest@0: ; MgC/ha*1000/10,000 = KgC/m2 forrest@0: forrest@0: factor_aboveground = 0.5 forrest@0: factor_unit_ob = 0.1 forrest@0: factor_unit_mod = 0.001 forrest@0: forrest@0: dataob = dataob * factor_aboveground * factor_unit_ob forrest@0: datamod0 = datamod0 * factor_unit_mod forrest@0: forrest@0: dataob@units = "KgC/m2" forrest@0: datamod0@units = "KgC/m2" forrest@0: forrest@0: dataob@long_name = "Amazon Biomass" forrest@0: datamod0@long_name = "Amazon Biomass" forrest@0: ;******************************************************** forrest@0: ; get subset of datamod0 that match dataob forrest@0: forrest@0: nlon = dimsizes(xo) forrest@0: nlat = dimsizes(yo) forrest@0: forrest@0: ind_lonL = ind(xm .eq. xo(0)) forrest@0: ind_lonR = ind(xm .eq. xo(nlon-1)) forrest@0: ind_latS = ind(ym .eq. yo(0)) forrest@0: ind_latN = ind(ym .eq. yo(nlat-1)) forrest@0: forrest@0: datamod = dataob forrest@0: datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR) forrest@0: forrest@0: area = dataob forrest@0: area(:,:) = area0(ind_latS:ind_latN,ind_lonL:ind_lonR) forrest@0: forrest@0: mask_amazon = dataob forrest@0: mask_amazon(:,:) = mask_amazon0(ind_latS:ind_latN,ind_lonL:ind_lonR) forrest@0: forrest@0: mask_amazon@units = "" forrest@0: ;******************************************************** forrest@0: ; sum over amazom_mask area: forrest@0: forrest@0: ; Peta g = 1.e15 g = 1.e12 Kg forrest@0: factor_unit = 1.e-12 forrest@0: forrest@0: ; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.) forrest@0: forrest@0: Sum_area = sum(area*mask_amazon)*factor_unit forrest@0: forrest@0: Sum_ob = sum(dataob*area*mask_amazon)*factor_unit forrest@0: Sum_mod = sum(datamod*area*mask_amazon)*factor_unit forrest@0: forrest@0: avg_ob = Sum_ob /Sum_area forrest@0: avg_mod = Sum_mod/Sum_area forrest@0: forrest@0: Sum_biomass_ob = sprintf("%.2f",Sum_ob ) forrest@0: Sum_biomass_mod = sprintf("%.2f",Sum_mod) forrest@0: forrest@0: ;---------------------------------------------------------------------- forrest@0: ; contour plot 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@gsnAddCyclic = False forrest@0: forrest@0: resg@gsnSpreadColors = True ; use full colormap forrest@0: resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals forrest@0: forrest@0: resg@mpMinLatF = -21.1 ; range to zoom in on forrest@0: resg@mpMaxLatF = 13.8 forrest@0: resg@mpMinLonF = 277.28 forrest@0: resg@mpMaxLonF = 326.43 forrest@0: ;------------------------------------------------------------------------ forrest@0: ; mask plot forrest@0: forrest@0: plot_name = "mask_ob" 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: ;----------------------------------------- forrest@0: ; plot area sum forrest@0: forrest@0: gRes = True forrest@0: gRes@txFontHeightF = 0.02 forrest@0: ; gRes@txAngleF = 90 forrest@0: forrest@0: area_sum_text = "(mask area = "+sprintf("%.2f", Sum_area)+"e12 m2)" forrest@0: forrest@0: gsn_text_ndc(wks,area_sum_text,0.50,0.80,gRes) forrest@0: ;----------------------------------------- forrest@0: forrest@0: resg@cnMinLevelValF = 0. ; Min level forrest@0: resg@cnMaxLevelValF = 1. ; Max level forrest@0: resg@cnLevelSpacingF = 0.1 ; interval forrest@0: forrest@0: resg@tiMainString = "Amazon Mask: grid = "+ model_grid forrest@0: forrest@0: plot = gsn_csm_contour_map_ce(wks,mask_amazon,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: ; contour ob forrest@0: forrest@0: resg@cnMinLevelValF = 0. ; Min level forrest@0: resg@cnMaxLevelValF = 30. ; Max level forrest@0: resg@cnLevelSpacingF = 2. ; interval forrest@0: forrest@0: plot_name = "global_ob" forrest@0: title = ob_name+" "+ model_grid 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,dataob,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: ; contour model forrest@0: forrest@0: resg@cnMinLevelValF = 0. ; Min level forrest@0: resg@cnMaxLevelValF = 30. ; Max level forrest@0: resg@cnLevelSpacingF = 2. ; interval forrest@0: forrest@0: plot_name = "global_model" forrest@0: title = "Model "+ model_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,datamod,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: ; contour model vs ob forrest@0: forrest@0: plot_name = "global_model_vs_ob" 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: delete (plot) 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: ;(d) compute correlation coef and M score forrest@0: forrest@0: uu = ndtooned(datamod) forrest@0: vv = ndtooned(dataob) forrest@0: forrest@0: good = ind(.not.ismissing(uu) .and. .not.ismissing(vv)) forrest@0: forrest@0: ug = uu(good) forrest@0: vg = vv(good) forrest@0: forrest@0: ccrG = esccr(ug,vg,0) forrest@0: forrest@0: score_max = 5.0 forrest@0: forrest@0: ; Miomass = (ccrG*ccrG)* score_max forrest@0: ; new eq forrest@0: bias = sum(abs(ug-vg)/(abs(ug)+abs(vg))) forrest@0: Mbiomass = (1. - (bias/dimsizes(ug)))*score_max forrest@0: M_biomass = sprintf("%.2f", Mbiomass) forrest@0: forrest@0: if (isvar("compare")) then forrest@0: system("sed -e '1,/M_biomass/s/M_biomass/"+M_biomass+"/' "+html_name2+" > "+html_new2+";"+ \ forrest@0: "mv -f "+html_new2+" "+html_name2) forrest@0: end if forrest@0: forrest@0: system("sed s#M_biomass#"+M_biomass+"# "+html_name+" > "+html_new+";"+ \ forrest@0: "mv -f "+html_new+" "+html_name) 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("%5.2f", ccrG)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) forrest@0: ;-------------------------------------------------------------------- forrest@0: forrest@0: ;(a) ob forrest@0: forrest@0: title = ob_name+" "+ model_grid forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(0) = gsn_csm_contour_map_ce(wks,dataob,resg) forrest@0: forrest@0: ;(b) 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,datamod,resg) forrest@0: forrest@0: ;(c) model-ob forrest@0: forrest@0: zz = datamod forrest@0: zz = datamod - dataob forrest@0: title = "Model_"+model_name+" - Observed" forrest@0: forrest@0: resg@cnMinLevelValF = -10. ; Min level forrest@0: resg@cnMaxLevelValF = 10. ; Max level forrest@0: resg@cnLevelSpacingF = 2. ; interval forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) forrest@0: forrest@0: pres = True ; panel plot mods desired forrest@0: ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around forrest@0: ; indiv. plots in panel 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: delete (zz) 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: resg@gsnFrame = True ; draw plot forrest@0: resg@gsnDraw = True ; advance frame forrest@0: ;------------------------------------------------------------------------ forrest@0: ; contour ob : masked forrest@0: forrest@0: resg@cnMinLevelValF = 0. ; Min level forrest@0: resg@cnMaxLevelValF = 30. ; Max level forrest@0: resg@cnLevelSpacingF = 2. ; interval forrest@0: forrest@0: plot_name = "global_mask_ob" forrest@0: title = ob_name+" "+ model_grid 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 average over mask region forrest@0: forrest@0: gRes = True forrest@0: gRes@txFontHeightF = 0.02 forrest@0: gRes@txAngleF = 0 forrest@0: forrest@0: area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_ob)+" Kg C/m2)" forrest@0: forrest@0: gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes) forrest@0: ;----------------------------------------- forrest@0: zo = dataob forrest@0: zo = dataob*mask_amazon forrest@0: zo = where((mask_amazon .le. 0.01), zo@_FillValue, zo) forrest@0: plot = gsn_csm_contour_map_ce(wks,zo,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: ; contour model: masked forrest@0: forrest@0: resg@cnMinLevelValF = 0. ; Min level forrest@0: resg@cnMaxLevelValF = 30. ; Max level forrest@0: resg@cnLevelSpacingF = 2. ; interval forrest@0: forrest@0: plot_name = "global_mask_model" forrest@0: title = "Model "+ model_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 average over mask region forrest@0: forrest@0: gRes = True forrest@0: gRes@txFontHeightF = 0.02 forrest@0: gRes@txAngleF = 0 forrest@0: forrest@0: area_avg_text = "(average over mask area = "+sprintf("%.2f", avg_mod)+" Kg C/m2)" forrest@0: forrest@0: gsn_text_ndc(wks,area_avg_text,0.50,0.81,gRes) forrest@0: ;----------------------------------------- forrest@0: zm = datamod forrest@0: zm = datamod*mask_amazon forrest@0: zm = where((mask_amazon .le. 0.01), zm@_FillValue, zm) forrest@0: plot = gsn_csm_contour_map_ce(wks,zm,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: ; contour model vs ob: masked forrest@0: forrest@0: plot_name = "global_mask_vs_ob" 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=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: ;(d) compute correlation coef and M score forrest@0: forrest@0: delete (good) forrest@0: delete (uu) forrest@0: delete (vv) forrest@0: delete (ug) forrest@0: delete (vg) forrest@0: forrest@0: score_max = 5. forrest@0: forrest@0: uu = ndtooned(zm) forrest@0: vv = ndtooned(zo) forrest@0: forrest@0: good = ind((uu .gt. 0.) .and. (vv .gt. 0.)) forrest@0: forrest@0: ug = uu(good) forrest@0: vg = vv(good) forrest@0: forrest@0: ccrG = esccr(ug,vg,0) forrest@0: forrest@0: ; Miomass = (ccrG*ccrG)*score_max forrest@0: ; new eq forrest@0: bias = sum(abs(ug-vg)/(abs(ug)+abs(vg))) forrest@0: Mbiomass2 = (1. - (bias/dimsizes(ug)))*score_max forrest@0: M_biomask = sprintf("%.2f", Mbiomass2) forrest@0: forrest@0: if (isvar("compare")) then forrest@0: system("sed -e '1,/M_biomask/s/M_biomask/"+M_biomask+"/' "+html_name2+" > "+html_new2+";"+ \ forrest@0: "mv -f "+html_new2+" "+html_name2) forrest@0: system("sed -e '1,/Sum_biomass_ob/s/Sum_biomass_ob/"+Sum_biomass_ob+"/' "+html_name2+" > "+html_new2+";"+ \ forrest@0: "mv -f "+html_new2+" "+html_name2) forrest@0: system("sed -e '1,/Sum_biomass_mod/s/Sum_biomass_mod/"+Sum_biomass_mod+"/' "+html_name2+" > "+html_new2+";"+ \ forrest@0: "mv -f "+html_new2+" "+html_name2) forrest@0: end if forrest@0: forrest@0: system("sed s#M_biomask#"+M_biomask+"# "+html_name+" > "+html_new+";"+ \ forrest@0: "mv -f "+html_new+" "+html_name) forrest@0: system("sed s#Sum_biomass_ob#"+Sum_biomass_ob+"# "+html_name+" > "+html_new+";"+ \ forrest@0: "mv -f "+html_new+" "+html_name) forrest@0: system("sed s#Sum_biomass_mod#"+Sum_biomass_mod+"# "+html_name+" > "+html_new+";"+ \ forrest@0: "mv -f "+html_new+" "+html_name) 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", ccrG)+")" forrest@0: forrest@0: gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes) forrest@0: ;-------------------------------------------------------------------- forrest@0: forrest@0: ;(a) ob forrest@0: forrest@0: title = ob_name+" "+ model_grid forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(0) = gsn_csm_contour_map_ce(wks,zo,resg) forrest@0: forrest@0: ;(b) 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,zm,resg) forrest@0: forrest@0: ;(c) model-ob forrest@0: forrest@0: zz = zo forrest@0: zz = zm - zo forrest@0: title = "Model_"+model_name+" - Observed" forrest@0: forrest@0: resg@cnMinLevelValF = -10. ; Min level forrest@0: resg@cnMaxLevelValF = 10. ; Max level forrest@0: resg@cnLevelSpacingF = 2. ; interval forrest@0: resg@tiMainString = title forrest@0: forrest@0: plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) forrest@0: forrest@0: pres = True ; panel plot mods desired forrest@0: ; pres@gsnPanelYWhiteSpacePercent = 5 ; increase white space around forrest@0: ; indiv. plots in panel 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: ;*************************************************************************** forrest@0: ; add total score and write to file forrest@0: ;*************************************************************************** forrest@0: M_total = Mbiomass + Mbiomass2 forrest@0: forrest@0: asciiwrite("M_save.biomass", M_total) forrest@0: forrest@0: ;*************************************************************************** forrest@0: ; output plots forrest@0: ;*************************************************************************** forrest@0: output_dir = model_name+"/biomass" forrest@0: forrest@0: system("mv *.png " + output_dir) forrest@0: ;*************************************************************************** forrest@0: end