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