1 | Subroutine trajet(sens, boolener, profssmaille,
2 | & boolssmaille,boolinterieur)
3 | IMPLICIT NONE
4 | c
5 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
6 | c suivi de rayon pour "cecile.f"
7 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
8 | c
9 | c.......................................................................
10 | c
11 | INCLUDE 'cecile.inc'
12 |
13 | include 'propradia.inc'
14 | include 'radiatif.inc'
15 |
16 | include 'entre.inc'
17 |
18 | integer hete
19 | logical boolener, profssmaille, debug,boolssmaille,boolinterieur
20 |
21 | c profssmaille est faux ==> pas de profil ss maille
22 | c.......................................................................
23 | c
24 | c debut du rayon = maille in
25 | c fin du rayon = condition d'arret = paroi noire
26 | c depot d'energie sur les iin traverses
27 | c
28 | c sens : si sens=1 on se propage vers les rayons croissants
29 | c si sens=2 on commence par se propager vers les rayons decroissants
30 | c w : poids pour les psi(in,iin) de iin=0 a n+1
31 | c b : poids energetiques pour les psi(in,iin)
32 | c
33 | c b_dt : poids pour les dpsi6dt(in,iin,itmp1)
34 | c w_dk : poids pour les dpsidk(in,iin,itmp1) de iin=0 a n+1 et pour tous
35 | c les itmp1 traverses
36 | c
37 | c iin_entree : indice de la premiere maille a traverser lorsque l'on
38 | c va vers les rayons croissants ... si "sens=1" alors
39 | c in_entree=in+1 ... si "sens=2" alors iin_entree=iin_min+1
40 | c
41 | c........................................................................
42 | c MISE EN GARDE :
43 | c si desprofils sous maille sont defini par un fonction f(x)
44 | c cela devient dangereux si vous utiliser rsup
45 | c il faut en tenir compte !!!!!!! (retrancher rsup)
46 | c
47 | c Le calcul des sensibilites en parcourant itmp1, peut etre reduit
48 | c en termes de boucles...car on n'est pas sensible aux mailles non
49 | c traversees: c'est un gain important en terme de temps de calcul
50 | c........................................................................
51 |
52 | INTEGER sens
53 | double precision rayon, profblack,proftemp_index.f" HREF="../src/proftemp_index.f.html">proftemp_index, proftemp
54 | double precision temperat
55 | double precision bi, bj, cg
56 | DOUBLE PRECISION w
57 | DOUBLE PRECISION w_dk (n_mx)
58 | DOUBLE PRECISION bi_dt (n_mx)
59 | DOUBLE PRECISION bj_dt (n_mx)
60 | double precision epocg(ngaz_mx)
61 | double precision phicg(ngaz_mx), cor(ngaz_mx), correct, cg_mix
62 | double precision distance, wcg
63 | INTEGER iin_entree,a,b,mmcte, iz
64 | c
65 | c.......................................................................
66 | c
67 | DOUBLE PRECISION tmp1,tmp2,tmp3,tmp4,tmp5,tmp6
68 | INTEGER itmp1,itmp2,itmp3,itmp4,itmp5
69 | logical choixhet, choixplus
70 | c
71 | c......................pour version CG..................................
72 | c
73 | ! variables interieures a la future version cd curtis godson
74 | integer test1, test2
75 | double precision tmp7, tmp8, tmp9
76 | DOUBLE PRECISION tdphi_j,tdphi_i,tphi,tphi_i,tphi_j,tk_i,tk_j
77 | DOUBLE PRECISION te,tb,ta,tdb_i,tdb_j,tda_i,tda_j,td2phi
78 | DOUBLE PRECISION tam12,tap12,tam32
79 | double precision tGa_i(ngaz_mx),tGa_j(ngaz_mx),tdGa_ij(ngaz_mx)
80 | DOUBLE PRECISION tGb_i(ngaz_mx),tGb_j(ngaz_mx),tdGb_ij(ngaz_mx)
81 | double precision tdGa_ji,toto,tGa, tGb, tGc,tGd
82 | double precision correction !rajouter cor(ngaz_mx) aussi
83 | c
84 | c.......................................................................
85 | c
86 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
87 | c
88 | c.......................................................................
89 | c mise a zero des poids etc...
90 | c.......................................................................
91 | choixhet = (hetero .le. 2) !sep ou cg sont a activer
92 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!! choixplus = (hetero .le. 0) !cg2
93 | debug=.False.
94 | ngaz = 3
95 | do igaz=1,ngaz
96 | cor(igaz)=1.d+0
97 | enddo
98 | correct = 1.d+0
99 | distance = 0.d+0
100 | mmcte = 110
101 | c rien a faire
102 | do iin=1,n
103 | w_dk (iin)= 0.d+0
104 | bi_dt (iin)=0.d+0
105 | bj_dt (iin)=0.d+0
106 | enddo
107 | iin = in
108 | c.......................................................................
109 | c pour la maille d'emission in: calcul des proprietes
110 | c.......................................................................
111 | c*****profil sous maille: sigma_1 deja tire
112 | !sigma_1 = 0.d+0
113 | if ((in.gt.0).and.(in.lt.(n+1))) then
114 | !volume (origine la sortie) + modif amaury 20 janv 2000
115 | if (sens.eq.2) then
116 | !sigma1 decroissant contre paroi int de la couronne
117 | rayon = dsqrt(delta**2 + ((long(in-1)/f)+(sigma_1/f))**2)
118 | else
119 | !sigma1 croissant contre paroi ext de la couronne
120 | rayon = dsqrt(delta**2 + ((long(in)/f)-(sigma_1/f))**2)
121 | endif
122 |
123 | else
124 | ! paroi
125 | sigma_1 = 0.d+0
126 | rayon=r(in)
127 | endif
128 | !------------CURTIS GODSON
129 | if (choixhet) then
130 | distance = distance + sigma_1
131 | do igaz=1,ngaz
132 | kbarphi_s(igaz) = kbarphi_s(igaz) +
133 | & kgbar(igaz,in)* phig(igaz,in) *sigma_1
134 | !if ((in.eq.25)) then
135 | ! print *,kbar_s(igaz),kgbar(igaz,in),sigma_1
136 | !endif
137 | kbar_s(igaz) = kbar_s(igaz) + kgbar(igaz,in)*sigma_1
138 | !A cette etape on a pas besoin de phicg
139 | !et le denominateur est souvent nul
140 | !phicg(igaz) = (kbarphi_s(igaz) / kbar_s(igaz))
141 | epocg(igaz) = kbar_s(igaz) !ATTENTION / distance
142 | ! if (kbar_s(igaz).lt.1.d-10) phicg(igaz)=phig(igaz,in)
143 | ! if (distance .eq. 0.d+0) epocg(igaz) = kgbar(igaz,in)
144 | !print *,igaz,kbar_s(igaz), phig(igaz,in), phicg(igaz)
145 | ! cor(igaz)=cg(in,iin,kgskgbar(igaz),phig(igaz,in),phicg(igaz)
146 | ! & ,kbar_s(igaz),kgbar(igaz,in),sigma_1 !)
147 | ! & ,kgbar(igaz,iin),0.d+0,phig(igaz,iin))
148 | !!print *,'EMISSION: ', igaz, phicg(igaz),kbar_s(igaz), in
149 | !!print *,'', cor(igaz)
150 | enddo
151 | ! if (choixplus) correct = cg_mix(epocg,phicg)
152 | endif
153 | !-------------------------
154 |
155 | c*****fonction: partie geom--> rien
156 | c*****fonction: partie energetique
157 | temperat= proftemp_index(in,rayon,boolinterieur)
158 | if (boolssmaille) temperat = proftemp(rayon, rsup)
159 | !print *, in, rayon-rsup, temperat
160 | !read *
161 | !call black(eta(ibande)*1.d+2, temp(in),emi)
162 | call black(eta(ibande)*1.d+2, temperat,emi)
163 | bi = (emi)/pi
164 |
165 | c pour papier milieux epais: mcep
166 | c bi = lblack(rayon)
167 | if (profssmaille) then
168 | call error(boolssmaille)
169 | bi = profblack(r(0), r(n+1), rayon, rsup)
170 | endif
171 | if (.not. boolener) then
172 | bi = 1.d+0
173 | bj= 0.d+0
174 | endif
175 | c*****pas d'attenuation ???
176 |
177 |
178 | c
179 | c.......................................................................
180 | c
181 | IF (sens.EQ.2) THEN
182 | c
183 | c.......................................................................
184 | c si sens=2
185 | c on fait tous les echanges avec les mailles plus petites
186 | c (jusqu'a iin_min)
187 | c si iin_min=0 c'est qu'on a rencontre la paroi interne, on fait l'echange
188 | c avec cette paroi et c'est fini
189 | c si iin_min>0 on repart ensuite vers les rayons croissants jusqu'a la
190 | c paroi externe (pour ca on retrouve l'algorithme de "sens=1")
191 | c.......................................................................
192 | c
193 | DO iin=in-1,iin_min+1,-1
194 |
195 | c****** profil sous maille iin: reception
196 | c -----{tirage de sigma_2 dans la maille iin}
197 | c
198 | call rand_uniforme(tmp2)
199 | sigma_2 =
200 | & - dlog(1.D+0
201 | & -tmp2*(1.D+0-dexp(-k(iin)*(long(iin)-long(iin-1)))))
202 | & /k(iin)
203 | if (k(iin).eq. 0.d+0) sigma_2 = 0.d+0
204 | !volume (origine a l'entree)
205 | cc rayon = dsqrt(delta**2 + ((long(iin-1)/f)+(sigma_2/f))**2)
206 | rayon = dsqrt(delta**2 + ((long(iin)/f)-(sigma_2/f))**2)
207 | ! print *, 'iin,rayon: ',iin , rayon
208 | ! read *
209 | !------------CURTIS GODSON
210 | if (choixhet) then
211 | tmp6 = distance + sigma_2
212 | do igaz=1,ngaz
213 | tmp4 = kbarphi_s(igaz) +
214 | & kgbar(igaz,iin)* phig(igaz,iin) *sigma_2
215 | !if ((in.eq.25).and.(iin.eq.123)) then
216 | ! print *,'**: ',kbar_s(igaz) + kgbar(igaz,iin)*sigma_2,
217 | c !& kbar_s(igaz),kgbar(igaz,in),sigma_2
218 | ! endif
219 | tmp5 = kbar_s(igaz) + kgbar(igaz,iin)*sigma_2
220 |
221 | phicg(igaz) = tmp4/ tmp5
222 | epocg(igaz) = tmp5 ! ATTENTION/ tmp6
223 | !if (tmp5.lt.1.d-10) phicg(igaz)=phig(igaz,in)
224 | !if (tmp6 .eq. 0.d+0) epocg(igaz) = kgbar(igaz,in)
225 | !print *,igaz,tmp5, phig(igaz,in), phicg(igaz)
226 | ! cor(igaz)=cg(in,iin,kgskgbar(igaz),phig(igaz,in),phicg(igaz),
227 | ! & tmp5,
228 | ! & kgbar(igaz,in),sigma_1 !)
229 | ! & ,kgbar(igaz,iin),sigma_2,phig(igaz,iin))
230 | !cor(igaz)=1
231 | !!print *, 'TRAV DEC',igaz, phicg(igaz),tmp5,in,iin
232 | !!print *,'', cor(igaz)
233 | enddo
234 |
235 | include 'curtis.inc'
236 | correct = wcg
237 | endif
238 | !-------------------------
239 |
240 | c****** calcul de la fonction psi(integrale) et ses derivees
241 | c -----{increment des poids optico-geom}
242 | tmp1=(long(iin)-long(iin-1))
243 | ! w= cor(1)*cor(2)*cor(3)*correct*
244 | w= correct*
245 | & poids*(1.D+0-dexp(-k(iin)*tmp1))
246 |
247 | c -----{increment des poids energetique: pas de cumul}
248 | temperat= proftemp_index(iin,rayon,boolinterieur)
249 | if (boolssmaille) temperat = proftemp(rayon, rsup)
250 | !print *, iin, rayon-rsup, temperat
251 | !read *
252 | call black (eta(ibande)*1.d+2, temperat,emj)
253 | bj = (emj)/pi
254 |
255 | c pour papier milieux epais: mcep
256 | c bi = lblack(rayon)
257 | if (profssmaille) then
258 | call error(boolssmaille)
259 | bj = profblack(r(0), r(n+1), rayon, rsup)
260 | endif
261 |
262 | c -----{calcul de psi...}
263 | if (debug.and.(in.eq.mmcte)) then
264 | print *, 'l1'
265 | print *, w, dexp(-k(iin)*tmp1), k(iin),poids
266 | endif
267 | if (.not. boolener) then
268 | bi = 1.d+0
269 | bj= 0.d+0
270 | endif
271 | call trajet_fonc(w,bi,bj)
272 |
273 | DO itmp1=in,iin,-1
274 | cccccccccccccccccccccccccccccccccccccccccccccc DO itmp1=1,n
275 | ! increment des poids optico-geom
276 | ! w_dk(itmp1)=cor(1)*cor(2)*cor(3)*correct*
277 | w_dk(itmp1)= correct*
278 | & poids*dpdksp(itmp1)*(1.D+0-dexp(-k(iin)*tmp1))
279 | ccc debug
280 | if (dexp(-k(iin)*tmp1).eq.1.d+0) then
281 | w_dk(itmp1) = 0.d+0
282 | endif
283 |
284 | c if (itmp1.eq.in)then
285 | c print *, sens, in, iin, w_dk(itmp1)
286 | c endif
287 | !! cas particulier de iin
288 | if (itmp1.eq.iin) then
289 | ! w_dk(iin)=cor(1)*cor(2)*cor(3)*correct*
290 | w_dk(iin)=correct*
291 | & poids*tmp1*dexp(-k(iin)*tmp1)
292 | endif
293 |
294 | ! increment poids energetique
295 | !!seule in et iin sont non nuls mais dans la boucle
296 | if (itmp1.eq.in) then
297 | !!temperat= proftemp_index(in,rayon)
298 | !!if (boolssmaille) temperat = proftemp(rayon, rsup)
299 | call black_dt (eta(ibande)*1.d+2, temperat,demi)
300 | ! print *, in, iin, ' ', eta(ibande)*1.d+2, temperat,demi
301 | ! read *
302 | bi_dt(itmp1)= (demi)/pi
303 | bj_dt(itmp1)= 0.d+0
304 | elseif (itmp1.eq.iin) then
305 | !!temperat= proftemp_index(iin,rayon)
306 | !!if (boolssmaille) temperat = proftemp(rayon, rsup)
307 | call black_dt (eta(ibande)*1.d+2, temperat,demj)
308 | bj_dt(itmp1)= (demj)/pi
309 | bi_dt(itmp1)=0.d+0
310 | else
311 | bj_dt(itmp1)= 0.d+0
312 | bi_dt(itmp1)= 0.d+0
313 | endif
314 |
315 |
316 | ! calcul de dpsi...
317 | call trajet_sens(itmp1,w,bi,bj,
318 | & w_dk(itmp1),bi_dt(itmp1),bj_dt(itmp1))
319 |
320 | ENDDO
321 |
322 |
323 | c
324 | c****** attenuation dans la maille iin}
325 | c
326 | tmp1=(long(iin)-long(iin-1))
327 | poids=poids*dexp(-k(iin)*tmp1)
328 | dpdksp(iin)=dpdksp(iin)-tmp1
329 | !------------CURTIS GODSON
330 | if (choixhet) then
331 | distance = distance + tmp1
332 | do igaz=1,ngaz
333 | kbarphi_s(igaz) = kbarphi_s(igaz) +
334 | & kgbar(igaz,iin)* phig(igaz,iin) *tmp1
335 | kbar_s(igaz) = kbar_s(igaz) + kgbar(igaz,iin)*tmp1
336 | phicg(igaz) = (kbarphi_s(igaz) / kbar_s(igaz))
337 | epocg(igaz) = kbar_s(igaz) ! ATTENTION / distance
338 | ! if (kbar_s(igaz).lt.1.d-10) phicg(igaz)=phig(igaz,in)
339 | ! if (distance .eq. 0.d+0) epocg(igaz) = kgbar(igaz,in)
340 | !print *,igaz,kbar_s(igaz), phig(igaz,in), phicg(igaz)
341 | ! cor(igaz)=cg(in,iin,kgskgbar(igaz),phig(igaz,in),phicg(igaz)
342 | ! & ,kbar_s(igaz),kgbar(igaz,in),sigma_1 !)
343 | ! & ,kgbar(igaz,iin),sigma_2,phig(igaz,iin))
344 | !!print *, 'TRAV DEC tot',igaz, phicg(igaz),kbar_s(igaz),in,iin
345 | !!print *,'', cor(igaz)
346 | !cor(igaz)=1
347 | !!!!!!!!!!!print *, igaz,phicg(igaz),kbarphi_s(igaz), kbar_s(igaz)
348 | !if(in.eq.3) print*, 'kbphi, kb', kbarphi_s(3), kbar_s(3)
349 | !if(in.eq.3) print*, '..............'
350 | !if(in.eq.3) print*, 'in=', iin
351 | !if(in.eq.3) print*, '..............'
352 | enddo
353 |
354 | ! if (choixplus) correct = cg_mix(epocg, phicg)
355 | endif
356 | !-------------------------
357 | c
358 | ENDDO
359 | c
360 | c ---{traitement particulier pour la maille iin_min= iin courant}
361 | c
362 | IF (iin_min .EQ.0) THEN
363 |
364 | c -----{paroi interne}
365 | c****** profil sous maille iin: reception
366 | ! pas d'épaisseur
367 | iin = iin_min
368 | rayon = r(iin)
369 | rayon = r(0)
370 | !------------CURTIS GODSON
371 | if (choixhet) then
372 | !rien a faire
373 | do igaz=1,ngaz
374 | tmp4 = kbarphi_s(igaz)
375 | tmp5 = kbar_s(igaz)
376 | phicg(igaz) = tmp4 /tmp5
377 | epocg(igaz) = tmp5
378 | !kbarphi_s(igaz) = kbarphi_s(igaz) +
379 | ! kgbar(igaz,in)* phig(igaz,in) *sigma_1
380 | ! kbar_s(igaz) = kbar_s(igaz) + kgbar(igaz,in)*sigma_1
381 | ! phicg(igaz) = (kbarphi_s(igaz) / kbar_s(igaz))
382 | ! cor(igaz)=cg(in,iin,kgskgbar(igaz),phig(igaz,in),phicg(igaz),
383 | ! & tmp5,
384 | ! & kgbar(igaz,in),sigma_1 !)
385 | ! & ,kgbar(igaz,iin),0.d+0,phig(igaz,iin))
386 | enddo
387 | !distance = distance +0
388 | include 'curtis.inc'
389 | correct = wcg
390 | endif
391 | !-------------------------
392 | c****** calcul de la fonction psi(integrale) et ses derivees
393 | c------{increment poids optico-geome}
394 | !w=poids*cor(1)*cor(2)*cor(3)*correct
395 | w=poids*correct
396 | c -----{increment des poids energetique: pas de cumul}
397 | temperat= proftemp_index(0,rayon,boolinterieur)
398 | if (boolssmaille) temperat = proftemp(rayon, rsup)
399 | !print *, 0, rayon-rsup, temperat
400 | !read *
401 | call black (eta(ibande)*1.d+2, temperat,emj)
402 | bj = (emj)/pi
403 |
404 | c pour papier milieux epais: mcep
405 | c bi = lblack(rayon)
406 | if (profssmaille) then
407 | call error(boolssmaille)
408 | bj = profblack(r(0), r(n+1), rayon, rsup)
409 | endif
410 |
411 | c -----{calcul de psi...}
412 | iin = 0
413 | if (debug.and.(in.eq.mmcte)) then
414 | print *,'l2'
415 | print *, w
416 | endif
417 | if (.not. boolener) then
418 | bi = 1.d+0
419 | bj= 0.d+0
420 | endif
421 | call trajet_fonc(w,bi,bj)
422 | DO itmp1=in,iin_min,-1
423 | ccccccccccccccccccccccccccccccccccccccccccccc DO itmp1=1,n
424 | ! increment poids optico-geom
425 | ! w_dk(itmp1)=cor(1)*cor(2)*cor(3)*correct*
426 | w_dk(itmp1)=correct*
427 | & poids*dpdksp(itmp1)
428 |
429 | ! increment poids energetique
430 | !!seule in et iin sont non nuls mais dans la boucle
431 | if (itmp1.eq.in) then
432 | !!temperat= proftemp_index(in,rayon)
433 | !!if (boolssmaille) temperat = proftemp(rayon, rsup)
434 | call black_dt (eta(ibande)*1.d+2, temperat,demi)
435 | ! print *, in, iin, ' ', eta(ibande)*1.d+2, temperat,demi
436 | ! read *
437 | bi_dt(itmp1)= (demi)/pi
438 | bj_dt(itmp1)= 0.d+0
439 | elseif (itmp1.eq.iin) then
440 | !!temperat= proftemp_index(iin,rayon)
441 | !!if (boolssmaille) temperat = proftemp(rayon, rsup)
442 | call black_dt (eta(ibande)*1.d+2, temperat,demj)
443 | bj_dt(itmp1)= (demj)/pi
444 | bi_dt(itmp1)=0.d+0
445 | else
446 | bj_dt(itmp1)= 0.d+0
447 | bi_dt(itmp1)= 0.d+0
448 | endif
449 |
450 |
451 | ! calcul de dpsi...
452 | call trajet_sens(itmp1,w,bi,bj,
453 | & w_dk(itmp1),bi_dt(itmp1),bj_dt(itmp1))
454 |
455 | ENDDO
456 | c****** attenuation dans une paroi: non
457 | GOTO 10002
458 | c
459 | ELSE
460 | iin = iin_min
461 | c****** profil sous maille iin: reception
462 | c -----{tirage de sigma_2 dans la maille iin_min}
463 | c
464 | call rand_uniforme(tmp2)
465 | sigma_2=
466 | & -dlog(1.D+0
467 | & -tmp2*(1.D+0-dexp(-k(iin_min)*(2.D+0*long(iin_min)))))
468 | & /k(iin_min)
469 | if (k(iin).eq. 0.d+0) sigma_2 = 0.d+0
470 | !sigma_2 = 0.d+0
471 | !volume (origine a l'entree)
472 | cc rayon = dsqrt(delta**2 + ((long(iin-1)/f)+(sigma_2/f))**2)
473 | rayon = dsqrt(delta**2 + ((long(iin_min)/f)
474 | & -(sigma_2/f))**2)
475 | !------------CURTIS GODSON
476 | if (choixhet) then
477 | tmp6 = distance + sigma_2
478 | do igaz=1,ngaz
479 | tmp4 = kbarphi_s(igaz) +
480 | & kgbar(igaz,iin)* phig(igaz,iin) *sigma_2
481 | tmp5 = kbar_s(igaz) + kgbar(igaz,iin)*sigma_2
482 |
483 | phicg(igaz) = tmp4 / tmp5
484 | epocg(igaz) = tmp5 ! ATTENTION / tmp6
485 | !if (tmp5.lt.1.d-10) phicg(igaz)=phig(igaz,in)
486 | !if (tmp6 .eq. 0.d+0) epocg(igaz) = kgbar(igaz,in)
487 | ! cor(igaz)=cg(in,iin,kgskgbar(igaz),phig(igaz,in),phicg(igaz),
488 | ! & tmp5,
489 | ! & kgbar(igaz,in),sigma_1 !)
490 | ! & ,kgbar(igaz,iin),sigma_2,phig(igaz,iin))
491 | !!print *, 'MIN VOL ',igaz, phicg(igaz),tmp5,in,iin
492 | !!print *,'', cor(igaz)
493 | !cor(igaz)=1
494 | enddo
495 | include 'curtis.inc'
496 | correct = wcg
497 |
498 | endif
499 | !-------------------------
500 |
501 | c****** calcul de la fonction psi(integrale) et ses derivees
502 | c -----{increment du poids optico-geom}
503 | c
504 | tmp1=(2.D+0*long(iin_min))
505 | ! w=cor(1)*cor(2)*cor(3)*correct*
506 | w=correct*
507 | & poids*(1.D+0-dexp(-k(iin_min)*tmp1))
508 |
509 | c -----{increment poids energetique}
510 | temperat= proftemp_index(iin_min,rayon,boolinterieur)
511 | if (boolssmaille) temperat = proftemp(rayon, rsup)
512 | !print *, iin_min, rayon-rsup, temperat
513 | !read *
514 | call black(eta(ibande)*1.d+2, temperat,emj)
515 | bj = (emj)/pi
516 |
517 | c pour papier milieux epais: mcep
518 | c bi = lblack(rayon)
519 | if (profssmaille) then
520 | call error(boolssmaille)
521 | bj = profblack(r(0), r(n+1), rayon, rsup)
522 | endif
523 |
524 | c -----{calcul de psi...}
525 | iin = iin_min
526 | if (debug.and.(in.eq.mmcte)) then
527 | print *, 'l3'
528 | print *, w
529 | endif
530 | if (.not. boolener) then
531 | bi = 1.d+0
532 | bj= 0.d+0
533 | endif
534 | call trajet_fonc(w,bi,bj)
535 |
536 | DO itmp1=in,iin_min,-1
537 | ccccccccccccccccccccccccccccccccccccccccccccccc DO itmp1=1,n
538 | ! increment poids optico-geom
539 | ! w_dk(itmp1)=cor(1)*cor(2)*cor(3)*correct*
540 | w_dk(itmp1)=correct*
541 | & poids*dpdksp(itmp1)*(1.D+0-dexp(-k(iin_min)*tmp1))
542 | ccc debug
543 | if (dexp(-k(iin_min)*tmp1).eq. 1.d+0) then
544 | w_dk(itmp1) = 0.d+0
545 | endif
546 | ccc if (itmp1.eq.in)then
547 | ccc print *, sens, in, iin, w_dk(itmp1)
548 | ccc endif
549 |
550 |
551 | if (itmp1.eq.iin_min) then
552 | ! w_dk(iin_min)=cor(1)*cor(2)*cor(3)*correct*
553 | w_dk(iin_min)=correct*
554 | & poids*tmp1*dexp(-k(iin_min)*tmp1)
555 | endif
556 |
557 | !increment poids energetique
558 | !!seule in et iin sont non nuls mais dans la boucle
559 | if (itmp1.eq.in) then
560 | !!temperat= proftemp_index(in,rayon)
561 | !!if (boolssmaille) temperat = proftemp(rayon, rsup)
562 | call black_dt (eta(ibande)*1.d+2, temperat,demi)
563 | bi_dt(itmp1)= (demi)/pi
564 | bj_dt(itmp1)= 0.d+0
565 | elseif (itmp1.eq.iin) then
566 | !!temperat= proftemp_index(iin,rayon)
567 | !!if (boolssmaille) temperat = proftemp(rayon, rsup)
568 | call black_dt (eta(ibande)*1.d+2, temperat,demj)
569 | bj_dt(itmp1)= (demj)/pi
570 | bi_dt(itmp1)=0.d+0
571 | else
572 | bj_dt(itmp1)= 0.d+0
573 | bi_dt(itmp1)= 0.d+0
574 | endif
575 |
576 |
577 | !calcul de dpsi...
578 | call trajet_sens(itmp1,w,bi,bj,
579 | & w_dk(itmp1),bi_dt(itmp1),bj_dt(itmp1))
580 | ENDDO
581 | c
582 | c****** attenuation dans la maille iin_min}
583 | c
584 | tmp1=(2.D+0*long(iin_min))
585 | poids=poids*dexp(-k(iin_min)*tmp1)
586 | dpdksp(iin_min)=dpdksp(iin_min)-tmp1
587 | !------------CURTIS GODSON
588 | if (choixhet) then
589 | distance = distance + tmp1
590 | do igaz=1,ngaz
591 | kbarphi_s(igaz) = kbarphi_s(igaz) +
592 | & kgbar(igaz,iin)* phig(igaz,iin) * tmp1
593 | kbar_s(igaz) = kbar_s(igaz) + kgbar(igaz,iin)* tmp1
594 |
595 | phicg(igaz) = (kbarphi_s(igaz) / kbar_s(igaz))
596 | epocg(igaz) = kbar_s(igaz) ! ATTENTION / distance
597 | ! if (kbar_s(igaz).lt. 1.d-10) phicg(igaz)=phig(igaz,in)
598 | ! if (distance .eq. 0.d+0) epocg(igaz) = kgbar(igaz,in)
599 | ! cor(igaz)=cg(in,iin,kgskgbar(igaz),phig(igaz,in),phicg(igaz),
600 | ! & kbar_s(igaz),kgbar(igaz,in),sigma_1 !)
601 | ! & ,kgbar(igaz,iin),sigma_2,phig(igaz,iin))
602 | !!print *, 'MIN VOL tot',igaz, phicg(igaz),kbar_s(igaz),in,iin
603 | !!print *,'', cor(igaz)
604 | !cor(igaz)=1
605 | enddo
606 |
607 | ! if (choixplus) correct = cg_mix(epocg, phicg)
608 | endif
609 | !-------------------------
610 | c
611 | iin_entree=iin_min+1
612 | ENDIF
613 | c
614 |
615 |
616 | ELSE
617 | ! sens=1 des le depart
618 | iin_entree=in+1
619 | ENDIF
620 |
621 | c.......................................................................
622 | c si sens=1
623 | c.......................................................................
624 | DO iin=iin_entree,n
625 | ! la condition ci-dessous elimine le calcul de psi(a,a) nul
626 | ! par definition.
627 | ! dans le cas general ca arrive tres rarement donc pas tres important
628 | IF (iin.NE.in) THEN
629 |
630 | c****** profil sous maille iin: reception
631 | c -------{tirage de sigma_2 dans la maille iin}
632 | c
633 | call rand_uniforme(tmp2)
634 | sigma_2=
635 | & -dlog(1.D+0
636 | & -tmp2*(1.D+0-dexp(-k(iin)*(long(iin)-long(iin-1)))))
637 | & /k(iin)
638 | if (k(iin).eq. 0.d+0) sigma_2 = 0.d+0
639 | !sigma_2 = 0.d+0
640 | !volume (origine a l'entree)
641 | rayon = dsqrt(delta**2 + ((long(iin-1)/f)+(sigma_2/f))**2)
642 | !------------CURTIS GODSON
643 | if (choixhet) then
644 | tmp6 = distance + sigma_2
645 | do igaz=1,ngaz
646 | tmp4 = kbarphi_s(igaz) +
647 | & kgbar(igaz,iin)* phig(igaz,iin) *sigma_2
648 | c !if ((in.eq.25).and.(iin.eq.123)) then
649 | c ! print *,'**: ',kbar_s(igaz) + kgbar(igaz,iin)*sigma_2,
650 | c !& kbar_s(igaz),kgbar(igaz,in),sigma_2
651 | c ! endif
652 | tmp5 = kbar_s(igaz) + kgbar(igaz,iin)*sigma_2
653 |
654 | phicg(igaz) = tmp4 / tmp5
655 | epocg(igaz) = tmp5 ! ATTENTION / tmp6
656 | !if (tmp5.lt.1.d-10) phicg(igaz)=phig(igaz,in)
657 | !if (tmp6 .eq. 0.d+0) epocg(igaz) = kgbar(igaz,in)
658 | ! cor(igaz)=cg(in,iin,kgskgbar(igaz),phig(igaz,in),phicg(igaz),
659 | ! & tmp5,
660 | ! & kgbar(igaz,in),sigma_1 !)
661 | ! & ,kgbar(igaz,iin),sigma_2,phig(igaz,iin))
662 | !!print *, 'CROIS',igaz, phicg(igaz),tmp5,in,iin
663 | !! print *,'', cor(igaz)
664 | !cor(igaz)=1
665 | enddo
666 | include 'curtis.inc'
667 | correct = wcg
668 | endif
669 | !-------------------------
670 |
671 | c****** calcul de la fonction psi(integrale) et ses derivees
672 | c -------{increment du poids optico-geom}
673 | c
674 | tmp1=(long(iin)-long(iin-1))
675 | ! w=cor(1)*cor(2)*cor(3)*correct*
676 | w=correct*
677 | & poids*(1.D+0-dexp(-k(iin)*tmp1))
678 |
679 | c -------{increment poids energetique}
680 | temperat= proftemp_index(iin,rayon,boolinterieur)
681 | if (boolssmaille) temperat = proftemp(rayon, rsup)
682 | !print *, iin, rayon-rsup, temperat
683 | !read *
684 | call black(eta(ibande)*1.d+2, temperat,emj)
685 | bj = (emj)/pi
686 |
687 | c pour papier milieux epais: mcep
688 | c bi = lblack(rayon)
689 | if (profssmaille) then
690 | call error(boolssmaille)
691 | bj = profblack(r(0), r(n+1), rayon, rsup)
692 | endif
693 |
694 | c -------{calcul de psi...}
695 | if (debug.and.(in.eq.mmcte)) then
696 | print *, 'l4'
697 | print *, w
698 | endif
699 | if (.not. boolener) then
700 | bi = 1.d+0
701 | bj= 0.d+0
702 | endif
703 | call trajet_fonc(w,bi,bj)
704 | if (sens.eq.2)then
705 | a = iin_min
706 | b = iin
707 | else
708 | a=in
709 | b=iin
710 | endif
711 | DO itmp1=a,b
712 | ccccccccccccccccccccccccccccccccccccccccccccccc DO itmp1=1,n
713 | !increment poids optico-geom
714 | ! w_dk(itmp1)=cor(1)*cor(2)*cor(3)*correct*
715 | w_dk(itmp1)=correct*
716 | & poids*dpdksp(itmp1)*(1.D+0-dexp(-k(iin)*tmp1))
717 | ccc debug
718 | if (dexp(-k(iin)*tmp1).eq.1.d+0) then
719 | w_dk(itmp1) = 0.d+0
720 | endif
721 | c if ((itmp1.eq.in).and.(in.eq.1)) then
722 | c print *, sens, in, iin, w_dk(itmp1)
723 | c print *, 'orig: ', poids,dpdksp(itmp1),
724 | c & (1.D+0-dexp(-k(iin)*tmp1))
725 | c read *
726 | c endif
727 |
728 |
729 | if (itmp1.eq.iin) then
730 | w_dk(iin)=
731 | & poids*tmp1*dexp(-k(iin)*tmp1)
732 | endif
733 |
734 | !increment poids energetique
735 | !!seule in et iin sont non nuls mais dans la boucle
736 | if (itmp1.eq.in) then
737 | !!temperat= proftemp_index(in,rayon)
738 | !!if (boolssmaille) temperat = proftemp(rayon, rsup)
739 | call black_dt (eta(ibande)*1.d+2, temperat,demi)
740 | bi_dt(itmp1)= (demi)/pi
741 | bj_dt(itmp1)= 0.d+0
742 | elseif (itmp1.eq.iin) then
743 | !!temperat= proftemp_index(iin,rayon)
744 | !!if (boolssmaille) temperat = proftemp(rayon, rsup)
745 | call black_dt (eta(ibande)*1.d+2, temperat,demj)
746 | bj_dt(itmp1)= (demj)/pi
747 | bi_dt(itmp1)=0.d+0
748 | else
749 | bj_dt(itmp1)= 0.d+0
750 | bi_dt(itmp1)= 0.d+0
751 | endif
752 |
753 |
754 | !calcul de dpsi...
755 | call trajet_sens(itmp1,w,bi,bj,
756 | & w_dk(itmp1),bi_dt(itmp1),bj_dt(itmp1))
757 | ENDDO
758 |
759 | c
760 | ENDIF
761 | c
762 | c******* attenuation dans la maille iin}
763 | c
764 | tmp1=(long(iin)-long(iin-1))
765 | poids=poids*dexp(-k(iin)*tmp1)
766 | dpdksp(iin)=dpdksp(iin)-tmp1
767 | !------------CURTIS GODSON
768 | if (choixhet) then
769 | distance = distance + tmp1
770 | do igaz=1,ngaz
771 | kbarphi_s(igaz) = kbarphi_s(igaz) +
772 | & kgbar(igaz,iin)* phig(igaz,iin) * tmp1
773 | kbar_s(igaz) = kbar_s(igaz) + kgbar(igaz,iin)* tmp1
774 | phicg(igaz) = (kbarphi_s(igaz) / kbar_s(igaz))
775 | epocg(igaz) = kbar_s(igaz) ! ATTENION / distance
776 | ! if (kbar_s(igaz).lt.1.d-10) phicg(igaz)=phig(igaz,in)
777 | ! if (distance .eq. 0.d+0) epocg(igaz) = kgbar(igaz,in)
778 | ! cor(igaz)=cg(in,iin,kgskgbar(igaz),phig(igaz,in),phicg(igaz)
779 | ! & ,kbar_s(igaz),kgbar(igaz,in),sigma_1 !)
780 | ! & ,kgbar(igaz,iin),sigma_2,phig(igaz,iin))
781 | !!print *, 'CROISS tot',igaz, phicg(igaz),kbar_s(igaz),in,iin
782 | !!print *,'', cor(igaz)
783 | !cor(igaz)=1
784 | enddo
785 |
786 | ! if (choixplus) correct = cg_mix(epocg, phicg)
787 | endif
788 | !-------------------------
789 | c
790 | ENDDO
791 | c
792 | c ---{paroi externe}
793 | c****** profil sous maille iin=n+1: reception
794 | !pas d'épaisseur
795 | iin = n+1
796 | rayon = r(iin)
797 | rayon = r(n+1)
798 | !------------CURTIS GODSON
799 | if (choixhet) then
800 | !rien a faire
801 | do igaz=1,ngaz
802 | tmp4 = kbarphi_s(igaz)
803 | tmp5 = kbar_s(igaz)
804 | phicg(igaz) = tmp4 /tmp5
805 | epocg(igaz) = tmp5
806 | !kbarphi_s(igaz) = kbarphi_s(igaz) +
807 | ! kgbar(igaz,in)* phig(igaz,in) *sigma_1
808 | ! kbar_s(igaz) = kbar_s(igaz) + kgbar(igaz,in)*sigma_1
809 | ! phicg(igaz) = (kbarphi_s(igaz) / kbar_s(igaz))
810 | ! cor(igaz)=cg(in,iin,kgskgbar(igaz),phig(igaz,in),phicg(igaz),
811 | ! & tmp5,
812 | ! & kgbar(igaz,in),sigma_1 !)
813 | ! & ,kgbar(igaz,iin),0.d+0,phig(igaz,iin))
814 | enddo
815 | include 'curtis.inc'
816 | correct = wcg
817 | endif
818 | !-------------------------
819 | c****** calcul de la fonction psi(integrale) et ses derivees
820 | IF (in.NE.n+1) THEN
821 | ! la condition ci-dessous elimine le calcul de psi(a,a) nul
822 | ! par definition.
823 | ! dans le cas general ca arrive tres rarement donc pas tres important
824 |
825 | c ------{increment poids optico-geom}
826 | ! w=poids*cor(1)*cor(2)*cor(3)*correct
827 | w=poids*correct
828 |
829 | c-------{increment poids energetique}
830 | temperat= proftemp_index(n+1,rayon,boolinterieur)
831 | if (boolssmaille) temperat = proftemp(rayon, rsup)
832 | !print *, n+1, rayon-rsup, temperat
833 | !read *
834 | call black(eta(ibande)*1.d+2, temperat,emj)
835 | bj = (emj)/pi
836 |
837 | c pour papier milieux epais: mcep
838 | c bi = lblack(rayon)
839 | if (profssmaille) then
840 | call error(boolssmaille)
841 | bj = profblack(r(0), r(n+1), rayon, rsup)
842 | endif
843 | c ------{calcul de psi...}
844 | iin = n+1
845 | if (debug.and.(in.eq.mmcte)) then
846 | print *, 'l5'
847 | print *, w
848 | endif
849 | if (.not. boolener) then
850 | bi = 1.d+0
851 | bj= 0.d+0
852 | endif
853 | call trajet_fonc(w,bi,bj)
854 | if (sens.eq.2)then
855 | a = iin_min
856 | b = n
857 | else
858 | a= in
859 | b= n
860 | endif
861 | DO itmp1=a,b
862 | ccccccccccccccccccccccccccccccccccccccccccccccc DO itmp1=1,n
863 | !increment poids optico-geom
864 | ! w_dk(itmp1)=cor(1)*cor(2)*cor(3)*correct*
865 | w_dk(itmp1)=correct*
866 | & poids*dpdksp(itmp1)
867 |
868 | !increment poids energetique
869 | !!seule in et iin sont non nuls mais dans la boucle
870 | if (itmp1.eq.in) then
871 | !!temperat= proftemp_index(in,rayon)
872 | !!if (boolssmaille) temperat = proftemp(rayon, rsup)
873 | call black_dt (eta(ibande)*1.d+2, temperat,demi)
874 | ! print *, in, iin, ' ', eta(ibande)*1.d+2, temperat,demi
875 | ! read *
876 | bi_dt(itmp1)= (demi)/pi
877 | bj_dt(itmp1)= 0.d+0
878 | elseif (itmp1.eq.iin) then
879 | !!temperat= proftemp_index(iin,rayon)
880 | !!if (boolssmaille) temperat = proftemp(rayon, rsup)
881 | call black_dt (eta(ibande)*1.d+2, temperat,demj)
882 | bj_dt(itmp1)= (demj)/pi
883 | bi_dt(itmp1)=0.d+0
884 | else
885 | bj_dt(itmp1)= 0.d+0
886 | bi_dt(itmp1)= 0.d+0
887 | endif
888 |
889 | ! calcul de dpsi...
890 | call trajet_sens(itmp1,w,bi,bj,
891 | & w_dk(itmp1),bi_dt(itmp1),bj_dt(itmp1))
892 |
893 |
894 | ENDDO
895 |
896 | ENDIF
897 | c****** attenuation dans une paroi: non
898 |
899 |
900 | 10002 CONTINUE
901 |
902 |
903 | c
904 | c.......................................................................
905 | c
906 | RETURN
907 | END
908 |
909 | c
910 | c----------------------------------------------------------------------
911 | c
912 | ! amaury nov 2002
913 | subroutine error(boolssmaille)
914 | implicit none
915 | logical boolssmaille
916 | if (.not. boolssmaille) then
917 | print *, 'ERROR: La description de la partie energetique'
918 | print *, ' en luminance ne peut pas se faire avec'
919 | print *, ' une forme CTE ou par INDEX'
920 | stop
921 | endif
922 | return
923 | end
trajet.f could be called by:
| mcecile.f | [archivage/code2000X_testCG] | - 753 - 855 - 969 - 1107 - 1230 |
| mcecile.f | [resultats/pt1_complet] | - 710 - 807 - 913 - 1044 - 1162 |
| mcecile.f | [src] | - 1015 - 1154 - 1283 - 1454 - 1612 |