!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS 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 FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module pmaxmin_mod

  implicit none
  private
  public :: pmaxming, pmaxmin, prt_maxmin_local

contains
! Parallelized utility routines for computing/printing
! max/min of an input array
!

  subroutine pmaxmin( qname, a, pmin, pmax, im, jt, fac )

#if defined( SPMD )
#define CPP_PRT_PREFIX  if(gid.eq.0)
    use mod_comm, only : gid, mp_reduce_max
#else
#define CPP_PRT_PREFIX
#endif
    use fv_arrays_mod, only: fv_array_check, fv_array_limits, fv_stack_push, fv_array_sync
    use mpp_mod, only: mpp_sync

    implicit none

    character*(*)  qname
    integer im, jt
    integer i, j
    real a(im,jt)

    real qmin(jt), qmax(jt)
    real pmax, pmin
    real fac                     ! multiplication factor
    real pm1(2)

    integer :: js, je
#ifdef use_shared_pointers
    pointer( p_qmin, qmin )
    pointer( p_qmax, qmax )
    call fv_stack_push( p_qmin, jt )
    call fv_stack_push( p_qmax, jt )
#endif

    call fv_array_check( LOC(a) )
    call fv_array_limits( 1, jt, js, je )
!$omp parallel do private(i, j, pmax, pmin)

    do j=js,je
       pmax = a(1,j)
       pmin = a(1,j)
       do i=2,im
          pmax = max(pmax, a(i,j))
          pmin = min(pmin, a(i,j))
       enddo
       qmax(j) = pmax
       qmin(j) = pmin
    enddo
    call fv_array_sync()
!
! Now find max/min of amax/amin
!
    pmax = qmax(1)
    pmin = qmin(1)
    do j=2,jt
       pmax = max(pmax, qmax(j))
       pmin = min(pmin, qmin(j))
    enddo


#if defined( SPMD )
    pm1(1) = pmax
    call mp_reduce_max(1, pm1)
    pmax=pm1(1)
    pm1(1) = -pmin
    call mp_reduce_max(1, pm1)
    pmin=-pm1(1)
#endif

    CPP_PRT_PREFIX write(6,*) qname, ' max = ', pmax*fac, ' min = ', pmin*fac

  end subroutine pmaxmin

  subroutine pmaxming(qname, a, im, jm, km,     &
       jfirst, jlast, ng_d, ng_s, fac)

    use fv_arrays_mod, only: fv_array_check, fv_array_sync, jsp, jep
    use mpp_mod, only: mpp_sync
#ifdef use_shared_pointers
   use fv_arrays_mod, only: fv_stack_push
#endif
    implicit none

    character*(*)  qname
    integer im, jm, km, jfirst, jlast, ng_d, ng_s
    real a(im,jfirst-ng_d:jlast+ng_s,km)
    real tmp(im, jfirst:jlast,km)
    real pmax, pmin, fac
    integer i, j, k
#ifdef use_shared_pointers
    pointer( p_tmp, tmp )
    call fv_stack_push( p_tmp, im*km*(jlast-jfirst+1) )
#endif
    call fv_array_check( LOC(a) )

!$omp parallel do private(i, j, k)
!balaji: note, parallelizing j not k... km may be 1.
    do k=1, km
       do j=jsp,jep
          do i=1, im
             tmp(i,j,k)=a(i,j,k)
          enddo
       enddo
    enddo
    call fv_array_sync()

    call pmaxmin(qname, tmp, pmin, pmax, im*(jlast-jfirst+1),km, fac)

  end subroutine pmaxming


! Routine for debug domain decomposed data

  subroutine prt_maxmin_local(gid, qname, q, im, jm, km, j1, j2)

    implicit none
    character*(*)  qname
    integer im, jm, km, j1, j2, gid
    real q(im,j1:j2,km)
    integer i,j,k,jp
    real qmax, qmin

    jp = max(1,j1)

    qmax = q(1,jp,1)
    qmin = qmax

    do k=1,km
       do j=max(2,j1), min(jm-1,j2)
          do i=1,im
             qmax = max(qmax, q(i,j,k))
             qmin = min(qmin, q(i,j,k))
          enddo
       enddo
    enddo
    write(*,*) 'GID=', gid, qname, ': Max=', qmax, ' Min=', qmin

  end subroutine prt_maxmin_local

end module pmaxmin_mod