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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: DIF_U_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_U_IS(DIF, A_M, M)
12     
13     ! Modules
14     !---------------------------------------------------------------------//
15           USE param, only: dimension_is, dimension_3, dimension_m
16           USE matrix, only: n, s, t, b
17           USE geometry, only: do_k
18           USE geometry, only: ody_n, ox_e, odz_t, axz_u, axy_u
19     
20           USE is, only: is_defined, is_plane
21           USE is, only: is_i_w, is_i_e, is_j_s, is_j_n, is_k_t, is_k_b
22     
23           USE fun_avg, only: avg_x_h, avg_y_h, avg_z_h
24     
25           USE functions, only: funijk, jp_of, kp_of
26           USE functions, only: east_of, north_of, top_of
27           USE functions, only: is_on_mype_plus2layers
28     
29           USE compar, only: dead_cell_at
30           USE compar, only: istart2, jstart2, kstart2
31           USE compar, only: iend2, jend2, kend2
32           IMPLICIT NONE
33     
34     ! Dummy arguments
35     !---------------------------------------------------------------------//
36     ! Gamma -- diffusion coefficient
37           DOUBLE PRECISION, INTENT(IN) :: Dif(DIMENSION_3)
38     ! Septadiagonal matrix A_m
39           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
40     ! Solids phase
41           INTEGER, INTENT(IN) :: M
42     
43     ! Local variables
44     !---------------------------------------------------------------------//
45     ! Internal surface
46           INTEGER :: L
47     ! Indices
48           INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK
49           INTEGER :: IJKE, IJKN, IJKT, IJPK, IJKP, IJKNE, IJKTE
50     ! Diffusion parameter
51           DOUBLE PRECISION :: D_f
52     !---------------------------------------------------------------------//
53     
54           DO L = 1, DIMENSION_IS
55              IF (IS_DEFINED(L)) THEN
56                 I1 = IS_I_W(L)
57                 I2 = IS_I_E(L)
58                 J1 = IS_J_S(L)
59                 J2 = IS_J_N(L)
60                 K1 = IS_K_B(L)
61                 K2 = IS_K_T(L)
62     
63     ! Limit I1, I2 and all to local processor first ghost layer
64                 IF(I1.LE.IEND2)   I1 = MAX(I1, ISTART2)
65                 IF(J1.LE.JEND2)   J1 = MAX(J1, JSTART2)
66                 IF(K1.LE.KEND2)   K1 = MAX(K1, KSTART2)
67                 IF(I2.GE.ISTART2) I2 = MIN(I2, IEND2)
68                 IF(J2.GE.JSTART2) J2 = MIN(J2, JEND2)
69                 IF(K2.GE.KSTART2) K2 = MIN(K2, KEND2)
70     
71                 IF (IS_PLANE(L) == 'N') THEN
72                    DO K = K1, K2
73                    DO J = J1, J2
74                    DO I = I1, I2
75                       IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
76                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
77                       IJK = FUNIJK(I,J,K)
78                       IJKE = EAST_OF(IJK)
79                       IJKN = NORTH_OF(IJK)
80                       IJKNE = EAST_OF(IJKN)
81                       IJPK = JP_OF(IJK)
82     
83                       D_F = AVG_X_H(AVG_Y_H(DIF(IJK),DIF(IJKN),J),&
84                                     AVG_Y_H(DIF(IJKE),DIF(IJKNE),J),I)*&
85                             ODY_N(J)*AXZ_U(IJK)
86     
87                       A_M(IJK,N,M) = A_M(IJK,N,M) - D_F
88                       A_M(IJPK,S,M) = A_M(IJPK,S,M) - D_F
89                    ENDDO
90                    ENDDO
91                    ENDDO
92     
93                 ELSEIF (IS_PLANE(L) == 'T') THEN
94                    IF (DO_K) THEN
95                       DO K = K1, K2
96                       DO J = J1, J2
97                       DO I = I1, I2
98                          IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
99                          IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
100                          IJK = FUNIJK(I,J,K)
101                          IJKE = EAST_OF(IJK)
102                          IJKT = TOP_OF(IJK)
103                          IJKTE = EAST_OF(IJKT)
104                          IJKP = KP_OF(IJK)
105     
106                          D_F = AVG_X_H(AVG_Z_H(DIF(IJK),DIF(IJKT),K),&
107                                        AVG_Z_H(DIF(IJKE),DIF(IJKTE),K),I)*&
108                                OX_E(I)*ODZ_T(K)*AXY_U(IJK)
109     
110                          A_M(IJK,T,M) = A_M(IJK,T,M) - D_F
111                          A_M(IJKP,B,M) = A_M(IJKP,B,M) - D_F
112                       ENDDO
113                       ENDDO
114                       ENDDO
115                    ENDIF
116                 ENDIF
117     
118              ENDIF   ! end if is_defined
119           ENDDO   ! end do dimension_is
120     
121     
122           RETURN
123           END SUBROUTINE DIF_U_IS
124     
125