ROMS
Loading...
Searching...
No Matches
stats.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE stats_mod
3!
4!git $Id$
5!================================================== Hernan G. Arango ===
6! Copyright (c) 2002-2025 The ROMS Group !
7! Licensed under a MIT/X style license !
8! See License_ROMS.md !
9!=======================================================================
10! !
11! This module contains several routines to compute the statistics of !
12! provided field array, F. The following information is computed: !
13! !
14! S % checksum bit sum hash !
15! S % count processed values count !
16! S % min minimum value !
17! S % max maximum value !
18! S % avg arithmetic mean !
19! S % rms root mean square !
20! !
21! Notice that to compute high order moments (standard deviation, !
22! variance, skewness, and kurosis) cannot be computed in parallel !
23! with a single call because we need to know the mean first. !
24! !
25! Routines: !
26! !
27! stats_2dfld Statistic information for 2D state field !
28! stats_3dfld Statistic information for 3D state field !
29! stats_4dfld Statistic information for 4D state field !
30! !
31!=======================================================================
32!
33 USE mod_param
34!
35 implicit none
36!
37 PUBLIC :: stats_2dfld
38 PUBLIC :: stats_3dfld
39 PUBLIC :: stats_4dfld
40!
41 CONTAINS
42!
43 SUBROUTINE stats_2dfld (ng, tile, model, gtype, S, &
44 & Extract_Flag, &
45 & LBi, UBi, LBj, UBj, &
46 & F, Fmask, debug)
47!
48!=======================================================================
49! !
50! This routine computes requested 2D-field statistics. !
51! !
52! On Input: !
53! !
54! ng Nested grid number (integer) !
55! tile Domain partition (integer) !
56! model Calling model identifier (integer) !
57! gtype Grid type (integer) !
58! Extract_Flag Extraction flag interpolation/decimation (integer) !
59! LBi I-dimension Lower bound (integer) !
60! UBi I-dimension Upper bound (integer) !
61! LBj J-dimension Lower bound (integer) !
62! UBj J-dimension Upper bound (integer) !
63! F 2D state field (2D tiled array) !
64! Fmask Land/sea mask (OPTIONAL, 2D tiled array) !
65! debug Switch to debug (OPTIONAL, logical) !
66! !
67! On Ouput: !
68! !
69! S Field statistics (structure) !
70! !
71!=======================================================================
72!
73 USE mod_param
74 USE mod_parallel
75 USE mod_iounits
76 USE mod_ncparam
77
78#ifdef DISTRIBUTE
79!
80 USE distribute_mod, ONLY : mp_reduce
81#endif
82 USE get_hash_mod, ONLY : get_hash
83!
84! Imported variable declarations.
85!
86 logical, intent(in), optional :: debug
87!
88 integer, intent(in) :: ng, tile, model, gtype
89 integer, intent(in) :: extract_flag
90 integer, intent(in) :: lbi, ubi, lbj, ubj
91!
92#ifdef ASSUMED_SHAPE
93 real(r8), intent(in) :: f(lbi:,lbj:)
94 real(r8), intent(in), optional :: fmask(lbi:,lbj:)
95#else
96 real(r8), intent(in) :: f(lbi:ubi,lbj:ubj)
97 real(r8), intent(in), optional :: fmask(lbi:ubi,lbj:ubj)
98#endif
99 TYPE(t_stats), intent(inout) :: s
100!
101! Local variable declarations.
102!
103 logical :: lprint
104!
105 integer :: imin, imax, jmin, jmax, npts, nsub
106 integer :: i, j, my_count
107
108 integer :: my_threadnum
109!
110 real(r8) :: fac
111 real(r8) :: my_max, my_min
112 real(r8) :: my_avg, my_rms
113
114 real(r8), allocatable :: cwrk(:)
115
116#ifdef DISTRIBUTE
117 real(r8), dimension(5) :: rbuffer
118 character (len=3), dimension(5) :: op_handle
119#endif
120!
121!-----------------------------------------------------------------------
122! Set tile indices.
123!-----------------------------------------------------------------------
124!
125! Determine I- and J-indices according to staggered C-grid type.
126#ifdef GRID_EXTRACT
127!
128! Here, the ELSE statements (Extract_Flag < 0) indicate the processing
129! of data already on the extracted grid, like geometry variables
130! requiring no extraction since they are read during initialization.
131#endif
132!
133 SELECT CASE (abs(gtype))
134 CASE (p2dvar)
135 IF (extract_flag.ge.0) THEN
136 imin=bounds(ng)%IstrP(tile)
137 imax=bounds(ng)%IendP(tile)
138 jmin=bounds(ng)%JstrP(tile)
139 jmax=bounds(ng)%JendP(tile)
140#ifdef GRID_EXTRACT
141 ELSE ! Already on extract grid
142 imin=xtr_bounds(ng)%IstrP(tile)
143 imax=xtr_bounds(ng)%IendP(tile)
144 jmin=xtr_bounds(ng)%JstrP(tile)
145 jmax=xtr_bounds(ng)%JendP(tile)
146#endif
147 END IF
148 CASE (r2dvar)
149 IF (extract_flag.ge.0) THEN
150 imin=bounds(ng)%IstrT(tile)
151 imax=bounds(ng)%IendT(tile)
152 jmin=bounds(ng)%JstrT(tile)
153 jmax=bounds(ng)%JendT(tile)
154#ifdef GRID_EXTRACT
155 ELSE ! Already on extract grid
156 imin=xtr_bounds(ng)%IstrT(tile)
157 imax=xtr_bounds(ng)%IendT(tile)
158 jmin=xtr_bounds(ng)%JstrT(tile)
159 jmax=xtr_bounds(ng)%JendT(tile)
160#endif
161 END IF
162 CASE (u2dvar)
163 IF (extract_flag.ge.0) THEN
164 imin=bounds(ng)%IstrP(tile)
165 imax=bounds(ng)%IendT(tile)
166 jmin=bounds(ng)%JstrT(tile)
167 jmax=bounds(ng)%JendT(tile)
168#ifdef GRID_EXTRACT
169 ELSE ! Already on extract grid
170 imin=xtr_bounds(ng)%IstrP(tile)
171 imax=xtr_bounds(ng)%IendT(tile)
172 jmin=xtr_bounds(ng)%JstrT(tile)
173 jmax=xtr_bounds(ng)%JendT(tile)
174#endif
175 END IF
176 CASE (v2dvar)
177 IF (extract_flag.ge.0) THEN
178 imin=bounds(ng)%IstrT(tile)
179 imax=bounds(ng)%IendT(tile)
180 jmin=bounds(ng)%JstrP(tile)
181 jmax=bounds(ng)%JendT(tile)
182#ifdef GRID_EXTRACT
183 ELSE ! Already on extract grid
184 imin=xtr_bounds(ng)%IstrT(tile)
185 imax=xtr_bounds(ng)%IendT(tile)
186 jmin=xtr_bounds(ng)%JstrP(tile)
187 jmax=xtr_bounds(ng)%JendT(tile)
188#endif
189 END IF
190 CASE DEFAULT
191 IF (extract_flag.ge.0) THEN
192 imin=bounds(ng)%IstrT(tile)
193 imax=bounds(ng)%IendT(tile)
194 jmin=bounds(ng)%JstrT(tile)
195 jmax=bounds(ng)%JendT(tile)
196#ifdef GRID_EXTRACT
197 ELSE ! Already on extract grid
198 imin=xtr_bounds(ng)%IstrT(tile)
199 imax=xtr_bounds(ng)%IendT(tile)
200 jmin=xtr_bounds(ng)%JstrT(tile)
201 jmax=xtr_bounds(ng)%JendT(tile)
202#endif
203 END IF
204 END SELECT
205!
206!-----------------------------------------------------------------------
207! Compute field statistics.
208!-----------------------------------------------------------------------
209!
210! Initialize local values.
211!
212 my_count=0
213 my_min= 1.0e+37_r8
214 my_max=-1.0e+37_r8
215 my_avg=0.0_r8
216 my_rms=0.0_r8
217!
218 IF (PRESENT(debug)) THEN
219 lprint=debug
220 ELSE
221 lprint=.false.
222 END IF
223!
224! Compute field mean, range, and processed values count for each tile.
225!
226 IF (PRESENT(fmask)) THEN
227 DO j=jmin,jmax
228 DO i=imin,imax
229 IF (fmask(i,j).gt.0.0_r8) THEN
230 my_count=my_count+1
231 my_min=min(my_min, f(i,j))
232 my_max=max(my_max, f(i,j))
233 my_avg=my_avg+f(i,j)
234 my_rms=my_rms+f(i,j)*f(i,j)
235 END IF
236 END DO
237 END DO
238 ELSE
239 DO j=jmin,jmax
240 DO i=imin,imax
241 my_count=my_count+1
242 my_min=min(my_min, f(i,j))
243 my_max=max(my_max, f(i,j))
244 my_avg=my_avg+f(i,j)
245 my_rms=my_rms+f(i,j)*f(i,j)
246 END DO
247 END DO
248 END IF
249!
250! Compute data checksum value.
251!
252 npts=(imax-imin+1)*(jmax-jmin+1)
253 IF (.not.allocated(cwrk)) allocate ( cwrk(npts) )
254 cwrk=pack(f(imin:imax, jmin:jmax), .true.)
255#ifdef DISTRIBUTE
256 CALL get_hash (cwrk, npts, s%checksum, .true.)
257#else
258 CALL get_hash (cwrk, npts, s%checksum)
259#endif
260 IF (allocated(cwrk)) deallocate (cwrk)
261!
262! Compute global field mean, range, and processed values count.
263!
264#ifdef DISTRIBUTE
265 nsub=1 ! distributed-memory
266#else
267 nsub=ntilex(ng)*ntilee(ng) ! tiled application
268#endif
269!$OMP CRITICAL (F_STATS)
270 s%count=s%count+my_count
271 s%min=min(s%min,my_min)
272 s%max=max(s%max,my_max)
273 s%avg=s%avg+my_avg
274 s%rms=s%rms+my_rms
275 IF (lprint) THEN
276 WRITE (stdout,10) &
277 & 'thread = ', my_threadnum(), &
278 & 'tile = ', tile, &
279 & 'tile_count = ', tile_count, &
280 & 'my_count = ', my_count, &
281 & 'S%count = ', s%count, &
282 & 'MINVAL = ', minval(f(imin:imax,jmin:jmax)), &
283 & 'MAXVAL = ', maxval(f(imin:imax,jmin:jmax)), &
284 & 'SUM = ', sum(f(imin:imax,jmin:jmax)), &
285 & 'MEAN = ', sum(f(imin:imax,jmin:jmax))/ &
286 & real(my_count,r8), &
287 & 'my_min = ', my_min, &
288 & 'my_max = ', my_max, &
289 & 'S%min = ', s%min, &
290 & 'S%max = ', s%max, &
291 & 'my_avg = ', my_avg, &
292 & 'S%avg = ', s%avg, &
293 & 'my_rms = ', my_rms, &
294 & 'S%rms = ', s%rms
295 10 FORMAT (10x,5(5x,a,i0),/, &
296 & 6(15x,a,1p,e15.8,0p,5x,a,1p,e15.8,0p,/))
297 FLUSH (stdout)
298 END IF
299 tile_count=tile_count+1
300 IF (tile_count.eq.nsub) THEN
301 tile_count=0
302#ifdef DISTRIBUTE
303 rbuffer(1)=real(s%count, r8)
304 rbuffer(2)=s%min
305 rbuffer(3)=s%max
306 rbuffer(4)=s%avg
307 rbuffer(5)=s%rms
308 op_handle(1)='SUM'
309 op_handle(2)='MIN'
310 op_handle(3)='MAX'
311 op_handle(4)='SUM'
312 op_handle(5)='SUM'
313 CALL mp_reduce (ng, model, 5, rbuffer, op_handle)
314 s%count=int(rbuffer(1))
315 s%min=rbuffer(2)
316 s%max=rbuffer(3)
317 s%avg=rbuffer(4)
318 s%rms=rbuffer(5)
319#endif
320!
321! Finalize computation of the mean and root mean square.
322!
323 IF (s%count.gt.0) THEN
324 fac=1.0_r8/real(s%count, r8)
325 s%avg=s%avg*fac
326 s%rms=sqrt(s%rms*fac)
327 ELSE
328 s%avg=1.0e+37_r8
329 s%rms=1.0e+37_r8
330 END IF
331 END IF
332!$OMP END CRITICAL (F_STATS)
333!$OMP BARRIER
334!
335 RETURN
336 END SUBROUTINE stats_2dfld
337!
338 SUBROUTINE stats_3dfld (ng, tile, model, gtype, S, &
339 & Extract_Flag, &
340 & LBi, UBi, LBj, UBj, LBk, UBk, &
341 & F, Fmask, debug)
342!
343!=======================================================================
344! !
345! This routine computes requested 3D-field statistics. !
346! !
347! On Input: !
348! !
349! ng Nested grid number (integer) !
350! tile Domain partition (integer) !
351! model Calling model identifier (integer) !
352! gtype Grid type (integer) !
353! Extract_Flag Extraction flag interpolation/decimation (integer) !
354! LBi I-dimension Lower bound (integer) !
355! UBi I-dimension Upper bound (integer) !
356! LBj J-dimension Lower bound (integer) !
357! UBj J-dimension Upper bound (integer) !
358! LBk K-dimension Lower bound (integer) !
359! UBk K-dimension Upper bound (integer) !
360! F 3D state field (3D tiled array) !
361! Fmask Land/sea mask (OPTIONAL, 2D tiled array) !
362! debug Switch to debug (OPTIONAL, logical) !
363! !
364! On Ouput: !
365! !
366! S Field statistics (structure) !
367! !
368!=======================================================================
369!
370 USE mod_param
371 USE mod_parallel
372 USE mod_iounits
373 USE mod_ncparam
374
375#ifdef DISTRIBUTE
376!
377 USE distribute_mod, ONLY : mp_reduce
378#endif
379 USE get_hash_mod, ONLY : get_hash
380!
381! Imported variable declarations.
382!
383 logical, intent(in), optional :: debug
384!
385 integer, intent(in) :: ng, tile, model, gtype
386 integer, intent(in) :: extract_flag
387 integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk
388!
389#ifdef ASSUMED_SHAPE
390 real(r8), intent(in) :: f(lbi:,lbj:,lbk:)
391 real(r8), intent(in), optional :: fmask(lbi:,lbj:)
392#else
393 real(r8), intent(in) :: f(lbi:ubi,lbj:ubj,lbk:ubk)
394 real(r8), intent(in), optional :: fmask(lbi:ubi,lbj:ubj)
395#endif
396 TYPE(t_stats), intent(inout) :: s
397!
398! Local variable declarations.
399!
400 logical :: lprint
401!
402 integer :: imin, imax, jmin, jmax, npts, nsub
403 integer :: i, j, k, my_count
404
405 integer :: my_threadnum
406!
407 real(r8) :: fac
408 real(r8) :: my_max, my_min
409 real(r8) :: my_avg, my_rms
410
411 real(r8), allocatable :: cwrk(:)
412
413#ifdef DISTRIBUTE
414 real(r8), dimension(5) :: rbuffer
415 character (len=3), dimension(5) :: op_handle
416#endif
417!
418!-----------------------------------------------------------------------
419! Set tile indices.
420!-----------------------------------------------------------------------
421!
422! Determine I- and J-indices according to staggered C-grid type.
423#ifdef GRID_EXTRACT
424!
425! Here, the ELSE statements (Extract_Flag < 0) indicate the processing
426! of data already on the extracted grid, like geometry variables
427! requiring no extraction since they are read during initialization.
428#endif
429!
430 SELECT CASE (abs(gtype))
431 CASE (p2dvar)
432 IF (extract_flag.ge.0) THEN
433 imin=bounds(ng)%IstrP(tile)
434 imax=bounds(ng)%IendP(tile)
435 jmin=bounds(ng)%JstrP(tile)
436 jmax=bounds(ng)%JendP(tile)
437#ifdef GRID_EXTRACT
438 ELSE ! Already on extract grid
439 imin=xtr_bounds(ng)%IstrP(tile)
440 imax=xtr_bounds(ng)%IendP(tile)
441 jmin=xtr_bounds(ng)%JstrP(tile)
442 jmax=xtr_bounds(ng)%JendP(tile)
443#endif
444 END IF
445 CASE (r2dvar)
446 IF (extract_flag.ge.0) THEN
447 imin=bounds(ng)%IstrT(tile)
448 imax=bounds(ng)%IendT(tile)
449 jmin=bounds(ng)%JstrT(tile)
450 jmax=bounds(ng)%JendT(tile)
451#ifdef GRID_EXTRACT
452 ELSE ! Already on extract grid
453 imin=xtr_bounds(ng)%IstrT(tile)
454 imax=xtr_bounds(ng)%IendT(tile)
455 jmin=xtr_bounds(ng)%JstrT(tile)
456 jmax=xtr_bounds(ng)%JendT(tile)
457#endif
458 END IF
459 CASE (u2dvar)
460 IF (extract_flag.ge.0) THEN
461 imin=bounds(ng)%IstrP(tile)
462 imax=bounds(ng)%IendT(tile)
463 jmin=bounds(ng)%JstrT(tile)
464 jmax=bounds(ng)%JendT(tile)
465#ifdef GRID_EXTRACT
466 ELSE ! Already on extract grid
467 imin=xtr_bounds(ng)%IstrP(tile)
468 imax=xtr_bounds(ng)%IendT(tile)
469 jmin=xtr_bounds(ng)%JstrT(tile)
470 jmax=xtr_bounds(ng)%JendT(tile)
471#endif
472 END IF
473 CASE (v2dvar)
474 IF (extract_flag.ge.0) THEN
475 imin=bounds(ng)%IstrT(tile)
476 imax=bounds(ng)%IendT(tile)
477 jmin=bounds(ng)%JstrP(tile)
478 jmax=bounds(ng)%JendT(tile)
479#ifdef GRID_EXTRACT
480 ELSE ! Already on extract grid
481 imin=xtr_bounds(ng)%IstrT(tile)
482 imax=xtr_bounds(ng)%IendT(tile)
483 jmin=xtr_bounds(ng)%JstrP(tile)
484 jmax=xtr_bounds(ng)%JendT(tile)
485#endif
486 END IF
487 CASE DEFAULT
488 IF (extract_flag.ge.0) THEN
489 imin=bounds(ng)%IstrT(tile)
490 imax=bounds(ng)%IendT(tile)
491 jmin=bounds(ng)%JstrT(tile)
492 jmax=bounds(ng)%JendT(tile)
493#ifdef GRID_EXTRACT
494 ELSE ! Already on extract grid
495 imin=xtr_bounds(ng)%IstrT(tile)
496 imax=xtr_bounds(ng)%IendT(tile)
497 jmin=xtr_bounds(ng)%JstrT(tile)
498 jmax=xtr_bounds(ng)%JendT(tile)
499#endif
500 END IF
501 END SELECT
502!
503!-----------------------------------------------------------------------
504! Compute field statistics.
505!-----------------------------------------------------------------------
506!
507! Initialize local values.
508!
509 my_count=0
510 my_min= 1.0e+37_r8
511 my_max=-1.0e+37_r8
512 my_avg=0.0_r8
513 my_rms=0.0_r8
514!
515 IF (PRESENT(debug)) THEN
516 lprint=debug
517 ELSE
518 lprint=.false.
519 END IF
520!
521! Compute field mean, range, and processed values count for each tile.
522!
523 IF (PRESENT(fmask)) THEN
524 DO k=lbk,ubk
525 DO j=jmin,jmax
526 DO i=imin,imax
527 IF (fmask(i,j).gt.0.0_r8) THEN
528 my_count=my_count+1
529 my_min=min(my_min, f(i,j,k))
530 my_max=max(my_max, f(i,j,k))
531 my_avg=my_avg+f(i,j,k)
532 my_rms=my_rms+f(i,j,k)*f(i,j,k)
533 END IF
534 END DO
535 END DO
536 END DO
537 ELSE
538 DO k=lbk,ubk
539 DO j=jmin,jmax
540 DO i=imin,imax
541 my_count=my_count+1
542 my_min=min(my_min, f(i,j,k))
543 my_max=max(my_max, f(i,j,k))
544 my_avg=my_avg+f(i,j,k)
545 my_rms=my_rms+f(i,j,k)*f(i,j,k)
546 END DO
547 END DO
548 END DO
549 END IF
550!
551! Compute data checksum value.
552!
553 npts=(imax-imin+1)*(jmax-jmin+1)*(ubk-lbk+1)
554 IF (.not.allocated(cwrk)) allocate ( cwrk(npts) )
555 cwrk=pack(f(imin:imax, jmin:jmax, lbk:ubk), .true.)
556#ifdef DISTRIBUTE
557 CALL get_hash (cwrk, npts, s%checksum, .true.)
558#else
559 CALL get_hash (cwrk, npts, s%checksum)
560#endif
561 IF (allocated(cwrk)) deallocate (cwrk)
562!
563! Compute global field mean, range, and processed values count.
564! Notice that the critical region F_STATS is the same as in
565! "stats_2dfld" but we are usinf different counter "thread_count"
566! to avoid interference and race conditions in shared-memory.
567!
568#ifdef DISTRIBUTE
569 nsub=1 ! distributed-memory
570#else
571 nsub=ntilex(ng)*ntilee(ng) ! tiled application
572#endif
573!$OMP CRITICAL (F_STATS)
574 s%count=s%count+my_count
575 s%min=min(s%min,my_min)
576 s%max=max(s%max,my_max)
577 s%avg=s%avg+my_avg
578 s%rms=s%rms+my_rms
579 IF (lprint) THEN
580 WRITE (stdout,10) &
581 & 'thread = ', my_threadnum(), &
582 & 'tile = ', tile, &
583 & 'tile_count = ', thread_count, &
584 & 'my_count = ', my_count, &
585 & 'S%count = ', s%count, &
586 & 'MINVAL = ', minval(f(imin:imax,jmin:jmax,:)), &
587 & 'MAXVAL = ', maxval(f(imin:imax,jmin:jmax,:)), &
588 & 'SUM = ', sum(f(imin:imax,jmin:jmax,:)), &
589 & 'MEAN = ', sum(f(imin:imax,jmin:jmax,:))/ &
590 & real(my_count,r8), &
591 & 'my_min = ', my_min, &
592 & 'my_max = ', my_max, &
593 & 'S%min = ', s%min, &
594 & 'S%max = ', s%max, &
595 & 'my_avg = ', my_avg, &
596 & 'S%avg = ', s%avg, &
597 & 'my_rms = ', my_rms, &
598 & 'S%rms = ', s%rms
599 10 FORMAT (10x,5(5x,a,i0),/, &
600 & 6(15x,a,1p,e15.8,0p,5x,a,1p,e15.8,0p,/))
601 FLUSH (stdout)
602 END IF
603 thread_count=thread_count+1
604 IF (thread_count.eq.nsub) THEN
605 thread_count=0
606#ifdef DISTRIBUTE
607 rbuffer(1)=real(s%count, r8)
608 rbuffer(2)=s%min
609 rbuffer(3)=s%max
610 rbuffer(4)=s%avg
611 rbuffer(5)=s%rms
612 op_handle(1)='SUM'
613 op_handle(2)='MIN'
614 op_handle(3)='MAX'
615 op_handle(4)='SUM'
616 op_handle(5)='SUM'
617 CALL mp_reduce (ng, model, 5, rbuffer, op_handle)
618 s%count=int(rbuffer(1))
619 s%min=rbuffer(2)
620 s%max=rbuffer(3)
621 s%avg=rbuffer(4)
622 s%rms=rbuffer(5)
623#endif
624!
625! Finalize computation of the mean and root mean square.
626!
627 IF (s%count.gt.0) THEN
628 fac=1.0_r8/real(s%count, r8)
629 s%avg=s%avg*fac
630 s%rms=sqrt(s%rms*fac)
631 ELSE
632 s%avg=1.0e+37_r8
633 s%rms=1.0e+37_r8
634 END IF
635 END IF
636!$OMP END CRITICAL (F_STATS)
637!$OMP BARRIER
638!
639 RETURN
640 END SUBROUTINE stats_3dfld
641!
642 SUBROUTINE stats_4dfld (ng, tile, model, gtype, S, &
643 & Extract_Flag, &
644 & LBi, UBi, LBj, UBj, LBk, UBk, LBt, UBt, &
645 & F, Fmask, debug)
646!
647!=======================================================================
648! !
649! This routine computes requested 4D-field statistics. !
650! !
651! On Input: !
652! !
653! ng Nested grid number (integer) !
654! tile Domain partition (integer) !
655! model Calling model identifier (integer) !
656! gtype Grid type (integer) !
657! Extract_Flag Extraction flag interpolation/decimation (integer) !
658! LBi I-dimension Lower bound (integer) !
659! UBi I-dimension Upper bound (integer) !
660! LBj J-dimension Lower bound (integer) !
661! UBj J-dimension Upper bound (integer) !
662! LBk K-dimension Lower bound (integer) !
663! UBk K-dimension Upper bound (integer) !
664! LBt Time-dimension Lower bound (integer) !
665! UBt Time-dimension Upper bound (integer) !
666! F 4D state field (4D tiled array) !
667! Fmask Land/sea mask (OPTIONAL, 2D tiled array) !
668! debug Switch to debug (OPTIONAL, logical) !
669! !
670! On Ouput: !
671! !
672! S Field statistics (structure) !
673! !
674!=======================================================================
675!
676 USE mod_param
677 USE mod_parallel
678 USE mod_iounits
679 USE mod_ncparam
680
681#ifdef DISTRIBUTE
682!
683 USE distribute_mod, ONLY : mp_reduce
684#endif
685 USE get_hash_mod, ONLY : get_hash
686!
687! Imported variable declarations.
688!
689 logical, intent(in), optional :: debug
690!
691 integer, intent(in) :: ng, tile, model, gtype
692 integer, intent(in) :: extract_flag
693 integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt
694!
695#ifdef ASSUMED_SHAPE
696 real(r8), intent(in) :: f(lbi:,lbj:,lbk:,lbt:)
697 real(r8), intent(in), optional :: fmask(lbi:,lbj:)
698#else
699 real(r8), intent(in) :: f(lbi:ubi,lbj:ubj,lbk:ubk)
700 real(r8), intent(in), optional :: fmask(lbi:ubi,lbj:ubj,lbt:ubt)
701#endif
702 TYPE(t_stats), intent(inout) :: s
703!
704! Local variable declarations.
705!
706 logical :: lprint
707!
708 integer :: imin, imax, jmin, jmax, npts, nsub
709 integer :: i, j, k, l, my_count
710
711 integer :: my_threadnum
712!
713 real(r8) :: fac
714 real(r8) :: my_max, my_min
715 real(r8) :: my_avg, my_rms
716
717 real(r8), allocatable :: cwrk(:)
718
719#ifdef DISTRIBUTE
720 real(r8), dimension(5) :: rbuffer
721 character (len=3), dimension(5) :: op_handle
722#endif
723!
724!-----------------------------------------------------------------------
725! Set tile indices.
726!-----------------------------------------------------------------------
727!
728! Determine I- and J-indices according to staggered C-grid type.
729#ifdef GRID_EXTRACT
730!
731! Here, the ELSE statements (Extract_Flag < 0) indicate the processing
732! of data already on the extracted grid, like geometry variables
733! requiring no extraction since they are read during initialization.
734#endif
735!
736 SELECT CASE (abs(gtype))
737 CASE (p2dvar)
738 IF (extract_flag.ge.0) THEN
739 imin=bounds(ng)%IstrP(tile)
740 imax=bounds(ng)%IendP(tile)
741 jmin=bounds(ng)%JstrP(tile)
742 jmax=bounds(ng)%JendP(tile)
743#ifdef GRID_EXTRACT
744 ELSE ! Already on extract grid
745 imin=xtr_bounds(ng)%IstrP(tile)
746 imax=xtr_bounds(ng)%IendP(tile)
747 jmin=xtr_bounds(ng)%JstrP(tile)
748 jmax=xtr_bounds(ng)%JendP(tile)
749#endif
750 END IF
751 CASE (r2dvar)
752 IF (extract_flag.ge.0) THEN
753 imin=bounds(ng)%IstrT(tile)
754 imax=bounds(ng)%IendT(tile)
755 jmin=bounds(ng)%JstrT(tile)
756 jmax=bounds(ng)%JendT(tile)
757#ifdef GRID_EXTRACT
758 ELSE ! Already on extract grid
759 imin=xtr_bounds(ng)%IstrT(tile)
760 imax=xtr_bounds(ng)%IendT(tile)
761 jmin=xtr_bounds(ng)%JstrT(tile)
762 jmax=xtr_bounds(ng)%JendT(tile)
763#endif
764 END IF
765 CASE (u2dvar)
766 IF (extract_flag.ge.0) THEN
767 imin=bounds(ng)%IstrP(tile)
768 imax=bounds(ng)%IendT(tile)
769 jmin=bounds(ng)%JstrT(tile)
770 jmax=bounds(ng)%JendT(tile)
771#ifdef GRID_EXTRACT
772 ELSE ! Already on extract grid
773 imin=xtr_bounds(ng)%IstrP(tile)
774 imax=xtr_bounds(ng)%IendT(tile)
775 jmin=xtr_bounds(ng)%JstrT(tile)
776 jmax=xtr_bounds(ng)%JendT(tile)
777#endif
778 END IF
779 CASE (v2dvar)
780 IF (extract_flag.ge.0) THEN
781 imin=bounds(ng)%IstrT(tile)
782 imax=bounds(ng)%IendT(tile)
783 jmin=bounds(ng)%JstrP(tile)
784 jmax=bounds(ng)%JendT(tile)
785#ifdef GRID_EXTRACT
786 ELSE ! Already on extract grid
787 imin=xtr_bounds(ng)%IstrT(tile)
788 imax=xtr_bounds(ng)%IendT(tile)
789 jmin=xtr_bounds(ng)%JstrP(tile)
790 jmax=xtr_bounds(ng)%JendT(tile)
791#endif
792 END IF
793 CASE DEFAULT
794 IF (extract_flag.ge.0) THEN
795 imin=bounds(ng)%IstrT(tile)
796 imax=bounds(ng)%IendT(tile)
797 jmin=bounds(ng)%JstrT(tile)
798 jmax=bounds(ng)%JendT(tile)
799#ifdef GRID_EXTRACT
800 ELSE ! Already on extract grid
801 imin=xtr_bounds(ng)%IstrT(tile)
802 imax=xtr_bounds(ng)%IendT(tile)
803 jmin=xtr_bounds(ng)%JstrT(tile)
804 jmax=xtr_bounds(ng)%JendT(tile)
805#endif
806 END IF
807 END SELECT
808!
809!-----------------------------------------------------------------------
810! Compute field statistics.
811!-----------------------------------------------------------------------
812!
813! Initialize local values.
814!
815 my_count=0
816 my_min= 1.0e+37_r8
817 my_max=-1.0e+37_r8
818 my_avg=0.0_r8
819 my_rms=0.0_r8
820!
821 IF (PRESENT(debug)) THEN
822 lprint=debug
823 ELSE
824 lprint=.false.
825 END IF
826!
827! Compute field mean, range, and processed values count for each tile.
828!
829 IF (PRESENT(fmask)) THEN
830 DO l=lbt,ubt
831 DO k=lbk,ubk
832 DO j=jmin,jmax
833 DO i=imin,imax
834 IF (fmask(i,j).gt.0.0_r8) THEN
835 my_count=my_count+1
836 my_min=min(my_min, f(i,j,k,l))
837 my_max=max(my_max, f(i,j,k,l))
838 my_avg=my_avg+f(i,j,k,l)
839 my_rms=my_rms+f(i,j,k,l)*f(i,j,k,l)
840 END IF
841 END DO
842 END DO
843 END DO
844 END DO
845 ELSE
846 DO l=lbt,ubt
847 DO k=lbk,ubk
848 DO j=jmin,jmax
849 DO i=imin,imax
850 my_count=my_count+1
851 my_min=min(my_min, f(i,j,k,l))
852 my_max=max(my_max, f(i,j,k,l))
853 my_avg=my_avg+f(i,j,k,l)
854 my_rms=my_rms+f(i,j,k,l)*f(i,j,k,l)
855 END DO
856 END DO
857 END DO
858 END DO
859 END IF
860!
861! Compute data checksum value.
862!
863 npts=(imax-imin+1)*(jmax-jmin+1)*(ubk-lbk+1)*(ubt-lbt+1)
864 IF (.not.allocated(cwrk)) allocate ( cwrk(npts) )
865 cwrk=pack(f(imin:imax, jmin:jmax, lbk:ubk, lbt:ubt), .true.)
866#ifdef DISTRIBUTE
867 CALL get_hash (cwrk, npts, s%checksum, .true.)
868#else
869 CALL get_hash (cwrk, npts, s%checksum)
870#endif
871 IF (allocated(cwrk)) deallocate (cwrk)
872!
873! Compute global field mean, range, and processed values count.
874! Notice that the critical region F_STATS is the same as in
875! "stats_2dfld" but we are usinf different counter "thread_count"
876! to avoid interference and race conditions in shared-memory.
877!
878#ifdef DISTRIBUTE
879 nsub=1 ! distributed-memory
880#else
881 nsub=ntilex(ng)*ntilee(ng) ! tiled application
882#endif
883!$OMP CRITICAL (F_STATS)
884 s%count=s%count+my_count
885 s%min=min(s%min,my_min)
886 s%max=max(s%max,my_max)
887 s%avg=s%avg+my_avg
888 s%rms=s%rms+my_rms
889 IF (lprint) THEN
890 WRITE (stdout,10) &
891 & 'thread = ', my_threadnum(), &
892 & 'tile = ', tile, &
893 & 'tile_count = ', thread_count, &
894 & 'my_count = ', my_count, &
895 & 'S%count = ', s%count, &
896 & 'MINVAL = ', minval(f(imin:imax,jmin:jmax,:,:)), &
897 & 'MAXVAL = ', maxval(f(imin:imax,jmin:jmax,:,:)), &
898 & 'SUM = ', sum(f(imin:imax,jmin:jmax,:,:)), &
899 & 'MEAN = ', sum(f(imin:imax,jmin:jmax,:,:))/ &
900 & real(my_count,r8), &
901 & 'my_min = ', my_min, &
902 & 'my_max = ', my_max, &
903 & 'S%min = ', s%min, &
904 & 'S%max = ', s%max, &
905 & 'my_avg = ', my_avg, &
906 & 'S%avg = ', s%avg, &
907 & 'my_rms = ', my_rms, &
908 & 'S%rms = ', s%rms
909 10 FORMAT (10x,5(5x,a,i0),/, &
910 & 6(15x,a,1p,e15.8,0p,5x,a,1p,e15.8,0p,/))
911 FLUSH (stdout)
912 END IF
913 thread_count=thread_count+1
914 IF (thread_count.eq.nsub) THEN
915 thread_count=0
916#ifdef DISTRIBUTE
917 rbuffer(1)=real(s%count, r8)
918 rbuffer(2)=s%min
919 rbuffer(3)=s%max
920 rbuffer(4)=s%avg
921 rbuffer(5)=s%rms
922 op_handle(1)='SUM'
923 op_handle(2)='MIN'
924 op_handle(3)='MAX'
925 op_handle(4)='SUM'
926 op_handle(5)='SUM'
927 CALL mp_reduce (ng, model, 5, rbuffer, op_handle)
928 s%count=int(rbuffer(1))
929 s%min=rbuffer(2)
930 s%max=rbuffer(3)
931 s%avg=rbuffer(4)
932 s%rms=rbuffer(5)
933#endif
934!
935! Finalize computation of the mean and root mean square.
936!
937 IF (s%count.gt.0) THEN
938 fac=1.0_r8/real(s%count, r8)
939 s%avg=s%avg*fac
940 s%rms=sqrt(s%rms*fac)
941 ELSE
942 s%avg=1.0e+37_r8
943 s%rms=1.0e+37_r8
944 END IF
945 END IF
946!$OMP END CRITICAL (F_STATS)
947!$OMP BARRIER
948!
949 RETURN
950 END SUBROUTINE stats_4dfld
951!
952 END MODULE stats_mod
distribute_mod::mp_reduce
Definition distribute.F:113
my_threadnum
integer function my_threadnum()
Definition mp_routines.F:121
distribute_mod
Definition distribute.F:3
get_hash_mod
Definition get_hash.F:5
get_hash_mod::get_hash
subroutine, public get_hash(a, asize, hash, lreduce)
Definition get_hash.F:72
mod_iounits
Definition mod_iounits.F:2
mod_iounits::stdout
integer stdout
Definition mod_iounits.F:195
mod_ncparam
Definition mod_ncparam.F:2
mod_parallel
Definition mod_parallel.F:2
mod_parallel::thread_count
integer thread_count
Definition mod_parallel.F:91
mod_parallel::tile_count
integer tile_count
Definition mod_parallel.F:89
mod_param
Definition mod_param.F:2
mod_param::bounds
type(t_bounds), dimension(:), allocatable bounds
Definition mod_param.F:232
mod_param::ntilex
integer, dimension(:), allocatable ntilex
Definition mod_param.F:685
mod_param::u2dvar
integer, parameter u2dvar
Definition mod_param.F:718
mod_param::ntilee
integer, dimension(:), allocatable ntilee
Definition mod_param.F:686
mod_param::p2dvar
integer, parameter p2dvar
Definition mod_param.F:716
mod_param::r2dvar
integer, parameter r2dvar
Definition mod_param.F:717
mod_param::v2dvar
integer, parameter v2dvar
Definition mod_param.F:719
stats_mod
Definition stats.F:2
stats_mod::stats_2dfld
subroutine, public stats_2dfld(ng, tile, model, gtype, s, extract_flag, lbi, ubi, lbj, ubj, f, fmask, debug)
Definition stats.F:47
stats_mod::stats_4dfld
subroutine, public stats_4dfld(ng, tile, model, gtype, s, extract_flag, lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt, f, fmask, debug)
Definition stats.F:646
stats_mod::stats_3dfld
subroutine, public stats_3dfld(ng, tile, model, gtype, s, extract_flag, lbi, ubi, lbj, ubj, lbk, ubk, f, fmask, debug)
Definition stats.F:342
mod_param::t_stats
Definition mod_param.F:420