35 public :: statistics_total
36 public :: statistics_detail
37 public :: statistics_horizontal_mean
38 public :: statistics_horizontal_min
39 public :: statistics_horizontal_max
41 public :: statistics_summation
42 public :: statistics_average
43 public :: statistics_variance
44 public :: statistics_stddev
56 interface statistics_total
59 end interface statistics_total
61 interface statistics_detail
64 end interface statistics_detail
66 interface statistics_horizontal_mean
69 end interface statistics_horizontal_mean
71 interface statistics_horizontal_max
74 end interface statistics_horizontal_max
76 interface statistics_horizontal_min
79 end interface statistics_horizontal_min
81 interface statistics_summation
83 module procedure statistics_summation_2d
84 module procedure statistics_summation_3d
85 end interface statistics_summation
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
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
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
119 logical,
private :: statistics_use_globalcomm = .false.
132 namelist / param_statistics / &
134 statistics_use_globalcomm
142 log_info(
"STATISTICS_setup",*)
'Setup'
148 log_info(
"STATISTICS_setup",*)
'Not found namelist. Default used.'
149 elseif( ierr > 0 )
then
150 log_error(
"STATISTICS_setup",*)
'Not appropriate names in namelist PARAM_STATISTICS. Check!'
153 log_nml(param_statistics)
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.'
160 log_info_cont(*)
'=> The total is calculated for the local domain.'
169 IA, IS, IE, JA, JS, JE, &
183 integer,
intent(in) :: IA, IS, IE
184 integer,
intent(in) :: JA, JS, JE
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
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
197 real(DP) :: sendbuf(2), recvbuf(2)
198 real(DP) :: sum_, mean_
200 logical :: suppress_, global_
207 if ( var(is,js) /= undef )
then
214 statval = statval + var(i,j) * area(i,j)
220 if ( .NOT. ( statval > -1.0_dp .OR. statval < 1.0_dp ) )
then
221 log_error(
"STATISTICS_total_2D",*)
'NaN is detected for ', trim(varname),
' in rank ',
prc_myrank
225 if (
present(log_suppress) )
then
226 suppress_ = log_suppress
231 if (
present(global) )
then
234 global_ = statistics_use_globalcomm
242 call mpi_allreduce( sendbuf(:), recvbuf(:), &
244 mpi_double_precision, &
246 prc_local_comm_world, &
250 if ( recvbuf(2) < eps )
then
255 mean_ = recvbuf(1) / recvbuf(2)
258 if ( .not. suppress_ )
then
259 log_info(
"STATISTICS_total_2D",
'(1x,A,A24,A,ES24.17)') &
260 '[', trim(varname),
'] MEAN(global) = ', mean_
263 if ( total < eps )
then
268 mean_ = statval / total
272 if ( .not. suppress_ )
then
273 log_info(
"STATISTICS_total_2D",
'(1x,A,A24,A,ES24.17)') &
274 '[', trim(varname),
'] MEAN(local) = ', mean_
278 if (
present(mean) ) mean = mean_
279 if (
present(sum ) ) sum = sum_
287 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
301 integer,
intent(in) :: KA, KS, KE
302 integer,
intent(in) :: IA, IS, IE
303 integer,
intent(in) :: JA, JS, JE
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
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
316 real(DP) :: sendbuf(2), recvbuf(2)
317 real(DP) :: mean_, sum_
319 logical :: suppress_, global_
328 if ( var(ks,is,js) /= undef )
then
339 statval = statval + var(k,i,j) * vol(k,i,j)
344 work = work + var(k,i,j) * vol(k,i,j)
346 statval = statval + work
353 if ( .NOT. ( statval > -1.0_dp .OR. statval < 1.0_dp ) )
then
354 log_error(
"STATISTICS_total_3D",*)
'NaN is detected for ', trim(varname),
' in rank ',
prc_myrank
358 if (
present(log_suppress) )
then
359 suppress_ = log_suppress
364 if (
present(global) )
then
367 global_ = statistics_use_globalcomm
375 call mpi_allreduce( sendbuf(:), recvbuf(:), &
377 mpi_double_precision, &
379 prc_local_comm_world, &
383 if ( recvbuf(2) < eps )
then
388 mean_ = recvbuf(1) / recvbuf(2)
391 if ( .not. suppress_ )
then
392 log_info(
"STATISTICS_total_3D",
'(1x,A,A24,A,ES24.17)') &
393 '[', trim(varname),
'] MEAN(global) = ', mean_
396 if ( total < eps )
then
401 mean_ = statval / total
405 if ( .not. suppress_ )
then
406 log_info(
"STATISTICS_total_3D",
'(1x,A,A24,A,ES24.17)') &
407 '[', trim(varname),
'] MEAN(local) = ', mean_
411 if (
present(mean) ) mean = mean_
412 if (
present(sum ) ) sum = sum_
420 IA, IS, IE, JA, JS, JE, &
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
431 real(DP) :: statval (2)
432 real(DP) :: allstatval(2)
447 if ( var(i,j) /= undef )
then
448 s1 = s1 + area(i,j) * var(i,j)
459 call mpi_allreduce( statval(:), &
462 mpi_double_precision, &
468 if ( allstatval(2) > 0.0_dp )
then
469 varmean = allstatval(1) / allstatval(2)
478 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
483 integer,
intent(in) :: KA, KS, KE
484 integer,
intent(in) :: IA, IS, IE
485 integer,
intent(in) :: JA, JS, JE
487 real(RP),
intent(in) :: var (KA,IA,JA)
488 real(RP),
intent(in) :: area( IA,JA)
489 real(RP),
intent(out) :: varmean(KA)
491 real(DP) :: statval (KS:KE,2)
492 real(DP) :: allstatval(KS:KE,2)
501 statval(:,:) = 0.0_dp
511 if ( var(k,i,j) /= undef )
then
513 statval(k,1) = statval(k,1) + area(i,j) * var(k,i,j)
516 statval(k,2) = statval(k,2) + area(i,j)
527 call mpi_allreduce( statval(:,:), &
530 mpi_double_precision, &
539 if ( allstatval(k,2) > 0.0_dp )
then
540 varmean(k) = allstatval(k,1) / allstatval(k,2)
565 IA, IS, IE, JA, JS, JE, &
573 integer,
intent(in) :: IA, IS, IE
574 integer,
intent(in) :: JA, JS, JE
576 real(RP),
intent(in) :: var(IA,JA)
577 real(RP),
intent(out) :: varmin
580 real(RP) :: allstatval
593 if ( var(i,j) /= undef .and. var(i,j) < statval )
then
602 call mpi_allreduce( statval, &
611 if ( allstatval < huge )
then
621 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
629 integer,
intent(in) :: KA, KS, KE
630 integer,
intent(in) :: IA, IS, IE
631 integer,
intent(in) :: JA, JS, JE
633 real(RP),
intent(in) :: var(KA,IA,JA)
634 real(RP),
intent(out) :: varmin(KA)
636 real(RP) :: statval (KA)
637 real(RP) :: allstatval(KA)
657 if ( var(k,i,j) /= undef )
then
659 statval(k) = min( var(k,i,j), statval(k) )
663 if ( var(k,i,j) /= undef .and. var(k,i,j) < statval(k) )
then
664 statval(k) = var(k,i,j)
675 call mpi_allreduce( statval(ks:ke), &
687 if ( allstatval(k) < huge )
then
688 varmin(k) = allstatval(k)
713 IA, IS, IE, JA, JS, JE, &
721 integer,
intent(in) :: IA, IS, IE
722 integer,
intent(in) :: JA, JS, JE
724 real(RP),
intent(in) :: var(IA,JA)
725 real(RP),
intent(out) :: varmax
728 real(RP) :: allstatval
741 if ( var(i,j) /= undef .and. var(i,j) > statval )
then
750 call mpi_allreduce( statval, &
759 if ( allstatval > - huge )
then
769 KA, KS, KE, IA, IS, IE, JA, JS, JE, &
777 integer,
intent(in) :: KA, KS, KE
778 integer,
intent(in) :: IA, IS, IE
779 integer,
intent(in) :: JA, JS, JE
781 real(RP),
intent(in) :: var(KA,IA,JA)
782 real(RP),
intent(out) :: varmax(KA)
784 real(RP) :: statval (KA)
785 real(RP) :: allstatval(KA)
805 if ( var(k,i,j) /= undef )
then
807 statval(k) = max( var(k,i,j), statval(k) )
811 if ( var(k,i,j) /= undef .and. var(k,i,j) > statval(k) )
then
812 statval(k) = var(k,i,j)
823 call mpi_allreduce( statval(ks:ke), &
835 if ( allstatval(k) > - huge )
then
836 varmax(k) = allstatval(k)
861 KA, KS, KE, IA, IS, IE, JA, JS, JE, VA, &
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
875 character(len=*),
intent(in) :: varname(VA)
876 real(RP),
intent(in) :: var(KA,IA,JA,VA)
878 logical,
intent(in),
optional :: local
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
893 do_globalcomm = statistics_use_globalcomm
894 if (
present(local) ) do_globalcomm = ( .not. local )
899 log_info(
"STATISTICS_detail_3D",*)
'Variable Statistics '
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
908 if ( var(k,i,j,v) > statval_l(v,1) )
then
909 statval_l( v,1) = var(k,i,j,v)
914 if ( var(k,i,j,v) < statval_l(v,2) )
then
915 statval_l( v,2) = var(k,i,j,v)
925 if ( do_globalcomm )
then
928 call mpi_allgather( statval_l(:,:), &
934 prc_local_comm_world, &
937 call mpi_allgather( statidx_l(:,:,:), &
943 prc_local_comm_world, &
949 allstatval(v,1) = statval(v,1,0)
950 allstatval(v,2) = statval(v,2,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)
957 if ( statval(v,2,p) < allstatval(v,2) )
then
958 allstatval(v,2) = statval(v,2,p)
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)),
')'
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),
',', &
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),
',', &
999 IA, IS, IE, JA, JS, JE, VA, &
1008 integer,
intent(in) :: IA, IS, IE
1009 integer,
intent(in) :: JA, JS, JE
1010 integer,
intent(in) :: VA
1012 character(len=*),
intent(in) :: varname(VA)
1013 real(RP),
intent(in) :: var(IA,JA,VA)
1015 logical,
intent(in),
optional :: local
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
1030 do_globalcomm = statistics_use_globalcomm
1031 if (
present(local) ) do_globalcomm = ( .not. local )
1036 log_info(
"STATISTICS_detail_2D",*)
'Variable Statistics '
1038 statval_l( v,:) = var(is,js,v)
1039 statidx_l(1,v,:) = is
1040 statidx_l(2,v,:) = js
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
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
1057 if ( do_globalcomm )
then
1060 call mpi_allgather( statval_l(:,:), &
1066 prc_local_comm_world, &
1069 call mpi_allgather( statidx_l(:,:,:), &
1075 prc_local_comm_world, &
1081 allstatval(v,1) = statval(v,1,0)
1082 allstatval(v,2) = statval(v,2,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)
1089 if ( statval(v,2,p) < allstatval(v,2) )
then
1090 allstatval(v,2) = statval(v,2,p)
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)),
')'
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),
')'
1139 integer,
intent(in) :: ia
1140 real(
rp),
intent(in) :: array(ia)
1142 real(
rp),
intent(in),
optional :: undef
1150 real(
rp) :: tmp_array(ia)
1155 if(
present( undef ) )
then
1156 if( any( abs( array(:) - undef ) > eps ) )
then
1157 where( abs( array(:) - undef ) > eps )
1158 tmp_array(:) = array(:)
1160 tmp_array(:) = 0.0_rp
1167 tmp_array(:) = array(:)
1191 function statistics_summation_2d( &
1198 integer ,
intent(in) :: ia
1199 integer ,
intent(in) :: ja
1200 real(
rp),
intent(in) :: array(ia,ja)
1202 real(
rp),
intent(in),
optional :: undef
1205 real(
rp) :: statistics_summation_2d
1210 real(
rp) :: tmp_array(ia,ja)
1215 if(
present( undef ) )
then
1216 if( any( abs( array(:,:) - undef ) > eps ) )
then
1217 where( abs( array(:,:) - undef ) > eps )
1218 tmp_array(:,:) = array(:,:)
1220 tmp_array(:,:) = 0.0_rp
1223 statistics_summation_2d = undef
1227 tmp_array(:,:) = array(:,:)
1230 statistics_summation_2d = 0.0_rp
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
1245 end function statistics_summation_2d
1253 function statistics_summation_3d( &
1254 KA, IA, JA, & ! (in)
1260 integer,
intent(in) :: ka
1261 integer,
intent(in) :: ia
1262 integer,
intent(in) :: ja
1263 real(
rp),
intent(in) :: array(ka,ia,ja)
1265 real(
rp),
intent(in),
optional :: undef
1268 real(
rp) :: statistics_summation_3d
1273 real(
rp) :: tmp_array(ka,ia,ja)
1278 if(
present( undef ) )
then
1279 if( any( abs( array(:,:,:) - undef ) > eps ) )
then
1280 where( abs( array(:,:,:) - undef ) > eps )
1281 tmp_array(:,:,:) = array(:,:,:)
1283 tmp_array(:,:,:) = 0.0_rp
1286 statistics_summation_3d = undef
1290 tmp_array(:,:,:) = array(:,:,:)
1293 statistics_summation_3d = 0.0_rp
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
1310 end function statistics_summation_3d
1318 function statistics_average_1d( &
1325 integer,
intent(in) :: ia
1326 real(
rp),
intent(in) :: array(ia)
1328 real(
rp),
intent(in),
optional :: undef
1331 real(
rp) :: statistics_average_1d
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 )
1339 statistics_average_1d = undef
1342 statistics_average_1d = statistics_summation( ia, array(:) ) &
1343 / real( ia, kind=
rp )
1347 end function statistics_average_1d
1355 function statistics_average_2d( &
1362 integer,
intent(in) :: ia
1363 integer,
intent(in) :: ja
1364 real(
rp),
intent(in) :: array(ia,ja)
1366 real(
rp),
intent(in),
optional :: undef
1369 real(
rp) :: statistics_average_2d
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 )
1377 statistics_average_2d = undef
1380 statistics_average_2d = statistics_summation( ia, ja, array(:,:) ) &
1381 / real( ia*ja, kind=
rp )
1385 end function statistics_average_2d
1393 function statistics_average_3d( &
1394 KA, IA, JA, & ! (in)
1400 integer,
intent(in) :: ka
1401 integer,
intent(in) :: ia
1402 integer,
intent(in) :: ja
1403 real(
rp),
intent(in) :: array(ka,ia,ja)
1405 real(
rp),
intent(in),
optional :: undef
1408 real(
rp) :: statistics_average_3d
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 )
1416 statistics_average_3d = undef
1419 statistics_average_3d = statistics_summation( ka, ia, ja, array(:,:,:) ) &
1420 / real( ka*ia*ja, kind=
rp )
1424 end function statistics_average_3d
1432 function statistics_variance_1d( &
1439 integer,
intent(in) :: ia
1440 real(
rp),
intent(in) :: array(ia)
1442 real(
rp),
intent(in),
optional :: undef
1445 real(
rp) :: statistics_variance_1d
1448 real(
rp) :: tmp_array(ia)
1452 if(
present( undef ) )
then
1453 if( count( abs( array(:) - undef ) > eps ) > 1 )
then
1454 tmp = statistics_average( ia, array(:), undef )
1456 where( abs( array(:) - undef ) > eps )
1457 tmp_array(:) = ( array(:) - tmp )**2
1459 tmp_array(:) = undef
1462 statistics_variance_1d = statistics_summation( ia, tmp_array(:), undef ) &
1463 / real( count( abs( array(:) - undef ) > eps ) - 1, kind=
rp )
1465 statistics_variance_1d = undef
1468 statistics_variance_1d = statistics_summation( ia, ( array(:) - statistics_average( ia, array(:) ) )**2 ) &
1469 / real( ia - 1, kind=
rp )
1473 end function statistics_variance_1d
1481 function statistics_variance_2d( &
1488 integer,
intent(in) :: ia
1489 integer,
intent(in) :: ja
1490 real(
rp),
intent(in) :: array(ia,ja)
1492 real(
rp),
intent(in),
optional :: undef
1495 real(
rp) :: statistics_variance_2d
1498 real(
rp) :: tmp_array(ia,ja)
1502 if(
present( undef ) )
then
1503 if( count( abs( array(:,:) - undef ) > eps ) > 1 )
then
1504 tmp = statistics_average( ia, ja, array(:,:), undef )
1506 where( abs( array(:,:) - undef ) > eps )
1507 tmp_array(:,:) = ( array(:,:) - tmp )**2
1509 tmp_array(:,:) = undef
1512 statistics_variance_2d = statistics_summation( ia, ja, tmp_array(:,:), undef ) &
1513 / real( count( abs( array(:,:) - undef ) > eps ) - 1, kind=
rp )
1515 statistics_variance_2d = undef
1518 statistics_variance_2d = statistics_summation( ia, ja, ( array(:,:) - statistics_average( ia, ja, array(:,:) ) )**2 ) &
1519 / real( ia*ja - 1, kind=
rp )
1523 end function statistics_variance_2d
1531 function statistics_variance_3d( &
1532 KA, IA, JA, & ! (in)
1538 integer,
intent(in) :: ka
1539 integer,
intent(in) :: ia
1540 integer,
intent(in) :: ja
1541 real(
rp),
intent(in) :: array(ka,ia,ja)
1543 real(
rp),
intent(in),
optional :: undef
1546 real(
rp) :: statistics_variance_3d
1549 real(
rp) :: tmp_array(ka,ia,ja)
1553 if(
present( undef ) )
then
1554 if( count( abs( array(:,:,:) - undef ) > eps ) > 1 )
then
1555 tmp = statistics_average( ka, ia, ja, array(:,:,:), undef )
1557 where( abs( array(:,:,:) - undef ) > eps )
1558 tmp_array(:,:,:) = ( array(:,:,:) - tmp )**2
1560 tmp_array(:,:,:) = undef
1563 statistics_variance_3d = statistics_summation( ka, ia, ja, tmp_array(:,:,:), undef ) &
1564 / real( count( abs( array(:,:,:) - undef ) > eps ) - 1, kind=
rp )
1566 statistics_variance_3d = undef
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 )
1574 end function statistics_variance_3d
1582 function statistics_stddev_1d( &
1589 integer,
intent(in) :: ia
1590 real(
rp),
intent(in) :: array(ia)
1592 real(
rp),
intent(in),
optional :: undef
1595 real(
rp) :: statistics_stddev_1d
1598 statistics_stddev_1d = sqrt( statistics_variance( ia, array(:), undef ) )
1601 end function statistics_stddev_1d
1609 function statistics_stddev_2d( &
1616 integer,
intent(in) :: ia
1617 integer,
intent(in) :: ja
1618 real(
rp),
intent(in) :: array(ia,ja)
1620 real(
rp),
intent(in),
optional :: undef
1623 real(
rp) :: statistics_stddev_2d
1626 statistics_stddev_2d = sqrt( statistics_variance( ia, ja, array(:,:), undef ) )
1629 end function statistics_stddev_2d
1637 function statistics_stddev_3d( &
1638 KA, IA, JA, & ! (in)
1644 integer,
intent(in) :: ka
1645 integer,
intent(in) :: ia
1646 integer,
intent(in) :: ja
1647 real(
rp),
intent(in) :: array(ka,ia,ja)
1649 real(
rp),
intent(in),
optional :: undef
1652 real(
rp) :: statistics_stddev_3d
1655 statistics_stddev_3d = sqrt( statistics_variance( ka, ia, ja, array(:,:,:), undef ) )
1658 end function statistics_stddev_3d
1674 integer,
intent(in) :: ia1, ia2
1675 real(
rp),
intent(in) :: array1(ia1)
1676 real(
rp),
intent(in) :: array2(ia2)
1678 real(
rp),
intent(in),
optional :: undef
1684 real(
rp) :: tmp_array1(ia1)
1685 real(
rp) :: tmp_array2(ia2)
1687 real(
rp),
allocatable :: fixed_array1(:)
1688 real(
rp),
allocatable :: fixed_array2(:)
1693 if(
present( undef ) )
then
1701 allocate( fixed_array1(n) )
1702 allocate( fixed_array2(n) )
1705 fixed_array1(:) = pack( tmp_array1(:), abs( tmp_array1(:) - undef ) > eps )
1706 fixed_array2(:) = pack( tmp_array2(:), abs( tmp_array2(:) - undef ) > eps )
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 )
1712 deallocate( fixed_array1 )
1713 deallocate( fixed_array2 )
1720 statistics_covariance = statistics_summation( n, ( array1(:) - statistics_average( n, array1(:) ) ) &
1721 * ( array2(:) - statistics_average( n, array2(:) ) ) ) &
1722 / real( n-1, kind=
rp )
1742 integer,
intent(in) :: ia1, ia2
1743 real(
rp),
intent(in) :: array1(ia1)
1744 real(
rp),
intent(in) :: array2(ia2)
1746 real(
rp),
intent(in),
optional :: undef
1752 if(
present( undef ) )
then
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
1760 / ( statistics_stddev( ia1, array1(:), undef ) &
1761 * statistics_stddev( ia2, array2(:), undef ) )
1768 / ( statistics_stddev( ia1, array1(:) ) &
1769 * statistics_stddev( ia2, array2(:) ) )
1789 integer,
intent(in) :: ia1, ia2
1790 real(
rp),
intent(in) :: array1(ia1)
1791 real(
rp),
intent(in) :: array2(ia2)
1793 real(
rp),
intent(in),
optional :: undef
1799 if(
present( undef ) )
then
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
1806 * statistics_stddev( ia1, array1(:), undef ) &
1807 / statistics_stddev( ia2, array2(:), undef )
1814 * statistics_stddev( ia1, array1(:) ) &
1815 / statistics_stddev( ia2, array2(:) )
1836 integer,
intent(in) :: ia1, ia2
1837 real(
rp),
intent(in) :: array1(ia1)
1838 real(
rp),
intent(in) :: array2(ia2)
1840 integer,
intent(in) :: lag
1842 real(
rp),
intent(in),
optional :: undef
1848 real(
rp) :: tmp(ia2)
1851 tmp(:) = eoshift( array2(:), lag, undef )
1865 IA1, IA2, IA3, & ! (in)
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)
1878 real(
rp),
intent(in),
optional :: undef
1889 if(
present( undef ) )
then
1894 if( corr12 == undef .or. &
1895 corr13 == undef .or. &
1896 corr23 == undef )
then
1899 else if( corr12 >= 1.0_rp .or. &
1900 corr13 >= 1.0_rp .or. &
1901 corr23 >= 1.0_rp )
then
1907 / sqrt( ( 1.0_rp - corr13**2 ) * ( 1.0_rp - corr23**2 ) )
1916 / sqrt( ( 1.0_rp - corr13**2 ) * ( 1.0_rp - corr23**2 ) )
1931 integer,
intent(in) :: ia
1932 real(
rp),
intent(in) :: array(ia)
1933 real(
rp),
intent(in) :: undef1
1934 real(
rp),
intent(in) :: undef2
1944 if( abs( array(i) - undef2 ) > eps )
then
1963 integer,
intent(in) :: ia1, ia2
1964 real(
rp),
intent(in) :: array1(ia1)
1965 real(
rp),
intent(in) :: array2(ia2)
1967 real(
rp),
intent(in),
optional :: undef
1977 if( ia1 /= ia2 )
then
1978 write(*,*)
"Error: different array size !!"
1983 if( abs( array2(i) - undef ) > eps )
then
2000 integer,
intent(in) :: ia1, ia2
2001 real(
rp),
intent(in) :: array1(ia1)
2002 real(
rp),
intent(in) :: array2(ia2)
2004 real(
rp),
intent(in),
optional :: undef
2014 if( ia1 /= ia2 )
then
2015 write(*,*)
"Error: different array size !!"
2019 if(
present( undef ) )
then
2022 if( abs( array1(i) - undef ) > eps .and. &
2023 abs( array2(i) - undef ) > eps )
then