source: palm/trunk/SOURCE/global_min_max.f90 @ 4892

!--------------------------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
! Public License as published by the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
! Public License for more details.
!
! You should have received a copy of the GNU General Public License along with PALM. If not, see
! <http://www.gnu.org/licenses/>.
!
! Copyright 1997-2021 Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
!
!
! Former revisions:
! -----------------
! $Id: global_min_max.f90 4855 2021-01-25 12:30:54Z suehring $
! bugfix for correct identification of indices of extreme values in case of non-cyclic
! boundary conditions
!
! 4828 2021-01-05 11:21:41Z Giersch
! preprocessor branch for ibm removed
!
! 4646 2020-08-24 16:02:40Z raasch
! file re-formatted to follow the PALM coding standard
!
! 4429 2020-02-27 15:24:30Z raasch
! bugfix: cpp-directives added for serial mode
!
! 4360 2020-01-07 11:25:50Z suehring
! OpenACC support added
!
! 4182 2019-08-22 15:20:23Z scharf
! Corrected "Former revisions" section
!
! 3655 2019-01-07 16:51:22Z knoop
! Corrected "Former revisions" section
!
! Revision 1.1  1997/07/24 11:14:03  raasch
! Initial revision
!
!
! Description:
! ------------
!> Determine the array minimum/maximum and the corresponding indices.
!--------------------------------------------------------------------------------------------------!
 SUBROUTINE global_min_max( i1, i2, j1, j2, k1, k2, ar, mode, offset, value, value_ijk, value1,    &
                            value1_ijk )


    USE control_parameters,                                                                        &
        ONLY:  bc_lr, bc_ns

    USE indices,                                                                                   &
        ONLY:  nbgp, ny, nx

    USE kinds

    USE pegrid

    IMPLICIT NONE

    CHARACTER (LEN=*) ::  mode  !<

    INTEGER(iwp) ::  i              !<
    INTEGER(iwp) ::  i1             !<
    INTEGER(iwp) ::  i2             !<
#if defined( __parallel )
    INTEGER(iwp) ::  id_fmax        !<
    INTEGER(iwp) ::  id_fmin        !<
#endif
    INTEGER(iwp) ::  j              !<
    INTEGER(iwp) ::  j1             !<
    INTEGER(iwp) ::  j2             !<
    INTEGER(iwp) ::  k              !<
    INTEGER(iwp) ::  k1             !<
    INTEGER(iwp) ::  k2             !<
    INTEGER(iwp) ::  value_ijk(3)   !<

    INTEGER(iwp), DIMENSION(3) ::  fmax_ijk    !<
    INTEGER(iwp), DIMENSION(3) ::  fmax_ijk_l  !<
    INTEGER(iwp), DIMENSION(3) ::  fmin_ijk    !<
    INTEGER(iwp), DIMENSION(3) ::  fmin_ijk_l  !<

    INTEGER(iwp), DIMENSION(3), OPTIONAL ::  value1_ijk  !<

    REAL(wp) ::  offset            !<
    REAL(wp) ::  value             !<
    REAL(wp), OPTIONAL ::  value1  !<

    REAL(wp), DIMENSION(2) ::  fmax    !<
    REAL(wp), DIMENSION(2) ::  fmax_l  !<
    REAL(wp), DIMENSION(2) ::  fmin    !<
    REAL(wp), DIMENSION(2) ::  fmin_l  !<

    REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) ::  ar  !<

#if defined( _OPENACC )
    INTEGER(iwp) ::  count_eq   !< counter for locations of maximum
    REAL(wp)     ::  red        !< scalar for reduction with OpenACC
#endif


!
!-- Determine array minimum
    IF ( mode == 'min'  .OR.  mode == 'minmax' )  THEN

!
!--    Determine the local minimum
       fmin_ijk_l = MINLOC( ar )
       fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - 1 ! MINLOC assumes lowerbound = 1
       fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - nbgp
       fmin_ijk_l(3) = k1 + fmin_ijk_l(3) - nbgp
       fmin_l(1)  = ar(fmin_ijk_l(1),fmin_ijk_l(2),fmin_ijk_l(3))

#if defined( __parallel )
       fmin_l(2)  = myid
       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
       CALL MPI_ALLREDUCE( fmin_l, fmin, 1, MPI_2REAL, MPI_MINLOC, comm2d, ierr )

