forrest@0
|
1 |
;********************************************************
|
forrest@0
|
2 |
;using model biome vlass
|
forrest@0
|
3 |
;
|
forrest@0
|
4 |
; required command line input parameters:
|
forrest@0
|
5 |
; ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
|
forrest@0
|
6 |
;
|
forrest@0
|
7 |
; histogram normalized by rain and compute correleration
|
forrest@0
|
8 |
;**************************************************************
|
forrest@0
|
9 |
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
|
forrest@0
|
10 |
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
|
forrest@0
|
11 |
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
|
forrest@0
|
12 |
;**************************************************************
|
forrest@0
|
13 |
procedure set_line(lines:string,nline:integer,newlines:string)
|
forrest@0
|
14 |
begin
|
forrest@0
|
15 |
; add line to ascci/html file
|
forrest@0
|
16 |
|
forrest@0
|
17 |
nnewlines = dimsizes(newlines)
|
forrest@0
|
18 |
if(nline+nnewlines-1.ge.dimsizes(lines))
|
forrest@0
|
19 |
print("set_line: bad index, not setting anything.")
|
forrest@0
|
20 |
return
|
forrest@0
|
21 |
end if
|
forrest@0
|
22 |
lines(nline:nline+nnewlines-1) = newlines
|
forrest@0
|
23 |
; print ("lines = " + lines(nline:nline+nnewlines-1))
|
forrest@0
|
24 |
nline = nline + nnewlines
|
forrest@0
|
25 |
return
|
forrest@0
|
26 |
end
|
forrest@0
|
27 |
;**************************************************************
|
forrest@0
|
28 |
; Main code.
|
forrest@0
|
29 |
begin
|
forrest@0
|
30 |
|
forrest@0
|
31 |
plot_type = "ps"
|
forrest@0
|
32 |
plot_type_new = "png"
|
forrest@0
|
33 |
|
forrest@0
|
34 |
;---------------------------------------------------
|
forrest@0
|
35 |
; model name and grid
|
forrest@0
|
36 |
|
forrest@0
|
37 |
model_grid = "T42"
|
forrest@0
|
38 |
|
forrest@0
|
39 |
model_name = "cn"
|
forrest@0
|
40 |
model_name1 = "i01.06cn"
|
forrest@0
|
41 |
model_name2 = "i01.10cn"
|
forrest@0
|
42 |
|
forrest@0
|
43 |
;---------------------------------------------------
|
forrest@0
|
44 |
; get biome data: model
|
forrest@0
|
45 |
|
forrest@0
|
46 |
biome_name_mod = "Model PFT Class"
|
forrest@0
|
47 |
|
forrest@0
|
48 |
dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
|
forrest@0
|
49 |
film = "class_pft_"+model_grid+".nc"
|
forrest@0
|
50 |
fm = addfile(dirm+film,"r")
|
forrest@0
|
51 |
|
forrest@0
|
52 |
classmod = fm->CLASS_PFT
|
forrest@0
|
53 |
|
forrest@0
|
54 |
delete (fm)
|
forrest@0
|
55 |
|
forrest@0
|
56 |
; model data has 17 land-type classes
|
forrest@0
|
57 |
|
forrest@0
|
58 |
nclass_mod = 17
|
forrest@0
|
59 |
|
forrest@0
|
60 |
;--------------------------------------------------
|
forrest@0
|
61 |
; get model data: landfrac and area
|
forrest@0
|
62 |
|
forrest@0
|
63 |
dirm = "/fis/cgd/cseg/people/jeff/surface_data/"
|
forrest@0
|
64 |
film = "lnd_T42.nc"
|
forrest@0
|
65 |
fm = addfile (dirm+film,"r")
|
forrest@0
|
66 |
|
forrest@0
|
67 |
landmask = fm->landmask
|
forrest@0
|
68 |
landfrac = fm->landfrac
|
forrest@0
|
69 |
area = fm->area
|
forrest@0
|
70 |
|
forrest@0
|
71 |
delete (fm)
|
forrest@0
|
72 |
|
forrest@0
|
73 |
; change area from km**2 to m**2
|
forrest@0
|
74 |
area = area * 1.e6
|
forrest@0
|
75 |
|
forrest@0
|
76 |
;----------------------------------------------------
|
forrest@0
|
77 |
; read data: time series, model
|
forrest@0
|
78 |
|
forrest@0
|
79 |
dirm = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
|
forrest@0
|
80 |
film = model_name2 + "_Fire_C_1979-2004_monthly.nc"
|
forrest@0
|
81 |
fm = addfile (dirm+film,"r")
|
forrest@0
|
82 |
|
forrest@0
|
83 |
data_mod = fm->COL_FIRE_CLOSS(18:25,:,:,:)
|
forrest@0
|
84 |
|
forrest@0
|
85 |
delete (fm)
|
forrest@0
|
86 |
|
forrest@0
|
87 |
; Units for these variables are:
|
forrest@0
|
88 |
; g C/m^2/s
|
forrest@0
|
89 |
|
forrest@0
|
90 |
; change unit to g C/m^2/month
|
forrest@0
|
91 |
|
forrest@0
|
92 |
nsec_per_month = 60*60*24*30
|
forrest@0
|
93 |
|
forrest@0
|
94 |
data_mod = data_mod * nsec_per_month
|
forrest@0
|
95 |
|
forrest@0
|
96 |
data_mod@unit = "gC/m2/month"
|
forrest@0
|
97 |
;----------------------------------------------------
|
forrest@0
|
98 |
; read data: time series, observed
|
forrest@0
|
99 |
|
forrest@0
|
100 |
dirm = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/"
|
forrest@0
|
101 |
film = "Fire_C_1997-2006_monthly_"+ model_grid+".nc"
|
forrest@0
|
102 |
fm = addfile (dirm+film,"r")
|
forrest@0
|
103 |
|
forrest@0
|
104 |
data_ob = fm->FIRE_C(0:7,:,:,:)
|
forrest@0
|
105 |
|
forrest@0
|
106 |
delete (fm)
|
forrest@0
|
107 |
|
forrest@0
|
108 |
ob_name = "GFEDv2"
|
forrest@0
|
109 |
|
forrest@0
|
110 |
; Units for these variables are:
|
forrest@0
|
111 |
; g C/m^2/month
|
forrest@0
|
112 |
|
forrest@0
|
113 |
;---------------------------------------------------
|
forrest@0
|
114 |
; take into account landfrac
|
forrest@0
|
115 |
|
forrest@0
|
116 |
data_mod = data_mod * conform(data_mod, landfrac, (/2,3/))
|
forrest@0
|
117 |
data_ob = data_ob * conform(data_ob, landfrac, (/2,3/))
|
forrest@0
|
118 |
|
forrest@0
|
119 |
;---------------------------------------------------
|
forrest@0
|
120 |
; get time-mean
|
forrest@0
|
121 |
|
forrest@0
|
122 |
x = dim_avg_Wrap(data_mod(lat|:,lon|:,month|:,year|:))
|
forrest@0
|
123 |
data_mod_m = dim_avg_Wrap( x(lat|:,lon|:,month|:))
|
forrest@0
|
124 |
delete (x)
|
forrest@0
|
125 |
|
forrest@0
|
126 |
x = dim_avg_Wrap( data_ob(lat|:,lon|:,month|:,year|:))
|
forrest@0
|
127 |
data_ob_m = dim_avg_Wrap( x(lat|:,lon|:,month|:))
|
forrest@0
|
128 |
delete (x)
|
forrest@0
|
129 |
|
forrest@0
|
130 |
; printVarSummary(y)
|
forrest@0
|
131 |
|
forrest@0
|
132 |
;----------------------------------------------------
|
forrest@0
|
133 |
; compute correlation coef
|
forrest@0
|
134 |
|
forrest@0
|
135 |
landmask_1d = ndtooned(landmask)
|
forrest@0
|
136 |
data_mod_1d = ndtooned(data_mod_m)
|
forrest@0
|
137 |
data_ob_1d = ndtooned(data_ob_m)
|
forrest@0
|
138 |
|
forrest@0
|
139 |
good = ind(landmask_1d .gt. 0.)
|
forrest@0
|
140 |
; print (dimsizes(good))
|
forrest@0
|
141 |
|
forrest@0
|
142 |
cc = esccr(data_mod_1d(good),data_ob_1d(good),0)
|
forrest@0
|
143 |
; print (cc)
|
forrest@0
|
144 |
|
forrest@0
|
145 |
delete (landmask_1d)
|
forrest@0
|
146 |
delete (data_mod_1d)
|
forrest@0
|
147 |
delete (data_ob_id)
|
forrest@0
|
148 |
|
forrest@0
|
149 |
;----------------------------------------------------
|
forrest@0
|
150 |
; compute M_global
|
forrest@0
|
151 |
|
forrest@0
|
152 |
score_max = 1.
|
forrest@0
|
153 |
|
forrest@0
|
154 |
Mscore = cc * cc * score_max
|
forrest@0
|
155 |
|
forrest@0
|
156 |
M_global = sprintf("%.2f", Mscore)
|
forrest@0
|
157 |
|
forrest@0
|
158 |
;----------------------------------------------------
|
forrest@0
|
159 |
; global res
|
forrest@0
|
160 |
|
forrest@0
|
161 |
resg = True ; Use plot options
|
forrest@0
|
162 |
resg@cnFillOn = True ; Turn on color fill
|
forrest@0
|
163 |
resg@gsnSpreadColors = True ; use full colormap
|
forrest@0
|
164 |
resg@cnLinesOn = False ; Turn off contourn lines
|
forrest@0
|
165 |
resg@mpFillOn = False ; Turn off map fill
|
forrest@0
|
166 |
resg@cnLevelSelectionMode = "ManualLevels" ; Manual contour invtervals
|
forrest@0
|
167 |
|
forrest@0
|
168 |
;----------------------------------------------------
|
forrest@0
|
169 |
; global contour: model vs ob
|
forrest@0
|
170 |
|
forrest@0
|
171 |
plot_name = "global_model_vs_ob"
|
forrest@0
|
172 |
|
forrest@0
|
173 |
wks = gsn_open_wks (plot_type,plot_name)
|
forrest@0
|
174 |
gsn_define_colormap(wks,"gui_default")
|
forrest@0
|
175 |
|
forrest@0
|
176 |
plot=new(3,graphic) ; create graphic array
|
forrest@0
|
177 |
|
forrest@0
|
178 |
resg@gsnFrame = False ; Do not draw plot
|
forrest@0
|
179 |
resg@gsnDraw = False ; Do not advance frame
|
forrest@0
|
180 |
|
forrest@0
|
181 |
;----------------------
|
forrest@0
|
182 |
; plot correlation coef
|
forrest@0
|
183 |
|
forrest@0
|
184 |
gRes = True
|
forrest@0
|
185 |
gRes@txFontHeightF = 0.02
|
forrest@0
|
186 |
gRes@txAngleF = 90
|
forrest@0
|
187 |
|
forrest@0
|
188 |
correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")"
|
forrest@0
|
189 |
|
forrest@0
|
190 |
gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
|
forrest@0
|
191 |
|
forrest@0
|
192 |
;-----------------------
|
forrest@0
|
193 |
; plot ob
|
forrest@0
|
194 |
|
forrest@0
|
195 |
data_ob_m = where(landmask .gt. 0., data_ob_m, data_ob_m@_FillValue)
|
forrest@0
|
196 |
|
forrest@0
|
197 |
title = ob_name
|
forrest@0
|
198 |
resg@tiMainString = title
|
forrest@0
|
199 |
|
forrest@0
|
200 |
resg@cnMinLevelValF = 1.
|
forrest@0
|
201 |
resg@cnMaxLevelValF = 10.
|
forrest@0
|
202 |
resg@cnLevelSpacingF = 1.
|
forrest@0
|
203 |
|
forrest@0
|
204 |
plot(0) = gsn_csm_contour_map_ce(wks,data_ob_m,resg)
|
forrest@0
|
205 |
|
forrest@0
|
206 |
;-----------------------
|
forrest@0
|
207 |
; plot model
|
forrest@0
|
208 |
|
forrest@0
|
209 |
data_mod_m = where(landmask .gt. 0., data_mod_m, data_mod_m@_FillValue)
|
forrest@0
|
210 |
|
forrest@0
|
211 |
title = "Model "+ model_name
|
forrest@0
|
212 |
resg@tiMainString = title
|
forrest@0
|
213 |
|
forrest@0
|
214 |
resg@cnMinLevelValF = 1.
|
forrest@0
|
215 |
resg@cnMaxLevelValF = 10.
|
forrest@0
|
216 |
resg@cnLevelSpacingF = 1.
|
forrest@0
|
217 |
|
forrest@0
|
218 |
plot(1) = gsn_csm_contour_map_ce(wks,data_mod_m,resg)
|
forrest@0
|
219 |
|
forrest@0
|
220 |
;-----------------------
|
forrest@0
|
221 |
; plot model-ob
|
forrest@0
|
222 |
|
forrest@0
|
223 |
resg@cnMinLevelValF = -8.
|
forrest@0
|
224 |
resg@cnMaxLevelValF = 2.
|
forrest@0
|
225 |
resg@cnLevelSpacingF = 1.
|
forrest@0
|
226 |
|
forrest@0
|
227 |
zz = data_ob_m
|
forrest@0
|
228 |
zz = data_mod_m - data_ob_m
|
forrest@0
|
229 |
title = "Model_"+model_name+" - Observed"
|
forrest@0
|
230 |
resg@tiMainString = title
|
forrest@0
|
231 |
|
forrest@0
|
232 |
plot(2) = gsn_csm_contour_map_ce(wks,zz,resg)
|
forrest@0
|
233 |
|
forrest@0
|
234 |
; plot panel
|
forrest@0
|
235 |
|
forrest@0
|
236 |
pres = True ; panel plot mods desired
|
forrest@0
|
237 |
pres@gsnMaximize = True ; fill the page
|
forrest@0
|
238 |
|
forrest@0
|
239 |
gsn_panel(wks,plot,(/3,1/),pres) ; create panel plot
|
forrest@0
|
240 |
|
forrest@0
|
241 |
exit
|
forrest@0
|
242 |
|
forrest@0
|
243 |
; system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new+";"+ \
|
forrest@0
|
244 |
; "rm "+plot_name+"."+plot_type)
|
forrest@0
|
245 |
system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)
|
forrest@0
|
246 |
|
forrest@0
|
247 |
clear (wks)
|
forrest@0
|
248 |
delete (plot)
|
forrest@0
|
249 |
|
forrest@0
|
250 |
delete (data_ob_m)
|
forrest@0
|
251 |
delete (data_mod_m)
|
forrest@0
|
252 |
delete (zz)
|
forrest@0
|
253 |
|
forrest@0
|
254 |
resg@gsnFrame = True ; Do advance frame
|
forrest@0
|
255 |
resg@gsnDraw = True ; Do draw plot
|
forrest@0
|
256 |
|
forrest@0
|
257 |
end
|
forrest@0
|
258 |
|