Merge pull request #12221 from mcgratta/master

FDS Source: Make NVEC (normal vector) a component of BOUNDARY_COORD

Source/ccib.f90

@@ -3495,10 +3495,11 @@ SUBROUTINE CC_COMPUTE_KRES(APPLY_TO_ESTIMATED_VARIABLES,NM)
                AF      = CUT_FACE(IFC2)%AREA(IFACE2)
                VELN = (1._EB-PRFCT)*CUT_FACE(IFC2)%VEL( IFACE2) + PRFCT *CUT_FACE(IFC2)%VELS(IFACE2)
                ! - to use velocity into gasphase, projected area.
+               BC => MESHES(NM)%BOUNDARY_COORD(CFACE(ICFA)%BC_INDEX)
                DO X1AXIS=IAXIS,KAXIS
-                 AUI  = ABS(CFACE(ICFA)%NVEC(X1AXIS)*AF)
+                 AUI  = ABS(BC%NVEC(X1AXIS)*AF)
                  ATOTV(X1AXIS) = ATOTV(X1AXIS) + AUI
-                 UVWAV(X1AXIS) = UVWAV(X1AXIS) - VELN*CFACE(ICFA)%NVEC(X1AXIS)*AUI
+                 UVWAV(X1AXIS) = UVWAV(X1AXIS) - VELN*BC%NVEC(X1AXIS)*AUI
                ENDDO
 
             END SELECT
@@ -3541,10 +3542,11 @@ SUBROUTINE CC_COMPUTE_KRES(APPLY_TO_ESTIMATED_VARIABLES,NM)
                AF      = CUT_FACE(IFC2)%AREA(IFACE2)
                VELN = (1._EB-PRFCT)*CUT_FACE(IFC2)%VEL( IFACE2) + PRFCT *CUT_FACE(IFC2)%VELS(IFACE2)
                ! - to use velocity into gasphase, projected area.
+               BC => MESHES(NM)%BOUNDARY_COORD(CFACE(ICFA)%BC_INDEX)
                DO X1AXIS=IAXIS,KAXIS
-                 AUI  = ABS(CFACE(ICFA)%NVEC(X1AXIS)*AF)
+                 AUI  = ABS(BC%NVEC(X1AXIS)*AF)
                  ATOTV(X1AXIS) = ATOTV(X1AXIS) + AUI
-                 UVWAV(X1AXIS) = UVWAV(X1AXIS) - VELN*CFACE(ICFA)%NVEC(X1AXIS)*AUI
+                 UVWAV(X1AXIS) = UVWAV(X1AXIS) - VELN*BC%NVEC(X1AXIS)*AUI
                ENDDO
 
             END SELECT
@@ -5844,12 +5846,13 @@ SUBROUTINE CC_COMPUTE_VISCOSITY(DT,NM)
 DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
 
    CFA => CFACE(ICF)
+   BC => BOUNDARY_COORD(CFA%BC_INDEX)
    B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
    B2 => BOUNDARY_PROP2(CFA%B2_INDEX)
    SF => SURFACE(CFA%SURF_INDEX)
 
    ! Surface Velocity:
-   NVEC = CFA%NVEC
+   NVEC = BC%NVEC
 
    ! right now VEL_T not defined for CFACEs
    TVEC1=(/ 0._EB,0._EB,0._EB/)
@@ -6825,7 +6828,7 @@ SUBROUTINE CFACE_THERMAL_GASVARS(ICF,B1)
             VVEL(VIND) = VVEL(VIND) + CUT_FACE(IND1)%INT_COEF(INPE)*CUT_FACE(IND1)%INT_FVARS(INT_VEL_IND,INPE)
          ENDDO
       ENDDO
-      V_TANG(IAXIS:KAXIS) = VVEL(IAXIS:KAXIS) - DOT_PRODUCT(VVEL,CFA%NVEC)*CFA%NVEC(IAXIS:KAXIS)
+      V_TANG(IAXIS:KAXIS) = VVEL(IAXIS:KAXIS) - DOT_PRODUCT(VVEL,BC%NVEC)*BC%NVEC(IAXIS:KAXIS)
       B1%U_TANG = SQRT(V_TANG(IAXIS)**2._EB+V_TANG(JAXIS)**2._EB+V_TANG(KAXIS)**2._EB)
 
    ELSE
