FDS Source : Fix issues with obst and geometry intersection.

Source/ccib.f90

@@ -9991,8 +9991,8 @@ SUBROUTINE CCENTHALPY_ADVECTION
 ! Now work with boundary faces:
 ! IAXIS faces:
 X1AXIS = IAXIS
-DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-   IW = CC_REGFACE_IAXIS_Z(IFACE)%IWC
+DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+   IW = CC_REGFACE_IAXIS_Z(IFACE)%IWC; IF(IW<1) CYCLE
    IF( WALL(IW)%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
        WALL(IW)%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
        WALL(IW)%BOUNDARY_TYPE==PERIODIC_BOUNDARY ) CYCLE
@@ -10031,8 +10031,8 @@ SUBROUTINE CCENTHALPY_ADVECTION
 
 ! JAXIS faces:
 X1AXIS = JAXIS
-DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-   IW = CC_REGFACE_JAXIS_Z(IFACE)%IWC
+DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+   IW = CC_REGFACE_JAXIS_Z(IFACE)%IWC; IF(IW<1) CYCLE
    IF( WALL(IW)%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
        WALL(IW)%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
        WALL(IW)%BOUNDARY_TYPE==PERIODIC_BOUNDARY ) CYCLE
@@ -10071,8 +10071,8 @@ SUBROUTINE CCENTHALPY_ADVECTION
 
 ! KAXIS faces:
 X1AXIS = KAXIS
-DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-   IW = CC_REGFACE_KAXIS_Z(IFACE)%IWC
+DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+   IW = CC_REGFACE_KAXIS_Z(IFACE)%IWC; IF(IW<1) CYCLE
    IF( WALL(IW)%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
        WALL(IW)%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
        WALL(IW)%BOUNDARY_TYPE==PERIODIC_BOUNDARY ) CYCLE
@@ -10215,8 +10215,8 @@ SUBROUTINE CCENTHALPY_ADVECTION
 ENDDO
 
 ! Finally Gasphase cut-faces:
-DO ICF = 1,MESHES(NM)%N_BBCUTFACE_MESH
-   IF ( CUT_FACE(ICF)%STATUS /= CC_GASPHASE ) CYCLE
+DO ICF = 1,MESHES(NM)%N_CUTFACE_MESH
+   IF ( CUT_FACE(ICF)%STATUS /= CC_GASPHASE .OR. MESHES(NM)%CUT_FACE(ICF)%IWC<1) CYCLE
    IW = CUT_FACE(ICF)%IWC
    IF( WALL(IW)%BOUNDARY_TYPE==NULL_BOUNDARY ) CYCLE
    WC => WALL(IW)
@@ -10557,8 +10557,8 @@ SUBROUTINE CC_DIFFUSIVE_HEAT_FLUXES
 
    ! IAXIS faces:
    X1AXIS = IAXIS
-   DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-      IW = CC_REGFACE_IAXIS_Z(IFACE)%IWC; WC => WALL(IW)
+   DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+      IW = CC_REGFACE_IAXIS_Z(IFACE)%IWC; IF(IW<1) CYCLE; WC => WALL(IW)
       IF( ANY(WC%BOUNDARY_TYPE==(/NULL_BOUNDARY,INTERPOLATED_BOUNDARY,PERIODIC_BOUNDARY/)) ) CYCLE
       I  = CC_REGFACE_IAXIS_Z(IFACE)%IJK(IAXIS)
       J  = CC_REGFACE_IAXIS_Z(IFACE)%IJK(JAXIS)
@@ -10588,8 +10588,8 @@ SUBROUTINE CC_DIFFUSIVE_HEAT_FLUXES
 
    ! JAXIS faces:
    X1AXIS = JAXIS
-   DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-      IW = CC_REGFACE_JAXIS_Z(IFACE)%IWC; WC => WALL(IW)
+   DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+      IW = CC_REGFACE_JAXIS_Z(IFACE)%IWC; IF(IW<1) CYCLE; WC => WALL(IW)
       IF( ANY(WC%BOUNDARY_TYPE==(/NULL_BOUNDARY,INTERPOLATED_BOUNDARY,PERIODIC_BOUNDARY/)) ) CYCLE
       I  = CC_REGFACE_JAXIS_Z(IFACE)%IJK(IAXIS)
       J  = CC_REGFACE_JAXIS_Z(IFACE)%IJK(JAXIS)
