- Timestamp:
- 12/21/06 10:23:42 (17 years ago)
- Location:
- trunk/NS3D_JMC
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NS3D_JMC/NS3D.F90
r12 r14 345 345 ! ****************************** 346 346 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 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_JMC/data.in
r13 r14 21 21 visc_________________________ 1.E-4 22 22 stratifie____________________ T 23 brunt_vaisala_frequency__ 0 .3923 brunt_vaisala_frequency__ 0 24 24 schmidt number___________ 1. 25 25 2omega_______________________ 0. … … 41 41 Bruit_blanc__________________ F 42 42 amplitude________________ 0.00001 43 44 *** gradient densite initiale 1D *** 45 Lecture fichier______________ T 46 Rhob_________________________ 1025.E3 43 47 44 48 *** forcage *** -
trunk/NS3D_JMC/parametres.F90
r13 r14 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 :: dRho_read 64 double precision, save :: Rhob 61 65 62 66 ! Bruit blanc et amplitude du bruit blanc … … 65 69 66 70 !introduction d'un forcage, appel de la subroutine forcage 67 logical, save :: forcage_sub 71 logical, save :: forcage_sub 68 72 double precision, save :: A_ics, freq 73 integer, save :: it0 69 74 70 75 ! nombre et instants où sont appeles les differentes fonctions de sorties … … 92 97 93 98 ! 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 95 100 double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), save :: vox2, voy2, voz2, rho2 96 101 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 31 31 32 32 open (unit=1,file='data.in') 33 34 33 35 34 36 read(1,*) … … 108 110 read(1,*) 109 111 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 110 120 ! pba : deb 111 121 ! forcage … … 118 128 119 129 !pba fin 120 121 130 122 131 ! sorties … … 247 256 write(*,aff_logical) 'Lecture fichier..............',base2D_read 248 257 258 249 259 write(*,*) 250 260 write(*,*) '---------------------------------------------------' … … 255 265 write(*,aff_logical) 'Bruit blanc..................',lnoise 256 266 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 258 276 write(*,*) 259 277 write(*,*)'---------------------------------------------------' … … 421 439 double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(in) :: fx,fy,fz 422 440 character(len=50) :: nomfic 423 character(len=1 0) :: sufx441 character(len=13) :: sufx 424 442 425 443 call time2string(it*dt + begin,sufx) … … 447 465 448 466 nomfic='velo_rho_vort.'//sufx 449 467 450 468 open (76,file=nomfic,form='unformatted') 451 469 write(76) 1, nxmax, nymax, nzmax, lx, ly, lz, dt, truncate, type_trunc, rtrunc_x, rtrunc_y, rtrunc_z, & … … 509 527 510 528 close(8) 511 529 512 530 end subroutine 513 531 … … 726 744 727 745 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 728 789 729 790 … … 864 925 865 926 double precision, intent(in) :: time 866 character(len=1 0), intent(out) :: string_time867 868 write(string_time,'(A2,I 4.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))) 869 930 end subroutine 870 931 … … 922 983 ! ************ 923 984 924 call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho )985 call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho,dRho_dz) 925 986 926 987 call cpu_time(tps_temp) … … 940 1001 ! Première étape du schéma d'ordre deux d'Adams-Bashforth 941 1002 ! ************************************************************************* 1003 942 1004 943 1005 implicit none … … 1295 1357 1296 1358 ! pba: deb 1297 ! ************************************************************************* 1359 ! ************************************************************************* 1298 1360 subroutine forcage(bx,by,bz) 1299 1361 ! ************************************************************************* … … 1301 1363 ! parametres du forcage entres ici provisoirement modifier 1302 1364 ! data.in etc pour entree plus propre 1303 ! ************************************************************************** 1304 1305 1306 implicit none 1307 1365 ! ************************************************************************** 1366 1367 1368 implicit none 1369 1308 1370 double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: bx, by, bz 1309 1371 double precision :: vbx0,vby0,vbz0,wbx0,wby0,wbz0 1310 1372 double precision :: ax0,ay0,az0,bx0,by0,bz0 1311 1373 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 1322 1388 i0=floor(0.5*nxmax) 1323 1389 j0=floor(0.5*nymax) 1324 1390 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) 1340 1411 ! Acc=A_ics*FF*GG*sin(arg) 1341 1412 !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 1344 1417 end do 1345 1418 end do 1346 end do 1419 1420 1347 1421 1348 1422 end subroutine 1349 1423 1350 1424 !pba : fin 1351 1352 1425 1353 1426 ! ************************************************************************* … … 1357 1430 ! ************************************************************************** 1358 1431 implicit none 1359 1360 1432 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 1362 1434 1363 1435 ! VERSION LINEARISEE … … 1402 1474 1403 1475 ! ************************************************************************* 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) 1407 1505 ! utile pour evaluer le terme non lineaire sur l'equation de la densite en economisant un champ de stockage... 1408 1506 ! ************************************************************************** 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 1414 1512 integer :: ix, iy, iz 1513 1514 1515 1415 1516 1416 1517 do iz=0,nzmax-1 1417 1518 do iy=0,nymax-1 1418 1519 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 1425 1526 end subroutine 1426 1527 … … 1488 1589 1489 1590 ! ************************************************************************* 1490 subroutine nonli(vox,voy,voz,rho,fx,fy,fz,frho )1591 subroutine nonli(vox,voy,voz,rho,fx,fy,fz,frho,dRho_dz) 1491 1592 ! ************************************************************************* 1492 1593 ! calcul des termes non linéaires f(t) et frho(t) en spectral 1493 1594 ! 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 1497 1599 double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(inout) :: vox, voy, voz, rho 1498 1600 double precision, dimension(0:dimx_tab-1,0:dimy_tab-1,0:dimz_tab-1), intent(out) :: fx, fy, fz, frho 1499 1601 1500 1602 ! necessaire au calcul d'une vitesse de translation 1501 1603 #ifdef calcul_ut … … 1506 1608 1507 1609 ! 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) 1508 1611 ! nécessaire pour le cas stratifié 1509 if (stratification) then1510 call cpu_time(tps_temp)1511 tps_copie = tps_copie - tps_temp1512 frho=voz1513 call cpu_time(tps_temp)1514 tps_copie = tps_copie + tps_temp1515 end if1612 !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 1516 1619 1517 1620 ! calcul de la vorticite de vo. Le resultat est placé dans f. … … 1552 1655 1553 1656 tps_bck_fft = tps_bck_fft - tps_temp 1554 1657 1658 1555 1659 call bck_fft(vox,fft_work) 1556 1660 call bck_fft(voy,fft_work) … … 1597 1701 !pba : deb 1598 1702 call forcage(fx,fy,fz) 1599 !pba : fin 1703 !pba : fin 1600 1704 call cpu_time(tps_temp) 1601 1705 tps_vecpro2 = tps_vecpro2 + tps_temp … … 1624 1728 call cpu_time(tps_temp) 1625 1729 tps_urho = tps_urho - tps_temp 1730 1731 1732 call vozdRho_dz(voz,dRho_dz,wdRho) 1626 1733 call urho(vox,voy,voz,rho) 1734 1627 1735 call cpu_time(tps_temp) 1628 1736 tps_urho = tps_urho + tps_temp 1629 1737 1630 ! on passe u*rho e n spectral1738 ! on passe u*rho et w*dRho_dz en spectral 1631 1739 call cpu_time(tps_temp) 1632 1740 tps_fwd_fft = tps_fwd_fft - tps_temp … … 1635 1743 call fwd_fft(voy,fft_work) 1636 1744 call fwd_fft(voz,fft_work) 1745 call fwd_fft(wdRho,fft_work) 1637 1746 1638 1747 call cpu_time(tps_temp) 1639 1748 tps_fwd_fft = tps_fwd_fft + tps_temp 1640 1749 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 ! 1642 1752 ! u*grad(rho) est calcule en prenant la divergence de u*rho parce que div(u) = 0 1643 1753 ! de plus on a ruse en stockant voz dans le champ frho pour gagner un champ en memoire 1644 1754 call cpu_time(tps_temp) 1645 1755 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) 1647 1757 call cpu_time(tps_temp) 1648 1758 tps_div_speciale = tps_div_speciale + tps_temp 1649 1759 1650 1760 ! 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 1652 1766 end if 1653 1767 … … 1768 1882 1769 1883 ! 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) 1771 1885 1772 1886 vox=vox2 … … 1797 1911 ! ************ 1798 1912 1799 call nonli(vox,voy,voz,rho,fx,fy,fz,frho )1913 call nonli(vox,voy,voz,rho,fx,fy,fz,frho,dRho_dz) 1800 1914 1801 1915 vox=vox2 … … 1817 1931 ! Etape 3 / Recalcul du terme f(V(t)) qui sera utile par la suite pour le schéma d'Adams-Bashforth 1818 1932 ! ************ 1819 call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho )1933 call nonli(vox2,voy2,voz2,rho2,fx,fy,fz,frho,dRho_dz) 1820 1934 1821 1935 ! les champs contiennent vo = V(t+dt) / vo2 = V(t) / f = f(V(t))
Note: See TracChangeset
for help on using the changeset viewer.