FireX: Merge with firemodels/master

Build/makefile

@@ -469,7 +469,7 @@ soot.o : prec.o mesh.o cons.o func.o turb.o
 pois.o : prec.o cons.o
 geom.o : gsmv.o func.o prec.o cons.o type.o mesh.o turb.o imkl.o data.o
 ccib.o : gsmv.o func.o prec.o cons.o type.o mesh.o turb.o geom.o
-radi.o : func.o prec.o cons.o type.o mesh.o rcal.o geom.o
+radi.o : func.o prec.o cons.o type.o mesh.o rcal.o geom.o ccib.o
 part.o : func.o prec.o cons.o type.o mesh.o devc.o data.o geom.o soot.o ccib.o
 vege.o : func.o prec.o cons.o type.o mesh.o geom.o
 ctrl.o : prec.o cons.o type.o mesh.o func.o
@@ -504,7 +504,7 @@ soot.obj : prec.obj mesh.obj cons.obj func.obj turb.obj
 pois.obj : prec.obj cons.obj
 geom.obj : gsmv.obj func.obj prec.obj cons.obj type.obj mesh.obj turb.obj imkl.obj data.obj
 ccib.obj : gsmv.obj func.obj prec.obj cons.obj type.obj mesh.obj turb.obj geom.obj
-radi.obj : func.obj prec.obj cons.obj type.obj mesh.obj rcal.obj geom.obj
+radi.obj : func.obj prec.obj cons.obj type.obj mesh.obj rcal.obj geom.obj ccib.obj
 part.obj : func.obj prec.obj cons.obj type.obj mesh.obj devc.obj data.obj geom.obj soot.obj ccib.obj
 vege.obj : func.obj prec.obj cons.obj type.obj mesh.obj geom.obj
 ctrl.obj : prec.obj cons.obj type.obj mesh.obj func.obj

Source/ccib.f90

@@ -132,7 +132,8 @@ MODULE CC_SCALARS
           GET_CFACE_OPEN_BC_COEF,GET_FN_DIVERGENCE_CUTCELL,GET_OPENBC_TANGENTIAL_CUTFACE_VEL,&
           GET_CUTCELL_DDDT,GET_H_CUTFACES,GET_H_MATRIX_CC,GET_H_GUARD_CUTCELL,GET_CRTCFCC_INT_STENCILS,GET_RCFACES_H, &
           GET_CC_MATRIXGRAPH_H,GET_CC_IROW,GET_CC_UNKH,GET_CUTCELL_HP, GET_LINKED_FV, GET_PRES_CFACE_BCS, &
-          GET_PRES_CFACE, GET_PRES_CFACE_TEST, GET_UVWGAS_CFACE, GET_MUDNS_CFACE, GET_BOUNDFACE_GEOM_INFO_H, &
+          GET_PRES_CFACE, GET_PRES_CFACE_TEST, GET_UVWGAS_CFACE, GET_CFACE_RAD_NVEC, GET_MUDNS_CFACE, &
+          GET_BOUNDFACE_GEOM_INFO_H, &
           GET_FH_FROM_PRHS_AND_BCS,GET_LINKED_VELOCITIES, &
           FINISH_CC, INIT_CUTCELL_DATA,MESH_CC_EXCHANGE,NUMBER_UNKH_CUTCELLS,&
           ROTATED_CUBE_ANN_SOLN,ROTATED_CUBE_VELOCITY_FLUX,ROTATED_CUBE_RHS_ZZ,&
@@ -154,11 +155,11 @@ MODULE CC_SCALARS
 SUBROUTINE INITIALIZE_BACK_CFACE_EXCHANGE
 
 INTEGER :: NM,NOM,IQ,NQ,IIO,JJO,KKO,NOM2,ICFACE,SURF_B,FRONT,OD_INDEX,N_RESOLVED,N_GLOBAL,IP,K,P,TOTAL_SEND,TOTAL_RECV,IERR
-REAL(EB) :: POS(3)
+REAL(EB) :: POS(3),NBACK(3)
 INTEGER,  ALLOCATABLE :: SEND_CNT(:),RECV_CNT(:),SDISPL(:),RDISPL(:),CNT3_S(:),CNT3_R(:),DSP3_S(:),DSP3_R(:),COUNTER(:)
 INTEGER,  ALLOCATABLE :: SEND_NOM(:),RECV_NOM(:),REPLY_IDX(:),REPLY_SURF(:),REPLY_DIM(:), &
                          REPLY_BACK(:),REPLY_BSURF(:),REPLY_BDIM(:),MAP_NM(:),MAP_FRONT(:)
