|
1 ; *********************************************** |
|
2 ; interpolate into model grids (T31) |
|
3 ; *********************************************** |
|
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
|
5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
|
6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" |
|
7 ;************************************************ |
|
8 begin |
|
9 |
|
10 ;************************************************ |
|
11 ; read in model data |
|
12 ;************************************************ |
|
13 ;fili2 = "b30.061n_1995-2004_ANN_climo_lnd.nc" |
|
14 ;model_grid = "T31" |
|
15 ;fili2 = "i01.03cn_1545-1569_ANN_climo.nc" |
|
16 ;model_grid = "T42" |
|
17 fili2 = "newcn05_ncep_1i_ANN_climo_lnd.nc" |
|
18 model_grid = "1.9" |
|
19 |
|
20 diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/" |
|
21 f = addfile (diri2+fili2,"r") |
|
22 |
|
23 lon = f->lon |
|
24 lat = f->lat |
|
25 nlon = dimsizes(lon) |
|
26 nlat = dimsizes(lat) |
|
27 |
|
28 ; print (xi) |
|
29 ; print (yi) |
|
30 ;************************************************ |
|
31 ; output data |
|
32 ;************************************************ |
|
33 diro = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/" |
|
34 filo = "amazon_biomass_"+model_grid+".nc" |
|
35 c = addfile(diro+filo,"c") |
|
36 ;************************************************ |
|
37 ; read in observed data |
|
38 ;************************************************ |
|
39 diri = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/" |
|
40 fili = "amazon_biomass.nc" |
|
41 g = addfile (diri+fili,"r") |
|
42 bi = g->BIOMASS |
|
43 xi = g->lon |
|
44 yi = g->lat |
|
45 ;************************************************ |
|
46 ; from class to value: |
|
47 ; Biomass map has 11 classes: (unit= Mg/ha) |
|
48 ; 0: (_FillValue) |
|
49 ; 1: 0-25 (12.5) |
|
50 ; 2: 25-50 (37.5) |
|
51 ; 3: 50-75 (62.5) |
|
52 ; 4: 75-100 (87.5) |
|
53 ; 5: 100-150 (125) |
|
54 ; 6: 150-200 (175) |
|
55 ; 7: 200-250 (225) |
|
56 ; 8: 250-300 (275) |
|
57 ; 9: 300-350 (325) |
|
58 ; 10: 350-400 (375) |
|
59 ; 11: >400 (425) |
|
60 ;-------------------------- |
|
61 |
|
62 bi@_FillValue = 1.e+36 |
|
63 |
|
64 bi = where(bi.eq.0., bi@_FillValue,bi) |
|
65 bi = where(bi.eq.1., 12.5,bi) |
|
66 bi = where(bi.eq.2., 37.5,bi) |
|
67 bi = where(bi.eq.3., 62.5,bi) |
|
68 bi = where(bi.eq.4., 87.5,bi) |
|
69 bi = where(bi.eq.5., 125.,bi) |
|
70 bi = where(bi.eq.6., 175.,bi) |
|
71 bi = where(bi.eq.7., 225.,bi) |
|
72 bi = where(bi.eq.8., 275.,bi) |
|
73 bi = where(bi.eq.9., 325.,bi) |
|
74 bi = where(bi.eq.10.,375.,bi) |
|
75 bi = where(bi.eq.11.,425.,bi) |
|
76 |
|
77 ;print("min/max = " + min(bi) + "/" + max(bi)) |
|
78 |
|
79 ;************************************************ |
|
80 ; Observed factor_aboveground = 0.5 |
|
81 ;************************************************ |
|
82 |
|
83 factor_aboveground = 0.5 |
|
84 ;bi = bi * factor_aboveground |
|
85 |
|
86 |
|
87 ;************************************************ |
|
88 ; find model grids that is inside observed grids |
|
89 ;************************************************ |
|
90 ind_lonL = min(ind(lon .ge. xi(0))) |
|
91 ind_lonR = max(ind(lon .le. xi(dimsizes(xi)-1))) |
|
92 |
|
93 print (ind_lonL) |
|
94 ;print (xi(0)) |
|
95 ;print (lon(ind_lonL)) |
|
96 |
|
97 print (ind_lonR) |
|
98 ;print (xi(dimsizes(xi)-1)) |
|
99 ;print (lon(ind_lonR)) |
|
100 |
|
101 ind_latS = min(ind(lat .ge. yi(0))) |
|
102 ind_latN = max(ind(lat .le. yi(dimsizes(yi)-1))) |
|
103 |
|
104 print (ind_latS) |
|
105 ;print (yi(0)) |
|
106 ;print (lat(ind_latS)) |
|
107 |
|
108 print (ind_latN) |
|
109 ;print (yi(dimsizes(yi)-1)) |
|
110 ;print (lat(ind_latN)) |
|
111 |
|
112 nlat_out = ind_latN - ind_latS + 1 |
|
113 nlon_out = ind_lonR - ind_lonL + 1 |
|
114 |
|
115 print (nlat_out) |
|
116 print (nlon_out) |
|
117 |
|
118 bo = new((/nlat_out,nlon_out/),float) |
|
119 lat_out = new((/nlat_out/),float) |
|
120 lon_out = new((/nlon_out/),float) |
|
121 |
|
122 do jj=0,nlat_out-1 |
|
123 j = ind_latS + jj |
|
124 lat_out(jj) = lat(j) |
|
125 if (j.eq.0 .or. j.eq.nlat-1) then |
|
126 if (j.eq.0) then |
|
127 LATS = -90. |
|
128 LATN = lat(j)+0.5*(lat(j+1)-lat(j)) |
|
129 end if |
|
130 if (j.eq.nlat-1) then |
|
131 LATS = lat(j)-0.5*(lat(j)-lat(j-1)) |
|
132 LATN = 90. |
|
133 end if |
|
134 else |
|
135 LATS = lat(j)-0.5*(lat(j)-lat(j-1)) |
|
136 LATN = lat(j)+0.5*(lat(j+1)-lat(j)) |
|
137 end if |
|
138 do ii=0,nlon_out-1 |
|
139 i = ind_lonL + ii |
|
140 if (jj.eq.0) then |
|
141 lon_out(ii) = lon(i) |
|
142 end if |
|
143 if (i.eq.0 .or. i.eq.nlon-1) then |
|
144 if (i.eq.0) then |
|
145 LONL = 0. |
|
146 LONR = lon(i)+0.5*(lon(i+1)-lon(i)) |
|
147 end if |
|
148 if (i.eq.nlon-1) then |
|
149 LONL = lon(i)-0.5*(lon(i)-lon(i-1)) |
|
150 LONR = 360. |
|
151 end if |
|
152 else |
|
153 LONL = lon(i)-0.5*(lon(i)-lon(i-1)) |
|
154 LONR = lon(i)+0.5*(lon(i+1)-lon(i)) |
|
155 end if |
|
156 |
|
157 print (LATS) |
|
158 print (LATN) |
|
159 print (LONL) |
|
160 print (LONR) |
|
161 |
|
162 bo(jj,ii) = avg(bi({LATS:LATN},{LONL:LONR})) |
|
163 end do |
|
164 end do |
|
165 |
|
166 lon_out!0 = "lon" |
|
167 lon_out@long_name = "longitude" |
|
168 lon_out@units = "degrees_east" |
|
169 lon_out&lon = lon_out |
|
170 |
|
171 lat_out!0 = "lat" |
|
172 lat_out@long_name = "latitude" |
|
173 lat_out@units = "degrees_north" |
|
174 lat_out&lat = lat_out |
|
175 |
|
176 bo!0 = "lat" |
|
177 bo!1 = "lon" |
|
178 bo&lat = lat_out |
|
179 bo&lon = lon_out |
|
180 bo@units = bi@units |
|
181 bo@long_name = bi@long_name |
|
182 ; bo@_FillValue = bi@_FillValue |
|
183 bo@_FillValue = 1.e+36 |
|
184 |
|
185 c->BIOMASS = bo |
|
186 c->lat = lat_out |
|
187 c->lon = lon_out |
|
188 end |