File: N:\mfix\model\GhdTheory\partial_elim_ghd.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 SUBROUTINE PARTIAL_ELIM_GHD_U(VAR_G, VAR_S, VXF, A_M, B_M)
23
24
25
26
27
28
29 USE param
30 USE param1
31 USE parallel
32 USE geometry
33 USE physprop
34 USE indices
35 USE run
36 USE compar
37 USE drag
38 USE fldvar
39 USE fun_avg
40 USE functions
41 IMPLICIT NONE
42
43
44
45
46
47
48
49
50 INTEGER IJK, IMJK, IJMK, IJKM, IPJK, IJPK, IJKP,LM,I,IJKE
51
52
53 DOUBLE PRECISION Var_g(DIMENSION_3)
54
55
56 DOUBLE PRECISION Var_s(DIMENSION_3, DIMENSION_M)
57
58
59 DOUBLE PRECISION VxF(DIMENSION_3, DIMENSION_M)
60
61
62 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
63
64
65 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
66
67
68 DOUBLE PRECISION SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
69 INTEGER L, M,LP
70
71
72 DOUBLE PRECISION a(0:DIMENSION_M), BB(0:DIMENSION_M), F(0:DIMENSION_M,0:DIMENSION_M),&
73 Saxf(0:DIMENSION_M)
74
75
76
77
78
79 DO IJK = ijkstart3, ijkend3
80 IF (FLOW_AT_E(IJK)) THEN
81 IMJK = IM_OF(IJK)
82 IJMK = JM_OF(IJK)
83 IPJK = IP_OF(IJK)
84 IJPK = JP_OF(IJK)
85 F = ZERO
86 DO M=0, MMAX
87 IF ((M==0 .OR. M==MMAX) .AND. MOMENTUM_X_EQ(M)) THEN
88 A(M) = A_M(IJK,0,M)
89 BB(M)= B_M(IJK,M)
90
91 if (m .ne. 0) then
92 IF (MOMENTUM_X_EQ(0)) F(M,0)=-VXF(IJK,M)
93 F(0,M) = F(M,0)
94 end if
95 DO L =1, MMAX
96 IF(L==0 .OR. L==MMAX) THEN
97 IF ((L .NE. M) .AND. (M .NE. 0) .AND. MOMENTUM_X_EQ(L)) THEN
98 LM = FUNLM(L,M)
99 IF (.NOT.IP_AT_E(IJK)) THEN
100 I = I_OF(IJK)
101 IJKE = EAST_OF(IJK)
102
103 F(M,L) = -AVG_X(F_SS(IJK,LM),F_SS(IJKE,LM),I)*VOL_U(IJK)
104 ENDIF
105 END IF
106 F(L,M)=F(M,L)
107 ENDIF
108 END DO
109
110 IF (M == 0 ) THEN
111 SAXF(M) = -(A_M(IJK,east,M)*VAR_G(IPJK)+A_M(IJK,west,M)*VAR_G(IMJK&
112 )+A_M(IJK,north,M)*VAR_G(IJPK)+A_M(IJK,south,M)*VAR_G(IJMK))
113 ELSE
114 SAXF(M) = -(A_M(IJK,east,M)*VAR_S(IPJK,M)+A_M(IJK,west,M)*VAR_S(&
115 IMJK,M)+A_M(IJK,north,M)*VAR_S(IJPK,M)+A_M(IJK,south,M)*VAR_S(&
116 IJMK,M))
117 ENDIF
118
119 IF (DO_K) THEN
120 IJKM = KM_OF(IJK)
121 IJKP = KP_OF(IJK)
122 IF (M ==0) THEN
123 SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_G(IJKP)+A_M(IJK,bottom,M)*&
124 VAR_G(IJKM))
125 ELSE
126 SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_S(IJKP,M)+A_M(IJK,bottom,M)*&
127 VAR_S(IJKM,M))
128 ENDIF
129 ENDIF
130 ENDIF
131
132 END DO
133
134 Do M=0,MMAX
135
136 IF ((M==0 .OR. M==MMAX) .AND. MOMENTUM_X_EQ(M)) THEN
137 SUM_A = ZERO
138 SUM_B = ZERO
139 do L =0,MMAX
140 IF(L==0 .OR. L==MMAX) THEN
141 SUM_A_LPRIME = ZERO
142 SUM_B_LPRIME = ZERO
143 DO LP=0,MMAX
144 IF(LP==0 .OR. LP==MMAX) THEN
145 IF ( LP .NE. M) THEN
146 SUM_A_LPRIME=SUM_A_LPRIME+F(L,LP)
147 IF (LP == 0) THEN
148 SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_G(IJK)
149 ELSE
150 SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_S(IJK,LP)
151 END IF
152 END IF
153 ENDIF
154 END DO
155 DEN = A(L) + SUM_A_LPRIME + F(L,M)
156 IF ( DEN .NE. ZERO) THEN
157 SUM_A = SUM_A + ((F(L,M)*(A(L)+SUM_A_LPRIME))/DEN)
158 SUM_B = SUM_B + F(L,M)*(SAXF(L) + BB(L)+SUM_B_LPRIME &
159 )/DEN
160 ENDIF
161 ENDIF
162 END DO
163 A_M(IJK,0,M) = SUM_A+A(M)
164 B_M(IJK,M) = SUM_B+BB(M)
165
166 DO L = 0, MMAX
167 IF(L==0 .OR. L==MMAX) THEN
168 IF(.NOT.MOMENTUM_X_EQ(L))THEN
169 IF( L /= M) THEN
170 IF(L == 0)THEN
171 A_M(IJK,0,M) = A_M(IJK,0,M) - VXF(IJK,M)
172 B_M(IJK,M) = B_M(IJK,M) - VXF(IJK,M) * VAR_G(IJK)
173 ELSE IF(M .NE. 0) THEN
174 LM = FUNLM(L,M)
175 IF (.NOT.IP_AT_E(IJK)) THEN
176 I = I_OF(IJK)
177 IJKE = EAST_OF(IJK)
178
179 F(M,L) = -AVG_X(F_SS(IJK,LM),F_SS(IJKE,LM),I)*VOL_U(IJK)
180 A_M(IJK,0,M) = A_M(IJK,0,M) + F(M,L)
181 B_M(IJK,M) = B_M(IJK,M) + F(M,L) * VAR_S(IJK, L)
182 ENDIF
183 ENDIF
184 ENDIF
185 ENDIF
186 ENDIF
187 ENDDO
188
189 ENDIF
190 END DO
191 ENDIF
192 END DO
193
194 RETURN
195 END SUBROUTINE PARTIAL_ELIM_GHD_U
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216 SUBROUTINE PARTIAL_ELIM_GHD_V(VAR_G, VAR_S, VXF, A_M, B_M)
217
218
219
220
221
222
223 USE param
224 USE param1
225 USE parallel
226 USE geometry
227 USE physprop
228 USE indices
229 USE run
230 USE compar
231 USE drag
232 USE fldvar
233 USE fun_avg
234 USE functions
235 IMPLICIT NONE
236
237
238
239
240
241
242
243
244 INTEGER IJK, IMJK, IJMK, IJKM, IPJK, IJPK, IJKP, IJKN,J,LM
245
246
247
248
249
250
251 DOUBLE PRECISION Var_g(DIMENSION_3)
252
253
254 DOUBLE PRECISION Var_s(DIMENSION_3, DIMENSION_M)
255
256
257 DOUBLE PRECISION VxF(DIMENSION_3, DIMENSION_M)
258
259
260 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
261
262
263 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
264
265
266 DOUBLE PRECISION SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
267 INTEGER L, M,LP
268
269
270 DOUBLE PRECISION a(0:DIMENSION_M), BB(0:DIMENSION_M), F(0:DIMENSION_M,0:DIMENSION_M),&
271 Saxf(0:DIMENSION_M)
272
273
274
275
276
277
278
279 DO IJK = ijkstart3, ijkend3
280 IF (FLOW_AT_N(IJK)) THEN
281 IMJK = IM_OF(IJK)
282 IJMK = JM_OF(IJK)
283 IPJK = IP_OF(IJK)
284 IJPK = JP_OF(IJK)
285 F = ZERO
286 DO M = 0, MMAX
287
288 IF ((M==0 .OR. M==MMAX) .AND. MOMENTUM_Y_EQ(M)) THEN
289 A(M)=A_M(IJK,0,M)
290 BB(M)=B_M(IJK,M)
291
292 if (m .ne. 0) then
293 IF (MOMENTUM_Y_EQ(0)) F(M,0)=-VXF(IJK,M)
294 F(0,M)=F(M,0)
295 end if
296 DO L =1, MMAX
297 IF(L==0 .OR. L==MMAX) THEN
298 IF ((L .NE. M) .AND. (M .NE. 0).AND. MOMENTUM_Y_EQ(L)) THEN
299 LM = FUNLM(L,M)
300 IF (.NOT.IP_AT_N(IJK)) THEN
301 J = J_OF(IJK)
302 IJKN = NORTH_OF(IJK)
303 F(M,L)=-AVG_Y(F_SS(IJK,LM),F_SS(IJKN,LM),J)*VOL_V(IJK)
304 END IF
305 END IF
306 F(L,M)=F(M,L)
307 ENDIF
308 END DO
309
310 IF (M == 0 ) THEN
311 SAXF(M) = -(A_M(IJK,east,M)*VAR_G(IPJK)+A_M(IJK,west,M)*VAR_G(IMJK&
312 )+A_M(IJK,north,M)*VAR_G(IJPK)+A_M(IJK,south,M)*VAR_G(IJMK))
313 ELSE
314 SAXF(M) = -(A_M(IJK,east,M)*VAR_S(IPJK,M)+A_M(IJK,west,M)*VAR_S(&
315 IMJK,M)+A_M(IJK,north,M)*VAR_S(IJPK,M)+A_M(IJK,south,M)*VAR_S(&
316 IJMK,M))
317 ENDIF
318
319 IF (DO_K) THEN
320 IJKM = KM_OF(IJK)
321 IJKP = KP_OF(IJK)
322 IF (M == 0) THEN
323 SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_G(IJKP)+A_M(IJK,bottom,M)*&
324 VAR_G(IJKM))
325 ELSE
326 SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_S(IJKP,M)+A_M(IJK,bottom,M)*&
327 VAR_S(IJKM,M))
328 ENDIF
329 ENDIF
330 ENDIF
331 END DO
332
333 Do M=0,MMAX
334 IF ((M==0 .OR. M==MMAX) .AND. MOMENTUM_Y_EQ(M)) THEN
335 SUM_A = ZERO
336 SUM_B = ZERO
337 do L =0,MMAX
338 IF(L==0 .OR. L==MMAX) THEN
339 SUM_A_LPRIME = ZERO
340 SUM_B_LPRIME = ZERO
341 DO LP=0,MMAX
342 IF(LP==0 .OR. LP==MMAX) THEN
343 IF ( LP .NE. M) THEN
344 SUM_A_LPRIME=SUM_A_LPRIME+F(L,LP)
345 IF (LP == 0) THEN
346 SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_G(IJK)
347 ELSE
348 SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_S(IJK,LP)
349 END IF
350 END IF
351 ENDIF
352 END DO
353 DEN = A(L) + SUM_A_LPRIME + F(L,M)
354 IF ( DEN .NE. ZERO) THEN
355 SUM_A = SUM_A + ((F(L,M)*(A(L)+SUM_A_LPRIME))/DEN)
356 SUM_B = SUM_B + F(L,M)*(SAXF(L) + BB(L)+SUM_B_LPRIME &
357 )/DEN
358 ENDIF
359 ENDIF
360 END DO
361 A_M(IJK,0,M)=SUM_A+A(M)
362 B_M(IJK,M) = SUM_B+BB(M)
363
364 DO L = 0, MMAX
365 IF(L==0 .OR. L==MMAX) THEN
366 IF(.NOT.MOMENTUM_Y_EQ(L))THEN
367 IF( L /= M) THEN
368 IF(L == 0)THEN
369 A_M(IJK,0,M) = A_M(IJK,0,M) - VXF(IJK,M)
370 B_M(IJK,M) = B_M(IJK,M) - VXF(IJK,M) * VAR_G(IJK)
371 ELSE IF(M .NE. 0) THEN
372 LM = FUNLM(L,M)
373 IF (.NOT.IP_AT_N(IJK)) THEN
374 J = J_OF(IJK)
375 IJKN = NORTH_OF(IJK)
376 F(M,L)=-AVG_Y(F_SS(IJK,LM),F_SS(IJKN,LM),J)*VOL_V(IJK)
377 A_M(IJK,0,M) = A_M(IJK,0,M) + F(M,L)
378 B_M(IJK,M) = B_M(IJK,M) + F(M,L) * VAR_S(IJK, L)
379 ENDIF
380 ENDIF
381 ENDIF
382 ENDIF
383 ENDIF
384 ENDDO
385
386 ENDIF
387 END DO
388 ENDIF
389 END DO
390 RETURN
391 END SUBROUTINE PARTIAL_ELIM_GHD_V
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412 SUBROUTINE PARTIAL_ELIM_GHD_W(VAR_G, VAR_S, VXF, A_M, B_M)
413
414
415
416
417
418
419 USE param
420 USE param1
421 USE parallel
422 USE geometry
423 USE physprop
424 USE indices
425 USE run
426 USE compar
427 USE drag
428 USE fldvar
429 USE fun_avg
430 USE functions
431 IMPLICIT NONE
432
433
434
435
436
437
438
439
440 INTEGER IJK, IMJK, IJMK, IJKM, IPJK, IJPK, IJKP,LM,K, IJKT
441
442
443 DOUBLE PRECISION Var_g(DIMENSION_3)
444
445
446 DOUBLE PRECISION Var_s(DIMENSION_3, DIMENSION_M)
447
448
449 DOUBLE PRECISION VxF(DIMENSION_3, DIMENSION_M)
450
451
452 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
453
454
455 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
456
457 DOUBLE PRECISION SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
458 INTEGER L, M,LP
459
460
461 DOUBLE PRECISION a(0:DIMENSION_M), BB(0:DIMENSION_M), F(0:DIMENSION_M,0:DIMENSION_M),&
462 Saxf(0:DIMENSION_M)
463
464
465
466
467
468 DO IJK = ijkstart3, ijkend3
469 IF (FLOW_AT_T(IJK)) THEN
470 IMJK = IM_OF(IJK)
471 IJMK = JM_OF(IJK)
472 IPJK = IP_OF(IJK)
473 IJPK = JP_OF(IJK)
474
475 F = ZERO
476 DO M=0, MMAX
477 IF ((M==0 .OR. M==MMAX) .AND. MOMENTUM_Z_EQ(M)) THEN
478 A(M)=A_M(IJK,0,M)
479 BB(M)=B_M(IJK,M)
480
481 if (m .ne. 0) then
482 IF (MOMENTUM_Z_EQ(0))F(M,0)=-VXF(IJK,M)
483 F(0,M)=F(M,0)
484 end if
485 DO L =1, MMAX
486 IF(L==0 .OR. L==MMAX) THEN
487 IF ((L .NE. M) .AND. (M .NE. 0) .AND. MOMENTUM_Z_EQ(L) ) THEN
488 LM = FUNLM(L,M)
489 IF (.NOT.IP_AT_T(IJK)) THEN
490 K = K_OF(IJK)
491 IJKT = TOP_OF(IJK)
492 F(M,L) = -AVG_Z(F_SS(IJK,LM),F_SS(IJKT,LM),K)*VOL_W(IJK)
493 ENDIF
494 END IF
495 F(L,M)=F(M,L)
496 ENDIF
497 END DO
498
499 IF (M == 0 ) THEN
500 SAXF(M) = -(A_M(IJK,east,M)*VAR_G(IPJK)+A_M(IJK,west,M)*VAR_G(IMJK&
501 )+A_M(IJK,north,M)*VAR_G(IJPK)+A_M(IJK,south,M)*VAR_G(IJMK))
502 ELSE
503 SAXF(M) = -(A_M(IJK,east,M)*VAR_S(IPJK,M)+A_M(IJK,west,M)*VAR_S(&
504 IMJK,M)+A_M(IJK,north,M)*VAR_S(IJPK,M)+A_M(IJK,south,M)*VAR_S(&
505 IJMK,M))
506 ENDIF
507
508 IF (DO_K) THEN
509 IJKM = KM_OF(IJK)
510 IJKP = KP_OF(IJK)
511 IF (M == 0) THEN
512 SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_G(IJKP)+A_M(IJK,bottom,M)*&
513 VAR_G(IJKM))
514 ELSE
515 SAXF(M) = SAXF(M) - (A_M(IJK,top,M)*VAR_S(IJKP,M)+A_M(IJK,bottom,M)*&
516 VAR_S(IJKM,M))
517 ENDIF
518 ENDIF
519 ENDIF
520 END DO
521
522 Do M=0,MMAX
523 IF ((M==0 .OR. M==MMAX) .AND. MOMENTUM_Z_EQ(M)) THEN
524 SUM_A = ZERO
525 SUM_B = ZERO
526 do L =0,MMAX
527 IF(L==0 .OR. L==MMAX) THEN
528 SUM_A_LPRIME = ZERO
529 SUM_B_LPRIME = ZERO
530 DO LP=0,MMAX
531 IF(LP==0 .OR. LP==MMAX) THEN
532 IF ( LP .NE. M) THEN
533 SUM_A_LPRIME=SUM_A_LPRIME+F(L,LP)
534 IF (LP == 0) THEN
535 SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_G(IJK)
536 ELSE
537 SUM_B_LPRIME=SUM_B_LPRIME+F(L,LP)*VAR_S(IJK,LP)
538 END IF
539 END IF
540 ENDIF
541 END DO
542
543 DEN = A(L) + SUM_A_LPRIME + F(L,M)
544 IF ( DEN .NE. ZERO) THEN
545 SUM_A = SUM_A + ((F(L,M)*(A(L)+SUM_A_LPRIME))/DEN)
546 SUM_B = SUM_B + F(L,M)*(SAXF(L) + BB(L)+SUM_B_LPRIME &
547 )/DEN
548 ENDIF
549 ENDIF
550 END DO
551 A_M(IJK,0,M)=SUM_A+A(M)
552 B_M(IJK,M) = SUM_B+BB(M)
553
554 DO L = 0, MMAX
555 IF(L==0 .OR. L==MMAX) THEN
556 IF(.NOT.MOMENTUM_Z_EQ(L))THEN
557 IF( L /= M) THEN
558 IF(L == 0)THEN
559 A_M(IJK,0,M) = A_M(IJK,0,M) - VXF(IJK,M)
560 B_M(IJK,M) = B_M(IJK,M) - VXF(IJK,M) * VAR_G(IJK)
561 ELSE IF(M .NE. 0) THEN
562 LM = FUNLM(L,M)
563 IF (.NOT.IP_AT_T(IJK)) THEN
564 K = K_OF(IJK)
565 IJKT = TOP_OF(IJK)
566 F(M,L) = -AVG_Z(F_SS(IJK,LM),F_SS(IJKT,LM),K)*VOL_W(IJK)
567 A_M(IJK,0,M) = A_M(IJK,0,M) + F(M,L)
568 B_M(IJK,M) = B_M(IJK,M) + F(M,L) * VAR_S(IJK, L)
569 ENDIF
570 ENDIF
571 ENDIF
572 ENDIF
573 ENDIF
574 ENDDO
575
576 ENDIF
577 END DO
578 ENDIF
579 END DO
580 RETURN
581 END SUBROUTINE PARTIAL_ELIM_GHD_W
582
583
584