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