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