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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: SOURCE_V_g(A_m, B_m, IER)                              C
4     !  Purpose: Determine contribution of cut-cell to source terms         C
5     !  for V_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_V_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: do_k, vol, vol_v
21           use geometry, only: axy, ayz, axz
22           use geometry, only: ayz_v, axz_v, axy_v
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, IJKN
45     ! Phase index
46           INTEGER :: M
47     ! Average volume fraction
48           DOUBLE PRECISION EPGA
49           INTEGER :: IM,JM,IP,JP,KM,KP
50           INTEGER :: IMJK,IPJK,IJMK,IJPK,IJKP,IJKM
51           INTEGER :: IJKC,IJKE,IJKNE,IJKW,IJKWN,IMJPK,IJPKM
52           INTEGER :: IJKT,IJKTN,IJKB,IJKBN
53           DOUBLE PRECISION :: Vn,Vs,Ve,Vw, Vt,Vb
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 :: VW_g
58           INTEGER :: BCV
59           INTEGER :: BCT
60     !  virtual (added) mass
61           DOUBLE PRECISION :: F_vir, ROP_MA
62           DOUBLE PRECISION :: Vsn, Vss, U_se, Usc, Usw, Vse, Vsw
63           DOUBLE PRECISION :: Wst, Wsb, Wsc, Vst, Vsb
64     ! Wall function
65           DOUBLE PRECISION :: W_F_Slip
66     !-----------------------------------------------
67     
68     
69           IF(CG_SAFE_MODE(4)==1) RETURN
70     
71           M = 0
72           IF (.NOT.MOMENTUM_Y_EQ(0)) RETURN
73     
74     !!!$omp  parallel do private( I, J, K, IJK, IJKN, ISV, Sdp, V0, Vpm, Vmt, Vbf, &
75     !!!$omp&  PGN, ROGA, MUGA, ROPGA, EPGA,VSH_n,VSH_s,VSH_e,VSH_w,&
76     !!!$omp&  VSH_p,Source_conv, SRT ) &
77     !!!$omp&  schedule(static)
78           DO IJK = ijkstart3, ijkend3
79              I = I_OF(IJK)
80              J = J_OF(IJK)
81              K = K_OF(IJK)
82              IJKN = NORTH_OF(IJK)
83              EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
84              IF (IP_AT_N(IJK)) THEN
85     ! do nothing
86     
87     ! dilute flow
88              ELSEIF (EPGA <= DIL_EP_S) THEN
89     ! do nothing
90     
91              ELSEIF(INTERIOR_CELL_AT(IJK)) THEN
92     
93                 BCV = BC_V_ID(IJK)
94     
95                 IF(BCV > 0 ) THEN
96                    BCT = BC_TYPE_ENUM(BCV)
97                 ELSE
98                    BCT = NONE
99                 ENDIF
100     
101                 SELECT CASE (BCT)
102                    CASE (CG_NSW)
103                       NOC_VG = .TRUE.
104                       VW_g = ZERO
105                       MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJKN)*MU_GT(IJKN))/(VOL(IJK) + VOL(IJKN))
106     
107                       IF(.NOT.K_EPSILON) THEN
108                          A_M(IJK,0,M) = A_M(IJK,0,M)  - MU_GT_CUT * Area_V_CUT(IJK)/DELH_V(IJK)
109                       ELSE
110                          CALL Wall_Function(IJK,IJK,ONE/DELH_V(IJK),W_F_Slip)
111                          A_M(IJK,0,M) = A_M(IJK,0,M)  - MU_GT_CUT * Area_V_CUT(IJK)*(ONE-W_F_Slip)/DELH_V(IJK)
112                       ENDIF
113                    CASE (CG_FSW)
114                       NOC_VG = .FALSE.
115                       VW_g = ZERO
116                    CASE(CG_PSW)
117                       IF(BC_HW_G(BCV)==UNDEFINED) THEN   ! same as NSW
118                          NOC_VG = .TRUE.
119                          VW_g = BC_VW_G(BCV)
120                          MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJKN)*MU_GT(IJKN))/(VOL(IJK) + VOL(IJKN))
121                          A_M(IJK,0,M) = A_M(IJK,0,M)  - MU_GT_CUT * Area_V_CUT(IJK)/DELH_V(IJK)
122                          B_M(IJK,M) = B_M(IJK,M) - MU_GT_CUT * VW_g * Area_V_CUT(IJK)/DELH_V(IJK)
123                       ELSEIF(BC_HW_G(BCV)==ZERO) THEN   ! same as FSW
124                          NOC_VG = .FALSE.
125                          VW_g = ZERO
126                       ELSE                              ! partial slip
127                          NOC_VG = .FALSE.
128                          VW_g = BC_VW_G(BCV)
129                          MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJKN)*MU_GT(IJKN))/(VOL(IJK) + VOL(IJKN))
130                          A_M(IJK,0,M) = A_M(IJK,0,M)  - MU_GT_CUT * Area_V_CUT(IJK)*(BC_HW_G(BCV))
131                          B_M(IJK,M) = B_M(IJK,M) - MU_GT_CUT * VW_g * Area_V_CUT(IJK)*(BC_HW_G(BCV))
132                       ENDIF
133                    CASE (NONE, CG_MI)
134                       NOC_VG = .FALSE.
135                 END SELECT
136     
137     
138                 IF(NOC_VG) THEN
139                    IMJK = IM_OF(IJK)
140                    IJMK = JM_OF(IJK)
141                    IJKM = KM_OF(IJK)
142                    IPJK = IP_OF(IJK)
143                    IJPK = JP_OF(IJK)
144                    IJKP = KP_OF(IJK)
145     
146                    Vn = Theta_Vn_bar(IJK)  * V_g(IJK)  + Theta_Vn(IJK)  * V_g(IJPK)
147                    Vs = Theta_Vn_bar(IJMK) * V_g(IJMK) + Theta_Vn(IJMK) * V_g(IJK)
148                    Ve = Theta_Ve_bar(IJK)  * V_g(IJK)  + Theta_Ve(IJK)  * V_g(IPJK)
149                    Vw = Theta_Ve_bar(IMJK) * V_g(IMJK) + Theta_Ve(IMJK) * V_g(IJK)
150     
151                    IJKN = NORTH_OF(IJK)
152                    IF (WALL_AT(IJK)) THEN
153                       IJKC = IJKN
154                    ELSE
155                       IJKC = IJK
156                    ENDIF
157                    JP = JP1(J)
158                    IJKE = EAST_OF(IJK)
159                    IJKNE = EAST_OF(IJKN)
160                    IM = IM1(I)
161                    IJKW = WEST_OF(IJK)
162                    IJKWN = NORTH_OF(IJKW)
163                    IMJPK = JP_OF(IMJK)
164                    KM = KM1(K)
165                    IJKT = TOP_OF(IJK)
166                    IJKTN = NORTH_OF(IJKT)
167                    IJKB = BOTTOM_OF(IJK)
168                    IJKBN = NORTH_OF(IJKB)
169     
170     
171                    MU_GT_E = AVG_Y_H(AVG_X_H(MU_GT(IJKC),MU_GT(IJKE),I),&
172                              AVG_X_H(MU_GT(IJKN),MU_GT(IJKNE),I),J)
173                    MU_GT_W = AVG_Y_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJKC),IM),&
174                              AVG_X_H(MU_GT(IJKWN),MU_GT(IJKN),IM),J)
175     
176                    MU_GT_N = MU_GT(IJKN)
177                    MU_GT_S = MU_GT(IJKC)
178     
179     
180                    B_NOC =     MU_GT_N * Axz_V(IJK)  * (Vn-VW_g) * NOC_V_N(IJK)   *2.0d0&
181                            -   MU_GT_S * Axz_V(IJMK) * (Vs-VW_g) * NOC_V_N(IJMK)  *2.0d0&
182                            +   MU_GT_E * Ayz_V(IJK)  * (Ve-VW_g) * NOC_V_E(IJK)   &
183                            -   MU_GT_W * Ayz_V(IMJK) * (Vw-VW_g) * NOC_V_E(IMJK)
184     
185                    IF(DO_K) THEN
186     
187                       Vt = Theta_Vt_bar(IJK)  * V_g(IJK)  + Theta_Vt(IJK)  * V_g(IJKP)
188                       Vb = Theta_Vt_bar(IJKM) * V_g(IJKM) + Theta_Vt(IJKM) * V_g(IJK)
189     
190                       MU_GT_T = AVG_Y_H(AVG_Z_H(MU_GT(IJKC),MU_GT(IJKT),K),&
191                                 AVG_Z_H(MU_GT(IJKN),MU_GT(IJKTN),K),J)
192     
193                       MU_GT_B = AVG_Y_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJKC),KM),&
194                                 AVG_Z_H(MU_GT(IJKBN),MU_GT(IJKN),KM),J)
195     
196                       B_NOC = B_NOC + MU_GT_T * Axy_V(IJK)  * (Vt-VW_g) * NOC_V_T(IJK)   &
197                                     - MU_GT_B * Axy_V(IJKM) * (Vb-VW_g) * NOC_V_T(IJKM)
198     
199                    ENDIF
200                       B_M(IJK,M) = B_M(IJK,M) + B_NOC
201     
202                 ENDIF
203     
204     
205                 IF(CUT_V_TREATMENT_AT(IJK)) THEN
206     
207     ! BEGIN VIRTUAL MASS SECTION (explicit terms)
208     ! adding transient term dVs/dt to virtual mass term
209                    F_vir = ZERO
210                    IF(Added_Mass) THEN
211                       F_vir = ( (V_s(IJK,M_AM) - V_sO(IJK,M_AM)) )*ODT*VOL_V(IJK)
212                       IM = I - 1
213                       JM = J - 1
214                       KM = K - 1
215                       IP = I + 1
216                       JP = J + 1
217                       KP = K + 1
218                       IMJK = FUNIJK(IM,J,K)
219                       IJMK = FUNIJK(I,JM,K)
220                       IPJK = FUNIJK(IP,J,K)
221                       IJPK = FUNIJK(I,JP,K)
222                       IJKP = FUNIJK(I,J,KP)
223                       IJKM = FUNIJK(I,J,KM)
224                       IMJPK = IM_OF(IJPK)
225                       IJKN = NORTH_OF(IJK)
226     
227     ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
228                       Vse = Theta_Ve_bar(IJK)  * V_s(IJK,M_AM)  + Theta_Ve(IJK)  * V_s(IPJK,M_AM)
229                       Vsw = Theta_Ve_bar(IMJK) * V_s(IMJK,M_AM) + Theta_Ve(IMJK) * V_s(IJK,M_AM)
230                       U_se =  Theta_V_ne(IJK)  * U_s(IJK,M_AM)  + Theta_V_se(IJK)  * U_s(IJPK,M_AM)
231                       Usw =  Theta_V_ne(IMJK) * U_s(IMJK,M_AM) + Theta_V_se(IMJK) * U_s(IMJPK,M_AM)
232                       Usc = (DELX_ve(IJK) * Usw + DELX_vw(IJK) * U_se) / (DELX_ve(IJK) + DELX_vw(IJK))
233                       Vsn = Theta_Vn_bar(IJK)  * V_s(IJK,M_AM)  + Theta_Vn(IJK)  * V_s(IJPK,M_AM)
234                       Vss = Theta_Vn_bar(IJMK) * V_s(IJMK,M_AM) + Theta_Vn(IJMK) * V_s(IJK,M_AM)
235     
236                       IF(DO_K) THEN
237                          IJPKM = KM_OF(IJPK)
238                          Vst = Theta_Vt_bar(IJK)  * V_s(IJK,M_AM)  + Theta_Vt(IJK)  * V_s(IJKP,M_AM)
239                          Vsb = Theta_Vt_bar(IJKM) * V_s(IJKM,M_AM) + Theta_Vt(IJKM) * V_s(IJK,M_AM)
240                          Wst = Theta_V_nt(IJK)  * W_s(IJK,M_AM)  + Theta_V_st(IJK)  * W_s(IJPK,M_AM)
241                          Wsb = Theta_V_nt(IJKM) * W_s(IJKM,M_AM) + Theta_V_st(IJKM) * W_s(IJPKM,M_AM)
242                          Wsc = (DELZ_vt(IJK) * Wsb + DELZ_vb(IJK) * Wst) / (DELZ_vt(IJK) + DELZ_vb(IJK))
243                          F_vir = F_vir +  Wsc* (Vst - Vsb)*AXY(IJK)
244                       ENDIF
245     
246     ! adding convective terms (U dV/dx + V dV/dy) to virtual mass; W dV/dz added above.
247                       F_vir = F_vir + V_s(IJK,M_AM)*(Vsn - Vss)*AXZ(IJK) + &
248                           Usc*(Vse - Vsw)*AYZ(IJK)
249     
250                       ROP_MA  = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M_AM)  + &
251                                  VOL(IJKN)*ROP_g(IJKN)*EP_s(IJKN,M_AM) )/(VOL(IJK) + VOL(IJKN))
252                       F_vir = F_vir * Cv * ROP_MA
253                       B_M(IJK,M) = B_M(IJK,M) - F_vir ! adding explicit-part of virtual mass force.
254                    ENDIF
255     ! END VIRTUAL MASS SECTION
256                 ENDIF
257     
258              ENDIF
259           ENDDO
260     
261           RETURN
262     
263         CONTAINS
264     
265           INCLUDE 'functions.inc'
266     
267           END SUBROUTINE CG_SOURCE_V_G
268     
269     
270     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
271     !                                                                      C
272     !  Module name: CG_SOURCE_V_g_BC(A_m, B_m, IER)                        C
273     !  Purpose: Determine contribution of cut-cell to source terms         C
274     !  for V_g momentum eq.                                                C
275     !                                                                      C
276     !                                                                      C
277     !  Author: Jeff Dietiker                              Date: 01-MAY-09  C
278     !                                                                      C
279     !                                                                      C
280     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
281           SUBROUTINE CG_SOURCE_V_G_BC(A_M, B_M)
282     
283     ! Modules 
284     !-----------------------------------------------
285           USE bc
286           USE compar
287           USE fldvar
288           USE fun_avg
289           use geometry, only: do_k, vol, vol_v
290           use geometry, only: ayz, axy, ayz_v, axy_v
291           use geometry, only: odx_e, ody_n
292           USE indices
293           USE param
294           USE param1
295           USE physprop
296           USE run
297           USE toleranc
298           USE visc_g
299           USE cutcell
300           USE quadric
301           IMPLICIT NONE
302     
303     ! Dummy arguments
304     !-----------------------------------------------
305     ! Septadiagonal matrix A_m
306           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
307     ! Vector b_m
308           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
309     
310     ! Local variables
311     !-----------------------------------------------
312     ! Indices
313           INTEGER :: I,  J, K, JM, IJK
314           INTEGER :: IM, KM, IJKS, IMJK
315     ! Solids phase
316           INTEGER :: M
317     ! Turbulent shear stress
318           DOUBLE PRECISION :: W_F_Slip
319           INTEGER :: BCV
320           INTEGER :: BCT
321     
322     !-----------------------------------------------
323     
324           IF(CG_SAFE_MODE(4)==1) RETURN
325     
326           M = 0
327     
328           DO IJK = ijkstart3, ijkend3
329     
330              BCV = BC_V_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_V_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                       IF(K_EPSILON) THEN
355     
356                          I = I_OF(IJK)
357                          J = J_OF(IJK)
358                          K = K_OF(IJK)
359     
360                          IM = I - 1
361                          JM = J - 1
362                          KM = K - 1
363     
364                       IF(DABS(NORMAL_V(IJK,2))/=ONE) THEN
365     
366                          IF (V_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
367                             CALL Wall_Function(IJK,EAST_OF(IJK),ODX_E(I),W_F_Slip)
368                             A_M(IJK,east,M) = W_F_Slip
369                             A_M(IJK,0,M) = -ONE
370                          ELSEIF (V_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
371                             CALL Wall_Function(IJK,WEST_OF(IJK),ODX_E(IM),W_F_Slip)
372                             A_M(IJK,west,M) = W_F_Slip
373                             A_M(IJK,0,M) = -ONE
374                          ELSEIF (V_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
375                             A_M(IJK,north,M) = ONE
376                          ELSEIF (V_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
377                             A_M(IJK,south,M) = ONE
378                          ELSEIF (V_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
379                             A_M(IJK,top,M) = ONE
380                          ELSEIF (V_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
381                             A_M(IJK,bottom,M) = ONE
382                          ENDIF
383     
384                       ENDIF
385     
386                       ENDIF
387     
388     
389                    ENDIF
390     
391                 CASE (CG_FSW)
392     
393                    IF(WALL_V_AT(IJK)) THEN
394     
395                       A_M(IJK,east,M) = ZERO
396                       A_M(IJK,west,M) = ZERO
397                       A_M(IJK,north,M) = ZERO
398                       A_M(IJK,south,M) = ZERO
399                       A_M(IJK,top,M) = ZERO
400                       A_M(IJK,bottom,M) = ZERO
401                       A_M(IJK,0,M) = -ONE
402     
403     !                  B_M(IJK,M) = - V_g(V_MASTER_OF(IJK))    ! Velocity of master node
404     
405                       B_M(IJK,M) = ZERO
406     
407                       IF(DABS(NORMAL_V(IJK,2))/=ONE) THEN
408     
409                          IF (V_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
410                             A_M(IJK,east,M) = ONE
411                          ELSEIF (V_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
412                             A_M(IJK,west,M) = ONE
413                          ELSEIF (V_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
414                             A_M(IJK,north,M) = ONE
415                          ELSEIF (V_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
416                             A_M(IJK,south,M) = ONE
417                          ELSEIF (V_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
418                             A_M(IJK,top,M) = ONE
419                          ELSEIF (V_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
420                             A_M(IJK,bottom,M) = ONE
421                          ENDIF
422     
423                       ENDIF
424     
425                    ENDIF
426     
427     
428                 CASE (CG_PSW)
429     
430                    IF(WALL_V_AT(IJK)) THEN
431     
432                       A_M(IJK,east,M) = ZERO
433                       A_M(IJK,west,M) = ZERO
434                       A_M(IJK,north,M) = ZERO
435                       A_M(IJK,south,M) = ZERO
436                       A_M(IJK,top,M) = ZERO
437                       A_M(IJK,bottom,M) = ZERO
438                       A_M(IJK,0,M) = -ONE
439     
440     
441                       IF(BC_HW_G(BCV)==UNDEFINED) THEN   ! same as NSW
442                          B_M(IJK,M) = -BC_VW_G(BCV)
443                       ELSEIF(BC_HW_G(BCV)==ZERO) THEN   ! same as FSW
444                          B_M(IJK,M) = ZERO
445     
446                          IF(DABS(NORMAL_V(IJK,2))/=ONE) THEN
447     
448                             IF (V_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
449                                A_M(IJK,east,M) = ONE
450                             ELSEIF (V_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
451                                A_M(IJK,west,M) = ONE
452                             ELSEIF (V_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
453                                A_M(IJK,north,M) = ONE
454                             ELSEIF (V_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
455                                A_M(IJK,south,M) = ONE
456                             ELSEIF (V_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
457                                A_M(IJK,top,M) = ONE
458                             ELSEIF (V_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
459                                A_M(IJK,bottom,M) = ONE
460                             ENDIF
461     
462                          ENDIF
463     
464                       ELSE                              ! partial slip
465     
466                          B_M(IJK,M) = ZERO
467     
468                          I = I_OF(IJK)
469                          J = J_OF(IJK)
470                          K = K_OF(IJK)
471     
472                          IM = I - 1
473                          JM = J - 1
474                          KM = K - 1
475     
476                          IMJK = FUNIJK(IM,J,K)
477     
478                          IF(DABS(NORMAL_V(IJK,2))/=ONE) THEN
479     
480                             IF (V_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
481     !                          A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ODX_E(I))
482     !                          A_M(IJK,east,M) = -(HALF*BC_HW_G(BCV)-ODX_E(I))
483     !                          B_M(IJK,M) = -BC_HW_G(BCV)*BC_VW_G(BCV)
484     
485                                A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ONEoDX_E_V(IJK))
486                                A_M(IJK,east,M) = -(HALF*BC_HW_G(BCV)-ONEoDX_E_V(IJK))
487                                B_M(IJK,M) = -BC_HW_G(BCV)*BC_VW_G(BCV)
488     
489                             ELSEIF (V_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
490     !                          A_M(IJK,west,M) = -(HALF*BC_HW_G(BCV)-ODX_E(IM))
491     !                          A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ODX_E(IM))
492     !                          B_M(IJK,M) = -BC_HW_G(BCV)*BC_VW_G(BCV)
493     
494                                A_M(IJK,west,M) = -(HALF*BC_HW_G(BCV)-ONEoDX_E_V(IMJK))
495                                A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ONEoDX_E_V(IMJK))
496                                B_M(IJK,M) = -BC_HW_G(BCV)*BC_VW_G(BCV)
497     
498     
499                             ELSEIF (V_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
500     
501                                   A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ODY_N(J))
502                                   A_M(IJK,north,M) = -(HALF*BC_HW_G(BCV)-ODY_N(J))
503                                   B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
504     
505     !                           print*,'vg master at north'
506                             ELSEIF (V_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
507     
508                                   A_M(IJK,south,M) = -(HALF*BC_HW_G(BCV)-ODY_N(JM))
509                                   A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ODY_N(JM))
510                                   B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
511     
512     !                           print*,'vg master at south'
513                             ELSEIF (V_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
514                                A_M(IJK,top,M) = ONE
515                             ELSEIF (V_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
516                                A_M(IJK,bottom,M) = ONE
517                             ENDIF
518     
519                          ENDIF
520                       ENDIF
521     
522                    ENDIF
523     
524     
525     
526                 CASE (CG_MI)
527     
528                    A_M(IJK,east,M) = ZERO
529                    A_M(IJK,west,M) = ZERO
530                    A_M(IJK,north,M) = ZERO
531                    A_M(IJK,south,M) = ZERO
532                    A_M(IJK,top,M) = ZERO
533                    A_M(IJK,bottom,M) = ZERO
534                    A_M(IJK,0,M) = -ONE
535     
536                    IF(BC_V_g(BCV)/=UNDEFINED) THEN
537                       B_M(IJK,M) = - BC_V_g(BCV)
538                    ELSE
539                       B_M(IJK,M) = - BC_VELMAG_g(BCV)*NORMAL_V(IJK,2)
540                    ENDIF
541     
542     
543                    IJKS = SOUTH_OF(IJK)
544                    IF(FLUID_AT(IJKS)) THEN
545     
546                       A_M(IJKS,east,M) = ZERO
547                       A_M(IJKS,west,M) = ZERO
548                       A_M(IJKS,north,M) = ZERO
549                       A_M(IJKS,south,M) = ZERO
550                       A_M(IJKS,top,M) = ZERO
551                       A_M(IJKS,bottom,M) = ZERO
552                       A_M(IJKS,0,M) = -ONE
553     
554                       IF(BC_V_g(BCV)/=UNDEFINED) THEN
555                          B_M(IJKS,M) = - BC_V_g(BCV)
556                       ELSE
557                          B_M(IJKS,M) = - BC_VELMAG_g(BCV)*NORMAL_V(IJK,2)
558                       ENDIF
559     
560     
561                    ENDIF
562     
563                 CASE (CG_PO)
564     
565                    A_M(IJK,east,M) = ZERO
566                    A_M(IJK,west,M) = ZERO
567                    A_M(IJK,north,M) = ZERO
568                    A_M(IJK,south,M) = ZERO
569                    A_M(IJK,top,M) = ZERO
570                    A_M(IJK,bottom,M) = ZERO
571                    A_M(IJK,0,M) = -ONE
572                    B_M(IJK,M) = ZERO
573     
574                    IJKS = SOUTH_OF(IJK)
575                    IF(FLUID_AT(IJKS)) THEN
576     
577                       A_M(IJK,south,M) = ONE
578                       A_M(IJK,0,M) = -ONE
579                       B_M(IJK,M) = ZERO
580     
581                    ENDIF
582     
583     
584              END SELECT
585     
586     
587              BCV = BC_ID(IJK)
588     
589              IF(BCV > 0 ) THEN
590                 BCT = BC_TYPE_ENUM(BCV)
591              ELSE
592                 BCT = NONE
593              ENDIF
594     
595              SELECT CASE (BCT)
596     
597                 CASE (CG_MI)
598     
599                    A_M(IJK,east,M) = ZERO
600                    A_M(IJK,west,M) = ZERO
601                    A_M(IJK,north,M) = ZERO
602                    A_M(IJK,south,M) = ZERO
603                    A_M(IJK,top,M) = ZERO
604                    A_M(IJK,bottom,M) = ZERO
605                    A_M(IJK,0,M) = -ONE
606     
607                    IF(BC_V_g(BCV)/=UNDEFINED) THEN
608                       B_M(IJK,M) = - BC_V_g(BCV)
609                    ELSE
610                       B_M(IJK,M) = - BC_VELMAG_g(BCV)*NORMAL_S(IJK,2)
611                    ENDIF
612     
613     
614                    IJKS = SOUTH_OF(IJK)
615                    IF(FLUID_AT(IJKS)) THEN
616     
617                       A_M(IJKS,east,M) = ZERO
618                       A_M(IJKS,west,M) = ZERO
619                       A_M(IJKS,north,M) = ZERO
620                       A_M(IJKS,south,M) = ZERO
621                       A_M(IJKS,top,M) = ZERO
622                       A_M(IJKS,bottom,M) = ZERO
623                       A_M(IJKS,0,M) = -ONE
624     
625                       IF(BC_V_g(BCV)/=UNDEFINED) THEN
626                          B_M(IJKS,M) = - BC_V_g(BCV)
627                       ELSE
628                          B_M(IJKS,M) = - BC_VELMAG_g(BCV)*NORMAL_S(IJK,2)
629                       ENDIF
630     
631     
632                    ENDIF
633     
634                 CASE (CG_PO)
635     
636                    A_M(IJK,east,M) = ZERO
637                    A_M(IJK,west,M) = ZERO
638                    A_M(IJK,north,M) = ZERO
639                    A_M(IJK,south,M) = ZERO
640                    A_M(IJK,top,M) = ZERO
641                    A_M(IJK,bottom,M) = ZERO
642                    A_M(IJK,0,M) = -ONE
643                    B_M(IJK,M) = ZERO
644     
645                    IJKS = SOUTH_OF(IJK)
646                    IF(FLUID_AT(IJKS)) THEN
647     
648                       A_M(IJK,south,M) = ONE
649                       A_M(IJK,0,M) = -ONE
650     
651                    ENDIF
652     
653              END SELECT
654     
655           ENDDO
656     
657           RETURN
658     
659         CONTAINS
660     
661           INCLUDE 'functions.inc'
662     
663           END SUBROUTINE CG_SOURCE_V_G_BC
664     
665