Changeset 9


Ignore:
Timestamp:
12/05/06 10:39:48 (17 years ago)
Author:
xlvlod
Message:

Ajout modifs Yannis

Location:
trunk/NS3D/SRC
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/NS3D/SRC/NS3D.F90

    r8 r9  
    345345    ! ****************************** 
    346346 
     347    if (Rho_read) then 
     348   
     349     write(*,*) "  -> lecture du fichier 'Rho.init'" 
     350     call readRho(Rhov) 
     351     
     352    end if 
     353 
     354 
    347355    if (reprise_run) then 
    348356        write(*,*) "Reprise d'un run" 
    … …  
    466474    ! Globalement ici, on rentre donc avec les valeurs à l'instant 1 et on en sort avec les valeurs à l'instant itmax 
    467475    do it=2,itmax 
    468  
     476        write(*,*) 'boucle it2, itmax'  
    469477        ! Shéma d'ordre deux d'Adams-Bashforth 
    470478        call adams_bashforth(vox,voy,voz,rho,fx,fy,fz,frho,vox2,voy2,voz2,rho2) 

  • trunk/NS3D/SRC/config.h

    r8 r9  
    1616#define DIMX 256 
    1717#define DIMY 256 
    18 #define DIMZ 256 
     18#define DIMZ 3 
    1919 
    2020! utile pour definir un mode "debugage" dans le code (pas d'initialisation du generateur de nombre aleatoire, 
    … …  
    3939! de l'architecture de calcul (NEC, processeur Intel...). 
    4040#define PPADX 0 
    41 #define PPADY 1 
    42 #define PPADZ 256 
     41#define PPADY 0 
     42#define PPADZ 0 
    4343 
    4444! choix de la fonction qui genere l'ecoulement de base. Actuellement est disponible : 

  • trunk/NS3D/SRC/data.in

    r8 r9  
    66ny___________________________   256 
    77nz___________________________   3 
    8 Lx___________________________   14.1322 
    9 Ly___________________________   60. 
     8Lx___________________________   3.E4 
     9Ly___________________________   25.E2 
    1010Lz___________________________   10. 
    11 Dt___________________________   0.01 
     11Dt___________________________   10. 
    1212begin________________________   0. 
    13 itmax________________________   15000 
     13itmax________________________   50000 
    1414Troncature_oui/non___________   T 
    1515    carree:1 / elliptique:2__   1 
    … …  
    1919 
    2020*** Parametres physiques *** 
    21 visc_________________________   5.E-6 
    22 stratifie____________________   F 
    23     brunt_vaisala_frequency__   10. 
     21visc_________________________   5.E-3 
     22stratifie____________________   T 
     23    brunt_vaisala_frequency__   5.E-3 
    2424    schmidt number___________   1. 
    25252omega_______________________   0. 
    2626 
    2727*** Type de simulation *** 
    28 En perturbation______________   T 
    29 Lineaire_____________________   T 
     28En perturbation______________   F 
     29Lineaire_____________________   F 
    3030 
    3131*** Reprise de run  *** 
    … …  
    3333 
    3434*** Etat de base (pour simulation en perturbation) *** 
    35 Call subroutine gen_base2D___   T 
     35Call subroutine gen_base2D___   F 
    3636Lecture fichier______________   F 
    3737 
    … …  
    3939Call subroutine gen_velo2D___   F 
    4040Lecture fichier______________   F 
    41 Bruit_blanc__________________   T 
     41Bruit_blanc__________________   F 
    4242    amplitude________________   0.00001 
    4343 
     44*** gradient densite initiale 1D *** 
     45Lecture fichier______________   T 
     46Rhob_________________________   1.E3 
     47 
     48*** forcage *** 
     49Call subroutine forcage______   T 
     50    A_ics____________________   15E2 
     51    freq_____________________   15.E-5 
     52      it0____________________   0 
     53 
    4454*** Sorties *** 
    45 periode out_physique1________   1 
    46 periode out_physique2________   2000 
    47 periode out_spectral1________   100 
     55periode out_physique1________   1000 
     56periode out_physique2________   1000 
     57periode out_spectral1________   1000 
    4858 
    4959*** Parametres additionnels... *** 
    … …  
    6171    Circulation______________   -6.28 
    6272    Rayon____________________   1. 
     73 

  • trunk/NS3D/SRC/parametres.F90

    r8 r9  
    5959     ! lecture de la vitesse initiale 2D en "unformatted" ou appel subroutine 
    6060    logical, save :: velo2D_sub, velo2D_read 
     61      
     62    ! lecture du profil de densite initiale 1D en "unformatted" 
     63    logical, save :: Rho_read 
     64    double precision, save :: Rhob 
    6165 
    6266    ! Bruit blanc et amplitude du bruit blanc 
    6367    logical, save :: lnoise 
    6468    double precision, save :: amplitude_noise 
     69 
     70    !introduction d'un forcage, appel de la subroutine forcage 
     71    logical, save :: forcage_sub    
     72    double precision, save :: A_ics, freq 
     73    integer, save :: it0   
    6574 
    6675    ! nombre et instants où sont appeles les differentes fonctions de sorties 
    … …  
    8897 
    8998    ! Vitesses, densité, termes non linéaire pour l'équation de la vitesse et de la densité 
    90     double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), save :: vox, voy, voz, rho 
     99    double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), save :: vox, voy, voz, rho, Rhov 
    91100    double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), save :: vox2, voy2, voz2, rho2 
    92101    double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), save :: fx, fy, fz, frho 

  • trunk/NS3D/SRC/sous_routines.F90

    r8 r9  
    3131 
    3232            open (unit=1,file='data.in') 
     33 
     34 
    3335 
    3436            read(1,*) 
    … …  
    107109            read(1,*) 
    108110            read(1,*) 
     111 
     112              
     113            !profil de densite  
     114            read(1,lecture_logical) bidon, Rho_read 
     115            read(1,lecture_reel) bidon, Rhob 
     116              
     117            read(1,*) 
     118            read(1,*) 
     119 
     120            ! pba : deb 
     121            ! forcage 
     122            read(1,lecture_logical) bidon, forcage_sub 
     123            read(1,lecture_reel) bidon, A_ics 
     124            read(1,lecture_reel) bidon, freq 
     125            read(1,lecture_entier) bidon, it0 
     126 
     127            read(1,*) 
     128            read(1,*) 
     129 
     130            !pba fin 
    109131 
    110132            ! sorties 
    … …  
    235257            write(*,aff_logical) 'Lecture fichier..............',base2D_read 
    236258 
     259 
    237260            write(*,*) 
    238261            write(*,*) '---------------------------------------------------' 
    … …  
    243266            write(*,aff_logical) 'Bruit blanc..................',lnoise 
    244267            write(*,aff_reel) '    amplitude................',amplitude_noise 
     268             
     269            write(*,*) 
     270            write(*,*) '---------------------------------------------------' 
     271            write(*,*) 'DENSITE INITIALE 1D' 
     272            write(*,*) '---------------------------------------------------' 
     273            write(*,aff_logical) 'Lecture fichier..............',Rho_read 
     274            write(*,aff_reel) 'densite de base..............',Rhob 
     275                
     276        
     277            write(*,*) 
     278            write(*,*)'---------------------------------------------------' 
     279            write(*,*)'FORCAGE EN VITESSE' 
     280            write(*,*)'---------------------------------------------------' 
     281            write(*,aff_logical)'call subroutine', forcage_sub 
     282            write(*,aff_reel)'Amplitude..................',A_ics 
     283            write(*,aff_reel)'Frequence..................',freq 
    245284 
    246285            write(*,*) 
    … …  
    401440            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx,fy,fz 
    402441            character(len=50) :: nomfic 
    403             character(len=10) :: sufx 
     442            character(len=13) :: sufx 
    404443 
    405444            call time2string(it*dt + begin,sufx) 
    … …  
    427466 
    428467            nomfic='velo_rho_vort.'//sufx 
    429  
     468             
    430469            open (76,file=nomfic,form='unformatted') 
    431470            write(76) 1, nxmax, nymax, nzmax, lx, ly, lz, dt, truncate, type_trunc, rtrunc_x, rtrunc_y, rtrunc_z, & 
    … …  
    443482 
    444483            close(76) 
    445  
     484     print *,'itfic:', it 
    446485        end subroutine 
    447486 
    … …  
    489528 
    490529            close(8) 
    491  
     530         
    492531        end subroutine 
    493532 
    … …  
    706745 
    707746        end subroutine 
     747 
     748 
     749! ******************************************************************************** 
     750        subroutine readRho(Rhov)  
     751!ajout d'une sous routine de lecture d'un gradient de densite moyen non constant 
     752! ******************************************************************************** 
     753!       Subroutine qui permet de lire un champ 2D en unformatted 
     754! ******************************************************************************** 
     755 
     756            implicit none 
     757 
     758            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: Rhov 
     759            double precision, dimension(0:nymax-1) :: Rhoi 
     760            integer :: flag, nxread, nyread 
     761            double precision :: lxread, lyread 
     762            integer :: ix,iy,iz 
     763            write(*,*) 'boucle lecture rho' 
     764            open (unit=88,file='Rho.init',status='old') 
     765            !read(88) flag, nxread, nyread, lxread, lyread 
     766            do iy=0,255 
     767            read(88,*) Rhoi(iy) 
     768            enddo 
     769             
     770            close(88) 
     771 
     772            ! on verifie que le format du fichier d'entree est correct 
     773            ! if (flag /= 1 .or. nyread /= nymax .or. lyread /= ly) then 
     774            !    write(*,*) "ERREUR: le fichier de reprise de run n'a pas le bon format" 
     775            !    write(*,*) "ARRET DU PROGRAMME..." 
     776            !    stop 
     777            !end if 
     778             
     779            do iz=0,nzmax-1 
     780                do iy=0,nymax-1 
     781                    do ix=0,nxmax-1 
     782                        Rhov(ix,iy,iz)=Rhoi(iy) 
     783                    end do 
     784                end do 
     785            end do 
     786            print*, Rhov 
     787 
     788        end subroutine 
     789 
     790 
    708791 
    709792 
    … …  
    844927 
    845928            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))) 
     929            character(len=13), intent(out) :: string_time 
     930 
     931            write(string_time,'(A2,I7.7,A1,I3.3)') 't=',int(time),'.',nint(1000*(time-int(time))) 
    849932      end subroutine 
    850933 
    … …  
    902985            ! ************ 
    903986 
    904             call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho) 
     987            call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho,Rhov) 
    905988 
    906989            call cpu_time(tps_temp) 
    … …  
    9201003!       Première étape du schéma d'ordre deux d'Adams-Bashforth 
    9211004! ************************************************************************* 
     1005 
    9221006 
    9231007            implicit none 
    … …  
    12741358        end subroutine 
    12751359 
    1276 ! ************************************************************************* 
    1277         subroutine urho(ax,ay,az,b) 
     1360! pba: deb 
     1361! ************************************************************************* 
     1362        subroutine forcage(bx,by,bz) 
     1363! ************************************************************************* 
     1364!       ajout d'un terme de forcage 
     1365!       parametres du forcage entres ici provisoirement modifier 
     1366!       data.in etc pour entree plus propre 
     1367! ************************************************************************** 
     1368 
     1369 
     1370            implicit none 
     1371 
     1372            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bx, by, bz 
     1373            double precision :: vbx0,vby0,vbz0,wbx0,wby0,wbz0 
     1374            double precision :: ax0,ay0,az0,bx0,by0,bz0 
     1375            double precision :: time, arg, sigma1, sigma2, FF, GG, acc 
     1376            integer :: ix, iy, iz, itime, i0, j0, k0    !itime a faire passer en parametre? 
     1377                ! itime=it0 
     1378 
     1379                ! itime=itime+1   ! penser a it0=it0+itmax(run precedent si restat 
     1380                ! print *,'temps:', itime  
     1381                !    A_ics=100.             !deplacement   a definir dans  data.in et parametres.F90 
     1382                !freq=2*pi/(3600*12.4)  !M2 a verifier idem 
     1383                arg=freq*it*Dt      !Dt est defini dans data.in ok 
     1384               ! print *,'arg:', arg  
     1385               ! print *,'A_ics:', A_ics  !acceleration 
     1386               
     1387 
     1388            ! localisation centre de la source 
     1389 
     1390              i0=floor(0.5*nxmax) 
     1391              j0=floor(0.5*nymax) 
     1392              k0=floor(0.5*nzmax) 
     1393              
     1394            ! enveloppe gaussienne 
     1395 
     1396!sigma1=4km, sigma2=100m     revoir les valeurs ci-dessous 
     1397              sigma1=20.   ! defini en unite de dx et dz 
     1398              sigma2=1.    !                    dy 
     1399 
     1400 
     1401                    do iz=0,nzmax-1 
     1402                        do iy=0,nymax-1 
     1403                            do ix=0,nxmax-1 
     1404 
     1405                                FF=exp(-((ix-i0)**2+(iy-j0)**2)/sigma1**2) 
     1406                                GG=exp(-(iz-k0)**2/sigma2**2) 
     1407 
     1408                                bx0=bx(ix,iy,iz) 
     1409                                by0=by(ix,iy,iz) 
     1410                                bz0=bz(ix,iy,iz) 
     1411 
     1412                                bx(ix,iy,iz)=bx0+A_ics*freq**2*FF*GG*sin(arg) 
     1413                                !  Acc=A_ics*FF*GG*sin(arg) 
     1414                                !print *,'acceleration:', Acc 
     1415                                by(ix,iy,iz)=by0 
     1416                                bz(ix,iy,iz)=bz0 
     1417 
     1418                            end do 
     1419                        end do 
     1420                    end do 
     1421 
     1422 
     1423 
     1424        end subroutine 
     1425 
     1426!pba : fin 
     1427 
     1428! ************************************************************************* 
     1429        subroutine urho(ax,ay,az,b,c) 
    12781430! ************************************************************************* 
    12791431!       term u.b pour l'equation de la densite (calcule dans l'espace physique) 
    12801432! ************************************************************************** 
    12811433            implicit none 
    1282  
     1434            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: c 
    12831435            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  
     1436            integer :: ix, iy, iz  
     1437            b=b+c; 
    12861438            ! VERSION LINEARISEE 
    12871439            if (lin) then 
    … …  
    13301482!       utile pour evaluer le terme non lineaire sur l'equation de la densite en economisant un champ de stockage... 
    13311483! ************************************************************************** 
    1332  
     1484        
    13331485            implicit none 
    13341486 
    … …  
    13361488            double precision ,dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bz 
    13371489            integer :: ix, iy, iz 
     1490             
     1491    !if (Rho_read) 
    13381492 
    13391493            do iz=0,nzmax-1 
    13401494                do iy=0,nymax-1 
    13411495                    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  
     1496                        bz(1,ix,iy,iz) = kx(ix)*ax(2,ix,iy,iz)+ky(iy)*ay(2,ix,iy,iz)+kz(iz)*az(2,ix,iy,iz) 
     1497                        bz(2,ix,iy,iz) = -(kx(ix)*ax(1,ix,iy,iz)+ky(iy)*ay(1,ix,iy,iz)+kz(iz)*az(1,ix,iy,iz)) 
     1498                    end do 
     1499                end do 
     1500            end do 
     1501 
     1502    !else 
     1503    !        do iz=0,nzmax-1 
     1504    !            do iy=0,nymax-1 
     1505    !                do ix=0,kxmax 
     1506    !                    bz(1,ix,iy,iz) = bz(1,ix,iy,iz) + (kx(ix)*ax(2,ix,iy,iz)+ky(iy)*ay(2,ix,iy,iz)+kz(iz)*az(2,ix,iy,iz)) 
     1507    !                    bz(2,ix,iy,iz) = bz(2,ix,iy,iz) - (kx(ix)*ax(1,ix,iy,iz)+ky(iy)*ay(1,ix,iy,iz)+kz(iz)*az(1,ix,iy,iz)) 
     1508    !                end do 
     1509    !            end do 
     1510    !        end do 
     1511  ! end if   
    13481512        end subroutine 
    13491513 
    … …  
    14111575 
    14121576! ************************************************************************* 
    1413         subroutine nonli(vox,voy,voz,rho,fx,fy,fz,frho) 
     1577        subroutine nonli(vox,voy,voz,rho,fx,fy,fz,frho,Rhov) 
    14141578! ************************************************************************* 
    14151579!       calcul des termes non linéaires f(t) et frho(t) en spectral 
    14161580!       en entrée v(t) et rho(t) en spectral 
    1417  
    1418             implicit none 
    1419  
     1581             
     1582            implicit none  
     1583 
     1584            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: Rhov 
    14201585            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho 
    14211586            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: fx, fy, fz, frho 
    1422  
     1587            
    14231588! necessaire au calcul d'une vitesse de translation 
    14241589#ifdef calcul_ut 
    … …  
    14751640 
    14761641            tps_bck_fft = tps_bck_fft - tps_temp 
    1477  
     1642             
     1643             
    14781644            call bck_fft(vox,fft_work) 
    14791645            call bck_fft(voy,fft_work) 
    … …  
    15181684            tps_vecpro2 = tps_vecpro2 - tps_temp 
    15191685            call vecpro2(vox,voy,voz,fx,fy,fz) 
     1686!pba : deb 
     1687            call forcage(fx,fy,fz) 
     1688!pba : fin 
    15201689            call cpu_time(tps_temp) 
    15211690            tps_vecpro2 = tps_vecpro2 + tps_temp 
    … …  
    15411710            if (stratification) then 
    15421711 
    1543                 ! on calcul u*rho dans l'espace physique 
     1712                ! on calcul u*rho dans l'espace physique (ou u*(rho+Rhov) si profil de densite initial) 
    15441713                call cpu_time(tps_temp) 
    15451714                tps_urho = tps_urho - tps_temp 
    1546                 call urho(vox,voy,voz,rho) 
     1715                 
     1716                            
     1717                 
     1718                call urho(vox,voy,voz,rho,Rhov) 
    15471719                call cpu_time(tps_temp) 
    15481720                tps_urho = tps_urho + tps_temp 
    15491721 
    1550                 ! on passe u*rho en spectral 
     1722                ! on passe u*rho (ou u*(rho+Rhov)) en spectral 
    15511723                call cpu_time(tps_temp) 
    15521724                tps_fwd_fft = tps_fwd_fft - tps_temp 
    … …  
    15591731                tps_fwd_fft = tps_fwd_fft + tps_temp 
    15601732 
    1561                 ! la fonction div_speciale renvoie frho = voz - u*grad(rho) 
     1733                ! la fonction div_speciale renvoie frho = voz - u*grad(rho) ou u*grad(rho+Rhov) si profil de densite initial 
     1734                ! 
    15621735                ! u*grad(rho) est calcule en prenant la divergence de u*rho parce que div(u) = 0 
    15631736                ! de plus on a ruse en stockant voz dans le champ frho pour gagner un champ en memoire 
    … …  
    15691742 
    15701743                ! ajout du terme de flottaison pour l'équation de la vitesse 
    1571                 fz=fz-xnsqr*rho 
     1744                !if (Rho_read) !on a lu un profil de densite initial 
     1745                fz=fz-9.81*rho/Rhob  
     1746                !else  
     1747                !fz=fz-xnsqr*rho  !si pas de profil de densité initial entré N est rescalé de manière à avoir dRho/dz=1          
     1748                  
    15721749            end if 
    15731750 
    … …  
    16881865 
    16891866            ! calcul des termes non linéaires f(V(t=0)) 
    1690             call nonli(vox,voy,voz,rho,fx,fy,fz,frho) 
     1867            call nonli(vox,voy,voz,rho,fx,fy,fz,frho,Rhov) 
    16911868 
    16921869            vox=vox2 
    … …  
    17171894            ! ************ 
    17181895 
    1719             call nonli(vox,voy,voz,rho,fx,fy,fz,frho) 
     1896            call nonli(vox,voy,voz,rho,fx,fy,fz,frho,Rhov) 
    17201897 
    17211898            vox=vox2 
    … …  
    17371914            ! Etape 3 / Recalcul du terme f(V(t)) qui sera utile par la suite pour le schéma d'Adams-Bashforth 
    17381915            ! ************ 
    1739             call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho) 
     1916            call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho,Rhov) 
    17401917 
    17411918            ! les champs contiennent vo = V(t+dt)  /  vo2 = V(t)  /   f = f(V(t)) 
Note: See TracChangeset for help on using the changeset viewer.