-REAL(EB), ALLOCATABLE :: SEND_POS(:),RECV_POS(:)
+REAL(EB), ALLOCATABLE :: SEND_POS(:),RECV_POS(:),SEND_NBACK(:),RECV_NBACK(:)
 
 N_RESOLVED = 0
 ALLOCATE(SEND_CNT(0:N_MPI_PROCESSES-1)) ; SEND_CNT = 0
@@ -169,10 +170,11 @@ SUBROUTINE INITIALIZE_BACK_CFACE_EXCHANGE
       IF (NQ==0) CYCLE
       IF (PROCESS(NOM)==MY_RANK) THEN  ! Owner is on this process: resolve directly.
          DO IQ=1,NQ
-            POS(1:3) = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_XYZ(IAXIS:KAXIS,IQ)
+            POS(1:3)   = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_XYZ(IAXIS:KAXIS,IQ)
+            NBACK(1:3) = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_NBACK(IAXIS:KAXIS,IQ)
             CALL SEARCH_OTHER_MESHES(POS(1),POS(2),POS(3),NOM2,IIO,JJO,KKO)
             IF (NOM2/=NOM) CYCLE
-            ICFACE = GET_BACK_CFACE_INDEX(NOM,IIO,JJO,KKO,POS,-1,-1)
+            ICFACE = GET_BACK_CFACE_INDEX(NOM,IIO,JJO,KKO,POS,-1,-1,NBACK)
             IF (ICFACE>0) THEN
                FRONT    = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_FRONT(IQ)
                OD_INDEX = MESHES(NM)%CFACE(FRONT)%OD_INDEX
@@ -206,10 +208,12 @@ SUBROUTINE INITIALIZE_BACK_CFACE_EXCHANGE
    TOTAL_SEND = SDISPL(N_MPI_PROCESSES-1) + SEND_CNT(N_MPI_PROCESSES-1)
    TOTAL_RECV = RDISPL(N_MPI_PROCESSES-1) + RECV_CNT(N_MPI_PROCESSES-1)
    CNT3_S = 3*SEND_CNT ; CNT3_R = 3*RECV_CNT ; DSP3_S = 3*SDISPL ; DSP3_R = 3*RDISPL
