Generalize support for downsample diagnostics

- Generalize _d2 diagnostics to arbitrary _dl
- The maximum level MAX_DSAMP_LEV needs to be specified.
  This limitation could later be enhanced by getting the desired levels of
  downsampling from MOM_parameter file, e.g., DOWNSAMPLE_DIAG_LEVELS = /2,3,4/

Author: nikizadehgfdl

config_src/infra/FMS2/MOM_diag_manager_infra.F90

@@ -10,6 +10,8 @@
 !! those APIs would be applied here).
 module MOM_diag_manager_infra
 
+#define MAX_DSAMP_LEV 3
+
 use, intrinsic :: iso_fortran_env, only : real64
 use diag_axis_mod,    only : fms_axis_init=>diag_axis_init
 use diag_axis_mod,    only : fms_get_diag_axis_name => get_diag_axis_name
@@ -116,10 +118,10 @@ integer function MOM_diag_axis_init(name, data, units, cart_name, long_name, MOM
       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
+    elseif (coarsening <= MAX_DSAMP_LEV) 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_d2, domain_position=position)
+              domain2=MOM_domain%mpp_domain_d(coarsening), domain_position=position)
     else
       call MOM_error(FATAL, "diag_axis_init called with an invalid value of coarsen.")
     endif

config_src/infra/FMS2/MOM_domain_infra.F90

@@ -34,6 +34,7 @@ module MOM_domain_infra
 
 ! The `group_pass_type` fields are never accessed, so we keep it as an FMS type
 use mpp_domains_mod, only : group_pass_type => mpp_group_update_type
+#define MAX_DSAMP_LEV 3
 
 implicit none ; private
 
@@ -133,7 +134,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 +1372,14 @@ 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
+  integer :: dl
   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)
+    do dl=2,MAX_DSAMP_LEV ; allocate(MOM_dom%mpp_domain_d(dl)) ; enddo
   endif
 
   MOM_dom%name = "MOM" ; if (present(domain_name)) MOM_dom%name = trim(domain_name)
@@ -1446,11 +1447,12 @@ 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)
+  do dl=2,MAX_DSAMP_LEV
+    !Downsample diagnostics calculations do not need halos. 
+    call clone_MD_to_d2D(MOM_dom, MOM_dom%mpp_domain_d(dl), coarsen=dl, halo_size=0, &
+                        domain_name="MOM_domain_d" // char(48+dl))
+  enddo
+
 end subroutine create_MOM_domain
 
 !> dealloc_MOM_domain deallocates memory associated with a pointer to a MOM_domain_type
@@ -1460,6 +1462,7 @@ subroutine deallocate_MOM_domain(MOM_domain, cursory)
   logical,  optional, intent(in) :: cursory    !< If true do not deallocate fields associated
                                                !! with the underlying infrastructure
   logical :: invasive  ! If true, deallocate fields associated with the underlying infrastructure
+  integer :: dl
 
   invasive = .true. ; if (present(cursory)) invasive = .not.cursory
 
@@ -1468,9 +1471,11 @@ 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
+      if (invasive) then
+        do dl=2,MAX_DSAMP_LEV ; call mpp_deallocate_domain(MOM_domain%mpp_domain_d(dl)); enddo
+      endif
+      deallocate(MOM_domain%mpp_domain_d)
     endif
     if (associated(MOM_domain%maskmap)) deallocate(MOM_domain%maskmap)
     deallocate(MOM_domain)
@@ -1567,7 +1572,7 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, domain
   integer, dimension(:), allocatable :: exnj ! The extents of the grid for each j-row of the layout.
                                              ! The sum of exni must equal MOM_dom%niglobal.
   integer :: qturns ! The number of quarter turns, restricted to the range of 0 to 3.
-  integer :: i, j, nl1, nl2
+  integer :: i, j, nl1, nl2, dl
   integer :: io_layout_in(2)
 
   qturns = 0
@@ -1582,7 +1587,7 @@ 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)
+    do dl=2,MAX_DSAMP_LEV ; allocate(MOM_dom%mpp_domain_d(dl)) ; enddo
   endif
 
 ! Save the extra data for creating other domains of different resolution that overlay this domain
