Skip to content

Commit 6539f69

Browse files
Hallberg-NOAAHallberg-NOAA
committed
Refactor btcalc to reduce duplicated code
 Refactored btcalc to include the application of open boundary conditions in the calculation of CS%frhatu and CS%frhatv within the same loops as the non-OBC calculations, and also to reduce the amount of duplicated code. This should increase the model efficiency with open boundary conditions, and have minimal performance impacts otherwise. All answers are bitwise identical and no public interfaces are changed.
1 parent c267ed6 commit 6539f69

File tree

1 file changed

+127
-160
lines changed

1 file changed

+127
-160
lines changed

‎src/core/MOM_barotropic.F90‎

Lines changed: 127 additions & 160 deletions
Original file line numberDiff line numberDiff line change
@@ -4279,11 +4279,8 @@ subroutine destroy_BT_OBC(BT_OBC)
42794279

42804280
end subroutine destroy_BT_OBC
42814281

4282-
!> btcalc calculates the barotropic velocities from the full velocity and
4283-
!! thickness fields, determines the fraction of the total water column in each
4284-
!! layer at velocity points, and determines a corrective fictitious mass source
4285-
!! that will drive the barotropic estimate of the free surface height toward the
4286-
!! baroclinic estimate.
4282+
!> btcalc determines the fraction of the total water column in each
4283+
!! layer at velocity points.
42874284
subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
42884285
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
42894286
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
@@ -4313,6 +4310,8 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
43134310
type(ocean_OBC_type), optional, pointer :: OBC !< Open boundary control structure.
43144311

43154312
! Local variables
4313+
real :: hatu(SZIB_(G),SZK_(GV)) ! The layer thicknesses interpolated to u points [H ~> m or kg m-2]
4314+
real :: hatv(SZI_(G),SZK_(GV)) ! The layer thicknesses interpolated to v points [H ~> m or kg m-2]
43164315
real :: hatutot(SZIB_(G)) ! The sum of the layer thicknesses interpolated to u points [H ~> m or kg m-2].
43174316
real :: hatvtot(SZI_(G)) ! The sum of the layer thicknesses interpolated to v points [H ~> m or kg m-2].
43184317
real :: Ihatutot(SZIB_(G)) ! Ihatutot is the inverse of hatutot [H-1 ~> m-1 or m2 kg-1].
@@ -4330,15 +4329,11 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
43304329
real :: D_shallow_v(SZIB_(G))! The height of the shallower of the adjacent bathymetric depths
43314330
! around a v-point (positive upward) [H ~> m or kg m-2]
43324331
real :: Z_to_H ! A local conversion factor [H Z-1 ~> nondim or kg m-3]
4333-
real :: htot ! The sum of the layer thicknesses [H ~> m or kg m-2].
4334-
real :: Ihtot ! The inverse of htot [H-1 ~> m-1 or m2 kg-1].
43354332

43364333
logical :: use_default, test_dflt
43374334
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, i, j, k
43384335
integer :: is_v, ie_v, Js_v, Je_v
43394336

4340-
! This section interpolates thicknesses onto u & v grid points with the
4341-
! second order accurate estimate h = 2*(h+ * h-)/(h+ + h-).
43424337
if (.not.CS%module_is_initialized) call MOM_error(FATAL, &
43434338
"btcalc: Module MOM_barotropic must be initialized before it is used.")
43444339

@@ -4362,191 +4357,163 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
43624357
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
43634358
h_neglect = GV%H_subroundoff
43644359

4365-
! This estimates the fractional thickness of each layer at the velocity
4366-
! points, using a harmonic mean estimate.
43674360

