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