File: /nfs/home/0/users/jenkins/mfix.git/model/GhdTheory/calc_d_ghd.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
32
33
34
35
36
37
38 SUBROUTINE CALC_D_ghd_E(A_M, VXF_GS, D_E)
39
40
41
42
43
44
45
46 USE param
47 USE param1
48 USE parallel
49 USE fldvar
50 USE geometry
51 USE indices
52 USE physprop
53 USE run
54 USE scales
55 USE compar
56 USE sendrecv
57 USE fun_avg
58 USE functions
59 IMPLICIT NONE
60
61
62
63
64
65
66
67
68 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
69
70
71 DOUBLE PRECISION VxF_gs(DIMENSION_3, DIMENSION_M)
72
73 DOUBLE PRECISION d_e(DIMENSION_3, 0:DIMENSION_M)
74
75
76 DOUBLE PRECISION EPGA
77
78 INTEGER M, L, I, IJK, IJKE
79
80 DOUBLE PRECISION EPSA(DIMENSION_M)
81
82 DOUBLE PRECISION other_ratio_1, FOA1
83
84 DOUBLE PRECISION SUM_VXF_GS
85
86 LOGICAL Pass1, Pass2
87
88 DOUBLE PRECISION numeratorxEP(DIMENSION_M)
89
90 DOUBLE PRECISION denominator(DIMENSION_M)
91
92 DOUBLE PRECISION other_denominator(DIMENSION_M)
93
94 DOUBLE PRECISION EPStmp
95
96
97 = .FALSE.
98 = .FALSE.
99 M = MMAX
100 if (MOMENTUM_X_EQ(0) .AND. MOMENTUM_X_EQ(M)) then
101 Pass1 = .TRUE.
102
103 elseif (MOMENTUM_X_EQ(M)) then
104 Pass2 = .TRUE.
105
106 endif
107
108
109 IF (Pass1) THEN
110
111
112
113
114 DO IJK = ijkstart3, ijkend3
115 IF (IP_AT_E(IJK) .OR. MFLOW_AT_E(IJK)) THEN
116 DO M= 0, MMAX
117 D_E(IJK,M) = ZERO
118 END DO
119 ELSE
120 I = I_OF(IJK)
121 IJKE = EAST_OF(IJK)
122 EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
123
124 SUM_VXF_GS = ZERO
125 EPSA(MMAX) = ZERO
126 DO M= 1, SMAX
127 EPSA(M) = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
128 EPSA(MMAX) = EPSA(MMAX) + EPSA(M)
129 SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)
130 END DO
131
132 M = MMAX
133 other_ratio_1 = ( VXF_GS(IJK,M)* (-A_M(IJK,0,M)) /&
134 ( (-A_M(IJK,0,M))+VXF_GS(IJK,M)+SMALL_NUMBER )&
135 )
136
137 IF (MODEL_B) THEN
138
139 if ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR. (other_ratio_1>SMALL_NUMBER) ) then
140 D_E(IJK,0) = P_SCALE*AYZ(IJK)/( (-A_M(IJK,0,0))+other_ratio_1 )
141 else
142 D_E(IJK,0) = ZERO
143 endif
144
145 = MMAX
146 if ( MOMENTUM_X_EQ(M) .AND. ( ((-A_M(IJK,0,M))>SMALL_NUMBER) .OR. &
147 (VXF_GS(IJK,M)>SMALL_NUMBER) ) ) then
148 D_E(IJK,M) = D_E(IJK,0)*(&
149 VXF_GS(IJK,M)/((-A_M(IJK,0,M))+VXF_GS(IJK,M))&
150 )
151 else
152 D_E(IJK,M) = ZERO
153 endif
154 ELSE
155 = ZERO
156 M = MMAX
157 FOA1 = FOA1 + (EPSA(M)*VXF_GS(IJK,M)/&
158 ( (-A_M(IJK,0,M))+VXF_GS(IJK,M)+SMALL_NUMBER )&
159 )
160 other_denominator(M) = VXF_GS(IJK,M)*( (-A_M(IJK,0,0))/&
161 ((-A_M(IJK,0,0))+SUM_VXF_GS+SMALL_NUMBER)&
162 )
163 numeratorxEP(M) = ZERO
164 denominator(M) = ZERO
165
166 if ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR. (other_ratio_1>SMALL_NUMBER) ) then
167 D_E(IJK,0) = P_SCALE*AYZ(IJK)*(EPGA+FOA1)/( (-A_M(IJK,0,0))+other_ratio_1 )
168 else
169 D_E(IJK,0) = ZERO
170 endif
171
172 = MMAX
173 if ( MOMENTUM_X_EQ(M) .AND. ( (-A_M(IJK,0,M)>SMALL_NUMBER) .OR. &
174 (other_denominator(M)>SMALL_NUMBER) ) ) then
175 D_E(IJK,M) = P_SCALE*AYZ(IJK)*(&
176 ( EPSA(M) + (VXF_GS(IJK,M)*EPGA/((-A_M(IJK,0,0))+SUM_VXF_GS+SMALL_NUMBER)) )/&
177 ( (-A_M(IJK,0,M))+other_denominator(M) )&
178 )
179 else
180 D_E(IJK,M) = ZERO
181 endif
182 ENDIF
183 ENDIF
184 ENDDO
185
186 ELSE IF (MOMENTUM_X_EQ(0)) THEN
187
188
189
190 DO IJK = ijkstart3, ijkend3
191 IF (IP_AT_E(IJK) .OR. MFLOW_AT_E(IJK)) THEN
192 (IJK,0) = ZERO
193 ELSE
194 I = I_OF(IJK)
195 IJKE = EAST_OF(IJK)
196 EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
197
198
199 SUM_VXF_GS = VXF_GS(IJK,MMAX)
200
201 IF ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR. (SUM_VXF_GS>SMALL_NUMBER) ) THEN
202 IF (MODEL_B) THEN
203 D_E(IJK,0) = P_SCALE*AYZ(IJK)/((-A_M(IJK,0,0))+SUM_VXF_GS)
204 ELSE
205 D_E(IJK,0) = P_SCALE*AYZ(IJK)*EPGA/((-A_M(IJK,0,0))+SUM_VXF_GS)
206 ENDIF
207 ELSE
208 D_E(IJK,0) = ZERO
209 ENDIF
210 ENDIF
211 END DO
212
213 ELSE IF (Pass2) THEN
214
215
216
217
218 DO IJK = ijkstart3, ijkend3
219 IF (IP_AT_E(IJK) .OR. MFLOW_AT_E(IJK) .OR. MODEL_B) THEN
220 DO M= 1, MMAX
221 D_E(IJK,M) = ZERO
222 END DO
223 ELSE
224 I = I_OF(IJK)
225 IJKE = EAST_OF(IJK)
226
227 M = MMAX
228 EPStmp = ZERO
229 DO L=1,SMAX
230 EPStmp = EPStmp+ AVG_X(EP_S(IJK,L),EP_S(IJKE,L),I)
231 ENDDO
232 EPSA(M) = EPStmp
233
234
235 = MMAX
236 if ( MOMENTUM_X_EQ(M) .AND. ( (-A_M(IJK,0,M)>SMALL_NUMBER) .OR. &
237 (VXF_GS(IJK,M)>SMALL_NUMBER) ) ) then
238 D_E(IJK,M) = P_SCALE*AYZ(IJK)*( EPSA(M) )/&
239 ( (-A_M(IJK,0,M))+VXF_GS(IJK,M) )
240 else
241 D_E(IJK,M) = ZERO
242 endif
243 ENDIF
244 END DO
245
246 ENDIF
247
248 RETURN
249 END SUBROUTINE CALC_D_ghd_E
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287 SUBROUTINE CALC_D_ghd_N(A_M, VXF_GS, D_N)
288
289
290
291
292
293
294
295 USE param
296 USE param1
297 USE parallel
298 USE fldvar
299 USE geometry
300 USE indices
301 USE physprop
302 USE run
303 USE scales
304 USE compar
305 USE sendrecv
306 USE fun_avg
307 USE functions
308 IMPLICIT NONE
309
310
311
312
313
314
315
316
317 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
318
319
320 DOUBLE PRECISION VxF_gs(DIMENSION_3, DIMENSION_M)
321
322 DOUBLE PRECISION d_n(DIMENSION_3, 0:DIMENSION_M)
323
324
325 DOUBLE PRECISION EPGA
326
327 INTEGER M, L, J, IJK, IJKN
328
329 DOUBLE PRECISION EPSA(DIMENSION_M)
330
331 DOUBLE PRECISION other_ratio_1, FOA1
332
333 DOUBLE PRECISION SUM_VXF_GS
334
335 LOGICAL Pass1, Pass2
336
337 DOUBLE PRECISION numeratorxEP(DIMENSION_M)
338
339 DOUBLE PRECISION denominator(DIMENSION_M)
340
341 DOUBLE PRECISION other_denominator(DIMENSION_M)
342
343 DOUBLE PRECISION EPStmp
344
345
346 = .FALSE.
347 = .FALSE.
348 M = MMAX
349 if (MOMENTUM_Y_EQ(0) .AND. MOMENTUM_Y_EQ(M)) then
350 Pass1 = .TRUE.
351
352 elseif (MOMENTUM_Y_EQ(M)) then
353 Pass2 = .TRUE.
354
355 endif
356
357
358 IF (Pass1) THEN
359
360
361
362
363 DO IJK = ijkstart3, ijkend3
364 IF (IP_AT_N(IJK) .OR. MFLOW_AT_N(IJK)) THEN
365 DO M= 0, MMAX
366 D_N(IJK,M) = ZERO
367 END DO
368 ELSE
369 J = J_OF(IJK)
370 IJKN = NORTH_OF(IJK)
371 EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
372
373 SUM_VXF_GS = ZERO
374 EPSA(MMAX) = ZERO
375 DO M= 1, SMAX
376 EPSA(M) = AVG_Y(EP_S(IJK,M),EP_S(IJKN,M),J)
377 EPSA(MMAX) = EPSA(MMAX) + EPSA(M)
378 SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)
379 END DO
380
381 M= MMAX
382 other_ratio_1 = ( VXF_GS(IJK,M)* (-A_M(IJK,0,M)) /&
383 ( (-A_M(IJK,0,M))+VXF_GS(IJK,M)+SMALL_NUMBER )&
384 )
385
386 IF (MODEL_B) THEN
387
388 if ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR. (other_ratio_1>SMALL_NUMBER) ) then
389 D_N(IJK,0) = P_SCALE*AXZ(IJK)/( (-A_M(IJK,0,0))+other_ratio_1 )
390 else
391 D_N(IJK,0) = ZERO
392 endif
393
394 = MMAX
395 if ( MOMENTUM_Y_EQ(M) .AND. ( ((-A_M(IJK,0,M))>SMALL_NUMBER) .OR. &
396 (VXF_GS(IJK,M)>SMALL_NUMBER) ) ) then
397 D_N(IJK,M) = D_N(IJK,0)*(&
398 VXF_GS(IJK,M)/((-A_M(IJK,0,M))+VXF_GS(IJK,M))&
399 )
400 else
401 D_N(IJK,M) = ZERO
402 endif
403 ELSE
404 = ZERO
405 M = MMAX
406 FOA1 = FOA1 + (EPSA(M)*VXF_GS(IJK,M)/&
407 ( (-A_M(IJK,0,M))+VXF_GS(IJK,M)+SMALL_NUMBER )&
408 )
409 other_denominator(M) = VXF_GS(IJK,M)*( (-A_M(IJK,0,0))/&
410 ((-A_M(IJK,0,0))+SUM_VXF_GS+SMALL_NUMBER)&
411 )
412 numeratorxEP(M) = ZERO
413 denominator(M) = ZERO
414
415 if ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR. (other_ratio_1>SMALL_NUMBER) ) then
416 D_N(IJK,0) = P_SCALE*AXZ(IJK)*(EPGA+FOA1)/( (-A_M(IJK,0,0))+other_ratio_1 )
417 else
418 D_N(IJK,0) = ZERO
419 endif
420
421 = MMAX
422 if ( MOMENTUM_Y_EQ(M) .AND. ( (-A_M(IJK,0,M)>SMALL_NUMBER) .OR. &
423 (other_denominator(M)>SMALL_NUMBER) ) ) then
424 D_N(IJK,M) = P_SCALE*AXZ(IJK)*(&
425 ( EPSA(M) + (VXF_GS(IJK,M)*EPGA/((-A_M(IJK,0,0))+SUM_VXF_GS+SMALL_NUMBER)) )/&
426 ( (-A_M(IJK,0,M))+other_denominator(M) )&
427 )
428 else
429 D_N(IJK,M) = ZERO
430 endif
431 ENDIF
432 ENDIF
433 END DO
434
435 ELSE IF (MOMENTUM_Y_EQ(0)) THEN
436
437
438
439 DO IJK = ijkstart3, ijkend3
440 IF (IP_AT_N(IJK) .OR. MFLOW_AT_N(IJK)) THEN
441 (IJK,0) = ZERO
442 ELSE
443 J = J_OF(IJK)
444 IJKN = NORTH_OF(IJK)
445 EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
446
447
448 SUM_VXF_GS = VXF_GS(IJK,MMAX)
449
450 IF ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR. (SUM_VXF_GS>SMALL_NUMBER) ) THEN
451 IF (MODEL_B) THEN
452 D_N(IJK,0) = P_SCALE*AXZ(IJK)/((-A_M(IJK,0,0))+SUM_VXF_GS)
453 ELSE
454 D_N(IJK,0) = P_SCALE*AXZ(IJK)*EPGA/((-A_M(IJK,0,0))+SUM_VXF_GS)
455 ENDIF
456 ELSE
457 D_N(IJK,0) = ZERO
458 ENDIF
459 ENDIF
460 END DO
461
462 ELSE IF (Pass2) THEN
463
464
465
466
467 DO IJK = ijkstart3, ijkend3
468 IF (IP_AT_N(IJK) .OR. MFLOW_AT_N(IJK) .OR. MODEL_B) THEN
469 DO M= 1, MMAX
470 D_N(IJK,M) = ZERO
471 END DO
472 ELSE
473 J = J_OF(IJK)
474 IJKN = NORTH_OF(IJK)
475
476 M = MMAX
477 EPStmp = ZERO
478 DO L=1,SMAX
479 EPStmp = EPStmp+ AVG_Y(EP_S(IJK,L),EP_S(IJKN,L),J)
480 ENDDO
481 EPSA(M) = EPStmp
482
483
484 = MMAX
485 if ( MOMENTUM_Y_EQ(M) .AND. ( (-A_M(IJK,0,M)>SMALL_NUMBER) .OR. &
486 (VXF_GS(IJK,M)>SMALL_NUMBER) ) ) then
487 D_N(IJK,M) = P_SCALE*AXZ(IJK)*( EPSA(M) )/&
488 ( (-A_M(IJK,0,M))+VXF_GS(IJK,M) )
489 else
490 D_N(IJK,M) = ZERO
491 endif
492 ENDIF
493 END DO
494
495 ENDIF
496
497 RETURN
498 END SUBROUTINE CALC_D_ghd_N
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536 SUBROUTINE CALC_D_ghd_T(A_M, VXF_GS, D_T)
537
538
539
540
541
542
543
544 USE param
545 USE param1
546 USE parallel
547 USE fldvar
548 USE geometry
549 USE indices
550 USE physprop
551 USE run
552 USE scales
553 USE compar
554 USE sendrecv
555 USE fun_avg
556 USE functions
557 IMPLICIT NONE
558
559
560
561
562
563
564
565
566 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
567
568
569 DOUBLE PRECISION VxF_gs(DIMENSION_3, DIMENSION_M)
570
571 DOUBLE PRECISION d_t(DIMENSION_3, 0:DIMENSION_M)
572
573
574 DOUBLE PRECISION EPGA
575
576 INTEGER M, L, K, IJK, IJKT
577
578 DOUBLE PRECISION EPSA(DIMENSION_M)
579
580 DOUBLE PRECISION other_ratio_1, FOA1
581
582 DOUBLE PRECISION SUM_VXF_GS
583
584 LOGICAL Pass1, Pass2
585
586 DOUBLE PRECISION numeratorxEP(DIMENSION_M)
587
588 DOUBLE PRECISION denominator(DIMENSION_M)
589
590 DOUBLE PRECISION other_denominator(DIMENSION_M)
591
592 DOUBLE PRECISION EPStmp
593
594
595 = .FALSE.
596 = .FALSE.
597 M = MMAX
598 if (MOMENTUM_Z_EQ(0) .AND. MOMENTUM_Z_EQ(M)) then
599 Pass1 = .TRUE.
600
601 elseif (MOMENTUM_Z_EQ(M)) then
602 Pass2 = .TRUE.
603
604 endif
605
606
607 IF (Pass1) THEN
608
609
610
611
612 DO IJK = ijkstart3, ijkend3
613 IF (IP_AT_T(IJK) .OR. MFLOW_AT_T(IJK)) THEN
614 DO M= 0, MMAX
615 D_T(IJK,M) = ZERO
616 END DO
617 ELSE
618 K = K_OF(IJK)
619 IJKT = TOP_OF(IJK)
620 EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
621
622 SUM_VXF_GS = ZERO
623 EPSA(MMAX) = ZERO
624 DO M= 1, SMAX
625 EPSA(M) = AVG_Z(EP_S(IJK,M),EP_S(IJKT,M),K)
626 EPSA(MMAX) = EPSA(MMAX) + EPSA(M)
627 SUM_VXF_GS = SUM_VXF_GS + VXF_GS(IJK,M)
628 END DO
629
630 M = MMAX
631 other_ratio_1 = ( VXF_GS(IJK,M)* (-A_M(IJK,0,M)) /&
632 ( (-A_M(IJK,0,M))+VXF_GS(IJK,M)+SMALL_NUMBER )&
633 )
634
635 IF (MODEL_B) THEN
636
637 if ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR. (other_ratio_1>SMALL_NUMBER) ) then
638 D_T(IJK,0) = P_SCALE*AXY(IJK)/( (-A_M(IJK,0,0))+other_ratio_1 )
639 else
640 D_T(IJK,0) = ZERO
641 endif
642
643 = MMAX
644 if ( MOMENTUM_Z_EQ(M) .AND. ( ((-A_M(IJK,0,M))>SMALL_NUMBER) .OR. &
645 (VXF_GS(IJK,M)>SMALL_NUMBER) ) ) then
646 D_T(IJK,M) = D_T(IJK,0)*(&
647 VXF_GS(IJK,M)/((-A_M(IJK,0,M))+VXF_GS(IJK,M))&
648 )
649 else
650 D_T(IJK,M) = ZERO
651 endif
652 ELSE
653 = ZERO
654 M = MMAX
655 FOA1 = FOA1 + (EPSA(M)*VXF_GS(IJK,M)/&
656 ( (-A_M(IJK,0,M))+VXF_GS(IJK,M)+SMALL_NUMBER )&
657 )
658 other_denominator(M) = VXF_GS(IJK,M)*( (-A_M(IJK,0,0))/&
659 ((-A_M(IJK,0,0))+SUM_VXF_GS+SMALL_NUMBER)&
660 )
661 numeratorxEP(M) = ZERO
662 denominator(M) = ZERO
663
664 if ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR. (other_ratio_1>SMALL_NUMBER) ) then
665 D_T(IJK,0) = P_SCALE*AXY(IJK)*(EPGA+FOA1)/( (-A_M(IJK,0,0))+other_ratio_1 )
666 else
667 D_T(IJK,0) = ZERO
668 endif
669
670 = MMAX
671 if ( MOMENTUM_Z_EQ(M) .AND. ( (-A_M(IJK,0,M)>SMALL_NUMBER) .OR. &
672 (other_denominator(M)>SMALL_NUMBER) ) ) then
673 D_T(IJK,M) = P_SCALE*AXY(IJK)*(&
674 ( EPSA(M) + (VXF_GS(IJK,M)*EPGA/((-A_M(IJK,0,0))+SUM_VXF_GS+SMALL_NUMBER)) )/&
675 ( (-A_M(IJK,0,M))+other_denominator(M) )&
676 )
677 else
678 D_T(IJK,M) = ZERO
679 endif
680 ENDIF
681 ENDIF
682 END DO
683
684 ELSE IF (MOMENTUM_Z_EQ(0)) THEN
685
686
687
688 DO IJK = ijkstart3, ijkend3
689 IF (IP_AT_T(IJK) .OR. MFLOW_AT_T(IJK)) THEN
690 (IJK,0) = ZERO
691 ELSE
692 K = K_OF(IJK)
693 IJKT = TOP_OF(IJK)
694 EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
695
696 SUM_VXF_GS = VXF_GS(IJK,MMAX)
697
698 IF ( ((-A_M(IJK,0,0))>SMALL_NUMBER) .OR. (SUM_VXF_GS>SMALL_NUMBER) ) THEN
699 IF (MODEL_B) THEN
700 D_T(IJK,0) = P_SCALE*AXY(IJK)/((-A_M(IJK,0,0))+SUM_VXF_GS)
701 ELSE
702 D_T(IJK,0) = P_SCALE*AXY(IJK)*EPGA/((-A_M(IJK,0,0))+SUM_VXF_GS)
703 ENDIF
704 ELSE
705 D_T(IJK,0) = ZERO
706 ENDIF
707 ENDIF
708 END DO
709
710 ELSE IF (Pass2) THEN
711
712
713
714
715 DO IJK = ijkstart3, ijkend3
716 IF (IP_AT_T(IJK) .OR. MFLOW_AT_T(IJK) .OR. MODEL_B) THEN
717 DO M= 1, MMAX
718 D_T(IJK,M) = ZERO
719 END DO
720 ELSE
721 K = K_OF(IJK)
722 IJKT = TOP_OF(IJK)
723
724 M = MMAX
725 EPStmp = ZERO
726 DO L=1,SMAX
727 EPStmp = EPStmp+ AVG_Z(EP_S(IJK,L),EP_S(IJKT,L),K)
728 ENDDO
729 EPSA(M) = EPStmp
730
731
732 = MMAX
733 if ( MOMENTUM_Z_EQ(M) .AND. ( (-A_M(IJK,0,M)>SMALL_NUMBER) .OR. &
734 (VXF_GS(IJK,M)>SMALL_NUMBER) ) ) then
735 D_T(IJK,M) = P_SCALE*AXY(IJK)*( EPSA(M) )/&
736 ( (-A_M(IJK,0,M))+VXF_GS(IJK,M) )
737 else
738 D_T(IJK,M) = ZERO
739 endif
740 ENDIF
741 END DO
742
743 ENDIF
744
745 RETURN
746 END SUBROUTINE CALC_D_ghd_T
747
748
749
750
751