1 | include(dom.inc)
2 |
3 | SUBROUTINE SPATIAL_SCHEME(I_SCHEME, Lb, Gi, G, k_abs, k_scat, Qr, &
4 | & Li, H, alpha)
5 |
6 | USE mod_pmm
7 | USE mod_slave
8 | #ifdef USEPALM
9 | USE palmlib
10 | #endif
11 |
12 | IMPLICIT NONE
13 |
14 | include 'dom_constants.h'
15 |
16 | ! IN
17 | DOM_INT :: I_SCHEME
18 | DOM_REAL :: alpha
19 | DOM_REAL,DIMENSION(is_ncells) :: Lb
20 | DOM_REAL,DIMENSION(is_ncells) :: k_abs, k_scat
21 | DOM_REAL,DIMENSION(is_ncells) :: Gi
22 | DOM_REAL,DIMENSION(is_nbfaces) :: Lo
23 |
24 |
25 | ! OUT
26 | DOM_REAL,DIMENSION(is_ncells) :: G
27 | DOM_REAL,DIMENSION(3,is_ncells) :: Qr
28 | DOM_REAL,DIMENSION(is_nfacesmax,is_ncells) :: Li
29 | DOM_REAL,DIMENSION(is_nbfaces) :: H
30 |
31 | ! LOCAL
32 | DOM_INT :: ierr, inode, ibnd
33 | DOM_INT :: idir,icell,i,j,k
34 | DOM_INT :: neighfaceid, neighcellid
35 | DOM_REAL :: nmin
36 | DOM_REAL :: Lpi
37 | DOM_REAL,DIMENSION(is_nfacesmax) :: Dm_ij
38 | DOM_REAL,DIMENSION(is_ncells) :: localG
39 | DOM_REAL,DIMENSION(is_nbfaces) :: localH
40 | DOM_REAL,DIMENSION(3,is_ncells) :: localQr
41 | DOM_REAL,DIMENSION(is_nfacesmax) :: Le
42 |
43 | ! ----------------------------------------------!
44 | ! Initialise vectors and rounding error catcher !
45 | ! ----------------------------------------------!
46 |
47 | H = 0.
48 | G = 0.
49 | Qr = 0.
50 | nmin = 0.
51 |
52 | localH = 0.
53 | localG = 0.
54 | localQr = 0.
55 |
56 | ! ----------------------------------!
57 | ! Looping over all slave directions !
58 | ! ----------------------------------!
59 |
60 | DO idir=1,is_ndir
61 |
62 | ! PRINT*, "dir -------", idir,is_dird, is_dirf
63 | ! IF (pmm_rank.eq.0) print*, "--- dir:", &
64 | ! & s_mu(idir),s_eta(idir),s_ksi(idir)
65 | ! IF (pmm_rank.eq.0) print*, " pathw:", is_pathway(:,idir)
66 |
67 | DO icell=1, is_ncells
68 |
69 | ibnd = 1
70 | j = is_pathway(icell,idir)
71 |
72 | ! IF (pmm_rank.eq.0) print*, " ++ Cell",icell,"/",is_ncells, &
73 | ! & ":",j
74 |
75 | ! --------------------------------------------!
76 | ! Calculate the direction.normal at each face !
77 | ! --------------------------------------------!
78 |
79 | DO i=1,is_nfcelt(j)
80 |
81 | Dm_ij(i) = s_norm(1,i,j)*s_mu(idir) + &
82 | & s_norm(2,i,j)*s_eta(idir) + &
83 | & s_norm(3,i,j)*s_ksi(idir)
84 | ENDDO
85 |
86 | ! ----------------------------------------------!
87 | ! Calculate intensity at exit faces of the cell !
88 | ! ----------------------------------------------!
89 |
90 | SELECT CASE (I_SCHEME)
91 |
92 | ! -----------------!
93 | ! MEAN FLUX SCHEME !
94 | ! -----------------!
95 | CASE(1)
96 | ! print*," (",pmm_rank,") Doing DMFS spatial discretization"
97 |
98 | CALL MFSCHEME(k_abs(j), k_scat(j), s_V(j), is_nfcelt(j), &
99 | & s_S(:,j), Dm_ij, Lb(j), Gi(j), Li(:,j), &
100 | & Le(:), Lpi, alpha, is_nfacesmax)
101 |
102 | ! -------------------!
103 | ! EXPONENTIAL SCHEME !
104 | ! -------------------!
105 | CASE(2)
106 | ! print*," (",pmm_rank,") Doing EXP spatial discretization"
107 |
108 | CALL EXPOSCHEME(k_abs(j), k_scat(j), s_V(j), &
109 | & is_nfcelt(j), s_S(:,j), Dm_ij, Lb(j), &
110 | & s_maxlen(icell,idir), s_ss(:,icell,idir), &
111 | & Gi(j), Li(:,j), Le(:), Lpi,is_nfacesmax)
112 |
113 | END SELECT
114 |
115 |
116 | localG(j) = localG(j) + Lpi*s_w(idir)
117 |
118 | localQr(1,j) = localQr(1,j) + Lpi*s_w(idir)*s_mu (idir)
119 | localQr(2,j) = localQr(2,j) + Lpi*s_w(idir)*s_eta(idir)
120 | localQr(3,j) = localQr(3,j) + Lpi*s_w(idir)*s_ksi(idir)
121 |
122 | ! IF (pmm_rank.eq.0) print*, " Ip =",Lpi
123 |
124 | ! -----------------------------------------------------------!
125 | ! Update incident intensities at the wall and neighbor cells !
126 | ! -----------------------------------------------------------!
127 |
128 | DO i=1,is_nfcelt(j)
129 |
130 | IF (Dm_ij(i).ge.nmin) THEN
131 |
132 | ! ----------------------------------!
133 | ! Detect if face is in the boundary !
134 | ! ----------------------------------!
135 |
136 | neighcellid = is_neighs(2*i-1,j)
137 |
138 | ! ----------------------------------------------!
139 | ! Apply luminance to the neighbor/boundary face !
140 | ! ----------------------------------------------!
141 |
142 | IF (neighcellid.ne.0) THEN
143 |
144 | neighfaceid = is_neighs(2*i,j)
145 | Li(neighfaceid, neighcellid) = Le(i)
146 |
147 | ELSE
148 |
149 | ! IF (pmm_rank.eq.0) print*, " face:", i
150 | CALL FIND_IBND(j,i,ibnd)
151 | localH(ibnd)= localH(ibnd) + Le(i)*s_w(idir)*Dm_ij(i)
152 |
153 | ENDIF
154 |
155 | ENDIF
156 |
157 | ENDDO
158 | ENDDO
159 | ENDDO
160 |
161 | ! -------------------------------------!
162 | ! Making sommation over ALL directions !
163 | ! -------------------------------------!
164 |
165 | #ifdef USEPALM
166 | CALL MPI_ALLREDUCE(localG, G, is_ncells, &
167 | & MPI_DOUBLE_PRECISION, MPI_SUM, PL_COMM_EXEC , &
168 | & ierr)
169 |
170 | CALL MPI_ALLREDUCE(localH, H, is_nbfaces, &
171 | & MPI_DOUBLE_PRECISION, MPI_SUM, PL_COMM_EXEC , &
172 | & ierr)
173 |
174 | CALL MPI_ALLREDUCE(localQr, Qr, 3*is_ncells, &
175 | & MPI_DOUBLE_PRECISION, MPI_SUM, PL_COMM_EXEC , &
176 | & ierr)
177 | #else
178 | CALL MPI_ALLREDUCE(localG, G, is_ncells, &
179 | & MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, &
180 | & ierr)
181 |
182 | CALL MPI_ALLREDUCE(localH, H, is_nbfaces, &
183 | & MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, &
184 | & ierr)
185 |
186 | CALL MPI_ALLREDUCE(localQr, Qr, 3*is_ncells, &
187 | & MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, &
188 | & ierr)
189 | #endif
190 |
191 | END SUBROUTINE SPATIAL_SCHEME
spatial_scheme.F could be called by: