1module gridBinning_mod
2 ! MODULE gridBinning_mod (prefix='gbi' category='4. Data Object transformations')
3 !
4 !:Purpose: To compute categorical mean and standard deviation for gridded data
5 ! contained in a gridStateVector or in an ensemble of
6 ! gridStateVectors (e.g. the respective mean over land and sea)
7 !
8 use midasMpi_mod
9 use ensembleStateVector_mod
10 use gridStateVector_mod
11 use gridStateVectorFileIO_mod
12 use utilities_mod
13 use horizontalCoord_mod
14 use timeCoord_mod
15 implicit none
16 save
17 private
18
19 ! public procedures
20 public :: struct_gbi
21 public :: gbi_setup, gbi_deallocate
22 public :: gbi_stdDev, gbi_mean
23
24 type :: struct_gbi
25 character(len=24) :: binningStrategy
26 type(struct_gsv) :: statevector_bin2d
27 integer :: numBins2d
28 end type struct_gbi
29
30 integer, external :: get_max_rss
31
32 ! Control parameter for the level of listing output
33 logical, parameter :: verbose = .true.
34
35 ! module interfaces
36 interface gbi_stdDev
37 module procedure gbi_stdDev_ens
38 end interface gbi_stdDev
39
40 interface gbi_mean
41 module procedure gbi_mean_gsv
42 end interface gbi_mean
43
44contains
45
46 !--------------------------------------------------------------------------
47 ! gbi_setup
48 !--------------------------------------------------------------------------
49 subroutine gbi_setup(gbi, binningStrategy, statevector_template, hco_coregrid, &
50 mpi_distribution_opt, writeBinsToFile_opt)
51 implicit none
52
53 ! Arguments:
54 type(struct_gbi), intent(inout) :: gbi
55 character(len=*), intent(in) :: binningStrategy
56 type(struct_gsv), intent(in) :: statevector_template
57 type(struct_hco), pointer, intent(in) :: hco_coregrid
58 character(len=*), optional, intent(in) :: mpi_distribution_opt
59 logical, optional, intent(in) :: writeBinsToFile_opt
60
61 ! Locals:
62 type(struct_gsv) :: statevector_landSeaTopo
63 integer :: myLonBeg, myLonEnd
64 integer :: myLatBeg, myLatEnd
65 integer :: latIndex, lonIndex
66 real(8) :: binCategory
67 logical :: allocHeightSfc
68 logical :: writeBinsToFile
69 logical :: skip
70 logical :: mpi_local
71 real(4), pointer :: bin2d(:,:,:)
72 real(4), pointer :: data2d(:,:,:)
73 character(len=24) :: mpi_distribution
74
75 !
76 !- 1. Allocate a 2D statevector that will contains the bin category for each
77 ! grid point
78 !
79 if (.not. gsv_isAllocated(statevector_template)) then
80 call utl_abort('gbi_setup: the input template statevector was not allocated')
81 end if
82
83 mpi_local = statevector_template%mpi_local
84 if (present(mpi_distribution_opt)) then
85 mpi_distribution = mpi_distribution_opt
86 if ( mpi_distribution == 'None') then
87 mpi_local = .false.
88 end if
89 else
90 mpi_distribution = 'Tiles'
91 end if
92
93 if (trim(binningStrategy) == 'landSeaTopo') then
94 allocHeightSfc = .true.
95 else
96 allocHeightSfc = .false.
97 end if
98
99 call gsv_allocate(gbi%statevector_bin2d, 1, statevector_template%hco, &
100 statevector_template%vco, varNames_opt=(/'BIN'/), &
101 mpi_distribution_opt=mpi_distribution, dataKind_opt=4, &
102 dateStamp_opt=statevector_template%dateStamp3d, &
103 mpi_local_opt=mpi_local, allocHeightSfc_opt=allocHeightSfc)
104
105 !
106 !- 2. Set the bins
107 !
108 call gsv_getField(gbi%statevector_bin2d,bin2d)
109
110 myLonBeg = gbi%statevector_bin2d%myLonBeg
111 myLonEnd = gbi%statevector_bin2d%myLonEnd
112 myLatBeg = gbi%statevector_bin2d%myLatBeg
113 myLatEnd = gbi%statevector_bin2d%myLatEnd
114
115 select case(trim(binningStrategy))
116 case ('GridPoint')
117 if (mmpi_myid == 0) then
118 write(*,*) 'gbi_setup : BIN_TYPE = No horizontal averaging '
119 write(*,*) ' > One bin per horizontal grid point'
120 end if
121
122 gbi%numBins2d = statevector_template%hco%ni * statevector_template%hco%nj
123 binCategory = 0
124 do latIndex = 1, statevector_template%hco%nj
125 do lonIndex = 1, statevector_template%hco%ni
126 binCategory = binCategory + 1
127 if (lonIndex >= myLonBeg .and. lonIndex <= myLatEnd .and. &
128 latIndex >= myLatBeg .and. latIndex <= myLatEnd ) then
129 bin2d(lonIndex,latIndex,1) = real(binCategory,4)
130 end if
131 end do
132 end do
133
134 case ('YrowBand')
135 if (mmpi_myid == 0) then
136 write(*,*)
137 write(*,*) 'gbi_setup : BIN_TYPE = One bin per Y row'
138 end if
139
140 gbi%numBins2d = statevector_template%hco%nj
141 binCategory = 0
142 do latIndex = 1, statevector_template%hco%nj
143 binCategory = binCategory + 1
144 if (latIndex >= myLatBeg .and. latIndex <= myLatEnd ) then
145 bin2d(:,latIndex,1) = real(binCategory,4)
146 end if
147 end do
148
149 case ('HorizontalMean','FullDomain')
150 if (mmpi_myid == 0) then
151 write(*,*)
152 write(*,*) 'gbi_setup : BIN_TYPE = Average over all horizontal points'
153 end if
154
155 gbi%numBins2d = 1
156 binCategory = 1
157 bin2d(:,:,1) = real(binCategory,4)
158
159 case ('landSeaTopo')
160 if (mmpi_myid == 0) then
161 write(*,*)
162 write(*,*) 'gbi_setup : BIN_TYPE = land/sea binning + storage of the topography field'
163 end if
164
165 call gsv_allocate(statevector_landSeaTopo, 1, statevector_template%hco, &
166 statevector_template%vco, varNames_opt=(/'MG'/), &
167 mpi_distribution_opt=mpi_distribution, dataKind_opt=4, &
168 dateStamp_opt=statevector_template%dateStamp3d, &
169 mpi_local_opt=mpi_local, hInterpolateDegree_opt='LINEAR', &
170 allocHeightSfc_opt=.true. )
171
172 call gio_readFromFile(statevector_landSeaTopo, './dataForGridBinning.fst', ' ', ' ', &
173 readHeightSfc_opt=.true.)
174
175 call gsv_getField(statevector_landSeaTopo,data2d)
176
177 gbi%numBins2d = 2 ! (land or sea)
178 do latIndex= myLatBeg, myLatEnd
179 do lonIndex = myLonBeg, myLonEnd
180
181 skip = .false.
182 if (.not. statevector_template%hco%global) then
183 if (lonIndex > hco_coregrid%ni .or. latIndex > hco_coregrid%nj) then
184 ! We are in the extension zone
185 skip = .true.
186 end if
187 end if
188
189 if (skip) then
190 bin2d(lonIndex,latIndex,1) = -1.0
191 else
192 if (data2D(lonIndex,latIndex,1) >= 0.05) then
193 bin2d(lonIndex,latIndex,1) = 1.0 ! land point
194 else
195 bin2d(lonIndex,latIndex,1) = 2.0 ! sea point
196 end if
197 end if
198
199 end do
200 end do
201
202 ! Store the topography in gbi%statevector_bin2d
203 call gsv_copyHeightSfc(statevector_landSeaTopo, & ! IN
204 gbi%statevector_bin2d ) ! OUT
205
206 call gsv_deallocate(statevector_landSeaTopo)
207
208 case default
209 call utl_abort('gbi_setup : Invalid Binning Strategy : '//trim(BinningStrategy))
210 end select
211
212 !
213 !- 3. Write the bins to a file (if desired)
214 !
215 if ( present(writeBinsToFile_opt) ) then
216 writeBinsToFile = writeBinsToFile_opt
217 else
218 writeBinsToFile = .false.
219 end if
220
221 if (writeBinsToFile .and. (mpi_local .or. mmpi_myid == 0) ) then
222 call gio_writeToFile(gbi%statevector_bin2d, './gridBinning.fst', & ! IN
223 'BINNING') ! IN
224 end if
225
226 end subroutine gbi_setup
227
228 !--------------------------------------------------------------------------
229 ! gbi_deallocate
230 !--------------------------------------------------------------------------
231 subroutine gbi_deallocate(gbi)
232 implicit none
233
234 ! Arguments:
235 type(struct_gbi), intent(inout) :: gbi
236
237 call gsv_deallocate(gbi%statevector_bin2d)
238
239 end subroutine gbi_deallocate
240
241 !--------------------------------------------------------------------------
242 ! gbi_mean_gsv
243 !--------------------------------------------------------------------------
244 subroutine gbi_mean_gsv(gbi, statevector_in, statevector_out)
245 implicit none
246
247 ! Arguments:
248 type(struct_gbi), intent(in) :: gbi
249 type(struct_gsv), intent(inout) :: statevector_in
250 type(struct_gsv), intent(inout) :: statevector_out
251
252 ! Locals:
253 integer :: myBinCount(gbi%numBins2d)
254 integer :: binCount(gbi%numBins2d)
255 real(8) :: myBinSum (gbi%numBins2d)
256 real(8) :: binMean (gbi%numBins2d)
257 real(8), pointer :: field4d(:,:,:,:)
258 real(8), pointer :: mean4d (:,:,:,:)
259 real(4), pointer :: bin2d(:,:,:)
260 integer :: myLonBeg, myLonEnd
261 integer :: myLatBeg, myLatEnd
262 integer :: latIndex, lonIndex, varLevIndex, stepIndex, binIndex
263 integer :: nVarLev, nStep, ier
264
265 if (mmpi_myid == 0 .and. verbose) then
266 write(*,*)
267 write(*,*) 'gbi_mean_gsv: Starting...'
268 end if
269
270 nVarLev = statevector_in%nk
271 nStep = statevector_in%numStep
272 myLonBeg = statevector_in%myLonBeg
273 myLonEnd = statevector_in%myLonEnd
274 myLatBeg = statevector_in%myLatBeg
275 myLatEnd = statevector_in%myLatEnd
276
277 call gsv_getField(statevector_in,field4d)
278 call gsv_getField(statevector_out,mean4d)
279
280 call gsv_getField(gbi%statevector_bin2d,bin2d)
281
282 do stepIndex = 1, nStep
283 do varLevIndex = 1, nVarLev
284
285 !- Sum the values per bin
286 myBinCount(:) = 0
287 myBinSum (:) = 0.d0
288 do latIndex= myLatBeg, myLatEnd
289 do lonIndex = myLonBeg, myLonEnd
290 binIndex = int(bin2d(lonIndex,latIndex,1))
291 if (binIndex /= -1) then
292 myBinCount(binIndex) = myBinCount(binIndex) + 1
293 myBinSum(binIndex) = myBinSum(binIndex) + field4d(lonIndex,latIndex,varLevIndex,stepIndex)
294 end if
295 end do
296 end do
297
298 !- Compute the mean per bin
299 call rpn_comm_allreduce(myBinCount,binCount,gbi%numBins2d,"MPI_INTEGER" ,"MPI_SUM","GRID",ier)
300 do binIndex = 1, gbi%numBins2d
301 binMean(binIndex) = myBinSum(binIndex)
302 call mmpi_allreduce_sumreal8scalar(binMean(binIndex),"GRID")
303 end do
304
305 do binIndex = 1, gbi%numBins2d
306 if (binCount(binIndex) /= 0 ) then
307 binMean(binIndex) = binMean(binIndex) / real(binCount(binIndex),8)
308 else
309 binMean(binIndex) = 0.d0
310 end if
311 end do
312
313 !- Distribute the bin mean values at each grid point
314 do latIndex= myLatBeg, myLatEnd
315 do lonIndex = myLonBeg, myLonEnd
316 binIndex = int(bin2d(lonIndex,latIndex,1))
317 if (binIndex /= -1) then
318 mean4d(lonIndex,latIndex,varLevIndex,stepIndex) = binMean(binIndex)
319 else
320 mean4d(lonIndex,latIndex,varLevIndex,stepIndex) = 0.d0
321 end if
322 end do
323 end do
324
325 end do
326 end do
327
328 if (mmpi_myid == 0 .and. verbose) then
329 write(*,*)
330 write(*,*) 'gbi_mean_gsv: Done!'
331 end if
332
333 end subroutine gbi_mean_gsv
334
335 !--------------------------------------------------------------------------
336 ! gbi_stdDev_ens
337 !--------------------------------------------------------------------------
338 subroutine gbi_stdDev_ens(gbi, ens, statevector)
339 implicit none
340
341 ! Arguments:
342 type(struct_gbi), intent(in) :: gbi
343 type(struct_ens), intent(inout) :: ens
344 type(struct_gsv), intent(inout) :: statevector
345
346 ! Locals:
347 integer :: myBinCount(gbi%numBins2d)
348 integer :: binCount (gbi%numBins2d)
349 real(8) :: myBinSum (gbi%numBins2d)
350 real(8) :: binStdDev (gbi%numBins2d)
351 real(4), pointer :: ptr4d_r4 (:,:,:,:)
352 real(8), pointer :: stdDev_r8(:,:,:,:)
353 real(4), pointer :: bin2d(:,:,:)
354 integer :: myLonBeg, myLonEnd
355 integer :: myLatBeg, myLatEnd
356 integer :: latIndex, lonIndex, varLevIndex, stepIndex, binIndex, memberIndex
357 integer :: nVarLev, nStep, nEns, ier
358 character(len=4), pointer :: varNamesList(:)
359
360 if (mmpi_myid == 0 .and. verbose) then
361 write(*,*)
362 write(*,*) 'gbi_stdDev_ens: Starting...'
363 end if
364
365 if (.not. gsv_isAllocated(statevector)) then
366 nullify(varNamesList)
367 call ens_varNamesList(varNamesList,ens)
368 call gsv_allocate(statevector, ens_getNumStep(ens), &
369 ens_getHco(ens), ens_getVco(ens), varNames_opt=varNamesList, &
370 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
371 dataKind_opt=8 )
372 end if
373
374 call ens_getLatLonBounds(ens, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
375
376 nVarLev = ens_getNumK(ens)
377 nStep = ens_getNumStep(ens)
378 nEns = ens_getNumMembers(ens)
379
380 call gsv_getField(statevector,stdDev_r8)
381 call gsv_getField(gbi%statevector_bin2d,bin2d)
382
383 do varLevIndex = 1, nVarLev
384
385 ptr4d_r4 => ens_getOneLev_r4(ens,varLevIndex)
386
387 do stepIndex = 1, nStep
388
389 !- Sum the square values per bin
390 myBinCount(:) = 0
391 myBinSum (:) = 0.d0
392 do latIndex = myLatBeg, myLatEnd
393 do lonIndex = myLonBeg, myLonEnd
394 binIndex = int(bin2d(lonIndex,latIndex,1))
395 if (binIndex /= -1) then
396 do memberIndex = 1, nEns
397 myBinCount(binIndex) = myBinCount(binIndex) + 1
398 myBinSum(binIndex) = myBinSum(binIndex) + real(ptr4d_r4(memberIndex,stepIndex,lonIndex,latIndex),8)**2
399 end do
400 end if
401 end do
402 end do
403
404 !- Compute the stdDev per bin
405 call rpn_comm_allreduce(myBinCount,binCount ,gbi%numBins2d,"MPI_INTEGER" ,"MPI_SUM","GRID",ier)
406 call rpn_comm_allreduce(myBinSum ,binStdDev,gbi%numBins2d,"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ier)
407
408 do binIndex = 1, gbi%numBins2d
409 if (binCount(binIndex) /= 0 ) then
410 binStdDev(binIndex) = sqrt(binStdDev(binIndex) / real(binCount(binIndex)-1,8))
411 else
412 binStdDev(binIndex) = 0.d0
413 end if
414 end do
415
416 !- Distribute the bin stdDev values at each grid point
417 do latIndex= myLatBeg, myLatEnd
418 do lonIndex = myLonBeg, myLonEnd
419 binIndex = int(bin2d(lonIndex,latIndex,1))
420 if (binIndex /= -1) then
421 stdDev_r8(lonIndex,latIndex,varLevIndex,stepIndex) = binStdDev(binIndex)
422 else
423 stdDev_r8(lonIndex,latIndex,varLevIndex,stepIndex) = 0.d0
424 end if
425 end do
426 end do
427
428 end do
429 end do
430
431 if (mmpi_myid == 0 .and. verbose) then
432 write(*,*)
433 write(*,*) 'gbi_stdDev_ens: Done!'
434 end if
435
436 end subroutine gbi_stdDev_ens
437
438end module gridBinning_mod
