mpp_util.inc

! -*-f90-*-
 
 
!***********************************************************************
!*                             Apache License 2.0
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* Licensed under the Apache License, Version 2.0 (the "License");
!* you may not use this file except in compliance with the License.
!* You may obtain a copy of the License at
!*
!*     http://www.apache.org/licenses/LICENSE-2.0
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
!* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
!* PARTICULAR PURPOSE. See the License for the specific language
!* governing permissions and limitations under the License.
!***********************************************************************
!> @file
!> @brief General utility functions for use in @ref mpp_mod
 
!> @addtogroup mpp_mod
!> @{
 
#if defined(use_libMPI)
#include <mpp_util_mpi.inc>
#else
#include <mpp_util_nocomm.inc>
#endif
 
  !> @brief This function returns the current standard fortran unit numbers for input.
  function stdin()
    integer :: stdin
    stdin = in_unit
    return
  end function stdin
 
  !> @brief This function returns the current  standard fortran unit numbers for output.
  function stdout()
    integer :: stdout
    stdout = out_unit
    if( pe.NE.root_pe )stdout = stdlog()
    return
  end function stdout
 
  !> @brief This function returns the current standard fortran unit numbers for error messages.
  function stderr()
    integer :: stderr
    stderr = err_unit
    return
  end function stderr
 
  !> @brief This function returns the current  standard fortran unit numbers for log messages.
  !!    Log messages, by convention, are written to the file <TT>logfile.out</TT>.
  function stdlog()
    integer :: stdlog
    logical :: opened
    character(len=11) :: this_pe
!$  logical           :: omp_in_parallel
!$  integer           :: omp_get_num_threads
!$  integer           :: errunit
 
 
!NOTES: We can not use mpp_error to handle the error because mpp_error
!       will call stdout and stdout will call stdlog for non-root-pe.
!       This will be a cicular call.
 
!$  if( omp_in_parallel() .and. (omp_get_num_threads() > 1) ) then
!$OMP single
!$      errunit = stderr()
!$      write( errunit,'(/a/)' ) 'FATAL: STDLOG: is called inside a OMP parallel region'
#ifdef use_libMPI
!$      call MPI_ABORT( MPI_COMM_WORLD, 1, error )
#else
!$      call ABORT()
#endif
!$OMP end single
!$  endif
 
    if( pe.EQ.root_pe )then
       write(this_pe,'(a,i6.6,a)') '.',pe,'.out'
       inquire( file=trim(configfile)//this_pe, opened=opened )
       if( opened )then
          FLUSH(log_unit)
       else
          open(newunit=log_unit, status='UNKNOWN', file=trim(configfile)//this_pe, position='APPEND', err=10 )
       end if
       stdlog = log_unit
    else
       inquire(unit=etc_unit, opened=opened )
       if( opened )then
          FLUSH(etc_unit)
       else
          open(newunit=etc_unit, status='UNKNOWN', file=trim(etcfile), position='APPEND', err=11 )
       end if
       stdlog = etc_unit
    end if
    return
10  call mpp_error( fatal, 'STDLOG: unable to open '//trim(configfile)//this_pe//'.' )
11  call mpp_error( fatal, 'STDLOG: unable to open '//trim(etcfile)//'.' )
  end function stdlog
 
  !#####################################################################
  subroutine mpp_init_logfile()
  integer :: p
  logical :: exist
  character(len=11) :: this_pe
  if( pe.EQ.root_pe )then
     do p=0,npes-1
       write(this_pe,'(a,i6.6,a)') '.',p,'.out'
       inquire( file=trim(configfile)//this_pe, exist=exist )
       if(exist)then
         open(newunit=log_unit, file=trim(configfile)//this_pe, status='REPLACE' )
         close(log_unit)
       endif
     end do
  end if
  end subroutine mpp_init_logfile
 
  !> Opens the warning log file, called during mpp_init
  subroutine mpp_init_warninglog()
  logical :: exist
  character(len=11) :: this_pe
  if( pe.EQ.root_pe )then
    write(this_pe,'(a,i6.6,a)') '.',pe,'.out'
    inquire( file=trim(warnfile)//this_pe, exist=exist )
    if(exist)then
      open(newunit=warn_unit, file=trim(warnfile)//this_pe, status='REPLACE' )
    else
      open(newunit=warn_unit, file=trim(warnfile)//this_pe, status='NEW' )
    endif
  end if
  end subroutine mpp_init_warninglog
 
  !> @brief This function returns unit number for the warning log
  !! if on the root pe, otherwise returns the etc_unit value (usually /dev/null)
  function warnlog()
    integer :: warnlog
    if(.not. module_is_initialized) call mpp_error(fatal, "mpp_mod: warnlog cannot be called before mpp_init")
    if(root_pe .eq. pe) then
      warnlog = warn_unit
    else
      warnlog = etc_unit
    endif
    return
  end function warnlog
 
  !#####################################################################
  subroutine mpp_set_warn_level(flag)
    integer, intent(in) :: flag
 
    if( flag.EQ.warning )then
       warnings_are_fatal = .false.
    else if( flag.EQ.fatal )then
       warnings_are_fatal = .true.
    else
       call mpp_error( fatal, 'MPP_SET_WARN_LEVEL: warning flag must be set to WARNING or FATAL.' )
    end if
    return
  end subroutine mpp_set_warn_level
 
  !#####################################################################
  function mpp_error_state()
    integer :: mpp_error_state
    mpp_error_state = error_state
    return
  end function mpp_error_state
 
!#####################################################################
!> @brief overloads to mpp_error_basic, support for error_mesg routine in FMS
subroutine mpp_error_mesg( routine, errormsg, errortype )
  character(len=*), intent(in) :: routine, errormsg
  integer,          intent(in) :: errortype
 
  call mpp_error( errortype, trim(routine)//': '//trim(errormsg) )
  return
end subroutine mpp_error_mesg
 
!#####################################################################
subroutine mpp_error_noargs()
  call mpp_error(fatal)
end subroutine mpp_error_noargs
 
!#####################################################################
subroutine mpp_error_is(errortype, errormsg1, mpp_ival, errormsg2)
  integer,          intent(in) :: errortype
  INTEGER,          intent(in) :: mpp_ival
  character(len=*), intent(in) :: errormsg1
  character(len=*),      intent(in), optional :: errormsg2
  call mpp_error( errortype, errormsg1, (/mpp_ival/), errormsg2)
end subroutine mpp_error_is
!#####################################################################
subroutine mpp_error_rs(errortype, errormsg1, mpp_rval, errormsg2)
  integer,          intent(in) :: errortype
  REAL,             intent(in) :: mpp_rval
  character(len=*), intent(in) :: errormsg1
  character(len=*),      intent(in), optional :: errormsg2
  call mpp_error( errortype, errormsg1, (/mpp_rval/), errormsg2)
end subroutine mpp_error_rs
!#####################################################################
subroutine mpp_error_ia(errortype, errormsg1, array, errormsg2)
  integer,               intent(in) :: errortype
  INTEGER, dimension(:), intent(in) :: array
  character(len=*),      intent(in) :: errormsg1
  character(len=*),      intent(in), optional :: errormsg2
  character(len=512) :: string
 
  string = errormsg1//trim(array_to_char(array))
  if(present(errormsg2)) string = trim(string)//errormsg2
  call mpp_error_basic( errortype, trim(string))
 
end subroutine mpp_error_ia
 
!#####################################################################
subroutine mpp_error_ra(errortype, errormsg1, array, errormsg2)
  integer,            intent(in) :: errortype
  REAL, dimension(:), intent(in) :: array
  character(len=*),      intent(in) :: errormsg1
  character(len=*),   intent(in), optional :: errormsg2
  character(len=512) :: string
 
  string = errormsg1//trim(array_to_char(array))
  if(present(errormsg2)) string = trim(string)//errormsg2
  call mpp_error_basic( errortype, trim(string))
 
end subroutine mpp_error_ra
 
!#####################################################################
#define _SUBNAME_ mpp_error_ia_ia
#define _ARRAY1TYPE_ integer
#define _ARRAY2TYPE_ integer
#include <mpp_error_a_a.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_ia_ra
#define _ARRAY1TYPE_ integer
#define _ARRAY2TYPE_ real
#include <mpp_error_a_a.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_ra_ia
#define _ARRAY1TYPE_ real
#define _ARRAY2TYPE_ integer
#include <mpp_error_a_a.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_ra_ra
#define _ARRAY1TYPE_ real
#define _ARRAY2TYPE_ real
#include <mpp_error_a_a.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_ia_is
#define _ARRAY1TYPE_ integer
#define _ARRAY2TYPE_ integer
#include <mpp_error_a_s.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_ia_rs
#define _ARRAY1TYPE_ integer
#define _ARRAY2TYPE_ real
#include <mpp_error_a_s.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_ra_is
#define _ARRAY1TYPE_ real
#define _ARRAY2TYPE_ integer
#include <mpp_error_a_s.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_ra_rs
#define _ARRAY1TYPE_ real
#define _ARRAY2TYPE_ real
#include <mpp_error_a_s.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_is_ia
#define _ARRAY1TYPE_ integer
#define _ARRAY2TYPE_ integer
#include <mpp_error_s_a.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_is_ra
#define _ARRAY1TYPE_ integer
#define _ARRAY2TYPE_ real
#include <mpp_error_s_a.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_rs_ia
#define _ARRAY1TYPE_ real
#define _ARRAY2TYPE_ integer
#include <mpp_error_s_a.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_rs_ra
#define _ARRAY1TYPE_ real
#define _ARRAY2TYPE_ real
#include <mpp_error_s_a.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_is_is
#define _ARRAY1TYPE_ integer
#define _ARRAY2TYPE_ integer
#include <mpp_error_s_s.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_is_rs
#define _ARRAY1TYPE_ integer
#define _ARRAY2TYPE_ real
#include <mpp_error_s_s.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_rs_is
#define _ARRAY1TYPE_ real
#define _ARRAY2TYPE_ integer
#include <mpp_error_s_s.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
#define _SUBNAME_ mpp_error_rs_rs
#define _ARRAY1TYPE_ real
#define _ARRAY2TYPE_ real
#include <mpp_error_s_s.fh>
#undef _SUBNAME_
#undef _ARRAY1TYPE_
#undef _ARRAY2TYPE_
!#####################################################################
function iarray_to_char(iarray) result(string)
integer, intent(in) :: iarray(:)
character(len=256) :: string
character(len=32)  :: chtmp
integer :: i, len_tmp, len_string
 
 string = ''
 do i=1,size(iarray)
   write(chtmp,'(i16)') iarray(i)
   chtmp = adjustl(chtmp)
   len_tmp = len_trim(chtmp)
   len_string  = len_trim(string)
   string(len_string+1:len_string+len_tmp) = trim(chtmp)
   string(len_string+len_tmp+1:len_string+len_tmp+1) = ','
 enddo
 len_string = len_trim(string)
 string(len_string:len_string) = ' ' ! remove trailing comma
 
end function iarray_to_char
!#####################################################################
function rarray_to_char(rarray) result(string)
real, intent(in) :: rarray(:)
character(len=256) :: string
character(len=32)  :: chtmp
integer :: i, len_tmp, len_string
 
 string = ''
 do i=1,size(rarray)
   write(chtmp,'(G16.9)') rarray(i)
   chtmp = adjustl(chtmp)
   len_tmp = len_trim(chtmp)
   len_string  = len_trim(string)
   string(len_string+1:len_string+len_tmp) = trim(chtmp)
   string(len_string+len_tmp+1:len_string+len_tmp+1) = ','
 enddo
 len_string = len_trim(string)
 string(len_string:len_string) = ' ' ! remove trailing comma
 
end function rarray_to_char
 
  !> @brief Returns processor ID.
  !!
  !> This returns the unique ID associated with a PE. This number runs
  !! between 0 and <TT>npes-1</TT>, where <TT>npes</TT> is the total
  !! processor count, returned by <TT>mpp_npes</TT>. For a uniprocessor
  !! application this will always return 0.
  function mpp_pe()
    integer :: mpp_pe
 
    if( .NOT.module_is_initialized )call mpp_error( fatal, 'MPP_PE: You must first call mpp_init.' )
    mpp_pe = pe
    return
  end function mpp_pe
 
  !#####################################################################
 
  !> @brief Returns processor count for current pelist
  !!
  !> This returns the number of PEs in the current pelist. For a uniprocessor application,
  !! it will always return 1.
  function mpp_npes()
    integer :: mpp_npes
 
    if( .NOT.module_is_initialized )call mpp_error( fatal, 'MPP_NPES: You must first call mpp_init.' )
    mpp_npes = size(peset(current_peset_num)%list(:))
    return
  end function mpp_npes
 
  !#####################################################################
  function mpp_root_pe()
    integer :: mpp_root_pe
 
    if( .NOT.module_is_initialized )call mpp_error( fatal, 'MPP_ROOT_PE: You must first call mpp_init.' )
    mpp_root_pe = root_pe
    return
  end function mpp_root_pe
 
  function mpp_comm()
    type(mpi_comm) :: mpp_comm
 
    if( .NOT.module_is_initialized )call mpp_error( fatal, 'mpp_comm: You must first call mpp_init.' )
    mpp_comm = peset(current_peset_num)%comm
  end function mpp_comm
 
  function mpp_commid()
    integer :: mpp_commid
 
    if( .NOT.module_is_initialized )call mpp_error( fatal, 'mpp_commID: You must first call mpp_init.' )
    call mpp_error(note, "mpp_commID() is deprecated. Please use mpp_comm() instead.")
    mpp_commid = peset(current_peset_num)%comm%mpi_val
  end function mpp_commid
 
  !#####################################################################
  subroutine mpp_set_root_pe(num)
    integer, intent(in) :: num
 
    if( .NOT.module_is_initialized )call mpp_error( fatal, 'MPP_SET_ROOT_PE: You must first call mpp_init.' )
    if( .NOT.(any(num.EQ.peset(current_peset_num)%list(:))) ) &
         call mpp_error( fatal, 'MPP_SET_ROOT_PE: you cannot set a root PE outside the current pelist.' )
    root_pe = num
    return
  end subroutine mpp_set_root_pe
 
  !> @brief Declare a pelist.
  !!
  !> This call is written specifically to accommodate a MPI restriction
  !! that requires a parent communicator to create a child communicator, In
  !! other words: a pelist cannot go off and declare a communicator, but
  !! every PE in the parent, including those not in pelist(:), must get
  !! together for the <TT>MPI_COMM_CREATE</TT> call. The parent is
  !! typically <TT>MPI_COMM_WORLD</TT>, though it could also be a subset
  !! that includes all PEs in <TT>pelist</TT>.
  !!
  !! This call implies synchronization across the PEs in the current
  !! pelist, of which <TT>pelist</TT> is a subset.
  subroutine mpp_declare_pelist_f08( pelist, name, comm )
    integer,                    intent(in) :: pelist(:) !> pelist you are declaring and storing within FMS
    character(len=*), intent(in), optional :: name      !> unique name for an input pelist
    type(mpi_comm), intent(out), optional  :: comm      !> mpi_comm communicator handle
    integer :: i
 
    if( .NOT.module_is_initialized )call mpp_error( fatal, 'MPP_DECLARE_PELIST: You must first call mpp_init.' )
    i = get_peset(pelist)
    write( peset(i)%name,'(a,i2.2)' ) 'PElist', i !default name
    if( PRESENT(name) ) peset(i)%name = name
    if( PRESENT(comm) ) then
      comm = peset(i)%comm
    endif
  end subroutine mpp_declare_pelist_f08
 
  subroutine mpp_declare_pelist_legacy( pelist, name, commID )
    integer,                    intent(in) :: pelist(:) !> pelist you are declaring and storing within FMS
    character(len=*), intent(in), optional :: name      !> unique name for an input pelist
    integer, intent(out)                   :: commID    !> integral MPI communicator handle
    type(mpi_comm) :: comm
 
    call mpp_declare_pelist_f08(pelist, name, comm)
    commid = comm%mpi_val
  end subroutine mpp_declare_pelist_legacy
 
  !#####################################################################
 
  !> @brief Set context pelist
  !!
  !! This call sets the value of the current pelist, which is the
  !! context for all subsequent "global" calls where the optional
  !! <TT>pelist</TT> argument is omitted. All the PEs that are to be in the
  !! current pelist must call it.
  !!
  !! In MPI, this call may hang unless <TT>pelist</TT> has been previous
  !! declared using @ref mpp_declare_pelist
  !!
  !! If the argument <TT>pelist</TT> is absent, the current pelist is
  !! set to the "world" pelist, of all PEs in the job.
  subroutine mpp_set_current_pelist( pelist, no_sync )
    !Once we branch off into a PE subset, we want subsequent "global" calls to
    !sync only across this subset. This is declared as the current pelist (peset(current_peset_num)%list)
    !when current_peset all pelist ops with no pelist should apply the current pelist.
    !also, we set the start PE in this pelist to be the root_pe.
    !unlike mpp_declare_pelist, this is called by the PEs in the pelist only
    !so if the PEset has not been previously declared, this will hang in MPI.
    !if pelist is omitted, we reset pelist to the world pelist.
    integer, intent(in), optional :: pelist(:)
    logical, intent(in), optional :: no_sync
 
    if( .NOT.module_is_initialized )call mpp_error( fatal, 'MPP_SET_CURRENT_PELIST: You must first call mpp_init.' )
    if( PRESENT(pelist) )then
       if( .NOT.any(pe.EQ.pelist) )call mpp_error( fatal, 'MPP_SET_CURRENT_PELIST: pe must be in pelist.' )
       current_peset_num = get_peset(pelist)
    else
       current_peset_num = world_peset_num
    end if
    call mpp_set_root_pe( minval(peset(current_peset_num)%list(:)) )
    if(.not.PRESENT(no_sync))call mpp_sync()  !this is called to make sure everyone in the current pelist is here.
    !      npes = mpp_npes()
    return
  end subroutine mpp_set_current_pelist
 
  !#####################################################################
  function mpp_get_current_pelist_name()
   ! Simply return the current pelist name
   character(len=len(peset(current_peset_num)%name)) :: mpp_get_current_pelist_name
 
   mpp_get_current_pelist_name = peset(current_peset_num)%name
  end function mpp_get_current_pelist_name
 
  !this is created for use by mpp_define_domains within a pelist
  !will be published but not publicized
  subroutine mpp_get_current_pelist_f08( pelist, name, comm )
    integer, intent(out)                    :: pelist(:) !> Array to copy the pelist into
    character(len=*), intent(out), optional :: name      !> Name of the pelist
    type(mpi_comm), intent(out), optional   :: comm      !> mpi_comm communicator handle
 
    if( size(pelist(:)).NE.size(peset(current_peset_num)%list(:)) ) &
         call mpp_error( fatal, 'MPP_GET_CURRENT_PELIST: size(pelist) is wrong.' )
    pelist(:) = peset(current_peset_num)%list(:)
    if( PRESENT(name) ) name = peset(current_peset_num)%name
    if( PRESENT(comm) ) then
      comm = peset(current_peset_num)%comm
    endif
  end subroutine mpp_get_current_pelist_f08
 
  subroutine mpp_get_current_pelist_legacy( pelist, name, commID )
    integer, intent(out)                    :: pelist(:) !> Array to copy the pelist into
    character(len=*), intent(out), optional :: name      !> Name of the pelist
    integer, intent(out)                    :: commID    !> Integral MPI communicator handle
    type(mpi_comm) :: comm
 
    call mpp_get_current_pelist_f08(pelist, name, comm)
    commid = comm%mpi_val
  end subroutine mpp_get_current_pelist_legacy
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                             !
  !                        PERFORMANCE PROFILING CALLS                          !
  !                                                                             !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!predefined clock granularities, but you can use any integer!using CLOCK_LOOP and above may distort coarser-grain measurements
  !> @brief Set the level of granularity of timing measurements.
  !!
  !> This routine and three other routines, mpp_clock_id, mpp_clock_begin(id),
  !!    and mpp_clock_end(id) may be used to time parallel code sections, and
  !!    extract parallel statistics. Clocks are identified by names, which
  !!    should be unique in the first 32 characters. The <TT>mpp_clock_id</TT>
  !!    call initializes a clock of a given name and returns an integer
  !!    <TT>id</TT>. This <TT>id</TT> can be used by subsequent
  !!    <TT>mpp_clock_begin</TT> and <TT>mpp_clock_end</TT> calls set around a
  !!    code section to be timed. Example:
  !!    <PRE>
  !!    integer :: id
  !!    id = mpp_clock_id( 'Atmosphere' )
  !!    call mpp_clock_begin(id)
  !!    call atmos_model()
  !!    call mpp_clock_end()
  !!    </PRE>
  !!     Two flags may be used to alter the behaviour of
  !!     <TT>mpp_clock</TT>. If the flag <TT>MPP_CLOCK_SYNC</TT> is turned on
  !!     by <TT>mpp_clock_id</TT>, the clock calls <TT>mpp_sync</TT> across all
  !!     the PEs in the current pelist at the top of the timed code section,
  !!     but allows each PE to complete the code section (and reach
  !!     <TT>mpp_clock_end</TT>) at different times. This allows us to measure
  !!     load imbalance for a given code section. Statistics are written to
  !!     <TT>stdout</TT> by <TT>mpp_exit</TT>.
  !!
  !!     The flag <TT>MPP_CLOCK_DETAILED</TT> may be turned on by
  !!     <TT>mpp_clock_id</TT> to get detailed communication
  !!     profiles. Communication events of the types <TT>SEND, RECV, BROADCAST,
  !!     REDUCE</TT> and <TT>WAIT</TT> are separately measured for data volume
  !!     and time. Statistics are written to <TT>stdout</TT> by
  !!     <TT>mpp_exit</TT>, and individual PE info is also written to the file
  !!     <TT>mpp_clock.out.####</TT> where <TT>####</TT> is the PE id given by
  !!     <TT>mpp_pe</TT>.
  !!
  !!     The flags <TT>MPP_CLOCK_SYNC</TT> and <TT>MPP_CLOCK_DETAILED</TT> are
  !!     integer parameters available by use association, and may be summed to
  !!     turn them both on.
  !!
  !!     While the nesting of clocks is allowed, please note that turning on
  !!     the non-optional flags on inner clocks has certain subtle issues.
  !!     Turning on <TT>MPP_CLOCK_SYNC</TT> on an inner
  !!     clock may distort outer clock measurements of load imbalance. Turning
  !!     on <TT>MPP_CLOCK_DETAILED</TT> will stop detailed measurements on its
  !!     outer clock, since only one detailed clock may be active at one time.
  !!     Also, detailed clocks only time a certain number of events per clock
  !!     (currently 40000) to conserve memory. If this array overflows, a
  !!     warning message is printed, and subsequent events for this clock are
  !!     not timed.
  !!
  !!     Timings are done using the <TT>f90</TT> standard
  !!     <TT>SYSTEM_CLOCK</TT> intrinsic.
  !!
  !!     The resolution of SYSTEM_CLOCK is often too coarse for use except
  !!     across large swaths of code. On SGI systems this is transparently
  !!     overloaded with a higher resolution clock made available in a
  !!     non-portable fortran interface made available by
  !!     <TT>nsclock.c</TT>. This approach will eventually be extended to other
  !!     platforms.
  !!
  !!     New behaviour added at the Havana release allows the user to embed
  !!     profiling calls at varying levels of granularity all over the code,
  !!     and for any particular run, set a threshold of granularity so that
  !!     finer-grained clocks become dormant.
  !!
  !!     The threshold granularity is held in the private module variable
  !!     <TT>clock_grain</TT>. This value may be modified by the call
  !!     <TT>mpp_clock_set_grain</TT>, and affect clocks initiated by
  !!     subsequent calls to <TT>mpp_clock_id</TT>. The value of
  !!     <TT>clock_grain</TT> is set to an arbitrarily large number initially.
  !!
  !!     Clocks initialized by <TT>mpp_clock_id</TT> can set a new optional
  !!     argument <TT>grain</TT> setting their granularity level. Clocks check
  !!     this level against the current value of <TT>clock_grain</TT>, and are
  !!     only triggered if they are <I>at or below ("coarser than")</I> the
  !!     threshold. Finer-grained clocks are dormant for that run.
  !!
  !!The following grain levels are pre-defined:
  !!
  !!<pre>
  !!
  !!
  !!  integer, parameter, public :: CLOCK_COMPONENT=1 !component level, e.g model, exchange
  !!  integer, parameter, public :: CLOCK_SUBCOMPONENT=11 !top level within a model component, e.g dynamics, physics
  !!  integer, parameter, public :: CLOCK_MODULE=21 !module level, e.g main subroutine of a physics module
  !!  integer, parameter, public :: CLOCK_ROUTINE=31 !level of individual subroutine or function
  !!  integer, parameter, public :: CLOCK_LOOP=41 !loops or blocks within a routine
  !!  integer, parameter, public :: CLOCK_INFRA=51 !infrastructure level, e.g halo update
  !!</pre>
  !!
  !!     Note that subsequent changes to <TT>clock_grain</TT> do not
  !!     change the status of already initiated clocks, and that if the
  !!     optional <TT>grain</TT> argument is absent, the clock is always
  !!     triggered. This guarantees backward compatibility.
  subroutine mpp_clock_set_grain( grain )
    integer, intent(in) :: grain
    !set the granularity of times: only clocks whose grain is lower than
    !clock_grain are triggered, finer-grained clocks are dormant.
    !clock_grain is initialized to CLOCK_LOOP, so all clocks above the loop level
    !are triggered if this is never called.
    if( .NOT.module_is_initialized )call mpp_error( fatal, 'MPP_CLOCK_SET_GRAIN: You must first call mpp_init.' )
 
    clock_grain = grain
    return
  end subroutine mpp_clock_set_grain
 
  !#####################################################################
  subroutine clock_init( id, name, flags, grain )
    integer,           intent(in) :: id
    character(len=*),  intent(in) :: name
    integer, intent(in), optional :: flags, grain
    integer                       :: i
 
    clocks(id)%name = name
    clocks(id)%hits = 0
    clocks(id)%tick = 0
    clocks(id)%total_ticks = 0
    clocks(id)%sync_on_begin = .false.
    clocks(id)%detailed      = .false.
    clocks(id)%peset_num = current_peset_num
    if( PRESENT(flags) )then
       if( btest(flags,0) )clocks(id)%sync_on_begin = .true.
       if( btest(flags,1) )clocks(id)%detailed      = .true.
    end if
    clocks(id)%grain = 0
    if( PRESENT(grain) )clocks(id)%grain = grain
    if( clocks(id)%detailed )then
       allocate( clocks(id)%events(max_event_types) )
       clocks(id)%events(event_allreduce)%name = 'ALLREDUCE'
       clocks(id)%events(event_broadcast)%name = 'BROADCAST'
       clocks(id)%events(event_recv)%name = 'RECV'
       clocks(id)%events(event_send)%name = 'SEND'
       clocks(id)%events(event_wait)%name = 'WAIT'
       do i=1,max_event_types
          clocks(id)%events(i)%ticks(:) = 0
          clocks(id)%events(i)%bytes(:) = 0
          clocks(id)%events(i)%calls = 0
       end do
       clock_summary(id)%name = name
       clock_summary(id)%event(event_allreduce)%name = 'ALLREDUCE'
       clock_summary(id)%event(event_broadcast)%name = 'BROADCAST'
       clock_summary(id)%event(event_recv)%name = 'RECV'
       clock_summary(id)%event(event_send)%name = 'SEND'
       clock_summary(id)%event(event_wait)%name = 'WAIT'
       do i=1,max_event_types
          clock_summary(id)%event(i)%msg_size_sums(:) = 0.0
          clock_summary(id)%event(i)%msg_time_sums(:) = 0.0
          clock_summary(id)%event(i)%total_data = 0.0
          clock_summary(id)%event(i)%total_time = 0.0
          clock_summary(id)%event(i)%msg_size_cnts(:) = 0
          clock_summary(id)%event(i)%total_cnts = 0
       end do
    end if
    return
  end subroutine clock_init
 
  !#####################################################################
  !> Return an ID for a new or existing clock
  function mpp_clock_id( name, flags, grain )
    integer                       :: mpp_clock_id
    character(len=*),  intent(in) :: name
    integer, intent(in), optional :: flags, grain
 
    if( .NOT.module_is_initialized )call mpp_error( fatal, 'MPP_CLOCK_ID: You must first call mpp_init.')
 
    !if grain is present, the clock is only triggered if it
    !is low ("coarse") enough: compared to clock_grain
    !finer-grained clocks are dormant.
    !if grain is absent, clock is triggered.
    if( PRESENT(grain) )then
       if( grain.GT.clock_grain )then
          mpp_clock_id = 0
          return
       end if
    end if
    mpp_clock_id = 1
 
    if( clock_num.EQ.0 )then  !first
       clock_num = mpp_clock_id
       call clock_init(mpp_clock_id,name,flags)
    else
       find_clock: do while( trim(name).NE.trim(clocks(mpp_clock_id)%name) )
          mpp_clock_id = mpp_clock_id + 1
          if( mpp_clock_id.GT.clock_num )then
             if( mpp_clock_id.GT.max_clocks )then
                call mpp_error( fatal, 'MPP_CLOCK_ID: too many clock requests, ' // &
                      'check your clock id request or increase MAX_CLOCKS.')
             else               !new clock: initialize
                clock_num = mpp_clock_id
                call clock_init(mpp_clock_id,name,flags,grain)
                exit find_clock
             end if
          end if
       end do find_clock
    endif
    return
  end function mpp_clock_id
 
  !#####################################################################
  subroutine mpp_clock_begin(id)
    integer, intent(in) :: id
 
    if( .NOT.module_is_initialized )call mpp_error( fatal, 'MPP_CLOCK_BEGIN: You must first call mpp_init.' )
    if( .not. mpp_record_timing_data)return
    if( id.EQ.0 )return
    if( id.LT.0 .OR. id.GT.clock_num )call mpp_error( fatal, 'MPP_CLOCK_BEGIN: invalid id.' )
 
!$OMP MASTER
    if( clocks(id)%peset_num.NE.current_peset_num ) &
         call mpp_error( fatal, 'MPP_CLOCK_BEGIN: cannot change pelist context of a clock.' )
    if( clocks(id)%is_on) call mpp_error(fatal, 'MPP_CLOCK_BEGIN: mpp_clock_begin is called again '// &
                'before calling mpp_clock_end for the clock '//trim(clocks(id)%name) )
    if( clocks(id)%sync_on_begin .OR. sync_all_clocks )then
       !do an untimed sync at the beginning of the clock
       !this puts all PEs in the current pelist on par, so that measurements begin together
       !ending time will be different, thus measuring load imbalance for this clock.
       call mpp_sync()
    end if
 
    if (debug) then
      num_clock_ids = num_clock_ids+1
      if(num_clock_ids > max_clocks)call mpp_error(fatal,'MPP_CLOCK_BEGIN: max num previous_clock exceeded.' )
      previous_clock(num_clock_ids) = current_clock
      current_clock = id
    endif
    call system_clock( clocks(id)%tick )
    clocks(id)%hits = clocks(id)%hits + 1
    clocks(id)%is_on = .true.
!$OMP END MASTER
    return
  end subroutine mpp_clock_begin
 
  !#####################################################################
  subroutine mpp_clock_end(id)
    integer, intent(in) :: id
    integer(i8_kind)  :: delta
    integer             :: errunit
 
    if( .NOT.module_is_initialized )call mpp_error( fatal, 'MPP_CLOCK_END: You must first call mpp_init.' )
    if( .not. mpp_record_timing_data)return
    if( id.EQ.0 )return
    if( id.LT.0 .OR. id.GT.clock_num )call mpp_error( fatal, 'MPP_CLOCK_BEGIN: invalid id.' )
!$OMP MASTER
    if( .NOT. clocks(id)%is_on) call mpp_error(fatal, 'MPP_CLOCK_END: mpp_clock_end is called '// &
                'before calling mpp_clock_begin for the clock '//trim(clocks(id)%name) )
 
    call system_clock(end_tick)
    if( clocks(id)%peset_num.NE.current_peset_num ) &
         call mpp_error( fatal, 'MPP_CLOCK_END: cannot change pelist context of a clock.' )
    delta = end_tick - clocks(id)%tick
    if( delta.LT.0 )then
       errunit = stderr()
       write( errunit,* )'pe, id, start_tick, end_tick, delta, max_ticks=', pe, id, clocks(id)%tick, end_tick, &
            &  delta, max_ticks
       delta = delta + max_ticks + 1
       call mpp_error( warning, 'MPP_CLOCK_END: Clock rollover, assumed single roll.' )
    end if
    clocks(id)%total_ticks = clocks(id)%total_ticks + delta
    if (debug) then
      if(num_clock_ids < 1) call mpp_error(note,'MPP_CLOCK_END: min num previous_clock < 1.' )
      current_clock = previous_clock(num_clock_ids)
      num_clock_ids = num_clock_ids-1
    endif
    clocks(id)%is_on = .false.
!$OMP END MASTER
    return
  end subroutine mpp_clock_end
 
 !#####################################################################
  subroutine mpp_record_time_start()
 
     mpp_record_timing_data = .true.
 
  end subroutine mpp_record_time_start
 
  !#####################################################################
  subroutine mpp_record_time_end()
 
     mpp_record_timing_data = .false.
 
  end subroutine mpp_record_time_end
 
 
  !#####################################################################
  subroutine increment_current_clock( event_id, bytes )
    integer,           intent(in) :: event_id
    integer, intent(in), optional :: bytes
    integer                       :: n
    integer(i8_kind)            :: delta
    integer                       :: errunit
 
    if( .not. mpp_record_timing_data )return
    if( .not.debug .or. (current_clock.EQ.0) )return
    if( current_clock.LT.0 .OR. current_clock.GT.clock_num )call mpp_error( fatal, &
       &  'MPP_CLOCK_BEGIN: invalid current_clock.' )
    if( .NOT.clocks(current_clock)%detailed )return
    call system_clock(end_tick)
    n = clocks(current_clock)%events(event_id)%calls + 1
 
    if( n.EQ.max_events )call mpp_error( warning, &
         'MPP_CLOCK: events exceed MAX_EVENTS, ignore detailed profiling data for clock '// &
         & trim(clocks(current_clock)%name) )
    if( n.GT.max_events )return
 
    clocks(current_clock)%events(event_id)%calls = n
    delta = end_tick - start_tick
    if( delta.LT.0 )then
       errunit = stderr()
       write( errunit,* )'pe, event_id, start_tick, end_tick, delta, max_ticks=', &
                           pe, event_id, start_tick, end_tick, delta, max_ticks
       delta = delta + max_ticks + 1
       call mpp_error( warning, 'MPP_CLOCK_END: Clock rollover, assumed single roll.' )
    end if
    clocks(current_clock)%events(event_id)%ticks(n) = delta
    if( PRESENT(bytes) )clocks(current_clock)%events(event_id)%bytes(n) = bytes
    return
  end subroutine increment_current_clock
 
  !#####################################################################
 
  subroutine dump_clock_summary()
 
    real              :: total_time,total_time_all,total_data
    real              :: msg_size,eff_bw,s
    integer           :: sd_unit, total_calls
    integer           :: j,k,ct, msg_cnt
    character(len=2)  :: u
    character(len=FMS_FILE_LEN) :: filename
    character(len=20),dimension(MAX_BINS),save :: bin
 
    data bin( 1)  /'  0   -    8    B:  '/
    data bin( 2)  /'  8   -   16    B:  '/
    data bin( 3)  /' 16   -   32    B:  '/
    data bin( 4)  /' 32   -   64    B:  '/
    data bin( 5)  /' 64   -  128    B:  '/
    data bin( 6)  /'128   -  256    B:  '/
    data bin( 7)  /'256   -  512    B:  '/
    data bin( 8)  /'512   - 1024    B:  '/
    data bin( 9)  /'  1.0 -    2.1 KB:  '/
    data bin(10)  /'  2.1 -    4.1 KB:  '/
    data bin(11)  /'  4.1 -    8.2 KB:  '/
    data bin(12)  /'  8.2 -   16.4 KB:  '/
    data bin(13)  /' 16.4 -   32.8 KB:  '/
    data bin(14)  /' 32.8 -   65.5 KB:  '/
    data bin(15)  /' 65.5 -  131.1 KB:  '/
    data bin(16)  /'131.1 -  262.1 KB:  '/
    data bin(17)  /'262.1 -  524.3 KB:  '/
    data bin(18)  /'524.3 - 1048.6 KB:  '/
    data bin(19)  /'  1.0 -    2.1 MB:  '/
    data bin(20)  /' >2.1          MB:  '/
 
    if( .NOT.any(clocks(1:clock_num)%detailed) )return
    write( filename,'(a,i6.6)' )'mpp_clock.out.', pe
 
    open(newunit=sd_unit,file=trim(filename),form='formatted')
 
    comm_type: do ct = 1,clock_num
 
       if( .NOT.clocks(ct)%detailed )cycle
       write(sd_unit,*) &
            clock_summary(ct)%name(1:15),' Communication Data for PE ',pe
 
       write(sd_unit,*) ' '
       write(sd_unit,*) ' '
 
       total_time_all = 0.0
       event_type: do k = 1,max_event_types-1
 
          if(clock_summary(ct)%event(k)%total_time == 0.0)cycle
 
          total_time = clock_summary(ct)%event(k)%total_time
          total_time_all = total_time_all + total_time
          total_data = clock_summary(ct)%event(k)%total_data
          total_calls = int(clock_summary(ct)%event(k)%total_cnts)
 
          write(sd_unit,1000) clock_summary(ct)%event(k)%name(1:9) // ':'
 
          write(sd_unit,1001) 'Total Data: ',total_data*1.0e-6, &
               'MB; Total Time: ', total_time, &
               'secs; Total Calls: ',total_calls
 
          write(sd_unit,*) ' '
          write(sd_unit,1002) '     Bin            Counts      Avg Size        Eff B/W'
          write(sd_unit,*) ' '
 
          bin_loop: do j=1,max_bins
 
             if(clock_summary(ct)%event(k)%msg_size_cnts(j)==0)cycle
 
             if(j<=8)then
                s = 1.0
                u = ' B'
             elseif(j<=18)then
                s = 1.0e-3
                u = 'KB'
             else
                s = 1.0e-6
                u = 'MB'
             endif
 
             msg_cnt = int(clock_summary(ct)%event(k)%msg_size_cnts(j))
             msg_size = &
                  s*(clock_summary(ct)%event(k)%msg_size_sums(j)/real(msg_cnt))
             eff_bw = (1.0e-6)*( clock_summary(ct)%event(k)%msg_size_sums(j) / &
                  clock_summary(ct)%event(k)%msg_time_sums(j) )
 
             write(sd_unit,1003) bin(j),msg_cnt,msg_size,u,eff_bw
 
          end do bin_loop
 
          write(sd_unit,*) ' '
          write(sd_unit,*) ' '
       end do event_type
 
       ! "Data-less" WAIT
 
       if(clock_summary(ct)%event(max_event_types)%total_time>0.0)then
 
          total_time = clock_summary(ct)%event(max_event_types)%total_time
          total_time_all = total_time_all + total_time
          total_calls = int(clock_summary(ct)%event(max_event_types)%total_cnts)
 
          write(sd_unit,1000) clock_summary(ct)%event(max_event_types)%name(1:9) // ':'
 
          write(sd_unit,1004) 'Total Calls: ',total_calls,'; Total Time: ', &
               total_time,'secs'
 
       endif
 
       write(sd_unit,*) ' '
       write(sd_unit,1005) 'Total communication time spent for ' // &
            clock_summary(ct)%name(1:9) // ': ',total_time_all,'secs'
       write(sd_unit,*) ' '
       write(sd_unit,*) ' '
       write(sd_unit,*) ' '
 
    end do comm_type
 
    close(sd_unit)
 
1000 format(a)
1001 format(a,f8.2,a,f8.2,a,i6)
1002 format(a)
1003 format(a,i6,'    ','  ',f9.1,a,'    ',f9.2,'MB/sec')
1004 format(a,i8,a,f9.2,a)
1005 format(a,f9.2,a)
    return
  end subroutine dump_clock_summary
 
  !#####################################################################
 
  integer function get_unit()
 
    integer,save :: i
    logical      :: l_open
 
    if (pe == root_pe) call mpp_error(warning, &
        'get_unit is deprecated and will be removed in a future release, please use the Fortran intrinsic newunit')
    do i=10,99
       inquire(unit=i,opened=l_open)
       if(.not.l_open)exit
    end do
 
    if(i==100)then
       call mpp_error(fatal,'Unable to get I/O unit')
    else
       get_unit = i
    endif
 
    return
  end function get_unit
 
  !#####################################################################
 
  subroutine sum_clock_data()
 
    integer :: i,j,k,ct,event_size,event_cnt
    real    :: msg_time
 
    clock_type: do ct=1,clock_num
       if( .NOT.clocks(ct)%detailed )cycle
       event_type: do j=1,max_event_types-1
          event_cnt = clocks(ct)%events(j)%calls
          event_summary: do i=1,event_cnt
 
             clock_summary(ct)%event(j)%total_cnts = &
                  clock_summary(ct)%event(j)%total_cnts + 1
 
             event_size = int(clocks(ct)%events(j)%bytes(i))
 
             k = find_bin(event_size)
 
             clock_summary(ct)%event(j)%msg_size_cnts(k) = &
                  clock_summary(ct)%event(j)%msg_size_cnts(k) + 1
 
             clock_summary(ct)%event(j)%msg_size_sums(k) = &
                  clock_summary(ct)%event(j)%msg_size_sums(k) &
                  + clocks(ct)%events(j)%bytes(i)
 
             clock_summary(ct)%event(j)%total_data = &
                  clock_summary(ct)%event(j)%total_data &
                  + clocks(ct)%events(j)%bytes(i)
 
             msg_time = clocks(ct)%events(j)%ticks(i)
             msg_time = tick_rate * real( clocks(ct)%events(j)%ticks(i) )
 
             clock_summary(ct)%event(j)%msg_time_sums(k) = &
                  clock_summary(ct)%event(j)%msg_time_sums(k) + msg_time
 
             clock_summary(ct)%event(j)%total_time = &
                  clock_summary(ct)%event(j)%total_time + msg_time
 
          end do event_summary
       end do event_type
 
       j = max_event_types ! WAITs
       ! "msg_size_cnts" doesn't really mean anything for WAIT
       ! but position will be used to store number of counts for now.
 
       event_cnt = clocks(ct)%events(j)%calls
       clock_summary(ct)%event(j)%msg_size_cnts(1) = event_cnt
       clock_summary(ct)%event(j)%total_cnts       = event_cnt
 
       msg_time = tick_rate * real( sum( clocks(ct)%events(j)%ticks(1:event_cnt) ) )
       clock_summary(ct)%event(j)%msg_time_sums(1) = &
            clock_summary(ct)%event(j)%msg_time_sums(1) + msg_time
 
       clock_summary(ct)%event(j)%total_time = clock_summary(ct)%event(j)%msg_time_sums(1)
 
    end do clock_type
 
    return
  contains
    integer function find_bin(event_size)
 
      integer,intent(in) :: event_size
      integer            :: k,msg_size
 
      msg_size = 8
      k = 1
      do while(event_size>msg_size .and. k<max_bins)
         k = k+1
         msg_size = msg_size*2
      end do
      find_bin = k
      return
    end function find_bin
 
  end subroutine sum_clock_data
 
  !#####################################################################
  !> This routine will double the size of peset and copy the original peset data
  !! into the expanded one. The maximum allowed to expand is PESET_MAX.
  subroutine expand_peset()
     integer :: old_peset_max,n
     type(communicator), allocatable :: peset_old(:)
 
     old_peset_max = current_peset_max
     if(old_peset_max .GE. peset_max) call mpp_error(fatal, &
         "mpp_mod(expand_peset): the number of peset reached PESET_MAX, increase PESET_MAX or contact developer")
 
     ! copy data to a tempoary data
     allocate(peset_old(0:old_peset_max))
     do n = 0, old_peset_max
        peset_old(n)%count      = peset(n)%count
        peset_old(n)%comm       = peset(n)%comm
        peset_old(n)%group      = peset(n)%group
        peset_old(n)%name       = peset(n)%name
        peset_old(n)%start      = peset(n)%start
        peset_old(n)%log2stride = peset(n)%log2stride
 
        if( ASSOCIATED(peset(n)%list) ) then
           allocate(peset_old(n)%list(size(peset(n)%list(:))) )
           peset_old(n)%list(:) = peset(n)%list(:)
           deallocate(peset(n)%list)
        endif
     enddo
     deallocate(peset)
 
     ! create the new peset
     current_peset_max = min(peset_max, 2*old_peset_max)
     allocate(peset(0:current_peset_max))
     peset(:)%count = -1
     peset(:)%comm = mpi_comm_null
     peset(:)%group = mpi_group_null
     peset(:)%start = -1
     peset(:)%log2stride = -1
     peset(:)%name  = " "
     do n = 0, old_peset_max
        peset(n)%count      = peset_old(n)%count
        peset(n)%comm       = peset_old(n)%comm
        peset(n)%group      = peset_old(n)%group
        peset(n)%name       = peset_old(n)%name
        peset(n)%start      = peset_old(n)%start
        peset(n)%log2stride = peset_old(n)%log2stride
 
        if( ASSOCIATED(peset_old(n)%list) ) then
           allocate(peset(n)%list(size(peset_old(n)%list(:))) )
           peset(n)%list(:) = peset_old(n)%list(:)
           deallocate(peset_old(n)%list)
        endif
     enddo
     deallocate(peset_old)
 
     call mpp_error(note, "mpp_mod(expand_peset): size of peset is expanded to ", current_peset_max)
 
  end subroutine expand_peset
  !#####################################################################
 
  function uppercase (cs)
    character(len=*), intent(in) :: cs
    character(len=len(cs)),target       :: uppercase
    integer                      :: k,tlen
    character, pointer :: ca
    integer, parameter :: co=iachar('A')-iachar('a') ! case offset
    !The transfer function truncates the string with xlf90_r
    tlen = len_trim(cs)
    if(tlen <= 0) then      ! catch IBM compiler bug
       uppercase = cs  ! simply return input blank string
    else
    uppercase = cs(1:tlen)
    do k=1, tlen
       ca => uppercase(k:k)
       if(ca >= "a" .and. ca <= "z") ca = achar(ichar(ca)+co)
    enddo
    endif
  end function uppercase
 
!#######################################################################
 
  function lowercase (cs)
    character(len=*), intent(in) :: cs
    character(len=len(cs)),target       :: lowercase
    integer, parameter :: co=iachar('a')-iachar('A') ! case offset
    integer                        :: k,tlen
    character, pointer :: ca
!  The transfer function truncates the string with xlf90_r
    tlen = len_trim(cs)
    if(tlen <= 0) then      ! catch IBM compiler bug
       lowercase = cs  ! simply return input blank string
    else
    lowercase = cs(1:tlen)
    do k=1, tlen
       ca => lowercase(k:k)
       if(ca >= "A" .and. ca <= "Z") ca = achar(ichar(ca)+co)
    enddo
    endif
  end function lowercase
 
 
  !#######################################################################
 
!-----------------------------------------------------------------------
!
! AUTHOR: Rusty Benson (rusty.benson@noaa.gov)
!
!
! THESE LINES MUST BE PRESENT IN MPP.F90
!
! ! public variable needed for reading an input nml file from an internal file
!   character(len=:), dimension(:), allocatable, public :: input_nml_file
!
 
!-----------------------------------------------------------------------
 
!> Reads an existing input nml file into a character array and broadcasts
!! it to the non-root mpi-tasks. This allows the use of reads from an
!! internal file for namelist settings (requires 2003 compliant compiler)
!!
!! read(input_nml_file, nml=<name_nml>, iostat=status)
!!
!!
  subroutine read_input_nml(pelist_name_in, alt_input_nml_path)
 
! Include variable "version" to be written to log file.
#include<file_version.h>
 
    character(len=*), intent(in), optional :: pelist_name_in
    character(len=*), intent(in), optional :: alt_input_nml_path
! private variables
    integer :: log_unit
    integer :: i
    integer, dimension(2) :: lines_and_length
    logical :: file_exist
    character(len=len(peset(current_peset_num)%name)) :: pelist_name
    character(len=FMS_PATH_LEN) :: filename
 
! check the status of input_nml_file
    if ( allocated(input_nml_file) ) then
      deallocate(input_nml_file)
    endif
 
! the following code is necessary for using alternate namelist files (nests, stretched grids, etc)
    if (PRESENT(pelist_name_in)) then
      ! test to make sure length of pelist_name_in is <= pelist_name
      if (len(pelist_name_in) > len(pelist_name)) then
        call mpp_error(fatal,  &
           "mpp_util.inc: read_input_nml optional argument pelist_name_in has size greater than local pelist_name")
      else
        pelist_name = pelist_name_in
      endif
    else
      pelist_name = mpp_get_current_pelist_name()
    endif
    filename='input_'//trim(pelist_name)//'.nml'
    inquire(file=filename, exist=file_exist)
    if (.not. file_exist ) then
       if (present(alt_input_nml_path)) then
          filename = alt_input_nml_path
       else
          filename = 'input.nml'
       end if
    endif
    lines_and_length = get_ascii_file_num_lines_and_length(filename)
    allocate(character(len=lines_and_length(2))::input_nml_file(lines_and_length(1)))
    call read_ascii_file(filename, lines_and_length(2), input_nml_file)
 
! write info logfile
    if (pe == root_pe) then
       log_unit = stdlog()
       write(log_unit,'(a)')  '========================================================================'
       write(log_unit,'(a)')  'READ_INPUT_NML: '//trim(version)
       write(log_unit,'(a)')  'READ_INPUT_NML: '//trim(filename)//' '
       do i = 1, lines_and_length(1)
          write(log_unit,*) trim(input_nml_file(i))
       enddo
    end if
  end subroutine read_input_nml
 
 
  !#######################################################################
  !z1l: This is extracted from read_ascii_file
  function get_ascii_file_num_lines(FILENAME, LENGTH, PELIST)
    character(len=*), intent(in) :: FILENAME
    integer, intent(in) :: LENGTH
    integer, intent(in), optional, dimension(:) :: PELIST
 
    integer :: num_lines, get_ascii_file_num_lines
    character(len=LENGTH) :: str_tmp
    character(len=5) :: text
    integer :: status, f_unit, from_pe
    logical :: file_exist
 
    if( read_ascii_file_on) then
       call mpp_error(fatal,  &
          "mpp_util.inc: get_ascii_file_num_lines is called again before calling read_ascii_file")
    endif
    read_ascii_file_on = .true.
 
    from_pe = root_pe
    get_ascii_file_num_lines = -1
    num_lines = -1
    if ( pe == root_pe ) then
       inquire(file=filename, exist=file_exist)
 
       if ( file_exist ) then
          open(newunit=f_unit, file=filename, action='READ', status='OLD', iostat=status)
 
          if ( status .ne. 0 ) then
             write (unit=text, fmt='(I5)') status
             call mpp_error(fatal, 'get_ascii_file_num_lines: Error opening file:' //trim(filename)// &
                            '.  (IOSTAT = '//trim(text)//')')
          else
             num_lines = 1
             do
                read (unit=f_unit, fmt='(A)', iostat=status) str_tmp
                if ( status .lt. 0 ) then
                   ! deprecate num_lines by 1 and ensure num_lines is at least 1
                   num_lines = max(num_lines - 1, 1)
                   exit
                endif
                if ( status .gt. 0 ) then
                   write (unit=text, fmt='(I5)') num_lines
                   call mpp_error(fatal, 'get_ascii_file_num_lines: Error reading line '//trim(text)// &
                        ' in file '//trim(filename)//'.')
                end if
                if ( len_trim(str_tmp) == length ) then
                   write(unit=text, fmt='(I5)') length
                   call mpp_error(fatal, 'get_ascii_file_num_lines: Length of output string ('//trim(text)//&
                                       & ' is too small. Increase the LENGTH value.')
                end if
                num_lines = num_lines + 1
             end do
             close(unit=f_unit)
          end if
       else
          call mpp_error(fatal, 'get_ascii_file_num_lines: File '//trim(filename)//' does not exist.')
       end if
    end if
 
    ! Broadcast number of lines
    call mpp_broadcast(num_lines, from_pe, pelist=pelist)
    get_ascii_file_num_lines = num_lines
 
  end function get_ascii_file_num_lines
 
  !#######################################################################
  !> @brief Function to determine the maximum line length and number of lines from an ascii file
  function get_ascii_file_num_lines_and_length(FILENAME, PELIST)
    character(len=*), intent(in) :: filename !< name of the file to be read
    integer, intent(in), optional, dimension(:) :: pelist !< optional pelist
 
    integer, dimension(2) :: get_ascii_file_num_lines_and_length !< number of lines (1) and
                                                                 !! max line length (2)
    integer :: num_lines, max_length
    integer, parameter :: length=1024
    character(len=LENGTH) :: str_tmp
    character(len=5) :: text
    integer :: status, f_unit, from_pe
    logical :: file_exist
 
    if( read_ascii_file_on) then
       call mpp_error(fatal,  &
          "mpp_util.inc: get_ascii_file_num_lines is called again before calling read_ascii_file")
    endif
    read_ascii_file_on = .true.
 
    from_pe = root_pe
    get_ascii_file_num_lines_and_length = -1
    num_lines = -1
    max_length = -1
    if ( pe == root_pe ) then
       inquire(file=filename, exist=file_exist)
 
       if ( file_exist ) then
          open(newunit=f_unit, file=filename, action='READ', status='OLD', iostat=status)
 
          if ( status .ne. 0 ) then
             write (unit=text, fmt='(I5)') status
             call mpp_error(fatal, 'get_ascii_file_num_lines: Error opening file:' //trim(filename)// &
                            '.  (IOSTAT = '//trim(text)//')')
          else
             num_lines = 1
             max_length = 1
             do
                read (unit=f_unit, fmt='(A)', iostat=status) str_tmp
                if ( status .lt. 0 ) then
                   ! deprecate num_lines by 1 and ensure num_lines is at least 1
                   num_lines = max(num_lines - 1, 1)
                   exit
                endif
                if ( status .gt. 0 ) then
                   write (unit=text, fmt='(I5)') num_lines
                   call mpp_error(fatal, 'get_ascii_file_num_lines: Error reading line '//trim(text)// &
                        ' in file '//trim(filename)//'.')
                end if
                if ( len_trim(str_tmp) == length) then
                   write(unit=text, fmt='(I5)') length
                   call mpp_error(fatal, 'get_ascii_file_num_lines: Length of output string ('//trim(text)//&
                                       & ' is too small. Increase the LENGTH value.')
                end if
                if (len_trim(str_tmp) > max_length) max_length = len_trim(str_tmp)
                num_lines = num_lines + 1
             end do
             close(unit=f_unit)
          end if
       else
          call mpp_error(fatal, 'get_ascii_file_num_lines: File '//trim(filename)//' does not exist.')
       end if
       max_length = max_length+1
    end if
 
    ! Broadcast number of lines
    call mpp_broadcast(num_lines, from_pe, pelist=pelist)
    call mpp_broadcast(max_length, from_pe, pelist=pelist)
    get_ascii_file_num_lines_and_length(1) = num_lines
    get_ascii_file_num_lines_and_length(2) = max_length
 
  end function get_ascii_file_num_lines_and_length
 
  !-----------------------------------------------------------------------
  !
  ! AUTHOR: Rusty Benson <rusty.benson@noaa.gov>,
  !         Seth Underwood <Seth.Underwood@noaa.gov>
  !
  !-----------------------------------------------------------------------
  ! subroutine READ_ASCII_FILE
  !
  !
  !> Reads any ascii file into a character array and broadcasts
  !! it to the non-root mpi-tasks.  Based off READ_INPUT_NML.
  !!
  !! Passed in 'Content' array, must be of the form:
  !! character(len=LENGTH), dimension(:), allocatable :: array_name
  !!
  !! Reads from this array must be done in a do loop over the number of
  !! lines, i.e.:
  !!
  !! do i=1, num_lines
  !!    read (UNIT=array_name(i), FMT=*) var1, var2, ...
  !! end do
  subroutine read_ascii_file(FILENAME, LENGTH, Content, PELIST)
    character(len=*),    intent(in)               :: FILENAME
    integer,             intent(in)               :: LENGTH
    character(len=*), intent(inout), dimension(:) :: Content
    integer, intent(in), optional,   dimension(:) :: PELIST
 
    ! Include variable "version" to be written to log file.
#include<file_version.h>
 
    character(len=5) :: text
    logical :: file_exist
    integer :: status, f_unit, log_unit
    integer :: from_pe
    integer :: pnum_lines, num_lines
    character(len=LENGTH) :: str_tmp !< Temporary variable to store line from file
 
    if( .NOT. read_ascii_file_on) then
       call mpp_error(fatal,  &
          "mpp_util.inc: get_ascii_file_num_lines needs to be called before calling read_ascii_file")
    endif
    read_ascii_file_on = .false.
 
    from_pe = root_pe
    num_lines = size(content(:))
 
    if ( pe == root_pe ) then
       ! write info logfile
       log_unit = stdlog()
       write(log_unit,'(a)')  '========================================================================'
       write(log_unit,'(a)')  'READ_ASCII_FILE: '//trim(version)
       write(log_unit,'(a)')  'READ_ASCII_FILE: File: '//trim(filename)
 
       inquire(file=filename, exist=file_exist)
 
       if ( file_exist ) then
          open(newunit=f_unit, file=filename, action='READ', status='OLD', iostat=status)
 
          if ( status .ne. 0 ) then
             write (unit=text, fmt='(I5)') status
             call mpp_error(fatal, 'READ_ASCII_FILE: Error opening file: '// &
                            & trim(filename)//'.  (IOSTAT = '//trim(text)//')')
          else
 
             if ( num_lines .gt. 0 ) then
                content(:) = ' '
 
                rewind(unit=f_unit, iostat=status)
                if ( status .ne. 0 ) then
                   write (unit=text, fmt='(I5)') status
                   call mpp_error(fatal, 'READ_ASCII_FILE: Unable to re-read file '//trim(filename)//'. (IOSTAT = '&
                        //trim(text)//'.')
                else
                   ! A second 'sanity' check on the file
                   pnum_lines = 1
 
                   do
                      read (unit=f_unit, fmt='(A)', iostat=status) str_tmp
 
                      if ( status .lt. 0 ) then
                         ! deprecate pnum_lines by 1 and ensure pnum_lines is at least 1
                         pnum_lines = max(pnum_lines - 1, 1)
                         exit
                      endif
                      if ( status .gt. 0 ) then
                         write (unit=text, fmt='(I5)') pnum_lines
                         call mpp_error(fatal, 'READ_ASCII_FILE: Error reading line '// &
                                        & trim(text)//' in file '//trim(filename)//'.')
                      end if
                      if(pnum_lines > num_lines) then
                         call mpp_error(fatal, 'READ_ASCII_FILE: number of lines in file '//trim(filename)// &
                                ' is greater than size(Content(:)). ')
                      end if
                      if ( len_trim(str_tmp) == length ) then
                         write(unit=text, fmt='(I5)') length
                         call mpp_error(fatal, 'READ_ASCII_FILE: Length of output string ('//trim(text)// &
                                             & ' is too small. Increase the LENGTH value.')
                      end if
                      content(pnum_lines) = str_tmp
                      pnum_lines = pnum_lines + 1
                   end do
                   if(num_lines .NE. pnum_lines) then
                      call mpp_error(fatal, 'READ_ASCII_FILE: number of lines in file '//trim(filename)// &
                          ' does not equal to size(Content(:)) ' )
                   end if
                end if
             end if
             close(unit=f_unit)
          end if
       else
          call mpp_error(fatal, 'READ_ASCII_FILE: File '//trim(filename)//' does not exist.')
       end if
    end if
 
    ! Broadcast character array
    call mpp_broadcast(content, length, from_pe, pelist=pelist)
 
  end subroutine read_ascii_file
 
  !> @brief Produce an inverse permutation. For example, transform [2, 3, 1, 4] to [3, 1, 2, 4].
  !!
  !! The purpose of this subroutine is to convert between (memory dimension) -> (logical dimension) and
  !! (logical dimension) -> (memory dimension) maps.
  !!
  !! @param [in] <x> The original permutation vector
  !! @param [out] <y> The inverted permutation vector
  subroutine inverse_permutation(x, y)
    integer, intent(in) :: x(:)
    integer, intent(out) :: y(size(x))
    integer :: i, n
 
    y = 0
 
    n = size(x)
    do i=1,n
      if (x(i).ge.1 .and. x(i).le.n) then
        y(x(i)) = i
      else
        block
          character(2) :: nstr
 
          write (nstr, "(I0)") n
          call mpp_error(fatal, "inverse_permutation: Invalid dimension map. &
                                 Values must be in the range from 1 to " // trim(nstr) // ".")
        end block
      endif
    enddo
 
    if (any(y.eq.0)) then
      call mpp_error(fatal, "inverse_permutation: Invalid dim_order. Values must be non-repeating.")
    endif
  end subroutine inverse_permutation
!> @}