feat: modern diag subaxis setup (NOAA-GFDL#1056)

diag_manager/Makefile.am

@@ -68,9 +68,11 @@ fms_diag_object_mod.$(FC_MODEXT): diag_data_mod.$(FC_MODEXT) fms_diag_file_objec
                                   fms_diag_time_utils_mod.$(FC_MODEXT) \
                                   fms_diag_buffer_mod.$(FC_MODEXT)
 fms_diag_field_object_mod.$(FC_MODEXT): diag_data_mod.$(FC_MODEXT) fms_diag_yaml_mod.$(FC_MODEXT) fms_diag_time_utils_mod.$(FC_MODEXT)
-fms_diag_file_object_mod.$(FC_MODEXT): diag_data_mod.$(FC_MODEXT) fms_diag_yaml_mod.$(FC_MODEXT) fms_diag_field_object_mod.$(FC_MODEXT) fms_diag_time_utils_mod.$(FC_MODEXT)
+fms_diag_file_object_mod.$(FC_MODEXT): diag_data_mod.$(FC_MODEXT) fms_diag_yaml_mod.$(FC_MODEXT) fms_diag_field_object_mod.$(FC_MODEXT) fms_diag_time_utils_mod.$(FC_MODEXT) \
+                                       fms_diag_axis_object_mod.$(FC_MODEXT)
 fms_diag_object_container_mod.$(FC_MODEXT): fms_diag_object_mod.$(FC_MODEXT) fms_diag_dlinked_list_mod.$(FC_MODEXT) fms_diag_time_utils_mod.$(FC_MODEXT)
-fms_diag_axis_object_mod.$(FC_MODEXT): diag_data_mod.$(FC_MODEXT) fms_diag_time_utils_mod.$(FC_MODEXT)
+fms_diag_axis_object_mod.$(FC_MODEXT): diag_data_mod.$(FC_MODEXT) fms_diag_time_utils_mod.$(FC_MODEXT) fms_diag_yaml_mod.$(FC_MODEXT) \
+                                       diag_grid_mod.$(FC_MODEXT)
 diag_manager_mod.$(FC_MODEXT): diag_axis_mod.$(FC_MODEXT) diag_data_mod.$(FC_MODEXT) diag_util_mod.$(FC_MODEXT) \
                                diag_output_mod.$(FC_MODEXT) diag_grid_mod.$(FC_MODEXT) diag_table_mod.$(FC_MODEXT) \
                                fms_diag_object_mod.$(FC_MODEXT) fms_diag_file_object_mod.$(FC_MODEXT) \

diag_manager/fms_diag_buffer.F90

@@ -27,7 +27,6 @@ module fms_diag_buffer_mod
 
 use platform_mod
 use iso_c_binding
-use fms_diag_axis_object_mod, only: diagDomain_t
 use time_manager_mod, only: time_type
 use mpp_mod, only: mpp_error, FATAL
 use diag_data_mod, only: DIAG_NULL, DIAG_NOT_REGISTERED, i4, i8, r4, r8

diag_manager/fms_diag_field_object.F90

@@ -330,6 +330,8 @@ subroutine set_diag_id(this , id)
  if (allocated(this%registered)) then
      if (this%registered) then
           call mpp_error("set_diag_id", "The variable"//this%varname//" is already registered", FATAL)
+     else
+       this%diag_id = id
      endif
  else
      this%diag_id = id

diag_manager/fms_diag_file_object.F90

@@ -29,13 +29,15 @@ module fms_diag_file_object_mod
                        get_instance_filename, open_file, close_file, get_mosaic_tile_file
 use diag_data_mod, only: DIAG_NULL, NO_DOMAIN, max_axes, SUB_REGIONAL, get_base_time, DIAG_NOT_REGISTERED, &
                          TWO_D_DOMAIN, UG_DOMAIN, prepend_date, DIAG_DAYS, VERY_LARGE_FILE_FREQ
-use time_manager_mod, only: time_type, operator(>), operator(/=), operator(==), get_date
+use time_manager_mod, only: time_type, operator(>), operator(/=), operator(==), get_date, &
+                            date_to_string
 use fms_diag_time_utils_mod, only: diag_time_inc, get_time_string
-use time_manager_mod, only: time_type, operator(/=), operator(==), date_to_string
-use fms_diag_yaml_mod, only: diag_yaml, diagYamlObject_type, diagYamlFiles_type
+use fms_diag_yaml_mod, only: diag_yaml, diagYamlObject_type, diagYamlFiles_type, subRegion_type
 use fms_diag_axis_object_mod, only: diagDomain_t, get_domain_and_domain_type, fmsDiagAxis_type, &
-                                    fmsDiagAxisContainer_type, DIAGDOMAIN2D_T, DIAGDOMAINUG_T
+                                    fmsDiagAxisContainer_type, DIAGDOMAIN2D_T, DIAGDOMAINUG_T, &
+                                    fmsDiagFullAxis_type, define_subaxis
 use mpp_mod, only: mpp_get_current_pelist, mpp_npes, mpp_root_pe, mpp_pe, mpp_error, FATAL, stdout
+
 implicit none
 private
 
@@ -95,8 +97,7 @@ module fms_diag_file_object_mod
  procedure, public :: get_file_freq
  procedure, public :: get_file_timeunit
  procedure, public :: get_file_unlimdim
-!! TODO get functions for sub region stuff
-! procedure, public :: get_file_sub_region
+ procedure, public :: get_file_sub_region
  procedure, public :: get_file_new_file_freq
  procedure, public :: get_file_new_file_freq_units
  procedure, public :: get_file_start_time
@@ -123,6 +124,7 @@ module fms_diag_file_object_mod
 type, extends (fmsDiagFile_type) :: subRegionalFile_type
   integer, dimension(:), allocatable :: sub_axis_ids !< Array of axis ids in the file
   logical :: write_on_this_pe !< Flag indicating if the subregion is on the current PE
+  logical :: is_subaxis_defined !< Flag indicating if the subaxes have already been defined
 end type subRegionalFile_type
 
 !> \brief A container for fmsDiagFile_type.  This is used to create the array of files
@@ -131,6 +133,8 @@ module fms_diag_file_object_mod
 
   contains
   procedure :: open_diag_file
+  procedure :: write_axis_metadata
+  procedure :: write_axis_data
 end type fmsDiagFileContainer_type
 
 !type(fmsDiagFile_type), dimension (:), allocatable, target :: FMS_diag_file !< The array of diag files
@@ -158,7 +162,9 @@ logical function fms_diag_files_object_init (files_array)
          type is (subRegionalFile_type)
            allocate(obj%sub_axis_ids(max_axes))
            obj%sub_axis_ids = diag_null
-           obj%write_on_this_pe = .true. !TODO this should be .false. probably
+           obj%write_on_this_pe = .false.
+           obj%is_subaxis_defined = .false.
+           obj%number_of_axis = 0
        end select
      else
        allocate(FmsDiagFile_type::files_array(i)%FMS_diag_file)
@@ -326,14 +332,13 @@ pure function get_file_unlimdim (this) result(res)
   res = this%diag_yaml_file%get_file_unlimdim()
 end function get_file_unlimdim
 
-!! TODO - get functions for sub region stuff
 !> \brief Returns a copy of file_sub_region from the yaml object
 !! \return Copy of file_sub_region
-!pure function get_file_sub_region (obj) result(res)
-! class(fmsDiagFile_type), intent(in) :: obj !< The file object
-! integer :: res
-!  res = obj%diag_yaml_file%get_file_sub_region()
-!end function get_file_sub_region
+function get_file_sub_region (obj) result(res)
+ class(fmsDiagFile_type), intent(in) :: obj !< The file object
+ type(subRegion_type) :: res
+  res = obj%diag_yaml_file%get_file_sub_region()
+end function get_file_sub_region
 
 !> \brief Returns a copy of file_new_file_freq from the yaml object
 !! \return Copy of file_new_file_freq
@@ -537,27 +542,51 @@ subroutine set_file_domain(this, domain, type_of_domain)
 end subroutine set_file_domain
 
 !> @brief Loops through a variable's axis_ids and adds them to the FMSDiagFile object if they don't exist
-subroutine add_axes(this, axis_ids)
-  class(fmsDiagFile_type), intent(inout)       :: this            !< The file object
-  integer,                 INTENT(in)          :: axis_ids(:)    !< Array of axes_ids
+subroutine add_axes(this, axis_ids, diag_axis, naxis)
+  class(fmsDiagFile_type),          intent(inout)       :: this          !< The file object
+  integer,                          INTENT(in)          :: axis_ids(:)   !< Array of axes_ids
+  class(fmsDiagAxisContainer_type), intent(inout)       :: diag_axis(:)  !< Diag_axis object
+  integer,                          intent(inout)       :: naxis         !< Number of axis that have been registered
 
   integer :: i, j !< For do loops
+  logical :: is_cube_sphere !< Flag indicating if the file's domain is a cubesphere
 
-  do i = 1, size(axis_ids)
-    do j = 1, this%number_of_axis
-      !> Check if the axis already exists, return
-      if (axis_ids(i) .eq. this%axis_ids(j)) return
-    enddo
-
-    !> If the axis does not exist add it to the list
-    this%number_of_axis = this%number_of_axis + 1
-    this%axis_ids(this%number_of_axis) = axis_ids(i)
+  is_cube_sphere = .false.
 
-    !> If this is a sub_regional file, set up the sub_axes
-    !> TO DO:
-    !!
-  enddo
+  select type(this)
+  type is (subRegionalFile_type)
+    if (.not. this%is_subaxis_defined) then
+      if (associated(this%domain)) then
+        if (this%domain%get_ntiles() .eq. 6) is_cube_sphere = .true.
+      endif
 
+      call define_subaxis(diag_axis, axis_ids, naxis, this%get_file_sub_region(), &
+        is_cube_sphere, this%write_on_this_pe)
+      this%is_subaxis_defined = .true.
+
+      !> add the axis to the list of axis in the file
+      if (this%write_on_this_pe) then
+        do i = 1, size(axis_ids)
+          this%number_of_axis = this%number_of_axis + 1 !< This is the current number of axis in the file
+          this%axis_ids(this%number_of_axis) = diag_axis(axis_ids(i))%axis%get_subaxes_id()
+        enddo
+      else
+        this%axis_ids = diag_null
+      endif
+    endif
+    return
+  type is (fmsDiagFile_type)
+    do i = 1, size(axis_ids)
+      do j = 1, this%number_of_axis
+        !> Check if the axis already exists, return
+        if (axis_ids(i) .eq. this%axis_ids(j)) return
+      enddo
+
+      !> If the axis does not exist add it to the list
+      this%number_of_axis = this%number_of_axis + 1
+      this%axis_ids(this%number_of_axis) = axis_ids(i)
+    enddo
+  end select
 end subroutine add_axes
 
 !> @brief adds the start time to the fileobj
@@ -592,7 +621,7 @@ subroutine add_start_time(this, start_time)
 subroutine dump_file_obj(this, unit_num)
   class(fmsDiagFile_type), intent(in) :: this !< the file object
   integer, intent(in) :: unit_num !< passed in from dump_diag_obj
-                                  !! will either be for new log file or stdout 
+                                  !! will either be for new log file or stdout
   write( unit_num, *) 'file id:', this%id
   write( unit_num, *) 'start time:', date_to_string(this%start_time)
   write( unit_num, *) 'last_output', date_to_string(this%last_output)
@@ -613,9 +642,11 @@ subroutine dump_file_obj(this, unit_num)
 end subroutine
 
 !< @brief Opens the diag_file if it is time to do so
-subroutine open_diag_file(this, time_step)
+subroutine open_diag_file(this, time_step, file_is_opened)
   class(fmsDiagFileContainer_type), intent(inout), target :: this            !< The file object
   TYPE(time_type),                  intent(in)            :: time_step       !< Current model step time
+  logical,                          intent(out)           :: file_is_opened  !< .true. if the file was opened in this
+                                                                             !! time
 
   class(fmsDiagFile_type), pointer     :: diag_file      !< Diag_file object to open
   class(diagDomain_t),     pointer     :: domain         !< The domain used in the file
@@ -642,6 +673,7 @@ subroutine open_diag_file(this, time_step)
   diag_file => this%FMS_diag_file
   domain => diag_file%domain
 
+  file_is_opened = .false.
   !< Go away if it is not time to open the file
   if (diag_file%next_open > time_step) return
 
@@ -751,9 +783,64 @@ subroutine open_diag_file(this, time_step)
     diag_file%next_open = diag_time_inc(diag_file%next_open, VERY_LARGE_FILE_FREQ, DIAG_DAYS)
   endif
 
-!TODO: closing the file here for now, just to see if it works
-  call close_file(diag_file%fileobj)
+  file_is_opened = .true.
+  domain => null()
+  diag_file => null()
 end subroutine open_diag_file
 
+!< @brief Writes the axis metadata for the file
+subroutine write_axis_metadata(this, diag_axis)
+  class(fmsDiagFileContainer_type), intent(inout), target :: this            !< The file object
+  class(fmsDiagAxisContainer_type), intent(in)            :: diag_axis(:)    !< Diag_axis object
+
+  class(fmsDiagFile_type), pointer     :: diag_file      !< Diag_file object to open
+  class(FmsNetcdfFile_t),  pointer     :: fileobj        !< The fileobj to write to
+  integer                              :: i              !< For do loops
+  integer                              :: j              !< diag_file%axis_ids(i) (for less typing)
+  integer                              :: parent_axis_id !< Id of the parent_axis
+
+  diag_file => this%FMS_diag_file
+  fileobj => diag_file%fileobj
+
+  do i = 1, diag_file%number_of_axis
+    j = diag_file%axis_ids(i)
+    parent_axis_id = diag_axis(j)%axis%get_parent_axis_id()
+    if (parent_axis_id .eq. DIAG_NULL) then
+      call diag_axis(j)%axis%write_axis_metadata(fileobj)
+    else
+      call diag_axis(j)%axis%write_axis_metadata(fileobj, diag_axis(parent_axis_id)%axis)
+    endif
+  enddo
+
+end subroutine write_axis_metadata
+
+!< @brief Writes the axis data for the file
+subroutine write_axis_data(this, diag_axis)
+  class(fmsDiagFileContainer_type), intent(inout), target :: this            !< The file object
+  class(fmsDiagAxisContainer_type), intent(in)            :: diag_axis(:)    !< Diag_axis object
+
+  class(fmsDiagFile_type), pointer     :: diag_file      !< Diag_file object to open
+  class(FmsNetcdfFile_t),  pointer     :: fileobj        !< The fileobj to write to
+  integer                              :: i              !< For do loops
+  integer                              :: j              !< diag_file%axis_ids(i) (for less typing)
+  integer                              :: parent_axis_id !< Id of the parent_axis
+
+  diag_file => this%FMS_diag_file
+  fileobj => diag_file%fileobj
+
+  do i = 1, diag_file%number_of_axis
+    j = diag_file%axis_ids(i)
+    parent_axis_id = diag_axis(j)%axis%get_parent_axis_id()
+    if (parent_axis_id .eq. DIAG_NULL) then
+      call diag_axis(j)%axis%write_axis_data(fileobj)
+    else
+      call diag_axis(j)%axis%write_axis_data(fileobj, diag_axis(parent_axis_id)%axis)
+    endif
+  enddo
+
+  !TODO: closing the file here for now, just to see if it works
+  call close_file(fileobj)
+end subroutine write_axis_data
+
 #endif
 end module fms_diag_file_object_mod

diag_manager/fms_diag_object.F90

@@ -211,15 +211,15 @@ integer function fms_register_diag_field_obj &
      fileptr => this%FMS_diag_files(file_ids(i))%FMS_diag_file
      call fileptr%add_field_id(fieldptr%get_id())
      call fileptr%set_file_domain(fieldptr%get_domain(), fieldptr%get_type_of_domain())
-     call fileptr%add_axes(axes)
+     call fileptr%add_axes(axes, this%diag_axis, this%registered_axis)
      call fileptr%add_start_time(init_time)
     enddo
   elseif (present(axes)) then !only axes present
     do i = 1, size(file_ids)
      fileptr => this%FMS_diag_files(file_ids(i))%FMS_diag_file
      call fileptr%add_field_id(fieldptr%get_id())
      call fileptr%set_file_domain(fieldptr%get_domain(), fieldptr%get_type_of_domain())
-     call fileptr%add_axes(axes)
+     call fileptr%add_axes(axes, this%diag_axis, this%registered_axis)
     enddo
   elseif (present(init_time)) then !only inti time present
     do i = 1, size(file_ids)
@@ -409,6 +409,7 @@ FUNCTION fms_diag_axis_init(this, axis_name, axis_data, units, cart_name, long_n
       & req=req, tile_count=tile_count, domain_position=domain_position)
 
     id = this%registered_axis
+    call axis%set_axis_id(id)
   end select
 #endif
 end function fms_diag_axis_init
@@ -426,9 +427,16 @@ subroutine fms_diag_send_complete(this, time_step)
 #else
   class(fmsDiagFileContainer_type), pointer :: diag_file !< Pointer to this%FMS_diag_files(i) (for convenience)
 
+  logical :: file_is_opened_this_time_step !< True if the file was opened in this time_step
+                                           !! If true the metadata will need to be written
+
   do i = 1, size(this%FMS_diag_files)
     diag_file => this%FMS_diag_files(i)
-    call diag_file%open_diag_file(time_step)
+    call diag_file%open_diag_file(time_step, file_is_opened_this_time_step)
+    if (file_is_opened_this_time_step) then
+      call diag_file%write_axis_metadata(this%diag_axis)
+      call diag_file%write_axis_data(this%diag_axis)
+    endif
   enddo
 #endif
 
@@ -619,7 +627,7 @@ subroutine dump_diag_obj( filename )
     write(unit_num, *) 'axes_initialized:', fms_diag_object%axes_initialized
     write(unit_num, *) 'Files:'
     if( fms_diag_object%files_initialized ) then
-      do i=1, SIZE(fms_diag_object%FMS_diag_files) 
+      do i=1, SIZE(fms_diag_object%FMS_diag_files)
         write(unit_num, *) 'File num:', i
         fileptr => fms_diag_object%FMS_diag_files(i)%FMS_diag_file
         call fileptr%dump_file_obj(unit_num)
@@ -628,7 +636,7 @@ subroutine dump_diag_obj( filename )
       write(unit_num, *) 'files not initialized'
     endif
     if( fms_diag_object%fields_initialized) then
-      do i=1, SIZE(fms_diag_object%FMS_diag_fields) 
+      do i=1, SIZE(fms_diag_object%FMS_diag_fields)
         write(unit_num, *) 'Field num:', i
         fieldptr => fms_diag_object%FMS_diag_fields(i)
         call fieldptr%dump_field_obj(unit_num)
@@ -642,5 +650,4 @@ subroutine dump_diag_obj( filename )
   call mpp_error( FATAL, "You can not use the modern diag manager without compiling with -Duse_yaml")
 #endif
 end subroutine
-
 end module fms_diag_object_mod

diag_manager/fms_diag_yaml.F90

@@ -47,7 +47,7 @@ module fms_diag_yaml_mod
 
 public :: diag_yaml
 public :: diag_yaml_object_init, diag_yaml_object_end
-public :: diagYamlObject_type, get_diag_yaml_obj
+public :: diagYamlObject_type, get_diag_yaml_obj, subRegion_type
 public :: diagYamlFiles_type, diagYamlFilesVar_type
 public :: get_num_unique_fields, find_diag_field, get_diag_fields_entries, get_diag_files_id
 public :: dump_diag_yaml_obj
@@ -1337,20 +1337,20 @@ end function get_diag_files_id
 
 !> Prints out values from diag_yaml object for debugging.
 !! Only writes on root.
-subroutine dump_diag_yaml_obj( filename ) 
+subroutine dump_diag_yaml_obj( filename )
   character(len=*), optional, intent(in)        :: filename !< optional name of logfile to write to, otherwise
                                                             !! prints to stdout
   type(diagyamlfilesvar_type), allocatable      :: fields(:)
   type(diagyamlfiles_type), allocatable         :: files(:)
   integer                                       :: i, unit_num
   if( present(filename)) then
-    open(newunit=unit_num, file=trim(filename), action='WRITE') 
+    open(newunit=unit_num, file=trim(filename), action='WRITE')
   else
-    unit_num = stdout() 
+    unit_num = stdout()
   endif
   !! TODO write to log
   if( mpp_pe() .eq. mpp_root_pe()) then
-    write(unit_num, *) '**********Dumping diag_yaml object**********' 
+    write(unit_num, *) '**********Dumping diag_yaml object**********'
     if( diag_yaml%has_diag_title())    write(unit_num, *) 'Title:', diag_yaml%diag_title
     if( diag_yaml%has_diag_basedate()) write(unit_num, *) 'basedate array:', diag_yaml%diag_basedate
     write(unit_num, *) 'FILES'
@@ -1359,7 +1359,7 @@ subroutine dump_diag_yaml_obj( filename )
     files = diag_yaml%get_diag_files()
     fields = diag_yaml%get_diag_fields()
     do i=1, SIZE(files)
-      write(unit_num, *) 'File: ', files(i)%get_file_fname() 
+      write(unit_num, *) 'File: ', files(i)%get_file_fname()
       if(files(i)%has_file_frequnit()) write(unit_num, *) 'file_frequnit:', files(i)%get_file_frequnit()
       if(files(i)%has_file_freq()) write(unit_num, *) 'freq:', files(i)%get_file_freq()
       if(files(i)%has_file_timeunit()) write(unit_num, *) 'timeunit:', files(i)%get_file_timeunit()
@@ -1378,7 +1378,7 @@ subroutine dump_diag_yaml_obj( filename )
     enddo
     write(unit_num, *) 'FIELDS'
     do i=1, SIZE(fields)
-      write(unit_num, *) 'Field: ', fields(i)%get_var_fname() 
+      write(unit_num, *) 'Field: ', fields(i)%get_var_fname()
       if(fields(i)%has_var_fname()) write(unit_num, *) 'fname:', fields(i)%get_var_fname()
       if(fields(i)%has_var_varname()) write(unit_num, *) 'varname:', fields(i)%get_var_varname()
       if(fields(i)%has_var_reduction()) write(unit_num, *) 'reduction:', fields(i)%get_var_reduction()