File: /nfs/home/0/users/jenkins/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, IER)
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     ! Error index
351           INTEGER, INTENT(INOUT) :: IER
352     !-----------------------------------------------
353     ! Local variables
354     !-----------------------------------------------
355     ! Indices
356           INTEGER :: IJK, IJKW
357     ! for cartesian grid:
358           INTEGER :: BCV
359           CHARACTER(LEN=9) :: BCT
360     !-----------------------------------------------
361     
362           IF(CG_SAFE_MODE(3)==1) RETURN
363     
364           DO IJK = ijkstart3, ijkend3
365     
366              BCV = BC_U_ID(IJK)
367              IF(BCV > 0 ) THEN
368                 BCT = BC_TYPE(BCV)
369              ELSE
370                 BCT = 'NONE'
371              ENDIF
372     
373              SELECT CASE (BCT)
374                 CASE ('CG_NSW')
375                    IF(WALL_U_AT(IJK)) THEN
376                       A_M(IJK,E,M) = ZERO
377                       A_M(IJK,W,M) = ZERO
378                       A_M(IJK,N,M) = ZERO
379                       A_M(IJK,S,M) = ZERO
380                       A_M(IJK,T,M) = ZERO
381                       A_M(IJK,B,M) = ZERO
382                       A_M(IJK,0,M) = -ONE
383                       B_M(IJK,M) = ZERO
384                    ENDIF
385     
386     
387                 CASE ('CG_FSW')
388                    IF(WALL_U_AT(IJK)) THEN
389                       A_M(IJK,E,M) = ZERO
390                       A_M(IJK,W,M) = ZERO
391                       A_M(IJK,N,M) = ZERO
392                       A_M(IJK,S,M) = ZERO
393                       A_M(IJK,T,M) = ZERO
394                       A_M(IJK,B,M) = ZERO
395                       A_M(IJK,0,M) = -ONE
396     !                  B_M(IJK,M) = - U_s(U_MASTER_OF(IJK),M)  ! Velocity of master node
397                       B_M(IJK,M) = ZERO
398                       IF(DABS(NORMAL_U(IJK,1))/=ONE) THEN
399                          IF (U_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
400                             A_M(IJK,E,M) = ONE
401                          ELSEIF (U_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
402                             A_M(IJK,W,M) = ONE
403                          ELSEIF (U_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
404                             A_M(IJK,N,M) = ONE
405                          ELSEIF (U_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
406                             A_M(IJK,S,M) = ONE
407                          ELSEIF (U_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
408                             A_M(IJK,T,M) = ONE
409                          ELSEIF (U_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
410                             A_M(IJK,B,M) = ONE
411                          ENDIF
412                       ENDIF
413                    ENDIF
414     
415     
416                 CASE ('CG_PSW')
417                    IF(WALL_U_AT(IJK)) THEN
418                       A_M(IJK,E,M) = ZERO
419                       A_M(IJK,W,M) = ZERO
420                       A_M(IJK,N,M) = ZERO
421                       A_M(IJK,S,M) = ZERO
422                       A_M(IJK,T,M) = ZERO
423                       A_M(IJK,B,M) = ZERO
424                       A_M(IJK,0,M) = -ONE
425     
426                       IF(BC_HW_S(BCV,M)==UNDEFINED) THEN   ! same as NSW
427                          B_M(IJK,M) = -BC_UW_S(BCV,M)
428                       ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN   ! same as FSW
429                          B_M(IJK,M) = ZERO
430                          IF(DABS(NORMAL_U(IJK,1))/=ONE) THEN
431                             IF (U_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
432                                A_M(IJK,E,M) = ONE
433                             ELSEIF (U_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
434                                A_M(IJK,W,M) = ONE
435                             ELSEIF (U_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
436                                A_M(IJK,N,M) = ONE
437                             ELSEIF (U_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
438                                A_M(IJK,S,M) = ONE
439                             ELSEIF (U_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
440                                A_M(IJK,T,M) = ONE
441                             ELSEIF (U_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
442                                A_M(IJK,B,M) = ONE
443                             ENDIF
444                          ENDIF
445                       ELSE                              ! partial slip  (WARNING:currently same as FSW)
446                          B_M(IJK,M) = ZERO
447                          IF(DABS(NORMAL_U(IJK,1))/=ONE) THEN
448                             IF (U_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
449                                A_M(IJK,E,M) = ONE
450                             ELSEIF (U_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
451                                A_M(IJK,W,M) = ONE
452                             ELSEIF (U_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
453                                A_M(IJK,N,M) = ONE
454                             ELSEIF (U_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
455                                A_M(IJK,S,M) = ONE
456                             ELSEIF (U_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
457                                A_M(IJK,T,M) = ONE
458                             ELSEIF (U_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
459                                A_M(IJK,B,M) = ONE
460                             ENDIF
461                          ENDIF
462                       ENDIF
463                    ENDIF
464     
465     
466                 CASE ('CG_MI')
467                    A_M(IJK,E,M) = ZERO
468                    A_M(IJK,W,M) = ZERO
469                    A_M(IJK,N,M) = ZERO
470                    A_M(IJK,S,M) = ZERO
471                    A_M(IJK,T,M) = ZERO
472                    A_M(IJK,B,M) = ZERO
473                    A_M(IJK,0,M) = -ONE
474     
475                    IF(BC_U_s(BCV,M)/=UNDEFINED) THEN
476                       B_M(IJK,M) = - BC_U_s(BCV,M)
477                    ELSE
478                       B_M(IJK,M) = - BC_VELMAG_s(BCV,M)*NORMAL_U(IJK,1)
479                    ENDIF
480     
481                    IJKW = WEST_OF(IJK)
482                    IF(FLUID_AT(IJKW)) THEN
483                       A_M(IJKW,E,M) = ZERO
484                       A_M(IJKW,W,M) = ZERO
485                       A_M(IJKW,N,M) = ZERO
486                       A_M(IJKW,S,M) = ZERO
487                       A_M(IJKW,T,M) = ZERO
488                       A_M(IJKW,B,M) = ZERO
489                       A_M(IJKW,0,M) = -ONE
490                       IF(BC_U_s(BCV,M)/=UNDEFINED) THEN
491                          B_M(IJKW,M) = - BC_U_s(BCV,M)
492                       ELSE
493                          B_M(IJKW,M) = - BC_VELMAG_s(BCV,M)*NORMAL_U(IJK,1)
494                       ENDIF
495                    ENDIF
496     
497     
498                 CASE ('CG_PO')
499                    A_M(IJK,E,M) = ZERO
500                    A_M(IJK,W,M) = ZERO
501                    A_M(IJK,N,M) = ZERO
502                    A_M(IJK,S,M) = ZERO
503                    A_M(IJK,T,M) = ZERO
504                    A_M(IJK,B,M) = ZERO
505                    A_M(IJK,0,M) = -ONE
506                    B_M(IJK,M) = ZERO
507                    IJKW = WEST_OF(IJK)
508                    IF(FLUID_AT(IJKW)) THEN
509                       A_M(IJK,W,M) = ONE
510                       A_M(IJK,0,M) = -ONE
511                    ENDIF
512     
513              END SELECT   ! end select bc_type
514     
515     
516     ! this is straight repeat of section of above code...?
517              BCV = BC_ID(IJK)
518              IF(BCV > 0 ) THEN
519                 BCT = BC_TYPE(BCV)
520              ELSE
521                 BCT = 'NONE'
522              ENDIF
523     
524     
525              SELECT CASE (BCT)
526                 CASE ('CG_MI')
527                    A_M(IJK,E,M) = ZERO
528                    A_M(IJK,W,M) = ZERO
529                    A_M(IJK,N,M) = ZERO
530                    A_M(IJK,S,M) = ZERO
531                    A_M(IJK,T,M) = ZERO
532                    A_M(IJK,B,M) = ZERO
533                    A_M(IJK,0,M) = -ONE
534     
535                    IF(BC_U_s(BCV,M)/=UNDEFINED) THEN
536                       B_M(IJK,M) = - BC_U_s(BCV,M)
537                    ELSE
538                       B_M(IJK,M) = - BC_VELMAG_s(BCV,M)*NORMAL_S(IJK,1)
539                    ENDIF
540     
541                    IJKW = WEST_OF(IJK)
542                    IF(FLUID_AT(IJKW)) THEN
543                       A_M(IJKW,E,M) = ZERO
544                       A_M(IJKW,W,M) = ZERO
545                       A_M(IJKW,N,M) = ZERO
546                       A_M(IJKW,S,M) = ZERO
547                       A_M(IJKW,T,M) = ZERO
548                       A_M(IJKW,B,M) = ZERO
549                       A_M(IJKW,0,M) = -ONE
550                       IF(BC_U_s(BCV,M)/=UNDEFINED) THEN
551                          B_M(IJKW,M) = - BC_U_s(BCV,M)
552                       ELSE
553                          B_M(IJKW,M) = - BC_VELMAG_s(BCV,M)*NORMAL_S(IJK,1)
554                       ENDIF
555                    ENDIF
556     
557     
558                 CASE ('CG_PO')
559                    A_M(IJK,E,M) = ZERO
560                    A_M(IJK,W,M) = ZERO
561                    A_M(IJK,N,M) = ZERO
562                    A_M(IJK,S,M) = ZERO
563                    A_M(IJK,T,M) = ZERO
564                    A_M(IJK,B,M) = ZERO
565                    A_M(IJK,0,M) = -ONE
566                    B_M(IJK,M) = ZERO
567                    IJKW = WEST_OF(IJK)
568                    IF(FLUID_AT(IJKW)) THEN
569                       A_M(IJK,W,M) = ONE
570                       A_M(IJK,0,M) = -ONE
571                    ENDIF
572     
573              END SELECT
574     
575           ENDDO
576     
577           RETURN
578           END SUBROUTINE CG_SOURCE_U_S_BC
579