MODERN_DIAG_MANAGER: fix error when setting zbounds for more than one variable

CMakeLists.txt

@@ -663,6 +663,7 @@ list(APPEND TEST_SRC_SINGLE_TARGET
   test_fms/diag_manager/test_diag_attribute_add.F90
   test_fms/diag_manager/check_new_file_freq.F90
   test_fms/diag_manager/test_zbounds_limits.F90
+  test_fms/diag_manager/test_multiple_zbounds.F90
   test_fms/drifters/test_cloud_interpolator.F90
   test_fms/drifters/test_drifters_io.F90
   test_fms/drifters/test_drifters_input.F90

diag_manager/fms_diag_axis_object.F90

@@ -40,7 +40,7 @@ module fms_diag_axis_object_mod
   use mpp_mod,         only:  FATAL, mpp_error, uppercase, mpp_pe, mpp_root_pe, stdout
   use fms2_io_mod,     only:  FmsNetcdfFile_t, FmsNetcdfDomainFile_t, FmsNetcdfUnstructuredDomainFile_t, &
                             & register_axis, register_field, register_variable_attribute, write_data
-  use fms_diag_yaml_mod, only: subRegion_type, diag_yaml, MAX_SUBAXES
+  use fms_diag_yaml_mod, only: subRegion_type, diag_yaml, MAX_SUBAXES, diagYamlFilesVar_type
   use diag_grid_mod,       only:  get_local_indices_cubesphere => get_local_indexes
   use axis_utils2_mod,   only: nearest_index
   implicit none
@@ -52,7 +52,7 @@ module fms_diag_axis_object_mod
           & DIAGDOMAIN2D_T, fmsDiagSubAxis_type, fmsDiagAxisContainer_type, fmsDiagFullAxis_type, DIAGDOMAINUG_T
   public :: define_new_axis, parse_compress_att, get_axis_id_from_name, define_diurnal_axis, &
           & fmsDiagDiurnalAxis_type, create_new_z_subaxis, is_parent_axis, define_new_subaxis_latlon, &
-          & define_new_subaxis_index
+          & define_new_subaxis_index, find_z_sub_axis_name
 
   !> @}
 
@@ -124,6 +124,7 @@ module fms_diag_axis_object_mod
       procedure :: get_starting_index
       procedure :: get_ending_index
       procedure :: get_compute_indices
+      procedure :: is_same_zbounds
   END TYPE fmsDiagSubAxis_type
 
   !> @brief Type to hold the diurnal axis
@@ -934,6 +935,19 @@ function get_compute_indices(this) result(indx)
     indx = this%compute_idx
   end function get_compute_indices
 
+  !> @brief Determines if the zbounds passed in are the same as those in the file
+  !! @return .True. if the zbounds are the same
+  function is_same_zbounds(this, zbounds) result(is_same)
+    class(fmsDiagSubAxis_type), intent(in) :: this       !< diag_sub_axis object
+    real(kind=r4_kind),         intent(in) :: zbounds(2) !< zbounds to compare with
+    logical :: is_same
+
+    is_same = .false.
+    if (zbounds(1) .eq. this%zbounds(1) .and. zbounds(2) .eq. this%zbounds(2)) then
+      is_same = .true.
+    endif
+  end function
+
   !> @brief Get the ntiles in a domain
   !> @return the number of tiles in a domain
   function get_ntiles(this) &
