Add more infrastructure for water tracers

CMakeLists.txt

@@ -94,6 +94,7 @@ endif()
 
 # Among other things, this handles the genf90 generation
 add_subdirectory(src)
+add_subdirectory(src/water_isotopes)
 
 file(GLOB FSOURCES "src/*.F90" "src/water_isotopes/*.F90" "RandNum/src/*.F90" "RandNum/src/*/*.F90")
 file(GLOB CSOURCES "src/*.c" "RandNum/src/*/*.c")
@@ -110,6 +111,7 @@ target_include_directories(csm_share PRIVATE ${CMAKE_BINARY_DIR})
 if(UNITTESTS)
    # need to turn the warning check off for pfunit
    set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -Wno-error ${CMAKE_Fortran_COMPILER_DIRECTIVE} -I${CMAKE_BINARY_DIR}/unittests/shr_assert_test/mod/assert/ ")
-   add_subdirectory(${CMAKE_SOURCE_DIR}/unit_test_stubs/util csm_share_stubs)
+   add_subdirectory(${CMAKE_SOURCE_DIR}/unit_test_stubs/util csm_share_stubs_util)
+   add_subdirectory(${CMAKE_SOURCE_DIR}/unit_test_stubs/gptl csm_share_stubs_gptl)
    add_subdirectory(${CMAKE_SOURCE_DIR}/test/unit ${CMAKE_BINARY_DIR}/unittests)
 endif()

src/CMakeLists.txt

@@ -26,7 +26,8 @@ list(APPEND share_sources
   shr_mpi_mod.F90
   shr_pio_mod.F90
   shr_wv_sat_mod.F90
-  m_MergeSorts.F90)
+  m_MergeSorts.F90
+  nuopc_shr_methods.F90)
 
 sourcelist_to_parent(share_sources)
 

src/water_isotopes/CMakeLists.txt

@@ -0,0 +1,4 @@
+list(APPEND share_sources
+  shr_wtracers_mod.F90)
+
+sourcelist_to_parent(share_sources)

src/water_isotopes/shr_wtracers_mod.F90

@@ -21,6 +21,7 @@ module shr_wtracers_mod
    use shr_log_mod       , only : s_logunit=>shr_log_Unit
    use shr_string_mod    , only : shr_string_listGetAllNames, shr_string_toUpper
    use shr_string_mod    , only : shr_string_withoutSuffix
+   use shr_infnan_mod    , only : shr_infnan_isnan
    use shr_sys_mod       , only : shr_sys_abort
    use nuopc_shr_methods , only : chkerr
    use NUOPC             , only : NUOPC_CompAttributeGet
@@ -35,15 +36,24 @@ module shr_wtracers_mod
    !--------------------------------------------------------------------------
 
    public :: shr_wtracers_init              ! initialize water tracer information
+   public :: shr_wtracers_init_directly_for_testing ! initialize water tracer information directly for the sake of unit testing
    public :: shr_wtracers_finalize          ! finalize water tracer information
+   public :: shr_wtracers_initialized       ! return true if this module has been initialized
    public :: shr_wtracers_is_wtracer_field  ! return true if the given field name is a water tracer field
+   public :: shr_wtracers_get_bulk_fieldname ! return the name of the equivalent bulk field corresponding to a water tracer field
    public :: shr_wtracers_present           ! return true if there are water tracers in this simulation
    public :: shr_wtracers_get_num_tracers   ! get number of water tracers in this simulation
    public :: shr_wtracers_get_name          ! get the name of a given tracer
    public :: shr_wtracers_get_species_type  ! get the species type value associated with a given tracer
    public :: shr_wtracers_get_species_name  ! get the species name associated with a given tracer
    public :: shr_wtracers_is_isotope        ! return true if a given tracer is an isotope
    public :: shr_wtracers_get_initial_ratio ! get the initial ratio for a given tracer
+   public :: shr_wtracers_check_tracer_ratios ! check tracer ratios against expectations
+
+   interface shr_wtracers_check_tracer_ratios
+      module procedure shr_wtracers_check_tracer_ratios_1d
+      module procedure shr_wtracers_check_tracer_ratios_2d
+   end interface shr_wtracers_check_tracer_ratios
 
    !--------------------------------------------------------------------------
    ! Private interfaces
@@ -127,6 +137,66 @@ subroutine shr_wtracers_init(driver, maintask, rc)
 
    end subroutine shr_wtracers_init
 
