Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
1 ;*************************************************
3 ;************************************************
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
6 ;************************************************
12 ;***********************************************
19 model_name = "10" + model
21 ;************************************************
22 ; read amazon mask data
23 ;************************************************
25 diri = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/"
26 fili = "amazon_mask_"+ model_grid + ".nc"
27 f = addfile (diri+fili,"r")
29 mask_amazon0 = f->mask_amazon
33 ;--------------------------------------------------
34 ; get model data: landfrac and area
36 dirs = "/fis/cgd/cseg/people/jeff/clamp_data/surface_model/"
37 fils = "lnd_"+ model_grid +".nc"
38 fs = addfile (dirs+fils,"r")
40 landfrac = fs->landfrac
45 ; change area from km**2 to m**2
48 ;-----------------------------
49 ; take into account landfrac
51 area0 = area0 * landfrac
53 ;--------------------------------------------
56 dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
57 film = "i01.10"+ model +"_1948-2004_ANN_climo.nc"
58 fm = addfile (dirm+film,"r")
60 if (model .eq. "cn") then
66 datamod0 = data1(0,:,:)
67 datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)
70 if (model .eq. "casa") then
76 datamod0 = data1(0,:,:)
77 datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)
79 ;----------------------------------------------------
85 ;************************************************
87 ;************************************************
88 ob_name = "LC15_Amazon_Biomass"
90 diro = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/"
91 filo = "amazon_biomass_"+model_grid+".nc"
92 fo = addfile (diro+filo,"r")
99 ;************************************************
100 ; Units for these variables are:
105 ; We want to convert these to KgC/m2
106 ; ha = 100m*100m = 10,000 m2
107 ; MgC/ha*1000/10,000 = KgC/m2
108 ; Peta g = 1.e15 g = 1.e12 Kg
110 factor_aboveground = 0.5
112 factor_unit_mod = 0.001
114 dataob = dataob * factor_aboveground * factor_unit_ob
115 datamod0 = datamod0 * factor_unit_mod
117 dataob@units = "KgC/m2"
118 datamod0@units = "KgC/m2"
120 dataob@long_name = "Amazon Biomass"
121 datamod0@long_name = "Amazon Biomass"
123 ;********************************************************
124 ; get subset of datamod that match dataob
129 ind_lonL = ind(xm .eq. xo(0))
130 ind_lonR = ind(xm .eq. xo(nlon-1))
131 ind_latS = ind(ym .eq. yo(0))
132 ind_latN = ind(ym .eq. yo(nlat-1))
135 datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)
138 area(:,:) = area0(ind_latS:ind_latN,ind_lonL:ind_lonR)
141 mask_amazon(:,:) = mask_amazon0(ind_latS:ind_latN,ind_lonL:ind_lonR)
143 ;********************************************************
144 ; sum over amazom_mask area:
146 ; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.)
148 area_sum = sum(area*mask_amazon)
150 area_sum_amazon_ob = sum(dataob*area*mask_amazon)
151 area_sum_amazon_mod = sum(datamod*area*mask_amazon)
153 ; Peta g = 1.e15 g = 1.e12 Kg
155 print (area_sum_amazon_ob)
156 print (area_sum_amazon_mod)
159 ;----------------------------------------------------------------------
163 resg = True ; Use plot options
164 resg@cnFillOn = True ; Turn on color fill
165 resg@gsnSpreadColors = True ; use full colormap
166 ; resg@cnFillMode = "RasterFill" ; Turn on raster color
167 ; resg@lbLabelAutoStride = True
168 resg@cnLinesOn = False ; Turn off contourn lines
169 resg@mpFillOn = False ; Turn off map fill
170 resg@gsnAddCyclic = False
172 resg@gsnSpreadColors = True ; use full colormap
173 resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
174 resg@cnMinLevelValF = 0. ; Min level
175 resg@cnMaxLevelValF = 30. ; Max level
176 resg@cnLevelSpacingF = 2. ; interval
178 resg@mpMinLatF = -21.1 ; range to zoom in on
179 resg@mpMaxLatF = 13.8
180 resg@mpMinLonF = 277.28
181 resg@mpMaxLonF = 326.43
182 ;------------------------------------------------------------------------
185 plot_name = "global_ob"
186 title = ob_name+" "+ model_grid
187 resg@tiMainString = title
189 wks = gsn_open_wks (plot_type,plot_name) ; open workstation
190 gsn_define_colormap(wks,"gui_default") ; choose colormap
192 plot = gsn_csm_contour_map_ce(wks,dataob,resg)