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

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