-   ALLOCATE(SEND_NOM(MAX(TOTAL_SEND,1)),SEND_POS(MAX(3*TOTAL_SEND,1)),MAP_NM(MAX(TOTAL_SEND,1)), &
+   ALLOCATE(SEND_NOM(MAX(TOTAL_SEND,1)),SEND_POS(MAX(3*TOTAL_SEND,1)),SEND_NBACK(MAX(3*TOTAL_SEND,1)), &
+            MAP_NM(MAX(TOTAL_SEND,1)), &
             MAP_FRONT(MAX(TOTAL_SEND,1)),REPLY_BACK(MAX(TOTAL_SEND,1)),REPLY_BSURF(MAX(TOTAL_SEND,1)), &
             REPLY_BDIM(MAX(TOTAL_SEND,1)))
-   ALLOCATE(RECV_NOM(MAX(TOTAL_RECV,1)),RECV_POS(MAX(3*TOTAL_RECV,1)),REPLY_IDX(MAX(TOTAL_RECV,1)), &
+   ALLOCATE(RECV_NOM(MAX(TOTAL_RECV,1)),RECV_POS(MAX(3*TOTAL_RECV,1)),RECV_NBACK(MAX(3*TOTAL_RECV,1)), &
+            REPLY_IDX(MAX(TOTAL_RECV,1)), &
             REPLY_SURF(MAX(TOTAL_RECV,1)),REPLY_DIM(MAX(TOTAL_RECV,1)))
    ALLOCATE(COUNTER(0:N_MPI_PROCESSES-1)) ; COUNTER = SDISPL
 
@@ -224,24 +228,28 @@ SUBROUTINE INITIALIZE_BACK_CFACE_EXCHANGE
          IP = PROCESS(NOM)
          DO IQ=1,NQ
             P = COUNTER(IP) + 1 ; COUNTER(IP) = P
-            SEND_NOM(P)         = NOM
-            SEND_POS(3*P-2:3*P) = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_XYZ(IAXIS:KAXIS,IQ)
-            MAP_NM(P)           = NM
-            MAP_FRONT(P)        = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_FRONT(IQ)
+            SEND_NOM(P)           = NOM
+            SEND_POS(3*P-2:3*P)   = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_XYZ(IAXIS:KAXIS,IQ)
+            SEND_NBACK(3*P-2:3*P) = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_NBACK(IAXIS:KAXIS,IQ)
+            MAP_NM(P)             = NM
+            MAP_FRONT(P)          = MESHES(NM)%OMESH(NOM)%CFACE_QUERY_FRONT(IQ)
          ENDDO
       ENDDO
    ENDDO
    CALL MPI_ALLTOALLV(SEND_NOM,SEND_CNT,SDISPL,MPI_INTEGER,RECV_NOM,RECV_CNT,RDISPL,MPI_INTEGER,MPI_COMM_WORLD,IERR)
    CALL MPI_ALLTOALLV(SEND_POS,CNT3_S,DSP3_S,MPI_DOUBLE_PRECISION,RECV_POS,CNT3_R,DSP3_R,MPI_DOUBLE_PRECISION,&
                       MPI_COMM_WORLD,IERR)
+   CALL MPI_ALLTOALLV(SEND_NBACK,CNT3_S,DSP3_S,MPI_DOUBLE_PRECISION,RECV_NBACK,CNT3_R,DSP3_R,MPI_DOUBLE_PRECISION,&
+                      MPI_COMM_WORLD,IERR)
 
    ! Owner side: resolve each received query against its local cut-face data.
    DO K=1,TOTAL_RECV
       NOM = RECV_NOM(K)
-      POS(1:3) = RECV_POS(3*K-2:3*K)
+      POS(1:3)   = RECV_POS(3*K-2:3*K)
+      NBACK(1:3) = RECV_NBACK(3*K-2:3*K)
       CALL SEARCH_OTHER_MESHES(POS(1),POS(2),POS(3),NOM2,IIO,JJO,KKO)
       ICFACE = 0 ; SURF_B = 0
-      IF (NOM2==NOM) ICFACE = GET_BACK_CFACE_INDEX(NOM,IIO,JJO,KKO,POS,-1,-1)
+      IF (NOM2==NOM) ICFACE = GET_BACK_CFACE_INDEX(NOM,IIO,JJO,KKO,POS,-1,-1,NBACK)
       IF (ICFACE>0) SURF_B = MESHES(NOM)%CFACE(ICFACE)%SURF_INDEX
       REPLY_IDX(K)  = ICFACE
       REPLY_SURF(K) = SURF_B
@@ -304,8 +312,8 @@ SUBROUTINE INITIALIZE_BACK_CFACE_EXCHANGE
    CALL MPI_ALLREDUCE(MPI_IN_PLACE,N_GLOBAL,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,IERR)
    IF (N_GLOBAL>0) CF_BACK_EXCHANGE_ACTIVE = .TRUE.
 
-   DEALLOCATE(RECV_CNT,SDISPL,RDISPL,CNT3_S,CNT3_R,DSP3_S,DSP3_R,COUNTER,SEND_NOM,SEND_POS,MAP_NM,MAP_FRONT, &
-              REPLY_BACK,REPLY_BSURF,REPLY_BDIM,RECV_NOM,RECV_POS,REPLY_IDX,REPLY_SURF,REPLY_DIM)
+   DEALLOCATE(RECV_CNT,SDISPL,RDISPL,CNT3_S,CNT3_R,DSP3_S,DSP3_R,COUNTER,SEND_NOM,SEND_POS,SEND_NBACK,MAP_NM,MAP_FRONT, &
+              REPLY_BACK,REPLY_BSURF,REPLY_BDIM,RECV_NOM,RECV_POS,RECV_NBACK,REPLY_IDX,REPLY_SURF,REPLY_DIM)
 
 ENDIF CROSS_PROCESS_IF
 DEALLOCATE(SEND_CNT)
@@ -622,7 +630,7 @@ SUBROUTINE GET_H_GUARD_CUTCELL(IPZ,HP)
                   ELSE;                  D_EXT = ABS(MESHES(NOM)%ZC(KKO) - XFN)
                   ENDIF
                END SELECT
-               ! Refinement seam only: match GET_H_MATRIX_CC sub-face selection (skip solid/no-unknown partners). 
+               ! Refinement seam only: match GET_H_MATRIX_CC sub-face selection (skip solid/no-unknown partners).
                IF (WC%BOUNDARY_TYPE == INTERPOLATED_BOUNDARY) THEN
                   IF (ALLOCATED(OM%MUNKH)) THEN; IF (OM%MUNKH(IIO,JJO,KKO) < 1) CYCLE; ENDIF
                   IF (MESHES(NOM)%FCVAR(IIO+ISHF(IAXIS),JJO+ISHF(JAXIS),KKO+ISHF(KAXIS),CC_FGSC,X1AXIS) == CC_SOLID) CYCLE
@@ -6756,6 +6764,20 @@ SUBROUTINE GET_UVWGAS_CFACE(U_CELL,V_CELL,W_CELL,IND1,IND2,UP,VP,WP,PREDFCT)
 RETURN
 END SUBROUTINE GET_UVWGAS_CFACE
 
+! ----------------------------------- GET_CFACE_RAD_NVEC ---------------------------------
+FUNCTION GET_CFACE_RAD_NVEC(CFA) RESULT(NVEC)
+
+TYPE(CFACE_TYPE), INTENT(IN) :: CFA
+REAL(EB) :: NVEC(IAXIS:KAXIS)
+INTEGER :: IBOD,IWSEL
+
+NVEC(IAXIS:KAXIS) = BOUNDARY_COORD(CFA%BC_INDEX)%NVEC(IAXIS:KAXIS)
+IBOD = CUT_FACE(CFA%CUT_FACE_IND1)%BODTRI(1,CFA%CUT_FACE_IND2)
+IWSEL= CUT_FACE(CFA%CUT_FACE_IND1)%BODTRI(2,CFA%CUT_FACE_IND2)
+IF (IBOD>=1 .AND. IWSEL>=1) NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL)
+
+END FUNCTION GET_CFACE_RAD_NVEC
+
 ! ----------------------------------- GET_PRES_CFACE ---------------------------------
 
 SUBROUTINE GET_PRES_CFACE(PRESS,IND1,IND2,CFA)
@@ -8936,7 +8958,7 @@ SUBROUTINE CCSPECIES_ADVECTION
 ! region and adds components to thermodynamic divergence.
 
 ! Local Variables:
-REAL(EB) :: RHO_Z_PV(-2:1), VELC, FN_ZZ, ZZ_GET_N
+REAL(EB) :: RHO_Z_PV(-2:1), VELC, FN_ZZ, FN_ZZ_USE, ZZ_GET_N
 REAL(EB), PARAMETER :: SGNFCT=1._EB
 INTEGER :: IOR, ICFA
 TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
@@ -9127,6 +9149,12 @@ SUBROUTINE CCSPECIES_ADVECTION
          CASE(INTERPOLATED_BOUNDARY)
             VELC = UVW_SAVE(IW)
       END SELECT
+      ! For a cut-cell against a mesh interface,flux-limited value built across the interface.
+      IF ((WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY) .AND. &
+          RCF%CELL_LIST(1,ISIDE+2)==CC_FTYPE_CFGAS) THEN
+         ICC = RCF%CELL_LIST(2,ISIDE+2); JCC = RCF%CELL_LIST(3,ISIDE+2)
+         CALL GET_CCIFACE_RHOZZ_LIM(X1AXIS,I,J,K,ISIDE+2,ICC,JCC,N,VELC,PRFCT,FN_ZZ)
+      ENDIF
       SELECT CASE(RCF%CELL_LIST(1,ISIDE+2))
       CASE(CC_FTYPE_RGGAS)
         SELECT CASE(X1AXIS)
@@ -9342,8 +9370,12 @@ SUBROUTINE CCSPECIES_ADVECTION
            RHOPV(ISIDE) = PRFCT *CUT_CELL(ICC)%RHO(JCC) + (1._EB-PRFCT) *CUT_CELL(ICC)%RHOS(JCC)
            ZZ_GET_N     = PRFCT*CUT_CELL(ICC)%ZZ(N,JCC) + (1._EB-PRFCT)*CUT_CELL(ICC)%ZZS(N,JCC)
            RHO_Z_PV(ISIDE) = RHOPV(ISIDE)*ZZ_GET_N
+           ! At a mesh interface use the flux-limited bar{rho*Z} across the interface.
+           FN_ZZ_USE = FN_ZZ
+           IF (WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY) &
+              CALL GET_CCIFACE_RHOZZ_LIM(X1AXIS,I,J,K,ISIDE+2,ICC,JCC,N,VELC,PRFCT,FN_ZZ_USE)
            CUT_CELL(ICC)%U_DOT_DEL_RHO_Z_VOL(N,JCC) = CUT_CELL(ICC)%U_DOT_DEL_RHO_Z_VOL(N,JCC) + &
-                                                      FCT*SGNFCT*(FN_ZZ-RHO_Z_PV(ISIDE))*VELC*AF
+                                                      FCT*SGNFCT*(FN_ZZ_USE-RHO_Z_PV(ISIDE))*VELC*AF
          ENDIF
       ENDDO ! IFACE
 
@@ -9473,7 +9505,7 @@ SUBROUTINE CCENTHALPY_ADVECTION
 
 
 ! Local Variables:
-REAL(EB) :: RHO_H_S_PV(-2:1), VELC, VELC2, FN_H_S, TMP_F_GAS
+REAL(EB) :: RHO_H_S_PV(-2:1), VELC, VELC2, FN_H_S, FN_H_S_USE, TMP_F_GAS
 INTEGER  :: IOR, ICFA
 TYPE(BOUNDARY_PROP1_TYPE), POINTER :: B1
 TYPE(CC_REGFACEZ_TYPE),  POINTER, DIMENSION(:) :: REGFACE_Z
@@ -10070,6 +10102,12 @@ SUBROUTINE CCENTHALPY_ADVECTION
    ZZ_GET(1:N_TRACKED_SPECIES) = B1%ZZ_F(1:N_TRACKED_SPECIES)
    CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,H_S,TMP_F_GAS)
    FN_H_S = B1%RHO_F*H_S ! bar{rho*hs}