@@ -10619,8 +10619,8 @@ SUBROUTINE CC_DIFFUSIVE_HEAT_FLUXES
 
    ! KAXIS faces:
    X1AXIS = KAXIS
-   DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-      IW = CC_REGFACE_KAXIS_Z(IFACE)%IWC; WC => WALL(IW)
+   DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+      IW = CC_REGFACE_KAXIS_Z(IFACE)%IWC; IF(IW<1) CYCLE; WC => WALL(IW)
       IF( ANY(WC%BOUNDARY_TYPE==(/NULL_BOUNDARY,INTERPOLATED_BOUNDARY,PERIODIC_BOUNDARY/)) ) CYCLE
       I  = CC_REGFACE_KAXIS_Z(IFACE)%IJK(IAXIS)
       J  = CC_REGFACE_KAXIS_Z(IFACE)%IJK(JAXIS)
@@ -10686,8 +10686,8 @@ SUBROUTINE CC_DIFFUSIVE_HEAT_FLUXES
 ENDDO
 
 ! GASPHASE cut-faces:
-DO ICF = 1,MESHES(NM)%N_BBCUTFACE_MESH
-   CF => CUT_FACE(ICF); IF ( CF%STATUS /= CC_GASPHASE ) CYCLE
+DO ICF = 1,MESHES(NM)%N_CUTFACE_MESH
+   CF => CUT_FACE(ICF); IF ( CF%STATUS /= CC_GASPHASE .OR. CF%IWC<1) CYCLE
    IW = CF%IWC; WC => WALL(IW)
    IF( ANY(WC%BOUNDARY_TYPE==(/NULL_BOUNDARY,INTERPOLATED_BOUNDARY,PERIODIC_BOUNDARY/)) ) CYCLE
    I  = CF%IJK(IAXIS); J = CF%IJK(JAXIS); K = CF%IJK(KAXIS); X1AXIS = CF%IJK(KAXIS+1)
@@ -11025,8 +11025,8 @@ SUBROUTINE CC_CONDUCTION_HEAT_FLUX
 ! Now do Boundary conditions for Conductive Heat Flux:
 ! IAXIS faces:
 X1AXIS = IAXIS
-DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-   IW = CC_REGFACE_IAXIS_Z(IFACE)%IWC
+DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+   IW = CC_REGFACE_IAXIS_Z(IFACE)%IWC; IF(IW<1) CYCLE
    IF( WALL(IW)%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
        WALL(IW)%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
        WALL(IW)%BOUNDARY_TYPE==PERIODIC_BOUNDARY ) CYCLE ! Already done on previous loops.
@@ -11043,8 +11043,8 @@ SUBROUTINE CC_CONDUCTION_HEAT_FLUX
 
 ! JAXIS faces:
 X1AXIS = JAXIS
-DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-   IW = CC_REGFACE_JAXIS_Z(IFACE)%IWC
+DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+   IW = CC_REGFACE_JAXIS_Z(IFACE)%IWC; IF(IW<1) CYCLE
    IF( WALL(IW)%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
        WALL(IW)%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
        WALL(IW)%BOUNDARY_TYPE==PERIODIC_BOUNDARY ) CYCLE ! Already done on previous loops.
@@ -11061,8 +11061,8 @@ SUBROUTINE CC_CONDUCTION_HEAT_FLUX
 
 ! KAXIS faces:
 X1AXIS = KAXIS
-DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-   IW = CC_REGFACE_KAXIS_Z(IFACE)%IWC
+DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+   IW = CC_REGFACE_KAXIS_Z(IFACE)%IWC; IF(IW<1) CYCLE
    IF( WALL(IW)%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
        WALL(IW)%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
        WALL(IW)%BOUNDARY_TYPE==PERIODIC_BOUNDARY ) CYCLE ! Already done on previous loops.
@@ -11117,8 +11117,8 @@ SUBROUTINE CC_CONDUCTION_HEAT_FLUX
 ENDDO
 
 ! GASPHASE cut-faces:
