Skip to content

Commit

Permalink
+Reproducing KPP_smooth_BLD when KPP%N_SMOOTH > 1
Browse files Browse the repository at this point in the history
  Revised KPP_smooth_BLD() to reproduce across processor count and layout when
USE_KPP is true and KPP%N_SMOOTH > 1.  The specific changes include adding a
variable with the total ocean depth before doing the iterations, doing a halo
update on this total ocean depth, marching in the working do-loop size with
successive iterations, and moving the code to calculate CS%kOBL into a separate
loop that is exercised after all of the iterations for the smoothing passes on
CS%OBLdepth.  This commit will change answers (so that they reproduce across
processor count and layout) when USE_KPP is true and KPP%N_SMOOTH >= 2, but it
gives bitwise identical answers when KPP%N_SMOOTH <= 1.
  • Loading branch information
Hallberg-NOAA committed Jul 12, 2024
1 parent 6ad1530 commit 920093a
Showing 1 changed file with 57 additions and 26 deletions.
83 changes: 57 additions & 26 deletions src/parameterizations/vertical/MOM_CVMix_KPP.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1342,47 +1342,60 @@ subroutine KPP_smooth_BLD(CS, G, GV, US, dz)
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: dz !< Layer thicknesses [Z ~> m]

! local
! local variables
real, dimension(SZI_(G),SZJ_(G)) :: OBLdepth_prev ! OBLdepth before s.th smoothing iteration [Z ~> m]
real, dimension(SZI_(G),SZJ_(G)) :: total_depth ! The total depth of the water column, adjusted
! for the minimum layer thickness [Z ~> m]
real, dimension( GV%ke ) :: cellHeight ! Cell center heights referenced to surface [Z ~> m]
! (negative in the ocean)
real, dimension( GV%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [Z ~> m]
! (negative in the ocean)
real :: wc, ww, we, wn, ws ! averaging weights for smoothing [nondim]
real :: dh ! The local thickness used for calculating interface positions [Z ~> m]
real :: h_cor(SZI_(G)) ! A cumulative correction arising from inflation of vanished layers [Z ~> m]
real :: hcorr ! A cumulative correction arising from inflation of vanished layers [Z ~> m]
integer :: i, j, k, s
integer :: i, j, k, s, halo

call cpu_clock_begin(id_clock_KPP_smoothing)

! Update halos
! Find the total water column thickness first, as it is reused for each smoothing pass.
total_depth(:,:) = 0.0

!$OMP parallel do default(shared) private(dh, h_cor)
do j = G%jsc, G%jec
h_cor(:) = 0.
do k=1,GV%ke
do i=G%isc,G%iec ; if (G%mask2dT(i,j) > 0.0) then
! This code replicates the interface height calculations below. It could be simpler, as shown below.
dh = dz(i,j,k) ! Nominal thickness to use for increment
dh = dh + h_cor(i) ! Take away the accumulated error (could temporarily make dh<0)
h_cor(i) = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0
dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness
total_depth(i,j) = total_depth(i,j) + dh
endif ; enddo
enddo
enddo
! A much simpler (but answer changing) version of the total_depth calculation would be
! do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
! total_depth(i,j) = total_depth(i,j) + dz(i,j,k)
! enddo ; enddo ; enddo

! Update halos once, then march inward for each iteration
if (CS%n_smooth > 1) call pass_var(total_depth, G%Domain, halo=CS%n_smooth, complete=.false.)
call pass_var(CS%OBLdepth, G%Domain, halo=CS%n_smooth)

if (CS%id_OBLdepth_original > 0) CS%OBLdepth_original = CS%OBLdepth
if (CS%id_OBLdepth_original > 0) CS%OBLdepth_original(:,:) = CS%OBLdepth(:,:)

do s=1,CS%n_smooth

OBLdepth_prev = CS%OBLdepth
OBLdepth_prev(:,:) = CS%OBLdepth(:,:)
halo = CS%n_smooth - s

! apply smoothing on OBL depth
!$OMP parallel do default(none) shared(G, GV, US, CS, dz, OBLdepth_prev) &
!$OMP private(wc, ww, we, wn, ws, dh, hcorr, cellHeight, iFaceHeight)
do j = G%jsc, G%jec
do i = G%isc, G%iec ; if (G%mask2dT(i,j) > 0.0) then

iFaceHeight(1) = 0.0 ! BBL is all relative to the surface
hcorr = 0.
do k=1,GV%ke

! cell center and cell bottom in meters (negative values in the ocean)
dh = dz(i,j,k) ! Nominal thickness to use for increment
dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0
dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness
cellHeight(k) = iFaceHeight(k) - 0.5 * dh
iFaceHeight(k+1) = iFaceHeight(k) - dh
enddo

!$OMP parallel do default(none) shared(G, GV, CS, OBLdepth_prev, total_depth, halo) &
!$OMP private(wc, ww, we, wn, ws)
do j = G%jsc-halo, G%jec+halo
do i = G%isc-halo, G%iec+halo ; if (G%mask2dT(i,j) > 0.0) then
! compute weights
ww = 0.125 * G%mask2dT(i-1,j)
we = 0.125 * G%mask2dT(i+1,j)
Expand All @@ -1400,19 +1413,37 @@ subroutine KPP_smooth_BLD(CS, G, GV, US, dz)
if (CS%deepen_only) CS%OBLdepth(i,j) = max(CS%OBLdepth(i,j), OBLdepth_prev(i,j))

! prevent OBL depths deeper than the bathymetric depth
CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(GV%ke+1) ) ! no deeper than bottom
CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )
CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), total_depth(i,j) ) ! no deeper than bottom
endif ; enddo
enddo

enddo ! s-loop

! Determine the fractional index of the bottom of the boundary layer.
!$OMP parallel do default(none) shared(G, GV, CS, dz) &
!$OMP private(dh, hcorr, cellHeight, iFaceHeight)
do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (G%mask2dT(i,j) > 0.0) then

iFaceHeight(1) = 0.0 ! BBL is all relative to the surface
hcorr = 0.
do k=1,GV%ke
! cell center and cell bottom in meters (negative values in the ocean)
dh = dz(i,j,k) ! Nominal thickness to use for increment
dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0
dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness
cellHeight(k) = iFaceHeight(k) - 0.5 * dh
iFaceHeight(k+1) = iFaceHeight(k) - dh
enddo

CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )
endif ; enddo ; enddo

call cpu_clock_end(id_clock_KPP_smoothing)

end subroutine KPP_smooth_BLD



!> Copies KPP surface boundary layer depth into BLD, in units of [Z ~> m] unless other units are specified.
subroutine KPP_get_BLD(CS, BLD, G, US, m_to_BLD_units)
type(KPP_CS), pointer :: CS !< Control structure for
Expand Down

0 comments on commit 920093a

Please sign in to comment.