1 | include(dom.inc)
2 |
3 | SUBROUTINE TAB_CASE(celldata, kabs, wquad, allkabs, DYH, DYC, &
4 | & DYCO, nYH, nYC, nYCO, DT, nT, nkabs, nelmts, &
5 | & nallbandes, all_WVNB, all_DWVNB, Lb, ieltd, &
6 | & ieltf)
7 |
8 | USE mod_inout
9 | USE mod_pmm
10 | #ifdef USEPALM
11 | USE palmlib
12 | #endif
13 |
14 | IMPLICIT NONE
15 |
16 | include 'dom_constants.h'
17 |
18 | DOM_INT :: nYH, nYC, nYCO, nT
19 | DOM_INT :: nkabs, nelmts, ierr
20 | DOM_INT :: i_bande,nallbandes
21 | DOM_REAL :: DYH, DYC, DYCO, DT
22 | DOM_REAL :: DYH2,DYC2,DYCO2, ave
23 | DOM_INT :: ieltd,ieltf,n_cl
24 | DOM_INT :: iq,ielt,i
25 | DOM_REAL, DIMENSION(8,nelmts) :: celldata
26 | DOM_INT, DIMENSION(8) :: bornes
27 | DOM_REAL, DIMENSION(nkabs,nelmts) :: kabs, local_kabs
28 | DOM_REAL, DIMENSION(nYH,nYC,nYCO,nT,nkabs) :: allkabs
29 | DOM_REAL, DIMENSION(16) :: d
30 | DOM_REAL, DIMENSION(4) :: dist
31 |
32 | DOM_REAL,DIMENSION(nallbandes) :: all_WVNB,all_DWVNB
33 | DOM_REAL,DIMENSION(nelmts) :: Lb
34 | DOM_REAL,DIMENSION(nelmts) :: K_SOOT, FSK_SOOT
35 | DOM_REAL :: blae, planck
36 | DOM_REAL :: WVNB_SI,DWVNB_SI
37 | DOM_REAL,DIMENSION(nallbandes) :: F
38 |
39 | DOM_REAL, PARAMETER :: x_min=0.0
40 | DOM_REAL, PARAMETER :: x_max=1.0
41 | DOM_REAL, dimension(nkabs) :: x_pts,wquad
42 |
43 | CALL gauleg(x_min,x_max,x_pts,wquad,nkabs)
44 |
45 | local_kabs = 0.
46 | kabs = 0.
47 | FSK_SOOT = 0.
48 |
49 | DYH2 = DYH*(DT/DYH)
50 | DYC2 = DYC*(DT/DYC)
51 | DYCO2 = DYCO*(DT/DYCO)
52 |
53 | ! ---------------------------------------------!
54 | ! Non-homogeneous system or homogeneous system !
55 | ! ---------------------------------------------!
56 |
57 | IF (homosyst.eq.'NO') THEN
58 | n_cl=ieltf
59 | ELSEIF (homosyst.eq.'YES') THEN
60 | n_cl=ieltd
61 | ENDIF
62 |
63 | DO ielt=ieltd,n_cl
64 |
65 | ! ---------------------------------!
66 | ! Mean soot absorption calculation !
67 | ! ---------------------------------!
68 |
69 | DO i_bande=1,nallbandes
70 |
71 | WVNB_SI = 100.*all_WVNB(i_bande)
72 | DWVNB_SI= 100.*all_DWVNB(i_bande)
73 |
74 | F(i_bande)=blae(WVNB_SI,celldata(1,ielt))/(pi*Lb(ielt))* &
75 | & DWVNB_SI
76 |
77 | K_SOOT(ielt)=WVNB_SI*5.5*celldata(8,ielt)
78 |
79 | FSK_SOOT(ielt)=FSK_SOOT(ielt)+F(i_bande)*K_SOOT(ielt)
80 |
81 | ENDDO
82 |
83 | FSK_SOOT(ielt) = FSK_SOOT(ielt) / SUM(F)
84 |
85 | ! ------------------------------------------------!
86 | ! Calculationn of the closest points in the table !
87 | ! ------------------------------------------------!
88 |
89 | IF (celldata(1,ielt).lt.300.d0) THEN
90 | bornes(1) = 1
91 | dist(1) = 0.d0
92 | ELSEIF (celldata(1,ielt).gt.2900.d0) THEN
93 | bornes(1) = int((2900.d0 - 300.d0)/DT) +1
94 | dist(1) = 0
95 | ELSE
96 | bornes(1) = int((celldata(1,ielt)-300.d0)/DT)+1
97 | dist(1) = (celldata(1,ielt)-300.d0) - ((bornes(1)-1)*1.0)*DT
98 | ENDIF
99 |
100 | IF ((modulo(celldata(1,ielt)-300.d0,DT).eq.0.).or. &
101 | & (celldata(1,ielt).gt.2900.d0)) THEN
102 | bornes(2) = bornes(1)
103 | ELSE
104 | bornes(2) = bornes(1)+1
105 | ENDIF
106 |
107 | ! ------------------------------------------!
108 | ! Test on species if no absorption products !
109 | ! ------------------------------------------!
110 |
111 | IF (celldata(3,ielt).le.DYH) THEN
112 | bornes(3)=1
113 | dist(2)=0.d0
114 | ELSE
115 | bornes(3) = int(celldata(3,ielt)/DYH)+1
116 | dist(2) = (celldata(3,ielt)-((bornes(3)-1)*1.0)*DYH)*(DT/DYH)
117 | ENDIF
118 |
119 | IF (celldata(4,ielt).le.DYC) THEN
120 | bornes(4)=1
121 | dist(3)=0.d0
122 | ELSE
123 | bornes(4) = int(celldata(4,ielt)/DYC)+1
124 | dist(3) = (celldata(4,ielt)-((bornes(4)-1)*1.0)*DYC)*(DT/DYC)
125 | ENDIF
126 |
127 | IF (celldata(5,ielt).le.DYCO) THEN
128 | bornes(5)=1
129 | dist(4)=0.d0
130 | ELSE
131 | bornes(5) = int(celldata(5,ielt)/DYCO)+1
132 | dist(4) =(celldata(5,ielt)-((bornes(5)-1)*1.0)*DYCO)*(DT/DYCO)
133 | ENDIF
134 |
135 | ! bornes(3) = int(celldata(3,ielt)/DYH)+1
136 | ! bornes(4) = int(celldata(4,ielt)/DYC)+1
137 | ! bornes(5) = int(celldata(5,ielt)/DYCO)+1
138 |
139 | IF( modulo(celldata(3,ielt),DYH).eq.0.) THEN
140 | bornes(6) = bornes(3)
141 | ELSE
142 | bornes(6) = bornes(3)+1
143 | ENDIF
144 |
145 | IF( modulo(celldata(4,ielt),DYC).eq.0.) THEN
146 | bornes(7) = bornes(4)
147 | ELSE
148 | bornes(7) = bornes(4)+1
149 | ENDIF
150 |
151 | IF( modulo(celldata(5,ielt),DYCO).eq.0.) THEN
152 | bornes(8) = bornes(5)
153 | ELSE
154 | bornes(8) = bornes(5)+1
155 | ENDIF
156 |
157 | ! ------------------------------------------------!
158 | ! Calculationn of the distances for interpolation !
159 | ! ------------------------------------------------!
160 | ! dist(2) = (celldata(3,ielt) - ((bornes(3)-1)*1.0)*DYH)*(DT/DYH)
161 | ! dist(3) = (celldata(4,ielt) - ((bornes(4)-1)*1.0)*DYC)*(DT/DYC)
162 | ! dist(4) = (celldata(5,ielt)- ((bornes(5)-1)*1.0) &
163 | ! & *DYCO)*(DT/DYCO)
164 |
165 | d(1) = sqrt(dist(1)**2 + dist (2)**2 + &
166 | & dist(3)**2 + 0.04*(dist(4)**2))
167 |
168 | d(2) = sqrt((DT - dist(1))**2 + dist(2)**2 + &
169 | & dist(3)**2 + 0.04*(dist(4)**2))
170 |
171 | d(3) = sqrt(dist(1)**2 + (DYH2 - dist(2))**2 + &
172 | & dist(3)**2 + 0.04*(dist(4)**2))
173 |
174 | d(4) = sqrt(dist(1)**2 + dist(2)**2 + &
175 | & (DYC2 - dist(3))**2 + 0.04*(dist(4)**2))
176 |
177 | d(5) = sqrt(dist(1)**2 + dist(2)**2 + &
178 | & dist(3)**2 + 0.04*(DYCO2 - dist(4))**2)
179 |
180 | d(6) = sqrt((DT - dist(1))**2 + &
181 | & (DYH2 - dist(2))**2 + dist(3)**2 + 0.04*(dist(4)**2))
182 |
183 | d(7) = sqrt((DT - dist(1))**2 + dist(2)**2 + &
184 | & (DYC2 - dist(3))**2 + 0.04*(dist(4)**2))
185 |
186 | d(8) = sqrt((DT - dist(1))**2 + dist(2)**2 + &
187 | & dist(3)**2 + 0.04*(DYCO2 - dist(4))**2)
188 |
189 | d(9) = sqrt(dist(1)**2 + (DYH2 - dist(2))**2 + &
190 | & (DYC2 - dist(3))**2 + 0.04*(dist(4)**2))
191 |
192 | d(10) = sqrt(dist(1)**2 + (DYH2 - dist(2))**2 + &
193 | & dist(3)**2 + 0.04*(DYCO2 - dist(4))**2)
194 |
195 | d(11) = sqrt(dist(1)**2 + dist(2)**2 + &
196 | & (DYC2 - dist(3))**2 + 0.04*(DYCO2 - dist(4))**2)
197 |
198 | d(16) = sqrt((DT - dist(1))**2 + (DYH2- dist(2))**2 + (DYC2 - &
199 | & dist(3))**2 + 0.04*(DYCO2 - dist(4))**2)
200 |
201 | d(12) = sqrt((DT - dist(1))**2 + (DYH2 - dist(2))**2 + &
202 | & (DYC2 - dist(3))**2 + 0.04*(dist(4)**2))
203 |
204 | d(13) = sqrt((DT - dist(1))**2 + (DYH2 - dist(2))**2 + &
205 | & dist(3)**2 + 0.04*(DYCO2 - dist(4))**2)
206 |
207 | d(14) = sqrt((DT - dist(1))**2 + dist(2)**2 + &
208 | & (DYC2 - dist(3))**2 + 0.04*(DYCO2 - dist(4))**2)
209 |
210 | d(15) = sqrt(dist(1)**2 + (DYH2 - dist(2))**2 + &
211 | & (DYC2 - dist(3))**2 + 0.04*(DYCO2 - dist(4))**2)
212 |
213 | ave = SUM(d(:))/16.d0
214 |
215 | ! -------------------------------------------------------!
216 | ! Interpolation of the tabulated absorption coefficients !
217 | ! -------------------------------------------------------!
218 |
219 | DO iq=1,nkabs
220 | local_kabs(iq,ielt) = &
221 | & allkabs(bornes(3),bornes(4),bornes(5), &
222 | & bornes(1),iq)*(ave/d(1)) + &
223 | & allkabs(bornes(3),bornes(4),bornes(5), &
224 | & bornes(2),iq)*(ave/d(2)) + &
225 | & allkabs(bornes(6),bornes(4),bornes(5), &
226 | & bornes(1),iq)*(ave/d(3)) + &
227 | & allkabs(bornes(3),bornes(7),bornes(5), &
228 | & bornes(1),iq)*(ave/d(4)) + &
229 | & allkabs(bornes(3),bornes(4),bornes(8), &
230 | & bornes(1),iq)*(ave/d(5)) + &
231 | & allkabs(bornes(6),bornes(4),bornes(5), &
232 | & bornes(2),iq)*(ave/d(6)) + &
233 | & allkabs(bornes(3),bornes(7),bornes(5), &
234 | & bornes(2),iq)*(ave/d(7)) + &
235 | & allkabs(bornes(3),bornes(4),bornes(8), &
236 | & bornes(2),iq)*(ave/d(8)) + &
237 | & allkabs(bornes(6),bornes(7),bornes(5), &
238 | & bornes(1),iq)*(ave/d(9)) + &
239 | & allkabs(bornes(6),bornes(5),bornes(8), &
240 | & bornes(1),iq)*(ave/d(10)) + &
241 | & allkabs(bornes(3),bornes(7),bornes(8), &
242 | & bornes(1),iq)*(ave/d(11)) + &
243 | & allkabs(bornes(6),bornes(7),bornes(5), &
244 | & bornes(2),iq)*(ave/d(12)) + &
245 | & allkabs(bornes(6),bornes(4),bornes(8), &
246 | & bornes(2),iq)*(ave/d(13)) + &
247 | & allkabs(bornes(3),bornes(7),bornes(8), &
248 | & bornes(2),iq)*(ave/d(14)) + &
249 | & allkabs(bornes(6),bornes(7),bornes(8), &
250 | & bornes(1),iq)*(ave/d(15)) + &
251 | & allkabs(bornes(6),bornes(7),bornes(8), &
252 | & bornes(2),iq)*(ave/d(16))
253 | ENDDO
254 |
255 | IF (celldata(1,ielt).lt.300.d0) THEN
256 | local_kabs(:,ielt) = 0.0d0
257 | ELSE
258 | local_kabs(:,ielt) = local_kabs(:,ielt)/16.d0 + FSK_SOOT(ielt)
259 | ENDIF
260 |
261 | ENDDO
262 |
263 | IF (homosyst.eq.'YES') THEN
264 | DO i=ieltd+1,ieltf
265 | local_kabs(:,i) = local_kabs(:,ieltd)
266 | ENDDO
267 | ENDIF
268 |
269 | ! ------------------------------------------------!
270 | ! sending all_k_abs....change to a pmm subroutine !
271 | ! ------------------------------------------------!
272 |
273 | #ifdef USEPALM
274 | CALL MPI_ALLREDUCE(local_kabs, kabs, nkabs*nelmts, &
275 | & MPI_DOUBLE_PRECISION, MPI_SUM, PL_COMM_EXEC ,&
276 | & ierr)
277 | #else
278 | CALL MPI_ALLREDUCE(local_kabs, kabs, nkabs*nelmts, &
279 | & MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD,&
280 | & ierr)
281 | #endif
282 |
283 | END SUBROUTINE TAB_CASE