@@ -15940,11 +15943,12 @@ SUBROUTINE GET_BOUND_VEL(X1AXIS,INBFC_CFCEN,XYZ_PP,VELX1)
 
 CFA => CFACE(ICF)
 B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
+BC => BOUNDARY_COORD(CFA%BC_INDEX)
 ! Velocity into Gas Region, component along X1AXIS:
 IF (PREDICTOR) THEN
-   VELX1 = -B1%U_NORMAL * CFACE(ICF)%NVEC(X1AXIS)
+   VELX1 = -B1%U_NORMAL * BC%NVEC(X1AXIS)
 ELSE
-   VELX1 = -B1%U_NORMAL_S* CFACE(ICF)%NVEC(X1AXIS)
+   VELX1 = -B1%U_NORMAL_S* BC%NVEC(X1AXIS)
 ENDIF
 
 RETURN
@@ -16583,7 +16587,7 @@ SUBROUTINE GET_CRTCFCC_INT_STENCILS
          N_CVAR_START = NPE_LIST_START
          DO IFACE=1,CUT_FACE(ICF)%NFACE
             ICF1 = CUT_FACE(ICF)%CFACE_INDEX(IFACE)
-            DV(IAXIS:KAXIS) = CFACE(ICF1)%NVEC(IAXIS:KAXIS)
+            DV(IAXIS:KAXIS) = MESHES(NM)%BOUNDARY_COORD(CFACE(ICF1)%BC_INDEX)%NVEC(IAXIS:KAXIS)
             CALL GET_DELN(1.001_EB,DELN,DXCELL(I),DYCELL(J),DZCELL(K),NVEC=DV,CLOSE_PT=.TRUE.)
             ! Origin in boundary point XYZ_PP(IAXIS:KAXIS):
             XYZ_PP(IAXIS:KAXIS)         = CUT_FACE(ICF)%XYZCEN(IAXIS:KAXIS,IFACE)
@@ -16611,7 +16615,7 @@ SUBROUTINE GET_CRTCFCC_INT_STENCILS
             XYZ_PP(IAXIS:KAXIS)         = CUT_FACE(ICF)%XYZCEN(IAXIS:KAXIS,IFACE)
             FOUND_INBFC(1:3) = (/ CC_FTYPE_CFINB, ICF, IFACE /)
             ICF1 = CUT_FACE(ICF)%CFACE_INDEX(IFACE)
-            DV(IAXIS:KAXIS) = CFACE(ICF1)%NVEC(IAXIS:KAXIS)
+            DV(IAXIS:KAXIS) = MESHES(NM)%BOUNDARY_COORD(CFACE(ICF1)%BC_INDEX)%NVEC(IAXIS:KAXIS)
             INT_NOUT(IAXIS:KAXIS,IFACE) = DV(IAXIS:KAXIS)
             INT_XN(0:INT_N_EXT_PTS) = 0._EB
             ! Initialize interpolation coefficients along the normal probe direction DV

Source/dump.f90

@@ -5265,7 +5265,7 @@ SUBROUTINE DUMP_CFACES_GEOM(FUNIT,STIME)
      IND1 = CFACE(ICF)%CUT_FACE_IND1
      IND2 = CFACE(ICF)%CUT_FACE_IND2
      XLOC(1:3) = (/ CFACE(ICF)%X, CFACE(ICF)%Y, CFACE(ICF)%Z /)
-     XLOC(4:6) = XLOC(1:3) + CUT_FACE(IND1)%INT_XN(1,IND2)*CFACE(ICF)%NVEC(IAXIS:KAXIS)
+     XLOC(4:6) = XLOC(1:3) + CUT_FACE(IND1)%INT_XN(1,IND2)*BOUNDARY_COORD(CFACE(ICF)%BC_INDEX)%NVEC(IAXIS:KAXIS)
      VERTS(6*(ICF-INTERNAL_CFACE_CELLS_LB-1)+1:6*(ICF-INTERNAL_CFACE_CELLS_LB)) = REAL(XLOC(1:6),FB)
    ENDDO
    WRITE(FUNIT) VERTS(1:6*N_INTERNAL_CFACE_CELLS)