43684361
!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_u,CS,h_neglect,h,use_default,G,GV) &
4369-
!$OMP private(hatutot,Ihatutot,e_u,D_shallow_u,h_arith,h_harm,wt_arith,Z_to_H,htot,Ihtot)
4362+
!$OMP private(hatu,hatutot,Ihatutot,e_u,D_shallow_u,h_arith,h_harm,wt_arith,Z_to_H)
43704363
do j=js,je
4364+
do I=is-1,ie ; hatutot(I) = 0.0 ; enddo
4365+
43714366
if (present(h_u)) then
4372-
do I=is-1,ie ; hatutot(I) = h_u(I,j,1) ; enddo
4373-
do k=2,nz ; do I=is-1,ie
4374-
hatutot(I) = hatutot(I) + h_u(I,j,k)
4367+
do k=1,nz ; do I=is-1,ie
4368+
hatu(I,k) = h_u(I,j,k)
4369+
hatutot(I) = hatutot(I) + hatu(I,k)
43754370
enddo ; enddo
4376-
do I=is-1,ie ; Ihatutot(I) = G%mask2dCu(I,j) / (hatutot(I) + h_neglect) ; enddo
4371+
elseif (CS%hvel_scheme == ARITHMETIC) then
43774372
do k=1,nz ; do I=is-1,ie
4378-
CS%frhatu(I,j,k) = h_u(I,j,k) * Ihatutot(I)
4373+
hatu(I,k) = 0.5 * (h(i+1,j,k) + h(i,j,k))
4374+
hatutot(I) = hatutot(I) + hatu(I,k)
43794375
enddo ; enddo
4380-
else
4381-
if (CS%hvel_scheme == ARITHMETIC) then
4382-
do I=is-1,ie
4383-
CS%frhatu(I,j,1) = 0.5 * (h(i+1,j,1) + h(i,j,1))
4384-
hatutot(I) = CS%frhatu(I,j,1)
4385-
enddo
4386-
do k=2,nz ; do I=is-1,ie
4387-
CS%frhatu(I,j,k) = 0.5 * (h(i+1,j,k) + h(i,j,k))
4388-
hatutot(I) = hatutot(I) + CS%frhatu(I,j,k)
4389-
enddo ; enddo
4390-
elseif (CS%hvel_scheme == HYBRID .or. use_default) then
4391-
Z_to_H = GV%Z_to_H ; if (.not.GV%Boussinesq) Z_to_H = GV%RZ_to_H * CS%Rho_BT_lin
4392-
do I=is-1,ie
4393-
e_u(I,nz+1) = -0.5 * Z_to_H * (G%bathyT(i+1,j) + G%bathyT(i,j))
4394-
D_shallow_u(I) = -Z_to_H * min(G%bathyT(i+1,j), G%bathyT(i,j))
4395-
hatutot(I) = 0.0
4396-
enddo
4397-
do k=nz,1,-1 ; do I=is-1,ie
4398-
e_u(I,K) = e_u(I,K+1) + 0.5 * (h(i+1,j,k) + h(i,j,k))
4399-
h_arith = 0.5 * (h(i+1,j,k) + h(i,j,k))
4400-
if (e_u(I,K+1) >= D_shallow_u(I)) then
4401-
CS%frhatu(I,j,k) = h_arith
4376+
elseif (CS%hvel_scheme == HYBRID .or. use_default) then
4377+
Z_to_H = GV%Z_to_H ; if (.not.GV%Boussinesq) Z_to_H = GV%RZ_to_H * CS%Rho_BT_lin
4378+
do I=is-1,ie
4379+
e_u(I,nz+1) = -0.5 * Z_to_H * (G%bathyT(i+1,j) + G%bathyT(i,j))
4380+
D_shallow_u(I) = -Z_to_H * min(G%bathyT(i+1,j), G%bathyT(i,j))
4381+
enddo
4382+
do k=nz,1,-1 ; do I=is-1,ie
4383+
e_u(I,K) = e_u(I,K+1) + 0.5 * (h(i+1,j,k) + h(i,j,k))
4384+
h_arith = 0.5 * (h(i+1,j,k) + h(i,j,k))
4385+
if (e_u(I,K+1) >= D_shallow_u(I)) then
4386+
hatu(I,k) = h_arith
4387+
else
4388+
h_harm = (h(i+1,j,k) * h(i,j,k)) / (h_arith + h_neglect)
4389+
if (e_u(I,K) <= D_shallow_u(I)) then
4390+
hatu(I,k) = h_harm
44024391
else
4403-
h_harm = (h(i+1,j,k) * h(i,j,k)) / (h_arith + h_neglect)
4404-
if (e_u(I,K) <= D_shallow_u(I)) then
4405-
CS%frhatu(I,j,k) = h_harm
4406-
else
4407-
wt_arith = (e_u(I,K) - D_shallow_u(I)) / (h_arith + h_neglect)
4408-
CS%frhatu(I,j,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
4409-
endif
4392+
wt_arith = (e_u(I,K) - D_shallow_u(I)) / (h_arith + h_neglect)
4393+
hatu(I,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
44104394
endif
4411-
hatutot(I) = hatutot(I) + CS%frhatu(I,j,k)
4412-
enddo ; enddo
4413-
elseif (CS%hvel_scheme == HARMONIC) then
4414-
do I=is-1,ie
4415-
CS%frhatu(I,j,1) = 2.0*(h(i+1,j,1) * h(i,j,1)) / &
4416-
((h(i+1,j,1) + h(i,j,1)) + h_neglect)
4417-
hatutot(I) = CS%frhatu(I,j,1)
4418-
enddo
4419-
do k=2,nz ; do I=is-1,ie
4420-
CS%frhatu(I,j,k) = 2.0*(h(i+1,j,k) * h(i,j,k)) / &
4421-
((h(i+1,j,k) + h(i,j,k)) + h_neglect)
4422-
hatutot(I) = hatutot(I) + CS%frhatu(I,j,k)
4423-
enddo ; enddo
4424-
endif
4425-
do I=is-1,ie ; Ihatutot(I) = G%mask2dCu(I,j) / (hatutot(I) + h_neglect) ; enddo
4395+
endif
4396+
hatutot(I) = hatutot(I) + hatu(I,k)
4397+
enddo ; enddo
4398+
elseif (CS%hvel_scheme == HARMONIC) then
4399+
! Interpolates thicknesses onto u grid points with the
4400+
! second order accurate estimate h = 2*(h+ * h-)/(h+ + h-).
44264401
do k=1,nz ; do I=is-1,ie
4427-
CS%frhatu(I,j,k) = CS%frhatu(I,j,k) * Ihatutot(I)
4402+
hatu(I,k) = 2.0*(h(i+1,j,k) * h(i,j,k)) / &
4403+
((h(i+1,j,k) + h(i,j,k)) + h_neglect)
4404+
hatutot(I) = hatutot(I) + hatu(I,k)
44284405
enddo ; enddo
44294406
endif
4407+
44304408
if (CS%BT_OBC%u_OBCs_on_PE) then
4431-
if (((j >= CS%BT_OBC%js_u_E_obc) .and. (j <= CS%BT_OBC%je_u_E_obc)) .or. &
4432-
((j >= CS%BT_OBC%js_u_W_obc) .and. (j <= CS%BT_OBC%je_u_W_obc))) then
4433-
do I=is-1,ie
4409+
! Reset velocity point thicknesses and their sums at OBC points
4410+
if ((j >= CS%BT_OBC%js_u_E_obc) .and. (j <= CS%BT_OBC%je_u_E_obc)) then
4411+
do I = max(is-1,CS%BT_OBC%Is_u_E_obc), min(ie,CS%BT_OBC%Ie_u_E_obc)
44344412
if (CS%BT_OBC%u_OBC_type(I,j) > 0) then ! Eastern boundary condition
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
4413+
hatutot(I) = 0.0
4414+
do k=1,nz
4415+
hatu(I,k) = h(i,j,k)
4416+
hatutot(I) = hatutot(I) + hatu(I,k)
4417+
enddo
44394418
endif
4419+
enddo
4420+
endif
4421+
if ((j >= CS%BT_OBC%js_u_W_obc) .and. (j <= CS%BT_OBC%je_u_W_obc)) then
4422+
do I = max(is-1,CS%BT_OBC%Is_u_W_obc), min(ie,CS%BT_OBC%Ie_u_W_obc)
44404423
if (CS%BT_OBC%u_OBC_type(I,j) < 0) then ! Western boundary condition
4441-
htot = h(i+1,j,1)
4442-
do k=2,nz ; htot = htot + h(i+1,j,k) ; enddo
4443-
Ihtot = G%mask2dCu(I,j) / (htot + h_neglect)
4444-
do k=1,nz ; CS%frhatu(I,j,k) = h(i+1,j,k) * Ihtot ; enddo
4424+
hatutot(I) = 0.0
4425+
do k=1,nz
4426+
hatu(I,k) = h(i+1,j,k)
4427+
hatutot(I) = hatutot(I) + hatu(I,k)
4428+
enddo
44454429
endif
44464430
enddo
44474431
endif
44484432
endif
4433+
4434+
! Determine the fractional thickness of each layer at the velocity points.
4435+
do I=is-1,ie ; Ihatutot(I) = G%mask2dCu(I,j) / (hatutot(I) + h_neglect) ; enddo
4436+
do k=1,nz ; do I=is-1,ie
4437+
CS%frhatu(I,j,k) = hatu(I,k) * Ihatutot(I)
4438+
enddo ; enddo
44494439
enddo
44504440

44514441
!$OMP parallel do default(none) shared(is,ie,js,je,nz,CS,G,GV,h_v,h_neglect,h,use_default) &
4452-
!$OMP private(hatvtot,Ihatvtot,e_v,D_shallow_v,h_arith,h_harm,wt_arith,Z_to_H)
4442+
!$OMP private(hatv,hatvtot,Ihatvtot,e_v,D_shallow_v,h_arith,h_harm,wt_arith,Z_to_H)
44534443
do J=js-1,je
4444+
do i=is,ie ; hatvtot(i) = 0.0 ; enddo
44544445
if (present(h_v)) then
4455-
do i=is,ie ; hatvtot(i) = h_v(i,J,1) ; enddo
4456-
do k=2,nz ; do i=is,ie
4457-
hatvtot(i) = hatvtot(i) + h_v(i,J,k)
4446+
do k=1,nz ; do i=is,ie
4447+
hatv(i,k) = h_v(i,J,k)
4448+
hatvtot(i) = hatvtot(i) + hatv(i,k)
44584449
enddo ; enddo
4459-
do i=is,ie ; Ihatvtot(i) = G%mask2dCv(i,J) / (hatvtot(i) + h_neglect) ; enddo
4450+
elseif (CS%hvel_scheme == ARITHMETIC) then
44604451
do k=1,nz ; do i=is,ie
4461-
CS%frhatv(i,J,k) = h_v(i,J,k) * Ihatvtot(i)
4452+
hatv(i,k) = 0.5 * (h(i,j+1,k) + h(i,j,k))
4453+
hatvtot(i) = hatvtot(i) + hatv(i,k)
44624454
enddo ; enddo
4463-
else
4464-
if (CS%hvel_scheme == ARITHMETIC) then
4465-
do i=is,ie
4466-
CS%frhatv(i,J,1) = 0.5 * (h(i,j+1,1) + h(i,j,1))
4467-
hatvtot(i) = CS%frhatv(i,J,1)
4468-
enddo
4469-
do k=2,nz ; do i=is,ie
4470-
CS%frhatv(i,J,k) = 0.5 * (h(i,j+1,k) + h(i,j,k))
4471-
hatvtot(i) = hatvtot(i) + CS%frhatv(i,J,k)
4472-
enddo ; enddo
4473-
elseif (CS%hvel_scheme == HYBRID .or. use_default) then
4474-
Z_to_H = GV%Z_to_H ; if (.not.GV%Boussinesq) Z_to_H = GV%RZ_to_H * CS%Rho_BT_lin
4475-
do i=is,ie
4476-
e_v(i,nz+1) = -0.5 * Z_to_H * (G%bathyT(i,j+1) + G%bathyT(i,j))
4477-
D_shallow_v(I) = -Z_to_H * min(G%bathyT(i,j+1), G%bathyT(i,j))
4478-
hatvtot(I) = 0.0
4479-
enddo
4480-
do k=nz,1,-1 ; do i=is,ie
4481-
e_v(i,K) = e_v(i,K+1) + 0.5 * (h(i,j+1,k) + h(i,j,k))
4482-
h_arith = 0.5 * (h(i,j+1,k) + h(i,j,k))
4483-
if (e_v(i,K+1) >= D_shallow_v(i)) then
4484-
CS%frhatv(i,J,k) = h_arith
4455+
elseif (CS%hvel_scheme == HYBRID .or. use_default) then
4456+
Z_to_H = GV%Z_to_H ; if (.not.GV%Boussinesq) Z_to_H = GV%RZ_to_H * CS%Rho_BT_lin
4457+
do i=is,ie
4458+
e_v(i,nz+1) = -0.5 * Z_to_H * (G%bathyT(i,j+1) + G%bathyT(i,j))
4459+
D_shallow_v(I) = -Z_to_H * min(G%bathyT(i,j+1), G%bathyT(i,j))
4460+
enddo
4461+
do k=nz,1,-1 ; do i=is,ie
4462+
e_v(i,K) = e_v(i,K+1) + 0.5 * (h(i,j+1,k) + h(i,j,k))
4463+
h_arith = 0.5 * (h(i,j+1,k) + h(i,j,k))
4464+
if (e_v(i,K+1) >= D_shallow_v(i)) then
4465+
hatv(i,k) = h_arith
4466+
else
4467+
h_harm = (h(i,j+1,k) * h(i,j,k)) / (h_arith + h_neglect)
4468+
if (e_v(i,K) <= D_shallow_v(i)) then
4469+
hatv(i,k) = h_harm
44854470
else
4486-
h_harm = (h(i,j+1,k) * h(i,j,k)) / (h_arith + h_neglect)
4487-
if (e_v(i,K) <= D_shallow_v(i)) then
4488-
CS%frhatv(i,J,k) = h_harm
4489-
else
4490-
wt_arith = (e_v(i,K) - D_shallow_v(i)) / (h_arith + h_neglect)
4491-
CS%frhatv(i,J,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
4492-
endif
4471+
wt_arith = (e_v(i,K) - D_shallow_v(i)) / (h_arith + h_neglect)
4472+
hatv(i,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
44934473
endif
4494-
hatvtot(i) = hatvtot(i) + CS%frhatv(i,J,k)
4495-
enddo ; enddo
4496-
elseif (CS%hvel_scheme == HARMONIC) then
4497-
do i=is,ie
4498-
CS%frhatv(i,J,1) = 2.0*(h(i,j+1,1) * h(i,j,1)) / &
4499-
((h(i,j+1,1) + h(i,j,1)) + h_neglect)
4500-
hatvtot(i) = CS%frhatv(i,J,1)
4501-
enddo
4502-
do k=2,nz ; do i=is,ie
4503-
CS%frhatv(i,J,k) = 2.0*(h(i,j+1,k) * h(i,j,k)) / &
4504-
((h(i,j+1,k) + h(i,j,k)) + h_neglect)
4505-
hatvtot(i) = hatvtot(i) + CS%frhatv(i,J,k)
4506-
enddo ; enddo
4507-
endif
4508-
do i=is,ie ; Ihatvtot(i) = G%mask2dCv(i,J) / (hatvtot(i) + h_neglect) ; enddo
4474+
endif
4475+
hatvtot(i) = hatvtot(i) + hatv(i,k)
4476+
enddo ; enddo
4477+
elseif (CS%hvel_scheme == HARMONIC) then
45094478
do k=1,nz ; do i=is,ie
4510-
CS%frhatv(i,J,k) = CS%frhatv(i,J,k) * Ihatvtot(i)
4479+
hatv(i,k) = 2.0*(h(i,j+1,k) * h(i,j,k)) / &
4480+
((h(i,j+1,k) + h(i,j,k)) + h_neglect)
4481+
hatvtot(i) = hatvtot(i) + hatv(i,k)
45114482
enddo ; enddo
45124483
endif
4513-
enddo
45144484

4515-
if (CS%BT_OBC%v_OBCs_on_PE) then
4516-
! Northern boundary conditions
4517-
is_v = max(is, CS%BT_OBC%is_v_N_obc) ; ie_v = min(ie, CS%BT_OBC%ie_v_N_obc)
4518-
Js_v = max(js-1, CS%BT_OBC%Js_v_N_obc) ; Je_v = min(je, CS%BT_OBC%Je_v_N_obc)
4519-
do J=Js_v,Je_v
4520-
do i=is_v,ie_v ; hatvtot(i) = h(i,j,1) ; enddo
4521-
do k=2,nz ; do i=is_v,ie_v
4522-
hatvtot(i) = hatvtot(i) + h(i,j,k)
4523-
enddo ; enddo
4524-
do i=is_v,ie_v
4525-
Ihatvtot(i) = G%mask2dCv(i,J) / (hatvtot(i) + h_neglect)
4526-
enddo
4527-
do k=1,nz ; do i=is_v,ie_v
4528-
if (CS%BT_OBC%v_OBC_type(i,J) > 0) & ! Northern boundary condition
4529-
CS%frhatv(i,J,k) = h(i,j,k) * Ihatvtot(i)
4530-
enddo ; enddo
4531-
enddo
4485+
if (CS%BT_OBC%v_OBCs_on_PE) then
4486+
! Reset v-velocity point thicknesses and their sums at OBC points
4487+
if ((J >= CS%BT_OBC%Js_v_N_obc) .and. (J <= CS%BT_OBC%Je_v_N_obc)) then
4488+
do i = max(is,CS%BT_OBC%is_v_N_obc), min(ie,CS%BT_OBC%ie_v_N_obc)
4489+
if (CS%BT_OBC%v_OBC_type(i,J) > 0) then ! Northern boundary condition
4490+
hatvtot(i) = 0.0
4491+
do k=1,nz
4492+
hatv(i,k) = h(i,j,k)
4493+
hatvtot(i) = hatvtot(i) + hatv(i,k)
4494+
enddo
4495+
endif
4496+
enddo
4497+
endif
4498+
if ((J >= CS%BT_OBC%Js_v_S_obc) .and. (J <= CS%BT_OBC%Je_v_S_obc)) then
4499+
do i = max(is,CS%BT_OBC%is_v_S_obc), min(ie,CS%BT_OBC%ie_v_S_obc)
4500+
if (CS%BT_OBC%v_OBC_type(i,J) < 0) then ! Southern boundary condition
4501+
hatvtot(i) = 0.0
4502+
do k=1,nz
4503+
hatv(i,k) = h(i,j+1,k)
4504+
hatvtot(i) = hatvtot(i) + hatv(i,k)
4505+
enddo
4506+
endif
4507+
enddo
4508+
endif
4509+
endif
45324510

4533-
! Southern boundary conditions
4534-
is_v = max(is, CS%BT_OBC%is_v_S_obc) ; ie_v = min(ie, CS%BT_OBC%ie_v_S_obc)
4535-
Js_v = max(js-1, CS%BT_OBC%Js_v_S_obc) ; Je_v = min(je, CS%BT_OBC%Je_v_S_obc)
4536-
do J=Js_v,Je_v
4537-
do i=is_v,ie_v ; hatvtot(i) = h(i,j+1,1) ; enddo
4538-
do k=2,nz ; do i=is_v,ie_v
4539-
hatvtot(i) = hatvtot(i) + h(i,j+1,k)
4540-
enddo ; enddo
4541-
do i=is_v,ie_v
4542-
Ihatvtot(i) = G%mask2dCv(i,J) / (hatvtot(i) + h_neglect)
4543-
enddo
4544-
do k=1,nz ; do i=is_v,ie_v
4545-
if (CS%BT_OBC%v_OBC_type(i,J) < 0) & ! Southern boundary condition
4546-
CS%frhatv(i,J,k) = h(i,j+1,k) * Ihatvtot(i)
4547-
enddo ; enddo
4548-
enddo
4549-
endif
4511+
! Determine the fractional thickness of each layer at the velocity points.
4512+
do i=is,ie ; Ihatvtot(i) = G%mask2dCv(i,J) / (hatvtot(i) + h_neglect) ; enddo
4513+
do k=1,nz ; do i=is,ie
4514+
CS%frhatv(i,J,k) = hatv(i,k) * Ihatvtot(i)
4515+
enddo ; enddo
4516+
enddo
45504517

45514518
if (CS%debug) then
45524519
call uvchksum("btcalc frhat[uv]", CS%frhatu, CS%frhatv, G%HI, &
@@ -4556,7 +4523,7 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
45564523
call uvchksum("btcalc h_[uv]", h_u, h_v, G%HI, haloshift=0, &
45574524
symmetric=.true., omit_corners=.true., unscale=GV%H_to_MKS, &
45584525
scalar_pair=.true.)
4559-
call hchksum(h, "btcalc h",G%HI, haloshift=1, unscale=GV%H_to_MKS)
4526+
call hchksum(h, "btcalc h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
45604527
endif
45614528

45624529
end subroutine btcalc

0 commit comments

Comments
 (0)