Changeset 14


Ignore:
Timestamp:
12/21/06 10:23:42 (17 years ago)
Author:
xlvlod
Message:

NS3D: modification Yannis

Location:
trunk/NS3D_JMC
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/NS3D_JMC/NS3D.F90

    r12 r14  
    345345    ! ****************************** 
    346346 
     347    if (dRho_read) then 
     348   
     349     write(*,*) "  -> lecture du fichier 'dRho_dz.init'" 
     350     call read_dRho(dRho_dz) 
     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_JMC/data.in

    r13 r14  
    2121visc_________________________   1.E-4 
    2222stratifie____________________   T 
    23     brunt_vaisala_frequency__   0.39 
     23    brunt_vaisala_frequency__   0 
    2424    schmidt number___________   1. 
    25252omega_______________________   0. 
    … …  
    4141Bruit_blanc__________________   F 
    4242    amplitude________________   0.00001 
     43 
     44*** gradient densite initiale 1D *** 
     45Lecture fichier______________   T 
     46Rhob_________________________   1025.E3 
    4347 
    4448*** forcage *** 

  • trunk/NS3D_JMC/parametres.F90

    r13 r14  
    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 :: dRho_read 
     64    double precision, save :: Rhob 
    6165 
    6266    ! Bruit blanc et amplitude du bruit blanc 
    … …  
    6569 
    6670    !introduction d'un forcage, appel de la subroutine forcage 
    67     logical, save :: forcage_sub 
     71    logical, save :: forcage_sub    
    6872    double precision, save :: A_ics, freq 
     73    integer, save :: it0   
    6974 
    7075    ! nombre et instants où sont appeles les differentes fonctions de sorties 
    … …  
    9297 
    9398    ! Vitesses, densité, termes non linéaire pour l'équation de la vitesse et de la densité 
    94     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, dRho_dz, wdRho 
    95100    double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), save :: vox2, voy2, voz2, rho2 
    96101    double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), save :: fx, fy, fz, frho 

  • trunk/NS3D_JMC/sous_routines.F90

    r13 r14  
    3131 
    3232            open (unit=1,file='data.in') 
     33 
     34 
    3335 
    3436            read(1,*) 
    … …  
    108110            read(1,*) 
    109111 
     112              
     113            !profil de densite  
     114            read(1,lecture_logical) bidon, dRho_read 
     115            read(1,lecture_reel) bidon, Rhob 
     116              
     117            read(1,*) 
     118            read(1,*) 
     119 
    110120            ! pba : deb 
    111121            ! forcage 
    … …  
    118128 
    119129            !pba fin 
    120  
    121130 
    122131            ! sorties 
    … …  
    247256            write(*,aff_logical) 'Lecture fichier..............',base2D_read 
    248257 
     258 
    249259            write(*,*) 
    250260            write(*,*) '---------------------------------------------------' 
    … …  
    255265            write(*,aff_logical) 'Bruit blanc..................',lnoise 
    256266            write(*,aff_reel) '    amplitude................',amplitude_noise 
    257  
     267             
     268            write(*,*) 
     269            write(*,*) '---------------------------------------------------' 
     270            write(*,*) 'DENSITE INITIALE 1D' 
     271            write(*,*) '---------------------------------------------------' 
     272            write(*,aff_logical) 'Lecture fichier..............',dRho_read 
     273            write(*,aff_reel) 'densite de base..............',Rhob 
     274                
     275        
    258276            write(*,*) 
    259277            write(*,*)'---------------------------------------------------' 
    … …  
    421439            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx,fy,fz 
    422440            character(len=50) :: nomfic 
    423             character(len=10) :: sufx 
     441            character(len=13) :: sufx 
    424442 
    425443            call time2string(it*dt + begin,sufx) 
    … …  
    447465 
    448466            nomfic='velo_rho_vort.'//sufx 
    449  
     467             
    450468            open (76,file=nomfic,form='unformatted') 
    451469            write(76) 1, nxmax, nymax, nzmax, lx, ly, lz, dt, truncate, type_trunc, rtrunc_x, rtrunc_y, rtrunc_z, & 
    … …  
    509527 
    510528            close(8) 
    511  
     529         
    512530        end subroutine 
    513531 
    … …  
    726744 
    727745        end subroutine 
     746 
     747 
     748! ******************************************************************************** 
     749        subroutine read_dRho(dRho_dz)  
     750!ajout d'une sous routine de lecture d'un gradient de densite moyen non constant 
     751! ******************************************************************************** 
     752!       Subroutine qui permet de lire un champ 2D en unformatted 
     753! ******************************************************************************** 
     754 
     755            implicit none 
     756 
     757            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: dRho_dz 
     758            double precision, dimension(0:nymax-1) :: dRho_dzi 
     759            integer :: flag, nxread, nyread 
     760            double precision :: lxread, lyread 
     761            integer :: ix,iy,iz 
     762            write(*,*) 'boucle lecture rho' 
     763            open (unit=88,file='dRho_dz.init',status='old') 
     764            !read(88) flag, nxread, nyread, lxread, lyread 
     765            do iy=0,nzmax-1 
     766            read(88,*) dRho_dzi(iy) 
     767            enddo 
     768             
     769            close(88) 
     770 
     771            ! on verifie que le format du fichier d'entree est correct 
     772            ! if (flag /= 1 .or. nyread /= nymax .or. lyread /= ly) then 
     773            !    write(*,*) "ERREUR: le fichier de reprise de run n'a pas le bon format" 
     774            !    write(*,*) "ARRET DU PROGRAMME..." 
     775            !    stop 
     776            !end if 
     777             
     778            do iz=0,nzmax-1 
     779                do iy=0,nymax-1 
     780                    do ix=0,nxmax-1 
     781                        dRho_dz(ix,iy,iz)=dRho_dzi(iy) 
     782                    end do 
     783                end do 
     784            end do 
     785 
     786        end subroutine 
     787 
     788 
    728789 
    729790 
    … …  
    864925 
    865926            double precision, intent(in) :: time 
    866             character(len=10), intent(out) :: string_time 
    867  
    868             write(string_time,'(A2,I4.4,A1,I3.3)') 't=',int(time),'.',nint(1000*(time-int(time))) 
     927            character(len=13), intent(out) :: string_time 
     928 
     929            write(string_time,'(A2,I7.7,A1,I3.3)') 't=',int(time),'.',nint(1000*(time-int(time))) 
    869930      end subroutine 
    870931 
    … …  
    922983            ! ************ 
    923984 
    924             call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho) 
     985            call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho,dRho_dz) 
    925986 
    926987            call cpu_time(tps_temp) 
    … …  
    9401001!       Première étape du schéma d'ordre deux d'Adams-Bashforth 
    9411002! ************************************************************************* 
     1003 
    9421004 
    9431005            implicit none 
    … …  
    12951357 
    12961358! pba: deb 
    1297 ! *************************************************************************        
     1359! ************************************************************************* 
    12981360        subroutine forcage(bx,by,bz) 
    12991361! ************************************************************************* 
    … …  
    13011363!       parametres du forcage entres ici provisoirement modifier 
    13021364!       data.in etc pour entree plus propre 
    1303 ! **************************************************************************                           
    1304                                  
    1305                                  
    1306             implicit none        
    1307                                  
     1365! ************************************************************************** 
     1366 
     1367 
     1368            implicit none 
     1369 
    13081370            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bx, by, bz 
    13091371            double precision :: vbx0,vby0,vbz0,wbx0,wby0,wbz0 
    13101372            double precision :: ax0,ay0,az0,bx0,by0,bz0 
    13111373            double precision :: time, arg, sigma1, sigma2, FF, GG, acc 
    1312             integer :: ix, iy, iz, itime, i0, j0, k0    !itime a faire passer en parametre?            
    1313             ! itime=it0 
    1314             ! itime=itime+1   ! penser a it0=it0+itmax(run precedent si restat                                 ! print *,'temps:', itime  
    1315             !    A_ics=100.             !deplacement   a definir dans  data.in et parametres.F90 
    1316             !freq=2*pi/(3600*12.4)  !M2 a verifier idem 
    1317              arg=freq*it*Dt      !Dt est defini dans data.in ok 
    1318             ! print *,'arg:', arg  
    1319             ! print *,'A_ics:', A_ics  !acceleration 
    1320                   
    1321             ! localisation centre de la source       
     1374            integer :: ix, iy, iz, itime, i0, j0, k0    !itime a faire passer en parametre? 
     1375                ! itime=it0 
     1376 
     1377                ! itime=itime+1   ! penser a it0=it0+itmax(run precedent si restat 
     1378                ! print *,'temps:', itime  
     1379                !    A_ics=100.             !deplacement   a definir dans  data.in et parametres.F90 
     1380                !freq=2*pi/(3600*12.4)  !M2 a verifier idem 
     1381                arg=freq*it*Dt      !Dt est defini dans data.in ok 
     1382               ! print *,'arg:', arg  
     1383               ! print *,'A_ics:', A_ics  !acceleration 
     1384               
     1385 
     1386            ! localisation centre de la source 
     1387 
    13221388              i0=floor(0.5*nxmax) 
    13231389              j0=floor(0.5*nymax) 
    13241390              k0=floor(0.5*nzmax) 
    1325                                  
    1326             ! enveloppe gaussienne               
    1327             !sigma1=4km, sigma2=100m     revoir les valeurs ci-dessous 
    1328             sigma1=20.   ! defini en unite de dx et dz 
    1329             sigma2=1.    !                    dy 
    1330                                       
    1331                do iz=0,nzmax-1 
    1332                    do iy=0,nymax-1 
    1333                         do ix=0,nxmax-1 
    1334                             FF=exp(-((ix-i0)**2+(iy-j0)**2)/sigma1**2)             
    1335                             GG=exp(-(iz-k0)**2/sigma2**2) 
    1336                             bx0=bx(ix,iy,iz) 
    1337                             by0=by(ix,iy,iz) 
    1338                             bz0=bz(ix,iy,iz) 
    1339                             bx(ix,iy,iz)=bx0+A_ics*freq**2*FF*GG*sin(arg) 
     1391              
     1392            ! enveloppe gaussienne 
     1393 
     1394!sigma1=4km, sigma2=100m     revoir les valeurs ci-dessous 
     1395              sigma1=20.   ! defini en unite de dx et dz 
     1396              sigma2=1.    !                    dy 
     1397 
     1398 
     1399                    do iz=0,nzmax-1 
     1400                        do iy=0,nymax-1 
     1401                            do ix=0,nxmax-1 
     1402 
     1403                                FF=exp(-((ix-i0)**2+(iy-j0)**2)/sigma1**2) 
     1404                                GG=exp(-(iz-k0)**2/sigma2**2) 
     1405 
     1406                                bx0=bx(ix,iy,iz) 
     1407                                by0=by(ix,iy,iz) 
     1408                                bz0=bz(ix,iy,iz) 
     1409 
     1410                                bx(ix,iy,iz)=bx0+A_ics*freq**2*FF*GG*sin(arg) 
    13401411                                !  Acc=A_ics*FF*GG*sin(arg) 
    13411412                                !print *,'acceleration:', Acc 
    1342                             by(ix,iy,iz)=by0 
    1343                             bz(ix,iy,iz)=bz0 
     1413                                by(ix,iy,iz)=by0 
     1414                                bz(ix,iy,iz)=bz0 
     1415 
     1416                            end do 
    13441417                        end do 
    13451418                    end do 
    1346                end do 
     1419 
     1420 
    13471421 
    13481422        end subroutine 
    13491423 
    13501424!pba : fin 
    1351  
    13521425 
    13531426! ************************************************************************* 
    … …  
    13571430! ************************************************************************** 
    13581431            implicit none 
    1359  
    13601432            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: ax, ay, az, b 
    1361             integer :: ix, iy, iz 
     1433            integer :: ix, iy, iz  
    13621434 
    13631435            ! VERSION LINEARISEE 
    … …  
    14021474 
    14031475! ************************************************************************* 
    1404         subroutine div_speciale(ax,ay,az,bz) 
    1405 ! ************************************************************************* 
    1406 !       Cette sous-routine renvoie dans le champ bz le terme bz - div(ax,ay,az) 
     1476        subroutine vozdRho_dz(a,b,c) 
     1477! ************************************************************************* 
     1478!       term c=a.b pour l'equation de la densite (avec a=voz  b=dRho_dz) (calcule dans l'espace physique) 
     1479! ************************************************************************** 
     1480            implicit none 
     1481            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: a, b 
     1482            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: c  
     1483            integer :: ix, iy, iz  
     1484 
     1485            
     1486                    do iz=0,nzmax-1 
     1487                        do iy=0,nymax-1 
     1488                            do ix=0,nxmax-1 
     1489                                c(ix,iy,iz)=b(ix,iy,iz)*a(ix,iy,iz) 
     1490                                 
     1491                            end do 
     1492                        end do 
     1493                    end do 
     1494            
     1495 
     1496         end subroutine 
     1497 
     1498 
     1499 
     1500 
     1501! ************************************************************************* 
     1502        subroutine div_speciale(ax,ay,az,bz1,bz2) 
     1503! ************************************************************************* 
     1504!       Cette sous-routine renvoie dans le champ bz2 le terme bz1 - div(ax,ay,az) 
    14071505!       utile pour evaluer le terme non lineaire sur l'equation de la densite en economisant un champ de stockage... 
    14081506! ************************************************************************** 
    1409  
    1410             implicit none 
    1411  
    1412             double precision ,dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: ax, ay, az 
    1413             double precision ,dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bz 
     1507        
     1508            implicit none 
     1509 
     1510            double precision ,dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: ax, ay, az, bz1 
     1511            double precision ,dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bz2 
    14141512            integer :: ix, iy, iz 
     1513             
     1514  
     1515                         
    14151516 
    14161517            do iz=0,nzmax-1 
    14171518                do iy=0,nymax-1 
    14181519                    do ix=0,kxmax 
    1419                         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)) 
    1420                         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)) 
    1421                     end do 
    1422                 end do 
    1423             end do 
    1424  
     1520                       bz2(1,ix,iy,iz) = - bz1(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)) 
     1521                       bz2(2,ix,iy,iz) = - bz1(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)) 
     1522                    end do 
     1523                end do 
     1524            end do 
     1525    
    14251526        end subroutine 
    14261527 
    … …  
    14881589 
    14891590! ************************************************************************* 
    1490         subroutine nonli(vox,voy,voz,rho,fx,fy,fz,frho) 
     1591        subroutine nonli(vox,voy,voz,rho,fx,fy,fz,frho,dRho_dz) 
    14911592! ************************************************************************* 
    14921593!       calcul des termes non linéaires f(t) et frho(t) en spectral 
    14931594!       en entrée v(t) et rho(t) en spectral 
    1494  
    1495             implicit none 
    1496  
     1595             
     1596            implicit none  
     1597 
     1598            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: dRho_dz 
    14971599            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho 
    14981600            double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: fx, fy, fz, frho 
    1499  
     1601            
    15001602! necessaire au calcul d'une vitesse de translation 
    15011603#ifdef calcul_ut 
    … …  
    15061608 
    15071609            ! sauvegarde de voz en sphérique dans le champ frho (ruse pour gagner de la place en memoire) 
     1610            !(cette ruse ne marche plus si on rentre un dRho_dz) 
    15081611            ! nécessaire pour le cas stratifié 
    1509             if (stratification) then 
    1510                 call cpu_time(tps_temp) 
    1511                 tps_copie = tps_copie - tps_temp 
    1512                 frho=voz 
    1513                 call cpu_time(tps_temp) 
    1514                 tps_copie = tps_copie + tps_temp 
    1515             end if 
     1612            !if (stratification) then 
     1613            !    call cpu_time(tps_temp) 
     1614            !    tps_copie = tps_copie - tps_temp 
     1615            !    frho=voz 
     1616            !    call cpu_time(tps_temp) 
     1617            !    tps_copie = tps_copie + tps_temp 
     1618            !end if 
    15161619 
    15171620            ! calcul de la vorticite de vo. Le resultat est placé dans f. 
    … …  
    15521655 
    15531656            tps_bck_fft = tps_bck_fft - tps_temp 
    1554  
     1657             
     1658             
    15551659            call bck_fft(vox,fft_work) 
    15561660            call bck_fft(voy,fft_work) 
    … …  
    15971701!pba : deb 
    15981702            call forcage(fx,fy,fz) 
    1599 !pba : fin                                                             
     1703!pba : fin 
    16001704            call cpu_time(tps_temp) 
    16011705            tps_vecpro2 = tps_vecpro2 + tps_temp 
    … …  
    16241728                call cpu_time(tps_temp) 
    16251729                tps_urho = tps_urho - tps_temp 
     1730                 
     1731                            
     1732                call vozdRho_dz(voz,dRho_dz,wdRho) 
    16261733                call urho(vox,voy,voz,rho) 
     1734                  
    16271735                call cpu_time(tps_temp) 
    16281736                tps_urho = tps_urho + tps_temp 
    16291737 
    1630                 ! on passe u*rho en spectral 
     1738                ! on passe u*rho et w*dRho_dz en spectral 
    16311739                call cpu_time(tps_temp) 
    16321740                tps_fwd_fft = tps_fwd_fft - tps_temp 
    … …  
    16351743                call fwd_fft(voy,fft_work) 
    16361744                call fwd_fft(voz,fft_work) 
     1745                call fwd_fft(wdRho,fft_work) 
    16371746 
    16381747                call cpu_time(tps_temp) 
    16391748                tps_fwd_fft = tps_fwd_fft + tps_temp 
    16401749 
    1641                 ! la fonction div_speciale renvoie frho = voz - u*grad(rho) 
     1750                ! la fonction div_speciale renvoie frho = voz*dRho_dz - u*grad(rho) (drho_dz= gradient de densité moyen rentré  au départ)  
     1751                ! 
    16421752                ! u*grad(rho) est calcule en prenant la divergence de u*rho parce que div(u) = 0 
    16431753                ! de plus on a ruse en stockant voz dans le champ frho pour gagner un champ en memoire 
    16441754                call cpu_time(tps_temp) 
    16451755                tps_div_speciale = tps_div_speciale - tps_temp 
    1646                 call div_speciale(vox,voy,voz,frho) 
     1756                call div_speciale(vox,voy,voz,wdRho,frho) 
    16471757                call cpu_time(tps_temp) 
    16481758                tps_div_speciale = tps_div_speciale + tps_temp 
    16491759 
    16501760                ! ajout du terme de flottaison pour l'équation de la vitesse 
    1651                 fz=fz-xnsqr*rho 
     1761                !if (Rho_read) !on a lu un profil de densite initial 
     1762                fz=fz-9.81*rho/Rhob  
     1763                !else  
     1764                !fz=fz-xnsqr*rho  !si pas de profil de densité initial entré N est rescalé de manière à avoir dRho/dz=1          
     1765                  
    16521766            end if 
    16531767 
    … …  
    17681882 
    17691883            ! calcul des termes non linéaires f(V(t=0)) 
    1770             call nonli(vox,voy,voz,rho,fx,fy,fz,frho) 
     1884            call nonli(vox,voy,voz,rho,fx,fy,fz,frho,dRho_dz) 
    17711885 
    17721886            vox=vox2 
    … …  
    17971911            ! ************ 
    17981912 
    1799             call nonli(vox,voy,voz,rho,fx,fy,fz,frho) 
     1913            call nonli(vox,voy,voz,rho,fx,fy,fz,frho,dRho_dz) 
    18001914 
    18011915            vox=vox2 
    … …  
    18171931            ! Etape 3 / Recalcul du terme f(V(t)) qui sera utile par la suite pour le schéma d'Adams-Bashforth 
    18181932            ! ************ 
    1819             call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho) 
     1933            call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho,dRho_dz) 
    18201934 
    18211935            ! les champs contiennent vo = V(t+dt)  /  vo2 = V(t)  /   f = f(V(t)) 
Note: See TracChangeset for help on using the changeset viewer.