@@ -9001,7 +9001,7 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(NM,INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL
 
          PVEC = ( PR1 - (PR2-PR1)*Z1**2 / (Z2**2-Z1**2) ) * NVEC ! surface normal pressure force
       ELSEIF (PRESENT(OPT_CFACE_INDEX)) THEN
-         NVEC = CFA%NVEC
+         NVEC = BC%NVEC
          ! find cut-cell adjacent to CFACE
          IND1 = CFA%CUT_FACE_IND1
          IND2 = CFA%CUT_FACE_IND2
@@ -9038,7 +9038,7 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(NM,INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL
          W_CELL = 0.5_EB*(W(IIG,JJG,KKG-1)+W(IIG,JJG,KKG))
          MU_WALL = MU_DNS(IIG,JJG,KKG)
       ELSEIF (PRESENT(OPT_CFACE_INDEX)) THEN
-         NVEC = CFA%NVEC
+         NVEC = BC%NVEC
          ! right now VEL_T not defined for CFACEs
          TVEC1=(/ 0._EB,0._EB,0._EB/)
          TVEC2=(/ 0._EB,0._EB,0._EB/)
@@ -9093,7 +9093,7 @@ REAL(EB) FUNCTION SOLID_PHASE_OUTPUT(NM,INDX,Y_INDEX,Z_INDEX,PART_INDEX,OPT_WALL
          KKG = BC%KKG
          PVEC = RHO(IIG,JJG,KKG)*H(IIG,JJG,KKG) * NVEC ! surface normal pressure force
       ELSEIF (PRESENT(OPT_CFACE_INDEX)) THEN
-         NVEC = CFA%NVEC
+         NVEC = BC%NVEC
          ! find cut-cell adjacent to CFACE
          IND1 = CFA%CUT_FACE_IND1
          IND2 = CFA%CUT_FACE_IND2

Source/func.f90

@@ -1638,9 +1638,6 @@ SUBROUTINE PACK_CFACE(NM,OS,CFA,SURF_INDEX,RC,IC,LC,UNPACK_IT,COUNT_ONLY)
 ! Assign and initialize reals
 
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%AREA,UNPACK_IT)
-RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%NVEC(1),UNPACK_IT)
-RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%NVEC(2),UNPACK_IT)
-RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%NVEC(3),UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%DUNDT,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%RSUM_G,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),CFA%TMP_G,UNPACK_IT)
@@ -1701,6 +1698,9 @@ SUBROUTINE PACK_BOUNDARY_COORD(NM,IC,RC,OS,BC_INDEX,UNPACK_IT,COUNT_ONLY)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%Y2,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%Z1,UNPACK_IT)
 RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%Z2,UNPACK_IT)
+RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%NVEC(1),UNPACK_IT)
+RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%NVEC(2),UNPACK_IT)
+RC=RC+1 ; IF (.NOT.COUNT_ONLY) CALL EQUATE(OS%REALS(RC),BC%NVEC(3),UNPACK_IT)
 
 END SUBROUTINE PACK_BOUNDARY_COORD
 

Source/geom.f90

@@ -669,7 +669,7 @@ SUBROUTINE POINT_IN_CFACE(NM,XP,YP,ZP,CFACE_INDEX,IN_CFACE)
 ! Normal, max normal component, define plane X2AXIS,X3AXIS to do search:
 VERT_CUTFACE = SIZE(MESHES(NM)%CUT_FACE(INBFC)%CFELEM, DIM=1); ALLOCATE(CFELEM(1:VERT_CUTFACE))
 CFELEM(1:VERT_CUTFACE)  = MESHES(NM)%CUT_FACE(INBFC)%CFELEM(1:VERT_CUTFACE,INBFC_LOC)
-NVEC(IAXIS:KAXIS) => CFACE(CFACE_INDEX)%NVEC(IAXIS:KAXIS)
+NVEC(IAXIS:KAXIS) => MESHES(NM)%BOUNDARY_COORD(CFACE(CFACE_INDEX)%BC_INDEX)%NVEC(IAXIS:KAXIS)
 
 ! Plane equation for INBOUNDARY cut-face plane:
 ! Location of first point in cf polygon is P0:
