Skip to content

Commit 5e4f97b

Browse files
kshedstromkshedstrom
authored
Fix interior obcs (#814)
* Adding a flag to turn off the bug fix (mostly) * Fix OBC check if not associated.
1 parent 9e7cfe9 commit 5e4f97b

File tree

6 files changed

+115
-19
lines changed

6 files changed

+115
-19
lines changed

‎src/core/MOM_barotropic.F90‎

Lines changed: 58 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -156,6 +156,7 @@ module MOM_barotropic
156156

157157
real, allocatable :: frhatu1(:,:,:) !< Predictor step values of frhatu stored for diagnostics [nondim]
158158
real, allocatable :: frhatv1(:,:,:) !< Predictor step values of frhatv stored for diagnostics [nondim]
159+
real, allocatable :: IareaT_OBCmask(:,:) !< If non-zero, work on given points [L-2 ~> m-2].
159160

160161
type(BT_OBC_type) :: BT_OBC !< A structure with all of this modules fields
161162
!! for applying open boundary conditions.
@@ -290,6 +291,8 @@ module MOM_barotropic
290291
!! consistent with tidal self-attraction and loading
291292
!! used within the barotropic solver
292293
logical :: wt_uv_bug = .true. !< If true, recover a bug that wt_[uv] that is not normalized.
294+
logical :: exterior_OBC_bug = .true. !< If true, recover a bug with boundary conditions
295+
!! inside the domain.
293296
type(time_type), pointer :: Time => NULL() !< A pointer to the ocean models clock.
294297
type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to regulate
295298
!! the timing of diagnostic output.
@@ -1961,7 +1964,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
19611964
enddo ; enddo
19621965
!$OMP do
19631966
do j=jsv-1,jev+1 ; do i=isv-1,iev+1
1964-
eta_pred(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + CS%IareaT(i,j) * &
1967+
eta_pred(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + CS%IareaT_OBCmask(i,j) * &
19651968
((uhbt_int(I-1,j) - uhbt_int(I,j)) + (vhbt_int(i,J-1) - vhbt_int(i,J)))
19661969
enddo ; enddo
19671970
elseif (use_BT_cont) then
@@ -1975,13 +1978,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
19751978
enddo ; enddo
19761979
!$OMP do
19771980
do j=jsv-1,jev+1 ; do i=isv-1,iev+1
1978-
eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * &
1981+
eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT_OBCmask(i,j)) * &
19791982
((uhbt(I-1,j) - uhbt(I,j)) + (vhbt(i,J-1) - vhbt(i,J)))
19801983
enddo ; enddo
19811984
else
19821985
!$OMP do
19831986
do j=jsv-1,jev+1 ; do i=isv-1,iev+1
1984-
eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * &
1987+
eta_pred(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT_OBCmask(i,j)) * &
19851988
(((Datu(I-1,j)*ubt(I-1,j) + uhbt0(I-1,j)) - &
19861989
(Datu(I,j)*ubt(I,j) + uhbt0(I,j))) + &
19871990
((Datv(i,J-1)*vbt(i,J-1) + vhbt0(i,J-1)) - &
@@ -2488,14 +2491,14 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
24882491
if (integral_BT_cont) then
24892492
!$OMP do
24902493
do j=jsv,jev ; do i=isv,iev
2491-
eta(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + CS%IareaT(i,j) * &
2494+
eta(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + CS%IareaT_OBCmask(i,j) * &
24922495
((uhbt_int(I-1,j) - uhbt_int(I,j)) + (vhbt_int(i,J-1) - vhbt_int(i,J)))
24932496
eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n)
24942497
enddo ; enddo
24952498
else
24962499
!$OMP do
24972500
do j=jsv,jev ; do i=isv,iev
2498-
eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * &
2501+
eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT_OBCmask(i,j)) * &
24992502
((uhbt(I-1,j) - uhbt(I,j)) + (vhbt(i,J-1) - vhbt(i,J)))
25002503
eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n)
25012504
enddo ; enddo
@@ -3283,7 +3286,7 @@ subroutine set_up_BT_OBC(OBC, eta, SpV_avg, BT_OBC, BT_Domain, G, GV, US, CS, MS
32833286
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
32843287
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
32853288
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3286-
type(barotropic_CS), intent(in) :: CS !< Barotropic control structure
3289+
type(barotropic_CS), intent(inout) :: CS !< Barotropic control structure
32873290
integer, intent(in) :: halo !< The extra halo size to use here.
32883291
logical, intent(in) :: use_BT_cont !< If true, use the BT_cont_types to calculate
32893292
!! transports.
@@ -3339,6 +3342,7 @@ subroutine set_up_BT_OBC(OBC, eta, SpV_avg, BT_OBC, BT_Domain, G, GV, US, CS, MS
33393342
allocate(BT_OBC%vhbt(isdw:iedw,jsdw-1:jedw), source=0.0)
33403343
allocate(BT_OBC%vbt_outer(isdw:iedw,jsdw-1:jedw), source=0.0)
33413344
allocate(BT_OBC%SSH_outer_v(isdw:iedw,jsdw-1:jedw), source=0.0)
3345+
33423346
BT_OBC%is_alloced = .true.
33433347
call create_group_pass(BT_OBC%pass_uv, BT_OBC%ubt_outer, BT_OBC%vbt_outer, BT_Domain)
33443348
call create_group_pass(BT_OBC%pass_uhvh, BT_OBC%uhbt, BT_OBC%vhbt, BT_Domain)
@@ -3481,6 +3485,7 @@ subroutine destroy_BT_OBC(BT_OBC)
34813485
deallocate(BT_OBC%vhbt)
34823486
deallocate(BT_OBC%vbt_outer)
34833487
deallocate(BT_OBC%SSH_outer_v)
3488+
34843489
BT_OBC%is_alloced = .false.
34853490
endif
34863491
end subroutine destroy_BT_OBC
@@ -4473,7 +4478,7 @@ end subroutine bt_mass_source
44734478
!! barotropic calculation and initializes any barotropic fields that have not
44744479
!! already been initialized.
44754480
subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
4476-
restart_CS, calc_dtbt, BT_cont, SAL_CSp, HA_CSp)
4481+
restart_CS, calc_dtbt, BT_cont, OBC, SAL_CSp, HA_CSp)
44774482
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
44784483
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
44794484
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
@@ -4494,7 +4499,8 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
44944499
type(BT_cont_type), pointer :: BT_cont !< A structure with elements that describe the
44954500
!! effective open face areas as a function of
44964501
!! barotropic flow.
4497-
type(SAL_CS), target, optional :: SAL_CSp !< A pointer to the control structure of the
4502+
type(ocean_OBC_type), pointer :: OBC !< The open boundary condition structure.
4503+
type(SAL_CS), target, optional :: SAL_CSp !< A pointer to the control structure of the
44984504
!! SAL module.
44994505
type(harmonic_analysis_CS), target, optional :: HA_CSp !< A pointer to the control structure of the
45004506
!! harmonic analysis module
@@ -4692,6 +4698,10 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
46924698
"If true, recover a bug in barotropic solver that uses an unnormalized weight "//&
46934699
"function for vertical averages of baroclinic velocity and forcing. Default "//&
46944700
"of this flag is set by VISC_REM_BUG.", default=visc_rem_bug)
4701+
call get_param(param_file, mdl, "EXTERIOR_OBC_BUG", CS%exterior_OBC_bug, &
4702+
"If true, recover a bug in barotropic solver and other routines when "//&
4703+
"boundary contitions interior to the domain are used.", &
4704+
default=.true., do_not_log=.true.)
46954705
call get_param(param_file, mdl, "TIDES", use_tides, &
46964706
"If true, apply tidal momentum forcing.", default=.false.)
46974707
if (use_tides .and. present(HA_CSp)) CS%HA_CSp => HA_CSp
@@ -4934,11 +4944,49 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
49344944
ALLOC_(CS%IdyCv(CS%isdw:CS%iedw,CS%jsdw-1:CS%jedw)) ; CS%IdyCv(:,:) = 0.0
49354945
ALLOC_(CS%dy_Cu(CS%isdw-1:CS%iedw,CS%jsdw:CS%jedw)) ; CS%dy_Cu(:,:) = 0.0
49364946
ALLOC_(CS%dx_Cv(CS%isdw:CS%iedw,CS%jsdw-1:CS%jedw)) ; CS%dx_Cv(:,:) = 0.0
4947+
allocate(CS%IareaT_OBCmask(isdw:iedw,jsdw:jedw), source=0.0)
49374948
do j=G%jsd,G%jed ; do i=G%isd,G%ied
49384949
CS%IareaT(i,j) = G%IareaT(i,j)
49394950
CS%bathyT(i,j) = G%bathyT(i,j)
4951+
CS%IareaT_OBCmask(i,j) = CS%IareaT(i,j)
49404952
enddo ; enddo
49414953

4954+
! Update IareaT_OBCmask so that nothing changes outside of the OBC (problem for interior OBCs only)
4955+
if (associated(OBC) .and. (CS%exterior_OBC_bug .eqv. .false.)) then
4956+
if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then
4957+
do i=isd,ied
4958+
do j=jsd,jed
4959+
if (OBC%segnum_u(I-1,j) /= OBC_NONE) then
4960+
if (OBC%segment(OBC%segnum_u(I-1,j))%direction == OBC_DIRECTION_E) then
4961+
CS%IareaT_OBCmask(i,j) = 0.0
4962+
endif
4963+
endif
4964+
if (OBC%segnum_u(I,j) /= OBC_NONE) then
4965+
if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then
4966+
CS%IareaT_OBCmask(i,j) = 0.0
4967+
endif
4968+
endif
4969+
enddo
4970+
enddo
4971+
endif
4972+
if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then
4973+
do j=jsd,jed
4974+
do i=isd,ied
4975+
if (OBC%segnum_v(i,J-1) /= OBC_NONE) then
4976+
if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then
4977+
CS%IareaT_OBCmask(i,j) = 0.0
4978+
endif
4979+
endif
4980+
if (OBC%segnum_v(i,J) /= OBC_NONE) then
4981+
if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then
4982+
CS%IareaT_OBCmask(i,j) = 0.0
4983+
endif
4984+
endif
4985+
enddo
4986+
enddo
4987+
endif
4988+
endif
4989+
49424990
! Note: G%IdxCu & G%IdyCv may be valid for a smaller extent than CS%IdxCu & CS%IdyCv, even without
49434991
! wide halos.
49444992
do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB
@@ -4949,6 +4997,7 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
49494997
enddo ; enddo
49504998
call create_group_pass(pass_static_data, CS%IareaT, CS%BT_domain, To_All)
49514999
call create_group_pass(pass_static_data, CS%bathyT, CS%BT_domain, To_All)
5000+
call create_group_pass(pass_static_data, CS%IareaT_OBCmask, CS%BT_domain, To_All)
49525001
call create_group_pass(pass_static_data, CS%IdxCu, CS%IdyCv, CS%BT_domain, To_All+Scalar_Pair)
49535002
call create_group_pass(pass_static_data, CS%dy_Cu, CS%dx_Cv, CS%BT_domain, To_All+Scalar_Pair)
49545003
call do_group_pass(pass_static_data, CS%BT_domain)
@@ -5302,6 +5351,7 @@ subroutine barotropic_end(CS)
53025351

53035352
if (allocated(CS%frhatu1)) deallocate(CS%frhatu1)
53045353
if (allocated(CS%frhatv1)) deallocate(CS%frhatv1)
5354+
if (allocated(CS%IareaT_OBCmask)) deallocate(CS%IareaT_OBCmask)
53055355
call deallocate_MOM_domain(CS%BT_domain)
53065356

53075357
! Allocated in restart registration, prior to timestep initialization

‎src/core/MOM_boundary_update.F90‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -145,9 +145,6 @@ subroutine update_OBC_data(OBC, G, GV, US, tv, h, CS, Time)
145145
type(update_OBC_CS), pointer :: CS !< Control structure for OBCs
146146
type(time_type), intent(in) :: Time !< Model time
147147

148-
! Something here... with CS%file_OBC_CSp?
149-
! if (CS%use_files) &
150-
! call update_OBC_segment_data(G, GV, OBC, tv, h, Time)
151148
if (CS%use_tidal_bay) &
152149
call tidal_bay_set_OBC_data(OBC, CS%tidal_bay_OBC, G, GV, US, h, Time)
153150
if (CS%use_Kelvin) &

‎src/core/MOM_dynamics_split_RK2.F90‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1554,7 +1554,8 @@ subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, p
15541554
do j=js,je ; do i=is,ie ; eta(i,j) = CS%eta(i,j) ; enddo ; enddo
15551555

15561556
call barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, &
1557-
CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, CS%SAL_CSp, HA_CSp)
1557+
CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, &
1558+
CS%OBC, CS%SAL_CSp, HA_CSp)
15581559

