File: RELATIVE:/../../../mfix.git/model/cartesian_grid/CG_source_u_s.f

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