da_create_bins.inc
References to this file elsewhere.
1 subroutine da_create_bins(ni, nj, nk, bin_type, num_bins, num_bins2d, bin, bin2d, &
2 lat_min, lat_max, binwidth_lat, hgt_min, hgt_max, binwidth_hgt, latitude, height)
3
4 !----------------------------------------------------------------------------
5 !
6 ! Purpose: To create the bins for calculation of statistics.
7 !
8 ! Input:
9 ! ni, nj, nk Dimensions
10 ! bin_type 0: No binning;
11 ! 1: bins for X-direction mean;
12 ! 2: bins for each of binwidth_lat/binwidth_hgt.
13 ! 3: bins for each of binwidth_lat/nk.
14 ! 4: num_hor_bins horizontal bins /nk.
15 ! 5: Average over all horizontal points (nk bins for 3D fields)
16 ! 6: Average over all points (only 1 bin).
17 ! Optional for bin_type = 2:
18 ! lat_min, lat_max Minimum/maximum latitude ranges for bin_type = 2
19 ! binwidth_lat interval between bins for latitude in degree for bin_type = 2
20 ! binwidth_hgt interval between bins for height in meter for bin_type = 2
21 ! num_bins_hgt the number of height bins for bin_type = 2
22 ! latitude 3d field latitude in degree for bin_type = 2
23 ! height 3d field height in meter for bin_type = 2
24 !
25 ! Output:
26 ! num_bins,num_bins2d ---- the number of bins for 3d and 2d fields
27 ! bin , bin2d ---- Assigned bin to a gridpoints for 3d and 2d fields
28 !
29 !----------------------------------------------------------------------------
30
31 implicit none
32
33 integer, intent(in) :: ni, nj, nk ! Dimensions read in.
34 integer, intent(in) :: bin_type ! Type of bin to average over
35 integer, intent(out) :: num_bins ! Number of bins.
36 integer, intent(out) :: num_bins2d ! Number of bins for 2D fields
37 integer, intent(out) :: bin(1:ni,1:nj,1:nk) ! Bin at each 3D point
38 integer, intent(out) :: bin2d(1:ni,1:nj) ! Bin at each 2D point
39
40 real, optional, intent(in) :: lat_min, lat_max ! Used if bin_type = 2 (deg)
41 real, optional, intent(in) :: binwidth_lat ! Used if bin_type = 2 (deg)
42 real, optional, intent(in) :: hgt_min, hgt_max ! Used if bin_type = 2 (deg)
43 real, optional, intent(in) :: binwidth_hgt ! Used if bin_type = 2 (m).
44 real, optional, intent(in) :: latitude(1:ni,1:nj) ! Latitude (degrees).
45 real, optional, intent(in) :: height(1:ni,1:nj,1:nk) ! Height (m).
46
47 integer :: b, i, j, k ! Loop counters.
48 integer :: count ! Counter
49 integer :: num_bins_lat ! Used if bin_type = 2.
50 integer :: num_bins_hgt ! Used if bin_type = 2.
51 integer :: bin_lat ! Latitude bin.
52 integer :: bin_hgt ! Height bin.
53 integer :: num_bins_i, num_bins_j ! Used if bin_type = 4.
54 integer :: nii, njj ! Used if bin_type = 4.
55 integer :: bin_i(1:ni), bin_j(1:nj) ! Used if bin_type = 4.
56 real, allocatable :: binstart_lat(:) ! Used if bin_type = 2 (deg)
57 real, allocatable :: binstart_hgt(:) ! Used if bin_type = 2 (deg)
58
59 if (trace_use) call da_trace_entry("da_create_bins")
60
61 if (bin_type == 0) then ! No averaging in space
62
63 num_bins = nk * nj * ni
64 num_bins2d = nj * ni ! Equals number of horizontal points.
65
66 count = 1
67 do k = 1, nk
68 do j = 1, nj
69 do i = 1, ni
70 bin(i,j,k) = count
71 count = count + 1
72 end do
73 end do
74 end do
75 bin2d(:,:) = bin(:,:,1)
76
77 else if (bin_type == 1) then ! Average over x-direction.
78
79 num_bins = nj * nk
80 num_bins2d = nj
81
82 count = 1
83 do k = 1, nk
84 do j = 1, nj
85 bin(1:ni,j,k) = count
86 count = count + 1
87 end do
88 end do
89 bin2d(:,:) = bin(:,:,1)
90
91 else if (bin_type == 2) then ! Global latitude/height bins:
92
93 ! Setup latitude bins:
94 write(unit=stdout,fmt='(/a,f12.5)')' Minimum latitude = ', lat_min
95 write(unit=stdout,fmt='(a,f12.5)')' Maximum latitude = ', lat_max
96 write(unit=stdout,fmt='(a,f12.5)') &
97 ' Latitude bin width = ', binwidth_lat
98 num_bins_lat = nint((lat_max - lat_min) / binwidth_lat)
99 write(unit=stdout,fmt='(a,i8)') &
100 ' Number of latitude bins = ', num_bins_lat
101
102 allocate(binstart_lat(1:num_bins_lat))
103 do b = 1, num_bins_lat ! Assume south to north (as in WRF).
104 binstart_lat(b) = lat_min + real(b-1) * binwidth_lat
105 end do
106
107 ! Setup height bins:
108 write(unit=stdout,fmt='(/a,f12.5)')' Minimum height = ', hgt_min
109 write(unit=stdout,fmt='(a,f12.5)')' Maximum height = ', hgt_max
110 write(unit=stdout,fmt='(a,f12.5)')' Height bin width = ', binwidth_hgt
111 num_bins_hgt = nint((hgt_max - hgt_min) / binwidth_hgt)
112 write(unit=stdout,fmt='(a,i8)') &
113 ' Number of height bins = ', num_bins_hgt
114
115 allocate(binstart_hgt(1:num_bins_hgt))
116 do b = 1, num_bins_hgt
117 binstart_hgt(b) = hgt_min + real(b-1) * binwidth_hgt
118 end do
119
120 num_bins = num_bins_lat * num_bins_hgt
121 num_bins2d = num_bins_lat
122
123 ! Select height bins:
124 do j = 1, nj
125 do i = 1, ni
126 do k = 1, nk
127 if (height(i,j,k) < binstart_hgt(1)) then
128 bin_hgt = 1 ! In first bin.
129 else if (height(i,j,k) >= binstart_hgt(num_bins_hgt)) then
130 bin_hgt = num_bins_hgt ! In final bin.
131 else
132 do b = 1, num_bins_hgt-1
133 if (height(i,j,k) >= binstart_hgt(b) .and. &
134 height(i,j,k) < binstart_hgt(b+1)) then
135 bin_hgt = b
136 exit
137 end if
138 end do
139 end if
140
141 ! Select latitude bin that point falls in:
142 if (k == 1) then
143 do b = 1, num_bins_lat-1
144 if (latitude(i,j) >= binstart_lat(b) .and. &
145 latitude(i,j) < binstart_lat(b+1)) then
146 bin_lat = b
147 exit
148 end if
149 end do
150 if (latitude(i,j) >= binstart_lat(num_bins_lat)) then
151 ! In final bin.
152 bin_lat = num_bins_lat
153 end if
154 bin2d(i,j) = bin_lat
155 end if
156 bin(i,j,k) = bin_lat + num_bins_lat * (bin_hgt - 1)
157 end do
158 end do
159 end do
160
161 deallocate(binstart_lat)
162 deallocate(binstart_hgt)
163
164 else if (bin_type == 3) then ! Latitude/nk bins:
165
166 ! Setup latitude bins:
167 write(unit=stdout,fmt='(/a,f12.5)')' Minimum latitude = ', lat_min
168 write(unit=stdout,fmt='(a,f12.5)')' Maximum latitude = ', lat_max
169 write(unit=stdout,fmt='(a,f12.5)')' Latitude bin width = ',binwidth_lat
170 num_bins_lat = nint((lat_max - lat_min) / binwidth_lat)
171 write(unit=stdout,fmt='(a,i8)') &
172 ' Number of latitude bins = ', num_bins_lat
173
174 allocate(binstart_lat(1:num_bins_lat))
175 do b = 1, num_bins_lat ! Assume south to north (as in WRF).
176 binstart_lat(b) = lat_min + real(b-1) * binwidth_lat
177 end do
178
179 num_bins = num_bins_lat * nk
180 num_bins2d = num_bins_lat
181
182 ! Select bins:
183 do j = 1, nj
184 do i = 1, ni
185 do k = 1, nk
186 ! Select latitude bin that point falls in:
187 if (k == 1) then
188 do b = 1, num_bins_lat-1
189 if (latitude(i,j) >= binstart_lat(b) .and. &
190 latitude(i,j) < binstart_lat(b+1)) then
191 bin_lat = b
192 exit
193 end if
194 end do
195 if (latitude(i,j) >= binstart_lat(num_bins_lat)) then
196 ! In final bin.
197 bin_lat = num_bins_lat
198 end if
199 bin2d(i,j) = bin_lat
200 end if
201 bin(i,j,k) = bin_lat + num_bins_lat * (k - 1)
202 end do
203 end do
204 end do
205
206 deallocate(binstart_lat)
207
208 else if (bin_type == 4) then ! binwidth_lat/nk bins:
209
210 ! Setup horizontal bins:
211 write(unit=stdout,fmt='(/a,f12.5)') &
212 ' Number of grid-cells to average over = ', binwidth_lat
213 ! use binwidth_lat, but actually an integer number of points.
214
215 num_bins_j = int(real(nj) / real(binwidth_lat))
216 njj = int(binwidth_lat) * num_bins_j
217 do j = 1, njj
218 bin_j(j) = 1 + int(real(j-1) / binwidth_lat)
219 end do
220 if (nj > njj) bin_j(njj+1:nj) = bin_j(njj)
221
222 num_bins_i = int(real(ni) / real(binwidth_lat))
223 nii = int(binwidth_lat) * num_bins_i
224 do i = 1, nii
225 bin_i(i) = 1 + int(real(i-1) / binwidth_lat)
226 end do
227 if (ni > nii) bin_i(nii+1:ni) = bin_i(nii)
228
229 num_bins2d = num_bins_i * num_bins_j
230 num_bins = num_bins2d * nk
231
232 do j = 1, nj
233 do i = 1, ni
234 bin2d(i,j) = bin_i(i) + (bin_j(j) - 1) * num_bins_i
235 do k = 1, nk
236 bin(i,j,k) = bin2d(i,j) + (k - 1) * num_bins2d
237 end do
238 end do
239 end do
240
241 else if (bin_type == 5) then ! Average over all horizontal points.
242
243 num_bins = nk
244 num_bins2d = 1
245
246 do k = 1, nk
247 bin(:,:,k) = k
248 end do
249 bin2d(:,:) = 1
250
251 else if (bin_type == 6) then ! Average over all points.
252
253 num_bins = 1
254 num_bins2d = 1
255 bin(:,:,:) = 1
256 bin2d(:,:) = 1
257 end if
258
259 if (trace_use) call da_trace_exit("da_create_bins")
260
261 end subroutine da_create_bins
262
263