1 | #include "config.h" |
---|
2 | |
---|
3 | module sous_routines |
---|
4 | |
---|
5 | use parametres |
---|
6 | use fft |
---|
7 | |
---|
8 | implicit none |
---|
9 | |
---|
10 | contains |
---|
11 | |
---|
12 | !********************************************************************************** |
---|
13 | !********************************************************************************** |
---|
14 | !********** SOUS ROUTINES D'ENTREE SORTIES ********** |
---|
15 | !********************************************************************************** |
---|
16 | !********************************************************************************** |
---|
17 | |
---|
18 | ! *************************************************************************** |
---|
19 | subroutine enter_param |
---|
20 | ! *************************************************************************** |
---|
21 | ! lecture des tous les parametres qui sont dans le fichier 'data.in' |
---|
22 | ! *************************************************************************** |
---|
23 | |
---|
24 | character(len=29) :: bidon |
---|
25 | character(len=100) :: lecture_reel, lecture_entier, lecture_logical |
---|
26 | integer :: i, nx,ny,nz |
---|
27 | |
---|
28 | lecture_reel = '(A29,E14.5)' |
---|
29 | lecture_entier = '(A29,I10)' |
---|
30 | lecture_logical = '(A29,L12)' |
---|
31 | |
---|
32 | open (unit=1,file='data.in') |
---|
33 | |
---|
34 | |
---|
35 | |
---|
36 | read(1,*) |
---|
37 | read(1,*) |
---|
38 | read(1,*) |
---|
39 | read(1,*) |
---|
40 | |
---|
41 | ! dimension du probleme (ces valeurs ne sont en fait pas reellement utilisees). |
---|
42 | ! les dimensions du probleme sont definies dans le fichier parametres.F90 par le preprocesseur |
---|
43 | ! en fonction des valeurs indiquees dans le fichier config.h |
---|
44 | read(1,lecture_entier) bidon,nx |
---|
45 | read(1,lecture_entier) bidon,ny |
---|
46 | read(1,lecture_entier) bidon,nz |
---|
47 | |
---|
48 | ! taille suivant x |
---|
49 | read(1,lecture_reel) bidon,lx |
---|
50 | read(1,lecture_reel) bidon,ly |
---|
51 | read(1,lecture_reel) bidon,lz |
---|
52 | |
---|
53 | ! Pas de temps et instant initial |
---|
54 | read(1,lecture_reel) bidon,dt |
---|
55 | read(1,lecture_reel) bidon,begin |
---|
56 | |
---|
57 | ! nombre d'itérations |
---|
58 | read(1,lecture_entier) bidon,itmax |
---|
59 | |
---|
60 | ! Troncature |
---|
61 | read(1,lecture_logical) bidon,truncate |
---|
62 | read(1,lecture_entier) bidon,type_trunc |
---|
63 | read(1,lecture_reel) bidon,rtrunc_x |
---|
64 | read(1,lecture_reel) bidon,rtrunc_y |
---|
65 | read(1,lecture_reel) bidon,rtrunc_z |
---|
66 | |
---|
67 | read(1,*) |
---|
68 | read(1,*) |
---|
69 | |
---|
70 | ! Viscosité |
---|
71 | read(1,lecture_reel) bidon,nu |
---|
72 | |
---|
73 | ! stratification, fréquence de Brunt Vaisala, vorticité planétaire et nombre de Schmidt |
---|
74 | read(1,lecture_logical) bidon,stratification |
---|
75 | read(1,lecture_reel) bidon,xns |
---|
76 | xnsqr=xns*xns |
---|
77 | read(1,lecture_reel) bidon,schm |
---|
78 | read(1,lecture_reel) bidon,omega2 |
---|
79 | |
---|
80 | read(1,*) |
---|
81 | read(1,*) |
---|
82 | |
---|
83 | ! type de run |
---|
84 | read(1,lecture_logical) bidon,perturbe |
---|
85 | read(1,lecture_logical) bidon,lin |
---|
86 | |
---|
87 | read(1,*) |
---|
88 | read(1,*) |
---|
89 | |
---|
90 | ! reprise de run |
---|
91 | read(1,lecture_logical) bidon,reprise_run |
---|
92 | |
---|
93 | read(1,*) |
---|
94 | read(1,*) |
---|
95 | |
---|
96 | ! etat de base 2D |
---|
97 | read(1,lecture_logical) bidon, base2D_sub |
---|
98 | read(1,lecture_logical) bidon, base2D_read |
---|
99 | |
---|
100 | read(1,*) |
---|
101 | read(1,*) |
---|
102 | |
---|
103 | ! vitesse initiale 2D |
---|
104 | read(1,lecture_logical) bidon, velo2D_sub |
---|
105 | read(1,lecture_logical) bidon, velo2D_read |
---|
106 | read(1,lecture_logical) bidon,lnoise |
---|
107 | read(1,lecture_reel) bidon,amplitude_noise |
---|
108 | |
---|
109 | read(1,*) |
---|
110 | read(1,*) |
---|
111 | |
---|
112 | |
---|
113 | !profil de densite |
---|
114 | read(1,lecture_logical) bidon, Rho_read |
---|
115 | read(1,lecture_reel) bidon, Rhob |
---|
116 | |
---|
117 | read(1,*) |
---|
118 | read(1,*) |
---|
119 | |
---|
120 | ! pba : deb |
---|
121 | ! forcage |
---|
122 | read(1,lecture_logical) bidon, forcage_sub |
---|
123 | read(1,lecture_reel) bidon, A_ics |
---|
124 | read(1,lecture_reel) bidon, freq |
---|
125 | read(1,lecture_entier) bidon, it0 |
---|
126 | |
---|
127 | read(1,*) |
---|
128 | read(1,*) |
---|
129 | |
---|
130 | !pba fin |
---|
131 | |
---|
132 | ! sorties |
---|
133 | read(1,lecture_entier) bidon,periode_out_physique1 |
---|
134 | read(1,lecture_entier) bidon,periode_out_physique2 |
---|
135 | read(1,lecture_entier) bidon,periode_out_spectral1 |
---|
136 | |
---|
137 | ! parametres additionnels |
---|
138 | read(1,*) |
---|
139 | read(1,*) |
---|
140 | |
---|
141 | ! Tourbillons gaussiens |
---|
142 | read(1,*) |
---|
143 | read(1,*) |
---|
144 | |
---|
145 | read(1,lecture_entier) bidon,nb_tourbillons |
---|
146 | do i=1,nb_tourbillons |
---|
147 | read(1,*) |
---|
148 | read(1,lecture_reel) bidon,xt(i) |
---|
149 | read(1,lecture_reel) bidon,yt(i) |
---|
150 | read(1,lecture_reel) bidon,gamma(i) |
---|
151 | read(1,lecture_reel) bidon,at(i) |
---|
152 | end do |
---|
153 | |
---|
154 | close(1) |
---|
155 | |
---|
156 | end subroutine |
---|
157 | |
---|
158 | ! *************************************************************************** |
---|
159 | subroutine out_param |
---|
160 | ! *************************************************************************** |
---|
161 | ! affichage des parametres du run |
---|
162 | ! *************************************************************************** |
---|
163 | |
---|
164 | implicit none |
---|
165 | |
---|
166 | character(len=100) :: aff_reel, aff_entier, aff_logical |
---|
167 | integer :: i |
---|
168 | |
---|
169 | ! descripteur des formats de sortie |
---|
170 | aff_reel = '(A29,TR3,G12.5)' |
---|
171 | aff_entier = '(A29,TR4,I7)' |
---|
172 | aff_logical = '(A29,TR4,L1)' |
---|
173 | |
---|
174 | write(*,*) |
---|
175 | write(*,*) '---------------------------------------------------' |
---|
176 | write(*,*) ' PROGRAMME NS3D' |
---|
177 | write(*,*) '---------------------------------------------------' |
---|
178 | write(*,*) |
---|
179 | write(*,*) |
---|
180 | ! Mode de calcul |
---|
181 | if (lin) then |
---|
182 | write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' |
---|
183 | write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' |
---|
184 | write(*,*) '@@@ @@@' |
---|
185 | write(*,*) ' VERSION PERTURBATIVE / LINEARISEE ' |
---|
186 | write(*,*) '@@@ @@@' |
---|
187 | write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' |
---|
188 | write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' |
---|
189 | else |
---|
190 | if (perturbe) then |
---|
191 | write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' |
---|
192 | write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' |
---|
193 | write(*,*) '@@@ @@@' |
---|
194 | write(*,*) ' VERSION PERTURBATIVE / NON LINEARISEE ' |
---|
195 | write(*,*) '@@@ @@@' |
---|
196 | write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' |
---|
197 | write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' |
---|
198 | else |
---|
199 | write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' |
---|
200 | write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' |
---|
201 | write(*,*) '@@@ @@@' |
---|
202 | write(*,*) ' VERSION NON-PERTURBATIVE ' |
---|
203 | write(*,*) '@@@ @@@' |
---|
204 | write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' |
---|
205 | write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' |
---|
206 | end if |
---|
207 | end if |
---|
208 | |
---|
209 | write(*,*) |
---|
210 | write(*,*) |
---|
211 | write(*,*) '---------------------------------------------------' |
---|
212 | write(*,*) 'PARAMETRES DE DISCRETISATION' |
---|
213 | write(*,*) '---------------------------------------------------' |
---|
214 | write(*,aff_entier) 'nx...........................',nxmax |
---|
215 | write(*,aff_entier) 'ny...........................',nymax |
---|
216 | write(*,aff_entier) 'nz...........................',nzmax |
---|
217 | write(*,aff_reel) 'Lx...........................',lx |
---|
218 | write(*,aff_reel) 'Ly...........................',ly |
---|
219 | write(*,aff_reel) 'Lz...........................',lz |
---|
220 | write(*,aff_reel) 'Dt...........................',dt |
---|
221 | write(*,aff_reel) 'begin........................',begin |
---|
222 | write(*,aff_entier) 'itmax........................',itmax |
---|
223 | write(*,aff_logical) 'Troncature...................',truncate |
---|
224 | write(*,aff_entier) ' carree:1 / elliptique:2..',type_trunc |
---|
225 | write(*,aff_reel) ' rayon troncature x.......',rtrunc_x |
---|
226 | write(*,aff_reel) ' rayon troncature y.......',rtrunc_y |
---|
227 | write(*,aff_reel) ' rayon troncature z........',rtrunc_z |
---|
228 | |
---|
229 | write(*,*) |
---|
230 | write(*,*) '---------------------------------------------------' |
---|
231 | write(*,*) 'PARAMETRES PHYSIQUES' |
---|
232 | write(*,*) '---------------------------------------------------' |
---|
233 | write(*,aff_reel) 'viscosite....................',nu |
---|
234 | write(*,aff_logical) 'stratification...............',stratification |
---|
235 | write(*,aff_reel) ' brunt_vaisala_frequency..',xns |
---|
236 | write(*,aff_reel) ' Schmidt number...........',schm |
---|
237 | write(*,aff_reel) 'Omega2.......................',omega2 |
---|
238 | |
---|
239 | write(*,*) |
---|
240 | write(*,*) '---------------------------------------------------' |
---|
241 | write(*,*) 'TYPE DE RUN' |
---|
242 | write(*,*) '---------------------------------------------------' |
---|
243 | write(*,aff_logical) 'en perturbation..............',perturbe |
---|
244 | write(*,aff_logical) 'lineaire.....................',lin |
---|
245 | |
---|
246 | write(*,*) |
---|
247 | write(*,*) '---------------------------------------------------' |
---|
248 | write(*,*) 'REPRISE DE RUN' |
---|
249 | write(*,*) '---------------------------------------------------' |
---|
250 | write(*,aff_logical) 'reprise de run...............',reprise_run |
---|
251 | |
---|
252 | write(*,*) |
---|
253 | write(*,*) '---------------------------------------------------' |
---|
254 | write(*,*) 'ETAT DE BASE 2D' |
---|
255 | write(*,*) '---------------------------------------------------' |
---|
256 | write(*,aff_logical) 'call subroutine..............',base2D_sub |
---|
257 | write(*,aff_logical) 'Lecture fichier..............',base2D_read |
---|
258 | |
---|
259 | |
---|
260 | write(*,*) |
---|
261 | write(*,*) '---------------------------------------------------' |
---|
262 | write(*,*) 'VITESSE INITIALE 2D' |
---|
263 | write(*,*) '---------------------------------------------------' |
---|
264 | write(*,aff_logical) 'call subroutine..............',velo2D_sub |
---|
265 | write(*,aff_logical) 'Lecture fichier..............',velo2D_read |
---|
266 | write(*,aff_logical) 'Bruit blanc..................',lnoise |
---|
267 | write(*,aff_reel) ' amplitude................',amplitude_noise |
---|
268 | |
---|
269 | write(*,*) |
---|
270 | write(*,*) '---------------------------------------------------' |
---|
271 | write(*,*) 'DENSITE INITIALE 1D' |
---|
272 | write(*,*) '---------------------------------------------------' |
---|
273 | write(*,aff_logical) 'Lecture fichier..............',Rho_read |
---|
274 | write(*,aff_reel) 'densite de base..............',Rhob |
---|
275 | |
---|
276 | |
---|
277 | write(*,*) |
---|
278 | write(*,*)'---------------------------------------------------' |
---|
279 | write(*,*)'FORCAGE EN VITESSE' |
---|
280 | write(*,*)'---------------------------------------------------' |
---|
281 | write(*,aff_logical)'call subroutine', forcage_sub |
---|
282 | write(*,aff_reel)'Amplitude..................',A_ics |
---|
283 | write(*,aff_reel)'Frequence..................',freq |
---|
284 | |
---|
285 | write(*,*) |
---|
286 | write(*,*) '---------------------------------------------------' |
---|
287 | write(*,*) 'SORTIES' |
---|
288 | write(*,*) '---------------------------------------------------' |
---|
289 | write(*,aff_entier) 'periode out_physique1.........',periode_out_physique1 |
---|
290 | write(*,aff_entier) 'periode out_physique2........',periode_out_physique2 |
---|
291 | write(*,aff_entier) 'periode out_spectral1........',periode_out_spectral1 |
---|
292 | |
---|
293 | write(*,*) |
---|
294 | write(*,*) '---------------------------------------------------' |
---|
295 | write(*,*) 'PARAMETRES ADDITIONNELS' |
---|
296 | write(*,*) '---------------------------------------------------' |
---|
297 | write(*,*) |
---|
298 | write(*,*) '** Tourbillons gaussiens **' |
---|
299 | write(*,aff_entier) 'nb tourbillons...............',nb_tourbillons |
---|
300 | do i=1,nb_tourbillons |
---|
301 | write(*,aff_entier) 'Tourbillon...................',i |
---|
302 | write(*,aff_reel) ' xt.......................',xt(i) |
---|
303 | write(*,aff_reel) ' xt.......................',yt(i) |
---|
304 | write(*,aff_reel) ' gamma....................',gamma(i) |
---|
305 | write(*,aff_reel) ' at.......................',at(i) |
---|
306 | end do |
---|
307 | |
---|
308 | write(*,*) |
---|
309 | write(*,*) |
---|
310 | |
---|
311 | end subroutine |
---|
312 | |
---|
313 | |
---|
314 | |
---|
315 | ! ********************************************************************** |
---|
316 | subroutine out_physique1(vox,voy,voz) |
---|
317 | ! ********************************************************************** |
---|
318 | ! Sortie effectuee de façon periodique et qui prend en argument le champ de vitesse |
---|
319 | ! dans l'espace physique. |
---|
320 | ! Ici on affiche le taux de croissance des perturbations ainsi que le temps cpu ecoule et |
---|
321 | ! une estimation du temps restant |
---|
322 | ! ********************************************************************** |
---|
323 | |
---|
324 | implicit none |
---|
325 | |
---|
326 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: vox, voy, voz |
---|
327 | double precision :: vquad, growthrate |
---|
328 | double precision, save :: tps_ecoule = 0, tps_ecoule_prec_iteration = 0 |
---|
329 | |
---|
330 | |
---|
331 | ! on calcule le temps cpu ecoule depuis la rentree dans la boucle et une estimation du temps |
---|
332 | ! restant (on se base sur le temps necessaire a la derniere iteration de la boucle) |
---|
333 | call cpu_time(tps_temp) |
---|
334 | tps_ecoule = tps_total + tps_temp |
---|
335 | |
---|
336 | write(*,*) |
---|
337 | write(*,'(A4,TR2,I6,TR3,A6,TR2,F10.3)') 'it =', it, 'time =', it*dt + begin |
---|
338 | write(*,'(A37,TR2,I10,TR3,I10)') ' temps cpu en sec (ecoule/restant) =', int(tps_ecoule), & |
---|
339 | int((tps_ecoule-tps_ecoule_prec_iteration)/periode_out_physique1*(itmax-it)) |
---|
340 | |
---|
341 | tps_ecoule_prec_iteration = tps_ecoule |
---|
342 | |
---|
343 | |
---|
344 | ! on calcule la vitesse quadratique moyenne des perturbations |
---|
345 | vquad = sum(vox(0:nxmax-1, 0:nymax-1, 0:nzmax-1)**2 + voy(0:nxmax-1, 0:nymax-1, 0:nzmax-1)**2 + & |
---|
346 | voz(0:nxmax-1, 0:nymax-1, 0:nzmax-1)**2) |
---|
347 | |
---|
348 | write(*,'(A32,TR2,G23.10)') ' vitesse quadratique moyenne =',vquad |
---|
349 | |
---|
350 | ! on calcule le taux de croissance des perturbations entre deux appels à la fonction out_physique1 |
---|
351 | if (vquadanc>epsilo) then |
---|
352 | growthrate=0.5*log(vquad/vquadanc)/(periode_out_physique1*dt) |
---|
353 | else |
---|
354 | growthrate=0. |
---|
355 | endif |
---|
356 | |
---|
357 | write(*,'(A34,TR2,G23.10)') ' taux de croissance instantane =', growthrate |
---|
358 | |
---|
359 | vquadanc = vquad |
---|
360 | |
---|
361 | end subroutine |
---|
362 | |
---|
363 | ! ************************************************************************* |
---|
364 | subroutine out_physique2(vox,voy,voz,rho,fx,fy,fz) |
---|
365 | ! ************************************************************************* |
---|
366 | ! Sortie effectuee de façon periodique et qui prend en argument le champ de vitesse |
---|
367 | ! la densite et la vorticite dans l'espace physique. |
---|
368 | ! ************************************************************************** |
---|
369 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: vox,voy,voz,rho |
---|
370 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx,fy,fz |
---|
371 | |
---|
372 | ! là on appelle simplement, une fonction qui écrit les champs dans un fichier. |
---|
373 | call out_champs(vox,voy,voz,rho,fx,fy,fz) |
---|
374 | |
---|
375 | end subroutine |
---|
376 | |
---|
377 | |
---|
378 | ! ************************************************************************* |
---|
379 | subroutine out_spectral1(ax,ay,az) |
---|
380 | ! ************************************************************************* |
---|
381 | ! Fonctions qui sauvegarde l'energie quadratique moyenné suivant ky dans l'espace spectral |
---|
382 | ! On fait d'abord une moyenne suivant ky et on enregistre l'energie pour chaque mode kx et kz |
---|
383 | ! ************************************************************************* |
---|
384 | |
---|
385 | implicit none |
---|
386 | |
---|
387 | double precision, dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: ax, ay, az |
---|
388 | double precision, dimension(0:kxmax, 0:nzmax-1) :: energie_quadratique_mode |
---|
389 | integer :: ix, iz |
---|
390 | |
---|
391 | |
---|
392 | do iz=0,nzmax-1 |
---|
393 | do ix=0,kxmax |
---|
394 | energie_quadratique_mode(ix, iz) = sum(ax(1,ix,0:nymax-1,iz)**2+ax(2,ix,0:nymax-1,iz)**2+ & |
---|
395 | ay(1,ix,0:nymax-1,iz)**2+ ay(2,ix,0:nymax-1,iz)**2+ & |
---|
396 | az(1,ix,0:nymax-1,iz)**2+az(2,ix,0:nymax-1,iz)**2)/nymax |
---|
397 | end do |
---|
398 | end do |
---|
399 | |
---|
400 | |
---|
401 | ! on l'écrit dans un fichier |
---|
402 | |
---|
403 | ! si c'est le premier appel à la fonction out_spectral1, on cree le fichier de sortie et on ecrit toutes les |
---|
404 | ! metadonnees. |
---|
405 | ! On pourrait eventuellement tester si le fichier existe deja inquire(file="energie_quadratique_mode", EXIST=existe) |
---|
406 | ! Cela permettrait de completer un fichier de sortie genere lors d'une precedente simulation... |
---|
407 | |
---|
408 | ! S'il ne s'agit pas du premier appel, on ecrit simplement l'instant et les valeurs de l'energie_quadratique_mode |
---|
409 | ! dans le fichier de sortie |
---|
410 | |
---|
411 | ! premier appel a periode_out_spectral1: on cree le fichier et on ecrit les metadonnees |
---|
412 | if (it == periode_out_spectral1) then |
---|
413 | open(unit=76, file='energie_quadratique_mode',form='unformatted',status="replace",action="write") |
---|
414 | write(76) 1, nxmax, nymax, nzmax, lx, ly, lz, dt, truncate, type_trunc, rtrunc_x, rtrunc_y, rtrunc_z, & |
---|
415 | nu, stratification, xns, schm, omega2, perturbe, lin |
---|
416 | close(76) |
---|
417 | end if |
---|
418 | |
---|
419 | ! ce n'est pas le premier appel : on ajoute a la fin du fichier existant l'instant et les valeurs de l'energie |
---|
420 | ! Si on ajoute l'argument action="write" dans la foncion write avec le compilateur Absoft Pro Fortran 8.2 sous Mac |
---|
421 | ! l'enregistrement ne se fait pas correctement alors que ca marche sans probleme avec le compilateur Intel Fortran |
---|
422 | ! ou le compilateur Fortran du NEC de l'IDRIS, bizarre, bizarre... |
---|
423 | open(unit=76, file='energie_quadratique_mode',form='unformatted',status="old",position="append") |
---|
424 | write(76) it*dt + begin |
---|
425 | write(76) energie_quadratique_mode |
---|
426 | close(76) |
---|
427 | |
---|
428 | end subroutine |
---|
429 | |
---|
430 | |
---|
431 | ! ************************************************************************* |
---|
432 | subroutine out_champs(vox,voy,voz,rho,fx,fy,fz) |
---|
433 | ! ************************************************************************* |
---|
434 | ! Assure la sortie des champs dans des fichiers |
---|
435 | ! ************************************************************************** |
---|
436 | |
---|
437 | implicit none |
---|
438 | |
---|
439 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: vox,voy,voz,rho |
---|
440 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx,fy,fz |
---|
441 | character(len=50) :: nomfic |
---|
442 | character(len=13) :: sufx |
---|
443 | |
---|
444 | call time2string(it*dt + begin,sufx) |
---|
445 | |
---|
446 | ! on ecrit la vitesse, la densite et la vorticite dans un fichier |
---|
447 | |
---|
448 | ! on écrit un certain nombre de metadonnees pour aider a identifier le fichier : |
---|
449 | ! c'est utile pour eviter bien des erreurs lors des entrees-sorties et faciliter le post-traitement... |
---|
450 | ! Cela permet notamment de s'assurer qu'on utilise le bon codage des données binaires (little- |
---|
451 | ! endian ou big-endian) et que les tableaux ont bien les dimensions requises. |
---|
452 | |
---|
453 | ! On écrit ensuite les sorties tableau par tableau dans le fichier final. |
---|
454 | |
---|
455 | ! Attention, Fortran ajoute au début et à la fin de chaque enregistrement, à chaque |
---|
456 | ! enregistrement binaire effectué avec l'instruction "write(,form='unformatted')", un entier |
---|
457 | ! correspondant a la taille de l'enregitrement en octets. Il faut faire attention |
---|
458 | ! et en tenir compte lors du traitement ulterieur éventuel des donnees (dans Matlab par ex) |
---|
459 | |
---|
460 | ! Le format final du fichier de sortie est donc (i et j sont des entiers correspondant |
---|
461 | ! à la taille de l'enregistrement qui respectivement suit et précède) : |
---|
462 | |
---|
463 | ! i,1,nxmax,nymax,nzmax,lx,ly,lz,truncate,type_trunc,rtrunc_x,rtrunc_y,rtrunc_z,nu,stratification, |
---|
464 | ! xns,schm,omega2,perturbe,lin,j,i,vox,j,i,voy,j,i,voz,j,i,rho,j,i,fx,j,i,fy,j,i,fz,j, |
---|
465 | ! i,vbx,vby,vbz,j,i,wbx,wby,wbz,j |
---|
466 | |
---|
467 | nomfic='velo_rho_vort.'//sufx |
---|
468 | |
---|
469 | open (76,file=nomfic,form='unformatted') |
---|
470 | write(76) 1, nxmax, nymax, nzmax, lx, ly, lz, dt, truncate, type_trunc, rtrunc_x, rtrunc_y, rtrunc_z, & |
---|
471 | nu, stratification, xns, schm, omega2, perturbe, lin |
---|
472 | write(76) it*dt + begin |
---|
473 | write(76) vox(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
474 | write(76) voy(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
475 | write(76) voz(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
476 | write(76) rho(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
477 | write(76) fx(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
478 | write(76) fy(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
479 | write(76) fz(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
480 | write(76) vbx,vby,vbz |
---|
481 | write(76) wbx,wby,wbz |
---|
482 | |
---|
483 | close(76) |
---|
484 | print *,'itfic:', it |
---|
485 | end subroutine |
---|
486 | |
---|
487 | ! *********************************************************** |
---|
488 | subroutine read_run(vox,voy,voz,rho,fx,fy,fz,vbx,vby,vbz,wbx,wby,wbz,begin) |
---|
489 | ! *********************************************************** |
---|
490 | ! Lecture des champs d'un run precedent dans le fichier nom_fichier |
---|
491 | ! *********************************************************** |
---|
492 | |
---|
493 | implicit none |
---|
494 | |
---|
495 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1),intent(out) :: vox,voy,voz,rho,fx,fy,fz |
---|
496 | double precision, dimension(0:nxmax-1,0:nymax-1), intent(out) :: vbx,vby,vbz,wbx,wby,wbz |
---|
497 | double precision, intent(out) :: begin |
---|
498 | integer :: flag, nxread, nyread, nzread |
---|
499 | double precision :: lxread, lyread, lzread |
---|
500 | double precision :: dbl_poub |
---|
501 | integer :: int_poub |
---|
502 | logical :: log_poub |
---|
503 | |
---|
504 | open (unit=8,file='velocity.init',form='unformatted') |
---|
505 | |
---|
506 | ! on lit les metadonnees |
---|
507 | read(8) flag,nxread,nyread,nzread,lxread,lyread,lzread,dbl_poub,log_poub,int_poub,dbl_poub,dbl_poub,dbl_poub, & |
---|
508 | dbl_poub,log_poub,dbl_poub,dbl_poub,dbl_poub,log_poub,log_poub |
---|
509 | read(8) begin |
---|
510 | |
---|
511 | ! on verifie que le format du fichier d'entree est correct. On ne verifie que les dimensions du probleme car ce sont |
---|
512 | ! les seuls parametres qui ne peuvent vraiment pas etre changes entre deux simulations |
---|
513 | if (flag /= 1 .or. nxread /= nxmax .or. nyread /= nymax .or. nzread /= nzmax .or. lxread /= lx .or. lyread /= ly) then |
---|
514 | write(*,*) "ERREUR: le fichier de reprise de run n'a pas le bon format" |
---|
515 | write(*,*) "ARRET DU PROGRAMME..." |
---|
516 | stop |
---|
517 | end if |
---|
518 | |
---|
519 | read(8) vox(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
520 | read(8) voy(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
521 | read(8) voz(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
522 | read(8) rho(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
523 | read(8) fx(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
524 | read(8) fy(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
525 | read(8) fz(0:nxmax-1,0:nymax-1,0:nzmax-1) |
---|
526 | read(8) vbx,vby,vbz |
---|
527 | read(8) wbx,wby,wbz |
---|
528 | |
---|
529 | close(8) |
---|
530 | |
---|
531 | end subroutine |
---|
532 | |
---|
533 | |
---|
534 | ! on choisit quelle fonction de generation de l'etat de base 2D doit etre utilisee |
---|
535 | #if choix_base2D==1 |
---|
536 | ! *********************************************************** |
---|
537 | subroutine gen_base2D(vbx,vby,vbz,wbx,wby,wbz,vox2,voy2,voz2) |
---|
538 | ! *********************************************************** |
---|
539 | ! Generation de l'ecoulement de base |
---|
540 | ! Ecoulement de Stuart |
---|
541 | ! *********************************************************** |
---|
542 | |
---|
543 | implicit none |
---|
544 | |
---|
545 | double precision, dimension(0:nxmax-1,0:nymax-1), intent(out) :: vbx,vby,vbz,wbx,wby,wbz |
---|
546 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: vox2, voy2, voz2 |
---|
547 | integer :: ix, iy |
---|
548 | double precision ::rho_stuart |
---|
549 | |
---|
550 | vbx = 0. |
---|
551 | vby = 0. |
---|
552 | vbz = 0. |
---|
553 | wbx = 0. |
---|
554 | wby = 0. |
---|
555 | wbz = 0. |
---|
556 | vox2 = 0. |
---|
557 | voy2 = 0. |
---|
558 | voz2 = 0. |
---|
559 | |
---|
560 | ! on définit ici le parametre de concentration des tourbillons de stuart |
---|
561 | rho_stuart = 0.25 |
---|
562 | |
---|
563 | do iy=0, nymax-1 |
---|
564 | do ix = 0,nxmax-1 |
---|
565 | vbx(ix,iy) = sinh(2*(yy(iy)-ly/2))/(cosh(2* (yy(iy)-ly/2) ) - rho_stuart*cos( 2*xx(ix) ) ) |
---|
566 | vby(ix,iy) = - rho_stuart * sin( 2*xx(ix) )/( cosh( 2*( yy(iy)-ly/2 ) ) - rho_stuart*cos(2*xx(ix))) |
---|
567 | wbz(ix,iy) = 2*(rho_stuart**2-1)/(cosh(2*(yy(iy)-ly/2)) - rho_stuart*cos(2*xx(ix)))**2 |
---|
568 | end do |
---|
569 | end do |
---|
570 | |
---|
571 | end subroutine |
---|
572 | |
---|
573 | |
---|
574 | #elif choix_base2D==2 |
---|
575 | ! *********************************************************** |
---|
576 | subroutine gen_base2D(vbx,vby,vbz,wbx,wby,wbz,vox2,voy2,voz2) |
---|
577 | ! *********************************************************** |
---|
578 | ! Generation de l'ecoulement de base |
---|
579 | ! Couche de mélange avec un profil vbx=tanh(y) |
---|
580 | ! *********************************************************** |
---|
581 | |
---|
582 | implicit none |
---|
583 | |
---|
584 | double precision, dimension(0:nxmax-1,0:nymax-1), intent(out) :: vbx,vby,vbz,wbx,wby,wbz |
---|
585 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: vox2, voy2, voz2 |
---|
586 | integer :: ix |
---|
587 | |
---|
588 | vbx = 0. |
---|
589 | vby = 0. |
---|
590 | vbz = 0. |
---|
591 | wbx = 0. |
---|
592 | wby = 0. |
---|
593 | wbz = 0. |
---|
594 | vox2 = 0. |
---|
595 | voy2 = 0. |
---|
596 | voz2 = 0. |
---|
597 | |
---|
598 | do ix = 0,nxmax-1 |
---|
599 | vbx(ix,:) = tanh(yy-ly/2) |
---|
600 | wbz(ix,:) = -(1-tanh(yy-ly/2)**2) |
---|
601 | end do |
---|
602 | |
---|
603 | end subroutine |
---|
604 | |
---|
605 | |
---|
606 | #elif choix_base2D==3 |
---|
607 | ! *********************************************************** |
---|
608 | subroutine gen_base2D(vbx,vby,vbz,wbx,wby,wbz,vox2,voy2,voz2) |
---|
609 | ! *********************************************************** |
---|
610 | ! Generation de l'ecoulement de base |
---|
611 | ! Superposition de tourbillons dont les caracteristiques |
---|
612 | ! sont donnees dans le fichier de parametres |
---|
613 | ! *********************************************************** |
---|
614 | |
---|
615 | implicit none |
---|
616 | |
---|
617 | double precision, dimension(0:nxmax-1,0:nymax-1), intent(out) :: vbx,vby,vbz,wbx,wby,wbz |
---|
618 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: vox2, voy2, voz2 |
---|
619 | double precision k2 |
---|
620 | integer num_tourbillon |
---|
621 | integer ix, iy, iz |
---|
622 | |
---|
623 | vbx = 0. |
---|
624 | vby = 0. |
---|
625 | vbz = 0. |
---|
626 | wbx = 0. |
---|
627 | wby = 0. |
---|
628 | wbz = 0. |
---|
629 | voz2 = 0. |
---|
630 | |
---|
631 | ! on determine la vorticite verticale de l'ecoulement |
---|
632 | do num_tourbillon=1,nb_tourbillons |
---|
633 | do iz=0,nzmax-1 |
---|
634 | do iy=0,nymax-1 |
---|
635 | do ix=0,nxmax-1 |
---|
636 | voz2(ix,iy,iz)=voz2(ix,iy,iz)+wz_gau((xx(ix)-xt(num_tourbillon)),(yy(iy)-yt(num_tourbillon)), & |
---|
637 | gamma(num_tourbillon),at(num_tourbillon)) |
---|
638 | end do |
---|
639 | end do |
---|
640 | end do |
---|
641 | end do |
---|
642 | |
---|
643 | ! initialisation de la vorticite verticale |
---|
644 | wbz = voz2(0:nxmax-1,0:nymax-1,0) |
---|
645 | |
---|
646 | ! inversion du laplacien 2D pour trouver la fct de courant 2D*1D |
---|
647 | call fwd_fft(voz2,fft_work) |
---|
648 | |
---|
649 | do iz=0,nzmax-1 |
---|
650 | do iy=0,nymax-1 |
---|
651 | do ix=0,2*kxmax+1 |
---|
652 | k2=max(ca2x(ix)+k2y(iy),epsilo) |
---|
653 | voz2(ix,iy,iz)= voz2(ix,iy,iz)/k2 |
---|
654 | end do |
---|
655 | end do |
---|
656 | end do |
---|
657 | |
---|
658 | ! vitesses en fct de la fct de courant 2D*1D |
---|
659 | do iz=0,nzmax-1 |
---|
660 | do iy=0,nymax-1 |
---|
661 | do ix=0,kxmax |
---|
662 | vox2(2*ix,iy,iz) = -ky(iy)*voz2(2*ix+1,iy,iz) |
---|
663 | vox2(2*ix+1,iy,iz) = ky(iy)*voz2(2*ix,iy,iz) |
---|
664 | voy2(2*ix,iy,iz) = kx(ix)*voz2(2*ix+1,iy,iz) |
---|
665 | voy2(2*ix+1,iy,iz) = -kx(ix)*voz2(2*ix,iy,iz) |
---|
666 | end do |
---|
667 | end do |
---|
668 | end do |
---|
669 | |
---|
670 | ! passage dans l'espace physique |
---|
671 | call bck_fft(vox2,fft_work) |
---|
672 | call bck_fft(voy2,fft_work) |
---|
673 | |
---|
674 | vbx = vox2(0:nxmax-1,0:nymax-1,0) |
---|
675 | vby = voy2(0:nxmax-1,0:nymax-1,0) |
---|
676 | |
---|
677 | end subroutine |
---|
678 | #endif |
---|
679 | |
---|
680 | |
---|
681 | ! ******************************************************************************** |
---|
682 | subroutine read_base2D(vbx,vby,vbz,wbx,wby,wbz) |
---|
683 | ! ******************************************************************************** |
---|
684 | ! Subroutine qui permet de lire un champ de base en formatted ou unformatted |
---|
685 | ! ******************************************************************************** |
---|
686 | |
---|
687 | implicit none |
---|
688 | |
---|
689 | double precision, dimension(0:nxmax-1,0:nymax-1), intent(out) :: vbx, vby, vbz, wbx, wby, wbz |
---|
690 | integer :: flag, nxread, nyread |
---|
691 | double precision :: lxread, lyread |
---|
692 | |
---|
693 | open (unit=88,file='base2D.init',form='unformatted') |
---|
694 | read(88) flag, nxread, nyread, lxread, lyread |
---|
695 | read(88) vbx, vby, vbz |
---|
696 | read(88) wbx, wby, wbz |
---|
697 | close(88) |
---|
698 | |
---|
699 | ! on verifie que le format du fichier d'entree est correct |
---|
700 | if (flag /= 1 .or. nxread /= nxmax .or. nyread /= nymax .or. lxread /= lx .or. lyread /= ly) then |
---|
701 | write(*,*) "ERREUR: le fichier de reprise de run n'a pas le bon format" |
---|
702 | write(*,*) "ARRET DU PROGRAMME..." |
---|
703 | stop |
---|
704 | end if |
---|
705 | |
---|
706 | end subroutine |
---|
707 | |
---|
708 | |
---|
709 | |
---|
710 | ! ******************************************************************************** |
---|
711 | subroutine read_velo2D(vox, voy, voz) |
---|
712 | ! ******************************************************************************** |
---|
713 | ! Subroutine qui permet de lire un champ 2D en unformatted |
---|
714 | ! ******************************************************************************** |
---|
715 | |
---|
716 | implicit none |
---|
717 | |
---|
718 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: vox, voy, voz |
---|
719 | double precision, dimension(0:nxmax-1,0:nymax-1) :: vi2Dx, vi2Dy, vi2Dz |
---|
720 | integer :: flag, nxread, nyread |
---|
721 | double precision :: lxread, lyread |
---|
722 | integer :: ix,iy,iz |
---|
723 | |
---|
724 | open (unit=88,file='velo2D.init',form='unformatted') |
---|
725 | read(88) flag, nxread, nyread, lxread, lyread |
---|
726 | read(88) vi2Dx,vi2Dy,vi2Dz |
---|
727 | close(88) |
---|
728 | |
---|
729 | ! on verifie que le format du fichier d'entree est correct |
---|
730 | if (flag /= 1 .or. nxread /= nxmax .or. nyread /= nymax .or. lxread /= lx .or. lyread /= ly) then |
---|
731 | write(*,*) "ERREUR: le fichier de reprise de run n'a pas le bon format" |
---|
732 | write(*,*) "ARRET DU PROGRAMME..." |
---|
733 | stop |
---|
734 | end if |
---|
735 | |
---|
736 | do iz=0,nzmax-1 |
---|
737 | do iy=0,nymax-1 |
---|
738 | do ix=0,nxmax-1 |
---|
739 | vox(ix,iy,iz)=vi2Dx(ix,iy) |
---|
740 | voy(ix,iy,iz)=vi2Dy(ix,iy) |
---|
741 | voz(ix,iy,iz)=vi2Dz(ix,iy) |
---|
742 | end do |
---|
743 | end do |
---|
744 | end do |
---|
745 | |
---|
746 | end subroutine |
---|
747 | |
---|
748 | |
---|
749 | ! ******************************************************************************** |
---|
750 | subroutine readRho(Rhov) |
---|
751 | !ajout d'une sous routine de lecture d'un gradient de densite moyen non constant |
---|
752 | ! ******************************************************************************** |
---|
753 | ! Subroutine qui permet de lire un champ 2D en unformatted |
---|
754 | ! ******************************************************************************** |
---|
755 | |
---|
756 | implicit none |
---|
757 | |
---|
758 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: Rhov |
---|
759 | double precision, dimension(0:nymax-1) :: Rhoi |
---|
760 | integer :: flag, nxread, nyread |
---|
761 | double precision :: lxread, lyread |
---|
762 | integer :: ix,iy,iz |
---|
763 | write(*,*) 'boucle lecture rho' |
---|
764 | open (unit=88,file='Rho.init',status='old') |
---|
765 | !read(88) flag, nxread, nyread, lxread, lyread |
---|
766 | do iy=0,255 |
---|
767 | read(88,*) Rhoi(iy) |
---|
768 | enddo |
---|
769 | |
---|
770 | close(88) |
---|
771 | |
---|
772 | ! on verifie que le format du fichier d'entree est correct |
---|
773 | ! if (flag /= 1 .or. nyread /= nymax .or. lyread /= ly) then |
---|
774 | ! write(*,*) "ERREUR: le fichier de reprise de run n'a pas le bon format" |
---|
775 | ! write(*,*) "ARRET DU PROGRAMME..." |
---|
776 | ! stop |
---|
777 | !end if |
---|
778 | |
---|
779 | do iz=0,nzmax-1 |
---|
780 | do iy=0,nymax-1 |
---|
781 | do ix=0,nxmax-1 |
---|
782 | Rhov(ix,iy,iz)=Rhoi(iy) |
---|
783 | end do |
---|
784 | end do |
---|
785 | end do |
---|
786 | print*, Rhov |
---|
787 | |
---|
788 | end subroutine |
---|
789 | |
---|
790 | |
---|
791 | |
---|
792 | |
---|
793 | ! ************************************************************************************ |
---|
794 | subroutine gen_velo2D(vox, voy, voz) |
---|
795 | ! ************************************************************************************ |
---|
796 | ! Generation d'un champ de vitesse correspondant a un dipole fait de |
---|
797 | ! tourbillons de Lamb-oseen. |
---|
798 | ! ************************************************************************************ |
---|
799 | |
---|
800 | implicit none |
---|
801 | |
---|
802 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: vox,voy,voz |
---|
803 | double precision k2 |
---|
804 | integer num_tourbillon |
---|
805 | integer ix, iy, iz |
---|
806 | |
---|
807 | vox = 0. |
---|
808 | voy = 0. |
---|
809 | voz = 0. |
---|
810 | |
---|
811 | ! on determine la vorticite verticale de l'ecoulement |
---|
812 | do num_tourbillon=1,nb_tourbillons |
---|
813 | do iz=0,nzmax-1 |
---|
814 | do iy=0,nymax-1 |
---|
815 | do ix=0,nxmax-1 |
---|
816 | voz(ix,iy,iz)=voz(ix,iy,iz)+wz_gau((xx(ix)-xt(num_tourbillon)),(yy(iy)-yt(num_tourbillon)), & |
---|
817 | gamma(num_tourbillon),at(num_tourbillon)) |
---|
818 | end do |
---|
819 | end do |
---|
820 | end do |
---|
821 | end do |
---|
822 | |
---|
823 | ! inversion du laplacien 2D pour trouver la fct de courant 2D*1D |
---|
824 | call fwd_fft(voz,fft_work) |
---|
825 | |
---|
826 | do iz=0,nzmax-1 |
---|
827 | do iy=0,nymax-1 |
---|
828 | do ix=0,2*kxmax+1 |
---|
829 | k2=max(ca2x(ix)+k2y(iy),epsilo) |
---|
830 | voz(ix,iy,iz)= voz(ix,iy,iz)/k2 |
---|
831 | end do |
---|
832 | end do |
---|
833 | end do |
---|
834 | |
---|
835 | ! vitesses en fct de la fct de courant 2D*1D |
---|
836 | do iz=0,nzmax-1 |
---|
837 | do iy=0,nymax-1 |
---|
838 | do ix=0,kxmax |
---|
839 | vox(2*ix,iy,iz) = -ky(iy)*voz(2*ix+1,iy,iz) |
---|
840 | vox(2*ix+1,iy,iz) = ky(iy)*voz(2*ix,iy,iz) |
---|
841 | voy(2*ix,iy,iz) = kx(ix)*voz(2*ix+1,iy,iz) |
---|
842 | voy(2*ix+1,iy,iz) = -kx(ix)*voz(2*ix,iy,iz) |
---|
843 | end do |
---|
844 | end do |
---|
845 | end do |
---|
846 | |
---|
847 | ! passage dans l'espace physique |
---|
848 | call bck_fft(vox,fft_work) |
---|
849 | call bck_fft(voy,fft_work) |
---|
850 | voz = 0 |
---|
851 | |
---|
852 | end subroutine |
---|
853 | |
---|
854 | ! *************************************************************** |
---|
855 | function wz_gau(x,y,gamma,a) |
---|
856 | ! *************************************************************** |
---|
857 | |
---|
858 | implicit none |
---|
859 | |
---|
860 | double precision, intent(in) :: x, y, gamma,a |
---|
861 | double precision :: wz_gau |
---|
862 | double precision :: r2 |
---|
863 | |
---|
864 | r2 = x**2 + y**2 |
---|
865 | wz_gau =gamma/(pi*a**2)*exp(-r2/a**2) |
---|
866 | |
---|
867 | end function |
---|
868 | |
---|
869 | |
---|
870 | ! ************************************************************************* |
---|
871 | subroutine white_noise(vox,voy,voz) |
---|
872 | ! ********************************************************* |
---|
873 | |
---|
874 | implicit none |
---|
875 | |
---|
876 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz |
---|
877 | double precision :: rnd_number, moy_x, moy_y, moy_z |
---|
878 | integer :: ix, iy, iz |
---|
879 | |
---|
880 | ! on initialise le generateur de nombre aleatoire avec une nouvelle "graine" |
---|
881 | ! ne semble pas fontionner sur toutes les architectures (bizarre, bizarre...) |
---|
882 | |
---|
883 | ! on desactive l'initialisation du generateur de nombre aleatoire en mode debugage de façon |
---|
884 | ! a pouvoir comparer plus facilement les resultats de runs differents. |
---|
885 | #ifndef debug |
---|
886 | call random_seed |
---|
887 | #endif |
---|
888 | |
---|
889 | ! on determine la moyenne initiale des champs vox, voy, voz |
---|
890 | moy_x = sum(vox(0:nxmax-1,0:nymax-1,0:nzmax-1))/(nxmax*nymax*nzmax) |
---|
891 | moy_y = sum(voy(0:nxmax-1,0:nymax-1,0:nzmax-1))/(nxmax*nymax*nzmax) |
---|
892 | moy_z = sum(voz(0:nxmax-1,0:nymax-1,0:nzmax-1))/(nxmax*nymax*nzmax) |
---|
893 | |
---|
894 | ! on ajoute le bruit blanc dans l'espace physique |
---|
895 | do iz=0,nzmax-1 |
---|
896 | do iy=0,nymax-1 |
---|
897 | do ix=0,nxmax-1 |
---|
898 | call random_number(rnd_number) |
---|
899 | vox(ix,iy,iz)=vox(ix,iy,iz)+amplitude_noise*rnd_number |
---|
900 | call random_number(rnd_number) |
---|
901 | voy(ix,iy,iz)=voy(ix,iy,iz)+amplitude_noise*rnd_number |
---|
902 | call random_number(rnd_number) |
---|
903 | voz(ix,iy,iz)=voz(ix,iy,iz)+amplitude_noise*rnd_number |
---|
904 | end do |
---|
905 | end do |
---|
906 | end do |
---|
907 | |
---|
908 | ! on fait en sorte que la moyenne du bruit que l'on a ajouté soit nulle |
---|
909 | ! (on utilise pour cela la difference des moyennes des champs vo avant et après ajout du bruit...) |
---|
910 | moy_x = moy_x - sum(vox(0:nxmax-1,0:nymax-1,0:nzmax-1))/(nxmax*nymax*nzmax) |
---|
911 | moy_y = moy_y - sum(voy(0:nxmax-1,0:nymax-1,0:nzmax-1))/(nxmax*nymax*nzmax) |
---|
912 | moy_z = moy_z - sum(voz(0:nxmax-1,0:nymax-1,0:nzmax-1))/(nxmax*nymax*nzmax) |
---|
913 | |
---|
914 | vox = vox + moy_x |
---|
915 | voy = voy + moy_y |
---|
916 | voz = voz + moy_z |
---|
917 | |
---|
918 | end subroutine |
---|
919 | |
---|
920 | ! ************************************************************************* |
---|
921 | subroutine time2string(time,string_time) |
---|
922 | ! ************************************************************************* |
---|
923 | ! CONVERSION TIME TO STRING t=0000.000 |
---|
924 | ! ************************************************************************* |
---|
925 | |
---|
926 | implicit none |
---|
927 | |
---|
928 | double precision, intent(in) :: time |
---|
929 | character(len=13), intent(out) :: string_time |
---|
930 | |
---|
931 | write(string_time,'(A2,I7.7,A1,I3.3)') 't=',int(time),'.',nint(1000*(time-int(time))) |
---|
932 | end subroutine |
---|
933 | |
---|
934 | |
---|
935 | !********************************************************************************** |
---|
936 | !********************************************************************************** |
---|
937 | !********** FIN DES ENTREE SORTIES ***************** |
---|
938 | !********************************************************************************** |
---|
939 | !********************************************************************************** |
---|
940 | |
---|
941 | ! ********************************************************************************* |
---|
942 | subroutine adams_bashforth(vox,voy,voz,rho,fx,fy,fz,frho,vox2,voy2,voz2,rho2) |
---|
943 | ! ********************************************************************************* |
---|
944 | |
---|
945 | ! coeur du programme avance temporelle sur itmax |
---|
946 | ! Time discretisation: second order Adams-Bashforth: |
---|
947 | ! step t+dt: input fields are : V(k,t) and F(k,t-dt) |
---|
948 | ! output fields are: V(k,t+dt) and F(k,t) |
---|
949 | |
---|
950 | implicit none |
---|
951 | |
---|
952 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho |
---|
953 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: fx, fy, fz, frho |
---|
954 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox2, voy2, voz2, rho2 |
---|
955 | |
---|
956 | call cpu_time(tps_temp) |
---|
957 | tps_copie = tps_copie - tps_temp |
---|
958 | vox2 = vox |
---|
959 | voy2 = voy |
---|
960 | voz2 = voz |
---|
961 | if (stratification) then |
---|
962 | rho2 = rho |
---|
963 | end if |
---|
964 | call cpu_time(tps_temp) |
---|
965 | tps_copie = tps_copie + tps_temp |
---|
966 | |
---|
967 | |
---|
968 | ! les champs contiennent vo = V(t) / vo2 = V(t) / f = f(V(t-dt)) |
---|
969 | |
---|
970 | ! ************ |
---|
971 | ! Etape 1 / V(t*,k)=E*V(t,k)-(dt/2)*E**2*F(t-dt,k) |
---|
972 | ! avec E=exp(-nu*k2*dt) des constantes precalculees |
---|
973 | ! ************ |
---|
974 | call cpu_time(tps_temp) |
---|
975 | tps_adams_step1 = tps_adams_step1 - tps_temp |
---|
976 | call adams_step1(vox,voy,voz,rho,fx,fy,fz,frho) |
---|
977 | call cpu_time(tps_temp) |
---|
978 | tps_adams_step1 = tps_adams_step1 + tps_temp |
---|
979 | |
---|
980 | ! les champs contiennent vo = V(t*) / vo2 = V(t) / f = f(V(t-dt)) |
---|
981 | |
---|
982 | ! ************ |
---|
983 | ! Etape 2 / V(t+dt) = V(t*)+(3.*dt/2.)*E*F(V(t)) |
---|
984 | ! avec E=exp(-nu*k2*dt) des constantes precalculees |
---|
985 | ! ************ |
---|
986 | |
---|
987 | call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho,Rhov) |
---|
988 | |
---|
989 | call cpu_time(tps_temp) |
---|
990 | tps_adams_step2 = tps_adams_step2 - tps_temp |
---|
991 | call adams_step2(vox,voy,voz,rho,fx,fy,fz,frho) |
---|
992 | call cpu_time(tps_temp) |
---|
993 | tps_adams_step2 = tps_adams_step2 + tps_temp |
---|
994 | |
---|
995 | ! les champs contiennent vo = V(t+dt) / vo2 = ## / f = f(V(t)) |
---|
996 | |
---|
997 | end subroutine |
---|
998 | |
---|
999 | |
---|
1000 | ! ************************************************************************* |
---|
1001 | subroutine adams_step1(vox,voy,voz,rho,fx,fy,fz,frho) |
---|
1002 | ! ************************************************************************* |
---|
1003 | ! Première étape du schéma d'ordre deux d'Adams-Bashforth |
---|
1004 | ! ************************************************************************* |
---|
1005 | |
---|
1006 | |
---|
1007 | implicit none |
---|
1008 | |
---|
1009 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho |
---|
1010 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx, fy, fz, frho |
---|
1011 | double precision :: ee,ees,e2,e2s |
---|
1012 | double precision :: temp_dt |
---|
1013 | integer :: ix, iy, iz |
---|
1014 | |
---|
1015 | temp_dt=dt/2. |
---|
1016 | |
---|
1017 | |
---|
1018 | do iz=0,nzmax-1 |
---|
1019 | do iy=0,nymax-1 |
---|
1020 | do ix=0,2*kxmax+1 |
---|
1021 | #ifdef test_troncature |
---|
1022 | if (trunctab(ix,iy,iz)) then |
---|
1023 | #endif |
---|
1024 | ee=eex(ix)*eey(iy)*eez(iz) |
---|
1025 | ees=exs(ix)*eeys(iy)*eezs(iz) |
---|
1026 | e2=ee*ee |
---|
1027 | e2s=ees*ees |
---|
1028 | vox(ix,iy,iz)=ee*vox(ix,iy,iz)-temp_dt*e2*fx(ix,iy,iz) |
---|
1029 | voy(ix,iy,iz)=ee*voy(ix,iy,iz)-temp_dt*e2*fy(ix,iy,iz) |
---|
1030 | voz(ix,iy,iz)=ee*voz(ix,iy,iz)-temp_dt*e2*fz(ix,iy,iz) |
---|
1031 | rho(ix,iy,iz)=ees*rho(ix,iy,iz)-temp_dt*e2s*frho(ix,iy,iz) |
---|
1032 | #ifdef test_troncature |
---|
1033 | end if |
---|
1034 | #endif |
---|
1035 | end do |
---|
1036 | end do |
---|
1037 | end do |
---|
1038 | |
---|
1039 | |
---|
1040 | end subroutine |
---|
1041 | |
---|
1042 | |
---|
1043 | ! ************************************************************************* |
---|
1044 | subroutine adams_step2(vox,voy,voz,rho,fx,fy,fz,frho) |
---|
1045 | ! ************************************************************************* |
---|
1046 | ! Deuxième étape du schéma d'ordre deux d'Adams-Bashforth |
---|
1047 | ! ************************************************************************* |
---|
1048 | |
---|
1049 | implicit none |
---|
1050 | |
---|
1051 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho |
---|
1052 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx, fy, fz, frho |
---|
1053 | double precision :: ee,ees |
---|
1054 | double precision :: temp_dt |
---|
1055 | integer :: ix, iy, iz |
---|
1056 | |
---|
1057 | temp_dt=3.*dt/2. |
---|
1058 | |
---|
1059 | |
---|
1060 | do iz=0,nzmax-1 |
---|
1061 | do iy=0,nymax-1 |
---|
1062 | do ix=0,2*kxmax+1 |
---|
1063 | #ifdef test_troncature |
---|
1064 | if (trunctab(ix,iy,iz)) then |
---|
1065 | #endif |
---|
1066 | ee=eex(ix)*eey(iy)*eez(iz) |
---|
1067 | ees=exs(ix)*eeys(iy)*eezs(iz) |
---|
1068 | vox(ix,iy,iz)=vox(ix,iy,iz)+ee*temp_dt*fx(ix,iy,iz) |
---|
1069 | voy(ix,iy,iz)=voy(ix,iy,iz)+ee*temp_dt*fy(ix,iy,iz) |
---|
1070 | voz(ix,iy,iz)=voz(ix,iy,iz)+ee*temp_dt*fz(ix,iy,iz) |
---|
1071 | rho(ix,iy,iz)=rho(ix,iy,iz)+ees*temp_dt*frho(ix,iy,iz) |
---|
1072 | #ifdef test_troncature |
---|
1073 | end if |
---|
1074 | #endif |
---|
1075 | end do |
---|
1076 | end do |
---|
1077 | end do |
---|
1078 | |
---|
1079 | |
---|
1080 | end subroutine |
---|
1081 | |
---|
1082 | |
---|
1083 | |
---|
1084 | ! ************************************************************************* |
---|
1085 | subroutine runge_step2(vox,voy,voz,rho,fx,fy,fz,frho) |
---|
1086 | ! ************************************************************************* |
---|
1087 | ! schema de type leap-frog centre a t+dt/2. |
---|
1088 | ! Utilise comme etape du subroutine raccro |
---|
1089 | ! qui effectue un pas Runge-Kutta a l'ordre deux pour redemarrer. |
---|
1090 | ! ( v(n+1/2)+1/2 - v(n+1/2)-1/2 )/(2.dt/2) = f( v(n+1/2) ) . |
---|
1091 | ! ************************************************************************** |
---|
1092 | |
---|
1093 | implicit none |
---|
1094 | |
---|
1095 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho |
---|
1096 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx, fy, fz, frho |
---|
1097 | double precision :: k2,ee,a,ees,as |
---|
1098 | integer :: ix, iy, iz |
---|
1099 | |
---|
1100 | do iz=0,nzmax-1 |
---|
1101 | do iy=0,nymax-1 |
---|
1102 | do ix=0,2*kxmax+1 |
---|
1103 | k2=max(ca2x(ix)+k2y(iy)+k2z(iz),epsilo) |
---|
1104 | ee=eex(ix)*eey(iy)*eez(iz) |
---|
1105 | a=(1.-ee)/(nu*k2) |
---|
1106 | ees=exs(ix)*eeys(iy)*eezs(iz) |
---|
1107 | as=(1.-ees)/(nu/schm*k2) |
---|
1108 | |
---|
1109 | vox(ix,iy,iz)=vox(ix,iy,iz)*ee+fx(ix,iy,iz)*a |
---|
1110 | voy(ix,iy,iz)=voy(ix,iy,iz)*ee+fy(ix,iy,iz)*a |
---|
1111 | voz(ix,iy,iz)=voz(ix,iy,iz)*ee+fz(ix,iy,iz)*a |
---|
1112 | rho(ix,iy,iz)=rho(ix,iy,iz)*ees+frho(ix,iy,iz)*as |
---|
1113 | end do |
---|
1114 | end do |
---|
1115 | end do |
---|
1116 | |
---|
1117 | end subroutine |
---|
1118 | |
---|
1119 | |
---|
1120 | ! ************************************************************************* |
---|
1121 | subroutine runge_step1(vox,voy,voz,rho,fx,fy,fz,frho) |
---|
1122 | ! ************************************************************************* |
---|
1123 | ! schema de type Euler. |
---|
1124 | ! Utilisé comme étape du subroutine raccro |
---|
1125 | ! qui effectue un pas Runge-Kutta a l'ordre deux pour redemarrer. |
---|
1126 | ! ( v(n+1/2)-v(n))/(dt/2) = f( v(n) ) |
---|
1127 | ! ************************************************************************* |
---|
1128 | |
---|
1129 | implicit none |
---|
1130 | |
---|
1131 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho |
---|
1132 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx, fy, fz, frho |
---|
1133 | double precision :: k2,ee,a,ees,as |
---|
1134 | integer :: ix, iy, iz |
---|
1135 | |
---|
1136 | do iz=0,nzmax-1 |
---|
1137 | do iy=0,nymax-1 |
---|
1138 | do ix=0,2*kxmax+1 |
---|
1139 | k2=max(ca2x(ix)+k2y(iy)+k2z(iz),epsilo) |
---|
1140 | ee=eex(ix)*eey(iy)*eez(iz) |
---|
1141 | a=(1.-ee)/(nu*k2) |
---|
1142 | ees=exs(ix)*eeys(iy)*eezs(iz) |
---|
1143 | as=(1.-ees)/(nu/schm*k2) |
---|
1144 | |
---|
1145 | vox(ix,iy,iz)=vox(ix,iy,iz)*ee+fx(ix,iy,iz)*a |
---|
1146 | voy(ix,iy,iz)=voy(ix,iy,iz)*ee+fy(ix,iy,iz)*a |
---|
1147 | voz(ix,iy,iz)=voz(ix,iy,iz)*ee+fz(ix,iy,iz)*a |
---|
1148 | rho(ix,iy,iz)=rho(ix,iy,iz)*ees+frho(ix,iy,iz)*as |
---|
1149 | end do |
---|
1150 | end do |
---|
1151 | end do |
---|
1152 | |
---|
1153 | end subroutine |
---|
1154 | |
---|
1155 | ! ************************************************************************* |
---|
1156 | subroutine rotv(ax,ay,az,bx,by,bz) |
---|
1157 | ! ************************************************************************* |
---|
1158 | ! calcul du rotationnel du champ (ax,ay,az). |
---|
1159 | ! ************************************************************************* |
---|
1160 | |
---|
1161 | implicit none |
---|
1162 | |
---|
1163 | double precision, dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: ax, ay, az |
---|
1164 | double precision, dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: bx, by, bz |
---|
1165 | integer :: ix, iy, iz |
---|
1166 | |
---|
1167 | |
---|
1168 | do iz=0,nzmax-1 |
---|
1169 | do iy=0,nymax-1 |
---|
1170 | do ix=0,kxmax |
---|
1171 | #ifndef mode_kz |
---|
1172 | bx(1,ix,iy,iz)=kz(iz)*ay(2,ix,iy,iz)-ky(iy)*az(2,ix,iy,iz) |
---|
1173 | bx(2,ix,iy,iz)=ky(iy)*az(1,ix,iy,iz)-kz(iz)*ay(1,ix,iy,iz) |
---|
1174 | by(1,ix,iy,iz)=kx(ix)*az(2,ix,iy,iz)-kz(iz)*ax(2,ix,iy,iz) |
---|
1175 | by(2,ix,iy,iz)=kz(iz)*ax(1,ix,iy,iz)-kx(ix)*az(1,ix,iy,iz) |
---|
1176 | #else |
---|
1177 | bx(1,ix,iy,iz)=-kz(iz)*ay(1,ix,iy,iz)-ky(iy)*az(2,ix,iy,iz) |
---|
1178 | bx(2,ix,iy,iz)=ky(iy)*az(1,ix,iy,iz)-kz(iz)*ay(2,ix,iy,iz) |
---|
1179 | by(1,ix,iy,iz)=kx(ix)*az(2,ix,iy,iz)+kz(iz)*ax(1,ix,iy,iz) |
---|
1180 | by(2,ix,iy,iz)=kz(iz)*ax(2,ix,iy,iz)-kx(ix)*az(1,ix,iy,iz) |
---|
1181 | #endif |
---|
1182 | |
---|
1183 | bz(1,ix,iy,iz)=ky(iy)*ax(2,ix,iy,iz)-kx(ix)*ay(2,ix,iy,iz) |
---|
1184 | bz(2,ix,iy,iz)=kx(ix)*ay(1,ix,iy,iz)-ky(iy)*ax(1,ix,iy,iz) |
---|
1185 | end do |
---|
1186 | end do |
---|
1187 | end do |
---|
1188 | |
---|
1189 | end subroutine |
---|
1190 | |
---|
1191 | |
---|
1192 | ! ************************************************************************* |
---|
1193 | subroutine trunc0(ax,ay,az,b) |
---|
1194 | ! ************************************************************************* |
---|
1195 | ! classical desaliasing in the spectral space. |
---|
1196 | ! ************************************************************************** |
---|
1197 | |
---|
1198 | implicit none |
---|
1199 | |
---|
1200 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: ax, ay, az, b |
---|
1201 | integer :: ix, iy, iz |
---|
1202 | |
---|
1203 | do iz=0,nzmax-1 |
---|
1204 | do iy=0,nymax-1 |
---|
1205 | do ix=0,2*kxmax+1 |
---|
1206 | if(.not. trunctab(ix,iy,iz)) then |
---|
1207 | ax(ix,iy,iz) = 0 |
---|
1208 | ay(ix,iy,iz) = 0 |
---|
1209 | az(ix,iy,iz) = 0 |
---|
1210 | b(ix,iy,iz) = 0 |
---|
1211 | end if |
---|
1212 | end do |
---|
1213 | end do |
---|
1214 | end do |
---|
1215 | |
---|
1216 | end subroutine |
---|
1217 | |
---|
1218 | ! ************************************************************************* |
---|
1219 | subroutine init_trunctab(trunctab) |
---|
1220 | ! ************************************************************************* |
---|
1221 | ! Génération du tableau utilisé lors de la troncature |
---|
1222 | ! ************************************************************************** |
---|
1223 | |
---|
1224 | implicit none |
---|
1225 | |
---|
1226 | logical, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: trunctab |
---|
1227 | double precision :: kx2max,ky2max,kz2max,dks |
---|
1228 | integer :: ix, iy, iz |
---|
1229 | |
---|
1230 | ! rayon du cercle: kmax*rtrunc |
---|
1231 | kx2max=(rtrunc_x*nxmax/2*2*pi/lx)**2 |
---|
1232 | ky2max=(rtrunc_y*nymax/2*2*pi/ly)**2 |
---|
1233 | |
---|
1234 | #ifndef mode_kz |
---|
1235 | kz2max=(rtrunc_z*nzmax/2*2*pi/lz)**2 |
---|
1236 | #else |
---|
1237 | kz2max=rtrunc_z**2 |
---|
1238 | #endif |
---|
1239 | |
---|
1240 | trunctab = .false. |
---|
1241 | do iz=0,nzmax-1 |
---|
1242 | do iy=0,nymax-1 |
---|
1243 | do ix=0,2*kxmax+1 |
---|
1244 | |
---|
1245 | if (type_trunc==1) then |
---|
1246 | dks = max(ca2x(ix)/kx2max, k2y(iy)/ky2max, k2z(iz)/kz2max) |
---|
1247 | else |
---|
1248 | dks = ca2x(ix)/kx2max + k2y(iy)/ky2max + k2z(iz)/kz2max |
---|
1249 | end if |
---|
1250 | |
---|
1251 | ! on rajoute un epsilo pour eviter que numériquement, un terme soit tronqué alors qu'il |
---|
1252 | ! ne devrait pas l'être. |
---|
1253 | if (dks <= 1. + epsilo) then |
---|
1254 | trunctab(ix,iy,iz)=.true. |
---|
1255 | end if |
---|
1256 | end do |
---|
1257 | end do |
---|
1258 | end do |
---|
1259 | |
---|
1260 | end subroutine |
---|
1261 | |
---|
1262 | |
---|
1263 | ! ************************************************************************* |
---|
1264 | subroutine vecpro2(ax,ay,az,bx,by,bz) |
---|
1265 | ! ************************************************************************* |
---|
1266 | ! produits vectoriels dans l'espace physique pour le code perturbatif. |
---|
1267 | ! inclut les termes avec l'ecoulement de base et coriolis |
---|
1268 | ! ************************************************************************** |
---|
1269 | |
---|
1270 | implicit none |
---|
1271 | |
---|
1272 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bx, by, bz |
---|
1273 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: ax, ay, az |
---|
1274 | double precision :: vbx0,vby0,vbz0,wbx0,wby0,wbz0 |
---|
1275 | double precision :: ax0,ay0,az0,bx0,by0,bz0 |
---|
1276 | integer :: ix, iy, iz |
---|
1277 | |
---|
1278 | ! VERSION LINEARISEE |
---|
1279 | if (lin) then |
---|
1280 | |
---|
1281 | do iz=0,nzmax-1 |
---|
1282 | do iy=0,nymax-1 |
---|
1283 | do ix=0,nxmax-1 |
---|
1284 | |
---|
1285 | vbx0=vbx(ix,iy) |
---|
1286 | vby0=vby(ix,iy) |
---|
1287 | vbz0=vbz(ix,iy) |
---|
1288 | wbx0=wbx(ix,iy) |
---|
1289 | wby0=wby(ix,iy) |
---|
1290 | wbz0=wbz(ix,iy)+omega2 |
---|
1291 | |
---|
1292 | ax0=ax(ix,iy,iz) |
---|
1293 | ay0=ay(ix,iy,iz) |
---|
1294 | az0=az(ix,iy,iz) |
---|
1295 | bx0=bx(ix,iy,iz) |
---|
1296 | by0=by(ix,iy,iz) |
---|
1297 | bz0=bz(ix,iy,iz) |
---|
1298 | |
---|
1299 | bx(ix,iy,iz)=vby0*bz0-vbz0*by0 + ay0*wbz0-az0*wby0 |
---|
1300 | by(ix,iy,iz)=vbz0*bx0-vbx0*bz0 + az0*wbx0-ax0*wbz0 |
---|
1301 | bz(ix,iy,iz)=vbx0*by0-vby0*bx0 + ax0*wby0-ay0*wbx0 |
---|
1302 | end do |
---|
1303 | end do |
---|
1304 | end do |
---|
1305 | |
---|
1306 | ! NON LINEAIRE |
---|
1307 | else |
---|
1308 | ! calcul en nonlineaire perturbatif : |
---|
1309 | if (perturbe) then |
---|
1310 | |
---|
1311 | do iz=0,nzmax-1 |
---|
1312 | do iy=0,nymax-1 |
---|
1313 | do ix=0,nxmax-1 |
---|
1314 | vbx0=vbx(ix,iy) |
---|
1315 | vby0=vby(ix,iy) |
---|
1316 | vbz0=vbz(ix,iy) |
---|
1317 | wbx0=wbx(ix,iy) |
---|
1318 | wby0=wby(ix,iy) |
---|
1319 | wbz0=wbz(ix,iy)+omega2 |
---|
1320 | |
---|
1321 | ax0=ax(ix,iy,iz) |
---|
1322 | ay0=ay(ix,iy,iz) |
---|
1323 | az0=az(ix,iy,iz) |
---|
1324 | bx0=bx(ix,iy,iz) |
---|
1325 | by0=by(ix,iy,iz) |
---|
1326 | bz0=bz(ix,iy,iz) |
---|
1327 | |
---|
1328 | bx(ix,iy,iz)=vby0*bz0-vbz0*by0 + ay0*(wbz0+bz0) - az0*(wby0+by0) |
---|
1329 | by(ix,iy,iz)=vbz0*bx0-vbx0*bz0 + az0*(wbx0+bx0) - ax0*(wbz0+bz0) |
---|
1330 | bz(ix,iy,iz)=vbx0*by0-vby0*bx0 + ax0*(wby0+by0) - ay0*(wbx0+bx0) |
---|
1331 | end do |
---|
1332 | end do |
---|
1333 | end do |
---|
1334 | |
---|
1335 | ! calcul en non perturbatif : |
---|
1336 | else |
---|
1337 | |
---|
1338 | do iz=0,nzmax-1 |
---|
1339 | do iy=0,nymax-1 |
---|
1340 | do ix=0,nxmax-1 |
---|
1341 | ax0=ax(ix,iy,iz) |
---|
1342 | ay0=ay(ix,iy,iz) |
---|
1343 | az0=az(ix,iy,iz) |
---|
1344 | bx0=bx(ix,iy,iz) |
---|
1345 | by0=by(ix,iy,iz) |
---|
1346 | bz0=bz(ix,iy,iz) |
---|
1347 | |
---|
1348 | bx(ix,iy,iz)=ay0*(omega2+bz0)-az0*by0 |
---|
1349 | by(ix,iy,iz)=az0*bx0-ax0*(omega2+bz0) |
---|
1350 | bz(ix,iy,iz)=ax0*by0-ay0*bx0 |
---|
1351 | end do |
---|
1352 | end do |
---|
1353 | end do |
---|
1354 | |
---|
1355 | end if |
---|
1356 | end if |
---|
1357 | |
---|
1358 | end subroutine |
---|
1359 | |
---|
1360 | ! pba: deb |
---|
1361 | ! ************************************************************************* |
---|
1362 | subroutine forcage(bx,by,bz) |
---|
1363 | ! ************************************************************************* |
---|
1364 | ! ajout d'un terme de forcage |
---|
1365 | ! parametres du forcage entres ici provisoirement modifier |
---|
1366 | ! data.in etc pour entree plus propre |
---|
1367 | ! ************************************************************************** |
---|
1368 | |
---|
1369 | |
---|
1370 | implicit none |
---|
1371 | |
---|
1372 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bx, by, bz |
---|
1373 | double precision :: vbx0,vby0,vbz0,wbx0,wby0,wbz0 |
---|
1374 | double precision :: ax0,ay0,az0,bx0,by0,bz0 |
---|
1375 | double precision :: time, arg, sigma1, sigma2, FF, GG, acc |
---|
1376 | integer :: ix, iy, iz, itime, i0, j0, k0 !itime a faire passer en parametre? |
---|
1377 | ! itime=it0 |
---|
1378 | |
---|
1379 | ! itime=itime+1 ! penser a it0=it0+itmax(run precedent si restat |
---|
1380 | ! print *,'temps:', itime |
---|
1381 | ! A_ics=100. !deplacement a definir dans data.in et parametres.F90 |
---|
1382 | !freq=2*pi/(3600*12.4) !M2 a verifier idem |
---|
1383 | arg=freq*it*Dt !Dt est defini dans data.in ok |
---|
1384 | ! print *,'arg:', arg |
---|
1385 | ! print *,'A_ics:', A_ics !acceleration |
---|
1386 | |
---|
1387 | |
---|
1388 | ! localisation centre de la source |
---|
1389 | |
---|
1390 | i0=floor(0.5*nxmax) |
---|
1391 | j0=floor(0.5*nymax) |
---|
1392 | k0=floor(0.5*nzmax) |
---|
1393 | |
---|
1394 | ! enveloppe gaussienne |
---|
1395 | |
---|
1396 | !sigma1=4km, sigma2=100m revoir les valeurs ci-dessous |
---|
1397 | sigma1=20. ! defini en unite de dx et dz |
---|
1398 | sigma2=1. ! dy |
---|
1399 | |
---|
1400 | |
---|
1401 | do iz=0,nzmax-1 |
---|
1402 | do iy=0,nymax-1 |
---|
1403 | do ix=0,nxmax-1 |
---|
1404 | |
---|
1405 | FF=exp(-((ix-i0)**2+(iy-j0)**2)/sigma1**2) |
---|
1406 | GG=exp(-(iz-k0)**2/sigma2**2) |
---|
1407 | |
---|
1408 | bx0=bx(ix,iy,iz) |
---|
1409 | by0=by(ix,iy,iz) |
---|
1410 | bz0=bz(ix,iy,iz) |
---|
1411 | |
---|
1412 | bx(ix,iy,iz)=bx0+A_ics*freq**2*FF*GG*sin(arg) |
---|
1413 | ! Acc=A_ics*FF*GG*sin(arg) |
---|
1414 | !print *,'acceleration:', Acc |
---|
1415 | by(ix,iy,iz)=by0 |
---|
1416 | bz(ix,iy,iz)=bz0 |
---|
1417 | |
---|
1418 | end do |
---|
1419 | end do |
---|
1420 | end do |
---|
1421 | |
---|
1422 | |
---|
1423 | |
---|
1424 | end subroutine |
---|
1425 | |
---|
1426 | !pba : fin |
---|
1427 | |
---|
1428 | ! ************************************************************************* |
---|
1429 | subroutine urho(ax,ay,az,b,c) |
---|
1430 | ! ************************************************************************* |
---|
1431 | ! term u.b pour l'equation de la densite (calcule dans l'espace physique) |
---|
1432 | ! ************************************************************************** |
---|
1433 | implicit none |
---|
1434 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: c |
---|
1435 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: ax, ay, az, b |
---|
1436 | integer :: ix, iy, iz |
---|
1437 | b=b+c; |
---|
1438 | ! VERSION LINEARISEE |
---|
1439 | if (lin) then |
---|
1440 | do iz=0,nzmax-1 |
---|
1441 | do iy=0,nymax-1 |
---|
1442 | do ix=0,nxmax-1 |
---|
1443 | ax(ix,iy,iz)=b(ix,iy,iz)*vbx(ix,iy) |
---|
1444 | ay(ix,iy,iz)=b(ix,iy,iz)*vby(ix,iy) |
---|
1445 | az(ix,iy,iz)=b(ix,iy,iz)*vbz(ix,iy) |
---|
1446 | end do |
---|
1447 | end do |
---|
1448 | end do |
---|
1449 | ! NON LINEAIRE |
---|
1450 | else |
---|
1451 | ! calcul en nonlineaire perturbatif : |
---|
1452 | if (perturbe) then |
---|
1453 | do iz=0,nzmax-1 |
---|
1454 | do iy=0,nymax-1 |
---|
1455 | do ix=0,nxmax-1 |
---|
1456 | ax(ix,iy,iz)=b(ix,iy,iz)*(ax(ix,iy,iz)+vbx(ix,iy)) |
---|
1457 | ay(ix,iy,iz)=b(ix,iy,iz)*(ay(ix,iy,iz)+vby(ix,iy)) |
---|
1458 | az(ix,iy,iz)=b(ix,iy,iz)*(az(ix,iy,iz)+vbz(ix,iy)) |
---|
1459 | end do |
---|
1460 | end do |
---|
1461 | end do |
---|
1462 | ! calcul en non perturbatif : |
---|
1463 | else |
---|
1464 | do iz=0,nzmax-1 |
---|
1465 | do iy=0,nymax-1 |
---|
1466 | do ix=0,nxmax-1 |
---|
1467 | ax(ix,iy,iz)=b(ix,iy,iz)*ax(ix,iy,iz) |
---|
1468 | ay(ix,iy,iz)=b(ix,iy,iz)*ay(ix,iy,iz) |
---|
1469 | az(ix,iy,iz)=b(ix,iy,iz)*az(ix,iy,iz) |
---|
1470 | end do |
---|
1471 | end do |
---|
1472 | end do |
---|
1473 | end if |
---|
1474 | end if |
---|
1475 | |
---|
1476 | end subroutine |
---|
1477 | |
---|
1478 | ! ************************************************************************* |
---|
1479 | subroutine div_speciale(ax,ay,az,bz) |
---|
1480 | ! ************************************************************************* |
---|
1481 | ! Cette sous-routine renvoie dans le champ bz le terme bz - div(ax,ay,az) |
---|
1482 | ! utile pour evaluer le terme non lineaire sur l'equation de la densite en economisant un champ de stockage... |
---|
1483 | ! ************************************************************************** |
---|
1484 | |
---|
1485 | implicit none |
---|
1486 | |
---|
1487 | double precision ,dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: ax, ay, az |
---|
1488 | double precision ,dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bz |
---|
1489 | integer :: ix, iy, iz |
---|
1490 | |
---|
1491 | !if (Rho_read) |
---|
1492 | |
---|
1493 | do iz=0,nzmax-1 |
---|
1494 | do iy=0,nymax-1 |
---|
1495 | do ix=0,kxmax |
---|
1496 | bz(1,ix,iy,iz) = kx(ix)*ax(2,ix,iy,iz)+ky(iy)*ay(2,ix,iy,iz)+kz(iz)*az(2,ix,iy,iz) |
---|
1497 | bz(2,ix,iy,iz) = -(kx(ix)*ax(1,ix,iy,iz)+ky(iy)*ay(1,ix,iy,iz)+kz(iz)*az(1,ix,iy,iz)) |
---|
1498 | end do |
---|
1499 | end do |
---|
1500 | end do |
---|
1501 | |
---|
1502 | !else |
---|
1503 | ! do iz=0,nzmax-1 |
---|
1504 | ! do iy=0,nymax-1 |
---|
1505 | ! do ix=0,kxmax |
---|
1506 | ! bz(1,ix,iy,iz) = bz(1,ix,iy,iz) + (kx(ix)*ax(2,ix,iy,iz)+ky(iy)*ay(2,ix,iy,iz)+kz(iz)*az(2,ix,iy,iz)) |
---|
1507 | ! bz(2,ix,iy,iz) = bz(2,ix,iy,iz) - (kx(ix)*ax(1,ix,iy,iz)+ky(iy)*ay(1,ix,iy,iz)+kz(iz)*az(1,ix,iy,iz)) |
---|
1508 | ! end do |
---|
1509 | ! end do |
---|
1510 | ! end do |
---|
1511 | ! end if |
---|
1512 | end subroutine |
---|
1513 | |
---|
1514 | ! ************************************************************************** |
---|
1515 | subroutine init_k(kx,k2x,ky,k2y,kz,k2z,cax,ca2x) |
---|
1516 | ! ************************************************************************** |
---|
1517 | ! Définition des vecteurs d'ondes |
---|
1518 | ! ************************************************************************** |
---|
1519 | |
---|
1520 | implicit none |
---|
1521 | |
---|
1522 | double precision, dimension(0:kxmax), intent(out) :: kx, k2x |
---|
1523 | double precision, dimension(0:nymax-1), intent(out) :: ky, k2y |
---|
1524 | double precision, dimension(0:nzmax-1), intent(out) :: kz, k2z |
---|
1525 | double precision, dimension(0:2*kxmax+1), intent(out) :: cax, ca2x |
---|
1526 | |
---|
1527 | #ifndef mode_kz |
---|
1528 | integer :: ix, iy, iz |
---|
1529 | #else |
---|
1530 | integer :: ix, iy |
---|
1531 | #endif |
---|
1532 | |
---|
1533 | do ix=0,kxmax |
---|
1534 | kx(ix)=ix*2.*pi/lx |
---|
1535 | end do |
---|
1536 | |
---|
1537 | do iy=0,nymax/2 |
---|
1538 | ky(iy)=iy*2.*pi/ly |
---|
1539 | end do |
---|
1540 | |
---|
1541 | do iy=nymax/2+1,nymax-1 |
---|
1542 | ky(iy)=(iy-nymax)*2.*pi/ly |
---|
1543 | end do |
---|
1544 | |
---|
1545 | #ifndef mode_kz |
---|
1546 | do iz=0,nzmax/2 |
---|
1547 | kz(iz)=iz*2.*pi/lz |
---|
1548 | end do |
---|
1549 | |
---|
1550 | do iz=nzmax/2+1,nzmax-1 |
---|
1551 | kz(iz)=(iz-nzmax)*2.*pi/lz |
---|
1552 | end do |
---|
1553 | #else |
---|
1554 | if (lz==0) then |
---|
1555 | kz = 0 |
---|
1556 | else |
---|
1557 | kz = 2.*pi/lz |
---|
1558 | end if |
---|
1559 | #endif |
---|
1560 | |
---|
1561 | k2x=kx*kx |
---|
1562 | k2y=ky*ky |
---|
1563 | k2z=kz*kz |
---|
1564 | |
---|
1565 | ! définition alternative des vecteurs d'onde suivant x |
---|
1566 | ! utile pour un certains nombre de sous-routines (améliore la vectorisation sur le NEC SX-5 de l'Idris...) |
---|
1567 | do ix=0,kxmax |
---|
1568 | cax(2*ix)=ix*2.*pi/lx |
---|
1569 | cax(2*ix+1)=ix*2.*pi/lx |
---|
1570 | end do |
---|
1571 | |
---|
1572 | ca2x=cax*cax |
---|
1573 | |
---|
1574 | end subroutine |
---|
1575 | |
---|
1576 | ! ************************************************************************* |
---|
1577 | subroutine nonli(vox,voy,voz,rho,fx,fy,fz,frho,Rhov) |
---|
1578 | ! ************************************************************************* |
---|
1579 | ! calcul des termes non linéaires f(t) et frho(t) en spectral |
---|
1580 | ! en entrée v(t) et rho(t) en spectral |
---|
1581 | |
---|
1582 | implicit none |
---|
1583 | |
---|
1584 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: Rhov |
---|
1585 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho |
---|
1586 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: fx, fy, fz, frho |
---|
1587 | |
---|
1588 | ! necessaire au calcul d'une vitesse de translation |
---|
1589 | #ifdef calcul_ut |
---|
1590 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1) :: dwz_dx, dwz_dy |
---|
1591 | integer :: ix, iy, iz |
---|
1592 | double precision :: u0, erreur |
---|
1593 | #endif |
---|
1594 | |
---|
1595 | ! sauvegarde de voz en sphérique dans le champ frho (ruse pour gagner de la place en memoire) |
---|
1596 | ! nécessaire pour le cas stratifié |
---|
1597 | if (stratification) then |
---|
1598 | call cpu_time(tps_temp) |
---|
1599 | tps_copie = tps_copie - tps_temp |
---|
1600 | frho=voz |
---|
1601 | call cpu_time(tps_temp) |
---|
1602 | tps_copie = tps_copie + tps_temp |
---|
1603 | end if |
---|
1604 | |
---|
1605 | ! calcul de la vorticite de vo. Le resultat est placé dans f. |
---|
1606 | call cpu_time(tps_temp) |
---|
1607 | tps_rotv = tps_rotv - tps_temp |
---|
1608 | call rotv(vox,voy,voz,fx,fy,fz) |
---|
1609 | call cpu_time(tps_temp) |
---|
1610 | tps_rotv = tps_rotv + tps_temp |
---|
1611 | |
---|
1612 | #ifdef calcul_ut |
---|
1613 | ! on determine la derivee de la vorticite verticale par rapport a x et par rapport a y |
---|
1614 | do iz=0, nzmax-1 |
---|
1615 | do iy=0, nymax-1 |
---|
1616 | do ix=0, kxmax |
---|
1617 | dwz_dx(2*ix,iy,iz) = -kx(ix)*fz(2*ix+1,iy,iz) |
---|
1618 | dwz_dx(2*ix+1,iy,iz) = kx(ix)*fz(2*ix,iy,iz) |
---|
1619 | |
---|
1620 | dwz_dy(2*ix,iy,iz) = -ky(iy)*fz(2*ix+1,iy,iz) |
---|
1621 | dwz_dy(2*ix+1,iy,iz) = ky(iy)*fz(2*ix,iy,iz) |
---|
1622 | end do |
---|
1623 | end do |
---|
1624 | end do |
---|
1625 | #endif |
---|
1626 | |
---|
1627 | ! on appelle ici éventuellement la fonction de sortie dans l'espace spectral out_spectral1 |
---|
1628 | call cpu_time(tps_temp) |
---|
1629 | tps_out = tps_out - tps_temp |
---|
1630 | |
---|
1631 | if (periode_out_spectral1 /= 0 .and. mod(it, periode_out_spectral1) == 0 .and. it /= 0) then |
---|
1632 | call out_spectral1(vox, voy, voz) |
---|
1633 | end if |
---|
1634 | |
---|
1635 | call cpu_time(tps_temp) |
---|
1636 | tps_out = tps_out + tps_temp |
---|
1637 | |
---|
1638 | ! passage dans l'espace physique de la vitesse, de la vorticite et de la densite |
---|
1639 | call cpu_time(tps_temp) |
---|
1640 | |
---|
1641 | tps_bck_fft = tps_bck_fft - tps_temp |
---|
1642 | |
---|
1643 | |
---|
1644 | call bck_fft(vox,fft_work) |
---|
1645 | call bck_fft(voy,fft_work) |
---|
1646 | call bck_fft(voz,fft_work) |
---|
1647 | call bck_fft(fx,fft_work) |
---|
1648 | call bck_fft(fy,fft_work) |
---|
1649 | call bck_fft(fz,fft_work) |
---|
1650 | |
---|
1651 | if (stratification) then |
---|
1652 | call bck_fft(rho,fft_work) |
---|
1653 | end if |
---|
1654 | |
---|
1655 | #ifdef calcul_ut |
---|
1656 | ! on passe dans l'espace physique les termes dwz_dx et dwz_dy |
---|
1657 | call bck_fft(dwz_dx,fft_work) |
---|
1658 | call bck_fft(dwz_dy,fft_work) |
---|
1659 | #endif |
---|
1660 | |
---|
1661 | call cpu_time(tps_temp) |
---|
1662 | tps_bck_fft = tps_bck_fft + tps_temp |
---|
1663 | |
---|
1664 | ! on appelle eventuellement les sorties periodiques dans l'espace physique |
---|
1665 | ! il y a deux sorties periodiques dans l'espace physique qui peuvent avoir des periodes |
---|
1666 | ! différentes |
---|
1667 | call cpu_time(tps_temp) |
---|
1668 | tps_out = tps_out - tps_temp |
---|
1669 | |
---|
1670 | if (periode_out_physique1 /= 0 .and. mod(it, periode_out_physique1) == 0 .and. it /= 0) then |
---|
1671 | call out_physique1(vox,voy,voz) |
---|
1672 | end if |
---|
1673 | |
---|
1674 | if (periode_out_physique2 /= 0 .and. mod(it, periode_out_physique2) == 0 .and. it /= 0) then |
---|
1675 | call out_physique2(vox,voy,voz,rho,fx,fy,fz) |
---|
1676 | end if |
---|
1677 | |
---|
1678 | call cpu_time(tps_temp) |
---|
1679 | tps_out = tps_out + tps_temp |
---|
1680 | |
---|
1681 | |
---|
1682 | ! calcul de V x W dans l'espace physique |
---|
1683 | call cpu_time(tps_temp) |
---|
1684 | tps_vecpro2 = tps_vecpro2 - tps_temp |
---|
1685 | call vecpro2(vox,voy,voz,fx,fy,fz) |
---|
1686 | !pba : deb |
---|
1687 | call forcage(fx,fy,fz) |
---|
1688 | !pba : fin |
---|
1689 | call cpu_time(tps_temp) |
---|
1690 | tps_vecpro2 = tps_vecpro2 + tps_temp |
---|
1691 | |
---|
1692 | #ifdef calcul_ut |
---|
1693 | ! on calcule une estimation de u0 |
---|
1694 | u0 = sum((vox(0:nxmax-1,0:nymax-1,0:nzmax-1)*dwz_dx(0:nxmax-1,0:nymax-1,0:nzmax-1)+voy(0:nxmax-1,0:nymax-1, & |
---|
1695 | 0:nzmax-1)*dwz_dy(0:nxmax-1,0:nymax-1,0:nzmax-1))*dwz_dx(0:nxmax-1,0:nymax-1,0:nzmax-1)) |
---|
1696 | |
---|
1697 | u0 = u0/max(sum(dwz_dx(0:nxmax-1,0:nymax-1,0:nzmax-1)**2), epsilo) |
---|
1698 | |
---|
1699 | ! on estime l'erreur que l'on commet |
---|
1700 | erreur = sum(((vox(0:nxmax-1,0:nymax-1,0:nzmax-1)-u0)*dwz_dx(0:nxmax-1,0:nymax-1,0:nzmax-1) + & |
---|
1701 | voy(0:nxmax-1,0:nymax-1,0:nzmax-1)*dwz_dy(0:nxmax-1,0:nymax-1,0:nzmax-1))**2)/(nxmax*nymax*nzmax) |
---|
1702 | |
---|
1703 | ! on imprime les resultats |
---|
1704 | write(*,'(A21,TR2,G23.10)') ' Estimation de ut =', u0 |
---|
1705 | write(*,'(A34,TR2,G23.10)') ' Estimation de l erreur sur ut =', erreur |
---|
1706 | |
---|
1707 | #endif |
---|
1708 | |
---|
1709 | ! calcul du terme non linéaire pour l'équation de la densité |
---|
1710 | if (stratification) then |
---|
1711 | |
---|
1712 | ! on calcul u*rho dans l'espace physique (ou u*(rho+Rhov) si profil de densite initial) |
---|
1713 | call cpu_time(tps_temp) |
---|
1714 | tps_urho = tps_urho - tps_temp |
---|
1715 | |
---|
1716 | |
---|
1717 | |
---|
1718 | call urho(vox,voy,voz,rho,Rhov) |
---|
1719 | call cpu_time(tps_temp) |
---|
1720 | tps_urho = tps_urho + tps_temp |
---|
1721 | |
---|
1722 | ! on passe u*rho (ou u*(rho+Rhov)) en spectral |
---|
1723 | call cpu_time(tps_temp) |
---|
1724 | tps_fwd_fft = tps_fwd_fft - tps_temp |
---|
1725 | |
---|
1726 | call fwd_fft(vox,fft_work) |
---|
1727 | call fwd_fft(voy,fft_work) |
---|
1728 | call fwd_fft(voz,fft_work) |
---|
1729 | |
---|
1730 | call cpu_time(tps_temp) |
---|
1731 | tps_fwd_fft = tps_fwd_fft + tps_temp |
---|
1732 | |
---|
1733 | ! la fonction div_speciale renvoie frho = voz - u*grad(rho) ou u*grad(rho+Rhov) si profil de densite initial |
---|
1734 | ! |
---|
1735 | ! u*grad(rho) est calcule en prenant la divergence de u*rho parce que div(u) = 0 |
---|
1736 | ! de plus on a ruse en stockant voz dans le champ frho pour gagner un champ en memoire |
---|
1737 | call cpu_time(tps_temp) |
---|
1738 | tps_div_speciale = tps_div_speciale - tps_temp |
---|
1739 | call div_speciale(vox,voy,voz,frho) |
---|
1740 | call cpu_time(tps_temp) |
---|
1741 | tps_div_speciale = tps_div_speciale + tps_temp |
---|
1742 | |
---|
1743 | ! ajout du terme de flottaison pour l'équation de la vitesse |
---|
1744 | !if (Rho_read) !on a lu un profil de densite initial |
---|
1745 | fz=fz-9.81*rho/Rhob |
---|
1746 | !else |
---|
1747 | !fz=fz-xnsqr*rho !si pas de profil de densité initial entré N est rescalé de manière à avoir dRho/dz=1 |
---|
1748 | |
---|
1749 | end if |
---|
1750 | |
---|
1751 | ! Retour dans l'espace spectral de f |
---|
1752 | call cpu_time(tps_temp) |
---|
1753 | tps_fwd_fft = tps_fwd_fft - tps_temp |
---|
1754 | |
---|
1755 | call fwd_fft(fx,fft_work) |
---|
1756 | call fwd_fft(fy,fft_work) |
---|
1757 | call fwd_fft(fz,fft_work) |
---|
1758 | |
---|
1759 | call cpu_time(tps_temp) |
---|
1760 | tps_fwd_fft = tps_fwd_fft + tps_temp |
---|
1761 | |
---|
1762 | end subroutine |
---|
1763 | |
---|
1764 | |
---|
1765 | ! ************************************************************************** |
---|
1766 | subroutine projection(ax,ay,az) |
---|
1767 | ! ************************************************************************** |
---|
1768 | ! calcul de la projection : Pij( V x W )(k) dans l'espace spectral. |
---|
1769 | ! ************************************************************************** |
---|
1770 | |
---|
1771 | implicit none |
---|
1772 | |
---|
1773 | double precision, dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: ax, ay, az |
---|
1774 | double precision inv_k2, s1, s2 |
---|
1775 | integer :: ix, iy, iz |
---|
1776 | |
---|
1777 | |
---|
1778 | do iz=0,nzmax-1 |
---|
1779 | do iy=0,nymax-1 |
---|
1780 | do ix=0,kxmax |
---|
1781 | #ifdef test_troncature |
---|
1782 | if (trunctab(ix,iy,iz)) then |
---|
1783 | #endif |
---|
1784 | ! On utilise epsilo pour éviter que le module de k ne soit nul |
---|
1785 | inv_k2=1./max(k2x(ix)+k2y(iy)+k2z(iz),epsilo) |
---|
1786 | |
---|
1787 | #ifndef mode_kz |
---|
1788 | s1 =(kx(ix)*ax(1,ix,iy,iz)+ky(iy)*ay(1,ix,iy,iz)+kz(iz)*az(1,ix,iy,iz))*inv_k2 |
---|
1789 | s2 =(kx(ix)*ax(2,ix,iy,iz)+ky(iy)*ay(2,ix,iy,iz)+kz(iz)*az(2,ix,iy,iz))*inv_k2 |
---|
1790 | #else |
---|
1791 | s1 =(kx(ix)*ax(1,ix,iy,iz)+ky(iy)*ay(1,ix,iy,iz)-kz(iz)*az(2,ix,iy,iz))*inv_k2 |
---|
1792 | s2 =(kx(ix)*ax(2,ix,iy,iz)+ky(iy)*ay(2,ix,iy,iz)+kz(iz)*az(1,ix,iy,iz))*inv_k2 |
---|
1793 | #endif |
---|
1794 | |
---|
1795 | ax(1,ix,iy,iz)= ax(1,ix,iy,iz)-s1*kx(ix) |
---|
1796 | ax(2,ix,iy,iz)= ax(2,ix,iy,iz)-s2*kx(ix) |
---|
1797 | ay(1,ix,iy,iz)= ay(1,ix,iy,iz)-s1*ky(iy) |
---|
1798 | ay(2,ix,iy,iz)= ay(2,ix,iy,iz)-s2*ky(iy) |
---|
1799 | |
---|
1800 | #ifndef mode_kz |
---|
1801 | az(1,ix,iy,iz)= az(1,ix,iy,iz)-s1*kz(iz) |
---|
1802 | az(2,ix,iy,iz)= az(2,ix,iy,iz)-s2*kz(iz) |
---|
1803 | #else |
---|
1804 | az(1,ix,iy,iz)= az(1,ix,iy,iz)-s2*kz(iz) |
---|
1805 | az(2,ix,iy,iz)= az(2,ix,iy,iz)+s1*kz(iz) |
---|
1806 | #endif |
---|
1807 | |
---|
1808 | |
---|
1809 | #ifdef test_troncature |
---|
1810 | end if |
---|
1811 | #endif |
---|
1812 | end do |
---|
1813 | end do |
---|
1814 | end do |
---|
1815 | |
---|
1816 | |
---|
1817 | end subroutine |
---|
1818 | |
---|
1819 | ! *************************************************************************** |
---|
1820 | subroutine precal2(dtl) |
---|
1821 | ! *************************************************************************** |
---|
1822 | ! Précalcul des termes exponentiels pour l'intégration en temps |
---|
1823 | ! *************************************************************************** |
---|
1824 | |
---|
1825 | implicit none |
---|
1826 | |
---|
1827 | double precision, intent(in) :: dtl |
---|
1828 | |
---|
1829 | eex=exp(-nu*ca2x*dtl) |
---|
1830 | exs=exp(-nu/schm*ca2x*dtl) |
---|
1831 | |
---|
1832 | eey=exp(-nu*k2y*dtl) |
---|
1833 | eeys=exp(-nu/schm*k2y*dtl) |
---|
1834 | |
---|
1835 | eez=exp(-nu*k2z*dtl) |
---|
1836 | eezs=exp(-nu/schm*k2z*dtl) |
---|
1837 | |
---|
1838 | end subroutine |
---|
1839 | |
---|
1840 | ! ********************************************************************** |
---|
1841 | subroutine raccro(vox,voy,voz,rho,fx,fy,fz,frho,vox2,voy2,voz2,rho2) |
---|
1842 | ! ********************************************************************** |
---|
1843 | ! Cette subroutine permet à partir du champs v(t) de déterminer |
---|
1844 | ! avec un schéma Runge-Kutta à l'ordre deux les champs v(t+dt) et f(t) |
---|
1845 | ! ********************************************************************** |
---|
1846 | |
---|
1847 | implicit none |
---|
1848 | |
---|
1849 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho |
---|
1850 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: fx, fy, fz, frho |
---|
1851 | double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox2, voy2, voz2, rho2 |
---|
1852 | |
---|
1853 | vox2=vox |
---|
1854 | voy2=voy |
---|
1855 | voz2=voz |
---|
1856 | if (stratification) then |
---|
1857 | rho2 = rho |
---|
1858 | end if |
---|
1859 | |
---|
1860 | ! les champs contiennent vo = V(t) / vo2 = V(t) / f = ## |
---|
1861 | |
---|
1862 | ! ************ |
---|
1863 | ! Etape 1 / Schema Euler avant : (V(t+dt/2)-V(t))/(dt/2) = f(V(t)) |
---|
1864 | ! ************ |
---|
1865 | |
---|
1866 | ! calcul des termes non linéaires f(V(t=0)) |
---|
1867 | call nonli(vox,voy,voz,rho,fx,fy,fz,frho,Rhov) |
---|
1868 | |
---|
1869 | vox=vox2 |
---|
1870 | voy=voy2 |
---|
1871 | voz=voz2 |
---|
1872 | if (stratification) then |
---|
1873 | rho = rho2 |
---|
1874 | end if |
---|
1875 | |
---|
1876 | ! les champs contiennent vo = V(t) / vo2 = V(t) / f = f(V(t)) |
---|
1877 | |
---|
1878 | ! precalcul des termes exponentiels associes au shéma Euler avant |
---|
1879 | call precal2(dt/2.) |
---|
1880 | |
---|
1881 | ! shema Euler avant proprement dit |
---|
1882 | call runge_step1(vox,voy,voz,rho,fx,fy,fz,frho) |
---|
1883 | |
---|
1884 | ! on projette et on tronque |
---|
1885 | call projection(vox,voy,voz) |
---|
1886 | if(truncate) then |
---|
1887 | call trunc0(vox,voy,voz,rho) |
---|
1888 | endif |
---|
1889 | |
---|
1890 | ! les champs contiennent vo = V(t+dt/2) / vo2 = V(t) / f = f(V(t)) |
---|
1891 | |
---|
1892 | ! ************ |
---|
1893 | ! Etape 2 / Schema leap-frog centre en t+dt/2 : ( V(t+dt)-V(t) )/(dt) = f( V(t+dt/2) ) |
---|
1894 | ! ************ |
---|
1895 | |
---|
1896 | call nonli(vox,voy,voz,rho,fx,fy,fz,frho,Rhov) |
---|
1897 | |
---|
1898 | vox=vox2 |
---|
1899 | voy=voy2 |
---|
1900 | voz=voz2 |
---|
1901 | rho=rho2 |
---|
1902 | |
---|
1903 | ! les champs contiennent vo = V(t) / vo2 = V(t) / f = f(V(t+dt/2)) |
---|
1904 | |
---|
1905 | ! precalcul des termes exponentiels associes au shéma Euler avant |
---|
1906 | call precal2(dt) |
---|
1907 | |
---|
1908 | ! shema leap-frog proprement dit |
---|
1909 | call runge_step2(vox,voy,voz,rho,fx,fy,fz,frho) |
---|
1910 | |
---|
1911 | ! les champs contiennent vo = V(t+dt) / vo2 = V(t) / f = f(V(t+dt/2)) |
---|
1912 | |
---|
1913 | ! ************ |
---|
1914 | ! Etape 3 / Recalcul du terme f(V(t)) qui sera utile par la suite pour le schéma d'Adams-Bashforth |
---|
1915 | ! ************ |
---|
1916 | call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho,Rhov) |
---|
1917 | |
---|
1918 | ! les champs contiennent vo = V(t+dt) / vo2 = V(t) / f = f(V(t)) |
---|
1919 | |
---|
1920 | end subroutine |
---|
1921 | |
---|
1922 | ! *************************************************************************** |
---|
1923 | subroutine init_coord(xx,yy,zz) |
---|
1924 | ! *************************************************************************** |
---|
1925 | ! Initialisations de tableaux de longueur |
---|
1926 | ! *************************************************************************** |
---|
1927 | |
---|
1928 | implicit none |
---|
1929 | |
---|
1930 | double precision, dimension(0:nxmax-1), intent(out) :: xx |
---|
1931 | double precision, dimension(0:nymax-1), intent(out) :: yy |
---|
1932 | double precision, dimension(0:nzmax-1), intent(out) :: zz |
---|
1933 | integer :: i |
---|
1934 | |
---|
1935 | do i=0,nxmax-1 |
---|
1936 | xx(i)=i*lx/nxmax |
---|
1937 | end do |
---|
1938 | |
---|
1939 | do i=0,nymax-1 |
---|
1940 | yy(i)=i*ly/nymax |
---|
1941 | end do |
---|
1942 | |
---|
1943 | do i=0,nzmax-1 |
---|
1944 | zz(i)=i*lz/nzmax |
---|
1945 | end do |
---|
1946 | |
---|
1947 | end subroutine |
---|
1948 | |
---|
1949 | end module |
---|