File: N:\mfix\model\cartesian_grid\CG_source_w_g.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CG_SOURCE_W_g(A_m, B_m, IER)                           C
4     !  Purpose: Determine contribution of cut-cell to source terms         C
5     !  for W_g momentum eq.                                                C
6     !                                                                      C
7     !                                                                      C
8     !  Author: Jeff Dietiker                              Date: 01-MAY-09  C
9     !                                                                      C
10     !                                                                      C
11     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
12           SUBROUTINE CG_SOURCE_W_G(A_M, B_M)
13     
14     ! Modules
15     !-----------------------------------------------
16           USE bc
17           USE compar
18           USE fldvar
19           USE fun_avg
20           use geometry, only: vol, vol_w
21           use geometry, only: axy, axz, ayz
22           use geometry, only: axy_w, axz_w, ayz_w
23           USE indices
24           USE param
25           USE param1
26           USE physprop
27           USE run
28           USE toleranc
29           USE visc_g
30           USE cutcell
31           USE quadric
32           IMPLICIT NONE
33     
34     ! Dummy arguments
35     !-----------------------------------------------
36     ! Septadiagonal matrix A_m
37           DOUBLE PRECISION :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
38     ! Vector b_m
39           DOUBLE PRECISION :: B_m(DIMENSION_3, 0:DIMENSION_M)
40     
41     ! Local variables
42     !-----------------------------------------------
43     ! Indices
44           INTEGER :: I, J, K, IJK, IJKT, IMJK, IJKP, IMJKP
45           INTEGER :: IJKE, IJKW, IJKTE, IM, IPJK
46     ! Phase index
47           INTEGER :: M
48     ! Average volume fraction
49           DOUBLE PRECISION :: EPGA
50           INTEGER :: JM,IP,JP,IJMK,IJPK
51           INTEGER :: IJKC,IJKN,IJKNE,IJKS,IJKSE,IPJMK,IJKM,KM,KP,IJMKP
52           INTEGER :: IJKTN,IJKWT,IJKST
53           DOUBLE PRECISION :: We,Ww,Wn,Ws,Wt,Wb
54           DOUBLE PRECISION :: B_NOC
55           DOUBLE PRECISION :: MU_GT_E,MU_GT_W,MU_GT_N,MU_GT_S
56           DOUBLE PRECISION :: MU_GT_T,MU_GT_B,MU_GT_CUT
57           DOUBLE PRECISION :: WW_g
58           INTEGER :: BCV
59           INTEGER :: BCT
60     ! virtual (added) mass
61           DOUBLE PRECISION :: F_vir, ROP_MA
62           DOUBLE PRECISION :: U_se, Usw, Wse, Wsw, Wsn, Wss
63           DOUBLE PRECISION :: Wst, Wsb, Usc,Vsc,Vsn,Vss
64     ! Wall function
65           DOUBLE PRECISION :: W_F_Slip
66     !-----------------------------------------------
67           IF(CG_SAFE_MODE(5)==1) RETURN
68     
69           M = 0
70           IF (.NOT.MOMENTUM_Z_EQ(0)) RETURN
71     !!!$omp  parallel do private( I, J, K, IJK, IJKT, ISV, Sdp, V0, Vpm, Vmt, Vbf, &
72     !!!$omp&  PGT, ROGA, IMJK, IJKP, IMJKP, IJKW, IJKTE, IJKTW, IM, IPJK,  &
73     !!!$omp&  CTE, CTW, SXZB, EPMUOX, VXZA, VXZB, UGT, VCOA, VCOB, IJKE,&
74     !!!$omp&  MUGA, ROPGA, EPGA, LINE)
75           DO IJK = ijkstart3, ijkend3
76              I = I_OF(IJK)
77              J = J_OF(IJK)
78              K = K_OF(IJK)
79              IJKT = TOP_OF(IJK)
80              EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
81              IF (IP_AT_T(IJK)) THEN
82     ! do nothing
83     
84     ! dilute flow
85              ELSE IF (EPGA <= DIL_EP_S) THEN
86     ! do nothing
87     
88              ELSEIF(INTERIOR_CELL_AT(IJK)) THEN
89                 BCV = BC_W_ID(IJK)
90     
91                 IF(BCV > 0 ) THEN
92                    BCT = BC_TYPE_ENUM(BCV)
93                 ELSE
94                    BCT = NONE
95                 ENDIF
96     
97                 SELECT CASE (BCT)
98                    CASE (CG_NSW)
99                       NOC_WG = .TRUE.
100                       WW_g = ZERO
101                       MU_GT_CUT =  (VOL(IJK)*MU_GT(IJK) + VOL(IJKT)*MU_GT(IJKT))/(VOL(IJK) + VOL(IJKT))
102     
103                       IF(.NOT.K_EPSILON) THEN
104                          A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_W_CUT(IJK)/DELH_W(IJK)
105                       ELSE
106                          CALL Wall_Function(IJK,IJK,ONE/DELH_W(IJK),W_F_Slip)
107                          A_M(IJK,0,M) = A_M(IJK,0,M)  - MU_GT_CUT * Area_W_CUT(IJK)*(ONE-W_F_Slip)/DELH_W(IJK)
108                       ENDIF
109     
110                    CASE (CG_FSW)
111                       NOC_WG = .FALSE.
112                       WW_g = ZERO
113                    CASE(CG_PSW)
114                       IF(BC_HW_G(BCV)==UNDEFINED) THEN   ! same as NSW
115                          NOC_WG = .TRUE.
116                          WW_g = BC_WW_G(BCV)
117                          MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJKT)*MU_GT(IJKT))/(VOL(IJK) + VOL(IJKT))
118                          A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_W_CUT(IJK)/DELH_W(IJK)
119                          B_M(IJK,M) = B_M(IJK,M) - MU_GT_CUT * WW_g * Area_W_CUT(IJK)/DELH_W(IJK)
120                       ELSEIF(BC_HW_G(BCV)==ZERO) THEN   ! same as FSW
121                          NOC_WG = .FALSE.
122                          WW_g = ZERO
123                       ELSE                              ! partial slip
124                          NOC_WG = .FALSE.
125                          WW_g = BC_WW_G(BCV)
126                          MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJKT)*MU_GT(IJKT))/(VOL(IJK) + VOL(IJKT))
127                          A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_W_CUT(IJK)*(BC_HW_G(BCV))
128                          B_M(IJK,M) = B_M(IJK,M) - MU_GT_CUT * WW_g * Area_W_CUT(IJK)*(BC_HW_G(BCV))
129     
130                       ENDIF
131                    CASE (NONE,CG_MI)
132                       NOC_WG = .FALSE.
133                 END SELECT
134     
135                 IF(NOC_WG) THEN
136     
137                    IMJK = IM_OF(IJK)
138                    IJMK = JM_OF(IJK)
139                    IJKM = KM_OF(IJK)
140                    IPJK = IP_OF(IJK)
141                    IJPK = JP_OF(IJK)
142                    IJKP = KP_OF(IJK)
143     
144                    We = Theta_We_bar(IJK)  * W_g(IJK)  + Theta_We(IJK)  * W_g(IPJK)
145                    Ww = Theta_We_bar(IMJK) * W_g(IMJK) + Theta_We(IMJK) * W_g(IJK)
146     
147                    Wn = Theta_Wn_bar(IJK)  * W_g(IJK)  + Theta_Wn(IJK)  * W_g(IJPK)
148                    Ws = Theta_Wn_bar(IJMK) * W_g(IJMK) + Theta_Wn(IJMK) * W_g(IJK)
149     
150                    Wt = Theta_Wt_bar(IJK)  * W_g(IJK)  + Theta_Wt(IJK)  * W_g(IJKP)
151                    Wb = Theta_Wt_bar(IJKM) * W_g(IJKM) + Theta_Wt(IJKM) * W_g(IJK)
152     
153                    IJKE = EAST_OF(IJK)
154     
155                    ijkt = top_of(ijk)
156     
157                    IF (WALL_AT(IJK)) THEN
158                       IJKC = IJKT
159                    ELSE
160                       IJKC = IJK
161                    ENDIF
162     
163                    IP = IP1(I)
164                    IM = IM1(I)
165                    IJKN = NORTH_OF(IJK)
166                    IJKNE = EAST_OF(IJKN)
167     
168                    JM = JM1(J)
169                    IPJMK = IP_OF(IJMK)
170                    IJKS = SOUTH_OF(IJK)
171                    IJKSE = EAST_OF(IJKS)
172     
173                    KP = KP1(K)
174                    IJKT = TOP_OF(IJK)
175                    IJKE = EAST_OF(IJK)
176                    IJKP = KP_OF(IJK)
177                    IJKTN = NORTH_OF(IJKT)
178                    IJKTE = EAST_OF(IJKT)
179                    IJKW = WEST_OF(IJK)
180                    IJKWT = TOP_OF(IJKW)
181                    IJKS = SOUTH_OF(IJK)
182                    IJKST = TOP_OF(IJKS)
183     
184                    MU_GT_E = AVG_Z_H(AVG_X_H(MU_GT(IJKC),MU_GT(IJKE),I),&
185                              AVG_X_H(MU_GT(IJKT),MU_GT(IJKTE),I),K)
186     
187                    MU_GT_W = AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJKC),IM),&
188                              AVG_X_H(MU_GT(IJKWT),MU_GT(IJKT),IM),K)
189     
190                    MU_GT_N = AVG_Z_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),&
191                              AVG_Y_H(MU_GT(IJKT),MU_GT(IJKTN),J),K)
192     
193                    MU_GT_S = AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),&
194                              AVG_Y_H(MU_GT(IJKST),MU_GT(IJKT),JM),K)
195     
196                    MU_GT_T = MU_GT(IJKT)
197                    MU_GT_B = MU_GT(IJKC)
198     
199                    B_NOC =     MU_GT_E * Ayz_W(IJK)  * (We-WW_g) * NOC_W_E(IJK)  &
200                            -   MU_GT_W * Ayz_W(IMJK) * (Ww-WW_g) * NOC_W_E(IMJK) &
201                            +   MU_GT_N * Axz_W(IJK)  * (Wn-WW_g) * NOC_W_N(IJK)  &
202                            -   MU_GT_S * Axz_W(IJMK) * (Ws-WW_g) * NOC_W_N(IJMK) &
203                            +   MU_GT_T * Axy_W(IJK)  * (Wt-WW_g) * NOC_W_T(IJK)  *2.0d0&
204                            -   MU_GT_B * Axy_W(IJKM) * (Wb-WW_g) * NOC_W_T(IJKM) *2.0D0
205     
206                    B_M(IJK,M) = B_M(IJK,M)   +  B_NOC
207                 ENDIF
208     
209                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
210     
211     ! BEGIN VIRTUAL MASS SECTION (explicit terms)
212     ! adding transient term  dWs/dt to virtual mass term
213                    F_vir = ZERO
214                    IF(Added_Mass) THEN
215     
216                       F_vir = ( (W_s(IJK,M_AM) - W_sO(IJK,M_AM)) )*ODT*VOL_W(IJK)
217                       I = I_OF(IJK)
218                       J = J_OF(IJK)
219                       K = K_OF(IJK)
220                       IM = I - 1
221                       JM = J - 1
222                       KM = K - 1
223                       IP = I + 1
224                       JP = J + 1
225                       KP = K + 1
226                       IMJK = FUNIJK(IM,J,K)
227                       IJMK = FUNIJK(I,JM,K)
228                       IPJK = FUNIJK(IP,J,K)
229                       IJPK = FUNIJK(I,JP,K)
230                       IJKP = FUNIJK(I,J,KP)
231                       IJKM = FUNIJK(I,J,KM)
232                       IMJKP = KP_OF(IMJK)
233                       IJMKP = KP_OF(IJMK)
234                       IJKE = EAST_OF(IJK)
235     
236     ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
237                       Wse = Theta_We_bar(IJK) * W_s(IJK,M_AM) + Theta_We(IJK) * W_s(IPJK,M_AM)
238                       Wsw = Theta_We_bar(IMJK) * W_s(IMJK,M_AM) + Theta_We(IMJK) * W_s(IJK,M_AM)
239                       U_se = Theta_W_te(IJK) * U_s(IJK,M_AM) + Theta_W_be(IJK) * U_s(IJKP,M_AM)
240                       Usw = Theta_W_te(IMJK) * U_s(IMJK,M_AM) + Theta_W_be(IMJK) * U_s(IMJKP,M_AM)
241                       Usc = (DELX_we(IJK) * Usw + DELX_ww(IJK) * U_se) / (DELX_we(IJK) + DELX_ww(IJK))
242                       Wsn = Theta_Wn_bar(IJK) * W_s(IJK,M_AM) + Theta_Wn(IJK) * W_s(IJPK,M_AM)
243                       Wss = Theta_Wn_bar(IJMK) * W_s(IJMK,M_AM) + Theta_Wn(IJMK) * W_s(IJK,M_AM)
244                       Vsn =  Theta_W_tn(IJK)  * V_s(IJK,M_AM)  + Theta_W_bn(IJK)  * V_s(IJKP,M_AM)
245                       Vss =  Theta_W_tn(IJMK) * V_s(IJMK,M_AM) + Theta_W_bn(IJMK) * V_s(IJMKP,M_AM)
246                       Vsc = (DELY_wn(IJK) * Vss + DELY_ws(IJK) * Vsn) / (DELY_wn(IJK) + DELY_ws(IJK))
247                       Wst = Theta_Wt_bar(IJK)  * W_s(IJK,M_AM)  + Theta_Wt(IJK)  * W_s(IJKP,M_AM)
248                       Wsb = Theta_Wt_bar(IMJK) * W_s(IMJK,M_AM) + Theta_Wt(IMJK) * W_s(IMJKP,M_AM)
249     
250     ! adding convective terms (U dW/dx + V dW/dy + W dW/dz) to virtual mass
251                       F_vir = F_vir +  Usc * (Wse - Wsw)*AYZ(IJK)    + &
252                                        Vsc * (Wsn - Wss)*AXZ(IJK)  + &
253                                        W_s(IJK,M_AM)*(Wst - Wsb)*AXY(IJK)
254     
255     
256                       ROP_MA = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M_AM) + &
257                                 VOL(IJKT)*ROP_g(IJKT)*EP_s(IJKT,M_AM))/(VOL(IJK) + VOL(IJKT))
258     
259                       F_vir = F_vir * Cv * ROP_MA
260                       B_M(IJK,M) = B_M(IJK,M) - F_vir ! explicit part of virtual mass force
261     
262                    ENDIF
263     ! END VIRTUAL MASS SECTION
264     
265                 ENDIF
266              ENDIF
267           ENDDO
268     
269           RETURN
270     
271         CONTAINS
272     
273           INCLUDE 'functions.inc'
274     
275           END SUBROUTINE CG_SOURCE_W_G
276     
277     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
278     !                                                                      C
279     !  Module name: CG_SOURCE_W_g_BC(A_m, B_m, IER)                        C
280     !  Purpose: Determine contribution of cut-cell to source terms         C
281     !  for W_g momentum eq.                                                C
282     !                                                                      C
283     !                                                                      C
284     !  Author: Jeff Dietiker                              Date: 01-MAY-09  C
285     !                                                                      C
286     !                                                                      C
287     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
288           SUBROUTINE CG_SOURCE_W_G_BC(A_M, B_M)
289     
290     ! Modules
291     !-----------------------------------------------
292           USE bc
293           USE compar
294           USE fldvar
295           USE fun_avg
296           USE indices
297           USE param
298           USE param1
299           USE physprop
300           USE run
301           USE toleranc
302           USE visc_g
303           USE cutcell
304           USE quadric
305           IMPLICIT NONE
306     
307     ! Dummy arguments
308     !-----------------------------------------------
309     ! Septadiagonal matrix A_m
310           DOUBLE PRECISION :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
311     ! Vector b_m
312           DOUBLE PRECISION :: B_m(DIMENSION_3, 0:DIMENSION_M)
313     
314     ! Local variables
315     !-----------------------------------------------
316     ! Indices
317           INTEGER :: IJK, IJKB
318     ! Solids phase
319           INTEGER :: M
320           INTEGER :: BCV
321           INTEGER :: BCT
322     
323     !-----------------------------------------------
324     
325           IF(CG_SAFE_MODE(5)==1) RETURN
326           M = 0
327     
328           DO IJK = ijkstart3, ijkend3
329     
330              BCV = BC_W_ID(IJK)
331     
332              IF(BCV > 0 ) THEN
333                 BCT = BC_TYPE_ENUM(BCV)
334              ELSE
335                 BCT = NONE
336              ENDIF
337     
338              SELECT CASE (BCT)
339     
340                 CASE (CG_NSW)
341     
342                    IF(WALL_W_AT(IJK)) THEN
343     
344                       A_M(IJK,east,M) = ZERO
345                       A_M(IJK,west,M) = ZERO
346                       A_M(IJK,north,M) = ZERO
347                       A_M(IJK,south,M) = ZERO
348                       A_M(IJK,top,M) = ZERO
349                       A_M(IJK,bottom,M) = ZERO
350                       A_M(IJK,0,M) = -ONE
351     
352                       B_M(IJK,M) = ZERO
353     
354                    ENDIF
355     
356                 CASE (CG_FSW)
357     
358                    IF(WALL_W_AT(IJK)) THEN
359     
360                       A_M(IJK,east,M) = ZERO
361                       A_M(IJK,west,M) = ZERO
362                       A_M(IJK,north,M) = ZERO
363                       A_M(IJK,south,M) = ZERO
364                       A_M(IJK,top,M) = ZERO
365                       A_M(IJK,bottom,M) = ZERO
366                       A_M(IJK,0,M) = -ONE
367     
368     !                  B_M(IJK,M) = - W_g(W_MASTER_OF(IJK))  ! Velocity of master node
369     
370                       B_M(IJK,M) = ZERO
371     
372                       IF(DABS(NORMAL_W(IJK,3))/=ONE) THEN
373     
374                          IF (W_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
375                             A_M(IJK,east,M) = ONE
376                          ELSEIF (W_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
377                             A_M(IJK,west,M) = ONE
378                          ELSEIF (W_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
379                             A_M(IJK,north,M) = ONE
380                          ELSEIF (W_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
381                             A_M(IJK,south,M) = ONE
382                          ELSEIF (W_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
383                             A_M(IJK,top,M) = ONE
384                          ELSEIF (W_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
385                             A_M(IJK,bottom,M) = ONE
386                          ENDIF
387     
388                       ENDIF
389     
390                    ENDIF
391     
392                 CASE (CG_PSW)
393     
394                    IF(WALL_W_AT(IJK)) THEN
395     
396                       A_M(IJK,east,M) = ZERO
397                       A_M(IJK,west,M) = ZERO
398                       A_M(IJK,north,M) = ZERO
399                       A_M(IJK,south,M) = ZERO
400                       A_M(IJK,top,M) = ZERO
401                       A_M(IJK,bottom,M) = ZERO
402                       A_M(IJK,0,M) = -ONE
403     
404     
405                       IF(BC_HW_G(BCV)==UNDEFINED) THEN   ! same as NSW
406                          B_M(IJK,M) = -BC_WW_G(BCV)
407                       ELSEIF(BC_HW_G(BCV)==ZERO) THEN   ! same as FSW
408                          B_M(IJK,M) = ZERO
409     
410                          IF(DABS(NORMAL_W(IJK,3))/=ONE) THEN
411     
412                             IF (W_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
413                                A_M(IJK,east,M) = ONE
414                             ELSEIF (W_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
415                                A_M(IJK,west,M) = ONE
416                             ELSEIF (W_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
417                                A_M(IJK,north,M) = ONE
418                             ELSEIF (W_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
419                                A_M(IJK,south,M) = ONE
420                             ELSEIF (W_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
421                                A_M(IJK,top,M) = ONE
422                             ELSEIF (W_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
423                                A_M(IJK,bottom,M) = ONE
424                             ENDIF
425     
426                          ENDIF
427     
428                       ELSE                              ! partial slip
429     
430                       ENDIF
431     
432                    ENDIF
433     
434     
435                 CASE (CG_MI)
436     
437                    A_M(IJK,east,M) = ZERO
438                    A_M(IJK,west,M) = ZERO
439                    A_M(IJK,north,M) = ZERO
440                    A_M(IJK,south,M) = ZERO
441                    A_M(IJK,top,M) = ZERO
442                    A_M(IJK,bottom,M) = ZERO
443                    A_M(IJK,0,M) = -ONE
444     
445                    IF(BC_W_g(BCV)/=UNDEFINED) THEN
446                       B_M(IJK,M) = - BC_W_g(BCV)
447                    ELSE
448                       B_M(IJK,M) = - BC_VELMAG_g(BCV)*NORMAL_W(IJK,3)
449                    ENDIF
450     
451     
452                    IJKB = BOTTOM_OF(IJK)
453                    IF(FLUID_AT(IJKB)) THEN
454     
455                       A_M(IJKB,east,M) = ZERO
456                       A_M(IJKB,west,M) = ZERO
457                       A_M(IJKB,north,M) = ZERO
458                       A_M(IJKB,south,M) = ZERO
459                       A_M(IJKB,top,M) = ZERO
460                       A_M(IJKB,bottom,M) = ZERO
461                       A_M(IJKB,0,M) = -ONE
462     
463                       IF(BC_W_g(BCV)/=UNDEFINED) THEN
464                          B_M(IJKB,M) = - BC_W_g(BCV)
465                       ELSE
466                          B_M(IJKB,M) = - BC_VELMAG_g(BCV)*NORMAL_W(IJK,3)
467                       ENDIF
468     
469     
470                    ENDIF
471     
472                 CASE (CG_PO)
473     
474                    A_M(IJK,east,M) = ZERO
475                    A_M(IJK,west,M) = ZERO
476                    A_M(IJK,north,M) = ZERO
477                    A_M(IJK,south,M) = ZERO
478                    A_M(IJK,top,M) = ZERO
479                    A_M(IJK,bottom,M) = ZERO
480                    A_M(IJK,0,M) = -ONE
481                    B_M(IJK,M) = ZERO
482     
483                    IJKB = BOTTOM_OF(IJK)
484                    IF(FLUID_AT(IJKB)) THEN
485     
486                       A_M(IJK,bottom,M) = ONE
487                       A_M(IJK,0,M) = -ONE
488     
489                    ENDIF
490     
491              END SELECT
492     
493              BCV = BC_ID(IJK)
494     
495              IF(BCV > 0 ) THEN
496                 BCT = BC_TYPE_ENUM(BCV)
497              ELSE
498                 BCT = NONE
499              ENDIF
500     
501              SELECT CASE (BCT)
502     
503                 CASE (CG_MI)
504     
505                    A_M(IJK,east,M) = ZERO
506                    A_M(IJK,west,M) = ZERO
507                    A_M(IJK,north,M) = ZERO
508                    A_M(IJK,south,M) = ZERO
509                    A_M(IJK,top,M) = ZERO
510                    A_M(IJK,bottom,M) = ZERO
511                    A_M(IJK,0,M) = -ONE
512     
513                    IF(BC_W_g(BCV)/=UNDEFINED) THEN
514                       B_M(IJK,M) = - BC_W_g(BCV)
515                    ELSE
516                       B_M(IJK,M) = - BC_VELMAG_g(BCV)*NORMAL_S(IJK,3)
517                    ENDIF
518     
519     
520                    IJKB = BOTTOM_OF(IJK)
521                    IF(FLUID_AT(IJKB)) THEN
522     
523                       A_M(IJKB,east,M) = ZERO
524                       A_M(IJKB,west,M) = ZERO
525                       A_M(IJKB,north,M) = ZERO
526                       A_M(IJKB,south,M) = ZERO
527                       A_M(IJKB,top,M) = ZERO
528                       A_M(IJKB,bottom,M) = ZERO
529                       A_M(IJKB,0,M) = -ONE
530     
531                       IF(BC_W_g(BCV)/=UNDEFINED) THEN
532                          B_M(IJKB,M) = - BC_W_g(BCV)
533                       ELSE
534                          B_M(IJKB,M) = - BC_VELMAG_g(BCV)*NORMAL_S(IJK,3)
535                       ENDIF
536     
537     
538                    ENDIF
539     
540                 CASE (CG_PO)
541     
542                    A_M(IJK,east,M) = ZERO
543                    A_M(IJK,west,M) = ZERO
544                    A_M(IJK,north,M) = ZERO
545                    A_M(IJK,south,M) = ZERO
546                    A_M(IJK,top,M) = ZERO
547                    A_M(IJK,bottom,M) = ZERO
548                    A_M(IJK,0,M) = -ONE
549                    B_M(IJK,M) = ZERO
550     
551                    IJKB = BOTTOM_OF(IJK)
552                    IF(FLUID_AT(IJKB)) THEN
553     
554                       A_M(IJK,bottom,M) = ONE
555                       A_M(IJK,0,M) = -ONE
556     
557                    ENDIF
558     
559              END SELECT
560     
561     
562           ENDDO
563     
564           RETURN
565     
566     
567         CONTAINS
568     
569           INCLUDE 'functions.inc'
570     
571           END SUBROUTINE CG_SOURCE_W_G_BC
572