+   ! At a mesh interface with a cut-cell, use the flux-limited bar{rho*h_s} built across the interface.
+   IF ((WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY) .AND. &
+       RC_FACE(IFACE)%CELL_LIST(1,ISIDE+2)==CC_FTYPE_CFGAS) THEN
+      ICC = RC_FACE(IFACE)%CELL_LIST(2,ISIDE+2); JCC = RC_FACE(IFACE)%CELL_LIST(3,ISIDE+2)
+      CALL GET_CCIFACE_RHOHS_LIM(X1AXIS,I,J,K,ISIDE+2,ICC,JCC,VELC,PRFCT,FN_H_S)
+   ENDIF
    ! Finally add to Div:
    SELECT CASE(RC_FACE(IFACE)%CELL_LIST(1,ISIDE+2))
    CASE(CC_FTYPE_RGGAS) ! Regular cell
@@ -10133,7 +10171,13 @@ SUBROUTINE CCENTHALPY_ADVECTION
                                  (1._EB-PRFCT)*CUT_CELL(ICC)%ZZS(1:N_TRACKED_SPECIES,JCC)
          CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,H_S,TMPV(ISIDE))
          RHO_H_S_PV(ISIDE) = RHOPV(ISIDE)*H_S
-         CUT_CELL(ICC)%DVOL(JCC)=CUT_CELL(ICC)%DVOL(JCC)+FCT*(-1._EB)*(FN_H_S-RHO_H_S_PV(ISIDE))*VELC*AF ! +ve or -ve dot
+         ! At a mesh interface use the flux-limited bar{rho*h_s} built across the interface.
+         FN_H_S_USE = FN_H_S
+         IF (WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY) &
+            CALL GET_CCIFACE_RHOHS_LIM(CUT_FACE(ICF)%IJK(KAXIS+1), &
+                 CUT_FACE(ICF)%IJK(IAXIS),CUT_FACE(ICF)%IJK(JAXIS),CUT_FACE(ICF)%IJK(KAXIS), &
+                 ISIDE+2,ICC,JCC,VELC,PRFCT,FN_H_S_USE)
+         CUT_CELL(ICC)%DVOL(JCC)=CUT_CELL(ICC)%DVOL(JCC)+FCT*(-1._EB)*(FN_H_S_USE-RHO_H_S_PV(ISIDE))*VELC*AF ! +ve or -ve dot
       END SELECT
    ENDDO ! IFACE
 ENDDO ! ICF
