From 502d710fa8b67efea2eadfdd02432b7901c3fba0 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 29 Jul 2023 05:31:31 -0400 Subject: [PATCH] +Add segment%dZtot Add the new element dZtot to the OBC_segment_type to hold the total vertical extent of the water column, and use thickness_to_dz in update_OBC_segment_data to convert the layer thicknesses to the vertical layer extents used to set dZtot. This change leads to the cancellation of 6 unit conversion factors. All answers are bitwise identical in Boussinesq mode, but they will change and become less dependent on the Boussinesq reference density in some Boussinesq cases with certain types of open boundary condition. --- src/core/MOM_open_boundary.F90 | 31 +++++++++++++++++++++++-------- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index fda119768a..1e6c057e63 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -14,6 +14,7 @@ module MOM_open_boundary use MOM_file_parser, only : get_param, log_version, param_file_type, log_param use MOM_grid, only : ocean_grid_type, hor_index_type use MOM_dyn_horgrid, only : dyn_horgrid_type +use MOM_interface_heights, only : thickness_to_dz use MOM_io, only : slasher, field_size, SINGLE_FILE use MOM_io, only : vardesc, query_vardesc, var_desc use MOM_restart, only : register_restart_field, register_restart_pair @@ -189,6 +190,7 @@ module MOM_open_boundary real, allocatable :: Cg(:,:) !< The external gravity wave speed [L T-1 ~> m s-1] !! at OBC-points. real, allocatable :: Htot(:,:) !< The total column thickness [H ~> m or kg m-2] at OBC-points. + real, allocatable :: dZtot(:,:) !< The total column vertical extent [Z ~> m] at OBC-points. real, allocatable :: h(:,:,:) !< The cell thickness [H ~> m or kg m-2] at OBC-points. real, allocatable :: normal_vel(:,:,:) !< The layer velocity normal to the OB !! segment [L T-1 ~> m s-1]. @@ -352,7 +354,7 @@ module MOM_open_boundary real, allocatable :: tres_y(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] logical :: debug !< If true, write verbose checksums for debugging purposes. real :: silly_h !< A silly value of thickness outside of the domain that can be used to test - !! the independence of the OBCs to this external data [H ~> m or kg m-2]. + !! the independence of the OBCs to this external data [Z ~> m]. real :: silly_u !< A silly value of velocity outside of the domain that can be used to test !! the independence of the OBCs to this external data [L T-1 ~> m s-1]. logical :: ramp = .false. !< If True, ramp from zero to the external values for SSH. @@ -3595,6 +3597,7 @@ subroutine allocate_OBC_segment_data(OBC, segment) ! If these are just Flather, change update_OBC_segment_data accordingly allocate(segment%Cg(IsdB:IedB,jsd:jed), source=0.0) allocate(segment%Htot(IsdB:IedB,jsd:jed), source=0.0) + allocate(segment%dZtot(IsdB:IedB,jsd:jed), source=0.0) allocate(segment%h(IsdB:IedB,jsd:jed,OBC%ke), source=0.0) allocate(segment%SSH(IsdB:IedB,jsd:jed), source=0.0) if (segment%radiation) & @@ -3630,6 +3633,7 @@ subroutine allocate_OBC_segment_data(OBC, segment) ! If these are just Flather, change update_OBC_segment_data accordingly allocate(segment%Cg(isd:ied,JsdB:JedB), source=0.0) allocate(segment%Htot(isd:ied,JsdB:JedB), source=0.0) + allocate(segment%dZtot(isd:ied,JsdB:JedB), source=0.0) allocate(segment%h(isd:ied,JsdB:JedB,OBC%ke), source=0.0) allocate(segment%SSH(isd:ied,JsdB:JedB), source=0.0) if (segment%radiation) & @@ -3671,6 +3675,7 @@ subroutine deallocate_OBC_segment_data(segment) if (allocated(segment%Cg)) deallocate(segment%Cg) if (allocated(segment%Htot)) deallocate(segment%Htot) + if (allocated(segment%dZtot)) deallocate(segment%dZtot) if (allocated(segment%h)) deallocate(segment%h) if (allocated(segment%SSH)) deallocate(segment%SSH) if (allocated(segment%rx_norm_rad)) deallocate(segment%rx_norm_rad) @@ -3753,7 +3758,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 do n = 1, OBC%number_of_segments do k = 1, GV%ke @@ -3804,6 +3809,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) integer :: ni_buf, nj_buf ! Number of filled values in tmp_buffer integer :: is_obc, ie_obc, js_obc, je_obc ! segment indices within local domain 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] integer :: is_obc2, js_obc2 @@ -3837,6 +3843,11 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10 endif + if (OBC%number_of_segments >= 1) then + call thickness_to_dz(h, tv, dz, G, GV, US) + call pass_var(dz, G%Domain) + endif + do n = 1, OBC%number_of_segments segment => OBC%segment(n) @@ -3869,11 +3880,13 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) I=segment%HI%IsdB do j=segment%HI%jsd,segment%HI%jed segment%Htot(I,j) = 0.0 + segment%dZtot(I,j) = 0.0 do k=1,GV%ke segment%h(I,j,k) = h(i+ishift,j,k) segment%Htot(I,j) = segment%Htot(I,j) + segment%h(I,j,k) + segment%dZtot(I,j) = segment%dZtot(I,j) + dz(i+ishift,j,k) enddo - segment%Cg(I,j) = sqrt(GV%g_prime(1)*segment%Htot(I,j)*GV%H_to_Z) + segment%Cg(I,j) = sqrt(GV%g_prime(1) * segment%dZtot(I,j)) enddo else! (segment%direction == OBC_DIRECTION_N .or. segment%direction == OBC_DIRECTION_S) allocate(normal_trans_bt(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB), source=0.0) @@ -3881,11 +3894,13 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) J=segment%HI%JsdB do i=segment%HI%isd,segment%HI%ied segment%Htot(i,J) = 0.0 + segment%dZtot(i,J) = 0.0 do k=1,GV%ke segment%h(i,J,k) = h(i,j+jshift,k) segment%Htot(i,J) = segment%Htot(i,J) + segment%h(i,J,k) + segment%dZtot(i,J) = segment%dZtot(i,J) + dz(i,j+jshift,k) enddo - segment%Cg(i,J) = sqrt(GV%g_prime(1)*segment%Htot(i,J)*GV%H_to_Z) + segment%Cg(i,J) = sqrt(GV%g_prime(1) * segment%dZtot(i,J)) enddo endif @@ -5638,8 +5653,8 @@ subroutine adjustSegmentEtaToFitBathymetry(G, GV, US, segment,fld) ! The normal slope at the boundary is zero by a ! previous call to open_boundary_impose_normal_slope do k=nz+1,1,-1 - if (-eta(i,j,k) > segment%Htot(i,j)*GV%H_to_Z + hTolerance) then - eta(i,j,k) = -segment%Htot(i,j)*GV%H_to_Z + if (-eta(i,j,k) > segment%dZtot(i,j) + hTolerance) then + eta(i,j,k) = -segment%dZtot(i,j) contractions = contractions + 1 endif enddo @@ -5657,10 +5672,10 @@ subroutine adjustSegmentEtaToFitBathymetry(G, GV, US, segment,fld) ! The whole column is dilated to accommodate deeper topography than ! the bathymetry would indicate. - if (-eta(i,j,nz+1) < (segment%Htot(i,j) * GV%H_to_Z) - hTolerance) then + if (-eta(i,j,nz+1) < segment%dZtot(i,j) - hTolerance) then dilations = dilations + 1 ! expand bottom-most cell only - eta(i,j,nz+1) = -(segment%Htot(i,j) * GV%H_to_Z) + eta(i,j,nz+1) = -segment%dZtot(i,j) segment%field(fld)%dz_src(i,j,nz)= eta(i,j,nz)-eta(i,j,nz+1) ! if (eta(i,j,1) <= eta(i,j,nz+1)) then ! do k=1,nz ; segment%field(fld)%dz_src(i,j,k) = (eta(i,j,1) + G%bathyT(i,j)) / real(nz) ; enddo