Generalize support for downsample diagnostics

- Generalize _d2 diagnostics to arbitrary _dl
- The maximum level MAX_DSAMP_LEV needs to be specified. Default is 4.
- This branch is based on the SPREAHI  code base hash 2d121dc

Author: nikizadehgfdl

config_src/infra/FMS2/MOM_diag_manager_infra.F90

@@ -1,3 +1,7 @@
+! This file is part of MOM6, the Modular Ocean Model version 6.
+! See the LICENSE file for licensing information.
+! SPDX-License-Identifier: Apache-2.0
+
 !> A wrapper for the FMS diag_manager routines. This module should be the
 !! only MOM6 module which imports the FMS shared infrastructure for
 !! diagnostics. Pass through interfaces are being documented
@@ -6,21 +10,24 @@
 !! those APIs would be applied here).
 module MOM_diag_manager_infra
 
-! This file is part of MOM6. See LICENSE.md for the license.
+#define MAX_DSAMP_LEV 4
 
+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
 use diag_axis_mod,    only : EAST, NORTH
 use diag_data_mod,    only : null_axis_id
 use diag_manager_mod, only : fms_diag_manager_init => diag_manager_init
 use diag_manager_mod, only : fms_diag_manager_end => diag_manager_end
+use diag_manager_mod, only : diag_send_complete
+use diag_manager_mod, only : diag_manager_set_time_end
 use diag_manager_mod, only : send_data_fms => send_data
 use diag_manager_mod, only : fms_diag_field_add_attribute => diag_field_add_attribute
 use diag_manager_mod, only : DIAG_FIELD_NOT_FOUND
 use diag_manager_mod, only : register_diag_field_fms => register_diag_field
 use diag_manager_mod, only : register_static_field_fms => register_static_field
 use diag_manager_mod, only : get_diag_field_id_fms => get_diag_field_id
-use MOM_time_manager, only : time_type
+use MOM_time_manager, only : time_type, set_time
 use MOM_domain_infra, only : MOM_domain_type
 use MOM_error_infra,  only : MOM_error => MOM_err, FATAL, WARNING
 
@@ -57,6 +64,8 @@ module MOM_diag_manager_infra
 public MOM_diag_manager_init
 public MOM_diag_manager_end
 public send_data_infra
+public diag_send_complete_infra
+public diag_manager_set_time_end_infra
 public MOM_diag_field_add_attribute
 public register_diag_field_infra
 public register_static_field_infra
@@ -109,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
@@ -357,7 +366,7 @@ end function send_data_infra_3d
 logical function send_data_infra_2d_r8(diag_field_id, field, is_in, ie_in, js_in, je_in, &
                                        time, mask, rmask, weight, err_msg)
   integer,                           intent(in) :: diag_field_id !< The diagnostic manager identifier for this field
-  real(kind=8), dimension(:,:),      intent(in) :: field !< A 2-d array of values being recorded
+  real(kind=real64), dimension(:,:), intent(in) :: field !< A 2-d array of values being recorded
   integer,                 optional, intent(in) :: is_in !< The starting i-index for the data being recorded
   integer,                 optional, intent(in) :: ie_in !< The end i-index for the data being recorded
   integer,                 optional, intent(in) :: js_in !< The starting j-index for the data being recorded
@@ -380,7 +389,7 @@ end function send_data_infra_2d_r8
 logical function send_data_infra_3d_r8(diag_field_id, field, is_in, ie_in, js_in, je_in, ks_in, ke_in, &
                                     time, mask, rmask, weight, err_msg)
   integer,                             intent(in) :: diag_field_id !< The diagnostic manager identifier for this field
-  real(kind=8), dimension(:,:,:),      intent(in) :: field !< A rank 1 array of floating point values being recorded
+  real(kind=real64), dimension(:,:,:), intent(in) :: field !< A rank 1 array of floating point values being recorded
   integer,                   optional, intent(in) :: is_in !< The starting i-index for the data being recorded
   integer,                   optional, intent(in) :: ie_in !< The end i-index for the data being recorded
   integer,                   optional, intent(in) :: js_in !< The starting j-index for the data being recorded
@@ -451,4 +460,19 @@ subroutine MOM_diag_field_add_attribute_i1d(diag_field_id, att_name, att_value)
 
 end subroutine MOM_diag_field_add_attribute_i1d
 