-DO ICF = 1,MESHES(NM)%N_BBCUTFACE_MESH
-   IF ( CUT_FACE(ICF)%STATUS /= CC_GASPHASE ) CYCLE
+DO ICF = 1,MESHES(NM)%N_CUTFACE_MESH
+   IF ( CUT_FACE(ICF)%STATUS /= CC_GASPHASE .OR. MESHES(NM)%CUT_FACE(ICF)%IWC<1) CYCLE
    IW = MESHES(NM)%CUT_FACE(ICF)%IWC
    IF( WALL(IW)%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
        WALL(IW)%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
@@ -11527,8 +11527,8 @@ SUBROUTINE CC_DIFFUSIVE_MASS_FLUXES(NM)
    ! done in WALL_LOOP_2 of DIVERGENCE_PART_1 (divg.f90):
    ! IAXIS faces:
    X1AXIS = IAXIS
-   DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-      IW = CC_REGFACE_IAXIS_Z(IFACE)%IWC
+   DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+      IW = CC_REGFACE_IAXIS_Z(IFACE)%IWC; IF(IW<1) CYCLE
       IF( WALL(IW)%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
           WALL(IW)%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
           WALL(IW)%BOUNDARY_TYPE==PERIODIC_BOUNDARY ) CYCLE ! Already done on previous loops.
@@ -11539,8 +11539,8 @@ SUBROUTINE CC_DIFFUSIVE_MASS_FLUXES(NM)
 
    ! JAXIS faces:
    X1AXIS = JAXIS
-   DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-      IW = CC_REGFACE_JAXIS_Z(IFACE)%IWC
+   DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+      IW = CC_REGFACE_JAXIS_Z(IFACE)%IWC; IF(IW<1) CYCLE
       IF( WALL(IW)%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
           WALL(IW)%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
           WALL(IW)%BOUNDARY_TYPE==PERIODIC_BOUNDARY ) CYCLE
@@ -11551,8 +11551,8 @@ SUBROUTINE CC_DIFFUSIVE_MASS_FLUXES(NM)
 
    ! KAXIS faces:
    X1AXIS = KAXIS
-   DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-      IW = CC_REGFACE_KAXIS_Z(IFACE)%IWC
+   DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+      IW = CC_REGFACE_KAXIS_Z(IFACE)%IWC; IF(IW<1) CYCLE
       IF( WALL(IW)%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
           WALL(IW)%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
           WALL(IW)%BOUNDARY_TYPE==PERIODIC_BOUNDARY ) CYCLE
@@ -11579,8 +11579,8 @@ SUBROUTINE CC_DIFFUSIVE_MASS_FLUXES(NM)
 
    ! GASPHASE cut-faces:
    ! In case of Cut Faces and OPEN boundaries redefine the location of the guard cells with atmospheric conditions:
-   DO ICF = 1,MESHES(NM)%N_BBCUTFACE_MESH
-      IF ( CUT_FACE(ICF)%STATUS /= CC_GASPHASE ) CYCLE
+   DO ICF = 1,MESHES(NM)%N_CUTFACE_MESH
+      IF ( CUT_FACE(ICF)%STATUS /= CC_GASPHASE .OR. MESHES(NM)%CUT_FACE(ICF)%IWC<1) CYCLE
       IW = MESHES(NM)%CUT_FACE(ICF)%IWC
       IF( WALL(IW)%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
           WALL(IW)%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
@@ -13630,9 +13630,8 @@ SUBROUTINE GET_ADVDIFFVECTOR_SCALAR_3D(N)
    ! IAXIS faces:
    X1AXIS = IAXIS
    REGFACE_Z => CC_REGFACE_IAXIS_Z
-   DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-      IW = REGFACE_Z(IFACE)%IWC
-      WC=>WALL(IW)
+   DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+      IW = REGFACE_Z(IFACE)%IWC; IF(IW<1) CYCLE; WC => WALL(IW)
       IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
           WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
           WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY ) CYCLE
@@ -13653,9 +13652,8 @@ SUBROUTINE GET_ADVDIFFVECTOR_SCALAR_3D(N)
    ! JAXIS faces:
    X1AXIS = JAXIS
    REGFACE_Z => CC_REGFACE_JAXIS_Z
-   DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-      IW = REGFACE_Z(IFACE)%IWC
-      WC=>WALL(IW)
+   DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+      IW = REGFACE_Z(IFACE)%IWC; IF(IW<1) CYCLE; WC => WALL(IW)
       IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
           WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
           WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY ) CYCLE