@@ -1704,7 +1709,11 @@ 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)
+  do dl=2,MAX_DSAMP_LEV
+    !Downsample diagnostics calculations do not need halos. 
+    call clone_MD_to_d2D(MOM_dom, MOM_dom%mpp_domain_d(dl), coarsen=dl, halo_size=0, &
+                        domain_name="MOM_domain_d" // char(48+dl))
+  enddo
 
 end subroutine clone_MD_to_MD
 
@@ -1842,12 +1851,12 @@ subroutine get_domain_extent_MD(Domain, isc, iec, jsc, jec, isd, ied, jsd, jed,
     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_)
+  elseif (coarsen_lev <= MAX_DSAMP_LEV) then
+    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_)
   else
     call MOM_error(FATAL, "get_domain_extent called with an unsupported level of coarsening.")
   endif

src/core/MOM_grid.F90

@@ -15,6 +15,7 @@ module MOM_grid
 implicit none ; private
 
 #include <MOM_memory.h>
+#define MAX_DSAMP_LEV 3
 
 public MOM_grid_init, MOM_grid_end, set_derived_metrics, set_first_direction
 public isPointInCell, hor_index_type, get_global_grid_size
@@ -29,7 +30,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) :: HId(2:MAX_DSAMP_LEV) !< Horizontal index ranges for downsampling, level 2 to MAX_DSAMP_LEV.
 
   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
@@ -237,7 +238,7 @@ subroutine MOM_grid_init(G, param_file, US, HI, global_indexing, bathymetry_at_v
   integer :: isd, ied, jsd, jed
   integer :: IsdB, IedB, JsdB, JedB
   integer :: ied_max, jed_max
-  integer :: niblock, njblock, nihalo, njhalo, nblocks, n, i, j
+  integer :: niblock, njblock, nihalo, njhalo, nblocks, n, i, j, dl
   logical :: local_indexing  ! If false use global index values instead of having
                              ! the data domain on each processor start at 1.
   ! This include declares and sets the variable "version".
@@ -402,22 +403,25 @@ 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)
+  ! Initialize the global grid extents for all levels of diagnosics coarsening.
+  do dl=2,MAX_DSAMP_LEV
+    call get_domain_extent(G%Domain, G%HId(dl)%isc, G%HId(dl)%iec, G%HId(dl)%jsc, G%HId(dl)%jec, &
+                          G%HId(dl)%isd, G%HId(dl)%ied, G%HId(dl)%jsd, G%HId(dl)%jed, &
+                          G%HId(dl)%isg, G%HId(dl)%ieg, G%HId(dl)%jsg, G%HId(dl)%jeg, coarsen=dl)
 
-  ! 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
+    ! Set array sizes for fields that are discretized at tracer cell boundaries.
+    G%HId(dl)%IscB = G%HId(dl)%isc ; G%HId(dl)%JscB = G%HId(dl)%jsc
+    G%HId(dl)%IsdB = G%HId(dl)%isd ; G%HId(dl)%JsdB = G%HId(dl)%jsd
+    G%HId(dl)%IsgB = G%HId(dl)%isg ; G%HId(dl)%JsgB = G%HId(dl)%jsg
+    if (G%symmetric) then
+      G%HId(dl)%IscB = G%HId(dl)%isc-1 ; G%HId(dl)%JscB = G%HId(dl)%jsc-1
+      G%HId(dl)%IsdB = G%HId(dl)%isd-1 ; G%HId(dl)%JsdB = G%HId(dl)%jsd-1
+      G%HId(dl)%IsgB = G%HId(dl)%isg-1 ; G%HId(dl)%JsgB = G%HId(dl)%jsg-1
+    endif
+    G%HId(dl)%IecB = G%HId(dl)%iec ; G%HId(dl)%JecB = G%HId(dl)%jec
+    G%HId(dl)%IedB = G%HId(dl)%ied ; G%HId(dl)%JedB = G%HId(dl)%jed
+    G%HId(dl)%IegB = G%HId(dl)%ieg ; G%HId(dl)%JegB = G%HId(dl)%jeg
+  enddo
 
 end subroutine MOM_grid_init