File: N:\mfix\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     !                                                                      C
10     !                                                                      C
11     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
12           SUBROUTINE CG_SOURCE_U_G(A_M, B_M)
13     
14     ! Modules
15     !-----------------------------------------------
16           USE bc
17           USE compar
18           USE fldvar
19           USE fun_avg
20           use geometry, only: do_k, vol, vol_u
21           use geometry, only: axy, ayz, axz, ayz_u, axz_u, axy_u
22           USE indices
23           USE param
24           USE param1
25           USE physprop
26           USE run
27           USE toleranc
28           USE visc_g
29           USE cutcell
30           USE quadric
31           IMPLICIT NONE
32     
33     ! Dummy arguments
34     !-----------------------------------------------
35     ! Septadiagonal matrix A_m
36           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
37     ! Vector b_m
38           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
39     
40     ! Local variables
41     !-----------------------------------------------
42     ! Indices
43           INTEGER          I, IJK, IJKE, IPJK, IJKM, IPJKM
44     ! Phase index
45           INTEGER          M
46     ! Average volume fraction
47           DOUBLE PRECISION EPGA
48     !
49           INTEGER :: J,K,IM,JM,IP,JP,KM,KP
50           INTEGER :: IMJK,IJMK,IJPK,IJKC,IJKN
51           INTEGER :: IJKNE,IJKS,IJKSE,IPJMK,IJKP
52           INTEGER :: IJKT,IJKTE,IJKB,IJKBE
53           DOUBLE PRECISION :: Ue,Uw,Un,Us,Ut,Ub
54           DOUBLE PRECISION :: B_NOC
55           DOUBLE PRECISION :: MU_GT_E,MU_GT_W,MU_GT_N,MU_GT_S
56           DOUBLE PRECISION :: MU_GT_T,MU_GT_B,MU_GT_CUT
57           DOUBLE PRECISION :: UW_g
58           INTEGER :: BCV
59           INTEGER(4) :: BCT
60     
61     ! virtual (added) mass
62           DOUBLE PRECISION :: F_vir, ROP_MA
63           DOUBLE PRECISION :: U_se, Usw, Vsn, Vss, Vsc, Usn, Uss
64           DOUBLE PRECISION :: Wsb, Wst, Wsc, Usb, Ust
65     ! Wall function
66           DOUBLE PRECISION :: W_F_Slip
67     !-----------------------------------------------
68     
69           IF(CG_SAFE_MODE(3)==1) RETURN
70     
71           M = 0
72           IF (.NOT.MOMENTUM_X_EQ(0)) RETURN
73     !
74     !!!$omp    parallel do private(I, IJK, IJKE, IJKM, IPJK, IPJKM,     &
75     !!!$omp&                  ISV, Sdp, V0, Vpm, Vmt, Vbf,              &
76     !!!$omp&                  Vcf, EPMUGA, VTZA, WGE, PGE, ROGA,        &
77     !!!$omp&                  MUGA, ROPGA, EPGA )
78           DO IJK = ijkstart3, ijkend3
79              I = I_OF(IJK)
80              IJKE = EAST_OF(IJK)
81              IJKM = KM_OF(IJK)
82              IPJK = IP_OF(IJK)
83              IPJKM = IP_OF(IJKM)
84              EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
85              IF (IP_AT_E(IJK)) THEN
86     ! do nothing
87     
88     ! dilute flow
89              ELSE IF (EPGA <= DIL_EP_S) THEN
90     ! do nothing
91     
92              ELSEIF(INTERIOR_CELL_AT(IJK)) THEN
93                 BCV = BC_U_ID(IJK)
94     
95                 IF(BCV > 0 ) THEN
96                    BCT = BC_TYPE_ENUM(BCV)
97                 ELSE
98                    BCT = NONE
99                 ENDIF
100     
101                 SELECT CASE (BCT)
102                    CASE (CG_NSW)
103                       NOC_UG = .TRUE.
104                       UW_g = ZERO
105                       MU_GT_CUT =  (VOL(IJK)*MU_GT(IJK) + VOL(IPJK)*MU_GT(IJKE))/(VOL(IJK) + VOL(IPJK))
106     
107                       IF(.NOT.K_EPSILON) THEN
108                          A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_U_CUT(IJK)/DELH_U(IJK)
109                       ELSE
110                          CALL Wall_Function(IJK,IJK,ONE/DELH_U(IJK),W_F_Slip)
111                          A_M(IJK,0,M) = A_M(IJK,0,M)  - MU_GT_CUT * Area_U_CUT(IJK)*(ONE-W_F_Slip)/DELH_U(IJK)
112                       ENDIF
113                    CASE (CG_FSW)
114                       NOC_UG = .FALSE.
115                       UW_g = ZERO
116                    CASE(CG_PSW)
117                       IF(BC_HW_G(BCV)==UNDEFINED) THEN   ! same as NSW
118                          NOC_UG = .TRUE.
119                          UW_g = BC_UW_G(BCV)
120                          MU_GT_CUT =  (VOL(IJK)*MU_GT(IJK) + VOL(IPJK)*MU_GT(IJKE))/(VOL(IJK) + VOL(IPJK))
121                          A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_U_CUT(IJK)/DELH_U(IJK)
122                          B_M(IJK,M) = B_M(IJK,M) - MU_GT_CUT * UW_g * Area_U_CUT(IJK)/DELH_U(IJK)
123                       ELSEIF(BC_HW_G(BCV)==ZERO) THEN   ! same as FSW
124                          NOC_UG = .FALSE.
125                          UW_g = ZERO
126                       ELSE                              ! partial slip
127                          NOC_UG = .FALSE.
128                          UW_g = BC_UW_G(BCV)
129                          MU_GT_CUT =  (VOL(IJK)*MU_GT(IJK) + VOL(IPJK)*MU_GT(IJKE))/(VOL(IJK) + VOL(IPJK))
130                          A_M(IJK,0,M) = A_M(IJK,0,M) - MU_GT_CUT * Area_U_CUT(IJK)*(BC_HW_G(BCV))
131                          B_M(IJK,M) = B_M(IJK,M) - MU_GT_CUT * UW_g * Area_U_CUT(IJK)*(BC_HW_G(BCV))
132                       ENDIF
133                    CASE (NONE, CG_MI)
134                       NOC_UG = .FALSE.
135                 END SELECT
136     
137                 IF(NOC_UG) THEN
138     
139                    I = I_OF(IJK)
140                    J = J_OF(IJK)
141                    K = K_OF(IJK)
142                    IMJK = IM_OF(IJK)
143                    IJMK = JM_OF(IJK)
144                    IJKM = KM_OF(IJK)
145                    IPJK = IP_OF(IJK)
146                    IJPK = JP_OF(IJK)
147                    IJKP = KP_OF(IJK)
148     
149                    Ue = Theta_Ue_bar(IJK)  * U_g(IJK)  + Theta_Ue(IJK)  * U_g(IPJK)
150                    Uw = Theta_Ue_bar(IMJK) * U_g(IMJK) + Theta_Ue(IMJK) * U_g(IJK)
151                    Un = Theta_Un_bar(IJK)  * U_g(IJK)  + Theta_Un(IJK)  * U_g(IJPK)
152                    Us = Theta_Un_bar(IJMK) * U_g(IJMK) + Theta_Un(IJMK) * U_g(IJK)
153     
154                    IJKE = EAST_OF(IJK)
155                    IF (WALL_AT(IJK)) THEN
156                       IJKC = IJKE
157                    ELSE
158                       IJKC = IJK
159                    ENDIF
160                    IP = IP1(I)
161                    IJKN = NORTH_OF(IJK)
162                    IJKNE = EAST_OF(IJKN)
163                    JM = JM1(J)
164                    IPJMK = IP_OF(IJMK)
165                    IJKS = SOUTH_OF(IJK)
166                    IJKSE = EAST_OF(IJKS)
167                    KM = KM1(K)
168                    IJKT = TOP_OF(IJK)
169                    IJKTE = EAST_OF(IJKT)
170                    IJKB = BOTTOM_OF(IJK)
171                    IJKBE = EAST_OF(IJKB)
172     
173                    MU_GT_E = MU_GT(IJKE)
174                    MU_GT_W = MU_GT(IJKC)
175                    MU_GT_N = AVG_X_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),&
176                                      AVG_Y_H(MU_GT(IJKE),MU_GT(IJKNE),J),I)
177                    MU_GT_S = AVG_X_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),&
178                                      AVG_Y_H(MU_GT(IJKSE),MU_GT(IJKE),JM),I)
179     
180                    B_NOC =     MU_GT_E * Ayz_U(IJK)  * (Ue-UW_g) * NOC_U_E(IJK)  *2.0d0&
181                            -   MU_GT_W * Ayz_U(IMJK) * (Uw-UW_g) * NOC_U_E(IMJK) *2.0d0&
182                            +   MU_GT_N * Axz_U(IJK)  * (Un-UW_g) * NOC_U_N(IJK)  &
183                            -   MU_GT_S * Axz_U(IJMK) * (Us-UW_g) * NOC_U_N(IJMK)
184     
185                    IF(DO_K) THEN
186                       Ut = Theta_Ut_bar(IJK)  * U_g(IJK)  + Theta_Ut(IJK)  * U_g(IJKP)
187                       Ub = Theta_Ut_bar(IJKM) * U_g(IJKM) + Theta_Ut(IJKM) * U_g(IJK)
188     
189                       MU_GT_T = AVG_X_H(AVG_Z_H(MU_GT(IJKC),MU_GT(IJKT),K),&
190                                 AVG_Z_H(MU_GT(IJKE),MU_GT(IJKTE),K),I)
191     
192                       MU_GT_B = AVG_X_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJKC),KM),&
193                                 AVG_Z_H(MU_GT(IJKBE),MU_GT(IJKE),KM),I)
194     
195                       B_NOC = B_NOC  +   MU_GT_T * Axy_U(IJK)  * (Ut-UW_g) * NOC_U_T(IJK)  &
196                                      -   MU_GT_B * Axy_U(IJKM) * (Ub-UW_g) * NOC_U_T(IJKM)
197                    ENDIF
198     
199                    B_M(IJK,M) = B_M(IJK,M)   +  B_NOC
200     
201     
202                 ENDIF
203     
204                 IF(CUT_U_TREATMENT_AT(IJK)) THEN
205     
206     ! BEGIN VIRTUAL MASS SECTION (explicit terms)
207     ! adding transient term  dUs/dt to virtual mass term
208                    F_vir = ZERO
209                    IF(Added_Mass) THEN
210     
211                       F_vir = ( (U_s(IJK,M_AM) - U_sO(IJK,M_AM)) )*ODT*VOL_U(IJK)
212                       J = J_OF(IJK)
213                       K = K_OF(IJK)
214                       IM = I - 1
215                       JM = J - 1
216                       KM = K - 1
217                       IP = I + 1
218                       JP = J + 1
219                       KP = K + 1
220                       IMJK = FUNIJK(IM,J,K)
221                       IJMK = FUNIJK(I,JM,K)
222                       IPJK = FUNIJK(IP,J,K)
223                       IJPK = FUNIJK(I,JP,K)
224                       IJKP = FUNIJK(I,J,KP)
225                       IJKM = FUNIJK(I,J,KM)
226                       IPJMK = IP_OF(IJMK)
227                       IJKE = EAST_OF(IJK)
228     
229     ! defining gas-particles velocity at momentum cell faces (or scalar cell center)
230                       U_se = Theta_Ue_bar(IJK)  * U_s(IJK,M_AM)  + Theta_Ue(IJK)  * U_s(IPJK,M_AM)
231                       Usw = Theta_Ue_bar(IMJK) * U_s(IMJK,M_AM) + Theta_Ue(IMJK) * U_s(IJK,M_AM)
232                       Usn = Theta_Un_bar(IJK)  * U_s(IJK,M_AM)  + Theta_Un(IJK)  * U_s(IJPK,M_AM)
233                       Uss = Theta_Un_bar(IJMK) * U_s(IJMK,M_AM) + Theta_Un(IJMK) * U_s(IJK,M_AM)
234                       Vsn =  Theta_U_ne(IJK)  * V_s(IJK,M_AM)  + Theta_U_nw(IJK)  * V_s(IPJK,M_AM)
235                       Vss =  Theta_U_ne(IJMK) * V_s(IJMK,M_AM) + Theta_U_nw(IJMK) * V_s(IPJMK,M_AM)
236                       Vsc = (DELY_un(IJK) * Vss + DELY_us(IJK) * Vsn) / (DELY_un(IJK) + DELY_us(IJK))
237     
238                       IF(DO_K) THEN
239                          IPJKM = IP_OF(IJKM)
240                          Ust = Theta_Ut_bar(IJK)  * U_s(IJK,M_AM)  + Theta_Ut(IJK)  * U_s(IJKP,M_AM)
241                          Usb = Theta_Ut_bar(IJKM) * U_s(IJKM,M_AM) + Theta_Ut(IJKM) * U_s(IJK,M_AM)
242                          Wst = Theta_U_te(IJK)  * W_s(IJK,M_AM)  + Theta_U_tw(IJK)  * W_s(IPJK,M_AM)
243                          Wsb = Theta_U_te(IJKM) * W_s(IJKM,M_AM) + Theta_U_tw(IJKM) * W_s(IPJKM,M_AM)
244                          Wsc = (DELZ_ut(IJK) * Wsb + DELZ_ub(IJK) * Wst) / (DELZ_ut(IJK) + DELZ_ub(IJK))
245                          F_vir = F_vir +  Wsc* (Ust - Usb)*AXY(IJK)
246                       ENDIF
247     
248     ! adding convective terms (U dU/dx + V dU/dy + W dU/dz) to virtual mass
249                       F_vir = F_vir + U_s(IJK,M_AM)*(U_se - Usw)*AYZ(IJK) + &
250                                   Vsc * (Usn - Uss)*AXZ(IJK)
251     
252                       ROP_MA = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M_AM) + &
253                                 VOL(IPJK)*ROP_g(IJKE)*EP_s(IJKE,M_AM) )/(VOL(IJK) + VOL(IPJK))
254     
255                       F_vir = F_vir * Cv * ROP_MA
256     
257                       B_M(IJK,M) = B_M(IJK,M) - F_vir ! explicit part of virtual mass force
258     
259                    ENDIF
260     ! END VIRTUAL MASS SECTION
261     
262                 ENDIF
263              ENDIF
264     
265           END DO
266     
267     
268           RETURN
269     
270         CONTAINS
271     
272           INCLUDE 'functions.inc'
273     
274           END SUBROUTINE CG_SOURCE_U_G
275     
276     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
277     !                                                                      C
278     !  Module name: CG_SOURCE_U_g_BC(A_m, B_m, IER)                        C
279     !  Purpose: Determine contribution of cut-cell to source terms         C
280     !  for U_g momentum eq.                                                C
281     !                                                                      C
282     !                                                                      C
283     !  Author: Jeff Dietiker                              Date: 01-MAY-09  C
284     !                                                                      C
285     !                                                                      C
286     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
287     
288           SUBROUTINE CG_SOURCE_U_G_BC(A_M, B_M)
289     
290     ! Modules
291     !-----------------------------------------------
292           USE bc
293           USE compar
294           USE fldvar
295           USE fun_avg
296           use geometry, only: odx_e, ody_n
297           USE indices
298           USE param
299           USE param1
300           USE physprop
301           USE run
302           USE toleranc
303           USE visc_g
304           USE cutcell
305           USE quadric
306           IMPLICIT NONE
307     
308     ! Dummy arguments
309     !-----------------------------------------------
310     ! Septadiagonal matrix A_m
311           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
312     ! Vector b_m
313           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
314     
315     ! Local variables
316     !-----------------------------------------------
317     ! Indices
318           INTEGER :: I,  J, K, IM, IJK
319           INTEGER :: JM, KM, IJKW, IMJK, IJMK,IJKM
320     ! Solids phase
321           INTEGER :: M
322     ! Turbulent shear stress
323           DOUBLE PRECISION :: W_F_Slip
324           INTEGER :: BCV
325           INTEGER :: BCT
326     !-----------------------------------------------
327     
328           IF(CG_SAFE_MODE(3)==1) RETURN
329     
330           M = 0
331     
332           DO IJK = ijkstart3, ijkend3
333     
334              BCV = BC_U_ID(IJK)
335     
336              IF(BCV > 0 ) THEN
337                 BCT = BC_TYPE_ENUM(BCV)
338              ELSE
339                 BCT = NONE
340              ENDIF
341     
342     
343              SELECT CASE (BCT)
344     
345                 CASE (CG_NSW)
346     
347                    IF(WALL_U_AT(IJK)) THEN
348     
349                       A_M(IJK,east,M) = ZERO
350                       A_M(IJK,west,M) = ZERO
351                       A_M(IJK,north,M) = ZERO
352                       A_M(IJK,south,M) = ZERO
353                       A_M(IJK,top,M) = ZERO
354                       A_M(IJK,bottom,M) = ZERO
355                       A_M(IJK,0,M) = -ONE
356     
357                       B_M(IJK,M) = ZERO
358     
359                       IF(K_EPSILON) THEN
360     
361                          I = I_OF(IJK)
362                          J = J_OF(IJK)
363                          K = K_OF(IJK)
364     
365                          IM = I - 1
366                          JM = J - 1
367                          KM = K - 1
368     
369                          IF(DABS(NORMAL_U(IJK,1))/=ONE) THEN
370     
371                             IF (U_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
372                                CALL Wall_Function(IJK,EAST_OF(IJK),ODX_E(I),W_F_Slip)
373                                A_M(IJK,east,M) = W_F_Slip
374                                A_M(IJK,0,M) = -ONE
375                             ELSEIF (U_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
376                                CALL Wall_Function(IJK,WEST_OF(IJK),ODX_E(IM),W_F_Slip)
377                                A_M(IJK,west,M) = W_F_Slip
378                                A_M(IJK,0,M) = -ONE
379                             ELSEIF (U_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
380                                CALL Wall_Function(IJK,NORTH_OF(IJK),ODY_N(J),W_F_Slip)
381                                A_M(IJK,north,M) = W_F_Slip
382                                A_M(IJK,0,M) = -ONE
383                             ELSEIF (U_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
384                                CALL Wall_Function(IJK,SOUTH_OF(IJK),ODY_N(JM),W_F_Slip)
385                                A_M(IJK,south,M) = W_F_Slip
386                                A_M(IJK,0,M) = -ONE
387                             ELSEIF (U_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
388                                A_M(IJK,top,M) = ONE
389                             ELSEIF (U_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
390                                A_M(IJK,bottom,M) = ONE
391                             ENDIF
392                          ENDIF
393     
394                       ENDIF
395     
396                    ENDIF
397     
398                 CASE (CG_FSW)
399     
400                   IF(WALL_U_AT(IJK)) THEN
401     
402                       A_M(IJK,east,M) = ZERO
403                       A_M(IJK,west,M) = ZERO
404                       A_M(IJK,north,M) = ZERO
405                       A_M(IJK,south,M) = ZERO
406                       A_M(IJK,top,M) = ZERO
407                       A_M(IJK,bottom,M) = ZERO
408                       A_M(IJK,0,M) = -ONE
409     
410     !                  B_M(IJK,M) = - U_g(U_MASTER_OF(IJK))  ! Velocity of master node
411     
412                       B_M(IJK,M) = ZERO
413     
414                       IF(DABS(NORMAL_U(IJK,1))/=ONE) THEN
415     
416                          IF (U_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
417                             A_M(IJK,east,M) = ONE
418                          ELSEIF (U_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
419                             A_M(IJK,west,M) = ONE
420                          ELSEIF (U_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
421                             A_M(IJK,north,M) = ONE
422                          ELSEIF (U_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
423                             A_M(IJK,south,M) = ONE
424                          ELSEIF (U_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
425                             A_M(IJK,top,M) = ONE
426                          ELSEIF (U_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
427                             A_M(IJK,bottom,M) = ONE
428                          ENDIF
429     
430                       ENDIF
431     
432                    ENDIF
433     
434     
435                 CASE (CG_PSW)
436     
437                    IF(WALL_U_AT(IJK)) THEN
438     
439                       A_M(IJK,east,M) = ZERO
440                       A_M(IJK,west,M) = ZERO
441                       A_M(IJK,north,M) = ZERO
442                       A_M(IJK,south,M) = ZERO
443                       A_M(IJK,top,M) = ZERO
444                       A_M(IJK,bottom,M) = ZERO
445                       A_M(IJK,0,M) = -ONE
446     
447     
448                       IF(BC_HW_G(BCV)==UNDEFINED) THEN   ! same as NSW
449                          B_M(IJK,M) = -BC_UW_G(BCV)
450                       ELSEIF(BC_HW_G(BCV)==ZERO) THEN   ! same as FSW
451                          B_M(IJK,M) = ZERO
452     
453                          IF(DABS(NORMAL_U(IJK,1))/=ONE) THEN
454     
455                             IF (U_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
456                                A_M(IJK,east,M) = ONE
457                             ELSEIF (U_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
458                                A_M(IJK,west,M) = ONE
459                             ELSEIF (U_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
460                                A_M(IJK,north,M) = ONE
461                             ELSEIF (U_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
462                                A_M(IJK,south,M) = ONE
463                             ELSEIF (U_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
464                                A_M(IJK,top,M) = ONE
465                             ELSEIF (U_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
466                                A_M(IJK,bottom,M) = ONE
467                             ENDIF
468     
469                          ENDIF
470     
471                       ELSE                              ! partial slip
472     
473                          B_M(IJK,M) = ZERO
474     
475                          I = I_OF(IJK)
476                          J = J_OF(IJK)
477                          K = K_OF(IJK)
478     
479                          IM = I - 1
480                          JM = J - 1
481                          KM = K - 1
482     
483                          IMJK = FUNIJK(IM,J,K)
484                          IJMK = FUNIJK(I,JM,K)
485                          IJKM = FUNIJK(I,J,KM)
486     
487     
488                          IF(DABS(NORMAL_U(IJK,1))/=ONE) THEN
489     
490                             IF (U_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
491                                A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ONEoDX_E_U(IJK))
492                                A_M(IJK,east,M) = -(HALF*BC_HW_G(BCV)-ONEoDX_E_U(IJK))
493                                B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
494     
495     !                           A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ODX_E(I))
496     !                           A_M(IJK,east,M) = -(HALF*BC_HW_G(BCV)-ODX_E(I))
497     !                           B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
498     !                           print*,'ug master at east'
499                             ELSEIF (U_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
500                                A_M(IJK,west,M) = -(HALF*BC_HW_G(BCV)-ONEoDX_E_U(IMJK))
501                                A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ONEoDX_E_U(IMJK))
502                                B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
503     
504     !                           A_M(IJK,west,M) = -(HALF*BC_HW_G(BCV)-ODX_E(IM))
505     !                           A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ODX_E(IM))
506     !                           B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
507     !                           print*,'ug master at west'
508                             ELSEIF (U_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
509     
510                                A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ONEoDY_N_U(IJK))
511                                A_M(IJK,north,M) = -(HALF*BC_HW_G(BCV)-ONEoDY_N_U(IJK))
512                                B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
513     
514     !                           A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ODY_N(J))
515     !                           A_M(IJK,north,M) = -(HALF*BC_HW_G(BCV)-ODY_N(J))
516     !                           B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
517                             ELSEIF (U_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
518                                A_M(IJK,south,M) = -(HALF*BC_HW_G(BCV)-ONEoDY_N_U(IJMK))
519                                A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ONEoDY_N_U(IJMK))
520                                B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
521     
522     !                           A_M(IJK,south,M) = -(HALF*BC_HW_G(BCV)-ODY_N(JM))
523     !                           A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ODY_N(JM))
524     !                           B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
525                             ELSEIF (U_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
526                                   A_M(IJK,0,M)=-(HALF*BC_HW_G(BCV)+ONEoDZ_T_U(IJK))
527                                   A_M(IJK,top,M)=-(HALF*BC_HW_G(BCV)-ONEoDZ_T_U(IJK))
528                                   B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
529     
530     !                              A_M(IJK,0,M)=-(HALF*BC_HW_G(L)+ODZ_T(K)*OX_E(I))
531     !                              A_M(IJK,top,M)=-(HALF*BC_HW_G(L)-ODZ_T(K)*OX_E(I))
532     !                              B_M(IJK,M) = -BC_HW_G(L)*BC_UW_G(L)
533                             ELSEIF (U_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
534                                   A_M(IJK,bottom,M) = -(HALF*BC_HW_G(BCV)-ONEoDZ_T_U(IJKM))
535                                   A_M(IJK,0,M) = -(HALF*BC_HW_G(BCV)+ONEoDZ_T_U(IJKM))
536                                   B_M(IJK,M) = -BC_HW_G(BCV)*BC_UW_G(BCV)
537     
538     !                              A_M(IJK,bottom,M) = -(HALF*BC_HW_G(L)-ODZ_T(KM)*OX_E(I&
539     !                                 ))
540     !                              A_M(IJK,0,M) = -(HALF*BC_HW_G(L)+ODZ_T(KM)*OX_E(I&
541     !                                 ))
542     !                              B_M(IJK,M) = -BC_HW_G(L)*BC_UW_G(L)
543                             ENDIF
544     
545                          ENDIF
546     
547                       ENDIF
548     
549                    ENDIF
550     
551                 CASE (CG_MI)
552     
553                    A_M(IJK,east,M) = ZERO
554                    A_M(IJK,west,M) = ZERO
555                    A_M(IJK,north,M) = ZERO
556                    A_M(IJK,south,M) = ZERO
557                    A_M(IJK,top,M) = ZERO
558                    A_M(IJK,bottom,M) = ZERO
559                    A_M(IJK,0,M) = -ONE
560     
561                    IF(BC_U_g(BCV)/=UNDEFINED) THEN
562                       B_M(IJK,M) = - BC_U_g(BCV)
563                    ELSE
564                       B_M(IJK,M) = - BC_VELMAG_g(BCV)*NORMAL_U(IJK,1)
565                    ENDIF
566     
567     
568                    IJKW = WEST_OF(IJK)
569                    IF(FLUID_AT(IJKW)) THEN
570     
571                       A_M(IJKW,east,M) = ZERO
572                       A_M(IJKW,west,M) = ZERO
573                       A_M(IJKW,north,M) = ZERO
574                       A_M(IJKW,south,M) = ZERO
575                       A_M(IJKW,top,M) = ZERO
576                       A_M(IJKW,bottom,M) = ZERO
577                       A_M(IJKW,0,M) = -ONE
578     
579                       IF(BC_U_g(BCV)/=UNDEFINED) THEN
580                          B_M(IJKW,M) = - BC_U_g(BCV)
581                       ELSE
582                          B_M(IJKW,M) = - BC_VELMAG_g(BCV)*NORMAL_U(IJK,1)
583                       ENDIF
584     
585     
586                    ENDIF
587     
588                 CASE (CG_PO)
589     
590                    A_M(IJK,east,M) = ZERO
591                    A_M(IJK,west,M) = ZERO
592                    A_M(IJK,north,M) = ZERO
593                    A_M(IJK,south,M) = ZERO
594                    A_M(IJK,top,M) = ZERO
595                    A_M(IJK,bottom,M) = ZERO
596                    A_M(IJK,0,M) = -ONE
597                    B_M(IJK,M) = ZERO
598     
599                    IJKW = WEST_OF(IJK)
600                    IF(FLUID_AT(IJKW)) THEN
601     
602                       A_M(IJK,west,M) = ONE
603                       A_M(IJK,0,M) = -ONE
604     
605                    ENDIF
606     
607              END SELECT
608     
609              BCV = BC_ID(IJK)
610     
611              IF(BCV > 0 ) THEN
612                 BCT = BC_TYPE_ENUM(BCV)
613              ELSE
614                 BCT = NONE
615              ENDIF
616     
617              SELECT CASE (BCT)
618     
619                 CASE (CG_MI)
620     
621                    A_M(IJK,east,M) = ZERO
622                    A_M(IJK,west,M) = ZERO
623                    A_M(IJK,north,M) = ZERO
624                    A_M(IJK,south,M) = ZERO
625                    A_M(IJK,top,M) = ZERO
626                    A_M(IJK,bottom,M) = ZERO
627                    A_M(IJK,0,M) = -ONE
628     
629                    IF(BC_U_g(BCV)/=UNDEFINED) THEN
630                       B_M(IJK,M) = - BC_U_g(BCV)
631                    ELSE
632                       B_M(IJK,M) = - BC_VELMAG_g(BCV)*NORMAL_S(IJK,1)
633                    ENDIF
634     
635     
636                    IJKW = WEST_OF(IJK)
637                    IF(FLUID_AT(IJKW)) THEN
638     
639                       A_M(IJKW,east,M) = ZERO
640                       A_M(IJKW,west,M) = ZERO
641                       A_M(IJKW,north,M) = ZERO
642                       A_M(IJKW,south,M) = ZERO
643                       A_M(IJKW,top,M) = ZERO
644                       A_M(IJKW,bottom,M) = ZERO
645                       A_M(IJKW,0,M) = -ONE
646     
647                       IF(BC_U_g(BCV)/=UNDEFINED) THEN
648                          B_M(IJKW,M) = - BC_U_g(BCV)
649                       ELSE
650                          B_M(IJKW,M) = - BC_VELMAG_g(BCV)*NORMAL_S(IJK,1)
651                       ENDIF
652     
653     
654                    ENDIF
655     
656                 CASE (CG_PO)
657     
658                    A_M(IJK,east,M) = ZERO
659                    A_M(IJK,west,M) = ZERO
660                    A_M(IJK,north,M) = ZERO
661                    A_M(IJK,south,M) = ZERO
662                    A_M(IJK,top,M) = ZERO
663                    A_M(IJK,bottom,M) = ZERO
664                    A_M(IJK,0,M) = -ONE
665                    B_M(IJK,M) = ZERO
666     
667     
668                    IJKW = WEST_OF(IJK)
669                    IF(FLUID_AT(IJKW)) THEN
670     
671                       A_M(IJK,west,M) = ONE
672                       A_M(IJK,0,M) = -ONE
673     
674                    ENDIF
675     
676              END SELECT
677     
678           ENDDO
679     
680     
681           RETURN
682     
683     
684         CONTAINS
685     
686           INCLUDE 'functions.inc'
687     
688           END SUBROUTINE CG_SOURCE_U_G_BC
689