Skip to content

Commit

Permalink
(*)non-Boussinesq Z unit update_OBC_segment_data
Browse files Browse the repository at this point in the history 
  Revised the update_OBC_segment_data code to keep the z-space input data in
height units ([Z ~> m]) rather than rescaling them to thickness units.  This
change means that the non-Boussinesq open boundary condition calculations avoid
using the Boussinesq reference density.  Associated with this change, 5 internal
variables in update_OBC_segment_data were renamed to reflect that they are layer
vertical extents instead of thicknesses.  Answers are bitwise identical in any
Boussinesq cases.  However, answers will change in any non-Boussinesq cases that
use the ALE_sponges.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Jan 24, 2024
1 parent d7d126a commit 9008c83
Showing 1 changed file with 50 additions and 51 deletions.
101 changes: 50 additions & 51 deletions src/core/MOM_open_boundary.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ module MOM_open_boundary
!! have the same units as the field they are being applied to?
integer :: nk_src !< Number of vertical levels in the source data
real, allocatable :: dz_src(:,:,:) !< vertical grid cell spacing of the incoming segment
!! data, set in [Z ~> m] then scaled to [H ~> m or kg m-2]
!! data in [Z ~> m].
real, allocatable :: buffer_dst(:,:,:) !< buffer src data remapped to the target vertical grid.
!! The values for tracers should have the same units as the field
!! they are being applied to?
Expand Down Expand Up @@ -3791,7 +3791,7 @@ subroutine open_boundary_test_extern_h(G, GV, OBC, h)

if (.not. associated(OBC)) return

silly_h = GV%Z_to_H * OBC%silly_h
silly_h = GV%Z_to_H * OBC%silly_h ! This rescaling is here because GV was initialized after OBC.

do n = 1, OBC%number_of_segments
do k = 1, GV%ke
Expand Down Expand Up @@ -3844,18 +3844,18 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
integer :: ishift, jshift ! offsets for staggered locations
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Distance between the interfaces around a layer [Z ~> m]
real, dimension(:,:,:), allocatable, target :: tmp_buffer ! A buffer for input data [various units]
real, dimension(:), allocatable :: h_stack ! Thicknesses at corner points [H ~> m or kg m-2]
real, dimension(:), allocatable :: dz_stack ! Distance between the interfaces at corner points [Z ~> m]
integer :: is_obc2, js_obc2
integer :: i_seg_offset, j_seg_offset
real :: net_H_src ! Total thickness of the incoming flow in the source field [H ~> m or kg m-2]
real :: net_H_int ! Total thickness of the incoming flow in the model [H ~> m or kg m-2]
real :: net_dz_src ! Total vertical extent of the incoming flow in the source field [Z ~> m]
real :: net_dz_int ! Total vertical extent of the incoming flow in the model [Z ~> m]
real :: scl_fac ! A scaling factor to compensate for differences in total thicknesses [nondim]
real :: tidal_vel ! Interpolated tidal velocity at the OBC points [L T-1 ~> m s-1]
real :: tidal_elev ! Interpolated tidal elevation at the OBC points [Z ~> m]
real, allocatable :: normal_trans_bt(:,:) ! barotropic transport [H L2 T-1 ~> m3 s-1]
integer :: turns ! Number of index quarter turns
real :: time_delta ! Time since tidal reference date [T ~> s]
real :: h_neglect, h_neglect_edge ! Small thicknesses [H ~> m or kg m-2]
real :: dz_neglect, dz_neglect_edge ! Small thicknesses [Z ~> m]

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
Expand All @@ -3869,11 +3869,13 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
if (OBC%add_tide_constituents) time_delta = US%s_to_T * time_type_to_real(Time - OBC%time_ref)

