diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index d55389620d..dd209fdb9f 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -363,9 +363,7 @@ module MOM_barotropic OBC_mask_v => NULL() integer, dimension(:,:), pointer :: & OBC_direction_u => NULL(), & - OBC_direction_v => NULL(), & - OBC_kind_u => NULL(), & - OBC_kind_v => NULL() + OBC_direction_v => NULL() real, dimension(:,:), pointer :: & Cg_u => NULL(), & ! The external wave speed at u-points, in m s-1. Cg_v => NULL(), & ! The external wave speed at u-points, in m s-1. @@ -2376,11 +2374,11 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, if (associated(BT_OBC%OBC_mask_u)) then do j=js,je ; do I=is-1,ie ; if (BT_OBC%OBC_mask_u(I,j)) then - if (BT_OBC%OBC_kind_u(I,j) == OBC_SIMPLE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified) then uhbt(I,j) = BT_OBC%uhbt(I,j) ubt(I,j) = BT_OBC%ubt_outer(I,j) vel_trans = ubt(I,j) - elseif (BT_OBC%OBC_kind_u(I,j) == OBC_FLATHER) then + elseif (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then if (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) then cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL u_inlet = cfl*ubt_old(I-1,j) + (1.0-cfl)*ubt_old(I,j) ! Valid for cfl<1 @@ -2422,7 +2420,7 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, endif endif - if (BT_OBC%OBC_kind_u(I,j) /= OBC_SIMPLE) then + if (.not. OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified) then if (use_BT_cont) then uhbt(I,j) = find_uhbt(vel_trans,BTCL_u(I,j)) + uhbt0(I,j) else @@ -2436,11 +2434,11 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, if (associated(BT_OBC%OBC_mask_v)) then do J=js-1,je ; do i=is,ie ; if (BT_OBC%OBC_mask_v(i,J)) then - if (BT_OBC%OBC_kind_v(i,J) == OBC_SIMPLE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%specified) then vhbt(i,J) = BT_OBC%vhbt(i,J) vbt(i,J) = BT_OBC%vbt_outer(i,J) vel_trans = vbt(i,J) - elseif (BT_OBC%OBC_kind_v(i,J) == OBC_FLATHER) then + elseif (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then if (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) then cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(I,j) ! CFL v_inlet = cfl*vbt_old(i,J-1) + (1.0-cfl)*vbt_old(i,J) ! Valid for cfl<1 @@ -2490,7 +2488,7 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, endif endif - if (BT_OBC%OBC_kind_v(i,J) /= OBC_SIMPLE) then + if (.not. OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%specified) then if (use_BT_cont) then vhbt(i,J) = find_vhbt(vel_trans,BTCL_v(i,J)) + vhbt0(i,J) else @@ -2538,50 +2536,54 @@ subroutine apply_eta_OBCs(OBC, eta, ubt, vbt, BT_OBC, G, MS, halo, dtbt) if ((OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west) .and. & associated(BT_OBC%OBC_mask_u)) then - do j=js,je ; do I=is-1,ie ; if (BT_OBC%OBC_kind_u(I,j) == OBC_FLATHER) then - if (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) then - cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL - u_inlet = cfl*ubt(I-1,j) + (1.0-cfl)*ubt(I,j) ! Valid for cfl <1 -! h_in = 2.0*cfl*eta(i,j) + (1.0-2.0*cfl)*eta(i+1,j) ! external - h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i-1,j)) ! internal - - H_u = BT_OBC%H_u(I,j) - eta(i+1,j) = 2.0 * 0.5*((BT_OBC%eta_outer_u(I,j)+h_in) + & - (H_u/BT_OBC%Cg_u(I,j))*(u_inlet-BT_OBC%ubt_outer(I,j))) - eta(i,j) - elseif (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W) then - cfl = dtbt*BT_OBC%Cg_u(I,j)*G%IdxCu(I,j) ! CFL - u_inlet = cfl*ubt(I+1,j) + (1.0-cfl)*ubt(I,j) ! Valid for cfl <1 -! h_in = 2.0*cfl*eta(i+1,j) + (1.0-2.0*cfl)*eta(i,j) ! external - h_in = eta(i+1,j) + (0.5-cfl)*(eta(i+1,j)-eta(i+2,j)) ! internal - - H_u = BT_OBC%H_u(I,j) - eta(i,j) = 2.0 * 0.5*((BT_OBC%eta_outer_u(I,j)+h_in) + & - (H_u/BT_OBC%Cg_u(I,j))*(BT_OBC%ubt_outer(I,j)-u_inlet)) - eta(i+1,j) + do j=js,je ; do I=is-1,ie ; if (OBC%OBC_segment_u(I,j) /= OBC_NONE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then + if (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) then + cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL + u_inlet = cfl*ubt(I-1,j) + (1.0-cfl)*ubt(I,j) ! Valid for cfl <1 +! h_in = 2.0*cfl*eta(i,j) + (1.0-2.0*cfl)*eta(i+1,j) ! external + h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i-1,j)) ! internal + + H_u = BT_OBC%H_u(I,j) + eta(i+1,j) = 2.0 * 0.5*((BT_OBC%eta_outer_u(I,j)+h_in) + & + (H_u/BT_OBC%Cg_u(I,j))*(u_inlet-BT_OBC%ubt_outer(I,j))) - eta(i,j) + elseif (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W) then + cfl = dtbt*BT_OBC%Cg_u(I,j)*G%IdxCu(I,j) ! CFL + u_inlet = cfl*ubt(I+1,j) + (1.0-cfl)*ubt(I,j) ! Valid for cfl <1 +! h_in = 2.0*cfl*eta(i+1,j) + (1.0-2.0*cfl)*eta(i,j) ! external + h_in = eta(i+1,j) + (0.5-cfl)*(eta(i+1,j)-eta(i+2,j)) ! internal + + H_u = BT_OBC%H_u(I,j) + eta(i,j) = 2.0 * 0.5*((BT_OBC%eta_outer_u(I,j)+h_in) + & + (H_u/BT_OBC%Cg_u(I,j))*(BT_OBC%ubt_outer(I,j)-u_inlet)) - eta(i+1,j) + endif endif endif ; enddo ; enddo endif if ((OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south) .and. & associated(BT_OBC%OBC_mask_v)) then - do J=js-1,je ; do i=is,ie ; if (BT_OBC%OBC_kind_v(i,J) == OBC_FLATHER) then - if (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) then - cfl = dtbt*BT_OBC%Cg_v(i,J)*G%IdyCv(i,J) ! CFL - v_inlet = cfl*vbt(i,J-1) + (1.0-cfl)*vbt(i,J) ! Valid for cfl <1 -! h_in = 2.0*cfl*eta(i,j) + (1.0-2.0*cfl)*eta(i,j+1) ! external - h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i,j-1)) ! internal - - H_v = BT_OBC%H_v(i,J) - eta(i,j+1) = 2.0 * 0.5*((BT_OBC%eta_outer_v(i,J)+h_in) + & - (H_v/BT_OBC%Cg_v(i,J))*(v_inlet-BT_OBC%vbt_outer(i,J))) - eta(i,j) - elseif (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) then - cfl = dtbt*BT_OBC%Cg_v(i,J)*G%IdyCv(i,J) ! CFL - v_inlet = cfl*vbt(i,J+1) + (1.0-cfl)*vbt(i,J) ! Valid for cfl <1 -! h_in = 2.0*cfl*eta(i,j+1) + (1.0-2.0*cfl)*eta(i,j) ! external - h_in = eta(i,j+1) + (0.5-cfl)*(eta(i,j+1)-eta(i,j+2)) ! internal - - H_v = BT_OBC%H_v(i,J) - eta(i,j) = 2.0 * 0.5*((BT_OBC%eta_outer_v(i,J)+h_in) + & - (H_v/BT_OBC%Cg_v(i,J))*(BT_OBC%vbt_outer(i,J)-v_inlet)) - eta(i,j+1) + do J=js-1,je ; do i=is,ie ; if (OBC%OBC_segment_v(i,J) /= OBC_NONE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then + if (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) then + cfl = dtbt*BT_OBC%Cg_v(i,J)*G%IdyCv(i,J) ! CFL + v_inlet = cfl*vbt(i,J-1) + (1.0-cfl)*vbt(i,J) ! Valid for cfl <1 +! h_in = 2.0*cfl*eta(i,j) + (1.0-2.0*cfl)*eta(i,j+1) ! external + h_in = eta(i,j) + (0.5-cfl)*(eta(i,j)-eta(i,j-1)) ! internal + + H_v = BT_OBC%H_v(i,J) + eta(i,j+1) = 2.0 * 0.5*((BT_OBC%eta_outer_v(i,J)+h_in) + & + (H_v/BT_OBC%Cg_v(i,J))*(v_inlet-BT_OBC%vbt_outer(i,J))) - eta(i,j) + elseif (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) then + cfl = dtbt*BT_OBC%Cg_v(i,J)*G%IdyCv(i,J) ! CFL + v_inlet = cfl*vbt(i,J+1) + (1.0-cfl)*vbt(i,J) ! Valid for cfl <1 +! h_in = 2.0*cfl*eta(i,j+1) + (1.0-2.0*cfl)*eta(i,j) ! external + h_in = eta(i,j+1) + (0.5-cfl)*(eta(i,j+1)-eta(i,j+2)) ! internal + + H_v = BT_OBC%H_v(i,J) + eta(i,j) = 2.0 * 0.5*((BT_OBC%eta_outer_v(i,J)+h_in) + & + (H_v/BT_OBC%Cg_v(i,J))*(BT_OBC%vbt_outer(i,J)-v_inlet)) - eta(i,j+1) + endif endif endif ; enddo ; enddo endif @@ -2640,7 +2642,6 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D allocate(BT_OBC%ubt_outer(isdw-1:iedw,jsdw:jedw)) ; BT_OBC%ubt_outer(:,:) = 0.0 allocate(BT_OBC%eta_outer_u(isdw-1:iedw,jsdw:jedw)) ; BT_OBC%eta_outer_u(:,:) = 0.0 allocate(BT_OBC%OBC_mask_u(isdw-1:iedw,jsdw:jedw)) ; BT_OBC%OBC_mask_u(:,:)=.false. - allocate(BT_OBC%OBC_kind_u(isdw-1:iedw,jsdw:jedw)) ; BT_OBC%OBC_kind_u(:,:)=OBC_NONE allocate(BT_OBC%OBC_direction_u(isdw-1:iedw,jsdw:jedw)); BT_OBC%OBC_direction_u(:,:)=OBC_NONE allocate(BT_OBC%Cg_v(isdw:iedw,jsdw-1:jedw)) ; BT_OBC%Cg_v(:,:) = 0.0 @@ -2649,13 +2650,11 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D allocate(BT_OBC%vbt_outer(isdw:iedw,jsdw-1:jedw)) ; BT_OBC%vbt_outer(:,:) = 0.0 allocate(BT_OBC%eta_outer_v(isdw:iedw,jsdw-1:jedw)) ; BT_OBC%eta_outer_v(:,:)=0.0 allocate(BT_OBC%OBC_mask_v(isdw:iedw,jsdw-1:jedw)) ; BT_OBC%OBC_mask_v(:,:)=.false. - allocate(BT_OBC%OBC_kind_v(isdw-1:iedw,jsdw:jedw)) ; BT_OBC%OBC_kind_v(:,:)=OBC_NONE allocate(BT_OBC%OBC_direction_v(isdw-1:iedw,jsdw:jedw)); BT_OBC%OBC_direction_v(:,:)=OBC_NONE if (associated(OBC%OBC_mask_u)) then do j=js-1,je+1 ; do I=is-1,ie BT_OBC%OBC_mask_u(I,j) = OBC%OBC_mask_u(I,j) - BT_OBC%OBC_kind_u(I,j) = OBC%OBC_kind_u(I,j) BT_OBC%OBC_direction_u(I,j) = OBC%OBC_direction_u(I,j) enddo ; enddo if (OBC%apply_OBC_u) then @@ -2664,7 +2663,7 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D enddo ; enddo ; enddo endif do j=js,je ; do I=is-1,ie ; if (OBC%OBC_mask_u(I,j)) then - if (OBC%OBC_kind_u(I,j) == OBC_SIMPLE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified) then if (use_BT_cont) then BT_OBC%ubt_outer(I,j) = uhbt_to_ubt(BT_OBC%uhbt(I,j),BTCL_u(I,j)) else @@ -2691,7 +2690,6 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D if (associated(OBC%OBC_mask_v)) then do J=js-1,je ; do i=is-1,ie+1 BT_OBC%OBC_mask_v(i,J) = OBC%OBC_mask_v(i,J) - BT_OBC%OBC_kind_v(i,J) = OBC%OBC_kind_v(i,J) BT_OBC%OBC_direction_v(i,J) = OBC%OBC_direction_v(i,J) enddo ; enddo if (OBC%apply_OBC_v) then @@ -2701,7 +2699,7 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D endif do J=js-1,je ; do i=is,ie ; if (OBC%OBC_mask_v(i,J)) then - if (OBC%OBC_kind_v(i,J) == OBC_SIMPLE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%specified) then if (use_BT_cont) then BT_OBC%vbt_outer(i,J) = vhbt_to_vbt(BT_OBC%vhbt(i,J),BTCL_v(i,J)) else @@ -2741,7 +2739,6 @@ subroutine destroy_BT_OBC(BT_OBC) type(BT_OBC_type), intent(inout) :: BT_OBC if (associated(BT_OBC%OBC_mask_u)) deallocate(BT_OBC%OBC_mask_u) - if (associated(BT_OBC%OBC_kind_u)) deallocate(BT_OBC%OBC_kind_u) if (associated(BT_OBC%OBC_direction_u)) deallocate(BT_OBC%OBC_direction_u) deallocate(BT_OBC%Cg_u) deallocate(BT_OBC%H_u) @@ -2750,7 +2747,6 @@ subroutine destroy_BT_OBC(BT_OBC) deallocate(BT_OBC%eta_outer_u) if (associated(BT_OBC%OBC_mask_v)) deallocate(BT_OBC%OBC_mask_v) - if (associated(BT_OBC%OBC_kind_v)) deallocate(BT_OBC%OBC_kind_v) if (associated(BT_OBC%OBC_direction_v)) deallocate(BT_OBC%OBC_direction_v) deallocate(BT_OBC%Cg_v) deallocate(BT_OBC%H_v) diff --git a/src/core/MOM_boundary_update.F90 b/src/core/MOM_boundary_update.F90 new file mode 100644 index 0000000000..8ad29f2cbc --- /dev/null +++ b/src/core/MOM_boundary_update.F90 @@ -0,0 +1,88 @@ +! This file is part of MOM6. See LICENSE.md for the license. +!> Controls where open boundary conditions are applied +module MOM_boundary_update + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE +use MOM_diag_mediator, only : diag_ctrl, time_type +use MOM_domains, only : pass_var, pass_vector +use MOM_domains, only : To_All, SCALAR_PAIR, CGRID_NE +use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING +use MOM_file_parser, only : get_param, log_version, param_file_type, log_param +use MOM_grid, only : ocean_grid_type +use MOM_dyn_horgrid, only : dyn_horgrid_type +use MOM_io, only : EAST_FACE, NORTH_FACE +use MOM_io, only : slasher, read_data +use MOM_open_boundary, only : ocean_obc_type +use MOM_tracer_registry, only : add_tracer_OBC_values, tracer_registry_type +use MOM_variables, only : thermo_var_ptrs +use tidal_bay_initialization, only : tidal_bay_set_OBC_data + +implicit none ; private + +#include + +public update_OBC_data + +integer :: id_clock_pass + +character(len=40) :: mod = "MOM_boundary_update" ! This module's name. +! This include declares and sets the variable "version". +#include "version_variable.h" + +contains + +!> Calls appropriate routine to update the open boundary conditions. +subroutine update_OBC_data(OBC, G, h, Time) + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< layer thickness + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure + type(time_type), intent(in) :: Time !< Model time + ! Local variables + logical :: read_OBC_eta = .false. + logical :: read_OBC_uv = .false. + logical :: read_OBC_TS = .false. + integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz + integer :: isd_off, jsd_off + integer :: IsdB, IedB, JsdB, JedB + character(len=40) :: mod = "set_Flather_Bdry_Conds" ! This subroutine's name. + character(len=200) :: filename, OBC_file, inputdir ! Strings for file/path + + real :: temp_u(G%domain%niglobal+1,G%domain%njglobal) + real :: temp_v(G%domain%niglobal,G%domain%njglobal+1) + + real, pointer, dimension(:,:,:) :: & + OBC_T_u => NULL(), & ! These arrays should be allocated and set to + OBC_T_v => NULL(), & ! specify the values of T and S that should come + OBC_S_u => NULL(), & + OBC_S_v => NULL() + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + + if (OBC%OBC_values_config == "tidal_bay") then + call tidal_bay_set_OBC_data(OBC, G, h, Time) + endif + +end subroutine update_OBC_data + +!> \namespace mom_boundary_update +!! This module updates the open boundary arrays when time-varying. +!! It caused a circular dependency with the tidal_bay setup when +!! MOM_open_boundary. +!! +!! A small fragment of the grid is shown below: +!! +!! j+1 x ^ x ^ x At x: q, CoriolisBu +!! j+1 > o > o > At ^: v, tauy +!! j x ^ x ^ x At >: u, taux +!! j > o > o > At o: h, bathyT, buoy, tr, T, S, Rml, ustar +!! j-1 x ^ x ^ x +!! i-1 i i+1 At x & ^: +!! i i+1 At > & o: +!! +!! The boundaries always run through q grid points (x). + +end module MOM_boundary_update diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 44dd7c9f23..91a544cd03 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -24,7 +24,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_SIMPLE, OBC_FLATHER +use MOM_open_boundary, only : ocean_OBC_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 @@ -156,7 +156,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, apply_OBC_u_flather_east = .false. ; apply_OBC_u_flather_west = .false. apply_OBC_v_flather_north = .false. ; apply_OBC_v_flather_south = .false. - if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_pe) then + if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then apply_OBC_u_flather_east = OBC%apply_OBC_u_flather_east apply_OBC_u_flather_west = OBC%apply_OBC_u_flather_west apply_OBC_v_flather_north = OBC%apply_OBC_v_flather_north @@ -193,20 +193,18 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, if (apply_OBC_u_flather_east .or. apply_OBC_u_flather_west) then do k=1,nz ; do j=LB%jsh,LB%jeh do I=LB%ish,LB%ieh+1 - if (OBC%OBC_kind_u(I-1,j) == OBC_FLATHER .and. (OBC%OBC_direction_u(I-1,j) == OBC_DIRECTION_E)) & - h(i,j,k) = h_input(i-1,j,k) + if (OBC%OBC_segment_u(I-1,j) /= OBC_NONE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I-1,j))%radiation & + .and. (OBC%OBC_direction_u(I-1,j) == OBC_DIRECTION_E)) & + h(i,j,k) = h_input(i-1,j,k) + endif enddo do i=LB%ish-1,LB%ieh - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W)) & - h(i,j,k) = h_input(i+1,j,k) - enddo - enddo - do J=LB%jsh-1,LB%jeh - do i=LB%ish-1,LB%ieh+1 - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_E)) & - v(i,J,k) = v(i-1,J,k) - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_W)) & - v(i,J,k) = v(i+1,J,k) + if (OBC%OBC_segment_u(I,j) /= OBC_NONE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation & + .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W)) & + h(i,j,k) = h_input(i+1,j,k) + endif enddo enddo ; enddo endif @@ -229,18 +227,18 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, if (apply_OBC_v_flather_north .or. apply_OBC_v_flather_south) then do k=1,nz do J=LB%jsh,LB%jeh+1 ; do i=LB%ish-1,LB%ieh+1 - if (OBC%OBC_kind_v(i,J-1) == OBC_FLATHER .and. (OBC%OBC_direction_v(i,J-1) == OBC_DIRECTION_N)) & - h(i,j,k) = h_input(i,j-1,k) + if (OBC%OBC_segment_v(i,J-1) /= OBC_NONE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J-1))%radiation & + .and. (OBC%OBC_direction_v(i,J-1) == OBC_DIRECTION_N)) & + h(i,j,k) = h_input(i,j-1,k) + endif enddo ; enddo do J=LB%jsh-1,LB%jeh ; do i=LB%ish-1,LB%ieh+1 - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S)) & - h(i,j,k) = h_input(i,j+1,k) - enddo ; enddo - do j=LB%jsh-1,LB%jeh+1 ; do I=LB%ish-1,LB%ieh - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_N)) & - u(I,j,k) = u(I,j-1,k) - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_S)) & - u(I,j,k) = u(I,j+1,k) + if (OBC%OBC_segment_v(i,J) /= OBC_NONE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation & + .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S)) & + h(i,j,k) = h_input(i,j+1,k) + endif enddo ; enddo enddo endif @@ -262,18 +260,18 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, if (apply_OBC_v_flather_north .or. apply_OBC_v_flather_south) then do k=1,nz do J=LB%jsh,LB%jeh+1 ; do i=LB%ish-1,LB%ieh+1 - if (OBC%OBC_kind_v(i,J-1) == OBC_FLATHER .and. (OBC%OBC_direction_v(i,J-1) == OBC_DIRECTION_N)) & - h(i,j,k) = h_input(i,j-1,k) + if (OBC%OBC_segment_v(i,J-1) /= OBC_NONE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J-1))%radiation & + .and. (OBC%OBC_direction_v(i,J-1) == OBC_DIRECTION_N)) & + h(i,j,k) = h_input(i,j-1,k) + endif enddo ; enddo do J=LB%jsh-1,LB%jeh ; do i=LB%ish-1,LB%ieh+1 - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S)) & - h(i,j,k) = h_input(i,j+1,k) - enddo ; enddo - do j=LB%jsh-1,LB%jeh+1 ; do I=LB%ish-1,LB%ieh - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_N)) & - u(I,j,k) = u(I,j-1,k) - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_S)) & - u(I,j,k) = u(I,j+1,k) + if (OBC%OBC_segment_v(i,J) /= OBC_NONE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation & + .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S)) & + h(i,j,k) = h_input(i,j+1,k) + endif enddo ; enddo enddo endif @@ -296,20 +294,18 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, if (apply_OBC_u_flather_east .or. apply_OBC_u_flather_west) then do k=1,nz ; do j=LB%jsh,LB%jeh do I=LB%ish,LB%ieh+1 - if (OBC%OBC_kind_u(I-1,j) == OBC_FLATHER .and. (OBC%OBC_direction_u(I-1,j) == OBC_DIRECTION_E)) & - h(i,j,k) = h_input(i-1,j,k) + if (OBC%OBC_segment_u(I-1,j) /= OBC_NONE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I-1,j))%radiation & + .and. (OBC%OBC_direction_u(I-1,j) == OBC_DIRECTION_E)) & + h(i,j,k) = h_input(i-1,j,k) + endif enddo do i=LB%ish-1,LB%ieh - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W)) & - h(i,j,k) = h_input(i+1,j,k) - enddo - enddo - do J=LB%jsh-1,LB%jeh - do i=LB%ish-1,LB%ieh+1 - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_E)) & - v(i,J,k) = v(i-1,J,k) - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER .and. (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_W)) & - v(i,J,k) = v(i+1,J,k) + if (OBC%OBC_segment_u(I,j) /= OBC_NONE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation & + .and. (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W)) & + h(i,j,k) = h_input(i+1,j,k) + endif enddo enddo ; enddo endif @@ -432,7 +428,8 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & uh(:,j,k), duhdu(:,k), visc_rem(:,k), & dt, G, j, ish, ieh, do_i, CS%vol_CFL) if (apply_OBC_u) then ; do I=ish-1,ieh - if (OBC%OBC_mask_u(I,j) .and. (OBC%OBC_kind_u(I,j) == OBC_SIMPLE)) & + if (OBC%OBC_mask_u(I,j) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified)) & uh(I,j,k) = OBC%uh(I,j,k) enddo ; endif enddo @@ -524,11 +521,11 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & if (present(uhbt) .or. do_aux .or. set_BT_cont) then if (apply_OBC_u) then ; do I=ish-1,ieh do_i(I) = .not.(OBC%OBC_mask_u(I,j) .and. & - (OBC%OBC_kind_u(I,j) == OBC_SIMPLE)) + (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified)) if (.not.do_i(I)) any_simple_OBC = .true. enddo ; else if (apply_OBC_flather) then ; do I=ish-1,ieh do_i(I) = .not.(OBC%OBC_mask_u(I,j) .and. & - (OBC%OBC_kind_u(I,j) == OBC_FLATHER) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) .and. & ((OBC%OBC_direction_u(I,j) == OBC_DIRECTION_N) .or. & (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_S))) enddo ; else ; do I=ish-1,ieh @@ -544,7 +541,8 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & if (present(u_cor)) then ; do k=1,nz do I=ish-1,ieh ; u_cor(I,j,k) = u(I,j,k) + du(I) * visc_rem(I,k) ; enddo if (apply_OBC_u) then ; do I=ish-1,ieh - if (OBC%OBC_mask_u(I,j) .and. (OBC%OBC_kind_u(I,j) == OBC_SIMPLE)) & + if (OBC%OBC_mask_u(I,j) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified)) & u_cor(I,j,k) = OBC%u(I,j,k) enddo ; endif enddo ; endif ! u-corrected @@ -559,7 +557,8 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & do k=1,nz do I=ish-1,ieh ; u_cor_aux(I,j,k) = u(I,j,k) + du(I) * visc_rem(I,k) ; enddo if (apply_OBC_u) then ; do I=ish-1,ieh - if (OBC%OBC_mask_u(I,j) .and. (OBC%OBC_kind_u(I,j) == OBC_SIMPLE)) & + if (OBC%OBC_mask_u(I,j) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified)) & u_cor_aux(I,j,k) = OBC%u(I,j,k) enddo ; endif enddo @@ -572,7 +571,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & if (any_simple_OBC) then do I=ish-1,ieh do_i(I) = (OBC%OBC_mask_u(I,j) .and. & - (OBC%OBC_kind_u(I,j) == OBC_SIMPLE)) + (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified)) if (do_i(I)) BT_cont%Fa_u_W0(I,j) = GV%H_subroundoff*G%dy_Cu(I,j) enddo do k=1,nz ; do I=ish-1,ieh ; if (do_i(I)) then @@ -1149,7 +1148,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & use_visc_rem = present(visc_rem_v) apply_OBC_v = .false. ; set_BT_cont = .false. ; apply_OBC_flather = .false. if (present(BT_cont)) set_BT_cont = (associated(BT_cont)) - if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_pe) then + if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then apply_OBC_v = OBC%apply_OBC_v apply_OBC_flather = OBC%apply_OBC_u_flather_east .or. & OBC%apply_OBC_u_flather_west .or. & @@ -1199,7 +1198,8 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & vh(:,J,k), dvhdv(:,k), visc_rem(:,k), & dt, G, J, ish, ieh, do_i, CS%vol_CFL) if (apply_OBC_v) then ; do i=ish,ieh - if (OBC%OBC_mask_v(i,J) .and. (OBC%OBC_kind_v(i,J) == OBC_SIMPLE)) & + if (OBC%OBC_mask_v(i,J) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%specified)) & vh(i,J,k) = OBC%vh(i,J,k) enddo ; endif enddo ! k-loop @@ -1287,11 +1287,11 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & if (present(vhbt) .or. do_aux .or. set_BT_cont) then if (apply_OBC_v) then ; do i=ish,ieh do_i(i) = .not.(OBC%OBC_mask_v(i,J) .and. & - (OBC%OBC_kind_v(i,J) == OBC_SIMPLE)) + (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%specified)) if (.not.do_i(i)) any_simple_OBC = .true. enddo ; else if (apply_OBC_flather) then ; do i=ish,ieh do_i(i) = .not.(OBC%OBC_mask_v(i,J) .and. & - (OBC%OBC_kind_v(i,J) == OBC_FLATHER) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) .and. & ((OBC%OBC_direction_v(i,J) == OBC_DIRECTION_E) .or. & (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_W))) enddo ; else ; do i=ish,ieh @@ -1307,7 +1307,8 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & if (present(v_cor)) then ; do k=1,nz do i=ish,ieh ; v_cor(i,J,k) = v(i,J,k) + dv(i) * visc_rem(i,k) ; enddo if (apply_OBC_v) then ; do i=ish,ieh - if (OBC%OBC_mask_v(i,J) .and. (OBC%OBC_kind_v(i,J) == OBC_SIMPLE)) & + if (OBC%OBC_mask_v(i,J) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%specified)) & v_cor(i,J,k) = OBC%v(i,J,k) enddo ; endif enddo ; endif ! v-corrected @@ -1321,7 +1322,8 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & do k=1,nz do i=ish,ieh ; v_cor_aux(i,J,k) = v(i,J,k) + dv(i) * visc_rem(i,k) ; enddo if (apply_OBC_v) then ; do i=ish,ieh - if (OBC%OBC_mask_v(i,J) .and. (OBC%OBC_kind_v(i,J) == OBC_SIMPLE)) & + if (OBC%OBC_mask_v(i,J) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%specified)) & v_cor_aux(i,J,k) = OBC%v(i,J,k) enddo ; endif enddo @@ -1334,7 +1336,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & if (any_simple_OBC) then do i=ish,ieh do_i(i) = (OBC%OBC_mask_v(i,J) .and. & - (OBC%OBC_kind_v(i,J) == OBC_SIMPLE)) + (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%specified)) if (do_i(i)) BT_cont%Fa_v_S0(i,J) = GV%H_subroundoff*G%dx_Cv(I,j) enddo do k=1,nz ; do i=ish,ieh ; if (do_i(i)) then diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index bce371cb19..4f4be21585 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -66,7 +66,6 @@ module MOM_dynamics_legacy_split !********+*********+*********+*********+*********+*********+*********+** -use MOM_open_boundary, only : ocean_OBC_type use MOM_variables, only : vertvisc_type, thermo_var_ptrs use MOM_variables, only : BT_cont_type, alloc_bt_cont_type, dealloc_bt_cont_type use MOM_variables, only : accel_diag_ptrs, ocean_internal_state, cont_diag_ptrs @@ -110,7 +109,9 @@ module MOM_dynamics_legacy_split use MOM_interface_heights, only : find_eta use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_MEKE_types, only : MEKE_type +use MOM_open_boundary, only : ocean_OBC_type use MOM_open_boundary, only : Radiation_Open_Bdry_Conds +use MOM_boundary_update, only : update_OBC_data use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_remnant @@ -495,6 +496,10 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & call cpu_clock_end(id_clock_pass) endif + if (associated(CS%OBC)) then; if (CS%OBC%update_OBC) then + call update_OBC_data(CS%OBC, G, h, Time_local) + endif; endif + ! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av call cpu_clock_begin(id_clock_Cor) call CorAdCalc(u_av, v_av, h_av, uh, vh, CS%CAu, CS%CAv, CS%ADp, G, GV, & @@ -784,6 +789,10 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & call cpu_clock_end(id_clock_pass) endif + if (associated(CS%OBC)) then; if (CS%OBC%update_OBC) then + call update_OBC_data(CS%OBC, G, h, Time_local) + endif; endif + if (BT_cont_BT_thick) then call cpu_clock_begin(id_clock_pass) call pass_vector(CS%BT_cont%h_u, CS%BT_cont%h_v, G%Domain, & diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 4bec4fec97..f4a895ded4 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -47,6 +47,7 @@ module MOM_dynamics_split_RK2 use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type use MOM_open_boundary, only : Radiation_Open_Bdry_Conds +use MOM_boundary_update, only : update_OBC_data use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS @@ -431,6 +432,10 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & call disable_averaging(CS%diag) if (showCallTree) call callTree_wayPoint("done with PressureForce (step_MOM_dyn_split_RK2)") + if (associated(CS%OBC)) then; if (CS%OBC%update_OBC) then + call update_OBC_data(CS%OBC, G, h, Time_local) + endif; endif + if (G%nonblocking_updates) then call cpu_clock_begin(id_clock_pass) call start_group_pass(CS%pass_eta_PF_eta, G%Domain) @@ -681,6 +686,10 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & if (showCallTree) call callTree_wayPoint("done with PressureForce[hp=(1-b).h+b.h] (step_MOM_dyn_split_RK2)") endif + if (associated(CS%OBC)) then; if (CS%OBC%update_OBC) then + call update_OBC_data(CS%OBC, G, h, Time_local) + endif; endif + if (G%nonblocking_updates) then call cpu_clock_begin(id_clock_pass) call complete_group_pass(CS%pass_av_uvh, G%Domain) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index fc392cbb88..b1587764e0 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -104,6 +104,7 @@ module MOM_dynamics_unsplit use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type use MOM_open_boundary, only : Radiation_Open_Bdry_Conds +use MOM_boundary_update, only : update_OBC_data use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS @@ -325,6 +326,10 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, & CS%PressureForce_CSp, CS%ALE_CSp, p_surf) call cpu_clock_end(id_clock_pres) + if (associated(CS%OBC)) then; if (CS%OBC%update_OBC) then + call update_OBC_data(CS%OBC, G, h, Time_local) + endif; endif + ! up = u + dt_pred * (PFu + CAu) call cpu_clock_begin(id_clock_mom_update) do k=1,nz ; do j=js,je ; do I=Isq,Ieq @@ -408,6 +413,10 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, & CS%PressureForce_CSp, CS%ALE_CSp, p_surf) call cpu_clock_end(id_clock_pres) + if (associated(CS%OBC)) then; if (CS%OBC%update_OBC) then + call update_OBC_data(CS%OBC, G, h, Time_local) + endif; endif + ! upp = u + dt/2 * ( PFu + CAu ) call cpu_clock_begin(id_clock_mom_update) do k=1,nz ; do j=js,je ; do I=Isq,Ieq @@ -482,6 +491,10 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, & CS%PressureForce_CSp, CS%ALE_CSp, p_surf) call cpu_clock_end(id_clock_pres) + if (associated(CS%OBC)) then; if (CS%OBC%update_OBC) then + call update_OBC_data(CS%OBC, G, h, Time_local) + endif; endif + ! u = u + dt * ( PFu + CAu ) do k=1,nz ; do j=js,je ; do I=Isq,Ieq u(i,j,k) = G%mask2dCu(i,j) * (u(i,j,k) + dt * & diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 8e228a9189..b87cdbe52b 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -102,6 +102,7 @@ module MOM_dynamics_unsplit_RK2 use MOM_MEKE_types, only : MEKE_type use MOM_open_boundary, only : ocean_OBC_type use MOM_open_boundary, only : Radiation_Open_Bdry_Conds +use MOM_boundary_update, only : update_OBC_data use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_CS use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_CS @@ -320,6 +321,10 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, call pass_vector(CS%CAu, CS%CAv, G%Domain) call cpu_clock_end(id_clock_pass) + if (associated(CS%OBC)) then; if (CS%OBC%update_OBC) then + call update_OBC_data(CS%OBC, G, h_in, Time_local) + endif; endif + ! up+[n-1/2] = u[n-1] + dt_pred * (PFu + CAu) call cpu_clock_begin(id_clock_mom_update) do k=1,nz ; do j=js,je ; do I=Isq,Ieq diff --git a/src/core/MOM_legacy_barotropic.F90 b/src/core/MOM_legacy_barotropic.F90 index cbf1e94a60..4cd5b21ccc 100644 --- a/src/core/MOM_legacy_barotropic.F90 +++ b/src/core/MOM_legacy_barotropic.F90 @@ -353,9 +353,7 @@ module MOM_legacy_barotropic OBC_mask_v => NULL() integer, dimension(:,:), pointer :: & OBC_direction_u => NULL(), & - OBC_direction_v => NULL(), & - OBC_kind_u => NULL(), & - OBC_kind_v => NULL() + OBC_direction_v => NULL() real, dimension(:,:), pointer :: & Cg_u => NULL(), & ! The external wave speed at u-points, in m s-1. Cg_v => NULL(), & ! The external wave speed at u-points, in m s-1. @@ -2234,11 +2232,11 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, if (associated(BT_OBC%OBC_mask_u)) then do j=js,je ; do I=is-1,ie ; if (BT_OBC%OBC_mask_u(I,j)) then - if (BT_OBC%OBC_kind_u(I,j) == OBC_SIMPLE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified) then uhbt(I,j) = BT_OBC%uhbt(I,j) ubt(I,j) = BT_OBC%ubt_outer(I,j) vel_trans = ubt(I,j) - elseif (BT_OBC%OBC_kind_u(I,j) == OBC_FLATHER) then + elseif (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then if (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) then cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL u_inlet = cfl*ubt_old(I-1,j) + (1.0-cfl)*ubt_old(I,j) ! Valid for cfl<1 @@ -2280,7 +2278,7 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, endif endif - if (BT_OBC%OBC_kind_u(I,j) /= OBC_SIMPLE) then + if (.not. OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified) then if (use_BT_cont) then uhbt(I,j) = find_uhbt(vel_trans,BTCL_u(I,j)) + uhbt0(I,j) else @@ -2294,11 +2292,11 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, if (associated(BT_OBC%OBC_mask_v)) then do J=js-1,je ; do i=is,ie ; if (BT_OBC%OBC_mask_v(i,J)) then - if (BT_OBC%OBC_kind_v(i,J) == OBC_SIMPLE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%specified) then vhbt(i,J) = BT_OBC%vhbt(i,J) vbt(i,J) = BT_OBC%vbt_outer(i,J) vel_trans = vbt(i,J) - elseif (BT_OBC%OBC_kind_v(i,J) == OBC_FLATHER) then + elseif (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then if (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) then cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(I,j) ! CFL v_inlet = cfl*vbt_old(i,J-1) + (1.0-cfl)*vbt_old(i,J) ! Valid for cfl<1 @@ -2348,7 +2346,7 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, endif endif - if (BT_OBC%OBC_kind_v(i,J) /= OBC_SIMPLE) then + if (OBC%OBC_segment_v(i,J) /= OBC_SIMPLE) then if (use_BT_cont) then vhbt(i,J) = find_vhbt(vel_trans,BTCL_v(i,J)) + vhbt0(i,J) else @@ -2397,7 +2395,7 @@ subroutine apply_eta_OBCs(OBC, eta, ubt, vbt, BT_OBC, G, MS, halo, dtbt) if ((OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west) .and. & associated(BT_OBC%OBC_mask_u)) then do j=js,je ; do I=is-1,ie ; if (BT_OBC%OBC_mask_u(I,j)) then - if (BT_OBC%OBC_kind_u(I,j) == OBC_FLATHER) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then if (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) then cfl = dtbt * BT_OBC%Cg_u(I,j) * G%IdxCu(I,j) ! CFL u_inlet = cfl*ubt(I-1,j) + (1.0-cfl)*ubt(I,j) ! Valid for cfl <1 @@ -2424,7 +2422,7 @@ subroutine apply_eta_OBCs(OBC, eta, ubt, vbt, BT_OBC, G, MS, halo, dtbt) if ((OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south) .and. & associated(BT_OBC%OBC_mask_v)) then do J=js-1,je ; do i=is,ie ; if (BT_OBC%OBC_mask_v(i,J)) then - if (BT_OBC%OBC_kind_v(i,J) == OBC_FLATHER) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then if (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) then cfl = dtbt*BT_OBC%Cg_v(i,J)*G%IdyCv(i,J) ! CFL v_inlet = cfl*vbt(i,J-1) + (1.0-cfl)*vbt(i,J) ! Valid for cfl <1 @@ -2502,7 +2500,6 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D allocate(BT_OBC%ubt_outer(isdw-1:iedw,jsdw:jedw)) ; BT_OBC%ubt_outer(:,:) = 0.0 allocate(BT_OBC%eta_outer_u(isdw-1:iedw,jsdw:jedw)) ; BT_OBC%eta_outer_u(:,:) = 0.0 allocate(BT_OBC%OBC_mask_u(isdw-1:iedw,jsdw:jedw)) ; BT_OBC%OBC_mask_u(:,:)=.false. - allocate(BT_OBC%OBC_kind_u(isdw-1:iedw,jsdw:jedw)) ; BT_OBC%OBC_kind_u(:,:)=OBC_NONE allocate(BT_OBC%OBC_direction_u(isdw-1:iedw,jsdw:jedw)); BT_OBC%OBC_direction_u(:,:)=OBC_NONE allocate(BT_OBC%Cg_v(isdw:iedw,jsdw-1:jedw)) ; BT_OBC%Cg_v(:,:) = 0.0 @@ -2511,13 +2508,11 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D allocate(BT_OBC%vbt_outer(isdw:iedw,jsdw-1:jedw)) ; BT_OBC%vbt_outer(:,:) = 0.0 allocate(BT_OBC%eta_outer_v(isdw:iedw,jsdw-1:jedw)) ; BT_OBC%eta_outer_v(:,:)=0.0 allocate(BT_OBC%OBC_mask_v(isdw:iedw,jsdw-1:jedw)) ; BT_OBC%OBC_mask_v(:,:)=.false. - allocate(BT_OBC%OBC_kind_v(isdw-1:iedw,jsdw:jedw)) ; BT_OBC%OBC_kind_v(:,:)=OBC_NONE allocate(BT_OBC%OBC_direction_v(isdw-1:iedw,jsdw:jedw)); BT_OBC%OBC_direction_v(:,:)=OBC_NONE if (associated(OBC%OBC_mask_u)) then do j=js-1,je+1 ; do I=is-1,ie BT_OBC%OBC_mask_u(I,j) = OBC%OBC_mask_u(I,j) - BT_OBC%OBC_kind_u(I,j) = OBC%OBC_kind_u(I,j) BT_OBC%OBC_direction_u(I,j) = OBC%OBC_direction_u(I,j) enddo ; enddo if (OBC%apply_OBC_u) then @@ -2526,7 +2521,7 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D enddo ; enddo ; enddo endif do j=js,je ; do I=is-1,ie ; if (OBC%OBC_mask_u(I,j)) then - if (OBC%OBC_kind_u(I,j) == OBC_SIMPLE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified) then if (use_BT_cont) then BT_OBC%ubt_outer(I,j) = uhbt_to_ubt(BT_OBC%uhbt(I,j),BTCL_u(I,j)) else @@ -2553,7 +2548,6 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D if (associated(OBC%OBC_mask_v)) then do J=js-1,je ; do i=is-1,ie+1 BT_OBC%OBC_mask_v(i,J) = OBC%OBC_mask_v(i,J) - BT_OBC%OBC_kind_v(i,J) = OBC%OBC_kind_v(i,J) BT_OBC%OBC_direction_v(i,J) = OBC%OBC_direction_v(i,J) enddo ; enddo if (OBC%apply_OBC_v) then @@ -2563,7 +2557,7 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D endif do J=js-1,je ; do i=is,ie ; if (OBC%OBC_mask_v(i,J)) then - if (OBC%OBC_kind_v(i,J) == OBC_SIMPLE) then + if (OBC%OBC_segment_v(i,J) == OBC_SIMPLE) then if (use_BT_cont) then BT_OBC%vbt_outer(i,J) = vhbt_to_vbt(BT_OBC%vhbt(i,J),BTCL_v(i,J)) else @@ -2603,7 +2597,6 @@ subroutine destroy_BT_OBC(BT_OBC) type(BT_OBC_type), intent(inout) :: BT_OBC if (associated(BT_OBC%OBC_mask_u)) deallocate(BT_OBC%OBC_mask_u) - if (associated(BT_OBC%OBC_kind_u)) deallocate(BT_OBC%OBC_kind_u) if (associated(BT_OBC%OBC_direction_u)) deallocate(BT_OBC%OBC_direction_u) deallocate(BT_OBC%Cg_u) deallocate(BT_OBC%H_u) @@ -2612,7 +2605,6 @@ subroutine destroy_BT_OBC(BT_OBC) deallocate(BT_OBC%eta_outer_u) if (associated(BT_OBC%OBC_mask_v)) deallocate(BT_OBC%OBC_mask_v) - if (associated(BT_OBC%OBC_kind_v)) deallocate(BT_OBC%OBC_kind_v) if (associated(BT_OBC%OBC_direction_v)) deallocate(BT_OBC%OBC_direction_v) deallocate(BT_OBC%Cg_v) deallocate(BT_OBC%H_v) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index d1bd34d6ad..31cddec1ad 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -14,6 +14,8 @@ module MOM_open_boundary use MOM_dyn_horgrid, only : dyn_horgrid_type use MOM_io, only : EAST_FACE, NORTH_FACE use MOM_io, only : slasher, read_data +use MOM_obsolete_params, only : obsolete_logical, obsolete_int, obsolete_real, obsolete_char +use MOM_string_functions, only : extract_word, remove_spaces use MOM_tracer_registry, only : add_tracer_OBC_values, tracer_registry_type use MOM_variables, only : thermo_var_ptrs @@ -28,18 +30,32 @@ module MOM_open_boundary public open_boundary_impose_normal_slope public open_boundary_impose_land_mask public Radiation_Open_Bdry_Conds -public set_Flather_positions public set_Flather_data integer, parameter, public :: OBC_NONE = 0, OBC_SIMPLE = 1, OBC_WALL = 2 integer, parameter, public :: OBC_FLATHER = 3 +integer, parameter, public :: OBC_RADIATION2D = 4 integer, parameter, public :: OBC_DIRECTION_N = 100 !< Indicates the boundary is an effective northern boundary integer, parameter, public :: OBC_DIRECTION_S = 200 !< Indicates the boundary is an effective southern boundary integer, parameter, public :: OBC_DIRECTION_E = 300 !< Indicates the boundary is an effective eastern boundary integer, parameter, public :: OBC_DIRECTION_W = 400 !< Indicates the boundary is an effective western boundary +!> Open boundary segment type - we'll have one for each open segment +!! to describe that segment. +type, public :: OBC_segment_type + logical :: radiation !< Radiation boundary. + logical :: radiation2D !< Oblique waves supported at radiation boundary. + logical :: nudged !< Optional supplement to radiation boundary. + logical :: specified !< Boundary fixed to external value. + logical :: gradient !< Zero gradient at boundary. + integer :: direction !< Boundary faces one of the four directions. + real :: Tnudge_in !< Nudging timescale on inflow. + real :: Tnudge_out !< Nudging timescale on outflow. +end type OBC_segment_type + !> Open-boundary data type, public :: ocean_OBC_type + integer :: number_of_segments = 0 !< The number of open-boundary segments. logical :: apply_OBC_u_flather_east = .false. !< True if any zonal velocity points in the !! local domain use east-facing Flather OBCs. logical :: apply_OBC_u_flather_west = .false. !< True if any zonal velocity points in the @@ -57,14 +73,18 @@ module MOM_open_boundary ! points, and can be OBC_NONE, OBC_SIMPLE, OBC_WALL, or OBC_FLATHER. Generally these ! should be consistent with OBC_mask_[uv], with OBC_mask_[uv] .false. for OBC_kind_[uv] = NONE ! and true for all other values. - integer, pointer, dimension(:,:) :: & - OBC_kind_u => NULL(), & !< Type of OBC at u-points. - OBC_kind_v => NULL() !< Type of OBC at v-points. ! These arrays indicate the outward-pointing orientation of the open boundary and will be set to ! one of OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_DIRECTION_E or OBC_DIRECTION_W. integer, pointer, dimension(:,:) :: & OBC_direction_u => NULL(), & !< Orientation of OBC at u-points. OBC_direction_v => NULL() !< Orientation of OBC at v-points. + ! Properties of the segments used. + type(OBC_segment_type), pointer, dimension(:) :: & + OBC_segment_list => NULL() !< List of segment objects. + ! Which segment object describes the current point. + integer, pointer, dimension(:,:) :: & + OBC_segment_u => NULL(), & !< Segment number of u-points. + OBC_segment_v => NULL() !< Segment number of v-points. ! The following apply at points with OBC_kind_[uv] = OBC_FLATHER. real, pointer, dimension(:,:,:) :: & rx_old_u => NULL(), & !< The rx_old_u value for radiation coeff for u-velocity in x-direction @@ -101,7 +121,9 @@ module MOM_open_boundary real :: rx_max !< The maximum magnitude of the baroclinic radiation !! velocity (or speed of characteristics), in m s-1. The !! default value is 10 m s-1. - logical :: this_pe !< Is there an open boundary on this tile? + logical :: OBC_pe !< Is there an open boundary on this tile? + logical :: update_OBC = .false. !< Is the open boundary info going to get updated? + character(len=200) :: OBC_values_config end type ocean_OBC_type integer :: id_clock_pass @@ -117,35 +139,92 @@ subroutine open_boundary_config(G, param_file, OBC) type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure type(param_file_type), intent(in) :: param_file !< Parameter file handle type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + ! Local variables + integer :: l ! For looping over segments + character(len=15) :: segment_param_str ! The run-time parameter name for each segment + character(len=100) :: segment_str ! The contents (rhs) for parameter "segment_param_str" allocate(OBC) call log_version(param_file, mod, version, "Controls where open boundaries are located, what "//& "kind of boundary condition to impose, and what data to apply, if any.") - call get_param(param_file, mod, "APPLY_OBC_U", OBC%apply_OBC_u, & - "If true, open boundary conditions may be set at some \n"//& - "u-points, with the configuration controlled by OBC_CONFIG", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_V", OBC%apply_OBC_v, & - "If true, open boundary conditions may be set at some \n"//& - "v-points, with the configuration controlled by OBC_CONFIG", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_EAST", OBC%apply_OBC_u_flather_east, & - "Apply a Flather open boundary condition on the eastern\n"//& - "side of the global domain", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_WEST", OBC%apply_OBC_u_flather_west, & - "Apply a Flather open boundary condition on the western\n"//& - "side of the global domain", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_NORTH", OBC%apply_OBC_v_flather_north, & - "Apply a Flather open boundary condition on the northern\n"//& - "side of the global domain", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_SOUTH", OBC%apply_OBC_v_flather_south, & - "Apply a Flather open boundary condition on the southern\n"//& - "side of the global domain", & - default=.false.) + call get_param(param_file, mod, "OBC_NUMBER_OF_SEGMENTS", OBC%number_of_segments, & + "The number of open boundary segments.", & + default=0) + if (OBC%number_of_segments == 0) then + ! For the interim, if segments are not set we'll interpret the old-style run-time OBC parameters as before + call get_param(param_file, mod, "APPLY_OBC_U", OBC%apply_OBC_u, & + "If true, open boundary conditions may be set at some \n"//& + "u-points, with the configuration controlled by OBC_CONFIG", & + default=.false.) + call get_param(param_file, mod, "APPLY_OBC_V", OBC%apply_OBC_v, & + "If true, open boundary conditions may be set at some \n"//& + "v-points, with the configuration controlled by OBC_CONFIG", & + default=.false.) + call get_param(param_file, mod, "OBC_CONFIG", OBC%OBC_values_config, & + "If set, open boundary configuration string", & + default="file") + call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_EAST", OBC%apply_OBC_u_flather_east, & + "Apply a Flather open boundary condition on the eastern\n"//& + "side of the global domain", & + default=.false.) + call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_WEST", OBC%apply_OBC_u_flather_west, & + "Apply a Flather open boundary condition on the western\n"//& + "side of the global domain", & + default=.false.) + call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_NORTH", OBC%apply_OBC_v_flather_north, & + "Apply a Flather open boundary condition on the northern\n"//& + "side of the global domain", & + default=.false.) + call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_SOUTH", OBC%apply_OBC_v_flather_south, & + "Apply a Flather open boundary condition on the southern\n"//& + "side of the global domain", & + default=.false.) + else + ! When specifying segments, we will consider the old-style run-time OBC parameters obsolete + ! Note: once we have finished transitioning to segments we can permanently obsolete these parameters + ! rather than conditionally. + call obsolete_logical(param_file, "APPLY_OBC_U", hint="APPLY_OBC_U cannot be used when using OBC_SEGMENTS") + call obsolete_logical(param_file, "APPLY_OBC_V", hint="APPLY_OBC_V cannot be used when using OBC_SEGMENTS") + call obsolete_logical(param_file, "APPLY_OBC_U_FLATHER_EAST", hint="APPLY_OBC_U_FLATHER_EAST cannot be used when using OBC_SEGMENTS") + call obsolete_logical(param_file, "APPLY_OBC_U_FLATHER_WEST", hint="APPLY_OBC_U_FLATHER_WEST cannot be used when using OBC_SEGMENTS") + call obsolete_logical(param_file, "APPLY_OBC_V_FLATHER_NORTH", hint="APPLY_OBC_V_FLATHER_NORTH cannot be used when using OBC_SEGMENTS") + call obsolete_logical(param_file, "APPLY_OBC_V_FLATHER_SOUTH", hint="APPLY_OBC_V_FLATHER_SOUTH cannot be used when using OBC_SEGMENTS") + ! Allocate everything + allocate(OBC%OBC_segment_list(0:OBC%number_of_segments)) + do l=0,OBC%number_of_segments + OBC%OBC_segment_list(l)%radiation = .false. + OBC%OBC_segment_list(l)%radiation2D = .false. + OBC%OBC_segment_list(l)%nudged = .false. + OBC%OBC_segment_list(l)%specified = .false. + OBC%OBC_segment_list(l)%gradient = .false. + OBC%OBC_segment_list(l)%direction = OBC_NONE + OBC%OBC_segment_list(l)%Tnudge_in = 0.0 + OBC%OBC_segment_list(l)%Tnudge_out = 0.0 + enddo + allocate(OBC%OBC_mask_u(G%IsdB:G%IedB,G%jsd:G%jed)) ; OBC%OBC_mask_u(:,:) = .false. + allocate(OBC%OBC_direction_u(G%IsdB:G%IedB,G%jsd:G%jed)) ; OBC%OBC_direction_u(:,:) = OBC_NONE + allocate(OBC%OBC_segment_u(G%IsdB:G%IedB,G%jsd:G%jed)) ; OBC%OBC_segment_u(:,:) = OBC_NONE + allocate(OBC%OBC_mask_v(G%isd:G%ied,G%JsdB:G%JedB)) ; OBC%OBC_mask_v(:,:) = .false. + allocate(OBC%OBC_direction_v(G%isd:G%ied,G%JsdB:G%JedB)) ; OBC%OBC_direction_v(:,:) = OBC_NONE + allocate(OBC%OBC_segment_v(G%isd:G%ied,G%JsdB:G%JedB)) ; OBC%OBC_segment_v(:,:) = OBC_NONE + + do l = 1, OBC%number_of_segments + write(segment_param_str(1:15),"('OBC_SEGMENT_',i3.3)") l + call get_param(param_file, mod, segment_param_str, segment_str, & + "Documentation needs to be dynamic?????", & + fail_if_missing=.true.) + segment_str = remove_spaces(segment_str) + if (segment_str(1:2) == 'I=') then + call setup_u_point_obc(OBC, G, segment_str, l) + elseif (segment_str(1:2) == 'J=') then + call setup_v_point_obc(OBC, G, segment_str, l) + else + call MOM_error(FATAL, "MOM_open_boundary.F90, open_boundary_config: "//& + "Unable to interpret "//segment_param_str//" = "//trim(segment_str)) + endif + enddo + endif ! Safety check if ((OBC%apply_OBC_u_flather_west .or. OBC%apply_OBC_v_flather_south) .and. & @@ -163,6 +242,294 @@ subroutine open_boundary_config(G, param_file, OBC) end subroutine open_boundary_config +!> Parse an OBC_SEGMENT_%%% string starting with "I=" and configure placement and type of OBC accordingly +subroutine setup_u_point_obc(OBC, G, segment_str, l_seg) + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure + character(len=*), intent(in) :: segment_str !< A string in form of "I=%,J=%:%,string" + integer, intent(in) :: l_seg !< which segment is this? + ! Local variables + integer :: I_obc, Js_obc, Je_obc ! Position of segment in global index space + integer :: j, this_kind + character(len=32) :: action_str + + ! This returns the global indices for the segment + call parse_segment_str(G%ieg, G%jeg, segment_str, I_obc, Js_obc, Je_obc, action_str ) + I_obc = I_obc - G%idg_offset ! Convert to local tile indices on this tile + Js_obc = Js_obc - G%jdg_offset ! Convert to local tile indices on this tile + Je_obc = Je_obc - G%jdg_offset ! Convert to local tile indices on this tile + +! if (Je_obc>Js_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_E +! if (Je_obcJs_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_E + if (Je_obcJs_obc) OBC%apply_OBC_u_flather_east = .true. ! This line will not be needed soon - AJA + if (Je_obcJs_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_E + if (Je_obcJs_obc) OBC%apply_OBC_u_flather_east = .true. ! This line will not be needed soon - AJA + if (Je_obcG%HI%IedB) return ! Boundary is not on tile + if (max(Js_obc,Je_obc)G%HI%JedB) return ! Segment is not on tile + + ! These four lines extend the open boundary into the halo region of tiles on the edge of the physical + ! domain. They are used to reproduce the checksums of the circle_obcs test case and will be removed + ! in the fullness of time. -AJA +! These were causing grief in the supercritical problem. - KSH +! if (Js_obc == G%HI%JscB) Js_obc = G%HI%jsd-1 +! if (Js_obc == G%HI%JecB) Js_obc = G%HI%jed +! if (Je_obc == G%HI%JscB) Je_obc = G%HI%jsd-1 +! if (Je_obc == G%HI%JecB) Je_obc = G%HI%jed + + do j=G%HI%jsd, G%HI%jed + if (j>min(Js_obc,Je_obc) .and. j<=max(Js_obc,Je_obc)) then + OBC%OBC_mask_u(I_obc,j) = .true. + OBC%OBC_segment_u(I_obc,j) = l_seg + if (Je_obc>Js_obc) then ! East is outward + if (this_kind == OBC_FLATHER) then + OBC%OBC_direction_u(I_obc,j) = OBC_DIRECTION_E ! We only use direction for Flather (maybe) + ! Set v points outside segment + OBC%OBC_mask_v(i_obc+1,J) = .true. + if (OBC%OBC_segment_v(i_obc+1,J) == OBC_NONE) then + OBC%OBC_direction_v(i_obc+1,J) = OBC_DIRECTION_E + OBC%OBC_segment_v(i_obc+1,J) = l_seg + endif + OBC%OBC_mask_v(i_obc+1,J-1) = .true. + if (OBC%OBC_segment_v(i_obc+1,J-1) == OBC_NONE) then + OBC%OBC_direction_v(i_obc+1,J-1) = OBC_DIRECTION_E + OBC%OBC_segment_v(i_obc+1,J-1) = l_seg + endif + endif + else ! West is outward + if (this_kind == OBC_FLATHER) then + OBC%OBC_direction_u(I_obc,j) = OBC_DIRECTION_W ! We only use direction for Flather + ! Set v points outside segment + OBC%OBC_mask_v(i_obc,J) = .true. + if (OBC%OBC_segment_v(i_obc,J) == OBC_NONE) then + OBC%OBC_direction_v(i_obc,J) = OBC_DIRECTION_W + OBC%OBC_segment_v(i_obc,J) = l_seg + endif + OBC%OBC_mask_v(i_obc,J-1) = .true. + if (OBC%OBC_segment_v(i_obc,J-1) == OBC_NONE) then + OBC%OBC_direction_v(i_obc,J-1) = OBC_DIRECTION_W + OBC%OBC_segment_v(i_obc,J-1) = l_seg + endif + endif + endif + endif + enddo + +end subroutine setup_u_point_obc + +!> Parse an OBC_SEGMENT_%%% string starting with "J=" and configure placement and type of OBC accordingly +subroutine setup_v_point_obc(OBC, G, segment_str, l_seg) + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure + character(len=*), intent(in) :: segment_str !< A string in form of "J=%,I=%:%,string" + integer, intent(in) :: l_seg !< which segment is this? + ! Local variables + integer :: J_obc, Is_obc, Ie_obc ! Position of segment in global index space + integer :: i, this_kind + character(len=32) :: action_str + + ! This returns the global indices for the segment + call parse_segment_str(G%ieg, G%jeg, segment_str, J_obc, Is_obc, Ie_obc, action_str ) + J_obc = J_obc - G%jdg_offset ! Convert to local tile indices on this tile + Is_obc = Is_obc - G%idg_offset ! Convert to local tile indices on this tile + Ie_obc = Ie_obc - G%idg_offset ! Convert to local tile indices on this tile + + if (trim(action_str) == 'FLATHER') then + this_kind = OBC_FLATHER + if (Ie_obc>Is_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_S + if (Ie_obcIs_obc) OBC%apply_OBC_v_flather_south = .true. ! This line will not be needed soon - AJA + if (Ie_obcIs_obc) OBC%OBC_segment_list(l_seg)%direction = OBC_DIRECTION_S + if (Ie_obcIs_obc) OBC%apply_OBC_v_flather_south = .true. ! This line will not be needed soon - AJA + if (Ie_obcG%HI%JedB) return ! Boundary is not on tile + if (max(Is_obc,Ie_obc)G%HI%IedB) return ! Segment is not on tile + + ! These four lines extend the open boundary into the halo region of tiles on the edge of the physical + ! domain. They are used to reproduce the checksums of the circle_obcs test case and will be removed + ! in the fullness of time. -AJA +! These cause trouble with DOME +! if (Is_obc == G%HI%IscB) Is_obc = G%HI%isd-1 +! if (Is_obc == G%HI%IecB) Is_obc = G%HI%ied +! if (Ie_obc == G%HI%IscB) Ie_obc = G%HI%isd-1 +! if (Ie_obc == G%HI%IecB) Ie_obc = G%HI%ied + + do i=G%HI%isd, G%HI%ied + if (i>min(Is_obc,Ie_obc) .and. i<=max(Is_obc,Ie_obc)) then + OBC%OBC_mask_v(i,J_obc) = .true. + OBC%OBC_segment_v(i,J_obc) = l_seg + if (Is_obc>Ie_obc) then ! North is outward + if (this_kind == OBC_FLATHER) then + OBC%OBC_direction_v(i,J_obc) = OBC_DIRECTION_N ! We only use direction for Flather + ! Set u points outside segment + OBC%OBC_mask_u(I,j_obc+1) = .true. + if (OBC%OBC_segment_u(I,j_obc+1) == OBC_NONE) then + OBC%OBC_direction_u(I,j_obc+1) = OBC_DIRECTION_N + OBC%OBC_segment_u(I,j_obc+1) = l_seg + endif + OBC%OBC_mask_u(I-1,j_obc+1) = .true. + if (OBC%OBC_segment_u(I-1,j_obc+1) == OBC_NONE) then + OBC%OBC_direction_u(I-1,j_obc+1) = OBC_DIRECTION_N + OBC%OBC_segment_u(I-1,j_obc+1) = l_seg + endif + endif + else ! South is outward + if (this_kind == OBC_FLATHER) then + OBC%OBC_direction_v(i,J_obc) = OBC_DIRECTION_S ! We only use direction for Flather + ! Set u points outside segment + OBC%OBC_mask_u(I,j_obc) = .true. + if (OBC%OBC_segment_u(I,j_obc) == OBC_NONE) then + OBC%OBC_direction_u(I,j_obc) = OBC_DIRECTION_S + OBC%OBC_segment_u(I,j_obc) = l_seg + endif + OBC%OBC_mask_u(I-1,j_obc) = .true. + if (OBC%OBC_segment_u(I-1,j_obc) == OBC_NONE) then + OBC%OBC_direction_u(I-1,j_obc) = OBC_DIRECTION_S + OBC%OBC_segment_u(I-1,j_obc) = l_seg + endif + endif + endif + endif + enddo + +end subroutine setup_v_point_obc + +!> Parse an OBC_SEGMENT_%%% string +subroutine parse_segment_str(ni_global, nj_global, segment_str, l, m, n, action_str ) + integer, intent(in) :: ni_global !< Number of h-points in zonal direction + integer, intent(in) :: nj_global !< Number of h-points in meridional direction + character(len=*), intent(in) :: segment_str !< A string in form of "I=l,J=m:n,string" or "J=l,I=m,n,string" + integer, intent(out) :: l !< The value of I=l, if segment_str begins with I=l, or the value of J=l + integer, intent(out) :: m !< The value of J=m, if segment_str begins with I=, or the value of I=m + integer, intent(out) :: n !< The value of J=n, if segment_str begins with I=, or the value of I=n + character(len=*), intent(out) :: action_str !< The "string" part of segment_str + ! Local variables + character(len=24) :: word1, word2, m_word, n_word !< Words delineated by commas in a string in form of "I=%,J=%:%,string" + integer :: l_max !< Either ni_global or nj_global, depending on whether segment_str begins with "I=" or "J=" + integer :: mn_max !< Either nj_global or ni_global, depending on whether segment_str begins with "I=" or "J=" + + ! Process first word which will started with either 'I=' or 'J=' + word1 = extract_word(segment_str,',',1) + word2 = extract_word(segment_str,',',2) + if (word1(1:2)=='I=') then + l_max = ni_global + mn_max = nj_global + if (.not. (word2(1:2)=='J=')) call MOM_error(FATAL, "MOM_open_boundary.F90, parse_segment_str: "//& + "Second word of string '"//trim(segment_str)//"' must start with 'J='.") + elseif (word1(1:2)=='J=') then ! Note that the file_parser uniformaly expands "=" to " = " + l_max = nj_global + mn_max = ni_global + if (.not. (word2(1:2)=='I=')) call MOM_error(FATAL, "MOM_open_boundary.F90, parse_segment_str: "//& + "Second word of string '"//trim(segment_str)//"' must start with 'I='.") + else + call MOM_error(FATAL, "MOM_open_boundary.F90, parse_segment_str"//& + "String '"//segment_str//"' must start with 'I=' or 'J='.") + endif + + ! Read l + l = interpret_int_expr( word1(3:24), l_max ) + if (l<0 .or. l>l_max) then + call MOM_error(FATAL, "MOM_open_boundary.F90, parse_segment_str: "//& + "First value from string '"//trim(segment_str)//"' is outside of the physical domain.") + endif + + ! Read m + m_word = extract_word(word2(3:24),':',1) + m = interpret_int_expr( m_word, mn_max ) + if (m<-1 .or. m>mn_max+1) then + call MOM_error(FATAL, "MOM_open_boundary.F90, parse_segment_str: "//& + "Beginning of range in string '"//trim(segment_str)//"' is outside of the physical domain.") + endif + + ! Read m + n_word = extract_word(word2(3:24),':',2) + n = interpret_int_expr( n_word, mn_max ) + if (n<-1 .or. n>mn_max+1) then + call MOM_error(FATAL, "MOM_open_boundary.F90, parse_segment_str: "//& + "End of range in string '"//trim(segment_str)//"' is outside of the physical domain.") + endif + + if (abs(n-m)==0) then + call MOM_error(FATAL, "MOM_open_boundary.F90, parse_segment_str: "//& + "Range in string '"//trim(segment_str)//"' must span one cell.") + endif + + ! Type of open boundary condition + action_str = extract_word(segment_str,',',3) + + contains + + ! Returns integer value interpreted from string in form of %I, N or N-%I + integer function interpret_int_expr(string, imax) + character(len=*), intent(in) :: string !< Integer in form or %I, N or N-%I + integer, intent(in) :: imax !< Value to replace 'N' with + ! Local variables + integer slen + + slen = len_trim(string) + if (slen==0) call MOM_error(FATAL, "MOM_open_boundary.F90, parse_segment_str"//& + "Parsed string was empty!") + if (len_trim(string)==1 .and. string(1:1)=='N') then + interpret_int_expr = imax + elseif (string(1:1)=='N') then + read(string(3:slen),*,err=911) interpret_int_expr + if (string(2:2)=='-') then + interpret_int_expr = imax - interpret_int_expr + elseif (string(2:2)=='+') then + interpret_int_expr = imax + interpret_int_expr + else + goto 911 + endif + else + read(string(1:slen),*,err=911) interpret_int_expr + endif + return + 911 call MOM_error(FATAL, "MOM_open_boundary.F90, parse_segment_str"//& + "Problem reading value from string '"//trim(string)//"'.") + end function interpret_int_expr +end subroutine parse_segment_str + !> Initialize open boundary control structure subroutine open_boundary_init(G, param_file, OBC) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure @@ -219,8 +586,9 @@ subroutine open_boundary_dealloc(OBC) if (.not. associated(OBC)) return if (associated(OBC%OBC_mask_u)) deallocate(OBC%OBC_mask_u) if (associated(OBC%OBC_mask_v)) deallocate(OBC%OBC_mask_v) - if (associated(OBC%OBC_kind_u)) deallocate(OBC%OBC_kind_u) - if (associated(OBC%OBC_kind_v)) deallocate(OBC%OBC_kind_v) + if (associated(OBC%OBC_segment_list)) deallocate(OBC%OBC_segment_list) + if (associated(OBC%OBC_segment_u)) deallocate(OBC%OBC_segment_u) + if (associated(OBC%OBC_segment_v)) deallocate(OBC%OBC_segment_v) if (associated(OBC%rx_old_u)) deallocate(OBC%rx_old_u) if (associated(OBC%ry_old_v)) deallocate(OBC%ry_old_v) if (associated(OBC%rx_old_h)) deallocate(OBC%rx_old_h) @@ -252,17 +620,29 @@ subroutine open_boundary_impose_normal_slope(OBC, G, depth) if (.not.associated(OBC)) return - if (associated(OBC%OBC_direction_u)) then + if (associated(OBC%OBC_segment_u)) then do j=G%jsd,G%jed ; do I=G%isd,G%ied-1 - if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) depth(i+1,j) = depth(i,j) - if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W) depth(i,j) = depth(i+1,j) + if (OBC%OBC_segment_u(I,j) /= OBC_NONE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%direction == & + OBC_DIRECTION_E) depth(i+1,j) = depth(i,j) + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%direction == & + OBC_DIRECTION_W) depth(i,j) = depth(i+1,j) + endif +! if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) depth(i+1,j) = depth(i,j) +! if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W) depth(i,j) = depth(i+1,j) enddo ; enddo endif - if (associated(OBC%OBC_kind_v)) then + if (associated(OBC%OBC_segment_v)) then do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied - if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) depth(i,j+1) = depth(i,j) - if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) depth(i,j) = depth(i,j+1) + if (OBC%OBC_segment_v(i,J) /= OBC_NONE) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%direction == & + OBC_DIRECTION_N) depth(i,j+1) = depth(i,j) + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%direction == & + OBC_DIRECTION_S) depth(i,j) = depth(i,j+1) + endif +! if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) depth(i,j+1) = depth(i,j) +! if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) depth(i,j) = depth(i,j+1) enddo ; enddo endif @@ -278,22 +658,24 @@ subroutine open_boundary_impose_land_mask(OBC, G) if (.not.associated(OBC)) return - if (associated(OBC%OBC_kind_u)) then - do j=G%jsd,G%jed ; do I=G%isd,G%ied-1 - if (G%mask2dCu(I,j) == 0) then - OBC%OBC_kind_u(I,j) = OBC_NONE - OBC%OBC_direction_u(I,j) = OBC_NONE - OBC%OBC_mask_u(I,j) = .false. + if (associated(OBC%OBC_segment_u)) then + do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB + if (G%mask2dCu(I,j) == 0 .and. (OBC%OBC_segment_u(I,j) /= OBC_NONE)) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then + OBC%OBC_direction_u(I,j) = OBC_NONE + OBC%OBC_mask_u(I,j) = .false. + endif endif enddo ; enddo endif - if (associated(OBC%OBC_kind_v)) then - do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied - if (G%mask2dCv(i,J) == 0) then - OBC%OBC_kind_v(i,J) = OBC_NONE - OBC%OBC_direction_v(i,J) = OBC_NONE - OBC%OBC_mask_v(i,J) = .false. + if (associated(OBC%OBC_segment_v)) then + do J=G%JsdB,G%JedB ; do i=G%isd,G%ied + if (G%mask2dCv(i,J) == 0 .and. (OBC%OBC_segment_v(i,J) /= OBC_NONE)) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then + OBC%OBC_direction_v(i,J) = OBC_NONE + OBC%OBC_mask_v(i,J) = .false. + endif endif enddo ; enddo endif @@ -307,9 +689,6 @@ subroutine open_boundary_impose_land_mask(OBC, G) ! bathymetry inside the boundary was do shallow and flagged as land. if (OBC%OBC_mask_u(I,j)) any_U = .true. enddo ; enddo -! if (.not. any_U) then -! deallocate(OBC%OBC_mask_u) -! endif endif any_V = .false. @@ -317,14 +696,10 @@ subroutine open_boundary_impose_land_mask(OBC, G) do J=G%JsdB,G%JedB ; do i=G%isd,G%ied if (OBC%OBC_mask_v(i,J)) any_V = .true. enddo ; enddo -! if (.not. any_V) then -! deallocate(OBC%OBC_mask_v) -! endif endif -! if (.not.(any_U .or. any_V)) call open_boundary_dealloc(OBC) - OBC%this_pe = .true. - if (.not.(any_U .or. any_V)) OBC%this_pe = .false. + OBC%OBC_pe = .true. + if (.not.(any_U .or. any_V)) OBC%OBC_pe = .false. end subroutine open_boundary_impose_land_mask @@ -340,9 +715,12 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_new real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old ! Local variables - real :: dhdt, dhdx, gamma_u, gamma_h, gamma_v + real, dimension(SZI_(G),SZJ_(G)) :: grad + real :: dhdt, dhdx, dhdy, gamma_u, gamma_h, gamma_v + real :: cff, Cx, Cy real :: rx_max, ry_max ! coefficients for radiation real :: rx_new, rx_avg ! coefficients for radiation + real, parameter :: eps = 1.0e-20 integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -355,9 +733,29 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & gamma_u = OBC%gamma_uv ; gamma_v = OBC%gamma_uv ; gamma_h = OBC%gamma_h rx_max = OBC%rx_max ; ry_max = OBC%rx_max - if (OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west) then - do k=1,nz ; do j=js,je ; do I=is-1,ie ; if (OBC%OBC_mask_u(I,j)) then - if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) then + do k=1,nz ; do j=js,je ; do I=is-1,ie ; if (OBC%OBC_mask_u(I,j)) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%direction == OBC_DIRECTION_E) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation2D) then + grad(I,J) = u_old(I,j+1,k) - u_old(I,j,k) + grad(I,J-1) = u_old(I,j,k) - u_old(I,j-1,k) + grad(I-1,J) = u_old(I-1,j+1,k) - u_old(I-1,j,k) + grad(I-1,J-1) = u_old(I-1,j,k) - u_old(I-1,j-1,k) + dhdt = u_old(I-1,j,k)-u_new(I-1,j,k) !old-new + dhdx = u_new(I-1,j,k)-u_new(I-2,j,k) !in new time backward sasha for I-1 + if (dhdt*dhdx < 0.0) dhdt = 0.0 + if (dhdt*(grad(I-1,J) + grad(I-1,J-1)) > 0.0) then + dhdy = grad(I-1,J-1) + else + dhdy = grad(I-1,J) + endif + cff = max(dhdx*dhdx + dhdy*dhdy, eps) + Cx = dhdt*dhdx + Cy = min(cff,max(dhdt*dhdy,-cff)) + ! Turning off radiation2D part + Cy = 0 + u_new(I,j,k) = ((cff*u_old(I,j,k) + Cx*u_new(I-1,j,k)) - & + (max(Cy,0.0)*grad(I,J-1) - min(Cy,0.0)*grad(I,J))) / (cff + Cx) + elseif (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then dhdt = u_old(I-1,j,k)-u_new(I-1,j,k) !old-new dhdx = u_new(I-1,j,k)-u_new(I-2,j,k) !in new time backward sasha for I-1 rx_new = 0.0 @@ -374,7 +772,30 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & ! OBC%rx_old_h(I,j,k) = rx_avg ! h_new(I+1,j,k) = (h_old(I+1,j,k) + rx_avg*h_new(I,j,k)) / (1.0+rx_avg) !original endif - if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W) then + endif + + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%direction == OBC_DIRECTION_W) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation2D) then + grad(I,J) = u_old(I,j+1,k) - u_old(I,j,k) + grad(I,J-1) = u_old(I,j,k) - u_old(I,j-1,k) + grad(I+1,J) = u_old(I+1,j+1,k) - u_old(I+1,j,k) + grad(I+1,J-1) = u_old(I+1,j,k) - u_old(I+1,j-1,k) + dhdt = u_old(I+1,j,k)-u_new(I+1,j,k) !old-new + dhdx = u_new(I+1,j,k)-u_new(I+2,j,k) !in new time backward sasha for I+1 + if (dhdt*dhdx < 0.0) dhdt = 0.0 + if (dhdt*(grad(I+1,J) + grad(I+1,J-1)) > 0.0) then + dhdy = grad(I+1,J-1) + else + dhdy = grad(I+1,J) + endif + cff = max(dhdx*dhdx + dhdy*dhdy, eps) + Cx = dhdt*dhdx + Cy = min(cff,max(dhdt*dhdy,-cff)) + ! Turning off radiation2D part + Cy = 0 + u_new(I,j,k) = ((cff*u_old(I,j,k) + Cx*u_new(I+1,j,k)) - & + (max(Cy,0.0)*grad(I,J-1) - min(Cy,0.0)*grad(I,J))) / (cff + Cx) + elseif (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then dhdt = u_old(I+1,j,k)-u_new(I+1,j,k) !old-new dhdx = u_new(I+1,j,k)-u_new(I+2,j,k) !in new time backward sasha for I+1 rx_new = 0.0 @@ -391,16 +812,88 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & ! OBC%rx_old_h(I,j,k) = rx_avg ! h_new(I,j,k) = (h_old(I,j,k) + rx_avg*h_new(I+1,j,k)) / (1.0+rx_avg) !original endif - endif ; enddo ; enddo ; enddo - endif + endif - if (OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south) then - do k=1,nz ; do J=js-1,je ; do i=is,ie ; if (OBC%OBC_mask_v(i,J)) then - if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%direction == OBC_DIRECTION_N) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%gradient) then + u_new(I,j,k) = u_new(I,j-1,k) + elseif (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then + grad(i,j) = u_old(I,j,k) - u_old(I-1,j,k) + grad(i,j-1) = u_old(I,j-1,k) - u_old(I-1,j-1,k) + grad(i+1,j) = u_old(I+1,j,k) - u_old(I,j,k) + grad(i+1,j-1) = u_old(I+1,j-1,k) - u_old(I,j-1,k) + dhdt = u_old(I,j-1,k)-u_new(I,j-1,k) !old-new + dhdy = u_new(I,j-1,k)-u_new(I,j-2,k) !in new time backward sasha for I+1 + if (dhdt*dhdy < 0.0) dhdt = 0.0 + if (dhdt*(grad(i,j-1) + grad(i+1,j-1)) > 0.0) then + dhdx = grad(i,j-1) + else + dhdx = grad(i+1,j-1) + endif + cff = max(dhdx*dhdx + dhdy*dhdy, eps) + Cx = 0.0 + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation2D) & + Cx = min(cff, max(dhdt*dhdx, -cff)) + Cy = dhdt*dhdy + u_new(I,j,k) = ((cff*u_old(I,j,k) + Cy*u_new(I,j-1,k)) - & + (max(Cx, 0.0)*grad(i,j) - min(Cx, 0.0)*grad(i+1,j)))/(cff + Cy) + endif + endif + + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%direction == OBC_DIRECTION_S) then + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%gradient) then + u_new(I,j,k) = u_new(I,j+1,k) + elseif (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) then + grad(i,j) = u_old(I,j,k) - u_old(I-1,j,k) + grad(i,j+1) = u_old(I,j+1,k) - u_old(I-1,j+1,k) + grad(i+1,j) = u_old(I+1,j,k) - u_old(I,j,k) + grad(i+1,j+1) = u_old(I+1,j+1,k) - u_old(I,j+1,k) + dhdt = u_old(I,j+1,k)-u_new(I,j+1,k) !old-new + dhdy = u_new(I,j+1,k)-u_new(I,j+2,k) !in new time backward sasha for I+1 + if (dhdt*dhdy < 0.0) dhdt = 0.0 + if (dhdt*(grad(i,j+1) + grad(i+1,j+1)) > 0.0) then + dhdx = grad(i,j+1) + else + dhdx = grad(i+1,j+1) + endif + cff = max(dhdx*dhdx + dhdy*dhdy, eps) + Cx = 0.0 + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation2D) & + Cx = min(cff, max(dhdt*dhdx, -cff)) + Cy = dhdt*dhdy + u_new(I,j,k) = ((cff*u_old(I,j,k) + Cy*u_new(I,j+1,k)) - & + (max(Cx, 0.0)*grad(i,j) - min(Cx, 0.0)*grad(i+1,j)))/(cff + Cy) + endif + endif + endif ; enddo ; enddo ; enddo + + do k=1,nz ; do J=js-1,je ; do i=is,ie ; if (OBC%OBC_mask_v(i,J)) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%direction == OBC_DIRECTION_N) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation2D) then + grad(I,J) = v_old(i+1,J,k) - v_old(i,J,k) + grad(I-1,J) = v_old(i,J,k) - v_old(i-1,J,k) + grad(I,J-1) = v_old(i+1,J-1,k) - v_old(i,J-1,k) + grad(I-1,J-1) = v_old(i,J-1,k) - v_old(i-1,J-1,k) + dhdt = v_old(i,J-1,k)-v_new(i,J-1,k) !old-new + dhdy = v_new(i,J-1,k)-v_new(i,J-2,k) !in new time backward sasha for J-1 + if (dhdt*dhdy < 0.0) dhdt = 0.0 + if (dhdt*(grad(I,J-1) + grad(I-1,J-1)) > 0.0) then + dhdx = grad(I-1,J-1) + else + dhdx = grad(I,J-1) + endif + cff = max(dhdx*dhdx + dhdy*dhdy, eps) + Cy = dhdt*dhdy + Cx = min(cff,max(dhdt*dhdx,-cff)) + ! Turning off radiation2D part + Cx = 0 + v_new(i,J,k) = ((cff*v_old(i,J,k) + Cy*v_new(i,J-1,k)) - & + (max(Cx,0.0)*grad(I-1,J) - min(Cx,0.0)*grad(I,J))) / (cff + Cy) + elseif (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then dhdt = v_old(i,J-1,k)-v_new(i,J-1,k) !old-new - dhdx = v_new(i,J-1,k)-v_new(i,J-2,k) !in new time backward sasha for J-1 + dhdy = v_new(i,J-1,k)-v_new(i,J-2,k) !in new time backward sasha for J-1 rx_new = 0.0 - if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max) + if (dhdt*dhdy > 0.0) rx_new = min( (dhdt/dhdy), rx_max) rx_avg = (1.0-gamma_v)*OBC%ry_old_v(i,J,k) + gamma_v*rx_new OBC%ry_old_v(i,J,k) = rx_avg v_new(i,J,k) = (v_old(i,J,k) + rx_avg*v_new(i,J-1,k)) / (1.0+rx_avg) @@ -413,12 +906,34 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & ! OBC%ry_old_h(i,J,k) = rx_avg ! h_new(i,J+1,k) = (h_old(i,J+1,k) + rx_avg*h_new(i,J,k)) / (1.0+rx_avg) !original endif + endif - if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%direction == OBC_DIRECTION_S) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation2D) then + grad(I,J) = v_old(i+1,J,k) - v_old(i,J,k) + grad(I-1,J) = v_old(i,J,k) - v_old(i-1,J,k) + grad(I,J+1) = v_old(i+1,J+1,k) - v_old(i,J+1,k) + grad(I-1,J+1) = v_old(i,J+1,k) - v_old(i-1,J+1,k) + dhdt = v_old(i,J+1,k)-v_new(i,J+1,k) !old-new + dhdy = v_new(i,J+1,k)-v_new(i,J+2,k) !in new time backward sasha for J+1 + if (dhdt*dhdy < 0.0) dhdt = 0.0 + if (dhdt*(grad(I,J+1) + grad(I-1,J+1)) > 0.0) then + dhdx = grad(I-1,J+1) + else + dhdx = grad(I,J+1) + endif + cff = max(dhdx*dhdx + dhdy*dhdy, eps) + Cy = dhdt*dhdy + Cx = min(cff,max(dhdt*dhdx,-cff)) + ! Turning off radiation2D part + Cx = 0 + v_new(i,J,k) = ((cff*v_old(i,J,k) + Cy*v_new(i,J+1,k)) - & + (max(Cx,0.0)*grad(I-1,J) - min(Cx,0.0)*grad(I,J))) / (cff + Cy) + elseif (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then dhdt = v_old(i,J+1,k)-v_new(i,J+1,k) !old-new - dhdx = v_new(i,J+1,k)-v_new(i,J+2,k) !in new time backward sasha for J+1 + dhdy = v_new(i,J+1,k)-v_new(i,J+2,k) !in new time backward sasha for J+1 rx_new = 0.0 - if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max) + if (dhdt*dhdy > 0.0) rx_new = min( (dhdt/dhdy), rx_max) rx_avg = (1.0-gamma_v)*OBC%ry_old_v(i,J,k) + gamma_v*rx_new OBC%ry_old_v(i,J,k) = rx_avg v_new(i,J,k) = (v_old(i,J,k) + rx_avg*v_new(i,J+1,k)) / (1.0+rx_avg) @@ -431,9 +946,59 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & ! OBC%ry_old_h(i,J,k) = rx_avg ! h_new(i,J,k) = (h_old(i,J,k) + rx_avg*h_new(i,J+1,k)) / (1.0+rx_avg) !original endif + endif - endif ; enddo ; enddo ; enddo - endif + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%direction == OBC_DIRECTION_E) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%gradient) then + v_new(i,J,k) = v_new(i-1,J,k) + elseif (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then + grad(i,j) = v_old(i,J,k) - v_old(i,J-1,k) + grad(i,j+1) = v_old(i,J+1,k) - v_old(i,J,k) + grad(i-1,j) = v_old(i-1,J,k) - v_old(i-1,J-1,k) + grad(i-1,j+1) = v_old(i-1,J+1,k) - v_old(i-1,J,k) + dhdt = v_old(i-1,J,k)-v_new(i-1,J,k) !old-new + dhdx = v_new(i-1,J,k)-v_new(i-2,J,k) !in new time backward sasha for I+1 + if (dhdt*dhdx < 0.0) dhdt = 0.0 + if (dhdt*(grad(i-1,j) + grad(i-1,j+1)) > 0.0) then + dhdy = grad(i-1,j) + else + dhdy = grad(i-1,j+1) + endif + cff = max(dhdx*dhdx + dhdy*dhdy, eps) + Cx = dhdt*dhdx + Cy = 0.0 + if (OBC%OBC_segment_list(OBC%OBC_segment_v(I,j))%radiation2D) & + Cy = min(cff, max(dhdt*dhdy, -cff)) + v_new(i,J,k) = ((cff*v_old(i,J,k) + Cx*v_new(i-1,J,k)) - & + (max(Cy, 0.0)*grad(i,j) - min(Cy, 0.0)*grad(i,j+1)))/(cff + Cx) + endif + endif + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%direction == OBC_DIRECTION_W) then + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%gradient) then + v_new(i,J,k) = v_new(i+1,J,k) + elseif (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) then + grad(i,j) = v_old(i,J,k) - v_old(i,J-1,k) + grad(i+1,j) = v_old(i+1,J,k) - v_old(i+1,J-1,k) + grad(i,j+1) = v_old(i,J+1,k) - v_old(i,J,k) + grad(i+1,j+1) = v_old(i+1,J+1,k) - v_old(i+1,J,k) + dhdt = v_old(i+1,J,k)-v_new(i+1,J,k) !old-new + dhdx = v_new(i+1,J,k)-v_new(i+2,J,k) !in new time backward sasha for I+1 + if (dhdt*dhdx < 0.0) dhdt = 0.0 + if (dhdt*(grad(i+1,j) + grad(i+1,j+1)) > 0.0) then + dhdy = grad(i+1,j) + else + dhdy = grad(i+1,j+1) + endif + cff = max(dhdx*dhdx + dhdy*dhdy, eps) + Cx = dhdt*dhdx + Cy = 0.0 + if (OBC%OBC_segment_list(OBC%OBC_segment_v(I,j))%radiation2D) & + Cy = min(cff, max(dhdt*dhdy, -cff)) + v_new(i,J,k) = ((cff*v_old(i,J,k) + Cx*v_new(i+1,J,k)) - & + (max(Cy, 0.0)*grad(i,j) - min(Cy, 0.0)*grad(i+1,j)))/(cff + Cx) + endif + endif + endif ; enddo ; enddo ; enddo call cpu_clock_begin(id_clock_pass) call pass_vector(u_new, v_new, G%Domain) @@ -442,135 +1007,6 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & end subroutine Radiation_Open_Bdry_Conds -!> Sets the domain boundaries as Flather open boundaries using the original -!! Flather run-time logicals -subroutine set_Flather_positions(G, OBC) - type(dyn_horgrid_type), intent(inout) :: G - type(ocean_OBC_type), pointer :: OBC - ! Local variables - integer :: east_boundary, west_boundary, north_boundary, south_boundary - integer :: i, j - - if (.not.associated(OBC%OBC_mask_u)) then - allocate(OBC%OBC_mask_u(G%IsdB:G%IedB,G%jsd:G%jed)) ; OBC%OBC_mask_u(:,:) = .false. - endif - if (.not.associated(OBC%OBC_direction_u)) then - allocate(OBC%OBC_direction_u(G%IsdB:G%IedB,G%jsd:G%jed)) ; OBC%OBC_direction_u(:,:) = OBC_NONE - endif - if (.not.associated(OBC%OBC_kind_u)) then - allocate(OBC%OBC_kind_u(G%IsdB:G%IedB,G%jsd:G%jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE - endif - if (.not.associated(OBC%OBC_mask_v)) then - allocate(OBC%OBC_mask_v(G%isd:G%ied,G%JsdB:G%JedB)) ; OBC%OBC_mask_v(:,:) = .false. - endif - if (.not.associated(OBC%OBC_direction_v)) then - allocate(OBC%OBC_direction_v(G%isd:G%ied,G%JsdB:G%JedB)) ; OBC%OBC_direction_v(:,:) = OBC_NONE - endif - if (.not.associated(OBC%OBC_kind_v)) then - allocate(OBC%OBC_kind_v(G%isd:G%ied,G%JsdB:G%JedB)) ; OBC%OBC_kind_v(:,:) = OBC_NONE - endif - - ! This code should be modified to allow OBCs to be applied anywhere. - - if (G%symmetric) then - east_boundary = G%ieg - west_boundary = G%isg-1 - north_boundary = G%jeg - south_boundary = G%jsg-1 - else - ! I am not entirely sure that this works properly. -RWH - east_boundary = G%ieg-1 - west_boundary = G%isg - north_boundary = G%jeg-1 - south_boundary = G%jsg - endif - - if (OBC%apply_OBC_u_flather_east) then - ! Determine where u points are applied at east side - do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB - if ((I+G%idg_offset) == east_boundary) then !eastern side - OBC%OBC_mask_u(I,j) = .true. - OBC%OBC_direction_u(I,j) = OBC_DIRECTION_E - OBC%OBC_kind_u(I,j) = OBC_FLATHER - OBC%OBC_mask_v(i+1,J) = .true. - if (OBC%OBC_direction_v(i+1,J) == OBC_NONE) then - OBC%OBC_direction_v(i+1,J) = OBC_DIRECTION_E - OBC%OBC_kind_v(i+1,J) = OBC_FLATHER - endif - OBC%OBC_mask_v(i+1,J-1) = .true. - if (OBC%OBC_direction_v(i+1,J-1) == OBC_NONE) then - OBC%OBC_direction_v(i+1,J-1) = OBC_DIRECTION_E - OBC%OBC_kind_v(i+1,J-1) = OBC_FLATHER - endif - endif - enddo ; enddo - endif - - if (OBC%apply_OBC_u_flather_west) then - ! Determine where u points are applied at west side - do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB - if ((I+G%idg_offset) == west_boundary) then !western side - OBC%OBC_mask_u(I,j) = .true. - OBC%OBC_direction_u(I,j) = OBC_DIRECTION_W - OBC%OBC_kind_u(I,j) = OBC_FLATHER - OBC%OBC_mask_v(i,J) = .true. - if (OBC%OBC_direction_v(i,J) == OBC_NONE) then - OBC%OBC_direction_v(i,J) = OBC_DIRECTION_W - OBC%OBC_kind_v(i,J) = OBC_FLATHER - endif - OBC%OBC_mask_v(i,J-1) = .true. - if (OBC%OBC_direction_v(i,J-1) == OBC_NONE) then - OBC%OBC_direction_v(i,J-1) = OBC_DIRECTION_W - OBC%OBC_kind_v(i,J-1) = OBC_FLATHER - endif - endif - enddo ; enddo - endif - - if (OBC%apply_OBC_v_flather_north) then - ! Determine where v points are applied at north side - do J=G%JsdB,G%JedB ; do i=G%isd,G%ied - if ((J+G%jdg_offset) == north_boundary) then !northern side - OBC%OBC_mask_v(i,J) = .true. - OBC%OBC_direction_v(i,J) = OBC_DIRECTION_N - OBC%OBC_kind_v(i,J) = OBC_FLATHER - OBC%OBC_mask_u(I,j+1) = .true. - if (OBC%OBC_direction_u(I,j+1) == OBC_NONE) then - OBC%OBC_direction_u(I,j+1) = OBC_DIRECTION_N - OBC%OBC_kind_u(I,j+1) = OBC_FLATHER - endif - OBC%OBC_mask_u(I-1,j+1) = .true. - if (OBC%OBC_direction_u(I-1,j+1) == OBC_NONE) then - OBC%OBC_direction_u(I-1,j+1) = OBC_DIRECTION_N - OBC%OBC_kind_u(I-1,j+1) = OBC_FLATHER - endif - endif - enddo ; enddo - endif - - if (OBC%apply_OBC_v_flather_south) then - ! Determine where v points are applied at south side - do J=G%JsdB,G%JedB ; do i=G%isd,G%ied - if ((J+G%jdg_offset) == south_boundary) then !southern side - OBC%OBC_mask_v(i,J) = .true. - OBC%OBC_direction_v(i,J) = OBC_DIRECTION_S - OBC%OBC_kind_v(i,J) = OBC_FLATHER - OBC%OBC_mask_u(I,j) = .true. - if (OBC%OBC_direction_u(I,j) == OBC_NONE) then - OBC%OBC_direction_u(I,j) = OBC_DIRECTION_S - OBC%OBC_kind_u(I,j) = OBC_FLATHER - endif - OBC%OBC_mask_u(I-1,j) = .true. - if (OBC%OBC_direction_u(I-1,j) == OBC_NONE) then - OBC%OBC_direction_u(I-1,j) = OBC_DIRECTION_S - OBC%OBC_kind_u(I-1,j) = OBC_FLATHER - endif - endif - enddo ; enddo - endif - -end subroutine set_Flather_positions - !> Sets the initial definitions of the characteristic open boundary conditions. !! \author Mehmet Ilicak subroutine set_Flather_data(OBC, tv, h, G, PF, tracer_Reg) diff --git a/src/framework/MOM_string_functions.F90 b/src/framework/MOM_string_functions.F90 index 7e02bb4788..e5e5df681b 100644 --- a/src/framework/MOM_string_functions.F90 +++ b/src/framework/MOM_string_functions.F90 @@ -36,6 +36,8 @@ module MOM_string_functions public left_real, left_reals public string_functions_unit_tests public extractWord +public extract_word +public remove_spaces public slasher contains @@ -206,44 +208,77 @@ function isFormattedFloatEqualTo(str, val) 987 return end function isFormattedFloatEqualTo -function extractWord(string,n) -! Returns string corresponding to the nth word in the argument -! or "" if the string is not long enough. Both spaces and commas -! are interpretted as separators. - character(len=*), intent(in) :: string - integer, intent(in) :: n - character(len=120) :: extractWord +!> Returns the string corresponding to the nth word in the argument +!! or "" if the string is not long enough. Both spaces and commas +!! are interpreted as separators. +character(len=120) function extractWord(string, n) + character(len=*), intent(in) :: string + integer, intent(in) :: n + + extractWord = extract_word(string, ' ,', n) + +end function extractWord + +!> Returns the string corresponding to the nth word in the argument +!! or "" if the string is not long enough. Words are delineated +!! by the mandatory separators argument. +character(len=120) function extract_word(string, separators, n) + character(len=*), intent(in) :: string !< String to scan + character(len=*), intent(in) :: separators !< Characters to use for delineation + integer, intent(in) :: n !< Number of word to extract ! Local variables integer :: ns, i, b, e, nw logical :: lastCharIsSeperator - extractWord = '' + extract_word = '' lastCharIsSeperator = .true. ns = len_trim(string) i = 0; b=0; e=0; nw=0; do while (i Returns string with all spaces removed. +character(len=120) function remove_spaces(string) + character(len=*), intent(in) :: string !< String to scan + ! Local variables + integer :: ns, i, o + logical :: lastCharIsSeperator + lastCharIsSeperator = .true. + ns = len_trim(string) + i = 0; o = 0 + do while (i OBC_mask_u - allocate(OBC%u(IsdB:IedB,jsd:jed,nz)) ; OBC%u(:,:,:) = 0.0 - allocate(OBC%uh(IsdB:IedB,jsd:jed,nz)) ; OBC%uh(:,:,:) = 0.0 - allocate(OBC%OBC_kind_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE - do j=jsd,jed ; do I=IsdB,IedB - if (OBC%OBC_mask_u(I,j)) OBC%OBC_kind_u(I,j) = OBC_SIMPLE - enddo ; enddo - endif - if (apply_OBC_v) then - OBC%apply_OBC_v = .true. - OBC%OBC_mask_v => OBC_mask_v - allocate(OBC%v(isd:ied,JsdB:JedB,nz)) ; OBC%v(:,:,:) = 0.0 - allocate(OBC%vh(isd:ied,JsdB:JedB,nz)) ; OBC%vh(:,:,:) = 0.0 - allocate(OBC%OBC_kind_v(isd:ied,JsdB:JedB)) ; OBC%OBC_kind_v(:,:) = OBC_NONE - do J=JsdB,JedB ; do i=isd,ied - if (OBC%OBC_mask_v(i,J)) OBC%OBC_kind_v(i,J) = OBC_SIMPLE - enddo ; enddo - endif +! Shouldn't be needed now, right??? -ksh +! if (apply_OBC_u) then +! OBC%apply_OBC_u = .true. +! OBC%OBC_mask_u => OBC_mask_u +! allocate(OBC%u(IsdB:IedB,jsd:jed,nz)) ; OBC%u(:,:,:) = 0.0 +! allocate(OBC%uh(IsdB:IedB,jsd:jed,nz)) ; OBC%uh(:,:,:) = 0.0 +! allocate(OBC%OBC_kind_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE +! do j=jsd,jed ; do I=IsdB,IedB +! if (OBC%OBC_mask_u(I,j)) OBC%OBC_kind_u(I,j) = OBC_SIMPLE +! enddo ; enddo +! endif +! if (apply_OBC_v) then +! OBC%apply_OBC_v = .true. +! OBC%OBC_mask_v => OBC_mask_v +! allocate(OBC%v(isd:ied,JsdB:JedB,nz)) ; OBC%v(:,:,:) = 0.0 +! allocate(OBC%vh(isd:ied,JsdB:JedB,nz)) ; OBC%vh(:,:,:) = 0.0 +! allocate(OBC%OBC_kind_v(isd:ied,JsdB:JedB)) ; OBC%OBC_kind_v(:,:) = OBC_NONE +! do J=JsdB,JedB ; do i=isd,ied +! if (OBC%OBC_mask_v(i,J)) OBC%OBC_kind_v(i,J) = OBC_SIMPLE +! enddo ; enddo +! endif if (apply_OBC_v) then do k=1,nz ; do J=Jsd,Jed ; do i=isd,ied diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 78939b1f75..2442fbf0bc 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -298,7 +298,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, h_neglect = GV%H_subroundoff h_neglect3 = h_neglect**3 - if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_pe) then + if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then apply_OBC = OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west .or. & OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south endif ; endif ; endif @@ -586,7 +586,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, G%IareaCu(I,j)) / (0.5*(h(i+1,j,k) + h(i,j,k)) + h_neglect) if (apply_OBC) then ; if (OBC%OBC_mask_u(I,j)) then - if ((OBC%OBC_kind_u(I,j) == OBC_FLATHER)) diffu(I,j,k) = 0.0 + if (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%radiation) diffu(I,j,k) = 0.0 endif ; endif enddo ; enddo @@ -598,7 +598,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, CS%DX2h(i,j+1)*str_xx(i,j+1))) * & G%IareaCv(i,J)) / (0.5*(h(i,j+1,k) + h(i,j,k)) + h_neglect) if (apply_OBC) then ; if (OBC%OBC_mask_v(i,J)) then - if ((OBC%OBC_kind_v(i,J) == OBC_FLATHER)) diffv(I,j,k) = 0.0 + if (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%radiation) diffv(I,j,k) = 0.0 endif ; endif enddo ; enddo diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 613672d46a..8e0299feee 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -388,16 +388,18 @@ subroutine vertvisc(u, v, h, fluxes, visc, dt, OBC, ADp, CDp, G, GV, CS, & call vertvisc_limit_vel(u, v, h, ADp, CDp, fluxes, visc, dt, G, GV, CS) ! Here the velocities associated with open boundary conditions are applied. - if (associated(OBC)) then ; if (OBC%this_pe) then + if (associated(OBC)) then ; if (OBC%OBC_pe) then if (OBC%apply_OBC_u) then do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq - if (OBC%OBC_mask_u(I,j) .and. (OBC%OBC_kind_u(I,j) == OBC_SIMPLE)) & + if (OBC%OBC_mask_u(I,j) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified)) & u(I,j,k) = OBC%u(I,j,k) enddo ; enddo ; enddo endif if (OBC%apply_OBC_v) then do k=1,nz ; do J=Jsq,Jeq ; do i=G%isc,G%iec - if (OBC%OBC_mask_v(i,J) .and. (OBC%OBC_kind_v(i,J) == OBC_SIMPLE)) & + if (OBC%OBC_mask_v(i,J) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%specified)) & v(i,J,k) = OBC%v(i,J,k) enddo ; enddo ; enddo endif diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 9686dc7eb2..5c9d81f902 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -470,7 +470,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & enddo ; enddo endif ! usePPM - if (associated(OBC)) then ; if (OBC%this_pe) then ; if (OBC%apply_OBC_u) then + if (associated(OBC)) then ; if (OBC%OBC_pe) then ; if (OBC%apply_OBC_u) then do_any_i = .false. do I=is-1,ie do_i(I) = .false. @@ -730,7 +730,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & enddo ; enddo endif ! usePPM - if (associated(OBC)) then ; if (OBC%this_pe) then ; if (OBC%apply_OBC_v) then + if (associated(OBC)) then ; if (OBC%OBC_pe) then ; if (OBC%apply_OBC_v) then do_any_i = .false. do i=is,ie do_i(i) = .false. diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index df686d546e..f5c33cedc1 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -38,7 +38,6 @@ module DOME_initialization public DOME_initialize_topography public DOME_initialize_thickness public DOME_initialize_sponges -public DOME_set_OBC_positions public DOME_set_OBC_data contains @@ -233,36 +232,6 @@ subroutine DOME_initialize_sponges(G, GV, tv, PF, CSp) end subroutine DOME_initialize_sponges -!> Set the positions of the open boundary needed for the DOME experiment. -subroutine DOME_set_OBC_positions(G, param_file, OBC) - type(dyn_horgrid_type), intent(in) :: G !< Grid structure. - type(param_file_type), intent(in) :: param_file !< Parameter file handle. - type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure. - ! Local variables - character(len=40) :: mod = "DOME_set_OBC_positions" ! This subroutine's name. - integer :: i, j - - if (.not.associated(OBC)) call MOM_error(FATAL, & - "DOME_initialization, DOME_set_OBC_positions: OBC type was not allocated!") - - if (OBC%apply_OBC_u) then - ! Set where u points are determined by OBCs. - !allocate(OBC_mask_u(IsdB:IedB,jsd:jed)) ; OBC_mask_u(:,:) = .false. - call MOM_error(FATAL,"DOME_initialization, DOME_set_OBC_positions: "//& - "APPLY_OBC_U=True is not coded for the DOME experiment") - endif - if (OBC%apply_OBC_v) then - ! Set where v points are determined by OBCs. - allocate(OBC%OBC_mask_v(G%isd:G%ied,G%JsdB:G%JedB)) ; OBC%OBC_mask_v(:,:) = .false. - do J=G%JsdB,G%JedB ; do i=G%isd,G%ied - if ((G%geoLonCv(i,J) > 1000.0) .and. (G%geoLonCv(i,J) < 1100.0) .and. & - (abs(G%geoLatCv(i,J) - G%gridLatB(G%JegB)) < 0.1)) then - OBC%OBC_mask_v(i,J) = .true. - endif - enddo ; enddo - endif -end subroutine DOME_set_OBC_positions - !> This subroutine sets the properties of flow at open boundary conditions. !! This particular example is for the DOME inflow describe in Legg et al. 2006. subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg) @@ -319,20 +288,10 @@ subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg) if (OBC%apply_OBC_u) then allocate(OBC%u(IsdB:IedB,jsd:jed,nz)) ; OBC%u(:,:,:) = 0.0 allocate(OBC%uh(IsdB:IedB,jsd:jed,nz)) ; OBC%uh(:,:,:) = 0.0 - allocate(OBC%OBC_kind_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE - allocate(OBC%OBC_direction_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_direction_u(:,:) = OBC_NONE - do j=jsd,jed ; do I=IsdB,IedB - if (OBC%OBC_mask_u(I,j)) OBC%OBC_kind_u(I,j) = OBC_SIMPLE - enddo ; enddo endif if (OBC%apply_OBC_v) then allocate(OBC%v(isd:ied,JsdB:JedB,nz)) ; OBC%v(:,:,:) = 0.0 allocate(OBC%vh(isd:ied,JsdB:JedB,nz)) ; OBC%vh(:,:,:) = 0.0 - allocate(OBC%OBC_kind_v(isd:ied,JsdB:JedB)) ; OBC%OBC_kind_v(:,:) = OBC_NONE - allocate(OBC%OBC_direction_v(isd:ied,JsdB:JedB)) ; OBC%OBC_direction_v(:,:) = OBC_NONE - do J=JsdB,JedB ; do i=isd,ied - if (OBC%OBC_mask_v(i,J)) OBC%OBC_kind_v(i,J) = OBC_SIMPLE - enddo ; enddo endif if (OBC%apply_OBC_v) then diff --git a/src/user/soliton_initialization.F90 b/src/user/soliton_initialization.F90 new file mode 100644 index 0000000000..deca8cdc47 --- /dev/null +++ b/src/user/soliton_initialization.F90 @@ -0,0 +1,116 @@ +!> Initial conditions for the Equatorial Rossby soliton test (Boyd). +module soliton_initialization + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, is_root_pe +use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_get_input, only : directories +use MOM_grid, only : ocean_grid_type +use MOM_io, only : close_file, fieldtype, file_exists +use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE +use MOM_io, only : write_field, slasher, vardesc +use MOM_variables, only : thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type +use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type +use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE +use regrid_consts, only : REGRIDDING_LAYER, REGRIDDING_ZSTAR +use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SIGMA + +implicit none ; private + +#include + +! Private (module-wise) parameters +character(len=40) :: mod = "soliton_initialization" !< This module's name. + +public soliton_initialize_thickness +public soliton_initialize_velocity + +contains + +!> Initialization of thicknesses in Equatorial Rossby soliton test +subroutine soliton_initialize_thickness(h, G) + type(ocean_grid_type), intent(in) :: G !< Grid structure + real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: h !< Thickness + + integer :: i, j, k, is, ie, js, je, nz + real :: x, y, x0, y0 + real :: val1, val2, val3, val4 + character(len=40) :: verticalCoordinate + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + + call MOM_mesg("soliton_initialization.F90, soliton_initialize_thickness: setting thickness") + + x0 = 2.0*G%len_lon/3.0 + y0 = 0.0 + val1 = 0.395 + val2 = 0.771*(val1*val1) + + do j = G%jsc,G%jec ; do i = G%isc,G%iec + do k = 1, nz + x = G%geoLonT(i,j)-x0 + y = G%geoLatT(i,j)-y0 + val3 = exp(-val1*x) + val4 = val2*((2.0*val3/(1.0+(val3*val3)))**2) + h(i,j,k) = 0.25*val4*(6.0*y*y+3.0)* & + exp(-0.5*y*y) + enddo + end do ; end do + +end subroutine soliton_initialize_thickness + + +!> Initialization of u and v in the equatorial Rossby soliton test +subroutine soliton_initialize_velocity(u, v, h, G) + type(ocean_grid_type), intent(in) :: G !< Grid structure + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: u !< i-component of velocity [m/s] + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: v !< j-component of velocity [m/s] + real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: h !< Thickness [H] + + real :: x, y, x0, y0 + real :: val1, val2, val3, val4 + integer :: i, j, k, is, ie, js, je, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + + x0 = 2.0*G%len_lon/3.0 + y0 = 0.0 + val1 = 0.395 + val2 = 0.771*(val1*val1) + + v(:,:,:) = 0.0 + u(:,:,:) = 0.0 + + do j = G%jsc,G%jec ; do I = G%isc-1,G%iec+1 + do k = 1, nz + x = 0.5*(G%geoLonT(i+1,j)+G%geoLonT(i,j))-x0 + y = 0.5*(G%geoLatT(i+1,j)+G%geoLatT(i,j))-y0 + val3 = exp(-val1*x) + val4 = val2*((2.0*val3/(1.0+(val3*val3)))**2) + u(I,j,k) = 0.25*val4*(6.0*y*y-9.0)* & + exp(-0.5*y*y) + enddo + enddo ; enddo + do j = G%jsc-1,G%jec+1 ; do I = G%isc,G%iec + do k = 1, nz + x = 0.5*(G%geoLonT(i,j+1)+G%geoLonT(i,j))-x0 + y = 0.5*(G%geoLatT(i,j+1)+G%geoLatT(i,j))-y0 + val3 = exp(-val1*x) + val4 = val2*((2.0*val3/(1.0+(val3*val3)))**2) + v(i,J,k) = 2.0*val4*y*(-2.0*val1*tanh(val1*x))* & + exp(-0.5*y*y) + enddo + enddo ; enddo + +end subroutine soliton_initialize_velocity + + +!> \namespace soliton_initialization +!! +!! \section section_soliton Description of the equatorial Rossby soliton initial +!! conditions +!! + +end module soliton_initialization diff --git a/src/user/supercritical_initialization.F90 b/src/user/supercritical_initialization.F90 new file mode 100644 index 0000000000..abb39e3ad4 --- /dev/null +++ b/src/user/supercritical_initialization.F90 @@ -0,0 +1,152 @@ +module supercritical_initialization +!*********************************************************************** +!* GNU General Public License * +!* This file is a part of MOM. * +!* * +!* MOM is free software; you can redistribute it and/or modify it and * +!* are expected to follow the terms of the GNU General Public License * +!* as published by the Free Software Foundation; either version 2 of * +!* the License, or (at your option) any later version. * +!* * +!* MOM is distributed in the hope that it will be useful, but WITHOUT * +!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * +!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * +!* License for more details. * +!* * +!* For the full text of the GNU General Public License, * +!* write to: Free Software Foundation, Inc., * +!* 675 Mass Ave, Cambridge, MA 02139, USA. * +!* or see: http://www.gnu.org/licenses/gpl.html * +!*********************************************************************** + +use MOM_dyn_horgrid, only : dyn_horgrid_type +use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, 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, OBC_SIMPLE +use MOM_open_boundary, only : open_boundary_query +use MOM_verticalGrid, only : verticalGrid_type +use MOM_time_manager, only : time_type, set_time, time_type_to_real + +implicit none ; private + +#include + +public supercritical_initialize_topography +public supercritical_initialize_velocity +public supercritical_set_OBC_data + +contains + +! ----------------------------------------------------------------------------- +!> This subroutine sets up the supercritical topography +subroutine supercritical_initialize_topography(D, G, param_file, max_depth) + type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type + real, dimension(G%isd:G%ied,G%jsd:G%jed), & + intent(out) :: D !< Ocean bottom depth in m + type(param_file_type), intent(in) :: param_file !< Parameter file structure + real, intent(in) :: max_depth !< Maximum depth of model in m + + real :: min_depth ! The minimum and maximum depths in m. + real :: PI +! This include declares and sets the variable "version". +#include "version_variable.h" + character(len=40) :: mod = "supercritical_initialize_topography" ! This subroutine's name. + integer :: i, j, is, ie, js, je, isd, ied, jsd, jed + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + PI = 4.0*atan(1.0) ; + + call MOM_mesg(" supercritical_initialization.F90, supercritical_initialize_topography: setting topography", 5) + + call log_version(param_file, mod, version, "") + call get_param(param_file, mod, "MINIMUM_DEPTH", min_depth, & + "The minimum depth of the ocean.", units="m", default=0.0) + + do j=js,je ; do i=is,ie + D(i,j)=max_depth + if ((G%geoLonT(i,j) > 10.0).AND. & + (atan2(G%geoLatT(i,j),G%geoLonT(i,j)-10.0) < 8.95*PI/180.)) then + D(i,j)=0.5*min_depth + endif + + if (D(i,j) > max_depth) D(i,j) = max_depth + if (D(i,j) < min_depth) D(i,j) = 0.5*min_depth + enddo ; enddo + +end subroutine supercritical_initialize_topography +! ----------------------------------------------------------------------------- +!> Initialization of u and v in the supercritical test +subroutine supercritical_initialize_velocity(u, v, h, G) + type(ocean_grid_type), intent(in) :: G !< Grid structure + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: u !< i-component of velocity [m/s] + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: v !< j-component of velocity [m/s] + real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: h !< Thickness [H] + + real :: y ! Non-dimensional coordinate across channel, 0..pi + integer :: i, j, k, is, ie, js, je, nz + character(len=40) :: verticalCoordinate + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + + v(:,:,:) = 0.0 + + do j = G%jsc,G%jec ; do I = G%isc-1,G%iec+1 + do k = 1, nz + u(I,j,k) = 8.57 * G%mask2dCu(I,j) ! Thermal wind starting at base of ML + enddo + enddo ; enddo + +end subroutine supercritical_initialize_velocity +! ----------------------------------------------------------------------------- +!> This subroutine sets the properties of flow at open boundary conditions. +subroutine supercritical_set_OBC_data(OBC, G) + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies + !! whether, where, and what open boundary + !! conditions are used. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + + ! The following variables are used to set up the transport in the TIDAL_BAY example. + character(len=40) :: mod = "supercritical_set_OBC_data" ! This subroutine's name. + integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz + integer :: IsdB, IedB, JsdB, JedB + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + + if (.not.associated(OBC)) return + + if (OBC%apply_OBC_u) then + allocate(OBC%u(IsdB:IedB,jsd:jed,nz)) ; OBC%u(:,:,:) = 0.0 + allocate(OBC%uh(IsdB:IedB,jsd:jed,nz)) ; OBC%uh(:,:,:) = 0.0 + endif + if (OBC%apply_OBC_v) then + allocate(OBC%v(isd:ied,JsdB:JedB,nz)) ; OBC%v(:,:,:) = 0.0 + allocate(OBC%vh(isd:ied,JsdB:JedB,nz)) ; OBC%vh(:,:,:) = 0.0 + endif + + do k=1,nz + do j=jsd,jed ; do I=IsdB,IedB + if (OBC%OBC_mask_u(I,j) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_u(I,j))%specified)) then + OBC%u(I,j,k) = 8.57 + OBC%uh(I,j,k) = 8.57 + endif + enddo ; enddo + do J=JsdB,JedB ; do i=isd,ied + if (OBC%OBC_mask_v(i,J) .and. & + (OBC%OBC_segment_list(OBC%OBC_segment_v(i,J))%specified)) then + OBC%v(i,J,k) = 0.0 + OBC%vh(i,J,k) = 0.0 + endif + enddo ; enddo + enddo + +end subroutine supercritical_set_OBC_data + +!> \class supercritical_initialization +!! +!! The module configures the model for the "supercritical" experiment. +!! https://marine.rutgers.edu/po/index.php?model=test-problems&title=supercritical +end module supercritical_initialization diff --git a/src/user/tidal_bay_initialization.F90 b/src/user/tidal_bay_initialization.F90 new file mode 100644 index 0000000000..8974987af6 --- /dev/null +++ b/src/user/tidal_bay_initialization.F90 @@ -0,0 +1,99 @@ +module tidal_bay_initialization +!*********************************************************************** +!* GNU General Public License * +!* This file is a part of MOM. * +!* * +!* MOM is free software; you can redistribute it and/or modify it and * +!* are expected to follow the terms of the GNU General Public License * +!* as published by the Free Software Foundation; either version 2 of * +!* the License, or (at your option) any later version. * +!* * +!* MOM is distributed in the hope that it will be useful, but WITHOUT * +!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * +!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * +!* License for more details. * +!* * +!* For the full text of the GNU General Public License, * +!* write to: Free Software Foundation, Inc., * +!* 675 Mass Ave, Cambridge, MA 02139, USA. * +!* or see: http://www.gnu.org/licenses/gpl.html * +!*********************************************************************** + +use MOM_dyn_horgrid, only : dyn_horgrid_type +use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, 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, OBC_SIMPLE +use MOM_open_boundary, only : open_boundary_query +use MOM_verticalGrid, only : verticalGrid_type +use MOM_time_manager, only : time_type, set_time, time_type_to_real + +implicit none ; private + +#include + +public tidal_bay_set_OBC_data + +contains + +!> This subroutine sets the properties of flow at open boundary conditions. +subroutine tidal_bay_set_OBC_data(OBC, G, h, Time) + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies + !! whether, where, and what open boundary + !! conditions are used. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< layer thickness. + type(time_type), intent(in) :: Time !< model time. + + logical :: apply_OBC_u, apply_OBC_v + ! The following variables are used to set up the transport in the tidal_bay example. + real :: time_sec, cff, cff2, tide_flow + real :: my_area, my_flux + real :: PI + character(len=40) :: mod = "tidal_bay_set_OBC_data" ! This subroutine's name. + integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz + integer :: IsdB, IedB, JsdB, JedB + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + + PI = 4.0*atan(1.0) ; + + if (.not.associated(OBC)) return + + time_sec = time_type_to_real(Time) + cff = 0.1*sin(2.0*PI*time_sec/(12.0*3600.0)) + tide_flow = 3.0e6 + my_area=0.0 + my_flux=0.0 + do j=jsd,jed ; do I=IsdB,IedB + if (OBC%OBC_mask_u(I,j)) then + do k=1,nz + cff2 = h(I,j,k)*G%dyCu(I,j) + my_area = my_area + cff2 + enddo + endif + enddo ; enddo + my_flux = -tide_flow*SIN(2.0*PI*time_sec/(12.0*3600.0)) + + do j=jsd,jed ; do I=IsdB,IedB + if (OBC%OBC_mask_u(I,j)) then + OBC%eta_outer_u(I,j) = cff + OBC%ubt_outer(I,j) = my_flux/my_area + endif + enddo ; enddo + do J=JsdB,JedB ; do i=isd,ied + if (OBC%OBC_mask_v(i,J)) then + OBC%eta_outer_v(i,J) = cff + OBC%vbt_outer(i,J) = 0.0 + endif + enddo ; enddo + +end subroutine tidal_bay_set_OBC_data + +!> \class tidal_bay_Initialization +!! +!! The module configures the model for the "tidal_bay" experiment. +!! tidal_bay = Tidally resonant bay from Zygmunt Kowalik's class on tides. +end module tidal_bay_initialization diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90 index c0868ef543..d030e5ec08 100644 --- a/src/user/user_initialization.F90 +++ b/src/user/user_initialization.F90 @@ -41,7 +41,7 @@ module user_initialization public USER_set_coord, USER_initialize_topography, USER_initialize_thickness public USER_initialize_velocity, USER_init_temperature_salinity public USER_init_mixed_layer_density, USER_initialize_sponges -public USER_set_OBC_positions, USER_set_OBC_data, USER_set_rotation +public USER_set_OBC_data, USER_set_rotation logical :: first_call = .true. @@ -199,23 +199,6 @@ subroutine USER_initialize_sponges(G, use_temperature, tv, param_file, CSp, h) end subroutine USER_initialize_sponges -!> This subroutine sets the location of open boundaries. -subroutine USER_set_OBC_positions(G, param_file, OBC) - type(dyn_horgrid_type), intent(in) :: G !< The ocean's grid structure. - type(param_file_type), intent(in) :: param_file !< A structure indicating the - !! open file to parse for model - !! parameter values. - type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies - !! whether, where, and what open boundary - !! conditions are used. -! call MOM_error(FATAL, & -! "USER_initialization.F90, USER_set_OBC_positions: " // & -! "Unmodified user routine called - you must edit the routine to use it") - - if (first_call) call write_user_log(param_file) - -end subroutine USER_set_OBC_positions - !> This subroutine sets the properties of flow at open boundary conditions. subroutine USER_set_OBC_data(OBC, tv, G, param_file, tr_Reg) type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies