!***********************************************************************************
MODULE para
IMPLICIT NONE
!_______________________________________________________________________________
INTEGER          :: N=512!Anzahl
DOUBLE PRECISION :: length=101.0d0!Länge des Systems
DOUBLE COMPLEX   :: imi=(0.0d0,1.0d0) !Imaginäre Einheit
DOUBLE PRECISION :: pi=DACOS(-1.0d0)
DOUBLE PRECISION :: deltat=0.001!0.0001!1e-4!deltat=0.05d0!Zeitschritt
INTEGER          :: screen_every=20!für dislin-Ausgabe
INTEGER          :: n_iter=1000!Anzahl der Zeitschritte
DOUBLE PRECISION :: alpha=1.0d0!2.0d0!0.5d0!0.0d0
DOUBLE PRECISION :: beta=2.0d0!-1.0d0!-2.0d0!-1.5d0!-4.0d0
DOUBLE PRECISION :: amp=0.01
END MODULE
!***********************************************************************************
!
!***********************************************************************************
PROGRAM GinzburgLandau
USE para
IMPLICIT NONE
!_______________________________________________________________________________
#include'fftw3.f'
!_______________________________________________________________________________
INTEGER                                       i,k
INTEGER*8                                     plan_hin,plan_her,plan_dummy
DOUBLE PRECISION,DIMENSION(:)  ,ALLOCATABLE:: dummy2_feld
DOUBLE PRECISION,DIMENSION(:)  ,ALLOCATABLE:: kvec
DOUBLE COMPLEX  ,Dimension(:)  ,ALLOCATABLE:: in,in2,psi,out,k1,k2,k3,k4
DOUBLE COMPLEX  ,Dimension(:)  ,ALLOCATABLE:: dummy_c_feld, dummy_feld
!_______________________________________________________________________________
ALLOCATE(in(N),in2(N),psi(N),out(N))
ALLOCATE(dummy_c_feld(N),dummy_feld(N),k1(N),k2(N),k3(N),k4(N))
ALLOCATE(kvec(N),dummy2_feld(N))
!_______________________________________________________________________________
!Programm zur Berechnung der Swift-Hohenberg-Gleichung
!_______________________________________________________________________________
!Initialisierung des Feldes
in=0

Call initial(in)

!_______________________________________________________________________________
! Wellenzahlen initialisieren
CALL kvec_init(kvec)
!_______________________________________________________________________________
! Dislin initialisieren
CALL METAFL ('XWIN')    ! Dislin Initialisierung
CALL PAGE(2500,2500)
CALL SCRMOD('REVERS')
CALL SCLMOD('FULL')
CALL DISINI()

!_______________________________________________________________________________
! Plaene initialisieren
CALL dfftw_plan_dft_1d(plan_hin,N,in,out,FFTW_FORWARD,FFTW_ESTIMATE)
CALL dfftw_plan_dft_1d(plan_her,N,out,in,FFTW_BACKWARD,FFTW_ESTIMATE)
CALL dfftw_plan_dft_1d(plan_dummy,N,dummy_c_feld,dummy_feld,FFTW_BACKWARD,FFTW_ESTIMATE)
!_______________________________________________________________________________
! in den Fourierraum
CALL dfftw_execute_dft(plan_hin,in,psi)
psi=psi/(DBLE(N))**0.5d0
!_______________________________________________________________________________
DO i=0,n_iter

   IF(MOD(i,screen_every).EQ.0) THEN
      write(*,*) i*deltat
     dummy_c_feld=psi
     CALL dfftw_execute_dft(plan_dummy,dummy_c_feld,dummy_feld)
     dummy_feld=(dummy_feld)/(DBLE(N))**0.5d0
     DO k=1,N
     dummy2_feld(k)=REAL(dummy_feld(k))
     END DO

     CALL display(dummy2_feld,DIMAG(dummy_feld))

   END IF

   
