diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index ba9e2941ec..23a0c5712c 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3,7 +3,8 @@ module MOM ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_variables, only : vertvisc_type, ocean_OBC_type +use MOM_variables, only : vertvisc_type +use MOM_open_boundary, only : ocean_OBC_type ! A Structure with pointers to forcing fields to drive MOM; ! all fluxes are positive downward. @@ -98,7 +99,7 @@ module MOM use MOM_mixed_layer_restrat, only : mixedlayer_restrat_register_restarts use MOM_neutral_diffusion, only : neutral_diffusion_CS, neutral_diffusion_diag_init use MOM_obsolete_diagnostics, only : register_obsolete_diagnostics -use MOM_open_boundary, only : Radiation_Open_Bdry_Conds, open_boundary_init +use MOM_open_boundary, only : Radiation_Open_Bdry_Conds use MOM_PressureForce, only : PressureForce, PressureForce_init, PressureForce_CS use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_init use MOM_set_visc, only : set_visc_register_restarts, set_visc_CS @@ -1821,7 +1822,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call callTree_waypoint("restart registration complete (initialize_MOM)") call cpu_clock_begin(id_clock_MOM_init) - call MOM_initialize_fixed(dG, param_file, write_geom_files, dirs%output_directory) + call MOM_initialize_fixed(dG, CS%OBC, param_file, write_geom_files, dirs%output_directory) call callTree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)") call MOM_initialize_coord(GV, param_file, write_geom_files, & dirs%output_directory, CS%tv, dG%max_depth) @@ -1952,7 +1953,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in) call wave_speed_init(Time, G, param_file, diag, CS%wave_speed_CSp) call VarMix_init(Time, G, param_file, diag, CS%VarMix, CS%wave_speed_CSp) call set_visc_init(Time, G, GV, param_file, diag, CS%visc, CS%set_visc_CSp) - if (CS%split) then allocate(eta(SZI_(G),SZJ_(G))) ; eta(:,:) = 0.0 if (CS%legacy_split) then diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index ad640b28c0..3fb030c3ce 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -105,13 +105,13 @@ module MOM_barotropic use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_io, only : vardesc, var_desc +use MOM_open_boundary, only : ocean_OBC_type, OBC_SIMPLE, OBC_NONE, OBC_FLATHER +use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W +use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS use MOM_tidal_forcing, only : tidal_forcing_sensitivity, tidal_forcing_CS use MOM_time_manager, only : time_type, set_time, operator(+), operator(-) use MOM_variables, only : BT_cont_type, alloc_bt_cont_type -use MOM_variables, only : ocean_OBC_type, OBC_SIMPLE, OBC_NONE -use MOM_variables, only : OBC_FLATHER_E, OBC_FLATHER_W -use MOM_variables, only : OBC_FLATHER_N, OBC_FLATHER_S use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -362,6 +362,8 @@ module MOM_barotropic OBC_mask_u => NULL(), & OBC_mask_v => NULL() integer, dimension(:,:), pointer :: & + OBC_direction_u => NULL(), & + OBC_direction_v => NULL(), & OBC_kind_u => NULL(), & OBC_kind_v => NULL() real, dimension(:,:), pointer :: & @@ -2378,44 +2380,46 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, 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_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 - ! 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) - vel_prev = ubt(I,j) - ubt(I,j) = 0.5*((u_inlet + BT_OBC%ubt_outer(I,j)) + & - (BT_OBC%Cg_u(I,j)/H_u) * (h_in-BT_OBC%eta_outer_u(I,j))) - - vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(I,j) - elseif (BT_OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) 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 - ! 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) - vel_prev = ubt(I,j) - ubt(I,j) = 0.5*((u_inlet+BT_OBC%ubt_outer(I,j)) + & - (BT_OBC%Cg_u(I,j)/H_u) * (BT_OBC%eta_outer_u(I,j)-h_in)) - - vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(I,j) - elseif (BT_OBC%OBC_kind_u(I,j) == OBC_FLATHER_N) then - if ((vbt(i,J-1)+vbt(i+1,J-1)) > 0.0) then - ubt(I,j) = 2.0*ubt(I,j-1)-ubt(I,j-2) - else - ubt(I,j) = BT_OBC%ubt_outer(I,j) - endif - vel_trans = ubt(I,j) - elseif (BT_OBC%OBC_kind_u(I,j) == OBC_FLATHER_S) then - if ((vbt(i,J)+vbt(i+1,J)) > 0.0) then - ubt(I,j) = 2.0*ubt(I,j+1)-ubt(I,j+2) - else - ubt(I,j) = BT_OBC%ubt_outer(I,j) + elseif (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_old(I-1,j) + (1.0-cfl)*ubt_old(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) + vel_prev = ubt(I,j) + ubt(I,j) = 0.5*((u_inlet + BT_OBC%ubt_outer(I,j)) + & + (BT_OBC%Cg_u(I,j)/H_u) * (h_in-BT_OBC%eta_outer_u(I,j))) + + vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(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_old(I+1,j) + (1.0-cfl)*ubt_old(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) + vel_prev = ubt(I,j) + ubt(I,j) = 0.5*((u_inlet+BT_OBC%ubt_outer(I,j)) + & + (BT_OBC%Cg_u(I,j)/H_u) * (BT_OBC%eta_outer_u(I,j)-h_in)) + + vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(I,j) + elseif (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_N) then + if ((vbt(i,J-1)+vbt(i+1,J-1)) > 0.0) then + ubt(I,j) = 2.0*ubt(I,j-1)-ubt(I,j-2) + else + ubt(I,j) = BT_OBC%ubt_outer(I,j) + endif + vel_trans = ubt(I,j) + elseif (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_S) then + if ((vbt(i,J)+vbt(i+1,J)) > 0.0) then + ubt(I,j) = 2.0*ubt(I,j+1)-ubt(I,j+2) + else + ubt(I,j) = BT_OBC%ubt_outer(I,j) + endif + vel_trans = ubt(I,j) endif - vel_trans = ubt(I,j) endif if (BT_OBC%OBC_kind_u(I,j) /= OBC_SIMPLE) then @@ -2436,52 +2440,54 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, 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_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 - ! 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) - vel_prev = vbt(i,J) - vbt(i,J) = 0.5*((v_inlet+BT_OBC%vbt_outer(i,J)) + & - (BT_OBC%Cg_v(i,J)/H_v) * (h_in-BT_OBC%eta_outer_v(i,J))) - - vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,J) - elseif (BT_OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) 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 - ! 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) - vel_prev = vbt(i,J) - vbt(i,J) = 0.5*((v_inlet+BT_OBC%vbt_outer(i,J)) + & - (BT_OBC%Cg_v(i,J)/H_v) * (BT_OBC%eta_outer_v(i,J)-h_in)) - - vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,J) - elseif (BT_OBC%OBC_kind_v(i,J) == OBC_FLATHER_E) then - if ((ubt(I-1,j)+ubt(I-1,j+1)) > 0.0) then - vbt(i,J) = 2.0*vbt(i-1,J)-vbt(i-2,J) - else - vbt(i,J) = BT_OBC%vbt_outer(i,J) - endif - vel_trans = vbt(i,J) + elseif (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_old(i,J-1) + (1.0-cfl)*vbt_old(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) + vel_prev = vbt(i,J) + vbt(i,J) = 0.5*((v_inlet+BT_OBC%vbt_outer(i,J)) + & + (BT_OBC%Cg_v(i,J)/H_v) * (h_in-BT_OBC%eta_outer_v(i,J))) + + vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(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_old(i,J+1) + (1.0-cfl)*vbt_old(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) + vel_prev = vbt(i,J) + vbt(i,J) = 0.5*((v_inlet+BT_OBC%vbt_outer(i,J)) + & + (BT_OBC%Cg_v(i,J)/H_v) * (BT_OBC%eta_outer_v(i,J)-h_in)) + + vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,J) + elseif (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_E) then + if ((ubt(I-1,j)+ubt(I-1,j+1)) > 0.0) then + vbt(i,J) = 2.0*vbt(i-1,J)-vbt(i-2,J) + else + vbt(i,J) = BT_OBC%vbt_outer(i,J) + endif + vel_trans = vbt(i,J) !!!!!!!!!!!!!!!!!!! CLAMPED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(i,J) ! ! vbt(i,J) = (vbt(i-1,J) + CFL*vbt(i,J)) / (1.0 + CFL) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - elseif (BT_OBC%OBC_kind_v(i,J) == OBC_FLATHER_W) then - if ((ubt(I,j)+ubt(I,j+1)) < 0.0) then - vbt(i,J) = 2.0*vbt(i+1,J)-vbt(i+2,J) - else - vbt(i,J) = BT_OBC%vbt_outer(i,J) - endif - vel_trans = vbt(i,J) + elseif (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_W) then + if ((ubt(I,j)+ubt(I,j+1)) < 0.0) then + vbt(i,J) = 2.0*vbt(i+1,J)-vbt(i+2,J) + else + vbt(i,J) = BT_OBC%vbt_outer(i,J) + endif + vel_trans = vbt(i,J) !!!!!!!!!!!!!!!!!! CLAMPED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(i,J) ! ! vbt(i,J) = (vbt(i-1,J) + CFL*vbt(i,J)) / (1.0 + CFL) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + endif endif if (BT_OBC%OBC_kind_v(i,J) /= OBC_SIMPLE) then @@ -2532,8 +2538,8 @@ 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_E) 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 @@ -2542,7 +2548,7 @@ subroutine apply_eta_OBCs(OBC, eta, ubt, vbt, BT_OBC, G, MS, halo, dtbt) 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_kind_u(I,j) == OBC_FLATHER_W) then + 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 @@ -2557,8 +2563,8 @@ 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_N) 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 @@ -2567,7 +2573,7 @@ subroutine apply_eta_OBCs(OBC, eta, ubt, vbt, BT_OBC, G, MS, halo, dtbt) 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_kind_v(i,J) == OBC_FLATHER_S) then + 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 @@ -2635,6 +2641,7 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D 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 allocate(BT_OBC%H_v(isdw:iedw,jsdw-1:jedw)) ; BT_OBC%H_v(:,:) = 0.0 @@ -2643,11 +2650,13 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D 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 do k=1,nz ; do j=js,je ; do I=is-1,ie @@ -2683,6 +2692,7 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D 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 do k=1,nz ; do J=js-1,je ; do i=is,ie @@ -2732,6 +2742,7 @@ subroutine destroy_BT_OBC(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) deallocate(BT_OBC%uhbt) @@ -2740,6 +2751,7 @@ subroutine destroy_BT_OBC(BT_OBC) 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) deallocate(BT_OBC%vhbt) diff --git a/src/core/MOM_continuity.F90 b/src/core/MOM_continuity.F90 index 3c378e02f1..912dde87a7 100644 --- a/src/core/MOM_continuity.F90 +++ b/src/core/MOM_continuity.F90 @@ -47,7 +47,8 @@ module MOM_continuity use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_string_functions, only : uppercase use MOM_grid, only : ocean_grid_type -use MOM_variables, only : ocean_OBC_type, BT_cont_type +use MOM_open_boundary, only : ocean_OBC_type +use MOM_variables, only : BT_cont_type use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -78,8 +79,8 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, & uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont) type(ocean_grid_type), intent(inout) :: G type(verticalGrid_type), intent(in) :: GV - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hin real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: uh diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index b479e65682..6529de4474 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -19,37 +19,14 @@ module MOM_continuity_PPM !* or see: http://www.gnu.org/licenses/gpl.html * !*********************************************************************** -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg and Alistair Adcroft, September 2006 - . * -!* * -!* This program contains the subroutine that advects layer * -!* thickness. The scheme here uses a Piecewise-Parabolic method with * -!* a positive definite limiter. * -!* * -!* Macros written all in capital letters are defined in MOM_memory.h. * -!* * -!* A small fragment of the grid is shown below: * -!* * -!* j+1 x ^ x ^ x At x: q * -!* j+1 > o > o > At ^: v, vh * -!* j x ^ x ^ x At >: u, uh * -!* j > o > o > At o: h, hin * -!* 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). * -!* * -!********+*********+*********+*********+*********+*********+*********+** - use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : time_type, diag_ctrl 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_variables, only : ocean_OBC_type, BT_cont_type, OBC_SIMPLE -use MOM_variables, only : OBC_FLATHER_E, OBC_FLATHER_W, OBC_FLATHER_N, OBC_FLATHER_S +use MOM_open_boundary, only : ocean_OBC_type, OBC_SIMPLE, OBC_FLATHER +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 implicit none ; private @@ -104,69 +81,57 @@ module MOM_continuity_PPM contains +!> This subroutine time steps the layer thicknesses, using a monotonically +!! limit, directionally split PPM scheme, based on Lin (1994). In the following +!! documentation, H is used for the units of thickness (usually m or kg m-2.) subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC, & visc_rem_u, visc_rem_v, u_cor, v_cor, & uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont) - type(ocean_grid_type), intent(inout) :: G - type(continuity_PPM_CS), pointer :: CS - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hin - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: uh - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: vh - real, intent(in) :: dt - type(verticalGrid_type), intent(in) :: GV - real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt - real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: vhbt - type(ocean_OBC_type), pointer, optional :: OBC - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in), optional :: visc_rem_u - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in), optional :: visc_rem_v - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor - real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt_aux - real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: vhbt_aux - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor_aux - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor_aux - type(BT_cont_type), pointer, optional :: BT_cont -! This subroutine time steps the layer thicknesses, using a monotonically -! limit, directionally split PPM scheme, based on Lin (1994). In the following -! documentation, H is used for the units of thickness (usually m or kg m-2.) - -! Arguments: u - Zonal velocity, in m s-1. -! (in) v - Meridional velocity, in m s-1. -! (in) hin - Initial layer thickness, in H. -! (out) h - Final layer thickness, in H. -! (out) uh - Volume flux through zonal faces = u*h*dy, H m2 s-1. -! (out) vh - Volume flux through meridional faces = v*h*dx, -! in H m2 s-1. -! (in) dt - Time increment in s. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) CS - The control structure returned by a previous call to -! continuity_PPM_init. -! (in, opt) uhbt - The summed volume flux through zonal faces, H m2 s-1. -! (in, opt) vhbt - The summed volume flux through meridional faces, H m2 s-1. -! (in, opt) OBC - This open boundary condition type specifies whether, where, -! and what open boundary conditions are used. -! (in, opt) visc_rem_u - Both the fraction of the momentum originally in a -! (in, opt) visc_rem_v - layer that remains after a time-step of viscosity, -! and the fraction of a time-step's worth of a -! barotropic acceleration that a layer experiences -! after viscosity is applied, in the zonal (_u) and -! meridional (_v) directions. Nondimensional between -! 0 (at the bottom) and 1 (far above the bottom). -! (out, opt) u_cor - The zonal velocities that give uhbt as the depth- -! integrated transport, in m s-1. -! (out, opt) v_cor - The meridional velocities that give vhbt as the -! depth-integrated transport, in m s-1. -! (in, opt) uhbt_aux - A second set of summed volume fluxes through zonal -! (in, opt) vhbt_aux - and meridional faces, both in H m2 s-1. -! (out, opt) u_cor_aux - The zonal and meridional velocities that give uhbt_aux -! (out, opt) v_cor_aux - and vhbt_aux as the depth-integrated transports, -! both in m s-1. -! (out, opt) BT_cont - A structure with elements that describe the effective -! open face areas as a function of barotropic flow. + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. + type(continuity_PPM_CS), pointer :: CS !< Module's control structure. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u !< Zonal velocity, in m s-1. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v !< Meridional velocity, in m s-1. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hin !< Initial layer thickness, in H. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Final layer thickness, in H. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: uh !< Zonal volume flux, + !! u*h*dy, H m2 s-1. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: vh !< Meridional volume flux, + !! v*h*dx, H m2 s-1. + real, intent(in) :: dt !< Time increment in s. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt !< + !! The summed volume flux through zonal faces, H m2 s-1. + real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: vhbt !< + !! The summed volume flux through meridional faces, H m2 s-1. + type(ocean_OBC_type), pointer, optional :: OBC !< + !! This open boundary condition type specifies whether, where, + !! and what open boundary conditions are used. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in), optional :: visc_rem_u !< + !! Both the fraction of the momentum originally in a layer that remains after a time-step + !! of viscosity, and the fraction of a time-step's worth of a barotropic acceleration that + !! a layer experiences after viscosity is applied, in the zonal (_u) + !! direction. Nondimensional between 0 (at the bottom) and 1 (far above the bottom). + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in), optional :: visc_rem_v !< + !! Same as above for meridional direction. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor !< + !! The zonal velocities that give uhbt as the depth- + !! integrated transport, in m s-1. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor !< + !! The meridional velocities that give vhbt as the + !! depth-integrated transport, in m s-1. + real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt_aux !< + !! A second set of summed volume fluxes through zonal faces, in H m2 s-1. + real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: vhbt_aux !< + !! A second set of summed volume fluxes through meridional faces, in H m2 s-1. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor_aux !< + !! The zonal velocities that give uhbt_aux as the + !! depth-integrated transports, in m s-1. + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out), optional :: v_cor_aux !< + !! The meridional velocities that give vhbt_aux as the + !! depth-integrated transports, in m s-1. + type(BT_cont_type), pointer, optional :: BT_cont !< + !! A structure with elements that describe the effective + !! open face areas as a function of barotropic flow. real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & h_input ! Left and right face thicknesses, in H. @@ -191,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 (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_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 @@ -201,7 +166,7 @@ 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 .or. & apply_OBC_v_flather_north .or. apply_OBC_v_flather_south) & h_input(:,:,:) = hin(:,:,:) - endif ; endif + endif ; endif ; endif if (present(visc_rem_u) .neqv. present(visc_rem_v)) call MOM_error(FATAL, & "MOM_continuity_PPM: Either both visc_rem_u and visc_rem_v or neither"// & @@ -228,13 +193,21 @@ 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_mask_u(I-1,j) .and. (OBC%OBC_kind_u(I-1,j) == OBC_FLATHER_E)) & + 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) enddo do i=LB%ish-1,LB%ieh - if (OBC%OBC_mask_u(I,j) .and. (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W)) & + 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) + enddo enddo ; enddo endif LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec @@ -256,13 +229,19 @@ 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_mask_v(i,J-1) .and. (OBC%OBC_kind_v(i,J-1) == OBC_FLATHER_N)) & + 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) enddo ; enddo do J=LB%jsh-1,LB%jeh ; do i=LB%ish-1,LB%ieh+1 - if (OBC%OBC_mask_v(i,J) .and. (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S)) & + 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) + enddo ; enddo enddo endif else ! .not. x_first @@ -283,13 +262,19 @@ 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_mask_v(i,J-1) .and. (OBC%OBC_kind_v(i,J-1) == OBC_FLATHER_N)) & + 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) enddo ; enddo do J=LB%jsh-1,LB%jeh ; do i=LB%ish-1,LB%ieh+1 - if (OBC%OBC_mask_v(i,J) .and. (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S)) & + 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) + enddo ; enddo enddo endif @@ -311,62 +296,65 @@ 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_mask_u(I-1,j) .and. (OBC%OBC_kind_u(I-1,j) == OBC_FLATHER_E)) & + 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) enddo do i=LB%ish-1,LB%ieh - if (OBC%OBC_mask_u(I,j) .and. (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W)) & + 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) + enddo enddo ; enddo endif endif end subroutine continuity_PPM +!> This subroutine calculates the mass or volume fluxes through the zonal +!! faces, and other related quantities. subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & visc_rem_u, u_cor, uhbt_aux, u_cor_aux, BT_cont) - type(ocean_grid_type), intent(inout) :: G - type(verticalGrid_type), intent(in) :: GV - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_in - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: uh - real, intent(in) :: dt - type(continuity_PPM_CS), pointer :: CS - type(loop_bounds_type), intent(in) :: LB - type(ocean_OBC_type), pointer, optional :: OBC - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in), optional :: visc_rem_u - real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt, uhbt_aux - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor, u_cor_aux - type(BT_cont_type), pointer, optional :: BT_cont -! This subroutine calculates the mass or volume fluxes through the zonal -! faces, and other related quantities. -! Arguments: u - Zonal velocity, in m s-1. -! (in) h_in - Layer thickness used to calculate the fluxes, in H. -! (out) uh - Volume flux through zonal faces = u*h*dy, H m2 s-1. -! (in) dt - Time increment in s. -! (in) G - The ocean's grid structure. -! (in) CS - The control structure returned by a previous call to -! continuity_PPM_init. -! (in) LB - A structure with the active loop bounds. -! (in, opt) uhbt - The summed volume flux through zonal faces, H m2 s-1. -! (in, opt) OBC - This open boundary condition type specifies whether, where, -! and what open boundary conditions are used. -! (in, opt) visc_rem_u - Both the fraction of the momentum originally in a -! layer that remains after a time-step of viscosity, -! and the fraction of a time-step's worth of a -! barotropic acceleration that a layer experiences -! after viscosity is applied. Nondimensional between -! 0 (at the bottom) and 1 (far above the bottom). -! (out, opt) u_cor - The zonal velocitiess (u with a barotropic correction) -! that give uhbt as the depth-integrated transport, m s-1. -! (in, opt) uhbt_aux - A second set of summed volume fluxes through zonal -! faces, in H m2 s-1. -! (out, opt) u_cor_aux - The zonal velocities (u with a barotropic correction) -! that give uhbt_aux as the depth-integrated transports, -! in m s-1. -! (out, opt) BT_cont - A structure with elements that describe the effective -! open face areas as a function of barotropic flow. + type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Zonal velocity, in m s-1. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_in !< Layer thickness used to + !! calculate fluxes, in H. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: uh !< Volume flux through zonal + !! faces = u*h*dy, H m2 s-1. + real, intent(in) :: dt !< Time increment in s. + type(continuity_PPM_CS), pointer :: CS !< This module's control structure. + type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure. + type(ocean_OBC_type), pointer, optional :: OBC !< + !! This open boundary condition type specifies whether, where, + !! and what open boundary conditions are used. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in), optional :: visc_rem_u !< + !! Both the fraction of the momentum originally in a + !! layer that remains after a time-step of viscosity, + !! and the fraction of a time-step's worth of a + !! barotropic acceleration that a layer experiences + !! after viscosity is applied. Nondimensional between + !! 0 (at the bottom) and 1 (far above the bottom). + real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt !< + !! The summed volume flux through zonal faces, H m2 s-1. + real, dimension(SZIB_(G),SZJ_(G)), intent(in), optional :: uhbt_aux !< + !! A second set of summed volume fluxes through zonal + !! faces, in H m2 s-1. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor !< + !! The zonal velocitiess (u with a barotropic correction) + !! that give uhbt as the depth-integrated transport, m s-1. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out), optional :: u_cor_aux !< + !! The zonal velocities (u with a barotropic correction) + !! that give uhbt_aux as the depth-integrated transports, in m s-1. + type(BT_cont_type), pointer, optional :: BT_cont !< + !! A structure with elements that describe the effective + !! open face areas as a function of barotropic flow. real, dimension(SZIB_(G),SZK_(G)) :: & duhdu ! Partial derivative of uh with u, in H m. @@ -390,13 +378,18 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & real :: dx_E, dx_W ! Effective x-grid spacings to the east and west, in m. integer :: i, j, k, ish, ieh, jsh, jeh, nz logical :: do_aux, apply_OBC_u, use_visc_rem, set_BT_cont, any_simple_OBC + logical :: apply_OBC_flather do_aux = (present(uhbt_aux) .and. present(u_cor_aux)) use_visc_rem = present(visc_rem_u) - apply_OBC_u = .false. ; set_BT_cont = .false. + apply_OBC_u = .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 apply_OBC_u = OBC%apply_OBC_u + apply_OBC_flather = 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 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = G%ke @@ -529,12 +522,17 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & any_simple_OBC = .false. if (present(uhbt) .or. do_aux .or. set_BT_cont) then - if (.not.apply_OBC_u) then ; do I=ish-1,ieh - do_i(I) = .true. - enddo ; else ; do I=ish-1,ieh + 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)) 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_direction_u(I,j) == OBC_DIRECTION_N) .or. & + (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_S))) + enddo ; else ; do I=ish-1,ieh + do_i(I) = .true. enddo ; endif endif @@ -606,35 +604,31 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, & end subroutine zonal_mass_flux +!> This subroutines evaluates the zonal mass or volume fluxes in a layer. subroutine zonal_flux_layer(u, h, hL, hR, uh, duhdu, visc_rem, dt, G, j, & ish, ieh, do_i, vol_CFL) - type(ocean_grid_type), intent(inout) :: G - real, dimension(SZIB_(G)), intent(in) :: u, visc_rem - real, dimension(SZI_(G)), intent(in) :: h, hL, hR - real, dimension(SZIB_(G)), intent(inout) :: uh, duhdu - real, intent(in) :: dt - integer, intent(in) :: j, ish, ieh - logical, dimension(SZIB_(G)), intent(in) :: do_i - logical, intent(in) :: vol_CFL -! This subroutines evaluates the zonal mass or volume fluxes in a layer. -! -! Arguments: u - Zonal velocity, in m s-1. -! (in) h - Layer thickness used to calculate the fluxes, in H. -! (in) hL, hR - Left- and right- thicknesses in the reconstruction, in H. -! (out) uh - The zonal mass or volume transport, in H m2 s-1. -! (out) duhdu - The partial derivative of uh with u, in H m. -! (in) dt - Time increment in s. -! (in) G - The ocean's grid structure. -! (in) visc_rem - Both the fraction of the momentum originally in a -! layer that remains after a time-step of viscosity, -! and the fraction of a time-step's worth of a -! barotropic acceleration that a layer experiences -! after viscosity is applied. Nondimensional between -! 0 (at the bottom) and 1 (far above the bottom). -! (in) j, ish, ieh - The index range to work on. -! (in) do_i - A logical flag indiciating which I values to work on. -! (in) vol_CFL - If true, rescale the ratio of face areas to the cell -! areas when estimating the CFL number. + type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. + real, dimension(SZIB_(G)), intent(in) :: u !< Zonal velocity, in m s-1. + real, dimension(SZIB_(G)), intent(in) :: visc_rem !< Both the fraction of the + !! momentum originally in a layer that remains after a time-step + !! of viscosity, and the fraction of a time-step's worth of a barotropic + !! acceleration that a layer experiences after viscosity is applied. + !! Nondimensional between 0 (at the bottom) and 1 (far above the bottom). + real, dimension(SZI_(G)), intent(in) :: h !< Layer thickness, in H. + real, dimension(SZI_(G)), intent(in) :: hL !< Left thickness, in H. + real, dimension(SZI_(G)), intent(in) :: hR !< Right thickness, in H. + real, dimension(SZIB_(G)), intent(inout) :: uh !< Zonal mass or volume + !! transport, in H m2 s-1. + real, dimension(SZIB_(G)), intent(inout) :: duhdu !< Partial derivative of uh + !! with u, in H m. + real, intent(in) :: dt !< Time increment in s. + integer, intent(in) :: j !< Spatial index. + integer, intent(in) :: ish !< Start of index range. + integer, intent(in) :: ieh !< End of index range. + logical, dimension(SZIB_(G)), intent(in) :: do_i !< Which i values to work on. + logical, intent(in) :: vol_CFL !< If true, rescale the + !! ratio of face areas to the cell areas when estimating the CFL number. + real :: CFL ! The CFL number based on the local velocity and grid spacing, ND. real :: curv_3 ! A measure of the thickness curvature over a grid length, ! with the same units as h_in. @@ -666,19 +660,21 @@ subroutine zonal_flux_layer(u, h, hL, hR, uh, duhdu, visc_rem, dt, G, j, & end subroutine zonal_flux_layer +!> This subroutines sets the effective interface thickness at each zonal +!! velocity point. subroutine zonal_face_thickness(u, h, hL, hR, h_u, dt, G, LB, vol_CFL, & marginal, visc_rem_u) type(ocean_grid_type), intent(inout) :: G real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h, hL, hR + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hL + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: hR real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: h_u real, intent(in) :: dt type(loop_bounds_type), intent(in) :: LB logical, intent(in) :: vol_CFL logical, intent(in) :: marginal real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in), optional :: visc_rem_u -! This subroutines sets the effective interface thickness at each zonal -! velocity point. ! ! Arguments: u - Zonal velocity, in m s-1. ! (in) h - Layer thickness used to calculate the fluxes, in H. @@ -1147,14 +1143,19 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & real :: dy_N, dy_S ! Effective y-grid spacings to the north and south, in m. integer :: i, j, k, ish, ieh, jsh, jeh, nz logical :: do_aux, apply_OBC_v, use_visc_rem, set_BT_cont, any_simple_OBC + logical :: apply_OBC_flather do_aux = (present(vhbt_aux) .and. present(v_cor_aux)) use_visc_rem = present(visc_rem_v) - apply_OBC_v = .false. ; set_BT_cont = .false. + 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 (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_pe) then apply_OBC_v = OBC%apply_OBC_v - endif ; endif + apply_OBC_flather = 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 ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = G%ke CFL_dt = CS%CFL_limit_adjust / dt @@ -1284,12 +1285,17 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, & any_simple_OBC = .false. if (present(vhbt) .or. do_aux .or. set_BT_cont) then - if (.not.apply_OBC_v) then ; do i=ish,ieh - do_i(i) = .true. - enddo ; else ; do i=ish,ieh + 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)) - if (.not.do_i(I)) any_simple_OBC = .true. + 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_direction_v(i,J) == OBC_DIRECTION_E) .or. & + (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_W))) + enddo ; else ; do i=ish,ieh + do_i(i) = .true. enddo ; endif endif @@ -2213,4 +2219,28 @@ subroutine continuity_PPM_end(CS) deallocate(CS) end subroutine continuity_PPM_end +!********+*********+*********+*********+*********+*********+*********+** +!* * +!* By Robert Hallberg and Alistair Adcroft, September 2006 - . * +!* * +!* This program contains the subroutine that advects layer * +!* thickness. The scheme here uses a Piecewise-Parabolic method with * +!* a positive definite limiter. * +!* * +!* Macros written all in capital letters are defined in MOM_memory.h. * +!* * +!* A small fragment of the grid is shown below: * +!* * +!* j+1 x ^ x ^ x At x: q * +!* j+1 > o > o > At ^: v, vh * +!* j x ^ x ^ x At >: u, uh * +!* j > o > o > At o: h, hin * +!* 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_continuity_PPM diff --git a/src/core/MOM_dynamics_legacy_split.F90 b/src/core/MOM_dynamics_legacy_split.F90 index afb455a2a9..cb9f1d62ff 100644 --- a/src/core/MOM_dynamics_legacy_split.F90 +++ b/src/core/MOM_dynamics_legacy_split.F90 @@ -66,7 +66,8 @@ module MOM_dynamics_legacy_split !********+*********+*********+*********+*********+*********+*********+** -use MOM_variables, only : vertvisc_type, ocean_OBC_type, thermo_var_ptrs +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 use MOM_forcing_type, only : forcing @@ -109,8 +110,7 @@ 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 : Radiation_Open_Bdry_Conds, open_boundary_init -use MOM_open_boundary, only : open_boundary_CS +use MOM_open_boundary, only : Radiation_Open_Bdry_Conds 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 @@ -240,7 +240,6 @@ module MOM_dynamics_legacy_split type(legacy_barotropic_CS), pointer :: barotropic_CSp => NULL() type(vertvisc_CS), pointer :: vertvisc_CSp => NULL() type(set_visc_CS), pointer :: set_visc_CSp => NULL() - type(open_boundary_CS), pointer :: open_boundary_CSp => NULL() type(ocean_OBC_type), pointer :: OBC => NULL() ! A pointer to an open boundary ! condition type that specifies whether, where, and what open boundary ! conditions are used. If no open BCs are used, this pointer stays @@ -738,7 +737,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & if (associated(CS%OBC)) then call Radiation_Open_Bdry_Conds(CS%OBC, u_av, u_old_rad_OBC, v_av, & - v_old_rad_OBC, hp, h_old_rad_OBC, G, CS%open_boundary_CSp) + v_old_rad_OBC, hp, h_old_rad_OBC, G) endif ! h_av = (h + hp)/2 @@ -997,7 +996,7 @@ subroutine step_MOM_dyn_legacy_split(u, v, h, tv, visc, & if (associated(CS%OBC)) then call Radiation_Open_Bdry_Conds(CS%OBC, u, u_old_rad_OBC, v, & - v_old_rad_OBC, h, h_old_rad_OBC, G, CS%open_boundary_CSp) + v_old_rad_OBC, h, h_old_rad_OBC, G) endif ! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in. @@ -1053,8 +1052,8 @@ end subroutine step_MOM_dyn_legacy_split subroutine adjustments_dyn_legacy_split(u, v, h, dt, G, GV, CS) type(ocean_grid_type), intent(inout) :: G type(verticalGrid_type), intent(in) :: GV - real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u - real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h real, intent(in) :: dt type(MOM_dyn_legacy_split_CS), pointer :: CS @@ -1379,10 +1378,7 @@ subroutine initialize_dyn_legacy_split(u, v, h, uh, vh, eta, Time, G, GV, param_ CS%set_visc_CSp => setVisc_CSp if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp - if (associated(OBC)) then - CS%OBC => OBC - call open_boundary_init(Time, G, param_file, diag, CS%open_boundary_CSp) - endif + if (associated(OBC)) CS%OBC => OBC if (.not. query_initialized(CS%eta,"sfc",restart_CS)) then ! Estimate eta based on the layer thicknesses - h. With the Boussinesq diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 2e112bf274..2e2173a502 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -3,7 +3,7 @@ module MOM_dynamics_split_RK2 ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_variables, only : vertvisc_type, ocean_OBC_type, thermo_var_ptrs +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 use MOM_forcing_type, only : forcing @@ -45,8 +45,8 @@ module MOM_dynamics_split_RK2 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 : Radiation_Open_Bdry_Conds, open_boundary_init -use MOM_open_boundary, only : open_boundary_CS +use MOM_open_boundary, only : ocean_OBC_type +use MOM_open_boundary, only : Radiation_Open_Bdry_Conds 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 @@ -160,7 +160,6 @@ module MOM_dynamics_split_RK2 type(barotropic_CS), pointer :: barotropic_CSp => NULL() type(vertvisc_CS), pointer :: vertvisc_CSp => NULL() type(set_visc_CS), pointer :: set_visc_CSp => NULL() - type(open_boundary_CS), pointer :: open_boundary_CSp => NULL() type(tidal_forcing_CS), pointer :: tides_CSp => NULL() type(ocean_OBC_type), pointer :: OBC => NULL() !< A pointer to an open boundary @@ -639,7 +638,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & if (associated(CS%OBC)) then call Radiation_Open_Bdry_Conds(CS%OBC, u_av, u_old_rad_OBC, v_av, & - v_old_rad_OBC, hp, h_old_rad_OBC, G, CS%open_boundary_CSp) + v_old_rad_OBC, hp, h_old_rad_OBC, G) endif ! h_av = (h + hp)/2 @@ -851,7 +850,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, & if (associated(CS%OBC)) then call Radiation_Open_Bdry_Conds(CS%OBC, u, u_old_rad_OBC, v, & - v_old_rad_OBC, h, h_old_rad_OBC, G, CS%open_boundary_CSp) + v_old_rad_OBC, h, h_old_rad_OBC, G) endif ! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in. @@ -1140,10 +1139,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, param_fil (LEN_TRIM(dirs%input_filename) == 1)) ) if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp - if (associated(OBC)) then - CS%OBC => OBC - call open_boundary_init(Time, G, param_file, diag, CS%open_boundary_CSp) - endif + if (associated(OBC)) CS%OBC => OBC if (.not. query_initialized(CS%eta,"sfc",restart_CS)) then ! Estimate eta based on the layer thicknesses - h. With the Boussinesq diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 7b4a4d4960..fc392cbb88 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -68,7 +68,7 @@ module MOM_dynamics_unsplit !********+*********+*********+*********+*********+*********+*********+** -use MOM_variables, only : vertvisc_type, ocean_OBC_type, thermo_var_ptrs +use MOM_variables, only : vertvisc_type, thermo_var_ptrs use MOM_variables, only : accel_diag_ptrs, ocean_internal_state, cont_diag_ptrs use MOM_forcing_type, only : forcing use MOM_checksum_packages, only : MOM_thermo_chksum, MOM_state_chksum, MOM_accel_chksum @@ -102,8 +102,8 @@ module MOM_dynamics_unsplit 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 : Radiation_Open_Bdry_Conds, open_boundary_init -use MOM_open_boundary, only : open_boundary_CS +use MOM_open_boundary, only : ocean_OBC_type +use MOM_open_boundary, only : Radiation_Open_Bdry_Conds 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 @@ -155,7 +155,6 @@ module MOM_dynamics_unsplit type(PressureForce_CS), pointer :: PressureForce_CSp => NULL() type(vertvisc_CS), pointer :: vertvisc_CSp => NULL() type(set_visc_CS), pointer :: set_visc_CSp => NULL() - type(open_boundary_CS), pointer :: open_boundary_CSp => NULL() type(ocean_OBC_type), pointer :: OBC => NULL() ! A pointer to an open boundary ! condition type that specifies whether, where, and what open boundary ! conditions are used. If no open BCs are used, this pointer stays @@ -679,10 +678,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, param_file, diag, CS, & CS%set_visc_CSp => setVisc_CSp if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp - if (associated(OBC)) then - CS%OBC => OBC - call open_boundary_init(Time, G, param_file, diag, CS%open_boundary_CSp) - endif + if (associated(OBC)) CS%OBC => OBC flux_units = get_flux_units(GV) CS%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, Time, & diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 72e4f56613..8e228a9189 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -65,7 +65,7 @@ module MOM_dynamics_unsplit_RK2 !* * !********+*********+*********+*********+*********+*********+*********+** -use MOM_variables, only : vertvisc_type, ocean_OBC_type, thermo_var_ptrs +use MOM_variables, only : vertvisc_type, thermo_var_ptrs use MOM_variables, only : ocean_internal_state, accel_diag_ptrs, cont_diag_ptrs use MOM_forcing_type, only : forcing use MOM_checksum_packages, only : MOM_thermo_chksum, MOM_state_chksum, MOM_accel_chksum @@ -100,8 +100,8 @@ module MOM_dynamics_unsplit_RK2 use MOM_hor_visc, only : horizontal_viscosity, hor_visc_init, hor_visc_CS use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_MEKE_types, only : MEKE_type -use MOM_open_boundary, only : Radiation_Open_Bdry_Conds, open_boundary_init -use MOM_open_boundary, only : open_boundary_CS +use MOM_open_boundary, only : ocean_OBC_type +use MOM_open_boundary, only : Radiation_Open_Bdry_Conds 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 @@ -161,7 +161,6 @@ module MOM_dynamics_unsplit_RK2 type(PressureForce_CS), pointer :: PressureForce_CSp => NULL() type(vertvisc_CS), pointer :: vertvisc_CSp => NULL() type(set_visc_CS), pointer :: set_visc_CSp => NULL() - type(open_boundary_CS), pointer :: open_boundary_CSp => NULL() type(ocean_OBC_type), pointer :: OBC => NULL() ! A pointer to an open boundary ! condition type that specifies whether, where, and what open boundary ! conditions are used. If no open BCs are used, this pointer stays @@ -643,10 +642,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, param_file, diag, CS CS%set_visc_CSp => setVisc_CSp if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp - if (associated(OBC)) then - CS%OBC => OBC - call open_boundary_init(Time, G, param_file, diag, CS%open_boundary_CSp) - endif + if (associated(OBC)) CS%OBC => OBC flux_units = get_flux_units(GV) CS%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, Time, & diff --git a/src/core/MOM_legacy_barotropic.F90 b/src/core/MOM_legacy_barotropic.F90 index 0b15fb6f3f..d8232d1b77 100644 --- a/src/core/MOM_legacy_barotropic.F90 +++ b/src/core/MOM_legacy_barotropic.F90 @@ -107,13 +107,13 @@ module MOM_legacy_barotropic use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_io, only : vardesc, var_desc +use MOM_open_boundary, only : ocean_OBC_type, OBC_SIMPLE, OBC_NONE, OBC_FLATHER +use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W +use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS use MOM_tidal_forcing, only : tidal_forcing_sensitivity, tidal_forcing_CS use MOM_time_manager, only : time_type, set_time, operator(+), operator(-) use MOM_variables, only : BT_cont_type, alloc_bt_cont_type -use MOM_variables, only : ocean_OBC_type, OBC_SIMPLE, OBC_NONE -use MOM_variables, only : OBC_FLATHER_E, OBC_FLATHER_W -use MOM_variables, only : OBC_FLATHER_N, OBC_FLATHER_S use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -352,6 +352,8 @@ module MOM_legacy_barotropic OBC_mask_u => NULL(), & OBC_mask_v => NULL() integer, dimension(:,:), pointer :: & + OBC_direction_u => NULL(), & + OBC_direction_v => NULL(), & OBC_kind_u => NULL(), & OBC_kind_v => NULL() real, dimension(:,:), pointer :: & @@ -2236,44 +2238,46 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, 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_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 - ! 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 + elseif (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_old(I-1,j) + (1.0-cfl)*ubt_old(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) + vel_prev = ubt(I,j) + ubt(I,j) = 0.5*((u_inlet + BT_OBC%ubt_outer(I,j)) + & + (BT_OBC%Cg_u(I,j)/H_u) * (h_in-BT_OBC%eta_outer_u(I,j))) - H_u = BT_OBC%H_u(I,j) - vel_prev = ubt(I,j) - ubt(I,j) = 0.5*((u_inlet + BT_OBC%ubt_outer(I,j)) + & - (BT_OBC%Cg_u(I,j)/H_u) * (h_in-BT_OBC%eta_outer_u(I,j))) - - vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(I,j) - elseif (BT_OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) 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 - ! 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 + vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(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_old(I+1,j) + (1.0-cfl)*ubt_old(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) - vel_prev = ubt(I,j) - ubt(I,j) = 0.5*((u_inlet+BT_OBC%ubt_outer(I,j)) + & - (BT_OBC%Cg_u(I,j)/H_u) * (BT_OBC%eta_outer_u(I,j)-h_in)) + H_u = BT_OBC%H_u(I,j) + vel_prev = ubt(I,j) + ubt(I,j) = 0.5*((u_inlet+BT_OBC%ubt_outer(I,j)) + & + (BT_OBC%Cg_u(I,j)/H_u) * (BT_OBC%eta_outer_u(I,j)-h_in)) - vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(I,j) - elseif (BT_OBC%OBC_kind_u(I,j) == OBC_FLATHER_N) then - if ((vbt(i,J-1)+vbt(i+1,J-1)) > 0.0) then - ubt(I,j) = 2.0*ubt(I,j-1)-ubt(I,j-2) - else - ubt(I,j) = BT_OBC%ubt_outer(I,j) - endif - vel_trans = ubt(I,j) - elseif (BT_OBC%OBC_kind_u(I,j) == OBC_FLATHER_S) then - if ((vbt(i,J)+vbt(i+1,J)) > 0.0) then - ubt(I,j) = 2.0*ubt(I,j+1)-ubt(I,j+2) - else - ubt(I,j) = BT_OBC%ubt_outer(I,j) + vel_trans = (1.0-bebt)*vel_prev + bebt*ubt(I,j) + elseif (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_N) then + if ((vbt(i,J-1)+vbt(i+1,J-1)) > 0.0) then + ubt(I,j) = 2.0*ubt(I,j-1)-ubt(I,j-2) + else + ubt(I,j) = BT_OBC%ubt_outer(I,j) + endif + vel_trans = ubt(I,j) + elseif (BT_OBC%OBC_direction_u(I,j) == OBC_DIRECTION_S) then + if ((vbt(i,J)+vbt(i+1,J)) > 0.0) then + ubt(I,j) = 2.0*ubt(I,j+1)-ubt(I,j+2) + else + ubt(I,j) = BT_OBC%ubt_outer(I,j) + endif + vel_trans = ubt(I,j) endif - vel_trans = ubt(I,j) endif if (BT_OBC%OBC_kind_u(I,j) /= OBC_SIMPLE) then @@ -2294,52 +2298,54 @@ subroutine apply_velocity_OBCs(OBC, ubt, vbt, uhbt, vhbt, ubt_trans, vbt_trans, 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_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 - ! 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) - vel_prev = vbt(i,J) - vbt(i,J) = 0.5*((v_inlet+BT_OBC%vbt_outer(i,J)) + & - (BT_OBC%Cg_v(i,J)/H_v) * (h_in-BT_OBC%eta_outer_v(i,J))) + elseif (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_old(i,J-1) + (1.0-cfl)*vbt_old(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) + vel_prev = vbt(i,J) + vbt(i,J) = 0.5*((v_inlet+BT_OBC%vbt_outer(i,J)) + & + (BT_OBC%Cg_v(i,J)/H_v) * (h_in-BT_OBC%eta_outer_v(i,J))) - vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,J) - elseif (BT_OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) 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 - ! 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 + vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(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_old(i,J+1) + (1.0-cfl)*vbt_old(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) - vel_prev = vbt(i,J) - vbt(i,J) = 0.5*((v_inlet+BT_OBC%vbt_outer(i,J)) + & - (BT_OBC%Cg_v(i,J)/H_v) * (BT_OBC%eta_outer_v(i,J)-h_in)) + H_v = BT_OBC%H_v(i,J) + vel_prev = vbt(i,J) + vbt(i,J) = 0.5*((v_inlet+BT_OBC%vbt_outer(i,J)) + & + (BT_OBC%Cg_v(i,J)/H_v) * (BT_OBC%eta_outer_v(i,J)-h_in)) - vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,J) - elseif (BT_OBC%OBC_kind_v(i,J) == OBC_FLATHER_E) then - if ((ubt(I-1,j)+ubt(I-1,j+1)) > 0.0) then - vbt(i,J) = 2.0*vbt(i-1,J)-vbt(i-2,J) - else - vbt(i,J) = BT_OBC%vbt_outer(i,J) - endif - vel_trans = vbt(i,J) + vel_trans = (1.0-bebt)*vel_prev + bebt*vbt(i,J) + elseif (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_E) then + if ((ubt(I-1,j)+ubt(I-1,j+1)) > 0.0) then + vbt(i,J) = 2.0*vbt(i-1,J)-vbt(i-2,J) + else + vbt(i,J) = BT_OBC%vbt_outer(i,J) + endif + vel_trans = vbt(i,J) !!!!!!!!!!!!!!!!!!! CLAMPED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(i,J) ! ! vbt(i,J) = (vbt(i-1,J) + CFL*vbt(i,J)) / (1.0 + CFL) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - elseif (BT_OBC%OBC_kind_v(i,J) == OBC_FLATHER_W) then - if ((ubt(I,j)+ubt(I,j+1)) < 0.0) then - vbt(i,J) = 2.0*vbt(i+1,J)-vbt(i+2,J) - else - vbt(i,J) = BT_OBC%vbt_outer(i,J) - endif - vel_trans = vbt(i,J) + elseif (BT_OBC%OBC_direction_v(i,J) == OBC_DIRECTION_W) then + if ((ubt(I,j)+ubt(I,j+1)) < 0.0) then + vbt(i,J) = 2.0*vbt(i+1,J)-vbt(i+2,J) + else + vbt(i,J) = BT_OBC%vbt_outer(i,J) + endif + vel_trans = vbt(i,J) !!!!!!!!!!!!!!!!!! CLAMPED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! cfl = dtbt * BT_OBC%Cg_v(i,J) * G%IdyCv(i,J) ! ! vbt(i,J) = (vbt(i-1,J) + CFL*vbt(i,J)) / (1.0 + CFL) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + endif endif if (BT_OBC%OBC_kind_v(i,J) /= OBC_SIMPLE) then @@ -2391,24 +2397,26 @@ 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_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_kind_u(I,j) == OBC_FLATHER_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) + 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) + endif endif endif ; enddo ; enddo endif @@ -2416,24 +2424,26 @@ 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_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_kind_v(i,J) == OBC_FLATHER_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) + 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) + endif endif endif ; enddo ; enddo endif @@ -2493,6 +2503,7 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D 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 allocate(BT_OBC%H_v(isdw:iedw,jsdw-1:jedw)) ; BT_OBC%H_v(:,:) = 0.0 @@ -2501,11 +2512,13 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D 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 do k=1,nz ; do j=js,je ; do I=is-1,ie @@ -2541,6 +2554,7 @@ subroutine set_up_BT_OBC(OBC, eta, BT_OBC, G, GV, MS, halo, use_BT_cont, Datu, D 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 do k=1,nz ; do J=js-1,je ; do i=is,ie @@ -2590,6 +2604,7 @@ subroutine destroy_BT_OBC(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) deallocate(BT_OBC%uhbt) @@ -2598,6 +2613,7 @@ subroutine destroy_BT_OBC(BT_OBC) 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) deallocate(BT_OBC%vhbt) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index c3b93add73..d1bd34d6ad 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -1,90 +1,345 @@ +! This file is part of MOM6. See LICENSE.md for the license. +!> Controls where open boundary conditions are applied module MOM_open_boundary -!*********************************************************************** -!* 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 * -!*********************************************************************** - -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Mehmet Ilicak and Robert Hallberg, 2010 * -!* * -!* This module implements some aspects of internal open boundary * -!* conditions in MOM. * -!* * -!* 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). * -!* * -!********+*********+*********+*********+*********+*********+*********+** + +! 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 +use MOM_file_parser, only : get_param, log_version, param_file_type, log_param use MOM_grid, only : ocean_grid_type -use MOM_variables, only : ocean_OBC_type, OBC_NONE, OBC_SIMPLE -use MOM_variables, only : OBC_FLATHER_E, OBC_FLATHER_W, OBC_FLATHER_N, OBC_FLATHER_S +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_tracer_registry, only : add_tracer_OBC_values, tracer_registry_type +use MOM_variables, only : thermo_var_ptrs implicit none ; private #include -public Radiation_Open_Bdry_Conds, open_boundary_init, open_boundary_end - -type, public :: open_boundary_CS ; private - real :: gamma_uv ! The relative weighting for the baroclinic radiation - ! velocities (or speed of characteristics) at the - ! new time level (1) or the running mean (0) for velocities. - ! Valid values range from 0 to 1, with a default of 0.3. - real :: gamma_h ! The relative weighting for the baroclinic radiation - ! velocities (or speed of characteristics) at the - ! new time level (1) or the running mean (0) for thicknesses. - ! Valid values range from 0 to 1, with a default of 0.2. - 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. -end type open_boundary_CS +public open_boundary_config +public open_boundary_init +public open_boundary_query +public open_boundary_end +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_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 data +type, public :: ocean_OBC_type + 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 + !! local domain use west-facing Flather OBCs. + logical :: apply_OBC_v_flather_north = .false. !< True if any zonal velocity points in the + !! local domain use north-facing Flather OBCs. + logical :: apply_OBC_v_flather_south = .false. !< True if any zonal velocity points in the + !! local domain use south-facing Flather OBCs. + logical :: apply_OBC_u = .false. !< True if any zonal velocity points in to local domain use OBCs. + logical :: apply_OBC_v = .false. !< True if any meridional velocity points in to local domain use OBCs. + logical, pointer, dimension(:,:) :: & + OBC_mask_u => NULL(), & !< True at zonal velocity points that have prescribed OBCs. + OBC_mask_v => NULL() !< True at meridional velocity points that have prescribed OBCs. + ! These arrays indicate the kind of open boundary conditions that are to be applied at the u and v + ! 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. + ! 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 + ry_old_v => NULL(), & !< The ry_old_v value for radiation coeff for v-velocity in y-direction + rx_old_h => NULL(), & !< The rx_old_h value for radiation coeff for layer thickness h in x-direction + ry_old_h => NULL() !< The ry_old_h value for radiation coeff for layer thickness h in y-direction + + ! The following can be used to specify the outer-domain values of the + ! surface height and barotropic velocity. If these are not allocated, the + ! default with Flather boundary conditions is the same as if they were + ! filled with zeros. With simple OBCs, these should not be allocated. + real, pointer, dimension(:,:) :: & + ubt_outer => NULL(), & !< The u-velocity in the outer domain, in m s-1. + vbt_outer => NULL(), & !< The v-velocity in the outer domain, in m s-1. + eta_outer_u => NULL(), & !< The SSH anomaly in the outer domain, in m or kg m-2. + eta_outer_v => NULL() !< The SSH anomaly in the outer domain, in m or kg m-2. + + ! The following apply at points with OBC_kind_[uv] = OBC_SIMPLE. + real, pointer, dimension(:,:,:) :: & + u => NULL(), & !< The prescribed values of the zonal velocity (u) at OBC points. + v => NULL(), & !< The prescribed values of the meridional velocity (v) at OBC points. + uh => NULL(), & !< The prescribed values of the zonal volume transport (uh) at OBC points. + vh => NULL() !< The prescribed values of the meridional volume transport (vh) at OBC points. + + ! The following parameters are used in the baroclinic radiation code: + real :: gamma_uv !< The relative weighting for the baroclinic radiation + !! velocities (or speed of characteristics) at the + !! new time level (1) or the running mean (0) for velocities. + !! Valid values range from 0 to 1, with a default of 0.3. + real :: gamma_h !< The relative weighting for the baroclinic radiation + !! velocities (or speed of characteristics) at the + !! new time level (1) or the running mean (0) for thicknesses. + !! Valid values range from 0 to 1, with a default of 0.2. + 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? +end type ocean_OBC_type integer :: id_clock_pass +character(len=40) :: mod = "MOM_open_boundary" ! This module's name. +! This include declares and sets the variable "version". +#include "version_variable.h" + contains +!> Enables OBC module and reads configuration parameters +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 + + 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.) + + ! Safety check + if ((OBC%apply_OBC_u_flather_west .or. OBC%apply_OBC_v_flather_south) .and. & + .not.G%symmetric ) call MOM_error(FATAL, & + "MOM_open_boundary, open_boundary_config: "//& + "Symmetric memory must be used when APPLY_OBC_U_FLATHER_WEST "//& + "or APPLY_OBC_V_FLATHER_SOUTH is true.") + + if (.not.(OBC%apply_OBC_u .or. OBC%apply_OBC_v .or. & + OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south .or. & + OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west)) then + ! No open boundaries have been requested + call open_boundary_dealloc(OBC) + endif + +end subroutine open_boundary_config + +!> Initialize open boundary control structure +subroutine open_boundary_init(G, param_file, OBC) + type(ocean_grid_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 + + if (.not.associated(OBC)) return + + if ( OBC%apply_OBC_v_flather_north .or. OBC%apply_OBC_v_flather_south .or. & + OBC%apply_OBC_u_flather_east .or. OBC%apply_OBC_u_flather_west ) then + call get_param(param_file, mod, "OBC_RADIATION_MAX", OBC%rx_max, & + "The maximum magnitude of the baroclinic radiation \n"//& + "velocity (or speed of characteristics). This is only \n"//& + "used if one of the APPLY_OBC_[UV]_FLATHER_... is true.", & + units="m s-1", default=10.0) + call get_param(param_file, mod, "OBC_RAD_VEL_WT", OBC%gamma_uv, & + "The relative weighting for the baroclinic radiation \n"//& + "velocities (or speed of characteristics) at the new \n"//& + "time level (1) or the running mean (0) for velocities. \n"//& + "Valid values range from 0 to 1. This is only used if \n"//& + "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & + units="nondim", default=0.3) + call get_param(param_file, mod, "OBC_RAD_THICK_WT", OBC%gamma_h, & + "The relative weighting for the baroclinic radiation \n"//& + "velocities (or speed of characteristics) at the new \n"//& + "time level (1) or the running mean (0) for thicknesses. \n"//& + "Valid values range from 0 to 1. This is only used if \n"//& + "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & + units="nondim", default=0.2) + endif + + id_clock_pass = cpu_clock_id('(Ocean OBC halo updates)', grain=CLOCK_ROUTINE) + +end subroutine open_boundary_init + +!> Query the state of open boundary module configuration +logical function open_boundary_query(OBC, apply_orig_OBCs, apply_orig_Flather) + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + logical, optional, intent(in) :: apply_orig_OBCs !< If present, returns True if APPLY_OBC_U/V was set + logical, optional, intent(in) :: apply_orig_Flather !< If present, returns True if APPLY_OBC_*_FLATHER_* was set + open_boundary_query = .false. + if (.not. associated(OBC)) return + if (present(apply_orig_OBCs)) open_boundary_query = OBC%apply_OBC_u .or. OBC%apply_OBC_v + if (present(apply_orig_Flather)) open_boundary_query = OBC%apply_OBC_v_flather_north .or. & + OBC%apply_OBC_v_flather_south .or. & + OBC%apply_OBC_u_flather_east .or. & + OBC%apply_OBC_u_flather_west +end function open_boundary_query + +!> Deallocate open boundary data +subroutine open_boundary_dealloc(OBC) + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + 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%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) + if (associated(OBC%ry_old_h)) deallocate(OBC%ry_old_h) + if (associated(OBC%ubt_outer)) deallocate(OBC%ubt_outer) + if (associated(OBC%vbt_outer)) deallocate(OBC%vbt_outer) + if (associated(OBC%eta_outer_u)) deallocate(OBC%eta_outer_u) + if (associated(OBC%eta_outer_v)) deallocate(OBC%eta_outer_v) + if (associated(OBC%u)) deallocate(OBC%u) + if (associated(OBC%v)) deallocate(OBC%v) + if (associated(OBC%uh)) deallocate(OBC%uh) + if (associated(OBC%vh)) deallocate(OBC%vh) + deallocate(OBC) +end subroutine open_boundary_dealloc + +!> Close open boundary data +subroutine open_boundary_end(OBC) + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + call open_boundary_dealloc(OBC) +end subroutine open_boundary_end + +!> Sets the slope of bathymetry normal to an open bounndary to zero. +subroutine open_boundary_impose_normal_slope(OBC, G, depth) + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure + real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: depth !< Bathymetry at h-points + ! Local variables + integer :: i, j + + if (.not.associated(OBC)) return + + if (associated(OBC%OBC_direction_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) + enddo ; enddo + endif + + if (associated(OBC%OBC_kind_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) + enddo ; enddo + endif + +end subroutine open_boundary_impose_normal_slope + +!> Reconcile masks and open boundaries, deallocate OBC on PEs where it is not needed +subroutine open_boundary_impose_land_mask(OBC, G) + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure + type(dyn_horgrid_type), intent(in) :: G !< Ocean grid structure + ! Local variables + integer :: i, j + logical :: any_U, any_V + + 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. + 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. + endif + enddo ; enddo + endif + + any_U = .false. + if (associated(OBC%OBC_mask_u)) then + do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB + ! G%mask2du will be open wherever bathymetry allows it. + ! Bathymetry outside of the open boundary was adjusted to match + ! the bathymetry inside so these points will be open unless the + ! 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. + if (associated(OBC%OBC_mask_v)) then + 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. + +end subroutine open_boundary_impose_land_mask + +!> Diagnose radiation conditions at open boundaries subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & - h_new, h_old, G, CS) - type(ocean_grid_type), intent(inout) :: G - type(ocean_OBC_type), pointer :: OBC + h_new, h_old, G) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(ocean_OBC_type), pointer :: OBC !< Open boundary control structure real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u_new real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u_old real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v_new real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: 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 - type(open_boundary_CS), pointer :: CS - + ! Local variables real :: dhdt, dhdx, gamma_u, gamma_h, gamma_v real :: rx_max, ry_max ! coefficients for radiation real :: rx_new, rx_avg ! coefficients for radiation @@ -96,15 +351,13 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & if (.not.(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)) & return - if (.not.associated(CS)) call MOM_error(FATAL, & - "MOM_open_boundary: Module must be initialized before it is used.") - gamma_u = CS%gamma_uv ; gamma_v = CS%gamma_uv ; gamma_h = CS%gamma_h - rx_max = CS%rx_max ; ry_max = CS%rx_max + 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_kind_u(I,j) == OBC_FLATHER_E) then + if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) 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 @@ -121,7 +374,7 @@ 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_kind_u(I,j) == OBC_FLATHER_W) then + if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W) 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 @@ -143,7 +396,7 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & 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_kind_v(i,J) == OBC_FLATHER_N) then + if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) 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 rx_new = 0.0 @@ -161,7 +414,7 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & ! 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 - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) then + if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) 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 rx_new = 0.0 @@ -189,77 +442,350 @@ subroutine Radiation_Open_Bdry_Conds(OBC, u_new, u_old, v_new, v_old, & end subroutine Radiation_Open_Bdry_Conds -subroutine open_boundary_init(Time, G, param_file, diag, CS) - type(time_type), target, intent(in) :: Time - type(ocean_grid_type), intent(in) :: G - type(param_file_type), intent(in) :: param_file - type(diag_ctrl), target, intent(inout) :: diag - type(open_boundary_CS), pointer :: CS -! Arguments: Time - The current model time. -! (in) G - The ocean's grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) diag - A structure that is used to regulate diagnostic output. -! (in/out) CS - A pointer that is set to point to the control structure -! for this module -! This include declares and sets the variable "version". -#include "version_variable.h" - character(len=40) :: mod = "MOM_open_boundary" ! This module's name. - logical :: flather_east, flather_west, flather_north, flather_south +!> 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 (associated(CS)) then - call MOM_error(WARNING, "continuity_init called with associated control structure.") - return + 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 - call log_version(param_file, mod, version, "") - call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_EAST", flather_east, & - "If true, some zonal velocity points use Flather open \n"//& - "boundary conditions on the east side of the ocean.", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_U_FLATHER_WEST", flather_west, & - "If true, some zonal velocity points use Flather open \n"//& - "boundary conditions on the west side of the ocean.", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_NORTH", flather_north, & - "If true, some meridional velocity points use Flather \n"//& - "open boundary conditions on the north side of the ocean.", & - default=.false.) - call get_param(param_file, mod, "APPLY_OBC_V_FLATHER_SOUTH", flather_south, & - "If true, some meridional velocity points use Flather \n"//& - "open boundary conditions on the north side of the ocean.", & + ! 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) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure + type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure + real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h !< Thickness + type(param_file_type), intent(in) :: PF !< Parameter file handle + type(tracer_registry_type), pointer :: tracer_Reg !< Tracer registry + ! 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 + + call get_param(PF, mod, "READ_OBC_UV", read_OBC_uv, & + "If true, read the values for the velocity open boundary \n"//& + "conditions from the file specified by OBC_FILE.", & default=.false.) - if (.not.(flather_east .or. flather_west .or. flather_north .or. & - flather_south)) return - - allocate(CS) - call get_param(param_file, mod, "OBC_RADIATION_MAX", CS%rx_max, & - "The maximum magnitude of the baroclinic radiation \n"//& - "velocity (or speed of characteristics). This is only \n"//& - "used if one of the APPLY_OBC_[UV]_FLATHER_... is true.", & - units="m s-1", default=10.0) - call get_param(param_file, mod, "OBC_RAD_VEL_WT", CS%gamma_uv, & - "The relative weighting for the baroclinic radiation \n"//& - "velocities (or speed of characteristics) at the new \n"//& - "time level (1) or the running mean (0) for velocities. \n"//& - "Valid values range from 0 to 1. This is only used if \n"//& - "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & - units="nondim", default=0.3) - call get_param(param_file, mod, "OBC_RAD_THICK_WT", CS%gamma_h, & - "The relative weighting for the baroclinic radiation \n"//& - "velocities (or speed of characteristics) at the new \n"//& - "time level (1) or the running mean (0) for thicknesses. \n"//& - "Valid values range from 0 to 1. This is only used if \n"//& - "one of the APPLY_OBC_[UV]_FLATHER_... is true.", & - units="nondim", default=0.2) + call get_param(PF, mod, "READ_OBC_ETA", read_OBC_eta, & + "If true, read the values for the sea surface height \n"//& + "open boundary conditions from the file specified by \n"//& + "OBC_FILE.", default=.false.) + call get_param(PF, mod, "READ_OBC_TS", read_OBC_TS, & + "If true, read the values for the temperature and \n"//& + "salinity open boundary conditions from the file \n"//& + "specified by OBC_FILE.", default=.false.) + if (read_OBC_uv .or. read_OBC_eta .or. read_OBC_TS) then + call get_param(PF, mod, "OBC_FILE", OBC_file, & + "The file from which the appropriate open boundary \n"//& + "condition values are read.", default="MOM_OBC_FILE.nc") + call get_param(PF, mod, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + filename = trim(inputdir)//trim(OBC_file) + call log_param(PF, mod, "INPUTDIR/OBC_FILE", filename) + endif - id_clock_pass = cpu_clock_id('(Ocean OBC halo updates)', grain=CLOCK_ROUTINE) + if (.not.associated(OBC%vbt_outer)) then + allocate(OBC%vbt_outer(isd:ied,JsdB:JedB)) ; OBC%vbt_outer(:,:) = 0.0 + endif -end subroutine open_boundary_init + if (.not.associated(OBC%ubt_outer)) then + allocate(OBC%ubt_outer(IsdB:IedB,jsd:jed)) ; OBC%ubt_outer(:,:) = 0.0 + endif -subroutine open_boundary_end(CS) - type(open_boundary_CS), pointer :: CS - deallocate(CS) -end subroutine open_boundary_end + if (.not.associated(OBC%eta_outer_u)) then + allocate(OBC%eta_outer_u(IsdB:IedB,jsd:jed)) ; OBC%eta_outer_u(:,:) = 0.0 + endif + + if (.not.associated(OBC%eta_outer_v)) then + allocate(OBC%eta_outer_v(isd:ied,JsdB:JedB)) ; OBC%eta_outer_v(:,:) = 0.0 + endif + + if (read_OBC_uv) then + call read_data(filename, 'ubt', OBC%ubt_outer, & + domain=G%Domain%mpp_domain, position=EAST_FACE) + call read_data(filename, 'vbt', OBC%vbt_outer, & + domain=G%Domain%mpp_domain, position=NORTH_FACE) + endif + + if (read_OBC_eta) then + call read_data(filename, 'eta_outer_u', OBC%eta_outer_u, & + domain=G%Domain%mpp_domain, position=EAST_FACE) + call read_data(filename, 'eta_outer_v', OBC%eta_outer_v, & + domain=G%Domain%mpp_domain, position=NORTH_FACE) + endif + + call pass_vector(OBC%eta_outer_u,OBC%eta_outer_v,G%Domain, To_All+SCALAR_PAIR, CGRID_NE) + call pass_vector(OBC%ubt_outer,OBC%vbt_outer,G%Domain) + + ! Define radiation coefficients r[xy]_old_[uvh] as needed. For now, there are + ! no radiation conditions applied to the thicknesses, since the thicknesses + ! might not be physically motivated. Instead, sponges should be used to + ! enforce the near-boundary layer structure. + if (OBC%apply_OBC_u_flather_west .or. OBC%apply_OBC_u_flather_east) then + allocate(OBC%rx_old_u(IsdB:IedB,jsd:jed,nz)) ; OBC%rx_old_u(:,:,:) = 0.0 + ! allocate(OBC%rx_old_h(Isd:Ied,jsd:jed,nz)) ; OBC%rx_old_h(:,:,:) = 0.0 + endif + if (OBC%apply_OBC_v_flather_south .or. OBC%apply_OBC_v_flather_north) then + allocate(OBC%ry_old_v(isd:ied,JsdB:JedB,nz)) ; OBC%ry_old_v(:,:,:) = 0.0 + ! allocate(OBC%ry_old_h(isd:ied,Jsd:Jed,nz)) ; OBC%ry_old_h(:,:,:) = 0.0 + endif + + + if (associated(tv%T)) then + allocate(OBC_T_u(IsdB:IedB,jsd:jed,nz)) ; OBC_T_u(:,:,:) = 0.0 + allocate(OBC_S_u(IsdB:IedB,jsd:jed,nz)) ; OBC_S_u(:,:,:) = 0.0 + allocate(OBC_T_v(isd:ied,JsdB:JedB,nz)) ; OBC_T_v(:,:,:) = 0.0 + allocate(OBC_S_v(isd:ied,JsdB:JedB,nz)) ; OBC_S_v(:,:,:) = 0.0 + + if (read_OBC_TS) then + call read_data(filename, 'OBC_T_u', OBC_T_u, & + domain=G%Domain%mpp_domain, position=EAST_FACE) + call read_data(filename, 'OBC_S_u', OBC_S_u, & + domain=G%Domain%mpp_domain, position=EAST_FACE) + + call read_data(filename, 'OBC_T_v', OBC_T_v, & + domain=G%Domain%mpp_domain, position=NORTH_FACE) + call read_data(filename, 'OBC_S_v', OBC_S_v, & + domain=G%Domain%mpp_domain, position=NORTH_FACE) + else + call pass_var(tv%T, G%Domain) + call pass_var(tv%S, G%Domain) + 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 + OBC_T_u(I,j,k) = tv%T(i,j,k) + OBC_S_u(I,j,k) = tv%S(i,j,k) + elseif (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W) then + OBC_T_u(I,j,k) = tv%T(i+1,j,k) + OBC_S_u(I,j,k) = tv%S(i+1,j,k) + elseif (G%mask2dT(i,j) + G%mask2dT(i+1,j) > 0) then + OBC_T_u(I,j,k) = (G%mask2dT(i,j)*tv%T(i,j,k) + G%mask2dT(i+1,j)*tv%T(i+1,j,k)) / & + (G%mask2dT(i,j) + G%mask2dT(i+1,j)) + OBC_S_u(I,j,k) = (G%mask2dT(i,j)*tv%S(i,j,k) + G%mask2dT(i+1,j)*tv%S(i+1,j,k)) / & + (G%mask2dT(i,j) + G%mask2dT(i+1,j)) + else ! This probably shouldn't happen or maybe it doesn't matter? + OBC_T_u(I,j,k) = 0.5*(tv%T(i,j,k)+tv%T(i+1,j,k)) + OBC_S_u(I,j,k) = 0.5*(tv%S(i,j,k)+tv%S(i+1,j,k)) + endif + else + OBC_T_u(I,j,k) = 0.5*(tv%T(i,j,k)+tv%T(i+1,j,k)) + OBC_S_u(I,j,k) = 0.5*(tv%S(i,j,k)+tv%S(i+1,j,k)) + 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_direction_v(i,J) == OBC_DIRECTION_N) then + OBC_T_v(i,J,k) = tv%T(i,j,k) + OBC_S_v(i,J,k) = tv%S(i,j,k) + elseif (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) then + OBC_T_v(i,J,k) = tv%T(i,j+1,k) + OBC_S_v(i,J,k) = tv%S(i,j+1,k) + elseif (G%mask2dT(i,j) + G%mask2dT(i,j+1) > 0) then + OBC_T_v(i,J,k) = (G%mask2dT(i,j)*tv%T(i,j,k) + G%mask2dT(i,j+1)*tv%T(i,j+1,k)) / & + (G%mask2dT(i,j) + G%mask2dT(i,j+1)) + OBC_S_v(i,J,k) = (G%mask2dT(i,j)*tv%S(i,j,k) + G%mask2dT(i,j+1)*tv%S(i,j+1,k)) / & + (G%mask2dT(i,j) + G%mask2dT(i,j+1)) + else ! This probably shouldn't happen or maybe it doesn't matter? + OBC_T_v(i,J,k) = 0.5*(tv%T(i,j,k)+tv%T(i,j+1,k)) + OBC_S_v(i,J,k) = 0.5*(tv%S(i,j,k)+tv%S(i,j+1,k)) + endif + else + OBC_T_v(i,J,k) = 0.5*(tv%T(i,j,k)+tv%T(i,j+1,k)) + OBC_S_v(i,J,k) = 0.5*(tv%S(i,j,k)+tv%S(i,j+1,k)) + endif + enddo; enddo ; enddo + endif + + call pass_vector(OBC_T_u, OBC_T_v, G%Domain, To_All+SCALAR_PAIR, CGRID_NE) + call pass_vector(OBC_S_u, OBC_S_v, G%Domain, To_All+SCALAR_PAIR, CGRID_NE) + + call add_tracer_OBC_values("T", tracer_Reg, OBC_in_u=OBC_T_u, & + OBC_in_v=OBC_T_v) + call add_tracer_OBC_values("S", tracer_Reg, OBC_in_u=OBC_S_u, & + OBC_in_v=OBC_S_v) + do k=1,nz ; do j=jsd,jed ; do I=isd,ied-1 + if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) then + tv%T(i+1,j,k) = tv%T(i,j,k) ; tv%S(i+1,j,k) = tv%S(i,j,k) + elseif (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W) then + tv%T(i,j,k) = tv%T(i+1,j,k) ; tv%S(i,j,k) = tv%S(i+1,j,k) + endif + enddo ; enddo ; enddo + do k=1,nz ; do J=jsd,jed-1 ; do i=isd,ied + if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) then + tv%T(i,j+1,k) = tv%T(i,j,k) ; tv%S(i,j+1,k) = tv%S(i,j,k) + elseif (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) then + tv%T(i,j,k) = tv%T(i,j+1,k) ; tv%S(i,j,k) = tv%S(i,j+1,k) + endif + enddo ; enddo ; enddo + endif + + do k=1,nz ; do j=jsd,jed ; do I=isd,ied-1 + if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E) h(i+1,j,k) = h(i,j,k) + if (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W) h(i,j,k) = h(i+1,j,k) + enddo ; enddo ; enddo + do k=1,nz ; do J=jsd,jed-1 ; do i=isd,ied + if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N) h(i,j+1,k) = h(i,j,k) + if (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S) h(i,j,k) = h(i,j+1,k) + enddo ; enddo ; enddo + +end subroutine set_Flather_data + +!> \namespace mom_open_boundary +!! This module implements some aspects of internal open boundary +!! conditions in MOM. +!! +!! 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_open_boundary diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 42a619c9b9..fa1b2da980 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -235,59 +235,6 @@ module MOM_variables ! at the interfaces between each layer, in m2 s-2. end type vertvisc_type -integer, parameter, public :: OBC_NONE = 0, OBC_SIMPLE = 1, OBC_WALL = 2 -integer, parameter, public :: OBC_FLATHER_E = 4, OBC_FLATHER_W = 5 -integer, parameter, public :: OBC_FLATHER_N = 6, OBC_FLATHER_S = 7 - -type, public :: ocean_OBC_type -! This structure is used to apply specified open boundary conditions. - logical :: apply_OBC_u_flather_east = .false. ! If true, some zonal velocity - logical :: apply_OBC_u_flather_west = .false. ! points in the local domain use flather open - ! boundary conditions. - logical :: apply_OBC_v_flather_north = .false. ! If true, some meridional velocity - logical :: apply_OBC_v_flather_south = .false. ! points in the local domain use flather open - ! boundary conditions. - logical :: apply_OBC_u = .false. ! If true, some zonal or meridional velocity - logical :: apply_OBC_v = .false. ! points in the local domain use open - ! boundary conditions. - logical, pointer, dimension(:,:) :: & - OBC_mask_u => NULL(), & ! These arrays are true at zonal or meridional - OBC_mask_v => NULL() ! velocity points that have prescribed open boundary - ! conditions. - integer, pointer, dimension(:,:) :: & - OBC_kind_u => NULL(), & ! These arrays indicate the kind of open boundary - OBC_kind_v => NULL() ! conditions that are to be applied at the u and v - ! points, and can be OBC_NONE, OBC_SIMPLE, OBC_WALL, - ! or one of OBC_FLATHER_[EWNS]. 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. - ! The following apply at points with OBC_kind_[uv] = OBC_FLATHER_x. - real, pointer, dimension(:,:,:) :: & - rx_old_u => NULL(), & ! The rx_old_u value for radiation coeff for u-velocity in x-direction - ry_old_v => NULL(), & ! The ry_old_v value for radiation coeff for v-velocity in y-direction - rx_old_h => NULL(), & ! The rx_old_h value for radiation coeff for layer thickness h in x-direction - ry_old_h => NULL() ! The ry_old_h value for radiation coeff for layer thickness h in y-direction - - ! The following can be used to specify the outer-domain values of the - ! surface height and barotropic velocity. If these are not allocated, the - ! default with Flather boundary conditions is the same as if they were - ! filled with zeros. With simple OBCs, these should not be allocated. - real, pointer, dimension(:,:) :: & - ubt_outer => NULL(), & ! The u-velocity in the outer domain, in m s-1. - vbt_outer => NULL(), & ! The v-velocity in the outer domain, in m s-1. - eta_outer_u => NULL(), & ! The sea surface height anomaly or water column - eta_outer_v => NULL() ! mass anomaly in the outer domain in m or kg m-2. - - ! The following apply at points with OBC_kind_[uv] = OBC_SIMPLE. - real, pointer, dimension(:,:,:) :: & - u => NULL(), & ! The prescribed values of the zonal (u) or meridional (v) - v => NULL(), & ! velocities at OBC points, in m s-1. - uh => NULL(), & ! The prescribed values of the zonal (uh) or meridional (vh) - vh => NULL() ! volume transports at OBC points, in m3 s-1. -end type ocean_OBC_type - - type, public :: BT_cont_type real, pointer, dimension(:,:) :: & FA_u_EE => NULL(), & ! The FA_u_XX variables are the effective open face diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90 index 252c4778c1..50635b5d6a 100644 --- a/src/initialization/MOM_fixed_initialization.F90 +++ b/src/initialization/MOM_fixed_initialization.F90 @@ -18,9 +18,13 @@ module MOM_fixed_initialization use MOM_io, only : slasher, vardesc, write_field, var_desc use MOM_io, only : EAST_FACE, NORTH_FACE use MOM_grid_initialize, only : initialize_masks, set_grid_metrics +use MOM_open_boundary, only : ocean_OBC_type +use MOM_open_boundary, only : open_boundary_config, open_boundary_query +use MOM_open_boundary, only : set_Flather_positions, open_boundary_impose_normal_slope +use MOM_open_boundary, only : open_boundary_impose_land_mask use MOM_string_functions, only : uppercase -use user_initialization, only : user_initialize_topography -use DOME_initialization, only : DOME_initialize_topography +use user_initialization, only : user_initialize_topography, USER_set_OBC_positions +use DOME_initialization, only : DOME_initialize_topography, DOME_set_OBC_positions use ISOMIP_initialization, only : ISOMIP_initialize_topography use benchmark_initialization, only : benchmark_initialize_topography use DOME2d_initialization, only : DOME2d_initialize_topography @@ -41,8 +45,9 @@ module MOM_fixed_initialization ! ----------------------------------------------------------------------------- !> MOM_initialize_fixed sets up time-invariant quantities related to MOM6's !! horizontal grid, bathymetry, and the Coriolis parameter. -subroutine MOM_initialize_fixed(G, PF, write_geom, output_dir) +subroutine MOM_initialize_fixed(G, OBC, PF, write_geom, output_dir) type(dyn_horgrid_type), intent(inout) :: G !< The ocean's grid structure. + type(ocean_OBC_type), pointer :: OBC !< Open boundary structure. type(param_file_type), intent(in) :: PF !< A structure indicating the open file !! to parse for model parameter values. logical, intent(in) :: write_geom !< If true, write grid geometry files. @@ -76,9 +81,40 @@ subroutine MOM_initialize_fixed(G, PF, write_geom, output_dir) ! masks, and Coriolis parameter. ! ==================================================================== -! This call sets seamasks that prohibit flow over any point with ! -! a bottom that is shallower than min_depth from PF. ! +! Determine the position of any open boundaries + call open_boundary_config(G, PF, OBC) + if (open_boundary_query(OBC, apply_orig_OBCs=.true.)) then + call get_param(PF, mod, "OBC_CONFIG", config, & + "A string that sets how the open boundary conditions are \n"//& + " configured: \n"//& + " \t DOME - use a slope and channel configuration for the \n"//& + " \t\t DOME sill-overflow test case. \n"//& + " \t USER - call a user modified routine.", default="file", & + fail_if_missing=.true.) + select case ( trim(config) ) + case ("none") + case ("DOME") ; call DOME_set_OBC_positions(G, PF, OBC) + case ("USER") ; call user_set_OBC_positions(G, PF, OBC) + case default ; call MOM_error(FATAL, "MOM_initialize_fixed: "// & + "The open boundary positions specified by OBC_CONFIG="//& + trim(config)//" have not been fully implemented.") + end select + elseif (open_boundary_query(OBC, apply_orig_Flather=.true.)) then + call set_Flather_positions(G, OBC) + endif + + ! To initialize masks, the bathymetry in halo regions must be filled in + call pass_var(G%bathyT, G%Domain) + + ! Make bathymetry consistent with open boundaries + call open_boundary_impose_normal_slope(OBC, G, G%bathyT) + + ! This call sets masks that prohibit flow over any point interpreted as land call initialize_masks(G, PF) + + ! Make OBC mask consistent with land mask, deallocate OBC on PEs where it is not needed + call open_boundary_impose_land_mask(OBC, G) + if (debug) then call hchksum(G%bathyT, 'MOM_initialize_fixed: depth ', G%HI, haloshift=1) call hchksum(G%mask2dT, 'MOM_initialize_fixed: mask2dT ', G%HI) diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 61a6f5ed66..c9505d3eb7 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -1347,8 +1347,6 @@ subroutine initialize_masks(G, PF) real :: Dmin, min_depth, mask_depth integer :: i, j - logical :: apply_OBC_u_flather_east, apply_OBC_u_flather_west - logical :: apply_OBC_v_flather_north, apply_OBC_v_flather_south character(len=40) :: mod = "MOM_grid_init initialize_masks" call callTree_enter("initialize_masks(), MOM_grid_initialize.F90") @@ -1362,65 +1360,13 @@ subroutine initialize_masks(G, PF) "The depth below which to mask points as land points, for which all\n"//& "fluxes are zeroed out. MASKING_DEPTH is ignored if negative.", & units="m", default=-9999.0) - call get_param(PF, mod, "APPLY_OBC_U_FLATHER_EAST", 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(PF, mod, "APPLY_OBC_U_FLATHER_WEST", 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(PF, mod, "APPLY_OBC_V_FLATHER_NORTH", 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(PF, mod, "APPLY_OBC_V_FLATHER_SOUTH", apply_OBC_v_flather_south,& - "Apply a Flather open boundary condition on the southern \n"//& - "side of the global domain", default=.false.) - - if ((apply_OBC_u_flather_west .or. apply_OBC_v_flather_south) .and. & - .not.G%symmetric ) & - call MOM_error(FATAL, "Symmetric memory must be used when "//& - "APPLY_OBC_U_FLATHER_WEST or APPLY_OBC_V_FLATHER_SOUTH is true.") Dmin = min_depth if (mask_depth>=0.) Dmin = mask_depth - call pass_var(G%bathyT, G%Domain) G%mask2dCu(:,:) = 0.0 ; G%mask2dCv(:,:) = 0.0 ; G%mask2dBu(:,:) = 0.0 - ! Extrapolate the bottom depths at any points that are subject to Flather - ! open boundary conditions. This should be generalized for Flather OBCs - ! that are not necessarily at the edges of the domain. - if (apply_OBC_u_flather_west) then - do j=G%jsd,G%jed ; do I=G%isd+1,G%ied - if ((I+G%idg_offset) == G%isg) then - G%bathyT(i-1,j) = G%bathyT(i,j) - endif - enddo; enddo - endif - - if (apply_OBC_u_flather_east) then - do j=G%jsd,G%jed ; do I=G%isd,G%ied-1 - if ((i+G%idg_offset) == G%ieg) then - G%bathyT(i+1,j) = G%bathyT(i,j) - endif - enddo; enddo - endif - - if (apply_OBC_v_flather_north) then - do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied - if ((j+G%jdg_offset) == G%jeg) then - G%bathyT(i,j+1) = G%bathyT(i,j) - endif - enddo; enddo - endif - - if (apply_OBC_v_flather_south) then - do J=G%jsd+1,G%jed ; do i=G%isd,G%ied - if ((J+G%jdg_offset) == G%jsg) then - G%bathyT(i,j-1) = G%bathyT(i,j) - endif - enddo; enddo - endif - + ! Construct the h-point or T-point mask do j=G%jsd,G%jed ; do i=G%isd,G%ied if (G%bathyT(i,j) <= Dmin) then G%mask2dT(i,j) = 0.0 diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 5f37446db6..7dbe4d2843 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -21,6 +21,9 @@ module MOM_state_initialization use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE, MULTIPLE use MOM_io, only : slasher, vardesc, write_field use MOM_io, only : EAST_FACE, NORTH_FACE +use MOM_open_boundary, only : ocean_OBC_type, open_boundary_init +use MOM_open_boundary, only : OBC_NONE, OBC_SIMPLE +use MOM_open_boundary, only : open_boundary_query, set_Flather_data, set_Flather_positions use MOM_grid_initialize, only : initialize_masks, set_grid_metrics use MOM_restart, only : restore_state, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, set_up_sponge_ML_density @@ -30,18 +33,17 @@ module MOM_state_initialization use MOM_string_functions, only : uppercase use MOM_time_manager, only : time_type, set_time use MOM_tracer_registry, only : add_tracer_OBC_values, tracer_registry_type -use MOM_variables, only : thermo_var_ptrs, ocean_OBC_type -use MOM_variables, only : OBC_NONE, OBC_SIMPLE, OBC_FLATHER_E, OBC_FLATHER_W -use MOM_variables, only : OBC_FLATHER_N, OBC_FLATHER_S +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : setVerticalGridAxes, verticalGrid_type use MOM_ALE, only : pressure_gradient_plm use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type use MOM_EOS, only : int_specific_vol_dp use user_initialization, only : user_initialize_thickness, user_initialize_velocity use user_initialization, only : user_init_temperature_salinity -use user_initialization, only : user_set_Open_Bdry_Conds, user_initialize_sponges +use user_initialization, only : user_set_OBC_positions, user_set_OBC_data +use user_initialization, only : user_initialize_sponges use DOME_initialization, only : DOME_initialize_thickness -use DOME_initialization, only : DOME_set_Open_Bdry_Conds +use DOME_initialization, only : DOME_set_OBC_positions, DOME_set_OBC_data use DOME_initialization, only : DOME_initialize_sponges use ISOMIP_initialization, only : ISOMIP_initialize_thickness use ISOMIP_initialization, only : ISOMIP_initialize_sponges @@ -141,9 +143,6 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & logical :: trim_ic_for_p_surf ! If true, remove the mass that would be displaced ! by a large surface pressure, such as with an ice sheet. logical :: Analytic_FV_PGF, obsol_test - logical :: apply_OBC_u, apply_OBC_v - logical :: apply_OBC_u_flather_east, apply_OBC_u_flather_west - logical :: apply_OBC_v_flather_north, apply_OBC_v_flather_south logical :: convert type(EOS_type), pointer :: eos => NULL() logical :: debug ! indicates whether to write debugging output @@ -422,48 +421,29 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, PF, dirs, & end select endif -! This subroutine call sets optional open boundary conditions. - call get_param(PF, mod, "APPLY_OBC_U", 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(PF, mod, "APPLY_OBC_V", 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.) - if (apply_OBC_u .or. apply_OBC_v) then - call get_param(PF, mod, "OBC_CONFIG", config, & - "A string that sets how the open boundary conditions are \n"//& - " configured: \n"//& - " \t DOME - use a slope and channel configuration for the \n"//& - " \t\t DOME sill-overflow test case. \n"//& - " \t USER - call a user modified routine.", default="file", & - fail_if_missing=.true.) + ! Reads OBC parameters not pertaining to the location of the boundaries + call open_boundary_init(G, PF, OBC) + + ! This is the legacy approach to turning on open boundaries + if (open_boundary_query(OBC, apply_orig_OBCs=.true.)) then + call get_param(PF, mod, "OBC_CONFIG", config, fail_if_missing=.true., do_not_log=.true.) if (trim(config) == "DOME") then - call DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, PF, tracer_Reg) + call DOME_set_OBC_data(OBC, tv, G, GV, PF, tracer_Reg) elseif (trim(config) == "USER") then - call user_set_Open_Bdry_Conds(OBC, tv, G, PF, tracer_Reg) + call user_set_OBC_data(OBC, tv, G, PF, tracer_Reg) else call MOM_error(FATAL, "The open boundary conditions specified by "//& "OBC_CONFIG = "//trim(config)//" have not been fully implemented.") call set_Open_Bdry_Conds(OBC, tv, G, GV, PF, tracer_Reg) endif + elseif (open_boundary_query(OBC, apply_orig_Flather=.true.)) then + call set_Flather_data(OBC, tv, h, G, PF, tracer_Reg) endif - - call get_param(PF, mod, "APPLY_OBC_U_FLATHER_EAST", 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(PF, mod, "APPLY_OBC_U_FLATHER_WEST", 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(PF, mod, "APPLY_OBC_V_FLATHER_NORTH", 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(PF, mod, "APPLY_OBC_V_FLATHER_SOUTH", apply_OBC_v_flather_south,& - "Apply a Flather open boundary condition on the southern \n"//& - "side of the global domain", default=.false.) - if (apply_OBC_u_flather_east .or. apply_OBC_u_flather_west .or. apply_OBC_v_flather_north .or. apply_OBC_v_flather_south) then - call set_Flather_Bdry_Conds(OBC, tv, h, G, PF, tracer_Reg) + if (debug.and.associated(OBC)) then + call hchksum(G%mask2dT, 'MOM_initialize_state: mask2dT ', G%HI) + call uchksum(G%mask2dCu, 'MOM_initialize_state: mask2dCu ', G%HI) + call vchksum(G%mask2dCv, 'MOM_initialize_state: mask2dCv ', G%HI) + call qchksum(G%mask2dBu, 'MOM_initialize_state: mask2dBu ', G%HI) endif call callTree_leave('MOM_initialize_state()') @@ -1695,361 +1675,6 @@ subroutine set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tracer_Reg) end subroutine set_Open_Bdry_Conds ! ----------------------------------------------------------------------------- -! ----------------------------------------------------------------------------- -subroutine set_Flather_Bdry_Conds(OBC, tv, h, G, PF, tracer_Reg) - type(ocean_grid_type), intent(inout) :: G - type(ocean_OBC_type), pointer :: OBC - type(thermo_var_ptrs), intent(inout) :: tv - real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(inout) :: h - type(param_file_type), intent(in) :: PF - type(tracer_registry_type), pointer :: tracer_Reg -! This subroutine sets the initial definitions of the characteristic open boundary -! conditions. Written by Mehmet Ilicak - -! Arguments: OBC - This open boundary condition type specifies whether, where, -! and what open boundary conditions are used. -! (out) tv - A structure containing pointers to any available -! thermodynamic fields, including potential temperature and -! salinity or mixed layer density. Absent fields have NULL ptrs. -! (in) G - The ocean's grid structure. -! (in) PF - A structure indicating the open file to parse for -! model parameter values. - - logical :: any_OBC ! Set to true if any points in this subdomain use - ! open boundary conditions. - - logical :: apply_OBC_u_flather_east = .false., apply_OBC_u_flather_west = .false. - logical :: apply_OBC_v_flather_north = .false., apply_OBC_v_flather_south = .false. - 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 - integer :: east_boundary, west_boundary, north_boundary, south_boundary - 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 - - call get_param(PF, mod, "APPLY_OBC_U_FLATHER_EAST", 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(PF, mod, "APPLY_OBC_U_FLATHER_WEST", 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(PF, mod, "APPLY_OBC_V_FLATHER_NORTH", 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(PF, mod, "APPLY_OBC_V_FLATHER_SOUTH", apply_OBC_v_flather_south,& - "Apply a Flather open boundary condition on the southern \n"//& - "side of the global domain", default=.false.) - - if (.not.(apply_OBC_u_flather_east .or. apply_OBC_u_flather_west .or. & - apply_OBC_v_flather_north .or. apply_OBC_v_flather_south)) return - - if (.not.associated(OBC)) allocate(OBC) - - OBC%apply_OBC_u_flather_east = apply_OBC_u_flather_east - OBC%apply_OBC_u_flather_west = apply_OBC_u_flather_west - OBC%apply_OBC_v_flather_north = apply_OBC_v_flather_north - OBC%apply_OBC_v_flather_south = apply_OBC_v_flather_south - - call get_param(PF, mod, "READ_OBC_UV", read_OBC_uv, & - "If true, read the values for the velocity open boundary \n"//& - "conditions from the file specified by OBC_FILE.", & - default=.false.) - call get_param(PF, mod, "READ_OBC_ETA", read_OBC_eta, & - "If true, read the values for the sea surface height \n"//& - "open boundary conditions from the file specified by \n"//& - "OBC_FILE.", default=.false.) - call get_param(PF, mod, "READ_OBC_TS", read_OBC_TS, & - "If true, read the values for the temperature and \n"//& - "salinity open boundary conditions from the file \n"//& - "specified by OBC_FILE.", default=.false.) - if (read_OBC_uv .or. read_OBC_eta .or. read_OBC_TS) then - call get_param(PF, mod, "OBC_FILE", OBC_file, & - "The file from which the appropriate open boundary \n"//& - "condition values are read.", default="MOM_OBC_FILE.nc") - call get_param(PF, mod, "INPUTDIR", inputdir, default=".") - inputdir = slasher(inputdir) - filename = trim(inputdir)//trim(OBC_file) - call log_param(PF, mod, "INPUTDIR/OBC_FILE", filename) - endif - - 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 (.not.associated(OBC%OBC_mask_u)) then - allocate(OBC%OBC_mask_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_mask_u(:,:) = .false. - endif - if (.not.associated(OBC%OBC_kind_u)) then - allocate(OBC%OBC_kind_u(IsdB:IedB,jsd:jed)) ; OBC%OBC_kind_u(:,:) = OBC_NONE - endif - if (.not.associated(OBC%OBC_mask_v)) then - allocate(OBC%OBC_mask_v(isd:ied,JsdB:JedB)) ; OBC%OBC_mask_v(:,:) = .false. - endif - if (.not.associated(OBC%OBC_kind_v)) then - allocate(OBC%OBC_kind_v(isd:ied,JsdB:JedB)) ; OBC%OBC_kind_v(:,:) = OBC_NONE - endif - - if (.not.associated(OBC%vbt_outer)) then - allocate(OBC%vbt_outer(isd:ied,JsdB:JedB)) ; OBC%vbt_outer(:,:) = 0.0 - endif - - if (.not.associated(OBC%ubt_outer)) then - allocate(OBC%ubt_outer(IsdB:IedB,jsd:jed)) ; OBC%ubt_outer(:,:) = 0.0 - endif - - if (.not.associated(OBC%eta_outer_u)) then - allocate(OBC%eta_outer_u(IsdB:IedB,jsd:jed)) ; OBC%eta_outer_u(:,:) = 0.0 - endif - - if (.not.associated(OBC%eta_outer_v)) then - allocate(OBC%eta_outer_v(isd:ied,JsdB:JedB)) ; OBC%eta_outer_v(:,:) = 0.0 - endif - - if (read_OBC_uv) then - call read_data(filename, 'ubt', OBC%ubt_outer, & - domain=G%Domain%mpp_domain, position=EAST_FACE) - call read_data(filename, 'vbt', OBC%vbt_outer, & - domain=G%Domain%mpp_domain, position=NORTH_FACE) - endif - - if (read_OBC_eta) then - call read_data(filename, 'eta_outer_u', OBC%eta_outer_u, & - domain=G%Domain%mpp_domain, position=EAST_FACE) - call read_data(filename, 'eta_outer_v', OBC%eta_outer_v, & - domain=G%Domain%mpp_domain, position=NORTH_FACE) - endif - - call pass_vector(OBC%eta_outer_u,OBC%eta_outer_v,G%Domain, To_All+SCALAR_PAIR, CGRID_NE) - call pass_vector(OBC%ubt_outer,OBC%vbt_outer,G%Domain) - - ! This code should be modified to allow OBCs to be applied anywhere. - - if (apply_OBC_u_flather_east) then - ! Determine where u points are applied at east side - do j=jsd,jed ; do I=IsdB,IedB - if ((I+G%idg_offset) == east_boundary) then !eastern side - if (G%mask2dCu(I,j) > 0.50) then - OBC%OBC_mask_u(I,j) = .true. - OBC%OBC_kind_u(I,j) = OBC_FLATHER_E - if (G%mask2dCv(i+1,J) > 0.50) then - OBC%OBC_mask_v(i+1,J) = .true. - if (OBC%OBC_kind_v(i+1,J) == OBC_NONE) OBC%OBC_kind_v(i+1,J) = OBC_FLATHER_E - endif - if (G%mask2dCv(i+1,J-1) > 0.50) then - OBC%OBC_mask_v(i+1,J-1) = .true. - if (OBC%OBC_kind_v(i+1,J-1) == OBC_NONE) OBC%OBC_kind_v(i+1,J-1) = OBC_FLATHER_E - endif - endif - endif - enddo ; enddo - endif - - if (apply_OBC_u_flather_west) then - ! Determine where u points are applied at west side - do j=jsd,jed ; do I=IsdB,IedB - if ((I+G%idg_offset) == west_boundary) then !western side - if (G%mask2dCu(I,j) > 0.50) then - OBC%OBC_mask_u(I,j) = .true. - OBC%OBC_kind_u(I,j) = OBC_FLATHER_W - if (G%mask2dCv(i,J) > 0.50) then - OBC%OBC_mask_v(i,J) = .true. - if (OBC%OBC_kind_v(i,J) == OBC_NONE) OBC%OBC_kind_v(i,J) = OBC_FLATHER_W - endif - if (G%mask2dCv(i,J-1) > 0.50) then - OBC%OBC_mask_v(i,J-1) = .true. - if (OBC%OBC_kind_v(i,J-1) == OBC_NONE) OBC%OBC_kind_v(i,J-1) = OBC_FLATHER_W - endif - endif - endif - enddo ; enddo - endif - - - if (apply_OBC_v_flather_north) then - ! Determine where v points are applied at north side - do J=JsdB,JedB ; do i=isd,ied - if ((J+G%jdg_offset) == north_boundary) then !northern side - if (G%mask2dCv(i,J) > 0.50) then - OBC%OBC_mask_v(i,J) = .true. - OBC%OBC_kind_v(i,J) = OBC_FLATHER_N - if (G%mask2dCu(I,j+1) > 0.50) then - OBC%OBC_mask_u(I,j+1) = .true. - if (OBC%OBC_kind_u(I,j+1) == OBC_NONE) OBC%OBC_kind_u(I,j+1) = OBC_FLATHER_N - endif - if (G%mask2dCu(I-1,j+1) > 0.50) then - OBC%OBC_mask_u(I-1,j+1) = .true. - if (OBC%OBC_kind_u(I-1,j+1) == OBC_NONE) OBC%OBC_kind_u(I-1,j+1) = OBC_FLATHER_N - endif - endif - endif - enddo ; enddo - endif - - if (apply_OBC_v_flather_south) then - ! Determine where v points are applied at south side - do J=JsdB,JedB ; do i=isd,ied - if ((J+G%jdg_offset) == south_boundary) then !southern side - if (G%mask2dCv(i,J) > 0.50) then - OBC%OBC_mask_v(i,J) = .true. - OBC%OBC_kind_v(i,J) = OBC_FLATHER_S - if (G%mask2dCu(I,j) > 0.50) then - OBC%OBC_mask_u(I,j) = .true. - if (OBC%OBC_kind_u(I,j) == OBC_NONE) OBC%OBC_kind_u(I,j) = OBC_FLATHER_S - endif - if (G%mask2dCu(I-1,j) > 0.50) then - OBC%OBC_mask_u(I-1,j) = .true. - if (OBC%OBC_kind_u(I-1,j) == OBC_NONE) OBC%OBC_kind_u(I-1,j) = OBC_FLATHER_S - endif - endif - endif - enddo ; enddo - endif - - ! If there are no OBC points on this PE, there is no reason to keep the OBC - ! type, and it could be deallocated. - - - ! Define radiation coefficients r[xy]_old_[uvh] as needed. For now, there are - ! no radiation conditions applied to the thicknesses, since the thicknesses - ! might not be physically motivated. Instead, sponges should be used to - ! enforce the near-boundary layer structure. - if (apply_OBC_u_flather_west .or. apply_OBC_u_flather_east) then - allocate(OBC%rx_old_u(IsdB:IedB,jsd:jed,nz)) ; OBC%rx_old_u(:,:,:) = 0.0 - ! allocate(OBC%rx_old_h(Isd:Ied,jsd:jed,nz)) ; OBC%rx_old_h(:,:,:) = 0.0 - endif - if (apply_OBC_v_flather_south .or. apply_OBC_v_flather_north) then - allocate(OBC%ry_old_v(isd:ied,JsdB:JedB,nz)) ; OBC%ry_old_v(:,:,:) = 0.0 - ! allocate(OBC%ry_old_h(isd:ied,Jsd:Jed,nz)) ; OBC%ry_old_h(:,:,:) = 0.0 - endif - - - if (associated(tv%T)) then - allocate(OBC_T_u(IsdB:IedB,jsd:jed,nz)) ; OBC_T_u(:,:,:) = 0.0 - allocate(OBC_S_u(IsdB:IedB,jsd:jed,nz)) ; OBC_S_u(:,:,:) = 0.0 - allocate(OBC_T_v(isd:ied,JsdB:JedB,nz)) ; OBC_T_v(:,:,:) = 0.0 - allocate(OBC_S_v(isd:ied,JsdB:JedB,nz)) ; OBC_S_v(:,:,:) = 0.0 - - if (read_OBC_TS) then - call read_data(filename, 'OBC_T_u', OBC_T_u, & - domain=G%Domain%mpp_domain, position=EAST_FACE) - call read_data(filename, 'OBC_S_u', OBC_S_u, & - domain=G%Domain%mpp_domain, position=EAST_FACE) - - call read_data(filename, 'OBC_T_v', OBC_T_v, & - domain=G%Domain%mpp_domain, position=NORTH_FACE) - call read_data(filename, 'OBC_S_v', OBC_S_v, & - domain=G%Domain%mpp_domain, position=NORTH_FACE) - else - call pass_var(tv%T, G%Domain) - call pass_var(tv%S, G%Domain) - do k=1,nz ; 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_FLATHER_E) then - OBC_T_u(I,j,k) = tv%T(i,j,k) - OBC_S_u(I,j,k) = tv%S(i,j,k) - elseif (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) then - OBC_T_u(I,j,k) = tv%T(i+1,j,k) - OBC_S_u(I,j,k) = tv%S(i+1,j,k) - elseif (G%mask2dT(i,j) + G%mask2dT(i+1,j) > 0) then - OBC_T_u(I,j,k) = (G%mask2dT(i,j)*tv%T(i,j,k) + G%mask2dT(i+1,j)*tv%T(i+1,j,k)) / & - (G%mask2dT(i,j) + G%mask2dT(i+1,j)) - OBC_S_u(I,j,k) = (G%mask2dT(i,j)*tv%S(i,j,k) + G%mask2dT(i+1,j)*tv%S(i+1,j,k)) / & - (G%mask2dT(i,j) + G%mask2dT(i+1,j)) - else ! This probably shouldn't happen or maybe it doesn't matter? - OBC_T_u(I,j,k) = 0.5*(tv%T(i,j,k)+tv%T(i+1,j,k)) - OBC_S_u(I,j,k) = 0.5*(tv%S(i,j,k)+tv%S(i+1,j,k)) - endif - else - OBC_T_u(I,j,k) = 0.5*(tv%T(i,j,k)+tv%T(i+1,j,k)) - OBC_S_u(I,j,k) = 0.5*(tv%S(i,j,k)+tv%S(i+1,j,k)) - 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_kind_v(i,J) == OBC_FLATHER_N) then - OBC_T_v(i,J,k) = tv%T(i,j,k) - OBC_S_v(i,J,k) = tv%S(i,j,k) - elseif (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) then - OBC_T_v(i,J,k) = tv%T(i,j+1,k) - OBC_S_v(i,J,k) = tv%S(i,j+1,k) - elseif (G%mask2dT(i,j) + G%mask2dT(i,j+1) > 0) then - OBC_T_v(i,J,k) = (G%mask2dT(i,j)*tv%T(i,j,k) + G%mask2dT(i,j+1)*tv%T(i,j+1,k)) / & - (G%mask2dT(i,j) + G%mask2dT(i,j+1)) - OBC_S_v(i,J,k) = (G%mask2dT(i,j)*tv%S(i,j,k) + G%mask2dT(i,j+1)*tv%S(i,j+1,k)) / & - (G%mask2dT(i,j) + G%mask2dT(i,j+1)) - else ! This probably shouldn't happen or maybe it doesn't matter? - OBC_T_v(i,J,k) = 0.5*(tv%T(i,j,k)+tv%T(i,j+1,k)) - OBC_S_v(i,J,k) = 0.5*(tv%S(i,j,k)+tv%S(i,j+1,k)) - endif - else - OBC_T_v(i,J,k) = 0.5*(tv%T(i,j,k)+tv%T(i,j+1,k)) - OBC_S_v(i,J,k) = 0.5*(tv%S(i,j,k)+tv%S(i,j+1,k)) - endif - enddo; enddo ; enddo - endif - - call pass_vector(OBC_T_u, OBC_T_v, G%Domain, To_All+SCALAR_PAIR, CGRID_NE) - call pass_vector(OBC_S_u, OBC_S_v, G%Domain, To_All+SCALAR_PAIR, CGRID_NE) - - call add_tracer_OBC_values("T", tracer_Reg, OBC_in_u=OBC_T_u, & - OBC_in_v=OBC_T_v) - call add_tracer_OBC_values("S", tracer_Reg, OBC_in_u=OBC_S_u, & - OBC_in_v=OBC_S_v) - do k=1,nz ; do j=js,je ; do I=is-1,ie - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) then - tv%T(i+1,j,k) = tv%T(i,j,k) ; tv%S(i+1,j,k) = tv%S(i,j,k) - elseif (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) then - tv%T(i,j,k) = tv%T(i+1,j,k) ; tv%S(i,j,k) = tv%S(i+1,j,k) - endif - enddo ; enddo ; enddo - do k=1,nz ; do J=js-1,je ; do i=is,ie - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) then - tv%T(i,j+1,k) = tv%T(i,j,k) ; tv%S(i,j+1,k) = tv%S(i,j,k) - elseif (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) then - tv%T(i,j,k) = tv%T(i,j+1,k) ; tv%S(i,j,k) = tv%S(i,j+1,k) - endif - enddo ; enddo ; enddo - endif - - do k=1,nz ; do j=js-1,je+1 ; do I=is-1,ie+1 - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E) h(i+1,j,k) = h(i,j,k) - if (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W) h(i,j,k) = h(i+1,j,k) - enddo ; enddo ; enddo - do k=1,nz ; do J=js-1,je+1 ; do i=is-1,ie+1 - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N) h(i,j+1,k) = h(i,j,k) - if (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S) h(i,j,k) = h(i,j+1,k) - enddo ; enddo ; enddo - -end subroutine set_Flather_Bdry_Conds -! ----------------------------------------------------------------------------- - ! ----------------------------------------------------------------------------- subroutine set_velocity_depth_max(G) type(ocean_grid_type), intent(inout) :: G diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 007804d62e..78939b1f75 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -93,8 +93,7 @@ module MOM_hor_visc use MOM_grid, only : ocean_grid_type use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_MEKE_types, only : MEKE_type -use MOM_variables, only : ocean_OBC_type, OBC_FLATHER_E, OBC_FLATHER_W -use MOM_variables, only : OBC_FLATHER_N, OBC_FLATHER_S +use MOM_open_boundary, only : ocean_OBC_type, OBC_FLATHER use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -299,10 +298,10 @@ 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 (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_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 ; endif ; endif if (.not.associated(CS)) call MOM_error(FATAL, & "MOM_hor_visc: Module must be initialized before it is used.") @@ -587,8 +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_E) .or. & - (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W)) diffu(I,j,k) = 0.0 + if ((OBC%OBC_kind_u(I,j) == OBC_FLATHER)) diffu(I,j,k) = 0.0 endif ; endif enddo ; enddo @@ -600,8 +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_N) .or. & - (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S)) diffv(I,j,k) = 0.0 + if ((OBC%OBC_kind_v(i,J) == OBC_FLATHER)) diffv(I,j,k) = 0.0 endif ; endif enddo ; enddo diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index ef46c1f8e0..92586a922f 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -64,7 +64,7 @@ module MOM_set_visc use MOM_io, only : vardesc, var_desc use MOM_restart, only : register_restart_field, MOM_restart_CS use MOM_variables, only : thermo_var_ptrs -use MOM_variables, only : vertvisc_type, ocean_OBC_type +use MOM_variables, only : vertvisc_type use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density, calculate_density_derivs diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index ffe469a694..cca1f7807a 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -11,12 +11,13 @@ module MOM_vert_friction use MOM_forcing_type, only : forcing use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type +use MOM_open_boundary, only : ocean_OBC_type, OBC_SIMPLE use MOM_PointAccel, only : write_u_accel, write_v_accel, PointAccel_init use MOM_PointAccel, only : PointAccel_CS use MOM_time_manager, only : time_type, time_type_to_real, operator(-) use MOM_variables, only : thermo_var_ptrs, vertvisc_type use MOM_variables, only : cont_diag_ptrs, accel_diag_ptrs -use MOM_variables, only : ocean_internal_state, ocean_OBC_type, OBC_SIMPLE +use MOM_variables, only : ocean_internal_state use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -387,7 +388,7 @@ 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 (associated(OBC)) then ; if (OBC%this_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)) & @@ -400,7 +401,7 @@ subroutine vertvisc(u, v, h, fluxes, visc, dt, OBC, ADp, CDp, G, GV, CS, & v(i,J,k) = OBC%v(i,J,k) enddo ; enddo ; enddo endif - endif + endif ; endif ! Offer diagnostic fields for averaging. if (CS%id_du_dt_visc > 0) & call post_data(CS%id_du_dt_visc, ADp%du_dt_visc, CS%diag) diff --git a/src/tracer/DOME_tracer.F90 b/src/tracer/DOME_tracer.F90 index e7dcb22836..993b67eea9 100644 --- a/src/tracer/DOME_tracer.F90 +++ b/src/tracer/DOME_tracer.F90 @@ -62,13 +62,14 @@ module DOME_tracer use MOM_hor_index, only : hor_index_type use MOM_grid, only : ocean_grid_type use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : register_restart_field, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS use MOM_time_manager, only : time_type, get_time use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values use MOM_tracer_registry, only : tracer_vertdiff -use MOM_variables, only : surface, ocean_OBC_type +use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type use coupler_util, only : set_coupler_values, ind_csurf diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index e63037286e..6b826f40d5 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -30,7 +30,8 @@ module ISOMIP_tracer use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values use MOM_tracer_registry, only : tracer_vertdiff -use MOM_variables, only : surface, ocean_OBC_type +use MOM_variables, only : surface +use MOM_open_boundary, only : ocean_OBC_type use MOM_verticalGrid, only : verticalGrid_type use coupler_util, only : set_coupler_values, ind_csurf diff --git a/src/tracer/MOM_OCMIP2_CFC.F90 b/src/tracer/MOM_OCMIP2_CFC.F90 index 41181f3ab8..8ad94af4c0 100644 --- a/src/tracer/MOM_OCMIP2_CFC.F90 +++ b/src/tracer/MOM_OCMIP2_CFC.F90 @@ -73,6 +73,7 @@ module MOM_OCMIP2_CFC use MOM_hor_index, only : hor_index_type use MOM_grid, only : ocean_grid_type use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS use MOM_time_manager, only : time_type, get_time @@ -80,7 +81,7 @@ module MOM_OCMIP2_CFC use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values use MOM_tracer_registry, only : tracer_vertdiff use MOM_tracer_Z_init, only : tracer_Z_init -use MOM_variables, only : surface, ocean_OBC_type +use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type use coupler_util, only : extract_coupler_values, set_coupler_values diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 26ed11fe75..9686dc7eb2 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -14,9 +14,9 @@ module MOM_tracer_advect use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, 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_DIRECTION_E +use MOM_open_boundary, only : OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chksum -use MOM_variables, only : ocean_OBC_type, OBC_FLATHER_E -use MOM_variables, only : OBC_FLATHER_W, OBC_FLATHER_N, OBC_FLATHER_S use MOM_verticalGrid, only : verticalGrid_type implicit none ; private @@ -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%apply_OBC_u) then + if (associated(OBC)) then ; if (OBC%this_pe) then ; if (OBC%apply_OBC_u) then do_any_i = .false. do I=is-1,ie do_i(I) = .false. @@ -478,9 +478,9 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & ! Tracer fluxes are set to prescribed values only for inflows ! from masked areas. if (((uhr(I,j,k) > 0.0) .and. ((G%mask2dT(i,j) < 0.5) .or. & - (OBC%OBC_kind_u(I,j) == OBC_FLATHER_W))) .or. & + (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_W))) .or. & ((uhr(I,j,k) < 0.0) .and. ((G%mask2dT(i+1,j) < 0.5) .or. & - (OBC%OBC_kind_u(I,j) == OBC_FLATHER_E))) ) then + (OBC%OBC_direction_u(I,j) == OBC_DIRECTION_E))) ) then do_i(I) = .true. ; do_any_i = .true. uhh(I) = uhr(I,j,k) endif @@ -491,7 +491,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & flux_x(I,m) = uhh(I)*Tr(m)%OBC_in_u(I,j,k) else ; flux_x(I,m) = uhh(I)*Tr(m)%OBC_inflow_conc ; endif endif ; enddo ; enddo ; endif - endif ; endif + endif ; endif ; endif ! Calculate new tracer concentration in each cell after accounting ! for the i-direction fluxes. @@ -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%apply_OBC_v) then + if (associated(OBC)) then ; if (OBC%this_pe) then ; if (OBC%apply_OBC_v) then do_any_i = .false. do i=is,ie do_i(i) = .false. @@ -738,9 +738,9 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & ! Tracer fluxes are set to prescribed values only for inflows ! from masked areas. if (((vhr(i,J,k) > 0.0) .and. ((G%mask2dT(i,j) < 0.5) .or. & - (OBC%OBC_kind_v(i,J) == OBC_FLATHER_S))) .or. & + (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_S))) .or. & ((vhr(i,J,k) < 0.0) .and. ((G%mask2dT(i,j+1) < 0.5) .or. & - (OBC%OBC_kind_v(i,J) == OBC_FLATHER_N))) ) then + (OBC%OBC_direction_v(i,J) == OBC_DIRECTION_N))) ) then do_i(i) = .true. ; do_any_i = .true. vhh(i,J) = vhr(i,J,k) endif @@ -751,7 +751,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & flux_y(i,m,J) = vhh(i,J)*Tr(m)%OBC_in_v(i,J,k) else ; flux_y(i,m,J) = vhh(i,J)*Tr(m)%OBC_inflow_conc ; endif endif ; enddo ; enddo ; endif - endif ; endif + endif ; endif ; endif else ! not domore_v. do i=is,ie ; vhh(i,J) = 0.0 ; enddo do m=1,ntr ; do i=is,ie ; flux_y(i,m,J) = 0.0 ; enddo ; enddo diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index f8299cf0a5..3d43765d9d 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -36,11 +36,12 @@ module MOM_tracer_flow_control use MOM_forcing_type, only : forcing, optics_type use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type +use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : MOM_restart_CS use MOM_sponge, only : sponge_CS use MOM_ALE_sponge, only : ALE_sponge_CS use MOM_tracer_registry, only : tracer_registry_type -use MOM_variables, only : surface, ocean_OBC_type, thermo_var_ptrs +use MOM_variables, only : surface, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type #include diff --git a/src/tracer/MOM_tracer_hor_diff.F90 b/src/tracer/MOM_tracer_hor_diff.F90 index ab4353b76b..91c4fdf1f7 100644 --- a/src/tracer/MOM_tracer_hor_diff.F90 +++ b/src/tracer/MOM_tracer_hor_diff.F90 @@ -22,8 +22,7 @@ module MOM_tracer_hor_diff use MOM_neutral_diffusion, only : neutral_diffusion_CS use MOM_neutral_diffusion, only : neutral_diffusion_calc_coeffs, neutral_diffusion use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chksum -use MOM_variables, only : ocean_OBC_type, thermo_var_ptrs, OBC_FLATHER_E -use MOM_variables, only : OBC_FLATHER_W, OBC_FLATHER_N, OBC_FLATHER_S +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type implicit none ; private diff --git a/src/tracer/advection_test_tracer.F90 b/src/tracer/advection_test_tracer.F90 index a490299e76..c5eae04fe2 100644 --- a/src/tracer/advection_test_tracer.F90 +++ b/src/tracer/advection_test_tracer.F90 @@ -62,13 +62,14 @@ module advection_test_tracer use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : register_restart_field, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS use MOM_time_manager, only : time_type, get_time use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values use MOM_tracer_registry, only : tracer_vertdiff -use MOM_variables, only : surface, ocean_OBC_type +use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type use coupler_util, only : set_coupler_values, ind_csurf diff --git a/src/tracer/dye_example.F90 b/src/tracer/dye_example.F90 index 489abebb32..4c4d7cc6bb 100644 --- a/src/tracer/dye_example.F90 +++ b/src/tracer/dye_example.F90 @@ -63,6 +63,7 @@ module regional_dyes use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS use MOM_time_manager, only : time_type, get_time @@ -70,7 +71,7 @@ module regional_dyes use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values use MOM_tracer_registry, only : tracer_vertdiff use MOM_tracer_Z_init, only : tracer_Z_init -use MOM_variables, only : surface, ocean_OBC_type +use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type use coupler_util, only : set_coupler_values, ind_csurf diff --git a/src/tracer/ideal_age_example.F90 b/src/tracer/ideal_age_example.F90 index 9f7be83e04..04bfa4532d 100644 --- a/src/tracer/ideal_age_example.F90 +++ b/src/tracer/ideal_age_example.F90 @@ -63,6 +63,7 @@ module ideal_age_example use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS use MOM_time_manager, only : time_type, get_time @@ -70,7 +71,7 @@ module ideal_age_example use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values use MOM_tracer_registry, only : tracer_vertdiff use MOM_tracer_Z_init, only : tracer_Z_init -use MOM_variables, only : surface, ocean_OBC_type +use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type use coupler_util, only : set_coupler_values, ind_csurf diff --git a/src/tracer/oil_tracer.F90 b/src/tracer/oil_tracer.F90 index 72b90bd135..d2725cdb5d 100644 --- a/src/tracer/oil_tracer.F90 +++ b/src/tracer/oil_tracer.F90 @@ -63,6 +63,7 @@ module oil_tracer use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS use MOM_time_manager, only : time_type, get_time @@ -70,7 +71,7 @@ module oil_tracer use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values use MOM_tracer_registry, only : tracer_vertdiff use MOM_tracer_Z_init, only : tracer_Z_init -use MOM_variables, only : surface, ocean_OBC_type +use MOM_variables, only : surface use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use coupler_util, only : set_coupler_values, ind_csurf diff --git a/src/tracer/tracer_example.F90 b/src/tracer/tracer_example.F90 index 3b63b32056..8b5218276e 100644 --- a/src/tracer/tracer_example.F90 +++ b/src/tracer/tracer_example.F90 @@ -59,12 +59,13 @@ module USER_tracer_example use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc +use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : register_restart_field, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS use MOM_time_manager, only : time_type, get_time use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values -use MOM_variables, only : surface, ocean_OBC_type +use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type use coupler_util, only : set_coupler_values, ind_csurf diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index 36c3996bc4..fab4e2f2dc 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -10,7 +10,7 @@ module DOME2d_initialization use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE use MOM_io, only : write_field, slasher, vardesc use MOM_sponge, only : sponge_CS, set_up_sponge_field, initialize_sponge -use MOM_variables, only : thermo_var_ptrs, ocean_OBC_type +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 diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 3d5e2c87a5..df686d546e 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -25,8 +25,9 @@ module DOME_initialization 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_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_SIMPLE use MOM_tracer_registry, only : tracer_registry_type, add_tracer_OBC_values -use MOM_variables, only : thermo_var_ptrs, ocean_OBC_type, OBC_NONE, OBC_SIMPLE +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type @@ -37,7 +38,8 @@ module DOME_initialization public DOME_initialize_topography public DOME_initialize_thickness public DOME_initialize_sponges -public DOME_set_Open_Bdry_Conds +public DOME_set_OBC_positions +public DOME_set_OBC_data contains @@ -230,12 +232,40 @@ subroutine DOME_initialize_sponges(G, GV, tv, PF, CSp) endif 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_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) +subroutine DOME_set_OBC_data(OBC, tv, G, GV, param_file, tr_Reg) type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies !! whether, where, and what open boundary !! conditions are used. @@ -249,12 +279,6 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) !! to parse for model parameter values. type(tracer_registry_type), pointer :: tr_Reg !< Tracer registry. - logical :: any_OBC ! Set to true if any points in this subdomain use - ! open boundary conditions. - logical, pointer, dimension(:,:) :: & - OBC_mask_u => NULL(), & ! These arrays are true at zonal or meridional - OBC_mask_v => NULL() ! velocity points that have prescribed open boundary - ! conditions. 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 @@ -276,7 +300,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) ! thickness D_edge, in the same units as lat. real :: Ri_trans ! The shear Richardson number in the transition ! region of the specified shear profile. - character(len=40) :: mod = "DOME_set_Open_Bdry_Conds" ! This subroutine's name. + character(len=40) :: mod = "DOME_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 @@ -289,75 +313,29 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) Ri_trans = 1.0/3.0 ! The shear Richardson number in the transition region ! region of the specified shear profile. - call get_param(param_file, mod, "APPLY_OBC_U", 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", 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.) - - if (apply_OBC_u) then - ! Determine where u points are applied. - allocate(OBC_mask_u(IsdB:IedB,jsd:jed)) ; OBC_mask_u(:,:) = .false. - any_OBC = .false. - do j=jsd,jed ; do I=IsdB,IedB - ! if (SOME_TEST_FOR_U_OPEN_BCS) then - ! OBC_mask_u(I,j) = .true. ; any_OBC = .true. - ! endif - enddo ; enddo - if (.not.any_OBC) then - ! This processor has no u points at which open boundary conditions are - ! to be applied. - apply_OBC_u = .false. - deallocate(OBC_mask_u) - endif - endif - if (apply_OBC_v) then - ! Determine where v points are applied. - allocate(OBC_mask_v(isd:ied,JsdB:JedB)) ; OBC_mask_v(:,:) = .false. - any_OBC = .false. - do J=JsdB,JedB ; do i=isd,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_mask_v(i,J) = .true. ; any_OBC = .true. - endif - enddo ; enddo - if (.not.any_OBC) then - ! This processor has no v points at which open boundary conditions are - ! to be applied. - apply_OBC_v = .false. - deallocate(OBC_mask_v) - endif - endif - - if (.not.(apply_OBC_u .or. apply_OBC_v)) return + if (.not.associated(OBC)) return + if (.not.(OBC%apply_OBC_u .or. OBC%apply_OBC_v)) return - if (.not.associated(OBC)) allocate(OBC) - - if (apply_OBC_u) then - OBC%apply_OBC_u = .true. - OBC%OBC_mask_u => OBC_mask_u + 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 (apply_OBC_v) then - OBC%apply_OBC_v = .true. - OBC%OBC_mask_v => OBC_mask_v + 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 (apply_OBC_v) then + if (OBC%apply_OBC_v) then g_prime_tot = (GV%g_Earth/GV%Rho0)*2.0 Def_Rad = sqrt(D_edge*g_prime_tot) / (1.0e-4*1000.0) tr_0 = (-D_edge*sqrt(D_edge*g_prime_tot)*0.5e3*Def_Rad) * GV%m_to_H @@ -380,7 +358,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) if (k == nz) tr_k = tr_k + tr_0 * (2.0/(Ri_trans*(2.0+Ri_trans))) * & log((2.0+Ri_trans)/(2.0-Ri_trans)) do J=JsdB,JedB ; do i=isd,ied - if (OBC_mask_v(i,J)) then + if (OBC%OBC_mask_v(i,J)) then ! This needs to be unneccesarily complicated without symmetric memory. lon_im1 = 2.0*G%geoLonCv(i,J) - G%geoLonBu(I,J) ! if (isd > IsdB) lon_im1 = G%geoLonBu(I-1,J) @@ -394,9 +372,9 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) enddo endif - if (apply_OBC_u) then + if (OBC%apply_OBC_u) then do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB - if (OBC_mask_u(I,j)) then + if (OBC%OBC_mask_u(I,j)) then ! An appropriate expression for the zonal inflow velocities and ! transports should go here. OBC%uh(I,j,k) = 0.0 * GV%m_to_H ; OBC%u(I,j,k) = 0.0 @@ -408,7 +386,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) ! The inflow values of temperature and salinity also need to be set here if ! these variables are used. The following code is just a naive example. - if (apply_OBC_u .or. apply_OBC_v) then + if (OBC%apply_OBC_u .or. OBC%apply_OBC_v) then if (associated(tv%S)) then ! In this example, all S inflows have values of 35 psu. call add_tracer_OBC_values("S", tr_Reg, OBC_inflow=35.0) @@ -428,13 +406,13 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) do k=1,nz ; T0(k) = T0(k) + (GV%Rlay(k)-rho_guess(k)) / drho_dT(k) ; enddo enddo - if (apply_OBC_u) then + if (OBC%apply_OBC_u) then allocate(OBC_T_u(IsdB:IedB,jsd:jed,nz)) do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB OBC_T_u(I,j,k) = T0(k) enddo ; enddo ; enddo endif - if (apply_OBC_v) then + if (OBC%apply_OBC_v) then allocate(OBC_T_v(isd:ied,JsdB:JedB,nz)) do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied OBC_T_v(i,J,k) = T0(k) @@ -445,8 +423,7 @@ subroutine DOME_set_Open_Bdry_Conds(OBC, tv, G, GV, param_file, tr_Reg) endif endif -end subroutine DOME_set_Open_Bdry_Conds -! ----------------------------------------------------------------------------- +end subroutine DOME_set_OBC_data !> \class DOME_initialization !! diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index c47671a2cc..4a21b31324 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -28,7 +28,7 @@ module ISOMIP_initialization 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, ocean_OBC_type +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 diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90 index 16970cb1a0..3f85c2ee06 100644 --- a/src/user/Phillips_initialization.F90 +++ b/src/user/Phillips_initialization.F90 @@ -28,10 +28,8 @@ module Phillips_initialization use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE use MOM_io, only : write_field, slasher use MOM_sponge, only : set_up_sponge_field, initialize_sponge, sponge_CS -use MOM_tracer_registry, only : tracer_registry_type, add_tracer_OBC_values +use MOM_tracer_registry, only : tracer_registry_type use MOM_variables, only : thermo_var_ptrs -use MOM_variables, only : ocean_OBC_type, OBC_NONE, OBC_SIMPLE -use MOM_variables, only : OBC_FLATHER_E, OBC_FLATHER_W, OBC_FLATHER_N, OBC_FLATHER_S use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type diff --git a/src/user/Rossby_front_2d_initialization.F90 b/src/user/Rossby_front_2d_initialization.F90 index 69a22ca757..29b3e7adf5 100644 --- a/src/user/Rossby_front_2d_initialization.F90 +++ b/src/user/Rossby_front_2d_initialization.F90 @@ -10,7 +10,7 @@ module Rossby_front_2d_initialization 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, ocean_OBC_type +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 diff --git a/src/user/adjustment_initialization.F90 b/src/user/adjustment_initialization.F90 index d8b2414df6..2831a5a47f 100644 --- a/src/user/adjustment_initialization.F90 +++ b/src/user/adjustment_initialization.F90 @@ -26,7 +26,7 @@ module adjustment_initialization 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 -use MOM_variables, only : thermo_var_ptrs, ocean_OBC_type +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 diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90 index e7953f86d1..6415db7a13 100644 --- a/src/user/benchmark_initialization.F90 +++ b/src/user/benchmark_initialization.F90 @@ -25,9 +25,8 @@ module benchmark_initialization 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_tracer_registry, only : tracer_registry_type, add_tracer_OBC_values +use MOM_tracer_registry, only : tracer_registry_type use MOM_variables, only : thermo_var_ptrs -use MOM_variables, only : ocean_OBC_type, OBC_NONE, OBC_SIMPLE use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type diff --git a/src/user/circle_obcs_initialization.F90 b/src/user/circle_obcs_initialization.F90 index 1c1502eaff..80e326674f 100644 --- a/src/user/circle_obcs_initialization.F90 +++ b/src/user/circle_obcs_initialization.F90 @@ -24,8 +24,8 @@ module circle_obcs_initialization 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_tracer_registry, only : tracer_registry_type, add_tracer_OBC_values -use MOM_variables, only : thermo_var_ptrs, ocean_OBC_type +use MOM_tracer_registry, only : tracer_registry_type +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type diff --git a/src/user/external_gwave_initialization.F90 b/src/user/external_gwave_initialization.F90 index 5ce6145a7a..cc1fe5f226 100644 --- a/src/user/external_gwave_initialization.F90 +++ b/src/user/external_gwave_initialization.F90 @@ -23,7 +23,7 @@ module external_gwave_initialization 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_tracer_registry, only : tracer_registry_type, add_tracer_OBC_values +use MOM_tracer_registry, only : tracer_registry_type use MOM_variables, only : thermo_var_ptrs implicit none ; private diff --git a/src/user/lock_exchange_initialization.F90 b/src/user/lock_exchange_initialization.F90 index 222e45dc98..71822a84c2 100644 --- a/src/user/lock_exchange_initialization.F90 +++ b/src/user/lock_exchange_initialization.F90 @@ -23,7 +23,7 @@ module lock_exchange_initialization 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_tracer_registry, only : tracer_registry_type, add_tracer_OBC_values +use MOM_tracer_registry, only : tracer_registry_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type diff --git a/src/user/seamount_initialization.F90 b/src/user/seamount_initialization.F90 index 79768bca5e..3fa5cb91c7 100644 --- a/src/user/seamount_initialization.F90 +++ b/src/user/seamount_initialization.F90 @@ -29,8 +29,8 @@ module seamount_initialization use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE use MOM_io, only : write_field, slasher, vardesc use MOM_sponge, only : set_up_sponge_field, initialize_sponge, sponge_CS -use MOM_tracer_registry, only : tracer_registry_type, add_tracer_OBC_values -use MOM_variables, only : thermo_var_ptrs, ocean_OBC_type +use MOM_tracer_registry, only : tracer_registry_type +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 diff --git a/src/user/sloshing_initialization.F90 b/src/user/sloshing_initialization.F90 index a24c8f845a..34bc7edc46 100644 --- a/src/user/sloshing_initialization.F90 +++ b/src/user/sloshing_initialization.F90 @@ -29,8 +29,8 @@ module sloshing_initialization use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE use MOM_io, only : write_field, slasher use MOM_sponge, only : set_up_sponge_field, initialize_sponge, sponge_CS -use MOM_tracer_registry, only : tracer_registry_type, add_tracer_OBC_values -use MOM_variables, only : thermo_var_ptrs, ocean_OBC_type +use MOM_tracer_registry, only : tracer_registry_type +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90 index a96b820c65..c0868ef543 100644 --- a/src/user/user_initialization.F90 +++ b/src/user/user_initialization.F90 @@ -27,10 +27,11 @@ module user_initialization 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 +use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_SIMPLE, OBC_FLATHER +use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S use MOM_sponge, only : set_up_sponge_field, initialize_sponge, sponge_CS use MOM_tracer_registry, only : tracer_registry_type, add_tracer_OBC_values -use MOM_variables, only : thermo_var_ptrs, ocean_OBC_type, OBC_NONE, OBC_SIMPLE -use MOM_variables, only : OBC_FLATHER_E, OBC_FLATHER_W, OBC_FLATHER_N, OBC_FLATHER_S +use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type implicit none ; private @@ -40,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_Open_Bdry_Conds, USER_set_rotation +public USER_set_OBC_positions, USER_set_OBC_data, USER_set_rotation logical :: first_call = .true. @@ -198,8 +199,25 @@ 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_Open_Bdry_Conds(OBC, tv, G, param_file, tr_Reg) +subroutine USER_set_OBC_data(OBC, tv, G, param_file, tr_Reg) type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies !! whether, where, and what open boundary !! conditions are used. @@ -213,12 +231,12 @@ subroutine USER_set_Open_Bdry_Conds(OBC, tv, G, param_file, tr_Reg) !! parameter values. type(tracer_registry_type), pointer :: tr_Reg !< Tracer registry. ! call MOM_error(FATAL, & -! "USER_initialization.F90, USER_set_Open_Bdry_Conds: " // & +! "USER_initialization.F90, USER_set_OBC_data: " // & ! "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_Open_Bdry_Conds +end subroutine USER_set_OBC_data subroutine USER_set_rotation(G, param_file) type(ocean_grid_type), intent(inout) :: G