- Timestamp:
- 12/05/06 10:39:48 (17 years ago)
- Location:
- trunk/NS3D/SRC
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NS3D/SRC/NS3D.F90
r8 r9 345 345 ! ****************************** 346 346 347 if (Rho_read) then 348 349 write(*,*) " -> lecture du fichier 'Rho.init'" 350 call readRho(Rhov) 351 352 end if 353 354 347 355 if (reprise_run) then 348 356 write(*,*) "Reprise d'un run" … … 466 474 ! Globalement ici, on rentre donc avec les valeurs à l'instant 1 et on en sort avec les valeurs à l'instant itmax 467 475 do it=2,itmax 468 476 write(*,*) 'boucle it2, itmax' 469 477 ! Shéma d'ordre deux d'Adams-Bashforth 470 478 call adams_bashforth(vox,voy,voz,rho,fx,fy,fz,frho,vox2,voy2,voz2,rho2) -
trunk/NS3D/SRC/config.h
r8 r9 16 16 #define DIMX 256 17 17 #define DIMY 256 18 #define DIMZ 25618 #define DIMZ 3 19 19 20 20 ! utile pour definir un mode "debugage" dans le code (pas d'initialisation du generateur de nombre aleatoire, … … 39 39 ! de l'architecture de calcul (NEC, processeur Intel...). 40 40 #define PPADX 0 41 #define PPADY 142 #define PPADZ 25641 #define PPADY 0 42 #define PPADZ 0 43 43 44 44 ! choix de la fonction qui genere l'ecoulement de base. Actuellement est disponible : -
trunk/NS3D/SRC/data.in
r8 r9 6 6 ny___________________________ 256 7 7 nz___________________________ 3 8 Lx___________________________ 14.13229 Ly___________________________ 60.8 Lx___________________________ 3.E4 9 Ly___________________________ 25.E2 10 10 Lz___________________________ 10. 11 Dt___________________________ 0.0111 Dt___________________________ 10. 12 12 begin________________________ 0. 13 itmax________________________ 1500013 itmax________________________ 50000 14 14 Troncature_oui/non___________ T 15 15 carree:1 / elliptique:2__ 1 … … 19 19 20 20 *** Parametres physiques *** 21 visc_________________________ 5.E- 622 stratifie____________________ F23 brunt_vaisala_frequency__ 10.21 visc_________________________ 5.E-3 22 stratifie____________________ T 23 brunt_vaisala_frequency__ 5.E-3 24 24 schmidt number___________ 1. 25 25 2omega_______________________ 0. 26 26 27 27 *** Type de simulation *** 28 En perturbation______________ T29 Lineaire_____________________ T28 En perturbation______________ F 29 Lineaire_____________________ F 30 30 31 31 *** Reprise de run *** … … 33 33 34 34 *** Etat de base (pour simulation en perturbation) *** 35 Call subroutine gen_base2D___ T35 Call subroutine gen_base2D___ F 36 36 Lecture fichier______________ F 37 37 … … 39 39 Call subroutine gen_velo2D___ F 40 40 Lecture fichier______________ F 41 Bruit_blanc__________________ T41 Bruit_blanc__________________ F 42 42 amplitude________________ 0.00001 43 43 44 *** gradient densite initiale 1D *** 45 Lecture fichier______________ T 46 Rhob_________________________ 1.E3 47 48 *** forcage *** 49 Call subroutine forcage______ T 50 A_ics____________________ 15E2 51 freq_____________________ 15.E-5 52 it0____________________ 0 53 44 54 *** Sorties *** 45 periode out_physique1________ 1 46 periode out_physique2________ 200047 periode out_spectral1________ 100 55 periode out_physique1________ 1000 56 periode out_physique2________ 1000 57 periode out_spectral1________ 1000 48 58 49 59 *** Parametres additionnels... *** … … 61 71 Circulation______________ -6.28 62 72 Rayon____________________ 1. 73 -
trunk/NS3D/SRC/parametres.F90
r8 r9 59 59 ! lecture de la vitesse initiale 2D en "unformatted" ou appel subroutine 60 60 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 61 65 62 66 ! Bruit blanc et amplitude du bruit blanc 63 67 logical, save :: lnoise 64 68 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 65 74 66 75 ! nombre et instants où sont appeles les differentes fonctions de sorties … … 88 97 89 98 ! 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 91 100 double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), save :: vox2, voy2, voz2, rho2 92 101 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 31 31 32 32 open (unit=1,file='data.in') 33 34 33 35 34 36 read(1,*) … … 107 109 read(1,*) 108 110 read(1,*) 111 112 113 !profil de densite 114 read(1,lecture_logical) bidon, Rho_read 115 read(1,lecture_reel) bidon, Rhob 116 117 read(1,*) 118 read(1,*) 119 120 ! pba : deb 121 ! forcage 122 read(1,lecture_logical) bidon, forcage_sub 123 read(1,lecture_reel) bidon, A_ics 124 read(1,lecture_reel) bidon, freq 125 read(1,lecture_entier) bidon, it0 126 127 read(1,*) 128 read(1,*) 129 130 !pba fin 109 131 110 132 ! sorties … … 235 257 write(*,aff_logical) 'Lecture fichier..............',base2D_read 236 258 259 237 260 write(*,*) 238 261 write(*,*) '---------------------------------------------------' … … 243 266 write(*,aff_logical) 'Bruit blanc..................',lnoise 244 267 write(*,aff_reel) ' amplitude................',amplitude_noise 268 269 write(*,*) 270 write(*,*) '---------------------------------------------------' 271 write(*,*) 'DENSITE INITIALE 1D' 272 write(*,*) '---------------------------------------------------' 273 write(*,aff_logical) 'Lecture fichier..............',Rho_read 274 write(*,aff_reel) 'densite de base..............',Rhob 275 276 277 write(*,*) 278 write(*,*)'---------------------------------------------------' 279 write(*,*)'FORCAGE EN VITESSE' 280 write(*,*)'---------------------------------------------------' 281 write(*,aff_logical)'call subroutine', forcage_sub 282 write(*,aff_reel)'Amplitude..................',A_ics 283 write(*,aff_reel)'Frequence..................',freq 245 284 246 285 write(*,*) … … 401 440 double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx,fy,fz 402 441 character(len=50) :: nomfic 403 character(len=1 0) :: sufx442 character(len=13) :: sufx 404 443 405 444 call time2string(it*dt + begin,sufx) … … 427 466 428 467 nomfic='velo_rho_vort.'//sufx 429 468 430 469 open (76,file=nomfic,form='unformatted') 431 470 write(76) 1, nxmax, nymax, nzmax, lx, ly, lz, dt, truncate, type_trunc, rtrunc_x, rtrunc_y, rtrunc_z, & … … 443 482 444 483 close(76) 445 484 print *,'itfic:', it 446 485 end subroutine 447 486 … … 489 528 490 529 close(8) 491 530 492 531 end subroutine 493 532 … … 706 745 707 746 end subroutine 747 748 749 ! ******************************************************************************** 750 subroutine readRho(Rhov) 751 !ajout d'une sous routine de lecture d'un gradient de densite moyen non constant 752 ! ******************************************************************************** 753 ! Subroutine qui permet de lire un champ 2D en unformatted 754 ! ******************************************************************************** 755 756 implicit none 757 758 double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: Rhov 759 double precision, dimension(0:nymax-1) :: Rhoi 760 integer :: flag, nxread, nyread 761 double precision :: lxread, lyread 762 integer :: ix,iy,iz 763 write(*,*) 'boucle lecture rho' 764 open (unit=88,file='Rho.init',status='old') 765 !read(88) flag, nxread, nyread, lxread, lyread 766 do iy=0,255 767 read(88,*) Rhoi(iy) 768 enddo 769 770 close(88) 771 772 ! on verifie que le format du fichier d'entree est correct 773 ! if (flag /= 1 .or. nyread /= nymax .or. lyread /= ly) then 774 ! write(*,*) "ERREUR: le fichier de reprise de run n'a pas le bon format" 775 ! write(*,*) "ARRET DU PROGRAMME..." 776 ! stop 777 !end if 778 779 do iz=0,nzmax-1 780 do iy=0,nymax-1 781 do ix=0,nxmax-1 782 Rhov(ix,iy,iz)=Rhoi(iy) 783 end do 784 end do 785 end do 786 print*, Rhov 787 788 end subroutine 789 790 708 791 709 792 … … 844 927 845 928 double precision, intent(in) :: time 846 character(len=1 0), intent(out) :: string_time847 848 write(string_time,'(A2,I 4.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))) 849 932 end subroutine 850 933 … … 902 985 ! ************ 903 986 904 call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho )987 call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho,Rhov) 905 988 906 989 call cpu_time(tps_temp) … … 920 1003 ! Première étape du schéma d'ordre deux d'Adams-Bashforth 921 1004 ! ************************************************************************* 1005 922 1006 923 1007 implicit none … … 1274 1358 end subroutine 1275 1359 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) 1278 1430 ! ************************************************************************* 1279 1431 ! term u.b pour l'equation de la densite (calcule dans l'espace physique) 1280 1432 ! ************************************************************************** 1281 1433 implicit none 1282 1434 double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: c 1283 1435 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; 1286 1438 ! VERSION LINEARISEE 1287 1439 if (lin) then … … 1330 1482 ! utile pour evaluer le terme non lineaire sur l'equation de la densite en economisant un champ de stockage... 1331 1483 ! ************************************************************************** 1332 1484 1333 1485 implicit none 1334 1486 … … 1336 1488 double precision ,dimension(2,0:dimx_tab/2-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bz 1337 1489 integer :: ix, iy, iz 1490 1491 !if (Rho_read) 1338 1492 1339 1493 do iz=0,nzmax-1 1340 1494 do iy=0,nymax-1 1341 1495 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 1348 1512 end subroutine 1349 1513 … … 1411 1575 1412 1576 ! ************************************************************************* 1413 subroutine nonli(vox,voy,voz,rho,fx,fy,fz,frho )1577 subroutine nonli(vox,voy,voz,rho,fx,fy,fz,frho,Rhov) 1414 1578 ! ************************************************************************* 1415 1579 ! calcul des termes non linéaires f(t) et frho(t) en spectral 1416 1580 ! 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 1420 1585 double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho 1421 1586 double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: fx, fy, fz, frho 1422 1587 1423 1588 ! necessaire au calcul d'une vitesse de translation 1424 1589 #ifdef calcul_ut … … 1475 1640 1476 1641 tps_bck_fft = tps_bck_fft - tps_temp 1477 1642 1643 1478 1644 call bck_fft(vox,fft_work) 1479 1645 call bck_fft(voy,fft_work) … … 1518 1684 tps_vecpro2 = tps_vecpro2 - tps_temp 1519 1685 call vecpro2(vox,voy,voz,fx,fy,fz) 1686 !pba : deb 1687 call forcage(fx,fy,fz) 1688 !pba : fin 1520 1689 call cpu_time(tps_temp) 1521 1690 tps_vecpro2 = tps_vecpro2 + tps_temp … … 1541 1710 if (stratification) then 1542 1711 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) 1544 1713 call cpu_time(tps_temp) 1545 1714 tps_urho = tps_urho - tps_temp 1546 call urho(vox,voy,voz,rho) 1715 1716 1717 1718 call urho(vox,voy,voz,rho,Rhov) 1547 1719 call cpu_time(tps_temp) 1548 1720 tps_urho = tps_urho + tps_temp 1549 1721 1550 ! on passe u*rho en spectral1722 ! on passe u*rho (ou u*(rho+Rhov)) en spectral 1551 1723 call cpu_time(tps_temp) 1552 1724 tps_fwd_fft = tps_fwd_fft - tps_temp … … 1559 1731 tps_fwd_fft = tps_fwd_fft + tps_temp 1560 1732 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 ! 1562 1735 ! u*grad(rho) est calcule en prenant la divergence de u*rho parce que div(u) = 0 1563 1736 ! de plus on a ruse en stockant voz dans le champ frho pour gagner un champ en memoire … … 1569 1742 1570 1743 ! 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 1572 1749 end if 1573 1750 … … 1688 1865 1689 1866 ! 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) 1691 1868 1692 1869 vox=vox2 … … 1717 1894 ! ************ 1718 1895 1719 call nonli(vox,voy,voz,rho,fx,fy,fz,frho )1896 call nonli(vox,voy,voz,rho,fx,fy,fz,frho,Rhov) 1720 1897 1721 1898 vox=vox2 … … 1737 1914 ! Etape 3 / Recalcul du terme f(V(t)) qui sera utile par la suite pour le schéma d'Adams-Bashforth 1738 1915 ! ************ 1739 call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho )1916 call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho,Rhov) 1740 1917 1741 1918 ! les champs contiennent vo = V(t+dt) / vo2 = V(t) / f = f(V(t))
Note: See TracChangeset
for help on using the changeset viewer.