File: N:\mfix\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
16           USE geometry, only: do_k
17           USE geometry, only: ody_n, ox_e, odz_t, axz_u, axy_u
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, jp_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           USE compar, only: istart2, jstart2, kstart2
30           USE compar, only: iend2, jend2, kend2
31           IMPLICIT NONE
32     
33     ! Dummy arguments
34     !---------------------------------------------------------------------//
35     ! Gamma -- diffusion coefficient
36           DOUBLE PRECISION, INTENT(IN) :: Dif(DIMENSION_3)
37     ! Septadiagonal matrix A_m
38           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
39     ! Solids phase
40           INTEGER, INTENT(IN) :: M
41     
42     ! Local variables
43     !---------------------------------------------------------------------//
44     ! Internal surface
45           INTEGER :: L
46     ! Indices
47           INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK
48           INTEGER :: IJKE, IJKN, IJKT, IJPK, IJKP, IJKNE, IJKTE
49     ! Diffusion parameter
50           DOUBLE PRECISION :: D_f
51     !---------------------------------------------------------------------//
52     
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     ! Limit I1, I2 and all to local processor first ghost layer
63                 IF(I1.LE.IEND2)   I1 = MAX(I1, ISTART2)
64                 IF(J1.LE.JEND2)   J1 = MAX(J1, JSTART2)
65                 IF(K1.LE.KEND2)   K1 = MAX(K1, KSTART2)
66                 IF(I2.GE.ISTART2) I2 = MIN(I2, IEND2)
67                 IF(J2.GE.JSTART2) J2 = MIN(J2, JEND2)
68                 IF(K2.GE.KSTART2) K2 = MIN(K2, KEND2)
69     
70                 IF (IS_PLANE(L) == 'N') THEN
71                    DO K = K1, K2
72                    DO J = J1, J2
73                    DO I = I1, I2
74                       IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
75                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
76                       IJK = FUNIJK(I,J,K)
77                       IJKE = EAST_OF(IJK)
78                       IJKN = NORTH_OF(IJK)
79                       IJKNE = EAST_OF(IJKN)
80                       IJPK = JP_OF(IJK)
81     
82                       D_F = AVG_X_H(AVG_Y_H(DIF(IJK),DIF(IJKN),J),&
83                                     AVG_Y_H(DIF(IJKE),DIF(IJKNE),J),I)*&
84                             ODY_N(J)*AXZ_U(IJK)
85     
86                       A_M(IJK,north,M) = A_M(IJK,north,M) - D_F
87                       A_M(IJPK,south,M) = A_M(IJPK,south,M) - D_F
88                    ENDDO
89                    ENDDO
90                    ENDDO
91     
92                 ELSEIF (IS_PLANE(L) == 'T') THEN
93                    IF (DO_K) THEN
94                       DO K = K1, K2
95                       DO J = J1, J2
96                       DO I = I1, I2
97                          IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
98                          IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
99                          IJK = FUNIJK(I,J,K)
100                          IJKE = EAST_OF(IJK)
101                          IJKT = TOP_OF(IJK)
102                          IJKTE = EAST_OF(IJKT)
103                          IJKP = KP_OF(IJK)
104     
105                          D_F = AVG_X_H(AVG_Z_H(DIF(IJK),DIF(IJKT),K),&
106                                        AVG_Z_H(DIF(IJKE),DIF(IJKTE),K),I)*&
107                                OX_E(I)*ODZ_T(K)*AXY_U(IJK)
108     
109                          A_M(IJK,top,M) = A_M(IJK,top,M) - D_F
110                          A_M(IJKP,bottom,M) = A_M(IJKP,bottom,M) - D_F
111                       ENDDO
112                       ENDDO
113                       ENDDO
114                    ENDIF
115                 ENDIF
116     
117              ENDIF   ! end if is_defined
118           ENDDO   ! end do dimension_is
119     
120     
121           RETURN
122           END SUBROUTINE DIF_U_IS
123     
124