|
1 ;--------------------------------------------------------------- |
|
2 ; convert from 1x1 to T31/T42 |
|
3 ; |
|
4 ; use model T31/T42 lat and lon |
|
5 ; output model area, too. |
|
6 ; |
|
7 ;------------------------------------------------------------- |
|
8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" |
|
9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" |
|
10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" |
|
11 undef("copyatt") |
|
12 load "/fis/cgd/cseg/people/jeff/ncl/remap_areaave.ncl" |
|
13 begin |
|
14 ;************************************************ |
|
15 ; output data |
|
16 ;************************************************ |
|
17 diro = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/" |
|
18 filo = "Fire_C_1997-2006_monthly_1.9.nc" |
|
19 c = addfile(diro+filo,"c") |
|
20 filedimdef(c,"year",-1,True) |
|
21 |
|
22 ;****************************************************** |
|
23 ; input grid: 1deg x 1deg, N->S, 180W-180E |
|
24 ;****************************************************** |
|
25 diri = "/fis/cgd/cseg/people/jeff/fire_data/ob/GFEDv2_C/" |
|
26 fili = "Fire_C_1997-2006_monthly_1x1.nc" |
|
27 g = addfile (diri+fili,"r") |
|
28 bi = g->FIRE_C |
|
29 xi = g->lon |
|
30 yi = g->lat |
|
31 |
|
32 dsizes = dimsizes(bi) |
|
33 nyear = dsizes(0) |
|
34 nmonth = dsizes(1) |
|
35 nlatx = dsizes(2) |
|
36 nlonx = dsizes(3) |
|
37 |
|
38 ;=================================================== |
|
39 ; change from 180W-180E, 90N-90S to 0-360E, 90S-90N |
|
40 ;=================================================== |
|
41 |
|
42 yi = (/ yi(::-1) /) |
|
43 bi = (/ bi(:,:,::-1,:) /) |
|
44 |
|
45 b2 = bi |
|
46 x2 = xi |
|
47 |
|
48 do i= 0,nlonx-1 |
|
49 if (i .lt. 180) then |
|
50 p = i + 180 |
|
51 xi(p) = x2(i) + 360. |
|
52 else |
|
53 p = i - 180 |
|
54 xi(p) = x2(i) |
|
55 end if |
|
56 bi(:,:,:,p)= b2(:,:,:,i) |
|
57 end do |
|
58 |
|
59 bi&lat = yi |
|
60 bi&lon = xi |
|
61 |
|
62 print (xi) |
|
63 print (yi) |
|
64 |
|
65 ;================================================== |
|
66 ; edges |
|
67 ;================================================== |
|
68 |
|
69 xi_edge = fspan(0.0, 360.0, nlonx+1) |
|
70 xi_edge@units = "degrees_east" |
|
71 yi_edge = fspan(-90.0, 90.0, nlatx+1) |
|
72 yi_edge@units = "degrees_north" |
|
73 |
|
74 ;================================================== |
|
75 ; area |
|
76 ;================================================== |
|
77 |
|
78 deg2rad = 4.0*atan(1.0)/180.0 |
|
79 earth_rad = 6.37122e6 |
|
80 dx = earth_rad * deg2rad * (xi_edge(1)-xi_edge(0)) |
|
81 dsin = sin(deg2rad*yi_edge(1:180)) - sin(deg2rad*yi_edge(0:179)) |
|
82 area_in = new((/ nlatx, nlonx /), float) |
|
83 area_in = conform(area_in, earth_rad*dx*dsin,0) |
|
84 area_in@units = "m^2" |
|
85 |
|
86 ;************************************************************** |
|
87 ; output grid: (0-360), (90S-90N) |
|
88 ;************************************************************** |
|
89 ; interpolate into model T31/T42, using model lat and lon |
|
90 |
|
91 diri = "/fis/cgd/cseg/people/jeff/surface_data/" |
|
92 fili = "lnd_1.9.nc" |
|
93 b = addfile(diri+fili,"r") |
|
94 |
|
95 yo = b->lat |
|
96 xo = b->lon |
|
97 nlat = dimsizes(yo) |
|
98 nlon = dimsizes(xo) |
|
99 |
|
100 xo_edge = fspan(0.0, 360.0, nlon+1) - 0.5*360.0/nlon |
|
101 xo_edge@units = "degrees_east" |
|
102 |
|
103 yo_edge = new(nlat+1, float) |
|
104 gau_info = gaus(nlat/2) |
|
105 yo_edge(0) = -90.0 |
|
106 yo_edge(1:nlat-1) = doubletofloat(0.5*(gau_info(0:nlat-2,0)+gau_info(1:nlat-1,0))) |
|
107 yo_edge(nlat) = 90.0 |
|
108 yo_edge@units = "degrees_north" |
|
109 |
|
110 ; yo = 0.5*(yo_edge(0:nlat-1) + yo_edge(1:nlat)) |
|
111 ; yo@units = yo_edge@units |
|
112 ; xo = 0.5*(xo_edge(0:nlon-1) + xo_edge(1:nlon)) |
|
113 ; xo@units = xo_edge@units |
|
114 |
|
115 ;=================================================== |
|
116 ; output |
|
117 ;=================================================== |
|
118 |
|
119 w = new((/nyear,nmonth,nlat,nlon/),float) |
|
120 |
|
121 do m = 0,nyear-1 |
|
122 do n = 0,nmonth-1 |
|
123 w(m,n,:,:) = remap_areaave(xi_edge, yi_edge, bi(m,n,:,:), xo_edge, yo_edge) |
|
124 end do |
|
125 end do |
|
126 |
|
127 w!0 = "year" |
|
128 w!1 = "month" |
|
129 w!2 = "lat" |
|
130 w&lat = yo |
|
131 w!3 = "lon" |
|
132 w&lon = xo |
|
133 w@units = bi@units |
|
134 w@long_name = bi@long_name |
|
135 w@_FillValue = bi@_FillValue |
|
136 |
|
137 c->FIRE_C = w |
|
138 c->date = g->date |
|
139 c->area = b->area |
|
140 |
|
141 end |