Cleaned up version of gmac_addhistogramandscalarcoord

src/ALE/MOM_remapping.F90

@@ -84,6 +84,7 @@ module MOM_remapping
 public initialize_remapping, end_remapping, remapping_set_param, extract_member_remapping_CS
 public remapping_unit_tests, build_reconstructions_1d, average_value_ppoly
 public interpolate_column, reintegrate_column, dzFromH1H2
+public histogram_column
 
 ! The following are private parameter constants
 integer, parameter  :: REMAPPING_PCM        = 0 !< O(h^1) remapping scheme
@@ -1385,6 +1386,32 @@ subroutine reintegrate_column(nsrc, h_src, uh_src, ndest, h_dest, uh_dest)
 
 end subroutine reintegrate_column
 
+!> Conservatively histogram data, uh_dest, weights mapping source grid to target grid
+subroutine histogram_column(nsrc, uh_src, ndest, uh_dest, histogram_weights)
+  integer,                intent(in)    :: nsrc    !< Number of source cells
+  real, dimension(nsrc),  intent(in)    :: uh_src  !< Values at source cell interfaces [A H]
+  integer,                intent(in)    :: ndest   !< Number of destination cells
+  real, dimension(ndest), intent(inout) :: uh_dest !< Interpolated value at destination cell interfaces [A H]
+  real, dimension(nsrc,ndest), intent(in) :: histogram_weights !< Weights mapping source to destination grid
+
+  ! Local variables
+  integer :: k_dest, k_src
+  real, dimension(nsrc) :: weights
+  real :: uh_dest_sum
+
+  uh_dest(:) = 0.0
+  do k_dest = 1,ndest ! Loop through destination grid layers
+    weights = histogram_weights(:,k_dest)
+    uh_dest_sum = 0.0
+    do k_src = 1,nsrc ! Loop through source grid
+      ! cumulatively sum field multiplied by the weight associated with each grid
+      uh_dest_sum = uh_dest_sum + uh_src(k_src)*weights(k_src)
+    enddo
+    uh_dest(k_dest) = uh_dest_sum
+  enddo
+
+end subroutine histogram_column
+
 !> Returns the average value of a reconstruction within a single source cell, i0,
 !! between the non-dimensional positions xa and xb (xa<=xb) with dimensional
 !! separation dh.

src/ALE/coord_rho.F90

@@ -7,6 +7,7 @@ module coord_rho
 use MOM_remapping,     only : remapping_CS, remapping_core_h
 use MOM_EOS,           only : EOS_type, calculate_density
 use regrid_interp,     only : interp_CS_type, build_and_interpolate_grid, DEGREE_MAX
+use regrid_interp,     only : build_histogram_weights
 
 implicit none ; private
 
@@ -29,6 +30,9 @@ module coord_rho
   !> Nominal density of interfaces [R ~> kg m-3]
   real, allocatable, dimension(:) :: target_density
 
+  !> If true, use a histogramming procedure for extensive diagnostics (rather than the reintegrate procedure)
+  logical :: histogram_extensive_diags = .false.
+
   !> Interpolation control structure
   type(interp_CS_type) :: interp_CS
 end type rho_CS
@@ -38,11 +42,12 @@ module coord_rho
 contains
 
 !> Initialise a rho_CS with pointers to parameters
-subroutine init_coord_rho(CS, nk, ref_pressure, target_density, interp_CS)
+subroutine init_coord_rho(CS, nk, ref_pressure, target_density, histogram_extensive_diags, interp_CS)
   type(rho_CS),         pointer    :: CS !< Unassociated pointer to hold the control structure
   integer,              intent(in) :: nk !< Number of layers in the grid
   real,                 intent(in) :: ref_pressure !< Coordinate reference pressure [R L2 T-2 ~> Pa]
   real, dimension(:),   intent(in) :: target_density !< Nominal density of interfaces [R ~> kg m-3]
