Merge pull request #16406 from ericvmueller/master

FDS Source: bug fix to properly sum particle drag FVX/Y/Z_D values at…

Author: ericvmueller

Source/main.f90

@@ -737,7 +737,7 @@ PROGRAM FDS
 
       COMPUTE_WALL_BC_LOOP_A: DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
          CALL WALL_BC(T,DT,NM)
-         IF (PARTICLE_DRAG) CALL PARTICLE_MOMENTUM_TRANSFER(DT,NM)
+         IF (PARTICLE_DRAG) CALL PARTICLE_MOMENTUM_TRANSFER(NM)
          CALL DIVERGENCE_PART_1(T,DT,NM)
       ENDDO COMPUTE_WALL_BC_LOOP_A
 
@@ -887,9 +887,15 @@ PROGRAM FDS
    COMPUTE_WALL_BC_2A: DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
       CALL PARTICLE_MASS_ENERGY_TRANSFER(T,DT,NM)
       CALL MOVE_PARTICLES(T,DT,NM)
-      IF (PARTICLE_DRAG) CALL PARTICLE_MOMENTUM_TRANSFER(DT,NM)
    ENDDO COMPUTE_WALL_BC_2A
 
+   IF (PARTICLE_DRAG) THEN
+      IF (NMESHES>1) CALL MESH_EXCHANGE(20)
+      DO NM=LOWER_MESH_INDEX,UPPER_MESH_INDEX
+         CALL PARTICLE_MOMENTUM_TRANSFER(NM)
+      ENDDO
+   ENDIF
+
    ! Exchange the number of particles sent from mesh to mesh (7), and if non-zero, exchange particles (11)
 
    IF (MPI_PARTICLE_EXCHANGE) THEN
@@ -2086,6 +2092,14 @@ SUBROUTINE INITIALIZE_MESH_EXCHANGE_1(NM)
    OM%FVX = 0._EB
    OM%FVY = 0._EB
    OM%FVZ = 0._EB
