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