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