Skip to content

Commit

Permalink
(*)Set loop bounds in calculate_density calls
Browse files Browse the repository at this point in the history 
  Explicitly set loop bounds in some calculate_density calls with halos that are
only set around velocity points to avoid errors with non-symmetric memory.  All
answers are bitwise identical in test cases, and this should fix a problem that
was detected by the automated testing.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Apr 15, 2020
1 parent 87054b4 commit 77b6b74
Show file tree
Hide file tree
Showing 4 changed files with 57 additions and 48 deletions.
46 changes: 25 additions & 21 deletions src/core/MOM_PressureForce_Montgomery.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,11 +124,12 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
real :: alpha_Lay(SZK_(G)) ! The specific volume of each layer [R-1 ~> m3 kg-1].
real :: dalpha_int(SZK_(G)+1) ! The change in specific volume across each
! interface [R-1 ~> m3 kg-1].
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts
integer :: i, j, k
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
nkmb=GV%nk_rho_varies
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2

use_p_atm = .false.
if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
Expand Down Expand Up @@ -227,8 +228,8 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
do k=1,nkmb ; do i=Isq,Ieq+1
tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
enddo ; enddo
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, &
tv%eqn_of_state, US=US)
do k=nkmb+1,nz ; do i=Isq,Ieq+1
if (GV%Rlay(k) < Rho_cv_BL(i)) then
tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
Expand All @@ -244,8 +245,8 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
endif
!$OMP parallel do default(shared) private(rho_in_situ)
do k=1,nz ; do j=Jsq,Jeq+1
call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_in_situ, start, npts, &
tv%eqn_of_state, US=US)
do i=Isq,Ieq+1 ; alpha_star(i,j,k) = 1.0 / rho_in_situ(i) ; enddo
enddo ; enddo
endif ! use_EOS
Expand Down Expand Up @@ -408,12 +409,13 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
! gradient terms are to be split into
! barotropic and baroclinic pieces.
type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts
integer :: i, j, k

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
nkmb=GV%nk_rho_varies
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2

use_p_atm = .false.
if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
Expand Down Expand Up @@ -482,8 +484,8 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
do k=1,nkmb ; do i=Isq,Ieq+1
tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
enddo ; enddo
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, &
tv%eqn_of_state, US=US)

do k=nkmb+1,nz ; do i=Isq,Ieq+1
if (GV%Rlay(k) < Rho_cv_BL(i)) then
Expand All @@ -503,8 +505,8 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
! will come down with a fatal error if there is any compressibility.
!$OMP parallel do default(shared)
do k=1,nz+1 ; do j=Jsq,Jeq+1
call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), start, npts, &
tv%eqn_of_state, US=US)
do i=Isq,Ieq+1 ; rho_star(i,j,k) = G_Rho0*rho_star(i,j,k) ; enddo
enddo ; enddo
endif ! use_EOS
Expand Down Expand Up @@ -630,9 +632,10 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
! an equation of state.
real :: z_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [Z ~> m].
integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k
integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k, start, npts

Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke
start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2

Rho0xG = Rho0 * GV%g_Earth
G_Rho0 = GV%g_Earth / GV%Rho0
Expand All @@ -659,8 +662,8 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
Ihtot(i) = GV%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect)
press(i) = -Rho0xG*e(i,j,1)
enddo
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, start, npts, &
tv%eqn_of_state, US=US)
do i=Isq,Ieq+1
pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i)) * GV%H_to_Z
enddo
Expand All @@ -670,8 +673,8 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
T_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
S_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
enddo
call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, start, npts, &
tv%eqn_of_state, US=US)
do i=Isq,Ieq+1
pbce(i,j,k) = pbce(i,j,k-1) + G_Rho0 * &
((e(i,j,K) - e(i,j,nz+1)) * Ihtot(i)) * &
Expand Down Expand Up @@ -730,9 +733,10 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
real :: dp_neglect ! A thickness that is so small it is usually lost
! in roundoff and can be neglected [R L2 T-2 ~> Pa].
logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k
integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k, start, npts

Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke
start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2

