Skip to content

Commit

Permalink
Use the 'dom=' interface to calculate_density
Browse files Browse the repository at this point in the history 
  Use the new variant of calculate_density with the 'dom' argument or no array
extent argument in calls in 15 files.  All answers are bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Apr 17, 2020
1 parent ab9d8e9 commit acf23a4
Show file tree
Hide file tree
Showing 15 changed files with 123 additions and 103 deletions.
6 changes: 2 additions & 4 deletions src/ALE/MOM_regridding.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1889,8 +1889,7 @@ subroutine convective_adjustment(G, GV, h, tv)
do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1

! Compute densities within current water column
call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, &
densities, 1, GV%ke, tv%eqn_of_state )
call calculate_density( tv%T(i,j,:), tv%S(i,j,:), p_col, densities, tv%eqn_of_state )

! Repeat restratification until complete
do
Expand All @@ -1909,8 +1908,7 @@ subroutine convective_adjustment(G, GV, h, tv)
tv%S(i,j,k) = S1 ; tv%S(i,j,k+1) = S0
h(i,j,k) = h1 ; h(i,j,k+1) = h0
! Recompute densities at levels k and k+1
call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), &
densities(k), tv%eqn_of_state )
call calculate_density( tv%T(i,j,k), tv%S(i,j,k), p_col(k), densities(k), tv%eqn_of_state)
call calculate_density( tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), &
densities(k+1), tv%eqn_of_state )
stratified = .false.
Expand Down
10 changes: 5 additions & 5 deletions src/ALE/coord_adapt.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,7 +156,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
0.5 * (tInt(i,j,2:nz) + tInt(i,j-1,2:nz)), &
0.5 * (sInt(i,j,2:nz) + sInt(i,j-1,2:nz)), &
0.5 * (zInt(i,j,2:nz) + zInt(i,j-1,2:nz)) * GV%H_to_Pa, &
alpha, beta, 2, nz - 1, tv%eqn_of_state)
alpha, beta, tv%eqn_of_state, dom=(/2,nz/))

del2sigma(2:nz) = del2sigma(2:nz) + &
(alpha(2:nz) * (tInt(i,j-1,2:nz) - tInt(i,j,2:nz)) + &
Expand All @@ -168,7 +168,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
0.5 * (tInt(i,j,2:nz) + tInt(i,j+1,2:nz)), &
0.5 * (sInt(i,j,2:nz) + sInt(i,j+1,2:nz)), &
0.5 * (zInt(i,j,2:nz) + zInt(i,j+1,2:nz)) * GV%H_to_Pa, &
alpha, beta, 2, nz - 1, tv%eqn_of_state)
alpha, beta, tv%eqn_of_state, dom=(/2,nz/))

del2sigma(2:nz) = del2sigma(2:nz) + &
(alpha(2:nz) * (tInt(i,j+1,2:nz) - tInt(i,j,2:nz)) + &
Expand All @@ -180,7 +180,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
0.5 * (tInt(i,j,2:nz) + tInt(i-1,j,2:nz)), &
0.5 * (sInt(i,j,2:nz) + sInt(i-1,j,2:nz)), &
0.5 * (zInt(i,j,2:nz) + zInt(i-1,j,2:nz)) * GV%H_to_Pa, &
alpha, beta, 2, nz - 1, tv%eqn_of_state)
alpha, beta, tv%eqn_of_state, dom=(/2,nz/))

del2sigma(2:nz) = del2sigma(2:nz) + &
(alpha(2:nz) * (tInt(i-1,j,2:nz) - tInt(i,j,2:nz)) + &
Expand All @@ -192,7 +192,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
0.5 * (tInt(i,j,2:nz) + tInt(i+1,j,2:nz)), &
0.5 * (sInt(i,j,2:nz) + sInt(i+1,j,2:nz)), &
0.5 * (zInt(i,j,2:nz) + zInt(i+1,j,2:nz)) * GV%H_to_Pa, &
alpha, beta, 2, nz - 1, tv%eqn_of_state)
alpha, beta, tv%eqn_of_state, dom=(/2,nz/))

del2sigma(2:nz) = del2sigma(2:nz) + &
(alpha(2:nz) * (tInt(i+1,j,2:nz) - tInt(i,j,2:nz)) + &
Expand All @@ -206,7 +206,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext)
! a positive curvature means we're too light relative to adjacent columns,
! so del2sigma needs to be positive too (push the interface deeper)
call calculate_density_derivs(tInt(i,j,:), sInt(i,j,:), zInt(i,j,:) * GV%H_to_Pa, &
alpha, beta, 1, nz + 1, tv%eqn_of_state)
alpha, beta, tv%eqn_of_state, dom=(/1,nz/)) !### This should be (/1,nz+1/) - see 25 lines below.
do K = 2, nz
! TODO make lower bound here configurable
del2sigma(K) = del2sigma(K) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / &
Expand Down
54 changes: 29 additions & 25 deletions src/core/MOM_PressureForce_Montgomery.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,12 +124,13 @@ 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, start, npts
integer, dimension(2) :: EOSdom ! The computational domain for the equation of state
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
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
EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1)

