#include "config.h" module sous_routines use parametres use fft implicit none contains !********************************************************************************** !********************************************************************************** !********** SOUS ROUTINES D'ENTREE SORTIES ********** !********************************************************************************** !********************************************************************************** ! *************************************************************************** subroutine enter_param ! *************************************************************************** ! lecture des tous les parametres qui sont dans le fichier 'data.in' ! *************************************************************************** character(len=29) :: bidon character(len=100) :: lecture_reel, lecture_entier, lecture_logical integer :: i, nx,ny,nz lecture_reel = '(A29,E14.5)' lecture_entier = '(A29,I10)' lecture_logical = '(A29,L12)' open (unit=1,file='data.in') read(1,*) read(1,*) read(1,*) read(1,*) ! dimension du probleme (ces valeurs ne sont en fait pas reellement utilisees). ! les dimensions du probleme sont definies dans le fichier parametres.F90 par le preprocesseur ! en fonction des valeurs indiquees dans le fichier config.h read(1,lecture_entier) bidon,nx read(1,lecture_entier) bidon,ny read(1,lecture_entier) bidon,nz ! taille suivant x read(1,lecture_reel) bidon,lx read(1,lecture_reel) bidon,ly read(1,lecture_reel) bidon,lz ! Pas de temps et instant initial read(1,lecture_reel) bidon,dt read(1,lecture_reel) bidon,begin ! nombre d'itérations read(1,lecture_entier) bidon,itmax ! Troncature read(1,lecture_logical) bidon,truncate read(1,lecture_entier) bidon,type_trunc read(1,lecture_reel) bidon,rtrunc_x read(1,lecture_reel) bidon,rtrunc_y read(1,lecture_reel) bidon,rtrunc_z read(1,*) read(1,*) ! Viscosité read(1,lecture_reel) bidon,nu ! stratification, fréquence de Brunt Vaisala, vorticité planétaire et nombre de Schmidt read(1,lecture_logical) bidon,stratification read(1,lecture_reel) bidon,xns xnsqr=xns*xns read(1,lecture_reel) bidon,schm read(1,lecture_reel) bidon,omega2 read(1,*) read(1,*) ! type de run read(1,lecture_logical) bidon,perturbe read(1,lecture_logical) bidon,lin read(1,*) read(1,*) ! reprise de run read(1,lecture_logical) bidon,reprise_run read(1,*) read(1,*) ! etat de base 2D read(1,lecture_logical) bidon, base2D_sub read(1,lecture_logical) bidon, base2D_read read(1,*) read(1,*) ! vitesse initiale 2D read(1,lecture_logical) bidon, velo2D_sub read(1,lecture_logical) bidon, velo2D_read read(1,lecture_logical) bidon,lnoise read(1,lecture_reel) bidon,amplitude_noise read(1,*) read(1,*) !profil de densite read(1,lecture_logical) bidon, Rho_read read(1,lecture_reel) bidon, Rhob read(1,*) read(1,*) ! pba : deb ! forcage read(1,lecture_logical) bidon, forcage_sub read(1,lecture_reel) bidon, A_ics read(1,lecture_reel) bidon, freq read(1,lecture_entier) bidon, it0 read(1,*) read(1,*) !pba fin ! sorties read(1,lecture_entier) bidon,periode_out_physique1 read(1,lecture_entier) bidon,periode_out_physique2 read(1,lecture_entier) bidon,periode_out_spectral1 ! parametres additionnels read(1,*) read(1,*) ! Tourbillons gaussiens read(1,*) read(1,*) read(1,lecture_entier) bidon,nb_tourbillons do i=1,nb_tourbillons read(1,*) read(1,lecture_reel) bidon,xt(i) read(1,lecture_reel) bidon,yt(i) read(1,lecture_reel) bidon,gamma(i) read(1,lecture_reel) bidon,at(i) end do close(1) end subroutine ! *************************************************************************** subroutine out_param ! *************************************************************************** ! affichage des parametres du run ! *************************************************************************** implicit none character(len=100) :: aff_reel, aff_entier, aff_logical integer :: i ! descripteur des formats de sortie aff_reel = '(A29,TR3,G12.5)' aff_entier = '(A29,TR4,I7)' aff_logical = '(A29,TR4,L1)' write(*,*) write(*,*) '---------------------------------------------------' write(*,*) ' PROGRAMME NS3D' write(*,*) '---------------------------------------------------' write(*,*) write(*,*) ! Mode de calcul if (lin) then write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' write(*,*) '@@@ @@@' write(*,*) ' VERSION PERTURBATIVE / LINEARISEE ' write(*,*) '@@@ @@@' write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' else if (perturbe) then write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' write(*,*) '@@@ @@@' write(*,*) ' VERSION PERTURBATIVE / NON LINEARISEE ' write(*,*) '@@@ @@@' write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' else write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' write(*,*) '@@@ @@@' write(*,*) ' VERSION NON-PERTURBATIVE ' write(*,*) '@@@ @@@' write(*,*) '@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@' write(*,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' end if end if write(*,*) write(*,*) write(*,*) '---------------------------------------------------' write(*,*) 'PARAMETRES DE DISCRETISATION' write(*,*) '---------------------------------------------------' write(*,aff_entier) 'nx...........................',nxmax write(*,aff_entier) 'ny...........................',nymax write(*,aff_entier) 'nz...........................',nzmax write(*,aff_reel) 'Lx...........................',lx write(*,aff_reel) 'Ly...........................',ly write(*,aff_reel) 'Lz...........................',lz write(*,aff_reel) 'Dt...........................',dt write(*,aff_reel) 'begin........................',begin write(*,aff_entier) 'itmax........................',itmax write(*,aff_logical) 'Troncature...................',truncate write(*,aff_entier) ' carree:1 / elliptique:2..',type_trunc write(*,aff_reel) ' rayon troncature x.......',rtrunc_x write(*,aff_reel) ' rayon troncature y.......',rtrunc_y write(*,aff_reel) ' rayon troncature z........',rtrunc_z write(*,*) write(*,*) '---------------------------------------------------' write(*,*) 'PARAMETRES PHYSIQUES' write(*,*) '---------------------------------------------------' write(*,aff_reel) 'viscosite....................',nu write(*,aff_logical) 'stratification...............',stratification write(*,aff_reel) ' brunt_vaisala_frequency..',xns write(*,aff_reel) ' Schmidt number...........',schm write(*,aff_reel) 'Omega2.......................',omega2 write(*,*) write(*,*) '---------------------------------------------------' write(*,*) 'TYPE DE RUN' write(*,*) '---------------------------------------------------' write(*,aff_logical) 'en perturbation..............',perturbe write(*,aff_logical) 'lineaire.....................',lin write(*,*) write(*,*) '---------------------------------------------------' write(*,*) 'REPRISE DE RUN' write(*,*) '---------------------------------------------------' write(*,aff_logical) 'reprise de run...............',reprise_run write(*,*) write(*,*) '---------------------------------------------------' write(*,*) 'ETAT DE BASE 2D' write(*,*) '---------------------------------------------------' write(*,aff_logical) 'call subroutine..............',base2D_sub write(*,aff_logical) 'Lecture fichier..............',base2D_read write(*,*) write(*,*) '---------------------------------------------------' write(*,*) 'VITESSE INITIALE 2D' write(*,*) '---------------------------------------------------' write(*,aff_logical) 'call subroutine..............',velo2D_sub write(*,aff_logical) 'Lecture fichier..............',velo2D_read write(*,aff_logical) 'Bruit blanc..................',lnoise write(*,aff_reel) ' amplitude................',amplitude_noise write(*,*) write(*,*) '---------------------------------------------------' write(*,*) 'DENSITE INITIALE 1D' write(*,*) '---------------------------------------------------' write(*,aff_logical) 'Lecture fichier..............',Rho_read write(*,aff_reel) 'densite de base..............',Rhob write(*,*) write(*,*)'---------------------------------------------------' write(*,*)'FORCAGE EN VITESSE' write(*,*)'---------------------------------------------------' write(*,aff_logical)'call subroutine', forcage_sub write(*,aff_reel)'Amplitude..................',A_ics write(*,aff_reel)'Frequence..................',freq write(*,*) write(*,*) '---------------------------------------------------' write(*,*) 'SORTIES' write(*,*) '---------------------------------------------------' write(*,aff_entier) 'periode out_physique1.........',periode_out_physique1 write(*,aff_entier) 'periode out_physique2........',periode_out_physique2 write(*,aff_entier) 'periode out_spectral1........',periode_out_spectral1 write(*,*) write(*,*) '---------------------------------------------------' write(*,*) 'PARAMETRES ADDITIONNELS' write(*,*) '---------------------------------------------------' write(*,*) write(*,*) '** Tourbillons gaussiens **' write(*,aff_entier) 'nb tourbillons...............',nb_tourbillons do i=1,nb_tourbillons write(*,aff_entier) 'Tourbillon...................',i write(*,aff_reel) ' xt.......................',xt(i) write(*,aff_reel) ' xt.......................',yt(i) write(*,aff_reel) ' gamma....................',gamma(i) write(*,aff_reel) ' at.......................',at(i) end do write(*,*) write(*,*) end subroutine ! ********************************************************************** subroutine out_physique1(vox,voy,voz) ! ********************************************************************** ! Sortie effectuee de façon periodique et qui prend en argument le champ de vitesse ! dans l'espace physique. ! Ici on affiche le taux de croissance des perturbations ainsi que le temps cpu ecoule et ! une estimation du temps restant ! ********************************************************************** implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: vox, voy, voz double precision :: vquad, growthrate double precision, save :: tps_ecoule = 0, tps_ecoule_prec_iteration = 0 ! on calcule le temps cpu ecoule depuis la rentree dans la boucle et une estimation du temps ! restant (on se base sur le temps necessaire a la derniere iteration de la boucle) call cpu_time(tps_temp) tps_ecoule = tps_total + tps_temp write(*,*) write(*,'(A4,TR2,I6,TR3,A6,TR2,F10.3)') 'it =', it, 'time =', it*dt + begin write(*,'(A37,TR2,I10,TR3,I10)') ' temps cpu en sec (ecoule/restant) =', int(tps_ecoule), & int((tps_ecoule-tps_ecoule_prec_iteration)/periode_out_physique1*(itmax-it)) tps_ecoule_prec_iteration = tps_ecoule ! on calcule la vitesse quadratique moyenne des perturbations vquad = sum(vox(0:nxmax-1, 0:nymax-1, 0:nzmax-1)**2 + voy(0:nxmax-1, 0:nymax-1, 0:nzmax-1)**2 + & voz(0:nxmax-1, 0:nymax-1, 0:nzmax-1)**2) write(*,'(A32,TR2,G23.10)') ' vitesse quadratique moyenne =',vquad ! on calcule le taux de croissance des perturbations entre deux appels à la fonction out_physique1 if (vquadanc>epsilo) then growthrate=0.5*log(vquad/vquadanc)/(periode_out_physique1*dt) else growthrate=0. endif write(*,'(A34,TR2,G23.10)') ' taux de croissance instantane =', growthrate vquadanc = vquad end subroutine ! ************************************************************************* subroutine out_physique2(vox,voy,voz,rho,fx,fy,fz) ! ************************************************************************* ! Sortie effectuee de façon periodique et qui prend en argument le champ de vitesse ! la densite et la vorticite dans l'espace physique. ! ************************************************************************** double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: vox,voy,voz,rho double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx,fy,fz ! là on appelle simplement, une fonction qui écrit les champs dans un fichier. call out_champs(vox,voy,voz,rho,fx,fy,fz) end subroutine ! ************************************************************************* subroutine out_spectral1(ax,ay,az) ! ************************************************************************* ! Fonctions qui sauvegarde l'energie quadratique moyenné suivant ky dans l'espace spectral ! On fait d'abord une moyenne suivant ky et on enregistre l'energie pour chaque mode kx et kz ! ************************************************************************* implicit none double precision, dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: ax, ay, az double precision, dimension(0:kxmax, 0:nzmax-1) :: energie_quadratique_mode integer :: ix, iz do iz=0,nzmax-1 do ix=0,kxmax energie_quadratique_mode(ix, iz) = sum(ax(1,ix,0:nymax-1,iz)**2+ax(2,ix,0:nymax-1,iz)**2+ & ay(1,ix,0:nymax-1,iz)**2+ ay(2,ix,0:nymax-1,iz)**2+ & az(1,ix,0:nymax-1,iz)**2+az(2,ix,0:nymax-1,iz)**2)/nymax end do end do ! on l'écrit dans un fichier ! si c'est le premier appel à la fonction out_spectral1, on cree le fichier de sortie et on ecrit toutes les ! metadonnees. ! On pourrait eventuellement tester si le fichier existe deja inquire(file="energie_quadratique_mode", EXIST=existe) ! Cela permettrait de completer un fichier de sortie genere lors d'une precedente simulation... ! S'il ne s'agit pas du premier appel, on ecrit simplement l'instant et les valeurs de l'energie_quadratique_mode ! dans le fichier de sortie ! premier appel a periode_out_spectral1: on cree le fichier et on ecrit les metadonnees if (it == periode_out_spectral1) then open(unit=76, file='energie_quadratique_mode',form='unformatted',status="replace",action="write") write(76) 1, nxmax, nymax, nzmax, lx, ly, lz, dt, truncate, type_trunc, rtrunc_x, rtrunc_y, rtrunc_z, & nu, stratification, xns, schm, omega2, perturbe, lin close(76) end if ! ce n'est pas le premier appel : on ajoute a la fin du fichier existant l'instant et les valeurs de l'energie ! Si on ajoute l'argument action="write" dans la foncion write avec le compilateur Absoft Pro Fortran 8.2 sous Mac ! l'enregistrement ne se fait pas correctement alors que ca marche sans probleme avec le compilateur Intel Fortran ! ou le compilateur Fortran du NEC de l'IDRIS, bizarre, bizarre... open(unit=76, file='energie_quadratique_mode',form='unformatted',status="old",position="append") write(76) it*dt + begin write(76) energie_quadratique_mode close(76) end subroutine ! ************************************************************************* subroutine out_champs(vox,voy,voz,rho,fx,fy,fz) ! ************************************************************************* ! Assure la sortie des champs dans des fichiers ! ************************************************************************** implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: vox,voy,voz,rho double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx,fy,fz character(len=50) :: nomfic character(len=13) :: sufx call time2string(it*dt + begin,sufx) ! on ecrit la vitesse, la densite et la vorticite dans un fichier ! on écrit un certain nombre de metadonnees pour aider a identifier le fichier : ! c'est utile pour eviter bien des erreurs lors des entrees-sorties et faciliter le post-traitement... ! Cela permet notamment de s'assurer qu'on utilise le bon codage des données binaires (little- ! endian ou big-endian) et que les tableaux ont bien les dimensions requises. ! On écrit ensuite les sorties tableau par tableau dans le fichier final. ! Attention, Fortran ajoute au début et à la fin de chaque enregistrement, à chaque ! enregistrement binaire effectué avec l'instruction "write(,form='unformatted')", un entier ! correspondant a la taille de l'enregitrement en octets. Il faut faire attention ! et en tenir compte lors du traitement ulterieur éventuel des donnees (dans Matlab par ex) ! Le format final du fichier de sortie est donc (i et j sont des entiers correspondant ! à la taille de l'enregistrement qui respectivement suit et précède) : ! i,1,nxmax,nymax,nzmax,lx,ly,lz,truncate,type_trunc,rtrunc_x,rtrunc_y,rtrunc_z,nu,stratification, ! 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, ! i,vbx,vby,vbz,j,i,wbx,wby,wbz,j nomfic='velo_rho_vort.'//sufx open (76,file=nomfic,form='unformatted') write(76) 1, nxmax, nymax, nzmax, lx, ly, lz, dt, truncate, type_trunc, rtrunc_x, rtrunc_y, rtrunc_z, & nu, stratification, xns, schm, omega2, perturbe, lin write(76) it*dt + begin write(76) vox(0:nxmax-1,0:nymax-1,0:nzmax-1) write(76) voy(0:nxmax-1,0:nymax-1,0:nzmax-1) write(76) voz(0:nxmax-1,0:nymax-1,0:nzmax-1) write(76) rho(0:nxmax-1,0:nymax-1,0:nzmax-1) write(76) fx(0:nxmax-1,0:nymax-1,0:nzmax-1) write(76) fy(0:nxmax-1,0:nymax-1,0:nzmax-1) write(76) fz(0:nxmax-1,0:nymax-1,0:nzmax-1) write(76) vbx,vby,vbz write(76) wbx,wby,wbz close(76) print *,'itfic:', it end subroutine ! *********************************************************** subroutine read_run(vox,voy,voz,rho,fx,fy,fz,vbx,vby,vbz,wbx,wby,wbz,begin) ! *********************************************************** ! Lecture des champs d'un run precedent dans le fichier nom_fichier ! *********************************************************** implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1),intent(out) :: vox,voy,voz,rho,fx,fy,fz double precision, dimension(0:nxmax-1,0:nymax-1), intent(out) :: vbx,vby,vbz,wbx,wby,wbz double precision, intent(out) :: begin integer :: flag, nxread, nyread, nzread double precision :: lxread, lyread, lzread double precision :: dbl_poub integer :: int_poub logical :: log_poub open (unit=8,file='velocity.init',form='unformatted') ! on lit les metadonnees read(8) flag,nxread,nyread,nzread,lxread,lyread,lzread,dbl_poub,log_poub,int_poub,dbl_poub,dbl_poub,dbl_poub, & dbl_poub,log_poub,dbl_poub,dbl_poub,dbl_poub,log_poub,log_poub read(8) begin ! on verifie que le format du fichier d'entree est correct. On ne verifie que les dimensions du probleme car ce sont ! les seuls parametres qui ne peuvent vraiment pas etre changes entre deux simulations if (flag /= 1 .or. nxread /= nxmax .or. nyread /= nymax .or. nzread /= nzmax .or. lxread /= lx .or. lyread /= ly) then write(*,*) "ERREUR: le fichier de reprise de run n'a pas le bon format" write(*,*) "ARRET DU PROGRAMME..." stop end if read(8) vox(0:nxmax-1,0:nymax-1,0:nzmax-1) read(8) voy(0:nxmax-1,0:nymax-1,0:nzmax-1) read(8) voz(0:nxmax-1,0:nymax-1,0:nzmax-1) read(8) rho(0:nxmax-1,0:nymax-1,0:nzmax-1) read(8) fx(0:nxmax-1,0:nymax-1,0:nzmax-1) read(8) fy(0:nxmax-1,0:nymax-1,0:nzmax-1) read(8) fz(0:nxmax-1,0:nymax-1,0:nzmax-1) read(8) vbx,vby,vbz read(8) wbx,wby,wbz close(8) end subroutine ! on choisit quelle fonction de generation de l'etat de base 2D doit etre utilisee #if choix_base2D==1 ! *********************************************************** subroutine gen_base2D(vbx,vby,vbz,wbx,wby,wbz,vox2,voy2,voz2) ! *********************************************************** ! Generation de l'ecoulement de base ! Ecoulement de Stuart ! *********************************************************** implicit none double precision, dimension(0:nxmax-1,0:nymax-1), intent(out) :: vbx,vby,vbz,wbx,wby,wbz double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: vox2, voy2, voz2 integer :: ix, iy double precision ::rho_stuart vbx = 0. vby = 0. vbz = 0. wbx = 0. wby = 0. wbz = 0. vox2 = 0. voy2 = 0. voz2 = 0. ! on définit ici le parametre de concentration des tourbillons de stuart rho_stuart = 0.25 do iy=0, nymax-1 do ix = 0,nxmax-1 vbx(ix,iy) = sinh(2*(yy(iy)-ly/2))/(cosh(2* (yy(iy)-ly/2) ) - rho_stuart*cos( 2*xx(ix) ) ) vby(ix,iy) = - rho_stuart * sin( 2*xx(ix) )/( cosh( 2*( yy(iy)-ly/2 ) ) - rho_stuart*cos(2*xx(ix))) wbz(ix,iy) = 2*(rho_stuart**2-1)/(cosh(2*(yy(iy)-ly/2)) - rho_stuart*cos(2*xx(ix)))**2 end do end do end subroutine #elif choix_base2D==2 ! *********************************************************** subroutine gen_base2D(vbx,vby,vbz,wbx,wby,wbz,vox2,voy2,voz2) ! *********************************************************** ! Generation de l'ecoulement de base ! Couche de mélange avec un profil vbx=tanh(y) ! *********************************************************** implicit none double precision, dimension(0:nxmax-1,0:nymax-1), intent(out) :: vbx,vby,vbz,wbx,wby,wbz double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: vox2, voy2, voz2 integer :: ix vbx = 0. vby = 0. vbz = 0. wbx = 0. wby = 0. wbz = 0. vox2 = 0. voy2 = 0. voz2 = 0. do ix = 0,nxmax-1 vbx(ix,:) = tanh(yy-ly/2) wbz(ix,:) = -(1-tanh(yy-ly/2)**2) end do end subroutine #elif choix_base2D==3 ! *********************************************************** subroutine gen_base2D(vbx,vby,vbz,wbx,wby,wbz,vox2,voy2,voz2) ! *********************************************************** ! Generation de l'ecoulement de base ! Superposition de tourbillons dont les caracteristiques ! sont donnees dans le fichier de parametres ! *********************************************************** implicit none double precision, dimension(0:nxmax-1,0:nymax-1), intent(out) :: vbx,vby,vbz,wbx,wby,wbz double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: vox2, voy2, voz2 double precision k2 integer num_tourbillon integer ix, iy, iz vbx = 0. vby = 0. vbz = 0. wbx = 0. wby = 0. wbz = 0. voz2 = 0. ! on determine la vorticite verticale de l'ecoulement do num_tourbillon=1,nb_tourbillons do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,nxmax-1 voz2(ix,iy,iz)=voz2(ix,iy,iz)+wz_gau((xx(ix)-xt(num_tourbillon)),(yy(iy)-yt(num_tourbillon)), & gamma(num_tourbillon),at(num_tourbillon)) end do end do end do end do ! initialisation de la vorticite verticale wbz = voz2(0:nxmax-1,0:nymax-1,0) ! inversion du laplacien 2D pour trouver la fct de courant 2D*1D call fwd_fft(voz2,fft_work) do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,2*kxmax+1 k2=max(ca2x(ix)+k2y(iy),epsilo) voz2(ix,iy,iz)= voz2(ix,iy,iz)/k2 end do end do end do ! vitesses en fct de la fct de courant 2D*1D do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,kxmax vox2(2*ix,iy,iz) = -ky(iy)*voz2(2*ix+1,iy,iz) vox2(2*ix+1,iy,iz) = ky(iy)*voz2(2*ix,iy,iz) voy2(2*ix,iy,iz) = kx(ix)*voz2(2*ix+1,iy,iz) voy2(2*ix+1,iy,iz) = -kx(ix)*voz2(2*ix,iy,iz) end do end do end do ! passage dans l'espace physique call bck_fft(vox2,fft_work) call bck_fft(voy2,fft_work) vbx = vox2(0:nxmax-1,0:nymax-1,0) vby = voy2(0:nxmax-1,0:nymax-1,0) end subroutine #endif ! ******************************************************************************** subroutine read_base2D(vbx,vby,vbz,wbx,wby,wbz) ! ******************************************************************************** ! Subroutine qui permet de lire un champ de base en formatted ou unformatted ! ******************************************************************************** implicit none double precision, dimension(0:nxmax-1,0:nymax-1), intent(out) :: vbx, vby, vbz, wbx, wby, wbz integer :: flag, nxread, nyread double precision :: lxread, lyread open (unit=88,file='base2D.init',form='unformatted') read(88) flag, nxread, nyread, lxread, lyread read(88) vbx, vby, vbz read(88) wbx, wby, wbz close(88) ! on verifie que le format du fichier d'entree est correct if (flag /= 1 .or. nxread /= nxmax .or. nyread /= nymax .or. lxread /= lx .or. lyread /= ly) then write(*,*) "ERREUR: le fichier de reprise de run n'a pas le bon format" write(*,*) "ARRET DU PROGRAMME..." stop end if end subroutine ! ******************************************************************************** subroutine read_velo2D(vox, voy, voz) ! ******************************************************************************** ! Subroutine qui permet de lire un champ 2D en unformatted ! ******************************************************************************** implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: vox, voy, voz double precision, dimension(0:nxmax-1,0:nymax-1) :: vi2Dx, vi2Dy, vi2Dz integer :: flag, nxread, nyread double precision :: lxread, lyread integer :: ix,iy,iz open (unit=88,file='velo2D.init',form='unformatted') read(88) flag, nxread, nyread, lxread, lyread read(88) vi2Dx,vi2Dy,vi2Dz close(88) ! on verifie que le format du fichier d'entree est correct if (flag /= 1 .or. nxread /= nxmax .or. nyread /= nymax .or. lxread /= lx .or. lyread /= ly) then write(*,*) "ERREUR: le fichier de reprise de run n'a pas le bon format" write(*,*) "ARRET DU PROGRAMME..." stop end if do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,nxmax-1 vox(ix,iy,iz)=vi2Dx(ix,iy) voy(ix,iy,iz)=vi2Dy(ix,iy) voz(ix,iy,iz)=vi2Dz(ix,iy) end do end do end do end subroutine ! ******************************************************************************** subroutine readRho(Rhov) !ajout d'une sous routine de lecture d'un gradient de densite moyen non constant ! ******************************************************************************** ! Subroutine qui permet de lire un champ 2D en unformatted ! ******************************************************************************** implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: Rhov double precision, dimension(0:nymax-1) :: Rhoi integer :: flag, nxread, nyread double precision :: lxread, lyread integer :: ix,iy,iz write(*,*) 'boucle lecture rho' open (unit=88,file='Rho.init',status='old') !read(88) flag, nxread, nyread, lxread, lyread do iy=0,255 read(88,*) Rhoi(iy) enddo close(88) ! on verifie que le format du fichier d'entree est correct ! if (flag /= 1 .or. nyread /= nymax .or. lyread /= ly) then ! write(*,*) "ERREUR: le fichier de reprise de run n'a pas le bon format" ! write(*,*) "ARRET DU PROGRAMME..." ! stop !end if do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,nxmax-1 Rhov(ix,iy,iz)=Rhoi(iy) end do end do end do print*, Rhov end subroutine ! ************************************************************************************ subroutine gen_velo2D(vox, voy, voz) ! ************************************************************************************ ! Generation d'un champ de vitesse correspondant a un dipole fait de ! tourbillons de Lamb-oseen. ! ************************************************************************************ implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: vox,voy,voz double precision k2 integer num_tourbillon integer ix, iy, iz vox = 0. voy = 0. voz = 0. ! on determine la vorticite verticale de l'ecoulement do num_tourbillon=1,nb_tourbillons do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,nxmax-1 voz(ix,iy,iz)=voz(ix,iy,iz)+wz_gau((xx(ix)-xt(num_tourbillon)),(yy(iy)-yt(num_tourbillon)), & gamma(num_tourbillon),at(num_tourbillon)) end do end do end do end do ! inversion du laplacien 2D pour trouver la fct de courant 2D*1D call fwd_fft(voz,fft_work) do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,2*kxmax+1 k2=max(ca2x(ix)+k2y(iy),epsilo) voz(ix,iy,iz)= voz(ix,iy,iz)/k2 end do end do end do ! vitesses en fct de la fct de courant 2D*1D do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,kxmax vox(2*ix,iy,iz) = -ky(iy)*voz(2*ix+1,iy,iz) vox(2*ix+1,iy,iz) = ky(iy)*voz(2*ix,iy,iz) voy(2*ix,iy,iz) = kx(ix)*voz(2*ix+1,iy,iz) voy(2*ix+1,iy,iz) = -kx(ix)*voz(2*ix,iy,iz) end do end do end do ! passage dans l'espace physique call bck_fft(vox,fft_work) call bck_fft(voy,fft_work) voz = 0 end subroutine ! *************************************************************** function wz_gau(x,y,gamma,a) ! *************************************************************** implicit none double precision, intent(in) :: x, y, gamma,a double precision :: wz_gau double precision :: r2 r2 = x**2 + y**2 wz_gau =gamma/(pi*a**2)*exp(-r2/a**2) end function ! ************************************************************************* subroutine white_noise(vox,voy,voz) ! ********************************************************* implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz double precision :: rnd_number, moy_x, moy_y, moy_z integer :: ix, iy, iz ! on initialise le generateur de nombre aleatoire avec une nouvelle "graine" ! ne semble pas fontionner sur toutes les architectures (bizarre, bizarre...) ! on desactive l'initialisation du generateur de nombre aleatoire en mode debugage de façon ! a pouvoir comparer plus facilement les resultats de runs differents. #ifndef debug call random_seed #endif ! on determine la moyenne initiale des champs vox, voy, voz moy_x = sum(vox(0:nxmax-1,0:nymax-1,0:nzmax-1))/(nxmax*nymax*nzmax) moy_y = sum(voy(0:nxmax-1,0:nymax-1,0:nzmax-1))/(nxmax*nymax*nzmax) moy_z = sum(voz(0:nxmax-1,0:nymax-1,0:nzmax-1))/(nxmax*nymax*nzmax) ! on ajoute le bruit blanc dans l'espace physique do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,nxmax-1 call random_number(rnd_number) vox(ix,iy,iz)=vox(ix,iy,iz)+amplitude_noise*rnd_number call random_number(rnd_number) voy(ix,iy,iz)=voy(ix,iy,iz)+amplitude_noise*rnd_number call random_number(rnd_number) voz(ix,iy,iz)=voz(ix,iy,iz)+amplitude_noise*rnd_number end do end do end do ! on fait en sorte que la moyenne du bruit que l'on a ajouté soit nulle ! (on utilise pour cela la difference des moyennes des champs vo avant et après ajout du bruit...) moy_x = moy_x - sum(vox(0:nxmax-1,0:nymax-1,0:nzmax-1))/(nxmax*nymax*nzmax) moy_y = moy_y - sum(voy(0:nxmax-1,0:nymax-1,0:nzmax-1))/(nxmax*nymax*nzmax) moy_z = moy_z - sum(voz(0:nxmax-1,0:nymax-1,0:nzmax-1))/(nxmax*nymax*nzmax) vox = vox + moy_x voy = voy + moy_y voz = voz + moy_z end subroutine ! ************************************************************************* subroutine time2string(time,string_time) ! ************************************************************************* ! CONVERSION TIME TO STRING t=0000.000 ! ************************************************************************* implicit none double precision, intent(in) :: time character(len=13), intent(out) :: string_time write(string_time,'(A2,I7.7,A1,I3.3)') 't=',int(time),'.',nint(1000*(time-int(time))) end subroutine !********************************************************************************** !********************************************************************************** !********** FIN DES ENTREE SORTIES ***************** !********************************************************************************** !********************************************************************************** ! ********************************************************************************* subroutine adams_bashforth(vox,voy,voz,rho,fx,fy,fz,frho,vox2,voy2,voz2,rho2) ! ********************************************************************************* ! coeur du programme avance temporelle sur itmax ! Time discretisation: second order Adams-Bashforth: ! step t+dt: input fields are : V(k,t) and F(k,t-dt) ! output fields are: V(k,t+dt) and F(k,t) implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: fx, fy, fz, frho double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox2, voy2, voz2, rho2 call cpu_time(tps_temp) tps_copie = tps_copie - tps_temp vox2 = vox voy2 = voy voz2 = voz if (stratification) then rho2 = rho end if call cpu_time(tps_temp) tps_copie = tps_copie + tps_temp ! les champs contiennent vo = V(t) / vo2 = V(t) / f = f(V(t-dt)) ! ************ ! Etape 1 / V(t*,k)=E*V(t,k)-(dt/2)*E**2*F(t-dt,k) ! avec E=exp(-nu*k2*dt) des constantes precalculees ! ************ call cpu_time(tps_temp) tps_adams_step1 = tps_adams_step1 - tps_temp call adams_step1(vox,voy,voz,rho,fx,fy,fz,frho) call cpu_time(tps_temp) tps_adams_step1 = tps_adams_step1 + tps_temp ! les champs contiennent vo = V(t*) / vo2 = V(t) / f = f(V(t-dt)) ! ************ ! Etape 2 / V(t+dt) = V(t*)+(3.*dt/2.)*E*F(V(t)) ! avec E=exp(-nu*k2*dt) des constantes precalculees ! ************ call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho,Rhov) call cpu_time(tps_temp) tps_adams_step2 = tps_adams_step2 - tps_temp call adams_step2(vox,voy,voz,rho,fx,fy,fz,frho) call cpu_time(tps_temp) tps_adams_step2 = tps_adams_step2 + tps_temp ! les champs contiennent vo = V(t+dt) / vo2 = ## / f = f(V(t)) end subroutine ! ************************************************************************* subroutine adams_step1(vox,voy,voz,rho,fx,fy,fz,frho) ! ************************************************************************* ! Première étape du schéma d'ordre deux d'Adams-Bashforth ! ************************************************************************* implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx, fy, fz, frho double precision :: ee,ees,e2,e2s double precision :: temp_dt integer :: ix, iy, iz temp_dt=dt/2. do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,2*kxmax+1 #ifdef test_troncature if (trunctab(ix,iy,iz)) then #endif ee=eex(ix)*eey(iy)*eez(iz) ees=exs(ix)*eeys(iy)*eezs(iz) e2=ee*ee e2s=ees*ees vox(ix,iy,iz)=ee*vox(ix,iy,iz)-temp_dt*e2*fx(ix,iy,iz) voy(ix,iy,iz)=ee*voy(ix,iy,iz)-temp_dt*e2*fy(ix,iy,iz) voz(ix,iy,iz)=ee*voz(ix,iy,iz)-temp_dt*e2*fz(ix,iy,iz) rho(ix,iy,iz)=ees*rho(ix,iy,iz)-temp_dt*e2s*frho(ix,iy,iz) #ifdef test_troncature end if #endif end do end do end do end subroutine ! ************************************************************************* subroutine adams_step2(vox,voy,voz,rho,fx,fy,fz,frho) ! ************************************************************************* ! Deuxième étape du schéma d'ordre deux d'Adams-Bashforth ! ************************************************************************* implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx, fy, fz, frho double precision :: ee,ees double precision :: temp_dt integer :: ix, iy, iz temp_dt=3.*dt/2. do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,2*kxmax+1 #ifdef test_troncature if (trunctab(ix,iy,iz)) then #endif ee=eex(ix)*eey(iy)*eez(iz) ees=exs(ix)*eeys(iy)*eezs(iz) vox(ix,iy,iz)=vox(ix,iy,iz)+ee*temp_dt*fx(ix,iy,iz) voy(ix,iy,iz)=voy(ix,iy,iz)+ee*temp_dt*fy(ix,iy,iz) voz(ix,iy,iz)=voz(ix,iy,iz)+ee*temp_dt*fz(ix,iy,iz) rho(ix,iy,iz)=rho(ix,iy,iz)+ees*temp_dt*frho(ix,iy,iz) #ifdef test_troncature end if #endif end do end do end do end subroutine ! ************************************************************************* subroutine runge_step2(vox,voy,voz,rho,fx,fy,fz,frho) ! ************************************************************************* ! schema de type leap-frog centre a t+dt/2. ! Utilise comme etape du subroutine raccro ! qui effectue un pas Runge-Kutta a l'ordre deux pour redemarrer. ! ( v(n+1/2)+1/2 - v(n+1/2)-1/2 )/(2.dt/2) = f( v(n+1/2) ) . ! ************************************************************************** implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx, fy, fz, frho double precision :: k2,ee,a,ees,as integer :: ix, iy, iz do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,2*kxmax+1 k2=max(ca2x(ix)+k2y(iy)+k2z(iz),epsilo) ee=eex(ix)*eey(iy)*eez(iz) a=(1.-ee)/(nu*k2) ees=exs(ix)*eeys(iy)*eezs(iz) as=(1.-ees)/(nu/schm*k2) vox(ix,iy,iz)=vox(ix,iy,iz)*ee+fx(ix,iy,iz)*a voy(ix,iy,iz)=voy(ix,iy,iz)*ee+fy(ix,iy,iz)*a voz(ix,iy,iz)=voz(ix,iy,iz)*ee+fz(ix,iy,iz)*a rho(ix,iy,iz)=rho(ix,iy,iz)*ees+frho(ix,iy,iz)*as end do end do end do end subroutine ! ************************************************************************* subroutine runge_step1(vox,voy,voz,rho,fx,fy,fz,frho) ! ************************************************************************* ! schema de type Euler. ! Utilisé comme étape du subroutine raccro ! qui effectue un pas Runge-Kutta a l'ordre deux pour redemarrer. ! ( v(n+1/2)-v(n))/(dt/2) = f( v(n) ) ! ************************************************************************* implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx, fy, fz, frho double precision :: k2,ee,a,ees,as integer :: ix, iy, iz do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,2*kxmax+1 k2=max(ca2x(ix)+k2y(iy)+k2z(iz),epsilo) ee=eex(ix)*eey(iy)*eez(iz) a=(1.-ee)/(nu*k2) ees=exs(ix)*eeys(iy)*eezs(iz) as=(1.-ees)/(nu/schm*k2) vox(ix,iy,iz)=vox(ix,iy,iz)*ee+fx(ix,iy,iz)*a voy(ix,iy,iz)=voy(ix,iy,iz)*ee+fy(ix,iy,iz)*a voz(ix,iy,iz)=voz(ix,iy,iz)*ee+fz(ix,iy,iz)*a rho(ix,iy,iz)=rho(ix,iy,iz)*ees+frho(ix,iy,iz)*as end do end do end do end subroutine ! ************************************************************************* subroutine rotv(ax,ay,az,bx,by,bz) ! ************************************************************************* ! calcul du rotationnel du champ (ax,ay,az). ! ************************************************************************* implicit none double precision, dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: ax, ay, az double precision, dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: bx, by, bz integer :: ix, iy, iz do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,kxmax #ifndef mode_kz bx(1,ix,iy,iz)=kz(iz)*ay(2,ix,iy,iz)-ky(iy)*az(2,ix,iy,iz) bx(2,ix,iy,iz)=ky(iy)*az(1,ix,iy,iz)-kz(iz)*ay(1,ix,iy,iz) by(1,ix,iy,iz)=kx(ix)*az(2,ix,iy,iz)-kz(iz)*ax(2,ix,iy,iz) by(2,ix,iy,iz)=kz(iz)*ax(1,ix,iy,iz)-kx(ix)*az(1,ix,iy,iz) #else bx(1,ix,iy,iz)=-kz(iz)*ay(1,ix,iy,iz)-ky(iy)*az(2,ix,iy,iz) bx(2,ix,iy,iz)=ky(iy)*az(1,ix,iy,iz)-kz(iz)*ay(2,ix,iy,iz) by(1,ix,iy,iz)=kx(ix)*az(2,ix,iy,iz)+kz(iz)*ax(1,ix,iy,iz) by(2,ix,iy,iz)=kz(iz)*ax(2,ix,iy,iz)-kx(ix)*az(1,ix,iy,iz) #endif bz(1,ix,iy,iz)=ky(iy)*ax(2,ix,iy,iz)-kx(ix)*ay(2,ix,iy,iz) bz(2,ix,iy,iz)=kx(ix)*ay(1,ix,iy,iz)-ky(iy)*ax(1,ix,iy,iz) end do end do end do end subroutine ! ************************************************************************* subroutine trunc0(ax,ay,az,b) ! ************************************************************************* ! classical desaliasing in the spectral space. ! ************************************************************************** implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: ax, ay, az, b integer :: ix, iy, iz do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,2*kxmax+1 if(.not. trunctab(ix,iy,iz)) then ax(ix,iy,iz) = 0 ay(ix,iy,iz) = 0 az(ix,iy,iz) = 0 b(ix,iy,iz) = 0 end if end do end do end do end subroutine ! ************************************************************************* subroutine init_trunctab(trunctab) ! ************************************************************************* ! Génération du tableau utilisé lors de la troncature ! ************************************************************************** implicit none logical, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: trunctab double precision :: kx2max,ky2max,kz2max,dks integer :: ix, iy, iz ! rayon du cercle: kmax*rtrunc kx2max=(rtrunc_x*nxmax/2*2*pi/lx)**2 ky2max=(rtrunc_y*nymax/2*2*pi/ly)**2 #ifndef mode_kz kz2max=(rtrunc_z*nzmax/2*2*pi/lz)**2 #else kz2max=rtrunc_z**2 #endif trunctab = .false. do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,2*kxmax+1 if (type_trunc==1) then dks = max(ca2x(ix)/kx2max, k2y(iy)/ky2max, k2z(iz)/kz2max) else dks = ca2x(ix)/kx2max + k2y(iy)/ky2max + k2z(iz)/kz2max end if ! on rajoute un epsilo pour eviter que numériquement, un terme soit tronqué alors qu'il ! ne devrait pas l'être. if (dks <= 1. + epsilo) then trunctab(ix,iy,iz)=.true. end if end do end do end do end subroutine ! ************************************************************************* subroutine vecpro2(ax,ay,az,bx,by,bz) ! ************************************************************************* ! produits vectoriels dans l'espace physique pour le code perturbatif. ! inclut les termes avec l'ecoulement de base et coriolis ! ************************************************************************** implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bx, by, bz double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: ax, ay, az double precision :: vbx0,vby0,vbz0,wbx0,wby0,wbz0 double precision :: ax0,ay0,az0,bx0,by0,bz0 integer :: ix, iy, iz ! VERSION LINEARISEE if (lin) then do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,nxmax-1 vbx0=vbx(ix,iy) vby0=vby(ix,iy) vbz0=vbz(ix,iy) wbx0=wbx(ix,iy) wby0=wby(ix,iy) wbz0=wbz(ix,iy)+omega2 ax0=ax(ix,iy,iz) ay0=ay(ix,iy,iz) az0=az(ix,iy,iz) bx0=bx(ix,iy,iz) by0=by(ix,iy,iz) bz0=bz(ix,iy,iz) bx(ix,iy,iz)=vby0*bz0-vbz0*by0 + ay0*wbz0-az0*wby0 by(ix,iy,iz)=vbz0*bx0-vbx0*bz0 + az0*wbx0-ax0*wbz0 bz(ix,iy,iz)=vbx0*by0-vby0*bx0 + ax0*wby0-ay0*wbx0 end do end do end do ! NON LINEAIRE else ! calcul en nonlineaire perturbatif : if (perturbe) then do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,nxmax-1 vbx0=vbx(ix,iy) vby0=vby(ix,iy) vbz0=vbz(ix,iy) wbx0=wbx(ix,iy) wby0=wby(ix,iy) wbz0=wbz(ix,iy)+omega2 ax0=ax(ix,iy,iz) ay0=ay(ix,iy,iz) az0=az(ix,iy,iz) bx0=bx(ix,iy,iz) by0=by(ix,iy,iz) bz0=bz(ix,iy,iz) bx(ix,iy,iz)=vby0*bz0-vbz0*by0 + ay0*(wbz0+bz0) - az0*(wby0+by0) by(ix,iy,iz)=vbz0*bx0-vbx0*bz0 + az0*(wbx0+bx0) - ax0*(wbz0+bz0) bz(ix,iy,iz)=vbx0*by0-vby0*bx0 + ax0*(wby0+by0) - ay0*(wbx0+bx0) end do end do end do ! calcul en non perturbatif : else do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,nxmax-1 ax0=ax(ix,iy,iz) ay0=ay(ix,iy,iz) az0=az(ix,iy,iz) bx0=bx(ix,iy,iz) by0=by(ix,iy,iz) bz0=bz(ix,iy,iz) bx(ix,iy,iz)=ay0*(omega2+bz0)-az0*by0 by(ix,iy,iz)=az0*bx0-ax0*(omega2+bz0) bz(ix,iy,iz)=ax0*by0-ay0*bx0 end do end do end do end if end if end subroutine ! pba: deb ! ************************************************************************* subroutine forcage(bx,by,bz) ! ************************************************************************* ! ajout d'un terme de forcage ! parametres du forcage entres ici provisoirement modifier ! data.in etc pour entree plus propre ! ************************************************************************** implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bx, by, bz double precision :: vbx0,vby0,vbz0,wbx0,wby0,wbz0 double precision :: ax0,ay0,az0,bx0,by0,bz0 double precision :: time, arg, sigma1, sigma2, FF, GG, acc integer :: ix, iy, iz, itime, i0, j0, k0 !itime a faire passer en parametre? ! itime=it0 ! itime=itime+1 ! penser a it0=it0+itmax(run precedent si restat ! print *,'temps:', itime ! A_ics=100. !deplacement a definir dans data.in et parametres.F90 !freq=2*pi/(3600*12.4) !M2 a verifier idem arg=freq*it*Dt !Dt est defini dans data.in ok ! print *,'arg:', arg ! print *,'A_ics:', A_ics !acceleration ! localisation centre de la source i0=floor(0.5*nxmax) j0=floor(0.5*nymax) k0=floor(0.5*nzmax) ! enveloppe gaussienne !sigma1=4km, sigma2=100m revoir les valeurs ci-dessous sigma1=20. ! defini en unite de dx et dz sigma2=1. ! dy do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,nxmax-1 FF=exp(-((ix-i0)**2+(iy-j0)**2)/sigma1**2) GG=exp(-(iz-k0)**2/sigma2**2) bx0=bx(ix,iy,iz) by0=by(ix,iy,iz) bz0=bz(ix,iy,iz) bx(ix,iy,iz)=bx0+A_ics*freq**2*FF*GG*sin(arg) ! Acc=A_ics*FF*GG*sin(arg) !print *,'acceleration:', Acc by(ix,iy,iz)=by0 bz(ix,iy,iz)=bz0 end do end do end do end subroutine !pba : fin ! ************************************************************************* subroutine urho(ax,ay,az,b,c) ! ************************************************************************* ! term u.b pour l'equation de la densite (calcule dans l'espace physique) ! ************************************************************************** implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: c double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: ax, ay, az, b integer :: ix, iy, iz b=b+c; ! VERSION LINEARISEE if (lin) then do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,nxmax-1 ax(ix,iy,iz)=b(ix,iy,iz)*vbx(ix,iy) ay(ix,iy,iz)=b(ix,iy,iz)*vby(ix,iy) az(ix,iy,iz)=b(ix,iy,iz)*vbz(ix,iy) end do end do end do ! NON LINEAIRE else ! calcul en nonlineaire perturbatif : if (perturbe) then do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,nxmax-1 ax(ix,iy,iz)=b(ix,iy,iz)*(ax(ix,iy,iz)+vbx(ix,iy)) ay(ix,iy,iz)=b(ix,iy,iz)*(ay(ix,iy,iz)+vby(ix,iy)) az(ix,iy,iz)=b(ix,iy,iz)*(az(ix,iy,iz)+vbz(ix,iy)) end do end do end do ! calcul en non perturbatif : else do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,nxmax-1 ax(ix,iy,iz)=b(ix,iy,iz)*ax(ix,iy,iz) ay(ix,iy,iz)=b(ix,iy,iz)*ay(ix,iy,iz) az(ix,iy,iz)=b(ix,iy,iz)*az(ix,iy,iz) end do end do end do end if end if end subroutine ! ************************************************************************* subroutine div_speciale(ax,ay,az,bz) ! ************************************************************************* ! Cette sous-routine renvoie dans le champ bz le terme bz - div(ax,ay,az) ! utile pour evaluer le terme non lineaire sur l'equation de la densite en economisant un champ de stockage... ! ************************************************************************** implicit none double precision ,dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: ax, ay, az double precision ,dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bz integer :: ix, iy, iz !if (Rho_read) do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,kxmax 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) 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)) end do end do end do !else ! do iz=0,nzmax-1 ! do iy=0,nymax-1 ! do ix=0,kxmax ! 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)) ! 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)) ! end do ! end do ! end do ! end if end subroutine ! ************************************************************************** subroutine init_k(kx,k2x,ky,k2y,kz,k2z,cax,ca2x) ! ************************************************************************** ! Définition des vecteurs d'ondes ! ************************************************************************** implicit none double precision, dimension(0:kxmax), intent(out) :: kx, k2x double precision, dimension(0:nymax-1), intent(out) :: ky, k2y double precision, dimension(0:nzmax-1), intent(out) :: kz, k2z double precision, dimension(0:2*kxmax+1), intent(out) :: cax, ca2x #ifndef mode_kz integer :: ix, iy, iz #else integer :: ix, iy #endif do ix=0,kxmax kx(ix)=ix*2.*pi/lx end do do iy=0,nymax/2 ky(iy)=iy*2.*pi/ly end do do iy=nymax/2+1,nymax-1 ky(iy)=(iy-nymax)*2.*pi/ly end do #ifndef mode_kz do iz=0,nzmax/2 kz(iz)=iz*2.*pi/lz end do do iz=nzmax/2+1,nzmax-1 kz(iz)=(iz-nzmax)*2.*pi/lz end do #else if (lz==0) then kz = 0 else kz = 2.*pi/lz end if #endif k2x=kx*kx k2y=ky*ky k2z=kz*kz ! définition alternative des vecteurs d'onde suivant x ! utile pour un certains nombre de sous-routines (améliore la vectorisation sur le NEC SX-5 de l'Idris...) do ix=0,kxmax cax(2*ix)=ix*2.*pi/lx cax(2*ix+1)=ix*2.*pi/lx end do ca2x=cax*cax end subroutine ! ************************************************************************* subroutine nonli(vox,voy,voz,rho,fx,fy,fz,frho,Rhov) ! ************************************************************************* ! calcul des termes non linéaires f(t) et frho(t) en spectral ! en entrée v(t) et rho(t) en spectral implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: Rhov double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: fx, fy, fz, frho ! necessaire au calcul d'une vitesse de translation #ifdef calcul_ut double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1) :: dwz_dx, dwz_dy integer :: ix, iy, iz double precision :: u0, erreur #endif ! sauvegarde de voz en sphérique dans le champ frho (ruse pour gagner de la place en memoire) ! nécessaire pour le cas stratifié if (stratification) then call cpu_time(tps_temp) tps_copie = tps_copie - tps_temp frho=voz call cpu_time(tps_temp) tps_copie = tps_copie + tps_temp end if ! calcul de la vorticite de vo. Le resultat est placé dans f. call cpu_time(tps_temp) tps_rotv = tps_rotv - tps_temp call rotv(vox,voy,voz,fx,fy,fz) call cpu_time(tps_temp) tps_rotv = tps_rotv + tps_temp #ifdef calcul_ut ! on determine la derivee de la vorticite verticale par rapport a x et par rapport a y do iz=0, nzmax-1 do iy=0, nymax-1 do ix=0, kxmax dwz_dx(2*ix,iy,iz) = -kx(ix)*fz(2*ix+1,iy,iz) dwz_dx(2*ix+1,iy,iz) = kx(ix)*fz(2*ix,iy,iz) dwz_dy(2*ix,iy,iz) = -ky(iy)*fz(2*ix+1,iy,iz) dwz_dy(2*ix+1,iy,iz) = ky(iy)*fz(2*ix,iy,iz) end do end do end do #endif ! on appelle ici éventuellement la fonction de sortie dans l'espace spectral out_spectral1 call cpu_time(tps_temp) tps_out = tps_out - tps_temp if (periode_out_spectral1 /= 0 .and. mod(it, periode_out_spectral1) == 0 .and. it /= 0) then call out_spectral1(vox, voy, voz) end if call cpu_time(tps_temp) tps_out = tps_out + tps_temp ! passage dans l'espace physique de la vitesse, de la vorticite et de la densite call cpu_time(tps_temp) tps_bck_fft = tps_bck_fft - tps_temp call bck_fft(vox,fft_work) call bck_fft(voy,fft_work) call bck_fft(voz,fft_work) call bck_fft(fx,fft_work) call bck_fft(fy,fft_work) call bck_fft(fz,fft_work) if (stratification) then call bck_fft(rho,fft_work) end if #ifdef calcul_ut ! on passe dans l'espace physique les termes dwz_dx et dwz_dy call bck_fft(dwz_dx,fft_work) call bck_fft(dwz_dy,fft_work) #endif call cpu_time(tps_temp) tps_bck_fft = tps_bck_fft + tps_temp ! on appelle eventuellement les sorties periodiques dans l'espace physique ! il y a deux sorties periodiques dans l'espace physique qui peuvent avoir des periodes ! différentes call cpu_time(tps_temp) tps_out = tps_out - tps_temp if (periode_out_physique1 /= 0 .and. mod(it, periode_out_physique1) == 0 .and. it /= 0) then call out_physique1(vox,voy,voz) end if if (periode_out_physique2 /= 0 .and. mod(it, periode_out_physique2) == 0 .and. it /= 0) then call out_physique2(vox,voy,voz,rho,fx,fy,fz) end if call cpu_time(tps_temp) tps_out = tps_out + tps_temp ! calcul de V x W dans l'espace physique call cpu_time(tps_temp) tps_vecpro2 = tps_vecpro2 - tps_temp call vecpro2(vox,voy,voz,fx,fy,fz) !pba : deb call forcage(fx,fy,fz) !pba : fin call cpu_time(tps_temp) tps_vecpro2 = tps_vecpro2 + tps_temp #ifdef calcul_ut ! on calcule une estimation de u0 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, & 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)) u0 = u0/max(sum(dwz_dx(0:nxmax-1,0:nymax-1,0:nzmax-1)**2), epsilo) ! on estime l'erreur que l'on commet erreur = sum(((vox(0:nxmax-1,0:nymax-1,0:nzmax-1)-u0)*dwz_dx(0:nxmax-1,0:nymax-1,0:nzmax-1) + & 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) ! on imprime les resultats write(*,'(A21,TR2,G23.10)') ' Estimation de ut =', u0 write(*,'(A34,TR2,G23.10)') ' Estimation de l erreur sur ut =', erreur #endif ! calcul du terme non linéaire pour l'équation de la densité if (stratification) then ! on calcul u*rho dans l'espace physique (ou u*(rho+Rhov) si profil de densite initial) call cpu_time(tps_temp) tps_urho = tps_urho - tps_temp call urho(vox,voy,voz,rho,Rhov) call cpu_time(tps_temp) tps_urho = tps_urho + tps_temp ! on passe u*rho (ou u*(rho+Rhov)) en spectral call cpu_time(tps_temp) tps_fwd_fft = tps_fwd_fft - tps_temp call fwd_fft(vox,fft_work) call fwd_fft(voy,fft_work) call fwd_fft(voz,fft_work) call cpu_time(tps_temp) tps_fwd_fft = tps_fwd_fft + tps_temp ! la fonction div_speciale renvoie frho = voz - u*grad(rho) ou u*grad(rho+Rhov) si profil de densite initial ! ! u*grad(rho) est calcule en prenant la divergence de u*rho parce que div(u) = 0 ! de plus on a ruse en stockant voz dans le champ frho pour gagner un champ en memoire call cpu_time(tps_temp) tps_div_speciale = tps_div_speciale - tps_temp call div_speciale(vox,voy,voz,frho) call cpu_time(tps_temp) tps_div_speciale = tps_div_speciale + tps_temp ! ajout du terme de flottaison pour l'équation de la vitesse !if (Rho_read) !on a lu un profil de densite initial fz=fz-9.81*rho/Rhob !else !fz=fz-xnsqr*rho !si pas de profil de densité initial entré N est rescalé de manière à avoir dRho/dz=1 end if ! Retour dans l'espace spectral de f call cpu_time(tps_temp) tps_fwd_fft = tps_fwd_fft - tps_temp call fwd_fft(fx,fft_work) call fwd_fft(fy,fft_work) call fwd_fft(fz,fft_work) call cpu_time(tps_temp) tps_fwd_fft = tps_fwd_fft + tps_temp end subroutine ! ************************************************************************** subroutine projection(ax,ay,az) ! ************************************************************************** ! calcul de la projection : Pij( V x W )(k) dans l'espace spectral. ! ************************************************************************** implicit none double precision, dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: ax, ay, az double precision inv_k2, s1, s2 integer :: ix, iy, iz do iz=0,nzmax-1 do iy=0,nymax-1 do ix=0,kxmax #ifdef test_troncature if (trunctab(ix,iy,iz)) then #endif ! On utilise epsilo pour éviter que le module de k ne soit nul inv_k2=1./max(k2x(ix)+k2y(iy)+k2z(iz),epsilo) #ifndef mode_kz s1 =(kx(ix)*ax(1,ix,iy,iz)+ky(iy)*ay(1,ix,iy,iz)+kz(iz)*az(1,ix,iy,iz))*inv_k2 s2 =(kx(ix)*ax(2,ix,iy,iz)+ky(iy)*ay(2,ix,iy,iz)+kz(iz)*az(2,ix,iy,iz))*inv_k2 #else s1 =(kx(ix)*ax(1,ix,iy,iz)+ky(iy)*ay(1,ix,iy,iz)-kz(iz)*az(2,ix,iy,iz))*inv_k2 s2 =(kx(ix)*ax(2,ix,iy,iz)+ky(iy)*ay(2,ix,iy,iz)+kz(iz)*az(1,ix,iy,iz))*inv_k2 #endif ax(1,ix,iy,iz)= ax(1,ix,iy,iz)-s1*kx(ix) ax(2,ix,iy,iz)= ax(2,ix,iy,iz)-s2*kx(ix) ay(1,ix,iy,iz)= ay(1,ix,iy,iz)-s1*ky(iy) ay(2,ix,iy,iz)= ay(2,ix,iy,iz)-s2*ky(iy) #ifndef mode_kz az(1,ix,iy,iz)= az(1,ix,iy,iz)-s1*kz(iz) az(2,ix,iy,iz)= az(2,ix,iy,iz)-s2*kz(iz) #else az(1,ix,iy,iz)= az(1,ix,iy,iz)-s2*kz(iz) az(2,ix,iy,iz)= az(2,ix,iy,iz)+s1*kz(iz) #endif #ifdef test_troncature end if #endif end do end do end do end subroutine ! *************************************************************************** subroutine precal2(dtl) ! *************************************************************************** ! Précalcul des termes exponentiels pour l'intégration en temps ! *************************************************************************** implicit none double precision, intent(in) :: dtl eex=exp(-nu*ca2x*dtl) exs=exp(-nu/schm*ca2x*dtl) eey=exp(-nu*k2y*dtl) eeys=exp(-nu/schm*k2y*dtl) eez=exp(-nu*k2z*dtl) eezs=exp(-nu/schm*k2z*dtl) end subroutine ! ********************************************************************** subroutine raccro(vox,voy,voz,rho,fx,fy,fz,frho,vox2,voy2,voz2,rho2) ! ********************************************************************** ! Cette subroutine permet à partir du champs v(t) de déterminer ! avec un schéma Runge-Kutta à l'ordre deux les champs v(t+dt) et f(t) ! ********************************************************************** implicit none double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: fx, fy, fz, frho double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox2, voy2, voz2, rho2 vox2=vox voy2=voy voz2=voz if (stratification) then rho2 = rho end if ! les champs contiennent vo = V(t) / vo2 = V(t) / f = ## ! ************ ! Etape 1 / Schema Euler avant : (V(t+dt/2)-V(t))/(dt/2) = f(V(t)) ! ************ ! calcul des termes non linéaires f(V(t=0)) call nonli(vox,voy,voz,rho,fx,fy,fz,frho,Rhov) vox=vox2 voy=voy2 voz=voz2 if (stratification) then rho = rho2 end if ! les champs contiennent vo = V(t) / vo2 = V(t) / f = f(V(t)) ! precalcul des termes exponentiels associes au shéma Euler avant call precal2(dt/2.) ! shema Euler avant proprement dit call runge_step1(vox,voy,voz,rho,fx,fy,fz,frho) ! on projette et on tronque call projection(vox,voy,voz) if(truncate) then call trunc0(vox,voy,voz,rho) endif ! les champs contiennent vo = V(t+dt/2) / vo2 = V(t) / f = f(V(t)) ! ************ ! Etape 2 / Schema leap-frog centre en t+dt/2 : ( V(t+dt)-V(t) )/(dt) = f( V(t+dt/2) ) ! ************ call nonli(vox,voy,voz,rho,fx,fy,fz,frho,Rhov) vox=vox2 voy=voy2 voz=voz2 rho=rho2 ! les champs contiennent vo = V(t) / vo2 = V(t) / f = f(V(t+dt/2)) ! precalcul des termes exponentiels associes au shéma Euler avant call precal2(dt) ! shema leap-frog proprement dit call runge_step2(vox,voy,voz,rho,fx,fy,fz,frho) ! les champs contiennent vo = V(t+dt) / vo2 = V(t) / f = f(V(t+dt/2)) ! ************ ! Etape 3 / Recalcul du terme f(V(t)) qui sera utile par la suite pour le schéma d'Adams-Bashforth ! ************ call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho,Rhov) ! les champs contiennent vo = V(t+dt) / vo2 = V(t) / f = f(V(t)) end subroutine ! *************************************************************************** subroutine init_coord(xx,yy,zz) ! *************************************************************************** ! Initialisations de tableaux de longueur ! *************************************************************************** implicit none double precision, dimension(0:nxmax-1), intent(out) :: xx double precision, dimension(0:nymax-1), intent(out) :: yy double precision, dimension(0:nzmax-1), intent(out) :: zz integer :: i do i=0,nxmax-1 xx(i)=i*lx/nxmax end do do i=0,nymax-1 yy(i)=i*ly/nymax end do do i=0,nzmax-1 zz(i)=i*lz/nzmax end do end subroutine end module