SCALE-RM
scale_statistics.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 #include "scalelib.h"
12  !-----------------------------------------------------------------------------
13  !
14  !++ used modules
15  !
16  use mpi ! TODO: to replace functions in scale_comm module
17  use scale_precision
18  use scale_io
19  use scale_prof
20  use scale_const, only: &
21  eps => const_eps
22  use scale_prc, only: &
24 #ifdef _OPENACC
25  use openacc
26 #endif
27  !-----------------------------------------------------------------------------
28  implicit none
29  private
30  !-----------------------------------------------------------------------------
31  !
32  !++ Public procedure
33  !
34  public :: statistics_setup
35  public :: statistics_total
36  public :: statistics_detail
37  public :: statistics_horizontal_mean
38  public :: statistics_horizontal_min
39  public :: statistics_horizontal_max
40 
41  public :: statistics_summation
42  public :: statistics_average
43  public :: statistics_variance
44  public :: statistics_stddev
45 
46  public :: statistics_covariance
47  public :: statistics_correlation
48  public :: statistics_regression
51 
52  public :: statistics_undef_replace
53  public :: statistics_undef_embed
55 
56  interface statistics_total
57  module procedure statistics_total_2d
58  module procedure statistics_total_3d
59  end interface statistics_total
60 
61  interface statistics_detail
62  module procedure statistics_detail_2d
63  module procedure statistics_detail_3d
64  end interface statistics_detail
65 
66  interface statistics_horizontal_mean
67  module procedure statistics_horizontal_mean_2d
68  module procedure statistics_horizontal_mean_3d
69  end interface statistics_horizontal_mean
70 
71  interface statistics_horizontal_max
72  module procedure statistics_horizontal_max_2d
73  module procedure statistics_horizontal_max_3d
74  end interface statistics_horizontal_max
75 
76  interface statistics_horizontal_min
77  module procedure statistics_horizontal_min_2d
78  module procedure statistics_horizontal_min_3d
79  end interface statistics_horizontal_min
80 
81  interface statistics_summation
82  module procedure statistics_summation_1d
83  module procedure statistics_summation_2d
84  module procedure statistics_summation_3d
85  end interface statistics_summation
86 
87  interface statistics_average
88  module procedure statistics_average_1d
89  module procedure statistics_average_2d
90  module procedure statistics_average_3d
91  end interface statistics_average
92 
93  interface statistics_variance
94  module procedure statistics_variance_1d
95  module procedure statistics_variance_2d
96  module procedure statistics_variance_3d
97  end interface statistics_variance
98 
99  interface statistics_stddev
100  module procedure statistics_stddev_1d
101  module procedure statistics_stddev_2d
102  module procedure statistics_stddev_3d
103  end interface statistics_stddev
104 
105  !-----------------------------------------------------------------------------
106  !
107  !++ Public parameters & variables
108  !
109  logical, public :: statistics_checktotal = .false.
110 
111  !-----------------------------------------------------------------------------
112  !
113  !++ Private procedure
114  !
115  !-----------------------------------------------------------------------------
116  !
117  !++ Private parameters & variables
118  !
119  logical, private :: statistics_use_globalcomm = .false.
120 
121  !-----------------------------------------------------------------------------
122 contains
123  !-----------------------------------------------------------------------------
125  subroutine statistics_setup
126  use scale_prc, only: &
127  prc_abort
128  use scale_comm_cartesc, only: &
129  comm_setup
130  implicit none
131 
132  namelist / param_statistics / &
134  statistics_use_globalcomm
135 
136  integer :: ierr
137  !---------------------------------------------------------------------------
138 
139  call comm_setup
140 
141  log_newline
142  log_info("STATISTICS_setup",*) 'Setup'
143 
144  !--- read namelist
145  rewind(io_fid_conf)
146  read(io_fid_conf,nml=param_statistics,iostat=ierr)
147  if( ierr < 0 ) then !--- missing
148  log_info("STATISTICS_setup",*) 'Not found namelist. Default used.'
149  elseif( ierr > 0 ) then !--- fatal error
150  log_error("STATISTICS_setup",*) 'Not appropriate names in namelist PARAM_STATISTICS. Check!'
151  call prc_abort
152  endif
153  log_nml(param_statistics)
154 
155  log_newline
156  log_info("STATISTICS_setup",*) 'Caluculate total statistics for monitoring? : ', statistics_checktotal
157  if ( statistics_use_globalcomm ) then
158  log_info_cont(*) '=> The total is calculated for the global domain.'
159  else
160  log_info_cont(*) '=> The total is calculated for the local domain.'
161  endif
162 
163  return
164  end subroutine statistics_setup
165 
166  !-----------------------------------------------------------------------------
168  subroutine statistics_total_2d( &
169  IA, IS, IE, JA, JS, JE, &
170  var, varname, &
171  area, total, &
172  log_suppress, &
173  global, &
174  mean, sum )
175  use scale_prc, only: &
176  prc_myrank, &
177  prc_abort
178  use scale_const, only: &
179  eps => const_eps, &
180  undef => const_undef
181  implicit none
182 
183  integer, intent(in) :: IA, IS, IE
184  integer, intent(in) :: JA, JS, JE
185 
186  real(RP), intent(in) :: var(IA,JA)
187  character(len=*), intent(in) :: varname
188  real(RP), intent(in) :: area(IA,JA)
189  real(RP), intent(in) :: total
190 
191  logical, intent(in), optional :: log_suppress
192  logical, intent(in), optional :: global
193  real(RP), intent(out), optional :: mean
194  real(DP), intent(out), optional :: sum
195 
196  real(DP) :: statval
197  real(DP) :: sendbuf(2), recvbuf(2)
198  real(DP) :: sum_, mean_
199 
200  logical :: suppress_, global_
201  integer :: ierr
202  integer :: i, j
203  !---------------------------------------------------------------------------
204 
205  statval = 0.0_dp
206  !$acc update host(var(IS,JS)) if(acc_is_present(var))
207  if ( var(is,js) /= undef ) then
208  !$omp parallel do private(i,j) OMP_SCHEDULE_ collapse(2) reduction(+:statval)
209  !$acc kernels copyin(var, area) if(acc_is_present(var))
210  !$acc loop reduction(+:statval)
211  do j = js, je
212  !$acc loop reduction(+:statval)
213  do i = is, ie
214  statval = statval + var(i,j) * area(i,j)
215  end do
216  end do
217  !$acc end kernels
218  end if
219 
220  if ( .NOT. ( statval > -1.0_dp .OR. statval < 1.0_dp ) ) then ! must be NaN
221  log_error("STATISTICS_total_2D",*) 'NaN is detected for ', trim(varname), ' in rank ', prc_myrank
222  call prc_abort
223  endif
224 
225  if ( present(log_suppress) ) then
226  suppress_ = log_suppress
227  else
228  suppress_ = .false.
229  end if
230 
231  if ( present(global) ) then
232  global_ = global
233  else
234  global_ = statistics_use_globalcomm
235  end if
236 
237  if ( global_ ) then
238  call prof_rapstart('COMM_Allreduce', 2)
239  sendbuf(1) = statval
240  sendbuf(2) = total
241  ! All reduce
242  call mpi_allreduce( sendbuf(:), recvbuf(:), &
243  2, &
244  mpi_double_precision, &
245  mpi_sum, &
246  prc_local_comm_world, &
247  ierr )
248  call prof_rapend ('COMM_Allreduce', 2)
249 
250  if ( recvbuf(2) < eps ) then
251  sum_ = undef
252  mean_ = undef
253  else
254  sum_ = recvbuf(1)
255  mean_ = recvbuf(1) / recvbuf(2)
256  end if
257  ! statistics over the all node
258  if ( .not. suppress_ ) then ! if varname is empty, suppress output
259  log_info("STATISTICS_total_2D",'(1x,A,A24,A,ES24.17)') &
260  '[', trim(varname), '] MEAN(global) = ', mean_
261  endif
262  else
263  if ( total < eps ) then
264  sum_ = undef
265  mean_ = undef
266  else
267  sum_ = statval
268  mean_ = statval / total
269  end if
270 
271  ! statistics on each node
272  if ( .not. suppress_ ) then ! if varname is empty, suppress output
273  log_info("STATISTICS_total_2D",'(1x,A,A24,A,ES24.17)') &
274  '[', trim(varname), '] MEAN(local) = ', mean_
275  endif
276  endif
277 
278  if ( present(mean) ) mean = mean_
279  if ( present(sum ) ) sum = sum_
280 
281  return
282  end subroutine statistics_total_2d
283 
284  !-----------------------------------------------------------------------------
286  subroutine statistics_total_3d( &
287  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
288  var, varname, &
289  vol, total, &
290  log_suppress, &
291  global, &
292  mean, sum )
293  use scale_prc, only: &
294  prc_myrank, &
295  prc_abort
296  use scale_const, only: &
297  eps => const_eps , &
298  undef => const_undef
299  implicit none
300 
301  integer, intent(in) :: KA, KS, KE
302  integer, intent(in) :: IA, IS, IE
303  integer, intent(in) :: JA, JS, JE
304 
305  real(RP), intent(in) :: var(KA,IA,JA)
306  character(len=*), intent(in) :: varname
307  real(RP), intent(in) :: vol(KA,IA,JA)
308  real(RP), intent(in) :: total
309 
310  logical, intent(in), optional :: log_suppress
311  logical, intent(in), optional :: global
312  real(RP), intent(out), optional :: mean
313  real(DP), intent(out), optional :: sum
314 
315  real(DP) :: statval
316  real(DP) :: sendbuf(2), recvbuf(2)
317  real(DP) :: mean_, sum_
318 
319  logical :: suppress_, global_
320 
321  real(DP) :: work
322  integer :: ierr
323  integer :: k, i, j
324  !---------------------------------------------------------------------------
325 
326  statval = 0.0_dp
327  !$acc update host(var(KS,IS,JS)) if(acc_is_present(var))
328  if ( var(ks,is,js) /= undef ) then
329  !$omp parallel do OMP_SCHEDULE_ reduction(+:statval) &
330  !$omp private(work)
331  !$acc kernels copyin(var, vol) if(acc_is_present(var))
332  !$acc loop reduction(statval)
333  do j = je, js, -1
334  !$acc loop reduction(statval)
335  do i = ie, is, -1
336 #ifdef _OPENACC
337  !$acc loop reduction(statval)
338  do k = ke, ks, -1
339  statval = statval + var(k,i,j) * vol(k,i,j)
340  enddo
341 #else
342  work = 0.0_rp
343  do k = ke, ks, -1
344  work = work + var(k,i,j) * vol(k,i,j)
345  enddo
346  statval = statval + work
347 #endif
348  enddo
349  enddo
350  !$acc end kernels
351  end if
352 
353  if ( .NOT. ( statval > -1.0_dp .OR. statval < 1.0_dp ) ) then ! must be NaN
354  log_error("STATISTICS_total_3D",*) 'NaN is detected for ', trim(varname), ' in rank ', prc_myrank
355  call prc_abort
356  endif
357 
358  if ( present(log_suppress) ) then
359  suppress_ = log_suppress
360  else
361  suppress_ = .false.
362  end if
363 
364  if ( present(global) ) then
365  global_ = global
366  else
367  global_ = statistics_use_globalcomm
368  end if
369 
370  if ( global_ ) then
371  call prof_rapstart('COMM_Allreduce', 2)
372  sendbuf(1) = statval
373  sendbuf(2) = total
374  ! All reduce
375  call mpi_allreduce( sendbuf(:), recvbuf(:), &
376  2, &
377  mpi_double_precision, &
378  mpi_sum, &
379  prc_local_comm_world, &
380  ierr )
381  call prof_rapend ('COMM_Allreduce', 2)
382 
383  if ( recvbuf(2) < eps ) then
384  sum_ = undef
385  mean_ = undef
386  else
387  sum_ = recvbuf(1)
388  mean_ = recvbuf(1) / recvbuf(2)
389  end if
390  ! statistics over the all node
391  if ( .not. suppress_ ) then ! if varname is empty, suppress output
392  log_info("STATISTICS_total_3D",'(1x,A,A24,A,ES24.17)') &
393  '[', trim(varname), '] MEAN(global) = ', mean_
394  endif
395  else
396  if ( total < eps ) then
397  sum_ = undef
398  mean_ = undef
399  else
400  sum_ = statval
401  mean_ = statval / total
402  end if
403 
404  ! statistics on each node
405  if ( .not. suppress_ ) then ! if varname is empty, suppress output
406  log_info("STATISTICS_total_3D",'(1x,A,A24,A,ES24.17)') &
407  '[', trim(varname), '] MEAN(local) = ', mean_
408  endif
409  endif
410 
411  if ( present(mean) ) mean = mean_
412  if ( present(sum ) ) sum = sum_
413 
414  return
415  end subroutine statistics_total_3d
416 
417  !-----------------------------------------------------------------------------
419  subroutine statistics_horizontal_mean_2d( &
420  IA, IS, IE, JA, JS, JE, &
421  var, area, &
422  varmean )
423  use scale_const, only: &
424  undef => const_undef
425  integer, intent(in) :: IA, IS, IE
426  integer, intent(in) :: JA, JS, JE
427  real(RP), intent(in) :: var (IA,JA)
428  real(RP), intent(in) :: area(IA,JA)
429  real(RP), intent(out) :: varmean
430 
431  real(DP) :: statval (2)
432  real(DP) :: allstatval(2)
433  real(DP) :: s1, s2
434 
435  integer :: ierr
436  integer :: i, j
437  !---------------------------------------------------------------------------
438 
439  s1 = 0.0_dp
440  s2 = 0.0_dp
441  !$omp parallel do reduction(+:s1,s2)
442  !$acc kernels copyin(area, var) if(acc_is_present(var))
443  !$acc loop reduction(+:s1,s2)
444  do j = js, je
445  !$acc loop reduction(+:s1,s2)
446  do i = is, ie
447  if ( var(i,j) /= undef ) then
448  s1 = s1 + area(i,j) * var(i,j)
449  s2 = s2 + area(i,j)
450  endif
451  enddo
452  enddo
453  !$acc end kernels
454  statval(1) = s1
455  statval(2) = s2
456 
457  call prof_rapstart('COMM_Allreduce', 2)
458  ! All reduce
459  call mpi_allreduce( statval(:), &
460  allstatval(:), &
461  2, &
462  mpi_double_precision, &
463  mpi_sum, &
465  ierr )
466  call prof_rapend ('COMM_Allreduce', 2)
467 
468  if ( allstatval(2) > 0.0_dp ) then
469  varmean = allstatval(1) / allstatval(2)
470  else
471  varmean = undef
472  end if
473 
474  return
475  end subroutine statistics_horizontal_mean_2d
476 
477  subroutine statistics_horizontal_mean_3d( &
478  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
479  var, area, &
480  varmean )
481  use scale_const, only: &
482  undef => const_undef
483  integer, intent(in) :: KA, KS, KE
484  integer, intent(in) :: IA, IS, IE
485  integer, intent(in) :: JA, JS, JE
486 
487  real(RP), intent(in) :: var (KA,IA,JA)
488  real(RP), intent(in) :: area( IA,JA)
489  real(RP), intent(out) :: varmean(KA)
490 
491  real(DP) :: statval (KS:KE,2)
492  real(DP) :: allstatval(KS:KE,2)
493 
494  integer :: ierr
495  integer :: k, i, j
496  !---------------------------------------------------------------------------
497 
498  !$acc data copyin(var, area) copyout(varmean) create(statval, allstatval) if(acc_is_present(var))
499 
500  !$acc kernels if(acc_is_present(var))
501  statval(:,:) = 0.0_dp
502  !$acc end kernels
503 
504 ! !$omp parallel do reduction(+:statval)
505  !$acc kernels if(acc_is_present(var))
506  !$acc loop independent
507  do j = js, je
508  !$acc loop independent
509  do i = is, ie
510  do k = ks, ke
511  if ( var(k,i,j) /= undef ) then
512  !$acc atomic update
513  statval(k,1) = statval(k,1) + area(i,j) * var(k,i,j)
514  !$acc end atomic
515  !$acc atomic update
516  statval(k,2) = statval(k,2) + area(i,j)
517  !$acc end atomic
518  endif
519  enddo
520  enddo
521  enddo
522  !$acc end kernels
523 
524  call prof_rapstart('COMM_Allreduce', 2)
525  ! All reduce
526  !$acc host_data use_device(statval, allstatval) if(acc_is_present(var))
527  call mpi_allreduce( statval(:,:), &
528  allstatval(:,:), &
529  (ke-ks+1)*2, &
530  mpi_double_precision, &
531  mpi_sum, &
533  ierr )
534  !$acc end host_data
535  call prof_rapend ('COMM_Allreduce', 2)
536 
537  !$acc kernels if(acc_is_present(var))
538  do k = ks, ke
539  if ( allstatval(k,2) > 0.0_dp ) then
540  varmean(k) = allstatval(k,1) / allstatval(k,2)
541  else
542  varmean(k) = undef
543  end if
544  enddo
545  !$acc end kernels
546  !$acc kernels if(acc_is_present(var))
547  do k = 1, ks-1
548  varmean(k) = undef
549  end do
550  !$acc end kernels
551  !$acc kernels if(acc_is_present(var))
552  do k = ke+1, ka
553  varmean(k) = undef
554  end do
555  !$acc end kernels
556 
557  !$acc end data
558 
559  return
560  end subroutine statistics_horizontal_mean_3d
561 
562  !-----------------------------------------------------------------------------
564  subroutine statistics_horizontal_min_2d( &
565  IA, IS, IE, JA, JS, JE, &
566  var, &
567  varmin )
568  use scale_const, only: &
569  undef => const_undef, &
570  huge => const_huge
571  use scale_comm_cartesc, only: &
573  integer, intent(in) :: IA, IS, IE
574  integer, intent(in) :: JA, JS, JE
575 
576  real(RP), intent(in) :: var(IA,JA)
577  real(RP), intent(out) :: varmin
578 
579  real(RP) :: statval
580  real(RP) :: allstatval
581 
582  integer :: ierr
583  integer :: i, j
584  !---------------------------------------------------------------------------
585 
586  statval = huge
587  !$omp parallel do reduction(min:statval)
588  !$acc kernels copyin(var) if(acc_is_present(var))
589  !$acc loop reduction(min:statval)
590  do j = js, je
591  !$acc loop reduction(min:statval)
592  do i = is, ie
593  if ( var(i,j) /= undef .and. var(i,j) < statval ) then
594  statval = var(i,j)
595  endif
596  enddo
597  enddo
598  !$acc end kernels
599 
600  call prof_rapstart('COMM_Allreduce', 2)
601  ! All reduce
602  call mpi_allreduce( statval, &
603  allstatval, &
604  1, &
605  comm_datatype, &
606  mpi_min, &
608  ierr )
609  call prof_rapend ('COMM_Allreduce', 2)
610 
611  if ( allstatval < huge ) then
612  varmin = allstatval
613  else
614  varmin = undef
615  end if
616 
617  return
618  end subroutine statistics_horizontal_min_2d
619 
620  subroutine statistics_horizontal_min_3d( &
621  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
622  var, &
623  varmin )
624  use scale_const, only: &
625  undef => const_undef, &
626  huge => const_huge
627  use scale_comm_cartesc, only: &
629  integer, intent(in) :: KA, KS, KE
630  integer, intent(in) :: IA, IS, IE
631  integer, intent(in) :: JA, JS, JE
632 
633  real(RP), intent(in) :: var(KA,IA,JA)
634  real(RP), intent(out) :: varmin(KA)
635 
636  real(RP) :: statval (KA)
637  real(RP) :: allstatval(KA)
638 
639  integer :: ierr
640  integer :: k, i, j
641  !---------------------------------------------------------------------------
642 
643  !$acc data copyin(var) copyout(varmin) create(statval, allstatval) if(acc_is_present(var))
644 
645  !$acc kernels if(acc_is_present(var))
646  statval(:) = huge
647  !$acc end kernels
648 
649 ! !$omp parallel do reduction(min:statval)
650  !$acc kernels if(acc_is_present(var))
651  !$acc loop independent
652  do j = js, je
653  !$acc loop independent
654  do i = is, ie
655  do k = ks, ke
656 #ifdef _OPENACC
657  if ( var(k,i,j) /= undef )then
658  !$acc atomic update
659  statval(k) = min( var(k,i,j), statval(k) )
660  !$acc end atomic
661  endif
662 #else
663  if ( var(k,i,j) /= undef .and. var(k,i,j) < statval(k) ) then
664  statval(k) = var(k,i,j)
665  endif
666 #endif
667  enddo
668  enddo
669  enddo
670  !$acc end kernels
671 
672  call prof_rapstart('COMM_Allreduce', 2)
673  ! All reduce
674  !$acc host_data use_device(statval, allstatval) if(acc_is_present(var))
675  call mpi_allreduce( statval(ks:ke), &
676  allstatval(ks:ke), &
677  ke-ks+1, &
678  comm_datatype, &
679  mpi_min, &
681  ierr )
682  !$acc end host_data
683  call prof_rapend ('COMM_Allreduce', 2)
684 
685  !$acc kernels if(acc_is_present(var))
686  do k = ks, ke
687  if ( allstatval(k) < huge ) then
688  varmin(k) = allstatval(k)
689  else
690  varmin(k) = undef
691  end if
692  enddo
693  !$acc end kernels
694  !$acc kernels if(acc_is_present(var))
695  do k = 1, ks-1
696  varmin(k) = undef
697  end do
698  !$acc end kernels
699  !$acc kernels if(acc_is_present(var))
700  do k = ke+1, ka
701  varmin(k) = undef
702  end do
703  !$acc end kernels
704 
705  !$acc end data
706 
707  return
708  end subroutine statistics_horizontal_min_3d
709 
710  !-----------------------------------------------------------------------------
712  subroutine statistics_horizontal_max_2d( &
713  IA, IS, IE, JA, JS, JE, &
714  var, &
715  varmax )
716  use scale_const, only: &
717  undef => const_undef, &
718  huge => const_huge
719  use scale_comm_cartesc, only: &
721  integer, intent(in) :: IA, IS, IE
722  integer, intent(in) :: JA, JS, JE
723 
724  real(RP), intent(in) :: var(IA,JA)
725  real(RP), intent(out) :: varmax
726 
727  real(RP) :: statval
728  real(RP) :: allstatval
729 
730  integer :: ierr
731  integer :: i, j
732  !---------------------------------------------------------------------------
733 
734  statval = - huge
735  !$omp parallel do reduction(max:statval)
736  !$acc kernels copyin(var) if(acc_is_present(var))
737  !$acc loop reduction(max:statval)
738  do j = js, je
739  !$acc loop reduction(max:statval)
740  do i = is, ie
741  if ( var(i,j) /= undef .and. var(i,j) > statval ) then
742  statval = var(i,j)
743  endif
744  enddo
745  enddo
746  !$acc end kernels
747 
748  call prof_rapstart('COMM_Allreduce', 2)
749  ! All reduce
750  call mpi_allreduce( statval, &
751  allstatval, &
752  1, &
753  comm_datatype, &
754  mpi_max, &
756  ierr )
757  call prof_rapend ('COMM_Allreduce', 2)
758 
759  if ( allstatval > - huge ) then
760  varmax = allstatval
761  else
762  varmax = undef
763  end if
764 
765  return
766  end subroutine statistics_horizontal_max_2d
767 
768  subroutine statistics_horizontal_max_3d( &
769  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
770  var, &
771  varmax )
772  use scale_const, only: &
773  undef => const_undef, &
774  huge => const_huge
775  use scale_comm_cartesc, only: &
777  integer, intent(in) :: KA, KS, KE
778  integer, intent(in) :: IA, IS, IE
779  integer, intent(in) :: JA, JS, JE
780 
781  real(RP), intent(in) :: var(KA,IA,JA)
782  real(RP), intent(out) :: varmax(KA)
783 
784  real(RP) :: statval (KA)
785  real(RP) :: allstatval(KA)
786 
787  integer :: ierr
788  integer :: k, i, j
789  !---------------------------------------------------------------------------
790 
791  !$acc data copyin(var) copyout(varmax) create(statval, allstatval) if(acc_is_present(var))
792 
793  !$acc kernels if(acc_is_present(var))
794  statval(:) = - huge
795  !$acc end kernels
796 
797 ! !$omp parallel do reduction(max:statval)
798  !$acc kernels if(acc_is_present(var))
799  !$acc loop independent
800  do j = js, je
801  !$acc loop independent
802  do i = is, ie
803  do k = ks, ke
804 #ifdef _OPENACC
805  if ( var(k,i,j) /= undef ) then
806  !$acc atomic update
807  statval(k) = max( var(k,i,j), statval(k) )
808  !$acc end atomic
809  endif
810 #else
811  if ( var(k,i,j) /= undef .and. var(k,i,j) > statval(k) ) then
812  statval(k) = var(k,i,j)
813  endif
814 #endif
815  enddo
816  enddo
817  enddo
818  !$acc end kernels
819 
820  call prof_rapstart('COMM_Allreduce', 2)
821  ! All reduce
822  !$acc host_data use_device(statval, allstatval) if(acc_is_present(var))
823  call mpi_allreduce( statval(ks:ke), &
824  allstatval(ks:ke), &
825  ke-ks+1, &
826  comm_datatype, &
827  mpi_max, &
829  ierr )
830  !$acc end host_data
831  call prof_rapend ('COMM_Allreduce', 2)
832 
833  !$acc kernels if(acc_is_present(var))
834  do k = ks, ke
835  if ( allstatval(k) > - huge ) then
836  varmax(k) = allstatval(k)
837  else
838  varmax(k) = undef
839  end if
840  enddo
841  !$acc end kernels
842  !$acc kernels if(acc_is_present(var))
843  do k = 1, ks-1
844  varmax(k) = undef
845  end do
846  !$acc end kernels
847  !$acc kernels if(acc_is_present(var))
848  do k = ke+1, ka
849  varmax(k) = undef
850  end do
851  !$acc end kernels
852 
853  !$acc end data
854 
855  return
856  end subroutine statistics_horizontal_max_3d
857 
858  !-----------------------------------------------------------------------------
860  subroutine statistics_detail_3d( &
861  KA, KS, KE, IA, IS, IE, JA, JS, JE, VA, &
862  varname, var, &
863  local )
864  use scale_prc, only: &
865  prc_nprocs
866  use scale_comm_cartesc, only: &
868  implicit none
869 
870  integer, intent(in) :: KA, KS, KE
871  integer, intent(in) :: IA, IS, IE
872  integer, intent(in) :: JA, JS, JE
873  integer, intent(in) :: VA
874 
875  character(len=*), intent(in) :: varname(VA)
876  real(RP), intent(in) :: var(KA,IA,JA,VA)
877 
878  logical, intent(in), optional :: local
879 
880  real(RP) :: statval_l ( VA,2)
881  integer :: statidx_l (3,VA,2)
882  real(RP) :: statval ( VA,2,0:PRC_nprocs-1)
883  integer :: statidx (3,VA,2,0:PRC_nprocs-1)
884  real(RP) :: allstatval(VA,2)
885  integer :: allstatidx(VA,2)
886  logical :: do_globalcomm
887 
888  integer :: k, i, j
889  integer :: ierr
890  integer :: v, p
891  !---------------------------------------------------------------------------
892 
893  do_globalcomm = statistics_use_globalcomm
894  if ( present(local) ) do_globalcomm = ( .not. local )
895 
896  !$acc update host(var) if ( acc_is_present(var) )
897 
898  log_newline
899  log_info("STATISTICS_detail_3D",*) 'Variable Statistics '
900  do v = 1, va
901  statval_l( v,:) = var(ks,is,js,v)
902  statidx_l(1,v,:) = ks
903  statidx_l(2,v,:) = is
904  statidx_l(3,v,:) = js
905  do j = js, je
906  do i = is, ie
907  do k = ks, ke
908  if ( var(k,i,j,v) > statval_l(v,1) ) then
909  statval_l( v,1) = var(k,i,j,v)
910  statidx_l(1,v,1) = k
911  statidx_l(2,v,1) = i
912  statidx_l(3,v,1) = j
913  end if
914  if ( var(k,i,j,v) < statval_l(v,2) ) then
915  statval_l( v,2) = var(k,i,j,v)
916  statidx_l(1,v,2) = k
917  statidx_l(2,v,2) = i
918  statidx_l(3,v,2) = j
919  end if
920  end do
921  end do
922  end do
923  enddo
924 
925  if ( do_globalcomm ) then
926  call prof_rapstart('COMM_Bcast', 2)
927 
928  call mpi_allgather( statval_l(:,:), &
929  va*2, &
930  comm_datatype, &
931  statval(:,:,:), &
932  va*2, &
933  comm_datatype, &
934  prc_local_comm_world, &
935  ierr )
936 
937  call mpi_allgather( statidx_l(:,:,:), &
938  3*va*2, &
939  mpi_integer, &
940  statidx(:,:,:,:), &
941  3*va*2, &
942  mpi_integer, &
943  prc_local_comm_world, &
944  ierr )
945 
946  call prof_rapend ('COMM_Bcast', 2)
947 
948  do v = 1, va
949  allstatval(v,1) = statval(v,1,0)
950  allstatval(v,2) = statval(v,2,0)
951  allstatidx(v,:) = 0
952  do p = 1, prc_nprocs-1
953  if ( statval(v,1,p) > allstatval(v,1) ) then
954  allstatval(v,1) = statval(v,1,p)
955  allstatidx(v,1) = p
956  end if
957  if ( statval(v,2,p) < allstatval(v,2) ) then
958  allstatval(v,2) = statval(v,2,p)
959  allstatidx(v,2) = p
960  end if
961  end do
962  log_info_cont(*) '[', trim(varname(v)), ']'
963  log_info_cont('(1x,A,ES17.10,A,4(I5,A))') ' MAX =', &
964  allstatval(v,1), ' (rank=', &
965  allstatidx(v,1), '; ', &
966  statidx(1,v,1,allstatidx(v,1)),',', &
967  statidx(2,v,1,allstatidx(v,1)),',', &
968  statidx(3,v,1,allstatidx(v,1)),')'
969  log_info_cont('(1x,A,ES17.10,A,4(I5,A))') ' MIN =', &
970  allstatval(v,2), ' (rank=', &
971  allstatidx(v,2), '; ', &
972  statidx(1,v,2,allstatidx(v,2)),',', &
973  statidx(2,v,2,allstatidx(v,2)),',', &
974  statidx(3,v,2,allstatidx(v,2)),')'
975  enddo
976  else
977  ! statistics on each node
978  do v = 1, va
979  log_info_cont(*) '[', trim(varname(v)), ']'
980  log_info_cont('(1x,A,ES17.10,A,3(I5,A))') 'MAX = ', &
981  statval_l( v,1),' (', &
982  statidx_l(1,v,1),',', &
983  statidx_l(2,v,1),',', &
984  statidx_l(3,v,1),')'
985  log_info_cont('(1x,A,ES17.10,A,3(I5,A))') 'MIN = ', &
986  statval_l( v,2),' (', &
987  statidx_l(1,v,2),',', &
988  statidx_l(2,v,2),',', &
989  statidx_l(3,v,2),')'
990  enddo
991  endif
992 
993  log_newline
994 
995  return
996  end subroutine statistics_detail_3d
997 
998  subroutine statistics_detail_2d( &
999  IA, IS, IE, JA, JS, JE, VA, &
1000  varname, var, &
1001  local )
1002  use scale_prc, only: &
1003  prc_nprocs
1004  use scale_comm_cartesc, only: &
1006  implicit none
1007 
1008  integer, intent(in) :: IA, IS, IE
1009  integer, intent(in) :: JA, JS, JE
1010  integer, intent(in) :: VA
1011 
1012  character(len=*), intent(in) :: varname(VA)
1013  real(RP), intent(in) :: var(IA,JA,VA)
1014 
1015  logical, intent(in), optional :: local ! calc in local node
1016 
1017  real(RP) :: statval_l ( VA,2)
1018  integer :: statidx_l (2,VA,2)
1019  real(RP) :: statval ( VA,2,0:PRC_nprocs-1)
1020  integer :: statidx (2,VA,2,0:PRC_nprocs-1)
1021  real(RP) :: allstatval(VA,2)
1022  integer :: allstatidx(VA,2)
1023  logical :: do_globalcomm
1024 
1025  integer :: i, j
1026  integer :: ierr
1027  integer :: v, p
1028  !---------------------------------------------------------------------------
1029 
1030  do_globalcomm = statistics_use_globalcomm
1031  if ( present(local) ) do_globalcomm = ( .not. local )
1032 
1033  !$acc update host(var) if ( acc_is_present(var) )
1034 
1035  log_newline
1036  log_info("STATISTICS_detail_2D",*) 'Variable Statistics '
1037  do v = 1, va
1038  statval_l( v,:) = var(is,js,v)
1039  statidx_l(1,v,:) = is
1040  statidx_l(2,v,:) = js
1041  do j = js, je
1042  do i = is, ie
1043  if ( var(i,j,v) > statval_l(v,1) ) then
1044  statval_l( v,1) = var(i,j,v)
1045  statidx_l(1,v,1) = i
1046  statidx_l(2,v,1) = j
1047  end if
1048  if ( var(i,j,v) < statval_l(v,2) ) then
1049  statval_l( v,2) = var(i,j,v)
1050  statidx_l(1,v,2) = i
1051  statidx_l(2,v,2) = j
1052  end if
1053  end do
1054  end do
1055  enddo
1056 
1057  if ( do_globalcomm ) then
1058  call prof_rapstart('COMM_Bcast', 2)
1059 
1060  call mpi_allgather( statval_l(:,:), &
1061  va*2, &
1062  comm_datatype, &
1063  statval(:,:,:), &
1064  va*2, &
1065  comm_datatype, &
1066  prc_local_comm_world, &
1067  ierr )
1068 
1069  call mpi_allgather( statidx_l(:,:,:), &
1070  2*va*2, &
1071  mpi_integer, &
1072  statidx(:,:,:,:), &
1073  2*va*2, &
1074  mpi_integer, &
1075  prc_local_comm_world, &
1076  ierr )
1077 
1078  call prof_rapend ('COMM_Bcast', 2)
1079 
1080  do v = 1, va
1081  allstatval(v,1) = statval(v,1,0)
1082  allstatval(v,2) = statval(v,2,0)
1083  allstatidx(v,:) = 0
1084  do p = 1, prc_nprocs-1
1085  if ( statval(v,1,p) > allstatval(v,1) ) then
1086  allstatval(v,1) = statval(v,1,p)
1087  allstatidx(v,1) = p
1088  end if
1089  if ( statval(v,2,p) < allstatval(v,2) ) then
1090  allstatval(v,2) = statval(v,2,p)
1091  allstatidx(v,2) = p
1092  end if
1093  end do
1094  log_info_cont(*) '[', trim(varname(v)), ']'
1095  log_info_cont('(1x,A,ES17.10,A,3(I5,A))') ' MAX =', &
1096  allstatval(v,1), ' (rank=', &
1097  allstatidx(v,1), '; ', &
1098  statidx(1,v,1,allstatidx(v,1)),',', &
1099  statidx(2,v,1,allstatidx(v,1)),')'
1100  log_info_cont('(1x,A,ES17.10,A,3(I5,A))') ' MIN =', &
1101  allstatval(v,2), ' (rank=', &
1102  allstatidx(v,2), '; ', &
1103  statidx(1,v,2,allstatidx(v,2)),',', &
1104  statidx(2,v,2,allstatidx(v,2)),')'
1105  enddo
1106  else
1107  ! statistics on each node
1108  do v = 1, va
1109  log_info_cont(*) '[', trim(varname(v)), ']'
1110  log_info_cont('(1x,A,ES17.10,A,2(I5,A))') 'MAX = ', &
1111  statval_l( v,1),' (', &
1112  statidx_l(1,v,1),',', &
1113  statidx_l(2,v,1),')'
1114  log_info_cont('(1x,A,ES17.10,A,2(I5,A))') 'MIN = ', &
1115  statval_l( v,2),' (', &
1116  statidx_l(1,v,2),',', &
1117  statidx_l(2,v,2),')'
1118  enddo
1119  endif
1120 
1121  log_newline
1122 
1123  return
1124  end subroutine statistics_detail_2d
1125 
1126  !-----------------------------------------------------------------------------
1127  !
1128  ! Summation (1D):
1129  ! Compenstaed summation of ARRAY (1D)
1130  !
1131  !-----------------------------------------------------------------------------
1132  function statistics_summation_1d( &
1133  IA, & ! (in)
1134  array, & ! (in)
1135  undef ) ! (in)
1136  implicit none
1137 
1138  ! arguments
1139  integer, intent(in) :: ia
1140  real(rp), intent(in) :: array(ia)
1141 
1142  real(rp), intent(in), optional :: undef
1143 
1144  ! function result
1145  real(rp) :: statistics_summation_1d
1146 
1147  ! work
1148  integer :: i
1149 
1150  real(rp) :: tmp_array(ia)
1151  real(rp) :: tmp
1152  real(rp) :: res
1153  !---------------------------------------------------------------------------
1154 
1155  if( present( undef ) ) then
1156  if( any( abs( array(:) - undef ) > eps ) ) then
1157  where( abs( array(:) - undef ) > eps )
1158  tmp_array(:) = array(:)
1159  else where
1160  tmp_array(:) = 0.0_rp
1161  end where
1162  else
1163  statistics_summation_1d = undef
1164  return ! end function
1165  end if
1166  else
1167  tmp_array(:) = array(:)
1168  end if
1169 
1170  statistics_summation_1d = 0.0_rp
1171 
1172  tmp = 0.0_rp
1173  res = 0.0_rp
1174 
1175  ! Kahan's Compensated Summation (Kahan 1965)
1176  do i = 1, ia
1177  tmp = statistics_summation_1d + ( tmp_array(i) + res )
1178  res = ( tmp_array(i) + res ) - ( tmp - statistics_summation_1d )
1180  end do
1181 
1182  return
1183  end function statistics_summation_1d
1184 
1185  !-----------------------------------------------------------------------------
1186  !
1187  ! Summation (2D):
1188  ! Compenstaed summation of ARRAY (2D)
1189  !
1190  !-----------------------------------------------------------------------------
1191  function statistics_summation_2d( &
1192  IA, JA, & ! (in)
1193  array, & ! (in)
1194  undef ) ! (in)
1195  implicit none
1196 
1197  ! arguments
1198  integer , intent(in) :: ia
1199  integer , intent(in) :: ja
1200  real(rp), intent(in) :: array(ia,ja)
1201 
1202  real(rp), intent(in), optional :: undef
1203 
1204  ! function result
1205  real(rp) :: statistics_summation_2d
1206 
1207  ! work
1208  integer :: i, j
1209 
1210  real(rp) :: tmp_array(ia,ja)
1211  real(rp) :: tmp
1212  real(rp) :: res
1213  !---------------------------------------------------------------------------
1214 
1215  if( present( undef ) ) then
1216  if( any( abs( array(:,:) - undef ) > eps ) ) then
1217  where( abs( array(:,:) - undef ) > eps )
1218  tmp_array(:,:) = array(:,:)
1219  else where
1220  tmp_array(:,:) = 0.0_rp
1221  end where
1222  else
1223  statistics_summation_2d = undef
1224  return ! end function
1225  end if
1226  else
1227  tmp_array(:,:) = array(:,:)
1228  end if
1229 
1230  statistics_summation_2d = 0.0_rp
1231 
1232  tmp = 0.0_rp
1233  res = 0.0_rp
1234 
1235  ! Kahan's Compensated Summation (Kahan 1965)
1236  do j = 1, ja
1237  do i = 1, ia
1238  tmp = statistics_summation_2d + ( tmp_array(i,j) + res )
1239  res = ( tmp_array(i,j) + res ) - ( tmp - statistics_summation_2d )
1240  statistics_summation_2d = tmp
1241  end do
1242  end do
1243 
1244  return
1245  end function statistics_summation_2d
1246 
1247  !-----------------------------------------------------------------------------
1248  !
1249  ! Summation (3D):
1250  ! Compenstaed summation of ARRAY (3D)
1251  !
1252  !-----------------------------------------------------------------------------
1253  function statistics_summation_3d( &
1254  KA, IA, JA, & ! (in)
1255  array, & ! (in)
1256  undef ) ! (in)
1257  implicit none
1258 
1259  ! arguments
1260  integer, intent(in) :: ka
1261  integer, intent(in) :: ia
1262  integer, intent(in) :: ja
1263  real(rp), intent(in) :: array(ka,ia,ja)
1264 
1265  real(rp), intent(in), optional :: undef
1266 
1267  ! function result
1268  real(rp) :: statistics_summation_3d
1269 
1270  ! work
1271  integer :: i, j, k
1272 
1273  real(rp) :: tmp_array(ka,ia,ja)
1274  real(rp) :: tmp
1275  real(rp) :: res
1276  !---------------------------------------------------------------------------
1277 
1278  if( present( undef ) ) then
1279  if( any( abs( array(:,:,:) - undef ) > eps ) ) then
1280  where( abs( array(:,:,:) - undef ) > eps )
1281  tmp_array(:,:,:) = array(:,:,:)
1282  else where
1283  tmp_array(:,:,:) = 0.0_rp
1284  end where
1285  else
1286  statistics_summation_3d = undef
1287  return ! end function
1288  end if
1289  else
1290  tmp_array(:,:,:) = array(:,:,:)
1291  end if
1292 
1293  statistics_summation_3d = 0.0_rp
1294 
1295  tmp = 0.0_rp
1296  res = 0.0_rp
1297 
1298  ! Kahan's Compensated Summation (Kahan 1965)
1299  do j = 1, ja
1300  do i = 1, ia
1301  do k = 1, ka
1302  tmp = statistics_summation_3d + ( tmp_array(i,j,k) + res )
1303  res = ( tmp_array(i,j,k) + res ) - ( tmp - statistics_summation_3d )
1304  statistics_summation_3d = tmp
1305  end do
1306  end do
1307  end do
1308 
1309  return
1310  end function statistics_summation_3d
1311 
1312  !-----------------------------------------------------------------------------
1313  !
1314  ! Average (1D):
1315  ! average of ARRAY (1D)
1316  !
1317  !-----------------------------------------------------------------------------
1318  function statistics_average_1d( &
1319  IA, & ! (in)
1320  array, & ! (in)
1321  undef ) ! (in)
1322  implicit none
1323 
1324  ! arguments
1325  integer, intent(in) :: ia
1326  real(rp), intent(in) :: array(ia)
1327 
1328  real(rp), intent(in), optional :: undef
1329 
1330  ! function result
1331  real(rp) :: statistics_average_1d
1332  !---------------------------------------------------------------------------
1333 
1334  if( present( undef ) ) then
1335  if( any( abs( array(:) - undef ) > eps ) ) then
1336  statistics_average_1d = statistics_summation( ia, array(:), undef ) &
1337  / real( count( abs( array(:) - undef ) > eps ), kind=rp )
1338  else
1339  statistics_average_1d = undef
1340  end if
1341  else
1342  statistics_average_1d = statistics_summation( ia, array(:) ) &
1343  / real( ia, kind=rp )
1344  end if
1345 
1346  return
1347  end function statistics_average_1d
1348 
1349  !-----------------------------------------------------------------------------
1350  !
1351  ! Average (2D):
1352  ! average of ARRAY (2D)
1353  !
1354  !-----------------------------------------------------------------------------
1355  function statistics_average_2d( &
1356  IA, JA, & ! (in)
1357  array, & ! (in)
1358  undef ) ! (in)
1359  implicit none
1360 
1361  ! arguments
1362  integer, intent(in) :: ia
1363  integer, intent(in) :: ja
1364  real(rp), intent(in) :: array(ia,ja)
1365 
1366  real(rp), intent(in), optional :: undef
1367 
1368  ! function result
1369  real(rp) :: statistics_average_2d
1370  !---------------------------------------------------------------------------
1371 
1372  if( present( undef ) ) then
1373  if( any( abs( array(:,:) - undef ) > eps ) ) then
1374  statistics_average_2d = statistics_summation( ia, ja, array(:,:), undef ) &
1375  / real( count( abs( array(:,:) - undef ) > eps ), kind=rp )
1376  else
1377  statistics_average_2d = undef
1378  end if
1379  else
1380  statistics_average_2d = statistics_summation( ia, ja, array(:,:) ) &
1381  / real( ia*ja, kind=rp )
1382  end if
1383 
1384  return
1385  end function statistics_average_2d
1386 
1387  !-----------------------------------------------------------------------------
1388  !
1389  ! Average (3D):
1390  ! average of ARRAY (3D)
1391  !
1392  !-----------------------------------------------------------------------------
1393  function statistics_average_3d( &
1394  KA, IA, JA, & ! (in)
1395  array, & ! (in)
1396  undef ) ! (in)
1397  implicit none
1398 
1399  ! arguments
1400  integer, intent(in) :: ka
1401  integer, intent(in) :: ia
1402  integer, intent(in) :: ja
1403  real(rp), intent(in) :: array(ka,ia,ja)
1404 
1405  real(rp), intent(in), optional :: undef
1406 
1407  ! function result
1408  real(rp) :: statistics_average_3d
1409  !---------------------------------------------------------------------------
1410 
1411  if( present( undef ) ) then
1412  if( any( abs( array(:,:,:) - undef ) > eps ) ) then
1413  statistics_average_3d = statistics_summation( ka, ia, ja, array(:,:,:), undef ) &
1414  / real( count( abs( array(:,:,:) - undef ) > eps ), kind=rp )
1415  else
1416  statistics_average_3d = undef
1417  end if
1418  else
1419  statistics_average_3d = statistics_summation( ka, ia, ja, array(:,:,:) ) &
1420  / real( ka*ia*ja, kind=rp )
1421  end if
1422 
1423  return
1424  end function statistics_average_3d
1425 
1426  !-----------------------------------------------------------------------------
1427  !
1428  ! Variance (1D):
1429  ! variance of ARRAY (1D)
1430  !
1431  !-----------------------------------------------------------------------------
1432  function statistics_variance_1d( &
1433  IA, & ! (in)
1434  array, & ! (in)
1435  undef ) ! (in)
1436  implicit none
1437 
1438  ! arguments
1439  integer, intent(in) :: ia
1440  real(rp), intent(in) :: array(ia)
1441 
1442  real(rp), intent(in), optional :: undef
1443 
1444  ! function result
1445  real(rp) :: statistics_variance_1d
1446 
1447  ! work
1448  real(rp) :: tmp_array(ia)
1449  real(rp) :: tmp
1450  !---------------------------------------------------------------------------
1451 
1452  if( present( undef ) ) then
1453  if( count( abs( array(:) - undef ) > eps ) > 1 ) then
1454  tmp = statistics_average( ia, array(:), undef )
1455 
1456  where( abs( array(:) - undef ) > eps )
1457  tmp_array(:) = ( array(:) - tmp )**2
1458  else where
1459  tmp_array(:) = undef
1460  end where
1461 
1462  statistics_variance_1d = statistics_summation( ia, tmp_array(:), undef ) &
1463  / real( count( abs( array(:) - undef ) > eps ) - 1, kind=rp )
1464  else
1465  statistics_variance_1d = undef
1466  end if
1467  else
1468  statistics_variance_1d = statistics_summation( ia, ( array(:) - statistics_average( ia, array(:) ) )**2 ) &
1469  / real( ia - 1, kind=rp )
1470  end if
1471 
1472  return
1473  end function statistics_variance_1d
1474 
1475  !-----------------------------------------------------------------------------
1476  !
1477  ! Variance (2D):
1478  ! variance of ARRAY (2D)
1479  !
1480  !-----------------------------------------------------------------------------
1481  function statistics_variance_2d( &
1482  IA, JA, & ! (in)
1483  array, & ! (in)
1484  undef ) ! (in)
1485  implicit none
1486 
1487  ! arguments
1488  integer, intent(in) :: ia
1489  integer, intent(in) :: ja
1490  real(rp), intent(in) :: array(ia,ja)
1491 
1492  real(rp), intent(in), optional :: undef
1493 
1494  ! function result
1495  real(rp) :: statistics_variance_2d
1496 
1497  ! work
1498  real(rp) :: tmp_array(ia,ja)
1499  real(rp) :: tmp
1500  !---------------------------------------------------------------------------
1501 
1502  if( present( undef ) ) then
1503  if( count( abs( array(:,:) - undef ) > eps ) > 1 ) then
1504  tmp = statistics_average( ia, ja, array(:,:), undef )
1505 
1506  where( abs( array(:,:) - undef ) > eps )
1507  tmp_array(:,:) = ( array(:,:) - tmp )**2
1508  else where
1509  tmp_array(:,:) = undef
1510  end where
1511 
1512  statistics_variance_2d = statistics_summation( ia, ja, tmp_array(:,:), undef ) &
1513  / real( count( abs( array(:,:) - undef ) > eps ) - 1, kind=rp )
1514  else
1515  statistics_variance_2d = undef
1516  end if
1517  else
1518  statistics_variance_2d = statistics_summation( ia, ja, ( array(:,:) - statistics_average( ia, ja, array(:,:) ) )**2 ) &
1519  / real( ia*ja - 1, kind=rp )
1520  end if
1521 
1522  return
1523  end function statistics_variance_2d
1524 
1525  !-----------------------------------------------------------------------------
1526  !
1527  ! Variance (3D):
1528  ! variance of ARRAY (3D)
1529  !
1530  !-----------------------------------------------------------------------------
1531  function statistics_variance_3d( &
1532  KA, IA, JA, & ! (in)
1533  array, & ! (in)
1534  undef ) ! (in)
1535  implicit none
1536 
1537  ! arguments
1538  integer, intent(in) :: ka
1539  integer, intent(in) :: ia
1540  integer, intent(in) :: ja
1541  real(rp), intent(in) :: array(ka,ia,ja)
1542 
1543  real(rp), intent(in), optional :: undef
1544 
1545  ! function result
1546  real(rp) :: statistics_variance_3d
1547 
1548  ! work
1549  real(rp) :: tmp_array(ka,ia,ja)
1550  real(rp) :: tmp
1551  !---------------------------------------------------------------------------
1552 
1553  if( present( undef ) ) then
1554  if( count( abs( array(:,:,:) - undef ) > eps ) > 1 ) then
1555  tmp = statistics_average( ka, ia, ja, array(:,:,:), undef )
1556 
1557  where( abs( array(:,:,:) - undef ) > eps )
1558  tmp_array(:,:,:) = ( array(:,:,:) - tmp )**2
1559  else where
1560  tmp_array(:,:,:) = undef
1561  end where
1562 
1563  statistics_variance_3d = statistics_summation( ka, ia, ja, tmp_array(:,:,:), undef ) &
1564  / real( count( abs( array(:,:,:) - undef ) > eps ) - 1, kind=rp )
1565  else
1566  statistics_variance_3d = undef
1567  end if
1568  else
1569  statistics_variance_3d = statistics_summation( ka, ia, ja, ( array(:,:,:) - statistics_average( ka, ia, ja, array(:,:,:) ) )**2 ) &
1570  / real( ka*ia*ja - 1, kind=rp )
1571  end if
1572 
1573  return
1574  end function statistics_variance_3d
1575 
1576  !-----------------------------------------------------------------------------
1577  !
1578  ! Standard deviation (1D):
1579  ! standard deviation of ARRAY (1D)
1580  !
1581  !-----------------------------------------------------------------------------
1582  function statistics_stddev_1d( &
1583  IA, & ! (in)
1584  array, & ! (in)
1585  undef ) ! (in)
1586  implicit none
1587 
1588  ! arguments
1589  integer, intent(in) :: ia
1590  real(rp), intent(in) :: array(ia)
1591 
1592  real(rp), intent(in), optional :: undef
1593 
1594  ! function result
1595  real(rp) :: statistics_stddev_1d
1596  !---------------------------------------------------------------------------
1597 
1598  statistics_stddev_1d = sqrt( statistics_variance( ia, array(:), undef ) )
1599 
1600  return
1601  end function statistics_stddev_1d
1602 
1603  !-----------------------------------------------------------------------------
1604  !
1605  ! Standard deviation (2D):
1606  ! standard deviation of ARRAY (2D)
1607  !
1608  !-----------------------------------------------------------------------------
1609  function statistics_stddev_2d( &
1610  IA, JA, & ! in)
1611  array, & ! (in)
1612  undef ) ! (in)
1613  implicit none
1614 
1615  ! arguments
1616  integer, intent(in) :: ia
1617  integer, intent(in) :: ja
1618  real(rp), intent(in) :: array(ia,ja)
1619 
1620  real(rp), intent(in), optional :: undef
1621 
1622  ! function result
1623  real(rp) :: statistics_stddev_2d
1624  !---------------------------------------------------------------------------
1625 
1626  statistics_stddev_2d = sqrt( statistics_variance( ia, ja, array(:,:), undef ) )
1627 
1628  return
1629  end function statistics_stddev_2d
1630 
1631  !-----------------------------------------------------------------------------
1632  !
1633  ! Standard deviation (3D):
1634  ! standard deviation of ARRAY (3D)
1635  !
1636  !-----------------------------------------------------------------------------
1637  function statistics_stddev_3d( &
1638  KA, IA, JA, & ! (in)
1639  array, & ! (in)
1640  undef ) ! (in)
1641  implicit none
1642 
1643  ! arguments
1644  integer, intent(in) :: ka
1645  integer, intent(in) :: ia
1646  integer, intent(in) :: ja
1647  real(rp), intent(in) :: array(ka,ia,ja)
1648 
1649  real(rp), intent(in), optional :: undef
1650 
1651  ! function result
1652  real(rp) :: statistics_stddev_3d
1653  !---------------------------------------------------------------------------
1654 
1655  statistics_stddev_3d = sqrt( statistics_variance( ka, ia, ja, array(:,:,:), undef ) )
1656 
1657  return
1658  end function statistics_stddev_3d
1659 
1660  !-----------------------------------------------------------------------------
1661  !
1662  ! Covariance:
1663  ! covariance between ARRAY1 and ARRAY2
1664  !
1665  !-----------------------------------------------------------------------------
1666  function statistics_covariance( &
1667  IA1, IA2, & ! (in)
1668  array1, & ! (in)
1669  array2, & ! (in)
1670  undef ) ! (in)
1671  implicit none
1672 
1673  ! argument
1674  integer, intent(in) :: ia1, ia2
1675  real(rp), intent(in) :: array1(ia1)
1676  real(rp), intent(in) :: array2(ia2)
1677 
1678  real(rp), intent(in), optional :: undef
1679 
1680  ! function result
1681  real(rp) :: statistics_covariance
1682 
1683  ! works
1684  real(rp) :: tmp_array1(ia1)
1685  real(rp) :: tmp_array2(ia2)
1686 
1687  real(rp), allocatable :: fixed_array1(:)
1688  real(rp), allocatable :: fixed_array2(:)
1689 
1690  integer :: n
1691  !---------------------------------------------------------------------------
1692 
1693  if( present( undef ) ) then
1694  n = statistics_undef_arraysize( ia1, ia2, array1(:), array2(:), undef )
1695 
1696  if( n > 1 ) then
1697  ! embed undefined value each other
1698  tmp_array1 = statistics_undef_embed( ia1, ia2, array1(:), array2(:), undef )
1699  tmp_array2 = statistics_undef_embed( ia2, ia1, array2(:), array1(:), undef )
1700 
1701  allocate( fixed_array1(n) )
1702  allocate( fixed_array2(n) )
1703 
1704  ! generate arrays without undefined value
1705  fixed_array1(:) = pack( tmp_array1(:), abs( tmp_array1(:) - undef ) > eps )
1706  fixed_array2(:) = pack( tmp_array2(:), abs( tmp_array2(:) - undef ) > eps )
1707 
1708  statistics_covariance = statistics_summation( n, ( fixed_array1(:) - statistics_average( n, array1(:), undef ) ) &
1709  * ( fixed_array2(:) - statistics_average( n, array2(:), undef ) ) ) &
1710  / real( n-1, kind=rp )
1711 
1712  deallocate( fixed_array1 )
1713  deallocate( fixed_array2 )
1714  else
1715  statistics_covariance = undef
1716  end if
1717  else
1718  n = statistics_undef_arraysize( ia1, ia2, array1(:), array2(:) )
1719  ! simple covariance
1720  statistics_covariance = statistics_summation( n, ( array1(:) - statistics_average( n, array1(:) ) ) &
1721  * ( array2(:) - statistics_average( n, array2(:) ) ) ) &
1722  / real( n-1, kind=rp )
1723  end if
1724 
1725  return
1726  end function statistics_covariance
1727 
1728  !-----------------------------------------------------------------------------
1729  !
1730  ! Correlation coefficient:
1731  ! correlation coffiecient between ARRAY1 and ARRAY2
1732  !
1733  !-----------------------------------------------------------------------------
1734  function statistics_correlation( &
1735  IA1, IA2, & ! (in)
1736  array1, & ! (in)
1737  array2, & ! (in)
1738  undef ) ! (in)
1739  implicit none
1740 
1741  ! arguments
1742  integer, intent(in) :: ia1, ia2
1743  real(rp), intent(in) :: array1(ia1)
1744  real(rp), intent(in) :: array2(ia2)
1745 
1746  real(rp), intent(in), optional :: undef
1747 
1748  ! function result
1749  real(rp) :: statistics_correlation
1750  !---------------------------------------------------------------------------
1751 
1752  if( present( undef ) ) then
1753  if( abs( statistics_covariance( ia1, ia2, array1(:), array2(:), undef ) - undef ) > eps .and. &
1754  abs( statistics_stddev( ia1, array1(:), undef ) - undef ) > eps .and. &
1755  abs( statistics_stddev( ia2, array2(:), undef ) - undef ) > eps .and. &
1756  abs( statistics_stddev( ia1, array1(:), undef ) ) > 0.0_rp .and. &
1757  abs( statistics_stddev( ia2, array2(:), undef ) ) > 0.0_rp ) then
1758  ! correlation coefficient without undefined value
1759  statistics_correlation = statistics_covariance( ia1, ia2, array1(:), array2(:), undef ) &
1760  / ( statistics_stddev( ia1, array1(:), undef ) &
1761  * statistics_stddev( ia2, array2(:), undef ) )
1762  else
1763  statistics_correlation = undef
1764  end if
1765  else
1766  ! simple correlation coefficient
1767  statistics_correlation = statistics_covariance( ia1, ia2, array1(:), array2(:) ) &
1768  / ( statistics_stddev( ia1, array1(:) ) &
1769  * statistics_stddev( ia2, array2(:) ) )
1770  end if
1771 
1772  return
1773  end function statistics_correlation
1774 
1775  !-----------------------------------------------------------------------------
1776  !
1777  ! Regression coefficient:
1778  ! regression coffiecient of ARRAY1 to ARRAY2
1779  !
1780  !-----------------------------------------------------------------------------
1781  function statistics_regression( &
1782  IA1, IA2, & ! (in)
1783  array1, & ! (in)
1784  array2, & ! (in)
1785  undef ) ! (in)
1786  implicit none
1787 
1788  ! arguments
1789  integer, intent(in) :: ia1, ia2
1790  real(rp), intent(in) :: array1(ia1)
1791  real(rp), intent(in) :: array2(ia2)
1792 
1793  real(rp), intent(in), optional :: undef
1794 
1795  ! function result
1796  real :: statistics_regression
1797  !---------------------------------------------------------------------------
1798 
1799  if( present( undef ) ) then
1800  if( abs( statistics_correlation( ia1, ia2, array1(:), array2(:), undef ) - undef ) > eps .and. &
1801  abs( statistics_stddev( ia1, array1(:), undef ) - undef ) > eps .and. &
1802  abs( statistics_stddev( ia2, array2(:), undef ) - undef ) > eps .and. &
1803  abs( statistics_stddev( ia2, array2(:), undef ) ) > 0.0_rp ) then
1804  ! regression coefficient without undefined value
1805  statistics_regression = statistics_correlation( ia1, ia2, array1(:), array2(:), undef ) &
1806  * statistics_stddev( ia1, array1(:), undef ) &
1807  / statistics_stddev( ia2, array2(:), undef )
1808  else
1809  statistics_regression = undef
1810  end if
1811  else
1812  ! simple regression coefficient
1813  statistics_regression = statistics_correlation( ia1, ia2, array1(:), array2(:) ) &
1814  * statistics_stddev( ia1, array1(:) ) &
1815  / statistics_stddev( ia2, array2(:) )
1816  end if
1817 
1818  return
1819  end function statistics_regression
1820 
1821  !-----------------------------------------------------------------------------
1822  !
1823  ! Lag correlation coefficient:
1824  ! correlation between ARRAY1(t) and ARRAY2(t+LAG)
1825  !
1826  !-----------------------------------------------------------------------------
1827  function statistics_lag_correlation( &
1828  IA1, IA2, & ! (in)
1829  array1, & ! (in)
1830  array2, & ! (in)
1831  lag, & ! (in)
1832  undef ) ! (in)
1833  implicit none
1834 
1835  ! arguments
1836  integer, intent(in) :: ia1, ia2
1837  real(rp), intent(in) :: array1(ia1)
1838  real(rp), intent(in) :: array2(ia2)
1839 
1840  integer, intent(in) :: lag
1841 
1842  real(rp), intent(in), optional :: undef
1843 
1844  ! function result
1846 
1847  ! works
1848  real(rp) :: tmp(ia2)
1849  !---------------------------------------------------------------------------
1850 
1851  tmp(:) = eoshift( array2(:), lag, undef )
1852 
1853  statistics_lag_correlation = statistics_correlation( ia1, ia2, array1(:), tmp(:), undef )
1854 
1855  return
1856  end function statistics_lag_correlation
1857 
1858  !-----------------------------------------------------------------------------
1859  !
1860  ! Partial correlation coefficient:
1861  ! correlation between ARRAY1 and ARRAY2 excepting influence of ARRAY3
1862  !
1863  !-----------------------------------------------------------------------------
1864  function statistics_partial_correlation( &
1865  IA1, IA2, IA3, & ! (in)
1866  array1, & ! (in)
1867  array2, & ! (in)
1868  array3, & ! (in)
1869  undef ) ! (in)
1870  implicit none
1871 
1872  ! arguments
1873  integer, intent(in) :: ia1, ia2, ia3
1874  real(rp), intent(in) :: array1(ia1)
1875  real(rp), intent(in) :: array2(ia2)
1876  real(rp), intent(in) :: array3(ia3)
1877 
1878  real(rp), intent(in), optional :: undef
1879 
1880  ! function result
1882 
1883  ! works
1884  real(rp) :: corr12
1885  real(rp) :: corr13
1886  real(rp) :: corr23
1887  !---------------------------------------------------------------------------
1888 
1889  if( present( undef ) ) then
1890  corr12 = statistics_correlation( ia1, ia2, array1(:), array2(:), undef )
1891  corr13 = statistics_correlation( ia1, ia3, array1(:), array3(:), undef )
1892  corr23 = statistics_correlation( ia2, ia3, array2(:), array3(:), undef )
1893 
1894  if( corr12 == undef .or. &
1895  corr13 == undef .or. &
1896  corr23 == undef ) then
1897  ! return undefined value if correlation can't calculate
1899  else if( corr12 >= 1.0_rp .or. &
1900  corr13 >= 1.0_rp .or. &
1901  corr23 >= 1.0_rp ) then
1902  ! return undefined value if correlation is 1
1904  else
1905  ! partial correlation without undefined value
1906  statistics_partial_correlation = ( corr12 - ( corr13 * corr23 ) ) &
1907  / sqrt( ( 1.0_rp - corr13**2 ) * ( 1.0_rp - corr23**2 ) )
1908  endif
1909  else
1910  corr12 = statistics_correlation( ia1, ia2, array1(:), array2(:) )
1911  corr13 = statistics_correlation( ia1, ia3, array1(:), array3(:) )
1912  corr23 = statistics_correlation( ia2, ia3, array2(:), array3(:) )
1913 
1914  ! simple partial correlation
1915  statistics_partial_correlation = ( corr12 - ( corr13 * corr23 ) ) &
1916  / sqrt( ( 1.0_rp - corr13**2 ) * ( 1.0_rp - corr23**2 ) )
1917  end if
1918 
1919  return
1920  end function statistics_partial_correlation
1921 
1922  !-----------------------------------------------------------------------------
1923  function statistics_undef_replace( &
1924  IA, & ! (in)
1925  array, & ! (in)
1926  undef1, & ! (in)
1927  undef2 ) ! (in)
1928  implicit none
1929 
1930  ! arguments
1931  integer, intent(in) :: ia
1932  real(rp), intent(in) :: array(ia)
1933  real(rp), intent(in) :: undef1
1934  real(rp), intent(in) :: undef2
1935 
1936  ! function result
1937  real(rp) :: statistics_undef_replace(ia)
1938 
1939  ! works
1940  integer :: i
1941  !---------------------------------------------------------------------------
1942 
1943  do i = 1, ia
1944  if( abs( array(i) - undef2 ) > eps ) then
1945  statistics_undef_replace(i) = array(i)
1946  else
1947  statistics_undef_replace(i) = undef1
1948  end if
1949  end do
1950 
1951  return
1952  end function statistics_undef_replace
1953 
1954  !-----------------------------------------------------------------------------
1955  function statistics_undef_embed( &
1956  IA1, IA2, & ! (in)
1957  array1, & ! (in)
1958  array2, & ! (in)
1959  undef ) ! (in)
1960  implicit none
1961 
1962  ! arguments
1963  integer, intent(in) :: ia1, ia2
1964  real(rp), intent(in) :: array1(ia1)
1965  real(rp), intent(in) :: array2(ia2)
1966 
1967  real(rp), intent(in), optional :: undef
1968 
1969  ! function result
1970  real(rp) :: statistics_undef_embed(ia1)
1971 
1972  ! works
1973  integer :: i
1974  !---------------------------------------------------------------------------
1975 
1976  ! check the length of input arrays
1977  if( ia1 /= ia2 ) then
1978  write(*,*) "Error: different array size !!"
1979  stop
1980  end if
1981 
1982  do i = 1, ia1
1983  if( abs( array2(i) - undef ) > eps ) then
1984  statistics_undef_embed(i) = array1(i)
1985  else
1986  statistics_undef_embed(i) = undef
1987  end if
1988  end do
1989  end function statistics_undef_embed
1990 
1991  !-----------------------------------------------------------------------------
1992  function statistics_undef_arraysize( &
1993  IA1, IA2, & ! (in)
1994  array1, & ! (in)
1995  array2, & ! (in)
1996  undef ) ! (in)
1997  implicit none
1998 
1999  ! arguments
2000  integer, intent(in) :: ia1, ia2
2001  real(rp), intent(in) :: array1(ia1)
2002  real(rp), intent(in) :: array2(ia2)
2003 
2004  real(rp), intent(in), optional :: undef
2005 
2006  ! function result
2007  integer :: statistics_undef_arraysize
2008 
2009  ! works
2010  integer :: i, cnt
2011  !---------------------------------------------------------------------------
2012 
2013  ! check the length of input arrays
2014  if( ia1 /= ia2 ) then
2015  write(*,*) "Error: different array size !!"
2016  stop
2017  end if
2018 
2019  if( present( undef ) ) then
2020  cnt = 0
2021  do i = 1, ia1
2022  if( abs( array1(i) - undef ) > eps .and. &
2023  abs( array2(i) - undef ) > eps ) then
2024  cnt = cnt + 1
2025  end if
2026  end do
2027  ! array size excepting number of undefined values
2029  else
2030  ! simple array size
2032  end if
2033 
2034  return
2035  end function statistics_undef_arraysize
2036 
2037 end module scale_statistics
scale_comm_cartesc::comm_setup
subroutine, public comm_setup
Setup.
Definition: scale_comm_cartesC.F90:200
scale_statistics
module Statistics
Definition: scale_statistics.F90:11
scale_comm_cartesc::comm_datatype
integer, public comm_datatype
datatype of variable
Definition: scale_comm_cartesC.F90:105
scale_statistics::statistics_horizontal_min_2d
subroutine statistics_horizontal_min_2d(IA, IS, IE, JA, JS, JE, var, varmin)
Calc horizontal minimum value.
Definition: scale_statistics.F90:568
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_statistics::statistics_horizontal_max_2d
subroutine statistics_horizontal_max_2d(IA, IS, IE, JA, JS, JE, var, varmax)
Calc horizontal maximum value.
Definition: scale_statistics.F90:716
scale_statistics::statistics_partial_correlation
real(rp) function, public statistics_partial_correlation(IA1, IA2, IA3, ARRAY1, ARRAY2, ARRAY3, UNDEF)
Definition: scale_statistics.F90:1870
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_prc::prc_local_comm_world
integer, public prc_local_comm_world
local communicator
Definition: scale_prc.F90:89
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:35
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level, disable_barrier)
Start raptime.
Definition: scale_prof.F90:174
scale_prc::prc_myrank
integer, public prc_myrank
process num in local communicator
Definition: scale_prc.F90:91
scale_statistics::statistics_undef_replace
real(rp) function, dimension(ia), public statistics_undef_replace(IA, ARRAY, UNDEF1, UNDEF2)
Definition: scale_statistics.F90:1928
scale_statistics::statistics_correlation
real(rp) function, public statistics_correlation(IA1, IA2, ARRAY1, ARRAY2, UNDEF)
Definition: scale_statistics.F90:1739
scale_statistics::statistics_covariance
real(rp) function, public statistics_covariance(IA1, IA2, ARRAY1, ARRAY2, UNDEF)
Definition: scale_statistics.F90:1671
scale_statistics::statistics_total_3d
subroutine statistics_total_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, var, varname, vol, total, log_suppress, global, mean, sum)
Calc domain sum and volume-weighted mean.
Definition: scale_statistics.F90:293
scale_statistics::statistics_horizontal_mean_2d
subroutine statistics_horizontal_mean_2d(IA, IS, IE, JA, JS, JE, var, area, varmean)
Calc horizontal mean value.
Definition: scale_statistics.F90:423
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_statistics::statistics_horizontal_min_3d
subroutine statistics_horizontal_min_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, var, varmin)
Definition: scale_statistics.F90:624
scale_io
module STDIO
Definition: scale_io.F90:10
scale_statistics::statistics_undef_arraysize
integer function, public statistics_undef_arraysize(IA1, IA2, ARRAY1, ARRAY2, UNDEF)
Definition: scale_statistics.F90:1997
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_statistics::statistics_lag_correlation
real(rp) function, public statistics_lag_correlation(IA1, IA2, ARRAY1, ARRAY2, LAG, UNDEF)
Definition: scale_statistics.F90:1833
scale_statistics::statistics_summation_1d
real(rp) function statistics_summation_1d(IA, ARRAY, UNDEF)
Definition: scale_statistics.F90:1136
scale_const::const_huge
real(rp), public const_huge
huge number
Definition: scale_const.F90:37
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_statistics::statistics_regression
real function, public statistics_regression(IA1, IA2, ARRAY1, ARRAY2, UNDEF)
Definition: scale_statistics.F90:1786
scale_statistics::statistics_horizontal_mean_3d
subroutine statistics_horizontal_mean_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, var, area, varmean)
Definition: scale_statistics.F90:481
scale_prc::prc_nprocs
integer, public prc_nprocs
myrank in local communicator
Definition: scale_prc.F90:90
scale_statistics::statistics_detail_2d
subroutine statistics_detail_2d(IA, IS, IE, JA, JS, JE, VA, varname, var, local)
Definition: scale_statistics.F90:1002
scale_statistics::statistics_checktotal
logical, public statistics_checktotal
calc&report variable totals to logfile?
Definition: scale_statistics.F90:109
scale_comm_cartesc
module COMMUNICATION
Definition: scale_comm_cartesC.F90:11
scale_statistics::statistics_undef_embed
real(rp) function, dimension(ia1), public statistics_undef_embed(IA1, IA2, ARRAY1, ARRAY2, UNDEF)
Definition: scale_statistics.F90:1960
scale_statistics::statistics_total_2d
subroutine statistics_total_2d(IA, IS, IE, JA, JS, JE, var, varname, area, total, log_suppress, global, mean, sum)
Calc domain sum and area-weighted mean.
Definition: scale_statistics.F90:175
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level, disable_barrier)
Save raptime.
Definition: scale_prof.F90:246
scale_statistics::statistics_setup
subroutine, public statistics_setup
Setup.
Definition: scale_statistics.F90:126
scale_statistics::statistics_horizontal_max_3d
subroutine statistics_horizontal_max_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, var, varmax)
Definition: scale_statistics.F90:772
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_statistics::statistics_detail_3d
subroutine statistics_detail_3d(KA, KS, KE, IA, IS, IE, JA, JS, JE, VA, varname, var, local)
Search global maximum & minimum value.
Definition: scale_statistics.F90:864
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:57