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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: DIF_W_IS                                                C
4     !  Purpose: Remove diffusive fluxes across internal surfaces.          C
5     !  (Make user defined internal surfaces non-conducting)                C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 30-APR-97  C
8     !                                                                      C
9     !                                                                      C
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
11           SUBROUTINE DIF_W_IS(DIF, A_M, M)
12     
13     ! Modules
14     !---------------------------------------------------------------------//
15           USE param, only: dimension_is, dimension_3, dimension_m
16           USE matrix, only: e, w, n, s
17           USE geometry, only: odx_e, ody_n, axz_w, ayz_w
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, jp_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     ! Dummy arguments
31     !---------------------------------------------------------------------//
32     ! Gamma -- diffusion coefficient
33           DOUBLE PRECISION, INTENT(IN) :: Dif(DIMENSION_3)
34     ! Septadiagonal matrix A_m
35           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
36     ! Solids phase
37           INTEGER, INTENT(IN) :: M
38     
39     ! Local variables
40     !---------------------------------------------------------------------//
41     ! Diffusion parameter
42           DOUBLE PRECISION :: D_f
43     ! Internal surface
44           INTEGER :: L
45     ! Indices
46           INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK
47           INTEGER :: IJKE, IJKN, IJKT, IJPK, IPJK, IJKTN, IJKTE
48     !---------------------------------------------------------------------//
49     
50     
51           DO L = 1, DIMENSION_IS
52              IF (IS_DEFINED(L)) THEN
53                 I1 = IS_I_W(L)
54                 I2 = IS_I_E(L)
55                 J1 = IS_J_S(L)
56                 J2 = IS_J_N(L)
57                 K1 = IS_K_B(L)
58                 K2 = IS_K_T(L)
59     
60                 IF (IS_PLANE(L) == 'N') THEN
61                    DO K = K1, K2
62                    DO J = J1, J2
63                    DO I = I1, I2
64                       IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
65                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
66                       IJK = FUNIJK(I,J,K)
67                       IJKT = TOP_OF(IJK)
68                       IJKN = NORTH_OF(IJK)
69                       IJKTN = TOP_OF(IJKN)
70                       IJPK = JP_OF(IJK)
71     
72                       D_F = AVG_Z_H(AVG_Y_H(DIF(IJK),DIF(IJKN),J),&
73                                     AVG_Y_H(DIF(IJKT),DIF(IJKTN),J),K)*&
74                             ODY_N(J)*AXZ_W(IJK)
75     
76                       A_M(IJK,N,M) = A_M(IJK,N,M) - D_F
77                       A_M(IJPK,S,M) = A_M(IJPK,S,M) - D_F
78                    ENDDO
79                    ENDDO
80                    ENDDO
81     
82                 ELSEIF (IS_PLANE(L) == 'E') THEN
83                    DO K = K1, K2
84                    DO J = J1, J2
85                    DO I = I1, I2
86                       IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
87                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
88                       IJK = FUNIJK(I,J,K)
89                       IJKE = EAST_OF(IJK)
90                       IJKT = TOP_OF(IJK)
91                       IJKTE = EAST_OF(IJKT)
92                       IPJK = IP_OF(IJK)
93     
94                       D_F = AVG_Z_H(AVG_X_H(DIF(IJK),DIF(IJKE),I),&
95                                     AVG_X_H(DIF(IJKT),DIF(IJKTE),I),K)*&
96                             ODX_E(I)*AYZ_W(IJK)
97     
98                       A_M(IJK,E,M) = A_M(IJK,E,M) - D_F
99                       A_M(IPJK,W,M) = A_M(IPJK,W,M) - D_F
100                    ENDDO
101                    ENDDO
102                    ENDDO
103                 ENDIF
104     
105              ENDIF   ! end if is_defined
106           ENDDO   ! end do dimension_is
107           RETURN
108           END SUBROUTINE DIF_W_IS
109