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

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