!
!--    Determine the global minimum. Result stored on PE0.
       id_fmin = fmin(2)
       IF ( id_fmin /= 0 )  THEN
          IF ( myid == 0 )  THEN
             CALL MPI_RECV( fmin_ijk, 3, MPI_INTEGER, id_fmin, 0, comm2d, status, ierr )
          ELSEIF ( myid == id_fmin )  THEN
             CALL MPI_SEND( fmin_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
          ENDIF
       ELSE
          fmin_ijk = fmin_ijk_l
       ENDIF
!
!--    Send the indices of the just determined array minimum to other PEs
       CALL MPI_BCAST( fmin_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
#else
       fmin(1)  = fmin_l(1)
       fmin_ijk = fmin_ijk_l
#endif

    ENDIF

!
!-- Determine array maximum
    IF ( mode == 'max'  .OR.  mode == 'minmax' )  THEN

!
!--    Determine the local maximum
       fmax_ijk_l = MAXLOC( ar )
       fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - 1 ! MAXLOC assumes lowerbound = 1
       fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - nbgp
       fmax_ijk_l(3) = k1 + fmax_ijk_l(3) - nbgp
       fmax_l(1) = ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3))

#if defined( __parallel )
       fmax_l(2)  = myid
       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, ierr )

!
!--    Determine the global maximum. Result stored on PE0.
       id_fmax = fmax(2)
       IF ( id_fmax /= 0 )  THEN
          IF ( myid == 0 )  THEN
             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, status, ierr )
          ELSEIF ( myid == id_fmax )  THEN
             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
          ENDIF
       ELSE
          fmax_ijk = fmax_ijk_l
       ENDIF
!
!--    Send the indices of the just determined array maximum to other PEs
       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
#else
       fmax(1)  = fmax_l(1)
       fmax_ijk = fmax_ijk_l
#endif

    ENDIF

!
!-- Determine absolute array maximum
    IF ( mode == 'abs' )  THEN

#if defined( _OPENACC )
       red = 0.0_wp
       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
       !$ACC PRESENT(ar) COPY(red) REDUCTION(MAX: red)
       DO  k = k1, k2
          DO  j = j1, j2
             DO  i = i1, i2
                IF ( ABS( ar(i,j,k) ) > red )  THEN
                   red = ABS( ar(i,j,k) )
                ENDIF
             ENDDO
          ENDDO
       ENDDO
       fmax_l(1) = red

!
!--    Determine the maximum's position and count how often it is found.
       count_eq = 0
       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
       !$ACC PRESENT(ar) COPY(fmax_ijk_l(1:3), count_eq) &
       !$ACC REDUCTION(+:count_eq)
       DO  k = k1, k2
          DO  j = j1, j2
             DO  i = i1, i2
                IF ( ABS( ar(i,j,k) ) == red )  THEN
                   fmax_ijk_l(1) = i
                   fmax_ijk_l(2) = j
                   fmax_ijk_l(3) = k
                   count_eq = count_eq + 1
                ENDIF
             ENDDO
          ENDDO
       ENDDO

       IF ( count_eq == 1 ) THEN
!
!--       We found a single maximum element and correctly got its position. Transfer its value to
!--       handle the negative case correctly.
          !$ACC UPDATE HOST(ar(fmax_ijk_l(1):fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)))
       ELSE
!
!--       We found no maximum element (?) or multiple, so the position is not correct.
!--       Copy the full array to the host and determine the maximum sequentially...
          !$ACC UPDATE HOST(ar(i1:i2,j1:j2,k1:k2))
#endif

!
!--    Determine the local absolut maximum
       fmax_l(1)     = 0.0_wp
       fmax_ijk_l(1) =  i1
       fmax_ijk_l(2) =  j1
       fmax_ijk_l(3) =  k1
       DO  k = k1, k2
          DO  j = j1, j2
             DO  i = i1, i2
                IF ( ABS( ar(i,j,k) ) > fmax_l(1) )  THEN
                   fmax_l(1) = ABS( ar(i,j,k) )
                   fmax_ijk_l(1) = i
                   fmax_ijk_l(2) = j
                   fmax_ijk_l(3) = k
                ENDIF
             ENDDO
          ENDDO
       ENDDO

#if defined( _OPENACC )
!
!--    Close ELSE case from above
       ENDIF
#endif