15591560
if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. &
15601561
.not. query_initialized(CS%diffv, "diffv", restart_CS)) then

‎src/core/MOM_dynamics_split_RK2b.F90‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1471,7 +1471,7 @@ subroutine initialize_dyn_split_RK2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US,
14711471

14721472
call barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, &
14731473
CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, &
1474-
CS%SAL_CSp, HA_CSp)
1474+
CS%OBC, CS%SAL_CSp, HA_CSp)
14751475

14761476
flux_units = get_flux_units(GV)
14771477
thickness_units = get_thickness_units(GV)

‎src/core/MOM_open_boundary.F90‎

Lines changed: 19 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -393,6 +393,7 @@ module MOM_open_boundary
393393
logical :: om4_remap_via_sub_cells !< If true, use the OM4 remapping algorithm
394394
character(40) :: remappingScheme !< String selecting the vertical remapping scheme
395395
type(group_pass_type) :: pass_oblique !< Structure for group halo pass
396+
logical :: exterior_OBC_bug !< If true, use incorrect form of tracers exterior to OBCs.
396397
end type ocean_OBC_type
397398

398399
!> Control structure for open boundaries that read from files.
@@ -561,6 +562,10 @@ subroutine open_boundary_config(G, US, param_file, OBC)
561562
"A silly value of velocities used outside of open boundary "//&
562563
"conditions for debugging.", units="m/s", default=0.0, scale=US%m_s_to_L_T, &
563564
do_not_log=.not.OBC%debug, debuggingParam=.true.)
565+
call get_param(param_file, mdl, "EXTERIOR_OBC_BUG", OBC%exterior_OBC_bug, &
566+
"If true, recover a bug in barotropic solver and other routines when "//&
567+
"boundary contitions interior to the domain are used.", &
568+
default=.true.)
564569
reentrant_x = .false.
565570
call get_param(param_file, mdl, "REENTRANT_X", reentrant_x, default=.true.)
566571
reentrant_y = .false.
@@ -5061,7 +5066,9 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC)
50615066
integer :: i, j
50625067
integer :: l_seg
50635068
logical :: fatal_error = .False.
5064-
real :: min_depth ! The minimum depth for ocean points [Z ~> m]
5069+
real :: min_depth ! The minimum depth for ocean points [Z ~> m]
5070+
real :: mask_depth ! The masking depth for ocean points [Z ~> m]
5071+
real :: Dmask ! The depth for masking in the same units as G%bathyT [Z ~> m].
50655072
integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2
50665073
character(len=256) :: mesg ! Message for error messages.
50675074
real, allocatable, dimension(:,:) :: color, color2 ! For sorting inside from outside,
@@ -5071,6 +5078,12 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC)
50715078

