|
forrest@0
|
1 |
;*************************************************
|
|
forrest@0
|
2 |
; ce_1.ncl
|
|
forrest@0
|
3 |
;************************************************
|
|
forrest@0
|
4 |
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
|
|
forrest@0
|
5 |
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
|
|
forrest@0
|
6 |
;************************************************
|
|
forrest@0
|
7 |
begin
|
|
forrest@0
|
8 |
|
|
forrest@0
|
9 |
plot_type = "ps"
|
|
forrest@0
|
10 |
plot_type_new = "png"
|
|
forrest@0
|
11 |
|
|
forrest@0
|
12 |
;***********************************************
|
|
forrest@0
|
13 |
|
|
forrest@0
|
14 |
model = "cn"
|
|
forrest@0
|
15 |
;model = "casa"
|
|
forrest@0
|
16 |
|
|
forrest@0
|
17 |
model_grid = "T42"
|
|
forrest@0
|
18 |
|
|
forrest@0
|
19 |
model_name = "10" + model
|
|
forrest@0
|
20 |
|
|
forrest@0
|
21 |
;************************************************
|
|
forrest@0
|
22 |
; read amazon mask data
|
|
forrest@0
|
23 |
;************************************************
|
|
forrest@0
|
24 |
|
|
forrest@0
|
25 |
diri = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/"
|
|
forrest@0
|
26 |
fili = "amazon_mask_"+ model_grid + ".nc"
|
|
forrest@0
|
27 |
f = addfile (diri+fili,"r")
|
|
forrest@0
|
28 |
|
|
forrest@0
|
29 |
mask_amazon0 = f->mask_amazon
|
|
forrest@0
|
30 |
|
|
forrest@0
|
31 |
delete (f)
|
|
forrest@0
|
32 |
|
|
forrest@0
|
33 |
;--------------------------------------------------
|
|
forrest@0
|
34 |
; get model data: landfrac and area
|
|
forrest@0
|
35 |
|
|
forrest@0
|
36 |
dirs = "/fis/cgd/cseg/people/jeff/clamp_data/surface_model/"
|
|
forrest@0
|
37 |
fils = "lnd_"+ model_grid +".nc"
|
|
forrest@0
|
38 |
fs = addfile (dirs+fils,"r")
|
|
forrest@0
|
39 |
|
|
forrest@0
|
40 |
landfrac = fs->landfrac
|
|
forrest@0
|
41 |
area0 = fs->area
|
|
forrest@0
|
42 |
|
|
forrest@0
|
43 |
delete (fs)
|
|
forrest@0
|
44 |
|
|
forrest@0
|
45 |
; change area from km**2 to m**2
|
|
forrest@0
|
46 |
area0 = area0 * 1.e6
|
|
forrest@0
|
47 |
|
|
forrest@0
|
48 |
;-----------------------------
|
|
forrest@0
|
49 |
; take into account landfrac
|
|
forrest@0
|
50 |
|
|
forrest@0
|
51 |
area0 = area0 * landfrac
|
|
forrest@0
|
52 |
|
|
forrest@0
|
53 |
;--------------------------------------------
|
|
forrest@0
|
54 |
; read model data
|
|
forrest@0
|
55 |
|
|
forrest@0
|
56 |
dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
|
|
forrest@0
|
57 |
film = "i01.10"+ model +"_1948-2004_ANN_climo.nc"
|
|
forrest@0
|
58 |
fm = addfile (dirm+film,"r")
|
|
forrest@0
|
59 |
|
|
forrest@0
|
60 |
if (model .eq. "cn") then
|
|
forrest@0
|
61 |
; unit: gC/m2
|
|
forrest@0
|
62 |
data1 = fm->LIVESTEMC
|
|
forrest@0
|
63 |
data2 = fm->DEADSTEMC
|
|
forrest@0
|
64 |
data3 = fm->LEAFC
|
|
forrest@0
|
65 |
|
|
forrest@0
|
66 |
datamod0 = data1(0,:,:)
|
|
forrest@0
|
67 |
datamod0 = data1(0,:,:) + data2(0,:,:) + data3(0,:,:)
|
|
forrest@0
|
68 |
end if
|
|
forrest@0
|
69 |
|
|
forrest@0
|
70 |
if (model .eq. "casa") then
|
|
forrest@0
|
71 |
; unit: gC/m2
|
|
forrest@0
|
72 |
; factor_WOODC = 0.7
|
|
forrest@0
|
73 |
data1 = fm->WOODC
|
|
forrest@0
|
74 |
data2 = fm->LEAFC
|
|
forrest@0
|
75 |
|
|
forrest@0
|
76 |
datamod0 = data1(0,:,:)
|
|
forrest@0
|
77 |
datamod0 = data1(0,:,:)*factor_WOODC + data2(0,:,:)
|
|
forrest@0
|
78 |
end if
|
|
forrest@0
|
79 |
;----------------------------------------------------
|
|
forrest@0
|
80 |
|
|
forrest@0
|
81 |
xm = fm->lon
|
|
forrest@0
|
82 |
ym = fm->lat
|
|
forrest@0
|
83 |
|
|
forrest@0
|
84 |
delete (fm)
|
|
forrest@0
|
85 |
;************************************************
|
|
forrest@0
|
86 |
; read observed data
|
|
forrest@0
|
87 |
;************************************************
|
|
forrest@0
|
88 |
ob_name = "LC15_Amazon_Biomass"
|
|
forrest@0
|
89 |
|
|
forrest@0
|
90 |
diro = "/fis/cgd/cseg/people/jeff/clamp_data/observed/biomass/"
|
|
forrest@0
|
91 |
filo = "amazon_biomass_"+model_grid+".nc"
|
|
forrest@0
|
92 |
fo = addfile (diro+filo,"r")
|
|
forrest@0
|
93 |
|
|
forrest@0
|
94 |
dataob = fo->BIOMASS
|
|
forrest@0
|
95 |
xo = fo->lon
|
|
forrest@0
|
96 |
yo = fo->lat
|
|
forrest@0
|
97 |
|
|
forrest@0
|
98 |
delete (fo)
|
|
forrest@0
|
99 |
;************************************************
|
|
forrest@0
|
100 |
; Units for these variables are:
|
|
forrest@0
|
101 |
|
|
forrest@0
|
102 |
; dataob : MgC/ha
|
|
forrest@0
|
103 |
; datamod0 : gC/m2
|
|
forrest@0
|
104 |
|
|
forrest@0
|
105 |
; We want to convert these to KgC/m2
|
|
forrest@0
|
106 |
; ha = 100m*100m = 10,000 m2
|
|
forrest@0
|
107 |
; MgC/ha*1000/10,000 = KgC/m2
|
|
forrest@0
|
108 |
; Peta g = 1.e15 g = 1.e12 Kg
|
|
forrest@0
|
109 |
|
|
forrest@0
|
110 |
factor_aboveground = 0.5
|
|
forrest@0
|
111 |
factor_unit_ob = 0.1
|
|
forrest@0
|
112 |
factor_unit_mod = 0.001
|
|
forrest@0
|
113 |
|
|
forrest@0
|
114 |
dataob = dataob * factor_aboveground * factor_unit_ob
|
|
forrest@0
|
115 |
datamod0 = datamod0 * factor_unit_mod
|
|
forrest@0
|
116 |
|
|
forrest@0
|
117 |
dataob@units = "KgC/m2"
|
|
forrest@0
|
118 |
datamod0@units = "KgC/m2"
|
|
forrest@0
|
119 |
|
|
forrest@0
|
120 |
dataob@long_name = "Amazon Biomass"
|
|
forrest@0
|
121 |
datamod0@long_name = "Amazon Biomass"
|
|
forrest@0
|
122 |
|
|
forrest@0
|
123 |
;********************************************************
|
|
forrest@0
|
124 |
; get subset of datamod that match dataob
|
|
forrest@0
|
125 |
|
|
forrest@0
|
126 |
nlon = dimsizes(xo)
|
|
forrest@0
|
127 |
nlat = dimsizes(yo)
|
|
forrest@0
|
128 |
|
|
forrest@0
|
129 |
ind_lonL = ind(xm .eq. xo(0))
|
|
forrest@0
|
130 |
ind_lonR = ind(xm .eq. xo(nlon-1))
|
|
forrest@0
|
131 |
ind_latS = ind(ym .eq. yo(0))
|
|
forrest@0
|
132 |
ind_latN = ind(ym .eq. yo(nlat-1))
|
|
forrest@0
|
133 |
|
|
forrest@0
|
134 |
datamod = dataob
|
|
forrest@0
|
135 |
datamod(:,:) = datamod0(ind_latS:ind_latN,ind_lonL:ind_lonR)
|
|
forrest@0
|
136 |
|
|
forrest@0
|
137 |
area = dataob
|
|
forrest@0
|
138 |
area(:,:) = area0(ind_latS:ind_latN,ind_lonL:ind_lonR)
|
|
forrest@0
|
139 |
|
|
forrest@0
|
140 |
mask_amazon = dataob
|
|
forrest@0
|
141 |
mask_amazon(:,:) = mask_amazon0(ind_latS:ind_latN,ind_lonL:ind_lonR)
|
|
forrest@0
|
142 |
|
|
forrest@0
|
143 |
;********************************************************
|
|
forrest@0
|
144 |
; sum over amazom_mask area:
|
|
forrest@0
|
145 |
|
|
forrest@0
|
146 |
; mask_amazon = where(mask_amazon .ge. 0.5, mask_amazon ,0.)
|
|
forrest@0
|
147 |
|
|
forrest@0
|
148 |
area_sum = sum(area*mask_amazon)
|
|
forrest@0
|
149 |
|
|
forrest@0
|
150 |
area_sum_amazon_ob = sum(dataob*area*mask_amazon)
|
|
forrest@0
|
151 |
area_sum_amazon_mod = sum(datamod*area*mask_amazon)
|
|
forrest@0
|
152 |
|
|
forrest@0
|
153 |
; Peta g = 1.e15 g = 1.e12 Kg
|
|
forrest@0
|
154 |
|
|
forrest@0
|
155 |
print (area_sum_amazon_ob)
|
|
forrest@0
|
156 |
print (area_sum_amazon_mod)
|
|
forrest@0
|
157 |
print (area_sum)
|
|
forrest@0
|
158 |
exit
|
|
forrest@0
|
159 |
;----------------------------------------------------------------------
|
|
forrest@0
|
160 |
; global contour
|
|
forrest@0
|
161 |
|
|
forrest@0
|
162 |
;res
|
|
forrest@0
|
163 |
resg = True ; Use plot options
|
|
forrest@0
|
164 |
resg@cnFillOn = True ; Turn on color fill
|
|
forrest@0
|
165 |
resg@gsnSpreadColors = True ; use full colormap
|
|
forrest@0
|
166 |
; resg@cnFillMode = "RasterFill" ; Turn on raster color
|
|
forrest@0
|
167 |
; resg@lbLabelAutoStride = True
|
|
forrest@0
|
168 |
resg@cnLinesOn = False ; Turn off contourn lines
|
|
forrest@0
|
169 |
resg@mpFillOn = False ; Turn off map fill
|
|
forrest@0
|
170 |
resg@gsnAddCyclic = False
|
|
forrest@0
|
171 |
|
|
forrest@0
|
172 |
resg@gsnSpreadColors = True ; use full colormap
|
|
forrest@0
|
173 |
resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
|
|
forrest@0
|
174 |
resg@cnMinLevelValF = 0. ; Min level
|
|
forrest@0
|
175 |
resg@cnMaxLevelValF = 30. ; Max level
|
|
forrest@0
|
176 |
resg@cnLevelSpacingF = 2. ; interval
|
|
forrest@0
|
177 |
|
|
forrest@0
|
178 |
resg@mpMinLatF = -21.1 ; range to zoom in on
|
|
forrest@0
|
179 |
resg@mpMaxLatF = 13.8
|
|
forrest@0
|
180 |
resg@mpMinLonF = 277.28
|
|
forrest@0
|
181 |
resg@mpMaxLonF = 326.43
|
|
forrest@0
|
182 |
;------------------------------------------------------------------------
|
|
forrest@0
|
183 |
;global contour ob
|
|
forrest@0
|
184 |
|
|
forrest@0
|
185 |
plot_name = "global_ob"
|
|
forrest@0
|
186 |
title = ob_name+" "+ model_grid
|
|
forrest@0
|
187 |
resg@tiMainString = title
|
|
forrest@0
|
188 |
|
|
forrest@0
|
189 |
wks = gsn_open_wks (plot_type,plot_name) ; open workstation
|
|
forrest@0
|
190 |
gsn_define_colormap(wks,"gui_default") ; choose colormap
|
|
forrest@0
|
191 |
|
|
forrest@0
|
192 |
plot = gsn_csm_contour_map_ce(wks,dataob,resg)
|
|
forrest@0
|
193 |
frame(wks)
|
|
forrest@0
|
194 |
end |