1 | PROGRAM mosim
2 |
3 |
4 | IMPLICIT NONE
5 | c
6 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
7 | c Modele approche de rayonnement: methode pertubationnelle
8 | c dst 1er ordre en kgaz et ksuies
9 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
10 | c
11 | c.......................................................................
12 | c
13 | INCLUDE 'cecile.inc'
14 | c
15 | c.......................................................................
16 | c
17 | c fichier_rayons : nom du fichier qui contient les rayons externes
18 | c de chaque couronne
19 | c fichier_ks_suie : nom du fichier qui contient les coefficients
20 | c d'absorption pour la suie dans chaque couronne
21 | c fichier_kgbar_gaz : nom du fichier qui contient les coefficients
22 | c d'absorption moyen pour le gaz dans chaque couronne
23 | c fichier_n_tirage : nom du fichier qui contient les nombres de
24 | c tirages effectues depuis chaque couronne
25 | c
26 | c.......................................................................
27 | c
28 | c......declaration pour black
29 | c implicit double precision (a-h,o-z)
30 | double precision sigma,c0,h,cbol,rind,c,c1,c2,blanu
31 |
32 | c......fin declaration pour black
33 |
34 |
35 |
36 | CHARACTER*100 fichier_rayons
37 | CHARACTER*100 fichier_ks_suie
38 | CHARACTER*100 fichier_kgbar_gaz
39 | CHARACTER*100 fichier_n_tirage
40 | character*100 fichier_temp
41 | character*100 texte
42 | integer a, b
43 | c
44 | c.......................................................................
45 | c
46 | DOUBLE PRECISION tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8
47 | INTEGER itmp1,itmp2,itmp3,itmp4,itmp5
48 | c
49 | c.......................................................................
50 | c
51 | c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
52 | c
53 | c.......................................................................
54 | c definition de la configuration, des proprietes physiques
55 | c et des controles du MCM
56 | c.......................................................................
57 | c
58 | print *,'LECTURE de LA CONFIG'
59 | read (*,*) a
60 | OPEN(10,FILE='cecile.in',STATUS='old')
61 | c
62 | READ (10,*)
63 | READ (10,*)
64 | READ (10,*)
65 | READ (10,*) n
66 | READ (10,*)
67 | IF (n.GT.n_mx)
68 | & PRINT *,'*************** ALERTE n ***************'
69 | c
70 | READ (10,*)
71 | READ (10,'(a)') fichier_rayons
72 | READ (10,*)
73 | OPEN(20,FILE=fichier_rayons,STATUS='old')
74 | READ (20,*) r(0)
75 | DO in=1,n
76 | READ (20,*) r(in)
77 | ENDDO
78 | CLOSE(20)
79 | r(n+1)=r(n)
80 | c
81 | READ (10,*)
82 | READ (10,'(a)') fichier_ks_suie
83 | READ (10,*)
84 | OPEN(20,FILE=fichier_ks_suie,STATUS='old')
85 | DO in=1,n
86 | READ (20,*) ks(in)
87 | ENDDO
88 | CLOSE(20)
89 | c
90 | READ (10,*)
91 | READ (10,'(a)') fichier_kgbar_gaz
92 | READ (10,*)
93 | OPEN(20,FILE=fichier_kgbar_gaz,STATUS='old')
94 | DO in=1,n
95 | READ (20,*) kgbar(in)
96 | ENDDO
97 | CLOSE(20)
98 | c
99 | READ (10,*)
100 | READ (10,*) phig
101 | READ (10,*)
102 | c
103 | READ (10,*)
104 | READ (10,'(a)') fichier_n_tirage
105 | READ (10,*)
106 | OPEN(20,FILE=fichier_n_tirage,STATUS='old')
107 | DO in=0,n+1
108 | READ (20,*) ntir(in)
109 | ENDDO
110 | CLOSE(20)
111 | c
112 | read (10,*)
113 | read (10,'(a)') fichier_temp
114 | read (10,*)
115 | c je prefere mettre la lecture du fichier dont on a eu le nom apres
116 | CLOSE(10)
117 |
118 | c lecture du fichier de temperature en K
119 | open (20,FILE=fichier_temp, STATUS='old')
120 | do in=0,n+1
121 | read(20,*) temp(in)
122 | enddo
123 | close (20)
124 | c.....lecture des temperatures en celsius
125 | c......les parois noires
126 | c temp(0)=temp(1)
127 | c.......l`exterieure refroidie
128 | c temp(n+1)=800
129 |
130 |
131 |
132 |
133 |
134 |
135 | c
136 | c.......................................................................
137 | c mise a zero de toutes les variables de sortie
138 | c.......................................................................
139 | c
140 | DO in=1,n
141 | bilan(in)=0.D+0
142 | var_bilan(in)=0.D+0
143 | ENDDO
144 | c
145 | DO in=0,n+1
146 | DO iin=0,n+1
147 | psi(in,iin)=0.D+0
148 | psir(in,iin)=0.D+0
149 | var_psi(in,iin)=0.D+0
150 | DO itmp1=1,n
151 | dpsidk(in,iin,itmp1)=0.D+0
152 | var_dpsidk(in,iin,itmp1)=0.D+0
153 | dpsidkgbar(in,iin,itmp1)=0.D+0
154 | var_dpsidkgbar(in,iin,itmp1)=0.D+0
155 | dkgbar(itmp1)=0.D+0
156 | dk(itmp1)=0.D+0
157 | ENDDO
158 | ENDDO
159 | ENDDO
160 |
161 |
162 |
163 | print *,'LECTURECECILE_sol de ref'
164 | c
165 | c.......................................................................
166 | c lectures des sorties DANS ECRAN de cecile psi et les derivees premiere/k
167 | c.......................................................................
168 | c
169 | c OPEN(10,FILE='cecile.out')
170 | c DO in=0,n+1
171 | c WRITE(10,*) 'bilan et ecart type en ',in
172 | c WRITE(10,*) bilan(in),var_bilan(in)
173 | c ENDDO
174 | c WRITE(10,*)
175 | c DO in=0,n+1
176 | c DO iin=0,n+1
177 | c WRITE(10,*) 'psi et ecart type entre ',in,' et ',iin
178 | c WRITE(10,*) psi(in,iin),var_psi(in,iin)
179 | c ENDDO
180 | c ENDDO
181 | c WRITE(10,*)
182 | c WRITE(10,*) 'et les DERIVEES ordre1 psi/k et leur ecart type'
183 | c DO in=0,n+1
184 | c write(10,*) 'entre ',in
185 | c DO iin=0,n+1
186 | c write(10,*) ' ...et ',iin
187 | c DO itmp1=1,n
188 | c write(10,*) ' ?itmp1 ',itmp1
189 | c WRITE(10,*) 'suies ',dpsidk(in,iin,itmp1),var_dpsidk(in,iin,itmp1)
190 | c WRITE(10,*) 'gaz ',dpsidkgbar(in,iin,itmp1)
191 | c & ,var_dpsidkgbar(in,iin,itmp1)
192 | c ENDDO
193 | c ENDDO
194 | c ENDDO
195 | c CLOSE(10)
196 |
197 | c...........................................................................
198 | c ......................lectures non commentees pour traitement rapide sft 98
199 | c...........................................................................
200 | c
201 | c......les psi de reference entre zones en tableau 2-d
202 | open(20,file='cecileout.don')
203 | c write(20,*) 'in iin psi ecart_type(dpsikg1_etc dpsiks1_etc)'
204 | read(20,*) texte
205 | do in=0,n+1
206 | do iin=0,n+1
207 | c write(20,*) in, iin, psi(in,iin)
208 | read(20,*) a, b, psi(in,iin), var_psi(in,iin)
209 | enddo
210 | enddo
211 | c.......les dpsi par rapport KGAZ
212 | c write(20,*) 'DPSIKGaz1_etc'
213 | read(20,*) texte
214 | do in=0,n+1
215 | do iin=0,n+1
216 | c write(20,*) '***********',in,iin
217 | read(20,*) texte, a, b
218 | do itmp1=1,n
219 | c write(20,*) dpsidkgbar(in,iin,itmp1)
220 | read(20,*) dpsidkgbar(in,iin,itmp1)
221 | enddo
222 | enddo
223 | enddo
224 |
225 | c.......les dpsi par rapport KSUIE
226 | c write(20,*) 'DPSIKSuie1_etc'
227 | read(20,*) texte
228 | do in=0,n+1
229 | do iin=0,n+1
230 | c write(20,*) '***********',in,iin
231 | read(20,*) texte, a, b
232 | do itmp1=1,n
233 | c write(20,*) dpsidk(in,iin,itmp1)
234 | read(20,*) dpsidk(in,iin,itmp1)
235 | enddo
236 | enddo
237 | enddo
238 | close(20)
239 | c.......................................................................
240 | c.................SAISIE de DKGBAR et DK................................
241 | do itmp1=1,n
242 | print *,'dk(suie):numero...',itmp1
243 | read *,dk(itmp1)
244 | enddo
245 |
246 | do itmp1=1,n
247 | print *,'dkgbar:numero...',itmp1
248 | read *,dkgbar(itmp1)
249 | enddo
250 |
251 | c.......................................................................
252 | c...ci-dessus on a saisie les valeurs de S(pour 1m de long)*ksi(in,iin)
253 | c........maintenant calcul de la fraction denergie emise sur............
254 | c........la bande etroite choisie de 1000 a 2000 cm-1...................
255 | c........c'est une reproduction du prog black...........................
256 | c***calcul de l'emittance monochronatique du CN a nbre donde et temp donn
257 | sigma=5.6693D-8
258 | c pi=datan(1.D+0)*4.D+0
259 | c0=2.9979D+8
260 | h=6.6262D-34
261 | cbol=1.3806D-23
262 | rind=1.D+0
263 | c=c0/rind
264 | c1=h*(c**2)
265 | c2=h*c/cbol
266 |
267 | c...les temperatures sont dans temp:saisie fichier(K)/blanu(m-1) est ici
268 | blanu=150000
269 | do in=0,n+1
270 | c je reutilise la var temp !!!!tres mauvais
271 | temp(in)=2.D+0*pi*c1*blanu**3/(dexp(c2*blanu/(temp(in)+273))-1.D+0)
272 | c temp represente maintenant l`emittance!!!!!!!!!!
273 | enddo
274 | c....pour avoir les luminances qui doivent venir en facteur de cecile
275 | c....je dois divise par pi et multiplier par deltanu: la bande etroite
276 | c....le passage des nomdbre d'onde a nu se fait nu=c*nbr onde
277 | c....PB: indice optique pqs connues donc 1 on va a vitesse lum du vide
278 | c....grosse erreur!!!!!!!!!!!!?????????????
279 | c deltanu=3.D+8*100000
280 | deltanu=2500
281 |
282 |
283 |
284 |
285 |
286 |
287 |
288 | c........................................................................
289 | c........................................................................
290 | c........................................................................
291 | print *,'DEBUTdeBOUCLE_calcul psi_cas reel ET ecart_type'
292 | read (*,*) a
293 | c a='****************'
294 | c........................................................................
295 | c boucle sur chaque psi in iin
296 | OPEN(30,FILE='ecran2000')
297 | do in=0,n+1
298 | do iin=0,n+1
299 | c write(30,*) in,iin
300 | psir(in,iin)=psi(in,iin)
301 | c ci-dessus premier terme factoriel 0
302 | c boucle sur les derivees
303 | do itmp1=1,n
304 |
305 | psir(in,iin)=psir(in,iin)+(dpsidkgbar(in,iin,itmp1)
306 | & *dkgbar(itmp1)) +(dpsidk(in,iin,itmp1))*dk(itmp1)
307 |
308 | enddo
309 | psir(in,iin)=psir(in,iin)*(temp(in)-temp(iin))*deltanu/pi
310 | c meme traitement pour les ecart types
311 | var_psi(in,iin)=var_psi(in,iin)*(temp(in)-temp(iin))*deltanu/pi
312 | c pour que les parois ne soient pas nulles a cause division selective par volume surface
313 | psirm(in,iin)=psir(in,iin)
314 | if ((in.ne.0.).and.(in.ne.(n+1))) then
315 | psirm(in,iin)=psir(in,iin)/(1*abs(pi*((r(in)**2-r(in-1)**2))))
316 | var_psi(in,iin)=var_psi(in,iin)/(1*abs(pi*((r(in)**2-r(in-1)**2))))
317 | endif
318 | c puissance en W/m3 de la maille courante
319 | if ((iin.ne.0.).and.(iin.ne.(n+1))) then
320 | psirm(in,iin)=psirm(in,iin)/(abs(pi*((r(iin)**2-r(iin-1)**2))))
321 | var_psi(in,iin)=var_psi(in,iin)/(abs(pi*((r(iin)**2-r(iin-1)**2))))
322 | endif
323 | c puissance en W/m3.m2 de la maille incrementee
324 | c write(30,*) in, iin, psirm(in,iin)
325 |
326 | enddo
327 | c PB formatage tableau fixe a n+3 mais nombre indispensable en f77
328 | c write(30,'(11(f7.3))') in, (psirm(in,iin), iin=0,n+1)
329 | write(30,'(i3,10(e12.3))') in, (psirm(in,iin), iin=0,n+1)
330 | c write(30,*) in, (psirm(in,iin), iin=0,n+1)
331 | enddo
332 | write(30,*) '********************les ecarts types'
333 | do in=0,n+1
334 | write(30,*) in, (var_psi(in,iin), iin=0,n+1)
335 | enddo
336 |
337 |
338 |
339 |
340 |
341 |
342 |
343 | c
344 | c.......................................................................
345 | c
346 | END