@@ -13676,9 +13674,8 @@ SUBROUTINE GET_ADVDIFFVECTOR_SCALAR_3D(N)
    ! KAXIS faces:
    X1AXIS = KAXIS
    REGFACE_Z => CC_REGFACE_KAXIS_Z
-   DO IFACE=1,MESHES(NM)%CC_NBBREGFACE_Z(X1AXIS)
-      IW = REGFACE_Z(IFACE)%IWC
-      WC=>WALL(IW)
+   DO IFACE=1,MESHES(NM)%CC_NREGFACE_Z(X1AXIS)
+      IW = REGFACE_Z(IFACE)%IWC; IF(IW<1) CYCLE; WC => WALL(IW)
       IF (WC%BOUNDARY_TYPE==NULL_BOUNDARY         .OR. &
           WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. &
           WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY ) CYCLE
@@ -21565,7 +21562,7 @@ SUBROUTINE GET_GASPHASE_REGRCFACES_DATA
 LOGICAL :: INLIST
 
 INTEGER :: IW,JJ,KK,IIF,JJF,KKF,IOR,LOWHIGH_TEST,IIG,JJG,KKG
-INTEGER :: IBNDINT,IC
+INTEGER :: IBNDINT,IC,IC2
 LOGICAL :: FLGIN
 INTEGER, PARAMETER :: OZPOS=0, ICPOS=1, JCPOS=2, IFPOS=3
 INTEGER :: IERR
@@ -21687,6 +21684,15 @@ SUBROUTINE GET_GASPHASE_REGRCFACES_DATA
                   IF ( CCVAR(I+1,J,K,CC_UNKZ)>0 ) LOHIBUFF(HIGH_IND,IREG)= .TRUE.
                   IC = CELL_INDEX(I,J,K)
                   IJKBUFFER(IAXIS:KAXIS+1,IREG) = (/ I, J, K, CELL(IC)%WALL_INDEX( X1AXIS) /)
+                  IC2 =  CELL_INDEX(I+1,J,K)
+                  IF(CELL(IC)%WALL_INDEX( X1AXIS)>0 .AND. CELL(IC2)%WALL_INDEX(-X1AXIS)>0) THEN
+                     IREG = IREG + 1
+                     LOHIBUFF(LOW_IND,IREG) = LOHIBUFF(LOW_IND,IREG-1)
+                     LOHIBUFF(HIGH_IND,IREG)= LOHIBUFF(HIGH_IND,IREG-1)
+                     IJKBUFFER(IAXIS:KAXIS+1,IREG) = (/ I, J, K, CELL(IC2)%WALL_INDEX(-X1AXIS) /)
+                  ELSEIF(CELL(IC2)%WALL_INDEX(-X1AXIS)>0) THEN
+                     IJKBUFFER(IAXIS:KAXIS+1,IREG) = (/ I, J, K, CELL(IC2)%WALL_INDEX(-X1AXIS) /)
+                  ENDIF
                ENDDO
             ENDDO
          ENDDO
@@ -21760,6 +21766,15 @@ SUBROUTINE GET_GASPHASE_REGRCFACES_DATA
                   IF ( CCVAR(I,J+1,K,CC_UNKZ)>0 ) LOHIBUFF(HIGH_IND,IREG)= .TRUE.
                   IC   = CELL_INDEX(I,J,K)
                   IJKBUFFER(IAXIS:KAXIS+1,IREG) = (/ I, J, K, CELL(IC)%WALL_INDEX( X1AXIS) /)
+                  IC2 =  CELL_INDEX(I,J+1,K)
+                  IF(CELL(IC)%WALL_INDEX( X1AXIS)>0 .AND. CELL(IC2)%WALL_INDEX(-X1AXIS)>0) THEN
+                     IREG = IREG + 1
+                     LOHIBUFF(LOW_IND,IREG) = LOHIBUFF(LOW_IND,IREG-1)
+                     LOHIBUFF(HIGH_IND,IREG)= LOHIBUFF(HIGH_IND,IREG-1)
+                     IJKBUFFER(IAXIS:KAXIS+1,IREG) = (/ I, J, K, CELL(IC2)%WALL_INDEX(-X1AXIS) /)
+                  ELSEIF(CELL(IC2)%WALL_INDEX(-X1AXIS)>0) THEN
+                     IJKBUFFER(IAXIS:KAXIS+1,IREG) = (/ I, J, K, CELL(IC2)%WALL_INDEX(-X1AXIS) /)
+                  ENDIF
                ENDDO
             ENDDO
          ENDDO
@@ -21833,6 +21848,15 @@ SUBROUTINE GET_GASPHASE_REGRCFACES_DATA
                   IF ( CCVAR(I,J,K+1,CC_UNKZ)>0 ) LOHIBUFF(HIGH_IND,IREG)= .TRUE.
                   IC   = CELL_INDEX(I,J,K)
                   IJKBUFFER(IAXIS:KAXIS+1,IREG) = (/ I, J, K, CELL(IC)%WALL_INDEX( X1AXIS) /)
+                  IC2 =  CELL_INDEX(I,J,K+1)
+                  IF(CELL(IC)%WALL_INDEX( X1AXIS)>0 .AND. CELL(IC2)%WALL_INDEX(-X1AXIS)>0) THEN
+                     IREG = IREG + 1
+                     LOHIBUFF(LOW_IND,IREG) = LOHIBUFF(LOW_IND,IREG-1)
+                     LOHIBUFF(HIGH_IND,IREG)= LOHIBUFF(HIGH_IND,IREG-1)
+                     IJKBUFFER(IAXIS:KAXIS+1,IREG) = (/ I, J, K, CELL(IC2)%WALL_INDEX(-X1AXIS) /)
+                  ELSEIF(CELL(IC2)%WALL_INDEX(-X1AXIS)>0) THEN
+                     IJKBUFFER(IAXIS:KAXIS+1,IREG) = (/ I, J, K, CELL(IC2)%WALL_INDEX(-X1AXIS) /)
+                  ENDIF
                ENDDO
             ENDDO
          ENDDO

Source/divg.f90

@@ -771,7 +771,7 @@ SUBROUTINE DIVERGENCE_PART_1(T,DT,NM)
    ENDDO WALL_LOOP4
 
 
-   CFACE_LOOP: DO ICF=N_EXTERNAL_CFACE_CELLS+1,N_EXTERNAL_CFACE_CELLS+N_INTERNAL_CFACE_CELLS
+   CFACE_LOOP: 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)

Source/hvac.f90

@@ -2201,7 +2201,7 @@ SUBROUTINE HVAC_BC_IN(NM)
    CALL INITIALIZE_HVAC
 ENDDO WALL_LOOP
 
-CFACE_LOOP: DO ICF=N_EXTERNAL_CFACE_CELLS+1,N_EXTERNAL_CFACE_CELLS+N_INTERNAL_CFACE_CELLS
+CFACE_LOOP: DO ICF=INTERNAL_CFACE_CELLS_LB+1,INTERNAL_CFACE_CELLS_LB+N_INTERNAL_CFACE_CELLS
    CFA => CFACE(ICF)
    B1 => BOUNDARY_PROP1(CFA%B1_INDEX)
    BC => BOUNDARY_COORD(CFA%BC_INDEX)

Source/init.f90

@@ -426,7 +426,7 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
 
 ! Over-ride default ambient gas species mass fractions if defined in INIT
 
-M%IWORK1=0 
+M%IWORK1=0
 DO N=1,N_INIT
    IN => INITIALIZATION(N)
    IF ((IN%NODE_ID/='null')) CYCLE
@@ -699,7 +699,7 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_1(DT,NM)
       IF (OB%J1==OB%J2 .AND. ABS(OB%Y2-OB%Y1)>TWO_EPSILON_EB) THEN
          DO K=OB%K1+1,OB%K2
             IF (OB%I1>0) THEN
-               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1 
+               M%N_THIN_WALL_CELLS = M%N_THIN_WALL_CELLS + 1
                CALL INIT_THIN_WALL_CELL(NM,OB%I1,OB%J1,K,N,M%N_THIN_WALL_CELLS,-1,OB%SURF_INDEX(-1),3)
             ENDIF
             IF (OB%I2<M%IBAR) THEN