use_p_atm = .false.
if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
Expand Down Expand Up @@ -228,8 +229,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(:), start, npts, &
tv%eqn_of_state, US=US)
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), &
tv%eqn_of_state, US=US, dom=EOSdom)
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 @@ -245,8 +246,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, start, npts, &
tv%eqn_of_state, US=US)
call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_in_situ, &
tv%eqn_of_state, US=US, dom=EOSdom)
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 @@ -409,13 +410,14 @@ 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, start, npts
integer, dimension(2) :: EOSdom ! The computational domain for the equation of state
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
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
EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1)

use_p_atm = .false.
if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif
Expand Down Expand Up @@ -484,8 +486,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(:), start, npts, &
tv%eqn_of_state, US=US)
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), &
tv%eqn_of_state, US=US, dom=EOSdom)

do k=nkmb+1,nz ; do i=Isq,Ieq+1
if (GV%Rlay(k) < Rho_cv_BL(i)) then
Expand All @@ -505,8 +507,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), start, npts, &
tv%eqn_of_state, US=US)
call calculate_density(tv_tmp%T(:,j,k), tv_tmp%S(:,j,k), p_ref, rho_star(:,j,k), &
tv%eqn_of_state, US=US, dom=EOSdom)
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 @@ -632,10 +634,11 @@ 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, start, npts
integer, dimension(2) :: EOSdom ! The computational domain for the equation of state
integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k

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
EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1)

Rho0xG = Rho0 * GV%g_Earth
G_Rho0 = GV%g_Earth / GV%Rho0
Expand All @@ -662,8 +665,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, start, npts, &
tv%eqn_of_state, US=US)
call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, &
tv%eqn_of_state, US=US, dom=EOSdom)
do i=Isq,Ieq+1
pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i)) * GV%H_to_Z
enddo
Expand All @@ -673,8 +676,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, start, npts, &
tv%eqn_of_state, US=US)
call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, &
tv%eqn_of_state, US=US, dom=EOSdom)
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 @@ -733,10 +736,11 @@ 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, start, npts
integer, dimension(2) :: EOSdom ! The computational domain for the equation of state
integer :: Isq, Ieq, Jsq, Jeq, nz, i, j, k

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
EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1)

use_EOS = associated(tv%eqn_of_state)

Expand All @@ -759,8 +763,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, start, npts, &
tv%eqn_of_state, US=US)
call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, &
tv%eqn_of_state, US=US, dom=EOSdom)
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 @@ -770,9 +774,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, 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)
call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, tv%eqn_of_state, US=US, dom=EOSdom)
call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, &
tv%eqn_of_state, US=US, dom=EOSdom)
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
30 changes: 16 additions & 14 deletions src/core/MOM_PressureForce_analytic_FV.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,13 +176,14 @@ 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, start, npts
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
integer, dimension(2) :: EOSdom
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
EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1)

if (.not.associated(CS)) call MOM_error(FATAL, &
"MOM_PressureForce_AFV_nonBouss: Module must be initialized before it is used.")
Expand Down Expand Up @@ -228,8 +229,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(:), start, npts, &
tv%eqn_of_state, US=US)
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), &
tv%eqn_of_state, US=US, dom=EOSdom)
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 @@ -334,8 +335,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, start, npts, &
tv%eqn_of_state, US=US)
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, &
tv%eqn_of_state, US=US, dom=EOSdom)

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 @@ -504,13 +505,14 @@ 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, start, npts
integer, dimension(2) :: EOSdom ! The computational domain for the equation of state
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
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
EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1)

if (.not.associated(CS)) call MOM_error(FATAL, &
"MOM_PressureForce_AFV_Bouss: Module must be initialized before it is used.")
Expand Down Expand Up @@ -578,8 +580,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(:), start, npts, &
tv%eqn_of_state, US=US)
call calculate_density(tv%T(:,j,nkmb), tv%S(:,j,nkmb), p_ref, Rho_cv_BL(:), &
tv%eqn_of_state, US=US, dom=EOSdom)

do k=nkmb+1,nz ; do i=Isq,Ieq+1
if (GV%Rlay(k) < Rho_cv_BL(i)) then
Expand All @@ -601,11 +603,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, start, npts, &
tv%eqn_of_state, US=US)
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, &
tv%eqn_of_state, US=US, dom=EOSdom)
else
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)
call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, &
tv%eqn_of_state, US=US, dom=EOSdom)
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
Loading

0 comments on commit acf23a4

Please sign in to comment.