Bug avec un code de prolongement vers le haut

[Flo Flo]

Bonjour à tous

Je me permets de solliciter votre aide parce que j'ai un souci avec un programme en Fortran que j'ai essayé de faire en me basant sur un autre programme existant.

Je souhaite un code qui me calcule un prolongement vers le haut d'un champ potentiel (En utilisant des transformées de Fourier). Le code que j'ai écrit est le suivant:

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
      program upward

c Déclarations

      character name*70
      real rft(512,512), flo(512,512), fla(512,512)

    1 format (a)
    2 format (f11.6,x,f10.6,x,f10.2)

      pi = 3.1415926535
      tpi = 2. * pi
      rad = pi / 180.

      dimension grid(nx*ny),store(2*nx*ny),nn(2)      
      real kx, ky, k
      complex cgrid, cmplx
      data pi/3.14159265/
      index(i,j,ncol)=(j-1)*ncol+i
      nn(1) = ny
      nn(2) = nx
      ndim = 2
      dkx = 2.*pi/(nx*dx)
      dky = 2.*pi/(ny*dy)

c Lecture des valeurs

      print*,'Nom de la grille d''anomalies magnetiques?'
      read (5,1) name
      print*,'Nombre de colonnes de la grille, en 2**n?'
      read (5,*) ny
      ny = 2 ** ny
      print*,'Nombre de lignes de la grille, en 2**n?'
      read (5,*) nx
      nx = 2 ** nx
      print*,'Pas d''echantillonnage en x?'
      read (5,*) dx
      print*,'Pas d''échantillonnage en y?'
      read (5,*) dy
      print*,'Hauteur du prolongement vers le haut?'
      read (5,*) dz
      
      store = 2. * nx * ny

c Programme

      eu = tpi / float(nx)
      ev = tpi / float(ny)
      nx2 = nx / 2 + 1
      ny2 = ny / 2 + 1

      open (70, file = name)
      do 3 i = 1, nx
      do 3 j = 1, ny
      read (70,*) flo(nx-i+1,j), fla(nx-i+1,j), rft(nx-i+1,j)
    3 continue
      close (70)

      do 4 i = 1, nx
      do 4 j = 1, ny

    4 continue

      do 10 j=1,nx
      do 10 i=1,ny
      ij=index(i,j,ny)
      store(2*ij-1)=grid(ij)
   10 store(2*ij)=0.
      call fourn(store,nn,ndim,-1)
      do 20 j=1,nx
      do 20 i=1,ny
      ij=index(i,j,ny)
      call kvalue(i,j,nx,ny,dkx,dky,kx,ky)
      k=sqrt(kx**2+ky**2)
      cont=exp(-k*dz)
      cgrid=cmplx(store(2*ij-1),store(2*ij))*cont
      store(2*ij-1)=real(cgrid)
   20 store(2*ij)=aimag(cgrid)
      endif
    5 continue

      call fourn(store,nn,ndim,+1)
      do 30 j=1,nx
      do 30 i=1,ny
      ij=index(i,j,ny)
   30 grid(ij)=store(2*ij-1)/(nx*ny)

      print*,'Nom de la grille prolongée vers le haut?'
      read (5,1) name
      open (70, file = name)
      do 6 i = 1, nx
      do 6 j = 1, ny
      rft(i,j) = rft(i,j) / acof * 2.
    6 continue
      do 7 i = 1, nx
      do 7 j = 1, ny
      write (70,2) flo(nx-i+1,j), fla(nx-i+1,j), rft(nx-i+1,j)
    7 continue
      close (70)

      stop
      end

Bien entendu, comme tout programme qui se respecte, il bugge

Je précise que je suis complètement nul en programmation. J'ai essayé de faire ce que je pouvais mais il continue de m'afficher toute une gamme de messages d'erreurs quand je le compile et je ne comprends pas forcément ce que je dois faire pour que ça fonctionne.

Je vous mets les subroutines avec lesquelles il fonctionne:

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
subroutine fourn(data,nn,ndim,isign)
c
c  Replaces DATA by its NDIM-dimensional discrete Fourier transform, 
c  if ISIGN is input as 1.  NN is an integer array of length NDIM, 
c  containing the lengths of each dimension (number of complex values), 
c  which must all be powers of 2.  DATA is a real array of length twice 
c  the product of these lengths, in which the data are stored as in a 
c  multidimensional complex Fortran array.  If ISIGN is input as -1, 
c  DATA is replaced by its inverse transform times the product of the
c  lengths of all dimensions.  From Press, W.H., Flannery, B.P., 
c  Teukolsky, S.A., and Vetterling, W.T., 1986, Numerical Recipes, 
c  Cambridge Univ. Press, p. 451-453.
c
      real*8 wr,wi,wpr,wpi,wtemp,theta
      dimension nn(ndim),data(*)
      ntot=1
      do 11 iidim=1,ndim
      ntot=ntot*nn(iidim)
   11    continue
      nprev=1
      do 18 iidim=1,ndim
        n=nn(iidim)
        nrem=ntot/(n*nprev)
        ip1=2*nprev
        ip2=ip1*n
        ip3=ip2*nrem
        i2rev=1
        do 14 i2=1,ip2,ip1
          if(i2.lt.i2rev)then
            do 13 i1=i2,i2+ip1-2,2
              do 12 i3=i1,ip3,ip2
                i3rev=i2rev+i3-i2
                tempr=data(i3)
                tempi=data(i3+1)
                data(i3)=data(i3rev)
                data(i3+1)=data(i3rev+1)
                data(i3rev)=tempr
                data(i3rev+1)=tempi
   12         continue
   13       continue
          endif
          ibit=ip2/2
   1      if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then
            i2rev=i2rev-ibit
            ibit=ibit/2
          go to 1
          endif
          i2rev=i2rev+ibit
   14   continue
        ifp1=ip1
   2    if(ifp1.lt.ip2)then
          ifp2=2*ifp1
          theta=isign*6.28318530717959d0/(ifp2/ip1)
          wpr=-2.d0*dsin(0.5d0*theta)**2
          wpi=dsin(theta)
          wr=1.d0
          wi=0.d0
          do 17 i3=1,ifp1,ip1
            do 16 i1=i3,i3+ip1-2,2
              do 15 i2=i1,ip3,ifp2
                k1=i2
                k2=k1+ifp1
                tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1)
                tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2)
                data(k2)=data(k1)-tempr
                data(k2+1)=data(k1+1)-tempi
                data(k1)=data(k1)+tempr
                data(k1+1)=data(k1+1)+tempi
   15         continue
   16       continue
            wtemp=wr
            wr=wr*wpr-wi*wpi+wr
            wi=wi*wpr+wtemp*wpi+wi
   17     continue
          ifp1=ifp2
        go to 2
        endif
        nprev=n*nprev
   18 continue
      return
      end

