FDS Source: Pad arrays for RHO, ZZ and TMP

Source/ccib.f90

@@ -56,7 +56,8 @@ MODULE CC_SCALARS
 REAL(EB), ALLOCATABLE, DIMENSION(:)  :: F_Z, RZ_Z, RZ_ZS, P_0_CV, TMP_0_CV, RHO_0_CV, ZCEN_CV
 REAL(EB), ALLOCATABLE, DIMENSION(:,:):: F_Z0, RZ_Z0
 
-REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,Z_TEMP,F_TEMP
+REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP
+REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP
 
 ! Types of interpolation:
 INTEGER, PARAMETER :: CC_LINEAR_INTERPOLATION    = 1
@@ -11252,7 +11253,7 @@ SUBROUTINE CC_DENSITY(T,DT)
       ! Store RHO*ZZ values at step n:
       IF (.NOT.ALLOCATED(MESHES(NM)%RHO_ZZN)) ALLOCATE(MESHES(NM)%RHO_ZZN(0:IBP1,0:JBP1,0:KBP1,N_TOTAL_SCALARS))
       DO N=1,N_TOTAL_SCALARS
-         MESHES(NM)%RHO_ZZN(:,:,:,N) = RHO(:,:,:)*ZZ(:,:,:,N)
+         MESHES(NM)%RHO_ZZN(0:IBP1,0:JBP1,0:KBP1,N) = RHO(0:IBP1,0:JBP1,0:KBP1)*ZZ(0:IBP1,0:JBP1,0:KBP1,N)
       ENDDO
 
       ! Second cut-cells, these variables being filled are only used for exporting to slices and applying Boundary

Source/divg.f90

@@ -155,7 +155,15 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM)
 
       ! Store max diffusivity for stability check
 
-      IF (CHECK_VN) D_Z_MAX = MAX(D_Z_MAX,RHO_D/(RHOP+TWO_EPSILON_EB))
+      IF (CHECK_VN) THEN
+         DO K=0,KBP1
+            DO J=0,JBP1
+               DO I=0,IBP1
+                  D_Z_MAX(I,J,K) = MAX(D_Z_MAX(I,J,K),RHO_D(I,J,K)/(RHOP(I,J,K)+TWO_EPSILON_EB))
+               ENDDO
+            ENDDO
+         ENDDO
+      ENDIF
 
       ! Compute rho*D del Z
 
@@ -784,13 +792,14 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S)
 REAL(EB), POINTER, DIMENSION(:,:,:) :: FX_H_S,FY_H_S,FZ_H_S,RHO_H_S_P,U_DOT_DEL_RHO_H_S
 REAL(EB) :: UN,UN_P,TMP_F_GAS,DU_P,DU_M,DV_P,DV_M,DW_P,DW_M,DU,H_S
 REAL(EB), ALLOCATABLE, DIMENSION(:) :: ZZ_GET
-REAL(EB), DIMENSION(0:3,0:3,0:3) :: Z_TEMP,U_TEMP,F_TEMP
+REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP
+REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP
 INTEGER :: IC,I,J,K,IW
 TYPE(WALL_TYPE), POINTER :: WC
 TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
 TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
 
-RHO_H_S_P=>WORK1 ; RHO_H_S_P = 0._EB
+RHO_H_S_P=>WORK_PAD ; RHO_H_S_P = 0._EB
 FX_H_S=>WORK2    ; FX_H_S = 0._EB
 FY_H_S=>WORK3    ; FY_H_S = 0._EB
 FZ_H_S=>WORK4    ; FZ_H_S = 0._EB
@@ -801,9 +810,9 @@ SUBROUTINE ENTHALPY_ADVECTION(U_DOT_DEL_RHO_H_S)
 !$OMP PARALLEL PRIVATE(ZZ_GET, H_S)
 ALLOCATE(ZZ_GET(1:N_TRACKED_SPECIES))
 !$OMP DO SCHEDULE(static)
