Skip to content

Commit

Permalink
+(*)Add tmp_scale arguments to global_area_mean
Browse files Browse the repository at this point in the history 
  Added new tmp_scale optional arguments to global_area_mean, global_layer_mean,
global_area_mean_u, global_area_mean_v, global_area_integral, global_volume_mean
and global_mass_integral.  This argument allows the reproducing sums to work
with rescaled variables without undoing the scaling of the returned variable.
Also corrected the area_rescaling in global_area_mean_u and global_area_mean_v,
which were substantially changing answers when horizontal distances were being
rescaled (i.e., if L_RESCALE_POWER is not 0).  These routines had only been
added recently, so this change is only changing answers with HOMOGENIZE_FORCINGS
= True.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Feb 27, 2022
1 parent a468bee commit b3e80f7
Showing 1 changed file with 78 additions and 27 deletions.
105 changes: 78 additions & 27 deletions src/diagnostics/MOM_spatial_means.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,89 +25,121 @@ module MOM_spatial_means
contains

!> Return the global area mean of a variable. This uses reproducing sums.
function global_area_mean(var, G, scale)
function global_area_mean(var, G, scale, tmp_scale)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
real, dimension(SZI_(G), SZJ_(G)), intent(in) :: var !< The variable to average
real, optional, intent(in) :: scale !< A rescaling factor for the variable

real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming
real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
!! variable that is reversed in the return value
real :: global_area_mean

! Local variables
real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming
real :: scalefac ! An overall scaling factor for the areas and variable.
real :: temp_scale ! A temporary scaling factor.
integer :: i, j, is, ie, js, je
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec

scalefac = G%US%L_to_m**2 ; if (present(scale)) scalefac = G%US%L_to_m**2*scale
temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
scalefac = G%US%L_to_m**2*temp_scale ; if (present(scale)) scalefac = scalefac * scale

tmpForSumming(:,:) = 0.
do j=js,je ; do i=is,ie
tmpForSumming(i,j) = var(i,j) * (scalefac * G%areaT(i,j) * G%mask2dT(i,j))
enddo ; enddo

global_area_mean = reproducing_sum(tmpForSumming) * G%IareaT_global

if ((temp_scale /= 0.0) .and. (temp_scale /= 1.0)) &
global_area_mean = global_area_mean / temp_scale

end function global_area_mean

!> Return the global area mean of a variable. This uses reproducing sums.
function global_area_mean_v(var, G)
function global_area_mean_v(var, G, tmp_scale)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
real, dimension(SZI_(G), SZJB_(G)), intent(in) :: var !< The variable to average
real, dimension(SZI_(G),SZJB_(G)), intent(in) :: var !< The variable to average
real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
!! variable that is reversed in the return value

real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming
real :: global_area_mean_v

! Local variables
real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming
real :: scalefac ! An overall scaling factor for the areas and variable.
real :: temp_scale ! A temporary scaling factor
integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB

temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
scalefac = G%US%L_to_m**2*temp_scale

tmpForSumming(:,:) = 0.
do J=js,je ; do i=is,ie
tmpForSumming(i,j) = G%areaT(i,j) * (var(i,J) * G%mask2dCv(i,J) + &
var(i,J-1) * G%mask2dCv(i,J-1)) &
/ max(1.e-20,G%mask2dCv(i,J)+G%mask2dCv(i,J-1))
tmpForSumming(i,j) = G%areaT(i,j) * scalefac * &
(var(i,J) * G%mask2dCv(i,J) + var(i,J-1) * G%mask2dCv(i,J-1)) / &
max(1.e-20, G%mask2dCv(i,J)+G%mask2dCv(i,J-1))
enddo ; enddo
global_area_mean_v = reproducing_sum(tmpForSumming) * G%IareaT_global
if ((temp_scale /= 0.0) .and. (temp_scale /= 1.0)) &
global_area_mean_v = global_area_mean_v / temp_scale

end function global_area_mean_v