#ifdef euler

  CALL EulerImplizit(psi,psi,kvec,plan_hin,plan_her)

#else

 CALL ETD2(psi,psi,kvec,plan_hin,plan_her)

#endif

!kompileren mit: ifort ... -fpp -D euler

END DO

!_______________________________________________________________________________
!Zerstörung der Plaene
CALL dfftw_destroy_plan(plan_hin)
CALL dfftw_destroy_plan(plan_her)
CALL dfftw_destroy_plan(plan_dummy)
!_______________________________________________________________________________
CALL DISFIN()
!_______________________________________________________________________________
END PROGRAM
!*******************************************************************************
!
!
!
!*******************************************************************************
SUBROUTINE kvec_init(kvec)
USE para
IMPLICIT NONE
!_______________________________________________________________________________
INTEGER                       :: k
DOUBLE PRECISION, DIMENSION(N):: kvec
!_______________________________________________________________________________
DO k=1,N
       if(k.le.N/2+1) kvec(k)=2.0*pi/length*(k-1) 
       if(k.gt.N/2+1) kvec(k)=2.0*pi/length*(-N+k-1)
END DO 

!_______________________________________________________________________________
PRINT*,""
PRINT*,"  kvec initialisiert"
PRINT*,""
END SUBROUTINE
!*******************************************************************************
!
!
!
!*******************************************************************************
SUBROUTINE display(real_feld,imagfeld)
USE para
IMPLICIT NONE
!_______________________________________________________________________________
INTEGER :: i
REAL :: mini, maxi
REAL, DIMENSION(N) :: x_achse, y_achse,y2_achse,y3_achse
DOUBLE PRECISION, DIMENSION(N) :: real_feld,imagfeld
!_______________________________________________________________________________
! Fuer z-Skalierung Minimum und Maximum bestimmen
!maxi=MAXVAL(real_feld)
!mini=MINVAL(real_feld)
!_______________________________________________________________________________
! Graph formatieren, ausgeben und beenden

DO i=1,N
   x_achse(i)=REAL(i-1)
   y_achse(i)=REAL(real_feld(i))
   y2_achse(i)=REAL(imagfeld(i))
   y3_achse(i)=(REAL(real_feld(i)))**2+(REAL(imagfeld(i)))**2
END DO

maxi=1.2
mini=-1.2

CALL ERASE()
CALL AXSSCL('LIN','XY')
CALL NAME('x-axis','X')
CALL NAME('y-axis','Y')
CALL GRAF(0.0,REAL(N),0.0,100.0,mini,maxi ,mini,(maxi-mini)/10.0)
CALL CURVE(x_achse,y_achse,n)
CALL COLOR('Green')
CALL CURVE(x_achse,y2_achse,n)
CALL COLOR('Red')
CALL CURVE(x_achse,y3_achse,n)
CALL COLOR('WHITE')
CALL ENDGRF()
!_______________________________________________________________________________
END SUBROUTINE
!***********************************************************************************
!
!
!
!***********************************************************************************
SUBROUTINE EulerImplizit(psi,out,kvec,plan_hin,plan_her)
USE para
IMPLICIT NONE
!_______________________________________________________________________________
INTEGER                        :: k
INTEGER*8                         plan_hin,plan_her
DOUBLE PRECISION,DIMENSION(N)  :: kvec
DOUBLE COMPLEX  ,Dimension(N)  :: psi,out,psinichtlin
!_______________________________________________________________________________
CALL Nichtlinear(psi,psinichtlin,plan_hin,plan_her)

DO k=1,N

   out(k)=(psi(k)/deltat-psinichtlin(k))/(1.0d0/deltat-1.0d0+(1.0d0+imi*alpha)*(kvec(k))**2)