+!> Finishes the diag manager reduction methods as needed for the time_step
+subroutine diag_send_complete_infra ()
+  !! The time_step in the diag_send_complete call is a dummy argument, needed for backwards compatibility
+  !! It won't be used at all when diag_manager_nml::use_modern_diag=.true.
+  !! It won't have any impact when diag_manager_nml::use_modern_diag=.false.
+  call diag_send_complete (set_time(0))
+end subroutine diag_send_complete_infra
+
+!> Sets the time that the simulation ends in the diag manager
+subroutine diag_manager_set_time_end_infra(time)
+  type(time_type),           optional, intent(in) :: time  !< The time the simulation ends
+
+  call diag_manager_set_time_end(time)
+end subroutine diag_manager_set_time_end_infra
+
 end module MOM_diag_manager_infra

config_src/infra/FMS2/MOM_domain_infra.F90

@@ -1,8 +1,10 @@
+! This file is part of MOM6, the Modular Ocean Model version 6.
+! See the LICENSE file for licensing information.
+! SPDX-License-Identifier: Apache-2.0
+
 !> Describes the decomposed MOM domain and has routines for communications across PEs
 module MOM_domain_infra
 
-! This file is part of MOM6. See LICENSE.md for the license.
-
 use MOM_coms_infra,  only : PE_here, root_PE, num_PEs
 use MOM_cpu_clock_infra, only : cpu_clock_begin, cpu_clock_end
 use MOM_error_infra, only : MOM_error=>MOM_err, NOTE, WARNING, FATAL
@@ -32,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 4
 
 implicit none ; private
 
@@ -131,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.
@@ -1212,7 +1215,7 @@ subroutine redistribute_array_2d(Domain1, array1, Domain2, array2, complete)
   ! Local variables
   logical :: do_complete
 
-  do_complete=.true.;if (PRESENT(complete)) do_complete = complete
+  do_complete=.true. ; if (PRESENT(complete)) do_complete = complete
 
   call mpp_redistribute(Domain1, array1, Domain2, array2, do_complete)
 
@@ -1231,7 +1234,7 @@ subroutine redistribute_array_3d(Domain1, array1, Domain2, array2, complete)
   ! Local variables
   logical :: do_complete
 
-  do_complete=.true.;if (PRESENT(complete)) do_complete = complete
+  do_complete=.true. ; if (PRESENT(complete)) do_complete = complete
 
   call mpp_redistribute(Domain1, array1, Domain2, array2, do_complete)
 
@@ -1250,61 +1253,97 @@ subroutine redistribute_array_4d(Domain1, array1, Domain2, array2, complete)
   ! Local variables
   logical :: do_complete
 
-  do_complete=.true.;if (PRESENT(complete)) do_complete = complete
+  do_complete=.true. ; if (PRESENT(complete)) do_complete = complete
 
   call mpp_redistribute(Domain1, array1, Domain2, array2, do_complete)
 
 end subroutine redistribute_array_4d
 
 
 !> Rescale the values of a 4-D array in its computational domain by a constant factor
-subroutine rescale_comp_data_4d(domain, array, scale)
+subroutine rescale_comp_data_4d(domain, array, scale, zero_zeros)
   type(MOM_domain_type),    intent(in)    :: domain !< MOM domain from which to extract information
   real, dimension(:,:,:,:), intent(inout) :: array  !< The array which is having the data in its
                                                     !! computational domain rescaled
   real,                     intent(in)    :: scale  !< A scaling factor by which to multiply the
                                                     !! values in the computational domain of array
-  integer :: is, ie, js, je
+  logical,        optional, intent(in)    :: zero_zeros !< If present and true, convert negative zeros
+                                                    !! into ordinary signless zeros.
+  logical :: unsign_zeros ! If true, convert negative zeros into ordinary signless zeros.
+  integer :: is, ie, js, je, i, j, k, m
 
-  if (scale == 1.0) return
+  unsign_zeros = .false. ; if (present(zero_zeros)) unsign_zeros = zero_zeros
+
+  if ((scale == 1.0) .and. (.not.unsign_zeros)) return
 
   call get_simple_array_i_ind(domain, size(array,1), is, ie)
   call get_simple_array_j_ind(domain, size(array,2), js, je)
-  array(is:ie,js:je,:,:) = scale*array(is:ie,js:je,:,:)
+  if (scale /= 1.0) &
+    array(is:ie,js:je,:,:) = scale*array(is:ie,js:je,:,:)
+
+  if (unsign_zeros) then ! Convert negative zeros into zeros
+    do m=1,size(array,4) ; do k=1,size(array,3) ; do j=js,je ; do i=is,ie
+      if (array(i,j,k,m) == 0.0) array(i,j,k,m) = 0.0
+    enddo ; enddo ; enddo ; enddo
+  endif
 
 end subroutine rescale_comp_data_4d
 
 !> Rescale the values of a 3-D array in its computational domain by a constant factor
