Diag_mediator: Fix non-standard axis support

When a diagnostic is registered, its axis group is compared to a
standard set of groups along the model grid. If it matches one of them,
the diagnostic is associated with this axis group.

When it does not match any of the axes, the intention was to retain the
pointer to the input axis group. However, there was a bug where the
axis group was pointed to the input stack value within the function,
rather than the original axis group. This caused the diagnostic to be
associated with an axis group which was largely nonsense.

This patch fixes this bug by allocating the new axis group inside the
function to heap, and then copying the contents of the input axis group
to the new allocated group. This ensures a permanent reference to the
new axis group.

It also probably creates a minor memory leak, but we can sort that out
as part of the overall memory cleanup PR (currently in preparation).

There was literally only one instance of a new axis, namely the use of
angles axes in the internal tide diagnostics. This patch fixes any runs
which use these diagnostics.

src/framework/MOM_diag_mediator.F90

@@ -1051,6 +1051,7 @@ subroutine define_axes_group(diag_cs, handles, axes, nz, vertical_coordinate_num
   else
     axes%v_cell_method = ''
   endif
+
   if (present(nz)) axes%nz = nz
   if (present(vertical_coordinate_number)) axes%vertical_coordinate_number = vertical_coordinate_number
   if (present(is_h_point)) axes%is_h_point = is_h_point
@@ -1971,39 +1972,44 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time
   type(diag_ctrl), pointer :: diag_cs => NULL()
   type(axes_grp), pointer :: remap_axes => null()
   type(axes_grp), pointer :: axes => null()
-  type(axes_grp), pointer :: paxes => null()
+  type(axes_grp), pointer :: axes_d2 => null()
   integer :: dm_id, i, dl
   character(len=256) :: msg, cm_string
   character(len=256) :: new_module_name
   character(len=480) :: module_list, var_list
   integer :: num_modnm, num_varnm
   logical :: active
 
-  axes => axes_in
-  MOM_missing_value = axes%diag_cs%missing_value
-  if (present(missing_value)) MOM_missing_value = missing_value
-
-  diag_cs => axes%diag_cs
-  dm_id = -1
+  diag_cs => axes_in%diag_cs
 
+  ! Check if the axes match a standard grid axis.
+  ! If not, allocate the new axis and copy the contents.
   if (axes_in%id == diag_cs%axesTL%id) then
     axes => diag_cs%axesTL
   elseif (axes_in%id == diag_cs%axesBL%id) then
     axes => diag_cs%axesBL
-  elseif (axes_in%id == diag_cs%axesCuL%id ) then
+  elseif (axes_in%id == diag_cs%axesCuL%id) then
     axes => diag_cs%axesCuL
   elseif (axes_in%id == diag_cs%axesCvL%id) then
     axes => diag_cs%axesCvL
   elseif (axes_in%id == diag_cs%axesTi%id) then
     axes => diag_cs%axesTi
   elseif (axes_in%id == diag_cs%axesBi%id) then
     axes => diag_cs%axesBi
-  elseif (axes_in%id == diag_cs%axesCui%id ) then
+  elseif (axes_in%id == diag_cs%axesCui%id) then
     axes => diag_cs%axesCui
   elseif (axes_in%id == diag_cs%axesCvi%id) then
     axes => diag_cs%axesCvi
+  else
+    allocate(axes)
+    axes = axes_in
   endif
-  paxes => axes
+
+  MOM_missing_value = axes%diag_cs%missing_value
+  if (present(missing_value)) MOM_missing_value = missing_value
+
+  diag_cs => axes%diag_cs
+  dm_id = -1
 
   module_list = "{"//trim(module_name)
   num_modnm = 1
@@ -2092,40 +2098,40 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time
     new_module_name = trim(module_name)//'_d2'
 
     if (axes_in%rank == 3 .or. axes_in%rank == 2 ) then
-      axes => null()
+      axes_d2 => null()
       if (axes_in%id == diag_cs%axesTL%id) then
-        axes => diag_cs%dsamp(dl)%axesTL
+        axes_d2 => diag_cs%dsamp(dl)%axesTL
       elseif (axes_in%id == diag_cs%axesBL%id) then
-        axes => diag_cs%dsamp(dl)%axesBL
+        axes_d2 => diag_cs%dsamp(dl)%axesBL
       elseif (axes_in%id == diag_cs%axesCuL%id ) then
-        axes => diag_cs%dsamp(dl)%axesCuL
+        axes_d2 => diag_cs%dsamp(dl)%axesCuL
       elseif (axes_in%id == diag_cs%axesCvL%id) then
-        axes => diag_cs%dsamp(dl)%axesCvL
+        axes_d2 => diag_cs%dsamp(dl)%axesCvL
       elseif (axes_in%id == diag_cs%axesTi%id) then
-        axes => diag_cs%dsamp(dl)%axesTi
+        axes_d2 => diag_cs%dsamp(dl)%axesTi
       elseif (axes_in%id == diag_cs%axesBi%id) then
-        axes => diag_cs%dsamp(dl)%axesBi
+        axes_d2 => diag_cs%dsamp(dl)%axesBi
       elseif (axes_in%id == diag_cs%axesCui%id ) then
-        axes => diag_cs%dsamp(dl)%axesCui
+        axes_d2 => diag_cs%dsamp(dl)%axesCui
       elseif (axes_in%id == diag_cs%axesCvi%id) then
-        axes => diag_cs%dsamp(dl)%axesCvi
+        axes_d2 => diag_cs%dsamp(dl)%axesCvi
       elseif (axes_in%id == diag_cs%axesT1%id) then
-        axes => diag_cs%dsamp(dl)%axesT1
+        axes_d2 => diag_cs%dsamp(dl)%axesT1
       elseif (axes_in%id == diag_cs%axesB1%id) then
-        axes => diag_cs%dsamp(dl)%axesB1
+        axes_d2 => diag_cs%dsamp(dl)%axesB1
       elseif (axes_in%id == diag_cs%axesCu1%id ) then
-        axes => diag_cs%dsamp(dl)%axesCu1
+        axes_d2 => diag_cs%dsamp(dl)%axesCu1
       elseif (axes_in%id == diag_cs%axesCv1%id) then
-        axes => diag_cs%dsamp(dl)%axesCv1
+        axes_d2 => diag_cs%dsamp(dl)%axesCv1
       else
         !Niki: Should we worry about these, e.g., diag_to_Z_CS?
         call MOM_error(WARNING,"register_diag_field: Could not find a proper axes for " &
               //trim(new_module_name)//"-"//trim(field_name))
       endif
     endif
     ! Register the native diagnostic
-    if (associated(axes)) then
-       active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes, &
+    if (associated(axes_d2)) then
+       active = register_diag_field_expand_cmor(dm_id, new_module_name, field_name, axes_d2, &
                 init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, &
                 range=range, mask_variant=mask_variant, standard_name=standard_name, &
                 verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, &
@@ -2193,7 +2199,7 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time
   if (is_root_pe() .and. (diag_CS%available_diag_doc_unit > 0)) then
     msg = ''
     if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"'
-    call attach_cell_methods(-1, paxes, cm_string, cell_methods, &
+    call attach_cell_methods(-1, axes, cm_string, cell_methods, &
                              x_cell_method, y_cell_method, v_cell_method, &
                              v_extensive=v_extensive)
     module_list = trim(module_list)//"}"
@@ -2218,7 +2224,7 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name,
   integer,          intent(inout) :: dm_id !< The diag_mediator ID for this diagnostic group
   character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" or "ice_shelf_model"
   character(len=*), intent(in) :: field_name !< Name of the diagnostic field
-  type(axes_grp), target, intent(in) :: axes !< Container w/ up to 3 integer handles that indicates axes
+  type(axes_grp),   intent(in) :: axes !< Container w/ up to 3 integer handles that indicates axes
                                              !! for this field
   type(time_type),  intent(in) :: init_time !< Time at which a field is first available?
   character(len=*), optional, intent(in) :: long_name !< Long name of a field.