@@ -8516,13 +8516,13 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
       ! Normal to cut-face:
       V2(IAXIS:KAXIS) = CF%XYZVERT(IAXIS:KAXIS,CF%CFELEM(2,IFACE))-CF%XYZCEN(IAXIS:KAXIS,IFACE)
       V3(IAXIS:KAXIS) = CF%XYZVERT(IAXIS:KAXIS,CF%CFELEM(3,IFACE))-CF%XYZCEN(IAXIS:KAXIS,IFACE)
-      CALL CROSS_PRODUCT(CFA%NVEC(IAXIS:KAXIS),V2,V3)
-      IF(NORM2(CFA%NVEC)>TWO_EPSILON_EB .AND. CF%CFACE_ORIGIN(IFACE)==BLOCKED_SPLIT_CELL) THEN
-         CFA%NVEC(IAXIS:KAXIS) = CFA%NVEC(IAXIS:KAXIS)/NORM2(CFA%NVEC)
+      CALL CROSS_PRODUCT(BC%NVEC(IAXIS:KAXIS),V2,V3)
+      IF(NORM2(BC%NVEC)>TWO_EPSILON_EB .AND. CF%CFACE_ORIGIN(IFACE)==BLOCKED_SPLIT_CELL) THEN
+         BC%NVEC(IAXIS:KAXIS) = BC%NVEC(IAXIS:KAXIS)/NORM2(BC%NVEC)
       ELSE
          IBOD =CF%BODTRI(1,IFACE)
          IWSEL=CF%BODTRI(2,IFACE)
-         CFA%NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL)
+         BC%NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL)
       ENDIF
 
       ! Boundary CFACES processed are defined of type SOLID_BOUNDARY
@@ -8544,9 +8544,9 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
          WC_BC => M%BOUNDARY_COORD(WC%BC_INDEX)
          IOR = WC_BC%IOR
          SELECT CASE(ABS(IOR))
-            CASE(IAXIS); CFA%NVEC(IAXIS:KAXIS) = (/ REAL(SIGN(1,IOR),EB), 0._EB, 0._EB /)
-            CASE(JAXIS); CFA%NVEC(IAXIS:KAXIS) = (/ 0._EB, REAL(SIGN(1,IOR),EB), 0._EB /)
-            CASE(KAXIS); CFA%NVEC(IAXIS:KAXIS) = (/ 0._EB, 0._EB, REAL(SIGN(1,IOR),EB) /)
+            CASE(IAXIS); BC%NVEC(IAXIS:KAXIS) = (/ REAL(SIGN(1,IOR),EB), 0._EB, 0._EB /)
+            CASE(JAXIS); BC%NVEC(IAXIS:KAXIS) = (/ 0._EB, REAL(SIGN(1,IOR),EB), 0._EB /)
+            CASE(KAXIS); BC%NVEC(IAXIS:KAXIS) = (/ 0._EB, 0._EB, REAL(SIGN(1,IOR),EB) /)
          END SELECT
          BC%IOR = IOR
 

Source/init.f90

@@ -1140,7 +1140,7 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM)
          IW  = MAXLOC(M%CUT_FACE(ICF)%AREA(1:M%CUT_FACE(ICF)%NFACE),DIM=1)
          CFA => M%CFACE( M%CUT_FACE(ICF)%CFACE_INDEX(IW) )
          BC  => M%BOUNDARY_COORD(CFA%BC_INDEX)
-         IF (CFA%NVEC(KAXIS)>-TWO_EPSILON_EB .AND.  CFA%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN
+         IF (BC%NVEC(KAXIS)>-TWO_EPSILON_EB .AND.  CFA%BOUNDARY_TYPE==SOLID_BOUNDARY) THEN
             IF (BC%KKG > M%K_AGL_SLICE(BC%IIG,BC%JJG,NSLICE)) THEN
                M%K_AGL_SLICE(BC%IIG,BC%JJG,NSLICE) = MIN( M%KBAR , M%K_AGL_SLICE(BC%IIG,BC%JJG,NSLICE)+BC%KKG )
             ENDIF
@@ -3142,6 +3142,14 @@ SUBROUTINE INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,TT)
 BC%JJG = JJG
 BC%KKG = KKG
 BC%IOR = IOR
