forrest@0: ;************************************************* forrest@0: ; ce_1.ncl 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: ;************************************************ forrest@0: begin forrest@0: forrest@0: plot_type = "ps" forrest@0: plot_type_new = "png" forrest@0: forrest@0: ;*********************************************** forrest@0: forrest@0: model = "cn" forrest@0: ;model = "casa" forrest@0: forrest@0: model_grid = "T42" forrest@0: forrest@0: model_name = "10" + model forrest@0: forrest@0: ;************************************************ forrest@0: ; read amazon mask data forrest@0: ;************************************************ forrest@0: forrest@0: diri = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/" forrest@0: fili = "amazon_mask_"+ model_grid + ".nc" forrest@0: f = addfile (diri+fili,"r") forrest@0: forrest@0: mask_amazon0 = f->mask_amazon forrest@0: forrest@0: delete (f) forrest@0: forrest@0: ;-------------------------------------------------- forrest@0: ; get model data: landfrac and area forrest@0: forrest@0: dirs = "/fis/cgd/cseg/people/jeff/clamp_data/surface_model/" forrest@0: fils = "lnd_"+ model_grid +".nc" forrest@0: fs = addfile (dirs+fils,"r") forrest@0: forrest@0: landfrac = fs->landfrac forrest@0: area0 = fs->area forrest@0: forrest@0: delete (fs) 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: dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/" forrest@0: film = "i01.10"+ model +"_1948-2004_ANN_climo.nc" forrest@0: fm = addfile (dirm+film,"r") forrest@0: forrest@0: if (model .eq. "cn") then forrest@0: ; unit: gC/m2 forrest@0: data1 = fm->LIVESTEMC forrest@0: data2 = fm->DEADSTEMC forrest@0: data3 = fm->LEAFC forrest@0: forrest@0: datamod0 = data1(0,:,:) forrest@0: datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:) forrest@0: end if forrest@0: forrest@0: if (model .eq. "casa") then forrest@0: ; unit: gC/m2 forrest@0: ; factor_WOODC = 0.7 forrest@0: data1 = fm->WOODC forrest@0: data2 = fm->LEAFC forrest@0: forrest@0: datamod0 = data1(0,:,:) forrest@0: datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:) forrest@0: end if forrest@0: ;---------------------------------------------------- forrest@0: forrest@0: xm = fm->lon forrest@0: ym = fm->lat forrest@0: forrest@0: delete (fm) forrest@0: ;************************************************ forrest@0: ; read observed data forrest@0: ;************************************************ forrest@0: ob_name = "LC15_Amazon_Biomass" forrest@0: forrest@0: diro = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/" forrest@0: filo = "amazon_biomass_"+model_grid+".nc" forrest@0: fo = addfile (diro+filo,"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: ; Units for these variables are: forrest@0: forrest@0: ; dataob : MgC/ha forrest@0: ; datamod0 : gC/m2 forrest@0: 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: ; Peta g = 1.e15 g = 1.e12 Kg 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: ;******************************************************** forrest@0: ; get subset of datamod 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: ;******************************************************** forrest@0: ; sum over amazom_mask area: forrest@0: forrest@0: ; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.) forrest@0: forrest@0: area_sum = sum(area*mask_amazon) forrest@0: forrest@0: area_sum_amazon_ob = sum(dataob*area*mask_amazon) forrest@0: area_sum_amazon_mod = sum(datamod*area*mask_amazon) forrest@0: forrest@0: ; Peta g = 1.e15 g = 1.e12 Kg forrest@0: forrest@0: print (area_sum_amazon_ob) forrest@0: print (area_sum_amazon_mod) forrest@0: print (area_sum) forrest@0: exit forrest@0: ;---------------------------------------------------------------------- forrest@0: ; global contour forrest@0: forrest@0: ;res 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@cnFillMode = "RasterFill" ; Turn on raster color forrest@0: ; resg@lbLabelAutoStride = True 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: resg@cnMinLevelValF = 0. ; Min level forrest@0: resg@cnMaxLevelValF = 30. ; Max level forrest@0: resg@cnLevelSpacingF = 2. ; interval 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: ;global contour ob 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: frame(wks) forrest@0: end