MODULE para
IMPLICIT NONE

INTEGER         :: N=32
DOUBLE COMPLEX   ::imi=(0.0d0,1.0d0)
DOUBLE PRECISION :: pi=DACOS(-1.0d0)
DOUBLE PRECISION :: length=1.0
END MODULE

!_______________________________________________________________________________

PROGRAM Zweidim
!USE Dislin
USE para
Implicit none
!_______________________________________________________________________________
#include'fftw3.f'
!Befehlzum Kompilieren: ifort Burgers.f90 -fpp -lfftw3 -WARN -ldislin  
!_______________________________________________________________________________
INTEGER                              k,i
!Parameter N=32
INTEGER*8                            plan_hin,plan_her!
DOUBLE PRECISION,DIMENSION(:,:), ALLOCATABLE       :: in,in2,in3
DOUBLE PRECISION,DIMENSION(:,:,:), ALLOCATABLE  :: kvec
DOUBLE COMPLEX,Dimension(:,:) , ALLOCATABLE     :: us,out,ableitung,delta
!_______________________________________________________________________________
!Programm zur Berechnung der Ableitungen eines zweidimensionalen Feldes
!_______________________________________________________________________________
length=2.0d0*pi


!Initialisierung des Feldes
ALLOCATE(in(N,N),in2(N,N),in3(N,N),us(N/2+1,N),out(N/2+1,N),kvec(N/2+1,N,3),ableitung(N/2+1,N),delta(N/2+1,N))

open(10,File='Sinusfunktion.dat')
DO k=1,N
  DO i=1,N

    in(i,k)=DSIN(2.0*pi*(i-1)/DBLE(N))+DSIN(2.0*pi*(k-1)/DBLE(N))


  write(10,*) in(i,k)
  END DO

  write(10,*) ""

END DO
close(10)
!_______________________________________________________________________________
! Wellenzahlen initialisieren
CALL kvec_init(kvec)
!_______________________________________________________________________________
! Plaene initialisieren
CALL dfftw_plan_dft_r2c_2d(plan_hin,N,N,in,out,FFTW_ESTIMATE)
CALL dfftw_plan_dft_c2r_2d(plan_her,N,N,out,in,FFTW_ESTIMATE)
!_______________________________________________________________________________
! in den Fourierraum
CALL dfftw_execute_dft_r2c(plan_hin,in,us)
us=us/DBLE(N)

CALL ableit(kvec,us,ableitung)

CALL laplace(kvec,us,delta)

!_______________________________________________________________________________
! in den Ortsraum
CALL dfftw_execute_dft_c2r(plan_her,us,in)
in=in/DBLE(N)

open(14,File='Ergebnis.dat')

DO k=1,N
   DO i=1,N
	  write(14,*) in(i,k)
   END DO
write(14,*) ""
END DO
close(14)


CALL dfftw_execute_dft_c2r(plan_her,ableitung,in2)
in2=in2/DBLE(N)
open(15,File='Ableitung.dat')
DO k=1,N
   DO i=1,N

	  write(15,*) in2(i,k)
   END DO
write(15,*) ""
END DO
close(15)


CALL dfftw_execute_dft_c2r(plan_her,delta,in3)
in3=in3/DBLE(N)
open(16,File='Laplace.dat')
DO k=1,N
   DO i=1,N

	  write(16,*) in3(i,k)
   END DO
write(16,*) ""
END DO
close(16)




!_______________________________________________________________________________
!Zerstörung der Plaene
CALL dfftw_destroy_plan(plan_hin)
CALL dfftw_destroy_plan(plan_her)

DEALLOCATE(in,in2,in3,us,out,kvec,ableitung,delta)

END PROGRAM Zweidim

!_______________________________________________________________________________

SUBROUTINE kvec_init(kvec)
USE para
IMPLICIT NONE
INTEGER                               :: i,k
DOUBLE PRECISION, DIMENSION(N/2+1,N,3):: kvec


DO k=1,N
  DO i=1,N/2+1

    IF(k.le.N/2+1) kvec(i,k,1:2)=(/ 2.0*pi/length*(i-1),2.0*pi/length*(k-1) /)
    IF(k.gt.N/2+1) kvec(i,k,1:2)=(/ 2.0*pi/length*(i-1),2.0*pi/length*(-N+k-1) /)

  END DO 
END DO 

DO k=1,N
   DO i=1,N/2+1
    kvec(i,k,3)= kvec(i,k,1)**2+kvec(i,k,2)**2 
  END DO
END DO

PRINT*,""
PRINT*,"  kvec initialisiert"
PRINT*,""

END SUBROUTINE
!_______________________________________________________________________________

SUBROUTINE Ableit(kvec,us,ableitung)
USE para
IMPLICIT NONE
INTEGER                               :: i,k
DOUBLE PRECISION, DIMENSION(N/2+1,N,3):: kvec
DOUBLE COMPLEX  , DIMENSION(N/2+1,N)  :: us,ableitung



DO i=1,N/2+1
   DO k=1,N
    ableitung(i,k)=imi*kvec(i,k,1)*us(i,k)
  End DO 
End DO 

PRINT*,""
PRINT*,"  Ableitung ausgeführt"
PRINT*,""
END SUBROUTINE
!_______________________________________________________________________________
SUBROUTINE Laplace(kvec,us,delta)
USE para
IMPLICIT NONE
INTEGER :: i,k
DOUBLE PRECISION, DIMENSION(N/2+1,N,3):: kvec
DOUBLE COMPLEX  , DIMENSION(N/2+1,N)  :: us,delta

DO i=1,N/2+1
  DO k=1,N
    delta(i,k)=-kvec(i,k,3)*us(i,k)
  End DO 
End DO 

PRINT*,""
PRINT*,"  Laplace angewendet"
PRINT*,""
END SUBROUTINE

!_______________________________________________________________________________