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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CG_SOURCE_U_g(A_m, B_m, IER)                           C
4     !  Purpose: Determine contribution of cut-cell to source terms         C
5     !  for U_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_U_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 fun_avg
50           USE functions
51     
52     !=======================================================================
53     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
54     !=======================================================================
55           USE cutcell
56           USE quadric
57     !=======================================================================
58     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
59     !=======================================================================
60     
61           IMPLICIT NONE
62     !-----------------------------------------------
63     !   G l o b a l   P a r a m e t e r s
64     !-----------------------------------------------
65     !-----------------------------------------------
66     !   D u m m y   A r g u m e n t s
67     !-----------------------------------------------
68     !
69     !                      Indices
70           INTEGER          I, IJK, IJKE, IPJK, IJKM, IPJKM
71     !
72     !                      Phase index
73           INTEGER          M
74     !
75     !                      Average volume fraction
76           DOUBLE PRECISION EPGA
77     !
78     !                      Septadiagonal matrix A_m
79           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
80     !
81     !                      Vector b_m
82           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
83     !
84     !=======================================================================
85     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
86     !=======================================================================
87           INTEGER :: J,K,IM,JM,IP,JP,KM,KP
88           INTEGER :: IMJK,IJMK,IJPK,IJKC,IJKN,IJKNE,IJKS,IJKSE,IPJMK,IJKP
89           INTEGER :: IJKT,IJKTE,IJKB,IJKBE
90           DOUBLE PRECISION :: Ue,Uw,Un,Us,Ut,Ub
91           DOUBLE PRECISION :: B_NOC
92           DOUBLE PRECISION :: MU_GT_E,MU_GT_W,MU_GT_N,MU_GT_S,MU_GT_T,MU_GT_B,MU_GT_CUT
93           DOUBLE PRECISION :: UW_g
94           INTEGER :: BCV
95           CHARACTER(LEN=9) :: BCT
96     
97     !                       virtual (added) mass
98           DOUBLE PRECISION F_vir, ROP_MA, U_se, Usw, Vsn, Vss, Vsc, Usn, Uss, Wsb, Wst, Wsc, Usb, Ust
99     ! Wall function
100           DOUBLE PRECISION :: W_F_Slip
101     !=======================================================================
102     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
103     !=======================================================================
104     
105           IF(CG_SAFE_MODE(3)==1) RETURN
106     
107     
108     !
109           M = 0
110           IF (.NOT.MOMENTUM_X_EQ(0)) RETURN
111     !
112     !!!$omp    parallel do private(I, IJK, IJKE, IJKM, IPJK, IPJKM,     &
113     !!!$omp&                  ISV, Sdp, V0, Vpm, Vmt, Vbf,              &
114     !!!$omp&                  Vcf, EPMUGA, VTZA, WGE, PGE, ROGA,        &
115     !!!$omp&                  MUGA, ROPGA, EPGA )
116           DO IJK = ijkstart3, ijkend3
117              I = I_OF(IJK)
118              IJKE = EAST_OF(IJK)
119              IJKM = KM_OF(IJK)
120              IPJK = IP_OF(IJK)
121              IPJKM = IP_OF(IJKM)
122              EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
123              IF (IP_AT_E(IJK)) THEN
124     !
125     !        do nothing
126     !
127     !       dilute flow
128              ELSE IF (EPGA <= DIL_EP_S) THEN
129     !
130     !        do nothing
131     !
132     !        ELSE
133              ELSEIF(INTERIOR_CELL_AT(IJK)) THEN
134     !
135                 BCV = BC_U_ID(IJK)
136     
137                 IF(BCV > 0 ) THEN
138                    BCT = BC_TYPE(BCV)
139                 ELSE
140                    BCT = 'NONE'
141                 ENDIF
142     
143                 SELECT CASE (BCT)
144                    CASE ('CG_NSW')
145                       NOC_UG = .TRUE.
146                       UW_g = ZERO
147                       MU_GT_CUT =  (VOL(IJK)*MU_GT(IJK) + VOL(IPJK)*MU_GT(IJKE))/(VOL(IJK) + VOL(IPJK))
148     
149                       IF(.NOT.K_EPSILON) THEN
150                          A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_U_CUT(IJK)/DELH_U(IJK)
151                       ELSE
152                          CALL Wall_Function(IJK,IJK,ONE/DELH_U(IJK),W_F_Slip)
153                          A_M(IJK,0,M) = A_M(IJK,0,M)  - MU_GT_CUT * Area_U_CUT(IJK)*(ONE-W_F_Slip)/DELH_U(IJK)
154                       ENDIF
155                    CASE ('CG_FSW')
156                       NOC_UG = .FALSE.
157                       UW_g = ZERO
158                    CASE('CG_PSW')
159                       IF(BC_HW_G(BCV)==UNDEFINED) THEN   ! same as NSW
160                          NOC_UG = .TRUE.
161                          UW_g = BC_UW_G(BCV)
162                          MU_GT_CUT =  (VOL(IJK)*MU_GT(IJK) + VOL(IPJK)*MU_GT(IJKE))/(VOL(IJK) + VOL(IPJK))
163                          A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_U_CUT(IJK)/DELH_U(IJK)
164                          B_M(IJK,M) = B_M(IJK,M) - MU_GT_CUT * UW_g * Area_U_CUT(IJK)/DELH_U(IJK)
165                       ELSEIF(BC_HW_G(BCV)==ZERO) THEN   ! same as FSW
166                          NOC_UG = .FALSE.
167                          UW_g = ZERO
168                       ELSE                              ! partial slip
169                          NOC_UG = .FALSE.
170                          UW_g = BC_UW_G(BCV)
171                          MU_GT_CUT =  (VOL(IJK)*MU_GT(IJK) + VOL(IPJK)*MU_GT(IJKE))/(VOL(IJK) + VOL(IPJK))
172                          A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_U_CUT(IJK)*(BC_HW_G(BCV))
173                          B_M(IJK,M) = B_M(IJK,M) - MU_GT_CUT * UW_g * Area_U_CUT(IJK)*(BC_HW_G(BCV))
174                       ENDIF
175                    CASE ('NONE', 'CG_MI')
176                       NOC_UG = .FALSE.
177                 END SELECT
178     
179                 IF(NOC_UG) THEN
180     
181                    I = I_OF(IJK)
182                    J = J_OF(IJK)
183                    K = K_OF(IJK)
184     
185                    IMJK = IM_OF(IJK)
186                    IJMK = JM_OF(IJK)
187                    IJKM = KM_OF(IJK)
188                    IPJK = IP_OF(IJK)
189                    IJPK = JP_OF(IJK)
190                    IJKP = KP_OF(IJK)
191     
192                    Ue = Theta_Ue_bar(IJK)  * U_g(IJK)  + Theta_Ue(IJK)  * U_g(IPJK)
193                    Uw = Theta_Ue_bar(IMJK) * U_g(IMJK) + Theta_Ue(IMJK) * U_g(IJK)
194     
195                    Un = Theta_Un_bar(IJK)  * U_g(IJK)  + Theta_Un(IJK)  * U_g(IJPK)
196                    Us = Theta_Un_bar(IJMK) * U_g(IJMK) + Theta_Un(IJMK) * U_g(IJK)
197     
198     
199                    IJKE = EAST_OF(IJK)
200     
201                    IF (WALL_AT(IJK)) THEN
202                       IJKC = IJKE
203                    ELSE
204                       IJKC = IJK
205                    ENDIF
206                    IP = IP1(I)
207                    IJKN = NORTH_OF(IJK)
208                    IJKNE = EAST_OF(IJKN)
209     
210                    JM = JM1(J)
211                    IPJMK = IP_OF(IJMK)
212                    IJKS = SOUTH_OF(IJK)
213                    IJKSE = EAST_OF(IJKS)
214     
215                    KM = KM1(K)
216     
217                    IJKT = TOP_OF(IJK)
218                    IJKTE = EAST_OF(IJKT)
219     
220                    IJKB = BOTTOM_OF(IJK)
221                    IJKBE = EAST_OF(IJKB)
222     
223                    MU_GT_E = MU_GT(IJKE)
224                    MU_GT_W = MU_GT(IJKC)
225     
226                    MU_GT_N = AVG_X_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),&
227                                      AVG_Y_H(MU_GT(IJKE),MU_GT(IJKNE),J),I)
228     
229                    MU_GT_S = AVG_X_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),&
230                                      AVG_Y_H(MU_GT(IJKSE),MU_GT(IJKE),JM),I)
231     
232                    B_NOC =     MU_GT_E * Ayz_U(IJK)  * (Ue-UW_g) * NOC_U_E(IJK)  *2.0d0&
233                            -   MU_GT_W * Ayz_U(IMJK) * (Uw-UW_g) * NOC_U_E(IMJK) *2.0d0&
234                            +   MU_GT_N * Axz_U(IJK)  * (Un-UW_g) * NOC_U_N(IJK)  &
235                            -   MU_GT_S * Axz_U(IJMK) * (Us-UW_g) * NOC_U_N(IJMK)
236     
237                    IF(DO_K) THEN
238     
239                       Ut = Theta_Ut_bar(IJK)  * U_g(IJK)  + Theta_Ut(IJK)  * U_g(IJKP)
240                       Ub = Theta_Ut_bar(IJKM) * U_g(IJKM) + Theta_Ut(IJKM) * U_g(IJK)
241     
242                       MU_GT_T = AVG_X_H(AVG_Z_H(MU_GT(IJKC),MU_GT(IJKT),K),&
243                                 AVG_Z_H(MU_GT(IJKE),MU_GT(IJKTE),K),I)
244     
245                       MU_GT_B = AVG_X_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJKC),KM),&
246                                 AVG_Z_H(MU_GT(IJKBE),MU_GT(IJKE),KM),I)
247     
248                       B_NOC = B_NOC  +   MU_GT_T * Axy_U(IJK)  * (Ut-UW_g) * NOC_U_T(IJK)  &
249                                      -   MU_GT_B * Axy_U(IJKM) * (Ub-UW_g) * NOC_U_T(IJKM)
250                    ENDIF
251     
252                    B_M(IJK,M) = B_M(IJK,M)   +  B_NOC
253     
254     
255                 ENDIF
256     
257                 IF(CUT_U_TREATMENT_AT(IJK)) THEN
258     !
259     !!! BEGIN VIRTUAL MASS SECTION (explicit terms)
260     ! adding transient term  dUs/dt to virtual mass term
261                    F_vir = ZERO
262                    IF(Added_Mass) THEN
263     
264                       F_vir = ( (U_s(IJK,M_AM) - U_sO(IJK,M_AM)) )*ODT*VOL_U(IJK)
265     
266                       J = J_OF(IJK)
267                       K = K_OF(IJK)
268     
269                       IM = I - 1
270                       JM = J - 1
271                       KM = K - 1
272     
273                       IP = I + 1
274                       JP = J + 1
275                       KP = K + 1
276     
277                       IMJK = FUNIJK(IM,J,K)
278                       IJMK = FUNIJK(I,JM,K)
279                       IPJK = FUNIJK(IP,J,K)
280                       IJPK = FUNIJK(I,JP,K)
281                       IJKP = FUNIJK(I,J,KP)
282                       IJKM = FUNIJK(I,J,KM)
283     
284                       IPJMK = IP_OF(IJMK)
285     
286                       IJKE = EAST_OF(IJK)
287     !
288     ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
289                       U_se = Theta_Ue_bar(IJK)  * U_s(IJK,M_AM)  + Theta_Ue(IJK)  * U_s(IPJK,M_AM)
290                       Usw = Theta_Ue_bar(IMJK) * U_s(IMJK,M_AM) + Theta_Ue(IMJK) * U_s(IJK,M_AM)
291     
292                       Usn = Theta_Un_bar(IJK)  * U_s(IJK,M_AM)  + Theta_Un(IJK)  * U_s(IJPK,M_AM)
293                       Uss = Theta_Un_bar(IJMK) * U_s(IJMK,M_AM) + Theta_Un(IJMK) * U_s(IJK,M_AM)
294     
295                       Vsn =  Theta_U_ne(IJK)  * V_s(IJK,M_AM)  + Theta_U_nw(IJK)  * V_s(IPJK,M_AM)
296                       Vss =  Theta_U_ne(IJMK) * V_s(IJMK,M_AM) + Theta_U_nw(IJMK) * V_s(IPJMK,M_AM)
297     
298                       Vsc = (DELY_un(IJK) * Vss + DELY_us(IJK) * Vsn) / (DELY_un(IJK) + DELY_us(IJK))
299     
300                       IF(DO_K) THEN
301     
302                          IPJKM = IP_OF(IJKM)
303     
304                          Ust = Theta_Ut_bar(IJK)  * U_s(IJK,M_AM)  + Theta_Ut(IJK)  * U_s(IJKP,M_AM)
305                          Usb = Theta_Ut_bar(IJKM) * U_s(IJKM,M_AM) + Theta_Ut(IJKM) * U_s(IJK,M_AM)
306     
307                          Wst = Theta_U_te(IJK)  * W_s(IJK,M_AM)  + Theta_U_tw(IJK)  * W_s(IPJK,M_AM)
308                          Wsb = Theta_U_te(IJKM) * W_s(IJKM,M_AM) + Theta_U_tw(IJKM) * W_s(IPJKM,M_AM)
309                          Wsc = (DELZ_ut(IJK) * Wsb + DELZ_ub(IJK) * Wst) / (DELZ_ut(IJK) + DELZ_ub(IJK))
310     
311                          F_vir = F_vir +  Wsc* (Ust - Usb)*AXY(IJK)
312     
313                       ENDIF
314     !
315     ! adding convective terms (U dU/dx + V dU/dy + W dU/dz) to virtual mass
316                       F_vir = F_vir + U_s(IJK,M_AM)*(U_se - Usw)*AYZ(IJK) + &
317                                   Vsc * (Usn - Uss)*AXZ(IJK)
318     
319                       ROP_MA = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M_AM) + VOL(IPJK)*ROP_g(IJKE)*EP_s(IJKE,M_AM) )/(VOL(IJK) + VOL(IPJK))
320     
321                       F_vir = F_vir * Cv * ROP_MA
322     
323                       B_M(IJK,M) = B_M(IJK,M) - F_vir ! explicit part of virtual mass force
324     
325                    ENDIF
326     !
327     !!! END VIRTUAL MASS SECTION
328     
329                 ENDIF
330     
331              ENDIF
332     
333           END DO
334     
335     
336           RETURN
337           END SUBROUTINE CG_SOURCE_U_G
338     !
339     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
340     !                                                                      C
341     !  Module name: CG_SOURCE_U_g_BC(A_m, B_m, IER)                        C
342     !  Purpose: Determine contribution of cut-cell to source terms         C
343     !  for U_g momentum eq.                                                C
344     !                                                                      C
345     !                                                                      C
346     !  Author: Jeff Dietiker                              Date: 01-MAY-09  C
347     !  Reviewer:                                          Date:            C
348     !                                                                      C
349     !                                                                      C
350     !  Literature/Document References:                                     C
351     !                                                                      C
352     !  Variables referenced:                                               C
353     !  Variables modified:                                                 C
354     !                                                                      C
355     !  Local variables:                                                    C
356     !                                                                      C
357     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
358     !
359           SUBROUTINE CG_SOURCE_U_G_BC(A_M, B_M)
360     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
361     !...Switches: -xf
362     !
363     !  Include param.inc file to specify parameter values
364     !
365     !-----------------------------------------------
366     !   M o d u l e s
367     !-----------------------------------------------
368           USE param
369           USE param1
370           USE parallel
371           USE matrix
372           USE scales
373           USE constant
374           USE physprop
375           USE fldvar
376           USE visc_g
377           USE rxns
378           USE run
379           USE toleranc
380           USE geometry
381           USE indices
382           USE is
383           USE tau_g
384           USE bc
385           USE output
386           USE compar
387           USE fun_avg
388           USE functions
389     
390     !=======================================================================
391     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
392     !=======================================================================
393           USE cutcell
394           USE quadric
395     !=======================================================================
396     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
397     !=======================================================================
398     
399     
400           IMPLICIT NONE
401     !-----------------------------------------------
402     !   G l o b a l   P a r a m e t e r s
403     !-----------------------------------------------
404     !-----------------------------------------------
405     !   D u m m y   A r g u m e n t s
406     !-----------------------------------------------
407     !
408     !                      Indices
409           INTEGER          I,  J, K, IM, IJK,&
410                            JM, KM, IJKW, IMJK, IJMK,IJKM
411     !
412     !                      Solids phase
413           INTEGER          M
414     !
415     !                      Septadiagonal matrix A_m
416           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
417     !
418     !                      Vector b_m
419           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
420     !
421     !                      Turbulent shear stress
422           DOUBLE PRECISION  W_F_Slip
423     
424     !=======================================================================
425     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
426     !=======================================================================
427           INTEGER :: BCV
428           CHARACTER(LEN=9) :: BCT
429     
430           IF(CG_SAFE_MODE(3)==1) RETURN
431     
432           M = 0
433     
434           DO IJK = ijkstart3, ijkend3
435     
436              BCV = BC_U_ID(IJK)
437     
438              IF(BCV > 0 ) THEN
439                 BCT = BC_TYPE(BCV)
440              ELSE
441                 BCT = 'NONE'
442              ENDIF
443     
444     
445              SELECT CASE (BCT)
446     
447                 CASE ('CG_NSW')
448     
449                    IF(WALL_U_AT(IJK)) THEN
450     
451                       A_M(IJK,E,M) = ZERO
452                       A_M(IJK,W,M) = ZERO
453                       A_M(IJK,N,M) = ZERO
454                       A_M(IJK,S,M) = ZERO
455                       A_M(IJK,T,M) = ZERO
456                       A_M(IJK,B,M) = ZERO
457                       A_M(IJK,0,M) = -ONE
458     
459                       B_M(IJK,M) = ZERO
460     
461                       IF(K_EPSILON) THEN
462     
463                          I = I_OF(IJK)
464                          J = J_OF(IJK)
465                          K = K_OF(IJK)
466     
467                          IM = I - 1
468                          JM = J - 1
469                          KM = K - 1
470     
471                          IF(DABS(NORMAL_U(IJK,1))/=ONE) THEN
472     
473                             IF (U_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
474                                CALL Wall_Function(IJK,EAST_OF(IJK),ODX_E(I),W_F_Slip)
475                                A_M(IJK,E,M) = W_F_Slip
476                                A_M(IJK,0,M) = -ONE
477                             ELSEIF (U_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
478                                CALL Wall_Function(IJK,WEST_OF(IJK),ODX_E(IM),W_F_Slip)
479                                A_M(IJK,W,M) = W_F_Slip
480                                A_M(IJK,0,M) = -ONE
481                             ELSEIF (U_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
482                                CALL Wall_Function(IJK,NORTH_OF(IJK),ODY_N(J),W_F_Slip)
483                                A_M(IJK,N,M) = W_F_Slip
484                                A_M(IJK,0,M) = -ONE
485                             ELSEIF (U_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
486                                CALL Wall_Function(IJK,SOUTH_OF(IJK),ODY_N(JM),W_F_Slip)
487                                A_M(IJK,S,M) = W_F_Slip
488                                A_M(IJK,0,M) = -ONE
489                             ELSEIF (U_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
490                                A_M(IJK,T,M) = ONE
491                             ELSEIF (U_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
492                                A_M(IJK,B,M) = ONE
493                             ENDIF
494                          ENDIF
495     
496                       ENDIF
497     
498                    ENDIF
499     
500                 CASE ('CG_FSW')
501     
502                   IF(WALL_U_AT(IJK)) THEN
503     
504                       A_M(IJK,E,M) = ZERO
505                       A_M(IJK,W,M) = ZERO
506                       A_M(IJK,N,M) = ZERO
507                       A_M(IJK,S,M) = ZERO
508                       A_M(IJK,T,M) = ZERO
509                       A_M(IJK,B,M) = ZERO
510                       A_M(IJK,0,M) = -ONE
511     
512     !                  B_M(IJK,M) = - U_g(U_MASTER_OF(IJK))  ! Velocity of master node
513     
514                       B_M(IJK,M) = ZERO
515     
516                       IF(DABS(NORMAL_U(IJK,1))/=ONE) THEN
517     
518                          IF (U_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
519                             A_M(IJK,E,M) = ONE
520                          ELSEIF (U_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
521                             A_M(IJK,W,M) = ONE
522                          ELSEIF (U_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
523                             A_M(IJK,N,M) = ONE
524                          ELSEIF (U_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
525                             A_M(IJK,S,M) = ONE
526                          ELSEIF (U_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
527                             A_M(IJK,T,M) = ONE
528                          ELSEIF (U_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
529                             A_M(IJK,B,M) = ONE
530                          ENDIF
531     
532                       ENDIF
533     
534                    ENDIF
535     
536     
537                 CASE ('CG_PSW')
538     
539                    IF(WALL_U_AT(IJK)) THEN
540     
541                       A_M(IJK,E,M) = ZERO
542                       A_M(IJK,W,M) = ZERO
543                       A_M(IJK,N,M) = ZERO
544                       A_M(IJK,S,M) = ZERO
545                       A_M(IJK,T,M) = ZERO
546                       A_M(IJK,B,M) = ZERO
547                       A_M(IJK,0,M) = -ONE
548     
549     
550                       IF(BC_HW_G(BCV)==UNDEFINED) THEN   ! same as NSW
551                          B_M(IJK,M) = -BC_UW_G(BCV)
552                       ELSEIF(BC_HW_G(BCV)==ZERO) THEN   ! same as FSW
553                          B_M(IJK,M) = ZERO
554     
555                          IF(DABS(NORMAL_U(IJK,1))/=ONE) THEN
556     
557                             IF (U_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
558                                A_M(IJK,E,M) = ONE
559                             ELSEIF (U_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
560                                A_M(IJK,W,M) = ONE
561                             ELSEIF (U_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
562                                A_M(IJK,N,M) = ONE
563                             ELSEIF (U_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
564                                A_M(IJK,S,M) = ONE
565                             ELSEIF (U_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
566                                A_M(IJK,T,M) = ONE
567                             ELSEIF (U_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
568                                A_M(IJK,B,M) = ONE
569                             ENDIF
570     
571                          ENDIF
572     
573                       ELSE                              ! partial slip
574     
575                          B_M(IJK,M) = ZERO
576     
577                          I = I_OF(IJK)
578                          J = J_OF(IJK)
579                          K = K_OF(IJK)
580     
581                          IM = I - 1
582                          JM = J - 1
583                          KM = K - 1
584     
585                          IMJK = FUNIJK(IM,J,K)
586                          IJMK = FUNIJK(I,JM,K)
587                          IJKM = FUNIJK(I,J,KM)
588     
589     
590                          IF(DABS(NORMAL_U(IJK,1))/=ONE) THEN
591     
592                             IF (U_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
593                                A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ONEoDX_E_U(IJK))
594                                A_M(IJK,E,M) = -(HALF*BC_HW_G(BCV)-ONEoDX_E_U(IJK))
595                                B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
596     
597     !                           A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ODX_E(I))
598     !                           A_M(IJK,E,M) = -(HALF*BC_HW_G(BCV)-ODX_E(I))
599     !                           B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
600     !                           print*,'ug master at east'
601                             ELSEIF (U_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
602                                A_M(IJK,W,M) = -(HALF*BC_HW_G(BCV)-ONEoDX_E_U(IMJK))
603                                A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ONEoDX_E_U(IMJK))
604                                B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
605     
606     !                           A_M(IJK,W,M) = -(HALF*BC_HW_G(BCV)-ODX_E(IM))
607     !                           A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ODX_E(IM))
608     !                           B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
609     !                           print*,'ug master at west'
610                             ELSEIF (U_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
611     
612                                A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ONEoDY_N_U(IJK))
613                                A_M(IJK,N,M) = -(HALF*BC_HW_G(BCV)-ONEoDY_N_U(IJK))
614                                B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
615     
616     !                           A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ODY_N(J))
617     !                           A_M(IJK,N,M) = -(HALF*BC_HW_G(BCV)-ODY_N(J))
618     !                           B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
619                             ELSEIF (U_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
620                                A_M(IJK,S,M) = -(HALF*BC_HW_G(BCV)-ONEoDY_N_U(IJMK))
621                                A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ONEoDY_N_U(IJMK))
622                                B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
623     
624     !                           A_M(IJK,S,M) = -(HALF*BC_HW_G(BCV)-ODY_N(JM))
625     !                           A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ODY_N(JM))
626     !                           B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
627                             ELSEIF (U_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
628                                   A_M(IJK,0,M)=-(HALF*BC_HW_G(BCV)+ONEoDZ_T_U(IJK))
629                                   A_M(IJK,T,M)=-(HALF*BC_HW_G(BCV)-ONEoDZ_T_U(IJK))
630                                   B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
631     
632     !                              A_M(IJK,0,M)=-(HALF*BC_HW_G(L)+ODZ_T(K)*OX_E(I))
633     !                              A_M(IJK,T,M)=-(HALF*BC_HW_G(L)-ODZ_T(K)*OX_E(I))
634     !                              B_M(IJK,M) = -BC_HW_G(L)*BC_UW_G(L)
635                             ELSEIF (U_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
636                                   A_M(IJK,B,M) = -(HALF*BC_HW_G(BCV)-ONEoDZ_T_U(IJKM))
637                                   A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ONEoDZ_T_U(IJKM))
638                                   B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
639     
640     !                              A_M(IJK,B,M) = -(HALF*BC_HW_G(L)-ODZ_T(KM)*OX_E(I&
641     !                                 ))
642     !                              A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODZ_T(KM)*OX_E(I&
643     !                                 ))
644     !                              B_M(IJK,M) = -BC_HW_G(L)*BC_UW_G(L)
645                             ENDIF
646     
647                          ENDIF
648     
649                       ENDIF
650     
651                    ENDIF
652     
653                 CASE ('CG_MI')
654     
655                    A_M(IJK,E,M) = ZERO
656                    A_M(IJK,W,M) = ZERO
657                    A_M(IJK,N,M) = ZERO
658                    A_M(IJK,S,M) = ZERO
659                    A_M(IJK,T,M) = ZERO
660                    A_M(IJK,B,M) = ZERO
661                    A_M(IJK,0,M) = -ONE
662     
663                    IF(BC_U_g(BCV)/=UNDEFINED) THEN
664                       B_M(IJK,M) = - BC_U_g(BCV)
665                    ELSE
666                       B_M(IJK,M) = - BC_VELMAG_g(BCV)*NORMAL_U(IJK,1)
667                    ENDIF
668     
669     
670                    IJKW = WEST_OF(IJK)
671                    IF(FLUID_AT(IJKW)) THEN
672     
673                       A_M(IJKW,E,M) = ZERO
674                       A_M(IJKW,W,M) = ZERO
675                       A_M(IJKW,N,M) = ZERO
676                       A_M(IJKW,S,M) = ZERO
677                       A_M(IJKW,T,M) = ZERO
678                       A_M(IJKW,B,M) = ZERO
679                       A_M(IJKW,0,M) = -ONE
680     
681                       IF(BC_U_g(BCV)/=UNDEFINED) THEN
682                          B_M(IJKW,M) = - BC_U_g(BCV)
683                       ELSE
684                          B_M(IJKW,M) = - BC_VELMAG_g(BCV)*NORMAL_U(IJK,1)
685                       ENDIF
686     
687     
688                    ENDIF
689     
690                 CASE ('CG_PO')
691     
692                    A_M(IJK,E,M) = ZERO
693                    A_M(IJK,W,M) = ZERO
694                    A_M(IJK,N,M) = ZERO
695                    A_M(IJK,S,M) = ZERO
696                    A_M(IJK,T,M) = ZERO
697                    A_M(IJK,B,M) = ZERO
698                    A_M(IJK,0,M) = -ONE
699                    B_M(IJK,M) = ZERO
700     
701                    IJKW = WEST_OF(IJK)
702                    IF(FLUID_AT(IJKW)) THEN
703     
704                       A_M(IJK,W,M) = ONE
705                       A_M(IJK,0,M) = -ONE
706     
707                    ENDIF
708     
709              END SELECT
710     
711              BCV = BC_ID(IJK)
712     
713              IF(BCV > 0 ) THEN
714                 BCT = BC_TYPE(BCV)
715              ELSE
716                 BCT = 'NONE'
717              ENDIF
718     
719              SELECT CASE (BCT)
720     
721                 CASE ('CG_MI')
722     
723                    A_M(IJK,E,M) = ZERO
724                    A_M(IJK,W,M) = ZERO
725                    A_M(IJK,N,M) = ZERO
726                    A_M(IJK,S,M) = ZERO
727                    A_M(IJK,T,M) = ZERO
728                    A_M(IJK,B,M) = ZERO
729                    A_M(IJK,0,M) = -ONE
730     
731                    IF(BC_U_g(BCV)/=UNDEFINED) THEN
732                       B_M(IJK,M) = - BC_U_g(BCV)
733                    ELSE
734                       B_M(IJK,M) = - BC_VELMAG_g(BCV)*NORMAL_S(IJK,1)
735                    ENDIF
736     
737     
738                    IJKW = WEST_OF(IJK)
739                    IF(FLUID_AT(IJKW)) THEN
740     
741                       A_M(IJKW,E,M) = ZERO
742                       A_M(IJKW,W,M) = ZERO
743                       A_M(IJKW,N,M) = ZERO
744                       A_M(IJKW,S,M) = ZERO
745                       A_M(IJKW,T,M) = ZERO
746                       A_M(IJKW,B,M) = ZERO
747                       A_M(IJKW,0,M) = -ONE
748     
749                       IF(BC_U_g(BCV)/=UNDEFINED) THEN
750                          B_M(IJKW,M) = - BC_U_g(BCV)
751                       ELSE
752                          B_M(IJKW,M) = - BC_VELMAG_g(BCV)*NORMAL_S(IJK,1)
753                       ENDIF
754     
755     
756                    ENDIF
757     
758                 CASE ('CG_PO')
759     
760                    A_M(IJK,E,M) = ZERO
761                    A_M(IJK,W,M) = ZERO
762                    A_M(IJK,N,M) = ZERO
763                    A_M(IJK,S,M) = ZERO
764                    A_M(IJK,T,M) = ZERO
765                    A_M(IJK,B,M) = ZERO
766                    A_M(IJK,0,M) = -ONE
767                    B_M(IJK,M) = ZERO
768     
769     
770                    IJKW = WEST_OF(IJK)
771                    IF(FLUID_AT(IJKW)) THEN
772     
773                       A_M(IJK,W,M) = ONE
774                       A_M(IJK,0,M) = -ONE
775     
776                    ENDIF
777     
778              END SELECT
779     
780           ENDDO
781     
782     
783           RETURN
784     
785     !=======================================================================
786     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
787     !=======================================================================
788     
789           END SUBROUTINE CG_SOURCE_U_G_BC
790     
791     !// Comments on the modifications for DMP version implementation
792     !// 001 Include header file and common declarations for parallelization
793     !// 350 Changed do loop limits: 1,kmax2->kmin3,kmax3
794     !// 360 Check if i,j,k resides on current processor
795