File: N:\mfix\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     
81                 M = 1
82                 IF (MMAX > 0) THEN
83                   ROPSX(:MMAX) = IC_ROP_S(L,:MMAX)
84                   TSX(:MMAX) = IC_T_S(L,:MMAX)
85                   USX(:MMAX) = IC_U_S(L,:MMAX)
86                   VSX(:MMAX) = IC_V_S(L,:MMAX)
87                   WSX(:MMAX) = IC_W_S(L,:MMAX)
88                   M = MMAX + 1
89                 ENDIF
90     
91                 DO K = IC_K_B(L), IC_K_T(L)
92                 DO J = IC_J_S(L), IC_J_N(L)
93                 DO I = IC_I_W(L), IC_I_E(L)
94                    IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
95                    IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
96                    IJK = FUNIJK(I,J,K)
97     
98                    IF (.NOT.WALL_AT(IJK)) THEN
99                       IF (EPGX /= UNDEFINED) EP_G(IJK) = EPGX
100     
101                       IF (IC_TYPE(L) == 'PATCH') THEN
102                           IF (PGX /= UNDEFINED) P_G(IJK) = SCALE_PRESSURE(PGX)
103                       ELSE
104                          P_G(IJK) = merge(SCALE_PRESSURE(PGX), UNDEFINED,           &
105                             PGX /= UNDEFINED)
106                       ENDIF
107     
108                       IF (PSX /= UNDEFINED) P_STAR(IJK) = PSX
109                       IF (TGX /= UNDEFINED) T_G(IJK) = TGX
110                       IF (IC_L_SCALE(L) /= UNDEFINED) L_SCALE(IJK) =       &
111                           IC_L_SCALE(L)
112     
113                       IF (NMAX(0) > 0) THEN
114                          WHERE (IC_X_G(L,:NMAX(0)) /= UNDEFINED)           &
115                             X_G(IJK,:NMAX(0)) = IC_X_G(L,:NMAX(0))
116                       ENDIF
117     
118                       IF (NScalar > 0) THEN
119                          WHERE (IC_Scalar(L,:NScalar) /= UNDEFINED)        &
120                             Scalar(IJK,:NScalar) = IC_Scalar(L,:NScalar)
121                       ENDIF
122     
123                       IF (K_Epsilon) THEN
124                           IF (IC_K_Turb_G(L) /= UNDEFINED)                 &
125                               K_Turb_G(IJK) = IC_K_Turb_G(L)
126                           IF (IC_E_Turb_G(L) /= UNDEFINED)                 &
127                               E_Turb_G(IJK) = IC_E_Turb_G(L)
128                       ENDIF
129     
130                       IF (UGX /= UNDEFINED) U_G(IJK) = UGX
131                       IF (VGX /= UNDEFINED) V_G(IJK) = VGX
132                       IF (WGX /= UNDEFINED) W_G(IJK) = WGX
133     
134                       GAMA_RG(IJK) = IC_GAMA_RG(L)
135                       T_RG(IJK) = merge(IC_T_RG(L), ZERO,                  &
136                          IC_T_RG(L) /= UNDEFINED)
137     
138                       DO M = 1, MMAX
139                         IF (ROPSX(M) /= UNDEFINED) ROP_S(IJK,M) = ROPSX(M)
140                         IF (TSX(M) /= UNDEFINED) T_S(IJK,M) = TSX(M)
141                         IF (IC_THETA_M(L,M) /= UNDEFINED)                  &
142                            THETA_M(IJK,M) = IC_THETA_M(L,M)
143                         IF (USX(M) /= UNDEFINED) U_S(IJK,M) = USX(M)
144                         IF (VSX(M) /= UNDEFINED) V_S(IJK,M) = VSX(M)
145                         IF (WSX(M) /= UNDEFINED) W_S(IJK,M) = WSX(M)
146     
147                         GAMA_RS(IJK,M) = IC_GAMA_RS(L,M)
148                         T_RS(IJK,M) = merge(IC_T_RS(L,M),ZERO,             &
149                            IC_T_RS(L,M) /= UNDEFINED)
150     
151                         IF (NMAX(M) > 0) THEN
152                             WHERE (IC_X_S(L,M,:NMAX(M)) /= UNDEFINED)      &
153                                X_S(IJK,M,:NMAX(M)) = IC_X_S(L,M,:NMAX(M))
154                         ENDIF
155                       ENDDO
156     
157     ! for GHD theory to compute mixture IC of velocity and density
158                       IF(KT_TYPE_ENUM == GHD_2007) THEN
159                          ROP_S(IJK,MMAX) = ZERO
160                          U_S(IJK,MMAX) = ZERO
161                          V_S(IJK,MMAX) = ZERO
162                          W_S(IJK,MMAX) = ZERO
163                          THETA_M(IJK,MMAX) = ZERO
164                          nTOT = ZERO
165                          nM = ZERO
166                          DO M = 1, SMAX
167                             IF (ROPSX(M) /= UNDEFINED) THEN
168                                ROP_S(IJK,MMAX) = ROP_S(IJK,MMAX) + ROPSX(M)
169                                nM = ROPSX(M)*6d0 /                         &
170                                   (PI*D_p(IJK,M)**3*RO_S(IJK,M))
171                                nTOT = nTOT + nM
172                             ENDIF
173                             IF (IC_THETA_M(L,M) /= UNDEFINED)              &
174                                THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) +     &
175     
176                                nM*IC_THETA_M(L,M)
177                             IF(USX(M) /= UNDEFINED .AND.                   &
178                                ROPSX(M) /= UNDEFINED)  U_S(IJK,MMAX) =     &
179                                U_S(IJK,MMAX) + ROPSX(M)*USX(M)
180     
181                             IF(VSX(M) /= UNDEFINED .AND.                   &
182                                ROPSX(M) /= UNDEFINED) V_S(IJK,MMAX) =      &
183                                V_S(IJK,MMAX) +  ROPSX(M)*VSX(M)
184     
185                             IF(WSX(M) /= UNDEFINED .AND.                   &
186                                ROPSX(M) /= UNDEFINED) W_S(IJK,MMAX) =      &
187                                W_S(IJK,MMAX) +  ROPSX(M)*WSX(M)
188                          ENDDO
189     
190     ! If ropsTotal > 0 then RoN_T > 0
191                          IF(ROP_S(IJK,MMAX) > ZERO) THEN
192                             U_S(IJK,MMAX) = U_S(IJK,MMAX) / ROP_S(IJK,MMAX)
193                             V_S(IJK,MMAX) = V_S(IJK,MMAX) / ROP_S(IJK,MMAX)
194                             W_S(IJK,MMAX) = W_S(IJK,MMAX) / ROP_S(IJK,MMAX)
195                             THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) / nTOT
196     ! For initially empty bed:
197                          ELSE
198                             U_S(IJK,MMAX) = U_S(IJK,MMAX)
199                             V_S(IJK,MMAX) = V_S(IJK,MMAX)
200                             W_S(IJK,MMAX) = W_S(IJK,MMAX)
201     ! Set T > 0 in case Ti > 0
202                             DO M = 1, SMAX
203                                THETA_M(IJK,MMAX) = THETA_M(IJK,M)
204                             ENDDO
205                             IF(THETA_M(IJK,MMAX)==ZERO)                     &
206                                THETA_M(IJK,MMAX) = small_number
207                         ENDIF
208                       ENDIF
209     ! end of modifications for GHD theory
210                    ENDIF     ! Fluid at
211                 ENDDO   ! over i
212                 ENDDO   ! over j
213                 ENDDO   ! over k
214              ENDIF   ! if (ic_defined)
215           ENDDO   ! over dimension_ic
216     
217           CALL SEND_RECV(L_SCALE,2)
218     
219           RETURN
220           END SUBROUTINE SET_IC
221