File: /nfs/home/0/users/jenkins/mfix.git/model/GhdTheory/calc_nflux.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !                                                                      C
4     !  Module name: CALC_NFLUX(IER)                                        C
5     !  Purpose: Calculate the convection fluxes based on total species     C
6     !           number density for use with GHD theory only.               C
7     !           Adapted from routine calc_mflux.                           C
8     !                                                                      C
9     !  Author: S. Benyahia                                Date: 31-DEC-08  C
10     !  Reviewer:                                          Date:            C
11     !                                                                      C
12     !                                                                      C
13     !  Literature/Document References:                                     C
14     !      Christine Hrenya handnotes and my own derivation                C
15     !                                                                      C
16     !  Variables referenced:                                               C
17     !  Variables modified:                                                 C
18     !                                                                      C
19     !  Local variables:                                                    C
20     !                                                                      C
21     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
22     !
23           SUBROUTINE CALC_NFLUX()
24     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
25     !...Switches: -xf
26     !
27     !  Include param.inc file to specify parameter values
28     !
29     !-----------------------------------------------
30     !   M o d u l e s
31     !-----------------------------------------------
32           USE param
33           USE param1
34           USE parallel
35           USE fldvar
36           USE geometry
37           USE physprop
38           USE constant
39           USE indices
40           USE mflux
41           USE compar
42           USE functions
43           IMPLICIT NONE
44     !-----------------------------------------------
45     !   G l o b a l   P a r a m e t e r s
46     !-----------------------------------------------
47     !-----------------------------------------------
48     !   D u m m y   A r g u m e n t s
49     !-----------------------------------------------
50     !
51     !
52     !
53     !                      Total Number density
54           DOUBLE PRECISION Ni_E(DIMENSION_3), Ni_N(DIMENSION_3), Ni_T(DIMENSION_3)
55     !
56     !                      Indices
57           INTEGER          IJK, IMJK, IJMK, IJKM, L
58     !
59     !                      Particle mass
60           DOUBLE PRECISION Mp_L
61     !
62     !  First compute total number density at faces
63     !!!$omp  parallel do private( IJK, IMJK, IJMK, IJKM) &
64     !!!$omp&  schedule(static)
65           DO IJK = ijkstart3, ijkend3
66             Ni_E(IJK) = ZERO
67             Ni_N(IJK) = ZERO
68             Ni_T(IJK) = ZERO
69     !
70             IF (FLUID_AT(IJK)) THEN
71                DO L = 1, SMAX
72                  MP_L = RO_S(IJK,L) * PI*D_P(IJK,L)**3/6d0
73                  Ni_E(IJK) = Ni_E(IJK) + ROP_sE(IJK,L) / MP_L
74                  IMJK = IM_OF(IJK)
75                  IF (.NOT.FLUID_AT(IMJK)) THEN
76                     Ni_E(IMJK) = Ni_E(IMJK) + ROP_sE(IMJK,L) / MP_L
77                  ENDIF
78                  Ni_N(IJK) = Ni_N(IJK) + ROP_sN(IJK,L) / MP_L
79                  IJMK = JM_OF(IJK)
80                  IF (.NOT.FLUID_AT(IJMK)) THEN
81                    Ni_N(IJMK) = Ni_N(IJMK) + ROP_sN(IJMK,L) / MP_L
82                  ENDIF
83     !
84                  IF (DO_K) THEN
85                    Ni_T(IJK) = Ni_T(IJK) + ROP_sT(IJK,L) / MP_L
86                    IJKM = KM_OF(IJK)
87                      IF (.NOT.FLUID_AT(IJKM)) THEN
88                        Ni_T(IJKM) = Ni_T(IJKM) + ROP_sT(IJKM,L) / MP_L
89                      ENDIF
90                  ENDIF
91                ENDDO
92             ENDIF
93           ENDDO
94     
95     !
96     !  Interpolate the face value of density for calculating the convection fluxes
97     !!!$omp  parallel do private( IJK, IMJK, IJMK, IJKM) &
98     !!!$omp&  schedule(static)
99           DO IJK = ijkstart3, ijkend3
100     !
101              IF (FLUID_AT(IJK)) THEN
102     !
103     !         East face (i+1/2, j, k)
104                 Flux_nE(IJK) = Ni_E(IJK)*AYZ(IJK)*U_s(IJK,MMAX)
105     !
106     !         North face (i, j+1/2, k)
107                 Flux_nN(IJK) = Ni_N(IJK)*AXZ(IJK)*V_s(IJK,MMAX)
108     !
109     !         Top face (i, j, k+1/2)
110                 IF (DO_K) THEN
111                    Flux_nT(IJK) = Ni_T(IJK)*AXY(IJK)*W_s(IJK,MMAX)
112                 ENDIF
113     !
114     !         West face (i-1/2, j, k)
115                 IMJK = IM_OF(IJK)
116                 IF (.NOT.FLUID_AT(IMJK)) THEN
117                    Flux_nE(IMJK) = Ni_E(IMJK)*AYZ(IMJK)*U_s(IMJK,MMAX)
118                 ENDIF
119     !
120     !         South face (i, j-1/2, k)
121                 IJMK = JM_OF(IJK)
122                 IF (.NOT.FLUID_AT(IJMK)) THEN
123                   Flux_nN(IJMK) = Ni_N(IJMK)*AXZ(IJMK)*V_s(IJMK,MMAX)
124                 ENDIF
125     !
126     !         Bottom face (i, j, k-1/2)
127                 IF (DO_K) THEN
128                    IJKM = KM_OF(IJK)
129                    IF (.NOT.FLUID_AT(IJKM)) THEN
130                      Flux_nT(IJKM) = Ni_T(IJKM)*AXY(IJKM)*W_s(IJKM,MMAX)
131                    ENDIF
132                 ENDIF
133              ENDIF
134           END DO
135     
136           RETURN
137           END SUBROUTINE CALC_NFLUX
138     
139     !// Comments on the modifications for DMP version implementation
140     !// 001 Include header file and common declarations for parallelization
141     !// 350 Changed do loop limits: 1,ijkmax2-> ijkstart3, ijkend3
142