!> Return the global area mean of a variable on U grid. This uses reproducing sums.
function global_area_mean_u(var, G)
function global_area_mean_u(var, G, tmp_scale)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
real, dimension(SZIB_(G), SZJ_(G)), intent(in) :: var !< The variable to average
real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: var !< The variable to average
real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
!! variable that is reversed in the return value
real :: global_area_mean_u

real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming
real :: global_area_mean_u
real :: scalefac ! An overall scaling factor for the areas and variable.
real :: temp_scale ! A temporary scaling factor
integer :: i, j, is, ie, js, je, isB, ieB, jsB, jeB

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB

temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
scalefac = G%US%L_to_m**2*temp_scale

tmpForSumming(:,:) = 0.
do j=js,je ; do i=is,ie
tmpForSumming(i,j) = G%areaT(i,j) * (var(I,j) * G%mask2dCu(I,j) + &
var(I-1,j) * G%mask2dCu(I-1,j)) &
/ max(1.e-20,G%mask2dCu(I,j)+G%mask2dCu(I-1,j))
tmpForSumming(i,j) = G%areaT(i,j) * scalefac * &
(var(I,j) * G%mask2dCu(I,j) + var(I-1,j) * G%mask2dCu(I-1,j)) / &
max(1.e-20, G%mask2dCu(I,j)+G%mask2dCu(I-1,j))
enddo ; enddo
global_area_mean_u = reproducing_sum(tmpForSumming) * G%IareaT_global
if ((temp_scale /= 0.0) .and. (temp_scale /= 1.0)) &
global_area_mean_u = global_area_mean_u / temp_scale

end function global_area_mean_u

!> Return the global area integral of a variable, by default using the masked area from the
!! grid, but an alternate could be used instead. This uses reproducing sums.
function global_area_integral(var, G, scale, area)
function global_area_integral(var, G, scale, area, tmp_scale)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: var !< The variable to integrate
real, optional, intent(in) :: scale !< A rescaling factor for the variable
real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: area !< The alternate area to use, including
!! any required masking [L2 ~> m2].
!! any required masking [L2 ~> m2].
real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
!! variable that is reversed in the return value
real :: global_area_integral !< The returned area integral, usually in the units of var times [m2].

! Local variables
real, dimension(SZI_(G),SZJ_(G)) :: tmpForSumming
real :: scalefac ! An overall scaling factor for the areas and variable.
real :: temp_scale ! A temporary scaling factor.
integer :: i, j, is, ie, js, je
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec

scalefac = G%US%L_to_m**2 ; if (present(scale)) scalefac = G%US%L_to_m**2*scale
temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
scalefac = G%US%L_to_m**2*temp_scale ; if (present(scale)) scalefac = scalefac * scale

tmpForSumming(:,:) = 0.
if (present(area)) then
Expand All @@ -119,27 +151,35 @@ function global_area_integral(var, G, scale, area)
tmpForSumming(i,j) = var(i,j) * (scalefac * G%areaT(i,j) * G%mask2dT(i,j))
enddo ; enddo
endif

global_area_integral = reproducing_sum(tmpForSumming)

if ((temp_scale /= 0.0) .and. (temp_scale /= 1.0)) &
global_area_integral = global_area_integral / temp_scale

end function global_area_integral

!> Return the layerwise global thickness-weighted mean of a variable. This uses reproducing sums.
function global_layer_mean(var, h, G, GV, scale)
function global_layer_mean(var, h, G, GV, scale, tmp_scale)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: var !< The variable to average
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
real, optional, intent(in) :: scale !< A rescaling factor for the variable
real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
!! variable that is reversed in the return value
real, dimension(SZK_(GV)) :: global_layer_mean

real, dimension(G%isc:G%iec, G%jsc:G%jec, SZK_(GV)) :: tmpForSumming, weight
type(EFP_type), dimension(2*SZK_(GV)) :: laysums
real, dimension(SZK_(GV)) :: global_temp_scalar, global_weight_scalar
real :: temp_scale ! A temporary scaling factor
real :: scalefac ! A scaling factor for the variable.
integer :: i, j, k, is, ie, js, je, nz
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

