diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 7ee90d746a..059e515051 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2883,6 +2883,26 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & ! reservoirs are used. call open_boundary_register_restarts(HI, GV, US, CS%OBC, CS%tracer_Reg, & param_file, restart_CSp, use_temperature) + if (turns /= 0) then + if (CS%OBC%radiation_BCs_exist_globally) then + OBC_in%rx_normal => CS%OBC%rx_normal + OBC_in%ry_normal => CS%OBC%ry_normal + endif + if (CS%OBC%oblique_BCs_exist_globally) then + OBC_in%rx_oblique_u => CS%OBC%rx_oblique_u + OBC_in%ry_oblique_u => CS%OBC%ry_oblique_u + OBC_in%rx_oblique_v => CS%OBC%rx_oblique_v + OBC_in%ry_oblique_v => CS%OBC%ry_oblique_v + OBC_in%cff_normal_u => CS%OBC%cff_normal_u + OBC_in%cff_normal_v => CS%OBC%cff_normal_v + endif + if (any(CS%OBC%tracer_x_reservoirs_used)) then + OBC_in%tres_x => CS%OBC%tres_x + endif + if (any(CS%OBC%tracer_y_reservoirs_used)) then + OBC_in%tres_y => CS%OBC%tres_y + endif + endif endif if (present(waves_CSp)) then diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index af2beca1fb..0f57335322 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -742,6 +742,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB integer :: ioff, joff integer :: l_seg + real :: factor(SZI_(G),SZJ_(G)) ! If non-zero, work on given points. if (.not.CS%module_is_initialized) call MOM_error(FATAL, & "btstep: Module MOM_barotropic must be initialized before it is used.") @@ -2457,17 +2458,52 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, haloshift=iev-ie, unscale=US%L_to_m**2*GV%H_to_m) endif + do j=jsv,jev + do i=isv,iev + factor(i,j) = CS%IareaT(i,j) + enddo + enddo + + ! Update factor so that nothing changes outside of the OBC (problem for interior OBCs only) + if (associated(OBC)) then ; if (OBC%OBC_pe) then + do j=jsv,jev + if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then + do i=isv,iev-1 ; if (OBC%segnum_u(I,j) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then + factor(i+1,j) = 0.0 + elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + factor(i,j) = 0.0 + endif + endif ; enddo + endif + if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then + do i=isv,iev + if (OBC%segnum_v(i,J-1) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then + factor(i,j) = 0.0 + endif + endif + if (OBC%segnum_v(i,J) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + factor(i,j) = 0.0 + endif + endif + enddo + endif + enddo + endif ; endif + if (integral_BT_cont) then !$OMP do do j=jsv,jev ; do i=isv,iev - eta(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + CS%IareaT(i,j) * & + eta(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + factor(i,j) * & ((uhbt_int(I-1,j) - uhbt_int(I,j)) + (vhbt_int(i,J-1) - vhbt_int(i,J))) eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n) enddo ; enddo else !$OMP do do j=jsv,jev ; do i=isv,iev - eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * & + eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * factor(i,j)) * & ((uhbt(I-1,j) - uhbt(I,j)) + (vhbt(i,J-1) - vhbt(i,J))) eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n) enddo ; enddo diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index f89c8953ab..bf2f7cb2e9 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -352,24 +352,24 @@ module MOM_open_boundary type(remapping_CS), pointer :: remap_h_CS=> NULL() !< ALE remapping control structure for !! thickness-based fields on segments type(OBC_registry_type), pointer :: OBC_Reg => NULL() !< Registry type for boundaries - real, allocatable :: rx_normal(:,:,:) !< Array storage for normal phase speed for EW radiation OBCs in units of + real, pointer :: rx_normal(:,:,:) !< Array storage for normal phase speed for EW radiation OBCs in units of !! grid points per timestep [nondim] - real, allocatable :: ry_normal(:,:,:) !< Array storage for normal phase speed for NS radiation OBCs in units of + real, pointer :: ry_normal(:,:,:) !< Array storage for normal phase speed for NS radiation OBCs in units of !! grid points per timestep [nondim] - real, allocatable :: rx_oblique_u(:,:,:) !< X-direction oblique boundary condition radiation speeds squared + real, pointer :: rx_oblique_u(:,:,:) !< X-direction oblique boundary condition radiation speeds squared !! at u points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: ry_oblique_u(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared + real, pointer :: ry_oblique_u(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared !! at u points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: rx_oblique_v(:,:,:) !< X-direction oblique boundary condition radiation speeds squared + real, pointer :: rx_oblique_v(:,:,:) !< X-direction oblique boundary condition radiation speeds squared !! at v points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: ry_oblique_v(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared + real, pointer :: ry_oblique_v(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared !! at v points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: cff_normal_u(:,:,:) !< Denominator for normalizing EW oblique boundary condition radiation + real, pointer :: cff_normal_u(:,:,:) !< Denominator for normalizing EW oblique boundary condition radiation !! rates at u points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: cff_normal_v(:,:,:) !< Denominator for normalizing NS oblique boundary condition radiation + real, pointer :: cff_normal_v(:,:,:) !< Denominator for normalizing NS oblique boundary condition radiation !! rates at v points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: tres_x(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] - real, allocatable :: tres_y(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] + real, pointer :: tres_x(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] + real, pointer :: 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 [Z ~> m]. @@ -1948,15 +1948,15 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CS) call create_group_pass(OBC%pass_oblique, OBC%cff_normal_u, OBC%cff_normal_v, G%Domain, To_All+Scalar_Pair) call do_group_pass(OBC%pass_oblique, G%Domain) endif - if (allocated(OBC%tres_x) .and. allocated(OBC%tres_y)) then + if (associated(OBC%tres_x) .and. associated(OBC%tres_y)) then do m=1,OBC%ntr call pass_vector(OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%Domain, To_All+Scalar_Pair) enddo - elseif (allocated(OBC%tres_x)) then + elseif (associated(OBC%tres_x)) then do m=1,OBC%ntr call pass_var(OBC%tres_x(:,:,:,m), G%Domain, position=EAST_FACE) enddo - elseif (allocated(OBC%tres_y)) then + elseif (associated(OBC%tres_y)) then do m=1,OBC%ntr call pass_var(OBC%tres_y(:,:,:,m), G%Domain, position=NORTH_FACE) enddo @@ -2001,16 +2001,16 @@ subroutine open_boundary_dealloc(OBC) if (allocated(OBC%segment)) deallocate(OBC%segment) if (allocated(OBC%segnum_u)) deallocate(OBC%segnum_u) if (allocated(OBC%segnum_v)) deallocate(OBC%segnum_v) - if (allocated(OBC%rx_normal)) deallocate(OBC%rx_normal) - if (allocated(OBC%ry_normal)) deallocate(OBC%ry_normal) - if (allocated(OBC%rx_oblique_u)) deallocate(OBC%rx_oblique_u) - if (allocated(OBC%ry_oblique_u)) deallocate(OBC%ry_oblique_u) - if (allocated(OBC%rx_oblique_v)) deallocate(OBC%rx_oblique_v) - if (allocated(OBC%ry_oblique_v)) deallocate(OBC%ry_oblique_v) - if (allocated(OBC%cff_normal_u)) deallocate(OBC%cff_normal_u) - if (allocated(OBC%cff_normal_v)) deallocate(OBC%cff_normal_v) - if (allocated(OBC%tres_x)) deallocate(OBC%tres_x) - if (allocated(OBC%tres_y)) deallocate(OBC%tres_y) + if (associated(OBC%rx_normal)) nullify(OBC%rx_normal) + if (associated(OBC%ry_normal)) nullify(OBC%ry_normal) + if (associated(OBC%rx_oblique_u)) nullify(OBC%rx_oblique_u) + if (associated(OBC%ry_oblique_u)) nullify(OBC%ry_oblique_u) + if (associated(OBC%rx_oblique_v)) nullify(OBC%rx_oblique_v) + if (associated(OBC%ry_oblique_v)) nullify(OBC%ry_oblique_v) + if (associated(OBC%cff_normal_u)) nullify(OBC%cff_normal_u) + if (associated(OBC%cff_normal_v)) nullify(OBC%cff_normal_v) + if (associated(OBC%tres_x)) nullify(OBC%tres_x) + if (associated(OBC%tres_y)) nullify(OBC%tres_y) if (associated(OBC%remap_z_CS)) deallocate(OBC%remap_z_CS) if (associated(OBC%remap_h_CS)) deallocate(OBC%remap_h_CS) deallocate(OBC) @@ -3369,7 +3369,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, haloshift=0, symmetric=sym, unscale=1.0/US%L_T_to_m_s**2) endif if (OBC%ntr == 0) return - if (.not. allocated (OBC%tres_x) .or. .not. allocated (OBC%tres_y)) return + if (.not. associated (OBC%tres_x) .or. .not. associated (OBC%tres_y)) return do m=1,OBC%ntr write(var_num,'(I3.3)') m call uvchksum("radiation_OBCs: OBC%tres_[xy]_"//var_num, OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%HI, & @@ -5053,7 +5053,9 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) integer :: i, j integer :: l_seg logical :: fatal_error = .False. - real :: min_depth ! The minimum depth for ocean points [Z ~> m] + real :: min_depth ! The minimum depth for ocean points [Z ~> m] + real :: mask_depth ! The masking depth for ocean points [Z ~> m] + real :: Dmask ! The depth for masking in the same units as G%bathyT [Z ~> m]. integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2 character(len=256) :: mesg ! Message for error messages. real, allocatable, dimension(:,:) :: color, color2 ! For sorting inside from outside, @@ -5063,6 +5065,12 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & units="m", default=0.0, scale=US%m_to_Z, do_not_log=.true.) + call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, & + units="m", default=-9999.0, scale=US%m_to_Z, do_not_log=.true.) + + Dmask = mask_depth + if (mask_depth == -9999.0*US%m_to_Z) Dmask = min_depth + ! The reference depth on a dyn_horgrid is 0, otherwise would need: min_depth = min_depth - G%Z_ref allocate(color(G%isd:G%ied, G%jsd:G%jed), source=0.0) @@ -5153,7 +5161,7 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) &"the masking of the outside grid points.")') i, j call MOM_error(WARNING,"MOM mask_outside_OBCs: "//mesg, all_print=.true.) endif - if (color(i,j) == cout) G%bathyT(i,j) = min_depth + if (color(i,j) == cout) G%bathyT(i,j) = Dmask enddo ; enddo if (fatal_error) call MOM_error(FATAL, & "MOM_open_boundary: inconsistent OBC segments.") @@ -5481,7 +5489,7 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) ((1.0-a_out+a_in)*segment%tr_Reg%Tr(m)%tres(I,j,k)+ & ((u_L_out+a_out)*Reg%Tr(ntr_id)%t(I+ishift,j,k) - & (u_L_in+a_in)*segment%tr_Reg%Tr(m)%t(I,j,k))) - if (allocated(OBC%tres_x)) OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k) + if (associated(OBC%tres_x)) OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k) enddo ; endif enddo enddo @@ -5521,7 +5529,7 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) ((1.0-a_out+a_in)*segment%tr_Reg%Tr(m)%tres(i,J,k) + & ((v_L_out+a_out)*Reg%Tr(ntr_id)%t(i,J+jshift,k) - & (v_L_in+a_in)*segment%tr_Reg%Tr(m)%t(i,J,k))) - if (allocated(OBC%tres_y)) OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) + if (associated(OBC%tres_y)) OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) enddo ; endif enddo enddo @@ -5597,7 +5605,7 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) ! Update tracer concentrations segment%tr_Reg%Tr(m)%tres(I,j,:) = tr_column(:) - if (allocated(OBC%tres_x)) then ; do k=1,nz + if (associated(OBC%tres_x)) then ; do k=1,nz OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k) enddo ; endif @@ -5664,7 +5672,7 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) ! Update tracer concentrations segment%tr_Reg%Tr(m)%tres(i,J,:) = tr_column(:) - if (allocated(OBC%tres_y)) then ; do k=1,nz + if (associated(OBC%tres_y)) then ; do k=1,nz OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) enddo ; endif @@ -6109,6 +6117,14 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%field(n)%buffer_src) endif + if (allocated(segment_in%field(n)%buffer_dst)) then + call allocate_rotated_array(segment_in%field(n)%buffer_dst, & + lbound(segment_in%field(n)%buffer_dst), turns, & + segment%field(n)%buffer_dst) + call rotate_array(segment_in%field(n)%buffer_dst, turns, & + segment%field(n)%buffer_dst) + endif + segment%field(n)%nk_src = segment_in%field(n)%nk_src if (allocated(segment_in%field(n)%dz_src)) then @@ -6122,6 +6138,8 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%field(n)%value = segment_in%field(n)%value enddo + call rotate_array(segment_in%SSH, turns, segment%SSH) + segment%temp_segment_data_exists = segment_in%temp_segment_data_exists segment%salt_segment_data_exists = segment_in%salt_segment_data_exists end subroutine rotate_OBC_segment_data diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index e8fdc3a2f0..70fe412b7a 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -1070,13 +1070,18 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only) if (associated(OBC)) then ; if (OBC%OBC_pe) then if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then - do i=is,ie-1 ; if (OBC%segnum_v(i,J) /= OBC_NONE) then - if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then - do_i(i,j+1) = .false. - elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then - do_i(i,j) = .false. + do i=is,ie + if (OBC%segnum_v(i,J-1) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then + do_i(i,j) = .false. + endif endif - endif ; enddo + if (OBC%segnum_v(i,J) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + do_i(i,j) = .false. + endif + endif + enddo endif endif ; endif