use_EOS = associated(tv%eqn_of_state)

Expand All @@ -755,8 +759,8 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
else
!$OMP parallel do default(shared) private(T_int,S_int,dR_dT,dR_dS,rho_in_situ)
do j=Jsq,Jeq+1
call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, start, npts, &
tv%eqn_of_state, US=US)
do i=Isq,Ieq+1
C_htot(i,j) = dP_dH / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect)
pbce(i,j,nz) = dP_dH / (rho_in_situ(i))
Expand All @@ -766,9 +770,9 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star)
T_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1))
S_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1))
enddo
call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, G%HI, tv%eqn_of_state, US, halo=1)
call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, start, npts, tv%eqn_of_state, US=US)
call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, start, npts, &
tv%eqn_of_state, US=US)
do i=Isq,Ieq+1
pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,K+1)-p(i,j,1))*C_htot(i,j)) * &
((dR_dT(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + &
Expand Down
26 changes: 14 additions & 12 deletions src/core/MOM_PressureForce_analytic_FV.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,12 +176,13 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p
real :: H_to_RL2_T2 ! A factor to convert from thicknesss units (H) to pressure units [R L2 T-2 H-1 ~> Pa H-1].
! real :: oneatm = 101325.0 ! 1 atm in [Pa] = [kg m-1 s-2]
real, parameter :: C1_6 = 1.0/6.0
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts
integer :: i, j, k

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
nkmb=GV%nk_rho_varies
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2

if (.not.associated(CS)) call MOM_error(FATAL, &
"MOM_PressureForce_AFV_nonBouss: Module must be initialized before it is used.")
Expand Down Expand Up @@ -227,8 +228,8 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p
do k=1,nkmb ; do i=Isq,Ieq+1
tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
enddo ; enddo
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, &
tv%eqn_of_state, US=US)
do k=nkmb+1,nz ; do i=Isq,Ieq+1
if (GV%Rlay(k) < Rho_cv_BL(i)) then
tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
Expand Down Expand Up @@ -333,8 +334,8 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p
if (use_EOS) then
!$OMP parallel do default(shared) private(rho_in_situ)
do j=Jsq,Jeq+1
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, start, npts, &
tv%eqn_of_state, US=US)

do i=Isq,Ieq+1
dM(i,j) = (CS%GFS_scale - 1.0) * (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j))
Expand Down Expand Up @@ -503,12 +504,13 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.

real, parameter :: C1_6 = 1.0/6.0
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts
integer :: i, j, k

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
nkmb=GV%nk_rho_varies
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2

if (.not.associated(CS)) call MOM_error(FATAL, &
"MOM_PressureForce_AFV_Bouss: Module must be initialized before it is used.")
Expand Down Expand Up @@ -576,8 +578,8 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
do k=1,nkmb ; do i=Isq,Ieq+1
tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
enddo ; enddo
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, &
tv%eqn_of_state, US=US)

do k=nkmb+1,nz ; do i=Isq,Ieq+1
if (GV%Rlay(k) < Rho_cv_BL(i)) then
Expand All @@ -599,11 +601,11 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1
if (use_p_atm) then
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, start, npts, &
tv%eqn_of_state, US=US)
else
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, start, npts, &
tv%eqn_of_state, US=US)
endif
do i=Isq,Ieq+1
dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * e(i,j,1)
Expand Down
26 changes: 14 additions & 12 deletions src/core/MOM_PressureForce_blocked_AFV.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,13 +173,14 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm,
real :: H_to_RL2_T2 ! A factor to convert from thicknesss units (H) to pressure units [R L2 T-2 H-1 ~> Pa H-1].
! real :: oneatm = 101325.0 ! 1 atm in [Pa] = [kg m-1 s-2]
real, parameter :: C1_6 = 1.0/6.0
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts
integer :: is_bk, ie_bk, js_bk, je_bk, Isq_bk, Ieq_bk, Jsq_bk, Jeq_bk
integer :: i, j, k, n, ib, jb, ioff_bk, joff_bk

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
nkmb=GV%nk_rho_varies
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2

