Added utilities directory with array and regex utilities

CMakeLists.txt

@@ -236,12 +236,12 @@ if (BUILD_WITH_FARGPARSE)
    add_subdirectory (docs)
    add_subdirectory (benchmarks)
 endif()
-
 add_subdirectory (geom)
 add_subdirectory (vertical_grid)
 add_subdirectory (regridder_mgr)
 add_subdirectory (hconfig)
 add_subdirectory (hconfig_utils)
+add_subdirectory (utilities)
 
 if (PFUNIT_FOUND)
   include (add_pfunit_ctest)

utilities/CMakeLists.txt

@@ -0,0 +1,13 @@
+esma_set_this (OVERRIDE MAPL.utilities)
+
+set(srcs utilities.F90)
+
+esma_add_library(${this} SRCS ${srcs} DEPENDENCIES MAPL.shared TYPE SHARED)
+target_include_directories (${this} PUBLIC $<BUILD_INTERFACE:${MAPL_SOURCE_DIR}/include>)
+
+add_subdirectory(arrays)
+add_subdirectory(regex)
+
+if (PFUNIT_FOUND)
+  add_subdirectory(tests EXCLUDE_FROM_ALL)
+endif ()

utilities/arrays/AreaMean.F90

@@ -0,0 +1,149 @@
+#include "MAPL.h"
+
+module mapl3g_AreaMean
+
+   use mpi
+   use, intrinsic :: iso_fortran_env, only: real32, real64
+   use MAPL_Constants, only: MAPL_UNDEFINED_REAL, MAPL_UNDEFINED_REAL64
+   use MAPL_ErrorHandling, only: MAPL_Verify, MAPL_Return, MAPL_Assert
+
+   implicit none
+   private
+
+   public :: AreaMean
+
+   interface AreaMean
+      ! module procedure AreaMean_2d_r8_bitrep
+      module procedure AreaMean_2d_r8
+   end interface AreaMean
+
+contains
+
+   ! subroutine AreaMean_2d_r8_bitrep( qave, q, area, grid, bitreproducible, rc )
+   !    real(kind=real64), intent(  OUT) :: qave
+   !    real,              intent(IN   ) :: q(:,:)
+   !    real,              intent(IN   ) :: area(:,:)
+   !    type(ESMF_Grid),   intent(INout) :: grid
+   !    logical,           intent(IN   ) :: bitreproducible
+   !    integer, optional, intent(  OUT) :: rc
+
+   !    ! Log err vars
+   !    integer :: status
+   !    character(len=ESMF_MAXSTR), parameter :: Iam='AreaMeanBR'
+
+   !    ! Local vars
+   !    real(kind=real64)  :: qdum(2)
+   !    integer :: im,jm
+   !    integer :: DIMS(3)
+
+   !    integer :: i,j
+   !    logical :: amIRoot
+
+   !    real, allocatable :: qglobal(:,:)
+   !    real, allocatable :: aglobal(:,:)
+
+   !    type(ESMF_VM) :: vm
+
+   !    if (.not. bitreproducible) then
+   !       call AreaMean(qave, q, area, grid, rc=status )
+   !       _RETURN(STATUS)
+   !    end if
+
+   !    ! get VM (should get from the grid, but this is quicker)
+   !    call ESMF_VmGetCurrent(vm, rc=status)
+   !    _VERIFY(STATUS)
+
+   !    amIRoot = MAPL_AM_I_Root(vm)
+
+   !    if (amIRoot) then
+   !       call MAPL_GridGet(GRID, globalCellCountPerDim=DIMS, RC=STATUS)
+   !       im = DIMS(1) ! global grid dim
+   !       jm = DIMS(2) ! global grid dim
+   !    else
+   !       im = 0 ! dummy sizes
+   !       jm = 0 ! dummy sizes
+   !    end if
+
+   !    allocate(qglobal(im,jm), stat=status)
+   !    _VERIFY(STATUS)
+   !    allocate(aglobal(im,jm), stat=status)
+   !    _VERIFY(STATUS)
+
+   !    call ArrayGather(local_array=area, global_array=aglobal, grid=grid, rc=status)
+   !    _VERIFY(STATUS)
+   !    call ArrayGather(local_array=q,    global_array=qglobal, grid=grid, rc=status)
+   !    _VERIFY(STATUS)
+   !    qdum = 0.0d+0
+   !    ! do calculation on root
+   !    if (amIRoot) then
+   !       do j=1,jm
+   !          do i=1,im
+   !             if (qglobal(i,j) == MAPL_Undef) cycle ! exclude any undefs
+   !             qdum(1) = qdum(1) + qglobal(i,j)*aglobal(i,j)
+   !             qdum(2) = qdum(2) + aglobal(i,j)
+   !          enddo
+   !       end do
+
+   !       if (qdum(2) /= 0.0d+0) then
+
+   !          qave = qdum(1) / qdum(2)
+
+   !          !ALT: convert the the result to single precision
+   !          !     This technically is not needed here, but is is done to be in sync
+   !          !     with the parallel method below
+   !          !     qave = real(qave, kind=4)
+   !       else
+   !          qave = MAPL_Undef
+   !       end if
+   !    end if
+   !    deallocate(aglobal)
+   !    deallocate(qglobal)
+
+   !    call MAPL_CommsBcast(vm, DATA=qave, N=1, Root=MAPL_Root, RC=status)
+   !    _VERIFY(STATUS)
+
+   !    _RETURN(ESMF_SUCCESS)
+   ! end subroutine AreaMean_2d_r8_bitrep
+
+   function AreaMean_2d_r8(q, area, comm, rc) result(qmean)
+      real, intent(in) :: q(:,:)
+      real, intent(in) :: area(:,:)
+      integer, intent(in) :: comm
+      integer, optional, intent(out) :: rc
+      real(real64) :: qmean ! result
+
+      real(kind=real64) :: sum_local(2), sum_global(2)
+      integer, parameter :: TWO = 2
+      integer :: im, jm, i, j, status
+
+      _ASSERT(all(shape(q) == shape(area)), "q and area need to be of the same shape")
+
+      im = size(area, 1) ! local grid dim
+      jm = size(area, 2) ! local grid dim
+
+      ! compute local sum
+      sum_local = 0.0d+0
+      do j = 1, jm
+         do i = 1, im
+            if (q(i, j) == MAPL_UNDEFINED_REAL) cycle ! exclude any undefs
+            sum_local(1) = sum_local(1) + q(i, j) * area(i, j)
+            sum_local(2) = sum_local(2) + area(i,j)
+         enddo
+      end do
+
+      call MPI_AllReduce(sum_local, sum_global, TWO, MPI_DOUBLE, MPI_SUM, comm, status)
+      _VERIFY(status)
+
+      if (sum_global(2) /= 0.0d+0) then
+         qmean = sum_global(1) / sum_global(2)
+         !ALT: convert the the result to single precision to get rid of
+         !     numerical non-associativity in floating point numbers
+         !     qmean = real(qmean, kind=4)
+      else
+         qmean = MAPL_UNDEFINED_REAL64
+      end if
+
+      _RETURN(_SUCCESS)
+   end function AreaMean_2d_r8
+
+end module mapl3g_AreaMean

