1 | include(dom.inc)
2 |
3 | PROGRAM ray
4 |
5 | ! ==================================================================!
6 | ! !
7 | ! ray.F : Program that calculates the luminance that exits a !
8 | ! 1D gaz column of length "e" using different spectral !
9 | ! methods !
10 | ! !
11 | ! author : J.AMAYA (April 2008) !
12 | ! !
13 | ! ==================================================================!
14 |
15 | IMPLICIT NONE
16 |
17 | include 'dom_constants.h'
18 |
19 | DOM_REAL :: e,t_start,t_stop,time
20 | DOM_REAL :: Tw, Ka0, I0, L0, Ltot, Lw, attenuation
21 | DOM_REAL :: ls, kabsmoy, Twallout, Lwallout
22 | DOM_REAL :: Kn, L, Lq, somme, Lnus, Qw
23 | DOM_REAL :: blae, x, planck
24 | DOM_INT :: n_segments, ngbands, j, s, i, k, n
25 | DOM_INT :: nbands, nquad, nintbands
26 | CHARACTER*3 :: homosyst
27 | CHARACTER*15 :: mediumtype
28 |
29 | DOM_REAL, allocatable, dimension(:,:) :: celldata,tau
30 | DOM_REAL, allocatable, dimension(:) :: tempe, Lb, press
31 | DOM_REAL, allocatable, dimension(:) :: XH2O, XCO2, XCO
32 | DOM_REAL, allocatable, dimension(:) :: XO2, XN2, XSOOT
33 | DOM_REAL, allocatable, dimension(:) :: all_WVNB, all_DWVNB
34 | DOM_REAL, allocatable, dimension(:) :: deltanu
35 | DOM_REAL, allocatable, dimension(:) :: wquad
36 | DOM_REAL, allocatable, dimension(:,:,:) :: kabs
37 | DOM_REAL, allocatable, dimension(:,:) :: KCO,DCO,KC,DC,KH,DH
38 |
39 | ! Tabulation variables
40 | DOM_INT :: nYH, nYC, nYCO, nT
41 | DOM_REAL :: DYH, DYC, DYCO, DT
42 | DOM_REAL :: MH2O, MCO2, MCO
43 | DOM_REAL, ALLOCATABLE,DIMENSION(:,:,:,:,:) :: tabkabs
44 |
45 | DOM_INT ::cnttemp
46 | character*64 ::outfile,tabfile
47 | character*80 ::pathspec
48 |
49 | ! ------------------!
50 | ! Read choices file !
51 | ! ------------------!
52 |
53 | call CPU_TIME(t_start)
54 |
55 | WRITE(*,*) " >> Reading choices file"
56 | OPEN (FILE_CHCS , FILE='ray.choices', FORM='FORMATTED')
57 | READ(FILE_CHCS,*) pathspec
58 | READ(FILE_CHCS,*) mediumtype
59 | READ(FILE_CHCS,*) homosyst
60 | READ(FILE_CHCS,*) e
61 | READ(FILE_CHCS,*) n_segments
62 | READ(FILE_CHCS,*) nquad
63 | READ(FILE_CHCS,*) Tw
64 | READ(FILE_CHCS,*) Twallout
65 | CLOSE(FILE_CHCS)
66 |
67 | IF (mediumtype.eq.'TAB') THEN
68 | deltanu = 1.
69 | tabfile = pathspec(1:len_trim(pathspec))//'/table.dat'
70 | PRINT*, "Tabulated File : ",tabfile
71 |
72 | OPEN(12,FILE=tabfile,FORM='FORMATTED')
73 | READ(12,*) nquad
74 | READ(12,*) DT
75 | READ(12,*) DYH, MH2O
76 | READ(12,*) DYC, MCO2
77 | READ(12,*) DYCO, MCO
78 | CLOSE(12)
79 |
80 | nT=INT((2900-300)/DT)+1
81 | nYH=INT((MH2O-0.0)/DYH)+1
82 | nYC=INT((MCO2-0.0)/DYC)+1
83 | nYCO=INT((MCO-0.0)/DYCO)+1
84 | ALLOCATE(tabkabs(nYH,nYC,nYCO,nT,nquad))
85 | CALL readtab(tabfile,tabkabs,nYH,nYC,nYCO,nT,nquad)
86 | PRINT*, "Table for FSK model readed"
87 | ENDIF
88 |
89 | ! ---------------------!
90 | ! Initialize variables !
91 | ! ---------------------!
92 |
93 | WRITE(*,*) " >> Initializing data and allocating memory"
94 | ls = e / n_segments
95 | nbands = 371
96 |
97 | IF (ALLOCATED(celldata)) DEALLOCATE(celldata)
98 | IF (ALLOCATED(Lb)) DEALLOCATE(Lb)
99 | IF (ALLOCATED(press)) DEALLOCATE(press)
100 | IF (ALLOCATED(XH2O)) DEALLOCATE(XH2O)
101 | IF (ALLOCATED(XCO)) DEALLOCATE(XCO)
102 | IF (ALLOCATED(XCO2)) DEALLOCATE(XCO2)
103 | IF (ALLOCATED(XO2)) DEALLOCATE(XO2)
104 | IF (ALLOCATED(XN2)) DEALLOCATE(XN2)
105 | IF (ALLOCATED(XSOOT)) DEALLOCATE(XSOOT)
106 | IF (ALLOCATED(tempe)) DEALLOCATE(tempe)
107 | IF (ALLOCATED(all_WVNB)) DEALLOCATE(all_WVNB)
108 | IF (ALLOCATED(all_DWVNB)) DEALLOCATE(all_DWVNB)
109 | IF (ALLOCATED(KCO)) DEALLOCATE(KCO)
110 | IF (ALLOCATED(DCO )) DEALLOCATE(DCO)
111 | IF (ALLOCATED(KC)) DEALLOCATE(KC)
112 | IF (ALLOCATED(DC )) DEALLOCATE(DC)
113 | IF (ALLOCATED(KH)) DEALLOCATE(KH)
114 | IF (ALLOCATED(DH)) DEALLOCATE(DH)
115 | IF (ALLOCATED(wquad)) DEALLOCATE(wquad)
116 |
117 | allocate(celldata (n_segments,8))
118 | allocate(Lb (n_segments))
119 | allocate(press (n_segments))
120 | allocate(XH2O (n_segments))
121 | allocate(XCO2 (n_segments))
122 | allocate(XCO (n_segments))
123 | allocate(XO2 (n_segments))
124 | allocate(XN2 (n_segments))
125 | allocate(XSOOT (n_segments))
126 | allocate(tempe (n_segments))
127 | allocate(all_WVNB (nbands))
128 | allocate(all_DWVNB(nbands))
129 | allocate(KCO (14,nbands))
130 | allocate(DCO (14,nbands))
131 | allocate(KC (14,nbands))
132 | allocate(DC (14,nbands))
133 | allocate(KH (14,nbands))
134 | allocate(DH (14,nbands))
135 | allocate(wquad (nquad))
136 |
137 | ! --------------------!
138 | ! Temperature profile !
139 | ! --------------------!
140 |
141 | DO i = 1, n_segments
142 | x = i*ls
143 | include 'profile.inc'
144 | print*, i,":", tempe(i)
145 | ENDDO
146 |
147 | celldata(:,1) = tempe(:)
148 | celldata(:,2) = press(:)
149 | celldata(:,3) = XH2O(:)
150 | celldata(:,4) = XCO2(:)
151 | celldata(:,5) = XCO(:)
152 | celldata(:,6) = XO2(:)
153 | celldata(:,7) = XN2(:)
154 | celldata(:,8) = XSOOT(:)
155 |
156 | deallocate(press, XH2O, XCO2, XCO, XO2, XN2, XSOOT)
157 |
158 | ! ------------------------------!
159 | ! Read data from spectral files !
160 | ! ------------------------------!
161 |
162 | WRITE(*,*) " >> Reading spectral data"
163 | OPEN(UNIT=49,FILE='SPECTRAL/SNBWN')
164 |
165 | i = 1
166 | DO WHILE (i.le. nbands)
167 |
168 | READ(49,*) all_WVNB(i),all_DWVNB(i)
169 | IF (all_WVNB(i).lt.0.) i = nbands
170 | i = i + 1
171 |
172 | ENDDO
173 |
174 | CLOSE(49)
175 |
176 | call PARAM(KCO,KC,KH,DCO,DC,DH,pathspec)
177 |
178 | ! ------------------------!
179 | ! Select sommation limits !
180 | ! ------------------------!
181 |
182 | WRITE(*,*) " >> Selecting limits for the integration"
183 |
184 | IF (mediumtype.eq.'FSK') THEN
185 | ngbands = nbands
186 | nintbands = 1
187 | ELSEIF (mediumtype.eq.'CK') THEN
188 | ngbands = 1
189 | nintbands = nbands
190 | ELSEIF (mediumtype.eq.'MALKMUS') THEN
191 | ngbands = 1
192 | nintbands = nbands
193 | ELSEIF (mediumtype.eq.'TAB') THEN
194 | ngbands = nbands
195 | nintbands = 1
196 | ENDIF
197 |
198 | IF (ALLOCATED(kabs)) DEALLOCATE(kabs)
199 | IF (ALLOCATED(deltanu)) DEALLOCATE(deltanu)
200 |
201 | allocate(kabs (n_segments, nquad, nintbands))
202 | allocate(deltanu(nintbands))
203 |
204 | IF (mediumtype.eq.'FSK') THEN
205 | deltanu = 1.
206 | ELSEIF (mediumtype.eq.'CK') THEN
207 | deltanu(:) = all_DWVNB(:) * 100.
208 | ELSEIF (mediumtype.eq.'MALKMUS') THEN
209 | deltanu(:) = all_DWVNB(:) * 100.
210 | ELSE IF (mediumtype.eq.'TAB') THEN
211 | deltanu = 1.
212 | ENDIF
213 |
214 | ! ---------------------!
215 | ! Printing information !
216 | ! ---------------------!
217 |
218 | WRITE(*,*) " __// RAY v1.0 \\__ "
219 | WRITE(*,*)
220 | WRITE(*,*) " Information:"
221 | WRITE(*,*) " ------------"
222 | WRITE(*,*)
223 | WRITE(*,*) " Spectral model :", mediumtype
224 | WRITE(*,*) " Number of segments :", n_segments
225 | WRITE(*,*) " Quadrature points :", nquad
226 | WRITE(*,*) " Total number of bands:", nbands
227 | WRITE(*,*) " Total grouped bands :", ngbands
228 | WRITE(*,*) " Total integrations :", nintbands
229 | WRITE(*,*) " Gaz wall thickness :", e
230 | WRITE(*,*) " Gaz segment thickness:", ls
231 | WRITE(*,*) " Wall temperature :", Tw
232 | WRITE(*,*)
233 |
234 | IF (nquad.ge.10) THEN
235 | outfile=trim(mediumtype)//'_Nq='//CHAR(48+int(nquad/10)) &
236 | & //CHAR(48+modulo(nquad,10))//'_out.dat'
237 | ELSE
238 | outfile = trim(mediumtype)//'_Nq='//CHAR(48+nquad)//'_out.dat'
239 | ENDIF
240 |
241 | OPEN (UNIT=102 , FILE=outfile, FORM='FORMATTED')
242 |
243 | ! ------------------------!
244 | ! Calculate spectral data !
245 | ! ------------------------!
246 |
247 | WRITE(*,*) " >> Calculate absorption coefficients"
248 |
249 | IF (mediumtype.eq.'FSK') THEN
250 | Lb(:) = planck(celldata(:,1))
251 | call FSK_CASE(kabs(:,:,1), wquad, Lb, celldata, n_segments, &
252 | & all_WVNB, all_DWVNB, KCO,DCO,KC,DC,KH,DH, &
253 | & homosyst, ngbands, nbands, nquad)
254 |
255 | ELSEIF(mediumtype.eq.'TAB') THEN
256 | Lb(:) = planck(celldata(:,1))
257 | CALL TAB_CASE(celldata,kabs(:,:,1), wquad, tabkabs, DYH, &
258 | & DYC, DYCO, nYH, nYC, nYCO, DT, nT, nquad, &
259 | & n_segments, nbands, all_WVNB, all_DWVNB, Lb)
260 |
261 | ELSEIF (mediumtype.eq.'CK') THEN
262 | DO i = 1, nintbands
263 | ! WRITE(*,*) " Absorption coeff - Doing band: ",i,"/",nbands
264 | Lb(:) = blae(all_WVNB(i)*100.,celldata(:,1))/pi
265 | call CK_CASE (kabs(:,:,i), wquad, Lb, celldata, &
266 | & n_segments, all_WVNB, all_DWVNB, KCO, DCO, &
267 | & KC, DC, KH, DH, homosyst,ngbands,i, nbands, &
268 | & nquad)
269 | ENDDO
270 |
271 | ELSEIF (mediumtype.eq.'MALKMUS') THEN
272 | IF (ALLOCATED(tau)) DEALLOCATE(tau)
273 | ALLOCATE(tau(n_segments,nbands))
274 | DO i = 1, nintbands
275 | call MALKMUS_CASE(tau(:,i),ls,celldata,n_segments, &
276 | & all_WVNB, all_DWVNB,KCO,DCO,KC,DC,KH,DH, &
277 | & i,nbands)
278 | ENDDO
279 |
280 | ENDIF
281 |
282 | ! --------------------------------!
283 | ! Integrate over all wave lengths !
284 | ! --------------------------------!
285 |
286 | L = 0.
287 | k = 0
288 | kabsmoy = 0.
289 |
290 | WRITE(*,*) " >> Luminance calculation"
291 | DO i = 1, nintbands
292 |
293 | ! -------------------------------!
294 | ! Integrate over all quadratures !
295 | ! -------------------------------!
296 |
297 | WRITE(*,*) " Narrow band ",i,"/",nintbands
298 | Lq = 0.
299 |
300 | IF (mediumtype.eq.'MALKMUS') THEN
301 |
302 | somme = 0.
303 | DO s = 1, n_segments
304 | ! print*, " Segment: ", s,"/", n_segments
305 |
306 | Kn = 1.
307 | DO n = s+1, n_segments
308 | Kn = Kn * tau (n,i)
309 | ENDDO
310 |
311 | Lnus = blae(all_WVNB(i)*100.,celldata(s,1))/pi
312 | somme = somme + (1-tau(s,i))*Lnus*Kn
313 | ! print*, " - tau(s,i) :", tau(s,i)
314 | ! print*, " - Kn, Lnus :", Kn, Lnus
315 |
316 | ENDDO
317 |
318 | L = L + deltanu(i) * somme
319 | ! print*, " - L, deltanu, Lq : ", L, deltanu(i), Lq
320 |
321 | ELSE
322 |
323 | DO j = 1, nquad
324 |
325 | ! ----------------------------------------------!
326 | ! Emission for each segment of the optical path !
327 | ! ----------------------------------------------!
328 | ! print*, " Quadrature point: ", j,"/", nquad
329 |
330 | somme = 0.
331 | DO s = 1, n_segments
332 | ! print*, " Segment: ", s,"/", n_segments
333 |
334 | Kn = 0.
335 | DO n = s+1, n_segments
336 | Kn = Kn + kabs(n,j,i)*ls
337 | ! print*, " Next segment ", n,"/",n_segments
338 | ! print*, " - kabs(n,j,i):", kabs(n,j,i)
339 | ! print*, " - ls :", ls
340 | ENDDO
341 |
342 | IF (mediumtype.eq.'CK') THEN
343 | Lnus = blae(all_WVNB(i)*100.,celldata(s,1))/pi
344 | ELSE IF (mediumtype.eq.'FSK') THEN
345 | Lnus = planck(celldata(s,1))
346 | ENDIF
347 |
348 | somme = somme + (1-exp(-kabs(s,j,i)*ls))*Lnus*exp(-Kn)
349 |
350 | ! print*, " - kabs(s,j,i):",kabs(s,j,i)
351 | ! print*, " - Kn, Lnus :", Kn, Lnus
352 | ENDDO
353 |
354 | ! ---------------------------------------!
355 | ! Calculate absorption of wall intensity !
356 | ! ---------------------------------------!
357 |
358 | I0 = blae(all_WVNB(i)*100.,Tw)/pi
359 | Ka0 = 0.
360 |
361 | DO s = 1, n_segments
362 | Ka0 = Ka0 + kabs(s,j,i)*ls
363 | ENDDO
364 |
365 | L0 = I0 * exp( -Ka0 ) * wquad(j)
366 |
367 | ! -----------------------------------------!
368 | ! integrate this spectral quadrature point !
369 | ! -----------------------------------------!
370 |
371 | Lq = Lq + somme * wquad(j)
372 | kabsmoy = kabsmoy + wquad(j) * deltanu(i) * kabs(1,j,i) &
373 | & * blae(all_WVNB(i)*100.,celldata(1,1))/pi
374 |
375 | IF (mediumtype.eq.'FSK') k = k + 1
376 |
377 | ENDDO
378 |
379 | L = L + deltanu(i) * Lq
380 | Lw = Lw + deltanu(i) * L0
381 |
382 | ! print*, " - L, deltanu, Lq : ", L, deltanu(i), Lq
383 |
384 | IF (mediumtype.eq.'CK') k = k + 1
385 |
386 | ENDIF
387 |
388 | ENDDO
389 |
390 | kabsmoy = kabsmoy / planck(celldata(1,1))
391 | Ltot = Lw + L
392 | attenuation = Ltot / planck(Tw)
393 | Lwallout = planck(Twallout)
394 | Qw = Lwallout - Ltot
395 |
396 | WRITE(*,*) " Total wall luminance reaching the exit, Lw =",Lw
397 | WRITE(*,*) " Emited/self-absorbed luminance of the gas, L =",L
398 | WRITE(*,*) " Luminance at the exit, Ltot = Lw + L = ", Ltot
399 | WRITE(*,*) " Out wall total heat flux, Qw =", Qw
400 |
401 | WRITE(*,*)
402 | WRITE(*,*) " kabsmoy = ", kabsmoy
403 | WRITE(*,*)
404 |
405 | WRITE(102,*) celldata(n_segments,1), Lw, L, Ltot, &
406 | & attenuation, kabsmoy, Qw
407 |
408 |
409 | CLOSE(102)
410 | call CPU_TIME(t_stop)
411 |
412 | PRINT*,'Elapsed CPU time = ',t_stop-t_start
413 |
414 | outfile="TIME_"//outfile
415 | OPEN(103,file=outfile,form='formatted')
416 | WRITE(103,*) t_stop-t_start
417 | CLOSE(103)
418 |
419 | END PROGRAM ray