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