@@ -1126,7 +1126,7 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_2(NM)
 
 ENDDO WALL_LOOP_0
 
-CFACE_LOOP_0: DO ICF=1,M%N_EXTERNAL_CFACE_CELLS+M%N_INTERNAL_CFACE_CELLS
+CFACE_LOOP_0: DO ICF=1,M%N_EXTERNAL_CFACE_CELLS+M%N_INTWALL_CFACE_CELLS+M%N_INTERNAL_CFACE_CELLS
    CFA => M%CFACE(ICF)
    BC  => M%BOUNDARY_COORD(CFA%BC_INDEX)
    B1  => M%BOUNDARY_PROP1(CFA%B1_INDEX)
@@ -1921,7 +1921,7 @@ SUBROUTINE ADJUST_HT3D_ONE_D
    SELECT CASE(ABS(BC%IOR))
       CASE(1) ; ONE_D%LAYER_THICKNESS(1) = OB%X2 - OB%X1
       CASE(2) ; ONE_D%LAYER_THICKNESS(1) = OB%Y2 - OB%Y1
-      CASE(3) ; ONE_D%LAYER_THICKNESS(1) = OB%Z2 - OB%Z1 
+      CASE(3) ; ONE_D%LAYER_THICKNESS(1) = OB%Z2 - OB%Z1
    END SELECT
 ENDIF
 
@@ -2001,7 +2001,7 @@ SUBROUTINE INITIALIZE_HT3D_WALL_CELLS(NM)
 
 M => MESHES(NM)
 
-! Loop over all 3-D wall cells, and for each interior node, find wall or thin wall cells in the other two "alternate" 
+! Loop over all 3-D wall cells, and for each interior node, find wall or thin wall cells in the other two "alternate"
 ! coordinate directions
 
 PRIMARY_WALL_LOOP: DO IW=1,M%N_EXTERNAL_WALL_CELLS+M%N_INTERNAL_WALL_CELLS
@@ -2035,7 +2035,7 @@ SUBROUTINE INITIALIZE_HT3D_WALL_CELLS(NM)
 
    PRIMARY_NODE_LOOP: DO I=1,NWP
 
-      SELECT CASE(BC%IOR)  
+      SELECT CASE(BC%IOR)
          CASE( 1) ; X1=BC%X2-ONE_D%X(I)   ; X2=BC%X2-ONE_D%X(I-1)
          CASE(-1) ; X1=BC%X1+ONE_D%X(I-1) ; X2=BC%X1+ONE_D%X(I)
          CASE( 2) ; Y1=BC%Y2-ONE_D%X(I)   ; Y2=BC%Y2-ONE_D%X(I-1)
@@ -2135,7 +2135,7 @@ SUBROUTINE INITIALIZE_HT3D_WALL_CELLS(NM)
 
    PRIMARY_THIN_NODE_LOOP: DO I=1,NWP
 
-      SELECT CASE(BC%IOR)  
+      SELECT CASE(BC%IOR)
          CASE( 1) ; X1=BC%X2-ONE_D%X(I)   ; X2=BC%X2-ONE_D%X(I-1)
          CASE(-1) ; X1=BC%X1+ONE_D%X(I-1) ; X2=BC%X1+ONE_D%X(I)
          CASE( 2) ; Y1=BC%Y2-ONE_D%X(I)   ; Y2=BC%Y2-ONE_D%X(I-1)
@@ -2236,17 +2236,17 @@ SUBROUTINE SEARCH_FOR_ALTERNATE_WALL_CELLS(NOM,WALL_INDEX,THIN_WALL_INDEX)
 
 ALTERNATE_NODE_LOOP: DO I2=1,NWP2  ! Loop over nodes of alternate wall cell
    SELECT CASE(BC2%IOR)
-      CASE( 1) ; XX1=BC2%X2-ONE_D2%X(I2)   ; XX2=BC2%X2-ONE_D2%X(I2-1) 
+      CASE( 1) ; XX1=BC2%X2-ONE_D2%X(I2)   ; XX2=BC2%X2-ONE_D2%X(I2-1)
                  DXX = MIN(XX2,X2)-MAX(XX1,X1) ; IF (DXX<=0._EB) CYCLE ALTERNATE_NODE_LOOP
