Skip to content

Commit c673ce7

Browse files
Hallberg-NOAAHallberg-NOAA
committed
Add btloop_add_dyn_PF and refactor OBCs in btcalc
 Moved the code to calculate the dynamic pressure into the new routine btloop_add_dyn_PF to simplify btloop_find_PF. Also revised btcalc to use the BT_OBC_type to streamline the application of open boundary conditions in the calculation of frhatu and frhatv. All answers are bitwise identical and no public interfaces are changed.
1 parent a563fa1 commit c673ce7

File tree

1 file changed

+140
-114
lines changed

1 file changed

+140
-114
lines changed

‎src/core/MOM_barotropic.F90

Lines changed: 140 additions & 114 deletions
Original file line numberDiff line numberDiff line change
@@ -2456,17 +2456,8 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL
24562456
if (CS%dynamic_psurf .or. (.not.CS%BT_project_velocity)) then
24572457
! Estimate the change in the free surface height.
24582458
call btloop_eta_predictor(n, dtbt, ubt, vbt, eta, ubt_int, vbt_int, uhbt, vhbt, uhbt0, vhbt0, &
2459-
uhbt_int, vhbt_int, BTCL_u, BTCL_v, Datu, Datv, &
2460-
eta_IC, eta_src, eta_pred, isv, iev, jsv, jev, &
2461-
integral_BT_cont, use_BT_cont, G, US, CS)
2462-
endif
2463-
2464-
! Use the change in eta to determine an additional divergence damping due to the ice.
2465-
if (CS%dynamic_psurf) then
2466-
!$OMP do
2467-
do j=jsv-1,jev+1 ; do i=isv-1,iev+1
2468-
p_surf_dyn(i,j) = dyn_coef_eta(i,j) * (eta_pred(i,j) - eta(i,j))
2469-
enddo ; enddo
2459+
uhbt_int, vhbt_int, BTCL_u, BTCL_v, Datu, Datv, eta_IC, eta_src, eta_pred, &
2460+
isv, iev, jsv, jev, integral_BT_cont, use_BT_cont, G, US, CS)
24702461
endif
24712462

24722463
if (interp_eta_PF) then
@@ -2480,41 +2471,43 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL
24802471

24812472
v_first = (MOD(n+G%first_direction,2)==1)
24822473

2474+
! Determine the pressure force accelerations due to the updated eta anomalies.
24832475
if (CS%BT_project_velocity) then
24842476
call btloop_find_PF(PFu, PFv, isv, iev, jsv, jev, eta, eta_PF, &
2485-
gtot_N, gtot_S, gtot_E, gtot_W, p_surf_dyn, dgeo_de, &
2486-
find_etaav, wt_accel2(n), eta_sum, v_first, G, US, CS)
2477+
gtot_N, gtot_S, gtot_E, gtot_W, dgeo_de, find_etaav, &
2478+
wt_accel2(n), eta_sum, v_first, G, US, CS)
24872479
else
24882480
call btloop_find_PF(PFu, PFv, isv, iev, jsv, jev, eta_pred, eta_PF, &
2489-
gtot_N, gtot_S, gtot_E, gtot_W, p_surf_dyn, dgeo_de, &
2490-
find_etaav, wt_accel2(n), eta_sum, v_first, G, US, CS)
2481+
gtot_N, gtot_S, gtot_E, gtot_W, dgeo_de, find_etaav, &
2482+
wt_accel2(n), eta_sum, v_first, G, US, CS)
2483+
endif
2484+
2485+
! Use the change in eta to determine an additional divergence damping due to the ice strength.
2486+
if (CS%dynamic_psurf) then
2487+
call btloop_add_dyn_PF(PFu, PFv, eta_pred, eta, dyn_coef_eta, p_surf_dyn, &
2488+
isv, iev, jsv, jev, v_first, G, US, CS)
24912489
endif
24922490