+SELECT CASE(BC%IOR)
+   CASE( 1) ; BC%NVEC=(/ 1._EB, 0._EB, 0._EB/)
+   CASE(-1) ; BC%NVEC=(/-1._EB, 0._EB, 0._EB/)
+   CASE( 2) ; BC%NVEC=(/ 0._EB, 1._EB, 0._EB/)
+   CASE(-2) ; BC%NVEC=(/ 0._EB,-1._EB, 0._EB/)
+   CASE( 3) ; BC%NVEC=(/ 0._EB, 0._EB, 1._EB/)
+   CASE(-3) ; BC%NVEC=(/ 0._EB, 0._EB,-1._EB/)
+END SELECT
 BC%X = XW
 BC%Y = YW
 BC%Z = ZW
@@ -3510,6 +3518,14 @@ SUBROUTINE INIT_THIN_WALL_CELL(NM,I,J,K,OBST_INDEX,ITW,IOR,SURF_INDEX,IEC)
 BC%JJG = J
 BC%KKG = K
 BC%IOR = IOR
+SELECT CASE(BC%IOR)
+   CASE( 1) ; BC%NVEC=(/ 1._EB, 0._EB, 0._EB/)
+   CASE(-1) ; BC%NVEC=(/-1._EB, 0._EB, 0._EB/)
+   CASE( 2) ; BC%NVEC=(/ 0._EB, 1._EB, 0._EB/)
+   CASE(-2) ; BC%NVEC=(/ 0._EB,-1._EB, 0._EB/)
+   CASE( 3) ; BC%NVEC=(/ 0._EB, 0._EB, 1._EB/)
+   CASE(-3) ; BC%NVEC=(/ 0._EB, 0._EB,-1._EB/)
+END SELECT
 BC%X  = M%X(I)
 BC%Y  = M%Y(J)
 BC%Z  = M%Z(K)

Source/part.f90

@@ -825,17 +825,17 @@ SUBROUTINE PARTICLE_FACE_INSERT(WALL_INDEX,CFACE_INDEX)
             ENDIF
          ELSEIF (PRESENT(CFACE_INDEX)) THEN
             CALL RANDOM_CFACE_XYZ(NM,CFA,CFA_X,CFA_Y,CFA_Z)
-            BC%X = CFA_X + CFA%NVEC(1)*VENT_OFFSET*DX(IIG)
-            BC%Y = CFA_Y + CFA%NVEC(2)*VENT_OFFSET*DY(JJG)
-            BC%Z = CFA_Z + CFA%NVEC(3)*VENT_OFFSET*DZ(KKG)
+            BC%X = CFA_X + BC%NVEC(1)*VENT_OFFSET*DX(IIG)
+            BC%Y = CFA_Y + BC%NVEC(2)*VENT_OFFSET*DY(JJG)
+            BC%Z = CFA_Z + BC%NVEC(3)*VENT_OFFSET*DZ(KKG)
             IF (INSERT_TYPE==1 .AND. LP%EMBER) THEN
                CALL RANDOM_NUMBER(RN3)
                BC%Z = CFA_Z + SF%EMBER_GENERATION_HEIGHT(1) + &
                       (SF%EMBER_GENERATION_HEIGHT(2)-SF%EMBER_GENERATION_HEIGHT(1))*REAL(RN3,EB)
             ENDIF
-            LP%U = DOT_PRODUCT(CFA%NVEC,(/-B1%U_NORMAL,SF%VEL_T(1),SF%VEL_T(2)/))
-            LP%V = DOT_PRODUCT(CFA%NVEC,(/SF%VEL_T(1),-B1%U_NORMAL,SF%VEL_T(2)/))
-            LP%W = DOT_PRODUCT(CFA%NVEC,(/SF%VEL_T(1),SF%VEL_T(2),-B1%U_NORMAL/))
+            LP%U = DOT_PRODUCT(BC%NVEC,(/-B1%U_NORMAL,SF%VEL_T(1),SF%VEL_T(2)/))
+            LP%V = DOT_PRODUCT(BC%NVEC,(/SF%VEL_T(1),-B1%U_NORMAL,SF%VEL_T(2)/))
+            LP%W = DOT_PRODUCT(BC%NVEC,(/SF%VEL_T(1),SF%VEL_T(2),-B1%U_NORMAL/))
          ENDIF WALL_OR_CFACE_IF_2
 
          ! Update idicies in case offset puts location in a different cell
