Add capability to use multiple diagnostics downsampling levels in MOM6

- This commit presents an enhancement of diag downsampling workflow for MOM6
- It allows an arbitrary number of downsampling factors in the same run (if they are commensurate with grid processor layout)
- The number of downsampling levels should be specified in MOM override parameters
  E.g., #override NUM_DIAG_DOWNSAMP_LEV = 3
- The downampling factors should be specified
  E.g., #override DIAG_DOWNSAMP_LEVS = 2,4,5
- The required downsampled domains and axes are then generated at the initialization by MOM_diag_mediator
  No need for other modules to know anything about the downsampling scheme.
- No need for a #define MAX_DSAMP_LEVEL, everything is allocated at runtime
- Since everything is allocated at runtime users need to be congnisent of the number of levels they request.
- Any downsampling diagnostics can be put in diag_table , only the levels that requested in parameters are handled
- The overall dependence of downsampling levels on the grid shape and pe layout remains
  Only levels should be requested where NIGLOBAL/x_layout and NJGLOBAL/y_layout are divisible by downsampling factor

Author: nikizadehgfdl

config_src/infra/FMS2/MOM_diag_manager_infra.F90

@@ -111,17 +111,15 @@ integer function MOM_diag_axis_init(name, data, units, cart_name, long_name, MOM
   endif ; endif
 
   if (present(MOM_domain)) then
-    coarsening = 1 ; if (present(coarsen)) coarsening = coarsen
-    if (coarsening == 1) then
+    coarsening = 0 ; if (present(coarsen)) coarsening = coarsen
+    if (coarsening == 0) then
       MOM_diag_axis_init = fms_axis_init(name, data, units, cart_name, long_name=long_name, &
               direction=direction, set_name=set_name, edges=edges, &
               domain2=MOM_domain%mpp_domain, domain_position=position)
-    elseif (coarsening == 2) then
+    else
       MOM_diag_axis_init = fms_axis_init(name, data, units, cart_name, long_name=long_name, &
               direction=direction, set_name=set_name, edges=edges, &
-              domain2=MOM_domain%mpp_domain_d2, domain_position=position)
-    else
-      call MOM_error(FATAL, "diag_axis_init called with an invalid value of coarsen.")
+              domain2=MOM_domain%mpp_domain_d(coarsening), domain_position=position)
     endif
   else
     if (present(coarsen)) then ; if (coarsen /= 1) then

config_src/infra/FMS2/MOM_domain_infra.F90

@@ -133,7 +133,7 @@ module MOM_domain_infra
   character(len=64) :: name     !< The name of this domain
   type(domain2D), pointer :: mpp_domain => NULL() !< The FMS domain with halos
                                 !! on this processor, centered at h points.
-  type(domain2D), pointer :: mpp_domain_d2 => NULL() !< A coarse FMS domain with halos
+  type(domain2D), pointer :: mpp_domain_d(:) => NULL() !< A coarse FMS domain with halos
                                 !! on this processor, centered at h points.
   integer :: niglobal           !< The total horizontal i-domain size.
   integer :: njglobal           !< The total horizontal j-domain size.