24932491
if (v_first) then
24942492
! On odd-steps, update v first.
2495-
call btloop_update_v(dtbt, ubt, vbt, v_accel_bt, Cor_v, PFv, &
2496-
isv-1, iev+1, jsv-1, jev, f_4_v, &
2497-
bt_rem_v, BT_force_v, vbt_prev, Cor_ref_v, Rayleigh_v, &
2493+
call btloop_update_v(dtbt, ubt, vbt, v_accel_bt, Cor_v, PFv, isv-1, iev+1, jsv-1, jev, &
2494+
f_4_v, bt_rem_v, BT_force_v, vbt_prev, Cor_ref_v, Rayleigh_v, &
24982495
wt_accel(n), G, US, CS)
24992496

25002497
! Now update the zonal velocity.
2501-
call btloop_update_u(dtbt, ubt, vbt, u_accel_bt, Cor_u, PFu, &
2502-
isv-1, iev, jsv, jev, f_4_u, &
2503-
bt_rem_u, BT_force_u, ubt_prev, Cor_ref_u, Rayleigh_u, &
2498+
call btloop_update_u(dtbt, ubt, vbt, u_accel_bt, Cor_u, PFu, isv-1, iev, jsv, jev, &
2499+
f_4_u, bt_rem_u, BT_force_u, ubt_prev, Cor_ref_u, Rayleigh_u, &
25042500
wt_accel(n), G, US, CS)
25052501

25062502
else
25072503
! On even steps, update u first.
2508-
call btloop_update_u(dtbt, ubt, vbt, u_accel_bt, Cor_u, PFu, &
2509-
isv-1, iev, jsv-1, jev+1, f_4_u, &
2510-
bt_rem_u, BT_force_u, ubt_prev, Cor_ref_u, Rayleigh_u, &
2504+
call btloop_update_u(dtbt, ubt, vbt, u_accel_bt, Cor_u, PFu, isv-1, iev, jsv-1, jev+1, &
2505+
f_4_u, bt_rem_u, BT_force_u, ubt_prev, Cor_ref_u, Rayleigh_u, &
25112506
wt_accel(n), G, US, CS)
25122507
! Now update the meridional velocity.
2513-
call btloop_update_v(dtbt, ubt, vbt, v_accel_bt, Cor_v, PFv, &
2514-
isv, iev, jsv-1, jev, f_4_v, &
2515-
bt_rem_v, BT_force_v, vbt_prev, Cor_ref_v, Rayleigh_v, &
2516-
wt_accel(n), G, US, CS, &
2517-
Cor_bracket_bug=CS%use_old_coriolis_bracket_bug)
2508+
call btloop_update_v(dtbt, ubt, vbt, v_accel_bt, Cor_v, PFv, isv, iev, jsv-1, jev, &
2509+
f_4_v, bt_rem_v, BT_force_v, vbt_prev, Cor_ref_v, Rayleigh_v, &
2510+
wt_accel(n), G, US, CS, Cor_bracket_bug=CS%use_old_coriolis_bracket_bug)
25182511
endif
25192512

25202513
! Determine the transports based on the updated velocities when no OBCs are applied
@@ -2961,8 +2954,8 @@ subroutine btloop_eta_predictor(n, dtbt, ubt, vbt, eta, ubt_int, vbt_int, uhbt,
29612954
end subroutine btloop_eta_predictor
29622955

29632956
subroutine btloop_find_PF(PFu, PFv, isv, iev, jsv, jev, eta_PF_BT, eta_PF, &
2964-
gtot_N, gtot_S, gtot_E, gtot_W, p_surf_dyn, dgeo_de, &
2965-
find_etaav, wt_accel2_n, eta_sum, v_first, G, US, CS)
2957+
gtot_N, gtot_S, gtot_E, gtot_W, dgeo_de, find_etaav, &
2958+
wt_accel2_n, eta_sum, v_first, G, US, CS)
29662959
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
29672960
type(barotropic_CS), intent(inout) :: CS !< Barotropic control structure
29682961
real, dimension(SZIBW_(CS),SZJW_(CS)), intent(inout) :: &
@@ -2977,9 +2970,9 @@ subroutine btloop_find_PF(PFu, PFv, isv, iev, jsv, jev, eta_PF_BT, eta_PF, &
29772970
eta_PF_BT !< The eta array (either the SSH anomaly or column mass anomaly) that
29782971
!! determines the barotropic pressure force [H ~> m or kg m-2]
29792972
real, dimension(SZIW_(CS),SZJW_(CS)), intent(in) :: &
2980-
eta_PF !< A local copy of the 2-D eta field (either SSH anomaly or
2981-
!! column mass anomaly) that was used to calculate the input
2982-
!! pressure gradient accelerations [H ~> m or kg m-2].
2973+
eta_PF !< The input 2-D eta field (either SSH anomaly or column mass anomaly)
2974+
!! that was used to calculate the input pressure gradient
2975+
!! accelerations [H ~> m or kg m-2].
29832976
real, dimension(SZIW_(CS),SZJW_(CS)), intent(in) :: &
29842977
gtot_N !< The effective total reduced gravity used to relate free surface height
29852978
!! deviations to pressure forces (including GFS and baroclinic contributions)
@@ -3002,8 +2995,6 @@ subroutine btloop_find_PF(PFu, PFv, isv, iev, jsv, jev, eta_PF_BT, eta_PF, &
30022995
!! in the barotropic momentum equations half a grid-point to the west of a
30032996
!! thickness point [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
30042997
!! (See Hallberg, J Comp Phys 1997 for a discussion of gtot_E and gtot_W.)
3005-
real, dimension(SZIW_(CS),SZJW_(CS)), intent(in) :: &
3006-
p_surf_dyn !< A dynamic surface pressure under rigid ice [L2 T-2 ~> m2 s-2].
30072998
real, intent(in) :: dgeo_de !< The constant of proportionality between geopotential and
30082999
!! sea surface height [nondim]. It is of order 1, but for stability this
30093000
!! may be made larger than the physical problem would suggest.
@@ -3037,21 +3028,8 @@ subroutine btloop_find_PF(PFu, PFv, isv, iev, jsv, jev, eta_PF_BT, eta_PF, &
30373028
PFv(i,J) = (((eta_PF_BT(i,j)-eta_PF(i,j))*gtot_N(i,j)) - &
30383029
((eta_PF_BT(i,j+1)-eta_PF(i,j+1))*gtot_S(i,j+1))) * &
30393030
dgeo_de * CS%IdyCv(i,J)
3040-
enddo ; enddo
3041-
!$OMP end do nowait
3042-
3043-
if (CS%dynamic_psurf) then
3044-
!$OMP do schedule(static)
3045-
do j=js_u,je_u ; do I=isv-1,iev
3046-
PFu(I,j) = PFu(I,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * CS%IdxCu(I,j)
3047-
enddo ; enddo
3048-
!$OMP end do nowait
3049-
!$OMP do schedule(static)
3050-
do J=jsv-1,jev ; do i=is_v,ie_v
3051-
PFv(i,J) = PFv(i,J) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * CS%IdyCv(i,J)
3052-
enddo ; enddo
3053-
!$OMP end do nowait
3054-
endif
3031+
enddo ; enddo
3032+
!$OMP end do nowait
30553033

30563034
if (find_etaav .and. (abs(wt_accel2_n) > 0.0)) then
30573035
!$OMP do
@@ -3063,6 +3041,63 @@ subroutine btloop_find_PF(PFu, PFv, isv, iev, jsv, jev, eta_PF_BT, eta_PF, &
30633041

30643042
end subroutine btloop_find_PF
30653043

3044+
!> This routine adds a dynamic pressure force based on the temporal changes in the predicted value
3045+
!! of eta, perhaps as an effective divergence damping to emulate the rigidity of an ice-sheet.
3046+
subroutine btloop_add_dyn_PF(PFu, PFv, eta_pred, eta, dyn_coef_eta, p_surf_dyn, &
3047+
isv, iev, jsv, jev, v_first, G, US, CS)
3048+
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
3049+
type(barotropic_CS), intent(inout) :: CS !< Barotropic control structure
3050+
real, dimension(SZIBW_(CS),SZJW_(CS)), intent(inout) :: &
3051+
PFu !< The anomalous zonal pressure force acceleration [L T-2 ~> m s-2].
3052+
real, dimension(SZIW_(CS),SZJBW_(CS)), intent(inout) :: &
3053+
PFv !< The meridional pressure force acceleration [L T-2 ~> m s-2].
3054+
real, dimension(SZIW_(CS),SZJW_(CS)), intent(in) :: &
3055+
eta_pred !< The updated eta field (either SSH anomaly or column mass anomaly) that is
3056+
!! used to estimate the divergence that is to be damped [H ~> m or kg m-2].
3057+
real, dimension(SZIW_(CS),SZJW_(CS)), intent(in) :: &
3058+
eta !< The previous eta field (either SSH anomaly or column mass anomaly) that is
3059+
!! used to estimate the divergence that is to be damped [H ~> m or kg m-2].
3060+
real, dimension(SZIW_(CS),SZJW_(CS)), intent(in) :: &
3061+
dyn_coef_eta !< The coefficient relating the changes in eta to the dynamic surface pressure
3062+
!! under rigid ice [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
3063+
real, dimension(SZIW_(CS),SZJW_(CS)), intent(inout) :: &
3064+
p_surf_dyn !< A dynamic surface pressure under rigid ice [L2 T-2 ~> m2 s-2].
3065+
integer, intent(in) :: isv !< The starting i-index of eta being set in ths loop
3066+
integer, intent(in) :: iev !< The ending i-index of eta_pred being set in ths loop
3067+
integer, intent(in) :: jsv !< The starting j-index of eta_pred being set in ths loop
3068+
integer, intent(in) :: jev !< The ending j-index of eta_pred being set in ths loop
3069+
logical, intent(in) :: v_first !< If true, update the v-velocity first with the present loop iteration
3070+
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3071+
3072+
! Local variables
3073+
integer :: i, j, js_u, je_u, is_v, ie_v
3074+
3075+
! Ensure that the extra points used for the temporally staggered Coriolis terms are updated.
3076+
if (v_first) then
3077+
is_v = isv-1 ; ie_v = iev+1 ; js_u = jsv ; je_u = jev
3078+
else
3079+
is_v = isv ; ie_v = iev ; js_u = jsv-1 ; je_u = jev+1
3080+
endif
3081+
3082+
! Use the change in eta to estimate the flow divergence and dynamic pressure.
3083+
!$OMP do
3084+
do j=jsv-1,jev+1 ; do i=isv-1,iev+1
3085+
p_surf_dyn(i,j) = dyn_coef_eta(i,j) * (eta_pred(i,j) - eta(i,j))
3086+
enddo ; enddo
3087+
3088+
!$OMP do schedule(static)
3089+
do j=js_u,je_u ; do I=isv-1,iev
3090+
PFu(I,j) = PFu(I,j) + (p_surf_dyn(i,j) - p_surf_dyn(i+1,j)) * CS%IdxCu(I,j)
3091+
enddo ; enddo
3092+
!$OMP end do nowait
3093+
!$OMP do schedule(static)
3094+
do J=jsv-1,jev ; do i=is_v,ie_v
3095+
PFv(i,J) = PFv(i,J) + (p_surf_dyn(i,j) - p_surf_dyn(i,j+1)) * CS%IdyCv(i,J)
3096+
enddo ; enddo
3097+
!$OMP end do nowait
3098+
3099+
end subroutine btloop_add_dyn_PF
3100+
30663101
!> Update meridional velocity.
30673102
subroutine btloop_update_v(dtbt, ubt, vbt, v_accel_bt, &
30683103
Cor_v, PFv, is_v, ie_v, Js_v, Je_v, f_4_v, &
@@ -4229,9 +4264,9 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
42294264
real :: htot ! The sum of the layer thicknesses [H ~> m or kg m-2].
42304265
real :: Ihtot ! The inverse of htot [H-1 ~> m-1 or m2 kg-1].
42314266

4232-
logical :: use_default, test_dflt, apply_OBCs
4267+
logical :: use_default, test_dflt
42334268
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, i, j, k
4234-
integer :: iss, ies, n
4269+
integer :: is_v, ie_v, Js_v, Je_v
42354270

42364271
! This section interpolates thicknesses onto u & v grid points with the
42374272
! second order accurate estimate h = 2*(h+ * h-)/(h+ + h-).
@@ -4254,13 +4289,6 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
42544289
"btcalc: Inconsistent settings of optional arguments and hvel_scheme.")
42554290
endif
42564291

4257-
apply_OBCs = .false.
4258-
if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then
4259-
! Some open boundary condition points might be in this processor's symmetric
4260-
! computational domain.
4261-
apply_OBCs = (OBC%number_of_segments > 0)
4262-
endif ; endif ; endif
4263-
42644292
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
42654293
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
42664294
h_neglect = GV%H_subroundoff
@@ -4330,6 +4358,25 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
43304358
CS%frhatu(I,j,k) = CS%frhatu(I,j,k) * Ihatutot(I)
43314359
enddo ; enddo
43324360
endif
4361+
if (CS%BT_OBC%u_OBCs_on_PE) then
4362+
if (((j >= CS%BT_OBC%js_u_E_obc) .and. (j <= CS%BT_OBC%je_u_E_obc)) .or. &
4363+
((j >= CS%BT_OBC%js_u_W_obc) .and. (j <= CS%BT_OBC%je_u_W_obc))) then
4364+
do I=is-1,ie
4365+
if (CS%BT_OBC%u_OBC_type(I,j) > 0) then ! Eastern boundary condition
4366+
htot = h(i,j,1)
4367+
do k=2,nz ; htot = htot + h(i,j,k) ; enddo
4368+
Ihtot = G%mask2dCu(I,j) / (htot + h_neglect)
4369+
do k=1,nz ; CS%frhatu(I,j,k) = h(i,j,k) * Ihtot ; enddo
4370+
endif
4371+
if (CS%BT_OBC%u_OBC_type(I,j) < 0) then ! Western boundary condition
4372+
htot = h(i+1,j,1)
4373+
do k=2,nz ; htot = htot + h(i+1,j,k) ; enddo
4374+
Ihtot = G%mask2dCu(I,j) / (htot + h_neglect)
4375+
do k=1,nz ; CS%frhatu(I,j,k) = h(i+1,j,k) * Ihtot ; enddo
4376+
endif
4377+
enddo
4378+
endif
4379+
endif
43334380
enddo
43344381

43354382
!$OMP parallel do default(none) shared(is,ie,js,je,nz,CS,G,GV,h_v,h_neglect,h,use_default) &
@@ -4396,62 +4443,41 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
43964443
endif
43974444
enddo
43984445

4399-
if (apply_OBCs) then ; do n=1,OBC%number_of_segments ! Test for segment type?
4400-
if (.not. OBC%segment(n)%on_pe) cycle
4401-
if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
4402-
J = OBC%segment(n)%HI%JsdB
4403-
if ((J >= js-1) .and. (J <= je)) then
4404-
iss = max(is,OBC%segment(n)%HI%isd) ; ies = min(ie,OBC%segment(n)%HI%ied)
4405-
do i=iss,ies ; hatvtot(i) = h(i,j,1) ; enddo
4406-
do k=2,nz ; do i=iss,ies
4407-
hatvtot(i) = hatvtot(i) + h(i,j,k)
4408-
enddo ; enddo
4409-
do i=iss,ies
4410-
Ihatvtot(i) = G%mask2dCv(i,J) / (hatvtot(i) + h_neglect)
4411-
enddo
4412-
do k=1,nz ; do i=iss,ies
4446+
if (CS%BT_OBC%v_OBCs_on_PE) then
4447+
! Northern boundary conditions
4448+
is_v = max(is, CS%BT_OBC%is_v_N_obc) ; ie_v = min(ie, CS%BT_OBC%ie_v_N_obc)
4449+
Js_v = max(js-1, CS%BT_OBC%Js_v_N_obc) ; Je_v = min(je, CS%BT_OBC%Je_v_N_obc)
4450+
do J=Js_v,Je_v
4451+
do i=is_v,ie_v ; hatvtot(i) = h(i,j,1) ; enddo
4452+
do k=2,nz ; do i=is_v,ie_v
4453+
hatvtot(i) = hatvtot(i) + h(i,j,k)
4454+
enddo ; enddo
4455+
do i=is_v,ie_v
4456+
Ihatvtot(i) = G%mask2dCv(i,J) / (hatvtot(i) + h_neglect)
4457+
enddo
4458+
do k=1,nz ; do i=is_v,ie_v
4459+
if (CS%BT_OBC%v_OBC_type(i,J) > 0) & ! Northern boundary condition
44134460
CS%frhatv(i,J,k) = h(i,j,k) * Ihatvtot(i)
4414-
enddo ; enddo
4415-
endif
4416-
elseif (OBC%segment(n)%direction == OBC_DIRECTION_S) then
4417-
J = OBC%segment(n)%HI%JsdB
4418-
if ((J >= js-1) .and. (J <= je)) then
4419-
iss = max(is,OBC%segment(n)%HI%isd) ; ies = min(ie,OBC%segment(n)%HI%ied)
4420-
do i=iss,ies ; hatvtot(i) = h(i,j+1,1) ; enddo
4421-
do k=2,nz ; do i=iss,ies
4422-
hatvtot(i) = hatvtot(i) + h(i,j+1,k)
4423-
enddo ; enddo
4424-
do i=iss,ies
4425-
Ihatvtot(i) = G%mask2dCv(i,J) / (hatvtot(i) + h_neglect)
4426-
enddo
4427-
do k=1,nz ; do i=iss,ies
4461+
enddo ; enddo
4462+
enddo
4463+
4464+
! Southern boundary conditions
4465+
is_v = max(is, CS%BT_OBC%is_v_S_obc) ; ie_v = min(ie, CS%BT_OBC%ie_v_S_obc)
4466+
Js_v = max(js-1, CS%BT_OBC%Js_v_S_obc) ; Je_v = min(je, CS%BT_OBC%Je_v_S_obc)
4467+
do J=Js_v,Je_v
4468+
do i=is_v,ie_v ; hatvtot(i) = h(i,j+1,1) ; enddo
4469+
do k=2,nz ; do i=is_v,ie_v
4470+
hatvtot(i) = hatvtot(i) + h(i,j+1,k)
4471+
enddo ; enddo
4472+
do i=is_v,ie_v
4473+
Ihatvtot(i) = G%mask2dCv(i,J) / (hatvtot(i) + h_neglect)
4474+
enddo
4475+
do k=1,nz ; do i=is_v,ie_v
4476+
if (CS%BT_OBC%v_OBC_type(i,J) < 0) & ! Southern boundary condition
44284477
CS%frhatv(i,J,k) = h(i,j+1,k) * Ihatvtot(i)
4429-
enddo ; enddo
4430-
endif
4431-
elseif (OBC%segment(n)%direction == OBC_DIRECTION_E) then
4432-
I = OBC%segment(n)%HI%IsdB
4433-
if ((I >= is-1) .and. (I <= ie)) then
4434-
do j = max(js,OBC%segment(n)%HI%jsd), min(je,OBC%segment(n)%HI%jed)
4435-
htot = h(i,j,1)
4436-
do k=2,nz ; htot = htot + h(i,j,k) ; enddo
4437-
Ihtot = G%mask2dCu(I,j) / (htot + h_neglect)
4438-
do k=1,nz ; CS%frhatu(I,j,k) = h(i,j,k) * Ihtot ; enddo
4439-
enddo
4440-
endif
4441-
elseif (OBC%segment(n)%direction == OBC_DIRECTION_W) then
4442-
I = OBC%segment(n)%HI%IsdB
4443-
if ((I >= is-1) .and. (I <= ie)) then
4444-
do j = max(js,OBC%segment(n)%HI%jsd), min(je,OBC%segment(n)%HI%jed)
4445-
htot = h(i+1,j,1)
4446-
do k=2,nz ; htot = htot + h(i+1,j,k) ; enddo
4447-
Ihtot = G%mask2dCu(I,j) / (htot + h_neglect)
4448-
do k=1,nz ; CS%frhatu(I,j,k) = h(i+1,j,k) * Ihtot ; enddo
4449-
enddo
4450-
endif
4451-
else
4452-
call MOM_error(fatal, "btcalc encountered an OBC segment of indeterminate direction.")
4453-
endif
4454-
enddo ; endif
4478+
enddo ; enddo
4479+
enddo
4480+
endif
44554481

44564482
if (CS%debug) then
44574483
call uvchksum("btcalc frhat[uv]", CS%frhatu, CS%frhatv, G%HI, &

0 commit comments

Comments
 (0)