|
1 ;************************************************* |
|
2 ; ce_1.ncl |
|
3 ;************************************************ |
|
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
|
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
|
6 ;************************************************ |
|
7 begin |
|
8 |
|
9 plot_type = "ps" |
|
10 plot_type_new = "png" |
|
11 |
|
12 ;*********************************************** |
|
13 |
|
14 model = "cn" |
|
15 ;model = "casa" |
|
16 |
|
17 model_grid = "T42" |
|
18 |
|
19 model_name = "10" + model |
|
20 |
|
21 ;************************************************ |
|
22 ; read amazon mask data |
|
23 ;************************************************ |
|
24 |
|
25 diri = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/" |
|
26 fili = "amazon_mask_"+ model_grid + ".nc" |
|
27 f = addfile (diri+fili,"r") |
|
28 |
|
29 mask_amazon0 = f->mask_amazon |
|
30 |
|
31 delete (f) |
|
32 |
|
33 ;-------------------------------------------------- |
|
34 ; get model data: landfrac and area |
|
35 |
|
36 dirs = "/fis/cgd/cseg/people/jeff/clamp_data/surface_model/" |
|
37 fils = "lnd_"+ model_grid +".nc" |
|
38 fs = addfile (dirs+fils,"r") |
|
39 |
|
40 landfrac = fs->landfrac |
|
41 area0 = fs->area |
|
42 |
|
43 delete (fs) |
|
44 |
|
45 ; change area from km**2 to m**2 |
|
46 area0 = area0 * 1.e6 |
|
47 |
|
48 ;----------------------------- |
|
49 ; take into account landfrac |
|
50 |
|
51 area0 = area0 * landfrac |
|
52 |
|
53 ;-------------------------------------------- |
|
54 ; read model data |
|
55 |
|
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") |
|
59 |
|
60 if (model .eq. "cn") then |
|
61 ; unit: gC/m2 |
|
62 data1 = fm->LIVESTEMC |
|
63 data2 = fm->DEADSTEMC |
|
64 data3 = fm->LEAFC |
|
65 |
|
66 datamod0 = data1(0,:,:) |
|
67 datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:) |
|
68 end if |
|
69 |
|
70 if (model .eq. "casa") then |
|
71 ; unit: gC/m2 |
|
72 ; factor_WOODC = 0.7 |
|
73 data1 = fm->WOODC |
|
74 data2 = fm->LEAFC |
|
75 |
|
76 datamod0 = data1(0,:,:) |
|
77 datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:) |
|
78 end if |
|
79 ;---------------------------------------------------- |
|
80 |
|
81 xm = fm->lon |
|
82 ym = fm->lat |
|
83 |
|
84 delete (fm) |
|
85 ;************************************************ |
|
86 ; read observed data |
|
87 ;************************************************ |
|
88 ob_name = "LC15_Amazon_Biomass" |
|
89 |
|
90 diro = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/" |
|
91 filo = "amazon_biomass_"+model_grid+".nc" |
|
92 fo = addfile (diro+filo,"r") |
|
93 |
|
94 dataob = fo->BIOMASS |
|
95 xo = fo->lon |
|
96 yo = fo->lat |
|
97 |
|
98 delete (fo) |
|
99 ;************************************************ |
|
100 ; Units for these variables are: |
|
101 |
|
102 ; dataob : MgC/ha |
|
103 ; datamod0 : gC/m2 |
|
104 |
|
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 |
|
109 |
|
110 factor_aboveground = 0.5 |
|
111 factor_unit_ob = 0.1 |
|
112 factor_unit_mod = 0.001 |
|
113 |
|
114 dataob = dataob * factor_aboveground * factor_unit_ob |
|
115 datamod0 = datamod0 * factor_unit_mod |
|
116 |
|
117 dataob@units = "KgC/m2" |
|
118 datamod0@units = "KgC/m2" |
|
119 |
|
120 dataob@long_name = "Amazon Biomass" |
|
121 datamod0@long_name = "Amazon Biomass" |
|
122 |
|
123 ;******************************************************** |
|
124 ; get subset of datamod that match dataob |
|
125 |
|
126 nlon = dimsizes(xo) |
|
127 nlat = dimsizes(yo) |
|
128 |
|
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)) |
|
133 |
|
134 datamod = dataob |
|
135 datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR) |
|
136 |
|
137 area = dataob |
|
138 area(:,:) = area0(ind_latS:ind_latN,ind_lonL:ind_lonR) |
|
139 |
|
140 mask_amazon = dataob |
|
141 mask_amazon(:,:) = mask_amazon0(ind_latS:ind_latN,ind_lonL:ind_lonR) |
|
142 |
|
143 ;******************************************************** |
|
144 ; sum over amazom_mask area: |
|
145 |
|
146 ; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.) |
|
147 |
|
148 area_sum = sum(area*mask_amazon) |
|
149 |
|
150 area_sum_amazon_ob = sum(dataob*area*mask_amazon) |
|
151 area_sum_amazon_mod = sum(datamod*area*mask_amazon) |
|
152 |
|
153 ; Peta g = 1.e15 g = 1.e12 Kg |
|
154 |
|
155 print (area_sum_amazon_ob) |
|
156 print (area_sum_amazon_mod) |
|
157 print (area_sum) |
|
158 exit |
|
159 ;---------------------------------------------------------------------- |
|
160 ; global contour |
|
161 |
|
162 ;res |
|
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 |
|
171 |
|
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 |
|
177 |
|
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 ;------------------------------------------------------------------------ |
|
183 ;global contour ob |
|
184 |
|
185 plot_name = "global_ob" |
|
186 title = ob_name+" "+ model_grid |
|
187 resg@tiMainString = title |
|
188 |
|
189 wks = gsn_open_wks (plot_type,plot_name) ; open workstation |
|
190 gsn_define_colormap(wks,"gui_default") ; choose colormap |
|
191 |
|
192 plot = gsn_csm_contour_map_ce(wks,dataob,resg) |
|
193 frame(wks) |
|
194 end |