File: N:\mfix\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
15           USE geometry, only: do_k
16           USE geometry, only: odx_e, ox, odz_t, axy_v, ayz_v
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, kp_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     
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     ! Internal surface
42           INTEGER :: L
43     ! Indices
44           INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK
45           INTEGER :: IJKE, IJKN, IJKT, IPJK, IJKP, IJKNE, IJKTN
46     ! Diffusion parameter
47           DOUBLE PRECISION :: D_f
48     !---------------------------------------------------------------------//
49     
50     
51     ! Make user defined internal surfaces non-conducting
52           DO L = 1, DIMENSION_IS
53              IF (IS_DEFINED(L)) THEN
54                 I1 = IS_I_W(L)
55                 I2 = IS_I_E(L)
56                 J1 = IS_J_S(L)
57                 J2 = IS_J_N(L)
58                 K1 = IS_K_B(L)
59                 K2 = IS_K_T(L)
60     
61                 IF (IS_PLANE(L) == 'E') THEN
62                    DO K = K1, K2
63                    DO J = J1, J2
64                    DO I = I1, I2
65                       IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
66                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
67                       IJK = FUNIJK(I,J,K)
68                       IJKE = EAST_OF(IJK)
69                       IJKN = NORTH_OF(IJK)
70                       IJKNE = EAST_OF(IJKN)
71                       IPJK = IP_OF(IJK)
72     
73                       D_F = AVG_Y_H(AVG_X_H(DIF(IJK),DIF(IJKE),I),&
74                                     AVG_X_H(DIF(IJKN),DIF(IJKNE),I),J)*&
75                             ODX_E(I)*AYZ_V(IJK)
76     
77                       A_M(IJK,east,M) = A_M(IJK,east,M) - D_F
78                       A_M(IPJK,west,M) = A_M(IPJK,west,M) - D_F
79                    ENDDO
80                    ENDDO
81                    ENDDO
82     
83                 ELSEIF (IS_PLANE(L) == 'T') THEN
84                    IF (DO_K) THEN
85                       DO K = K1, K2
86                       DO J = J1, J2
87                       DO I = I1, I2
88                       IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
89                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
90                          IJK = FUNIJK(I,J,K)
91                          IJKN = NORTH_OF(IJK)
92                          IJKT = TOP_OF(IJK)
93                          IJKTN = NORTH_OF(IJKT)
94                          IJKP = KP_OF(IJK)
95     
96                          D_F = AVG_Y_H(AVG_Z_H(DIF(IJK),DIF(IJKT),K),&
97                                        AVG_Z_H(DIF(IJKN),DIF(IJKTN),K),J)*&
98                                OX(I)*ODZ_T(K)*AXY_V(IJK)
99     
100                          A_M(IJK,top,M) = A_M(IJK,top,M) - D_F
101                          A_M(IJKP,bottom,M) = A_M(IJKP,bottom,M) - D_F
102                       ENDDO
103                       ENDDO
104                       ENDDO
105                    ENDIF
106                 ENDIF
107     
108              ENDIF   ! end if is_defined
109           ENDDO   ! end do dimension_is
110           RETURN
111           END SUBROUTINE DIF_V_IS
112