File: N:\mfix\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
16           USE geometry, only: odx_e, ody_n, axz_w, ayz_w
17     
18           USE is, only: is_defined, is_plane
19           USE is, only: is_i_w, is_i_e, is_j_s, is_j_n, is_k_t, is_k_b
20     
21           USE fun_avg, only: avg_x_h, avg_y_h, avg_z_h
22     
23           USE functions, only: funijk, ip_of, jp_of
24           USE functions, only: east_of, north_of, top_of
25           USE functions, only: is_on_mype_plus2layers
26     
27           USE compar, only: dead_cell_at
28           IMPLICIT NONE
29     ! Dummy arguments
30     !---------------------------------------------------------------------//
31     ! Gamma -- diffusion coefficient
32           DOUBLE PRECISION, INTENT(IN) :: Dif(DIMENSION_3)
33     ! Septadiagonal matrix A_m
34           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
35     ! Solids phase
36           INTEGER, INTENT(IN) :: M
37     
38     ! Local variables
39     !---------------------------------------------------------------------//
40     ! Diffusion parameter
41           DOUBLE PRECISION :: D_f
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, IJPK, IPJK, IJKTN, IJKTE
47     !---------------------------------------------------------------------//
48     
49     
50           DO L = 1, DIMENSION_IS
51              IF (IS_DEFINED(L)) THEN
52                 I1 = IS_I_W(L)
53                 I2 = IS_I_E(L)
54                 J1 = IS_J_S(L)
55                 J2 = IS_J_N(L)
56                 K1 = IS_K_B(L)
57                 K2 = IS_K_T(L)
58     
59                 IF (IS_PLANE(L) == 'N') THEN
60                    DO K = K1, K2
61                    DO J = J1, J2
62                    DO I = I1, I2
63                       IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
64                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
65                       IJK = FUNIJK(I,J,K)
66                       IJKT = TOP_OF(IJK)
67                       IJKN = NORTH_OF(IJK)
68                       IJKTN = TOP_OF(IJKN)
69                       IJPK = JP_OF(IJK)
70     
71                       D_F = AVG_Z_H(AVG_Y_H(DIF(IJK),DIF(IJKN),J),&
72                                     AVG_Y_H(DIF(IJKT),DIF(IJKTN),J),K)*&
73                             ODY_N(J)*AXZ_W(IJK)
74     
75                       A_M(IJK,north,M) = A_M(IJK,north,M) - D_F
76                       A_M(IJPK,south,M) = A_M(IJPK,south,M) - D_F
77                    ENDDO
78                    ENDDO
79                    ENDDO
80     
81                 ELSEIF (IS_PLANE(L) == 'E') THEN
82                    DO K = K1, K2
83                    DO J = J1, J2
84                    DO I = I1, I2
85                       IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
86                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
87                       IJK = FUNIJK(I,J,K)
88                       IJKE = EAST_OF(IJK)
89                       IJKT = TOP_OF(IJK)
90                       IJKTE = EAST_OF(IJKT)
91                       IPJK = IP_OF(IJK)
92     
93                       D_F = AVG_Z_H(AVG_X_H(DIF(IJK),DIF(IJKE),I),&
94                                     AVG_X_H(DIF(IJKT),DIF(IJKTE),I),K)*&
95                             ODX_E(I)*AYZ_W(IJK)
96     
97                       A_M(IJK,east,M) = A_M(IJK,east,M) - D_F
98                       A_M(IPJK,west,M) = A_M(IPJK,west,M) - D_F
99                    ENDDO
100                    ENDDO
101                    ENDDO
102                 ENDIF
103     
104              ENDIF   ! end if is_defined
105           ENDDO   ! end do dimension_is
106           RETURN
107           END SUBROUTINE DIF_W_IS
108