1 | program DOMASIUM
2 |
3 | !##############################################################################
4 | ! ACQUISITION DES PARAMETRES COMMUNS
5 | !##############################################################################
6 |
7 | USE paramax
8 |
9 | IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N)
10 |
11 | CHARACTER (len=60) :: nomfich
12 | CHARACTER (len=15) :: meshstatus,datastatus,meshtype,quadtype,&
13 | mediumtype,spascheme
14 |
15 | INTEGER, DIMENSION (nmax,2*nfacemax) :: neighs
16 |
17 | INTEGER, DIMENSION (nmax) :: nfcelt, nptelt, pathway
18 |
19 | REAL(8) :: kbar, phi_f, DWVNB
20 |
21 | REAL(8),DIMENSION(2) :: SNBDATA
22 |
23 | REAL(8),DIMENSION(n_SNBmax) :: WV
24 |
25 | REAL(8),DIMENSION (nkabs) :: wkabs, ka_g_i, k_absmel,WVSOOT
26 |
27 | REAL(8),DIMENSION(14,n_SNBmax) :: KCO,DCO,KC,DC,KH,DH
28 |
29 | REAL(8),DIMENSION(nmax,8) :: celldata
30 |
31 | REAL(8),DIMENSION(nmax,nkabs) :: all_k_abs
32 |
33 | REAL(8),DIMENSION(n_SNBmax) :: all_WVNB,all_DWVNB
34 |
35 | REAL(8),DIMENSION(nquadmax) :: mu, eta, ksi, w
36 |
37 | REAL(8),DIMENSION(nfacemax) :: Dm_ij
38 |
39 | REAL(8),DIMENSION(nmax) :: Lb, Lpi, Spi, Gi, Gj, G, Gtot, T, &
40 | xc, yc, zc, V, k_scat, k_abs, Sr, Srtot, maxlen, Qin, Qintot, &
41 | wk, all_a1, all_a2, Gtotpar, Srtotpar,kabs_gray
42 |
43 | REAL(8),DIMENSION(nmax,3) :: Q_rtot
44 |
45 | REAL(8), DIMENSION (nmax,nfacemax) :: Lo, Lotot, Li, Le, L, H, Htot, &
46 | S, epsil, Tf, s_s, xfc, yfc, zfc, Htotpar, Lototpar
47 |
48 | REAL(8), DIMENSION (nmax,nfacemax,3) :: norm
49 |
50 |
51 | REAL :: etime ! Declare the type of etime()
52 | REAL, DIMENSION (2) :: elapsed ! For receiving user and system time
53 | REAL :: totaltime ! For receiving total time
54 |
55 |
56 | INTEGER :: finprog
57 |
58 | !##############################################################################
59 | ! DECLARATION DE L'ENVIRONNEMENT PARALLELE
60 | !##############################################################################
61 |
62 |
63 | !##############################################################################
64 | ! ACQUISITION DES 1eres DONNEES
65 | !##############################################################################
66 | ! Lecture du fichier d'entree INPUT :
67 | ! meshtype : GAMBIT ou AVBP
68 | ! nomfich : Chemin d'acces au fichier complet
69 | ! Meshstatus : NEW ou OLD ( OPTION utile pour CREATELT.f90 uniquement )
70 | ! Quadtype : Type de Quadrature (SNDOM, TNDOM ou FVM)
71 | ! Spascsheme : Type de Schema de derivation spatial (MFS ou EXPON)
72 | ! critconv : Critere de convergence (utile en cas de Reflexion et/ou Diffusion)
73 | ! ncells : Nombre de cellules composant le maillage
74 | !##############################################################################
75 |
76 |
77 | CALL READ_DATA(nomfich,meshstatus,datastatus,meshtype,quadtype,&
78 | mediumtype,spascheme,neighs,nfcelt,KCO,DCO,KC,DC,KH,DH,&
79 | xc, yc, zc, V, k_scat, kabs_gray,S, epsil, xfc, yfc, zfc, Tf,&
80 | norm,celldata,mu, eta, ksi, w,critconv,alpha,ncells,I_SCHEME,&
81 | ndir,nbounds,end_prog)
82 |
83 | !###################################################################!
84 | !###################################################################!
85 | ! !
86 | ! FIN DE CHARGEMENT DES DONNEES !
87 | ! !
88 | !###################################################################!
89 | !###################################################################!
90 |
91 |
92 | finprog=0
93 |
94 |
95 | !###################################################################!
96 | !###################################################################!
97 | ! !
98 | ! INITIALISATION DES MATRICES RESULTATS !
99 | ! !
100 | !###################################################################!
101 | !###################################################################!
102 |
103 |
104 | Gtot=0.
105 | Htot=0.
106 | Lotot=0.
107 | Srtot=0.
108 | Q_rtot=0.
109 |
110 | !###################################################################!
111 | !###################################################################!
112 | ! !
113 | ! SPECTRAL TREATMENT !
114 | ! !
115 | !###################################################################!
116 | !###################################################################!
117 | ! !
118 | ! GRAY CASE !
119 | ! !
120 | !###################################################################!
121 | !###################################################################!
122 |
123 | IF (mediumtype=='GRAY') THEN
124 |
125 | CALL GRAY_CASE(mediumtype,nfcelt,wkabs,celldata,all_k_abs,Lb, &
126 | kabs_gray,Lo, epsil, Tf, ncells,n_k_abs,DWVNB_SI)
127 |
128 |
129 | CALL BAND_INTEG(neighs, nfcelt, pathway, wkabs, all_k_abs, mu, &
130 | eta, ksi, w, Dm_ij, Lb, Gtot, V, k_scat, maxlen,&
131 | Srtot, Q_rtot, Lo, Lotot, Htot, S, epsil, s_s, norm,&
132 | n_k_abs, ncells, I_SCHEME, ndir, DWVNB_SI, alpha)
133 |
134 |
135 | ELSE
136 |
137 |
138 | !###################################################################!
139 | !###################################################################!
140 | ! !
141 | ! NON-GRAY CASE !
142 | ! !
143 | !###################################################################!
144 | !###################################################################!
145 |
146 |
147 | OPEN(UNIT=49,FILE='SPECTRAL/SNBWN')
148 |
149 | i_bande=1
150 | DO WHILE (i_bande .le. n_SNBmax)
151 |
152 | READ(49,*) all_WVNB(i_bande),all_DWVNB(i_bande)
153 | IF (all_WVNB(i_bande).LT.0.) THEN
154 | i_bande=n_SNBmax !arret de la boucle
155 | ENDIF
156 | i_bande=i_bande+1
157 | ENDDO
158 |
159 |
160 | CLOSE(49)
161 |
162 |
163 | DO i_bande=1,n_SNBmax
164 | ! DO i_bande=1,10
165 |
166 | WVNB=all_WVNB(i_bande)
167 | DWVNB=all_DWVNB(i_bande)
168 |
169 | PRINT*, 'Bande: ',i_bande,'/'
170 |
171 | WVSOOT=100.*WVNB*5.5*XSOOT
172 |
173 | all_k_abs=0.
174 | k_absmel=0.
175 |
176 | !###################################################################!
177 | !###################################################################!
178 | ! !
179 | ! CALCULS DES INDICES SPECTRAUX POUR CHAQUE ESPECES !
180 | ! !
181 | !###################################################################!
182 | !###################################################################!
183 |
184 | IF (WVNB>9300.0) THEN
185 | k_absmel=0.0
186 | wkabs=1.0/nkabs
187 | ELSE
188 | CALL FINDI(LICO,LICO2,LIH2O,ICO,ICO2,IH2O,C1,C2,WVNB,DWVNB)
189 | ENDIF
190 |
191 |
192 | !###################################################################!
193 |
194 | !###################################################################!
195 | ! !
196 | ! INITIALISATION DES MATRICES ET DU COMPTEUR D'ITERATION !
197 | ! !
198 | !###################################################################!
199 | !###################################################################!
200 |
201 |
202 | WVNB_SI=WVNB*100.
203 | DWVNB_SI=DWVNB*100.
204 |
205 | !###################################################################!
206 | !###################################################################!
207 | ! !
208 | ! CALCUL DE TOUTES LES LUMINANCES NOIRES ET DE TOUS {K_ij,w_ij} !
209 | ! DU MAILLAGE POUR LA BANDE ETROITE CONSIDEREE !
210 | ! !
211 | !###################################################################!
212 | !###################################################################!
213 |
214 | IF (datastatus=='EXTERNAL') THEN
215 |
216 |
217 | CALL EXTERN(nfcelt,wkabs,celldata,all_k_abs,kabs_gray,Lo, epsil,&
218 | ncells,WVNB_SI,n_k_abs,Tf,Lb)
219 |
220 | ELSE
221 |
222 | IF (homosyst=='NO') THEN
223 |
224 | CALL NO_HOMOSYST(nfcelt,SNBDATA,wkabs, k_absmel,KCO,DCO,KC,DC,KH,DH,&
225 | celldata,Lb,all_a1,all_a2,Lo,epsil,Tf,ncells,WVNB_SI,&
226 | WVNB,LICO,LICO2,LIH2O,ICO,ICO2,IH2O,all_k_abs,WVSOOT)
227 |
228 |
229 | ELSE IF (homosyst=='YES') THEN
230 |
231 | CALL YES_HOMOSYST(nfcelt,SNBDATA,wkabs,k_absmel,KCO,DCO,KC,DC,KH,DH,&
232 | celldata,Lb,Lo, epsil,Tf,WVSOOT,all_k_abs,ncells,&
233 | WVNB_SI,LICO,LICO2,LIH2O,ICO,ICO2,IH2O)
234 |
235 | ELSE
236 |
237 | finprog=1
238 | ENDIF
239 |
240 | n_k_abs=nkabs
241 |
242 | ENDIF !IF (datastatus=='EXTERNAL') THEN
243 |
244 |
245 | !###################################################################!
246 | !###################################################################!
247 | ! !
248 | ! INTEGRATION SUR LA BANDE ETROITE PAR SOMMATION SUR LES g_ij !
249 | ! !
250 | !###################################################################!
251 | !###################################################################!
252 |
253 | IF (finprog .ne. 1) THEN
254 |
255 | CALL BAND_INTEG(neighs, nfcelt, pathway, wkabs, all_k_abs, mu, &
256 | eta, ksi, w, Dm_ij, Lb, Gtot, V, k_scat, maxlen,&
257 | Srtot, Q_rtot, Lo, Lotot, Htot, S, epsil, s_s, norm,&
258 | n_k_abs, ncells, I_SCHEME, ndir, DWVNB_SI, alpha)
259 | ENDIF
260 |
261 | ENDDO
262 |
263 | ENDIF !IF (mediumtype=='GRAY') THEN
264 |
265 |
266 | IF (finprog .ne. 1) THEN
267 |
268 |
269 | totaltime = etime(elapsed)
270 |
271 | print*,'End total=',totaltime,' user=',elapsed(1),&
272 | ' system=',elapsed(2)
273 |
274 | WRITE(1,*)
275 | WRITE(1,*) 'That s all Folks !'
276 | WRITE(1,*)
277 | WRITE(1,*)
278 |
279 | CLOSE(1)
280 |
281 | CALL OUTPROCESSING(ncells,nbounds,xc,yc,zc,Gtot,Srtot,Htot,&
282 | Lotot,Q_rtot,Lb,epsil)
283 |
284 | ENDIF
285 |
286 | END program DOMASIUM
287 |
288 |
289 |
290 |
291 | !********************************************************************
292 | !********************************************************************
293 | !********************************************************************
294 | !********************************************************************
295 | !********************************************************************
296 | !********************************************************************
297 | !********************************************************************
298 | !********************************************************************
299 | !********************************************************************
300 | !********************************************************************
301 | !********************************************************************
302 | !********************************************************************
303 | !*** ****
304 | !*** **** * * **** **** *** * * ***** * * * ***** ****
305 | !*** * * * * * * * * * * * * * ** * * ****
306 | !*** *** * * **** **** * * * * * * * * * *** ****
307 | !*** * * * * * * * * * * * * * * ** * ****
308 | !*** **** *** **** * * *** *** * * * * ***** ****
309 | !*** ****
310 | !********************************************************************
311 | !********************************************************************
312 | !********************************************************************
313 | !********************************************************************
314 | !********************************************************************
315 | !********************************************************************
316 | !********************************************************************
317 | !********************************************************************
318 | !********************************************************************
319 | !********************************************************************
320 | !********************************************************************
321 | !********************************************************************
322 |
323 |
324 |
325 |
326 | SUBROUTINE MFSCHEME(k_abs,w_kabs,k_scat,V,nface,S,Di,Lb,Gi,Li,epsil,&
327 | Le,Lpi,alpha)
328 |
329 | !********************************************************************
330 | ! kabs --> Coeff. d'absorption gris !
331 | ! wkabs --> Poids associe au numero du coeff !
332 | ! kscat --> Coeff de diffusion isotrope !
333 | ! V --> Volume de la cellule !
334 | ! S --> Vecteur des 4 surfaces (ABC, ABD, ACD et BCD) !
335 | ! Di --> Vecteur des produits (normes * direction discrete)!
336 | ! Lb --> Luminace noire emise au centre de la cellule !
337 | ! Gi --> Luminance incidente a l'iteration precedente !
338 | ! Li --> Lminance aux faces d'entree !
339 | ! epsil --> Emissivite( -1 pour les faces non parietales ) !
340 | ! Le --> Luminance calculees aux faces de sortie !
341 | ! Lpi --> Lpi total au centre de la cellule !
342 | !********************************************************************
343 |
344 | USE paramax
345 |
346 | IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N)
347 |
348 | REAL(8) :: L_IN,AREA_IN,k_abs,w_kabs
349 | REAL(8) :: Lb, Gi, k_scat, V, ko, Lpi
350 | REAL(8), DIMENSION (nfacemax) :: Li, Le, Di, S, epsil
351 |
352 | L_IN=0.
353 | AREA_IN=0.
354 |
355 | DO i=1,nface
356 | IF (Di(i)<0.) THEN
357 | L_IN=L_IN+S(i)*Di(i)*Li(i)
358 | AREA_IN=AREA_IN+S(i)*Di(i)
359 | ENDIF
360 | ENDDO
361 |
362 | Y1=-L_IN+alpha*V*(k_abs*Lb+k_scat/&
363 | (4*pi)*Gi)
364 | Y2=(k_abs+kscat)*alpha*V-AREA_IN
365 | Lpi=Y1/Y2
366 |
367 | DO ici1=1,nface
368 | IF (Di(ici1)>0.) THEN
369 | Le(ici1)=(Lpi-(1.-alpha)*L_IN/AREA_IN)/alpha
370 | ! IF (Le(ici1)<0) THEN
371 | ! PRINT*,Le(ici1)
372 | ! ENDIF
373 | ENDIF
374 | ENDDO
375 |
376 | END SUBROUTINE MFSCHEME
377 |
378 |
379 | !********************************************************************
380 | !********************************************************************
381 | !********************************************************************
382 | !********************************************************************
383 | !********************************************************************
384 | !********************************************************************
385 |
386 |
387 | SUBROUTINE EXPOSCHEME(k_abs,wkabs,k_scat,V,S,Di,maxlen,s_s,&
388 | Lb,Gi,Li,epsil,Le,X,X3)
389 |
390 | !******************************************************************
391 | ! k_abs --> Coeff. d'absorption gris !
392 | ! nkabs --> Nombre de coeff gris ( le + svt = 1 ) !
393 | ! wkabs --> Poids associe au numero du coeff !
394 | ! k_scat --> Coeff de diffusion isotrope !
395 | ! maxlen --> Longueur carac. du phenomene d'absorption !
396 | ! ntype --> Type de la cellule v-a-v. du la dir. discrete !
397 | ! s_s --> Vecteur des Coeff. d'echange surface/surface !
398 | ! V --> Volume de la cellule !
399 | ! S --> Vecteur des 4 surfaces (ABC, ABD, ACD et BCD) !
400 | ! Di --> Vecteur des Dij (normales*direction discrete) !
401 | ! Lb --> Luminace noire emise au centre de la cellule !
402 | ! Gi --> Luminance incidente a l'iteration precedente !
403 | ! Li --> Lminance aux faces d'entree !
404 | ! epsil --> Emissivite ( -1 pour les faces non parietales ) !
405 | ! Le --> Luminance calculees aux faces de sortie !
406 | ! X --> Lpi au centre de la cellule !
407 | ! Y0 --> Terme Source par Emission et diffusion incidente!
408 | ! Y1 --> Terme ne dependant que de l'epaisseur optique !
409 | !******************************************************************
410 |
411 | USE paramax
412 |
413 | IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N)
414 |
415 | REAL(8) :: Lb, Gi, Lpi, V, maxlen, KSI,&
416 | k_abs, wkabs, k_scat, omega, tau
417 | REAL(8), DIMENSION (nfacemax) :: Li, Le, Di, S, epsil, s_s, X2
418 |
419 | X=0.
420 | beta=k_scat+k_abs
421 | omega=k_scat/beta
422 | tau=beta*maxlen
423 | Y0=(1.-omega)*Lb+omega*Gi/(4*pi)
424 | KSI=(2./tau)*(1.-(1./tau)*(1.-exp(-tau)))
425 | Y2=KSI*sum(Li*s_s)+Y0*(1-KSI)
426 | X0=Y2/(beta*V)
427 | X2=Li/(beta*V)
428 | X1=Y0
429 |
430 | ! Sommation sur les 4 faces !
431 | DO ici1=1,4
432 | IF (Di(ici1)>0.) THEN
433 | X=X+X0*S(ici1)*Di(ici1) ! Faces de sorties
434 | IF (epsil(ici1)==-1.) THEN
435 | Le(ici1)=Y2 ! Faces parietales
436 | ENDIF
437 | ELSEIF (Di(ici1)<0.) THEN
438 | X=X+X2(ici1)*S(ici1)*Di(ici1) ! Faces d'entrees
439 | ENDIF
440 | ENDDO
441 |
442 | X=X1-X
443 |
444 | ! IF (X<0.) THEN
445 | ! PRINT*,'Negative Lpi for K=',k_abs
446 | ! ENDIF
447 |
448 |
449 | END SUBROUTINE EXPOSCHEME
450 |
451 |
452 | !***************************************************************************************
453 | !***************************************************************************************
454 | !***************************************************************************************
455 | !***************************************************************************************
456 | !***************************************************************************************
457 | !***************************************************************************************
458 |
459 |
460 | SUBROUTINE OUTPROCESSING(ncells,nbounds,xc,yc,zc,G,Sr,H,Lo,Q_r,Lb,epsil)
461 |
462 | USE paramax
463 |
464 | IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N)
465 |
466 | REAL(8) , DIMENSION (nfacemax) :: xfc,yfc,zfc
467 | REAL(8) , DIMENSION (nmax) :: G,xc,yc,zc,Lb,Sr
468 | REAL(8) , DIMENSION (nmax,nfacemax) :: H,Lo,epsil
469 | REAL(8) , DIMENSION (nmax,3) :: Q_r
470 | INTEGER , DIMENSION (nmax) :: nfcelt
471 |
472 |
473 | OPEN(UNIT=10,FILE='OUTFILES/G.out',FORM='unformatted')
474 | OPEN(UNIT=20,FILE='OUTFILES/Sr.out',FORM='unformatted')
475 | OPEN(UNIT=30,FILE='OUTFILES/H.out',FORM='unformatted')
476 | OPEN(UNIT=31,FILE='OUTFILES/Qw.out',FORM='unformatted')
477 | OPEN(UNIT=32,FILE='OUTFILES/Qr.out',FORM='unformatted')
478 | OPEN(UNIT=40,FILE='INFILES/Centerfaces.in',FORM='unformatted')
479 |
480 | WRITE(10) ncells
481 | WRITE(20) ncells
482 | WRITE(30) nbounds
483 | WRITE(31) nbounds
484 |
485 | DO i=1,ncells
486 |
487 | READ(40) icell,nfcelt(icell),((xfc(l),yfc(l),zfc(l)),l=1,nfcelt(icell))
488 | WRITE(10) xc(icell),yc(icell),zc(icell),G(icell)
489 | WRITE(20) xc(icell),yc(icell),zc(icell),Sr(icell)
490 | WRITE(32) xc(icell),yc(icell),zc(icell),(Q_r(icell,l),l=1,3)
491 |
492 | DO j=1,nfcelt(icell)
493 |
494 | IF (epsil(icell,j)/=-1.) THEN
495 |
496 | WRITE(30) xfc(j),yfc(j),zfc(j),H(icell,j)
497 |
498 | Qw=epsil(icell,j)*H(icell,j)-pi*Lo(icell,j)
499 |
500 | WRITE(31) xfc(j),yfc(j),zfc(j),Qw
501 |
502 | ENDIF
503 |
504 | ENDDO
505 |
506 | ENDDO
507 |
508 | CLOSE(10)
509 | CLOSE(20)
510 | CLOSE(30)
511 | CLOSE(31)
512 | CLOSE(32)
513 | CLOSE(40)
514 |
515 | END SUBROUTINE OUTPROCESSING
516 |
517 | !*****************************************************************
518 | !*****************************************************************
519 | !*****************************************************************
520 | !*****************************************************************
521 | !*****************************************************************
522 | ! *
523 | ! THIS SUBROUTINE READS THE MODEL PARAMETERS *
524 | ! *
525 | !*****************************************************************
526 | !
527 |
528 |
529 | SUBROUTINE PARAM(KCO,KC,KH,DCO,DC,DH)
530 |
531 | USE paramax
532 |
533 | REAL(8),DIMENSION(14,n_SNBmax)::KCO,DCO,KC,DC,KH,DH
534 |
535 | OPEN(UNIT=1010,FILE='SPECTRAL/SNBH2O_K')
536 | OPEN(UNIT=1020,FILE='SPECTRAL/SNBH2O_Phi')
537 | OPEN(UNIT=1030,FILE='SPECTRAL/SNBCO2_K')
538 | OPEN(UNIT=1040,FILE='SPECTRAL/SNBCO2_Phi')
539 | OPEN(UNIT=1050,FILE='SPECTRAL/SNBCO_K')
540 | OPEN(UNIT=1060,FILE='SPECTRAL/SNBCO_Phi')
541 |
542 | !
543 | ! READING THE PARAMETERS
544 | !
545 | DO I=1,48
546 | READ(1050,*) (KCO(J,I),J=1,12)
547 | READ(1060,*) (DCO(J,I),J=1,12)
548 | END DO
549 | DO I=1,96
550 | READ(1030,*) (KC(J,I),J=1,14)
551 | READ(1040,*) (DC(J,I),J=1,14)
552 | END DO
553 | DO I=1,367
554 | READ(1010,*) (KH(J,I),J=1,14)
555 | READ(1020,*) (DH(J,I),J=1,14)
556 | END DO
557 |
558 | CLOSE(1010)
559 | CLOSE(1020)
560 | CLOSE(1030)
561 | CLOSE(1040)
562 | CLOSE(1050)
563 | CLOSE(1060)
564 |
565 |
566 | RETURN
567 |
568 | END SUBROUTINE PARAM
569 |
570 | !********************************************************************
571 | !********************************************************************
572 | !********************************************************************
573 | !********************************************************************
574 | !********************************************************************
575 | !********************************************************************
576 |
577 |
578 |
579 | SUBROUTINE FINDI(LICO,LICO2,LIH2O,ICO,ICO2,IH2O,C1,C2,WVNB,DWVNB)
580 |
581 |
582 | !*****************************************************************
583 | ! *
584 | ! THIS SUBROUTINE SEARCHS THE PARAMETER INDEXES *
585 | ! CORRESPONDING THE WAVE NUMBER 'WVNB' *
586 | ! *
587 | !*****************************************************************
588 | !
589 |
590 | USE paramax
591 |
592 | IMPLICIT REAL*8 (A-H,O-Z), INTEGER (I-N)
593 |
594 | LOGICAL LICO,LICO2,LIH2O
595 |
596 |
597 | !
598 | ! CALCULATION OF PLANCK LAW PARAMETERS
599 | !
600 | C1=(0.11909E-7)*WVNB**3*DWVNB
601 | C2=1.4388*WVNB
602 |
603 | !
604 | ! LOOKING FOR THE INDEXES
605 | !
606 | IF(WVNB.LT.450.0.OR.WVNB.GT.1200) GO TO 10
607 | ICO2=INT((WVNB-450.)/25.)+1
608 | GO TO 14
609 | 10 IF(WVNB.LT.1950.0.OR.WVNB.GT.2450.) GO TO 11
610 | ICO2=INT((WVNB-1950.)/25.)+32
611 | GO TO 14
612 | 11 IF(WVNB.LT.3300.0.OR.WVNB.GT.3800.) GO TO 12
613 | ICO2=INT((WVNB-3300.)/25.)+53
614 | GO TO 14
615 | 12 IF(WVNB.LT.4700.0.OR.WVNB.GT.5250.) GO TO 13
616 | ICO2=INT((WVNB-4700.)/25.)+74
617 | GO TO 14
618 | 13 ICO2=-9
619 | 14 CONTINUE
620 | IF(WVNB.GE.150.AND.WVNB.LE.9300.) THEN
621 | IH2O=INT((WVNB-150.)/25.)+1
622 | ELSE
623 | IH2O=-9
624 | ENDIF
625 | IF(WVNB.LT.1750.0.OR.WVNB.GT.2325.) GO TO 19
626 | ICO=INT((WVNB-1750.)/25.)+1
627 | GO TO 21
628 | 19 IF(WVNB.LT.3775.0.OR.WVNB.GT.4350.) GO TO 20
629 | ICO=INT((WVNB-3775.)/25.)+25
630 | GO TO 21
631 | 20 ICO=-9
632 | 21 CONTINUE
633 | LICO=ICO.GT.0
634 | LICO2=ICO2.GT.0
635 | LIH2O=IH2O.GT.0
636 | RETURN
637 | END SUBROUTINE FINDI
638 |
639 |
640 |
641 | !********************************************************************
642 | !********************************************************************
643 | !********************************************************************
644 | !********************************************************************
645 | !********************************************************************
646 | !********************************************************************
647 |
648 |
649 |
650 | SUBROUTINE KBARANDPHI(celldata,LICO,LICO2,LIH2O,ICO,ICO2,IH2O,&
651 | KCO,KC,KH,DCO,DC,DH,MIXDATA)
652 |
653 |
654 | !*******************************************************************!
655 | !*******************************************************************!
656 | ! !
657 | ! SUBROUTINE FOR THE CALCULATION OF Kappa moy AND Phi FROM THE !
658 | ! CELLS' TEMPERATURE, PRESSURE AND SPECIES MASS CONCENTRATIONS !
659 | ! !
660 | ! INPUT : T en K ---> celldata(1) !
661 | ! : P en atm ---> celldata(2) !
662 | ! : X_H2O ---> celldata(3) !
663 | ! : X_CO2 ---> celldata(4) !
664 | ! : X_CO ---> celldata(5) !
665 | ! : X_O2 ---> celldata(6) !
666 | ! : X_N2 ---> celldata(7) !
667 | ! : X_SOOT ---> celldata(8) !
668 | ! !
669 | ! OUTPUT : Kappa_moy_Mix en m-1 ---> MIXDATA(1) !
670 | ! Phi_Mix sans unite ---> MIXDATA(2) !
671 | ! !
672 | ! Kappa_moy_H2O en m-1 ---> SNBDATA(1,...) !
673 | ! Kappa_moy_CO2 en m-1 ---> SNBDATA(...,3,...) !
674 | ! Kappa_moy_CO en m-1 ---> SNBDATA(...,5,...) !
675 | ! Phi_H2O sans unite ---> SNBDATA(...,2,...) !
676 | ! Phi_CO2 sans unite ---> SNBDATA(...,4...) !
677 | ! Phi_CO sans unite ---> SNBDATA(...,6) !
678 | ! !
679 | ! !
680 | !*******************************************************************!
681 | !*******************************************************************!
682 |
683 | USE paramax
684 |
685 | IMPLICIT DOUBLE PRECISION (A-C,E-H,O-Z)
686 |
687 | REAL(8),DIMENSION(14,n_SNBmax)::KCO,DCO,KC,DC,KH,DH
688 | REAL(8), DIMENSION(8) :: celldata
689 | REAL(8), DIMENSION(6) :: SNBDATA
690 | REAL(8), DIMENSION(2) :: MIXDATA
691 | REAL(8), DIMENSION(3) :: BIGB
692 | LOGICAL LICO,LICO2,LIH2O
693 |
694 | CALL TMNO(celldata(1),RT,IT)
695 |
696 | RRT=RT
697 | IIT=IT
698 | T296=296./celldata(1)
699 | T273=273./celldata(1)
700 | T900=900./celldata(1)
701 | XN2=1.-celldata(5)-celldata(4)-celldata(3)
702 |
703 | IF(LICO) THEN
704 | GAM=0.07*celldata(4)+0.06*(celldata(5)+XN2+celldata(3))
705 | GAM=celldata(2)*GAM*SQRT(T273)
706 | SNBDATA(5)=KCO(IT,ICO)+RT*(KCO(IT+1,ICO)-KCO(IT,ICO))
707 | SNBDATA(5)=100.*SNBDATA(5)*celldata(2)*celldata(5)
708 | XDCO=DCO(IT,ICO)+RT*(DCO(IT+1,ICO)-DCO(IT,ICO))
709 | SNBDATA(6)=2.*GAM*XDCO
710 | BIGB(3)=SNBDATA(5)**2/SNBDATA(6)
711 | ENDIF
712 | IF(LICO2) THEN
713 | GAM=0.07*celldata(4)+0.058*XN2+0.15*celldata(3)
714 | IF(celldata(1).LE.900.) THEN
715 | GAM=celldata(2)*GAM*(T296)**0.7
716 | ELSE
717 | GAM=celldata(2)*GAM*0.45913*DSQRT(T900)
718 | ENDIF
719 | SNBDATA(3)=KC(IT,ICO2)+RT*(KC(IT+1,ICO2)-KC(IT,ICO2))
720 | SNBDATA(3)=100.*SNBDATA(3)*celldata(2)*celldata(4)
721 | XDCO2=DC(IT,ICO2)+RT*(DC(IT+1,ICO2)-DC(IT,ICO2))
722 | SNBDATA(4)=2.*GAM*XDCO2
723 | BIGB(2)=SNBDATA(3)**2/SNBDATA(4)
724 | ENDIF
725 | IF(LIH2O) THEN
726 | RATT=DSQRT(T296)
727 | GAM=0.066*(7.0*RATT*celldata(3)+1.2*(celldata(3)+XN2)&
728 | +1.5*celldata(4))*RATT
729 | GAM=celldata(2)*GAM
730 | SNBDATA(1)=KH(IT,IH2O)+RT*(KH(IT+1,IH2O)-KH(IT,IH2O))
731 | SNBDATA(1)=100.*SNBDATA(1)*celldata(2)*celldata(3)
732 | XDH2O=DH(IT,IH2O)+RT*(DH(IT+1,IH2O)-DH(IT,IH2O))
733 | SNBDATA(2)=2.*GAM*XDH2O
734 | BIGB(1)=SNBDATA(1)**2/SNBDATA(2)
735 | ENDIF
736 |
737 | ! Gaz mixture modelling !
738 |
739 | MIXDATA(1)=SNBDATA(1)+SNBDATA(3)+SNBDATA(5)
740 | MIXDATA(2)=MIXDATA(1)**2/(BIGB(1)+BIGB(2)+BIGB(3))
741 |
742 |
743 |
744 | RETURN
745 | END SUBROUTINE KBARANDPHI
746 |
747 |
748 |
749 | !********************************************************************
750 | !********************************************************************
751 | !********************************************************************
752 | !********************************************************************
753 | !********************************************************************
754 | !********************************************************************
755 |
756 |
757 |
758 | SUBROUTINE TMNO(TM,RT,IT)
759 |
760 |
761 | !***************************************************************
762 | !***************************************************************
763 | ! *
764 | ! THIS SUBROUTINE CALCULATES THE INTERPOLATION COEFFICIENTS *
765 | ! ( linked to the subroutine KBARANDPHI ) *
766 | ! *
767 | !***************************************************************
768 | !***************************************************************
769 |
770 |
771 | DOUBLE PRECISION TM,RT
772 |
773 | IF(TM.GT.300.0) THEN
774 | IF(TM.LT.2900.0) THEN
775 | RT=(TM-300.0)/200.0
776 | IT=INT(RT+1.0E-6)
777 | RT=RT-IT
778 | IT=IT+1
779 | ELSE
780 | RT=1.0
781 | IT=11
782 | ENDIF
783 | ELSE
784 | RT=0.0
785 | IT=1
786 | ENDIF
787 | RETURN
788 | END SUBROUTINE TMNO
789 |
790 |
791 |
792 | !********************************************************************
793 | !********************************************************************
794 | !********************************************************************
795 | !********************************************************************
796 | !********************************************************************
797 | !********************************************************************
798 |
799 |
800 |
801 | SUBROUTINE k_distributeur(gazdata,n,ka_g_i,w)
802 |
803 |
804 | !##############################################################
805 | !##############################################################
806 | !##############################################################
807 | !##############################################################
808 | !## ##
809 | !## SUBROUTINE EFFECTUANT LA QUADRATURE EN K-DISTRIBUTION ##
810 | !## ##
811 | !##############################################################
812 | !##############################################################
813 | !##############################################################
814 | !##############################################################
815 |
816 |
817 | IMPLICIT NONE
818 |
819 | !.....................................................................
820 | !.....................................................................
821 | INTEGER, PARAMETER :: n_mx = 100
822 | INTEGER n,i,jphi,jkbar
823 | !.....................................................................
824 |
825 | REAL, DIMENSION (n_mx) :: kdeg
826 | REAL(8), PARAMETER :: x_min=0.0
827 | REAL(8), PARAMETER :: x_max=1.0
828 | REAL(8), DIMENSION (n_mx) :: x,w,ka_g_i
829 | !.....................................................................
830 | REAL(8) :: phi,kbar,long,phi2,kbar2
831 | REAL(8), DIMENSION (2) :: gazdata
832 | REAL(8), DIMENSION (0:n_mx+1) :: g
833 | !.....................................................................
834 |
835 |
836 | kbar=gazdata(1)
837 | phi=gazdata(2)
838 |
839 | !.....................................................................
840 | !c Definition de la quadrature
841 | !.....................................................................
842 |
843 | CALL gauleg(x_min,x_max,x,w,n)
844 |
845 | !.....................................................................
846 | !c Parametrisation du tableau de quadrature
847 | !.....................................................................
848 |
849 | g(0)=0.
850 | DO i=1,n
851 | g(i)=dble(x(i))
852 | ENDDO
853 | g(n+1)=0.
854 | long=1.E+0
855 |
856 | !.....................................................................
857 | !c Calcul de la transmittance par k-distribution
858 | !.....................................................................
859 |
860 | DO i=1,n
861 |
862 | CALL COFG(phi,kbar,g(i),g(n),ka_g_i(i))
863 |
864 | ENDDO
865 |
866 | RETURN
867 | END SUBROUTINE k_distributeur
868 |
869 |
870 | !.....................................................................
871 | !.....................................................................
872 | !.....................................................................
873 |
874 |
875 | SUBROUTINE gauleg(x1,x2,x,w,n)
876 |
877 | INTEGER n
878 | REAL(8) :: x1,x2,x(n),w(n)
879 | DOUBLE PRECISION EPS
880 | PARAMETER (EPS=3.d-14)
881 | INTEGER i,j,m
882 | DOUBLE PRECISION p1,p2,p3,pp,xl,xm,z,z1
883 | m=(n+1)/2
884 | xm=0.5d0*(x2+x1)
885 | xl=0.5d0*(x2-x1)
886 | DO i=1,m
887 | z=cos(3.141592654d0*(i-.25d0)/(n+.5d0))
888 | 1 continue
889 | p1=1.d0
890 | p2=0.d0
891 | DO j=1,n
892 | p3=p2
893 | p2=p1
894 | p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j
895 | ENDDO
896 | pp=n*(z*p1-p2)/(z*z-1.d0)
897 | z1=z
898 |
899 | z=z1-p1/pp
900 | if(abs(z-z1).gt.EPS)goto 1
901 | x(i)=xm-xl*z
902 | x(n+1-i)=xm+xl*z
903 | w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
904 | w(n+1-i)=w(i)
905 | ENDDO
906 |
907 | ! x(1)=0.00000
908 | ! w(1)=0.04500
909 | ! x(2)=0.15541
910 | ! w(2)=0.24500
911 | ! x(3)=0.45000
912 | ! w(3)=0.32000
913 | ! x(4)=0.74459
914 | ! w(4)=0.24500
915 | ! x(5)=0.90000
916 | ! w(5)=0.05611
917 | ! x(6)=0.93551
918 | ! w(6)=0.05125
919 | ! x(7)=0.98449
920 | ! w(7)=0.03764
921 |
922 |
923 | return
924 |
925 |
926 | END SUBROUTINE gauleg
927 |
928 |
929 | !.....................................................................
930 | !.....................................................................
931 | !.....................................................................
932 |
933 | !.......................................................................
934 | !c
935 | !c inverse of the cumulative for the inverse transmission function
936 | !c by dichotomy
937 | !c
938 | !c in : * $\phi$ ---> phcofg
939 | !c * $\overline k$ ---> cbcofg
940 | !c * $g$ ---> gcofg
941 | !c
942 | !c out : * $k$ ---> dk
943 | !c
944 | !.......................................................................
945 | SUBROUTINE COFG(phcofg,cbcofg,gcofg,gmax,dk)
946 | !.......................................................................
947 | implicit double precision (a-h,o-z)
948 | !.......................................................................
949 | parameter (icofg1=100)
950 | !.......................................................................
951 | !c Recherche des bornes pour la dichitomie (dkinf,dksup)
952 | !.....................................................................
953 | !c PRINT *,'cbcofg dans cofg en entree',cbcofg
954 |
955 | dksup=cbcofg
956 | dkinf=0.D+0
957 |
958 | 10 CALL cdss(phcofg,cbcofg,0.D+0,dksup,gdek)
959 |
960 | IF (gdek.lt.gcofg) THEN
961 | dkinf=dksup
962 | dksup=dksup*10
963 |
964 | GO TO 10
965 | ENDIF
966 | !
967 | !.....................................................................
968 | !c Dichotomie pour obtenir dk par evaluation de gdek
969 | !.....................................................................
970 |
971 | ! iteration=0
972 | 20 dk=(dksup+dkinf)/2
973 |
974 | CALL cdss(phcofg,cbcofg,0.D+0,dk,gdek)
975 | ! iteration=iteration+1
976 | IF (gdek.lt.gcofg) THEN
977 | dkinf=dk
978 | ELSE
979 | dksup=dk
980 | ENDIF
981 |
982 | !.....................................................................
983 | err=dabs(gdek-gcofg)
984 | eps=gmax/1000000
985 | IF (err.gt.eps) THEN
986 | ! PRINT *,cbcofg,dk,eps,err
987 | GO TO 20
988 | ENDIF
989 |
990 | RETURN
991 | !.......................................................................
992 | END SUBROUTINE COFG
993 |
994 | !.....................................................................
995 | !.....................................................................
996 | !.....................................................................
997 |
998 | !.......................................................................
999 | !c
1000 | !c k-cumulative for surface-surface exchanges
1001 | !c
1002 | !c in : * $\phi$ ---> phcdss
1003 | !c * $\overline k$ ---> cbcdss
1004 | !c * $l$ ---> alcdss
1005 | !c * $k$ ---> ccdss
1006 | !c
1007 | !c out : * $g^{ss} (k;l)$ ---> cdcdss
1008 | !c
1009 | !.......................................................................
1010 | SUBROUTINE cdss(phcdss,cbcdss,alcdss,ccdss,cdcdss)
1011 | !.......................................................................
1012 | implicit double precision (a-h,o-z)
1013 | !.......................................................................
1014 |
1015 | cdss1=dsqrt(1.d0+2.d0/phcdss*cbcdss*alcdss)
1016 |
1017 | cdss2=phcdss*cdss1
1018 |
1019 | cdss3=cbcdss/cdss1
1020 |
1021 | call gicd(cdss2,cdss3,ccdss,cdcdss)
1022 |
1023 | !.......................................................................
1024 | return
1025 | END SUBROUTINE cdss
1026 |
1027 | !.....................................................................
1028 | !.....................................................................
1029 | !.....................................................................
1030 |
1031 | !.......................................................................
1032 | !c
1033 | !c inverse gaussian cumulative
1034 | !c
1035 | !c in : * $\phi$ ---> phgicd
1036 | !c * $\mu$ ---> umgicd
1037 | !c * $x$ ---> xgicd
1038 | !c
1039 | !c out : * $\int_0^x \Im (t ; \mu , \mu \phi ) dt$
1040 | !c ---> cdgicd
1041 | !c
1042 | !.......................................................................
1043 | subroutine gicd(phgicd,umgicd,xgicd,cdgicd)
1044 | !.......................................................................
1045 | implicit double precision (a-h,o-z)
1046 | !.......................................................................
1047 |
1048 | gicd1=xgicd/umgicd
1049 | gicd2=dsqrt(phgicd/gicd1)
1050 | gicd3=-gicd2*(1.d0-gicd1)
1051 | gicd4=-gicd2*(1.d0+gicd1)
1052 |
1053 | call stcd(gicd3,gicd5)
1054 |
1055 | call stcd(gicd4,gicd6)
1056 |
1057 | cdgicd=gicd5+exp(2.d0*phgicd)*gicd6
1058 |
1059 | !.......................................................................
1060 | return
1061 | end
1062 |
1063 | !.....................................................................
1064 | !.....................................................................
1065 | !.....................................................................
1066 |
1067 | !.......................................................................
1068 | !c
1069 | !c cumulative of the standard normal
1070 | !c
1071 | !c in : * $x$ ---> xstcd
1072 | !c
1073 | !c out : * $\Gamma (x)$
1074 | !c ---> cdstcd
1075 | !c
1076 | !.......................................................................
1077 | subroutine stcd(xstcd,cdstcd)
1078 | !.......................................................................
1079 | implicit double precision (a-h,o-z)
1080 | !.......................................................................
1081 | real stcd2,erfc
1082 | !.......................................................................
1083 | stcd1=-xstcd/dsqrt(2.d0)
1084 | stcd2=real(stcd1)
1085 | stcd2=erfc(stcd2)
1086 | cdstcd=dble(stcd2)/2.d0
1087 | !.......................................................................
1088 | return
1089 | end
1090 |
1091 | !.....................................................................
1092 | !.....................................................................
1093 | !.....................................................................
1094 |
1095 | SUBROUTINE GSER(GAMSER,A,X,GLN)
1096 | PARAMETER (ITMAX=100,EPS=3.E-7)
1097 | GLN=GAMMLNA>(A)
1098 | IF(X.LE.0.)THEN
1099 | IF(X.LT.0.)STOP
1100 | GAMSER=0.
1101 | RETURN
1102 | ENDIF
1103 | APA>=A
1104 | SUMAA>=1./A
1105 | DEL=SUMA
1106 | DO 11 N=1,ITMAX
1107 | AP=AP+1.
1108 | DEL=DEL*X/AP
1109 | SUMA=SUMA+DEL
1110 | IF(ABS(DEL).LT.ABS(SUMA)*EPS)GO TO 1
1111 | 11 CONTINUE
1112 | STOP 'A too large, ITMAX too small'
1113 | 1 GAMSER=SUMA*EXP(-XA>+A*LOG(X)-GLN)
1114 | RETURN
1115 | END
1116 |
1117 |
1118 | !.....................................................................
1119 | !.....................................................................
1120 | !.....................................................................
1121 |
1122 | SUBROUTINE GCF(GAMMCF,A,X,GLN)
1123 | PARAMETER (ITMAX=100,EPS=3.E-7)
1124 | GLN=GAMMLNA>(A)
1125 | GOLD=0.
1126 | A0=1.
1127 | A1=X
1128 | B0=0.
1129 | B1=1.
1130 | FAC=1.
1131 | DO 11 N=1,ITMAX
1132 | AN=FLOAT(N)
1133 | ANA=ANA>-A
1134 | A0=(A1+A0*ANA)*FAC
1135 | B0=(B1+B0*ANA)*FAC
1136 | ANF=AN*FAC
1137 | A1=X*A0+ANF*A1
1138 | B1=X*B0+ANF*B1
1139 | IF(A1.NE.0.)THEN
1140 | FAC=1./A1
1141 | G=B1*FAC
1142 | IF(ABS((G-GOLD)/G).LT.EPS)GO TO 1
1143 | GOLD=G
1144 | ENDIF
1145 | 11 CONTINUE
1146 | STOP 'A too large, ITMAX too small'
1147 | 1 GAMMCF=EXP(-XA>+A*ALOG(X)-GLN)*G
1148 | RETURN
1149 | END
1150 |
1151 |
1152 |
1153 | !***********************************************************************
1154 | SUBROUTINE READ_DATA(nomfich,meshstatus,datastatus,meshtype,quadtype,&
1155 | mediumtype,spascheme,neighs,nfcelt,KCO,DCO,KC,DC,KH,DH,&
1156 | xc, yc, zc, V, k_scat, kabs_gray,S, epsil, xfc, yfc, zfc, Tf,&
1157 | norm,celldata,mu, eta, ksi, w,critconv,alpha,ncells,I_SCHEME,&
1158 | ndir,nbounds,end_prog)
1159 |
1160 | !######################################################
1161 | ! ACQUISITION DES 1eres DONNEES
1162 | !#######################################################
1163 | ! Lecture du fichier d'entree INPUT :
1164 | ! Meshstatus :
1165 | ! Quadtype :
1166 | !######################################################
1167 |
1168 | USE paramax
1169 |
1170 |
1171 | IMPLICIT NONE
1172 |
1173 |
1174 | CHARACTER (len=60) :: nomfich
1175 | CHARACTER (len=15) :: meshstatus,datastatus,meshtype,quadtype,&
1176 | mediumtype,spascheme
1177 |
1178 | INTEGER, DIMENSION (nmax,2*nfacemax) :: neighs
1179 |
1180 | INTEGER, DIMENSION (nmax) :: nfcelt
1181 |
1182 |
1183 | REAL(8),DIMENSION(14,n_SNBmax) :: KCO,DCO,KC,DC,KH,DH
1184 |
1185 | REAL(8),DIMENSION(nmax) :: xc, yc, zc, V, k_scat, kabs_gray
1186 |
1187 | REAL(8), DIMENSION (nmax,nfacemax) :: S, epsil, xfc, yfc, zfc, Tf
1188 |
1189 | REAL(8), DIMENSION (nmax,nfacemax,3) :: norm
1190 |
1191 | REAL(8),DIMENSION(nmax,8) :: celldata
1192 |
1193 | REAL(8),DIMENSION(nquadmax) :: mu, eta, ksi, w
1194 |
1195 |
1196 | REAL*8 :: critconv
1197 | REAL*8 :: alpha
1198 |
1199 | INTEGER :: ncells
1200 |
1201 | INTEGER :: I_SCHEME
1202 | INTEGER :: ndir
1203 | INTEGER :: nbounds
1204 |
1205 |
1206 |
1207 | INTEGER :: i, j, k, i2, icell, nquad, nquad1, nquad2, nquad3
1208 |
1209 | INTEGER :: end_prog
1210 |
1211 |
1212 | end_prog=0
1213 |
1214 |
1215 | OPEN(UNIT=1,FILE='INPUT')
1216 |
1217 | READ(1,*)
1218 | READ(1,*) meshtype
1219 |
1220 | READ(1,*)
1221 | READ(1,*) nomfich
1222 |
1223 | READ(1,*)
1224 | READ(1,*) meshstatus
1225 |
1226 | READ(1,*)
1227 | READ(1,*) datastatus
1228 |
1229 | READ(1,*)
1230 | READ(1,*) quadtype
1231 |
1232 | READ(1,*)
1233 | IF (quadtype==SNDOM) THEN
1234 | READ(1,*) nquad
1235 | ELSEIF (quadtype==FVM) THEN
1236 | READ(1,*) nquad1,nquad2
1237 | ELSEIF (quadtype==TNDOM) THEN
1238 | READ(1,*) nquad3
1239 | ENDIF
1240 |
1241 | READ(1,*)
1242 | READ(1,*) spascheme
1243 |
1244 | READ(1,*)
1245 | READ(1,*) mediumtype
1246 |
1247 | READ(1,*)
1248 | READ(1,*) critconv
1249 |
1250 | READ(1,*)
1251 | READ(1,*) ncells
1252 |
1253 | CLOSE(1)
1254 |
1255 | OPEN(UNIT=1,FILE='OUTPUT',FORM='formatted')
1256 |
1257 | !####################################################################
1258 | !####################################################################
1259 | ! !
1260 | ! CREATION DES DONNEES DE CALCUL !
1261 | ! !
1262 | !####################################################################
1263 | !####################################################################
1264 |
1265 | !####################################################################
1266 | ! !
1267 | ! CHOIX DU SCHEMA DE DERIVATION SPATIALE !
1268 | ! !
1269 | !####################################################################
1270 |
1271 |
1272 | IF (spascheme==DMFS) THEN
1273 | alpha=0.5
1274 | I_SCHEME=1
1275 | ELSE IF (spascheme==SMFS) THEN
1276 | alpha=1.0
1277 | I_SCHEME=1
1278 | ELSE IF (spascheme==EXPON) THEN
1279 | alpha=0.0
1280 | I_SCHEME=2
1281 | ELSE
1282 | ! GO TO 9999
1283 | end_prog=1
1284 | ENDIF
1285 |
1286 |
1287 | !####################################################################
1288 | ! !
1289 | ! CREATION DES DONNEES DE QUADRATURE ANGULAIRE !
1290 | ! (Discretisation angulaire selon la quadrature choisie) !
1291 | ! !
1292 | !####################################################################
1293 |
1294 |
1295 |
1296 | OPEN(UNIT=5,FILE='INFILES/Quadrature.in',FORM='unformatted')
1297 |
1298 | IF (quadtype==SNDOM) THEN
1299 | ndir=nquad*(nquad+2)
1300 | IF (nquad==9) THEN !
1301 | ndir=ndir-3 ! Condition de lecture special pour la LC11 !
1302 | ENDIF !
1303 | ELSEIF (quadtype==FVM) THEN
1304 | ndir=8*nquad1*nquad2
1305 | ELSEIF (quadtype==TNDOM) THEN
1306 | ndir=8*nquad3**2
1307 | ENDIF
1308 | DO i=1,ndir
1309 | READ(5) mu(i),eta(i),ksi(i),w(i)
1310 | ENDDO
1311 |
1312 | CLOSE(5)
1313 |
1314 | !...................................................................!
1315 | ! Fin d'acquisition de quadrature !
1316 | !...................................................................!
1317 |
1318 |
1319 | !####################################################################
1320 | !####################################################################
1321 | ! !
1322 | ! CREATION DES DONNEES DE CALCUL !
1323 | ! !
1324 | !####################################################################
1325 | !####################################################################
1326 |
1327 | !####################################################################
1328 | ! !
1329 | ! CHOIX DU SCHEMA DE DERIVATION SPATIALE !
1330 | ! !
1331 | !####################################################################
1332 |
1333 |
1334 | IF (spascheme==DMFS) THEN
1335 | alpha=0.5
1336 | I_SCHEME=1
1337 | ELSE IF (spascheme==SMFS) THEN
1338 | alpha=1.0
1339 | I_SCHEME=1
1340 | ELSE IF (spascheme==EXPON) THEN
1341 | alpha=0.0
1342 | I_SCHEME=2
1343 | ELSE
1344 | ! GO TO 9999
1345 | end_prog=1
1346 | ENDIF
1347 |
1348 |
1349 | !####################################################################
1350 | ! !
1351 | ! CREATION DES DONNEES DE QUADRATURE ANGULAIRE !
1352 | ! (Discretisation angulaire selon la quadrature choisie) !
1353 | ! !
1354 | !####################################################################
1355 |
1356 |
1357 |
1358 | OPEN(UNIT=5,FILE='INFILES/Quadrature.in',FORM='unformatted')
1359 |
1360 | !####################################################################
1361 | ! !
1362 | ! Debut d'acquisition des donnees !
1363 | ! !
1364 | !####################################################################
1365 |
1366 | OPEN(UNIT=11,FILE='INFILES/Cell2cells.in',FORM='unformatted')
1367 | OPEN(UNIT=12,FILE='INFILES/Centerfaces.in',FORM='unformatted')
1368 | OPEN(UNIT=20,FILE='INFILES/Centercells.in',FORM='unformatted')
1369 | OPEN(UNIT=21,FILE='INFILES/Externaldata.in',FORM='formatted')
1370 | OPEN(UNIT=22,FILE='INFILES/Emissivities.in',FORM='unformatted')
1371 | OPEN(UNIT=24,FILE='INFILES/K_Scattering.in',FORM='unformatted')
1372 | OPEN(UNIT=25,FILE='INFILES/Properties.in',FORM='unformatted')
1373 | OPEN(UNIT=27,FILE='INFILES/CLProperties.in',FORM='unformatted')
1374 | OPEN(UNIT=28,FILE='INFILES/CLFaces.in',FORM='unformatted')
1375 | OPEN(UNIT=30,FILE='INFILES/Normals.in',FORM='unformatted')
1376 | OPEN(UNIT=40,FILE='INFILES/Volumesareas.in',FORM='unformatted')
1377 |
1378 | DO i=1,ncells
1379 |
1380 | READ(11) icell,nfcelt(icell),(neighs(icell,j),&
1381 | j=1,(2*nfcelt(icell)))
1382 |
1383 | READ(20) icell,xc(icell),yc(icell),zc(icell)
1384 |
1385 | READ(12) icell,nfcelt(icell),(xfc(icell,i2),yfc(icell,i2),&
1386 | zfc(icell,i2),i2=1,nfcelt(icell))
1387 |
1388 | READ(25) icell,(celldata(icell,j),j=1,8)
1389 |
1390 | READ(27) icell,nfcelt(icell),(Tf(icell,j),j=1,nfcelt(icell))
1391 |
1392 | READ(30) icell,nfcelt(icell),((norm(icell,j,k),k=1,3),&
1393 | j=1,nfcelt(icell))
1394 |
1395 | READ(40) icell,nfcelt(icell),V(icell),(S(icell,j),&
1396 | j=1,nfcelt(icell))
1397 |
1398 | READ(22) icell,nfcelt(icell),(epsil(icell,j), j=1,nfcelt(icell))
1399 |
1400 | READ(24) icell,k_scat(icell)
1401 |
1402 | IF (datastatus=='EXTERNAL') THEN
1403 | READ(21,*) icell,celldata(icell,1),kabs_gray(icell)
1404 | ENDIF
1405 |
1406 | ENDDO
1407 |
1408 | READ(28) nbounds
1409 |
1410 | CLOSE(10)
1411 | CLOSE(11)
1412 | CLOSE(12)
1413 | CLOSE(20)
1414 | CLOSE(21)
1415 | CLOSE(23)
1416 | CLOSE(24)
1417 | CLOSE(25)
1418 | CLOSE(27)
1419 | CLOSE(28)
1420 | CLOSE(30)
1421 | CLOSE(40)
1422 |
1423 |
1424 | CALL PARAM(KCO,KC,KH,DCO,DC,DH)
1425 |
1426 |
1427 | !###################################################################!
1428 | !###################################################################!
1429 | ! !
1430 | ! FIN DE CHARGEMENT DES DONNEES !
1431 | ! !
1432 | !###################################################################!
1433 | !###################################################################!
1434 |
1435 | RETURN
1436 | END
1437 |
1438 |
1439 | !***********************************************************************
1440 | SUBROUTINE GRAY_CASE(mediumtype,nfcelt,wkabs,celldata,all_k_abs,Lb, &
1441 | kabs_gray,Lo, epsil, Tf, ncells,n_k_abs,DWVNB_SI)
1442 |
1443 |
1444 | USE paramax
1445 |
1446 |
1447 | IMPLICIT NONE
1448 |
1449 | CHARACTER (len=15) :: mediumtype
1450 |
1451 | INTEGER, DIMENSION (nmax) :: nfcelt
1452 |
1453 | REAL(8),DIMENSION (nkabs) :: wkabs
1454 |
1455 | REAL(8),DIMENSION(nmax,8) :: celldata
1456 |
1457 | REAL(8),DIMENSION(nmax,nkabs) :: all_k_abs
1458 |
1459 | REAL(8),DIMENSION(nmax) :: Lb, kabs_gray
1460 |
1461 | REAL(8), DIMENSION (nmax,nfacemax) :: Lo, epsil, Tf
1462 |
1463 |
1464 |
1465 | INTEGER :: ielt,m,ncells,n_k_abs,DWVNB_SI
1466 |
1467 | REAL*8 :: planck
1468 |
1469 |
1470 | !###################################################################!
1471 | !###################################################################!
1472 | ! !
1473 | ! SPECTRAL TREATMENT !
1474 | ! !
1475 | !###################################################################!
1476 | !###################################################################!
1477 | ! !
1478 | ! GRAY CASE !
1479 | ! !
1480 | !###################################################################!
1481 | !###################################################################!
1482 |
1483 |
1484 | IF (homosyst=='NO') THEN
1485 |
1486 | DO ielt=1,ncells
1487 |
1488 | DO m=1,nfcelt(ielt)
1489 |
1490 | ! Luminances noires aux parois
1491 | Lo(ielt,m)=epsil(ielt,m)*planck(Tf(ielt,m))
1492 |
1493 | ENDDO
1494 |
1495 | ! Luminances noires dans le milieu
1496 | Lb(ielt)=planck(celldata(ielt,1))
1497 |
1498 | all_k_abs(ielt,:)=kabs_gray(ielt)
1499 |
1500 | ! PRINT*,ielt,celldata(ielt,1),kabs_gray(ielt),(Tf(ielt,m),m=1,nfcelt(ielt))
1501 |
1502 | ENDDO
1503 |
1504 | ELSE IF (homosyst=='YES') THEN
1505 |
1506 | DO ielt=1,ncells
1507 |
1508 | DO m=1,nfcelt(ielt)
1509 | ! Luminances noires aux parois
1510 | Lo(ielt,m)=epsil(ielt,m)*planck(Tf(ielt,m))
1511 |
1512 | ENDDO
1513 | ! Luminances noires dans le milieu
1514 | Lb(ielt)=planck(celldata(ielt,1))
1515 |
1516 | all_k_abs(ielt,:)=kabs
1517 |
1518 | ENDDO
1519 |
1520 | END IF
1521 |
1522 | n_k_abs=1
1523 | wkabs=1.0
1524 | DWVNB_SI=1.0
1525 |
1526 |
1527 |
1528 | RETURN
1529 | END
1530 |
1531 | !***********************************************************************
1532 |
1533 |
1534 | !***********************************************************************
1535 | SUBROUTINE EXTERN(nfcelt,wkabs,celldata,all_k_abs,kabs_gray,Lo, epsil,&
1536 | ncells,WVNB_SI,n_k_abs,Tf,Lb)
1537 |
1538 |
1539 | USE paramax
1540 |
1541 |
1542 | IMPLICIT NONE
1543 |
1544 | INTEGER, DIMENSION (nmax) :: nfcelt
1545 |
1546 | REAL(8),DIMENSION (nkabs) :: wkabs
1547 |
1548 | REAL(8),DIMENSION(nmax,8) :: celldata
1549 |
1550 | REAL(8),DIMENSION(nmax,nkabs) :: all_k_abs
1551 |
1552 | REAL(8),DIMENSION(nmax) :: Lb,kabs_gray
1553 |
1554 | REAL(8), DIMENSION (nmax,nfacemax) :: Lo, epsil, Tf
1555 |
1556 | INTEGER :: ncells
1557 |
1558 | REAL(8) :: blae
1559 | REAL(8) :: WVNB_SI
1560 | INTEGER :: n_k_abs
1561 |
1562 | INTEGER :: ielt,m
1563 |
1564 |
1565 | DO ielt=1,ncells
1566 |
1567 | DO m=1,nfcelt(ielt)
1568 |
1569 | ! Luminances noires aux parois
1570 | Lo(ielt,m)=epsil(ielt,m)/pi*blae(WVNB_SI,Tf(ielt,m))
1571 |
1572 | ENDDO
1573 |
1574 | ! Luminances noires dans le milieu
1575 | Lb(ielt)=blae(WVNB_SI,celldata(ielt,1))/pi
1576 |
1577 | all_k_abs(ielt,:)=1.E8*kabs_gray(ielt)*WVNB_SI
1578 |
1579 | ENDDO
1580 |
1581 | n_k_abs=1
1582 | wkabs=1.0
1583 |
1584 | RETURN
1585 | END
1586 |
1587 | !***********************************************************************
1588 |
1589 |
1590 | !***********************************************************************
1591 | SUBROUTINE NO_HOMOSYST(nfcelt,SNBDATA,wkabs, k_absmel,KCO,DCO,KC,DC,KH,DH,&
1592 | celldata,Lb,all_a1,all_a2,Lo,epsil,Tf,ncells,WVNB_SI,&
1593 | WVNB,LICO,LICO2,LIH2O,ICO,ICO2,IH2O,all_k_abs,WVSOOT)
1594 |
1595 |
1596 | USE paramax
1597 |
1598 |
1599 | IMPLICIT NONE
1600 |
1601 | INTEGER, DIMENSION (nmax) :: nfcelt
1602 |
1603 | REAL(8),DIMENSION(2) :: SNBDATA
1604 |
1605 | REAL(8),DIMENSION (nkabs) :: wkabs, k_absmel
1606 |
1607 | REAL(8),DIMENSION(14,n_SNBmax) :: KCO,DCO,KC,DC,KH,DH
1608 |
1609 | REAL(8),DIMENSION(nmax,8) :: celldata
1610 |
1611 | REAL(8),DIMENSION(nmax) :: Lb,all_a1,all_a2
1612 |
1613 | REAL(8), DIMENSION (nmax,nfacemax) :: Lo, epsil, Tf
1614 |
1615 | REAL(8),DIMENSION (nkabs) :: WVSOOT
1616 |
1617 | REAL(8),DIMENSION(nmax,nkabs) :: all_k_abs
1618 |
1619 | INTEGER :: ncells
1620 |
1621 | REAL(8) :: blae
1622 | REAL(8) :: WVNB_SI,WVNB
1623 |
1624 | LOGICAL :: LICO,LICO2,LIH2O
1625 | INTEGER :: ICO,ICO2,IH2O
1626 |
1627 |
1628 | INTEGER :: ielt,m
1629 |
1630 | DO ielt=1,ncells
1631 |
1632 | DO m=1,nfcelt(ielt)
1633 |
1634 | ! Luminances noires aux parois
1635 | Lo(ielt,m)=epsil(ielt,m)/pi*blae(WVNB_SI,Tf(ielt,m))
1636 |
1637 | ENDDO
1638 |
1639 | ! Luminances noires dans le milieu
1640 | Lb(ielt)=blae(WVNB_SI,celldata(ielt,1))/pi
1641 |
1642 | IF (WVNB<=9300.0) THEN
1643 |
1644 | CALL KBARANDPHI(celldata(ielt,:),LICO,LICO2,LIH2O,ICO,ICO2,IH2O,&
1645 | KCO,KC,KH,DCO,DC,DH,SNBDATA)
1646 |
1647 | CALL k_distributeur(SNBDATA,nkabs,k_absmel,wkabs)
1648 |
1649 | ENDIF
1650 |
1651 | all_k_abs(ielt,:)=k_absmel+WVSOOT
1652 |
1653 | all_a1(ielt)=sum(all_k_abs(ielt,:)*wkabs)
1654 | all_a2=4.*pi*all_a1
1655 |
1656 | ENDDO
1657 |
1658 |
1659 | RETURN
1660 | END
1661 |
1662 | !***********************************************************************
1663 |
1664 | !***********************************************************************
1665 | SUBROUTINE YES_HOMOSYST(nfcelt,SNBDATA,wkabs,k_absmel,KCO,DCO,KC,DC,KH,DH,&
1666 | celldata,Lb,Lo, epsil,Tf,WVSOOT,all_k_abs,ncells,&
1667 | WVNB_SI,LICO,LICO2,LIH2O,ICO,ICO2,IH2O)
1668 |
1669 |
1670 | USE paramax
1671 |
1672 |
1673 | IMPLICIT NONE
1674 |
1675 | INTEGER, DIMENSION (nmax) :: nfcelt
1676 |
1677 | REAL(8),DIMENSION(2) :: SNBDATA
1678 |
1679 | REAL(8),DIMENSION (nkabs) :: wkabs, k_absmel
1680 |
1681 | REAL(8),DIMENSION(14,n_SNBmax) :: KCO,DCO,KC,DC,KH,DH
1682 |
1683 | REAL(8),DIMENSION(nmax,8) :: celldata
1684 |
1685 | REAL(8),DIMENSION(nmax) :: Lb
1686 |
1687 | REAL(8), DIMENSION (nmax,nfacemax) :: Lo, epsil, Tf
1688 |
1689 | REAL(8),DIMENSION (nkabs) :: WVSOOT
1690 |
1691 | REAL(8),DIMENSION(nmax,nkabs) :: all_k_abs
1692 |
1693 | INTEGER :: ncells
1694 |
1695 | REAL(8) :: blae
1696 | REAL(8) :: WVNB_SI
1697 |
1698 | LOGICAL :: LICO,LICO2,LIH2O
1699 | INTEGER :: ICO,ICO2,IH2O
1700 |
1701 |
1702 | INTEGER :: ielt,m
1703 |
1704 | CALL KBARANDPHI(celldata(3206,:),LICO,LICO2,LIH2O,ICO,ICO2,IH2O,&
1705 | KCO,KC,KH,DCO,DC,DH,SNBDATA)
1706 |
1707 | CALL k_distributeur(SNBDATA,nkabs,k_absmel,wkabs)
1708 |
1709 | k_absmel=k_absmel+WVSOOT
1710 |
1711 | DO ielt=1,ncells
1712 |
1713 | DO m=1,nfcelt(ielt)
1714 |
1715 | ! Luminances noires aux parois
1716 | Lo(ielt,m)=epsil(ielt,m)/pi*blae(WVNB_SI,Tf(ielt,m))
1717 |
1718 | ENDDO
1719 |
1720 |
1721 | ! Luminances noires dans le milieu
1722 | Lb(ielt)=blae(WVNB_SI,celldata(ielt,1))/pi
1723 |
1724 | all_k_abs(ielt,:)=k_absmel
1725 |
1726 | ENDDO
1727 |
1728 | RETURN
1729 | END
1730 |
1731 | !***********************************************************************
1732 |
1733 |
1734 | !***********************************************************************
1735 | SUBROUTINE SCHEME_DMFS(neighs,nfcelt,wkabs,pathway,mu,eta,ksi,w,Dm_ij,Lb, Lpi,&
1736 | Gi,G,V,k_scat,k_abs,Q_rtot,Li,Le,H,S,epsil,norm,&
1737 | alpha,ndir,ncells,i_gij)
1738 |
1739 |
1740 | USE paramax
1741 |
1742 |
1743 | IMPLICIT NONE
1744 |
1745 |
1746 | INTEGER, DIMENSION (nmax,2*nfacemax) :: neighs
1747 |
1748 | INTEGER, DIMENSION (nmax) :: nfcelt, pathway
1749 |
1750 | REAL(8),DIMENSION (nkabs) :: wkabs
1751 |
1752 | REAL(8),DIMENSION(nquadmax) :: mu, eta, ksi, w
1753 |
1754 | REAL(8),DIMENSION(nfacemax) :: Dm_ij
1755 |
1756 | REAL(8),DIMENSION(nmax) :: Lb, Lpi, Gi, G, V, k_scat, k_abs
1757 |
1758 | REAL(8),DIMENSION(nmax,3) :: Q_rtot
1759 |
1760 | REAL(8), DIMENSION (nmax,nfacemax) :: Li, Le, H, S, epsil
1761 |
1762 | REAL(8), DIMENSION (nmax,nfacemax,3) :: norm
1763 |
1764 |
1765 | REAL(8) :: alpha
1766 | INTEGER :: ndir,ncells,i_gij
1767 |
1768 | INTEGER :: iquad,icell,i,j,n
1769 |
1770 |
1771 | DO iquad=1,ndir
1772 |
1773 | READ(200) (pathway(n), n=1,ncells)
1774 |
1775 | DO icell=1,ncells
1776 |
1777 | j=pathway(icell)
1778 |
1779 | DO i=1,nfcelt(j)
1780 | Dm_ij(i)=norm(j,i,1)*mu(iquad)+norm(j,i,2)*&
1781 | eta(iquad)+norm(j,i,3)*ksi(iquad)
1782 | ENDDO
1783 |
1784 | CALL MFSCHEME(k_abs(j),wkabs(i_gij),k_scat(j),&
1785 | V(j),nfcelt(j),S(j,:),Dm_ij,Lb(j),Gi(j),Li(j,:),&
1786 | epsil(j,:),Le(j,:),Lpi(j),alpha)
1787 |
1788 | G(j)=G(j)+Lpi(j)*w(iquad)
1789 | Q_rtot(j,1)=Q_rtot(j,1)+Lpi(j)*w(iquad)*mu(iquad)
1790 | Q_rtot(j,2)=Q_rtot(j,2)+Lpi(j)*w(iquad)*eta(iquad)
1791 | Q_rtot(j,3)=Q_rtot(j,3)+Lpi(j)*w(iquad)*ksi(iquad)
1792 |
1793 |
1794 |
1795 |
1796 | !***********************************************!
1797 | ! Mise a jour des faces de sorties de la cellule!
1798 | ! Calcul du flux incident Hw a la parois !
1799 | !***********************************************!
1800 |
1801 | DO i=1,nfcelt(j)
1802 | IF (Dm_ij(i)>=0.) THEN
1803 | IF (epsil(j,i)==-1.) THEN
1804 | Li(neighs(j,(2*i-1)),neighs(j,(2*i)))=&
1805 | Le(j,i)
1806 | ELSE
1807 | H(j,i)=H(j,i)+Le(j,i)*w(iquad)*Dm_ij(i)
1808 | ENDIF
1809 | ENDIF
1810 |
1811 | ENDDO
1812 |
1813 | !***********************************************!
1814 |
1815 | ENDDO
1816 |
1817 | ENDDO
1818 |
1819 |
1820 | RETURN
1821 | END
1822 |
1823 | !***********************************************************************
1824 |
1825 | !***********************************************************************
1826 | SUBROUTINE SCHEME_EXP(neighs,nfcelt, pathway,wkabs,mu, eta, ksi, w,Dm_ij,&
1827 | Lb, Lpi, Gi, G, V, k_scat, k_abs, maxlen,Li, Le, H,&
1828 | S, epsil, s_s,norm, ndir, ncells,i_gij)
1829 |
1830 |
1831 | USE paramax
1832 |
1833 |
1834 | IMPLICIT NONE
1835 |
1836 | INTEGER, DIMENSION (nmax,2*nfacemax) :: neighs
1837 |
1838 | INTEGER, DIMENSION (nmax) :: nfcelt, pathway
1839 |
1840 | REAL(8),DIMENSION (nkabs) :: wkabs
1841 |
1842 | REAL(8),DIMENSION(nquadmax) :: mu, eta, ksi, w
1843 |
1844 | REAL(8),DIMENSION(nfacemax) :: Dm_ij
1845 |
1846 | REAL(8),DIMENSION(nmax) :: Lb, Lpi, Gi, G, &
1847 | V, k_scat, k_abs, maxlen
1848 |
1849 | REAL(8), DIMENSION (nmax,nfacemax) :: Li, Le, H, &
1850 | S, epsil, s_s
1851 |
1852 | REAL(8), DIMENSION (nmax,nfacemax,3) :: norm
1853 |
1854 | INTEGER :: ndir, ncells,i_gij
1855 |
1856 | INTEGER :: iquad, n, k_cell, k_coef, icell, j, i
1857 |
1858 | s_s=0.
1859 |
1860 | OPEN(UNIT=240,FILE='INFILES/L_SPEC.in',&
1861 | FORM='unformatted')
1862 | OPEN(UNIT=260,FILE='INFILES/S_SPEC.in',&
1863 | FORM='unformatted')
1864 |
1865 | DO iquad=1,ndir
1866 |
1867 | READ(200) (pathway(n), n=1,ncells)
1868 | READ(240) (maxlen(k_cell), k_cell=1,ncells)
1869 | READ(260) ((s_s(k_cell,k_coef),k_cell=1,ncells),&
1870 | k_coef=1,4)
1871 |
1872 | DO icell=1,ncells
1873 |
1874 | j=pathway(icell)
1875 |
1876 | DO i=1,nfcelt(j)
1877 | Dm_ij(i)=norm(j,i,1)*mu(iquad)+norm(j,i,2)*&
1878 | eta(iquad)+norm(j,i,3)*ksi(iquad)
1879 | ENDDO
1880 |
1881 | CALL EXPOSCHEME(k_abs(j),wkabs(i_gij),k_scat(j),&
1882 | V(j),S(j,:),Dm_ij,maxlen(icell),s_s(icell,:),&
1883 | Lb(j),Gi(j),Li(j,:),epsil(j,:),Le(j,:),Lpi(j))
1884 |
1885 | G(j)=G(j)+Lpi(j)*w(iquad)
1886 |
1887 | !***********************************************!
1888 | ! Mise a jour des faces de sorties de la cellule!
1889 | ! Calcul du flux incident Hw a la parois !
1890 | !***********************************************!
1891 |
1892 | DO i=1,nfcelt(j)
1893 | IF (Dm_ij(i)>=0) THEN
1894 | IF (epsil(j,i)==-1.) THEN
1895 | Li(neighs(j,(2*i-1)),neighs(j,(2*i)))=&
1896 | Le(j,i)
1897 | ELSE
1898 | H(j,i)=H(j,i)+Le(j,i)*w(iquad)*Dm_ij(i)
1899 | ENDIF
1900 | ENDIF
1901 | ENDDO
1902 |
1903 | ENDDO
1904 |
1905 | ENDDO
1906 |
1907 | CLOSE(240)
1908 | CLOSE(260)
1909 |
1910 |
1911 | RETURN
1912 | END
1913 | !***********************************************************************
1914 |
1915 |
1916 | !***********************************************************************
1917 | SUBROUTINE BAND_INTEG(neighs, nfcelt, pathway, wkabs, all_k_abs, mu, &
1918 | eta, ksi, w, Dm_ij, Lb, Gtot, V, k_scat, maxlen,&
1919 | Srtot, Q_rtot, Lo, Lotot, Htot, S, epsil, s_s, norm,&
1920 | n_k_abs, ncells, I_SCHEME, ndir, DWVNB_SI, alpha)
1921 |
1922 |
1923 | USE paramax
1924 |
1925 | IMPLICIT NONE
1926 |
1927 |
1928 | INTEGER, DIMENSION (nmax,2*nfacemax) :: neighs
1929 |
1930 | INTEGER, DIMENSION (nmax) :: nfcelt, pathway
1931 |
1932 | REAL(8),DIMENSION (nkabs) :: wkabs
1933 |
1934 | REAL(8),DIMENSION(nmax,nkabs) :: all_k_abs
1935 |
1936 | REAL(8),DIMENSION(nquadmax) :: mu, eta, ksi, w
1937 |
1938 | REAL(8),DIMENSION(nfacemax) :: Dm_ij
1939 |
1940 | REAL(8),DIMENSION(nmax) :: Lb, Gtot, V, k_scat, maxlen, Srtot
1941 |
1942 | REAL(8),DIMENSION(nmax,3) :: Q_rtot
1943 |
1944 | REAL(8), DIMENSION (nmax,nfacemax) :: Lo, Lotot, Htot, &
1945 | S, epsil, s_s
1946 |
1947 | REAL(8), DIMENSION (nmax,nfacemax,3) :: norm
1948 |
1949 | INTEGER :: n_k_abs, ncells, I_SCHEME, ndir
1950 |
1951 | REAL(8) :: DWVNB_SI, alpha
1952 |
1953 | !***********************************************************************
1954 | REAL(8),DIMENSION(nmax) :: Lpi, Gi, Gj, G, k_abs, Sr
1955 |
1956 | REAL(8), DIMENSION (nmax,nfacemax) :: Li, Le, L, H
1957 |
1958 | INTEGER :: i_gij, n_iter, i, j
1959 |
1960 | REAL(8) :: aiminan
1961 |
1962 | DO i_gij=1,n_k_abs
1963 |
1964 | G=0.
1965 | Sr=0.
1966 | Lpi=0.
1967 | H=0.
1968 | Le=0.
1969 |
1970 | !###################################################################!
1971 | !###################################################################!
1972 | ! !
1973 | ! POINT DE MISE A JOUR PAR ITERATION (REFLECTION/DIFFUSION) !
1974 | ! !
1975 | !###################################################################!
1976 | !###################################################################!
1977 |
1978 | n_iter=-1 ! Compteur d'iteration initialise
1979 |
1980 | ! IF (n_iter==-1) GOTO 10000
1981 |
1982 |
1983 | n_iter=n_iter+1 ! Compteur d'iteration mis a jour
1984 | ! ( =0 pour le premier tour )
1985 |
1986 |
1987 | !*******************************************************************!
1988 | ! Introduction de la diffusion au niveau du Terme Source S !
1989 | ! Stockage des luminances au faces mailles et de la luminance !
1990 | ! calculees pour chaque maille a l'iteration precedente dans !
1991 | ! les Ij et Gi !
1992 | ! Mise a jour de la luminance aux faces des mailles: Ii=I calculees !
1993 | ! aux parois !
1994 | ! On remet I=Io et Q=Io pour calculer resp. la luminance emise aux !
1995 | ! parois avec reflexion et le flux aux parois !
1996 | ! Mise a zero des matrices G et Qaxes !
1997 | !*******************************************************************!
1998 |
1999 |
2000 | Gi=G
2001 |
2002 | DO i=1,ncells
2003 | DO j=1,nfcelt(i)
2004 | IF (epsil(i,j)/=-1.) THEN
2005 | Li(i,j)=Lo(i,j)+(1-epsil(i,j))/pi*H(i,j)
2006 | ELSE
2007 | Li(i,j)=Le(i,j)
2008 | ENDIF
2009 | ENDDO
2010 | ENDDO
2011 |
2012 | G=0.
2013 | H=0.
2014 | Q_rtot=0.
2015 |
2016 | !*******************************************************************!
2017 | ! Lecture des itineraires de resolutions !
2018 | !*******************************************************************!
2019 |
2020 | k_abs=all_k_abs(:,i_gij)
2021 |
2022 | OPEN(UNIT=200,FILE='INFILES/PROGRESS.in',FORM='unformatted')
2023 |
2024 |
2025 | SELECT CASE (I_SCHEME)
2026 |
2027 |
2028 | !*******************************************************************!
2029 | ! Procedure de MEAN FLUX SCHEME !
2030 | !*******************************************************************!
2031 |
2032 |
2033 | CASE(1)
2034 |
2035 | CALL SCHEME_DMFS(neighs,nfcelt,wkabs,pathway,mu,eta,ksi,w,Dm_ij,Lb, Lpi,&
2036 | Gi,G,V,k_scat,k_abs,Q_rtot,Li,Le,H,S,epsil,norm,&
2037 | alpha,ndir,ncells,i_gij)
2038 | !*******************************************************************!
2039 | ! Procedure de SCHEMA EXPONENTIEL !
2040 | !*******************************************************************!
2041 |
2042 |
2043 | CASE(2)
2044 |
2045 | CALL SCHEME_EXP(neighs,nfcelt, pathway,wkabs,mu, eta, ksi, w,Dm_ij,&
2046 | Lb, Lpi, Gi, G, V, k_scat, k_abs, maxlen,Li, Le, H,&
2047 | S, epsil, s_s,norm, ndir, ncells,i_gij)
2048 |
2049 |
2050 |
2051 | END SELECT
2052 |
2053 | CLOSE(200)
2054 |
2055 |
2056 | !********************************************************************
2057 | ! Module d'iteration avec reflexion *
2058 | !********************************************************************
2059 | ! *
2060 | ! Tant qu'il existera un element parietal dont la valeur sera *
2061 | ! modifiee a plus du pourcentage fixe par critconv la convergence *
2062 | ! ne sera pas satisfaite et la mise a jour des radiosites sera *
2063 | ! conserver pour une nouvelle estimation des champs de luminance *
2064 | ! *
2065 | !********************************************************************
2066 |
2067 |
2068 | DO i=1,ncells
2069 |
2070 | IF ((n_iter/=0).and.(G(i)/=0.)) THEN
2071 | Gj(i)=abs((G(i)-Gi(i))/G(i))
2072 | ELSE
2073 | Gj(i)=1.
2074 | ENDIF
2075 |
2076 | ENDDO
2077 |
2078 | aiminan=maxval(Gj)
2079 |
2080 |
2081 | Gtot=Gtot+G*wkabs(i_gij)*DWVNB_SI !(i_bande)
2082 | Htot=Htot+H*wkabs(i_gij)*DWVNB_SI
2083 | Lotot=Lotot+Lo*wkabs(i_gij)*DWVNB_SI
2084 | Sr=(4*pi*Lb-G)*all_k_abs(:,i_gij)*wkabs(i_gij)*DWVNB_SI
2085 | Srtot=Srtot+Sr
2086 |
2087 | !##############################################
2088 |
2089 | ENDDO
2090 |
2091 |
2092 |
2093 |
2094 |
2095 |
2096 | RETURN
2097 | END
2098 |
2099 | !***********************************************************************
2100 |
2101 |
2102 | !***********************************************************************
2103 |
2104 | !*************************************************************************
2105 | !*************************************************************************
2106 | !*************************************************************************
2107 | !*************************************************************************
2108 | !*************************************************************************
2109 | !*************************************************************************
2110 | !*************************************************************************
2111 | !*************************************************************************
2112 | !*************************************************************************
2113 | !*************************************************************************
2114 | !*************************************************************************
2115 | !*************************************************************************
2116 | !****** ******
2117 | !****** ***** *** * * *** ***** * *** * * **** ******
2118 | !****** * * * ** * * * * * * * ** * * ******
2119 | !****** *** * * * * * * * * * * * * * *** ******
2120 | !****** * * * * ** * * * * * * * ** * ******
2121 | !****** * *** * * *** * * *** * * **** ******
2122 | !****** ******
2123 | !*************************************************************************
2124 | !*************************************************************************
2125 | !*************************************************************************
2126 | !*************************************************************************
2127 | !*************************************************************************
2128 | !*************************************************************************
2129 | !*************************************************************************
2130 | !*************************************************************************
2131 | !*************************************************************************
2132 | !*************************************************************************
2133 | !*************************************************************************
2134 | !*************************************************************************
2135 |
2136 |
2137 | !c.......................................................................
2138 | !c
2139 | !c calcul de l'emittance monochromatique du corps noir a un nombre
2140 | !c d'onde et une temperature donnee (A.Delataillade, 2001)
2141 | !c
2142 | !c en entree : * le nombre d'onde blanu (m-1)
2143 | !c * la temperature blat (K)
2144 | !c
2145 | !!c en sortie : * l'emittance blae
2146 | !c
2147 | !c.......................................................................
2148 | REAL(8) FUNCTION BLAE(blanu,blat)
2149 | !c.......................................................................
2150 | IMPLICIT REAL*8 (a-h,o-z)
2151 | !c.......................................................................
2152 | REAL(8):: sigma,pi,c0,hp,cbol,rind,c,c1,c2
2153 | !c.......................................................................
2154 | sigma=5.6693D-8
2155 | pi=datan(1.D+0)*4.D+0
2156 | c0=2.9979D+8
2157 | hp=6.6262D-34
2158 | cbol=1.3806D-23
2159 | rind=1.D+0
2160 | c=c0/rind
2161 | c1=hp*(c**2)
2162 | c2=hp*c/cbol
2163 | !c.......................................................................
2164 | IF (blat .eq. 0.d+0) THEN
2165 | blae = 0d+0
2166 | ELSE
2167 | blae=2.D+0*pi*c1*blanu**3/(dexp(c2*blanu/blat)-1.D+0)
2168 | !cc blae=blat
2169 | ENDIF
2170 | !c les erreurs due a une temperature proche de zero ne sont pas evitees
2171 | !c le probleme ne s'est jamais pose: amaury le 1999 II 22
2172 | !c.......................................................................
2173 | return
2174 | END FUNCTION BLAE
2175 | !c.......................................................................
2176 |
2177 | !c.......................................................................
2178 | !c
2179 | !c calcul de l'emittance monochromatique du corps noir a un nombre
2180 | !c d'onde et une temperature donnee (C.Caliot, 2004)
2181 | !c
2182 | !c en entree : * le nombre d'onde blanu (m-1)
2183 | !c * la temperature blat (K)
2184 | !c
2185 | !!c en sortie : * l'emittance blae
2186 | !c
2187 | !c.......................................................................
2188 |
2189 | FUNCTION plancknu(wavnb,temp,lum)
2190 | ! wavnb en metre-1
2191 | real*8 pi,h,c,kb,plancknu
2192 | real*8 c1,c2,c3,lum,lamb,temp,wavnb
2193 | parameter(h=6.626176d-34,kb=1.380662d-23,c=2.998d8)
2194 | ! J.s J/K m/s
2195 | c3=2.d0*h*c*c
2196 | lum=(h*c*wavnb)/(kb*temp)
2197 | lum=exp(lum)
2198 | lum= lum-1.d0
2199 | lum= (c3*wavnb*wavnb*wavnb)/lum
2200 | plancknu=lum
2201 | ! W.sr-1.m-2.(m-1)-1
2202 | end FUNCTION plancknu
2203 | !
2204 |
2205 | !.....................................................................
2206 | !.....................................................................
2207 | !.....................................................................
2208 | !c.......................................................................
2209 | !c
2210 | !c calcul de l'emittance du corps noir a une temperature donnee
2211 | !c (D.J., 2004)
2212 | !c
2213 | !c en entree : * la temperature temp (K)
2214 | !c
2215 | !c
2216 | !!c en sortie : * l'emittance planck
2217 | !c
2218 | !c.......................................................................
2219 |
2220 | FUNCTION planck(temp)
2221 | !c.......................................................................
2222 | !c.......................................................................
2223 | USE paramax
2224 |
2225 | IMPLICIT REAL*8 (a-h,o-z)
2226 | !c.......................................................................
2227 | !c.......................................................................
2228 |
2229 | planck=sigma*(temp)**4/pi
2230 |
2231 |
2232 | end FUNCTION planck
2233 | !
2234 |
2235 | !.....................................................................
2236 | !.....................................................................
2237 | !.....................................................................
2238 |
2239 | FUNCTION ERFC(X)
2240 | IF(X.LT.0.)THEN
2241 | ERFC=1.+GAMMP(.5,X**2)
2242 | ELSE
2243 | ERFC=GAMMQ(.5,X**2)
2244 | ENDIF
2245 | RETURN
2246 | END
2247 |
2248 | !.....................................................................
2249 | !.....................................................................
2250 | !.....................................................................
2251 |
2252 |
2253 | FUNCTION GAMMP(A,X)
2254 | IF(XA>.LT.0..OR.A.LE.0.)STOP
2255 | IF(XA>.LT.A+1.)THEN
2256 | CALL GSER(GAMSERA>,A,X,GLN)
2257 | GAMMP=GAMSER
2258 | ELSE
2259 | CALL GCF(GAMMCFA>,A,X,GLN)
2260 | GAMMP=1.-GAMMCF
2261 | ENDIF
2262 | RETURN
2263 | END
2264 |
2265 | !.....................................................................
2266 | !.....................................................................
2267 | !.....................................................................
2268 |
2269 | FUNCTION GAMMQ(A,X)
2270 | IF(XA>.LT.0..OR.A.LE.0.)STOP
2271 | IF(XA>.LT.A+1.)THEN
2272 | CALL GSER(GAMSERA>,A,X,GLN)
2273 | GAMMQ=1.-GAMSER
2274 | ELSE
2275 | CALL GCF(GAMMCFA>,A,X,GLN)
2276 | GAMMQ=GAMMCF
2277 | ENDIF
2278 | RETURN
2279 | END
2280 |
2281 |
2282 | !.....................................................................
2283 | !.....................................................................
2284 | !.....................................................................
2285 |
2286 | FUNCTION gammln(xx)
2287 | REAL gammln,xx
2288 | INTEGER j
2289 | DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
2290 | SAVE cof,stp
2291 | DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,&
2292 | 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,&
2293 | -.5395239384953d-5,2.5066282746310005d0/
2294 | x=xx
2295 | y=x
2296 | tmp=x+5.5d0
2297 | tmp=(x+0.5d0)*log(tmp)-tmp
2298 | ser=1.000000000190015d0
2299 | do 11 j=1,6
2300 | y=y+1.d0
2301 | ser=ser+cof(j)/y
2302 | 11 continue
2303 | gammln=tmp+log(stp*ser/x)
2304 |
2305 | return
2306 | END