1 | Program cecile
2 |
3 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
4 | c****************** PROGRAMME CECILE **********************************************
5 | c PROCHAIN CHANTIERS: RESTART PRECOMPIL ECRITUREZIPPEE
6 | c prog principal
7 | c entree:cecile.in(don:ascii) sortie:cecile.out(son:ascii)
8 | c sortie ecran: ascii formate html
9 | c aucune interface graphique
10 | c aucune intervention utilisateur
11 | c pas de reprise possible calcul en cas d'interruption
12 | c diagostics en cours de programme
13 | c calcul par la methode de Monte Carlo des echanges entre zones cyl. ou planes 1-D
14 | c calcul sur le spectre IR de bandes etroites: decoupage de Soufiani et Taine
15 | c calcul avec gaz+suie + hétérogénéités
16 | c pas de gradient de pression p=1 atm
17 |
18 | c les variables "var" ont été supprimées par "cccc" car memoire vive trop petite
19 | c le 17 nov 1999 (sauf var_psi)
20 | c voir les fichiers: cecile.inc trajet_sens.f mcecile.f fich.f
21 |
22 | c redecouper le code comme dans cecile 2D avec notion de chemin a passer a trajet
23 | c se conformer a la notion d'objet dans l'ecriture des ss prog
24 | c penser aussi le futur decoupage pour la parallelisation
25 | c le decoupage des données entrée <> données sortie: subdiviser (ou) profil ss maill
26 | c repenser le nombre de tirage par maille avec la notion d'énergie ????
27 | c penser à la notion de CL radiative !!!!
28 |
29 | c integrer modif pour spectre des suies patrice
30 | c traitement des ecarts types differents si correlés ou pas !!!!
31 | c l'integrale du chemin optique sous exp ne peut pas glisser dans le lot des int stat!
32 | c sauf si on linearise exp
33 | c un autre modele de spectre (priorite) ! reflexion speculaire ! a voir...
34 | c auto-apprentissage des lois de proba ????????
35 | c Makefile plus intelligent qui detecte le systeme
36 | c sous CVS ? multiversion multi utilisateur
37 |
38 | c Thèse A. de Lataillade 1997-2001
39 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
40 |
41 |
42 |
43 |
44 |
45 |
46 | c***********************************************************************
47 | c declaration des modules
48 | c***********************************************************************
49 | c neant
50 |
51 |
52 |
53 |
54 |
55 |
56 |
57 | c***********************************************************************
58 | c declarations des var du programme
59 | c***********************************************************************
60 |
61 | implicit none
62 |
63 | c par fichier include par bloc
64 | INCLUDE 'cecile.inc'
65 |
66 | include 'propradia.inc'
67 | include 'propradiabis.inc'
68 |
69 | include 'radiatif.inc'
70 | include 'entre.inc'
71 |
72 |
73 | c par declarations detaillees
74 | c fichier_rayons : nom du fichier qui contient les rayons externes
75 | c de chaque couronne
76 | c fichier_ks_suie : nom du fichier qui contient les coefficients
77 | c d'absorption pour la suie dans chaque couronne
78 | c fichier_kgbar_gaz : nom du fichier qui contient les coefficients
79 | c d'absorption moyen pour le gaz dans chaque couronne
80 | c fichier_n_tirage : nom du fichier qui contient les nombres de
81 | c tirages effectues depuis chaque couronne
82 | c fichier_temp:nom du fichier qui contient les temperatures de chaque
83 | c maille
84 | character*100 fichier_n_tirage
85 | CHARACTER*100 fichier_rayons
86 | CHARACTER*100 fichier_rsup
87 | CHARACTER*100 fichier_pression_tot
88 | CHARACTER*100 fichier_temperature
89 | CHARACTER*100 fichier_fm_co
90 | CHARACTER*100 fichier_fm_co2
91 | CHARACTER*100 fichier_fm_h2o
92 | CHARACTER*100 fichier_fv
93 | CHARACTER*100 affichage
94 | CHARACTER*10 text
95 | CHARACTER*2 choixheterogeneite, ver
96 |
97 | c.......................................................................
98 | c boolspec = .True. on fait l'integration sur le spectre IR
99 | c .False. on fait las calculs sur une seule bande etroite
100 | c dans ce cas pour l'instant les 3 gaz prennent les
101 | c valeurs pour kbar et phi
102 | c boolprof = .True. on fait avec profil sous maille
103 |
104 |
105 | DOUBLE PRECISION tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8,tmp9
106 | double precision speckgbar, specphi, tmin, tmax,tmp20, tmp10
107 | double precision alphac, skl, moyco, moyco2, moyh2o, moyfv
108 | double precision proba, tran !observer stat des bandes
109 | double precision bilfin, kbartot !pour hypothese optiquementmince
110 | double precision totalbis
111 | double precision correc,invarg !inhmogeneites
112 | double precision taub12, taub1, taub2,kcg,phicg
113 | double precision kbl, phieq
114 | INTEGER itmp1,itmp2,itmp3,itmp4,itmp5,iseed
115 | integer i,j, total
116 | logical boolspec, boolprof, booltrie, boolener, boolplat
117 | logical opt_bande, vacuum, boolssmaille, boolinterieur
118 | logical boolgraine, flag0,flag1
119 | dimension invarg(ngaz_mx)
120 | dimension skl(nbande_mx),proba(nbande_mx),tran(nbande_mx)
121 | dimension kbl(ngaz_mx,nbande_mx),phieq(ngaz_mx,nbande_mx)
122 | dimension bilfin(0:n_mx+1)
123 | dimension correc(ngaz_mx)
124 | common/comsee/iseed
125 | !namelist: lecture de fichier par blocs identifiés
126 | namelist/meshformule/n,rsup,boolener
127 | namelist/meshprof/boolinterieur,boolssmaille
128 | namelist/spectral/boolspec,speckgbar,specphi,boolprof
129 | &,choixheterogeneite
130 | namelist/stat/booltrie,opt_bande
131 | namelist/convcrit/tmp1
132 | namelist/file/fichier_n_tirage,
133 | &fichier_rayons,fichier_pression_tot,
134 | &fichier_temperature,fichier_fm_co,fichier_fm_co2,fichier_fm_h2o,
135 | &fichier_fv
136 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
137 | c
138 | c.......................................................................
139 | c definition de la configuration,(des proprietes)et des controles du MCM
140 | c.......................................................................
141 | c
142 | c-----------ATTENTION : SI ON TIRE PLAT METTRE boolplat=.True.
143 | c on perd moins a tirer plat qd on est épais
144 | c qu'a tirer lambertien quand on est mince
145 | boolplat =.True.
146 | c----------- ci dessous: la varibale "boolinterieur" sert a specifier si on
147 | c utilise un profil sous maille: False==pas de profil sous maille
148 | c /!\ !!! variable suivante a False AUSSI
149 | c True == oui un profil, maniere de le calculer
150 | c voir la variable suivante
151 | boolinterieur = .True.
152 | c----------- boolssmaille si False == profil ss mille temp est par index
153 | c sinon il est fonctionnel (True)
154 | boolssmaille = .False.
155 | if (boolinterieur .eqv. .False.) boolssmaille = .False.
156 | c----------- CHOIX DE LA METHODE POUR INHOMOGENEITE: CG CK NO ?
157 | c NO = separabilite = CG avec la correction a f(phi)--> f(phiCG): 1 (dissymétrie)
158 | c CG = curtis-godson: (probleme des deux parois dans trajet) 2 (symétrie)
159 | c CK = ck methode: 3 (symétrie)
160 | choixheterogeneite = 'CG'
161 | !open(21,FILE='z.debug',STATUS='unknown')
162 | c----------- reutilisation d'un germe different a chaque fois
163 | c POUR COMMENCER TOUJOURS AVEC LE MEME GERME: COMMENTER LES 3 LIGNES SUIVANTES
164 | boolgraine= .False. !sans effet pour l'instant
165 | open(10,FILE='SNBiseed',STATUS='unknown')
166 | read(10,*) iseed
167 | close(10)
168 | iseed=889886264 !la graine reste fixe
169 | open(10,FILE='iseed.out',STATUS='unknown')
170 | write(10,*) 'graine de depard dans fichier SNBiseed: ', iseed
171 | close(10)
172 | ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
173 | call ftete
174 | call frouge
175 | write(*,*) 'Lecture des donnees'
176 | call fligne
177 | call freset
178 | write(*,*) 'iseed= ', iseed
179 | call fligne
180 | OPEN(9,FILE='cecile.in',STATUS='unknown')
181 | c
182 | READ (9,*) ver
183 | ! detection de la version du fichier cecile.in
184 | if ((ver .eq. 'v2').or.(ver .eq. 'V2')) then
185 | print*, 'cecile.in ver. 2'
186 | rewind(9) ! on revient au debut du fichier pour le scruter
187 | 10 continue
188 | read(9,meshformule,end=11,err=10)
189 | goto 10
190 | 11 continue
191 | rewind(9)
192 | 12 continue
193 | read(9,spectral,end=13,err=12)
194 | goto 12
195 | 13 continue
196 | rewind(9)
197 | 14 continue
198 | read(9,stat,end=15,err=14)
199 | goto 14
200 | 15 continue
201 | rewind(9)
202 | 16 continue
203 | read(9,file,end=17,err=16)
204 | goto 16
205 | 17 continue
206 | rewind(9)
207 | 18 continue
208 | read(9,meshprof,end=19,err=18)
209 | goto 18
210 | 19 continue
211 | rewind(9)
212 | else
213 | print*, 'cecile.in ver. 1'
214 | c.................................................................
215 | !debut de passage pour version 1 de cecile.in
216 | READ (9,*)
217 | READ (9,*)
218 | READ (9,*) n
219 | READ (9,*)
220 | c
221 | READ (9,*)
222 | READ (9,'(a)') fichier_n_tirage
223 | READ (9,*)
224 | c
225 | READ (9,*)
226 | READ (9,'(a)') fichier_rayons
227 | READ (9,*)
228 | c
229 | READ (9,*)
230 | READ (9,*) rsup
231 | READ (9,*)
232 | c
233 | READ (9,*)
234 | READ (9,'(a)') fichier_pression_tot
235 | READ (9,*)
236 | c
237 | READ (9,*)
238 | READ (9,'(a)') fichier_temperature
239 | READ (9,*)
240 | c
241 | READ (9,*)
242 | READ (9,'(a)') fichier_fm_co
243 | READ (9,*)
244 | c
245 | READ (9,*)
246 | READ (9,'(a)') fichier_fm_co2
247 | READ (9,*)
248 | c
249 | READ (9,*)
250 | READ (9,'(a)') fichier_fm_h2o
251 | READ (9,*)
252 | c
253 | READ (9,*)
254 | READ (9,'(a)') fichier_fv
255 | READ (9,*)
256 | c
257 | READ (9,*)
258 | READ (9,*) boolener
259 | READ (9,*)
260 | c
261 | READ (9,*)
262 | READ (9,*) boolspec
263 | READ (9,*)
264 | c
265 | READ (9,*)
266 | READ (9,*) speckgbar
267 | READ (9,*)
268 | c
269 | READ (9,*)
270 | READ (9,*) specphi
271 | READ (9,*)
272 | c
273 | READ (9,*)
274 | READ (9,*) boolprof
275 | READ (9,*)
276 | c
277 | READ (9,*)
278 | READ (9,*) booltrie
279 | READ (9,*)
280 | c
281 | READ (9,*)
282 | READ (9,*) opt_bande
283 | READ (9,*)
284 | c
285 | CLOSE(9)
286 | c fin du passage pour la version 1 de cecile.in
287 | endif
288 | c fin du IF pour gerer les deux versions de cecile.in
289 | c................................................................
290 | IF (n.GT.n_mx)
291 | & PRINT *,'*************** ALERTE n ***************'
292 | c
293 | OPEN(20,FILE=fichier_n_tirage,STATUS='old')
294 | DO in=0,n+1
295 | READ (20,*) ntir(in)
296 | ENDDO
297 | CLOSE(20)
298 | c
299 | OPEN(20,FILE=fichier_rayons,STATUS='old')
300 | READ (20,*) r(0)
301 | DO in=1,n
302 | READ (20,*) r(in)
303 | ! tester si rayons croissant
304 | if (r(in) .lt. r(in-1)) then
305 | print *, 'ERROR: rayons non-croissant'
306 | STOP
307 | endif
308 | ENDDO
309 | read (20,*) r(n+1)
310 | CLOSE(20)
311 | c r(n+1)=r(n)
312 | c
313 | OPEN(20,FILE=fichier_pression_tot,STATUS='old')
314 | read (20,*) ptot
315 | c DO in=1,n
316 | c READ (20,*) ks(in)
317 | c ENDDO
318 | c affectation arbitraire pression
319 | c do in=0,n+1
320 | c ptot(in)=1*1e5
321 | c enddo
322 | CLOSE(20)
323 | c
324 | OPEN(20,FILE=fichier_temperature,STATUS='old')
325 | DO in=0,n+1
326 | READ (20,*) temp(in)
327 | ENDDO
328 | c quand on part d'une paroi il faut tirer le kbar des gaz
329 | c cette valeur est tire aleatoirement mais il faut que la
330 | c variable declaree est une valeur: cf les 2 affectation ci-dessous
331 | c kgbar(0)=kgbar(1)
332 | c kgbar(n+1)=kgbar(n)
333 | CLOSE(20)
334 | c
335 | OPEN(20,FILE=fichier_fm_co,STATUS='old')
336 | DO in=1,n
337 | READ (20,*) fm(1,in)
338 | ENDDO
339 | fm(1,0)= 0.d+0
340 | fm(1,n+1)= 0.d+0
341 | close(20)
342 | c
343 | OPEN(20,FILE=fichier_fm_co2,STATUS='old')
344 | DO in=1,n
345 | READ (20,*) fm(2,in)
346 | ENDDO
347 | fm(2,0)= 0.d+0
348 | fm(2,n+1)= 0.d+0
349 | close(20)
350 | c
351 | OPEN(20,FILE=fichier_fm_h2o,STATUS='old')
352 | DO in=1,n
353 | READ (20,*) fm(3,in)
354 | ENDDO
355 | fm(3,0)= 0.d+0
356 | fm(3,n+1)= 0.d+0
357 | close(20)
358 | c
359 | OPEN(20,FILE=fichier_fv,STATUS='old')
360 | DO in=1,n
361 | READ (20,*) fv(in)
362 | ENDDO
363 | if (fv(n/2) .eq. 0.d+0) then
364 | write(*,*) 'SOOT is not used'
365 | call fligne
366 | endif
367 | fv(0)= 0.d+0
368 | fv(n+1)= 0.d+0
369 | CLOSE(20)
370 | c................................................................
371 |
372 | c boolprof ne concerne que la partie energetique qd on est sur un
373 | c bande étroite: False description en temperature, True en lumina
374 | c il faut donc imposer boolprof a False si boolspec est true
375 | c Amaury novembre 2002
376 | if (boolspec .eqv. .True.) boolprof = .False.
377 | ! print *, 'TEST PASSAGE'
378 | ! print *, n, rsup, fichier_n_tirage, fichier_rayons,
379 | ! & fichier_pression_tot,
380 | ! & fichier_temperature, fichier_fm_h2o, fichier_fm_co2,
381 | ! & fichier_fm_co, fichier_fv,
382 | ! & boolener, boolspec, speckgbar, specphi,boolprof,
383 | ! & booltrie, opt_bande, boolinterieur, boolssmaille
384 | ! print *, 'FIN TEST'
385 | ! read *
386 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
387 | c lecture de la structure en bande commune a tous les spectres
388 | c moleculaires
389 | open(unit=9, file='SNBWN')
390 |
391 | c lecture initial pour boucle conditionnelle
392 | ibande =1
393 | read(9,*) eta(ibande), delta_eta(ibande)
394 |
395 | c lecture boucle
396 | ibande=1
397 | do while (eta(ibande).ge.0)
398 | read(9,*) eta(ibande+1), delta_eta(ibande+1)
399 | ibande=ibande+1
400 | enddo
401 | c si on decouvre le nombre de bandes
402 | c nbande =ibande - 1
403 |
404 | close(9)
405 | c------------------------------------------------------------------
406 | c lecture des donnees de bandes de Taine: kbar et 1/delta
407 | call parambgaz
408 |
409 | c------------------------------------------------------------------
410 | c correspondance entre index de bandes et index des param de band
411 | call paramind
412 |
413 | c------------------------------------------------------------------
414 | c ck par méthode d'interpolation numerique a partir tabulations
415 | call init_ck
416 | c------------------------------------------------------------------
417 | write(*,*) 'fin lecture des donnees'
418 | call fligne
419 |
420 | call frouge
421 | print *, 'mises à zéro ET précalcul'
422 | call fligne
423 | call freset
424 | c.......................................................................
425 | c vrai pour les flags de warning etc...
426 | c.......................................................................
427 |
428 | flag0 = .True.
429 | flag1 = .True.
430 |
431 |
432 | c.......................................................................
433 | c mise a zero de toutes les variables de sortie
434 | c.......................................................................
435 | print *, 'sur l''espace:...'
436 | c phig=0.D+0
437 |
438 | print *, 'kbar et phig...'
439 | DO in=1,n
440 |
441 | cccc bilan(in)=0.D+0
442 | cccc var_bilan(in)=0.D+0
443 |
444 | ksbar(in)=0.D+0
445 | do itmp1=1,3
446 | kgbar(itmp1,in)=0.D+0
447 | phig(itmp1,in)=1.D-2
448 | enddo
449 |
450 | ENDDO
451 |
452 |
453 | print *,'min max moyenne...'
454 | tmin = temp(0)
455 | tmax = temp(0)
456 | moyco = 0.d+0
457 | moyco2=0.d+0
458 | moyh2o=0.d+0
459 | moyfv = 0.d+0
460 | print *, 'psi var_psi...'
461 | DO in=0,n+1
462 | moyco = moyco + fm(1,in)
463 | moyco2=moyco2 + fm(2,in)
464 | moyh2o=moyh2o + fm(3,in)
465 | moyfv=moyfv + fv(in)
466 | c ci-dessous pour traiter le plan parallele qd rsup grand
467 | c cette valeur doit etre retranchee a la fin des calculs
468 | r(in)=r(in) + rsup
469 | c pour la partie optimisation tirage de bande on veut connaitre
470 | c la Tmin et Tmax du systeme
471 | if (temp(in).gt.tmax) then
472 | tmax = temp(in)
473 | endif
474 | if (temp(in).lt.tmin) then
475 | tmin = temp(in)
476 | endif
477 |
478 | DO iin=0,n+1
479 | psi(in,iin)=0.D+0
480 | var_psi(in,iin)=0.D+0
481 | DO itmp1=1,n
482 | c dpsidk(in,iin,itmp1)=0.D+0
483 | c var_dpsidk(in,iin,itmp1)=0.D+0
484 | dpsi6dfv(in,iin,itmp1)=0.D+0
485 | cccc var_dpsi6dfv(in,iin,itmp1)=0.D+0
486 |
487 | cc dpsidkgbar(in,iin,itmp1)=0.D+0
488 | cc var_dpsidkgbar(in,iin,itmp1)=0.D+0
489 | dpsi6dfmco(in,iin,itmp1)=0.D+0
490 | cccc var_dpsi6dfmco(in,iin,itmp1)=0.D+0
491 | dpsi6dfmco2(in,iin,itmp1)=0.D+0
492 | cccc var_dpsi6dfmco2(in,iin,itmp1)=0.D+0
493 | dpsi6dfmh2o(in,iin,itmp1)=0.D+0
494 | cccc var_dpsi6dfmh2o(in,iin,itmp1)=0.D+0
495 | cc
496 | dpsi6dt(in,iin,itmp1)=0.D+0
497 | cccc var_dpsi6dt(in,iin,itmp1)=0.D+0
498 | ENDDO
499 | ENDDO
500 | ENDDO
501 | c............boucle pour calul temp au interface de 1 a n
502 | print *,'valeurs interfaciales des volumes...'
503 | call fligne
504 | tempinterf(0) = temp(0)
505 | ! print *, tempinterf(0)
506 | do i=1,n-1
507 | tmp1 = r(i+1)-r(i)
508 | tmp2 = r(i) - r(i-1)
509 | tempinterf(i)=(tmp1*temp(i)+tmp2*temp(i+1))/(r(i+1)-r(i-1))
510 | ! print *, tempinterf(i)
511 | enddo
512 | tempinterf(n) = temp(n+1)
513 | ! print *, tempinterf(n)
514 | ! read *
515 | c...les rcentres
516 | rc(0)=r(0)
517 | ! print *, rc(0)-rsup
518 | do i=1,n
519 | rc(i)=0.5*(r(i-1)+r(i))
520 | ! print *, rc(i)-rsup
521 | enddo
522 | rc(n+1)=r(n+1)
523 | ! print *, rc(n+1)-rsup
524 | ! read *
525 | c........................................................
526 | moyco = moyco / (n+2)
527 | moyco2=moyco2 / (n+2)
528 | moyh2o=moyh2o / (n+2)
529 | moyfv =moyfv /(n+2)
530 |
531 |
532 |
533 | c.........................................................
534 | print *, 'sur le spectre:...'
535 | print *, 'tran proba psi_par_bande...'
536 | print *, 'kbl phieq...'
537 | call fligne
538 | ngaz = 3
539 | do ibande=1,nbande
540 | tran(ibande) = 0.d+0
541 | proba(ibande)=0.d+0
542 | psi60bande(ibande)= 0.d+0
543 | !.........
544 | !calcul de phieq
545 | do igaz= 1,ngaz
546 | kbl(igaz,ibande) = 0.d+0
547 | phieq(igaz,ibande) = 0.d+0
548 | enddo
549 | do in=1,n
550 | iin=in
551 | call modbgaz(ptot,temp(in),fm(1,in),fm(2,in),fm(3,in),flag1)
552 | do igaz=1,ngaz
553 | kbl(igaz,ibande)=kbl(igaz,ibande)
554 | & + kgbar(igaz,in)*(r(in)-r(in-1)) *1.66
555 | phieq(igaz,ibande)= phieq(igaz,ibande)
556 | & + phig(igaz,in)*kgbar(igaz,in)*(r(in)-r(in-1))*1.66
557 | enddo
558 | enddo
559 | do igaz=1,ngaz
560 | phieq(igaz,ibande)=phieq(igaz,ibande)/kbl(igaz,ibande)
561 | if (kbl(igaz,ibande).lt.1.d-20) phieq(igaz,ibande)=1.d-2
562 | enddo
563 | enddo
564 | total = 0.+0
565 | totalbis = 0.d+0
566 |
567 | c
568 | c-------------------------------------------------------------------
569 | c
570 |
571 | print *
572 | print *,'SAISIE: reaffichage fichier pour verif'
573 | print *
574 | call fligne
575 |
576 | ! call cdss(0.280095026,1.d+0,0.d+0,89.3853589,tmp9)
577 | ! print *, 'tmp9',tmp9
578 |
579 |
580 | c print *,'donnees numero et nu bande'
581 | c do ibande=1,nbande
582 | c print *, ibande, eta(ibande), delta_eta(ibande)
583 |
584 | c enddo
585 |
586 |
587 |
588 |
589 | c--------------------------------------------------------------------
590 | c pre - calcul sur les variables d'entree
591 | c....................................................................
592 | c
593 | c on a choisi un spectre independant de l'espace car h20~cte
594 | c donc calcul une fois sur toute les bandes
595 |
596 | c call calcparamspect(ptot_spectre,t_spectre,fm_h2o_spectre
597 | c & ,fm_co2_spectre ,fm_co_spectre)
598 |
599 |
600 |
601 |
602 |
603 |
604 |
605 |
606 | c**********************************************************************
607 | c fin des saisies debut du calcul
608 | c**********************************************************************
609 | call frouge
610 | print *
611 | print *,'EN COURS DE PROGRAMME'
612 | print *
613 | call freset
614 | call fligne
615 | c
616 | c.......................................................................
617 | c boucle sur les mailles : pour chaque maille un EMCM qui calcule
618 | c tous les echanges avec les autres mailles d'un coup (parcours deterministe)
619 | c.......................................................................
620 | c
621 | if (rsup .gt. 0.d+0) then
622 | print*, 'FORMULE: géométrie plane (rsup=',rsup,')'
623 | else
624 | print*, 'FORMULE: géométrie cylindrique'
625 | endif
626 | call fligne
627 |
628 | if (boolener) then
629 | if (boolprof) then
630 | print*,'FORMULE: OPTICO-GEOM: SNBk-dist ENER:description en lum'
631 | call fligne
632 | else
633 | print*,'FORMULE: OPTICO-GEOM: SNBk-dist ENER:description en temp'
634 | call fligne
635 | endif
636 | else
637 | print*,'FORMULE: OPTICO-GEOM: SNBk-dist ENER:néant'
638 | call fligne
639 | endif
640 |
641 | if (boolplat) then
642 | print*,'WARNING: on tire plat en angle solide (pour cas mince,
643 | & changer dans fichier source) '
644 | call fligne
645 | else
646 | print*,'WARNING: on tire en lambert pour les angles'
647 | call fligne
648 | endif
649 |
650 | if (boolinterieur) then
651 | if (boolssmaille) then
652 | print*, 'WARNING: profil sous maille en temperature (ou lum),
653 | & ANALYTIQUE'
654 | call fligne
655 | else
656 | print*, 'WARNING: profil sous maille en temperature (ou lum),
657 | & INDEXATION'
658 | call fligne
659 | endif
660 | else
661 | print*, 'WARNING: profil de temperature (ou lum) constant,
662 | & maille homogène'
663 | call fligne
664 | endif
665 |
666 | if (choixheterogeneite .eq. 'NO' ) then
667 | print*, 'WARNING: methode SEPARABILITE + mélange gaz/suie ok'
668 | hetero=1
669 | call fligne
670 | endif
671 | if (choixheterogeneite .eq. 'CG' ) then
672 | print*, 'WARNING: methode CURTIS-GODSON + mélange gaz/suie ok'
673 | hetero=2
674 | call fligne
675 | endif
676 | if (choixheterogeneite .eq. 'CK' ) then
677 | print*, 'WARNING: methode Correlated K + mélange gaz/suie ok'
678 | hetero=3
679 | call fligne
680 | endif
681 | print *, 'ALGO: pre-processing NO'; call fligne
682 | print *, 'ALGO: parallélisation NO'; call fligne
683 |
684 | print*, 'boucle sur les mailles'
685 | call fligne
686 | DO 10001 in=0,n+1
687 | !if (in.eq.4) then
688 | ! iseed = 736263595
689 | !endif
690 | call fligne !mieux de decider avant d'ecrire
691 | print *,'Numero de la maille en cours:',in,'/',n+1,' typeM: statM:
692 | & estimé: écoulé: restant: '
693 | !temps estimé= nbr total operation * temps écoulé/nbre op écoulées
694 | !temps écoulé OK
695 |
696 | c calcul des taubar et DELTAb pour optimisation du tirage de la bande--------
697 | c modif pour prendre en compte l'ambiance
698 | if (in .eq. 0) then
699 | tmp5 = (r(n+1)-r(0))* 1.66
700 | !tmp6=temp(1)
701 | tmp6 = 1500
702 | tmp7=fm(1,1)
703 | tmp8=fm(2,1)
704 | tmp9=fm(3,1)
705 | ! si je prends l'ambiance que voit la paroi gauche
706 | tmp7=moyco
707 | tmp8=moyco2
708 | tmp9=moyh2o
709 | tmp10=moyfv
710 | elseif (in.eq.(n+1)) then
711 | tmp5 = (r(n+1)-r(0))* 1.66
712 | !tmp6=temp(n)
713 | tmp6 = 1500
714 | tmp7=fm(1,n)
715 | tmp8=fm(2,n)
716 | tmp9=fm(3,n)
717 | ! ambiance pareil avec paroi droite
718 | tmp7=moyco
719 | tmp8=moyco2
720 | tmp9=moyh2o
721 | tmp10=moyfv
722 | else
723 | tmp5 = (r(in)-r(in-1))* 1.66
724 | tmp6=temp(in)
725 | tmp7=fm(1,in)
726 | tmp8=fm(2,in)
727 | tmp9=fm(3,in)
728 | tmp10=fv(in)
729 | !ambiance vue des volumes
730 | !tmp7=moyco
731 | !tmp8 = moyco2
732 | !tmp9 = moyh2o
733 | !tmp10 = moyfv
734 | endif
735 |
736 | !recherche si la maille est vide de gaz et suies
737 | !eviter ainsi un bug dans tirage de bandes (div par zero)
738 | if ((tmp7.eq.0.d+0).and.(tmp8.eq.0.d+0).and.(tmp9.eq.0.d+0)
739 | & .and.(fv(in).eq.0.d+0)) then
740 | vacuum = .True.
741 | else
742 | vacuum = .False.
743 | endif
744 |
745 |
746 | do ibande=1, nbande
747 | ! taubar
748 | iin =in
749 | call modbsuie(moyfv)
750 | call modbgaz(ptot,tmp6,tmp7,tmp8,tmp9,flag1)
751 | do igaz = 1, 3
752 | c la ligne suivante n'a rien a voir ces pour initialiser pour ck
753 | correc(igaz) = 1.d+0
754 | c if (( kgbar(igaz,in)*(r(n+1)-r(0)) ).le.1.d-5) then
755 | if ( (kgbar(igaz,in)*tmp5) .le.1.d-6) then
756 | taub(igaz,ibande)=1.d+0
757 | else
758 | call trss(phig(igaz,in),kgbar(igaz,in),tmp5,
759 | & taub(igaz,ibande))
760 | endif
761 | enddo
762 | taubs(ibande)=dexp(-ksbar(iin)*tmp5)
763 |
764 | tmp20=taub(1,ibande)*taub(2,ibande)*taub(3,ibande)
765 | & *taubs(ibande)
766 | if (tmp20 .le. (1.d+0-1.d-6)) then
767 | mtaub(ibande)=1.d+0 - tmp20
768 | else
769 | !mince 1-exp(-kl)= kl
770 | mtaub(ibande)=(kgbar(1,in)+kgbar(2,in)+kgbar(3,in)
771 | & +ksbar(iin))*tmp5
772 | endif
773 | skl(ibande) =(kgbar(1,in)+kgbar(2,in)+kgbar(3,in)
774 | & +ksbar(iin))*tmp5
775 |
776 |
777 | ! partie DELTAb le max del'ecart
778 | call black (eta(ibande)*1.d+2, temp(in), tmp1)
779 | call black (eta(ibande)*1.d+2, tmin, tmp2)
780 | call black (eta(ibande)*1.d+2, tmax, tmp3)
781 | tmp4 = abs(tmp2-tmp1)
782 | if (tmp4 .lt. (abs(tmp3-tmp1))) then
783 | dbmax(ibande) = abs(tmp3-tmp1)
784 | else
785 | dbmax(ibande) = tmp4
786 | endif
787 | enddo
788 | c------------------------------------------------------------------
789 | c
790 | c.......................................................................
791 | c boucle sur les nombres de tirage (parcours deterministe)
792 | c.......................................................................
793 | c
794 | call fligne
795 | !print *,'boucle sur nbr de tirages:', ntir(in)
796 | print *, ntir(in),' tirages: '
797 | DO 10002 itir=1,ntir(in)
798 |
799 | !itmp2 est une variable qui va de zero à 100 %
800 | !on choisi d'imprimer un point tous les 10%
801 | itmp2 = aint( itir /(ntir(in)/4.) * 25)+1
802 | if (mod(itmp2,10) .ne. 0) flag0=.True.
803 | if ((mod(itmp2,10) .eq. 0) .and. flag0) then
804 | if ((mod(itmp2,50) .eq. 0) .and. flag0) then
805 | call system('printf "%c" M ') !repere 50%
806 | else
807 | call system('printf "%c" . ')
808 | endif
809 | flag0=.False.
810 | endif
811 | ! write(text,*) itir
812 | ! affichage = 'printf "nbr de tirages: %f \b " '
813 | ! & //text
814 | ! print *, affichage
815 | ! call system(affichage)
816 |
817 | c
818 | c.......................................................................
819 | c poids pour la premiere fois initiale
820 |
821 | poids=2.D+0*pi*r(in)
822 | kbarphi_s(1) = 0.d+0
823 | kbar_s(1) =0.d+0
824 | kbarphi_s(2) = 0.d+0
825 | kbar_s(2) =0.d+0
826 | kbarphi_s(3) = 0.d+0
827 | kbar_s(3) =0.d+0
828 | c
829 | c ponderer la sommation qd integrale sur bande etroite
830 | c le 99 II 18
831 | if (boolspec) then
832 | poids=poids * nbande
833 | endif
834 | c on genere une bande de maniere identique paroi volume
835 | call genere_b(1,opt_bande,vacuum,proba,tran,total)
836 | !print *, 'ibande=', ibande
837 | c
838 | DO itmp1=1,n
839 | dpdksp(itmp1)=0.D+0
840 | ENDDO
841 | c
842 | c ---{P1 - angle - Lambert}
843 | c
844 | c suite a pb cas tres mince il faut tirer en lambert ou plat - dec 1999
845 | c deux etapes 1/ savoir tirer sur pdf_adaptee 2/modif du poids
846 | c 1/ pdf_adapte = alpha * pdf_plat + alpha * lambert
847 | c * trouver alpha pour savoir si on tire avec plat ou lambert
848 | c -etablir la loi de ponderation alpha=alpha critique
849 | alphac= 1.d+0 / (1.d+0 + skl(ibande))
850 | c print *, alphac
851 | c -tirer alpha puis test bernouilli
852 | call rand_uniforme(tmp3)
853 | !tmp3 < alphac on tire plat sinon lambert
854 | c * tirer sur la loi (apha ou lambert)
855 | CALL rand_uniforme(tmp1)
856 |
857 | if (boolplat) then !si vrai on tire plat tout le temps
858 | !!!!!!!!!if (skl(ibande).le. 1.d+0) then !on tire plat tout le temps
859 | !!!!!!!!!if (tmp3.le.alphac) then ! plat ou Lambert(pas de modif)
860 | !on tire uniforme en angle solide
861 | !cela revient a remplacer le tmp1 de Lambert par
862 | tmp1 = 1.d+0 -(tmp1*tmp1)
863 | ! tmp1= cos devient 1-cos2=sin2
864 | endif
865 |
866 | CALL rand_uniforme(tmp2)
867 | tmp2=tmp2*pi/2.D+0
868 | f=1.D+0 / dsqrt(1.D+0-tmp1*dcos(tmp2)**2)
869 | sin_teta=dsin(tmp2)*sqrt(tmp1) * f
870 | delta=r(in)*sin_teta
871 |
872 |
873 | c 2/ modif du poids consecutif a pdf_adapte en angulaire: on ne modifie que si tire plat
874 | if (boolplat) then !si vrai on tire plat tout le temps
875 | !!!!!!!!!if (skl(ibande).le. 1.d+0) then !on tire plat tout le temps
876 | !!!!!!!!!if (tmp3.le.alphac) then !plat ou Lambert
877 | !on tire uniforme en angle solide
878 | ! cos teta = sqrt(1-tmp1)
879 | !si Lambert et plat: poids = poids / ((1.d+0-tmp3)+(tmp3/(2*dsqrt(1.d+0-tmp1))))
880 | poids=poids * 2 *dsqrt(1.d+0-tmp1)
881 | endif
882 |
883 | !!!!!!!!!pour tirer tout droit
884 | !delta = 0.d+0
885 | !f = 1.d+0
886 | c
887 | c on veut pi*DELTA B !!!
888 | poids=poids*pi
889 | c---------------------------------------------------------------------------
890 | c ---{identification d'un chemin a partir de la maille d'emission}
891 | c en 1-D cylindrique sans reflexions (parois noires)
892 | c la determination de la plus petite maille est suffisante
893 | c ---{determination de la plus petite maille traversee}
894 | c
895 | iin_min=in
896 | 1001 CONTINUE
897 | IF (iin_min.GT.0) THEN
898 | IF (r(iin_min-1).GT.delta) THEN
899 | iin_min=iin_min-1
900 | GOTO 1001
901 | ENDIF
902 | ENDIF
903 |
904 |
905 |
906 | c
907 | ccccccccccccccccccccccc ---{cas particulier de la surface interne}
908 | c
909 | IF (in.EQ.0) THEN
910 | !if(itir .eq. 5161) print *, 'CAS 1: Surface INTERNE'
911 | c
912 | c -----{calcul des long necessaires : surface interne}
913 | c
914 | DO iin=0,n
915 | long(iin)=dsqrt(r(iin)**2-delta**2) * f
916 | ENDDO
917 |
918 | c
919 | c---{tirage de la bande etroite} en fonction maille in amaury juin98
920 | c
921 | c call genere_b(1,ibande)
922 |
923 | c calcul des proprietes radiatives sur la bande: les gaz
924 | c kgbar et phi pour cette maille et kgbar pour toute les autres
925 | do iin =1 ,n
926 | call modbgaz(ptot,temp(iin),fm(1,iin),fm(2,iin),
927 | & fm(3,iin),flag1)
928 | call modbsuie(fv(iin))
929 | c write (*,*) kgbar(iin),'phig',phig
930 | c read (*,*)
931 | if (.not. boolspec) then
932 | do igaz=1,3
933 | kgbar(igaz,iin) = speckgbar
934 | phig(igaz,iin) = specphi
935 | if (iin.eq.1) phig(igaz,iin) = specphi*2.
936 | enddo
937 | kgbar(1,iin)=0.d+0
938 | kgbar(2,iin)=0.d+0
939 | !reste que l'eau
940 | ksbar(iin) = 0.d+0
941 | !suie supprime
942 | endif
943 | !richard
944 | DO igaz=1,3
945 | IF (kgbar(igaz,iin).LE.1.d-20)
946 | & phig(igaz,in)=phieq(igaz,ibande)
947 | ENDDO
948 | enddo
949 | ccccccccccccoucou 1
950 | c kgbar(3)= 0.01D-99
951 | c on emet d'une parois mais on genere k dans la maille
952 | c de gaz a cote: ici on est en in =0 et on prend de l'eau
953 | c phig =phi(3,ibande)
954 | do igaz=1,3
955 | kgbar(igaz,0)=kgbar(igaz,0+1)
956 | !phig(igaz,0)=phig(igaz,0+1)
957 | !richard
958 | phig(igaz,0)=phieq(igaz,ibande)
959 | ksbar(0)=ksbar(0+1)
960 | enddo
961 | c sinon propriete radiative si dessous
962 |
963 | c calcul des proprietes radiatives sur la bande: les suies
964 | c pour l'instant rien lecture dans fichier
965 |
966 |
967 |
968 | c
969 | c -----{tirage de kg}
970 | c
971 |
972 | ccccccccccccccccc if ( kgbar(3,0) .eq. 0.d+0) then
973 | ccccccccccccccc kgskgbar(3) = 1
974 |
975 | ccccccccccccc else
976 | cccccc print *, 'titi'
977 | c pour les paroi dans l'epaisseur on met dim systeme
978 | CALL genere_kgaz(1,r(n+1)-r(0))
979 | cccccc print *, 'fin de genere_kgaz'
980 | cccccc print*, 'kgskgbar(1),kgskgbar(2),kgskgbar(3)'
981 | ccccc print *, kgskgbar(1),kgskgbar(2),kgskgbar(3)
982 | c read *
983 | ccccccccccccc endif
984 | c write (*,*) 'kgskgbar(3)',kgskgbar(3)
985 | c read(*,*)
986 |
987 | c
988 | c -----{calcul des k necessaires : surface interne}
989 | c
990 | !calcul de l'invariant maille d'emission
991 | do igaz=1,3
992 | call cdss(phig(igaz,in),1.d+0,0.d+0,kgskgbar(igaz),invarg(igaz))
993 | ! print *, 'invar',igaz, invarg(igaz),phig(igaz,in),
994 | ! & phig(igaz,iin),kgskgbar(igaz)
995 | enddo
996 | !---------------------------------------
997 | DO iin=1,n
998 | !inhmoog--------
999 | do igaz=1,3
1000 | ! print *, 'in,iin/ITIR/igaz-----------',in,iin,itir ,igaz
1001 | if (hetero.eq.3) call cor_tab_ck(invarg(igaz),phig(igaz,in)
1002 | & ,phig(igaz,iin)
1003 | & ,correc(igaz))
1004 | enddo
1005 | !---------------
1006 | k(iin)=1.D+0 * ksbar(iin)
1007 | & +kgskgbar(3)*kgbar(3,iin)*correc(3)
1008 | & +kgskgbar(2)*kgbar(2,iin)*correc(2)
1009 | & +kgskgbar(1)*kgbar(1,iin)*correc(1)
1010 | c write(*,*) k(iin)
1011 | c read(*,*)
1012 | ENDDO
1013 | c
1014 | c on est sur la surface interne: sigma1 = 0 pas d'épaissur
1015 | CALL trajet(1,boolener, boolprof,boolssmaille,boolinterieur)
1016 | c
1017 | GOTO 10003
1018 | ENDIF
1019 | c
1020 | cccccccccccccccccccccccccccccccccccc ---{cas particulier de la surface externe}
1021 | c
1022 | IF (in.EQ.n+1) THEN
1023 | !if(itir .eq. 5161) print *, 'CAS 2: Surface EXTERNE'
1024 | c
1025 | c -----{calcul des long necessaires : rayon vers l'interieur}
1026 | c
1027 | IF (iin_min.EQ.0) THEN
1028 | DO iin=0,n
1029 | long(iin)=dsqrt(r(iin)**2-delta**2) * f
1030 | ENDDO
1031 | ELSE
1032 | DO iin=iin_min,n
1033 | long(iin)=dsqrt(r(iin)**2-delta**2) * f
1034 | ENDDO
1035 | ENDIF
1036 |
1037 |
1038 | c
1039 | c---{tirage de la bande etroite} amaury juin98
1040 | c
1041 | c call genere_b(memoire)
1042 | c
1043 | c call genere_b(1,ibande)
1044 |
1045 |
1046 | c calcul des proprietes radiatives sur la bande: les gaz
1047 | c kgbar et phi pour cette maille et kgbar pour toute les autres
1048 | do iin =1 ,n
1049 | call modbgaz(ptot,temp(iin),fm(1,iin),fm(2,iin),
1050 | & fm(3,iin),flag1)
1051 | call modbsuie(fv(iin))
1052 | if (.not. boolspec) then
1053 | do igaz=1,3
1054 | kgbar(igaz,iin) = speckgbar
1055 | phig(igaz,iin) = specphi
1056 | if (iin.eq.1) phig(igaz,iin) = specphi*2.
1057 | enddo
1058 | kgbar(1,iin)=0.d+0
1059 | kgbar(2,iin)=0.d+0
1060 | !reste que l'eau
1061 | ksbar(iin) = 0.d+0
1062 | !suie supprime
1063 | endif
1064 | !richard
1065 | DO igaz=1,3
1066 | IF (kgbar(igaz,iin).LE.1.d-20)
1067 | & phig(igaz,in)=phieq(igaz,ibande)
1068 | ENDDO
1069 | enddo
1070 | ccccccccccccoucou 1
1071 | c kgbar(3)= 0.01D-99
1072 | c on emet d'une parois mais on genere k avec loi gaz
1073 | c phig =phi(3,ibande) et on prend de l'eau
1074 | do igaz=1,3
1075 | kgbar(igaz,n+1)=kgbar(igaz,n)
1076 | !phig(igaz,n+1)=phig(igaz,n)
1077 | !richard
1078 | phig(igaz,n+1)=phieq(igaz,ibande)
1079 | ksbar(n+1)=ksbar(n)
1080 | enddo
1081 | c sinon propriete radiative si dessous
1082 |
1083 | c calcul des proprietes radiatives sur la bande: les suies
1084 | c pour l'instant rien lecture dans fichier
1085 |
1086 |
1087 | c
1088 | c -----{tirage de kg}
1089 | c
1090 |
1091 | cccccccccccccccc if ( kgbar(3,n+1) .eq. 0.d+0) then
1092 | cccccccccccccccc kgskgbar(3) = 1
1093 |
1094 | cccccccccccccccccc else
1095 | c pour les parois dans l'epaisseur on met dim systeme
1096 | CALL genere_kgaz(1,r(n+1)-r(0))
1097 | cccccccccccccccccccc endif
1098 | c
1099 | c -----{calcul des k necessaires : rayon vers l'interieur}
1100 | c
1101 | IF (iin_min.EQ.0) THEN
1102 | !calcul de l'invariant maille d'emission
1103 | do igaz=1,3
1104 | call cdss(phig(igaz,in),1.d+0,0.d+0,kgskgbar(igaz),
1105 | & invarg(igaz))
1106 | enddo
1107 | !------------------------------
1108 | DO iin=1,n
1109 | !inhmoog--------
1110 | do igaz=1,3
1111 | ! print *, 'in,iin/ITIR/igaz--',in,iin,itir ,igaz
1112 | ! print *,phig(igaz,in),phig(igaz,iin),
1113 | ! & kgskgbar(igaz),invarg(igaz)
1114 | if (hetero.eq.3) call cor_tab_ck(invarg(igaz),phig(igaz,in)
1115 | & ,phig(igaz,iin)
1116 | & ,correc(igaz))
1117 | enddo
1118 | !---------------
1119 | k(iin)=1D+0*ksbar(iin)
1120 | & +kgskgbar(3)*kgbar(3,iin)*correc(3)
1121 | & +kgskgbar(2)*kgbar(2,iin)*correc(2)
1122 | & +kgskgbar(1)*kgbar(1,iin)*correc(1)
1123 |
1124 | c print *, 'rayons vers interieur'
1125 | c print *, k(iin)
1126 | c read *
1127 | ENDDO
1128 | ELSE
1129 | !calcul de l'invariant maille d'emission
1130 | do igaz=1,3
1131 | call cdss(phig(igaz,in),1.d+0,0.d+0,kgskgbar(igaz),
1132 | & invarg(igaz))
1133 | enddo
1134 | !----------------------------------------
1135 | DO iin=iin_min,n
1136 | !inhmoog--------
1137 | do igaz=1,3
1138 | ! print *, 'in,iin/ITIR/igaz--',in,iin,itir ,igaz
1139 | ! print *,phig(igaz,in),phig(igaz,iin),
1140 | ! & kgskgbar(igaz),invarg(igaz)
1141 | if (hetero.eq.3) call cor_tab_ck(invarg(igaz),phig(igaz,in)
1142 | & ,phig(igaz,iin)
1143 | & ,correc(igaz))
1144 | enddo
1145 | !---------------
1146 | k(iin)=1.D+0 * ksbar(iin)
1147 | & +kgskgbar(3)*kgbar(3,iin)*correc(3)
1148 | & +kgskgbar(2)*kgbar(2,iin)*correc(2)
1149 | & +kgskgbar(1)*kgbar(1,iin)*correc(1)
1150 | ENDDO
1151 | ENDIF
1152 | c
1153 | c on est sur la surface externe sigma_1 = 0 pas d'épaisseur
1154 | CALL trajet(2, boolener, boolprof,boolssmaille,boolinterieur)
1155 | c
1156 | GOTO 10003
1157 | ENDIF
1158 | c
1159 | cccccccccccccccccccccccccccccccccccccccccccc ---{les mailles de volume}
1160 | c
1161 | c ---{P1 - position - exponentielle}
1162 | c
1163 | IF (iin_min.EQ.in) THEN
1164 | !if(itir .eq. 5161) print *, 'CAS 3: RASANT (croissant tjrs)'
1165 | c
1166 | c -----{calcul des long necessaires : rayon rasant}
1167 | c
1168 | DO iin=iin_min,n
1169 | long(iin)=dsqrt(r(iin)**2-delta**2) * f
1170 | ENDDO
1171 | !if(itir .eq. 5161) print *, 'delate, f', delta, f
1172 | c
1173 | c---{tirage de la bande etroite} amaury juin98
1174 | c
1175 | c call genere_b(memoire)
1176 | c
1177 | c call genere_b(1,ibande)
1178 |
1179 | c calcul des proprietes radiatives sur la bande: les gaz
1180 | c kgbar et phi pour cette maille et kgbar pour toute les autres
1181 | do iin =1 ,n
1182 | call modbgaz(ptot,temp(iin),fm(1,iin),fm(2,iin),
1183 | & fm(3,iin),flag1)
1184 | call modbsuie(fv(iin))
1185 | c write (*,*) kgbar(iin),'phig',phig
1186 | c read (*,*)
1187 | if (.not. boolspec) then
1188 | do igaz=1,3
1189 | kgbar(igaz,iin) = speckgbar
1190 | phig(igaz,iin) = specphi
1191 | if (iin.eq.1) phig(igaz,iin) = specphi*2.
1192 | enddo
1193 | kgbar(1,iin)=0.d+0
1194 | kgbar(2,iin)=0.d+0
1195 | !reste que l'eau
1196 | ksbar(iin) = 0.d+0
1197 | !suie supprime
1198 | endif
1199 | !richard
1200 | DO igaz=1,3
1201 | IF (kgbar(igaz,iin).LE.1.d-20)
1202 | & phig(igaz,in)=phieq(igaz,ibande)
1203 | ENDDO
1204 | enddo
1205 | ccccccccccccoucou 1
1206 | c kgbar(3)=0.01D-99
1207 | if ((in.eq.0).or.(in.eq.(n+1))) then
1208 | c on emet d'une parois mais on genere k avec loi gaz
1209 | ccccccccc phig =0.1
1210 | ccccccccc kgbar(in)=1.
1211 | endif
1212 | c sinon propriete radiative si dessous
1213 |
1214 | c calcul des proprietes radiatives sur la bande: les suies
1215 | c pour l'instant rien lecture dans fichier
1216 |
1217 |
1218 |
1219 | c
1220 | c -----{tirage de kg}
1221 | c
1222 |
1223 | tmp1=2.D+0*long(in)
1224 |
1225 | cccccccccccccccccccc if ( kgbar(3,in) .eq. 0.d+0) then
1226 | ccccccccccccccccc kgskgbar(3) = 1
1227 |
1228 | cccccccccccccccccc else
1229 | c if(itir .eq. 5161) print *, 'poids avt kgaz', poids
1230 | CALL genere_kgaz(2,tmp1)
1231 | c if(itir .eq. 5161) print *, 'poids ape kgaz', poids
1232 | ccccccccccccccccc endif
1233 | c
1234 | c -----{calcul des k necessaires : rayon rasant}
1235 | c
1236 | c print *,'rayons rasants'
1237 | !calcul de l'invariant maille d'emission
1238 | do igaz=1,3
1239 | call cdss(phig(igaz,in),1.d+0,0.d+0,kgskgbar(igaz)
1240 | & ,invarg(igaz))
1241 | enddo
1242 | !------------------------------------------
1243 | DO iin=iin_min,n
1244 | !inhmoog--------
1245 | do igaz=1,3
1246 | ! print *, 'in,iin/ITIR/igaz--',in,iin,itir ,igaz
1247 | ! print *,phig(igaz,in),phig(igaz,iin),
1248 | ! & kgskgbar(igaz),invarg(igaz)
1249 | if (hetero.eq.3) call cor_tab_ck(invarg(igaz),phig(igaz,in)
1250 | & ,phig(igaz,iin)
1251 | & ,correc(igaz))
1252 | enddo
1253 | !---------------
1254 | k(iin)=1.D+0 * ksbar(iin)
1255 | & +kgskgbar(3)*kgbar(3,iin)*correc(3)
1256 | & +kgskgbar(2)*kgbar(2,iin)*correc(2)
1257 | & +kgskgbar(1)*kgbar(1,iin)*correc(1)
1258 | c print *, k(iin)
1259 | c read *
1260 | ENDDO
1261 | c
1262 | c -----{tirage de sigma_1}
1263 | c
1264 | call rand_uniforme(tmp1)
1265 | tmp2=2.D+0*long(in)
1266 | sigma_1=-dlog(1.D+0-tmp1*(1.D+0-dexp(-k(in)*tmp2)))
1267 | & /k(in)
1268 | ccc debug
1269 | if (k(in).eq. 0.d+0) then
1270 | sigma_1 = 0.d+0
1271 | endif
1272 | c
1273 | dpdksp(in)=dpdksp(in)
1274 | & + tmp2*dexp(-k(in)*tmp2)/(1.D+0-dexp(-k(in)*tmp2))
1275 | c if(itir .eq. 5161) print *, 'poids', poids
1276 | poids=poids*(1.D+0-dexp(-k(in)*tmp2))
1277 | ccc debug2
1278 | if (dexp(-k(in)*tmp3).eq. 1.d+0) then
1279 | dpdksp(in)= tmp3 !!! a revoir un delire du 13 mars 2000
1280 | endif
1281 | c
1282 | c if(itir .eq. 5161) print *, 'poids', poids
1283 | CALL trajet(1, boolener, boolprof,boolssmaille,boolinterieur)
1284 | c
1285 | ELSE
1286 | c
1287 | prob1=0.5D+0
1288 | prob2=1.D+0 - prob1
1289 | call rand_uniforme(tmp1)
1290 |
1291 | c
1292 | c -----{P1 - tirage d'un point proche de Q}
1293 | IF (tmp1.LT.prob1) THEN
1294 | !if(itir .eq. 5161)
1295 | !& print *, 'CAS 4: TRAVERSANT PROCHE de Q = rasant(croissant)'
1296 | c
1297 | poids=poids/prob1
1298 | c
1299 | c -------{calcul des long necessaires : rayon non rasant vers l'exterieur}
1300 | c
1301 | DO iin=in-1,n
1302 | long(iin)=dsqrt(r(iin)**2-delta**2) * f
1303 | ENDDO
1304 |
1305 |
1306 |
1307 |
1308 | c
1309 | c---{tirage de la bande etroite} amaury juin98
1310 | c
1311 | c call genere_b(memoire)
1312 | c
1313 | c if (in.eq.110) then
1314 | c print *, poids
1315 | c read *
1316 | c endif
1317 | cccc call genere_b(1,ibande)
1318 | c if (in.eq.110) then
1319 | c print *, poids
1320 | c read *
1321 | c endif
1322 |
1323 | c calcul des proprietes radiatives sur la bande: les gaz
1324 | c kgbar et phi pour cette maille et kgbar pour toute les autres
1325 | do iin =1 ,n
1326 | call modbgaz(ptot,temp(iin),fm(1,iin),fm(2,iin),
1327 | & fm(3,iin),flag1)
1328 | call modbsuie(fv(iin))
1329 | c write (*,*) kgbar(iin),'phig',phig
1330 | c read (*,*)
1331 | if (.not. boolspec) then
1332 | do igaz=1,3
1333 | kgbar(igaz,iin) = speckgbar
1334 | phig(igaz,iin) = specphi
1335 | if (iin.eq.1) phig(igaz,iin) = specphi*2.
1336 | enddo
1337 | kgbar(1,iin)=0.d+0
1338 | kgbar(2,iin)=0.d+0
1339 | !reste que l'eau
1340 | ksbar(iin) = 0.d+0
1341 | !suie supprime
1342 | endif
1343 | !richard
1344 | DO igaz=1,3
1345 | IF (kgbar(igaz,iin).LE.1.d-20)
1346 | & phig(igaz,in)=phieq(igaz,ibande)
1347 | ENDDO
1348 | enddo
1349 | ccccccccccccoucou 1
1350 | c kgbar(3)= 0.01D-99
1351 | if ((in.eq.0).or.(in.eq.(n+1))) then
1352 | c on emet d'une parois mais on genere k avec loi gaz
1353 | cccccccccccc phig =0.1
1354 | cccccccccccccc kgbar(in)=1.
1355 | endif
1356 | c sinon propriete radiative si dessous
1357 |
1358 | c calcul des proprietes radiatives sur la bande: les suies
1359 | c pour l'instant rien lecture dans fichier
1360 |
1361 |
1362 | c
1363 | c -------{tirage de kg}
1364 | c
1365 |
1366 |
1367 |
1368 | tmp1=(long(in)-long(in-1))
1369 |
1370 | cccccccccccccccccccc if ( kgbar(3,in) .eq. 0.d+0) then
1371 | ccccccccccccccccc kgskgbar(3) = 1
1372 |
1373 | cccccccccccccccc else
1374 | CALL genere_kgaz(2,tmp1)
1375 | ccccccccccccccccccccc endif
1376 | c
1377 | c -------{calcul des k necessaires : rayon non rasant vers l'exterieur}
1378 | c
1379 | IF (in.EQ.1) THEN
1380 | c print *, 'label!1'
1381 | !calcul de l'invariant maille d'emission
1382 | do igaz=1,3
1383 | call cdss(phig(igaz,in),1.d+0,0.d+0,kgskgbar(igaz)
1384 | & ,invarg(igaz))
1385 | enddo
1386 | !-------------------------------------------
1387 | DO iin=1,n
1388 | !inhmoog--------
1389 | do igaz=1,3
1390 | ! print *, 'in,iin/ITIR/igaz--',in,iin,itir ,igaz
1391 | ! print *,phig(igaz,in),phig(igaz,iin),
1392 | ! & kgskgbar(igaz),invarg(igaz)
1393 | if (hetero.eq.3) call cor_tab_ck(invarg(igaz),phig(igaz,in)
1394 | & ,phig(igaz,iin)
1395 | & ,correc(igaz))
1396 | enddo
1397 | !---------------
1398 | k(iin)=1.D+0 * ksbar(iin)
1399 | & +kgskgbar(3)*kgbar(3,iin)*correc(3)
1400 | & +kgskgbar(2)*kgbar(2,iin)*correc(2)
1401 | & +kgskgbar(1)*kgbar(1,iin)*correc(1)
1402 | c print *, k(iin)
1403 | c read *
1404 | ENDDO
1405 | ELSE
1406 | c print *, 'label2'
1407 | !calcul de l'invariant maille d'emission
1408 | do igaz=1,3
1409 | call cdss(phig(igaz,in),1.d+0,0.d+0,kgskgbar(igaz)
1410 | & ,invarg(igaz))
1411 | enddo
1412 | !------------------------------------------------
1413 | DO iin=in-1,n
1414 | !inhmoog--------
1415 | do igaz=1,3
1416 | ! print *, 'in,iin/ITIR/igaz--',in,iin,itir ,igaz
1417 | ! print *,phig(igaz,in),phig(igaz,iin),
1418 | ! & kgskgbar(igaz),invarg(igaz)
1419 | if (hetero.eq.3) call cor_tab_ck(invarg(igaz),phig(igaz,in)
1420 | & ,phig(igaz,iin)
1421 | & ,correc(igaz))
1422 | enddo
1423 | !---------------
1424 | k(iin)=1.D+0 * ksbar(iin)
1425 | & +kgskgbar(3)*kgbar(3,iin)*correc(3)
1426 | & +kgskgbar(2)*kgbar(2,iin)*correc(2)
1427 | & +kgskgbar(1)*kgbar(1,iin)*correc(1)
1428 | c print *, k(iin)
1429 | c read *
1430 | ENDDO
1431 | ENDIF
1432 | c
1433 | c -------{tirage de sigma_1}
1434 | c
1435 | call rand_uniforme(tmp2)
1436 | tmp3=(long(in)-long(in-1))
1437 | sigma_1=
1438 | & -dlog(1.D+0-tmp2*(1.D+0-dexp(-k(in)*tmp3)))
1439 | & /k(in)
1440 | ccc debug
1441 | if (k(in).eq. 0.d+0) then
1442 | sigma_1 = 0.d+0
1443 | endif
1444 |
1445 | c
1446 | dpdksp(in)=dpdksp(in)
1447 | & + tmp3*dexp(-k(in)*tmp3)/(1.D+0-dexp(-k(in)*tmp3))
1448 | poids=poids*(1.D+0-dexp(-k(in)*tmp3))
1449 | ccc debug2
1450 | if (dexp(-k(in)*tmp3).eq. 1.d+0) then
1451 | dpdksp(in)= tmp3 !!! a revoir un delire du 13 mars 2000
1452 | endif
1453 | c
1454 | CALL trajet(1, boolener, boolprof,boolssmaille,boolinterieur)
1455 | c
1456 | c -----{P1 - tirage d'un point loin de Q}
1457 | c
1458 | ELSE
1459 | !if(itir .eq. 5161)
1460 | !& print *,'CAS 5: TRAVERSANT Loin de Q = traversant reel'
1461 | poids=poids/prob2
1462 | c
1463 | c -------{calcul des long necessaires : rayon non rasant vers l'interieur}
1464 | c
1465 | IF (iin_min.EQ.0) THEN
1466 | DO iin=0,in
1467 | long(iin)=dsqrt(r(iin)**2-delta**2) * f
1468 | ENDDO
1469 | ELSE
1470 | DO iin=iin_min,n
1471 | long(iin)=dsqrt(r(iin)**2-delta**2) * f
1472 | ENDDO
1473 | ENDIF
1474 |
1475 | c
1476 | c---{tirage de la bande etroite} amaury juin98
1477 | c
1478 | c call genere_b(memoire)
1479 | c
1480 | cccc call genere_b(1,ibande)
1481 |
1482 |
1483 | c calcul des proprietes radiatives sur la bande: les gaz
1484 | c kgbar et phi pour cette maille et kgbar pour toute les autres
1485 | do iin =1 ,n
1486 | call modbgaz(ptot,temp(iin),fm(1,iin),fm(2,iin),
1487 | & fm(3,iin),flag1)
1488 | call modbsuie(fv(iin))
1489 | c write (*,*) kgbar(iin),'phig',phig
1490 | c read (*,*)
1491 | if (.not. boolspec) then
1492 | do igaz=1,3
1493 | kgbar(igaz,iin) = speckgbar
1494 | phig(igaz,iin) = specphi
1495 | if (iin.eq.1) phig(igaz,iin) = specphi*2.
1496 | enddo
1497 | kgbar(1,iin)=0.d+0
1498 | kgbar(2,iin)=0.d+0
1499 | !reste que l'eau
1500 | ksbar(iin) = 0.d+0
1501 | !suie supprime
1502 | endif
1503 | !richard
1504 | DO igaz=1,3
1505 | IF (kgbar(igaz,iin).LE.1.d-20)
1506 | & phig(igaz,in)=phieq(igaz,ibande)
1507 | ENDDO
1508 | enddo
1509 | ccccccccccccoucou 1
1510 | c kgbar(3)= 0.01D-99
1511 | if ((in.eq.0).or.(in.eq.(n+1))) then
1512 | c on emet d'une parois mais on genere k avec loi gaz
1513 | cccccccc phig =0.1
1514 | cccccccccc kgbar(in)=1.
1515 | endif
1516 | c sinon propriete radiative si dessous
1517 |
1518 | c calcul des proprietes radiatives sur la bande: les suies
1519 | c pour l'instant rien lecture dans fichier
1520 |
1521 |
1522 |
1523 |
1524 | c
1525 | c -------{tirage de kg}
1526 | c
1527 |
1528 | tmp1=(long(in)-long(in-1))
1529 |
1530 | ccccccccccccccccccc if ( kgbar(3,in) .eq. 0.d+0) then
1531 | cccccccccccccccccccccc kgskgbar(3) = 1
1532 |
1533 | ccccccccccccccc else
1534 | CALL genere_kgaz(2,tmp1)
1535 | cccccccccccccccc endif
1536 | c
1537 | c -------{calcul des k necessaires : rayon non rasant vers l'interieur}
1538 | c
1539 | IF (iin_min.EQ.0) THEN
1540 | !calcul de l'invariant maille d'emission
1541 | do igaz=1,3
1542 | call cdss(phig(igaz,in),1.d+0,0.d+0,kgskgbar(igaz)
1543 | & ,invarg(igaz))
1544 | enddo
1545 | !----------------------------------------------
1546 | DO iin=1,in
1547 | !inhmoog--------
1548 | do igaz=1,3
1549 | ! print *, 'in,iin/ITIR/igaz--',in,iin,itir ,igaz
1550 | ! print *,phig(igaz,in),phig(igaz,iin),
1551 | ! & kgskgbar(igaz),invarg(igaz)
1552 | if (hetero.eq.3) call cor_tab_ck(invarg(igaz),phig(igaz,in)
1553 | & ,phig(igaz,iin)
1554 | & ,correc(igaz))
1555 | enddo
1556 | !---------------
1557 | k(iin)=1.D+0 * ksbar(iin)
1558 | & +kgskgbar(3)*kgbar(3,iin)*correc(3)
1559 | & +kgskgbar(2)*kgbar(2,iin)*correc(2)
1560 | & +kgskgbar(1)*kgbar(1,iin)*correc(1)
1561 | ENDDO
1562 | ELSE
1563 | !calcul de l'invariant maille d'emission
1564 | do igaz=1,3
1565 | call cdss(phig(igaz,in),1.d+0,0.d+0,kgskgbar(igaz)
1566 | & ,invarg(igaz))
1567 | enddo
1568 | !___________________________________________
1569 | DO iin=iin_min,n
1570 | !inhmoog--------
1571 | do igaz=1,3
1572 | ! print *, 'in,iin/ITIR/igaz--',in,iin,itir ,igaz
1573 | ! print *,phig(igaz,in),phig(igaz,iin)
1574 | ! & ,kgskgbar(igaz),invarg(igaz)
1575 | if (hetero.eq.3) call cor_tab_ck(invarg(igaz),phig(igaz,in)
1576 | & ,phig(igaz,iin)
1577 | & ,correc(igaz))
1578 | enddo
1579 | !---------------
1580 | k(iin)=1.D+0 * ksbar(iin)
1581 | & +kgskgbar(3)*kgbar(3,iin)*correc(3)
1582 | & +kgskgbar(2)*kgbar(2,iin)*correc(2)
1583 | & +kgskgbar(1)*kgbar(1,iin)*correc(1)
1584 | ENDDO
1585 | ENDIF
1586 | c
1587 | c -------{tirage de sigma_1_star}
1588 | c
1589 | call rand_uniforme(tmp2)
1590 | tmp3=(long(in)-long(in-1))
1591 | sigma_1_star=
1592 | & -dlog(1.D+0-tmp2*(1.D+0-dexp(-k(in)*tmp3)))
1593 | & /k(in)
1594 | !en fait la distinction sig1 etsig1_star est purement
1595 | !inutil qd on connait le chemin venir ici sens 2
1596 | sigma_1 = sigma_1_star
1597 | c
1598 | ccc debug
1599 | if (k(in).eq. 0.d+0) then
1600 | sigma_1 = 0.d+0
1601 | endif
1602 |
1603 |
1604 | dpdksp(in)=dpdksp(in)
1605 | & + tmp3*dexp(-k(in)*tmp3)/(1.D+0-dexp(-k(in)*tmp3))
1606 | poids=poids*(1.D+0-dexp(-k(in)*tmp3))
1607 | ccc debug2
1608 | if (dexp(-k(in)*tmp3).eq. 1.d+0) then
1609 | dpdksp(in)= tmp3 !!! a revoir un delire du 13 mars 2000
1610 | endif
1611 | c
1612 | CALL trajet(2, boolener, boolprof,boolssmaille,boolinterieur)
1613 | c
1614 | ENDIF
1615 | ENDIF
1616 | c
1617 | c.......................................................................
1618 | c
1619 | 10003 CONTINUE
1620 | c
1621 | c.......................................................................
1622 |
1623 | c if (in.eq.3) then
1624 | c print *, in, ksbar(in), k(in)
1625 | c read (*,*)
1626 | c endif
1627 |
1628 | c.......................................................................
1629 | c
1630 | 10002 CONTINUE
1631 | c
1632 | c.......................................................................
1633 | c
1634 | 10001 CONTINUE
1635 | c
1636 | c***********************************************************************
1637 | c.......................................................................
1638 | c post-traitement pour le Monte Carlo
1639 | c.......................................................................
1640 | c
1641 | c bandes etroite a deplacer
1642 | c rajouter par amaury le 20/04/98
1643 | c precedemment pour une maille (in) et un tirage on a calcule 'ksi' ?
1644 | c il faut le multiplier par (delta nu) * (L - L) pour avoir la valeur
1645 | c energetique de l'echange.
1646 | c on fait une boucle sur tous les elements en regard surface ou volume
1647 | c apres ce traitement les psi sont algebriques et les ecart-types tjrs pos.
1648 | c donc j'ai mis valeur absolue
1649 | cccccccccc eta(1) = eta(1) *1.d+2
1650 | cccccccccc do in=0,n+1
1651 | cccccccccc call black (eta(1), temp(in),emi)
1652 | cccccccccc do iin=0,n+1
1653 | cccccccccc call black (eta(1), temp(iin),emj)
1654 | ccccccccccc il faut diviser les emittances par pi pour avoir luminances
1655 | cccccccccc psi(in,iin)= psi(in,iin)*delta_eta(1)*(emi-emj)/pi
1656 | cccccccccc var_psi(in,iin)= var_psi(in,iin)*(delta_eta(1)*
1657 | cccccccccc & (emi-emj)/pi)**2
1658 | c print *, 'emi','emj',emi, emj
1659 | c read (*,*)
1660 | cccccccccc do itmp1=1,n
1661 | cccccccccc dpsidk(in,iin,itmp1)=
1662 | cccccccccc & dpsidk(in,iin,itmp1)*delta_eta(1)*(emi-emj)/pi
1663 | cccccccccc var_dpsidk(in,iin,itmp1)=
1664 | cccccccccc & var_dpsidk(in,iin,itmp1)*(delta_eta(1)*(emi-emj)/pi)**2
1665 | cccccccccc dpsidkgbar(in,iin,itmp1)=
1666 | cccccccccc & dpsidkgbar(in,iin,itmp1)*delta_eta(1)*(emi-emj)/pi
1667 | cccccccccc var_dpsidkgbar(in,iin,itmp1)=
1668 | cccccccccc & var_dpsidkgbar(in,iin,itmp1)*(delta_eta(1)*(emi-emj)/pi)**2
1669 |
1670 | c print *, in,iin,itmp1, var_psi(in,iin),
1671 | c & var_dpsidk(in,iin,itmp1),var_dpsidkgbar(in,iin,itmp1)
1672 | cccccccccc enddo
1673 |
1674 | cccccccccc enddo
1675 | cccccccccc enddo
1676 | call fligne
1677 | call frouge
1678 | print *, 'POST-TRAITEMENT SPECIAL MCM:'
1679 | call freset
1680 | c....................................................................
1681 | c post traitement moyenne macro sur nbre tirages
1682 | c...................................................................
1683 | c ne marche que si psi et var_psi positif ?
1684 |
1685 | DO in=0,n+1
1686 | c traitement des rayons pour le cas simulation 1D plan parallele
1687 | c on le fait car tous les calculs sont termines !!!!!!!!!!!
1688 | r(in)=r(in)- rsup
1689 | DO iin=0,n+1
1690 | psi(in,iin)=psi(in,iin)/dble(ntir(in))
1691 | var_psi(in,iin)=var_psi(in,iin)/dble(ntir(in))
1692 | var_psi(in,iin)=var_psi(in,iin)-psi(in,iin)**2
1693 | var_psi(in,iin)=dsqrt(var_psi(in,iin)/dble(ntir(in)))
1694 | DO itmp1=1,n
1695 |
1696 | c dpsidk(in,iin,itmp1)=
1697 | c & dpsidk(in,iin,itmp1)/dble(ntir(in))
1698 | c var_dpsidk(in,iin,itmp1)=
1699 | c & var_dpsidk(in,iin,itmp1)/dble(ntir(in))
1700 | c var_dpsidk(in,iin,itmp1)=
1701 | c & var_dpsidk(in,iin,itmp1)-dpsidk(in,iin,itmp1)**2
1702 | c var_dpsidk(in,iin,itmp1)=
1703 | c & dsqrt(var_dpsidk(in,iin,itmp1)/dble(ntir(in)))
1704 |
1705 | dpsi6dfv(in,iin,itmp1)=
1706 | & dpsi6dfv(in,iin,itmp1)/dble(ntir(in))
1707 | cccc var_dpsi6dfv(in,iin,itmp1)=
1708 | cccc & var_dpsi6dfv(in,iin,itmp1)/dble(ntir(in))
1709 | cccc var_dpsi6dfv(in,iin,itmp1)=
1710 | cccc & var_dpsi6dfv(in,iin,itmp1)-dpsi6dfv(in,iin,itmp1)**2
1711 | cccc var_dpsi6dfv(in,iin,itmp1)=
1712 | cccc & dsqrt(var_dpsi6dfv(in,iin,itmp1)/dble(ntir(in)))
1713 |
1714 |
1715 | c dpsidkgbar(in,iin,itmp1)=
1716 | c & dpsidkgbar(in,iin,itmp1)/dble(ntir(in))
1717 | c var_dpsidkgbar(in,iin,itmp1)=
1718 | c & var_dpsidkgbar(in,iin,itmp1)/dble(ntir(in))
1719 | c var_dpsidkgbar(in,iin,itmp1)=
1720 | c & var_dpsidkgbar(in,iin,itmp1) - dpsidkgbar(in,iin,itmp1)**2
1721 | c var_dpsidkgbar(in,iin,itmp1)=
1722 | c & dsqrt(var_dpsidkgbar(in,iin,itmp1)/dble(ntir(in)))
1723 |
1724 | dpsi6dfmco(in,iin,itmp1)=
1725 | & dpsi6dfmco(in,iin,itmp1)/dble(ntir(in))
1726 | cccc var_dpsi6dfmco(in,iin,itmp1)=
1727 | cccc & var_dpsi6dfmco(in,iin,itmp1)/dble(ntir(in))
1728 | cccc var_dpsi6dfmco(in,iin,itmp1)=
1729 | cccc & var_dpsi6dfmco(in,iin,itmp1)-dpsi6dfmco(in,iin,itmp1)**2
1730 | cccc var_dpsi6dfmco(in,iin,itmp1)=
1731 | cccc & dsqrt(var_dpsi6dfmco(in,iin,itmp1)/dble(ntir(in)))
1732 | c
1733 | dpsi6dfmco2(in,iin,itmp1)=
1734 | & dpsi6dfmco2(in,iin,itmp1)/dble(ntir(in))
1735 | cccc var_dpsi6dfmco2(in,iin,itmp1)=
1736 | cccc & var_dpsi6dfmco2(in,iin,itmp1)/dble(ntir(in))
1737 | cccc var_dpsi6dfmco2(in,iin,itmp1)=
1738 | cccc & var_dpsi6dfmco2(in,iin,itmp1)-dpsi6dfmco2(in,iin,itmp1)**2
1739 | cccc var_dpsi6dfmco2(in,iin,itmp1)=
1740 | cccc & dsqrt(var_dpsi6dfmco2(in,iin,itmp1)/dble(ntir(in)))
1741 | c
1742 | dpsi6dfmh2o(in,iin,itmp1)=
1743 | & dpsi6dfmh2o(in,iin,itmp1)/dble(ntir(in))
1744 | cccc var_dpsi6dfmh2o(in,iin,itmp1)=
1745 | cccc & var_dpsi6dfmh2o(in,iin,itmp1)/dble(ntir(in))
1746 | cccc var_dpsi6dfmh2o(in,iin,itmp1)=
1747 | cccc & var_dpsi6dfmh2o(in,iin,itmp1)-dpsi6dfmh2o(in,iin,itmp1)**2
1748 | cccc var_dpsi6dfmh2o(in,iin,itmp1)=
1749 | cccc & dsqrt(var_dpsi6dfmh2o(in,iin,itmp1)/dble(ntir(in)))
1750 |
1751 |
1752 |
1753 | dpsi6dt(in,iin,itmp1)=
1754 | & dpsi6dt(in,iin,itmp1)/dble(ntir(in))
1755 | cccc var_dpsi6dt(in,iin,itmp1)=
1756 | cccc & var_dpsi6dt(in,iin,itmp1)/dble(ntir(in))
1757 | cccc var_dpsi6dt(in,iin,itmp1)=
1758 | cccc & var_dpsi6dt(in,iin,itmp1)-dpsi6dt(in,iin,itmp1)**2
1759 | cccc var_dpsi6dt(in,iin,itmp1)=
1760 | cccc & dsqrt(var_dpsi6dt(in,iin,itmp1)/dble(ntir(in)))
1761 |
1762 |
1763 |
1764 |
1765 | ENDDO
1766 | ENDDO
1767 | ENDDO
1768 | c
1769 | cc call densite
1770 |
1771 | c...................................................................
1772 | c post traitement tri selon ecart type
1773 | c...................................................................
1774 | c avant le traitement par luminances les ksi ont meme norme
1775 | c et meme signe
1776 | if (booltrie) then
1777 | DO in=0,n+1
1778 | c traitement des rayons pour le cas simulation 1D plan parallele
1779 | c on le fait car tous les calculs sont termines !!!!!!!!!!!
1780 | c r(in)=r(in)- rsup
1781 | DO iin=0,n+1
1782 | if (var_psi(in,iin).le.var_psi(iin,in)) then
1783 | psi(iin,in)= - psi(in,iin)
1784 | var_psi(iin,in)=var_psi(in,iin)
1785 | c DO itmp1=1,n
1786 | c var_dpsidk(iin,in,itmp1)
1787 | c & = var_dpsidk(in,iin,itmp1)
1788 | c var_dpsidkgbar(iin,in,itmp1)
1789 | c & = var_dpsidk(in,iin,itmp1)
1790 | cc dpsi6dfv(iin,in,itmp1)
1791 | cc & = - dpsid6dfv(in,iin,itmp1)
1792 | cc var_dpsidk(iin,in,itmp1)
1793 | cc & = var_dpsidk(in,iin,itmp1)
1794 | cc var_dpsidk(iin,in,itmp1)
1795 | cc & = var_dpsidk(in,iin,itmp1)
1796 | c enddo
1797 | else
1798 | psi(in,iin) = - psi(iin,in)
1799 | var_psi(in,iin)=var_psi(iin,in)
1800 | do itmp1=1,n
1801 | c var_dpsidk(in,iin,itmp1)
1802 | c & = var_dpsidk(iin,in,itmp1)
1803 | c var_dpsidkgbar(in,iin,itmp1)
1804 | c & = var_dpsidk(iin,in,itmp1)
1805 | enddo
1806 | endif
1807 | ENDDO
1808 | ENDDO
1809 | endif
1810 |
1811 |
1812 |
1813 | c
1814 | c***********************************************************************
1815 | c.......................................................................
1816 | c sorties
1817 | c.......................................................................
1818 | c
1819 | call frouge
1820 | print *
1821 | print *,'AFFICHAGE FINAL DES RESULTATS:'
1822 | call freset
1823 | call fligne
1824 | c
1825 | c.......................................................................
1826 | c sorties ecran
1827 | c.......................................................................
1828 | c
1829 | print *,'SORTIES ecran: formatage ascii et en petites colonnes'
1830 | print *
1831 | call fligne
1832 |
1833 | print *
1834 | print *,'sorties brutes'
1835 | print *
1836 | call fligne
1837 | c......les psi et ecart type: entre deux elements
1838 | cc print *, 'in iin psi ecart_type'
1839 | cc do in=0,n+1
1840 | cc do iin=0,n+1
1841 | cc print *, in, iin, psi(in,iin), var_psi(in,iin)
1842 | cc enddo
1843 | cc enddo
1844 | c.......les DERIVEES de psi par rapport KGAZ: ordre 1 et ecart type
1845 | cc print *, 'DpsiDkgaz: (Dpsi/dk(itmp1))(in,iin)'
1846 | cc do in=0,n+1
1847 | cc do iin=0,n+1
1848 | cc print *,'in iin *******',in,iin
1849 | cc do itmp1=1,n
1850 | cc print *, dpsidkgbar(in,iin,itmp1)
1851 | cc & ,var_dpsidkgbar(in,iin,itmp1)
1852 | cc enddo
1853 | cc enddo
1854 | cc enddo
1855 |
1856 | c.......les DERIVEES de psi par rapport KSUIE
1857 | cc print *, 'DpsiDksuie1: (Dpsi/dk(itmp1))(in,iin)'
1858 | cc do in=0,n+1
1859 | cc do iin=0,n+1
1860 | cc print *, 'in iin *******',in,iin
1861 | cc do itmp1=1,n
1862 | cc print *, dpsidk(in,iin,itmp1)
1863 | cc & ,var_dpsidk(in,iin,itmp1)
1864 | cc enddo
1865 | cc enddo
1866 | cc enddo
1867 |
1868 |
1869 |
1870 | c.......sortie derivees
1871 | PRINT *
1872 | print *,'Sortie derivees'
1873 | call fligne
1874 | ccc print *
1875 | ccc DO in=0,n+1
1876 | ccc print *, 'bilan et ecart type en ',in,':', bilan(in)
1877 | cccc &,var_bilan(in)
1878 | ccc ENDDO
1879 |
1880 | close(20)
1881 |
1882 | PRINT *
1883 | PRINT *
1884 | PRINT *,'en W/m2'
1885 | call fligne
1886 | PRINT *
1887 | PRINT *
1888 | c
1889 | cc DO in=0,n+1
1890 | cc DO iin=0,n+1
1891 | cc PRINT *,'psi',in,iin
1892 | cc & ,psi(in,iin)/2.D+0/pi/r(0)
1893 | cc & ,var_psi(in,iin)/2.D+0/pi/r(0)
1894 | cc PRINT *,'dpsidk2',in,iin
1895 | cc & ,dpsidk(in,iin,2)/2.D+0/pi/r(0)
1896 | cc & ,var_dpsidk(in,iin,2)/2.D+0/pi/r(0)
1897 | cc PRINT *,'dpsidkgbar2',in,iin
1898 | cc & ,dpsidkgbar(in,iin,2)/2.D+0/pi/r(0)
1899 | cc & ,var_dpsidkgbar(in,iin,2)/2.D+0/pi/r(0)
1900 | cc ENDDO
1901 | cc ENDDO
1902 |
1903 |
1904 | c.......................................................................
1905 | c sortie fichiers
1906 | c.......................................................................
1907 | c 1/ les donnees d'entree sont a regrouper en colonnes si le
1908 | c visualisateur le necessite ==> tot.in ???
1909 | open(10,FILE='SNBiseed',STATUS='unknown')
1910 | write(10,*) iseed
1911 | close(10)
1912 |
1913 | c 2/ la description des fichiers de sortie est dans cecile.out
1914 | c (sorties brutes et derivees)
1915 | call fich
1916 | open (10, file = 'proba_bande.out', status='unknown')
1917 | open (20, file = 'stat_bande.out', status='unknown')
1918 | do itmp4 = 1,nbande
1919 | write(10,*) itmp4, proba(itmp4)
1920 | write(20,*) itmp4, 0.d+0, 0.d+0 !tran(itmp4)/total, psi60bande(itmp4)/ntir(60)
1921 | enddo
1922 | close(10)
1923 | close(20)
1924 | c 2bis/ sortie pour curtis-godson
1925 | ! print *, 'paroi g: ',kgbar(3,0), phig(3,0)
1926 | ! print *, 'volume 1: ',kgbar(3,1), phig(3,1)
1927 | ! print *, 'volume 2: ',kgbar(3,2), phig(3,2)
1928 | ! print *, 'paroi d: ',kgbar(3,3), phig(3,3)
1929 |
1930 | ! tmp1 = r(1)-r(0)
1931 | ! tmp2 = r(2)-r(1)
1932 | ! tmp3 = r(2)-r(0)
1933 |
1934 | ! CALL trss(phig(3,2),kgbar(3,2),tmp2,taub2)
1935 |
1936 | ! kcg = (kgbar(3,1) * tmp1 + kgbar(3,2) * tmp2) / (tmp1 + tmp2)
1937 | ! phicg = (kgbar(3,1)*tmp1*phig(3,1) + kgbar(3,2)*phig(3,2)*tmp2)
1938 | ! & / (tmp1* kgbar(3,1) + tmp2*kgbar(3,2))
1939 |
1940 | ! CALL trss(phicg,kcg,tmp3,taub1)
1941 | ! PRINT *, 'taub2 -taub1: ', taub2 -taub1
1942 | ! PRINT *, 'psi 1 3(W) : ', (psi(1,3)/(25.D+2*2.d+0*pi*rsup))/pi ,
1943 | ! & (var_psi(1,3)/(25.D+2*2.d+0*pi*rsup))/pi
1944 | ! PRINT *, 'psi 3 1(W) : ', (psi(3,1)/(25.D+2*2.d+0*pi*rsup))/pi ,
1945 | ! & (var_psi(3,1)/(25.D+2*2.d+0*pi*rsup))/pi
1946 | ! PRINT *, '---'
1947 | ! PRINT *, '1-taub2 : ', 1- taub2
1948 |
1949 | ! PRINT *, 'psi 2 3(W) : ', (psi(2,3)/(25.D+2*2.d+0*pi*rsup))/pi,
1950 | ! & (var_psi(2,3)/(25.D+2*2.d+0*pi*rsup))/pi
1951 | ! PRINT *, 'psi 3 2(W) : ', (psi(3,2)/(25.D+2*2.d+0*pi*rsup))/pi,
1952 | ! & (var_psi(3,2)/(25.D+2*2.d+0*pi*rsup))/pi
1953 |
1954 |
1955 |
1956 |
1957 | c 3/ sorties des intermediaires de calcul (supplement)
1958 |
1959 | c *********sorties pour une maille: le spectre
1960 | c *********necessite de refaire un calcul deterministe !!!
1961 | ccc iin = 1
1962 | ccc ibande = 200
1963 | ccc open (10, file = 'sup_5conf.out', status='unknown')
1964 | ccc write (10,*) iin
1965 | ccc write (10,*) temp(iin)
1966 | ccc write (10,*) fm(1,iin)
1967 | ccc write (10,*) fm(2,iin)
1968 | ccc write (10,*) fm(3,iin)
1969 | ccc write (10,*) fv(iin)
1970 | ccc close(10)
1971 |
1972 | ccc open (40, file = 'sup_5spt_co.out', status='unknown')
1973 | ccc open (20, file = 'sup_5spt_co2.out', status='unknown')
1974 | ccc open (30, file = 'sup_5spt_h2o.out', status='unknown')
1975 | ccc do ibande=1,nbande
1976 | ccc call modbgaz(ptot,temp(iin),fm(1,iin),fm(2,iin), fm(3,iin)
1977 | ccc & ,flag1)
1978 |
1979 | ccc write (40,*) ibande, eta(ibande), delta_eta(in),
1980 | ccc & kgbar(1,iin),
1981 | ccc & gamma(1,ibande),
1982 | ccc & dinv(1,ibande),
1983 | ccc & phig(1,iin)
1984 |
1985 | ccc write (20,*) ibande, eta(ibande), delta_eta(in),
1986 | ccc & kgbar(2,iin),
1987 | ccc & gamma(2,ibande),
1988 | ccc & dinv(2,ibande),
1989 | ccc & phig(2,iin)
1990 |
1991 | ccc write (30,*) ibande, eta(ibande), delta_eta(in),
1992 | ccc & kgbar(3,iin),
1993 | ccc & gamma(3,ibande),
1994 | ccc & dinv(3,ibande),
1995 | ccc & phig(3,iin)
1996 |
1997 | ccc enddo
1998 |
1999 | ccc close(40)
2000 | ccc close(20)
2001 | ccc close(30)
2002 | c ***********sorties pour statistiques des tirages
2003 | call fpied
2004 |
2005 |
2006 | c
2007 | END