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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CG_SOURCE_V_s                                           C
4     !  Purpose: Determine contribution of cut-cell to source terms         C
5     !           for V_s momentum eq.                                       C
6     !                                                                      C
7     !                                                                      C
8     !  Author: Jeff Dietiker                              Date: 01-MAY-09  C
9     !  Reviewer:                                          Date:            C
10     !                                                                      C
11     !  Revision Number:                                                    C
12     !                                                                      C
13     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
14     
15           SUBROUTINE CG_SOURCE_V_S(A_M, B_M, M)
16     
17     !-----------------------------------------------
18     ! Modules
19     !-----------------------------------------------
20           USE param, only: dimension_3, dimension_m
21           USE param1, only: zero, undefined
22     
23     ! number of solids phases
24           USE physprop, only: smax, mmax
25     ! virtual (added) mass coefficient Cv
26           USE physprop, only: cv
27     
28           USE bc
29           USE calc_gr_boundary
30           USE compar
31           USE cutcell
32           USE fldvar
33           USE fun_avg
34           USE geometry
35           USE indices
36           USE quadric
37           USE run
38           USE toleranc, only: dil_ep_s
39           USE visc_s, only: mu_s
40     
41           IMPLICIT NONE
42     !-----------------------------------------------
43     ! Dummy Arguments
44     !-----------------------------------------------
45     ! Solids phase index
46           INTEGER, INTENT(IN) :: M
47     ! Septadiagonal matrix A_m
48           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
49     ! Vector b_m
50           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
51     !-----------------------------------------------
52     ! Local variables
53     !-----------------------------------------------
54     ! Indices
55           INTEGER :: I, J, K, IJK,IMJK, IJMK, IJKM, IJKN
56     ! Solids phase index
57           INTEGER :: L
58     ! average volume fraction
59           DOUBLE PRECISION EPSA, EPStmp
60     ! virtual (added) mass
61           DOUBLE PRECISION :: F_vir, ROP_MA
62           DOUBLE PRECISION :: Vgn, Vgs, Uge, Ugw, Ugc, Vge, &
63                               Vgw, Wgt, Wgb, Wgc,Vgt, Vgb
64     !
65           DOUBLE PRECISION :: F_2
66     ! for cartesian grid
67           INTEGER :: IM,JM,IP,JP,KM,KP
68           INTEGER :: IPJK,IJPK,IJKP,IJKC,IJKE,IJKNE,IJKW,IJKWN,IMJPK,IJPKM
69           INTEGER :: IJKT,IJKTN,IJKB,IJKBN
70           DOUBLE PRECISION :: Vn,Vs,Ve,Vw, Vt,Vb
71           DOUBLE PRECISION :: B_NOC
72           DOUBLE PRECISION :: MU_S_E,MU_S_W,MU_S_N,MU_S_S,MU_S_T,MU_S_B,MU_S_CUT
73           DOUBLE PRECISION :: VW_s
74           INTEGER :: BCV
75           INTEGER :: BCT
76     !-----------------------------------------------
77     
78           IF(CG_SAFE_MODE(4)==1) RETURN
79     
80           IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
81              (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
82     
83              IF (MOMENTUM_Y_EQ(M)) THEN
84     
85     !!$omp  parallel do private(I, J, K, IJK, IJKN, &
86     !!$omp&                     EPSA,EPStmp,ROP_MA,F_vir) &
87     !!$omp&  schedule(static)
88                 DO IJK = ijkstart3, ijkend3
89     
90     ! Skip walls where some values are undefined.
91                    IF(WALL_AT(IJK)) cycle
92     
93                    I = I_OF(IJK)
94                    J = J_OF(IJK)
95                    K = K_OF(IJK)
96                    IMJK = IM_OF(IJK)
97                    IJMK = JM_OF(IJK)
98                    IJKM = KM_OF(IJK)
99                    IJKN = NORTH_OF(IJK)
100     
101                    IF (KT_TYPE_ENUM == GHD_2007) THEN   ! with ghd theory, m = mmax
102                       EPStmp = ZERO
103                       DO L = 1, SMAX
104                          EPStmp = EPStmp + AVG_Y(EP_S(IJK,L),EP_S(IJKN,L),J)
105                       ENDDO
106                       EPSA = EPStmp
107                    ELSE
108                       EPSA = AVG_Y(EP_S(IJK,M),EP_S(IJKN,M),J)
109                    ENDIF
110     
111     
112     ! Impermeable internal surface
113                    IF (IP_AT_N(IJK)) THEN   ! do nothing
114     
115     ! Semi-permeable internal surface
116                    ELSEIF (SIP_AT_N(IJK)) THEN   ! do nothing
117     
118     ! dilute flow
119                    ELSEIF (EPSA <= DIL_EP_S) THEN   ! do nothing
120     
121     
122                    ELSEIF(INTERIOR_CELL_AT(IJK)) THEN
123                       BCV = BC_V_ID(IJK)
124                       IF(BCV > 0 ) THEN
125                          BCT = BC_TYPE_ENUM(BCV)
126                       ELSE
127                          BCT = NONE
128                       ENDIF
129     
130                       SELECT CASE (BCT)
131                          CASE (CG_NSW)
132                             NOC_VS = .TRUE.
133                             VW_s = ZERO
134                             MU_S_CUT = (VOL(IJK)*MU_S(IJK,M) + VOL(IJKN)*MU_S(IJKN,M))/(VOL(IJK) + VOL(IJKN))
135                             A_M(IJK,0,M) = A_M(IJK,0,M)  - MU_S_CUT * Area_V_CUT(IJK)/DELH_V(IJK)
136                          CASE (CG_FSW)
137                             NOC_VS = .FALSE.
138                             VW_s = ZERO
139                          CASE(CG_PSW)
140                             IF(BC_JJ_PS(BCV)==1) THEN   ! Johnson-Jackson partial slip bc
141                                NOC_VS = .FALSE.
142                                VW_s = BC_VW_S(BCV,M)
143                                CALL CG_CALC_GRBDRY(IJK, 'V_MOMENTUM', M, BCV, F_2)
144                                A_M(IJK,0,M) = A_M(IJK,0,M) - Area_V_CUT(IJK)*F_2
145                                B_M(IJK,M) = B_M(IJK,M) - Area_V_CUT(IJK)*F_2*VW_s
146                             ELSEIF(BC_HW_S(BCV,M)==UNDEFINED) THEN   ! same as NSW
147                                NOC_VS = .TRUE.
148                                VW_s = BC_VW_S(BCV,M)
149                                MU_S_CUT = (VOL(IJK)*MU_S(IJK,M) + VOL(IJKN)*MU_S(IJKN,M))/(VOL(IJK) + VOL(IJKN))
150                                A_M(IJK,0,M) = A_M(IJK,0,M)  - MU_S_CUT * Area_V_CUT(IJK)/DELH_V(IJK)
151                                B_M(IJK,M) = B_M(IJK,M) - MU_S_CUT * VW_s * Area_V_CUT(IJK)/DELH_V(IJK)
152                             ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN   ! same as FSW
153                                NOC_VS = .FALSE.
154                                VW_s = ZERO
155                             ELSE                              ! partial slip
156                                NOC_VS = .FALSE.
157                                VW_s = BC_VW_S(BCV,M)
158                                MU_S_CUT = (VOL(IJK)*MU_S(IJK,M) + VOL(IJKN)*MU_S(IJKN,M))/(VOL(IJK) + VOL(IJKN))
159                                A_M(IJK,0,M) = A_M(IJK,0,M) - MU_S_CUT * Area_V_CUT(IJK)*(BC_HW_S(BCV,M))
160                                B_M(IJK,M) = B_M(IJK,M) - MU_S_CUT * VW_s * Area_V_CUT(IJK)*(BC_HW_S(BCV,M))
161                             ENDIF
162                          CASE (NONE, CG_MI)
163                             NOC_VS = .FALSE.
164     
165                       END SELECT
166     
167     
168                       IF(NOC_VS) THEN
169                          IMJK = IM_OF(IJK)
170                          IJMK = JM_OF(IJK)
171                          IJKM = KM_OF(IJK)
172                          IPJK = IP_OF(IJK)
173                          IJPK = JP_OF(IJK)
174                          IJKP = KP_OF(IJK)
175     
176                          Vn = Theta_Vn_bar(IJK)  * V_S(IJK,M)  + Theta_Vn(IJK)  * V_S(IJPK,M)
177                          Vs = Theta_Vn_bar(IJMK) * V_S(IJMK,M) + Theta_Vn(IJMK) * V_S(IJK,M)
178                          Ve = Theta_Ve_bar(IJK)  * V_S(IJK,M)  + Theta_Ve(IJK)  * V_S(IPJK,M)
179                          Vw = Theta_Ve_bar(IMJK) * V_S(IMJK,M) + Theta_Ve(IMJK) * V_S(IJK,M)
180     
181                          IJKN = NORTH_OF(IJK)
182                          IF (WALL_AT(IJK)) THEN
183                             IJKC = IJKN
184                          ELSE
185                             IJKC = IJK
186                          ENDIF
187     
188                          JP = JP1(J)
189                          IJKE = EAST_OF(IJK)
190                          IJKNE = EAST_OF(IJKN)
191                          IM = IM1(I)
192                          IJKW = WEST_OF(IJK)
193                          IJKWN = NORTH_OF(IJKW)
194                          IMJPK = JP_OF(IMJK)
195                          KM = KM1(K)
196                          IJKT = TOP_OF(IJK)
197                          IJKTN = NORTH_OF(IJKT)
198                          IJKB = BOTTOM_OF(IJK)
199                          IJKBN = NORTH_OF(IJKB)
200     
201                          MU_S_E = AVG_Y_H(AVG_X_H(MU_S(IJKC,M),MU_S(IJKE,M),I),&
202                                   AVG_X_H(MU_S(IJKN,M),MU_S(IJKNE,M),I),J)
203                          MU_S_W = AVG_Y_H(AVG_X_H(MU_S(IJKW,M),MU_S(IJKC,M),IM),&
204                                   AVG_X_H(MU_S(IJKWN,M),MU_S(IJKN,M),IM),J)
205                          MU_S_N = MU_S(IJKN,M)
206                          MU_S_S = MU_S(IJKC,M)
207     
208                          B_NOC =   MU_S_N * Axz_V(IJK)  * (Vn-VW_s) * NOC_V_N(IJK)   *2.0d0&
209                                  - MU_S_S * Axz_V(IJMK) * (Vs-VW_s) * NOC_V_N(IJMK)  *2.0d0&
210                                  + MU_S_E * Ayz_V(IJK)  * (Ve-VW_s) * NOC_V_E(IJK)   &
211                                  - MU_S_W * Ayz_V(IMJK) * (Vw-VW_s) * NOC_V_E(IMJK)
212     
213                          IF(DO_K) THEN
214     
215                             Vt = Theta_Vt_bar(IJK)  * V_S(IJK,M)  + Theta_Vt(IJK)  * V_S(IJKP,M)
216                             Vb = Theta_Vt_bar(IJKM) * V_S(IJKM,M) + Theta_Vt(IJKM) * V_S(IJK,M)
217     
218                             MU_S_T = AVG_Y_H(AVG_Z_H(MU_S(IJKC,M),MU_S(IJKT,M),K),&
219                                      AVG_Z_H(MU_S(IJKN,M),MU_S(IJKTN,M),K),J)
220                             MU_S_B = AVG_Y_H(AVG_Z_H(MU_S(IJKB,M),MU_S(IJKC,M),KM),&
221                                      AVG_Z_H(MU_S(IJKBN,M),MU_S(IJKN,M),KM),J)
222     
223                             B_NOC = B_NOC + MU_S_T * Axy_V(IJK)  * (Vt-VW_s) * NOC_V_T(IJK)   &
224                                           - MU_S_B * Axy_V(IJKM) * (Vb-VW_s) * NOC_V_T(IJKM)
225     
226                          ENDIF
227                          B_M(IJK,M) = B_M(IJK,M) + B_NOC
228                       ENDIF   ! end if noc_vs
229     
230     
231                       IF(CUT_V_TREATMENT_AT(IJK)) THEN
232     ! VIRTUAL MASS SECTION (explicit terms)
233     ! adding transient term dVg/dt to virtual mass term
234                          F_vir = ZERO
235                          IF(Added_Mass.AND. M==M_AM ) THEN
236                             F_vir = ( (V_g(IJK) - V_gO(IJK)) )*ODT*VOL_V(IJK)
237                             IM = I - 1
238                             JM = J - 1
239                             KM = K - 1
240                             IP = I + 1
241                             JP = J + 1
242                             KP = K + 1
243                             IMJK = FUNIJK(IM,J,K)
244                             IJMK = FUNIJK(I,JM,K)
245                             IPJK = FUNIJK(IP,J,K)
246                             IJPK = FUNIJK(I,JP,K)
247                             IJKP = FUNIJK(I,J,KP)
248                             IJKM = FUNIJK(I,J,KM)
249                             IMJPK = IM_OF(IJPK)
250                             IJKN = NORTH_OF(IJK)
251     
252     ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
253                             Vge = Theta_Ve_bar(IJK)  * V_g(IJK)  + Theta_Ve(IJK)  * V_g(IPJK)
254                             Vgw = Theta_Ve_bar(IMJK) * V_g(IMJK) + Theta_Ve(IMJK) * V_g(IJK)
255                             Uge =  Theta_V_ne(IJK)  * U_g(IJK)  + Theta_V_se(IJK)  * U_g(IJPK)
256                             Ugw =  Theta_V_ne(IMJK) * U_g(IMJK) + Theta_V_se(IMJK) * U_g(IMJPK)
257                             Ugc = (DELX_ve(IJK) * Ugw + DELX_vw(IJK) * Uge) / (DELX_ve(IJK) + DELX_vw(IJK))
258                             Vgn = Theta_Vn_bar(IJK)  * V_g(IJK)  + Theta_Vn(IJK)  * V_g(IJPK)
259                             Vgs = Theta_Vn_bar(IJMK) * V_g(IJMK) + Theta_Vn(IJMK) * V_g(IJK)
260                             IF(DO_K) THEN
261                                IJPKM = KM_OF(IJPK)
262                                Vgt = Theta_Vt_bar(IJK)  * V_g(IJK)  + Theta_Vt(IJK)  * V_g(IJKP)
263                                Vgb = Theta_Vt_bar(IJKM) * V_g(IJKM) + Theta_Vt(IJKM) * V_g(IJK)
264                                Wgt = Theta_V_nt(IJK)  * W_g(IJK)  + Theta_V_st(IJK)  * W_g(IJPK)
265                                Wgb = Theta_V_nt(IJKM) * W_g(IJKM) + Theta_V_st(IJKM) * W_g(IJPKM)
266                                Wgc = (DELZ_vt(IJK) * Wgb + DELZ_vb(IJK) * Wgt) / (DELZ_vt(IJK) + DELZ_vb(IJK))
267                                F_vir = F_vir +  Wgc* (Vgt - Vgb)*AXY(IJK)
268                             ENDIF
269     
270     ! adding convective terms (U dV/dx + V dV/dy) to virtual mass; W dV/dz added above.
271                             F_vir = F_vir + V_g(IJK)*(Vgn - Vgs)*AXZ(IJK) + &
272                                     Ugc*(Vge - Vgw)*AYZ(IJK)
273                             ROP_MA =  (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M) + &
274                                        VOL(IJKN)*ROP_g(IJKN)*EP_s(IJKN,M))/&
275                                       (VOL(IJK) + VOL(IJKN))
276                             F_vir = F_vir * Cv * ROP_MA
277                             B_M(IJK,M) = B_M(IJK,M) - F_vir ! adding explicit-part of virtual mass force.
278                          ENDIF   ! end if added_mass .and. m==m_am)
279                       ENDIF   ! end if cut_v_treatment_at
280     
281     
282                    ENDIF ! end if sip or ip or dilute flow branch
283                 ENDDO   ! end do ijk
284              ENDIF   ! end if momentum_y_eq
285           ENDIF   ! end if for GHD Theory
286     
287           RETURN
288     
289         CONTAINS
290     
291           INCLUDE 'functions.inc'
292     
293           END SUBROUTINE CG_SOURCE_V_S
294     
295     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
296     !                                                                      C
297     !  Subroutine: CG_SOURCE_V_s_BC                                        C
298     !  Purpose: Determine contribution of cut-cell to source terms         C
299     !           for V_s momentum eq.                                       C
300     !                                                                      C
301     !                                                                      C
302     !  Author: Jeff Dietiker                              Date: 01-MAY-09  C
303     !  Reviewer:                                          Date:            C
304     !                                                                      C
305     !                                                                      C
306     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
307     
308           SUBROUTINE CG_SOURCE_V_S_BC(A_M, B_M, M)
309     
310     !-----------------------------------------------
311     ! Modules
312     !-----------------------------------------------
313           USE param
314           USE param1, only: zero, one, undefined
315           USE fldvar
316           USE bc
317           USE compar
318           USE geometry
319           USE indices
320           USE cutcell
321           USE quadric
322           USE fun_avg
323     
324           IMPLICIT NONE
325     !-----------------------------------------------
326     ! Dummy arguments
327     !-----------------------------------------------
328     ! Solids phase
329           INTEGER, INTENT(IN) :: M
330     ! Septadiagonal matrix A_m
331           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
332     ! Vector b_m
333           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
334     !-----------------------------------------------
335     ! Local variables
336     !-----------------------------------------------
337     ! Indices
338           INTEGER :: IJK, IJKS
339     ! for cartesian grid:
340           INTEGER :: BCV
341           INTEGER :: BCT
342     !-----------------------------------------------
343     
344           IF(CG_SAFE_MODE(4)==1) RETURN
345     
346           DO IJK = ijkstart3, ijkend3
347     
348              BCV = BC_V_ID(IJK)
349              IF(BCV > 0 ) THEN
350                 BCT = BC_TYPE_ENUM(BCV)
351              ELSE
352                 BCT = NONE
353              ENDIF
354     
355              SELECT CASE (BCT)
356                 CASE (CG_NSW)
357                    IF(WALL_V_AT(IJK)) THEN
358                       A_M(IJK,east,M) = ZERO
359                       A_M(IJK,west,M) = ZERO
360                       A_M(IJK,north,M) = ZERO
361                       A_M(IJK,south,M) = ZERO
362                       A_M(IJK,top,M) = ZERO
363                       A_M(IJK,bottom,M) = ZERO
364                       A_M(IJK,0,M) = -ONE
365                       B_M(IJK,M) = ZERO
366                    ENDIF
367     
368     
369                 CASE (CG_FSW)
370                    IF(WALL_V_AT(IJK)) THEN
371                       A_M(IJK,east,M) = ZERO
372                       A_M(IJK,west,M) = ZERO
373                       A_M(IJK,north,M) = ZERO
374                       A_M(IJK,south,M) = ZERO
375                       A_M(IJK,top,M) = ZERO
376                       A_M(IJK,bottom,M) = ZERO
377                       A_M(IJK,0,M) = -ONE
378     !                  B_M(IJK,M) = - V_s(V_MASTER_OF(IJK),M)    ! Velocity of master node
379                       B_M(IJK,M) = ZERO
380                       IF(DABS(NORMAL_V(IJK,2))/=ONE) THEN
381                          IF (V_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
382                             A_M(IJK,east,M) = ONE
383                          ELSEIF (V_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
384                             A_M(IJK,west,M) = ONE
385                          ELSEIF (V_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
386                             A_M(IJK,north,M) = ONE
387                          ELSEIF (V_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
388                             A_M(IJK,south,M) = ONE
389                          ELSEIF (V_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
390                             A_M(IJK,top,M) = ONE
391                          ELSEIF (V_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
392                             A_M(IJK,bottom,M) = ONE
393                          ENDIF
394                       ENDIF
395                    ENDIF
396     
397     
398                 CASE (CG_PSW)
399                    IF(WALL_V_AT(IJK)) THEN
400                       A_M(IJK,east,M) = ZERO
401                       A_M(IJK,west,M) = ZERO
402                       A_M(IJK,north,M) = ZERO
403                       A_M(IJK,south,M) = ZERO
404                       A_M(IJK,top,M) = ZERO
405                       A_M(IJK,bottom,M) = ZERO
406                       A_M(IJK,0,M) = -ONE
407                       IF(BC_HW_S(BCV,M)==UNDEFINED) THEN   ! same as NSW
408                          B_M(IJK,M) = -BC_VW_S(BCV,M)
409                       ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN   ! same as FSW
410                          B_M(IJK,M) = ZERO
411                          IF(DABS(NORMAL_V(IJK,2))/=ONE) THEN
412                             IF (V_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
413                                A_M(IJK,east,M) = ONE
414                             ELSEIF (V_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
415                                A_M(IJK,west,M) = ONE
416                             ELSEIF (V_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
417                                A_M(IJK,north,M) = ONE
418                             ELSEIF (V_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
419                                A_M(IJK,south,M) = ONE
420                             ELSEIF (V_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
421                                A_M(IJK,top,M) = ONE
422                             ELSEIF (V_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
423                                A_M(IJK,bottom,M) = ONE
424                             ENDIF
425                          ENDIF
426                       ELSE                              ! partial slip  (WARNING:currently same as FSW)
427                          B_M(IJK,M) = ZERO
428                          IF(DABS(NORMAL_V(IJK,2))/=ONE) THEN
429                             IF (V_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
430                                A_M(IJK,east,M) = ONE
431                             ELSEIF (V_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
432                                A_M(IJK,west,M) = ONE
433                             ELSEIF (V_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
434                                A_M(IJK,north,M) = ONE
435                             ELSEIF (V_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
436                                A_M(IJK,south,M) = ONE
437                             ELSEIF (V_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
438                                A_M(IJK,top,M) = ONE
439                             ELSEIF (V_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
440                                A_M(IJK,bottom,M) = ONE
441                             ENDIF
442                          ENDIF
443                       ENDIF
444                    ENDIF
445     
446     
447                 CASE (CG_MI)
448                    A_M(IJK,east,M) = ZERO
449                    A_M(IJK,west,M) = ZERO
450                    A_M(IJK,north,M) = ZERO
451                    A_M(IJK,south,M) = ZERO
452                    A_M(IJK,top,M) = ZERO
453                    A_M(IJK,bottom,M) = ZERO
454                    A_M(IJK,0,M) = -ONE
455                    IF(BC_V_s(BCV,M)/=UNDEFINED) THEN
456                       B_M(IJK,M) = - BC_V_s(BCV,M)
457                    ELSE
458                       B_M(IJK,M) = - BC_VELMAG_s(BCV,M)*NORMAL_V(IJK,2)
459                    ENDIF
460                    IJKS = SOUTH_OF(IJK)
461                    IF(FLUID_AT(IJKS)) THEN
462                       A_M(IJKS,east,M) = ZERO
463                       A_M(IJKS,west,M) = ZERO
464                       A_M(IJKS,north,M) = ZERO
465                       A_M(IJKS,south,M) = ZERO
466                       A_M(IJKS,top,M) = ZERO
467                       A_M(IJKS,bottom,M) = ZERO
468                       A_M(IJKS,0,M) = -ONE
469                       IF(BC_V_s(BCV,M)/=UNDEFINED) THEN
470                          B_M(IJKS,M) = - BC_V_s(BCV,M)
471                       ELSE
472                          B_M(IJKS,M) = - BC_VELMAG_s(BCV,M)*NORMAL_V(IJK,2)
473                       ENDIF
474                    ENDIF
475     
476     
477                 CASE (CG_PO)
478                    A_M(IJK,east,M) = ZERO
479                    A_M(IJK,west,M) = ZERO
480                    A_M(IJK,north,M) = ZERO
481                    A_M(IJK,south,M) = ZERO
482                    A_M(IJK,top,M) = ZERO
483                    A_M(IJK,bottom,M) = ZERO
484                    A_M(IJK,0,M) = -ONE
485                    B_M(IJK,M) = ZERO
486     
487                    IJKS = SOUTH_OF(IJK)
488                    IF(FLUID_AT(IJKS)) THEN
489     
490                       A_M(IJK,south,M) = ONE
491                       A_M(IJK,0,M) = -ONE
492                       B_M(IJK,M) = ZERO
493                    ENDIF
494     
495              END SELECT   ! end select bc_type
496     
497     
498     ! this is straight repeat of section of above code...?
499              BCV = BC_ID(IJK)
500              IF(BCV > 0 ) THEN
501                 BCT = BC_TYPE_ENUM(BCV)
502              ELSE
503                 BCT = NONE
504              ENDIF
505     
506     
507              SELECT CASE (BCT)
508                 CASE (CG_MI)
509                    A_M(IJK,east,M) = ZERO
510                    A_M(IJK,west,M) = ZERO
511                    A_M(IJK,north,M) = ZERO
512                    A_M(IJK,south,M) = ZERO
513                    A_M(IJK,top,M) = ZERO
514                    A_M(IJK,bottom,M) = ZERO
515                    A_M(IJK,0,M) = -ONE
516                    IF(BC_V_s(BCV,M)/=UNDEFINED) THEN
517                       B_M(IJK,M) = - BC_V_s(BCV,M)
518                    ELSE
519                       B_M(IJK,M) = - BC_VELMAG_s(BCV,M)*NORMAL_S(IJK,2)
520                    ENDIF
521                    IJKS = SOUTH_OF(IJK)
522                    IF(FLUID_AT(IJKS)) THEN
523                       A_M(IJKS,east,M) = ZERO
524                       A_M(IJKS,west,M) = ZERO
525                       A_M(IJKS,north,M) = ZERO
526                       A_M(IJKS,south,M) = ZERO
527                       A_M(IJKS,top,M) = ZERO
528                       A_M(IJKS,bottom,M) = ZERO
529                       A_M(IJKS,0,M) = -ONE
530                       IF(BC_V_s(BCV,M)/=UNDEFINED) THEN
531                          B_M(IJKS,M) = - BC_V_s(BCV,M)
532                       ELSE
533                          B_M(IJKS,M) = - BC_VELMAG_s(BCV,M)*NORMAL_S(IJK,2)
534                       ENDIF
535                    ENDIF
536     
537     
538                 CASE (CG_PO)
539                    A_M(IJK,east,M) = ZERO
540                    A_M(IJK,west,M) = ZERO
541                    A_M(IJK,north,M) = ZERO
542                    A_M(IJK,south,M) = ZERO
543                    A_M(IJK,top,M) = ZERO
544                    A_M(IJK,bottom,M) = ZERO
545                    A_M(IJK,0,M) = -ONE
546                    B_M(IJK,M) = ZERO
547                    IJKS = SOUTH_OF(IJK)
548                    IF(FLUID_AT(IJKS)) THEN
549                       A_M(IJK,south,M) = ONE
550                       A_M(IJK,0,M) = -ONE
551                    ENDIF
552     
553              END SELECT
554     
555           ENDDO
556     
557           RETURN
558     
559         CONTAINS
560     
561           INCLUDE 'functions.inc'
562     
563           END SUBROUTINE CG_SOURCE_V_S_BC
564