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