diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index daf2d448ba..e316a068f5 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2006,6 +2006,31 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, & e_anom(i,j) = dgeo_de * (0.5 * (eta(i,j) + eta_in(i,j)) - eta_PF(i,j)) enddo ; enddo endif + if (apply_OBCs) then + if (CS%BT_OBC%apply_u_OBCs) then ! copy back the value for u-points on the boundary. +!GOMP parallel do default(none) shared(is,ie,js,je,ubt_sum_prev,ubt_sum,uhbt_sum_prev,& +!GOMP uhbt_sum,ubt_wtd_prev,ubt_wtd) + do j=js,je ; do I=is-1,ie + if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then + e_anom(i+1,j) = e_anom(i,j) + elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + e_anom(i,j) = e_anom(i+1,j) + endif + enddo ; enddo + endif + + if (CS%BT_OBC%apply_v_OBCs) then ! copy back the value for v-points on the boundary. +!GOMP parallel do default(none) shared(is,ie,js,je,vbt_sum_prev,vbt_sum,vhbt_sum_prev, & +!GOMP vhbt_sum,vbt_wtd_prev,vbt_wtd) + do J=js-1,je ; do I=is,ie + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then + e_anom(i,j+1) = e_anom(i,j) + elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + e_anom(i,j) = e_anom(i,j+1) + endif + enddo ; enddo + endif + endif ! It is possible that eta_out and eta_in are the same. do j=js,je ; do i=is,ie diff --git a/src/core/MOM_continuity.F90 b/src/core/MOM_continuity.F90 index 33b3526ee7..a2a50594cd 100644 --- a/src/core/MOM_continuity.F90 +++ b/src/core/MOM_continuity.F90 @@ -45,7 +45,7 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, & type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity, in m/s. real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity, in m/s. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hin !< Initial layer thickness, in m or kg/m2. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: hin !< Initial layer thickness, in m or kg/m2. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Final layer thickness, in m or kg/m2. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: uh !< Volume flux through zonal faces = !! u*h*dy, in m3/s. diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 34e5154585..31dc507800 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -8,7 +8,7 @@ module MOM_continuity_PPM use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE +use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type, OBC_NONE use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_variables, only : BT_cont_type use MOM_verticalGrid, only : verticalGrid_type @@ -79,7 +79,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, type(continuity_PPM_CS), pointer :: CS !< Module's control structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity, in m s-1. real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity, in m s-1. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hin !< Initial layer thickness, in H. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: hin !< Initial layer thickness, in H. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Final layer thickness, in H. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: uh !< Zonal volume flux, !! u*h*dy, H m2 s-1. @@ -286,7 +286,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity, in m s-1. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_in !< Layer thickness used to + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_in !< Layer thickness used to !! calculate fluxes, in H. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: uh !< Volume flux through zonal !! faces = u*h*dy, H m2 s-1. @@ -331,18 +331,21 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & real :: I_dt ! 1.0 / dt, in s-1. real :: du_lim ! The velocity change that give a relative CFL of 1, in m s-1. real :: dx_E, dx_W ! Effective x-grid spacings to the east and west, in m. - integer :: i, j, k, ish, ieh, jsh, jeh, nz + integer :: i, j, k, ish, ieh, jsh, jeh, n, nz logical :: do_aux, local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC - logical :: local_Flather_OBC, is_simple + logical :: local_Flather_OBC, local_open_BC, is_simple + type(OBC_segment_type), pointer :: segment do_aux = (present(uhbt_aux) .and. present(u_cor_aux)) use_visc_rem = present(visc_rem_u) local_specified_BC = .false. ; set_BT_cont = .false. ; local_Flather_OBC = .false. + local_open_BC = .false. if (present(BT_cont)) set_BT_cont = (associated(BT_cont)) if (present(OBC)) then ; if (associated(OBC)) then local_specified_BC = OBC%specified_u_BCs_exist_globally local_Flather_OBC = OBC%Flather_u_BCs_exist_globally .or. & OBC%Flather_v_BCs_exist_globally + local_open_BC = OBC%open_u_BCs_exist_globally endif ; endif ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = G%ke @@ -364,6 +367,27 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & endif do I=ish-1,ieh ; visc_rem(I,k) = 1.0 ; enddo enddo + if (local_open_BC) then + do n=1, OBC%number_of_segments + segment => OBC%segment(n) + if (.not. segment%on_pe .or. segment%specified) cycle + if (segment%direction == OBC_DIRECTION_E) then + I=segment%HI%IsdB + do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed + h_in(i+1,j,k) = h_in(i,j,k) + h_L(i+1,j,k) = h_in(i,j,k) + h_R(i+1,j,k) = h_in(i,j,k) + enddo ; enddo + elseif (segment%direction == OBC_DIRECTION_W) then + I=segment%HI%IsdB + do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed + h_in(i,j,k) = h_in(i+1,j,k) + h_L(i,j,k) = h_in(i+1,j,k) + h_R(i,j,k) = h_in(i+1,j,k) + enddo ; enddo + endif + enddo + endif call cpu_clock_end(id_clock_update) call cpu_clock_begin(id_clock_correct) @@ -385,10 +409,12 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & call zonal_flux_layer(u(:,j,k), h_in(:,j,k), h_L(:,j,k), h_R(:,j,k), & uh(:,j,k), duhdu(:,k), visc_rem(:,k), & dt, G, j, ish, ieh, do_I, CS%vol_CFL) - if (local_specified_BC) then ; do I=ish-1,ieh - if (OBC%segment(OBC%segnum_u(I,j))%specified) & - uh(I,j,k) = OBC%segment(OBC%segnum_u(I,j))%normal_trans(I,j,k) - enddo ; endif + if (local_specified_BC) then + do I=ish-1,ieh + if (OBC%segment(OBC%segnum_u(I,j))%specified) & + uh(I,j,k) = OBC%segment(OBC%segnum_u(I,j))%normal_trans(I,j,k) + enddo + endif enddo if ((.not.use_visc_rem).or.(.not.CS%use_visc_rem_max)) then ; do I=ish-1,ieh @@ -1027,7 +1053,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Meridional velocity, in m s-1. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_in !< Layer thickness used to + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_in !< Layer thickness used to !! calculate fluxes, in H. real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: vh !< Volume flux through meridional !! faces = v*h*dx, H m2 s-1. @@ -1078,18 +1104,21 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & real :: I_dt ! 1.0 / dt, in s-1. real :: dv_lim ! The velocity change that give a relative CFL of 1, in m s-1. real :: dy_N, dy_S ! Effective y-grid spacings to the north and south, in m. - integer :: i, j, k, ish, ieh, jsh, jeh, nz + integer :: i, j, k, ish, ieh, jsh, jeh, n, nz logical :: do_aux, local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC - logical :: local_Flather_OBC, is_simple + logical :: local_Flather_OBC, is_simple, local_open_BC + type(OBC_segment_type), pointer :: segment do_aux = (present(vhbt_aux) .and. present(v_cor_aux)) use_visc_rem = present(visc_rem_v) local_specified_BC = .false. ; set_BT_cont = .false. ; local_Flather_OBC = .false. + local_open_BC = .false. if (present(BT_cont)) set_BT_cont = (associated(BT_cont)) if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then local_specified_BC = OBC%specified_v_BCs_exist_globally local_Flather_OBC = OBC%Flather_u_BCs_exist_globally .or. & OBC%Flather_v_BCs_exist_globally + local_open_BC = OBC%open_v_BCs_exist_globally endif ; endif ; endif ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = G%ke @@ -1111,6 +1140,27 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & endif do i=ish,ieh ; visc_rem(i,k) = 1.0 ; enddo enddo + if (local_open_BC) then + do n=1, OBC%number_of_segments + segment => OBC%segment(n) + if (.not. segment%on_pe .or. segment%specified) cycle + if (segment%direction == OBC_DIRECTION_N) then + J=segment%HI%JsdB + do k=1,nz ; do i=segment%HI%isd,segment%HI%ied + h_in(i,j+1,k) = h_in(i,j,k) + h_L(i,j+1,k) = h_in(i,j,k) + h_R(i,j+1,k) = h_in(i,j,k) + enddo ; enddo + elseif (segment%direction == OBC_DIRECTION_S) then + J=segment%HI%JsdB + do k=1,nz ; do i=segment%HI%isd,segment%HI%ied + h_in(i,j,k) = h_in(i,j+1,k) + h_L(i,j,k) = h_in(i,j+1,k) + h_R(i,j,k) = h_in(i,j+1,k) + enddo ; enddo + endif + enddo + endif call cpu_clock_end(id_clock_update) call cpu_clock_begin(id_clock_correct) @@ -1134,10 +1184,12 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & call merid_flux_layer(v(:,J,k), h_in(:,:,k), h_L(:,:,k), h_R(:,:,k), & vh(:,J,k), dvhdv(:,k), visc_rem(:,k), & dt, G, J, ish, ieh, do_I, CS%vol_CFL) - if (local_specified_BC) then ; do i=ish,ieh - if (OBC%segment(OBC%segnum_v(i,J))%specified) & - vh(i,J,k) = OBC%segment(OBC%segnum_v(i,J))%normal_trans(i,J,k) - enddo ; endif + if (local_specified_BC) then + do i=ish,ieh + if (OBC%segment(OBC%segnum_v(i,J))%specified) & + vh(i,J,k) = OBC%segment(OBC%segnum_v(i,J))%normal_trans(i,J,k) + enddo + endif enddo ! k-loop if ((.not.use_visc_rem) .or. (.not.CS%use_visc_rem_max)) then ; do i=ish,ieh visc_rem_max(i) = 1.0 diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index e643ca9505..90a07bab1f 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -1062,7 +1062,7 @@ subroutine adjustments_dyn_legacy_split(u, v, h, dt, G, GV, CS) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< The zonal velocity, in m s-1 real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< The meridional velocity, in m s-1 - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thicknesses, in H (usually m or kg m-2) real, intent(in) :: dt !< The baroclinic dynamics time step, in s type(MOM_dyn_legacy_split_CS), pointer :: CS diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index ae0280dbb6..f3ec01aadf 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -74,6 +74,7 @@ module MOM_open_boundary logical :: oblique !< Oblique waves supported at radiation boundary. logical :: nudged !< Optional supplement to radiation boundary. logical :: specified !< Boundary fixed to external value. + logical :: open !< Boundary is open for continuity solver. logical :: gradient !< Zero gradient at boundary. logical :: values_needed !< Whether or not external OBC fields are needed. integer :: direction !< Boundary faces one of the four directions. @@ -279,6 +280,7 @@ subroutine open_boundary_config(G, param_file, OBC) OBC%segment(l)%oblique = .false. OBC%segment(l)%nudged = .false. OBC%segment(l)%specified = .false. + OBC%segment(l)%open = .false. OBC%segment(l)%gradient = .false. OBC%segment(l)%values_needed = .false. OBC%segment(l)%direction = OBC_NONE @@ -601,13 +603,16 @@ subroutine setup_u_point_obc(OBC, G, segment_str, l_seg) cycle elseif (trim(action_str(a_loop)) == 'FLATHER') then OBC%segment(l_seg)%Flather = .true. + OBC%segment(l_seg)%open = .true. OBC%Flather_u_BCs_exist_globally = .true. OBC%open_u_BCs_exist_globally = .true. elseif (trim(action_str(a_loop)) == 'ORLANSKI') then OBC%segment(l_seg)%radiation = .true. + OBC%segment(l_seg)%open = .true. OBC%open_u_BCs_exist_globally = .true. elseif (trim(action_str(a_loop)) == 'OBLIQUE') then OBC%segment(l_seg)%oblique = .true. + OBC%segment(l_seg)%open = .true. OBC%oblique_BCs_exist_globally = .true. OBC%open_u_BCs_exist_globally = .true. elseif (trim(action_str(a_loop)) == 'NUDGED') then @@ -617,6 +622,7 @@ subroutine setup_u_point_obc(OBC, G, segment_str, l_seg) OBC%nudged_u_BCs_exist_globally = .true. elseif (trim(action_str(a_loop)) == 'GRADIENT') then OBC%segment(l_seg)%gradient = .true. + OBC%segment(l_seg)%open = .true. OBC%open_u_BCs_exist_globally = .true. elseif (trim(action_str(a_loop)) == 'LEGACY') then OBC%segment(l_seg)%Flather = .true. @@ -700,13 +706,16 @@ subroutine setup_v_point_obc(OBC, G, segment_str, l_seg) cycle elseif (trim(action_str(a_loop)) == 'FLATHER') then OBC%segment(l_seg)%Flather = .true. + OBC%segment(l_seg)%open = .true. OBC%Flather_v_BCs_exist_globally = .true. OBC%open_v_BCs_exist_globally = .true. elseif (trim(action_str(a_loop)) == 'ORLANSKI') then OBC%segment(l_seg)%radiation = .true. + OBC%segment(l_seg)%open = .true. OBC%open_v_BCs_exist_globally = .true. elseif (trim(action_str(a_loop)) == 'OBLIQUE') then OBC%segment(l_seg)%oblique = .true. + OBC%segment(l_seg)%open = .true. OBC%oblique_BCs_exist_globally = .true. OBC%open_v_BCs_exist_globally = .true. elseif (trim(action_str(a_loop)) == 'NUDGED') then @@ -716,6 +725,7 @@ subroutine setup_v_point_obc(OBC, G, segment_str, l_seg) OBC%nudged_v_BCs_exist_globally = .true. elseif (trim(action_str(a_loop)) == 'GRADIENT') then OBC%segment(l_seg)%gradient = .true. + OBC%segment(l_seg)%open = .true. OBC%open_v_BCs_exist_globally = .true. elseif (trim(action_str(a_loop)) == 'LEGACY') then OBC%segment(l_seg)%radiation = .true. diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index 4714bd7a15..00491f43cb 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/src/user/Kelvin_initialization.F90 @@ -220,6 +220,9 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, h, Time) if (segment%direction == OBC_DIRECTION_E) cycle if (segment%direction == OBC_DIRECTION_N) cycle + ! This should be somewhere else... + segment%Tnudge_in = 1.0/(0.3*86400) + if (segment%direction == OBC_DIRECTION_W) then IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB jsd = segment%HI%jsd ; jed = segment%HI%jed