END DO	
!_______________________________________________________________________________
END SUBROUTINE
!***********************************************************************************
!
!
!
!***********************************************************************************
SUBROUTINE Nichtlinear(in,out,plan_hin,plan_her)
USE para
IMPLICIT NONE
!_______________________________________________________________________________
#include'fftw3.f'
!_______________________________________________________________________________
INTEGER                       :: k
INTEGER*8                        plan_hin,plan_her
DOUBLE COMPLEX  ,Dimension(N) :: psiquadratpsi,out,in,psiquadratpsirealraum,reserve,psirealraum
!_______________________________________________________________________________
DO k=1,N
      reserve(k)=in(k)
END DO
!___________________________________________________________
!Übergang in Realraum zur Berechnung der Nichtlinearität
CALL dfftw_execute_dft(plan_her,reserve,psirealraum)
psirealraum=psirealraum/(DBLE(N))**0.5d0
!_______________________________________________________________________________
DO k=1,N
      psiquadratpsirealraum(k)=(1.0d0+imi*beta)*psirealraum(k)*(DBLE(psirealraum(k))**2+DIMAG(psirealraum(k))**2)
END DO
!_______________________________________________________________________________
!Übergang zurück in den Fourier-Raum
CALL dfftw_execute_dft(plan_hin,psiquadratpsirealraum,psiquadratpsi)! 
psiquadratpsi=psiquadratpsi/(DBLE(N))**0.5d0
!_______________________________________________________________________________
DO k=1,N
      out(k)=psiquadratpsi(k)
END DO	
!_______________________________________________________________________________
END SUBROUTINE Nichtlinear
!***********************************************************************************
!
!
!
!***********************************************************************************
SUBROUTINE Initial(in)
USE para
IMPLICIT NONE
!_______________________________________________________________________________
INTEGER                    :: k
DOUBLE COMPLEX  ,DIMENSION(N):: in
DOUBLE PRECISION,DIMENSION(N):: realfeld,imagfeld
!_______________________________________________________________________________
DO k=1,N
      CALL RANDOM_NUMBER(realfeld(k))
      CALL RANDOM_NUMBER(imagfeld(k))
      in(k)=amp*((realfeld(k)-0.5)+imi*(imagfeld(k)-0.5))
    !in(k)=sqrt(1.0d0-(20.0d0*pi/length)**2)*exp(imi*20.0d0*pi/length*(REAL(k-1)/REAL(N)*(length)-50.0d0))+in(k)
    in(k)=1.0d0+in(k)
    !in(k)=        2.0d0/(exp((k*length/real(N)+length/4.0d0)**2)+exp(-(k*length/real(N)+length/4.0d0)**2))&
     !      +0.8d0*2.0d0/(exp((k*length/real(N)-length/4.0d0)**2)+exp(-(k*length/real(N)-length/4.0d0)**2))&
     !      +in(k)

END DO
!_______________________________________________________________________________
END SUBROUTINE Initial
!***********************************************************************************
!
!
!
!***********************************************************************************
SUBROUTINE ETD2(psi,out,kvec,plan_hin,plan_her)
USE para
IMPLICIT NONE
!_______________________________________________________________________________
INTEGER                        :: k
INTEGER*8                         plan_hin,plan_her
DOUBLE PRECISION,DIMENSION(N)  :: kvec
DOUBLE COMPLEX  ,Dimension(N)  :: psi,out,psinichtlin,psialt
DOUBLE PRECISION               :: q
!_______________________________________________________________________________
CALL Nichtlinear(psi,psinichtlin,plan_hin,plan_her)

DO k=1,N
   q=(1.0d0-(1.0d0+imi*alpha)*(kvec(k))**2)
   out(k)=psi(k)*exp(q*deltat)&
         +(-psinichtlin(k))*((1.0d0+deltat*q)*exp(q*deltat)-1.0d0-2.0d0*q*deltat)/(deltat*q**2)&
         +(-psialt(k))*(-exp(q*deltat)+1.0d0+deltat*q)/(deltat*q**2)

   psialt(k)=psinichtlin(k)

END DO	
!_______________________________________________________________________________
END SUBROUTINE
!***********************************************************************************
!