+   IF (PARTICLE_DRAG) THEN
+      ALLOCATE(OM%FVX_D(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
+      ALLOCATE(OM%FVY_D(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
+      ALLOCATE(OM%FVZ_D(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
+      OM%FVX_D = 0._EB
+      OM%FVY_D = 0._EB
+      OM%FVZ_D = 0._EB
+   ENDIF
    ALLOCATE(OM%KRES(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
    OM%KRES = 0._EB
    ALLOCATE(OM%Q(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX))
@@ -2814,6 +2828,7 @@ SUBROUTINE POST_RECEIVES(CODE)
             IJK_SIZE = (M3%I_MAX_R-M3%I_MIN_R+1)*(M3%J_MAX_R-M3%J_MIN_R+1)*(M3%K_MAX_R-M3%K_MIN_R+1)
             ALLOCATE(M3%REAL_RECV_PKG1(M3%NIC_R*(6+2*N_TOTAL_SCALARS)))
             ALLOCATE(M3%REAL_RECV_PKG3(IJK_SIZE*4))
+            ALLOCATE(M3%REAL_RECV_PKG4(M3%NIC_R))
             ALLOCATE(M3%REAL_RECV_PKG7(M3%NIC_R*3))
             ALLOCATE(M3%REAL_RECV_PKG8(M3%NIC_R*2))
 
@@ -2846,6 +2861,10 @@ SUBROUTINE POST_RECEIVES(CODE)
             CALL MPI_RECV_INIT(M3%REAL_RECV_PKG3(1),SIZE(M3%REAL_RECV_PKG3),MPI_DOUBLE_PRECISION,SNODE,NOM,MPI_COMM_WORLD,&
                                REQ3(N_REQ3),IERR)
 
+            N_REQ4 = N_REQ4 + 1
+            CALL MPI_RECV_INIT(M3%REAL_RECV_PKG4(1),SIZE(M3%REAL_RECV_PKG4),MPI_DOUBLE_PRECISION,SNODE,NOM,MPI_COMM_WORLD,&
+                               REQ4(N_REQ4),IERR)
+
             N_REQ7 = N_REQ7 + 1
             CALL MPI_RECV_INIT(M3%REAL_RECV_PKG7(1),SIZE(M3%REAL_RECV_PKG7),MPI_DOUBLE_PRECISION,SNODE,NOM,MPI_COMM_WORLD,&
                                REQ7(N_REQ7),IERR)
@@ -3001,6 +3020,7 @@ SUBROUTINE MESH_EXCHANGE(CODE)
             IJK_SIZE = (M3%I_MAX_S-M3%I_MIN_S+1)*(M3%J_MAX_S-M3%J_MIN_S+1)*(M3%K_MAX_S-M3%K_MIN_S+1)
             ALLOCATE(M3%REAL_SEND_PKG1(M3%NIC_S*(6+2*N_TOTAL_SCALARS)))
             ALLOCATE(M3%REAL_SEND_PKG3(IJK_SIZE*4))
+            ALLOCATE(M3%REAL_SEND_PKG4(M3%NIC_S))
             ALLOCATE(M3%REAL_SEND_PKG7(M3%NIC_S*3))
 
             IF (LEVEL_SET_MODE>0 .OR. TERRAIN_CASE) ALLOCATE(M3%REAL_SEND_PKG14(5*M3%NIC_S))
@@ -3036,6 +3056,10 @@ SUBROUTINE MESH_EXCHANGE(CODE)
             CALL MPI_SEND_INIT(M3%REAL_SEND_PKG3(1),SIZE(M3%REAL_SEND_PKG3),MPI_DOUBLE_PRECISION,RNODE,NM,MPI_COMM_WORLD,&
                                REQ3(N_REQ3),IERR)
 
+            N_REQ4 = N_REQ4 + 1
+            CALL MPI_SEND_INIT(M3%REAL_SEND_PKG4(1),SIZE(M3%REAL_SEND_PKG4),MPI_DOUBLE_PRECISION,RNODE,NM,MPI_COMM_WORLD,&
+                               REQ4(N_REQ4),IERR)
+
             N_REQ7 = N_REQ7 + 1
             CALL MPI_SEND_INIT(M3%REAL_SEND_PKG7(1),SIZE(M3%REAL_SEND_PKG7),MPI_DOUBLE_PRECISION,RNODE,NM,MPI_COMM_WORLD,&
                                REQ7(N_REQ7),IERR)
@@ -3217,6 +3241,28 @@ SUBROUTINE MESH_EXCHANGE(CODE)
          ENDIF
       ENDIF
 
+      ! Exchange particle drag at interpolated boundaries
+
+      IF (CODE==20 .AND. M3%NIC_S>0) THEN
+         IF (RNODE/=SNODE) THEN
+            PACK_REAL_SEND_PKG4: DO LL=1,M3%NIC_S
+               SELECT CASE(M3%IOR_S(LL))
+                  CASE(-1) ; M3%REAL_SEND_PKG4(LL) = M%FVX_D(M3%IIO_S(LL)-1,M3%JJO_S(LL)  ,M3%KKO_S(LL)  )
+                  CASE( 1) ; M3%REAL_SEND_PKG4(LL) = M%FVX_D(M3%IIO_S(LL)  ,M3%JJO_S(LL)  ,M3%KKO_S(LL)  )
+                  CASE(-2) ; M3%REAL_SEND_PKG4(LL) = M%FVY_D(M3%IIO_S(LL)  ,M3%JJO_S(LL)-1,M3%KKO_S(LL)  )
+                  CASE( 2) ; M3%REAL_SEND_PKG4(LL) = M%FVY_D(M3%IIO_S(LL)  ,M3%JJO_S(LL)  ,M3%KKO_S(LL)  )
+                  CASE(-3) ; M3%REAL_SEND_PKG4(LL) = M%FVZ_D(M3%IIO_S(LL)  ,M3%JJO_S(LL)  ,M3%KKO_S(LL)-1)
+                  CASE( 3) ; M3%REAL_SEND_PKG4(LL) = M%FVZ_D(M3%IIO_S(LL)  ,M3%JJO_S(LL)  ,M3%KKO_S(LL)  )
+               END SELECT
+            ENDDO PACK_REAL_SEND_PKG4
+         ELSE
+            M2=>MESHES(NOM)%OMESH(NM)
+            M2%FVX_D(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX) = M%FVX_D(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
+            M2%FVY_D(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX) = M%FVY_D(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
+            M2%FVZ_D(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX) = M%FVZ_D(IMIN:IMAX,JMIN:JMAX,KMIN:KMAX)
+         ENDIF
+      ENDIF
+
       ! Send pressure information at the end of the PREDICTOR (CODE=3) or CORRECTOR (CODE=6) stage of the time step
 
       IF ((CODE==3.OR.CODE==6) .AND. M3%NIC_S>0) THEN
@@ -3445,6 +3491,11 @@ SUBROUTINE MESH_EXCHANGE(CODE)
    CALL TIMEOUT('MPI exchange of pressure (CODE=5)',N_REQ7,REQ7(1:N_REQ7))
 ENDIF
 
+IF (N_MPI_PROCESSES>1 .AND. CODE==20 .AND. N_REQ4>0) THEN
+   CALL MPI_STARTALL(N_REQ4,REQ4(1:N_REQ4),IERR)
+   CALL TIMEOUT('MPI exchange of particle drag (CODE=20)',N_REQ4,REQ4(1:N_REQ4))
+ENDIF
+
 IF (N_MPI_PROCESSES>1 .AND. CODE==6 .AND. N_REQ6>0) THEN
    CALL MPI_STARTALL(N_REQ6,REQ6(1:N_REQ6),IERR)
    CALL TIMEOUT('MPI exchange of back wall info (CODE=6)',N_REQ6,REQ6(1:N_REQ6))
@@ -3547,6 +3598,21 @@ SUBROUTINE MESH_EXCHANGE(CODE)
          ENDDO UNPACK_REAL_RECV_PKG7
       ENDIF
 
+      ! Unpack particle drag at interpolated boundaries
+
+      IF (CODE==20 .AND. M2%NIC_R>0 .AND. RNODE/=SNODE) THEN
+         UNPACK_REAL_RECV_PKG4: DO LL=1,M2%NIC_R
+            SELECT CASE(M2%IOR_R(LL))
+               CASE(-1) ; M2%FVX_D(M2%IIO_R(LL)-1,M2%JJO_R(LL)  ,M2%KKO_R(LL)  ) = M2%REAL_RECV_PKG4(LL)
+               CASE( 1) ; M2%FVX_D(M2%IIO_R(LL)  ,M2%JJO_R(LL)  ,M2%KKO_R(LL)  ) = M2%REAL_RECV_PKG4(LL)
+               CASE(-2) ; M2%FVY_D(M2%IIO_R(LL)  ,M2%JJO_R(LL)-1,M2%KKO_R(LL)  ) = M2%REAL_RECV_PKG4(LL)
+               CASE( 2) ; M2%FVY_D(M2%IIO_R(LL)  ,M2%JJO_R(LL)  ,M2%KKO_R(LL)  ) = M2%REAL_RECV_PKG4(LL)
+               CASE(-3) ; M2%FVZ_D(M2%IIO_R(LL)  ,M2%JJO_R(LL)  ,M2%KKO_R(LL)-1) = M2%REAL_RECV_PKG4(LL)
+               CASE( 3) ; M2%FVZ_D(M2%IIO_R(LL)  ,M2%JJO_R(LL)  ,M2%KKO_R(LL)  ) = M2%REAL_RECV_PKG4(LL)
+            END SELECT
+         ENDDO UNPACK_REAL_RECV_PKG4
+      ENDIF
+
       ! Unpack pressure following PREDICTOR stage of time step
 
       IF ((CODE==3.OR.CODE==6) .AND. M2%NIC_R>0 .AND. RNODE/=SNODE) THEN

Source/part.f90

@@ -1855,10 +1855,6 @@ SUBROUTINE MOVE_PARTICLES(T,DT,NM)
 
 PART_UVWMAX = 0._EB
 
-IF (MESHES(NM)%NLP==0) RETURN
-
-! Set the CPU timer and point to the current mesh variables
-
 TNOW=CURRENT_TIME()
 CALL POINT_TO_MESH(NM)
 
@@ -1868,6 +1864,8 @@ SUBROUTINE MOVE_PARTICLES(T,DT,NM)
    FVZ_D = 0._EB
 ENDIF
 
+IF (MESHES(NM)%NLP==0) RETURN
+
 IF (CC_IBM) CALL CUTFACE_VELOCITIES(NM,U,V,W,CUTFACES=.TRUE.)
 
 ! Loop through all Lagrangian particles and move them one time step
@@ -2564,6 +2562,8 @@ SUBROUTINE MOVE_PARTICLES(T,DT,NM)
 
 IF (CC_IBM) CALL CUTFACE_VELOCITIES(NM,U,V,W,CUTFACES=.FALSE.)
 
+IF (PARTICLE_DRAG) CALL CLIP_PARTICLE_DRAG(DT,NM)
+
 ! Remove out-of-bounds particles
 
 CALL REMOVE_PARTICLES(T,NM)
@@ -4570,21 +4570,18 @@ END SUBROUTINE PARTICLE_RUNNING_AVERAGES
 END SUBROUTINE PARTICLE_MASS_ENERGY_TRANSFER
 
 
-!> \brief Add PARTICLE momentum as a force term in momentum equation
-!> \param DT Current time step (s)
+!> \brief Clip particle drag fluxes to plus/minus abs(u)/dt
+!> \param DT Time step size (s)
 !> \param NM Current mesh number
 
-SUBROUTINE PARTICLE_MOMENTUM_TRANSFER(DT,NM)
+SUBROUTINE CLIP_PARTICLE_DRAG(DT,NM)
 
-USE CC_SCALARS, ONLY: CUTFACE_VELOCITIES
 INTEGER, INTENT(IN) :: NM
 REAL(EB), INTENT(IN) :: DT
 REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW
 REAL(EB) :: RDT,UODT,VODT,WODT
 INTEGER :: I,J,K
 
-IF (MESHES(NM)%NLP==0) RETURN
-
 CALL POINT_TO_MESH(NM)
 
 RDT = 1._EB/DT
@@ -4599,25 +4596,172 @@ SUBROUTINE PARTICLE_MOMENTUM_TRANSFER(DT,NM)
    WW => WS
 ENDIF
 
-IF (CC_IBM) CALL CUTFACE_VELOCITIES(NM,UU,VV,WW,CUTFACES=.TRUE.)
-
-! Add summed particle accelerations to the momentum equation. Limit the value to plus/minus abs(u)/dt to prevent a sudden
-! change in gas direction.
+! Limit drag flux to plus/minus abs(u)/dt to prevent a sudden change in gas direction.
 
 DO K=0,KBAR
    DO J=0,JBAR
       DO I=0,IBAR
          UODT = ABS(UU(I,J,K)*RDT)
          VODT = ABS(VV(I,J,K)*RDT)
          WODT = ABS(WW(I,J,K)*RDT)
-         FVX(I,J,K) = FVX(I,J,K) + MIN(UODT,MAX(-UODT,FVX_D(I,J,K)))
-         FVY(I,J,K) = FVY(I,J,K) + MIN(VODT,MAX(-VODT,FVY_D(I,J,K)))
-         FVZ(I,J,K) = FVZ(I,J,K) + MIN(WODT,MAX(-WODT,FVZ_D(I,J,K)))
+         FVX_D(I,J,K) = MIN(UODT,MAX(-UODT,FVX_D(I,J,K)))
+         FVY_D(I,J,K) = MIN(VODT,MAX(-VODT,FVY_D(I,J,K)))
+         FVZ_D(I,J,K) = MIN(WODT,MAX(-WODT,FVZ_D(I,J,K)))
       ENDDO
    ENDDO
 ENDDO
 
-IF (CC_IBM) CALL CUTFACE_VELOCITIES(NM,UU,VV,WW,CUTFACES=.FALSE.)
+END SUBROUTINE CLIP_PARTICLE_DRAG
+
+
+!> \brief Add PARTICLE momentum as a force term in momentum equation
+!> \param DT Current time step (s)
+!> \param NM Current mesh number
+
+SUBROUTINE PARTICLE_MOMENTUM_TRANSFER(NM)
+
+INTEGER, INTENT(IN) :: NM
+INTEGER :: NOM,II,JJ,KK,IOR,IW,IIO,JJO,KKO
+REAL(EB) :: DA_OTHER,FVX_OTHER,FVY_OTHER,FVZ_OTHER
+TYPE (OMESH_TYPE), POINTER :: OM
+TYPE (MESH_TYPE), POINTER :: M2
+TYPE(WALL_TYPE), POINTER :: WC
+TYPE(EXTERNAL_WALL_TYPE), POINTER :: EWC
+TYPE(BOUNDARY_COORD_TYPE), POINTER :: BC
+
+CALL POINT_TO_MESH(NM)
+
+! Add area-weighted particle drag from neighboring meshes at interpolated boundaries.
+
+IF (NMESHES>1 .AND. CORRECTOR) THEN
+
+   EXTERNAL_WALL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS
+
+      WC=>WALL(IW)
+      EWC=>EXTERNAL_WALL(IW)
+      IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) CYCLE EXTERNAL_WALL_LOOP
+
+      BC => BOUNDARY_COORD(WC%BC_INDEX)
+      II  = BC%II
+      JJ  = BC%JJ
+      KK  = BC%KK
+      IOR = BC%IOR
+      NOM = EWC%NOM
+      OM => OMESH(NOM)
+      M2 => MESHES(NOM)
+
+      DA_OTHER = 0._EB
+
+      SELECT CASE(ABS(IOR))
+         CASE(1)
+            DO KKO=EWC%KKO_MIN,EWC%KKO_MAX
+               DO JJO=EWC%JJO_MIN,EWC%JJO_MAX
+                  DO IIO=EWC%IIO_MIN,EWC%IIO_MAX
+                     DA_OTHER = DA_OTHER + M2%DY(JJO)*M2%DZ(KKO)
+                  ENDDO
+               ENDDO
+            ENDDO
+         CASE(2)
+            DO KKO=EWC%KKO_MIN,EWC%KKO_MAX
+               DO JJO=EWC%JJO_MIN,EWC%JJO_MAX
+                  DO IIO=EWC%IIO_MIN,EWC%IIO_MAX
+                     DA_OTHER = DA_OTHER + M2%DX(IIO)*M2%DZ(KKO)
+                  ENDDO
+               ENDDO
+            ENDDO
+         CASE(3)
+            DO KKO=EWC%KKO_MIN,EWC%KKO_MAX
+               DO JJO=EWC%JJO_MIN,EWC%JJO_MAX
+                  DO IIO=EWC%IIO_MIN,EWC%IIO_MAX
+                     DA_OTHER = DA_OTHER + M2%DX(IIO)*M2%DY(JJO)
+                  ENDDO
+               ENDDO
+            ENDDO
+      END SELECT
+
+      SELECT CASE(IOR)
+
+         CASE( 1)
+
+            FVX_OTHER = 0._EB
+            DO KKO=EWC%KKO_MIN,EWC%KKO_MAX
+               DO JJO=EWC%JJO_MIN,EWC%JJO_MAX
+                  DO IIO=EWC%IIO_MIN,EWC%IIO_MAX
+                     FVX_OTHER = FVX_OTHER + OM%FVX_D(IIO,JJO,KKO)*M2%DY(JJO)*M2%DZ(KKO)/DA_OTHER
+                  ENDDO
+               ENDDO
+            ENDDO
+            FVX_D(0,JJ,KK) = FVX_D(0,JJ,KK) + FVX_OTHER
+
+         CASE(-1)
+
+            FVX_OTHER = 0._EB
+            DO KKO=EWC%KKO_MIN,EWC%KKO_MAX
+               DO JJO=EWC%JJO_MIN,EWC%JJO_MAX
+                  DO IIO=EWC%IIO_MIN,EWC%IIO_MAX
+                     FVX_OTHER = FVX_OTHER + OM%FVX_D(IIO-1,JJO,KKO)*M2%DY(JJO)*M2%DZ(KKO)/DA_OTHER
+                  ENDDO
+               ENDDO
+            ENDDO
+            FVX_D(IBAR,JJ,KK) = FVX_D(IBAR,JJ,KK) + FVX_OTHER
+
+         CASE( 2)
+
+            FVY_OTHER = 0._EB
+            DO KKO=EWC%KKO_MIN,EWC%KKO_MAX
+               DO JJO=EWC%JJO_MIN,EWC%JJO_MAX
+                  DO IIO=EWC%IIO_MIN,EWC%IIO_MAX
+                     FVY_OTHER = FVY_OTHER + OM%FVY_D(IIO,JJO,KKO)*M2%DX(IIO)*M2%DZ(KKO)/DA_OTHER
+                  ENDDO
+               ENDDO
+            ENDDO
+            FVY_D(II,0,KK) = FVY_D(II,0,KK) + FVY_OTHER
+
+         CASE(-2)
+
+            FVY_OTHER = 0._EB
+            DO KKO=EWC%KKO_MIN,EWC%KKO_MAX
+               DO JJO=EWC%JJO_MIN,EWC%JJO_MAX
+                  DO IIO=EWC%IIO_MIN,EWC%IIO_MAX
+                     FVY_OTHER = FVY_OTHER + OM%FVY_D(IIO,JJO-1,KKO)*M2%DX(IIO)*M2%DZ(KKO)/DA_OTHER
+                  ENDDO
+               ENDDO
+            ENDDO
+            FVY_D(II,JBAR,KK) = FVY_D(II,JBAR,KK) + FVY_OTHER
+
+         CASE( 3)
+
+            FVZ_OTHER = 0._EB
+            DO KKO=EWC%KKO_MIN,EWC%KKO_MAX
+               DO JJO=EWC%JJO_MIN,EWC%JJO_MAX
+                  DO IIO=EWC%IIO_MIN,EWC%IIO_MAX
+                     FVZ_OTHER = FVZ_OTHER + OM%FVZ_D(IIO,JJO,KKO)*M2%DX(IIO)*M2%DY(JJO)/DA_OTHER
+                  ENDDO
+               ENDDO
+            ENDDO
+            FVZ_D(II,JJ,0) = FVZ_D(II,JJ,0) + FVZ_OTHER
+
+         CASE(-3)
+
+            FVZ_OTHER = 0._EB
+            DO KKO=EWC%KKO_MIN,EWC%KKO_MAX
+               DO JJO=EWC%JJO_MIN,EWC%JJO_MAX
+                  DO IIO=EWC%IIO_MIN,EWC%IIO_MAX
+                     FVZ_OTHER = FVZ_OTHER + OM%FVZ_D(IIO,JJO,KKO-1)*M2%DX(IIO)*M2%DY(JJO)/DA_OTHER
+                  ENDDO
+               ENDDO
+            ENDDO
+            FVZ_D(II,JJ,KBAR) = FVZ_D(II,JJ,KBAR) + FVZ_OTHER
+
+      END SELECT
+
+   ENDDO EXTERNAL_WALL_LOOP
+
+ENDIF
+
+FVX = FVX + FVX_D
+FVY = FVY + FVY_D
+FVZ = FVZ + FVZ_D
 
 END SUBROUTINE PARTICLE_MOMENTUM_TRANSFER
 

Source/type.f90

@@ -1025,7 +1025,8 @@ MODULE TYPES
 TYPE (SURFACE_TYPE), DIMENSION(:), ALLOCATABLE, TARGET :: SURFACE
 
 TYPE OMESH_TYPE
-   REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) :: MU,RHO,RHOS,TMP,U,V,W,US,VS,WS,H,HS,FVX,FVY,FVZ,D,DS,KRES,IL_S,IL_R,IL_R_OLD,Q
+   REAL(EB), ALLOCATABLE, DIMENSION(:,:,:) ::&
+         MU,RHO,RHOS,TMP,U,V,W,US,VS,WS,H,HS,FVX,FVY,FVZ,FVX_D,FVY_D,FVZ_D,D,DS,KRES,IL_S,IL_R,IL_R_OLD,Q
    REAL(EB), ALLOCATABLE, DIMENSION(:,:,:,:) :: ZZ,ZZS
    INTEGER, ALLOCATABLE, DIMENSION(:) :: IIO_R,JJO_R,KKO_R,IOR_R,IIO_S,JJO_S,KKO_S,IOR_S
    INTEGER :: I_MIN_R=-10,I_MAX_R=-10,J_MIN_R=-10,J_MAX_R=-10,K_MIN_R=-10,K_MAX_R=-10,NIC_R=0, &