File: /nfs/home/0/users/jenkins/mfix.git/model/cartesian_grid/CG_source_w_s.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 SUBROUTINE CG_SOURCE_W_S(A_M, B_M, M)
16
17
18
19
20 USE param, only: dimension_3, dimension_m
21 USE param1, only: zero, undefined
22
23
24 USE physprop, only: smax, mmax
25
26 USE physprop, only: cv
27 USE visc_s, only: mu_s
28 USE fldvar
29 USE run
30
31 USE toleranc, only: dil_ep_s
32 USE bc
33
34 USE geometry
35 USE indices
36 USE compar
37
38 USE cutcell
39 USE quadric
40 USE fun_avg
41 USE functions
42
43 IMPLICIT NONE
44
45
46
47
48
49 INTEGER, INTENT(IN) :: M
50
51 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
52
53 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
54
55
56
57
58 INTEGER :: I, J, K, IJK, IJKT, IMJK, IJMK, IJKM, IJKP, IMJKP
59 INTEGER :: IJKE, IJKW, IJKTE, IM, IPJK
60
61 INTEGER :: L
62
63 DOUBLE PRECISION :: EPSA, EPStmp
64
65 DOUBLE PRECISION :: F_vir, ROP_MA
66 DOUBLE PRECISION :: Uge, Ugw, Wge, Wgw, Wgn, &
67 Wgs, Wgt, Wgb, Ugc, Vgc, Vgn, Vgs
68
69 DOUBLE PRECISION :: F_2
70
71 INTEGER :: JM,IP,JP,IJPK,IJKC,IJKN,IJKNE,IJKS,IJKSE,IPJMK,KM,KP,IJMKP
72 INTEGER :: IJKTN,IJKWT,IJKST
73 DOUBLE PRECISION :: We,Ww,Wn,Ws,Wt,Wb
74 DOUBLE PRECISION :: B_NOC
75 DOUBLE PRECISION :: MU_S_E,MU_S_W,MU_S_N,MU_S_S,MU_S_T,MU_S_B,MU_S_CUT
76 DOUBLE PRECISION :: WW_s
77 INTEGER :: BCV
78 CHARACTER(LEN=9) :: BCT
79
80
81 IF(CG_SAFE_MODE(5)==1) RETURN
82
83
84 IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
85 (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
86
87 IF (MOMENTUM_Z_EQ(M)) THEN
88
89
90
91
92
93
94 DO IJK = ijkstart3, ijkend3
95
96
97 IF(WALL_AT(IJK)) cycle
98
99 I = I_OF(IJK)
100 J = J_OF(IJK)
101 K = K_OF(IJK)
102 IMJK = IM_OF(IJK)
103 IJMK = JM_OF(IJK)
104 IJKM = KM_OF(IJK)
105 IJKT = TOP_OF(IJK)
106
107 IF (KT_TYPE_ENUM == GHD_2007) THEN
108 = ZERO
109 DO L = 1, SMAX
110 EPStmp = EPStmp + AVG_Z(EP_S(IJK,L),EP_S(IJKT,L),K)
111 ENDDO
112 EPSA = EPStmp
113 ELSE
114 EPSA = AVG_Z(EP_S(IJK,M),EP_S(IJKT,M),K)
115 ENDIF
116
117
118 IF (IP_AT_T(IJK)) THEN
119
120
121 ELSEIF (SIP_AT_T(IJK)) THEN
122
123
124 ELSEIF (EPSA <= DIL_EP_S) THEN
125
126
127 ELSEIF(INTERIOR_CELL_AT(IJK)) THEN
128 BCV = BC_W_ID(IJK)
129 IF(BCV > 0 ) THEN
130 BCT = BC_TYPE(BCV)
131 ELSE
132 BCT = 'NONE'
133 ENDIF
134
135 SELECT CASE (BCT)
136 CASE ('CG_NSW')
137 NOC_WS = .TRUE.
138 WW_s = ZERO
139 MU_S_CUT = (VOL(IJK)*MU_S(IJK,M) + VOL(IJKT)*MU_S(IJKT,M))/(VOL(IJK) + VOL(IJKT))
140 A_M(IJK,0,M) = A_M(IJK,0,M) - MU_S_CUT * Area_W_CUT(IJK)/DELH_W(IJK)
141 CASE ('CG_FSW')
142 NOC_WS = .FALSE.
143 WW_s = ZERO
144 CASE('CG_PSW')
145 IF(BC_JJ_PS(BCV)==1) THEN
146 = .FALSE.
147 WW_s = BC_WW_S(BCV,M)
148 CALL CG_CALC_GRBDRY(IJK, 'W_MOMENTUM', M, BCV, F_2)
149 A_M(IJK,0,M) = A_M(IJK,0,M) - Area_W_CUT(IJK)*F_2
150 B_M(IJK,M) = B_M(IJK,M) - Area_W_CUT(IJK)*F_2*WW_s
151 ELSEIF(BC_HW_S(BCV,M)==UNDEFINED) THEN
152 = .TRUE.
153 WW_s = BC_WW_S(BCV,M)
154 MU_S_CUT = (VOL(IJK)*MU_S(IJK,M) + VOL(IJKT)*MU_S(IJKT,M))/(VOL(IJK) + VOL(IJKT))
155 A_M(IJK,0,M) = A_M(IJK,0,M) - MU_S_CUT * Area_W_CUT(IJK)/DELH_W(IJK)
156 B_M(IJK,M) = B_M(IJK,M) - MU_S_CUT * WW_s * Area_W_CUT(IJK)/DELH_W(IJK)
157 ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN
158 = .FALSE.
159 WW_s = ZERO
160 ELSE
161 = .FALSE.
162 WW_s = BC_WW_S(BCV,M)
163 MU_S_CUT = (VOL(IJK)*MU_S(IJK,M) + VOL(IJKT)*MU_S(IJKT,M))/(VOL(IJK) + VOL(IJKT))
164 A_M(IJK,0,M) = A_M(IJK,0,M) - MU_S_CUT * Area_W_CUT(IJK)*(BC_HW_S(BCV,M))
165 B_M(IJK,M) = B_M(IJK,M) - MU_S_CUT * WW_s * Area_W_CUT(IJK)*(BC_HW_S(BCV,M))
166 ENDIF
167 CASE ('NONE')
168 NOC_WS = .FALSE.
169
170 END SELECT
171
172
173 IF(NOC_WS) THEN
174 IMJK = IM_OF(IJK)
175 IJMK = JM_OF(IJK)
176 IJKM = KM_OF(IJK)
177 IPJK = IP_OF(IJK)
178 IJPK = JP_OF(IJK)
179 IJKP = KP_OF(IJK)
180
181 We = Theta_We_bar(IJK) * W_S(IJK,M) + Theta_We(IJK) * W_S(IPJK,M)
182 Ww = Theta_We_bar(IMJK) * W_S(IMJK,M) + Theta_We(IMJK) * W_S(IJK,M)
183 Wn = Theta_Wn_bar(IJK) * W_S(IJK,M) + Theta_Wn(IJK) * W_S(IJPK,M)
184 Ws = Theta_Wn_bar(IJMK) * W_S(IJMK,M) + Theta_Wn(IJMK) * W_S(IJK,M)
185 Wt = Theta_Wt_bar(IJK) * W_S(IJK,M) + Theta_Wt(IJK) * W_S(IJKP,M)
186 Wb = Theta_Wt_bar(IJKM) * W_S(IJKM,M) + Theta_Wt(IJKM) * W_S(IJK,M)
187
188 IF (WALL_AT(IJK)) THEN
189 IJKC = IJKT
190 ELSE
191 IJKC = IJK
192 ENDIF
193
194 IJKE = EAST_OF(IJK)
195 ijkt = top_of(ijk)
196 IP = IP1(I)
197 IM = IM1(I)
198 IJKN = NORTH_OF(IJK)
199 IJKNE = EAST_OF(IJKN)
200 JM = JM1(J)
201 IPJMK = IP_OF(IJMK)
202 IJKS = SOUTH_OF(IJK)
203 IJKSE = EAST_OF(IJKS)
204 KP = KP1(K)
205 IJKT = TOP_OF(IJK)
206 IJKE = EAST_OF(IJK)
207 IJKP = KP_OF(IJK)
208 IJKTN = NORTH_OF(IJKT)
209 IJKTE = EAST_OF(IJKT)
210 IJKW = WEST_OF(IJK)
211 IJKWT = TOP_OF(IJKW)
212 IJKS = SOUTH_OF(IJK)
213 IJKST = TOP_OF(IJKS)
214
215 MU_S_E = AVG_Z_H(AVG_X_H(MU_S(IJKC,M),MU_S(IJKE,M),I),&
216 AVG_X_H(MU_S(IJKT,M),MU_S(IJKTE,M),I),K)
217 MU_S_W = AVG_Z_H(AVG_X_H(MU_S(IJKW,M),MU_S(IJKC,M),IM),&
218 AVG_X_H(MU_S(IJKWT,M),MU_S(IJKT,M),IM),K)
219 MU_S_N = AVG_Z_H(AVG_Y_H(MU_S(IJKC,M),MU_S(IJKN,M),J),&
220 AVG_Y_H(MU_S(IJKT,M),MU_S(IJKTN,M),J),K)
221 MU_S_S = AVG_Z_H(AVG_Y_H(MU_S(IJKS,M),MU_S(IJKC,M),JM),&
222 AVG_Y_H(MU_S(IJKST,M),MU_S(IJKT,M),JM),K)
223 MU_S_T = MU_S(IJKT,M)
224 MU_S_B = MU_S(IJKC,M)
225
226 B_NOC = MU_S_E * Ayz_W(IJK) * (We-WW_s) * NOC_W_E(IJK) &
227 - MU_S_W * Ayz_W(IMJK) * (Ww-WW_s) * NOC_W_E(IMJK) &
228 + MU_S_N * Axz_W(IJK) * (Wn-WW_s) * NOC_W_N(IJK) &
229 - MU_S_S * Axz_W(IJMK) * (Ws-WW_s) * NOC_W_N(IJMK) &
230 + MU_S_T * Axy_W(IJK) * (Wt-WW_s) * NOC_W_T(IJK) *2.0d0&
231 - MU_S_B * Axy_W(IJKM) * (Wb-WW_s) * NOC_W_T(IJKM) *2.0D0
232
233 B_M(IJK,M) = B_M(IJK,M) + B_NOC
234 ENDIF
235
236
237 IF(CUT_W_TREATMENT_AT(IJK)) THEN
238
239
240 = ZERO
241 IF(Added_Mass.AND. M == M_AM ) THEN
242 F_vir = ( (W_g(IJK) - W_gO(IJK)) )*ODT*VOL_W(IJK)
243 I = I_OF(IJK)
244 J = J_OF(IJK)
245 K = K_OF(IJK)
246 IM = I - 1
247 JM = J - 1
248 KM = K - 1
249 IP = I + 1
250 JP = J + 1
251 KP = K + 1
252 IMJK = FUNIJK(IM,J,K)
253 IJMK = FUNIJK(I,JM,K)
254 IPJK = FUNIJK(IP,J,K)
255 IJPK = FUNIJK(I,JP,K)
256 IJKP = FUNIJK(I,J,KP)
257 IJKM = FUNIJK(I,J,KM)
258 IMJKP = KP_OF(IMJK)
259 IJMKP = KP_OF(IJMK)
260 IJKE = EAST_OF(IJK)
261
262
263 = Theta_We_bar(IJK) * W_g(IJK) + Theta_We(IJK) * W_g(IPJK)
264 Wgw = Theta_We_bar(IMJK) * W_g(IMJK) + Theta_We(IMJK) * W_g(IJK)
265 Uge = Theta_W_te(IJK) * U_g(IJK) + Theta_W_be(IJK) * U_g(IJKP)
266 Ugw = Theta_W_te(IMJK) * U_g(IMJK) + Theta_W_be(IMJK) * U_g(IMJKP)
267 Ugc = (DELX_we(IJK) * Ugw + DELX_ww(IJK) * Uge) / (DELX_we(IJK) + DELX_ww(IJK))
268 Wgn = Theta_Wn_bar(IJK) * W_g(IJK) + Theta_Wn(IJK) * W_g(IJPK)
269 Wgs = Theta_Wn_bar(IJMK) * W_g(IJMK) + Theta_Wn(IJMK) * W_g(IJK)
270 Vgn = Theta_W_tn(IJK) * V_g(IJK) + Theta_W_bn(IJK) * V_g(IJKP)
271 Vgs = Theta_W_tn(IJMK) * V_g(IJMK) + Theta_W_bn(IJMK) * V_g(IJMKP)
272 Vgc = (DELY_wn(IJK) * Vgs + DELY_ws(IJK) * Vgn) / (DELY_wn(IJK) + DELY_ws(IJK))
273 Wgt = Theta_Wt_bar(IJK) * W_g(IJK) + Theta_Wt(IJK) * W_g(IJKP)
274 Wgb = Theta_Wt_bar(IMJK) * W_g(IMJK) + Theta_Wt(IMJK) * W_g(IMJKP)
275
276
277 = F_vir + Ugc * (Wge - Wgw)*AYZ(IJK) + &
278 Vgc * (Wgn - Wgs)*AXZ(IJK) + &
279 W_g(IJK)*(Wgt - Wgb)*AXY(IJK)
280 ROP_MA = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M) + &
281 VOL(IJKT)*ROP_g(IJKT)*EP_s(IJKT,M))/&
282 (VOL(IJK) + VOL(IJKT))
283 F_vir = F_vir * Cv * ROP_MA
284 B_M(IJK,M) = B_M(IJK,M) - F_vir
285 ENDIF
286 ENDIF
287
288
289 ENDIF
290 ENDDO
291 ENDIF
292 ENDIF
293
294 RETURN
295 END SUBROUTINE CG_SOURCE_W_S
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312 SUBROUTINE CG_SOURCE_W_S_BC(A_M, B_M, M)
313
314
315
316
317 USE param, only: dimension_3, dimension_m
318 USE param1, only: zero, one, undefined
319
320 USE matrix
321 USE fldvar
322 USE bc
323
324 USE compar
325 USE geometry
326 USE indices
327
328 USE cutcell
329 USE quadric
330 USE fun_avg
331 USE functions
332
333 IMPLICIT NONE
334
335
336
337
338 INTEGER, INTENT(IN) :: M
339
340 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
341
342 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
343
344
345
346
347 INTEGER :: IJK, IJKB
348
349 INTEGER :: BCV
350 CHARACTER(LEN=9) :: BCT
351
352
353 IF(CG_SAFE_MODE(5)==1) RETURN
354
355 DO IJK = ijkstart3, ijkend3
356
357 BCV = BC_W_ID(IJK)
358 IF(BCV > 0 ) THEN
359 BCT = BC_TYPE(BCV)
360 ELSE
361 BCT = 'NONE'
362 ENDIF
363
364 SELECT CASE (BCT)
365 CASE ('CG_NSW')
366 IF(WALL_W_AT(IJK)) THEN
367 A_M(IJK,E,M) = ZERO
368 A_M(IJK,W,M) = ZERO
369 A_M(IJK,N,M) = ZERO
370 A_M(IJK,S,M) = ZERO
371 A_M(IJK,T,M) = ZERO
372 A_M(IJK,B,M) = ZERO
373 A_M(IJK,0,M) = -ONE
374 B_M(IJK,M) = ZERO
375 ENDIF
376
377
378 CASE ('CG_FSW')
379 IF(WALL_W_AT(IJK)) THEN
380 A_M(IJK,E,M) = ZERO
381 A_M(IJK,W,M) = ZERO
382 A_M(IJK,N,M) = ZERO
383 A_M(IJK,S,M) = ZERO
384 A_M(IJK,T,M) = ZERO
385 A_M(IJK,B,M) = ZERO
386 A_M(IJK,0,M) = -ONE
387
388 (IJK,M) = ZERO
389 IF(DABS(NORMAL_W(IJK,3))/=ONE) THEN
390 IF (W_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
391 A_M(IJK,E,M) = ONE
392 ELSEIF (W_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
393 A_M(IJK,W,M) = ONE
394 ELSEIF (W_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
395 A_M(IJK,N,M) = ONE
396 ELSEIF (W_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
397 A_M(IJK,S,M) = ONE
398 ELSEIF (W_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
399 A_M(IJK,T,M) = ONE
400 ELSEIF (W_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
401 A_M(IJK,B,M) = ONE
402 ENDIF
403 ENDIF
404 ENDIF
405
406
407 CASE ('CG_PSW')
408 IF(WALL_W_AT(IJK)) THEN
409 A_M(IJK,E,M) = ZERO
410 A_M(IJK,W,M) = ZERO
411 A_M(IJK,N,M) = ZERO
412 A_M(IJK,S,M) = ZERO
413 A_M(IJK,T,M) = ZERO
414 A_M(IJK,B,M) = ZERO
415 A_M(IJK,0,M) = -ONE
416 IF(BC_HW_S(BCV,M)==UNDEFINED) THEN
417 (IJK,M) = -BC_WW_S(BCV,M)
418 ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN
419 (IJK,M) = ZERO
420 IF(DABS(NORMAL_W(IJK,3))/=ONE) THEN
421 IF (W_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
422 A_M(IJK,E,M) = ONE
423 ELSEIF (W_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
424 A_M(IJK,W,M) = ONE
425 ELSEIF (W_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
426 A_M(IJK,N,M) = ONE
427 ELSEIF (W_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
428 A_M(IJK,S,M) = ONE
429 ELSEIF (W_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
430 A_M(IJK,T,M) = ONE
431 ELSEIF (W_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
432 A_M(IJK,B,M) = ONE
433 ENDIF
434 ENDIF
435 ELSE
436 (IJK,M) = ZERO
437 IF(DABS(NORMAL_W(IJK,3))/=ONE) THEN
438 IF (W_MASTER_OF(IJK) == EAST_OF(IJK)) THEN
439 A_M(IJK,E,M) = ONE
440 ELSEIF (W_MASTER_OF(IJK) == WEST_OF(IJK)) THEN
441 A_M(IJK,W,M) = ONE
442 ELSEIF (W_MASTER_OF(IJK) == NORTH_OF(IJK)) THEN
443 A_M(IJK,N,M) = ONE
444 ELSEIF (W_MASTER_OF(IJK) == SOUTH_OF(IJK)) THEN
445 A_M(IJK,S,M) = ONE
446 ELSEIF (W_MASTER_OF(IJK) == TOP_OF(IJK)) THEN
447 A_M(IJK,T,M) = ONE
448 ELSEIF (W_MASTER_OF(IJK) == BOTTOM_OF(IJK)) THEN
449 A_M(IJK,B,M) = ONE
450 ENDIF
451 ENDIF
452 ENDIF
453 ENDIF
454
455
456 CASE ('CG_MI')
457 A_M(IJK,E,M) = ZERO
458 A_M(IJK,W,M) = ZERO
459 A_M(IJK,N,M) = ZERO
460 A_M(IJK,S,M) = ZERO
461 A_M(IJK,T,M) = ZERO
462 A_M(IJK,B,M) = ZERO
463 A_M(IJK,0,M) = -ONE
464 IF(BC_W_s(BCV,M)/=UNDEFINED) THEN
465 B_M(IJK,M) = - BC_W_s(BCV,M)
466 ELSE
467 B_M(IJK,M) = - BC_VELMAG_s(BCV,M)*NORMAL_W(IJK,3)
468 ENDIF
469 IJKB = BOTTOM_OF(IJK)
470 IF(FLUID_AT(IJKB)) THEN
471 A_M(IJKB,E,M) = ZERO
472 A_M(IJKB,W,M) = ZERO
473 A_M(IJKB,N,M) = ZERO
474 A_M(IJKB,S,M) = ZERO
475 A_M(IJKB,T,M) = ZERO
476 A_M(IJKB,B,M) = ZERO
477 A_M(IJKB,0,M) = -ONE
478 IF(BC_W_s(BCV,M)/=UNDEFINED) THEN
479 B_M(IJKB,M) = - BC_W_s(BCV,M)
480 ELSE
481 B_M(IJKB,M) = - BC_VELMAG_s(BCV,M)*NORMAL_W(IJK,3)
482 ENDIF
483 ENDIF
484
485
486 CASE ('CG_PO')
487 A_M(IJK,E,M) = ZERO
488 A_M(IJK,W,M) = ZERO
489 A_M(IJK,N,M) = ZERO
490 A_M(IJK,S,M) = ZERO
491 A_M(IJK,T,M) = ZERO
492 A_M(IJK,B,M) = ZERO
493 A_M(IJK,0,M) = -ONE
494 B_M(IJK,M) = ZERO
495 IJKB = BOTTOM_OF(IJK)
496 IF(FLUID_AT(IJKB)) THEN
497 A_M(IJK,B,M) = ONE
498 A_M(IJK,0,M) = -ONE
499 ENDIF
500 END SELECT
501
502
503
504 = BC_ID(IJK)
505 IF(BCV > 0 ) THEN
506 BCT = BC_TYPE(BCV)
507 ELSE
508 BCT = 'NONE'
509 ENDIF
510
511
512 SELECT CASE (BCT)
513 CASE ('CG_MI')
514 A_M(IJK,E,M) = ZERO
515 A_M(IJK,W,M) = ZERO
516 A_M(IJK,N,M) = ZERO
517 A_M(IJK,S,M) = ZERO
518 A_M(IJK,T,M) = ZERO
519 A_M(IJK,B,M) = ZERO
520 A_M(IJK,0,M) = -ONE
521 IF(BC_W_s(BCV,M)/=UNDEFINED) THEN
522 B_M(IJK,M) = - BC_W_s(BCV,M)
523 ELSE
524 B_M(IJK,M) = - BC_VELMAG_s(BCV,M)*NORMAL_S(IJK,3)
525 ENDIF
526 IJKB = BOTTOM_OF(IJK)
527 IF(FLUID_AT(IJKB)) THEN
528 A_M(IJKB,E,M) = ZERO
529 A_M(IJKB,W,M) = ZERO
530 A_M(IJKB,N,M) = ZERO
531 A_M(IJKB,S,M) = ZERO
532 A_M(IJKB,T,M) = ZERO
533 A_M(IJKB,B,M) = ZERO
534 A_M(IJKB,0,M) = -ONE
535 IF(BC_W_s(BCV,M)/=UNDEFINED) THEN
536 B_M(IJKB,M) = - BC_W_s(BCV,M)
537 ELSE
538 B_M(IJKB,M) = - BC_VELMAG_s(BCV,M)*NORMAL_S(IJK,3)
539 ENDIF
540 ENDIF
541
542
543 CASE ('CG_PO')
544 A_M(IJK,E,M) = ZERO
545 A_M(IJK,W,M) = ZERO
546 A_M(IJK,N,M) = ZERO
547 A_M(IJK,S,M) = ZERO
548 A_M(IJK,T,M) = ZERO
549 A_M(IJK,B,M) = ZERO
550 A_M(IJK,0,M) = -ONE
551 B_M(IJK,M) = ZERO
552 IJKB = BOTTOM_OF(IJK)
553 IF(FLUID_AT(IJKB)) THEN
554 A_M(IJK,B,M) = ONE
555 A_M(IJK,0,M) = -ONE
556 ENDIF
557
558 END SELECT
559
560 ENDDO
561
562 RETURN
563 END SUBROUTINE CG_SOURCE_W_S_BC
564