scalefac = 1.0 ; if (present(scale)) scalefac = scale
temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
scalefac = temp_scale ; if (present(scale)) scalefac = scale * temp_scale
tmpForSumming(:,:,:) = 0. ; weight(:,:,:) = 0.

do k=1,nz ; do j=js,je ; do i=is,ie
Expand All @@ -152,29 +192,33 @@ function global_layer_mean(var, h, G, GV, scale)
call EFP_sum_across_PEs(laysums, 2*nz)

do k=1,nz
global_layer_mean(k) = EFP_to_real(laysums(k)) / EFP_to_real(laysums(nz+k))
global_layer_mean(k) = EFP_to_real(laysums(k)) / (temp_scale * EFP_to_real(laysums(nz+k)))
enddo

end function global_layer_mean

!> Find the global thickness-weighted mean of a variable. This uses reproducing sums.
function global_volume_mean(var, h, G, GV, scale)
function global_volume_mean(var, h, G, GV, scale, tmp_scale)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: var !< The variable being averaged
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
real, optional, intent(in) :: scale !< A rescaling factor for the variable
real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
!! variable that is reversed in the return value
real :: global_volume_mean !< The thickness-weighted average of var

real :: temp_scale ! A temporary scaling factor
real :: scalefac ! A scaling factor for the variable.
real :: weight_here
real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming, sum_weight
integer :: i, j, k, is, ie, js, je, nz
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

scalefac = 1.0 ; if (present(scale)) scalefac = scale
temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
scalefac = temp_scale ; if (present(scale)) scalefac = temp_scale * scale
tmpForSumming(:,:) = 0. ; sum_weight(:,:) = 0.

do k=1,nz ; do j=js,je ; do i=is,ie
Expand All @@ -183,13 +227,13 @@ function global_volume_mean(var, h, G, GV, scale)
sum_weight(i,j) = sum_weight(i,j) + weight_here
enddo ; enddo ; enddo
global_volume_mean = (reproducing_sum(tmpForSumming)) / &
(reproducing_sum(sum_weight))
(temp_scale * reproducing_sum(sum_weight))

end function global_volume_mean


!> Find the global mass-weighted integral of a variable. This uses reproducing sums.
function global_mass_integral(h, G, GV, var, on_PE_only, scale)
function global_mass_integral(h, G, GV, var, on_PE_only, scale, tmp_scale)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
Expand All @@ -199,16 +243,20 @@ function global_mass_integral(h, G, GV, var, on_PE_only, scale)
logical, optional, intent(in) :: on_PE_only !< If present and true, the sum is only
!! done on the local PE, and it is _not_ order invariant.
real, optional, intent(in) :: scale !< A rescaling factor for the variable
real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
!! variable that is reversed in the return value
real :: global_mass_integral !< The mass-weighted integral of var (or 1) in
!! kg times the units of var

real, dimension(SZI_(G), SZJ_(G)) :: tmpForSumming
real :: scalefac ! An overall scaling factor for the areas and variable.
real :: temp_scale ! A temporary scaling factor.
logical :: global_sum
integer :: i, j, k, is, ie, js, je, nz
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

scalefac = G%US%L_to_m**2 ; if (present(scale)) scalefac = G%US%L_to_m**2*scale
temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
scalefac = G%US%L_to_m**2*temp_scale ; if (present(scale)) scalefac = scalefac * scale
tmpForSumming(:,:) = 0.0

if (present(var)) then
Expand All @@ -232,6 +280,9 @@ function global_mass_integral(h, G, GV, var, on_PE_only, scale)
enddo ; enddo
endif

if ((temp_scale /= 0.0) .and. (temp_scale /= 1.0)) &
global_mass_integral = global_mass_integral / temp_scale

end function global_mass_integral

!> Find the global mass-weighted order invariant integral of a variable in mks units,
Expand Down

0 comments on commit b3e80f7

Please sign in to comment.