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 : Sélectionner tout - Visualiser dans une fenêtre à part
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 : Sélectionner tout - Visualiser dans une fenêtre à part
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 : Sélectionner tout - Visualiser dans une fenêtre à part
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