+  logical,   intent(in)            :: histogram_extensive_diags !< Boolean to select how to deal with extensive diagnostics
   type(interp_CS_type), intent(in) :: interp_CS !< Controls for interpolation
 
   if (associated(CS)) call MOM_error(FATAL, "init_coord_rho: CS already associated!")
@@ -52,6 +57,7 @@ subroutine init_coord_rho(CS, nk, ref_pressure, target_density, interp_CS)
   CS%nk                = nk
   CS%ref_pressure      = ref_pressure
   CS%target_density(:) = target_density(:)
+  CS%histogram_extensive_diags = histogram_extensive_diags
   CS%interp_CS         = interp_CS
 
 end subroutine init_coord_rho
@@ -67,21 +73,23 @@ subroutine end_coord_rho(CS)
 end subroutine end_coord_rho
 
 !> This subroutine can be used to set the parameters for the coord_rho module
-subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, interp_CS, ref_pressure)
+subroutine set_rho_params(CS, min_thickness, integrate_downward_for_e, histogram_extensive_diags, interp_CS, ref_pressure)
   type(rho_CS),      pointer    :: CS !< Coordinate control structure
   real,    optional, intent(in) :: min_thickness !< Minimum allowed thickness [H ~> m or kg m-2]
   logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface
                                       !! positions from the top downward.  If false, integrate
                                       !! from the bottom upward, as does the rest of the model.
   real,    optional, intent(in) :: ref_pressure     !< The reference pressure for density-dependent
                                                     !! coordinates [R L2 T-2 ~> Pa]
+  logical,    optional, intent(in) :: histogram_extensive_diags     !< If true, use histogram approach to for outputing extensive diags
 
   type(interp_CS_type), optional, intent(in) :: interp_CS !< Controls for interpolation
 
   if (.not. associated(CS)) call MOM_error(FATAL, "set_rho_params: CS not associated")
 
   if (present(min_thickness)) CS%min_thickness = min_thickness
   if (present(integrate_downward_for_e)) CS%integrate_downward_for_e = integrate_downward_for_e
+  if (present(histogram_extensive_diags)) CS%histogram_extensive_diags = histogram_extensive_diags
   if (present(interp_CS)) CS%interp_CS = interp_CS
   if (present(ref_pressure)) CS%ref_pressure = ref_pressure
 end subroutine set_rho_params
@@ -91,7 +99,7 @@ end subroutine set_rho_params
 !! 1. Density profiles are calculated on the source grid.
 !! 2. Positions of target densities (for interfaces) are found by interpolation.
 subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, &
-                            z_rigid_top, eta_orig, h_neglect, h_neglect_edge)
+                            histogram_weights, z_rigid_top, eta_orig, h_neglect, h_neglect_edge)
   type(rho_CS),        intent(in)    :: CS !< coord_rho control structure
   integer,             intent(in)    :: nz !< Number of levels on source grid (i.e. length of  h, T, S)
   real,                intent(in)    :: depth !< Depth of ocean bottom (positive downward) [H ~> m or kg m-2]
@@ -109,6 +117,7 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, &
                                              !! of cell reconstructions [H ~> m or kg m-2]
   real,      optional, intent(in)    :: h_neglect_edge !< A negligibly small width for the purpose
                                              !! of edge value calculations [H ~> m or kg m-2]
+  real, optional, dimension(nz,CS%nk), intent (inout) :: histogram_weights ! Matrix of weights mapping source grid cells to target grid layers
 
   ! Local variables
   integer :: k, count_nonzero_layers
@@ -141,6 +150,12 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, &
                                     h_nv, xTmp, CS%target_density, CS%nk, h_new, &
                                     x1, h_neglect, h_neglect_edge)
 
+    if (CS%histogram_extensive_diags) then
+      ! Based on source column density profile, derive weights that map source to target grid
+      call build_histogram_weights(CS%interp_CS, densities, nz, h, CS%target_density, CS%nk, &
+                                    histogram_weights, h_neglect, h_neglect_edge)
+    endif
+
     ! Inflate vanished layers
     call old_inflate_layers_1d(CS%min_thickness, CS%nk, h_new)