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