File: /nfs/home/0/users/jenkins/mfix.git/model/set_ic.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SET_IC                                                  C
4     !  Purpose: This module sets all the initial conditions.               C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 21-JAN-92  C
7     !  Reviewer:M. Syamlal, S. Venkatesan, P. Nicoletti,  Date: 29-JAN-92  C
8     !           W. Rogers                                                  C
9     !                                                                      C
10     !  Revision Number: 1                                                  C
11     !  Purpose: Modifications for initializing solids pressure             C
12     !  Author: M. Syamlal                                 Date: 11-MAY-93  C
13     !                                                                      C
14     !                                                                      C
15     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
16     
17           SUBROUTINE SET_IC
18     
19     !-----------------------------------------------
20     ! Modules
21     !-----------------------------------------------
22           USE param
23           USE param1
24           USE geometry
25           USE constant
26           USE physprop
27           USE ic
28           USE fldvar
29           USE visc_g
30           USE indices
31           USE scales
32           USE energy
33           USE scalars
34           USE compar
35           USE run
36           USE sendrecv
37           USE solids_pressure
38           USE functions
39           IMPLICIT NONE
40     !-----------------------------------------------
41     ! Local variables
42     !-----------------------------------------------
43     ! indices
44           INTEGER :: I, J, K, IJK, M
45     ! Local index for initial condition
46           INTEGER :: L
47     ! Temporary variable for storing IC_EP_g
48           DOUBLE PRECISION :: EPGX
49     ! Temporary variable for storing IC_P_g
50           DOUBLE PRECISION :: PGX
51     ! Temporary variable for storing P_star
52           DOUBLE PRECISION :: PSX
53     ! Temporary variable for storing IC_T_g
54           DOUBLE PRECISION :: TGX
55     ! Temporary variable for storing IC_U_g
56           DOUBLE PRECISION :: UGX
57     ! Temporary variable for storing IC_V_g
58           DOUBLE PRECISION :: VGX
59     ! Temporary variable for storing IC_W_g
60           DOUBLE PRECISION :: WGX
61     ! Temporary variable for storing IC_ROP_s
62           DOUBLE PRECISION :: ROPSX (DIMENSION_M)
63     ! Temporary variable for storing IC_T_s
64           DOUBLE PRECISION :: TSX (DIMENSION_M)
65     ! Temporary variable for storing IC_U_s
66           DOUBLE PRECISION :: USX (DIMENSION_M)
67     ! Temporary variable for storing IC_V_s
68           DOUBLE PRECISION :: VSX (DIMENSION_M)
69     ! Temporary variable for storing IC_W_s
70           DOUBLE PRECISION :: WSX (DIMENSION_M)
71     ! number density for GHD theory
72           DOUBLE PRECISION :: nM, nTOT
73     !-----------------------------------------------
74     
75     !  Set the initial conditions.
76           DO L = 1, DIMENSION_IC
77              IF (IC_DEFINED(L)) THEN
78                 EPGX = IC_EP_G(L)
79                 PGX = IC_P_G(L)
80                 PSX = IC_P_STAR(L)
81                 IF (PSX==UNDEFINED .AND. IC_TYPE(L)/='PATCH') PSX = ZERO
82                 TGX = IC_T_G(L)
83                 UGX = IC_U_G(L)
84                 VGX = IC_V_G(L)
85                 WGX = IC_W_G(L)
86                 IF (IC_L_SCALE(L) /= UNDEFINED) RECALC_VISC_G = .TRUE.
87                 IF(K_Epsilon) RECALC_VISC_G = .TRUE.
88     
89                 M = 1
90                 IF (MMAX > 0) THEN
91                   ROPSX(:MMAX) = IC_ROP_S(L,:MMAX)
92                   TSX(:MMAX) = IC_T_S(L,:MMAX)
93                   USX(:MMAX) = IC_U_S(L,:MMAX)
94                   VSX(:MMAX) = IC_V_S(L,:MMAX)
95                   WSX(:MMAX) = IC_W_S(L,:MMAX)
96                   M = MMAX + 1
97                 ENDIF
98     
99                 DO K = IC_K_B(L), IC_K_T(L)
100                 DO J = IC_J_S(L), IC_J_N(L)
101                 DO I = IC_I_W(L), IC_I_E(L)
102                    IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
103                    IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
104                    IJK = FUNIJK(I,J,K)
105     
106                    IF (.NOT.WALL_AT(IJK)) THEN
107                       IF (EPGX /= UNDEFINED) EP_G(IJK) = EPGX
108                       IF (IC_TYPE(L) == 'PATCH') THEN
109                           IF (PGX /= UNDEFINED) P_G(IJK) = SCALE(PGX)
110                       ELSE
111                          P_G(IJK) = merge(SCALE(PGX), UNDEFINED,           &
112                             PGX /= UNDEFINED)
113                       ENDIF
114                       IF (PSX /= UNDEFINED) P_STAR(IJK) = PSX
115                       IF (TGX /= UNDEFINED) T_G(IJK) = TGX
116                       IF (IC_L_SCALE(L) /= UNDEFINED) L_SCALE(IJK) =       &
117                           IC_L_SCALE(L)
118     
119                       IF (NMAX(0) > 0) THEN
120                          WHERE (IC_X_G(L,:NMAX(0)) /= UNDEFINED)           &
121                             X_G(IJK,:NMAX(0)) = IC_X_G(L,:NMAX(0))
122                       ENDIF
123     
124                       IF (NScalar > 0) THEN
125                          WHERE (IC_Scalar(L,:NScalar) /= UNDEFINED)        &
126                             Scalar(IJK,:NScalar) = IC_Scalar(L,:NScalar)
127                       ENDIF
128     
129                       IF (K_Epsilon) THEN
130                           IF (IC_K_Turb_G(L) /= UNDEFINED)                 &
131                               K_Turb_G(IJK) = IC_K_Turb_G(L)
132                           IF (IC_E_Turb_G(L) /= UNDEFINED)                 &
133                               E_Turb_G(IJK) = IC_E_Turb_G(L)
134                       ENDIF
135     
136                       IF (UGX /= UNDEFINED) U_G(IJK) = UGX
137                       IF (VGX /= UNDEFINED) V_G(IJK) = VGX
138                       IF (WGX /= UNDEFINED) W_G(IJK) = WGX
139     
140                       GAMA_RG(IJK) = IC_GAMA_RG(L)
141                       T_RG(IJK) = merge(IC_T_RG(L), ZERO,                  &
142                          IC_T_RG(L) /= UNDEFINED)
143     
144                       DO M = 1, MMAX
145                         IF (ROPSX(M) /= UNDEFINED) ROP_S(IJK,M) = ROPSX(M)
146                         IF (TSX(M) /= UNDEFINED) T_S(IJK,M) = TSX(M)
147                         IF (IC_THETA_M(L,M) /= UNDEFINED)                  &
148                            THETA_M(IJK,M) = IC_THETA_M(L,M)
149                         IF (USX(M) /= UNDEFINED) U_S(IJK,M) = USX(M)
150                         IF (VSX(M) /= UNDEFINED) V_S(IJK,M) = VSX(M)
151                         IF (WSX(M) /= UNDEFINED) W_S(IJK,M) = WSX(M)
152     
153                         GAMA_RS(IJK,M) = IC_GAMA_RS(L,M)
154                         T_RS(IJK,M) = merge(IC_T_RS(L,M),ZERO,             &
155                            IC_T_RS(L,M) /= UNDEFINED)
156     
157                         IF (NMAX(M) > 0) THEN
158                             WHERE (IC_X_S(L,M,:NMAX(M)) /= UNDEFINED)      &
159                                X_S(IJK,M,:NMAX(M)) = IC_X_S(L,M,:NMAX(M))
160                         ENDIF
161                       ENDDO
162     
163     ! for GHD theory to compute mixture IC of velocity and density
164                       IF(KT_TYPE_ENUM == GHD_2007) THEN
165                          ROP_S(IJK,MMAX) = ZERO
166                          U_S(IJK,MMAX) = ZERO
167                          V_S(IJK,MMAX) = ZERO
168                          W_S(IJK,MMAX) = ZERO
169                          THETA_M(IJK,MMAX) = ZERO
170                          nTOT = ZERO
171                          nM = ZERO
172                          DO M = 1, SMAX
173                             IF (ROPSX(M) /= UNDEFINED) THEN
174                                ROP_S(IJK,MMAX) = ROP_S(IJK,MMAX) + ROPSX(M)
175                                nM = ROPSX(M)*6d0 /                         &
176                                   (PI*D_p(IJK,M)**3*RO_S(IJK,M))
177                                nTOT = nTOT + nM
178                             ENDIF
179                             IF (IC_THETA_M(L,M) /= UNDEFINED)              &
180                                THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) +     &
181     
182                                nM*IC_THETA_M(L,M)
183                             IF(USX(M) /= UNDEFINED .AND.                   &
184                                ROPSX(M) /= UNDEFINED)  U_S(IJK,MMAX) =     &
185                                U_S(IJK,MMAX) + ROPSX(M)*USX(M)
186     
187                             IF(VSX(M) /= UNDEFINED .AND.                   &
188                                ROPSX(M) /= UNDEFINED) V_S(IJK,MMAX) =      &
189                                V_S(IJK,MMAX) +  ROPSX(M)*VSX(M)
190     
191                             IF(WSX(M) /= UNDEFINED .AND.                   &
192                                ROPSX(M) /= UNDEFINED) W_S(IJK,MMAX) =      &
193                                W_S(IJK,MMAX) +  ROPSX(M)*WSX(M)
194                          ENDDO
195     
196     ! If ropsTotal > 0 then RoN_T > 0
197                          IF(ROP_S(IJK,MMAX) > ZERO) THEN
198                             U_S(IJK,MMAX) = U_S(IJK,MMAX) / ROP_S(IJK,MMAX)
199                             V_S(IJK,MMAX) = V_S(IJK,MMAX) / ROP_S(IJK,MMAX)
200                             W_S(IJK,MMAX) = W_S(IJK,MMAX) / ROP_S(IJK,MMAX)
201                             THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) / nTOT
202     ! For initially empty bed:
203                          ELSE
204                             U_S(IJK,MMAX) = U_S(IJK,MMAX)
205                             V_S(IJK,MMAX) = V_S(IJK,MMAX)
206                             W_S(IJK,MMAX) = W_S(IJK,MMAX)
207     ! Set T > 0 in case Ti > 0
208                             DO M = 1, SMAX
209                                THETA_M(IJK,MMAX) = THETA_M(IJK,M)
210                             ENDDO
211                             IF(THETA_M(IJK,MMAX)==ZERO)                     &
212                                THETA_M(IJK,MMAX) = small_number
213                         ENDIF
214                       ENDIF
215     ! end of modifications for GHD theory
216                    ENDIF     ! Fluid at
217                 ENDDO   ! over i
218                 ENDDO   ! over j
219                 ENDDO   ! over k
220              ENDIF   ! if (ic_defined)
221           ENDDO   ! over dimension_ic
222     
223           CALL SEND_RECV(L_SCALE,2)
224     
225           RETURN
226           END SUBROUTINE SET_IC
227