+   !-----------------------------------------------------------------------
+   subroutine shr_wtracers_init_directly_for_testing( &
+              water_tracer_names, water_tracer_species, water_tracer_initial_ratios, rc)
+      !
+      ! !DESCRIPTION:
+      ! Initialize water tracer information directly for the sake of unit testing
+      !
+      ! If there are any errors, an ESMF error code is returned in rc
+      !
+      ! !ARGUMENTS
+      character(len=*), intent(in) :: water_tracer_names(:)
+      character(len=*), intent(in) :: water_tracer_species(:)  ! expected to be uppercase
+      real(r8), intent(in) :: water_tracer_initial_ratios(:)
+      integer, intent(out) :: rc
+      !
+      ! !LOCAL VARIABLES
+      character(len=*), parameter :: subname='shr_wtracers_init_directly_for_testing'
+      !---------------------------------------------------------------
+
+      rc = ESMF_SUCCESS
+
+      if (water_tracers_initialized) then
+         call shr_log_error("Attempt to call "//subname//" multiple times", rc=rc)
+         return
+      end if
+
+      num_tracers = size(water_tracer_names)
+      if (size(water_tracer_species) /= num_tracers .or. &
+          size(water_tracer_initial_ratios) /= num_tracers) then
+         call shr_log_error(subname//": Input array sizes disagree", rc=rc)
+         return
+      end if
+
+      tracer_names = water_tracer_names
+      tracer_species_names = water_tracer_species
+      call shr_wtracers_set_species_types(rc=rc)
+      if (chkerr(rc,__LINE__,u_FILE_u)) return
+      tracer_initial_ratios = water_tracer_initial_ratios
+
+      water_tracers_initialized = .true.
+
+   end subroutine shr_wtracers_init_directly_for_testing
+
+   !-----------------------------------------------------------------------
+   function shr_wtracers_initialized()
+      !
+      ! !DESCRIPTION:
+      ! Return true if this module has been initialized
+      !
+      ! !ARGUMENTS
+      logical :: shr_wtracers_initialized  ! function result
+      !
+      ! !LOCAL VARIABLES:
+      character(len=*), parameter :: subname='shr_wtracers_initialized'
+      !-----------------------------------------------------------------------
+
+      shr_wtracers_initialized = water_tracers_initialized
+
+   end function shr_wtracers_initialized
+
    !-----------------------------------------------------------------------
    subroutine shr_wtracers_parse_attributes(driver, rc)
       !
@@ -388,6 +458,45 @@ function shr_wtracers_is_wtracer_field(fieldname)
       shr_wtracers_is_wtracer_field = is_tracer
    end function shr_wtracers_is_wtracer_field
 
+   !-----------------------------------------------------------------------
+   subroutine shr_wtracers_get_bulk_fieldname(fieldname, is_wtracer_field, bulk_fieldname)
+      !
+      ! !DESCRIPTION:
+      ! Return the name of the equivalent bulk field corresponding to a water tracer field
+      !
+      ! If fieldname is the name of a water tracer field, based on naming conventions,
+      ! then is_wtracer_field will be true and bulk_fieldname will hold the name of the
+      ! corresponding bulk field.
+      !
+      ! If fieldname is *not* the name of a water tracer field, then is_wtracer_field will
+      ! be false, and bulk_fieldname will be the same as fieldname.
+      !
+      ! Note that, unlike most other routines in this module, this function works even if
+      ! the data in this module has not been initialized (i.e., even if shr_wtracers_init
+      ! has not been called): it works simply based on naming conventions.
+      !
+      ! !ARGUMENTS
+      character(len=*), intent(in)  :: fieldname
+      logical         , intent(out) :: is_wtracer_field
+      character(len=*), intent(out) :: bulk_fieldname
+      !
+      ! !LOCAL VARIABLES:
+      integer :: localrc
+
+      character(len=*), parameter :: subname='shr_wtracers_get_bulk_fieldname'
+      !-----------------------------------------------------------------------
+
+      call shr_string_withoutSuffix( &
+           in_str = fieldname, &
+           suffix = WTRACERS_SUFFIX, &
+           has_suffix = is_wtracer_field, &
+           out_str = bulk_fieldname, &
+           rc = localrc)
+      if (localrc /= 0) then
+         call shr_sys_abort(subname//": ERROR in shr_string_withoutSuffix")
+      end if
+   end subroutine shr_wtracers_get_bulk_fieldname
+
    !-----------------------------------------------------------------------
    function shr_wtracers_present()
       !
@@ -531,6 +640,146 @@ function shr_wtracers_get_initial_ratio(tracer_num)
       shr_wtracers_get_initial_ratio = tracer_initial_ratios(tracer_num)
    end function shr_wtracers_get_initial_ratio
 
+   !-----------------------------------------------------------------------
+   subroutine shr_wtracers_check_tracer_ratios_1d(tracers, bulk, name, extra_dim_index)
+      !
+      ! !DESCRIPTION:
+      ! Check tracer ratios (tracer/bulk) against expectations
+      !
+      ! Aborts if any inconsistencies are found
+      !
+      ! Should only be called in simulations set up to maintain constant water tracer
+      ! ratios: in general, water tracers will deviate from their initial, fixed ratios,
+      ! and so it makes no sense to perform these checks since they will always fail.
+      !
+      ! !ARGUMENTS
+      real(r8), intent(in) :: tracers(:,:)  ! dimensioned [tracerNum, gridcell]
+      real(r8), intent(in) :: bulk(:)
+      character(len=*), intent(in) :: name  ! for diagnostic output
+      integer, intent(in), optional :: extra_dim_index  ! index of extra dimension (for error messages)
+      !
+      ! !LOCAL VARIABLES
+      integer :: n, i
+      logical :: arrays_equal
+      integer :: diff_tracer, diff_loc
+      real(r8) :: val_bulk, val_tracer
+
+      real(r8), parameter :: tolerance = 1.0e-7_r8
+
+      character(len=*), parameter :: subname='shr_wtracers_check_tracer_ratios_1d'
+      ! In some error messages, it makes more sense to print the generic name:
+      character(len=*), parameter :: subname_generic='shr_wtracers_check_tracer_ratios'
+      !-----------------------------------------------------------------------
+      if (.not. water_tracers_initialized) then
+         call shr_sys_abort(subname//" ERROR: water tracers not yet initialized")
+      end if
+      if (size(tracers, 1) /= num_tracers) then
+         call shr_sys_abort(subname//" ERROR: unexpected number of tracers")
+      end if
+      if (size(tracers, 2) /= size(bulk)) then
+         call shr_sys_abort(subname//" ERROR: inconsistent sizes for tracers and bulk")
+      end if
+
+      arrays_equal = .true.
+      tracer_loop: do n = 1, num_tracers
+         ! We may eventually want a mechanism to denote certain tracers as for-checking
+         ! and others not-for-checking (probably via another nuopc attribute parallel to
+         ! the other water tracers attributes). If we do, we could check that flag here
+         ! for tracer #n and go on to the next tracer if this is a not-for-checking
+         ! tracer.
+
+         do i = 1, size(bulk)
+            if (.not. shr_infnan_isnan(bulk(i)) .and. .not. shr_infnan_isnan(tracers(n,i))) then
+               ! neither value is nan: check error tolerance
+               val_bulk = bulk(i) * tracer_initial_ratios(n)
+               val_tracer = tracers(n,i)
+               if (val_bulk == 0.0_r8 .and. val_tracer == 0.0_r8) then
+                  ! trap special case where both are zero to avoid division by zero: values equal (do nothing)
+               else if (abs(val_bulk - val_tracer) / max(abs(val_bulk), abs(val_tracer)) > tolerance) then
+                  arrays_equal = .false.
+                  diff_tracer = n
+                  diff_loc = i
+                  exit tracer_loop
+               else
+                  ! error < tolerance: values considered equal (do nothing)
+               end if
+            else if (shr_infnan_isnan(bulk(i)) .and. shr_infnan_isnan(tracers(n,i))) then
+               ! both values are nan: values are considered equal (do nothing)
+            else
+               ! only one value is nan: not equal
+               arrays_equal = .false.
+               diff_tracer = n
+               diff_loc = i
+               exit tracer_loop
+            end if
+         end do
+      end do tracer_loop
+
+      if (.not. arrays_equal) then
+         write(s_logunit, '(A,A)') subname_generic, " ERROR: tracer does not agree with bulk water"
+         write(s_logunit, '(A,A)') "Variable: ", trim(name)
+         if (present(extra_dim_index)) then
+            write(s_logunit, '(A,I0)') "Extra dimension index: ", extra_dim_index
+         end if
+         write(s_logunit, '(A,I0,A,A)') "First difference found for tracer #", diff_tracer, &
+              ": ", tracer_names(diff_tracer)
+         write(s_logunit, '(A,I0)') "First difference at index: ", diff_loc
+         write(s_logunit, '(A, ES25.17)') "Bulk  : ", bulk(diff_loc)
+         write(s_logunit, '(A, ES25.17)') "Tracer: ", tracers(diff_tracer, diff_loc)
+         write(s_logunit, '(A, ES25.17)') "Expected ratio: ", tracer_initial_ratios(diff_tracer)
+         if (.not. shr_infnan_isnan(bulk(diff_loc))) then
+            write(s_logunit, '(A, ES25.17)') "Bulk*ratio: ", bulk(diff_loc) * tracer_initial_ratios(diff_tracer)
+         end if
+         call shr_sys_abort(subname_generic//" ERROR: tracer does not agree with bulk water")
+      end if
+
+   end subroutine shr_wtracers_check_tracer_ratios_1d
+
+   !-----------------------------------------------------------------------
+   subroutine shr_wtracers_check_tracer_ratios_2d(tracers, bulk, name)
+      !
+      ! !DESCRIPTION:
+      ! Check tracer ratios (tracer/bulk) against expectations for 2-d bulk arrays
+      !
+      ! This is a lightweight wrapper around shr_wtracers_check_tracer_ratios_1d that
+      ! loops over the first dimension of the bulk array.
+      !
+      ! Aborts if any inconsistencies are found
+      !
+      ! Should only be called in simulations set up to maintain constant water tracer
+      ! ratios: in general, water tracers will deviate from their initial, fixed ratios,
+      ! and so it makes no sense to perform these checks since they will always fail.
+      !
+      ! !ARGUMENTS
+      real(r8), intent(in) :: tracers(:,:,:)  ! dimensioned [ungriddedDim, tracerNum, gridcell]
+      real(r8), intent(in) :: bulk(:,:)       ! dimensioned [ungriddedDim, gridcell]
+      character(len=*), intent(in) :: name    ! for diagnostic output
+      !
+      ! !LOCAL VARIABLES
+      integer :: i
+
+      character(len=*), parameter :: subname='shr_wtracers_check_tracer_ratios_2d'
+      !-----------------------------------------------------------------------
+      if (.not. water_tracers_initialized) then
+         call shr_sys_abort(subname//" ERROR: water tracers not yet initialized")
+      end if
+      if (size(tracers, 1) /= size(bulk, 1)) then
+         call shr_sys_abort(subname//" ERROR: inconsistent size for tracers dim 1 and bulk dim 1")
+      end if
+      if (size(tracers, 2) /= num_tracers) then
+         call shr_sys_abort(subname//" ERROR: unexpected number of tracers")
+      end if
+      if (size(tracers, 3) /= size(bulk, 2)) then
+         call shr_sys_abort(subname//" ERROR: inconsistent size for tracers dim 3 and bulk dim 2")
+      end if
+
+      do i = 1, size(bulk, 1)
+         call shr_wtracers_check_tracer_ratios_1d(tracers(i,:,:), bulk(i,:), name, &
+              extra_dim_index=i)
+      end do
+
+   end subroutine shr_wtracers_check_tracer_ratios_2d
+
    !-----------------------------------------------------------------------
    subroutine shr_wtracers_finalize(rc)
       !

test/unit/CMakeLists.txt

@@ -18,8 +18,10 @@ add_subdirectory(dynamic_vector)
 
 add_subdirectory(shr_vmath_test)
 
+add_subdirectory(shr_wtracers_test)
+
 add_subdirectory(shr_wv_sat_test)
 
 add_subdirectory(shr_precip_test)
 
-add_subdirectory(shr_cal_test)
+add_subdirectory(shr_cal_test)

test/unit/shr_wtracers_test/CMakeLists.txt

@@ -0,0 +1,29 @@
+set (pf_sources
+  test_shr_wtracers.pf
+  )
+
+set(sources_needed
+  shr_wtracers_mod.F90
+  shr_assert_mod.F90
+  shr_infnan_mod.F90
+  shr_kind_mod.F90
+  shr_log_mod.F90
+  shr_strconvert_mod.F90
+  shr_string_mod.F90
+  shr_sys_mod.nompi_abortthrows.F90
+  shr_abort_mod.abortthrows.F90
+  shr_timer_mod.F90
+  nuopc_shr_methods.F90
+  gptl.F90)
+
+extract_sources("${sources_needed}" "${share_sources}" test_sources)
+
+add_pfunit_ctest(shr_wtracers
+  TEST_SOURCES "${pf_sources}"
+  OTHER_SOURCES "${test_sources}")
+
+declare_generated_dependencies(shr_wtracers "${share_genf90_sources}")
+
+target_link_libraries(shr_wtracers esmf)
+# The following adds all dependencies of ESMF, including PIO, NetCDF, etc.:
+target_link_libraries(shr_wtracers ESMF::ESMF)