@@ -1146,7 +1146,7 @@ SUBROUTINE INSERT_VOLUMETRIC_PARTICLES
                      IF (DIST<DIST_MIN) THEN
                         DIST_MIN = DIST
                         P_VECTOR_MIN = P_VECTOR
-                        NVEC_MIN = CFACE(CF%CFACE_INDEX(IFACE))%NVEC
+                        NVEC_MIN = BOUNDARY_COORD(CFACE(CF%CFACE_INDEX(IFACE))%BC_INDEX)%NVEC
                      ENDIF
                   ENDDO CFA_LOOP1
                   IF (DOT_PRODUCT(NVEC_MIN,P_VECTOR_MIN) > TWO_EPSILON_EB) CC_VALID=.TRUE.
@@ -1319,7 +1319,7 @@ SUBROUTINE INSERT_VOLUMETRIC_PARTICLES
                               IF (DIST<DIST_MIN) THEN
                                  DIST_MIN = DIST
                                  P_VECTOR_MIN = P_VECTOR
-                                 NVEC_MIN = CFACE(CF%CFACE_INDEX(IFACE))%NVEC
+                                 NVEC_MIN = BOUNDARY_COORD(CFACE(CF%CFACE_INDEX(IFACE))%BC_INDEX)%NVEC
                               ENDIF
                            ENDDO CFA_LOOP2
                            IF (DOT_PRODUCT(NVEC_MIN,P_VECTOR_MIN) > TWO_EPSILON_EB) EXIT RAND_LOCATION_LOOP
@@ -1994,16 +1994,16 @@ SUBROUTINE MOVE_PARTICLES(T,DT,NM)
                CFA2 => CFACE(ICF)
                CFA2_BC => BOUNDARY_COORD(CFA2%BC_INDEX)
 
-               DIST2 = DOT_PRODUCT(CFA%NVEC,GVEC/(NORM2(GVEC)+TWO_EPSILON_EB))
+               DIST2 = DOT_PRODUCT(CFA_BC%NVEC,GVEC/(NORM2(GVEC)+TWO_EPSILON_EB))
                IF(DIST2<-0.99_EB) THEN ! Only for slopes less than 8 degrees.
 
                   IN_CFACE=.TRUE.
                   ! Test for case of CFACE_INDEX switching to an ICF with similar slope, if so don't do anything.
-                  IF (.NOT.(LP%CFACE_INDEX/=ICF .AND. (DOT_PRODUCT(CFA%NVEC,CFA2%NVEC))>0.99_EB) ) THEN
+                  IF (.NOT.(LP%CFACE_INDEX/=ICF .AND. (DOT_PRODUCT(CFA_BC%NVEC,CFA2_BC%NVEC))>0.99_EB) ) THEN
 
-                     IF (LP%CFACE_INDEX/=ICF .AND. (DOT_PRODUCT(CFA%NVEC,CFA2%NVEC))<=0.99_EB) THEN
+                     IF (LP%CFACE_INDEX/=ICF .AND. (DOT_PRODUCT(CFA_BC%NVEC,CFA2_BC%NVEC))<=0.99_EB) THEN
                         ! Case of switching to different ICF with different slope:
-                        DIST2 = DOT_PRODUCT(CFA2%NVEC,GVEC/(NORM2(GVEC)+TWO_EPSILON_EB))
+                        DIST2 = DOT_PRODUCT(CFA2_BC%NVEC,GVEC/(NORM2(GVEC)+TWO_EPSILON_EB))
                         ! If ICF is almost vertical pointing in the direction of velocity, set creep velocity.
                         IF(ABS(DIST2)<0.01_EB .AND. (CFA2_BC%Z<CFA_BC%Z)) IN_CFACE=.FALSE.
                      ELSE
@@ -2058,9 +2058,9 @@ SUBROUTINE MOVE_PARTICLES(T,DT,NM)
                CFA => CFACE(ICF)
                CFA_BC => BOUNDARY_COORD(CFA%BC_INDEX)
                P_VECTOR = (/BC%X-CFA_BC%X, BC%Y-CFA_BC%Y, BC%Z-CFA_BC%Z/)
