File: RELATIVE:/../../../mfix.git/model/dif_v_is.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: DIF_V_IS                                                C
4     !  Purpose: Remove diffusive fluxes across internal surfaces.          C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 30-APR-97  C
7     !                                                                      C
8     !                                                                      C
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
10           SUBROUTINE DIF_V_IS(DIF, A_M, M)
11     
12     ! Modules
13     !---------------------------------------------------------------------//
14           USE param, only: dimension_is, dimension_3, dimension_m
15           USE matrix, only: e, w, t, b
16           USE geometry, only: do_k
17           USE geometry, only: odx_e, ox, odz_t, axy_v, ayz_v
18     
19           USE is, only: is_defined, is_plane
20           USE is, only: is_i_w, is_i_e, is_j_s, is_j_n, is_k_t, is_k_b
21     
22           USE fun_avg, only: avg_x_h, avg_y_h, avg_z_h
23     
24           USE functions, only: funijk, ip_of, kp_of
25           USE functions, only: east_of, north_of, top_of
26           USE functions, only: is_on_mype_plus2layers
27     
28           USE compar, only: dead_cell_at
29           IMPLICIT NONE
30     
31     ! Dummy arguments
32     !---------------------------------------------------------------------//
33     ! Gamma -- diffusion coefficient
34           DOUBLE PRECISION, INTENT(IN) :: Dif(DIMENSION_3)
35     ! Septadiagonal matrix A_m
36           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
37     ! Solids phase
38           INTEGER, INTENT(IN) :: M
39     
40     ! Local variables
41     !---------------------------------------------------------------------//
42     ! Internal surface
43           INTEGER :: L
44     ! Indices
45           INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK
46           INTEGER :: IJKE, IJKN, IJKT, IPJK, IJKP, IJKNE, IJKTN
47     ! Diffusion parameter
48           DOUBLE PRECISION :: D_f
49     !---------------------------------------------------------------------//
50     
51     
52     ! Make user defined internal surfaces non-conducting
53           DO L = 1, DIMENSION_IS
54              IF (IS_DEFINED(L)) THEN
55                 I1 = IS_I_W(L)
56                 I2 = IS_I_E(L)
57                 J1 = IS_J_S(L)
58                 J2 = IS_J_N(L)
59                 K1 = IS_K_B(L)
60                 K2 = IS_K_T(L)
61     
62                 IF (IS_PLANE(L) == 'E') THEN
63                    DO K = K1, K2
64                    DO J = J1, J2
65                    DO I = I1, I2
66                       IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
67                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
68                       IJK = FUNIJK(I,J,K)
69                       IJKE = EAST_OF(IJK)
70                       IJKN = NORTH_OF(IJK)
71                       IJKNE = EAST_OF(IJKN)
72                       IPJK = IP_OF(IJK)
73     
74                       D_F = AVG_Y_H(AVG_X_H(DIF(IJK),DIF(IJKE),I),&
75                                     AVG_X_H(DIF(IJKN),DIF(IJKNE),I),J)*&
76                             ODX_E(I)*AYZ_V(IJK)
77     
78                       A_M(IJK,E,M) = A_M(IJK,E,M) - D_F
79                       A_M(IPJK,W,M) = A_M(IPJK,W,M) - D_F
80                    ENDDO
81                    ENDDO
82                    ENDDO
83     
84                 ELSEIF (IS_PLANE(L) == 'T') THEN
85                    IF (DO_K) THEN
86                       DO K = K1, K2
87                       DO J = J1, J2
88                       DO I = I1, I2
89                       IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
90                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
91                          IJK = FUNIJK(I,J,K)
92                          IJKN = NORTH_OF(IJK)
93                          IJKT = TOP_OF(IJK)
94                          IJKTN = NORTH_OF(IJKT)
95                          IJKP = KP_OF(IJK)
96     
97                          D_F = AVG_Y_H(AVG_Z_H(DIF(IJK),DIF(IJKT),K),&
98                                        AVG_Z_H(DIF(IJKN),DIF(IJKTN),K),J)*&
99                                OX(I)*ODZ_T(K)*AXY_V(IJK)
100     
101                          A_M(IJK,T,M) = A_M(IJK,T,M) - D_F
102                          A_M(IJKP,B,M) = A_M(IJKP,B,M) - D_F
103                       ENDDO
104                       ENDDO
105                       ENDDO
106                    ENDIF
107                 ENDIF
108     
109              ENDIF   ! end if is_defined
110           ENDDO   ! end do dimension_is
111           RETURN
112           END SUBROUTINE DIF_V_IS
113