utilities/arrays/CMakeLists.txt

@@ -0,0 +1,4 @@
+target_sources(MAPL.utilities PRIVATE
+  MaxMin.F90
+  AreaMean.F90
+)

utilities/arrays/MaxMin.F90

@@ -0,0 +1,140 @@
+!------------------------------------------------------------------------------
+!               Global Modeling and Assimilation Office (GMAO)                !
+!                    Goddard Earth Observing System (GEOS)                    !
+!                                 MAPL Component                              !
+!------------------------------------------------------------------------------
+!>
+!### MODULE: `mapl3g_MaxMin`
+!
+! Author: GMAO SI-Team
+!
+! `mapl3g_MaxMin` --- 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.
+
+#include "MAPL.h"
+
+module mapl3g_MaxMin
+
+   use mpi
+   use, intrinsic :: iso_fortran_env, only: real32, real64
+
+   use MAPL_ErrorHandling, only: MAPL_Verify, MAPL_Assert, MAPL_Return
+
+   implicit none
+   private
+
+   public ::  MaxMin
+
+   interface MaxMin
+      ! real32
+      module procedure pmaxmin1d_r4
+      module procedure pmaxmin2d_r4
+      module procedure pmaxmin3d_r4
+      ! real64
+      module procedure pmaxmin1d_r8
+      module procedure pmaxmin2d_r8
+      module procedure pmaxmin3d_r8
+   end interface MaxMin
+
+contains
+
+   function pmaxmin1d_r4(p, comm, rc) result(pmaxmin)
+      real, intent(in) :: p(:) ! input array
+      integer, intent(in) :: comm
+      integer, optional, intent(out) :: rc
+      real :: pmaxmin(2) !  [pmax, pmin]
+
+      integer :: im, jt, status
+
+      im = size(p)
+      jt = 1
+      pmaxmin = pmaxmin2d_r4(reshape(p, [im, jt]), comm, _RC)
+
+      _RETURN(_SUCCESS)
+   end function pmaxmin1d_r4
+
+   function pmaxmin2d_r4(p, comm, rc) result(pmaxmin)
+      real, intent(in)  :: p(:,:) ! input array
+      integer, intent(in) :: comm
+      integer, optional, intent(out) :: rc
+      real :: pmaxmin(2) ! [pmax, pmin]
+
+      real :: pmax, pmin, pm_send(2), pm_recv(2)
+      integer, parameter :: TWO=2
+      logical :: has_nans
+      integer :: status
+
+      has_nans = any(p /= p)
+      _ASSERT(.not. has_nans, "input array contains NaNs")
+
+      pm_send = [maxval(p), -minval(p)]
+      call MPI_AllReduce(pm_send, pm_recv, TWO, MPI_REAL, MPI_MAX, comm, status)
+      _VERIFY(status)
+      pmaxmin = [pm_recv(1), -pm_recv(2)]
+
+      _RETURN(_SUCCESS)
+   end function pmaxmin2d_r4
+
+   function pmaxmin3d_r4(p, comm, rc) result(pmaxmin)
+      real, intent(in) :: p(:,:,:) ! input array
+      integer, intent(in) :: comm
+      integer, optional, intent(out) :: rc
+      real :: pmaxmin(2) ! [pmax, pmin]
+
+      integer :: im, jt, status
+
+      im = size(p, 1) * size(p,2)
+      jt = size(p, 3)
+      pmaxmin = pmaxmin2d_r4(reshape(p, [im, jt]), comm, _RC)
+
+      _RETURN(_SUCCESS)
+   end function pmaxmin3d_r4
+
+   function pmaxmin1d_r8(p, comm, rc) result(pmaxmin)
+      real(real64), intent(in) :: p(:) ! input array
+      integer, intent(in) :: comm
+      integer, optional, intent(out) :: rc
+      real(real64) :: pmaxmin(2) ! [pmax, pmin]
+
+      real(real32) :: pmaxmin_r4(2)
+      integer :: status
+
+      pmaxmin_r4 =  pmaxmin1d_r4(real(p, kind=real32), comm, _RC)
+      pmaxmin = pmaxmin_r4
+
+      _RETURN(_SUCCESS)
+   end function pmaxmin1d_r8
+
+   function pmaxmin2d_r8(p, comm, rc) result(pmaxmin)
+      real(real64), intent(in) :: p(:,:) ! input array
+      integer, intent(in) :: comm
+      integer, optional, intent(out) :: rc
+      real(real64) :: pmaxmin(2) ! [pmax, pmin]
+
+      real(real32) :: pmaxmin_r4(2)
+      integer :: status
+
+      pmaxmin_r4 = pmaxmin2d_r4(real(p, kind=real32), comm, _RC)
+      pmaxmin = pmaxmin_r4
+
+      _RETURN(_SUCCESS)
+   end function pmaxmin2d_r8
+
+   function pmaxmin3d_r8(p, comm, rc) result(pmaxmin)
+      real(real64), intent(in) :: p(:,:,:) ! input array
+      integer, intent(in) :: comm
+      integer, optional, intent(out) :: rc
+      real(real64) :: pmaxmin(2) ! [pmax, pmin]
+
+      real(real32) :: pmaxmin_r4(2)
+      integer :: status
+
+      pmaxmin_r4 =  pmaxmin3d_r4(real(p, kind=real32), comm, _RC)
+      pmaxmin = pmaxmin_r4
+
+      _RETURN(_SUCCESS)
+   end function pmaxmin3d_r8
+
+end module mapl3g_MaxMin

