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