@@ -13449,6 +13493,112 @@ SUBROUTINE GET_RHOZZVECTOR_SCALAR_3D(N)
 END SUBROUTINE GET_RHOZZVECTOR_SCALAR_3D
 
 
+! ----------------------------- GET_CCIFACE_RHOZZ_LIM --------------------------
+! Flux-limited bar{rho*Z_N} face value at an INTERPOLATED / PERIODIC-boundary gasphase cut-face.
+SUBROUTINE GET_CCIFACE_RHOZZ_LIM(X1AXIS,IFC,JFC,KFC,ILOC,ICC,JCC,N,VELC,PRFCT,FN_ZZ)
+
+INTEGER, INTENT(IN) :: X1AXIS,IFC,JFC,KFC,ILOC,ICC,JCC,N
+REAL(EB), INTENT(IN) :: VELC,PRFCT
+REAL(EB), INTENT(OUT):: FN_ZZ
+INTEGER :: DD,ILSD,ILOUT,II,JJ,KK,SHFT(3)
+REAL(EB) :: RZPV(-2:1)
+REAL(EB), POINTER, DIMENSION(:,:,:)   :: RHOP
+REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZP
+REAL(EB), TARGET  :: ZT(0:3,1,1),UT(1,1,1),FT(1,1,1)
+REAL(EB), POINTER, DIMENSION(:,:,:) :: ZTP,UTP,FTP
+
+IF (PRFCT > 0.5_EB) THEN; RHOP=>RHO ; ZZP=>ZZ   ! PRFCT=1 -> n-level fields
+ELSE                    ; RHOP=>RHOS; ZZP=>ZZS  ! PRFCT=0 -> star fields
+ENDIF
+
+SELECT CASE(X1AXIS)
+CASE(IAXIS); SHFT=(/1,0,0/)
+CASE(JAXIS); SHFT=(/0,1,0/)
+CASE(KAXIS); SHFT=(/0,0,1/)
+END SELECT
+
+! Cartesian rho*Z on the four stencil cells: cell(DD) = base_face_cell + (1+DD) along the face axis.
+DO DD=-2,1
+   II=IFC+(1+DD)*SHFT(1); JJ=JFC+(1+DD)*SHFT(2); KK=KFC+(1+DD)*SHFT(3)
+   RZPV(DD)=RHOP(II,JJ,KK)*ZZP(II,JJ,KK,N)
+ENDDO
+
+! Override the local immediate cut-cell (side ILOC) with the detailed cut-cell rho*Z:
+ILSD = ILOC-2  ! LOW_IND(1) -> -1 ; HIGH_IND(2) -> 0
+RZPV(ILSD) = (PRFCT*CUT_CELL(ICC)%RHO(JCC) +(1._EB-PRFCT)*CUT_CELL(ICC)%RHOS(JCC)) * &
+             (PRFCT*CUT_CELL(ICC)%ZZ(N,JCC)+(1._EB-PRFCT)*CUT_CELL(ICC)%ZZS(N,JCC))
+
+! Solid fall-back only for the local-outer stencil cell (the ghost-outer cell is remote gas):
+IF (ILOC==LOW_IND) THEN; ILOUT=-2; ELSE; ILOUT=1; ENDIF
+II=IFC+(1+ILOUT)*SHFT(1); JJ=JFC+(1+ILOUT)*SHFT(2); KK=KFC+(1+ILOUT)*SHFT(3)
+IF (II>=1.AND.II<=IBAR .AND. JJ>=1.AND.JJ<=JBAR .AND. KK>=1.AND.KK<=KBAR) THEN
+   IF (CELL(CELL_INDEX(II,JJ,KK))%SOLID .OR. CCVAR(II,JJ,KK,CC_CGSC)==CC_SOLID) RZPV(ILOUT)=RZPV(ILSD)
+ENDIF
+ZT(0:3,1,1)=RZPV(-2:1); UT(1,1,1)=VELC
+ZTP=>ZT; UTP=>UT; FTP=>FT
+CALL GET_SCALAR_FACE_VALUE(UTP,ZTP,FTP,1,1,1,1,1,1,1,I_FLUX_LIMITER)
+FN_ZZ=FT(1,1,1)
+
+RETURN
+END SUBROUTINE GET_CCIFACE_RHOZZ_LIM
+
+
+! ----------------------------- GET_CCIFACE_RHOHS_LIM --------------------------
+! Enthalpy counterpart of GET_CCIFACE_RHOZZ_LIM: flux-limited bar{rho*h_s} face value at an
+! INTERPOLATED / PERIODIC-boundary gasphase cut-face.
+SUBROUTINE GET_CCIFACE_RHOHS_LIM(X1AXIS,IFC,JFC,KFC,ILOC,ICC,JCC,VELC,PRFCT,FN_H_S)
+USE PHYSICAL_FUNCTIONS, ONLY: GET_SENSIBLE_ENTHALPY
+
+INTEGER, INTENT(IN) :: X1AXIS,IFC,JFC,KFC,ILOC,ICC,JCC
+REAL(EB), INTENT(IN) :: VELC,PRFCT
+REAL(EB), INTENT(OUT):: FN_H_S
+INTEGER :: DD,ILSD,ILOUT,II,JJ,KK,SHFT(3)
+REAL(EB) :: RHS_PV(-2:1),HS1,ZZ_GET(1:N_TRACKED_SPECIES)
+REAL(EB), POINTER, DIMENSION(:,:,:)   :: RHOP
+REAL(EB), POINTER, DIMENSION(:,:,:,:) :: ZZP
+REAL(EB), TARGET  :: ZT(0:3,1,1),UT(1,1,1),FT(1,1,1)
+REAL(EB), POINTER, DIMENSION(:,:,:) :: ZTP,UTP,FTP
+
+IF (PRFCT > 0.5_EB) THEN; RHOP=>RHO ; ZZP=>ZZ
+ELSE                    ; RHOP=>RHOS; ZZP=>ZZS
+ENDIF
+
+SELECT CASE(X1AXIS)
+CASE(IAXIS); SHFT=(/1,0,0/)
+CASE(JAXIS); SHFT=(/0,1,0/)
+CASE(KAXIS); SHFT=(/0,0,1/)
+END SELECT
+
+! Cartesian rho*h_s on the four stencil cells: cell(DD) = base_face_cell + (1+DD) along the axis.
+DO DD=-2,1
+   II=IFC+(1+DD)*SHFT(1); JJ=JFC+(1+DD)*SHFT(2); KK=KFC+(1+DD)*SHFT(3)
+   ZZ_GET(1:N_TRACKED_SPECIES)=ZZP(II,JJ,KK,1:N_TRACKED_SPECIES)
+   CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,HS1,TMP(II,JJ,KK))
+   RHS_PV(DD)=RHOP(II,JJ,KK)*HS1
+ENDDO
+
+! Override the local immediate cut-cell (side ILOC) with the detailed cut-cell rho*h_s:
+ILSD = ILOC-2
+ZZ_GET(1:N_TRACKED_SPECIES)=PRFCT*CUT_CELL(ICC)%ZZ(1:N_TRACKED_SPECIES,JCC)+ &
+                            (1._EB-PRFCT)*CUT_CELL(ICC)%ZZS(1:N_TRACKED_SPECIES,JCC)
+CALL GET_SENSIBLE_ENTHALPY(ZZ_GET,HS1,CUT_CELL(ICC)%TMP(JCC))
+RHS_PV(ILSD)=(PRFCT*CUT_CELL(ICC)%RHO(JCC)+(1._EB-PRFCT)*CUT_CELL(ICC)%RHOS(JCC))*HS1
+
+! Solid fall-back for the local-outer stencil cell only:
+IF (ILOC==LOW_IND) THEN; ILOUT=-2; ELSE; ILOUT=1; ENDIF
+II=IFC+(1+ILOUT)*SHFT(1); JJ=JFC+(1+ILOUT)*SHFT(2); KK=KFC+(1+ILOUT)*SHFT(3)
+IF (II>=1.AND.II<=IBAR .AND. JJ>=1.AND.JJ<=JBAR .AND. KK>=1.AND.KK<=KBAR) THEN
+   IF (CELL(CELL_INDEX(II,JJ,KK))%SOLID .OR. CCVAR(II,JJ,KK,CC_CGSC)==CC_SOLID) RHS_PV(ILOUT)=RHS_PV(ILSD)
+ENDIF
+ZT(0:3,1,1)=RHS_PV(-2:1); UT(1,1,1)=VELC
+ZTP=>ZT; UTP=>UT; FTP=>FT
+CALL GET_SCALAR_FACE_VALUE(UTP,ZTP,FTP,1,1,1,1,1,1,1,I_FLUX_LIMITER)
+FN_H_S=FT(1,1,1)
+
+RETURN
+END SUBROUTINE GET_CCIFACE_RHOHS_LIM
+
+
 ! ------------------------- GET_ADVDIFFVECTOR_SCALAR_3D -------------------------
 
 SUBROUTINE GET_ADVDIFFVECTOR_SCALAR_3D(N)
