scale_rm_statistics.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
#include "inc_openmp.h"
module scale_rm_statistics
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use mpi
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_process, only: &
     prc_local_comm_world
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: stat_setup
  public :: stat_total
  public :: stat_detail

  interface stat_total
     module procedure stat_total_2d
     module procedure stat_total_3d
  end interface stat_total

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  logical, public :: statistics_checktotal = .false. 

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  logical, private :: statistics_use_globalcomm = .false. 

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine stat_setup
    use scale_process, only: &
       prc_mpistop
    implicit none

    namelist / param_statistics / &
       statistics_checktotal, &
       statistics_use_globalcomm

    integer :: ierr
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[STATISTICS] / Categ[ATMOS-RM COMM] / Origin[SCALElib]'

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_statistics,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_STATISTICS. Check!'
       call prc_mpistop
    endif
    if( io_lnml ) write(io_fid_log,nml=param_statistics)

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '*** Caluculate statistics?                     : ', statistics_checktotal
    if( io_l ) write(io_fid_log,*) '*** Allow global communication for statistics? : ', statistics_use_globalcomm
    if ( statistics_use_globalcomm ) then
       if( io_l ) write(io_fid_log,*) '*** Global total is calculated using MPI_ALLreduce.'
    else
       if( io_l ) write(io_fid_log,*) '*** Local total is calculated in each process.'
    endif

    return
  end subroutine stat_setup

  !-----------------------------------------------------------------------------
  subroutine stat_total_2d( allstatval, var, varname )
    use scale_process, only: &
       prc_myrank, &
       prc_mpistop
    use scale_comm, only: &
       comm_datatype
    use scale_grid_real, only: &
       area => real_area
    implicit none

    real(RP),         intent(out) :: allstatval
    real(RP),         intent(in)  :: var(ia,ja)
    character(len=*), intent(in)  :: varname

    character(len=H_SHORT) :: varname_trim
    real(RP) :: statval

    integer :: ierr
    integer :: i, j
    !---------------------------------------------------------------------------

    varname_trim = trim(varname)

    statval = 0.0_rp
    !$omp parallel do private(i,j) OMP_SCHEDULE_ collapse(2) reduction(+:statval)
    do j = js, je
    do i = is, ie
       statval = statval + var(i,j) * area(i,j)
    enddo
    enddo

    if ( .NOT. ( statval > -1.0_rp .OR. statval < 1.0_rp ) ) then ! must be NaN
       write(*,*) 'xxx [STAT_total] NaN is detected for ', varname_trim, ' in rank ', prc_myrank
       call prc_mpistop
    endif

    if ( statistics_use_globalcomm ) then
       call prof_rapstart('COMM_Allreduce', 2)
       ! All reduce
       call mpi_allreduce( statval,              &
                           allstatval,           &
                           1,                    &
                           comm_datatype,        &
                           mpi_sum,              &
                           prc_local_comm_world, &
                           ierr                  )

       call prof_rapend  ('COMM_Allreduce', 2)

       ! statistics over the all node
       if ( varname_trim /= "" ) then ! if varname is empty, suppress output
          if( io_l ) write(io_fid_log,'(1x,A,A,A,ES24.17)') &
                     '[', varname_trim, '] SUM(global) =', allstatval
       endif
    else
       allstatval = statval

       ! statistics on each node
       if ( varname_trim /= "" ) then ! if varname is empty, suppress output
          if( io_l ) write(io_fid_log,'(1x,A,A,A,ES24.17)') &
                     '[', varname_trim, '] SUM(local)  =', statval
       endif
    endif

    return
  end subroutine stat_total_2d

  !-----------------------------------------------------------------------------
  subroutine stat_total_3d( allstatval, var, varname )
    use scale_process, only: &
       prc_myrank, &
       prc_mpistop
    use scale_comm, only: &
       comm_datatype
    use scale_grid_real, only: &
       vol  => real_vol
    implicit none

    real(RP),         intent(out) :: allstatval
    real(RP),         intent(in)  :: var(ka,ia,ja)
    character(len=*), intent(in)  :: varname

    character(len=H_SHORT) :: varname_trim
    real(RP) :: statval

    integer :: ierr
    integer :: k, i, j
    !---------------------------------------------------------------------------

    varname_trim = trim(varname)

    statval = 0.0_rp
    !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2) reduction(+:statval)
    do j = js, je
    do i = is, ie
    do k = ks, ke
       statval = statval + var(k,i,j) * vol(k,i,j)
    enddo
    enddo
    enddo

    if ( .NOT. ( statval > -1.0_rp .OR. statval < 1.0_rp ) ) then ! must be NaN
       write(*,*) 'xxx [STAT_total] NaN is detected for ', varname_trim, ' in rank ', prc_myrank
       call prc_mpistop
    endif

    if ( statistics_use_globalcomm ) then
       call prof_rapstart('COMM_Allreduce', 2)
       ! All reduce
       call mpi_allreduce( statval,              &
                           allstatval,           &
                           1,                    &
                           comm_datatype,        &
                           mpi_sum,              &
                           prc_local_comm_world, &
                           ierr                  )

       call prof_rapend  ('COMM_Allreduce', 2)

       ! statistics over the all node
       if ( varname_trim /= "" ) then ! if varname is empty, suppress output
          if( io_l ) write(io_fid_log,'(1x,A,A,A,ES24.17)') &
                     '[', varname_trim, '] SUM(global) =', allstatval
       endif
    else
       allstatval = statval

       ! statistics on each node
       if ( varname_trim /= "" ) then ! if varname is empty, suppress output
          if( io_l ) write(io_fid_log,'(1x,A,A,A,ES24.17)') &
                     '[', varname_trim, '] SUM(local)  =', statval
       endif
    endif

    return
  end subroutine stat_total_3d

  !-----------------------------------------------------------------------------
  subroutine stat_detail(var, varname, supress_globalcomm)
    use scale_process, only: &
       prc_nprocs,   &
       prc_myrank
    use scale_const, only: &
       const_undef8, &
       const_undef2
    use scale_comm, only: &
       comm_datatype
    implicit none

    real(RP),         intent(inout) :: var(:,:,:,:)
    character(len=*), intent(in)    :: varname(:)
    logical,          intent(in), optional :: supress_globalcomm

    logical , allocatable :: halomask  (:,:,:)
    real(RP), allocatable :: statval   (:,:,:)
    integer,  allocatable :: statidx   (:,:,:,:)
    real(RP), allocatable :: allstatval(:,:)
    integer,  allocatable :: allstatidx(:,:,:)
    integer               :: ksize, vsize
    logical               :: do_globalcomm

    integer :: ierr
    integer :: v, p
    !---------------------------------------------------------------------------

    do_globalcomm = statistics_use_globalcomm
    if ( present(supress_globalcomm) ) then
       if ( supress_globalcomm ) then
          do_globalcomm = .false.
       endif
    endif

    ksize = size(var(:,:,:,:),1)
    vsize = size(var(:,:,:,:),4)

    allocate( halomask(ksize,ia,ja) ); halomask(:,:,:) = .false.

    if ( ksize == ka ) then
       halomask(ks:ke,is:ie,js:je) = .true.
    else
       halomask(:,is:ie,js:je) = .true.
    endif

    allocate( statval(  vsize,2,0:prc_nprocs-1) ); statval(:,:,:)   = const_undef8
    allocate( statidx(3,vsize,2,0:prc_nprocs-1) ); statidx(:,:,:,:) = const_undef2

    allocate( allstatval(  vsize,2) ); allstatval(:,:)   = const_undef8
    allocate( allstatidx(1,vsize,2) ); allstatidx(:,:,:) = const_undef2

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '*** Variable Statistics ***'
    do v = 1, vsize
       statval(  v,1,prc_myrank) = maxval(var(:,:,:,v),mask=halomask)
       statval(  v,2,prc_myrank) = minval(var(:,:,:,v),mask=halomask)
       statidx(:,v,1,prc_myrank) = maxloc(var(:,:,:,v),mask=halomask)
       statidx(:,v,2,prc_myrank) = minloc(var(:,:,:,v),mask=halomask)
    enddo

    if ( do_globalcomm ) then
       call prof_rapstart('COMM_Bcast', 2)
       do p = 0, prc_nprocs-1

          call mpi_bcast( statval(1,1,p),       &
                          vsize*2,              &
                          comm_datatype,        &
                          p,                    &
                          prc_local_comm_world, &
                          ierr                  )

          call mpi_bcast( statidx(1,1,1,p),     &
                          3*vsize*2,            &
                          mpi_integer,          &
                          p,                    &
                          prc_local_comm_world, &
                          ierr                  )

       enddo
       call prof_rapend  ('COMM_Bcast', 2)

       do v = 1, vsize
          allstatval(v,1)   = maxval(statval(v,1,:))
          allstatval(v,2)   = minval(statval(v,2,:))
          allstatidx(:,v,1) = maxloc(statval(v,1,:))-1
          allstatidx(:,v,2) = minloc(statval(v,2,:))-1
          if( io_l ) write(io_fid_log,*) '[', trim(varname(v)), ']'
          if( io_l ) write(io_fid_log,'(1x,A,ES17.10,A,4(I5,A))') '  MAX =', &
                                                       allstatval(  v,1), '(', &
                                                       allstatidx(1,v,1), ',', &
                                         statidx(1,v,1,allstatidx(1,v,1)),',', &
                                         statidx(2,v,1,allstatidx(1,v,1)),',', &
                                         statidx(3,v,1,allstatidx(1,v,1)),')'
          if( io_l ) write(io_fid_log,'(1x,A,ES17.10,A,4(I5,A))') '  MIN =', &
                                                       allstatval(  v,2), '(', &
                                                       allstatidx(1,v,2), ',', &
                                         statidx(1,v,2,allstatidx(1,v,2)),',', &
                                         statidx(2,v,2,allstatidx(1,v,2)),',', &
                                         statidx(3,v,2,allstatidx(1,v,2)),')'
       enddo
    else
       ! statistics on each node
       do v = 1, vsize
          if( io_l ) write(io_fid_log,*) '*** [', trim(varname(v)), ']'
          if( io_l ) write(io_fid_log,'(1x,A,ES17.10,A,3(I5,A))') '*** MAX = ', &
                                                statval(  v,1,prc_myrank),'(', &
                                                statidx(1,v,1,prc_myrank),',', &
                                                statidx(2,v,1,prc_myrank),',', &
                                                statidx(3,v,1,prc_myrank),')'
          if( io_l ) write(io_fid_log,'(1x,A,ES17.10,A,3(I5,A))') '*** MIN = ', &
                                                statval(  v,2,prc_myrank),'(', &
                                                statidx(1,v,2,prc_myrank),',', &
                                                statidx(2,v,2,prc_myrank),',', &
                                                statidx(3,v,2,prc_myrank),')'
       enddo
    endif

    if( io_l ) write(io_fid_log,*)

    deallocate( halomask )

    deallocate( statval )
    deallocate( statidx )

    deallocate( allstatval )
    deallocate( allstatidx )

    return
  end subroutine stat_detail

end module scale_rm_statistics