File: /nfs/home/0/users/jenkins/mfix.git/model/cartesian_grid/CG_source_w_g.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CG_SOURCE_W_g(A_m, B_m, IER)                           C
4     !  Purpose: Determine contribution of cut-cell to source terms         C
5     !  for W_g momentum eq.                                                C
6     !                                                                      C
7     !                                                                      C
8     !  Author: Jeff Dietiker                              Date: 01-MAY-09  C
9     !  Reviewer:                                          Date:            C
10     !                                                                      C
11     !                                                                      C
12     !  Literature/Document References:                                     C
13     !                                                                      C
14     !  Variables referenced:                                               C
15     !  Variables modified:                                                 C
16     !                                                                      C
17     !  Local variables:                                                    C
18     !                                                                      C
19     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20     !
21           SUBROUTINE CG_SOURCE_W_G(A_M, B_M)
22     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
23     !...Switches: -xf
24     !
25     !  Include param.inc file to specify parameter values
26     !
27     !-----------------------------------------------
28     !   M o d u l e s
29     !-----------------------------------------------
30           USE param
31           USE param1
32           USE parallel
33           USE matrix
34           USE scales
35           USE constant
36           USE physprop
37           USE fldvar
38           USE visc_g
39           USE rxns
40           USE run
41           USE toleranc
42           USE geometry
43           USE indices
44           USE is
45           USE tau_g
46           USE bc
47           USE compar
48           USE sendrecv
49           USE ghdtheory
50           USE drag
51           USE fun_avg
52           USE functions
53     !=======================================================================
54     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
55     !=======================================================================
56           USE cutcell
57           USE quadric
58     !=======================================================================
59     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
60     !=======================================================================
61     
62           IMPLICIT NONE
63     !-----------------------------------------------
64     !   G l o b a l   P a r a m e t e r s
65     !-----------------------------------------------
66     !-----------------------------------------------
67     !   D u m m y   A r g u m e n t s
68     !-----------------------------------------------
69     !
70     !
71     !                      Error index
72           INTEGER          IER
73     !
74     !                      Indices
75           INTEGER          I, J, K, IJK, IJKT, IMJK, IJKP, IMJKP,&
76                            IJKE, IJKW, IJKTE, IM, IPJK
77     !
78     !                      Phase index
79           INTEGER          M
80     !
81     !                      Average volume fraction
82           DOUBLE PRECISION EPGA
83     !
84     !                      Septadiagonal matrix A_m
85           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
86     !
87     !                      Vector b_m
88           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
89     !
90     !=======================================================================
91     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
92     !=======================================================================
93           INTEGER :: JM,IP,JP,IJMK,IJPK,IJKC,IJKN,IJKNE,IJKS,IJKSE,IPJMK,IJKM,KM,KP,IJMKP
94           INTEGER :: IJKTN,IJKWT,IJKST
95           DOUBLE PRECISION :: We,Ww,Wn,Ws,Wt,Wb
96           DOUBLE PRECISION :: B_NOC
97           DOUBLE PRECISION :: MU_GT_E,MU_GT_W,MU_GT_N,MU_GT_S,MU_GT_T,MU_GT_B,MU_GT_CUT
98           DOUBLE PRECISION :: WW_g
99           INTEGER :: BCV
100           CHARACTER(LEN=9) :: BCT
101     !                       virtual (added) mass
102           DOUBLE PRECISION F_vir, ROP_MA, U_se, Usw, Wse, Wsw, Wsn, Wss, Wst, Wsb, Usc,Vsc,Vsn,Vss
103     ! Wall function
104           DOUBLE PRECISION :: W_F_Slip
105     !=======================================================================
106     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
107     !=======================================================================
108     !-----------------------------------------------
109     
110           IF(CG_SAFE_MODE(5)==1) RETURN
111     
112     !
113           M = 0
114           IF (.NOT.MOMENTUM_Z_EQ(0)) RETURN
115     !
116     !
117     !!!$omp  parallel do private( I, J, K, IJK, IJKT, ISV, Sdp, V0, Vpm, Vmt, Vbf, &
118     !!!$omp&  PGT, ROGA, IMJK, IJKP, IMJKP, IJKW, IJKTE, IJKTW, IM, IPJK,  &
119     !!!$omp&  CTE, CTW, SXZB, EPMUOX, VXZA, VXZB, UGT, VCOA, VCOB, IJKE,&
120     !!!$omp&  MUGA, ROPGA, EPGA, LINE)
121           DO IJK = ijkstart3, ijkend3
122              I = I_OF(IJK)
123              J = J_OF(IJK)
124              K = K_OF(IJK)
125              IJKT = TOP_OF(IJK)
126              EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
127              IF (IP_AT_T(IJK)) THEN
128     !
129     !        do nothing
130     !
131     !       dilute flow
132              ELSE IF (EPGA <= DIL_EP_S) THEN
133     !
134     !        do nothing
135     !
136     !         ELSE
137              ELSEIF(INTERIOR_CELL_AT(IJK)) THEN
138     !
139                 BCV = BC_W_ID(IJK)
140     
141                 IF(BCV > 0 ) THEN
142                    BCT = BC_TYPE(BCV)
143                 ELSE
144                    BCT = 'NONE'
145                 ENDIF
146     
147                 SELECT CASE (BCT)
148                    CASE ('CG_NSW')
149                       NOC_WG = .TRUE.
150                       WW_g = ZERO
151                       MU_GT_CUT =  (VOL(IJK)*MU_GT(IJK) + VOL(IJKT)*MU_GT(IJKT))/(VOL(IJK) + VOL(IJKT))
152     
153                       IF(.NOT.K_EPSILON) THEN
154                          A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_W_CUT(IJK)/DELH_W(IJK)
155                       ELSE
156                          CALL Wall_Function(IJK,IJK,ONE/DELH_W(IJK),W_F_Slip)
157                          A_M(IJK,0,M) = A_M(IJK,0,M)  - MU_GT_CUT * Area_W_CUT(IJK)*(ONE-W_F_Slip)/DELH_W(IJK)
158                       ENDIF
159     
160                    CASE ('CG_FSW')
161                       NOC_WG = .FALSE.
162                       WW_g = ZERO
163                    CASE('CG_PSW')
164                       IF(BC_HW_G(BCV)==UNDEFINED) THEN   ! same as NSW
165                          NOC_WG = .TRUE.
166                          WW_g = BC_WW_G(BCV)
167                          MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJKT)*MU_GT(IJKT))/(VOL(IJK) + VOL(IJKT))
168                          A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_W_CUT(IJK)/DELH_W(IJK)
169                          B_M(IJK,M) = B_M(IJK,M) - MU_GT_CUT * WW_g * Area_W_CUT(IJK)/DELH_W(IJK)
170                       ELSEIF(BC_HW_G(BCV)==ZERO) THEN   ! same as FSW
171                          NOC_WG = .FALSE.
172                          WW_g = ZERO
173                       ELSE                              ! partial slip
174                          NOC_WG = .FALSE.
175                          WW_g = BC_WW_G(BCV)
176                          MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJKT)*MU_GT(IJKT))/(VOL(IJK) + VOL(IJKT))
177                          A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_W_CUT(IJK)*(BC_HW_G(BCV))
178                          B_M(IJK,M) = B_M(IJK,M) - MU_GT_CUT * WW_g * Area_W_CUT(IJK)*(BC_HW_G(BCV))
179     
180                       ENDIF
181                    CASE ('NONE')
182                       NOC_WG = .FALSE.
183                 END SELECT
184     
185                 IF(NOC_WG) THEN
186     
187                    IMJK = IM_OF(IJK)
188                    IJMK = JM_OF(IJK)
189                    IJKM = KM_OF(IJK)
190                    IPJK = IP_OF(IJK)
191                    IJPK = JP_OF(IJK)
192                    IJKP = KP_OF(IJK)
193     
194                    We = Theta_We_bar(IJK)  * W_g(IJK)  + Theta_We(IJK)  * W_g(IPJK)
195                    Ww = Theta_We_bar(IMJK) * W_g(IMJK) + Theta_We(IMJK) * W_g(IJK)
196     
197                    Wn = Theta_Wn_bar(IJK)  * W_g(IJK)  + Theta_Wn(IJK)  * W_g(IJPK)
198                    Ws = Theta_Wn_bar(IJMK) * W_g(IJMK) + Theta_Wn(IJMK) * W_g(IJK)
199     
200                    Wt = Theta_Wt_bar(IJK)  * W_g(IJK)  + Theta_Wt(IJK)  * W_g(IJKP)
201                    Wb = Theta_Wt_bar(IJKM) * W_g(IJKM) + Theta_Wt(IJKM) * W_g(IJK)
202     
203                    IJKE = EAST_OF(IJK)
204     
205                    ijkt = top_of(ijk)
206     
207                    IF (WALL_AT(IJK)) THEN
208                       IJKC = IJKT
209                    ELSE
210                       IJKC = IJK
211                    ENDIF
212     
213                    IP = IP1(I)
214                    IM = IM1(I)
215                    IJKN = NORTH_OF(IJK)
216                    IJKNE = EAST_OF(IJKN)
217     
218                    JM = JM1(J)
219                    IPJMK = IP_OF(IJMK)
220                    IJKS = SOUTH_OF(IJK)
221                    IJKSE = EAST_OF(IJKS)
222     
223                    KP = KP1(K)
224                    IJKT = TOP_OF(IJK)
225                    IJKE = EAST_OF(IJK)
226                    IJKP = KP_OF(IJK)
227                    IJKTN = NORTH_OF(IJKT)
228                    IJKTE = EAST_OF(IJKT)
229                    IJKW = WEST_OF(IJK)
230                    IJKWT = TOP_OF(IJKW)
231                    IJKS = SOUTH_OF(IJK)
232                    IJKST = TOP_OF(IJKS)
233     
234                    MU_GT_E = AVG_Z_H(AVG_X_H(MU_GT(IJKC),MU_GT(IJKE),I),&
235                              AVG_X_H(MU_GT(IJKT),MU_GT(IJKTE),I),K)
236     
237                    MU_GT_W = AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJKC),IM),&
238                              AVG_X_H(MU_GT(IJKWT),MU_GT(IJKT),IM),K)
239     
240                    MU_GT_N = AVG_Z_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),&
241                              AVG_Y_H(MU_GT(IJKT),MU_GT(IJKTN),J),K)
242     
243                    MU_GT_S = AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),&
244                              AVG_Y_H(MU_GT(IJKST),MU_GT(IJKT),JM),K)
245     
246                    MU_GT_T = MU_GT(IJKT)
247                    MU_GT_B = MU_GT(IJKC)
248     
249                    B_NOC =     MU_GT_E * Ayz_W(IJK)  * (We-WW_g) * NOC_W_E(IJK)  &
250                            -   MU_GT_W * Ayz_W(IMJK) * (Ww-WW_g) * NOC_W_E(IMJK) &
251                            +   MU_GT_N * Axz_W(IJK)  * (Wn-WW_g) * NOC_W_N(IJK)  &
252                            -   MU_GT_S * Axz_W(IJMK) * (Ws-WW_g) * NOC_W_N(IJMK) &
253                            +   MU_GT_T * Axy_W(IJK)  * (Wt-WW_g) * NOC_W_T(IJK)  *2.0d0&
254                            -   MU_GT_B * Axy_W(IJKM) * (Wb-WW_g) * NOC_W_T(IJKM) *2.0D0
255     
256                    B_M(IJK,M) = B_M(IJK,M)   +  B_NOC
257                 ENDIF
258     
259                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
260     !
261     !!! BEGIN VIRTUAL MASS SECTION (explicit terms)
262     ! adding transient term  dWs/dt to virtual mass term
263                    F_vir = ZERO
264                    IF(Added_Mass) THEN
265     
266                       F_vir = ( (W_s(IJK,M_AM) - W_sO(IJK,M_AM)) )*ODT*VOL_W(IJK)
267     
268                       I = I_OF(IJK)
269                       J = J_OF(IJK)
270                       K = K_OF(IJK)
271     
272                       IM = I - 1
273                       JM = J - 1
274                       KM = K - 1
275     
276                       IP = I + 1
277                       JP = J + 1
278                       KP = K + 1
279     
280                       IMJK = FUNIJK(IM,J,K)
281                       IJMK = FUNIJK(I,JM,K)
282                       IPJK = FUNIJK(IP,J,K)
283                       IJPK = FUNIJK(I,JP,K)
284                       IJKP = FUNIJK(I,J,KP)
285                       IJKM = FUNIJK(I,J,KM)
286     
287                       IMJKP = KP_OF(IMJK)
288                       IJMKP = KP_OF(IJMK)
289     
290                       IJKE = EAST_OF(IJK)
291     !
292     ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
293     
294     
295                       Wse = Theta_We_bar(IJK) * W_s(IJK,M_AM) + Theta_We(IJK) * W_s(IPJK,M_AM)
296                       Wsw = Theta_We_bar(IMJK) * W_s(IMJK,M_AM) + Theta_We(IMJK) * W_s(IJK,M_AM)
297     
298                       U_se = Theta_W_te(IJK) * U_s(IJK,M_AM) + Theta_W_be(IJK) * U_s(IJKP,M_AM)
299                       Usw = Theta_W_te(IMJK) * U_s(IMJK,M_AM) + Theta_W_be(IMJK) * U_s(IMJKP,M_AM)
300     
301                       Usc = (DELX_we(IJK) * Usw + DELX_ww(IJK) * U_se) / (DELX_we(IJK) + DELX_ww(IJK))
302     
303     
304                       Wsn = Theta_Wn_bar(IJK) * W_s(IJK,M_AM) + Theta_Wn(IJK) * W_s(IJPK,M_AM)
305                       Wss = Theta_Wn_bar(IJMK) * W_s(IJMK,M_AM) + Theta_Wn(IJMK) * W_s(IJK,M_AM)
306     
307                       Vsn =  Theta_W_tn(IJK)  * V_s(IJK,M_AM)  + Theta_W_bn(IJK)  * V_s(IJKP,M_AM)
308                       Vss =  Theta_W_tn(IJMK) * V_s(IJMK,M_AM) + Theta_W_bn(IJMK) * V_s(IJMKP,M_AM)
309     
310                       Vsc = (DELY_wn(IJK) * Vss + DELY_ws(IJK) * Vsn) / (DELY_wn(IJK) + DELY_ws(IJK))
311     
312     
313                       Wst = Theta_Wt_bar(IJK)  * W_s(IJK,M_AM)  + Theta_Wt(IJK)  * W_s(IJKP,M_AM)
314                       Wsb = Theta_Wt_bar(IMJK) * W_s(IMJK,M_AM) + Theta_Wt(IMJK) * W_s(IMJKP,M_AM)
315     
316     !
317     ! adding convective terms (U dW/dx + V dW/dy + W dW/dz) to virtual mass
318     
319                       F_vir = F_vir +  Usc * (Wse - Wsw)*AYZ(IJK)    + &
320                                        Vsc * (Wsn - Wss)*AXZ(IJK)  + &
321                                        W_s(IJK,M_AM)*(Wst - Wsb)*AXY(IJK)
322     
323     
324                       ROP_MA = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M_AM) + VOL(IJKT)*ROP_g(IJKT)*EP_s(IJKT,M_AM))/(VOL(IJK) + VOL(IJKT))
325     
326                       F_vir = F_vir * Cv * ROP_MA
327     
328                       B_M(IJK,M) = B_M(IJK,M) - F_vir ! explicit part of virtual mass force
329     
330                    ENDIF
331     !
332     !!! END VIRTUAL MASS SECTION
333     
334                 ENDIF
335     
336              ENDIF
337           END DO
338     
339           RETURN
340           END SUBROUTINE CG_SOURCE_W_G
341     
342     !
343     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
344     !                                                                      C
345     !  Module name: CG_SOURCE_W_g_BC(A_m, B_m, IER)                        C
346     !  Purpose: Determine contribution of cut-cell to source terms         C
347     !  for W_g momentum eq.                                                C
348     !                                                                      C
349     !                                                                      C
350     !  Author: Jeff Dietiker                              Date: 01-MAY-09  C
351     !  Reviewer:                                          Date:            C
352     !                                                                      C
353     !                                                                      C
354     !  Literature/Document References:                                     C
355     !                                                                      C
356     !  Variables referenced:                                               C
357     !  Variables modified:                                                 C
358     !                                                                      C
359     !  Local variables:                                                    C
360     !                                                                      C
361     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
362     !
363           SUBROUTINE CG_SOURCE_W_G_BC(A_M, B_M)
364     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
365     !...Switches: -xf
366     !
367     !  Include param.inc file to specify parameter values
368     !
369     !-----------------------------------------------
370     !   M o d u l e s
371     !-----------------------------------------------
372           USE param
373           USE param1
374           USE parallel
375           USE matrix
376           USE scales
377           USE constant
378           USE physprop
379           USE fldvar
380           USE visc_g
381           USE rxns
382           USE run
383           USE toleranc
384           USE geometry
385           USE indices
386           USE is
387           USE tau_g
388           USE bc
389           USE output
390           USE compar
391           USE fun_avg
392           USE functions
393     
394     !=======================================================================
395     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
396     !=======================================================================
397           USE cutcell
398           USE quadric
399     !=======================================================================
400     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
401     !=======================================================================
402     
403     
404           IMPLICIT NONE
405     !-----------------------------------------------
406     !   G l o b a l   P a r a m e t e r s
407     !-----------------------------------------------
408     !-----------------------------------------------
409     !   D u m m y   A r g u m e n t s
410     !-----------------------------------------------
411     !
412     !                      Indices
413           INTEGER          IJK, IJKB
414     !
415     !                      Solids phase
416           INTEGER          M
417     !
418     !                      Septadiagonal matrix A_m
419           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
420     !
421     !                      Vector b_m
422           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
423     !
424     
425     !=======================================================================
426     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
427     !=======================================================================
428           INTEGER :: BCV
429           CHARACTER(LEN=9) :: BCT
430     
431     !-----------------------------------------------
432     
433           IF(CG_SAFE_MODE(5)==1) RETURN
434     
435     !
436           M = 0
437     
438           DO IJK = ijkstart3, ijkend3
439     
440              BCV = BC_W_ID(IJK)
441     
442              IF(BCV > 0 ) THEN
443                 BCT = BC_TYPE(BCV)
444              ELSE
445                 BCT = 'NONE'
446              ENDIF
447     
448              SELECT CASE (BCT)
449     
450                 CASE ('CG_NSW')
451     
452                    IF(WALL_W_AT(IJK)) THEN
453     
454                       A_M(IJK,E,M) = ZERO
455                       A_M(IJK,W,M) = ZERO
456                       A_M(IJK,N,M) = ZERO
457                       A_M(IJK,S,M) = ZERO
458                       A_M(IJK,T,M) = ZERO
459                       A_M(IJK,B,M) = ZERO
460                       A_M(IJK,0,M) = -ONE
461     
462                       B_M(IJK,M) = ZERO
463     
464                    ENDIF
465     
466                 CASE ('CG_FSW')
467     
468                    IF(WALL_W_AT(IJK)) THEN
469     
470                       A_M(IJK,E,M) = ZERO
471                       A_M(IJK,W,M) = ZERO
472                       A_M(IJK,N,M) = ZERO
473                       A_M(IJK,S,M) = ZERO
474                       A_M(IJK,T,M) = ZERO
475                       A_M(IJK,B,M) = ZERO
476                       A_M(IJK,0,M) = -ONE
477     
478     !                  B_M(IJK,M) = - W_g(W_MASTER_OF(IJK))  ! Velocity of master node
479     
480                       B_M(IJK,M) = ZERO
481     
482                       IF(DABS(NORMAL_W(IJK,3))/=ONE) THEN
483     
484                          IF (W_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
485                             A_M(IJK,E,M) = ONE
486                          ELSEIF (W_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
487                             A_M(IJK,W,M) = ONE
488                          ELSEIF (W_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
489                             A_M(IJK,N,M) = ONE
490                          ELSEIF (W_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
491                             A_M(IJK,S,M) = ONE
492                          ELSEIF (W_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
493                             A_M(IJK,T,M) = ONE
494                          ELSEIF (W_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
495                             A_M(IJK,B,M) = ONE
496                          ENDIF
497     
498                       ENDIF
499     
500                    ENDIF
501     
502                 CASE ('CG_PSW')
503     
504                    IF(WALL_W_AT(IJK)) THEN
505     
506                       A_M(IJK,E,M) = ZERO
507                       A_M(IJK,W,M) = ZERO
508                       A_M(IJK,N,M) = ZERO
509                       A_M(IJK,S,M) = ZERO
510                       A_M(IJK,T,M) = ZERO
511                       A_M(IJK,B,M) = ZERO
512                       A_M(IJK,0,M) = -ONE
513     
514     
515                       IF(BC_HW_G(BCV)==UNDEFINED) THEN   ! same as NSW
516                          B_M(IJK,M) = -BC_WW_G(BCV)
517                       ELSEIF(BC_HW_G(BCV)==ZERO) THEN   ! same as FSW
518                          B_M(IJK,M) = ZERO
519     
520                          IF(DABS(NORMAL_W(IJK,3))/=ONE) THEN
521     
522                             IF (W_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
523                                A_M(IJK,E,M) = ONE
524                             ELSEIF (W_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
525                                A_M(IJK,W,M) = ONE
526                             ELSEIF (W_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
527                                A_M(IJK,N,M) = ONE
528                             ELSEIF (W_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
529                                A_M(IJK,S,M) = ONE
530                             ELSEIF (W_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
531                                A_M(IJK,T,M) = ONE
532                             ELSEIF (W_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
533                                A_M(IJK,B,M) = ONE
534                             ENDIF
535     
536                          ENDIF
537     
538                       ELSE                              ! partial slip
539     
540     
541     
542     
543     
544                       ENDIF
545     
546                    ENDIF
547     
548     
549                 CASE ('CG_MI')
550     
551                    A_M(IJK,E,M) = ZERO
552                    A_M(IJK,W,M) = ZERO
553                    A_M(IJK,N,M) = ZERO
554                    A_M(IJK,S,M) = ZERO
555                    A_M(IJK,T,M) = ZERO
556                    A_M(IJK,B,M) = ZERO
557                    A_M(IJK,0,M) = -ONE
558     
559                    IF(BC_W_g(BCV)/=UNDEFINED) THEN
560                       B_M(IJK,M) = - BC_W_g(BCV)
561                    ELSE
562                       B_M(IJK,M) = - BC_VELMAG_g(BCV)*NORMAL_W(IJK,3)
563                    ENDIF
564     
565     
566                    IJKB = BOTTOM_OF(IJK)
567                    IF(FLUID_AT(IJKB)) THEN
568     
569                       A_M(IJKB,E,M) = ZERO
570                       A_M(IJKB,W,M) = ZERO
571                       A_M(IJKB,N,M) = ZERO
572                       A_M(IJKB,S,M) = ZERO
573                       A_M(IJKB,T,M) = ZERO
574                       A_M(IJKB,B,M) = ZERO
575                       A_M(IJKB,0,M) = -ONE
576     
577                       IF(BC_W_g(BCV)/=UNDEFINED) THEN
578                          B_M(IJKB,M) = - BC_W_g(BCV)
579                       ELSE
580                          B_M(IJKB,M) = - BC_VELMAG_g(BCV)*NORMAL_W(IJK,3)
581                       ENDIF
582     
583     
584                    ENDIF
585     
586                 CASE ('CG_PO')
587     
588                    A_M(IJK,E,M) = ZERO
589                    A_M(IJK,W,M) = ZERO
590                    A_M(IJK,N,M) = ZERO
591                    A_M(IJK,S,M) = ZERO
592                    A_M(IJK,T,M) = ZERO
593                    A_M(IJK,B,M) = ZERO
594                    A_M(IJK,0,M) = -ONE
595                    B_M(IJK,M) = ZERO
596     
597                    IJKB = BOTTOM_OF(IJK)
598                    IF(FLUID_AT(IJKB)) THEN
599     
600                       A_M(IJK,B,M) = ONE
601                       A_M(IJK,0,M) = -ONE
602     
603                    ENDIF
604     
605              END SELECT
606     
607              BCV = BC_ID(IJK)
608     
609              IF(BCV > 0 ) THEN
610                 BCT = BC_TYPE(BCV)
611              ELSE
612                 BCT = 'NONE'
613              ENDIF
614     
615              SELECT CASE (BCT)
616     
617                 CASE ('CG_MI')
618     
619                    A_M(IJK,E,M) = ZERO
620                    A_M(IJK,W,M) = ZERO
621                    A_M(IJK,N,M) = ZERO
622                    A_M(IJK,S,M) = ZERO
623                    A_M(IJK,T,M) = ZERO
624                    A_M(IJK,B,M) = ZERO
625                    A_M(IJK,0,M) = -ONE
626     
627                    IF(BC_W_g(BCV)/=UNDEFINED) THEN
628                       B_M(IJK,M) = - BC_W_g(BCV)
629                    ELSE
630                       B_M(IJK,M) = - BC_VELMAG_g(BCV)*NORMAL_S(IJK,3)
631                    ENDIF
632     
633     
634                    IJKB = BOTTOM_OF(IJK)
635                    IF(FLUID_AT(IJKB)) THEN
636     
637                       A_M(IJKB,E,M) = ZERO
638                       A_M(IJKB,W,M) = ZERO
639                       A_M(IJKB,N,M) = ZERO
640                       A_M(IJKB,S,M) = ZERO
641                       A_M(IJKB,T,M) = ZERO
642                       A_M(IJKB,B,M) = ZERO
643                       A_M(IJKB,0,M) = -ONE
644     
645                       IF(BC_W_g(BCV)/=UNDEFINED) THEN
646                          B_M(IJKB,M) = - BC_W_g(BCV)
647                       ELSE
648                          B_M(IJKB,M) = - BC_VELMAG_g(BCV)*NORMAL_S(IJK,3)
649                       ENDIF
650     
651     
652                    ENDIF
653     
654                 CASE ('CG_PO')
655     
656                    A_M(IJK,E,M) = ZERO
657                    A_M(IJK,W,M) = ZERO
658                    A_M(IJK,N,M) = ZERO
659                    A_M(IJK,S,M) = ZERO
660                    A_M(IJK,T,M) = ZERO
661                    A_M(IJK,B,M) = ZERO
662                    A_M(IJK,0,M) = -ONE
663                    B_M(IJK,M) = ZERO
664     
665                    IJKB = BOTTOM_OF(IJK)
666                    IF(FLUID_AT(IJKB)) THEN
667     
668                       A_M(IJK,B,M) = ONE
669                       A_M(IJK,0,M) = -ONE
670     
671                    ENDIF
672     
673              END SELECT
674     
675     
676           ENDDO
677     
678           RETURN
679     
680     !=======================================================================
681     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
682     !=======================================================================
683     
684           END SUBROUTINE CG_SOURCE_W_G_BC
685     
686     !// Comments on the modifications for DMP version implementation
687     !// 001 Include header file and common declarations for parallelization
688     !// 350 Changed do loop limits: 1,kmax2->kmin3,kmax3
689     !// 360 Check if i,j,k resides on current processor
690