-      CASE(-1) ; XX1=BC2%X1+ONE_D2%X(I2-1) ; XX2=BC2%X1+ONE_D2%X(I2) 
+      CASE(-1) ; XX1=BC2%X1+ONE_D2%X(I2-1) ; XX2=BC2%X1+ONE_D2%X(I2)
                  DXX = MIN(XX2,X2)-MAX(XX1,X1) ; IF (DXX<=0._EB) CYCLE ALTERNATE_NODE_LOOP
-      CASE( 2) ; YY1=BC2%Y2-ONE_D2%X(I2)   ; YY2=BC2%Y2-ONE_D2%X(I2-1) 
+      CASE( 2) ; YY1=BC2%Y2-ONE_D2%X(I2)   ; YY2=BC2%Y2-ONE_D2%X(I2-1)
                  DYY = MIN(YY2,Y2)-MAX(YY1,Y1) ; IF (DYY<=0._EB) CYCLE ALTERNATE_NODE_LOOP
-      CASE(-2) ; YY1=BC2%Y1+ONE_D2%X(I2-1) ; YY2=BC2%Y1+ONE_D2%X(I2) 
+      CASE(-2) ; YY1=BC2%Y1+ONE_D2%X(I2-1) ; YY2=BC2%Y1+ONE_D2%X(I2)
                  DYY = MIN(YY2,Y2)-MAX(YY1,Y1) ; IF (DYY<=0._EB) CYCLE ALTERNATE_NODE_LOOP
-      CASE( 3) ; ZZ1=BC2%Z2-ONE_D2%X(I2)   ; ZZ2=BC2%Z2-ONE_D2%X(I2-1) 
+      CASE( 3) ; ZZ1=BC2%Z2-ONE_D2%X(I2)   ; ZZ2=BC2%Z2-ONE_D2%X(I2-1)
                  DZZ = MIN(ZZ2,Z2)-MAX(ZZ1,Z1) ; IF (DZZ<=0._EB) CYCLE ALTERNATE_NODE_LOOP
-      CASE(-3) ; ZZ1=BC2%Z1+ONE_D2%X(I2-1) ; ZZ2=BC2%Z1+ONE_D2%X(I2) 
+      CASE(-3) ; ZZ1=BC2%Z1+ONE_D2%X(I2-1) ; ZZ2=BC2%Z1+ONE_D2%X(I2)
                  DZZ = MIN(ZZ2,Z2)-MAX(ZZ1,Z1) ; IF (DZZ<=0._EB) CYCLE ALTERNATE_NODE_LOOP
    END SELECT
    OVERLAP_VOLUME = DXX*DYY*DZZ
@@ -2309,11 +2309,11 @@ SUBROUTINE INITIALIZE_MESH_VARIABLES_3(NM)
 M => MESHES(NM)
 
 IF (M%N_CFACE_CELLS_DIM<1) RETURN
-   
+
 N_REALS    = 0
 N_INTEGERS = 0
 N_LOGICALS = 0
-DO ICF=1,M%N_EXTERNAL_CFACE_CELLS+M%N_INTERNAL_CFACE_CELLS
+DO ICF=1,M%N_EXTERNAL_CFACE_CELLS+M%N_INTWALL_CFACE_CELLS+M%N_INTERNAL_CFACE_CELLS
    CFA => M%CFACE(ICF)
    N_REALS    = MAX(N_REALS   ,CFA%N_REALS)
    N_INTEGERS = MAX(N_INTEGERS,CFA%N_INTEGERS)
@@ -3354,7 +3354,7 @@ SUBROUTINE INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,TT)
 
    B1%AREA_ADJUST = SF%AREA_MULTIPLIER * VT%INPUT_AREA/VT%FDS_AREA
    IF (B1%AREA_ADJUST<=TWO_EPSILON_EB) B1%AREA_ADJUST = 1._EB
-   
+
    IF (VT%CTRL_INDEX > 0) THEN
       IF (.NOT. CONTROL(VT%CTRL_INDEX)%CURRENT_STATE) T_ACTIVATE=1E10_EB
    ENDIF