if (OBC%remap_answer_date >= 20190101) then
h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff
elseif (GV%Boussinesq) then
h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10
dz_neglect = US%m_to_Z * 1.0e-30 ; dz_neglect_edge = US%m_to_Z * 1.0e-10
elseif (GV%semi_Boussinesq) then
dz_neglect = GV%kg_m2_to_H*GV%H_to_Z * 1.0e-30 ; dz_neglect_edge = GV%kg_m2_to_H*GV%H_to_Z * 1.0e-10
else
h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10
dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff
endif

if (OBC%number_of_segments >= 1) then
Expand Down Expand Up @@ -3937,7 +3939,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
enddo
endif

allocate(h_stack(GV%ke), source=0.0)
allocate(dz_stack(GV%ke), source=0.0)
do m = 1,segment%num_fields
!This field may not require a high frequency OBC segment update and might be allowed
!a less frequent update as set by the parameter update_OBC_period_max in MOM.F90.
Expand Down Expand Up @@ -4149,6 +4151,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
endif
endif

! The units of ...%dz_src are no longer changed from [Z ~> m] to [H ~> m or kg m-2] here.
call adjustSegmentEtaToFitBathymetry(G,GV,US,segment,m)

if (segment%is_E_or_W) then
Expand All @@ -4162,26 +4165,26 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
! Pretty sure we need to check for source/target grid consistency here
segment%field(m)%buffer_dst(I,J,:)=0.0 ! initialize remap destination buffer
if (G%mask2dCu(I,j)>0. .and. G%mask2dCu(I,j+1)>0.) then
h_stack(:) = 0.5*(h(i+ishift,j,:) + h(i+ishift,j+1,:))
dz_stack(:) = 0.5*(dz(i+ishift,j,:) + dz(i+ishift,j+1,:))
call remapping_core_h(OBC%remap_CS, &
segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), &
segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), &
segment%field(m)%buffer_src(I,J,:), &
GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), &
h_neglect, h_neglect_edge)
GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), &
dz_neglect, dz_neglect_edge)
elseif (G%mask2dCu(I,j)>0.) then
h_stack(:) = h(i+ishift,j,:)
dz_stack(:) = dz(i+ishift,j,:)
call remapping_core_h(OBC%remap_CS, &
segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), &
segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), &
segment%field(m)%buffer_src(I,J,:), &
GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), &
h_neglect, h_neglect_edge)
GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), &
dz_neglect, dz_neglect_edge)
elseif (G%mask2dCu(I,j+1)>0.) then
h_stack(:) = h(i+ishift,j+1,:)
dz_stack(:) = dz(i+ishift,j+1,:)
call remapping_core_h(OBC%remap_CS, &
segment%field(m)%nk_src,segment%field(m)%dz_src(I,j,:), &
segment%field(m)%nk_src, segment%field(m)%dz_src(I,j,:), &
segment%field(m)%buffer_src(I,J,:), &
GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), &
h_neglect, h_neglect_edge)
GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), &
dz_neglect, dz_neglect_edge)
endif
enddo
else
Expand All @@ -4190,14 +4193,14 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
! Pretty sure we need to check for source/target grid consistency here
segment%field(m)%buffer_dst(I,j,:)=0.0 ! initialize remap destination buffer
if (G%mask2dCu(I,j)>0.) then
net_H_src = sum( segment%field(m)%dz_src(I,j,:) )
net_H_int = sum( h(i+ishift,j,:) )
scl_fac = net_H_int / net_H_src
net_dz_src = sum( segment%field(m)%dz_src(I,j,:) )
net_dz_int = sum( dz(i+ishift,j,:) )
scl_fac = net_dz_int / net_dz_src
call remapping_core_h(OBC%remap_CS, &
segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(I,j,:), &
segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(I,j,:), &
segment%field(m)%buffer_src(I,j,:), &
GV%ke, h(i+ishift,j,:), segment%field(m)%buffer_dst(I,j,:), &
h_neglect, h_neglect_edge)
GV%ke, dz(i+ishift,j,:), segment%field(m)%buffer_dst(I,j,:), &
dz_neglect, dz_neglect_edge)
endif
enddo
endif
Expand All @@ -4212,26 +4215,26 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
if (G%mask2dCv(i,J)>0. .and. G%mask2dCv(i+1,J)>0.) then
! Using the h remapping approach
! Pretty sure we need to check for source/target grid consistency here
h_stack(:) = 0.5*(h(i,j+jshift,:) + h(i+1,j+jshift,:))
dz_stack(:) = 0.5*(dz(i,j+jshift,:) + dz(i+1,j+jshift,:))
call remapping_core_h(OBC%remap_CS, &
segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), &
segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), &
segment%field(m)%buffer_src(I,J,:), &
GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), &
h_neglect, h_neglect_edge)
GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), &
dz_neglect, dz_neglect_edge)
elseif (G%mask2dCv(i,J)>0.) then
h_stack(:) = h(i,j+jshift,:)
dz_stack(:) = dz(i,j+jshift,:)
call remapping_core_h(OBC%remap_CS, &
segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), &
segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), &
segment%field(m)%buffer_src(I,J,:), &
GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), &
h_neglect, h_neglect_edge)
GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), &
dz_neglect, dz_neglect_edge)
elseif (G%mask2dCv(i+1,J)>0.) then
h_stack(:) = h(i+1,j+jshift,:)
dz_stack(:) = dz(i+1,j+jshift,:)
call remapping_core_h(OBC%remap_CS, &
segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), &
segment%field(m)%nk_src, segment%field(m)%dz_src(I,J,:), &
segment%field(m)%buffer_src(I,J,:), &
GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), &
h_neglect, h_neglect_edge)
GV%ke, dz_stack, segment%field(m)%buffer_dst(I,J,:), &
dz_neglect, dz_neglect_edge)
endif
enddo
else
Expand All @@ -4240,14 +4243,14 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
! Pretty sure we need to check for source/target grid consistency here
segment%field(m)%buffer_dst(i,J,:)=0.0 ! initialize remap destination buffer
if (G%mask2dCv(i,J)>0.) then
net_H_src = sum( segment%field(m)%dz_src(i,J,:) )
net_H_int = sum( h(i,j+jshift,:) )
scl_fac = net_H_int / net_H_src
net_dz_src = sum( segment%field(m)%dz_src(i,J,:) )
net_dz_int = sum( dz(i,j+jshift,:) )
scl_fac = net_dz_int / net_dz_src
call remapping_core_h(OBC%remap_CS, &
segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(i,J,:), &
segment%field(m)%nk_src, scl_fac* segment%field(m)%dz_src(i,J,:), &
segment%field(m)%buffer_src(i,J,:), &
GV%ke, h(i,j+jshift,:), segment%field(m)%buffer_dst(i,J,:), &
h_neglect, h_neglect_edge)
GV%ke, dz(i,j+jshift,:), segment%field(m)%buffer_dst(i,J,:), &
dz_neglect, dz_neglect_edge)
endif
enddo
endif
Expand Down Expand Up @@ -4516,7 +4519,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
endif

enddo ! end field loop
deallocate(h_stack)
deallocate(dz_stack)
deallocate(normal_trans_bt)

enddo ! end segment loop
Expand Down Expand Up @@ -5750,10 +5753,6 @@ subroutine adjustSegmentEtaToFitBathymetry(G, GV, US, segment,fld)
! endif
!do k=nz,2,-1 ; eta(i,j,K) = eta(i,j,K+1) + segment%field(fld)%dz_src(i,j,k) ; enddo
endif
! Now convert thicknesses to units of H.
do k=1,nz
segment%field(fld)%dz_src(i,j,k) = segment%field(fld)%dz_src(i,j,k)*GV%Z_to_H
enddo
enddo ; enddo

! can not do communication call here since only PEs on the current segment are here
Expand Down

0 comments on commit 9008c83

Please sign in to comment.