Skip to content

Commit 98feea2

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 bed5ba7 commit 98feea2

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
@@ -4210,11 +4210,8 @@ subroutine destroy_BT_OBC(BT_OBC)
42104210

42114211
end subroutine destroy_BT_OBC
42124212

4213-
!> btcalc calculates the barotropic velocities from the full velocity and
4214-
!! thickness fields, determines the fraction of the total water column in each
4215-
!! layer at velocity points, and determines a corrective fictitious mass source
4216-
!! that will drive the barotropic estimate of the free surface height toward the
4217-
!! baroclinic estimate.
4213+
!> btcalc determines the fraction of the total water column in each
4214+
!! layer at velocity points.
42184215
subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
42194216
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
42204217
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
@@ -4244,6 +4241,8 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
42444241
type(ocean_OBC_type), optional, pointer :: OBC !< Open boundary control structure.
42454242

42464243
! Local variables
4244+
real :: hatu(SZIB_(G),SZK_(GV)) ! The layer thicknesses interpolated to u points [H ~> m or kg m-2]
4245+
real :: hatv(SZI_(G),SZK_(GV)) ! The layer thicknesses interpolated to v points [H ~> m or kg m-2]
42474246
real :: hatutot(SZIB_(G)) ! The sum of the layer thicknesses interpolated to u points [H ~> m or kg m-2].
42484247
real :: hatvtot(SZI_(G)) ! The sum of the layer thicknesses interpolated to v points [H ~> m or kg m-2].
42494248
real :: Ihatutot(SZIB_(G)) ! Ihatutot is the inverse of hatutot [H-1 ~> m-1 or m2 kg-1].
@@ -4261,15 +4260,11 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
42614260
real :: D_shallow_v(SZIB_(G))! The height of the shallower of the adjacent bathymetric depths
42624261
! around a v-point (positive upward) [H ~> m or kg m-2]
42634262
real :: Z_to_H ! A local conversion factor [H Z-1 ~> nondim or kg m-3]
4264-
real :: htot ! The sum of the layer thicknesses [H ~> m or kg m-2].
4265-
real :: Ihtot ! The inverse of htot [H-1 ~> m-1 or m2 kg-1].
42664263

42674264
logical :: use_default, test_dflt
42684265
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, i, j, k
42694266
integer :: is_v, ie_v, Js_v, Je_v
42704267

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

@@ -4293,191 +4288,163 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
42934288
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
42944289
h_neglect = GV%H_subroundoff
42954290

4296-
! This estimates the fractional thickness of each layer at the velocity
4297-
! points, using a harmonic mean estimate.
42984291

42994292
!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_u,CS,h_neglect,h,use_default,G,GV) &
4300-
!$OMP private(hatutot,Ihatutot,e_u,D_shallow_u,h_arith,h_harm,wt_arith,Z_to_H,htot,Ihtot)
4293+
!$OMP private(hatu,hatutot,Ihatutot,e_u,D_shallow_u,h_arith,h_harm,wt_arith,Z_to_H)
43014294
do j=js,je
4295+
do I=is-1,ie ; hatutot(I) = 0.0 ; enddo
4296+
43024297
if (present(h_u)) then
4303-
do I=is-1,ie ; hatutot(I) = h_u(I,j,1) ; enddo
4304-
do k=2,nz ; do I=is-1,ie
4305-
hatutot(I) = hatutot(I) + h_u(I,j,k)
4298+
do k=1,nz ; do I=is-1,ie
4299+
hatu(I,k) = h_u(I,j,k)
4300+
hatutot(I) = hatutot(I) + hatu(I,k)
43064301
enddo ; enddo
4307-
do I=is-1,ie ; Ihatutot(I) = G%mask2dCu(I,j) / (hatutot(I) + h_neglect) ; enddo
4302+
elseif (CS%hvel_scheme == ARITHMETIC) then
43084303
do k=1,nz ; do I=is-1,ie
4309-
CS%frhatu(I,j,k) = h_u(I,j,k) * Ihatutot(I)
4304+
hatu(I,k) = 0.5 * (h(i+1,j,k) + h(i,j,k))
4305+
hatutot(I) = hatutot(I) + hatu(I,k)
43104306
enddo ; enddo
4311-
else
4312-
if (CS%hvel_scheme == ARITHMETIC) then
4313-
do I=is-1,ie
4314-
CS%frhatu(I,j,1) = 0.5 * (h(i+1,j,1) + h(i,j,1))
4315-
hatutot(I) = CS%frhatu(I,j,1)
4316-
enddo
4317-
do k=2,nz ; do I=is-1,ie
4318-
CS%frhatu(I,j,k) = 0.5 * (h(i+1,j,k) + h(i,j,k))
4319-
hatutot(I) = hatutot(I) + CS%frhatu(I,j,k)
4320-
enddo ; enddo
4321-
elseif (CS%hvel_scheme == HYBRID .or. use_default) then
4322-
Z_to_H = GV%Z_to_H ; if (.not.GV%Boussinesq) Z_to_H = GV%RZ_to_H * CS%Rho_BT_lin
4323-
do I=is-1,ie
4324-
e_u(I,nz+1) = -0.5 * Z_to_H * (G%bathyT(i+1,j) + G%bathyT(i,j))
4325-
D_shallow_u(I) = -Z_to_H * min(G%bathyT(i+1,j), G%bathyT(i,j))
4326-
hatutot(I) = 0.0
4327-
enddo
4328-
do k=nz,1,-1 ; do I=is-1,ie
4329-
e_u(I,K) = e_u(I,K+1) + 0.5 * (h(i+1,j,k) + h(i,j,k))
4330-
h_arith = 0.5 * (h(i+1,j,k) + h(i,j,k))
4331-
if (e_u(I,K+1) >= D_shallow_u(I)) then
4332-
CS%frhatu(I,j,k) = h_arith
4307+
elseif (CS%hvel_scheme == HYBRID .or. use_default) then
4308+
Z_to_H = GV%Z_to_H ; if (.not.GV%Boussinesq) Z_to_H = GV%RZ_to_H * CS%Rho_BT_lin
4309+
do I=is-1,ie
4310+
e_u(I,nz+1) = -0.5 * Z_to_H * (G%bathyT(i+1,j) + G%bathyT(i,j))
4311+
D_shallow_u(I) = -Z_to_H * min(G%bathyT(i+1,j), G%bathyT(i,j))
4312+
enddo
4313+
do k=nz,1,-1 ; do I=is-1,ie
4314+
e_u(I,K) = e_u(I,K+1) + 0.5 * (h(i+1,j,k) + h(i,j,k))
4315+
h_arith = 0.5 * (h(i+1,j,k) + h(i,j,k))
4316+
if (e_u(I,K+1) >= D_shallow_u(I)) then
4317+
hatu(I,k) = h_arith
4318+
else
4319+
h_harm = (h(i+1,j,k) * h(i,j,k)) / (h_arith + h_neglect)
4320+
if (e_u(I,K) <= D_shallow_u(I)) then
4321+
hatu(I,k) = h_harm
43334322
else
4334-
h_harm = (h(i+1,j,k) * h(i,j,k)) / (h_arith + h_neglect)
4335-
if (e_u(I,K) <= D_shallow_u(I)) then
4336-
CS%frhatu(I,j,k) = h_harm
4337-
else
4338-
wt_arith = (e_u(I,K) - D_shallow_u(I)) / (h_arith + h_neglect)
4339-
CS%frhatu(I,j,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
4340-
endif
4323+
wt_arith = (e_u(I,K) - D_shallow_u(I)) / (h_arith + h_neglect)
4324+
hatu(I,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
43414325
endif
4342-
hatutot(I) = hatutot(I) + CS%frhatu(I,j,k)
4343-
enddo ; enddo
4344-
elseif (CS%hvel_scheme == HARMONIC) then
4345-
do I=is-1,ie
4346-
CS%frhatu(I,j,1) = 2.0*(h(i+1,j,1) * h(i,j,1)) / &
4347-
((h(i+1,j,1) + h(i,j,1)) + h_neglect)
4348-
hatutot(I) = CS%frhatu(I,j,1)
4349-
enddo
4350-
do k=2,nz ; do I=is-1,ie
4351-
CS%frhatu(I,j,k) = 2.0*(h(i+1,j,k) * h(i,j,k)) / &
4352-
((h(i+1,j,k) + h(i,j,k)) + h_neglect)
4353-
hatutot(I) = hatutot(I) + CS%frhatu(I,j,k)
4354-
enddo ; enddo
4355-
endif
4356-
do I=is-1,ie ; Ihatutot(I) = G%mask2dCu(I,j) / (hatutot(I) + h_neglect) ; enddo
4326+
endif
4327+
hatutot(I) = hatutot(I) + hatu(I,k)
4328+
enddo ; enddo
4329+
elseif (CS%hvel_scheme == HARMONIC) then
4330+
! Interpolates thicknesses onto u grid points with the
4331+
! second order accurate estimate h = 2*(h+ * h-)/(h+ + h-).
43574332
do k=1,nz ; do I=is-1,ie
4358-
CS%frhatu(I,j,k) = CS%frhatu(I,j,k) * Ihatutot(I)
4333+
hatu(I,k) = 2.0*(h(i+1,j,k) * h(i,j,k)) / &
4334+
((h(i+1,j,k) + h(i,j,k)) + h_neglect)
4335+
hatutot(I) = hatutot(I) + hatu(I,k)
43594336
enddo ; enddo
43604337
endif
4338+
43614339
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
4340+
! Reset velocity point thicknesses and their sums at OBC points
4341+
if ((j >= CS%BT_OBC%js_u_E_obc) .and. (j <= CS%BT_OBC%je_u_E_obc)) then
4342+
do I = max(is-1,CS%BT_OBC%Is_u_E_obc), min(ie,CS%BT_OBC%Ie_u_E_obc)
43654343
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
4344+
hatutot(I) = 0.0
4345+
do k=1,nz
4346+
hatu(I,k) = h(i,j,k)
4347+
hatutot(I) = hatutot(I) + hatu(I,k)
4348+
enddo
43704349
endif
4350+
enddo
4351+
endif
4352+
if ((j >= CS%BT_OBC%js_u_W_obc) .and. (j <= CS%BT_OBC%je_u_W_obc)) then
4353+
do I = max(is-1,CS%BT_OBC%Is_u_W_obc), min(ie,CS%BT_OBC%Ie_u_W_obc)
43714354
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
4355+
hatutot(I) = 0.0
4356+
do k=1,nz
4357+
hatu(I,k) = h(i+1,j,k)
4358+
hatutot(I) = hatutot(I) + hatu(I,k)
4359+
enddo
43764360
endif
43774361
enddo
43784362
endif
43794363
endif
4364+
4365+
! Determine the fractional thickness of each layer at the velocity points.
4366+
do I=is-1,ie ; Ihatutot(I) = G%mask2dCu(I,j) / (hatutot(I) + h_neglect) ; enddo
4367+
do k=1,nz ; do I=is-1,ie
4368+
CS%frhatu(I,j,k) = hatu(I,k) * Ihatutot(I)
4369+
enddo ; enddo
43804370
enddo
43814371

43824372
!$OMP parallel do default(none) shared(is,ie,js,je,nz,CS,G,GV,h_v,h_neglect,h,use_default) &
4383-
!$OMP private(hatvtot,Ihatvtot,e_v,D_shallow_v,h_arith,h_harm,wt_arith,Z_to_H)
4373+
!$OMP private(hatv,hatvtot,Ihatvtot,e_v,D_shallow_v,h_arith,h_harm,wt_arith,Z_to_H)
43844374
do J=js-1,je
4375+
do i=is,ie ; hatvtot(i) = 0.0 ; enddo
43854376
if (present(h_v)) then
4386-
do i=is,ie ; hatvtot(i) = h_v(i,J,1) ; enddo
4387-
do k=2,nz ; do i=is,ie
4388-
hatvtot(i) = hatvtot(i) + h_v(i,J,k)
4377+
do k=1,nz ; do i=is,ie
4378+
hatv(i,k) = h_v(i,J,k)
4379+
hatvtot(i) = hatvtot(i) + hatv(i,k)
43894380
enddo ; enddo
4390-
do i=is,ie ; Ihatvtot(i) = G%mask2dCv(i,J) / (hatvtot(i) + h_neglect) ; enddo
4381+
elseif (CS%hvel_scheme == ARITHMETIC) then
43914382
do k=1,nz ; do i=is,ie
4392-
CS%frhatv(i,J,k) = h_v(i,J,k) * Ihatvtot(i)
4383+
hatv(i,k) = 0.5 * (h(i,j+1,k) + h(i,j,k))
4384+
hatvtot(i) = hatvtot(i) + hatv(i,k)
43934385
enddo ; enddo
4394-
else
4395-
if (CS%hvel_scheme == ARITHMETIC) then
4396-
do i=is,ie
4397-
CS%frhatv(i,J,1) = 0.5 * (h(i,j+1,1) + h(i,j,1))
4398-
hatvtot(i) = CS%frhatv(i,J,1)
4399-
enddo
4400-
do k=2,nz ; do i=is,ie
4401-
CS%frhatv(i,J,k) = 0.5 * (h(i,j+1,k) + h(i,j,k))
4402-
hatvtot(i) = hatvtot(i) + CS%frhatv(i,J,k)
4403-
enddo ; enddo
4404-
elseif (CS%hvel_scheme == HYBRID .or. use_default) then
4405-
Z_to_H = GV%Z_to_H ; if (.not.GV%Boussinesq) Z_to_H = GV%RZ_to_H * CS%Rho_BT_lin
4406-
do i=is,ie
4407-
e_v(i,nz+1) = -0.5 * Z_to_H * (G%bathyT(i,j+1) + G%bathyT(i,j))
4408-
D_shallow_v(I) = -Z_to_H * min(G%bathyT(i,j+1), G%bathyT(i,j))
4409-
hatvtot(I) = 0.0
4410-
enddo
4411-
do k=nz,1,-1 ; do i=is,ie
4412-
e_v(i,K) = e_v(i,K+1) + 0.5 * (h(i,j+1,k) + h(i,j,k))
4413-
h_arith = 0.5 * (h(i,j+1,k) + h(i,j,k))
4414-
if (e_v(i,K+1) >= D_shallow_v(i)) then
4415-
CS%frhatv(i,J,k) = h_arith
4386+
elseif (CS%hvel_scheme == HYBRID .or. use_default) then
4387+
Z_to_H = GV%Z_to_H ; if (.not.GV%Boussinesq) Z_to_H = GV%RZ_to_H * CS%Rho_BT_lin
4388+
do i=is,ie
4389+
e_v(i,nz+1) = -0.5 * Z_to_H * (G%bathyT(i,j+1) + G%bathyT(i,j))
4390+
D_shallow_v(I) = -Z_to_H * min(G%bathyT(i,j+1), G%bathyT(i,j))
4391+
enddo
4392+
do k=nz,1,-1 ; do i=is,ie
4393+
e_v(i,K) = e_v(i,K+1) + 0.5 * (h(i,j+1,k) + h(i,j,k))
4394+
h_arith = 0.5 * (h(i,j+1,k) + h(i,j,k))
4395+
if (e_v(i,K+1) >= D_shallow_v(i)) then
4396+
hatv(i,k) = h_arith
4397+
else
4398+
h_harm = (h(i,j+1,k) * h(i,j,k)) / (h_arith + h_neglect)
4399+
if (e_v(i,K) <= D_shallow_v(i)) then
4400+
hatv(i,k) = h_harm
44164401
else
4417-
h_harm = (h(i,j+1,k) * h(i,j,k)) / (h_arith + h_neglect)
4418-
if (e_v(i,K) <= D_shallow_v(i)) then
4419-
CS%frhatv(i,J,k) = h_harm
4420-
else
4421-
wt_arith = (e_v(i,K) - D_shallow_v(i)) / (h_arith + h_neglect)
4422-
CS%frhatv(i,J,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
4423-
endif
4402+
wt_arith = (e_v(i,K) - D_shallow_v(i)) / (h_arith + h_neglect)
4403+
hatv(i,k) = wt_arith*h_arith + (1.0-wt_arith)*h_harm
44244404
endif
4425-
hatvtot(i) = hatvtot(i) + CS%frhatv(i,J,k)
4426-
enddo ; enddo
4427-
elseif (CS%hvel_scheme == HARMONIC) then
4428-
do i=is,ie
4429-
CS%frhatv(i,J,1) = 2.0*(h(i,j+1,1) * h(i,j,1)) / &
4430-
((h(i,j+1,1) + h(i,j,1)) + h_neglect)
4431-
hatvtot(i) = CS%frhatv(i,J,1)
4432-
enddo
4433-
do k=2,nz ; do i=is,ie
4434-
CS%frhatv(i,J,k) = 2.0*(h(i,j+1,k) * h(i,j,k)) / &
4435-
((h(i,j+1,k) + h(i,j,k)) + h_neglect)
4436-
hatvtot(i) = hatvtot(i) + CS%frhatv(i,J,k)
4437-
enddo ; enddo
4438-
endif
4439-
do i=is,ie ; Ihatvtot(i) = G%mask2dCv(i,J) / (hatvtot(i) + h_neglect) ; enddo
4405+
endif
4406+
hatvtot(i) = hatvtot(i) + hatv(i,k)
4407+
enddo ; enddo
4408+
elseif (CS%hvel_scheme == HARMONIC) then
44404409
do k=1,nz ; do i=is,ie
4441-
CS%frhatv(i,J,k) = CS%frhatv(i,J,k) * Ihatvtot(i)
4410+
hatv(i,k) = 2.0*(h(i,j+1,k) * h(i,j,k)) / &
4411+
((h(i,j+1,k) + h(i,j,k)) + h_neglect)
4412+
hatvtot(i) = hatvtot(i) + hatv(i,k)
44424413
enddo ; enddo
44434414
endif
4444-
enddo
44454415

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
4460-
CS%frhatv(i,J,k) = h(i,j,k) * Ihatvtot(i)
4461-
enddo ; enddo
4462-
enddo
4416+
if (CS%BT_OBC%v_OBCs_on_PE) then
4417+
! Reset v-velocity point thicknesses and their sums at OBC points
4418+
if ((J >= CS%BT_OBC%Js_v_N_obc) .and. (J <= CS%BT_OBC%Je_v_N_obc)) then
4419+
do i = max(is,CS%BT_OBC%is_v_N_obc), min(ie,CS%BT_OBC%ie_v_N_obc)
4420+
if (CS%BT_OBC%v_OBC_type(i,J) > 0) then ! Northern boundary condition
4421+
hatvtot(i) = 0.0
4422+
do k=1,nz
4423+
hatv(i,k) = h(i,j,k)
4424+
hatvtot(i) = hatvtot(i) + hatv(i,k)
4425+
enddo
4426+
endif
4427+
enddo
4428+
endif
4429+
if ((J >= CS%BT_OBC%Js_v_S_obc) .and. (J <= CS%BT_OBC%Je_v_S_obc)) then
4430+
do i = max(is,CS%BT_OBC%is_v_S_obc), min(ie,CS%BT_OBC%ie_v_S_obc)
4431+
if (CS%BT_OBC%v_OBC_type(i,J) < 0) then ! Southern boundary condition
4432+
hatvtot(i) = 0.0
4433+
do k=1,nz
4434+
hatv(i,k) = h(i,j+1,k)
4435+
hatvtot(i) = hatvtot(i) + hatv(i,k)
4436+
enddo
4437+
endif
4438+
enddo
4439+
endif
4440+
endif
44634441

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
4477-
CS%frhatv(i,J,k) = h(i,j+1,k) * Ihatvtot(i)
4478-
enddo ; enddo
4479-
enddo
4480-
endif
4442+
! Determine the fractional thickness of each layer at the velocity points.
4443+
do i=is,ie ; Ihatvtot(i) = G%mask2dCv(i,J) / (hatvtot(i) + h_neglect) ; enddo
4444+
do k=1,nz ; do i=is,ie
4445+
CS%frhatv(i,J,k) = hatv(i,k) * Ihatvtot(i)
4446+
enddo ; enddo
4447+
enddo
44814448

44824449
if (CS%debug) then
44834450
call uvchksum("btcalc frhat[uv]", CS%frhatu, CS%frhatv, G%HI, &
@@ -4487,7 +4454,7 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
44874454
call uvchksum("btcalc h_[uv]", h_u, h_v, G%HI, haloshift=0, &
44884455
symmetric=.true., omit_corners=.true., unscale=GV%H_to_MKS, &
44894456
scalar_pair=.true.)
4490-
call hchksum(h, "btcalc h",G%HI, haloshift=1, unscale=GV%H_to_MKS)
4457+
call hchksum(h, "btcalc h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
44914458
endif
44924459

44934460
end subroutine btcalc

0 commit comments

Comments
 (0)