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