1 | c.......................................................................
2 | c
3 | c programme de traitement des donnees psi --> densites et bilans
4 | c comme majorant standard: on additionne les erreurs
5 | c ici on est passer au carre, c'est un peu plus petit
6 | c.......................................................................
7 |
8 | program gnutrait2paw
9 |
10 |
11 | implicit none
12 | c utilisation des variables du common, pas de passage de var
13 | include 'cecile.inc'
14 | integer itmp1, i,j,l, m
15 | double precision bil ,var_bil, psid
16 | double precision var_bilpg, var_bilpd, var_bilv
17 | double precision bilpg,bilpd, bilv
18 | double precision tempo, cum1, cum2, psifin, toto
19 | dimension bil(0:n_mx+1), var_bil(0:n_mx+1), psid(0:n_mx,0:n_mx)
20 | dimension bilpg(0:n_mx+1),bilpd(0:n_mx+1),bilv(0:n_mx+1)
21 | dimension var_bilpg(0:n_mx+1),var_bilpd(0:n_mx+1)
22 | & ,var_bilv(0:n_mx+1),psifin(0:n_mx+1)
23 | logical booltrie
24 | character*1 choix
25 |
26 |
27 | c la description des fichiers de sortie est dans cecile.out
28 | c print *, 'SORTIES fichiers: cecile.out'
29 | c print *
30 | c open (50,file='cecile.out',status='unknown')
31 | c write(50,*) 'Fichier de description des fichiers sorties'
32 | c write(50,*)
33 | c write(50,*) 'cecile.ecr: redirection de l''affichage ecran si b.'
34 | c write(50,*)
35 | c write(50,*) 'sorties ascii --> .out '
36 | c write(50,*) 'ind.out: Indices imax=n jmax=n kmax=n'
37 | c write(50,*) 'psi.out: Puissance nette echangees en W ---
38 | c & Ecart type pne'
39 | c write(50,*) 'dpsidk.out: DERIV/kgaz --- ecart --- DERIV /ksuie ---
40 | c & ecart'
41 | c write(50,*)
42 | c write(50,*) 'sortie binaires --> .oub'
43 | c close(50)
44 |
45 |
46 |
47 | c ***********************************************************************
48 | c ***********************************************************************
49 | c Sorties ascii: ACCES SEQUENTIEL
50 | c Leur format (commentaires ...) doit s'adapter a celui utilise par le
51 | c logiciel de relecture.
52 | c Faire une boucle sur chaque indice sauf si tableau mieux
53 |
54 |
55 | c ************************************************************************
56 | c LECTURE FORMAT GNUPLOT <=> sur une ligne: param,..,function (un param varie/bloc)
57 | open(20,file='ind.out', status='unknown')
58 | read(20,*) n, rsup
59 | read(20,*) n, rsup
60 | read(20,*) n, rsup
61 | close(20)
62 |
63 | c..................................................................................
64 | c open(20,file='psi.out',status='unknown', recl=21*(6))
65 | open(20,file='psi.out',status='unknown')
66 | c write(20,*) 'psi --- ecart type'
67 | cum1 = 0.d+0
68 | cum2 = 0.d+0
69 | do i = 0,n+1
70 | do j=0,n+1
71 | c read(20,'(i4,i4,4(e22.13))') in, iin, r(in), r(iin)
72 | read(20,*) in, iin, r(i), r(j)
73 | & ,psi(i,j), var_psi(i,j)
74 | c
75 | cc if (i .eq. (n+1)) then
76 | cc psi(121,0)= 0.d+0
77 | cc var_psi(121,0)= 0.d+0
78 | cc cum1 = cum1 + psi(i, j)
79 | cc cum2 = cum2 + var_psi(i,j)
80 | cc endif
81 | c
82 | enddo
83 | read(20,*)
84 | enddo
85 | close (20)
86 | cc cum1 = cum1 / (2 * pi * r(121))
87 | cc cum2 = cum2 / (2 * pi * r(121))
88 | cc print *, 'psi', i, ' ',cum1
89 | cc print *, cum2
90 | cc read *
91 |
92 | c rendre artificiellemt mauvais les psi retenus en emission et ceux en absortion
93 | c tjours aussi bien : test pour trie apres
94 | cc do i=0, n
95 | cc do j=i+1, n
96 | cc psi(i,j)=psi(i,j)+ 100.d+0
97 | cc enddo
98 | cc enddo
99 |
100 | c special triche pour debug
101 | c do i=0,n+1
102 | c psi(111,i)= -psi(i,111)
103 | c var_psi(111,i)= var_psi(i,111)
104 | c psi(110,i)= -psi(i,110)
105 | c var_psi(110,i)= var_psi(i,110)
106 | c psi(1,i)= -psi(i,1)
107 | c var_psi(1,i)= var_psi(i,1)
108 | c psi(0,i)= -psi(i,0)
109 | c var_psi(0,i)= var_psi(i,0)
110 | c enddo
111 |
112 | c...................................................................................
113 | c open(20,file='dpsidk.out',status='unknown',recl=21*(7))
114 | open(20,file='dpsi.out',status='unknown')
115 | c write(20,*)'DERIv psi/kgaz --- Ecart type --- DERIV psi/ksuie ---'
116 | do i=0,n+1
117 | do j=0,n+1
118 | do l=1,n
119 | c read(20,'(3(i4),4(e22.13))') in, iin,
120 | c & itmp1, dpsidkgbar(i,j,l), var_dpsidkgbar(i,j,l),
121 | c & dpsidk(i,j,l), var_dpsidk(i,j,l)
122 |
123 | read(20,'(3(i4),2(e22.13))') in, iin,
124 | & itmp1, dpsi6dfv(i,j,l)
125 | cccc & , var_dpsi6dfv(i,j,l)
126 |
127 | enddo
128 | read(20,*)
129 | enddo
130 | enddo
131 | close(20)
132 | c......................................................................
133 | c......................................................................
134 | c le fichier dpsidfm a ete rajoute le 14 novembre 1999
135 | open(20,file='dpsidfm.out',
136 | & status='unknown')
137 | do in=0,n+1
138 | do iin=0,n+1
139 | do itmp1=1,n
140 | c write(20,'(3(i4),2(e22.13))') in, iin,
141 | read(20,'(6(e22.13))')
142 | & dpsi6dfmco(in,iin,itmp1),
143 | cccc & var_dpsi6dfmco(in,iin,itmp1),
144 | & dpsi6dfmco2(in,iin,itmp1),
145 | cccc & var_dpsi6dfmco2(in,iin,itmp1),
146 | & dpsi6dfmh2o(in,iin,itmp1)
147 | cccc & , var_dpsi6dfmh2o(in,iin,itmp1)
148 | enddo
149 | enddo
150 | enddo
151 | close(20)
152 |
153 | c.....................................................................................
154 | c open(20,file='dpsi.out',status='unknown',recl=21*(7))
155 | open(20,file='dpsidt.out'
156 | & ,status
157 | & ='unknown')
158 | c write(20,*)'DERIv psi/kgaz --- Ecart type --- DERIV psi/ksuie ---'
159 | do i=0,n+1
160 | do j=0,n+1
161 | do m=1,n
162 | c read(20,'(3(i4),2(e22.13))') in, iin,
163 | read(20,*) in, iin,
164 | c write(20,*) in, iin,
165 | & itmp1, dpsi6dt(i,j,m)
166 | cccc & , var_dpsi6dt(i,j,m)
167 | enddo
168 | c simple retour a la ligne
169 | read(20,*)
170 | enddo
171 | enddo
172 | close(20)
173 | c..........................
174 | ! open(20,file='psibilfin.out'
175 | ! & ,status
176 | ! & ='unknown')
177 | ! do i=0,n+1
178 | ! read(20,*) in, toto, psifin(i)
179 | ! enddo
180 | ! close(20)
181 | c......................................................................
182 | c......................................................................
183 | c sortie binaire: ACCES DIRECT (plus universel)
184 | c cas ou les sorties ascii ne sont pas bien mise en forme et ou il faut
185 | c tout lire rapidement
186 | c chaque ligne est un enregistrement qu'on peut lire d'un bloc
187 | c enregistrement depend de la machine: sun double precision ici
188 | c on peut lire les lignes dans le desordre
189 |
190 | c numero ligne ou on enregistre moins un
191 | ccc itmp5= 0
192 | ccc open(20,file='psi.oub',access='direct',recl=8
193 | ccc &*(n+1+1))
194 | ccc write(20,*) 'Matrice des psi'
195 | ccc do in=0,n+1
196 | ccc write(20,rec=in+1) (psi(in,iin), iin=0,n+1)
197 | ccc enddo
198 | c derniere valeur de rec
199 | ccc itmp5=in+1
200 |
201 | c write(20,*)'Matrice des ecarttypes de psi')
202 | ccc do in=0+itmp5,n+1+itmp5
203 | ccc write(20,rec=in+1) (var_psi(in,iin), iin=0,n+1)
204 | ccc enddo
205 |
206 | c
207 | c gerer bloc acces direct connitre numero de ligne
208 | c
209 | ccc close(20)
210 |
211 |
212 |
213 | c ************************************************************************
214 | c TRAITEMENT des DONNEES
215 | open(10, file='ind.oup', status='unknown')
216 | write(10,*) n, rsup
217 | write(10,*) n, rsup
218 | write(10,*) n, rsup
219 | close(10)
220 |
221 |
222 | c-------------------------------------------------------------------------
223 | c...................................................................
224 | c post traitement tri selon ecart type
225 | c...................................................................
226 | c avant le traitement par luminances les ksi ont meme norme
227 | c et meme signe
228 | write(*,*) 'Voulez-vous le trie par ecart-type (n/y)?'
229 | cccccccc read (*,'(a)') choix
230 | write(*,*) 'Oui par defaut'
231 | choix = 'y'
232 | if (choix.eq.'y') then
233 | booltrie = .True.
234 | else
235 | booltrie =.False.
236 | endif
237 |
238 | if (booltrie) then
239 | DO in=0,n+1
240 | DO iin=0,n+1
241 | if (var_psi(in,iin).le.var_psi(iin,in)) then
242 | psi(iin,in)= - psi(in,iin)
243 | var_psi(iin,in)=var_psi(in,iin)
244 | DO itmp1=1,n
245 | c var_dpsidk(iin,in,itmp1)
246 | c & = var_dpsidk(in,iin,itmp1)
247 | c var_dpsidkgbar(iin,in,itmp1)
248 | c & = var_dpsidk(in,iin,itmp1)
249 | dpsi6dfv(iin,in,itmp1)
250 | & = - dpsi6dfv(in,iin,itmp1)
251 | dpsi6dfmco(iin,in,itmp1)
252 | & = - dpsi6dfmco(in,iin,itmp1)
253 | dpsi6dfmco2(iin,in,itmp1)
254 | & = - dpsi6dfmco2(in,iin,itmp1)
255 | dpsi6dfmh2o(iin,in,itmp1)
256 | & = - dpsi6dfmh2o(in,iin,itmp1)
257 | dpsi6dt(iin,in,itmp1)
258 | & = - dpsi6dt(in,iin,itmp1)
259 | enddo
260 | else
261 | psi(in,iin) = - psi(iin,in)
262 | var_psi(in,iin)=var_psi(iin,in)
263 | do itmp1=1,n
264 | c var_dpsidk(in,iin,itmp1)
265 | c & = var_dpsidk(iin,in,itmp1)
266 | c var_dpsidkgbar(in,iin,itmp1)
267 | c & = var_dpsidk(iin,in,itmp1)
268 | dpsi6dfv(in,iin,itmp1)
269 | & = - dpsi6dfv(iin,in,itmp1)
270 | dpsi6dfmco(in,iin,itmp1)
271 | & = - dpsi6dfmco(iin,in,itmp1)
272 | dpsi6dfmco2(in,iin,itmp1)
273 | & = - dpsi6dfmco2(iin,in,itmp1)
274 | dpsi6dfmh2o(in,iin,itmp1)
275 | & = - dpsi6dfmh2o(iin,in,itmp1)
276 | dpsi6dt(in,iin,itmp1)
277 | & = - dpsi6dt(iin,in,itmp1)
278 | enddo
279 | endif
280 | ENDDO
281 | ENDDO
282 | endif
283 | c-------------------------------------------------------------------------
284 | !rajoute le 17 fev 2000 --- Amaury
285 | if (choix.eq.'y') then
286 | write(*,*) 'Vous avez choisi un tri par ecart type: ok,
287 | & la moitie de l"information est perdue...'
288 |
289 | else
290 | write(*, *) 'Sans tri il y a deux logiques:'
291 | write(*, *) ' logique d"emission: retenir les
292 | & echanges qui partent d"ou je suis --> e'
293 | write(*, *) ' logique en reception: retenir les
294 | & echanges qui arrivent ou je suis --> a'
295 | c read (*,'(a)') choix
296 | write (*,*) 'Emission par defaut'
297 | choix = 'e'
298 | if (choix.eq.'e') then
299 | write (*, *) ' bilan(i)= Sj psi(i,j)'
300 | else
301 | write (*,*) ' bilan(j)= Si psi(i,j)'
302 | !pour effectuer cette deuxieme etape
303 | ! 1/ il faut se rappeler on a des psi(emis, absorb)
304 | ! 2/ le deuxieme cas = deuxieme_som psi(emis, absorb)
305 | ! ou bien = permiere_som psi(absorb, emis)
306 | ! donc pour garder une meme suite au code
307 | ! je fais les modifs suivantes pour inverser les positions
308 | ! et avoir des psi dont premier indice = maille absorption
309 | do i=0,n+1
310 | do j=0,n+1
311 | tempo = psi(i,j)
312 | ! nvelle affectations
313 | psi(i,j) = - psi(j,i)
314 | psi(j,i) = - tempo
315 | enddo
316 | enddo
317 | endif
318 | endif
319 | c-------------------------------------------------------------------------
320 | cc write (*,*) psi(50,55), psi(55,50)
321 | cc write (*,*) dpsi6dfmco2(5,55,30), dpsi6dfmco2(55,5,30)
322 | cc write (*,*) dpsi6dfmh2o(5,55,30), dpsi6dfmh2o(55,5,30)
323 | c write (*,*) dpsi6dt(50,55,55), dpsi6dt(55,50,55)
324 | cc write (*,*) dpsi6dt(1,44,1), dpsi6dt(44,1,1)
325 | cc write(*,*)
326 | cc write (*,*) var_psi(50,55), var_psi(55,50)
327 | cc read *
328 |
329 |
330 | write(*,*) 'MENU TRAITEMENT DES DONNEES'
331 | write(*,*) '___________________________'
332 | write(*,*) 'Densite de puissance nette echangees: '
333 |
334 | !densite
335 | write(*,*) 'MAILLE Volume/Volume [W/m origine == W]'
336 | write(*,*) 'densite Volume/Volume [W/m6]'
337 | !initialisation
338 | do i=0,n+1
339 | r(i)=r(i)+rsup
340 | ! toutes les divisions utilises r(i) donc
341 | if (r(i).eq.0.d+0) then
342 | r(i)= 1.d-100
343 | endif
344 | do j=0,n+1
345 | psid(i,j)=psi(i,j)
346 | enddo
347 | enddo
348 |
349 | do i=1,n
350 | do j=1,n
351 | !division par le volume origine i (1 metre de haut)
352 | psid(i,j)=psid(i,j)/((pi*r(i)**2)-(pi*r(i-1)**2))*1
353 | !division par un metre du volume de hauteur infinie
354 | psid(i,j)=psid(i,j)/((pi*r(j)**2)-(pi*r(j-1)**2))*1
355 | enddo
356 | enddo
357 |
358 |
359 | write(*,*) 'MAILLE Volume/Surface [W/m origine == W]'
360 | write(*,*) 'densite Volume/Surface [W/m5]'
361 | ! pas de boucle sur les parois
362 | write(*,*) 'paroi i=0'
363 | i=0
364 | do j=1,n
365 | !division par surface maille i d'un metre de haut
366 | psid(i,j)=psid(i,j)/(2*pi*r(i)*1)
367 | psid(j,i)=psi(j,i)/(2*pi*r(i)*1)
368 | !division par le volume de hauteur infinie j
369 | psid(i,j)=psid(i,j)/((pi*r(j)**2-(pi*r(j-1)**2))*1)
370 | psid(j,i)=psid(j,i)/((pi*r(j)**2-(pi*r(j-1)**2))*1)
371 | enddo
372 | write(*,*) 'paroi i=',n+1
373 | i=n+1
374 | do j=1,n
375 | !division par surface maille i d'un metre de haut
376 | psid(i,j)=psid(i,j)/(2*pi*r(i)*1)
377 | psid(j,i)=psid(j,i)/(2*pi*r(i)*1)
378 | !division par le volume de hauteur infinie j
379 | psid(i,j)=psid(i,j)/((pi*r(j)**2-(pi*r(j-1)**2))*1)
380 | psid(j,i)=psid(j,i)/((pi*r(j)**2-(pi*r(j-1)**2))*1)
381 | enddo
382 |
383 |
384 | write(*,*) 'MAILLE Surface/Surface [W/m origine == W]'
385 | write(*,*) 'densite Surface/Surface [W/m4]'
386 | !Division par la suface de un metre de haut
387 | psid(0,n+1)=psid(0,n+1)/(2*pi*r(0)*1)
388 | psid(0,n+1)=psid(0,n+1)/(2*pi*r(n+1)*1)
389 |
390 | psid(n+1,0)=psid(n+1,0)/(2*pi*r(n+1)*1)
391 | psid(n+1,0)=psid(n+1,0)/(2*pi*r(0)*1)
392 |
393 |
394 | !______________________________________________________
395 | ! format paw
396 | write (*,*) 'les matrices d"echange sont limitees a vol/vol'
397 | open(10, file='psid.oup', status='unknown')
398 | open(11, file='psitout.oup', status='unknown')
399 | do i = 1,n
400 | do j=1,n
401 |
402 | c !enlever tout ce qui concerne les parois
403 | c if (.not.(i.eq.0 .or. i.eq.(n+1).or. j.eq.0 .or.
404 | c & j.eq.(n+1))) then
405 |
406 |
407 |
408 | c write(10,*) i,j, real(r(i)),real( r(j))
409 | c & ,real(psid(i,j)), real(var_psi(i,j))
410 | write(10,*) real(psid(i,j))
411 | write(11,*) real(psi(i,j))
412 |
413 | if ((i.eq.3).and.(j.eq.7)) then
414 | print *, '3,7: ',psid(i,j), real(psid(i,j))
415 | endif
416 | if ((i.eq.7).and.(j.eq.3)) then
417 | print *, '7,3: ',psid(i,j), real(psid(i,j))
418 | endif
419 |
420 |
421 |
422 |
423 |
424 |
425 |
426 |
427 | c endif
428 |
429 |
430 | enddo
431 | c write(10,*)
432 | enddo
433 | close (10)
434 | close (11)
435 |
436 |
437 | !______________________________________________________
438 |
439 |
440 |
441 |
442 | write(*,*) '______________________________'
443 | write(*,*) 'Bilan de puissance par maille: '
444 | !bilan
445 | write(*,*) 'BILAN en [W/m2 hauteur]'
446 | do i=0,n+1
447 | bil(i)=0.d+0
448 | bilpg(i)=0.d+0
449 | bilpd(i)=0.d+0
450 | bilv(i)= 0.d+0
451 | var_bil(i)=0.d+0
452 | var_bilv(i)=0.d+0
453 | var_bilpg(i)=0.d+0
454 | var_bilpd(i)=0.d+0
455 | do j=0,n+1
456 | !inversion indice
457 | !psi(i,j)=-psi(j,i)
458 | !var_psi(i,j)=var_psi(j,i)
459 | bil(i)=bil(i)+(-psi(i,j))/1
460 | c var_bil(i)=dsqrt(var_bil(i)**2+abs(var_psi(i,j))**2)
461 | c & /1.d+0 *1.d+100
462 | !var_bil(i)=dsqrt(var_bil(i)**2+abs(var_psi(i,j))**2)
463 | var_bil(i)=var_bil(i)+abs(var_psi(i,j))
464 |
465 | if (j.eq.0) then !.or.(j.eq.(n+1))) then
466 | bilpg(i)=bilpg(i)+(-psi(i,j))
467 | !var_bilpg(i)=dsqrt(var_bilpg(i)**2+abs(var_psi(i,j))**2)
468 | var_bilpg(i)=var_bilpg(i)+abs(var_psi(i,j))
469 |
470 | elseif (j.eq.(n+1) ) then
471 | bilpd(i)=bilpd(i)+(-psi(i,j))
472 | !var_bilpd(i)=dsqrt(var_bilpd(i)**2+abs(var_psi(i,j))**2)
473 | var_bilpd(i)=var_bilpd(i)+abs(var_psi(i,j))
474 |
475 | else
476 |
477 | bilv(i)=bilv(i)+(-psi(i,j))
478 | !var_bilv(i)=dsqrt(var_bilv(i)**2+abs(var_psi(i,j))**2)
479 | var_bilv(i)=var_bilv(i)+abs(var_psi(i,j))
480 | endif
481 | enddo
482 | enddo
483 |
484 |
485 | write(*,*) 'BILAN en [W/m2 hauteur /vol ou surf ]'
486 | do i=0,n+1
487 | if ((i.ne.0).and.(i.ne.(n+1))) then
488 |
489 |
490 | bil(i) = bil(i)/((pi*r(i)**2)-(pi*r(i-1)**2))*1
491 | bilpg(i) = bilpg(i)/((pi*r(i)**2)-(pi*r(i-1)**2))*1
492 | bilpd(i) = bilpd(i)/((pi*r(i)**2)-(pi*r(i-1)**2))*1
493 | bilv(i) = bilv(i)/((pi*r(i)**2)-(pi*r(i-1)**2))*1
494 |
495 | var_bil(i)=abs(var_bil(i))/((pi*r(i)**2)-
496 | & (pi*r(i-1)**2))*1
497 |
498 | var_bilpg(i)=abs(var_bilpg(i))/((pi*r(i)**2)-
499 | & (pi*r(i-1)**2))*1
500 | var_bilpd(i)=abs(var_bilpd(i))/((pi*r(i)**2)-
501 | & (pi*r(i-1)**2))*1
502 | var_bilv(i)=abs(var_bilv(i))/((pi*r(i)**2)-
503 | & (pi*r(i-1)**2))*1
504 |
505 |
506 |
507 | !bil(i) = bil(i) * 1.d-100
508 | !var_bil(i) = var_bil(i) * 1.d-100
509 |
510 | !if (i.eq.(n+1))then
511 | ! print *,'i ',i
512 | ! print *,bilpg(i)
513 | ! print *,bilpd(i)
514 | !!!print *,'psi paroi',psi(54,0), psi(54,111)
515 | ! print *, bilv(i)
516 | ! print *,bil(i)*1.d-100
517 | ! read *
518 | !endif
519 | else
520 |
521 |
522 | bil(i) = bil(i) /(2*pi*r(i) *1)
523 | bilpg(i) = bilpg(i)/(2*pi*r(i) *1)
524 | bilpd(i) = bilpd(i)/(2*pi*r(i) *1)
525 | bilv(i) = bilv(i)/(2*pi*r(i) *1)
526 |
527 | var_bil(i)=abs(var_bil(i))/(2*pi*r(i) *1)
528 | var_bilpg(i)=abs(var_bilpg(i))/(2*pi*r(i) *1)
529 | var_bilpd(i)=abs(var_bilpd(i))/(2*pi*r(i) *1)
530 | var_bilv(i)=abs(var_bilv(i))/(2*pi*r(i) *1)
531 |
532 | !bil(i) = bil(i) * 1.d-100
533 | !var_bil(i) = var_bil(i) * 1.d-100
534 |
535 | !! print *,'i ',i
536 | !! print *,bilpg(i)
537 | !! print *,bilpd(i)
538 | !!!!print *,'psi paroi',psi(54,0), psi(54,111)
539 | !! print *, bilv(i)
540 | !! print *,bil(i)
541 | !! read *
542 |
543 |
544 | endif
545 |
546 | enddo
547 | !bilan de puissance calcule en partant des parois
548 | !c'est le meme qu'avant qd il y a le trie
549 | c**************************************************************************
550 | c do j=0,n+1
551 | c bil(i)=0.d+0
552 | c bilpg(i)=0.d+0
553 | c bilpd(i)=0.d+0
554 | c bilv(i)= 0.d+0
555 | c var_bil(i)=0.d+0
556 | c var_bilv(i)=0.d+0
557 | c var_bilpg(i)=0.d+0
558 | c var_bilpd(i)=0.d+0
559 | c do i=0,n+1
560 | !inversion indice
561 | !psi(i,j)=-psi(j,i)
562 | !var_psi(i,j)=var_psi(j,i)
563 |
564 | c bil(j)=bil(j)+(-psi(i,j))/1
565 | c var_bil(j)=dsqrt(var_bil(j)**2+abs(var_psi(i,j))**2)/1
566 |
567 | c if (i.eq.0) then !.or.(j.eq.(n+1))) then
568 | c bilpg(j)=bilpg(j)+(-psi(i,j))
569 | c var_bilpg(j)=dsqrt(var_bilpg(j)**2+abs(var_psi(i,j))**2)
570 |
571 | c elseif (i.eq.(n+1) ) then
572 | c bilpd(j)=bilpd(j)+(-psi(i,j))
573 | c var_bilpd(j)=dsqrt(var_bilpd(j)**2+abs(var_psi(i,j))**2)
574 |
575 | c else
576 |
577 | c bilv(j)=bilv(j)+(-psi(i,j))
578 | c var_bilv(j)=dsqrt(var_bilv(j)**2+abs(var_psi(i,j))**2)
579 | c endif
580 | c enddo
581 | c enddo
582 |
583 |
584 | c write(*,*) 'BILAN en [W/m2 hauteur /vol ou surf ]'
585 | c do i=0,n+1
586 | c if ((i.ne.0).and.(i.ne.(n+1))) then
587 |
588 |
589 | c bil(i) = bil(i)/((pi*r(i)**2)-(pi*r(i-1)**2))*1
590 | c bilpg(i) = -bilpg(i)/((pi*r(i)**2)-(pi*r(i-1)**2))*1
591 | c bilpd(i) = -bilpd(i)/((pi*r(i)**2)-(pi*r(i-1)**2))*1
592 | c bilv(i) = -bilv(i)/((pi*r(i)**2)-(pi*r(i-1)**2))*1
593 |
594 | c var_bil(i)=abs(var_bil(i))/((pi*r(i)**2)-
595 | c & (pi*r(i-1)**2))*1
596 |
597 | c var_bilpg(i)=abs(var_bilpg(i))/((pi*r(i)**2)-
598 | c & (pi*r(i-1)**2))*1
599 | c var_bilpd(i)=abs(var_bilpd(i))/((pi*r(i)**2)-
600 | c & (pi*r(i-1)**2))*1
601 | c var_bilv(i)=abs(var_bilv(i))/((pi*r(i)**2)-
602 | c & (pi*r(i-1)**2))*1
603 |
604 |
605 |
606 | c !bil(i) = -bil(i) * 1.d-100
607 | c !var_bil(i) = var_bil(i) * 1.d-100
608 |
609 | !if (i.eq.(n+1))then
610 | ! print *,'i ',i
611 | ! print *,bilpg(i)
612 | ! print *,bilpd(i)
613 | !!!print *,'psi paroi',psi(54,0), psi(54,111)
614 | ! print *, bilv(i)
615 | ! print *,bil(i)*1.d-100
616 | ! read *
617 | !endif
618 | c else
619 |
620 |
621 | c bil(i) = bil(i) /(2*pi*r(i) *1)
622 | c bilpg(i) = -bilpg(i)/(2*pi*r(i) *1)
623 | c bilpd(i) = -bilpd(i)/(2*pi*r(i) *1)
624 | c bilv(i) = -bilv(i)/(2*pi*r(i) *1)
625 |
626 | c var_bil(i)=abs(var_bil(i))/(2*pi*r(i) *1)
627 | c var_bilpg(i)=abs(var_bilpg(i))/(2*pi*r(i) *1)
628 | c var_bilpd(i)=abs(var_bilpd(i))/(2*pi*r(i) *1)
629 | c var_bilv(i)=abs(var_bilv(i))/(2*pi*r(i) *1)
630 |
631 | c !bil(i) = -bil(i) * 1.d-100
632 | c !var_bil(i) = var_bil(i) * 1.d-100
633 |
634 | c print *,'i ',i
635 | c print *,bilpg(i)
636 | c print *,bilpd(i)
637 | c !!!!print *,'psi paroi',psi(54,0), psi(54,111)
638 | c print *, bilv(i)
639 | c print *,bil(i)
640 | c read *
641 | c
642 | c
643 | c endif
644 | c
645 | c enddo
646 |
647 | c**************************************************************************
648 |
649 | !______________________________________________________
650 | open(10, file='psibil.oup', status='unknown')
651 | write(*,*) 'Suppression des parois pour les bilans'
652 | c elles sont inscrite a l'ecran seulement
653 | c i=0
654 | c write(10,*) i, real(r(i)), real(bil(i))
655 | do i = 1,n
656 | write(10,'(i4,10(e22.13))') i,
657 | & real((r(i)+r(i-1)-2*rsup)/2),
658 | & real(bil(i) /(4 *pi* 100 * 25*1.d+2 ))
659 | & ,real (var_bil(i) /(4 *pi*100 * 25*1.d+2 ))
660 | & , real(bilpg(i)),real(bilpd(i)),real(bilv(i))
661 | & , real(var_bilpg(i)),real(var_bilpd(i)),real(var_bilv(i))
662 | !!!!!!!!!!!!! &, real (var_bil(i)* 10)
663 | !!!!!!!!!!!!!!!!! ! derniere colonne bilan en erg/s.cm3
664 | enddo
665 | c i=n+1
666 | c write(10,*) i, real(r(i)), real(bil(i))
667 | close (10)
668 |
669 |
670 | !______________________________________________________
671 | open(10, file='psibil_parois.oup', status='unknown')
672 | i = 0
673 | write(10,'(i4,10(e22.13))') i,
674 | & real((r(i)-rsup)),
675 | & real(bil(i))
676 | & ,real (var_bil(i))
677 | & , real(bilpg(i)),real(bilpd(i)),real(bilv(i))
678 | & , real(var_bilpg(i)),real(var_bilpd(i)),real(var_bilv(i))
679 | !!!!!!!!!!!!! &, real (var_bil(i)* 10)
680 | !!!!!!!!!!!!!!!!! ! derniere colonne bilan en erg/s.cm3
681 | i = n+1
682 | write(10,'(i4,10(e22.13))') i,
683 | & real((r(i)-rsup)),
684 | & real(bil(i))
685 | & ,real (var_bil(i))
686 | & , real(bilpg(i)),real(bilpd(i)),real(bilv(i))
687 | & , real(var_bilpg(i)),real(var_bilpd(i)),real(var_bilv(i))
688 | close (10)
689 |
690 |
691 | !______________________________________________________
692 |
693 | open(10, file='rayons.oup', status='unknown')
694 | write(10,*) real (r(0)-rsup)
695 | do i = 1,n
696 | write(10,*) real((r(i)+r(i-1)-2*rsup)/2)
697 | enddo
698 | write(10,*) real (r(n+1)-rsup)
699 | c i=n+1
700 | c write(10,*) i, real(r(i)), real(bil(i))
701 | close (10)
702 |
703 | !________________________________________________________
704 | ! volumes seuls
705 | ! open(10, file='psibilfin.oup', status='unknown')
706 | ! do i = 1,n
707 | ! write(10,*) i, real((r(i)+r(i-1)-2*rsup)/2), psifin(i)
708 | ! enddo
709 | ! close (10)
710 |
711 | end
712 |
713 |