-               TEST_POS = .FALSE.; IF (LP%CFACE_INDEX == 0) TEST_POS = DOT_PRODUCT(CFA%NVEC,P_VECTOR) > TWO_EPSILON_EB
+               TEST_POS = .FALSE.; IF (LP%CFACE_INDEX == 0) TEST_POS = DOT_PRODUCT(CFA_BC%NVEC,P_VECTOR) > TWO_EPSILON_EB
 
-               CFACE_ATTACH : IF (DOT_PRODUCT(CFA%NVEC,GVEC)>0._EB .OR. TEST_POS) THEN
+               CFACE_ATTACH : IF (DOT_PRODUCT(CFA_BC%NVEC,GVEC)>0._EB .OR. TEST_POS) THEN
 
                   ! Normal points down or particle in gas phase. Let particle move freely:
 
@@ -2071,8 +2071,8 @@ SUBROUTINE MOVE_PARTICLES(T,DT,NM)
 
 
                   IF(LPC%ADHERE_TO_SOLID) THEN
-                     CALL CROSS_PRODUCT(VEL_VECTOR_1,CFA%NVEC,GVEC)
-                     CALL CROSS_PRODUCT(VEL_VECTOR_2,VEL_VECTOR_1,CFA%NVEC)
+                     CALL CROSS_PRODUCT(VEL_VECTOR_1,CFA_BC%NVEC,GVEC)
+                     CALL CROSS_PRODUCT(VEL_VECTOR_2,VEL_VECTOR_1,CFA_BC%NVEC)
                      CFACE_SLOPE : IF (NORM2(VEL_VECTOR_2) > TWO_EPSILON_EB .AND. ABS(LP%W) > TWO_EPSILON_EB) THEN
                         ! The surface is tilted; particles go down slope:
                         VEL_VECTOR_1 = VEL_VECTOR_2/NORM2(VEL_VECTOR_2)
@@ -2092,7 +2092,7 @@ SUBROUTINE MOVE_PARTICLES(T,DT,NM)
                      ELSEIF (LP%CFACE_INDEX /= 0 .AND. ABS(LP%W)<TWO_EPSILON_EB) THEN CFACE_SLOPE
                         ! Particle moving in horizontal direction and assumed crossing into solid.
                         ! Bounce back on random direction, maintaining CFACE_INDEX:
-                        IF (DOT_PRODUCT( (/ LP%U, LP%V /) ,CFA%NVEC(IAXIS:JAXIS))<-TWO_EPSILON_EB)THEN
+                        IF (DOT_PRODUCT( (/ LP%U, LP%V /) ,CFA_BC%NVEC(IAXIS:JAXIS))<-TWO_EPSILON_EB)THEN
                            CALL RANDOM_NUMBER(RN)
                            DIST2_MIN = (1._EB-SIGN(1._EB,LP%V))*PI/2._EB
                            IF(ABS(LP%U) > TWO_EPSILON_EB) DIST2_MIN = ATAN2(LP%V,LP%U)
@@ -2117,12 +2117,12 @@ SUBROUTINE MOVE_PARTICLES(T,DT,NM)
 
                   PVEC_L = NORM2(P_VECTOR)
                   IF (PVEC_L>TWO_EPSILON_EB) THEN
-                     THETA = ACOS(DOT_PRODUCT(CFACE(ICF)%NVEC,P_VECTOR/PVEC_L))
+                     THETA = ACOS(DOT_PRODUCT(BOUNDARY_COORD(CFACE(ICF)%BC_INDEX)%NVEC,P_VECTOR/PVEC_L))
                      IF (THETA>PIO2) THEN
                         DELTA = PVEC_L*SIN(THETA-0.5_EB*PI)+TWO_EPSILON_EB
-                        BC%X = BC%X + DELTA*CFA%NVEC(1)
-                        BC%Y = BC%Y + DELTA*CFA%NVEC(2)
-                        BC%Z = BC%Z + DELTA*CFA%NVEC(3)
+                        BC%X = BC%X + DELTA*CFA_BC%NVEC(1)
+                        BC%Y = BC%Y + DELTA*CFA_BC%NVEC(2)
+                        BC%Z = BC%Z + DELTA*CFA_BC%NVEC(3)
                      ENDIF
                   ENDIF
                   LP%CFACE_INDEX = ICF