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