if (.not.associated(CS)) call MOM_error(FATAL, &
"MOM_PressureForce_AFV_nonBouss: Module must be initialized before it is used.")
Expand Down Expand Up @@ -223,8 +224,8 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm,
do k=1,nkmb ; do i=Isq,Ieq+1
tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
enddo ; enddo
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, &
tv%eqn_of_state, US=US)
do k=nkmb+1,nz ; do i=Isq,Ieq+1
if (GV%Rlay(k) < Rho_cv_BL(i)) then
tv_tmp%T(i,j,k) = tv%T(i,j,nkmb) ; tv_tmp%S(i,j,k) = tv%S(i,j,nkmb)
Expand Down Expand Up @@ -301,8 +302,8 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm,
if (use_EOS) then
!$OMP parallel do default(shared) private(rho_in_situ)
do j=Jsq,Jeq+1
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, start, npts, &
tv%eqn_of_state, US=US)

do i=Isq,Ieq+1
dM(i,j) = (CS%GFS_scale - 1.0) * (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j))
Expand Down Expand Up @@ -488,14 +489,15 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp,
type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.

real, parameter :: C1_6 = 1.0/6.0
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, start, npts
integer :: is_bk, ie_bk, js_bk, je_bk, Isq_bk, Ieq_bk, Jsq_bk, Jeq_bk
integer :: ioff_bk, joff_bk
integer :: i, j, k, n, ib, jb

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
nkmb=GV%nk_rho_varies
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2

if (.not.associated(CS)) call MOM_error(FATAL, &
"MOM_PressureForce_AFV_Bouss: Module must be initialized before it is used.")
Expand Down Expand Up @@ -563,8 +565,8 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp,
do k=1,nkmb ; do i=Isq,Ieq+1
tv_tmp%T(i,j,k) = tv%T(i,j,k) ; tv_tmp%S(i,j,k) = tv%S(i,j,k)
enddo ; enddo
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), start, npts, &
tv%eqn_of_state, US=US)

do k=nkmb+1,nz ; do i=Isq,Ieq+1
if (GV%Rlay(k) < Rho_cv_BL(i)) then
Expand All @@ -586,11 +588,11 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp,
!$OMP parallel do default(shared)
do j=Jsq,Jeq+1
if (use_p_atm) then
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, start, npts, &
tv%eqn_of_state, US=US)
else
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, start, npts, &
tv%eqn_of_state, US=US)
endif
do i=Isq,Ieq+1
dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * e(i,j,1)
Expand Down
7 changes: 4 additions & 3 deletions src/parameterizations/vertical/MOM_set_viscosity.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -273,12 +273,13 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize)
! accuracy of a single L(:) Newton iteration
logical :: use_L0, do_one_L_iter ! Control flags for L(:) Newton iteration
logical :: use_BBL_EOS, do_i(SZIB_(G))
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m, n, K2, nkmb, nkml
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m, n, K2, nkmb, nkml, start, npts
integer :: itt, maxitt=20
type(ocean_OBC_type), pointer :: OBC => NULL()

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
start = Isq - (G%isd-1) ; npts = G%iec - Isq + 2
nkmb = GV%nk_rho_varies ; nkml = GV%nkml
h_neglect = GV%H_subroundoff
Rho0x400_G = 400.0*(GV%Rho0 / (US%L_to_Z**2 * GV%g_Earth)) * GV%Z_to_H
Expand Down Expand Up @@ -315,8 +316,8 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, symmetrize)
do i=Isq,Ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
!$OMP parallel do default(shared)
do k=1,nkmb ; do j=Jsq,Jeq+1
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, Rml(:,j,k), G%HI, &
tv%eqn_of_state, US, halo=1)
call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, Rml(:,j,k), start, npts, &
tv%eqn_of_state, US=US)
enddo ; enddo
endif

Expand Down

0 comments on commit 77b6b74

Please sign in to comment.