@@ -1371,14 +1371,12 @@ subroutine create_MOM_domain(MOM_dom, n_global, n_halo, reentrant, tripolar_N, l
   integer, dimension(4) :: global_indices ! The lower and upper global i- and j-index bounds
   integer :: X_FLAGS  ! A combination of integers encoding the x-direction grid connectivity.
   integer :: Y_FLAGS  ! A combination of integers encoding the y-direction grid connectivity.
-  integer :: xhalo_d2, yhalo_d2
   character(len=200) :: mesg    ! A string for use in error messages
   logical :: mask_table_exists  ! Mask_table is present and the file it points to exists
 
   if (.not.associated(MOM_dom)) then
     allocate(MOM_dom)
     allocate(MOM_dom%mpp_domain)
-    allocate(MOM_dom%mpp_domain_d2)
   endif
 
   MOM_dom%name = "MOM" ; if (present(domain_name)) MOM_dom%name = trim(domain_name)
@@ -1446,11 +1444,6 @@ subroutine create_MOM_domain(MOM_dom, n_global, n_halo, reentrant, tripolar_N, l
 
   call clone_MD_to_d2D(MOM_dom, MOM_dom%mpp_domain)
 
-  !For downsampled domain, recommend a halo of 1 (or 0?) since we're not doing wide-stencil computations.
-  !But that does not work because the downsampled field would not have the correct size to pass the checks, e.g., we get
-  !error: downsample_diag_indices_get: peculiar size 28 in i-direction\ndoes not match one of 24 25 26 27
-  ! call clone_MD_to_d2D(MOM_dom, MOM_dom%mpp_domain_d2, halo_size=(MOM_dom%nihalo/2), coarsen=2)
-  call clone_MD_to_d2D(MOM_dom, MOM_dom%mpp_domain_d2, coarsen=2)
 end subroutine create_MOM_domain
 
 !> dealloc_MOM_domain deallocates memory associated with a pointer to a MOM_domain_type
@@ -1468,9 +1461,8 @@ subroutine deallocate_MOM_domain(MOM_domain, cursory)
       if (invasive) call mpp_deallocate_domain(MOM_domain%mpp_domain)
       deallocate(MOM_domain%mpp_domain)
     endif
-    if (associated(MOM_domain%mpp_domain_d2)) then
-      if (invasive) call mpp_deallocate_domain(MOM_domain%mpp_domain_d2)
-      deallocate(MOM_domain%mpp_domain_d2)
+    if (associated(MOM_domain%mpp_domain_d)) then
+      deallocate(MOM_domain%mpp_domain_d)
     endif
     if (associated(MOM_domain%maskmap)) deallocate(MOM_domain%maskmap)
     deallocate(MOM_domain)
@@ -1582,7 +1574,6 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, domain
   if (.not.associated(MOM_dom)) then
     allocate(MOM_dom)
     allocate(MOM_dom%mpp_domain)
-    allocate(MOM_dom%mpp_domain_d2)
   endif
 
 ! Save the extra data for creating other domains of different resolution that overlay this domain
@@ -1704,7 +1695,6 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, domain
   endif
 
   call clone_MD_to_d2D(MOM_dom, MOM_dom%mpp_domain, xextent=exni, yextent=exnj)
-  call clone_MD_to_d2D(MOM_dom, MOM_dom%mpp_domain_d2, domain_name=MOM_dom%name, coarsen=2)
 
 end subroutine clone_MD_to_MD
 
@@ -1836,20 +1826,18 @@ subroutine get_domain_extent_MD(Domain, isc, iec, jsc, jec, isd, ied, jsd, jed,
   local = .true. ; if (present(local_indexing)) local = local_indexing
   ind_off = 0 ; if (present(index_offset)) ind_off = index_offset
 
-  coarsen_lev = 1 ; if (present(coarsen)) coarsen_lev = coarsen
+  coarsen_lev = 0 ; if (present(coarsen)) coarsen_lev = coarsen
 
-  if (coarsen_lev == 1) then
+  if (coarsen_lev == 0) then
     call mpp_get_compute_domain(Domain%mpp_domain, isc, iec, jsc, jec)
     call mpp_get_data_domain(Domain%mpp_domain, isd, ied, jsd, jed)
     call mpp_get_global_domain(Domain%mpp_domain, isg_, ieg_, jsg_, jeg_)
-  elseif (coarsen_lev == 2) then
-    if (.not.associated(Domain%mpp_domain_d2)) call MOM_error(FATAL, &
-            "get_domain_extent called with coarsen=2, but Domain%mpp_domain_d2 is not associated.")
-    call mpp_get_compute_domain(Domain%mpp_domain_d2, isc, iec, jsc, jec)
-    call mpp_get_data_domain(Domain%mpp_domain_d2, isd, ied, jsd, jed)
-    call mpp_get_global_domain(Domain%mpp_domain_d2, isg_, ieg_, jsg_, jeg_)
   else
-    call MOM_error(FATAL, "get_domain_extent called with an unsupported level of coarsening.")
+    if (.not.associated(Domain%mpp_domain_d)) call MOM_error(FATAL, &
+            "get_domain_extent called with coarsen_lev, but Domain%mpp_domain_d(coarsen_lev) is not associated.")
+    call mpp_get_compute_domain(Domain%mpp_domain_d(coarsen_lev), isc, iec, jsc, jec)
+    call mpp_get_data_domain(Domain%mpp_domain_d(coarsen_lev), isd, ied, jsd, jed)
+    call mpp_get_global_domain(Domain%mpp_domain_d(coarsen_lev), isg_, ieg_, jsg_, jeg_)
   endif
 
   if (local) then

src/core/MOM_grid.F90

@@ -29,7 +29,7 @@ module MOM_grid
   type(MOM_domain_type), pointer :: Domain => NULL() !< Ocean model domain
   type(MOM_domain_type), pointer :: Domain_aux => NULL() !< A non-symmetric auxiliary domain type.
   type(hor_index_type) :: HI !< Horizontal index ranges
-  type(hor_index_type) :: HId2 !< Horizontal index ranges for level-2-downsampling
+  type(hor_index_type), allocatable :: HId(:) !< Horizontal index ranges for downsampling
 
   integer :: isc !< The start i-index of cell centers within the computational domain
   integer :: iec !< The end i-index of cell centers within the computational domain
@@ -402,23 +402,6 @@ subroutine MOM_grid_init(G, param_file, US, HI, global_indexing, bathymetry_at_v
   if ( G%block(nblocks)%jed+G%block(nblocks)%jdg_offset > G%HI%jed + G%HI%jdg_offset ) &
         call MOM_error(FATAL, "MOM_grid_init: G%jed_bk > G%jed")
 
-  call get_domain_extent(G%Domain, G%HId2%isc, G%HId2%iec, G%HId2%jsc, G%HId2%jec, &
-                         G%HId2%isd, G%HId2%ied, G%HId2%jsd, G%HId2%jed, &
-                         G%HId2%isg, G%HId2%ieg, G%HId2%jsg, G%HId2%jeg, coarsen=2)
-
-  ! Set array sizes for fields that are discretized at tracer cell boundaries.
-  G%HId2%IscB = G%HId2%isc ; G%HId2%JscB = G%HId2%jsc
-  G%HId2%IsdB = G%HId2%isd ; G%HId2%JsdB = G%HId2%jsd
-  G%HId2%IsgB = G%HId2%isg ; G%HId2%JsgB = G%HId2%jsg
-  if (G%symmetric) then
-    G%HId2%IscB = G%HId2%isc-1 ; G%HId2%JscB = G%HId2%jsc-1
-    G%HId2%IsdB = G%HId2%isd-1 ; G%HId2%JsdB = G%HId2%jsd-1
-    G%HId2%IsgB = G%HId2%isg-1 ; G%HId2%JsgB = G%HId2%jsg-1
-  endif
-  G%HId2%IecB = G%HId2%iec ; G%HId2%JecB = G%HId2%jec
-  G%HId2%IedB = G%HId2%ied ; G%HId2%JedB = G%HId2%jed
-  G%HId2%IegB = G%HId2%ieg ; G%HId2%JegB = G%HId2%jeg
-
 end subroutine MOM_grid_init
 
 !> set_derived_metrics calculates metric terms that are derived from other metrics.