@@ -1451,52 +1465,57 @@ subroutine create_new_z_subaxis(zbounds, var_axis_ids, diag_axis, naxis, file_ax
                                             !! axis
     integer           :: i                  !< For do loops
     integer           :: subaxis_id         !< The id of the new z subaxis
-    logical           :: axis_found         !< Flag that indicated if the zsubaxis already exists
+    integer           :: parent_axis_id     !< Id of parent axis id
+    integer           :: zaxis_index        !< Index of the z axis (i.e 3 if the variable is x,y,z)
+    type(fmsDiagFullAxis_type), pointer :: parent_axis !< Pointer to the parent axis
+
+    parent_axis_id = diag_null
+    zaxis_index = diag_null
+
+    !< Determine which axis is the z axis:
+    do i = 1, size(var_axis_ids)
+      select type (parent_axis => diag_axis(var_axis_ids(i))%axis)
+      type is (fmsDiagFullAxis_type)
+        if (parent_axis%cart_name .eq. "Z") then
+          parent_axis_id = var_axis_ids(i)
+          zaxis_index = i
+        endif
+      end select
+    enddo
 
     !< Determine if the axis was already created
-    axis_found = .false.
     do i = 1, nfile_axis
       select type (axis => diag_axis(file_axis_id(i))%axis)
       type is (fmsDiagSubAxis_type)
+        if (axis%parent_axis_id .ne. parent_axis_id) cycle
         if (axis%zbounds(1) .eq. zbounds(1) .and. axis%zbounds(2) .eq. zbounds(2)) then
-          axis_found = .true.
-          subaxis_id = file_axis_id(i)
-          exit
+          var_axis_ids(zaxis_index) = file_axis_id(i)
+          return
         endif
       end select
     enddo
 
-    !< Determine which of the variable's axis is the zaxis!
-    do i = 1, size(var_axis_ids)
-      select type (parent_axis => diag_axis(var_axis_ids(i))%axis)
-      type is (fmsDiagFullAxis_type)
-        if (parent_axis%cart_name .eq. "Z") then
-          !< If the axis was previously defined set the var_axis_ids and leave
-          if (axis_found) then
-            var_axis_ids(i) = subaxis_id
-            return
-          endif
-          zaxis_data => parent_axis%axis_data
-
-          select type(zaxis_data)
-          type is (real(kind=r4_kind))
-            !TODO need to include the conversion to "real" because nearest_index doesn't take r4s and r8s
-            subaxis_indices(1) = nearest_index(real(zbounds(1)), real(zaxis_data))
-            subaxis_indices(2) = nearest_index(real(zbounds(2)), real(zaxis_data))
-          type is (real(kind=r8_kind))
-            subaxis_indices(1) = nearest_index(real(zbounds(1)), real(zaxis_data))
-            subaxis_indices(2) = nearest_index(real(zbounds(2)), real(zaxis_data))
-          end select
-
-          nz_subaxis = nz_subaxis + 1
-          call define_new_axis(diag_axis, parent_axis, naxis, parent_axis%axis_id, &
+    select type (axis => diag_axis(parent_axis_id)%axis)
+    type is (fmsDiagFullAxis_type)
+      zaxis_data => axis%axis_data
+      parent_axis => axis
+    end select
+
+    select type(zaxis_data)
+      type is (real(kind=r4_kind))
+        !TODO need to include the conversion to "real" because nearest_index doesn't take r4s and r8s
+        subaxis_indices(1) = nearest_index(real(zbounds(1)), real(zaxis_data))
+        subaxis_indices(2) = nearest_index(real(zbounds(2)), real(zaxis_data))
+      type is (real(kind=r8_kind))
+        subaxis_indices(1) = nearest_index(real(zbounds(1)), real(zaxis_data))
+        subaxis_indices(2) = nearest_index(real(zbounds(2)), real(zaxis_data))
+    end select
+
+    nz_subaxis = nz_subaxis + 1
+    call define_new_axis(diag_axis, parent_axis, naxis, parent_axis%axis_id, &
                         &subaxis_indices(1), subaxis_indices(2), (/lbound(zaxis_data,1), ubound(zaxis_data,1)/), &
                         &new_axis_id=subaxis_id, zbounds=zbounds, nz_subaxis=nz_subaxis)
-          var_axis_ids(i) = subaxis_id
-          return
-        endif
-      end select
-    enddo
+    var_axis_ids(zaxis_index) = subaxis_id
 
   end subroutine
 
@@ -1517,6 +1536,32 @@ function is_parent_axis(axis_id, parent_axis_id, diag_axis) &
     end select
   end function is_parent_axis
 
+  !> @brief Determine the name of the z subaxis by matching the parent axis id and the zbounds
+  !! in the diag table yaml
+  subroutine find_z_sub_axis_name(dim_name, parent_axis_id, file_axis_id, field_yaml, diag_axis)
+    character(len=*),                intent(inout) :: dim_name         !< Name of z subaxis
+    integer,                         intent(in)    :: parent_axis_id   !< Axis id of the parent
+    integer,                         intent(in)    :: file_axis_id(:)  !< Axis ids of the file
+    type(diagYamlFilesVar_type),     intent(in)    :: field_yaml       !< Field info from diag_table yaml
+    class(fmsDiagAxisContainer_type),intent(in)    :: diag_axis(:)     !< Array of axis objections
+
+    integer :: id
+    integer :: i
+
+    do i = 1, size(file_axis_id)
+      id = file_axis_id(i)
+      select type (axis_ptr => diag_axis(id)%axis)
+      type is (fmsDiagSubAxis_type)
+        if (axis_ptr%parent_axis_id .eq. parent_axis_id) then
+          if (axis_ptr%is_same_zbounds(field_yaml%get_var_zbounds())) then
+            dim_name = axis_ptr%subaxis_name
+            return
+          endif
+        endif
+      end select
+    enddo
+    dim_name = ""
+  end subroutine
 #endif
 end module fms_diag_axis_object_mod
 !> @}

diag_manager/fms_diag_field_object.F90

@@ -36,7 +36,7 @@ module fms_diag_field_object_mod
 use fms_diag_yaml_mod, only:  diagYamlFilesVar_type, get_diag_fields_entries, get_diag_files_id, &
   & find_diag_field, get_num_unique_fields, diag_yaml
 use fms_diag_axis_object_mod, only: diagDomain_t, get_domain_and_domain_type, fmsDiagAxis_type, &
-  & fmsDiagAxisContainer_type, fmsDiagFullAxis_Type
+  & fmsDiagAxisContainer_type, fmsDiagFullAxis_Type, find_z_sub_axis_name
 use time_manager_mod, ONLY: time_type, get_date
 use fms2_io_mod, only: FmsNetcdfFile_t, FmsNetcdfDomainFile_t, FmsNetcdfUnstructuredDomainFile_t, register_field, &
                        register_variable_attribute
@@ -1162,14 +1162,16 @@ pure function get_longname_to_write(this, field_yaml) &
 end function get_longname_to_write
 
 !> @brief Determine the dimension names to use when registering the field to fms2_io
-subroutine get_dimnames(this, diag_axis, field_yaml, unlim_dimname, dimnames, is_regional)
+subroutine get_dimnames(this, diag_axis, field_yaml, unlim_dimname, dimnames, is_regional, &
+                        file_axis_ids)
   class (fmsDiagField_type),        target, intent(inout) :: this          !< diag field
   class(fmsDiagAxisContainer_type), target, intent(in)    :: diag_axis(:)  !< Diag_axis object
   type(diagYamlFilesVar_type),              intent(in)    :: field_yaml    !< Field info from diag_table yaml
   character(len=*),                         intent(in)    :: unlim_dimname !< The name of unlimited dimension
   character(len=120), allocatable,          intent(out)   :: dimnames(:)   !< Array of the dimension names
                                                                            !! for the field
   logical,                                  intent(in)    :: is_regional   !< Flag indicating if the field is regional
+  integer,                                  intent(in)    :: file_axis_ids(:) !< Ids of the file axis
 
   integer                                   :: i     !< For do loops
   integer                                   :: naxis !< Number of axis for the field
@@ -1193,7 +1195,7 @@ subroutine get_dimnames(this, diag_axis, field_yaml, unlim_dimname, dimnames, is
     do i = 1, size(this%axis_ids)
       axis_ptr => diag_axis(this%axis_ids(i))
       if (axis_ptr%axis%is_z_axis()) then
-        dimnames(i) = axis_ptr%axis%get_axis_name(is_regional)//"_sub01"
+        call find_z_sub_axis_name(dimnames(i), this%axis_ids(i), file_axis_ids, field_yaml, diag_axis)
       else
         dimnames(i) = axis_ptr%axis%get_axis_name(is_regional)
       endif
@@ -1238,7 +1240,7 @@ end subroutine register_field_wrap
 
 !> @brief Write the field's metadata to the file
 subroutine write_field_metadata(this, fms2io_fileobj, file_id, yaml_id, diag_axis, unlim_dimname, is_regional, &
-                                cell_measures, use_collective_writes)
+                                cell_measures, use_collective_writes, file_axis_ids)
   class (fmsDiagField_type), target, intent(inout) :: this          !< diag field
   class(FmsNetcdfFile_t),            INTENT(INOUT) :: fms2io_fileobj!< Fms2_io fileobj to write to
   integer,                           intent(in)    :: file_id       !< File id of the file to write to
@@ -1249,6 +1251,7 @@ subroutine write_field_metadata(this, fms2io_fileobj, file_id, yaml_id, diag_axi
   character(len=*),                  intent(in)    :: cell_measures !< The cell measures attribute to write
   logical,                           intent(in)    :: use_collective_writes !< True if using collective writes
                                                                             !! for this variable
+  integer,                           intent(in)    :: file_axis_ids(:) !< Ids of all of the axis in thje file
 
   type(diagYamlFilesVar_type), pointer     :: field_yaml  !< pointer to the yaml entry
   character(len=:),            allocatable :: var_name    !< Variable name
@@ -1267,7 +1270,7 @@ subroutine write_field_metadata(this, fms2io_fileobj, file_id, yaml_id, diag_axi
   var_name = field_yaml%get_var_outname()
 
   if (allocated(this%axis_ids)) then
-    call this%get_dimnames(diag_axis, field_yaml, unlim_dimname, dimnames, is_regional)
+    call this%get_dimnames(diag_axis, field_yaml, unlim_dimname, dimnames, is_regional, file_axis_ids)
 
     !! Collective writes are only used for 2D+ variables
     if ((use_collective_writes .and. size(this%axis_ids) >= 2) .or. field_yaml%has_chunksizes()) then

diag_manager/fms_diag_file_object.F90

@@ -1851,7 +1851,7 @@ subroutine write_field_metadata(this, diag_field, diag_axis)
 
     call field_ptr%write_field_metadata(fms2io_fileobj, diag_file%id, diag_file%yaml_ids(i), diag_axis, &
       this%FMS_diag_file%get_file_unlimdim(), is_regional, cell_measures, &
-      diag_file%is_using_collective_writes())
+      diag_file%is_using_collective_writes(), diag_file%axis_ids(1:diag_file%number_of_axis))
   enddo
 
   if (need_associated_files) &

test_fms/diag_manager/Makefile.am

@@ -35,7 +35,7 @@ check_PROGRAMS = test_diag_manager test_diag_manager_time \
  check_time_pow check_time_rms check_subregional test_cell_measures test_var_masks \
  check_var_masks test_multiple_send_data test_diag_out_yaml test_output_every_freq \
  test_dm_weights test_prepend_date test_ens_runs test_diag_multi_file test_diag_attribute_add \
- check_new_file_freq test_zbounds_limits
+ check_new_file_freq test_zbounds_limits test_multiple_zbounds
 
 # This is the source code for the test.
 test_output_every_freq_SOURCES = test_output_every_freq.F90
@@ -71,6 +71,7 @@ test_ens_runs_SOURCES = test_ens_runs.F90
 test_diag_attribute_add_SOURCES = test_diag_attribute_add.F90
 check_new_file_freq_SOURCES = check_new_file_freq.F90
 test_zbounds_limits_SOURCES = test_zbounds_limits.F90
+test_multiple_zbounds_SOURCES = test_multiple_zbounds.F90
 
 TEST_EXTENSIONS = .sh
 SH_LOG_DRIVER = env AM_TAP_AWK='$(AWK)' $(SHELL) \
@@ -81,7 +82,7 @@ TESTS = test_diag_manager2.sh test_time_none.sh test_time_min.sh test_time_max.s
         test_time_avg.sh test_time_pow.sh test_time_rms.sh test_time_diurnal.sh test_cell_measures.sh \
         test_subregional.sh test_var_masks.sh test_multiple_send_data.sh test_output_every_freq.sh \
         test_dm_weights.sh test_flush_nc_file.sh test_prepend_date.sh test_ens_runs.sh test_diag_multi_file.sh \
-        test_diag_attribute_add.sh test_time_none_netcdf4.sh test_zbounds_limits.sh
+        test_diag_attribute_add.sh test_time_none_netcdf4.sh test_zbounds_limits.sh test_multiple_zbounds.sh
 
 testing_utils.mod: testing_utils.$(OBJEXT)
 
@@ -91,7 +92,7 @@ EXTRA_DIST = test_diag_manager2.sh check_crashes.sh test_time_none.sh test_time_
              test_cell_measures.sh test_subregional.sh test_var_masks.sh test_multiple_send_data.sh \
              test_flush_nc_file.sh test_dm_weights.sh test_output_every_freq.sh test_prepend_date.sh \
              test_ens_runs.sh test_diag_multi_file.sh test_diag_attribute_add.sh test_time_none_netcdf4.sh \
-             test_zbounds_limits.sh
+             test_zbounds_limits.sh test_multiple_zbounds.sh
 
 if USING_YAML
   parser_skip=""