@@ -13460,7 +13610,7 @@ SUBROUTINE GET_ADVDIFFVECTOR_SCALAR_3D(N)
 REAL(EB):: PRFCT
 INTEGER :: X1AXIS,IFACE,IND(LOW_IND:HIGH_IND),IND_LOC(LOW_IND:HIGH_IND),ICF,IND1,IND2,IOR
 INTEGER :: LOCROW_1,LOCROW_2,ILOC,IROW,ICC,JCC,IW
-REAL(EB):: AF,VELC,RHO_Z,FN_ZZ,FCT
+REAL(EB):: AF,VELC,RHO_Z,FN_ZZ,FCT,FN_ZZ_USE
 REAL(EB), POINTER, DIMENSION(:,:,:)  :: RHOP,UP,VP,WP
 REAL(EB), POINTER, DIMENSION(:,:,:,:)::  ZZP
 TYPE(CC_REGFACEZ_TYPE),  POINTER, DIMENSION(:) :: REGFACE_Z
@@ -13620,6 +13770,12 @@ SUBROUTINE GET_ADVDIFFVECTOR_SCALAR_3D(N)
          CASE(INTERPOLATED_BOUNDARY)
             VELC = UVW_SAVE(IW)
       END SELECT
+      ! For a cut-cell against a mesh interface , use the flux-limited bar{rho*Z} built across the interface.
+      IF ((WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY) .AND. &
+          MESHES(NM)%RC_FACE(IFACE)%CELL_LIST(1,ILOC)==CC_FTYPE_CFGAS) THEN
+         ICC = MESHES(NM)%RC_FACE(IFACE)%CELL_LIST(2,ILOC); JCC = MESHES(NM)%RC_FACE(IFACE)%CELL_LIST(3,ILOC)
+         CALL GET_CCIFACE_RHOZZ_LIM(X1AXIS,I,J,K,ILOC,ICC,JCC,N,VELC,PRFCT,FN_ZZ)
+      ENDIF
 
       IROW=IND_LOC(ILOC)   ! Process Local Unknown number.
       FCT = REAL(3-2*ILOC,EB)
