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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CG_SOURCE_W_s                                           C
4     !  Purpose: Determine contribution of cut-cell to source terms         C
5     !           for W_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_W_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     !-----------------------------------------------
44     ! Dummy Arguments
45     !-----------------------------------------------
46     ! Solids phase index
47           INTEGER, INTENT(IN) :: M
48     ! Septadiagonal matrix A_m
49           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
50     ! Vector b_m
51           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
52     !-----------------------------------------------
53     ! Local variables
54     !-----------------------------------------------
55     ! Indices
56           INTEGER :: I, J, K, IJK, IJKT, IMJK, IJMK, IJKM, IJKP, IMJKP
57           INTEGER :: IJKE, IJKW, IJKTE, IM, IPJK
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 :: Uge, Ugw, Wge, Wgw, Wgn, &
65                               Wgs, Wgt, Wgb, Ugc, Vgc, Vgn, Vgs
66     !
67           DOUBLE PRECISION :: F_2
68     ! for cartesian grid:
69           INTEGER :: JM,IP,JP,IJPK,IJKC,IJKN,IJKNE,IJKS,IJKSE,IPJMK,KM,KP,IJMKP
70           INTEGER :: IJKTN,IJKWT,IJKST
71           DOUBLE PRECISION :: We,Ww,Wn,Ws,Wt,Wb
72           DOUBLE PRECISION :: B_NOC
73           DOUBLE PRECISION :: MU_S_E,MU_S_W,MU_S_N,MU_S_S,MU_S_T,MU_S_B,MU_S_CUT
74           DOUBLE PRECISION :: WW_s
75           INTEGER :: BCV
76           INTEGER :: BCT
77     !-----------------------------------------------
78     
79           IF(CG_SAFE_MODE(5)==1) RETURN
80     
81     
82           IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
83              (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
84     
85              IF (MOMENTUM_Z_EQ(M)) THEN
86     
87     !!$omp  parallel do &
88     !!$omp& private(IJK, I, J, K, IJKT, EPSA, EPStmp, &
89     !!$omp&         IMJK,IJKP,IMJKP, &
90     !!$omp&         IJKE,IJKW,IJKTE,IJKTW,IM,IPJK, &
91     !!$omp&         ROP_MA,F_vir)
92                 DO IJK = ijkstart3, ijkend3
93     
94     ! Skip walls where some values are undefined.
95                    IF(WALL_AT(IJK)) cycle
96     
97                    I = I_OF(IJK)
98                    J = J_OF(IJK)
99                    K = K_OF(IJK)
100                    IMJK = IM_OF(IJK)
101                    IJMK = JM_OF(IJK)
102                    IJKM = KM_OF(IJK)
103                    IJKT = TOP_OF(IJK)
104     
105                    IF (KT_TYPE_ENUM == GHD_2007) THEN   ! with ghd theory, m = mmax
106                       EPStmp = ZERO
107                       DO L = 1, SMAX
108                         EPStmp = EPStmp + AVG_Z(EP_S(IJK,L),EP_S(IJKT,L),K)
109                       ENDDO
110                       EPSA = EPStmp
111                     ELSE
112                       EPSA = AVG_Z(EP_S(IJK,M),EP_S(IJKT,M),K)
113                     ENDIF
114     
115     ! Impermeable internal surface
116                     IF (IP_AT_T(IJK)) THEN   ! do nothing
117     
118     ! Semi-permeable internal surface
119                     ELSEIF (SIP_AT_T(IJK)) THEN   ! do nothing
120     
121     ! dilute flow
122                     ELSEIF (EPSA <= DIL_EP_S) THEN   ! do nothing
123     
124     
125                     ELSEIF(INTERIOR_CELL_AT(IJK)) THEN
126                       BCV = BC_W_ID(IJK)
127                       IF(BCV > 0 ) THEN
128                          BCT = BC_TYPE_ENUM(BCV)
129                       ELSE
130                          BCT = NONE
131                       ENDIF
132     
133                       SELECT CASE (BCT)
134                          CASE (CG_NSW)
135                             NOC_WS = .TRUE.
136                             WW_s = ZERO
137                             MU_S_CUT = (VOL(IJK)*MU_S(IJK,M) + VOL(IJKT)*MU_S(IJKT,M))/(VOL(IJK) + VOL(IJKT))
138                             A_M(IJK,0,M) = A_M(IJK,0,M) - MU_S_CUT * Area_W_CUT(IJK)/DELH_W(IJK)
139                          CASE (CG_FSW)
140                             NOC_WS = .FALSE.
141                             WW_s = ZERO
142                          CASE(CG_PSW)
143                             IF(BC_JJ_PS(BCV)==1) THEN   ! Johnson-Jackson partial slip bc
144                                NOC_WS = .FALSE.
145                                WW_s = BC_WW_S(BCV,M)
146                                CALL CG_CALC_GRBDRY(IJK, 'W_MOMENTUM', M, BCV, F_2)
147                                A_M(IJK,0,M) = A_M(IJK,0,M) - Area_W_CUT(IJK)*F_2
148                                B_M(IJK,M) = B_M(IJK,M) - Area_W_CUT(IJK)*F_2*WW_s
149                             ELSEIF(BC_HW_S(BCV,M)==UNDEFINED) THEN   ! same as NSW
150                                NOC_WS = .TRUE.
151                                WW_s = BC_WW_S(BCV,M)
152                                MU_S_CUT = (VOL(IJK)*MU_S(IJK,M) + VOL(IJKT)*MU_S(IJKT,M))/(VOL(IJK) + VOL(IJKT))
153                                A_M(IJK,0,M) = A_M(IJK,0,M) - MU_S_CUT * Area_W_CUT(IJK)/DELH_W(IJK)
154                                B_M(IJK,M) = B_M(IJK,M) - MU_S_CUT * WW_s * Area_W_CUT(IJK)/DELH_W(IJK)
155                             ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN   ! same as FSW
156                                NOC_WS = .FALSE.
157                                WW_s = ZERO
158                             ELSE                              ! partial slip
159                                NOC_WS = .FALSE.
160                                WW_s = BC_WW_S(BCV,M)
161                                MU_S_CUT = (VOL(IJK)*MU_S(IJK,M) + VOL(IJKT)*MU_S(IJKT,M))/(VOL(IJK) + VOL(IJKT))
162                                A_M(IJK,0,M) = A_M(IJK,0,M) - MU_S_CUT * Area_W_CUT(IJK)*(BC_HW_S(BCV,M))
163                                B_M(IJK,M) = B_M(IJK,M) - MU_S_CUT * WW_s * Area_W_CUT(IJK)*(BC_HW_S(BCV,M))
164                             ENDIF
165                          CASE (NONE,CG_MI)
166                             NOC_WS = .FALSE.
167     
168                       END SELECT
169     
170     
171                       IF(NOC_WS) THEN
172                          IMJK = IM_OF(IJK)
173                          IJMK = JM_OF(IJK)
174                          IJKM = KM_OF(IJK)
175                          IPJK = IP_OF(IJK)
176                          IJPK = JP_OF(IJK)
177                          IJKP = KP_OF(IJK)
178     
179                          We = Theta_We_bar(IJK)  * W_S(IJK,M)  + Theta_We(IJK)  * W_S(IPJK,M)
180                          Ww = Theta_We_bar(IMJK) * W_S(IMJK,M) + Theta_We(IMJK) * W_S(IJK,M)
181                          Wn = Theta_Wn_bar(IJK)  * W_S(IJK,M)  + Theta_Wn(IJK)  * W_S(IJPK,M)
182                          Ws = Theta_Wn_bar(IJMK) * W_S(IJMK,M) + Theta_Wn(IJMK) * W_S(IJK,M)
183                          Wt = Theta_Wt_bar(IJK)  * W_S(IJK,M)  + Theta_Wt(IJK)  * W_S(IJKP,M)
184                          Wb = Theta_Wt_bar(IJKM) * W_S(IJKM,M) + Theta_Wt(IJKM) * W_S(IJK,M)
185     
186                          IF (WALL_AT(IJK)) THEN
187                             IJKC = IJKT
188                          ELSE
189                             IJKC = IJK
190                          ENDIF
191     
192                          IJKE = EAST_OF(IJK)
193                          ijkt = top_of(ijk)
194                          IP = IP1(I)
195                          IM = IM1(I)
196                          IJKN = NORTH_OF(IJK)
197                          IJKNE = EAST_OF(IJKN)
198                          JM = JM1(J)
199                          IPJMK = IP_OF(IJMK)
200                          IJKS = SOUTH_OF(IJK)
201                          IJKSE = EAST_OF(IJKS)
202                          KP = KP1(K)
203                          IJKT = TOP_OF(IJK)
204                          IJKE = EAST_OF(IJK)
205                          IJKP = KP_OF(IJK)
206                          IJKTN = NORTH_OF(IJKT)
207                          IJKTE = EAST_OF(IJKT)
208                          IJKW = WEST_OF(IJK)
209                          IJKWT = TOP_OF(IJKW)
210                          IJKS = SOUTH_OF(IJK)
211                          IJKST = TOP_OF(IJKS)
212     
213                          MU_S_E = AVG_Z_H(AVG_X_H(MU_S(IJKC,M),MU_S(IJKE,M),I),&
214                                   AVG_X_H(MU_S(IJKT,M),MU_S(IJKTE,M),I),K)
215                          MU_S_W = AVG_Z_H(AVG_X_H(MU_S(IJKW,M),MU_S(IJKC,M),IM),&
216                                   AVG_X_H(MU_S(IJKWT,M),MU_S(IJKT,M),IM),K)
217                          MU_S_N = AVG_Z_H(AVG_Y_H(MU_S(IJKC,M),MU_S(IJKN,M),J),&
218                                   AVG_Y_H(MU_S(IJKT,M),MU_S(IJKTN,M),J),K)
219                          MU_S_S = AVG_Z_H(AVG_Y_H(MU_S(IJKS,M),MU_S(IJKC,M),JM),&
220                                   AVG_Y_H(MU_S(IJKST,M),MU_S(IJKT,M),JM),K)
221                          MU_S_T = MU_S(IJKT,M)
222                          MU_S_B = MU_S(IJKC,M)
223     
224                          B_NOC =     MU_S_E * Ayz_W(IJK)  * (We-WW_s) * NOC_W_E(IJK)  &
225                                  -   MU_S_W * Ayz_W(IMJK) * (Ww-WW_s) * NOC_W_E(IMJK) &
226                                  +   MU_S_N * Axz_W(IJK)  * (Wn-WW_s) * NOC_W_N(IJK)  &
227                                  -   MU_S_S * Axz_W(IJMK) * (Ws-WW_s) * NOC_W_N(IJMK) &
228                                  +   MU_S_T * Axy_W(IJK)  * (Wt-WW_s) * NOC_W_T(IJK)  *2.0d0&
229                                  -   MU_S_B * Axy_W(IJKM) * (Wb-WW_s) * NOC_W_T(IJKM) *2.0D0
230     
231                          B_M(IJK,M) = B_M(IJK,M)   +  B_NOC
232                       ENDIF
233     
234     
235                       IF(CUT_W_TREATMENT_AT(IJK)) THEN
236     ! VIRTUAL MASS SECTION (explicit terms)
237     ! adding transient term  dWg/dt to virtual mass term
238                          F_vir = ZERO
239                          IF(Added_Mass.AND. M == M_AM ) THEN
240                             F_vir = ( (W_g(IJK) - W_gO(IJK)) )*ODT*VOL_W(IJK)
241                             I = I_OF(IJK)
242                             J = J_OF(IJK)
243                             K = K_OF(IJK)
244                             IM = I - 1
245                             JM = J - 1
246                             KM = K - 1
247                             IP = I + 1
248                             JP = J + 1
249                             KP = K + 1
250                             IMJK = FUNIJK(IM,J,K)
251                             IJMK = FUNIJK(I,JM,K)
252                             IPJK = FUNIJK(IP,J,K)
253                             IJPK = FUNIJK(I,JP,K)
254                             IJKP = FUNIJK(I,J,KP)
255                             IJKM = FUNIJK(I,J,KM)
256                             IMJKP = KP_OF(IMJK)
257                             IJMKP = KP_OF(IJMK)
258                             IJKE = EAST_OF(IJK)
259     
260     ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
261                             Wge = Theta_We_bar(IJK) * W_g(IJK) + Theta_We(IJK) * W_g(IPJK)
262                             Wgw = Theta_We_bar(IMJK) * W_g(IMJK) + Theta_We(IMJK) * W_g(IJK)
263                             Uge = Theta_W_te(IJK) * U_g(IJK) + Theta_W_be(IJK) * U_g(IJKP)
264                             Ugw = Theta_W_te(IMJK) * U_g(IMJK) + Theta_W_be(IMJK) * U_g(IMJKP)
265                             Ugc = (DELX_we(IJK) * Ugw + DELX_ww(IJK) * Uge) / (DELX_we(IJK) + DELX_ww(IJK))
266                             Wgn = Theta_Wn_bar(IJK) * W_g(IJK) + Theta_Wn(IJK) * W_g(IJPK)
267                             Wgs = Theta_Wn_bar(IJMK) * W_g(IJMK) + Theta_Wn(IJMK) * W_g(IJK)
268                             Vgn =  Theta_W_tn(IJK)  * V_g(IJK)  + Theta_W_bn(IJK)  * V_g(IJKP)
269                             Vgs =  Theta_W_tn(IJMK) * V_g(IJMK) + Theta_W_bn(IJMK) * V_g(IJMKP)
270                             Vgc = (DELY_wn(IJK) * Vgs + DELY_ws(IJK) * Vgn) / (DELY_wn(IJK) + DELY_ws(IJK))
271                             Wgt = Theta_Wt_bar(IJK)  * W_g(IJK)  + Theta_Wt(IJK)  * W_g(IJKP)
272                             Wgb = Theta_Wt_bar(IMJK) * W_g(IMJK) + Theta_Wt(IMJK) * W_g(IMJKP)
273     
274     ! adding convective terms (U dW/dx + V dW/dy + W dW/dz) to virtual mass
275                             F_vir = F_vir +  Ugc * (Wge - Wgw)*AYZ(IJK)    + &
276                                              Vgc * (Wgn - Wgs)*AXZ(IJK)  + &
277                                              W_g(IJK)*(Wgt - Wgb)*AXY(IJK)
278                             ROP_MA = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M) + &
279                                       VOL(IJKT)*ROP_g(IJKT)*EP_s(IJKT,M))/&
280                                      (VOL(IJK) + VOL(IJKT))
281                             F_vir = F_vir * Cv * ROP_MA
282                             B_M(IJK,M) = B_M(IJK,M) - F_vir ! explicit part of virtual mass force
283                          ENDIF   ! end if added_mass .and. m==m_am)
284                       ENDIF   ! end if cut_w_treatment_at
285     
286     
287                    ENDIF ! end if sip or ip or dilute flow branch
288                 ENDDO   ! end do ijk
289              ENDIF   ! end if momentum_z_eq
290           ENDIF   ! end if for GHD Theory
291     
292           RETURN
293     
294         CONTAINS
295     
296           INCLUDE 'functions.inc'
297     
298           END SUBROUTINE CG_SOURCE_W_S
299     
300     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
301     !                                                                      C
302     !  Subroutine name: CG_SOURCE_W_s_BC                                   C
303     !  Purpose: Determine contribution of cut-cell to source terms         C
304     !           for W_s momentum eq.                                       C
305     !                                                                      C
306     !                                                                      C
307     !  Author: Jeff Dietiker                              Date: 01-MAY-09  C
308     !  Reviewer:                                          Date:            C
309     !                                                                      C
310     !                                                                      C
311     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
312     
313           SUBROUTINE CG_SOURCE_W_S_BC(A_M, B_M, M)
314     
315     !-----------------------------------------------
316     ! Modules
317     !-----------------------------------------------
318           USE param
319           USE param1, only: zero, one, undefined
320     
321           USE fldvar
322           USE bc
323     
324           USE compar
325           USE geometry
326           USE indices
327     
328           USE cutcell
329           USE quadric
330           USE fun_avg
331     
332           IMPLICIT NONE
333     !-----------------------------------------------
334     ! Dummy arguments
335     !-----------------------------------------------
336     ! Solids phase
337           INTEGER, INTENT(IN) :: M
338     ! Septadiagonal matrix A_m
339           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
340     ! Vector b_m
341           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
342     !-----------------------------------------------
343     ! Local variables
344     !-----------------------------------------------
345     ! Indices
346           INTEGER :: IJK, IJKB
347     ! for cartesian grid:
348           INTEGER :: BCV
349           INTEGER :: BCT
350     !-----------------------------------------------
351     
352           IF(CG_SAFE_MODE(5)==1) RETURN
353     
354           DO IJK = ijkstart3, ijkend3
355     
356              BCV = BC_W_ID(IJK)
357              IF(BCV > 0 ) THEN
358                 BCT = BC_TYPE_ENUM(BCV)
359              ELSE
360                 BCT = NONE
361              ENDIF
362     
363              SELECT CASE (BCT)
364                 CASE (CG_NSW)
365                    IF(WALL_W_AT(IJK)) THEN
366                       A_M(IJK,east,M) = ZERO
367                       A_M(IJK,west,M) = ZERO
368                       A_M(IJK,north,M) = ZERO
369                       A_M(IJK,south,M) = ZERO
370                       A_M(IJK,top,M) = ZERO
371                       A_M(IJK,bottom,M) = ZERO
372                       A_M(IJK,0,M) = -ONE
373                       B_M(IJK,M) = ZERO
374                    ENDIF
375     
376     
377                 CASE (CG_FSW)
378                    IF(WALL_W_AT(IJK)) THEN
379                       A_M(IJK,east,M) = ZERO
380                       A_M(IJK,west,M) = ZERO
381                       A_M(IJK,north,M) = ZERO
382                       A_M(IJK,south,M) = ZERO
383                       A_M(IJK,top,M) = ZERO
384                       A_M(IJK,bottom,M) = ZERO
385                       A_M(IJK,0,M) = -ONE
386     !                  B_M(IJK,M) = - W_s(W_MASTER_OF(IJK),M)  ! Velocity of master node
387                       B_M(IJK,M) = ZERO
388                       IF(DABS(NORMAL_W(IJK,3))/=ONE) THEN
389                          IF (W_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
390                             A_M(IJK,east,M) = ONE
391                          ELSEIF (W_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
392                             A_M(IJK,west,M) = ONE
393                          ELSEIF (W_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
394                             A_M(IJK,north,M) = ONE
395                          ELSEIF (W_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
396                             A_M(IJK,south,M) = ONE
397                          ELSEIF (W_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
398                             A_M(IJK,top,M) = ONE
399                          ELSEIF (W_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
400                             A_M(IJK,bottom,M) = ONE
401                          ENDIF
402                       ENDIF
403                    ENDIF
404     
405     
406                 CASE (CG_PSW)
407                    IF(WALL_W_AT(IJK)) THEN
408                       A_M(IJK,east,M) = ZERO
409                       A_M(IJK,west,M) = ZERO
410                       A_M(IJK,north,M) = ZERO
411                       A_M(IJK,south,M) = ZERO
412                       A_M(IJK,top,M) = ZERO
413                       A_M(IJK,bottom,M) = ZERO
414                       A_M(IJK,0,M) = -ONE
415                       IF(BC_HW_S(BCV,M)==UNDEFINED) THEN   ! same as NSW
416                          B_M(IJK,M) = -BC_WW_S(BCV,M)
417                       ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN   ! same as FSW
418                          B_M(IJK,M) = ZERO
419                          IF(DABS(NORMAL_W(IJK,3))/=ONE) THEN
420                             IF (W_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
421                                A_M(IJK,east,M) = ONE
422                             ELSEIF (W_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
423                                A_M(IJK,west,M) = ONE
424                             ELSEIF (W_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
425                                A_M(IJK,north,M) = ONE
426                             ELSEIF (W_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
427                                A_M(IJK,south,M) = ONE
428                             ELSEIF (W_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
429                                A_M(IJK,top,M) = ONE
430                             ELSEIF (W_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
431                                A_M(IJK,bottom,M) = ONE
432                             ENDIF
433                          ENDIF
434                       ELSE                              ! partial slip  (WARNING:currently same as FSW)
435                          B_M(IJK,M) = ZERO
436                          IF(DABS(NORMAL_W(IJK,3))/=ONE) THEN
437                             IF (W_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
438                                A_M(IJK,east,M) = ONE
439                             ELSEIF (W_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
440                                A_M(IJK,west,M) = ONE
441                             ELSEIF (W_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
442                                A_M(IJK,north,M) = ONE
443                             ELSEIF (W_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
444                                A_M(IJK,south,M) = ONE
445                             ELSEIF (W_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
446                                A_M(IJK,top,M) = ONE
447                             ELSEIF (W_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
448                                A_M(IJK,bottom,M) = ONE
449                             ENDIF
450                          ENDIF
451                       ENDIF
452                    ENDIF
453     
454     
455                 CASE (CG_MI)
456                    A_M(IJK,east,M) = ZERO
457                    A_M(IJK,west,M) = ZERO
458                    A_M(IJK,north,M) = ZERO
459                    A_M(IJK,south,M) = ZERO
460                    A_M(IJK,top,M) = ZERO
461                    A_M(IJK,bottom,M) = ZERO
462                    A_M(IJK,0,M) = -ONE
463                    IF(BC_W_s(BCV,M)/=UNDEFINED) THEN
464                       B_M(IJK,M) = - BC_W_s(BCV,M)
465                    ELSE
466                       B_M(IJK,M) = - BC_VELMAG_s(BCV,M)*NORMAL_W(IJK,3)
467                    ENDIF
468                    IJKB = BOTTOM_OF(IJK)
469                    IF(FLUID_AT(IJKB)) THEN
470                       A_M(IJKB,east,M) = ZERO
471                       A_M(IJKB,west,M) = ZERO
472                       A_M(IJKB,north,M) = ZERO
473                       A_M(IJKB,south,M) = ZERO
474                       A_M(IJKB,top,M) = ZERO
475                       A_M(IJKB,bottom,M) = ZERO
476                       A_M(IJKB,0,M) = -ONE
477                       IF(BC_W_s(BCV,M)/=UNDEFINED) THEN
478                          B_M(IJKB,M) = - BC_W_s(BCV,M)
479                       ELSE
480                          B_M(IJKB,M) = - BC_VELMAG_s(BCV,M)*NORMAL_W(IJK,3)
481                       ENDIF
482                    ENDIF
483     
484     
485                 CASE (CG_PO)
486                    A_M(IJK,east,M) = ZERO
487                    A_M(IJK,west,M) = ZERO
488                    A_M(IJK,north,M) = ZERO
489                    A_M(IJK,south,M) = ZERO
490                    A_M(IJK,top,M) = ZERO
491                    A_M(IJK,bottom,M) = ZERO
492                    A_M(IJK,0,M) = -ONE
493                    B_M(IJK,M) = ZERO
494                    IJKB = BOTTOM_OF(IJK)
495                    IF(FLUID_AT(IJKB)) THEN
496                       A_M(IJK,bottom,M) = ONE
497                       A_M(IJK,0,M) = -ONE
498                    ENDIF
499              END SELECT
500     
501     
502     ! this is straight repeat of section of above code...?
503              BCV = BC_ID(IJK)
504              IF(BCV > 0 ) THEN
505                 BCT = BC_TYPE_ENUM(BCV)
506              ELSE
507                 BCT = NONE
508              ENDIF
509     
510     
511              SELECT CASE (BCT)
512                 CASE (CG_MI)
513                    A_M(IJK,east,M) = ZERO
514                    A_M(IJK,west,M) = ZERO
515                    A_M(IJK,north,M) = ZERO
516                    A_M(IJK,south,M) = ZERO
517                    A_M(IJK,top,M) = ZERO
518                    A_M(IJK,bottom,M) = ZERO
519                    A_M(IJK,0,M) = -ONE
520                    IF(BC_W_s(BCV,M)/=UNDEFINED) THEN
521                       B_M(IJK,M) = - BC_W_s(BCV,M)
522                    ELSE
523                       B_M(IJK,M) = - BC_VELMAG_s(BCV,M)*NORMAL_S(IJK,3)
524                    ENDIF
525                    IJKB = BOTTOM_OF(IJK)
526                    IF(FLUID_AT(IJKB)) THEN
527                       A_M(IJKB,east,M) = ZERO
528                       A_M(IJKB,west,M) = ZERO
529                       A_M(IJKB,north,M) = ZERO
530                       A_M(IJKB,south,M) = ZERO
531                       A_M(IJKB,top,M) = ZERO
532                       A_M(IJKB,bottom,M) = ZERO
533                       A_M(IJKB,0,M) = -ONE
534                       IF(BC_W_s(BCV,M)/=UNDEFINED) THEN
535                          B_M(IJKB,M) = - BC_W_s(BCV,M)
536                       ELSE
537                          B_M(IJKB,M) = - BC_VELMAG_s(BCV,M)*NORMAL_S(IJK,3)
538                       ENDIF
539                    ENDIF
540     
541     
542                 CASE (CG_PO)
543                    A_M(IJK,east,M) = ZERO
544                    A_M(IJK,west,M) = ZERO
545                    A_M(IJK,north,M) = ZERO
546                    A_M(IJK,south,M) = ZERO
547                    A_M(IJK,top,M) = ZERO
548                    A_M(IJK,bottom,M) = ZERO
549                    A_M(IJK,0,M) = -ONE
550                    B_M(IJK,M) = ZERO
551                    IJKB = BOTTOM_OF(IJK)
552                    IF(FLUID_AT(IJKB)) THEN
553                       A_M(IJK,bottom,M) = ONE
554                       A_M(IJK,0,M) = -ONE
555                    ENDIF
556     
557              END SELECT
558     
559           ENDDO
560     
561           RETURN
562     
563         CONTAINS
564     
565           INCLUDE 'functions.inc'
566     
567           END SUBROUTINE CG_SOURCE_W_S_BC
568