-subroutine rescale_comp_data_3d(domain, array, scale)
+subroutine rescale_comp_data_3d(domain, array, scale, zero_zeros)
   type(MOM_domain_type),  intent(in)    :: domain !< MOM domain from which to extract information
   real, dimension(:,:,:), intent(inout) :: array  !< The array which is having the data in its
                                                   !! computational domain rescaled
   real,                   intent(in)    :: scale  !< A scaling factor by which to multiply the
                                                   !! values in the computational domain of array
-  integer :: is, ie, js, je
+  logical,      optional, intent(in)    :: zero_zeros !< If present and true, convert negative zeros
+                                                  !! into ordinary signless zeros.
+  logical :: unsign_zeros ! If true, convert negative zeros into ordinary signless zeros.
+  integer :: is, ie, js, je, i, j, k
+
+  unsign_zeros = .false. ; if (present(zero_zeros)) unsign_zeros = zero_zeros
 
-  if (scale == 1.0) return
+  if ((scale == 1.0) .and. (.not.unsign_zeros)) return
 
   call get_simple_array_i_ind(domain, size(array,1), is, ie)
   call get_simple_array_j_ind(domain, size(array,2), js, je)
-  array(is:ie,js:je,:) = scale*array(is:ie,js:je,:)
+  if (scale /= 1.0) &
+    array(is:ie,js:je,:) = scale*array(is:ie,js:je,:)
+
+  if (unsign_zeros) then ! Convert negative zeros into zeros
+    do k=1,size(array,3) ; do j=js,je ; do i=is,ie
+      if (array(i,j,k) == 0.0) array(i,j,k) = 0.0
+    enddo ; enddo ; enddo
+  endif
 
 end subroutine rescale_comp_data_3d
 
 !> Rescale the values of a 2-D array in its computational domain by a constant factor
-subroutine rescale_comp_data_2d(domain, array, scale)
+subroutine rescale_comp_data_2d(domain, array, scale, zero_zeros)
   type(MOM_domain_type), intent(in)    :: domain !< MOM domain from which to extract information
   real, dimension(:,:),  intent(inout) :: array  !< The array which is having the data in its
                                                  !! computational domain rescaled
   real,                  intent(in)    :: scale  !< A scaling factor by which to multiply the
                                                  !! values in the computational domain of array
-  integer :: is, ie, js, je
+  logical,      optional, intent(in)   :: zero_zeros !< If present and true, convert negative zeros
+                                                  !! into ordinary signless zeros.
+  logical :: unsign_zeros ! If true, convert negative zeros into ordinary signless zeros.
+  integer :: is, ie, js, je, i, j
 
-  if (scale == 1.0) return
+  unsign_zeros = .false. ; if (present(zero_zeros)) unsign_zeros = zero_zeros
+
+  if ((scale == 1.0) .and. (.not.unsign_zeros)) return
 
   call get_simple_array_i_ind(domain, size(array,1), is, ie)
   call get_simple_array_j_ind(domain, size(array,2), js, je)
-  array(is:ie,js:je) = scale*array(is:ie,js:je)
+  if (scale /= 1.0) &
+    array(is:ie,js:je) = scale*array(is:ie,js:je)
+
+  if (unsign_zeros) then ! Convert negative zeros into zeros
+    do j=js,je ; do i=is,ie
+      if (array(i,j) == 0.0) array(i,j) = 0.0
+    enddo ; enddo
+  endif
 
 end subroutine rescale_comp_data_2d
 
@@ -1333,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)
@@ -1354,8 +1393,10 @@ subroutine create_MOM_domain(MOM_dom, n_global, n_halo, reentrant, tripolar_N, l
       "TRIPOLAR_N and REENTRANT_Y may not be used together.")
   endif
 
-  MOM_dom%nonblocking_updates = nonblocking
-  MOM_dom%thin_halo_updates = thin_halos
+  MOM_dom%nonblocking_updates = .false.
+  if (present(nonblocking)) MOM_dom%nonblocking_updates = nonblocking
+  MOM_dom%thin_halo_updates = .false.
+  if (present(thin_halos)) MOM_dom%thin_halo_updates = thin_halos
   MOM_dom%symmetric = .true. ; if (present(symmetric)) MOM_dom%symmetric = symmetric
   MOM_dom%niglobal = n_global(1) ; MOM_dom%njglobal = n_global(2)
   MOM_dom%nihalo = n_halo(1) ; MOM_dom%njhalo = n_halo(2)
@@ -1406,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
@@ -1420,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
 
@@ -1428,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)
@@ -1527,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
@@ -1542,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
@@ -1664,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
 
@@ -1802,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