diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 4688e646c5..3042449314 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -4,6 +4,7 @@ module MOM_barotropic ! This file is part of MOM6. See LICENSE.md for the license. use MOM_checksums, only : chksum0 +use MOM_coms, only : any_across_PEs use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_debugging, only : hchksum, uvchksum use MOM_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field @@ -265,12 +266,18 @@ module MOM_barotropic !! velocities. The streaming band-pass filter must be turned on. logical :: use_wide_halos !< If true, use wide halos and march in during the !! barotropic time stepping for efficiency. + integer :: min_stencil !< The minimum stencil width to use with the wide halo iterations. + !! A nonzero value may reflect the distribution of OBC faces or it + !! may be useful for debugging purposes. logical :: clip_velocity !< If true, limit any velocity components that are !! are large enough for a CFL number to exceed !! CFL_trunc. This should only be used as a !! desperate debugging measure. logical :: debug !< If true, write verbose checksums for debugging purposes. - logical :: debug_bt !< If true, write verbose checksums for debugging purposes. + logical :: debug_bt !< If true, write verbose checksums from within the barotropic + !! time-stepping loop for debugging purposes. + logical :: debug_wide_halos !< If true, write the checksums on the full wide halos. Otherwise + !! only the output for the final computational domain is written. real :: vel_underflow !< Velocity components smaller than vel_underflow !! are set to 0 [L T-1 ~> m s-1]. real :: maxvel !< Velocity components greater than maxvel are @@ -756,9 +763,9 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, interp_eta_PF = associated(eta_PF_start) ! Figure out the fullest arrays that could be updated. - stencil = 1 + stencil = max(1, CS%min_stencil) if ((.not.use_BT_cont) .and. CS%Nonlinear_continuity .and. & - (CS%Nonlin_cont_update_period > 0)) stencil = 2 + (CS%Nonlin_cont_update_period > 0)) stencil = max(2, CS%min_stencil) find_etaav = present(etaav) @@ -2395,6 +2402,7 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL integer :: stencil ! The stencil size of the algorithm, often 1 or 2. integer :: err_count ! A counter to limit the volume of error messages written to stdout. integer :: i, j, n, is, ie, js, je + integer :: debug_halo ! The halo size to use for debugging checksums integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -2403,9 +2411,10 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL err_count = 0 ! Figure out the fullest arrays that could be updated. - stencil = 1 + stencil = max(1, CS%min_stencil) if ((.not.use_BT_cont) .and. CS%Nonlinear_continuity .and. (CS%Nonlin_cont_update_period > 0)) & - stencil = 2 + stencil = max(2, CS%min_stencil) + num_cycles = 1 if (CS%use_wide_halos) & num_cycles = min((is-CS%isdw) / stencil, (js-CS%jsdw) / stencil) @@ -2508,6 +2517,14 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL jsv = jsv+stencil ; jev = jev-stencil endif + ! Store the previous velocities for time-filtered transports and OBCs. + do j=jsv,jev ; do I=isv-2,iev+1 + ubt_prev(I,j) = ubt(I,j) + enddo ; enddo + do J=jsv-2,jev+1 ; do i=isv,iev + vbt_prev(i,J) = vbt(i,J) + enddo ; enddo + if (integral_BT_cont) then !$OMP parallel do default(shared) do j=jsv-1,jev+1 ; do I=isv-2,iev+1 @@ -2563,22 +2580,22 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL if (v_first) then ! On odd-steps, update v first. call btloop_update_v(dtbt, ubt, vbt, v_accel_bt, Cor_v, PFv, isv-1, iev+1, jsv-1, jev, & - f_4_v, bt_rem_v, BT_force_v, vbt_prev, Cor_ref_v, Rayleigh_v, & + f_4_v, bt_rem_v, BT_force_v, Cor_ref_v, Rayleigh_v, & wt_accel(n), G, US, CS) ! Now update the zonal velocity. call btloop_update_u(dtbt, ubt, vbt, u_accel_bt, Cor_u, PFu, isv-1, iev, jsv, jev, & - f_4_u, bt_rem_u, BT_force_u, ubt_prev, Cor_ref_u, Rayleigh_u, & + f_4_u, bt_rem_u, BT_force_u, Cor_ref_u, Rayleigh_u, & wt_accel(n), G, US, CS) else ! On even steps, update u first. call btloop_update_u(dtbt, ubt, vbt, u_accel_bt, Cor_u, PFu, isv-1, iev, jsv-1, jev+1, & - f_4_u, bt_rem_u, BT_force_u, ubt_prev, Cor_ref_u, Rayleigh_u, & + f_4_u, bt_rem_u, BT_force_u, Cor_ref_u, Rayleigh_u, & wt_accel(n), G, US, CS) ! Now update the meridional velocity. call btloop_update_v(dtbt, ubt, vbt, v_accel_bt, Cor_v, PFv, isv, iev, jsv-1, jev, & - f_4_v, bt_rem_v, BT_force_v, vbt_prev, Cor_ref_v, Rayleigh_v, & + f_4_v, bt_rem_v, BT_force_v, Cor_ref_v, Rayleigh_v, & wt_accel(n), G, US, CS, Cor_bracket_bug=CS%use_old_coriolis_bracket_bug) endif @@ -2632,23 +2649,24 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL ! This might need to be moved outside of the OMP do loop directives. if (CS%debug_bt) then write(mesg,'("BT vel update ",I4)') n - call uvchksum(trim(mesg)//" PF[uv]", PFu, PFv, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_T_to_m_s*US%s_to_T) - call uvchksum(trim(mesg)//" Cor_[uv]", Cor_u, Cor_v, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_T_to_m_s*US%s_to_T) - call uvchksum(trim(mesg)//" BT_force_[uv]", BT_force_u, BT_force_v, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_T_to_m_s*US%s_to_T) - call uvchksum(trim(mesg)//" BT_rem_[uv]", BT_rem_u, BT_rem_v, CS%debug_BT_HI, & - haloshift=iev-ie, scalar_pair=.true.) - call uvchksum(trim(mesg)//" [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_T_to_m_s) - call uvchksum(trim(mesg)//" [uv]bt_trans", ubt_trans, vbt_trans, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_T_to_m_s) - call uvchksum(trim(mesg)//" [uv]hbt", uhbt, vhbt, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%s_to_T*US%L_to_m**2*GV%H_to_m) + debug_halo = 0 ; if (CS%debug_wide_halos) debug_halo = iev - ie + call uvchksum(trim(mesg)//" PF[uv]", PFu, PFv, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_T_to_m_s*US%s_to_T) + call uvchksum(trim(mesg)//" Cor_[uv]", Cor_u, Cor_v, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_T_to_m_s*US%s_to_T) + call uvchksum(trim(mesg)//" BT_force_[uv]", BT_force_u, BT_force_v, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_T_to_m_s*US%s_to_T) + call uvchksum(trim(mesg)//" BT_rem_[uv]", BT_rem_u, BT_rem_v, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., scalar_pair=.true.) + call uvchksum(trim(mesg)//" [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_T_to_m_s) + call uvchksum(trim(mesg)//" [uv]bt_trans", ubt_trans, vbt_trans, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_T_to_m_s) + call uvchksum(trim(mesg)//" [uv]hbt", uhbt, vhbt, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%s_to_T*US%L_to_m**2*GV%H_to_m) if (integral_BT_cont) & - call uvchksum(trim(mesg)//" [uv]hbt_int", uhbt_int, vhbt_int, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_to_m**2*GV%H_to_m) + call uvchksum(trim(mesg)//" [uv]hbt_int", uhbt_int, vhbt_int, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_to_m**2*GV%H_to_m) endif ! Apply open boundary condition considerations to revise the updated velocities and transports. @@ -2685,11 +2703,11 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL !$OMP end do nowait if (CS%debug_bt) then - call uvchksum("BT [uv]hbt just after OBC", uhbt, vhbt, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%s_to_T*US%L_to_m**2*GV%H_to_m) + call uvchksum("BT [uv]hbt just after OBC", uhbt, vhbt, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%s_to_T*US%L_to_m**2*GV%H_to_m) if (integral_BT_cont) & call uvchksum("BT [uv]hbt_int just after OBC", uhbt_int, vhbt_int, CS%debug_BT_HI, & - haloshift=iev-ie, unscale=US%L_to_m**2*GV%H_to_m) + haloshift=debug_halo, symmetric=.true., unscale=US%L_to_m**2*GV%H_to_m) endif ! Update eta in a corrector step using the barotropic continuity equation. @@ -2711,9 +2729,9 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL if (CS%debug_bt) then write(mesg,'("BT step ",I4)') n - call uvchksum(trim(mesg)//" [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_T_to_m_s) - call hchksum(eta, trim(mesg)//" eta", CS%debug_BT_HI, haloshift=iev-ie, unscale=GV%H_to_MKS) + call uvchksum(trim(mesg)//" [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_T_to_m_s) + call hchksum(eta, trim(mesg)//" eta", CS%debug_BT_HI, haloshift=debug_halo, unscale=GV%H_to_MKS) endif ! Issue warnings if there are unphysical values of the sea surface height or total water column mass. @@ -3190,7 +3208,7 @@ end subroutine btloop_add_dyn_PF !> Update meridional velocity. subroutine btloop_update_v(dtbt, ubt, vbt, v_accel_bt, & Cor_v, PFv, is_v, ie_v, Js_v, Je_v, f_4_v, & - bt_rem_v, BT_force_v, vbt_prev, Cor_ref_v, Rayleigh_v, & + bt_rem_v, BT_force_v, Cor_ref_v, Rayleigh_v, & wt_accel_n, G, US, CS, Cor_bracket_bug) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(barotropic_CS), intent(inout) :: CS !< Barotropic control structure @@ -3216,8 +3234,6 @@ subroutine btloop_update_v(dtbt, ubt, vbt, v_accel_bt, & integer, intent(in) :: ie_v !< The ending i-index of the range of v-point values to calculate integer, intent(in) :: Js_v !< The starting j-index of the range of v-point values to calculate integer, intent(in) :: Je_v !< The ending j-index of the range of v-point values to calculate - real, dimension(SZIW_(CS),SZJBW_(CS)), intent(inout) :: & - vbt_prev !< The previous velocity, stored for time-filtered transports and OBCs [L T-1 ~> m s-1] real, dimension(SZIW_(CS),SZJBW_(CS)), intent(in) :: & bt_rem_v !< The fraction of the barotropic meridional velocity that !! remains after a time step, the rest being lost to bottom @@ -3265,7 +3281,6 @@ subroutine btloop_update_v(dtbt, ubt, vbt, v_accel_bt, & !$OMP do schedule(static) ! This updates the v-velocity, except at OBC points. do J=Js_v,Je_v ; do i=is_v,ie_v - vbt_prev(i,J) = vbt(i,J) vbt(i,J) = bt_rem_v(i,J) * (vbt(i,J) + & dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J))) if (abs(vbt(i,J)) < CS%vel_underflow) vbt(i,J) = 0.0 @@ -3290,7 +3305,7 @@ end subroutine btloop_update_v !> Update zonal velocity. subroutine btloop_update_u(dtbt, ubt, vbt, u_accel_bt, & Cor_u, PFu, Is_u, Ie_u, js_u, je_u, f_4_u, & - bt_rem_u, BT_force_u, ubt_prev, Cor_ref_u, Rayleigh_u, & + bt_rem_u, BT_force_u, Cor_ref_u, Rayleigh_u, & wt_accel_n, G, US, CS) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(barotropic_CS), intent(inout) :: CS !< Barotropic control structure @@ -3324,8 +3339,6 @@ subroutine btloop_update_u(dtbt, ubt, vbt, u_accel_bt, & real, dimension(SZIBW_(CS),SZJW_(CS)), intent(in) :: & BT_force_u !< The vertical average of all of the v-accelerations that are !! not explicitly included in the barotropic equation [L T-2 ~> m s-2]. - real, dimension(SZIBW_(CS),SZJW_(CS)), intent(inout) :: & - ubt_prev !< The previous velocity, stored for time-filtered transports and OBCs [L T-1 ~> m s-1] real, dimension(SZIBW_(CS),SZJW_(CS)), intent(in) :: & Cor_ref_u !< The meridional barotropic Coriolis acceleration due !! to the reference velocities [L T-2 ~> m s-2]. @@ -3347,7 +3360,6 @@ subroutine btloop_update_u(dtbt, ubt, vbt, u_accel_bt, & ((f_4_u(3,I,j) * vbt(i,J)) + (f_4_u(2,I,j) * vbt(i+1,J-1)))) - & Cor_ref_u(I,j) - ubt_prev(I,j) = ubt(I,j) ubt(I,j) = bt_rem_u(I,j) * (ubt(I,j) + & dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j))) if (abs(ubt(I,j)) < CS%vel_underflow) ubt(I,j) = 0.0 @@ -4022,6 +4034,9 @@ subroutine initialize_BT_OBC(OBC, BT_OBC, G, CS) ! converted to real numbers to work with the MOM6 halo update code [nondim] real :: OBC_sign ! A sign encoding the direction of the OBC being used at a point [nondim] real :: OBC_type ! A real copy of the integer encoding the type of OBC being used at a point [nondim] + logical :: reversed_OBCs ! True of there any OBCs in the opposite halo on this PE, e.g. points + ! with a southern OBC in a northern halo. + logical :: any_reversed_OBCs integer :: i, j, isdw, iedw, jsdw, jedw integer :: l_seg, Flather_OBC_in_halo @@ -4121,6 +4136,15 @@ subroutine initialize_BT_OBC(OBC, BT_OBC, G, CS) BT_OBC%u_OBCs_on_PE = ((BT_OBC%Is_u_E_obc <= iedw) .or. (BT_OBC%Is_u_W_obc <= iedw)) BT_OBC%v_OBCs_on_PE = ((BT_OBC%is_v_N_obc <= iedw) .or. (BT_OBC%is_v_S_obc <= iedw)) + + ! Determine whether there are any OBCs in the opposite halo on any processors in the domain, e.g., + ! points with OBC_DIRECTION_S in a northern halo. + reversed_OBCs = (BT_OBC%u_OBCs_on_PE .and. ((BT_OBC%Is_u_E_obc <= G%isc-1) .or. (BT_OBC%Ie_u_W_obc >= G%iec))) .or. & + (BT_OBC%v_OBCs_on_PE .and. ((BT_OBC%Js_v_N_obc <= G%jsc-1) .or. (BT_OBC%Je_v_S_obc >= G%jec))) + any_reversed_OBCs = any_across_PEs(reversed_OBCs) + if (any_reversed_OBCs) call MOM_mesg("OBCs in an opposite halo require the use of a wider stencil.", 5) + if (any_reversed_OBCs) CS%min_stencil = max(CS%min_stencil, 2) + ! Allocate time-varying arrays that will be used for open boundary conditions. ! This pair is used with either Flather or specified OBCs. @@ -5428,6 +5452,11 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & call get_param(param_file, mdl, "BTHALO", bt_halo_sz, & "The minimum halo size for the barotropic solver.", default=0, & layoutParam=.true.) + call get_param(param_file, mdl, "BT_WIDE_HALO_MIN_STENCIL", CS%min_stencil, & + "The minimum stencil width to use with the wide halo iterations. "//& + "A nonzero value may be useful for debugging purposes, but at the "//& + "cost of reducing the efficiency gain from BT_USE_WIDE_HALOS.", & + default=0, layoutParam=.true., do_not_log=.not.CS%use_wide_halos) #ifdef STATIC_MEMORY_ if ((bt_halo_sz > 0) .and. (bt_halo_sz /= BTHALO_)) call MOM_error(FATAL, & "barotropic_init: Run-time values of BTHALO must agree with the "//& @@ -5646,6 +5675,12 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & "barotropic time-stepping loop. The data volume can be "//& "quite large if this is true.", default=CS%debug, & debuggingParam=.true.) + call get_param(param_file, mdl, "DEBUG_BT_WIDE_HALOS", CS%debug_wide_halos, & + "If true, write the checksums on the full wide halos. Otherwise only the "//& + "output for the final computational domain is written. This can be valuable "//& + "for debugging certain cases where the stencil used in the wide halo "//& + "iterations depends on which opoen boundary conditions are in the halos.", & + default=.true., do_not_log=.not.(CS%debug_bt.and.CS%use_wide_halos), debuggingParam=.true.) call get_param(param_file, mdl, "LINEARIZED_BT_CORIOLIS", CS%linearized_BT_PV, & "If true use the bottom depth instead of the total water column thickness "//&