@@ -13652,7 +13808,15 @@ SUBROUTINE GET_ADVDIFFVECTOR_SCALAR_3D(N)
          ! First (rho hs)_i,j,k:
          IF (CUT_FACE(ICF)%CELL_LIST(1,ILOC,IFACE) == CC_FTYPE_CFGAS) THEN
             IROW=IND_LOC(ILOC)   ! Process Local Unknown number.
-            F_Z(IROW) = F_Z(IROW) + FCT*FN_ZZ*VELC*AF
+            FN_ZZ_USE = FN_ZZ
+            ! Build the flux-limited bar{rho*Z} across the interface.
+            IF (WC%BOUNDARY_TYPE==INTERPOLATED_BOUNDARY .OR. WC%BOUNDARY_TYPE==PERIODIC_BOUNDARY) THEN
+               ICC = CUT_FACE(ICF)%CELL_LIST(2,ILOC,IFACE); JCC = CUT_FACE(ICF)%CELL_LIST(3,ILOC,IFACE)
+               CALL GET_CCIFACE_RHOZZ_LIM(CUT_FACE(ICF)%IJK(KAXIS+1), &
+                    CUT_FACE(ICF)%IJK(IAXIS),CUT_FACE(ICF)%IJK(JAXIS),CUT_FACE(ICF)%IJK(KAXIS), &
+                    ILOC,ICC,JCC,N,VELC,PRFCT,FN_ZZ_USE)
+            ENDIF
+            F_Z(IROW) = F_Z(IROW) + FCT*FN_ZZ_USE*VELC*AF
          ENDIF
       ENDDO ! IFACE