-DO K=0,KBP1
-   DO J=0,JBP1
-      DO I=0,IBP1
+DO K=-1,KBP1+1
+   DO J=-1,JBP1+1
+      DO I=-1,IBP1+1
          ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(I,J,K,1:N_TRACKED_SPECIES)
          CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,H_S,TMP(I,J,K))
          RHO_H_S_P(I,J,K) = RHOP(I,J,K)*H_S
@@ -964,7 +973,8 @@ SUBROUTINE SPECIES_ADVECTION_PART_1
 REAL(EB), POINTER, DIMENSION(:,:,:) :: RHO_Z_P,RHO_RMW
 REAL(EB), POINTER, DIMENSION(:,:,:,:) :: FX_ZZ,FY_ZZ,FZ_ZZ
 REAL(EB) :: ZZ_GET(1:N_TRACKED_SPECIES),MW_G
-REAL(EB), DIMENSION(0:3,0:3,0:3) :: Z_TEMP,U_TEMP,F_TEMP
+REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP
+REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP
 INTEGER :: I,J,K,IW,N !,II,JJ,KK,IOR,IC,IIG,JJG,KKG
 TYPE(WALL_TYPE), POINTER :: WC
 TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
@@ -973,16 +983,16 @@ SUBROUTINE SPECIES_ADVECTION_PART_1
 FX_ZZ=>SWORK1
 FY_ZZ=>SWORK2
 FZ_ZZ=>SWORK3
-RHO_Z_P=>WORK6
+RHO_Z_P=>WORK_PAD
 
 ! Species face values
 
 SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS
 
    !$OMP PARALLEL DO
-   DO K=0,KBP1
-      DO J=0,JBP1
-         DO I=0,IBP1
+   DO K=-1,KBP1+1
+      DO J=-1,JBP1+1
+         DO I=-1,IBP1+1
             RHO_Z_P(I,J,K) = RHOP(I,J,K)*ZZP(I,J,K,N)
          ENDDO
       ENDDO
@@ -1068,12 +1078,12 @@ SUBROUTINE SPECIES_ADVECTION_PART_1
 
    ! Repeat the above for DENSITY
 
-   RHO_RMW=>WORK6
+   RHO_RMW=>WORK_PAD
 
    !$OMP PARALLEL DO PRIVATE(ZZ_GET,MW_G) SCHEDULE(STATIC)
-   DO K=0,KBP1
-      DO J=0,JBP1
-         DO I=0,IBP1
+   DO K=-1,KBP1+1
+      DO J=-1,JBP1+1
+         DO I=-1,IBP1+1
             ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(I,J,K,1:N_TRACKED_SPECIES)
             CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_G)
             RHO_RMW(I,J,K) = RHOP(I,J,K)/MW_G
@@ -1189,7 +1199,6 @@ END SUBROUTINE SPECIES_ADVECTION_PART_1
 
 SUBROUTINE SPECIES_ADVECTION_PART_2(N,U_DOT_DEL_RHO_Z)
 
-USE MATH_FUNCTIONS, ONLY: GET_SCALAR_FACE_VALUE
 INTEGER, INTENT(IN) :: N
 REAL(EB), POINTER, DIMENSION(:,:,:) :: U_DOT_DEL_RHO_Z
 REAL(EB), POINTER, DIMENSION(:,:,:,:) :: FX_ZZ,FY_ZZ,FZ_ZZ

Source/fire.f90

@@ -1672,11 +1672,11 @@ SUBROUTINE CONDENSATION_EVAPORATION(DT,NM)
 CALL POINT_TO_MESH(NM)
 
 ZZ_INTERIM=> SWORK1
-ZZ_INTERIM(:,:,:,1:) = ZZ(:,:,:,1:) ! Lower bound may be zero for MW flux correction
+ZZ_INTERIM(0:IBP1,0:JBP1,0:KBP1,1:) = ZZ(0:IBP1,0:JBP1,0:KBP1,1:) ! Lower bound may be zero for MW flux correction
 RHO_INTERIM => WORK1
-RHO_INTERIM = RHO
+RHO_INTERIM(0:IBP1,0:JBP1,0:KBP1) = RHO(0:IBP1,0:JBP1,0:KBP1)
 TMP_INTERIM => WORK2