!
!--    Set a flag in case that the determined value is negative.
!--    A constant offset has to be subtracted in order to handle the special case i=0 correctly.
       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0_wp )  THEN
          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
       ENDIF

#if defined( __parallel )
       fmax_l(2)  = myid
       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, ierr )

!
!--    Determine the global absolut maximum. Result stored on PE0.
       id_fmax = fmax(2)
       IF ( id_fmax /= 0 )  THEN
          IF ( myid == 0 )  THEN
             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, status, ierr )
          ELSEIF ( myid == id_fmax )  THEN
             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
          ENDIF
       ELSE
          fmax_ijk = fmax_ijk_l
       ENDIF
!
!--    Send the indices of the just determined absolut maximum to other PEs
       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
#else
       fmax(1)  = fmax_l(1)
       fmax_ijk = fmax_ijk_l
#endif

    ENDIF

!
!-- Determine absolute maximum of ( array - offset )
    IF ( mode == 'absoff' )  THEN

!
!--    Determine the local absolut maximum
       fmax_l(1)     = 0.0_wp
       fmax_ijk_l(1) =  i1
       fmax_ijk_l(2) =  j1
       fmax_ijk_l(3) =  k1
       DO  k = k1, k2
          DO  j = j1, j2
!
!--          Attention: the lowest gridpoint is excluded here, because there is no advection at
!--          ---------- nzb=0 and mode 'absoff' is only used for calculating u,v extrema for
!--                     CFL-criteria.
             DO  i = i1+1, i2
                IF ( ABS( ar(i,j,k) - offset ) > fmax_l(1) )  THEN
                   fmax_l(1) = ABS( ar(i,j,k) - offset )
                   fmax_ijk_l(1) = i
                   fmax_ijk_l(2) = j
                   fmax_ijk_l(3) = k
                ENDIF
             ENDDO
          ENDDO
       ENDDO

!
!--    Set a flag in case that the determined value is negative.
!--    A constant offset has to be subtracted in order to handle the special case i=0 correctly.
       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0_wp )  THEN
          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
       ENDIF

#if defined( __parallel )
       fmax_l(2)  = myid
       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, ierr )

!
!--    Determine the global absolut maximum. Result stored on PE0.
       id_fmax = fmax(2)
       IF ( id_fmax /= 0 )  THEN
          IF ( myid == 0 )  THEN
             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, status, ierr )
          ELSEIF ( myid == id_fmax )  THEN
             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
          ENDIF
       ELSE
          fmax_ijk = fmax_ijk_l
       ENDIF
!
!--    Send the indices of the just determined absolut maximum to other PEs
       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
#else
       fmax(1)  = fmax_l(1)
       fmax_ijk = fmax_ijk_l
#endif

    ENDIF

!
!-- Determine output parameters
    SELECT CASE( mode )

       CASE( 'min' )

          value     = fmin(1)
          value_ijk = fmin_ijk

       CASE( 'max' )

          value     = fmax(1)
          value_ijk = fmax_ijk

       CASE( 'minmax' )

          value      = fmin(1)
          value_ijk  = fmin_ijk
          value1     = fmax(1)
          value1_ijk = fmax_ijk

       CASE( 'abs', 'absoff' )

          value     = fmax(1)
          value_ijk = fmax_ijk
          IF ( fmax_ijk(1) <= -10 )  THEN
!
!--          Index needs to be corrected because it has been modified above to indicate negative
!--          values
             value_ijk(1) = -value_ijk(1) - 10
!
!--          For this reason also change the sign of the quantity
             value        = -value
          ENDIF

    END SELECT

!
!-- Limit index values to the range 0..nx, 0..ny. Non-cyclic setups may have extrema at the
!-- outer borders, which should be correctly identified.
    IF ( bc_lr == 'cyclic' )  THEN
       IF ( value_ijk(3) < 0  ) value_ijk(3) = nx +1 + value_ijk(3)
       IF ( value_ijk(3) > nx ) value_ijk(3) = value_ijk(3) - (nx+1)
    ENDIF
    IF ( bc_ns == 'cyclic' )  THEN
       IF ( value_ijk(2) < 0  ) value_ijk(2) = ny +1 + value_ijk(2)
       IF ( value_ijk(2) > ny ) value_ijk(2) = value_ijk(2) - (ny+1)
    ENDIF

 END SUBROUTINE global_min_max