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