Et

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
subroutine kvalue(i,j,nx,ny,dkx,dky,kx,ky)
c  Subroutine KVALUE finds the wavenumber coordinates of one 
c  element of a rectangular grid from subroutine FOURN.  
c
c  Input parameters:
c    i  - index in the ky direction.
c    j  - index in the kx direction.
c    nx - dimension of grid in ky direction (a power of two).
c    ny - dimension of grid in kx direction (a power of two).
c    dkx - sample interval in the kx direction.
c    dky - sample interval in the ky direction.
c
c  Output parameters:
c    kx - the wavenumber coordinate in the kx direction.
c    ky - the wavenumber coordinate in the ky direction.
c
      real kx,ky
      nyqx=nx/2+1
      nyqy=ny/2+1
      if(j.le.nyqx)then
         kx=(j-1)*dkx
         else
            kx=(j-nx-1)*dkx
            end if
      if(i.le.nyqy)then
        ky=(i-1)*dky
        else
            ky=(i-ny-1)*dky
            end if
      return
      end

Voilà voilà.

Je vous remercie par avance pour toute réponse qui pourrait me sortir de ma galère actuelle

[Sylvain Bergeron]

Je n'ai pas le temps de garder à fond, mais je vois rapidement certaines erreurs :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
      program upward

c Déclarations

      character name*70
      real rft(512,512), flo(512,512), fla(512,512)

    1 format (a)
    2 format (f11.6,x,f10.6,x,f10.2)

      pi = 3.1415926535
      tpi = 2. * pi
      rad = pi / 180.

      dimension grid(nx*ny),store(2*nx*ny),nn(2)      
      real kx, ky, k
      complex cgrid, cmplx
      data pi/3.14159265/
...

Les lignes "1 format" à "rad = " doivent être déplacées juste après la ligne "data" car ce sont des lignes exécutables et en Fortran, les lignes exécutables doivent suivre les lignes de déclaration. De plus, pi est assigné 2 fois (pi = et data pi...).

Ensuite, tu déclares des tableaux (grid et store) dont tu ne connais pas la taille. Tu dois les délacer comme suit :

Code :

real, allocatable :: grid(:),store(:)

ou (au choix) :

Code :

real, allocatable, dimension(:) :: grid,store

Quand nx et ny seront connus (donc après les "read"), tu devras allouer les tableaux :

Code :

allocate (grid(nx*ny),store(2*nx*ny))

Finalement, comme ton programme est petit et que tu utilises 2 routines, je te conseille d'ajouter une section "contains" à la fin de ton programme, juste avant le "end" et y insérer les 2 routines (attention: les end à la fin des subroutines doivent être remplacés par des end subroutine) . Le compilateurs fera alors un meilleur diagnostic des erreurs :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
      close (70)

      stop

      contains

      subroutine fourn(data,nn,ndim,isign)
c
c  Replaces DATA by its NDIM-dimensional discrete Fourier transform, 
...
   18 continue
      return
      end subroutine

      subroutine kvalue(i,j,nx,ny,dkx,dky,kx,ky)
c  Subroutine KVALUE finds the wavenumber coordinates of one 
...
            end if
      return
      end subroutine

      end

[Ehouarn]

Bonjour,

Envoyé par Flo Flo

Bien entendu, comme tout programme qui se respecte, il bugge

Je précise que je suis complètement nul en programmation. J'ai essayé de faire ce que je pouvais mais il continue de m'afficher toute une gamme de messages d'erreurs quand je le compile et je ne comprends pas forcément ce que je dois faire pour que ça fonctionne.

Je ne peux que te recommander chaudement d'utiliser les options de compilation du compilateur qui permettent justement d'identifier les problèmes et erreurs "de base" (variable utilisée alors que non initialisée, tableau non alloué, dépassement d'indice dans un tableau, etc.).
Par exemple avec le compilateur gfortran: "-g -O2 -Wall -fbounds-check -ffpe-trap=invalid,zero,overflow".
Consulter leur doc pour les autres compilateurs.

Bonne continuation.