utilities/regex/CMakeLists.txt

@@ -0,0 +1,4 @@
+target_sources(MAPL.utilities PRIVATE
+  regex_F.c
+  regex_module.F90
+)

utilities/regex/regex_F.c

@@ -0,0 +1,51 @@
+#include <sys/types.h>
+#include <regex.h>
+#include <string.h>
+#include <stdlib.h>
+
+void C_regalloc(regex_t **preg_return) {
+  *preg_return = malloc(sizeof(**preg_return));
+}
+
+/* pattern must be NUL terminated. */
+void C_regcomp(regex_t *preg, const char *pattern,
+               const char *flags, int *status_return) {
+  int i, cflags=0;
+  for (i=0;flags[i];i++) {
+    switch (flags[i]) {
+      case 'i': cflags |= REG_ICASE; break;
+      case 'm': cflags |= REG_NEWLINE; break;
+      case 'x': cflags |= REG_EXTENDED; break;
+      case 'n': cflags |= REG_NOSUB; break;
+      case ' ': break;
+      default: *status_return=-2; return;
+    }
+  }
+  *status_return = regcomp(preg,pattern,cflags);
+}
+
+void C_regexec(const regex_t *preg, const char *string, int nmatch,
+               int matches[nmatch][2], const char *flags,
+               int *status_return) {
+  int i, eflags=0;
+  regmatch_t *pmatch;
+  for (i=0;flags[i];i++) {
+    switch (flags[i]) {
+      case 'b': eflags |= REG_NOTBOL; break;
+      case 'e': eflags |= REG_NOTEOL; break;
+      case ' ': break;
+      default: *status_return=-2; return;
+    }
+  }
+  if (nmatch>0 && sizeof(pmatch->rm_so)!=sizeof(matches[0][0])) {
+    pmatch = malloc(sizeof(regmatch_t)*nmatch);
+    *status_return = regexec(preg,string,nmatch,pmatch,eflags);
+    for (i=0;i<nmatch;i++) {
+      matches[i][0]=pmatch[i].rm_so;
+      matches[i][1]=pmatch[i].rm_eo;
+    }
+    free(pmatch);
+  } else {
+    *status_return = regexec(preg,string,nmatch,(regmatch_t*)&(matches[0][0]),eflags);
+  }
+}