50725079
call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
50735080
units="m", default=0.0, scale=US%m_to_Z, do_not_log=.true.)
5081+
call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, &
5082+
units="m", default=-9999.0, scale=US%m_to_Z, do_not_log=.true.)
5083+
5084+
Dmask = mask_depth
5085+
if (mask_depth == -9999.0*US%m_to_Z) Dmask = min_depth
5086+
50745087
! The reference depth on a dyn_horgrid is 0, otherwise would need: min_depth = min_depth - G%Z_ref
50755088

50765089
allocate(color(G%isd:G%ied, G%jsd:G%jed), source=0.0)
@@ -5161,7 +5174,7 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC)
51615174
&"the masking of the outside grid points.")') i, j
51625175
call MOM_error(WARNING,"MOM mask_outside_OBCs: "//mesg, all_print=.true.)
51635176
endif
5164-
if (color(i,j) == cout) G%bathyT(i,j) = min_depth
5177+
if (color(i,j) == cout) G%bathyT(i,j) = Dmask
51655178
enddo ; enddo
51665179
if (fatal_error) call MOM_error(FATAL, &
51675180
"MOM_open_boundary: inconsistent OBC segments.")
@@ -5463,8 +5476,8 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
54635476
if (G%mask2dT(I+ishift,j) == 0.0) cycle
54645477
! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep
54655478
do m=1,segment%tr_Reg%ntseg
5466-
ntr_id = segment%tr_reg%Tr(m)%ntr_index
5467-
fd_id = segment%tr_reg%Tr(m)%fd_index
5479+
ntr_id = segment%tr_Reg%Tr(m)%ntr_index
5480+
fd_id = segment%tr_Reg%Tr(m)%fd_index
54685481
if (fd_id == -1) then
54695482
resrv_lfac_out = 1.0
54705483
resrv_lfac_in = 1.0
@@ -5507,8 +5520,8 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
55075520
if (G%mask2dT(i,j+jshift) == 0.0) cycle
55085521
! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep
55095522
do m=1,segment%tr_Reg%ntseg
5510-
ntr_id = segment%tr_reg%Tr(m)%ntr_index
5511-
fd_id = segment%tr_reg%Tr(m)%fd_index
5523+
ntr_id = segment%tr_Reg%Tr(m)%ntr_index
5524+
fd_id = segment%tr_Reg%Tr(m)%fd_index
55125525
if (fd_id == -1) then
55135526
resrv_lfac_out = 1.0
55145527
resrv_lfac_in = 1.0

‎src/tracer/MOM_tracer_advect.F90‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -648,6 +648,21 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
648648
endif
649649
enddo
650650

651+
! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only)
652+
if (associated(OBC)) then
653+
if ((OBC%exterior_OBC_bug .eqv. .false.) .and. (OBC%OBC_pe)) then
654+
if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then
655+
do i=is,ie-1 ; if (OBC%segnum_u(I,j) /= OBC_NONE) then
656+
if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then
657+
do_i(i+1,j) = .false.
658+
elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then
659+
do_i(i,j) = .false.
660+
endif
661+
endif ; enddo
662+
endif
663+
endif
664+
endif
665+
651666
! update tracer concentration from i-flux and save some diagnostics
652667
do m=1,ntr
653668

@@ -1039,6 +1054,26 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
10391054
else ; do_i(i,j) = .false. ; endif
10401055
enddo
10411056

1057+
! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only)
1058+
if (associated(OBC)) then
1059+
if ((OBC%exterior_OBC_bug .eqv. .false.) .and. (OBC%OBC_pe)) then
1060+
if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then
1061+
do i=is,ie
1062+
if (OBC%segnum_v(i,J-1) /= OBC_NONE) then
1063+
if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then
1064+
do_i(i,j) = .false.
1065+
endif
1066+
endif
1067+
if (OBC%segnum_v(i,J) /= OBC_NONE) then
1068+
if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then
1069+
do_i(i,j) = .false.
1070+
endif
1071+
endif
1072+
enddo
1073+
endif
1074+
endif
1075+
endif
1076+
10421077
! update tracer and save some diagnostics
10431078
do m=1,ntr
10441079
do i=is,ie ; if (do_i(i,j)) then

0 commit comments

Comments
 (0)