File: /nfs/home/0/users/jenkins/mfix.git/model/source_u_s.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29 SUBROUTINE SOURCE_U_S(A_M, B_M, IER)
30
31
32
33
34 USE param
35 USE param1
36 USE parallel
37 USE matrix
38 USE scales
39 USE constant
40 USE physprop
41 USE fldvar
42 USE visc_s
43 USE rxns
44 USE run
45 USE toleranc
46 USE geometry
47 USE indices
48 USE is
49 USE tau_s
50 USE bc
51 USE compar
52 USE sendrecv
53 use kintheory
54 USE ghdtheory
55 USE drag
56 USE cutcell
57 USE quadric
58 USE mms
59 USE bodyforce
60 USE fun_avg
61 USE functions
62 IMPLICIT NONE
63
64
65
66
67 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
68
69
70 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
71
72 INTEGER, INTENT(INOUT) :: IER
73
74
75
76
77 INTEGER :: I, IJK,IMJK, IJMK, IJKE, IJKM, IPJK, IPJKM, &
78 J, K, IJPK, IPJMK, IJKP
79
80 INTEGER :: M, MM, L
81
82 INTEGER :: ISV
83
84 DOUBLE PRECISION :: PgE
85
86 DOUBLE PRECISION :: EPSA, EPStmp, epse, epsw, epsn, epss, &
87 epst, epsb, epsMix, epsMixE
88 DOUBLE PRECISION :: SUM_EPS_CP
89
90 DOUBLE PRECISION :: ROPSA
91
92 DOUBLE PRECISION :: dro1, dro2, droa
93
94 DOUBLE PRECISION :: wse, EPMUSA
95
96 DOUBLE PRECISION :: Sdp, Sdps
97
98 DOUBLE PRECISION :: V0, Vmt, Vbf, Vcf, Vtza, Vmttmp
99
100 DOUBLE PRECISION :: Ghd_drag, avgRop
101
102 DOUBLE PRECISION :: HYS_drag, avgDrag
103
104 DOUBLE PRECISION :: ROP_MA, Uge, Ugw, Vgw, Vge, Ugn,&
105 Ugs, Wgb, Wgt, Wge, Ugb, Ugt
106 DOUBLE PRECISION :: F_vir
107
108
109 DO M = 1, MMAX
110 IF(KT_TYPE_ENUM /= GHD_2007 .OR. &
111 (KT_TYPE_ENUM == GHD_2007 .AND. M==MMAX)) THEN
112
113 IF (MOMENTUM_X_EQ(M)) THEN
114
115
116
117
118
119
120
121
122
123
124
125 DO IJK = ijkstart3, ijkend3
126
127
128 IF(WALL_AT(IJK)) cycle
129
130 I = I_OF(IJK)
131 J = J_OF(IJK)
132 K = K_OF(IJK)
133 IJKE = EAST_OF(IJK)
134 IMJK = IM_OF(IJK)
135 IJMK = JM_OF(IJK)
136 IJKM = KM_OF(IJK)
137 IPJK = IP_OF(IJK)
138 IJPK = JP_OF(IJK)
139 IPJMK = IP_OF(IJMK)
140 IPJKM = IP_OF(IJKM)
141 IJKP = KP_OF(IJK)
142
143 IF (KT_TYPE_ENUM == GHD_2007) THEN
144
145 = ZERO
146 epsMix = ZERO
147 epsMixE= ZERO
148 DO L = 1, SMAX
149 EPStmp = EPStmp + AVG_X(EP_S(IJK,L),EP_S(IJKE,L),I)
150 epsMix = epsMix + EP_S(IJK,L)
151 = epsMixE + EP_S(IJKE,L)
152 IF(IP_AT_E(IJK)) THEN
153 U_S(IJK,L) = ZERO
154 ELSEIF(SIP_AT_E(IJK)) THEN
155 ISV = IS_ID_AT_E(IJK)
156 U_S(IJK,L) = IS_VEL_S(ISV,L)
157 ENDIF
158 ENDDO
159 EPSA = EPStmp
160 ELSE
161 EPSA = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
162 ENDIF
163
164
165 IF (IP_AT_E(IJK)) THEN
166 A_M(IJK,E,M) = ZERO
167 A_M(IJK,W,M) = ZERO
168 A_M(IJK,N,M) = ZERO
169 A_M(IJK,S,M) = ZERO
170 A_M(IJK,T,M) = ZERO
171 A_M(IJK,B,M) = ZERO
172 A_M(IJK,0,M) = -ONE
173 B_M(IJK,M) = ZERO
174
175
176 ELSEIF (SIP_AT_E(IJK)) THEN
177 A_M(IJK,E,M) = ZERO
178 A_M(IJK,W,M) = ZERO
179 A_M(IJK,N,M) = ZERO
180 A_M(IJK,S,M) = ZERO
181 A_M(IJK,T,M) = ZERO
182 A_M(IJK,B,M) = ZERO
183 A_M(IJK,0,M) = -ONE
184 ISV = IS_ID_AT_E(IJK)
185 B_M(IJK,M) = -IS_VEL_S(ISV,M)
186
187
188 ELSEIF (EPSA <= DIL_EP_S) THEN
189 A_M(IJK,E,M) = ZERO
190 A_M(IJK,W,M) = ZERO
191 A_M(IJK,N,M) = ZERO
192 A_M(IJK,S,M) = ZERO
193 A_M(IJK,T,M) = ZERO
194 A_M(IJK,B,M) = ZERO
195 A_M(IJK,0,M) = -ONE
196 B_M(IJK,M) = ZERO
197 IF (KT_TYPE_ENUM == GHD_2007) THEN
198 EPSw = ZERO
199 EPSe = ZERO
200 EPSn = ZERO
201 EPSs = ZERO
202 EPSt = ZERO
203 EPSb = ZERO
204 DO L = 1, SMAX
205 EPSw = EPSw + EP_S(WEST_OF(IJK),L)
206 EPSe = EPSe + EP_S(EAST_OF(IJK),L)
207 EPSn = EPSn + EP_S(NORTH_OF(IJK),L)
208 EPSs = EPSs + EP_S(SOUTH_OF(IJK),L)
209 IF(.NOT. NO_K) THEN
210 EPSt = EPSt + EP_S(TOP_OF(IJK),L)
211 EPSb = EPSb + EP_S(BOTTOM_OF(IJK),L)
212 ENDIF
213 ENDDO
214 ELSE
215 EPSw = EP_S(WEST_OF(IJK),M)
216 EPSe = EP_S(EAST_OF(IJK),M)
217 EPSn = EP_S(NORTH_OF(IJK),M)
218 EPSs = EP_S(SOUTH_OF(IJK),M)
219 IF(.NOT. NO_K) THEN
220 EPSt = EP_S(TOP_OF(IJK),M)
221 EPSb = EP_S(BOTTOM_OF(IJK),M)
222 ENDIF
223 ENDIF
224
225 IF (EPSw > DIL_EP_S .AND. .NOT.IS_AT_E(IMJK)) A_M(IJK,W,M) = ONE
226 IF (EPSe > DIL_EP_S .AND. .NOT.IS_AT_E(IJK)) A_M(IJK,E,M) = ONE
227 IF (EPSs > DIL_EP_S .AND. .NOT.IS_AT_N(IJMK)) A_M(IJK,S,M) = ONE
228 IF (EPSn > DIL_EP_S .AND. .NOT.IS_AT_N(IJK)) A_M(IJK,N,M) = ONE
229 IF(.NOT. NO_K) THEN
230 IF (EPSb > DIL_EP_S .AND. .NOT.IS_AT_T(IJKM)) A_M(IJK,B,M) = ONE
231 IF (EPSt > DIL_EP_S .AND. .NOT.IS_AT_T(IJK)) A_M(IJK,T,M) = ONE
232 ENDIF
233 IF((A_M(IJK,W,M)+A_M(IJK,E,M)+A_M(IJK,S,M)+A_M(IJK,N,M)+ &
234 A_M(IJK,B,M)+A_M(IJK,T,M)) == ZERO) THEN
235 B_M(IJK,M) = -U_S(IJK,M)
236 ELSE
237 A_M(IJK,0,M) = -(A_M(IJK,E,M)+A_M(IJK,W,M)+A_M(IJK,N,M)+ &
238 A_M(IJK,S,M)+A_M(IJK,T,M)+A_M(IJK,B,M))
239 ENDIF
240
241
242 ELSEIF (BLOCKED_U_CELL_AT(IJK)) THEN
243 A_M(IJK,E,M) = ZERO
244 A_M(IJK,W,M) = ZERO
245 A_M(IJK,N,M) = ZERO
246 A_M(IJK,S,M) = ZERO
247 A_M(IJK,T,M) = ZERO
248 A_M(IJK,B,M) = ZERO
249 A_M(IJK,0,M) = -ONE
250 B_M(IJK,M) = ZERO
251
252
253 ELSE
254
255
256
257 = P_G(IJKE)
258 IF (CYCLIC_X_PD) THEN
259
260
261
262 IF (IMAP(I_OF(IJK)).EQ.IMAX1) PGE = P_G(IJKE) - DELP_X
263 ENDIF
264 IF (MODEL_B) THEN
265 SDP = ZERO
266 ELSE
267 IF(.NOT.CUT_U_TREATMENT_AT(IJK)) THEN
268 SDP = -P_SCALE*EPSA*(PGE - P_G(IJK))*AYZ(IJK)
269 ELSE
270 SDP = -P_SCALE*EPSA*(PGE * A_UPG_E(IJK) - P_G(IJK) * A_UPG_W(IJK) )
271 ENDIF
272 ENDIF
273
274 IF (CLOSE_PACKED(M)) THEN
275 IF(SMAX > 1 .AND. KT_TYPE_ENUM /= GHD_2007) THEN
276 SUM_EPS_CP=0.0
277 DO MM=1,SMAX
278 IF (CLOSE_PACKED(MM))&
279 SUM_EPS_CP=SUM_EPS_CP+AVG_X(EP_S(IJK,MM),EP_S(IJKE,MM),I)
280 ENDDO
281 SUM_EPS_CP = Max(SUM_EPS_CP, small_number)
282 SDPS = -( (P_S(IJKE,M)-P_S(IJK,M))+(EPSA/SUM_EPS_CP)*&
283 (P_STAR(IJKE)-P_STAR(IJK)) )*AYZ(IJK)
284 ELSE
285 IF(.NOT.CUT_U_TREATMENT_AT(IJK)) THEN
286 SDPS =-((P_S(IJKE,M)-P_S(IJK,M))+(P_STAR(IJKE)-P_STAR(IJK)))*AYZ(IJK)
287 ELSE
288 SDPS =-((P_S(IJKE,M)* A_UPG_E(IJK)-P_S(IJK,M)* A_UPG_W(IJK)) &
289 +(P_STAR(IJKE)* A_UPG_E(IJK)-P_STAR(IJK)* A_UPG_W(IJK)))
290 ENDIF
291 ENDIF
292 ELSE
293 IF(.NOT.CUT_U_TREATMENT_AT(IJK)) THEN
294 SDPS = -(P_S(IJKE,M)-P_S(IJK,M))*AYZ(IJK)
295 ELSE
296 SDPS = - (P_S(IJKE,M) * A_UPG_E(IJK) - P_S(IJK,M) * A_UPG_W(IJK))
297 ENDIF
298 ENDIF
299
300
301 IF(.NOT.CUT_U_TREATMENT_AT(IJK)) THEN
302
303 = HALF * (VOL(IJK)*ROP_S(IJK,M) + &
304 VOL(IPJK)*ROP_S(IJKE,M))/VOL_U(IJK)
305
306 = HALF * (VOL(IJK)*ROP_SO(IJK,M) +&
307 VOL(IPJK)*ROP_SO(IJKE,M))*ODT/VOL_U(IJK)
308
309 IF(Added_Mass .AND. M==M_AM) THEN
310 ROP_MA = AVG_X(ROP_g(IJK)*EP_s(IJK,M),ROP_g(IJKE)*EP_s(IJKE,M),I)
311 V0 = V0 + Cv * ROP_MA * ODT
312 ENDIF
313 ELSE
314
315 = (VOL(IJK)*ROP_S(IJK,M) + &
316 VOL(IPJK)*ROP_S(IJKE,M))/(VOL(IJK) + VOL(IPJK))
317
318 = (VOL(IJK)*ROP_SO(IJK,M) + &
319 VOL(IPJK)*ROP_SO(IJKE,M))*ODT/&
320 (VOL(IJK) + VOL(IPJK))
321
322 IF(Added_Mass .AND. M==M_AM) THEN
323 ROP_MA = (VOL(IJK)*ROP_g(IJK)*EP_s(IJK,M) +&
324 VOL(IPJK)*ROP_g(IJKE)*EP_s(IJKE,M))/&
325 (VOL(IJK) + VOL(IPJK))
326 V0 = V0 + Cv * ROP_MA * ODT
327 ENDIF
328 ENDIF
329
330
331
332 = ZERO
333 IF(Added_Mass .AND. M==M_AM .AND.&
334 (.NOT.CUT_U_TREATMENT_AT(IJK))) THEN
335 F_vir = ( (U_G(IJK) - U_GO(IJK)) )*ODT*VOL_U(IJK)
336
337
338 = AVG_X_E(U_G(IMJK),U_G(IJK),I)
339 Uge = AVG_X_E(U_G(IJK),U_G(IPJK),IP1(I))
340 Vgw = AVG_Y_N(V_G(IJMK),V_G(IJK))
341 Vge = AVG_Y_N(V_G(IPJMK),V_G(IPJK))
342 Ugs = AVG_Y(U_G(IJMK),U_G(IJK),JM1(J))
343 Ugn = AVG_Y(U_G(IJK),U_G(IJPK),J)
344 IF(DO_K) THEN
345 Wgb = AVG_Z_T(W_g(IJKM),W_g(IJK))
346 Wgt = AVG_Z_T(W_g(IPJKM),W_g(IPJK))
347 Wge = AVG_X(Wgb,Wgt,I)
348 Ugb = AVG_Z(U_g(IJKM),U_g(IJK),KM1(K))
349 Ugt = AVG_Z(U_g(IJK),U_g(IJKP),K)
350 F_vir = F_vir + Wge*OX_E(I) * &
351 (Ugt - Ugb) *AXY(IJK)
352
353 IF (CYLINDRICAL) F_vir = F_vir - &
354 Wge**2*OX_E(I)
355 ENDIF
356
357 = F_vir + U_g(IJK)*(Uge - Ugw) *AYZ(IJK) + &
358 AVG_X(Vgw,Vge,I) * (Ugn - Ugs) *AXZ(IJK)
359 F_vir = F_vir * Cv * ROP_MA
360 ENDIF
361
362
363
364 IF (KT_TYPE_ENUM == GHD_2007) THEN
365 VMTtmp = ZERO
366 DO L = 1,SMAX
367 VMTtmp = VMTtmp + HALF*(VOL(IJK)*SUM_R_S(IJK,L) + &
368 VOL(IPJK)*SUM_R_S(IJKE,L))/VOL_U(IJK)
369 ENDDO
370 VMT = VMTtmp
371 ELSE
372 IF(.NOT.CUT_U_TREATMENT_AT(IJK)) THEN
373 VMT = HALF * (VOL(IJK)*SUM_R_S(IJK,M) + &
374 VOL(IPJK)*SUM_R_S(IJKE,M))/VOL_U(IJK)
375 ELSE
376 VMT = (VOL(IJK)*SUM_R_S(IJK,M) + &
377 VOL(IPJK)*SUM_R_S(IJKE,M))/(VOL(IJK) + VOL(IPJK))
378 ENDIF
379 ENDIF
380
381
382 IF (MODEL_B) THEN
383 IF (KT_TYPE_ENUM == GHD_2007) THEN
384 DRO1 = ROP_S(IJK,M) - RO_G(IJK) *epsMix
385 DRO2 = ROP_S(IJKE,M) - RO_G(IJKE)*epsMixE
386 DROA = AVG_X(DRO1,DRO2,I)
387 VBF = DROA*BFX_S(IJK,M)
388 ELSE
389 DRO1 = (RO_S(IJK,M)-RO_G(IJK))*EP_S(IJK,M)
390 DRO2 = (RO_S(IJK,M)-RO_G(IJKE))*EP_S(IJKE,M)
391 DROA = AVG_X(DRO1,DRO2,I)
392 VBF = DROA*BFX_S(IJK,M)
393 ENDIF
394 ELSE
395 = ROPSA*BFX_S(IJK,M)
396 ENDIF
397
398
399 = ZERO
400 IF (KT_TYPE_ENUM == GHD_2007) THEN
401 DO L = 1,SMAX
402 avgRop = AVG_X(ROP_S(IJK,L),ROP_S(IJKE,L),I)
403 if(avgRop > ZERO) Ghd_drag = Ghd_drag - &
404 AVG_X(F_GS(IJK,L),F_GS(IJKE,L),I) * &
405 JoiX(IJK,L) / avgRop
406 ENDDO
407 ENDIF
408
409
410 = ZERO
411 IF (DRAG_TYPE_ENUM .EQ. HYS .AND. &
412 KT_TYPE_ENUM /= GHD_2007) THEN
413 DO L = 1,MMAX
414 IF (L /= M) THEN
415 avgDrag = AVG_X(beta_ij(IJK,M,L),beta_ij(IJKE,M,L),I)
416 HYS_drag = HYS_drag - avgDrag * &
417 (U_g(ijk) - U_s(IJK,L))
418 ENDIF
419 ENDDO
420 ENDIF
421
422
423
424 IF (CYLINDRICAL) THEN
425
426 = AVG_X(HALF*(W_S(IJK,M)+W_S(IJKM,M)),&
427 HALF*(W_S(IPJK,M)+W_S(IPJKM,M)),I)
428 VCF = ROPSA*WSE**2*OX_E(I)
429 IF(Added_Mass .AND. M==M_AM) &
430 VCF = VCF + Cv*ROP_MA*WSE**2*OX_E(I)
431
432 = AVG_X(MU_S(IJK,M),MU_S(IJKE,M),I)
433 VTZA = 2.d0*EPMUSA*OX_E(I)*OX_E(I)
434 ELSE
435 VCF = ZERO
436 VTZA = ZERO
437 ENDIF
438
439
440 (IJK,0,M) = -(A_M(IJK,E,M)+A_M(IJK,W,M)+&
441 A_M(IJK,N,M)+A_M(IJK,S,M)+A_M(IJK,T,M)+&
442 A_M(IJK,B,M)+(V0+ZMAX(VMT)+VTZA)*VOL_U(IJK))
443
444 B_M(IJK,M) = B_M(IJK,M) - (SDP + SDPS + &
445 TAU_U_S(IJK,M) + F_vir + &
446 ( (V0+ZMAX((-VMT)))*U_SO(IJK,M)+&
447 VBF + VCF + HYS_drag)*VOL_U(IJK) )
448
449 IF(USE_MMS)B_M(IJK,M) = &
450 B_M(IJK,M) - MMS_U_S_SRC(IJK)*VOL_U(IJK)
451
452 IF (KT_TYPE_ENUM == IA_2005) THEN
453 B_M(IJK,M) = B_M(IJK,M) - KTMOM_U_S(IJK,M)
454 ELSEIF (KT_TYPE_ENUM == GHD_2007) THEN
455 B_M(IJK,M) = B_M(IJK,M) - Ghd_drag*VOL_U(IJK)
456 ENDIF
457
458 ENDIF
459 ENDDO
460
461
462
463 IF(CARTESIAN_GRID) CALL CG_SOURCE_U_S (A_M, B_M, M, IER)
464
465 CALL SOURCE_U_S_BC (A_M, B_M, M)
466 IF(CARTESIAN_GRID) CALL CG_SOURCE_U_S_BC (A_M, B_M, M, IER)
467
468 ENDIF
469 ENDIF
470 ENDDO
471
472 RETURN
473 END SUBROUTINE SOURCE_U_S
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494 SUBROUTINE SOURCE_U_S_BC(A_M, B_M, M)
495
496
497
498
499 USE param
500 USE param1
501 USE parallel
502 USE matrix
503 USE scales
504 USE constant
505 USE physprop
506 USE fldvar
507 USE visc_s
508 USE rxns
509 USE run
510 USE toleranc
511 USE geometry
512 USE indices
513 USE is
514 USE tau_s
515 USE bc
516 USE output
517 USE compar
518 USE fun_avg
519 USE functions
520 IMPLICIT NONE
521
522
523
524
525 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
526
527 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
528
529 INTEGER, INTENT(IN) :: M
530
531
532
533
534 INTEGER :: L
535
536 INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK,&
537 IM, JM, KM, IJKW, IMJK, IPJK, IP
538
539
540
541 IF (DO_K) THEN
542 K1 = 1
543 DO J1 = jmin3,jmax3
544 DO I1 = imin3, imax3
545 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
546 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE
547 = FUNIJK(I1,J1,K1)
548 IF (NS_WALL_AT(IJK)) THEN
549
550 (IJK,E,M) = ZERO
551 A_M(IJK,W,M) = ZERO
552 A_M(IJK,N,M) = ZERO
553 A_M(IJK,S,M) = ZERO
554 A_M(IJK,T,M) = -ONE
555 A_M(IJK,B,M) = ZERO
556 A_M(IJK,0,M) = -ONE
557 B_M(IJK,M) = ZERO
558 ELSEIF (FS_WALL_AT(IJK)) THEN
559
560 (IJK,E,M) = ZERO
561 A_M(IJK,W,M) = ZERO
562 A_M(IJK,N,M) = ZERO
563 A_M(IJK,S,M) = ZERO
564 A_M(IJK,T,M) = ONE
565 A_M(IJK,B,M) = ZERO
566 A_M(IJK,0,M) = -ONE
567 B_M(IJK,M) = ZERO
568 ENDIF
569 ENDDO
570 ENDDO
571
572 K1 = KMAX2
573 DO J1 = jmin3,jmax3
574 DO I1 = imin3, imax3
575 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
576 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE
577 = FUNIJK(I1,J1,K1)
578 IF (NS_WALL_AT(IJK)) THEN
579 A_M(IJK,E,M) = ZERO
580 A_M(IJK,W,M) = ZERO
581 A_M(IJK,N,M) = ZERO
582 A_M(IJK,S,M) = ZERO
583 A_M(IJK,T,M) = ZERO
584 A_M(IJK,B,M) = -ONE
585 A_M(IJK,0,M) = -ONE
586 B_M(IJK,M) = ZERO
587 ELSEIF (FS_WALL_AT(IJK)) THEN
588 A_M(IJK,E,M) = ZERO
589 A_M(IJK,W,M) = ZERO
590 A_M(IJK,N,M) = ZERO
591 A_M(IJK,S,M) = ZERO
592 A_M(IJK,T,M) = ZERO
593 A_M(IJK,B,M) = ONE
594 A_M(IJK,0,M) = -ONE
595 B_M(IJK,M) = ZERO
596 ENDIF
597 ENDDO
598 ENDDO
599 ENDIF
600
601 J1 = 1
602 DO K1 = kmin3, kmax3
603 DO I1 = imin3, imax3
604 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
605 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE
606 = FUNIJK(I1,J1,K1)
607 IF (NS_WALL_AT(IJK)) THEN
608 A_M(IJK,E,M) = ZERO
609 A_M(IJK,W,M) = ZERO
610 A_M(IJK,N,M) = -ONE
611 A_M(IJK,S,M) = ZERO
612 A_M(IJK,T,M) = ZERO
613 A_M(IJK,B,M) = ZERO
614 A_M(IJK,0,M) = -ONE
615 B_M(IJK,M) = ZERO
616 ELSEIF (FS_WALL_AT(IJK)) THEN
617 A_M(IJK,E,M) = ZERO
618 A_M(IJK,W,M) = ZERO
619 A_M(IJK,N,M) = ONE
620 A_M(IJK,S,M) = ZERO
621 A_M(IJK,T,M) = ZERO
622 A_M(IJK,B,M) = ZERO
623 A_M(IJK,0,M) = -ONE
624 B_M(IJK,M) = ZERO
625 ENDIF
626 ENDDO
627 ENDDO
628
629 J1 = JMAX2
630 DO K1 = kmin3, kmax3
631 DO I1 = imin3, imax3
632 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
633 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE
634 = FUNIJK(I1,J1,K1)
635 IF (NS_WALL_AT(IJK)) THEN
636 A_M(IJK,E,M) = ZERO
637 A_M(IJK,W,M) = ZERO
638 A_M(IJK,N,M) = ZERO
639 A_M(IJK,S,M) = -ONE
640 A_M(IJK,T,M) = ZERO
641 A_M(IJK,B,M) = ZERO
642 A_M(IJK,0,M) = -ONE
643 B_M(IJK,M) = ZERO
644 ELSEIF (FS_WALL_AT(IJK)) THEN
645 A_M(IJK,E,M) = ZERO
646 A_M(IJK,W,M) = ZERO
647 A_M(IJK,N,M) = ZERO
648 A_M(IJK,S,M) = ONE
649 A_M(IJK,T,M) = ZERO
650 A_M(IJK,B,M) = ZERO
651 A_M(IJK,0,M) = -ONE
652 B_M(IJK,M) = ZERO
653 ENDIF
654 ENDDO
655 ENDDO
656
657
658
659
660 DO L = 1, DIMENSION_BC
661 IF (BC_DEFINED(L)) THEN
662
663
664 IF (BC_TYPE(L) == 'NO_SLIP_WALL') THEN
665 I1 = BC_I_W(L)
666 I2 = BC_I_E(L)
667 J1 = BC_J_S(L)
668 J2 = BC_J_N(L)
669 K1 = BC_K_B(L)
670 K2 = BC_K_T(L)
671 IF (BC_JJ_PS(L) == 0) THEN
672 DO K = K1, K2
673 DO J = J1, J2
674 DO I = I1, I2
675 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
676 IF (DEAD_CELL_AT(I,J,K)) CYCLE
677 = FUNIJK(I,J,K)
678 IF (.NOT.WALL_AT(IJK)) CYCLE
679 (IJK,E,M) = ZERO
680 A_M(IJK,W,M) = ZERO
681 A_M(IJK,N,M) = ZERO
682 A_M(IJK,S,M) = ZERO
683 A_M(IJK,T,M) = ZERO
684 A_M(IJK,B,M) = ZERO
685 A_M(IJK,0,M) = -ONE
686 B_M(IJK,M) = ZERO
687 IF (FLUID_AT(NORTH_OF(IJK))) THEN
688 A_M(IJK,N,M) = -ONE
689 ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
690 A_M(IJK,S,M) = -ONE
691 ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
692 A_M(IJK,T,M) = -ONE
693 ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
694 A_M(IJK,B,M) = -ONE
695 ENDIF
696 ENDDO
697 ENDDO
698 ENDDO
699 ELSE
700 CALL JJ_BC_U_S (I1, I2, J1, J2, K1, K2, L, M, A_M, B_M)
701 ENDIF
702
703 ELSEIF (BC_TYPE(L) == 'FREE_SLIP_WALL') THEN
704 I1 = BC_I_W(L)
705 I2 = BC_I_E(L)
706 J1 = BC_J_S(L)
707 J2 = BC_J_N(L)
708 K1 = BC_K_B(L)
709 K2 = BC_K_T(L)
710 IF (BC_JJ_PS(L) == 0) THEN
711 DO K = K1, K2
712 DO J = J1, J2
713 DO I = I1, I2
714 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
715 IF (DEAD_CELL_AT(I,J,K)) CYCLE
716 = FUNIJK(I,J,K)
717 IF (.NOT.WALL_AT(IJK)) CYCLE
718 (IJK,E,M) = ZERO
719 A_M(IJK,W,M) = ZERO
720 A_M(IJK,N,M) = ZERO
721 A_M(IJK,S,M) = ZERO
722 A_M(IJK,T,M) = ZERO
723 A_M(IJK,B,M) = ZERO
724 A_M(IJK,0,M) = -ONE
725 B_M(IJK,M) = ZERO
726 IF (FLUID_AT(NORTH_OF(IJK))) THEN
727 A_M(IJK,N,M) = ONE
728 ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
729 A_M(IJK,S,M) = ONE
730 ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
731 A_M(IJK,T,M) = ONE
732 ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
733 A_M(IJK,B,M) = ONE
734 ENDIF
735 ENDDO
736 ENDDO
737 ENDDO
738 ELSE
739 CALL JJ_BC_U_S (I1, I2, J1, J2, K1, K2, L, M, A_M, B_M)
740 ENDIF
741
742 ELSEIF (BC_TYPE(L) == 'PAR_SLIP_WALL') THEN
743 I1 = BC_I_W(L)
744 I2 = BC_I_E(L)
745 J1 = BC_J_S(L)
746 J2 = BC_J_N(L)
747 K1 = BC_K_B(L)
748 K2 = BC_K_T(L)
749 IF (BC_JJ_PS(L) == 0) THEN
750 DO K = K1, K2
751 DO J = J1, J2
752 DO I = I1, I2
753 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
754 IF (DEAD_CELL_AT(I,J,K)) CYCLE
755 = FUNIJK(I,J,K)
756 IF (.NOT.WALL_AT(IJK)) CYCLE
757 = JM1(J)
758 KM = KM1(K)
759 A_M(IJK,E,M) = ZERO
760 A_M(IJK,W,M) = ZERO
761 A_M(IJK,N,M) = ZERO
762 A_M(IJK,S,M) = ZERO
763 A_M(IJK,T,M) = ZERO
764 A_M(IJK,B,M) = ZERO
765 A_M(IJK,0,M) = -ONE
766 B_M(IJK,M) = ZERO
767 IF (FLUID_AT(NORTH_OF(IJK))) THEN
768 IF (BC_HW_S(L,M) == UNDEFINED) THEN
769 A_M(IJK,N,M) = -HALF
770 A_M(IJK,0,M) = -HALF
771 B_M(IJK,M) = -BC_UW_S(L,M)
772 ELSE
773 A_M(IJK,0,M) = -(HALF*BC_HW_S(L,M)+ODY_N(J))
774 A_M(IJK,N,M) = -(HALF*BC_HW_S(L,M)-ODY_N(J))
775 B_M(IJK,M) = -BC_HW_S(L,M)*BC_UW_S(L,M)
776 ENDIF
777 ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
778 IF (BC_HW_S(L,M) == UNDEFINED) THEN
779 A_M(IJK,S,M) = -HALF
780 A_M(IJK,0,M) = -HALF
781 B_M(IJK,M) = -BC_UW_S(L,M)
782 ELSE
783 A_M(IJK,S,M) = -(HALF*BC_HW_S(L,M)-ODY_N(JM))
784 A_M(IJK,0,M) = -(HALF*BC_HW_S(L,M)+ODY_N(JM))
785 B_M(IJK,M) = -BC_HW_S(L,M)*BC_UW_S(L,M)
786 ENDIF
787 ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
788 IF (BC_HW_S(L,M) == UNDEFINED) THEN
789 A_M(IJK,T,M) = -HALF
790 A_M(IJK,0,M) = -HALF
791 B_M(IJK,M) = -BC_UW_S(L,M)
792 ELSE
793 A_M(IJK,0,M) = -(HALF*BC_HW_S(L,M)+ODZ_T(K)*&
794 OX_E(I))
795 A_M(IJK,T,M) = -(HALF*BC_HW_S(L,M)-ODZ_T(K)*&
796 OX_E(I))
797 B_M(IJK,M) = -BC_HW_S(L,M)*BC_UW_S(L,M)
798 ENDIF
799 ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
800 IF (BC_HW_S(L,M) == UNDEFINED) THEN
801 A_M(IJK,B,M) = -HALF
802 A_M(IJK,0,M) = -HALF
803 B_M(IJK,M) = -BC_UW_S(L,M)
804 ELSE
805 A_M(IJK,B,M) = -(HALF*BC_HW_S(L,M)-ODZ_T(KM)*&
806 OX_E(I))
807 A_M(IJK,0,M) = -(HALF*BC_HW_S(L,M)+ODZ_T(KM)*&
808 OX_E(I))
809 B_M(IJK,M) = -BC_HW_S(L,M)*BC_UW_S(L,M)
810 ENDIF
811 ENDIF
812 ENDDO
813 ENDDO
814 ENDDO
815 ELSE
816 CALL JJ_BC_U_S (I1, I2, J1, J2, K1, K2, L, M, A_M, B_M)
817 ENDIF
818
819
820 ELSEIF (BC_TYPE(L)=='P_INFLOW' .OR. BC_TYPE(L)=='P_OUTFLOW') THEN
821 IF (BC_PLANE(L) == 'W') THEN
822 I1 = BC_I_W(L)
823 I2 = BC_I_E(L)
824 J1 = BC_J_S(L)
825 J2 = BC_J_N(L)
826 K1 = BC_K_B(L)
827 K2 = BC_K_T(L)
828 DO K = K1, K2
829 DO J = J1, J2
830 DO I = I1, I2
831 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
832 IF (DEAD_CELL_AT(I,J,K)) CYCLE
833 = FUNIJK(I,J,K)
834 A_M(IJK,E,M) = ZERO
835 A_M(IJK,W,M) = ONE
836 A_M(IJK,N,M) = ZERO
837 A_M(IJK,S,M) = ZERO
838 A_M(IJK,T,M) = ZERO
839 A_M(IJK,B,M) = ZERO
840 A_M(IJK,0,M) = -ONE
841 B_M(IJK,M) = ZERO
842 ENDDO
843 ENDDO
844 ENDDO
845 ENDIF
846
847 ELSEIF (BC_TYPE(L) == 'OUTFLOW') THEN
848 IF (BC_PLANE(L) == 'W') THEN
849 I1 = BC_I_W(L)
850 I2 = BC_I_E(L)
851 J1 = BC_J_S(L)
852 J2 = BC_J_N(L)
853 K1 = BC_K_B(L)
854 K2 = BC_K_T(L)
855 DO K = K1, K2
856 DO J = J1, J2
857 DO I = I1, I2
858 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
859 IF (DEAD_CELL_AT(I,J,K)) CYCLE
860 = FUNIJK(I,J,K)
861 A_M(IJK,E,M) = ZERO
862 A_M(IJK,W,M) = ONE
863 A_M(IJK,N,M) = ZERO
864 A_M(IJK,S,M) = ZERO
865 A_M(IJK,T,M) = ZERO
866 A_M(IJK,B,M) = ZERO
867 A_M(IJK,0,M) = -ONE
868 B_M(IJK,M) = ZERO
869 IM = IM1(I)
870 IMJK = IM_OF(IJK)
871 A_M(IMJK,E,M) = ZERO
872 A_M(IMJK,W,M) = X_E(IM)/X_E(IM1(IM))
873 A_M(IMJK,N,M) = ZERO
874 A_M(IMJK,S,M) = ZERO
875 A_M(IMJK,T,M) = ZERO
876 A_M(IMJK,B,M) = ZERO
877 A_M(IMJK,0,M) = -ONE
878 B_M(IMJK,M) = ZERO
879 ENDDO
880 ENDDO
881 ENDDO
882 ELSEIF (BC_PLANE(L) == 'E') THEN
883 I1 = BC_I_W(L)
884 I2 = BC_I_E(L)
885 J1 = BC_J_S(L)
886 J2 = BC_J_N(L)
887 K1 = BC_K_B(L)
888 K2 = BC_K_T(L)
889 DO K = K1, K2
890 DO J = J1, J2
891 DO I = I1, I2
892 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
893 IF (DEAD_CELL_AT(I,J,K)) CYCLE
894 = FUNIJK(I,J,K)
895 IP = IP1(I)
896 IPJK = IP_OF(IJK)
897 A_M(IPJK,E,M) = X_E(IP)/X_E(I)
898 A_M(IPJK,W,M) = ZERO
899 A_M(IPJK,N,M) = ZERO
900 A_M(IPJK,S,M) = ZERO
901 A_M(IPJK,T,M) = ZERO
902 A_M(IPJK,B,M) = ZERO
903 A_M(IPJK,0,M) = -ONE
904 B_M(IPJK,M) = ZERO
905 ENDDO
906 ENDDO
907 ENDDO
908 ENDIF
909
910
911 ELSE
912 I1 = BC_I_W(L)
913 I2 = BC_I_E(L)
914 J1 = BC_J_S(L)
915 J2 = BC_J_N(L)
916 K1 = BC_K_B(L)
917 K2 = BC_K_T(L)
918 DO K = K1, K2
919 DO J = J1, J2
920 DO I = I1, I2
921 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
922 IF (DEAD_CELL_AT(I,J,K)) CYCLE
923
924 = FUNIJK(I,J,K)
925 A_M(IJK,E,M) = ZERO
926 A_M(IJK,W,M) = ZERO
927 A_M(IJK,N,M) = ZERO
928 A_M(IJK,S,M) = ZERO
929 A_M(IJK,T,M) = ZERO
930 A_M(IJK,B,M) = ZERO
931 A_M(IJK,0,M) = -ONE
932 B_M(IJK,M) = -U_S(IJK,M)
933 IF (BC_PLANE(L) == 'W') THEN
934 IJKW = WEST_OF(IJK)
935 A_M(IJKW,E,M) = ZERO
936 A_M(IJKW,W,M) = ZERO
937 A_M(IJKW,N,M) = ZERO
938 A_M(IJKW,S,M) = ZERO
939 A_M(IJKW,T,M) = ZERO
940 A_M(IJKW,B,M) = ZERO
941 A_M(IJKW,0,M) = -ONE
942 B_M(IJKW,M) = -U_S(IJKW,M)
943 ENDIF
944 ENDDO
945 ENDDO
946 ENDDO
947
948 ENDIF
949 ENDIF
950 ENDDO
951
952 RETURN
953 END SUBROUTINE SOURCE_U_S_BC
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969 SUBROUTINE JJ_BC_U_S(I1, I2, J1, J2, K1, K2, L, M, A_M, B_M)
970
971
972
973
974 USE param
975 USE param1
976 USE parallel
977 USE matrix
978 USE scales
979 USE constant
980 USE physprop
981 USE fldvar
982 USE visc_s
983 USE rxns
984 USE run
985 USE toleranc
986 USE geometry
987 USE indices
988 USE is
989 USE tau_s
990 USE bc
991 USE output
992 USE compar
993 USE functions
994 IMPLICIT NONE
995
996
997
998
999 INTEGER, INTENT(IN) :: L
1000
1001 INTEGER, INTENT(IN) :: I1, I2, J1, J2, K1, K2
1002
1003 INTEGER, INTENT(IN) :: M
1004
1005 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
1006
1007 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
1008
1009
1010
1011
1012 INTEGER :: I, J, K, IJK, JM, KM, IPJK
1013
1014 DOUBLE PRECISION :: hw, gw, cw
1015
1016
1017 DO K = K1, K2
1018 DO J = J1, J2
1019 DO I = I1, I2
1020 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
1021 IF (DEAD_CELL_AT(I,J,K)) CYCLE
1022
1023 = FUNIJK(I,J,K)
1024 IF (.NOT.WALL_AT(IJK)) CYCLE
1025 = JM1(J)
1026 KM = KM1(K)
1027 A_M(IJK,E,M) = ZERO
1028 A_M(IJK,W,M) = ZERO
1029 A_M(IJK,N,M) = ZERO
1030 A_M(IJK,S,M) = ZERO
1031 A_M(IJK,T,M) = ZERO
1032 A_M(IJK,B,M) = ZERO
1033 A_M(IJK,0,M) = -ONE
1034 B_M(IJK,M) = ZERO
1035
1036 IF (FLUID_AT(NORTH_OF(IJK))) THEN
1037 IPJK = IP_OF(NORTH_OF(IJK))
1038 IF (WALL_AT(IPJK)) CYCLE
1039 IF (EP_S(NORTH_OF(IJK),M) <= DIL_EP_S) THEN
1040 A_M(IJK,N,M) = ONE
1041 ELSE
1042 IF (BC_JJ_PS(L) == 1) THEN
1043
1044
1045 CALL CALC_GRBDRY (IJK, NORTH_OF(IJK), 'N', 'U',&
1046 M, L, GW, HW, CW)
1047 ELSEIF (BC_JJ_PS(L) == 2) THEN
1048
1049
1050 = 0D0
1051 HW = 1D0
1052 CW = BC_UW_S(L,M)
1053 ELSE
1054
1055
1056 = 1D0
1057 CW = 0D0
1058 HW = 0D0
1059 ENDIF
1060 A_M(IJK,N,M) = -(HALF*HW - ODY_N(J)*GW)
1061 A_M(IJK,0,M) = -(HALF*HW + ODY_N(J)*GW)
1062 B_M(IJK,M) = -CW
1063 ENDIF
1064
1065 ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
1066 IPJK = IP_OF(SOUTH_OF(IJK))
1067 IF (WALL_AT(IPJK)) CYCLE
1068 IF (EP_S(SOUTH_OF(IJK),M) <= DIL_EP_S) THEN
1069 A_M(IJK,S,M) = ONE
1070 ELSE
1071 IF (BC_JJ_PS(L) == 1) THEN
1072 CALL CALC_GRBDRY (IJK, SOUTH_OF(IJK), 'S', 'U',&
1073 M, L, GW, HW, CW)
1074 ELSEIF (BC_JJ_PS(L) == 2) THEN
1075 GW = 0D0
1076 HW = 1D0
1077 CW = BC_UW_S(L,M)
1078 ELSE
1079 GW = 1D0
1080 CW = 0D0
1081 HW = 0D0
1082 ENDIF
1083 A_M(IJK,S,M) = -(HALF*HW - ODY_N(JM)*GW)
1084 A_M(IJK,0,M) = -(HALF*HW + ODY_N(JM)*GW)
1085 B_M(IJK,M) = -CW
1086 ENDIF
1087
1088 ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
1089 IPJK = IP_OF(TOP_OF(IJK))
1090 IF (WALL_AT(IPJK)) CYCLE
1091 IF (EP_S(TOP_OF(IJK),M) <= DIL_EP_S) THEN
1092 A_M(IJK,T,M) = ONE
1093 ELSE
1094 IF (BC_JJ_PS(L) == 1) THEN
1095 CALL CALC_GRBDRY (IJK, TOP_OF(IJK), 'T', 'U',&
1096 M, L, GW, HW, CW)
1097 ELSEIF (BC_JJ_PS(L) == 2) THEN
1098 GW = 0D0
1099 HW = 1D0
1100 CW = BC_UW_S(L,M)
1101 ELSE
1102 GW = 1D0
1103 CW = 0D0
1104 HW = 0D0
1105 ENDIF
1106 A_M(IJK,T,M) = -(HALF*HW - ODZ_T(K)*OX_E(I)*GW)
1107 A_M(IJK,0,M) = -(HALF*HW + ODZ_T(K)*OX_E(I)*GW)
1108 B_M(IJK,M) = -CW
1109 ENDIF
1110
1111 ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
1112 IPJK = IP_OF(BOTTOM_OF(IJK))
1113 IF (WALL_AT(IPJK)) CYCLE
1114 IF (EP_S(BOTTOM_OF(IJK),M) <= DIL_EP_S) THEN
1115 A_M(IJK,B,M) = ONE
1116 ELSE
1117 IF (BC_JJ_PS(L) == 1) THEN
1118 CALL CALC_GRBDRY (IJK, BOTTOM_OF(IJK), 'B', 'U',&
1119 M,L, GW, HW, CW)
1120 ELSEIF (BC_JJ_PS(L) == 2) THEN
1121 GW = 0D0
1122 HW = 1D0
1123 CW = BC_UW_S(L,M)
1124 ELSE
1125 GW = 1D0
1126 CW = 0D0
1127 HW = 0D0
1128 ENDIF
1129 A_M(IJK,B,M) = -(HALF*HW - ODZ_T(KM)*OX_E(I)*GW)
1130 A_M(IJK,0,M) = -(HALF*HW + ODZ_T(KM)*OX_E(I)*GW)
1131 B_M(IJK,M) = -CW
1132 ENDIF
1133 ENDIF
1134
1135 ENDDO
1136 ENDDO
1137 ENDDO
1138
1139 RETURN
1140 END SUBROUTINE JJ_BC_U_S
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153 SUBROUTINE POINT_SOURCE_U_S(A_M, B_M)
1154
1155
1156
1157
1158 use compar
1159 use constant
1160 use fldvar
1161 use geometry
1162 use indices
1163 use physprop
1164 use ps
1165 use run
1166 use functions
1167 IMPLICIT NONE
1168
1169
1170
1171
1172 DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
1173
1174 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
1175
1176
1177
1178
1179 INTEGER :: IJK, I, J, K
1180 INTEGER :: PSV, M
1181
1182 INTEGER :: lIE, lIW
1183
1184 DOUBLE PRECISION :: pSource
1185
1186
1187 do M=1 , MMAX
1188
1189
1190
1191 PS_LP: do PSV = 1, DIMENSION_PS
1192 if(.NOT.PS_DEFINED(PSV)) cycle PS_LP
1193 if(abs(PS_U_s(PSV,M)) < small_number) cycle PS_LP
1194
1195 if(PS_U_s(PSV,M) < 0.0d0) then
1196 lIW = PS_I_W(PSV) - 1
1197 lIE = PS_I_E(PSV) - 1
1198 else
1199 lIW = PS_I_W(PSV)
1200 lIE = PS_I_E(PSV)
1201 endif
1202
1203 do k = PS_K_B(PSV), PS_K_T(PSV)
1204 do j = PS_J_S(PSV), PS_J_N(PSV)
1205 do i = lIW, lIE
1206
1207 if(.NOT.IS_ON_myPE_plus2layers(I,J,K)) cycle
1208 IF (DEAD_CELL_AT(I,J,K)) CYCLE
1209
1210 = funijk(i,j,k)
1211 if(.NOT.fluid_at(ijk)) cycle
1212
1213 if(A_M(IJK,0,M) == -ONE .AND. &
1214 B_M(IJK,M) == -U_s(IJK,M)) then
1215 B_M(IJK,M) = -PS_U_s(PSV,M) * PS_VEL_MAG_S(PSV,M)
1216 else
1217 pSource = PS_MASSFLOW_S(PSV,M) * &
1218 (VOL(IJK)/PS_VOLUME(PSV))
1219
1220 B_M(IJK,M) = B_M(IJK,M) - pSource * &
1221 PS_U_s(PSV,M) * PS_VEL_MAG_S(PSV,M)
1222 endif
1223
1224 enddo
1225 enddo
1226 enddo
1227
1228 enddo PS_LP
1229
1230 enddo
1231
1232
1233 RETURN
1234 END SUBROUTINE POINT_SOURCE_U_S
1235