diff --git a/CHANGELOG.md b/CHANGELOG.md
index 51e7b7f1b51..f28f4c2a765 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -86,6 +86,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 
+- Fixed the behavior of MAPL_MaxMin in presence of NaN
 - Fixed bug with return codes and macros in udunits2f
 
 ### Removed
diff --git a/base/MAPL_MaxMinMod.F90 b/base/MAPL_MaxMinMod.F90
index 8aa0bef53b4..e63f4586619 100644
--- a/base/MAPL_MaxMinMod.F90
+++ b/base/MAPL_MaxMinMod.F90
@@ -8,7 +8,7 @@
 !
 ! Author: GMAO SI-Team
 !
-! `MAPL_MaxMinMo` --- Global Max/Min of Arrays
+! `MAPL_MaxMinMod` --- Global Max/Min of Arrays
 !
 ! This module implements functions for calculating/printing out the global min/max
 ! of fortran arrays. Derived from GEOS-4 pmaxmin() functions.
@@ -20,6 +20,7 @@ module MAPL_MaxMinMod
 
       Use ESMF
       Use MAPL_CommsMod
+      Use mpi
 
       implicit None
 
@@ -66,21 +67,19 @@ subroutine pmaxmin2d_r4 ( qname, a, pmin_, pmax_, fac_ )
 
       real(ESMF_KIND_R4) :: pmax, pmin, fac
 
-      integer :: im, jt
-
       integer :: i, j, two=2
 
-      real, allocatable :: qmin(:), qmax(:)
       real pm1(2)
       real pm_res(2)
       type(ESMF_VM) :: vm
 
       character(len=16) :: name
+!NOTE: the current version does not trap error conditions returned in status
       integer :: status
-
-      im = size(a,1)
-      jt = size(a,2)
-      allocate(qmin(jt),qmax(jt))
+      integer :: comm
+      logical :: has_nans
+      logical :: has_nans_local
+      character(len=:), allocatable :: buf
 
       if ( present(fac_) ) then
          fac = fac_
@@ -90,25 +89,10 @@ subroutine pmaxmin2d_r4 ( qname, a, pmin_, pmax_, fac_ )
 
       call ESMF_VmGetCurrent(vm=vm, rc=status)
 
-      do j=1,jt
-         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
-!
-! 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
+      has_nans_local = any(a /= a)
+
+      pmin = minval(a, mask=(a==a))
+      pmax = maxval(a, mask=(a==a))
 
       pm1(1) = pmax
       pm1(2) = -pmin
@@ -118,13 +102,19 @@ subroutine pmaxmin2d_r4 ( qname, a, pmin_, pmax_, fac_ )
 
       if ( present(pmax_) ) pmax_ = pmax
       if ( present(pmin_) ) pmin_ = pmin
-      deallocate(qmax,qmin)
+
+      call ESMF_VmGet(VM, mpicommunicator=comm, rc=status)
+      call MPI_Reduce(has_nans_local, has_nans, 1, MPI_LOGICAL, MPI_LOR, 0, comm, status)
 
       if ( fac /= 0.0 ) then  ! trick to prevent printing
-         if ( MAPL_am_I_root() ) then
+         if ( MAPL_am_I_root(vm) ) then
             name = '            '
             name(1:len(qname)) = qname
-            write(*,*) name, ' max = ', pmax*fac, ' min = ', pmin*fac
+            buf = ""
+            if (has_nans) then
+               buf = " has NaN"
+            end if
+            write(*,*) name, ' max = ', pmax*fac, ' min = ', pmin*fac, buf
             return
          end if
       end if