-TMP_INTERIM = TMP
+TMP_INTERIM(0:IBP1,0:JBP1,0:KBP1) = TMP(0:IBP1,0:JBP1,0:KBP1)
 
 DO IW=1,N_EXTERNAL_WALL_CELLS+N_INTERNAL_WALL_CELLS
    WC => WALL(IW)

Source/func.f90

@@ -1293,7 +1293,7 @@ END FUNCTION F_B
 
 SUBROUTINE GET_SCALAR_FACE_VALUE(A,U,F,I1,I2,J1,J2,K1,K2,IOR,LIMITER)
 
-REAL(EB), INTENT(IN) :: A(0:,0:,0:),U(0:,0:,0:)
+REAL(EB), INTENT(IN) :: A(0:,0:,0:),U(-1:,-1:,-1:)
 REAL(EB), INTENT(OUT) :: F(0:,0:,0:)
 INTEGER, INTENT(IN) :: LIMITER,I1,I2,J1,J2,K1,K2,IOR
 REAL(EB) :: R,B,DU_UP,DU_LOC,V(-2:2)

Source/init.f90

@@ -520,9 +520,17 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
 JBAR =>M%JBAR
 KBAR =>M%KBAR
 
-ALLOCATE(M%RHO(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)   ; CALL ChkMemErr('INIT','RHO',IZERO)
-ALLOCATE(M%RHOS(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)  ; CALL ChkMemErr('INIT','RHOS',IZERO) ; M%RHOS = RHOA
-ALLOCATE(M%TMP(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)   ; CALL ChkMemErr('INIT','TMP',IZERO)
+! Special arrays that require an extra layer of ghost cells for ease in calculating flux limiters at interpolated boundaries
+
+ALLOCATE(M%TMP( -1:IBP1+1,-1:JBP1+1,-1:KBP1+1),STAT=IZERO) ; CALL ChkMemErr('INIT','TMP',IZERO) ; M%TMP = TMPA
+ALLOCATE(M%RHO( -1:IBP1+1,-1:JBP1+1,-1:KBP1+1),STAT=IZERO) ; CALL ChkMemErr('INIT','RHO',IZERO) ; M%RHO = RHOA
+ALLOCATE(M%RHOS(-1:IBP1+1,-1:JBP1+1,-1:KBP1+1),STAT=IZERO) ; CALL ChkMemErr('INIT','RHOS',IZERO) ; M%RHOS = RHOA
+ALLOCATE(M%ZZ(  -1:IBP1+1,-1:JBP1+1,-1:KBP1+1,N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','ZZ',IZERO)  ; M%ZZ=0._EB
+ALLOCATE(M%ZZS( -1:IBP1+1,-1:JBP1+1,-1:KBP1+1,N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','ZZS',IZERO) ; M%ZZS=0._EB
+ALLOCATE(M%WORK_PAD(-1:IBP1+1,-1:JBP1+1,-1:KBP1+1),STAT=IZERO) ; CALL ChkMemErr('INIT','IWORK1',IZERO)
+
+! Regular 3D arrays
+
 ALLOCATE(M%U(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)     ; CALL ChkMemErr('INIT','U',IZERO)
 ALLOCATE(M%V(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)     ; CALL ChkMemErr('INIT','V',IZERO)
 ALLOCATE(M%W(0:IBP1,0:JBP1,0:KBP1),STAT=IZERO)     ; CALL ChkMemErr('INIT','W',IZERO)
@@ -576,8 +584,6 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
 
 ! Allocate species arrays
 
-ALLOCATE( M%ZZ(0:IBP1,0:JBP1,0:KBP1,N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','ZZ',IZERO)  ; M%ZZ=0._EB
-ALLOCATE(M%ZZS(0:IBP1,0:JBP1,0:KBP1,N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','ZZS',IZERO) ; M%ZZS=0._EB
 ALLOCATE(M%DEL_RHO_D_DEL_Z(0:IBP1,0:JBP1,0:KBP1,N_TOTAL_SCALARS),STAT=IZERO) ; CALL ChkMemErr('INIT','DEL_RHO_D_DEL_Z',IZERO)
 M%DEL_RHO_D_DEL_Z = 0._EB
 

Source/mass.f90

@@ -24,7 +24,8 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM)
 
 INTEGER, INTENT(IN) :: NM
 REAL(EB) :: TNOW,MW_F,MW_G,ZZ_GET(1:N_TRACKED_SPECIES)
-REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,Z_TEMP,F_TEMP
+REAL(EB), DIMENSION(0:3,0:3,0:3) :: U_TEMP,F_TEMP
+REAL(EB), DIMENSION(-1:3,-1:3,-1:3) :: Z_TEMP
 REAL(EB), PARAMETER :: DUMMY=0._EB
 INTEGER  :: I,J,K,N,IOR,IW,IIG,JJG,KKG,II,JJ,KK,IC
 REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,RHOP
@@ -62,12 +63,12 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM)
 
 SPECIES_LOOP: DO N=1,N_TOTAL_SCALARS
 
-   RHO_Z_P=>WORK1
+   RHO_Z_P=>WORK_PAD
 
    !$OMP PARALLEL DO
-   DO K=0,KBP1
-      DO J=0,JBP1
-         DO I=0,IBP1
+   DO K=-1,KBP1+1
+      DO J=-1,JBP1+1
+         DO I=-1,IBP1+1
             RHO_Z_P(I,J,K) = RHOP(I,J,K)*ZZP(I,J,K,N)
          ENDDO
       ENDDO
@@ -182,12 +183,12 @@ SUBROUTINE MASS_FINITE_DIFFERENCES(NM)
 
    ! Repeat the above for DENSITY
 
-   RHO_RMW=>WORK1
+   RHO_RMW=>WORK_PAD
 
    !$OMP PARALLEL DO PRIVATE(ZZ_GET,MW_G) SCHEDULE(STATIC)
-   DO K=0,KBP1
-      DO J=0,JBP1
-         DO I=0,IBP1
+   DO K=-1,KBP1+1
+      DO J=-1,JBP1+1
+         DO I=-1,IBP1+1
             ZZ_GET(1:N_TRACKED_SPECIES) = ZZP(I,J,K,1:N_TRACKED_SPECIES)
             CALL GET_MOLECULAR_WEIGHT(ZZ_GET,MW_G)
             RHO_RMW(I,J,K) = RHOP(I,J,K)/MW_G
@@ -444,7 +445,9 @@ SUBROUTINE DENSITY(T,DT,NM)
 
    ! Add gas production source term
 
-   IF (ALLOCATED(MESHES(NM)%M_DOT_PPP)) ZZS = ZZS + DT*M_DOT_PPP
+   IF (ALLOCATED(MESHES(NM)%M_DOT_PPP)) ZZS(0:IBP1,0:JBP1,0:KBP1,1:N_TRACKED_SPECIES) = &
+                                        ZZS(0:IBP1,0:JBP1,0:KBP1,1:N_TRACKED_SPECIES) + &
+                                        DT*M_DOT_PPP(0:IBP1,0:JBP1,0:KBP1,1:N_TRACKED_SPECIES)
 
    ! Manufactured solution
 
@@ -605,7 +608,8 @@ SUBROUTINE DENSITY(T,DT,NM)
    ! Add gas production source term
 
    IF (ALLOCATED(MESHES(NM)%M_DOT_PPP)) THEN
-      ZZ = ZZ + 0.5_EB*DT*M_DOT_PPP
+      ZZ(0:IBP1,0:JBP1,0:KBP1,1:N_TRACKED_SPECIES) = ZZ(0:IBP1,0:JBP1,0:KBP1,1:N_TRACKED_SPECIES) + &
+                                                     0.5_EB*DT*M_DOT_PPP(0:IBP1,0:JBP1,0:KBP1,1:N_TRACKED_SPECIES)
       IF (.NOT. CC_IBM) THEN ! We will use these for Regular cells in cut-cell region in CC_DENSITY.
          M_DOT_PPP = 0._EB
          D_SOURCE  = 0._EB

Source/mesh.f90

@@ -116,6 +116,7 @@ MODULE MESH_VARIABLES
    REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: SWORK4
    REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: WORK1,WORK2,WORK3,WORK4,WORK5,WORK6,WORK7,WORK8,WORK9
    INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: IWORK1
+   REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: WORK_PAD
    REAL(EB),     ALLOCATABLE, DIMENSION(:,:,:) :: PWORK1,PWORK2,PWORK3,PWORK4
    COMPLEX(EB),  ALLOCATABLE, DIMENSION(:,:,:) :: PWORK5,PWORK6,PWORK7,PWORK8
    REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: TURB_WORK1,TURB_WORK2,TURB_WORK3,TURB_WORK4
@@ -371,6 +372,7 @@ MODULE MESH_POINTERS
 INTEGER, POINTER, DIMENSION(:,:,:) :: MUNKH
 REAL(EB), POINTER, DIMENSION(:,:,:) :: WORK1,WORK2,WORK3,WORK4,WORK5,WORK6,WORK7,WORK8,WORK9
 INTEGER, POINTER, DIMENSION(:,:,:) :: IWORK1
+REAL(EB), POINTER, DIMENSION(:,:,:) :: WORK_PAD
 
 REAL(EB),     POINTER, DIMENSION(:,:,:) :: PWORK1,PWORK2,PWORK3,PWORK4
 COMPLEX(EB),  POINTER, DIMENSION(:,:,:) :: PWORK5,PWORK6,PWORK7,PWORK8
@@ -617,6 +619,7 @@ SUBROUTINE POINT_TO_MESH(NM)
 WORK8=>M%WORK8
 WORK9=>M%WORK9
 IWORK1=>M%IWORK1
+WORK_PAD=>M%WORK_PAD
 PWORK1=>M%PWORK1
 PWORK2=>M%PWORK2
 PWORK3=>M%PWORK3

Source/part.f90

@@ -3630,9 +3630,9 @@ SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM)
       B2%WORK1 = B1%TMP_F
    ENDDO
 
-   RHO_INTERIM => WORK1 ; RHO_INTERIM = RHO
-   TMP_INTERIM => WORK2 ; TMP_INTERIM = TMP
-   ZZ_INTERIM  => SWORK1 ; ZZ_INTERIM(:,:,:,1:) = ZZ(:,:,:,1:)
+   RHO_INTERIM => WORK1 ; RHO_INTERIM(0:IBP1,0:JBP1,0:KBP1) = RHO(0:IBP1,0:JBP1,0:KBP1)
+   TMP_INTERIM => WORK2 ; TMP_INTERIM(0:IBP1,0:JBP1,0:KBP1) = TMP(0:IBP1,0:JBP1,0:KBP1)
+   ZZ_INTERIM  => SWORK1 ; ZZ_INTERIM(0:IBP1,0:JBP1,0:KBP1,1:) = ZZ(0:IBP1,0:JBP1,0:KBP1,1:)
 
 ENDIF
 

Source/pres.f90

@@ -5256,7 +5256,7 @@ SUBROUTINE PRESSURE_SOLVER_CHECK_RESIDUALS_U(NM)
 
 IF (ITERATE_BAROCLINIC_TERM) THEN
    P => WORK7
-   P = RHOP*(HP-KRES)
+   P(0:IBP1,0:JBP1,0:KBP1) = RHOP(0:IBP1,0:JBP1,0:KBP1)*(HP(0:IBP1,0:JBP1,0:KBP1)-KRES(0:IBP1,0:JBP1,0:KBP1))
    RESIDUAL => WORK8(1:IBAR,1:JBAR,1:KBAR); RESIDUAL = 0._EB
    IF(CC_IBM) CALL COMPUTE_LINKED_CUTFACE_BAROCLINIC(NM,HP,RHOP,P)
    !$OMP PARALLEL PRIVATE(I,J,K,RHSS,LHSS,IMFCT,JMFCT,KMFCT,IPFCT,JPFCT,KPFCT)

Source/turb.f90

@@ -443,7 +443,7 @@ SUBROUTINE VARDEN_DYNSMAG(NM)
 REAL(EB) :: TEMP_TERM,DUDY,DUDZ,DVDX,DVDZ,DWDX,DWDY,ONTHDIV
 INTEGER :: I,J,K
 
-REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,UP,VP,WP,RHOP,RHOPHAT
+REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,UP,VP,WP,RHOP,RHOPHAT,RHOPTMP
 REAL(EB), POINTER, DIMENSION(:,:,:) :: S11,S22,S33,S12,S13,S23,SS
 REAL(EB), POINTER, DIMENSION(:,:,:) :: SHAT11,SHAT22,SHAT33,SHAT12,SHAT13,SHAT23,SSHAT
 REAL(EB), POINTER, DIMENSION(:,:,:) :: BETA11,BETA22,BETA33,BETA12,BETA13,BETA23
@@ -620,7 +620,9 @@ SUBROUTINE VARDEN_DYNSMAG(NM)
 ! test filter the density
 
 RHOPHAT => WORK7
-CALL TEST_FILTER(RHOPHAT,RHOP)
+RHOPTMP => WORK1
+RHOPTMP(0:IBP1,0:JBP1,0:KBP1) = RHOP(0:IBP1,0:JBP1,0:KBP1)
+CALL TEST_FILTER(RHOPHAT,RHOPTMP)
 
 ! calculate the Mij tensor
 
@@ -750,12 +752,12 @@ SUBROUTINE CALC_VARDEN_LEONARD_TERM
 RUW => TURB_WORK8
 RVW => TURB_WORK9
 
-RUU = RHOP*UP*UP
-RVV = RHOP*VP*VP
-RWW = RHOP*WP*WP
-RUV = RHOP*UP*VP
-RUW = RHOP*UP*WP
-RVW = RHOP*VP*WP
+RUU(0:IBP1,0:JBP1,0:KBP1) = RHOP(0:IBP1,0:JBP1,0:KBP1)*UP(0:IBP1,0:JBP1,0:KBP1)*UP(0:IBP1,0:JBP1,0:KBP1)
+RVV(0:IBP1,0:JBP1,0:KBP1) = RHOP(0:IBP1,0:JBP1,0:KBP1)*VP(0:IBP1,0:JBP1,0:KBP1)*VP(0:IBP1,0:JBP1,0:KBP1)
+RWW(0:IBP1,0:JBP1,0:KBP1) = RHOP(0:IBP1,0:JBP1,0:KBP1)*WP(0:IBP1,0:JBP1,0:KBP1)*WP(0:IBP1,0:JBP1,0:KBP1)
+RUV(0:IBP1,0:JBP1,0:KBP1) = RHOP(0:IBP1,0:JBP1,0:KBP1)*UP(0:IBP1,0:JBP1,0:KBP1)*VP(0:IBP1,0:JBP1,0:KBP1)
+RUW(0:IBP1,0:JBP1,0:KBP1) = RHOP(0:IBP1,0:JBP1,0:KBP1)*UP(0:IBP1,0:JBP1,0:KBP1)*WP(0:IBP1,0:JBP1,0:KBP1)
+RVW(0:IBP1,0:JBP1,0:KBP1) = RHOP(0:IBP1,0:JBP1,0:KBP1)*VP(0:IBP1,0:JBP1,0:KBP1)*WP(0:IBP1,0:JBP1,0:KBP1)
 
 ! extrapolate to ghost cells
 
@@ -788,9 +790,9 @@ SUBROUTINE CALC_VARDEN_LEONARD_TERM
 RV => TURB_WORK2
 RW => TURB_WORK3
 
-RU = RHOP*UP
-RV = RHOP*VP
-RW = RHOP*WP
+RU(0:IBP1,0:JBP1,0:KBP1) = RHOP(0:IBP1,0:JBP1,0:KBP1)*UP(0:IBP1,0:JBP1,0:KBP1)
+RV(0:IBP1,0:JBP1,0:KBP1) = RHOP(0:IBP1,0:JBP1,0:KBP1)*VP(0:IBP1,0:JBP1,0:KBP1)
+RW(0:IBP1,0:JBP1,0:KBP1) = RHOP(0:IBP1,0:JBP1,0:KBP1)*WP(0:IBP1,0:JBP1,0:KBP1)
 
 ! extrapolate to ghost cells