!***********************************************************************************
MODULE para
IMPLICIT NONE
!_______________________________________________________________________________
INTEGER          :: N=128!Anzahl
DOUBLE PRECISION :: length=50.0d0!Länge des Systems
DOUBLE COMPLEX   :: imi=(0.0d0,1.0d0) !Imaginäre Einheit
DOUBLE PRECISION :: pi=DACOS(-1.0d0)
DOUBLE PRECISION :: eps=1.0d0!Vorfaktor für linearen psi-Term
DOUBLE PRECISION :: deltat=0.1!deltat=0.05d0!Zeitschritt
INTEGER          :: screen_every=2000!für dislin-Ausgabe
INTEGER          :: n_iter=1000000!Anzahl der Zeitschritte
DOUBLE PRECISION :: alpha=1.0d0
DOUBLE PRECISION :: beta=2.0d0
DOUBLE PRECISION :: amp=0.01
Logical          :: screen=.true.!.FALSE.
Integer          :: ps=2500
END MODULE
!***********************************************************************************
!
!***********************************************************************************
PROGRAM GinzburgLandau
USE para
IMPLICIT NONE
!_______________________________________________________________________________
#include'fftw3.f'
!_______________________________________________________________________________
INTEGER                                         i
INTEGER*8                                       plan_hin,plan_her,plan_dummy
DOUBLE PRECISION,DIMENSION(:,:)  ,ALLOCATABLE:: dummy2_feld,dummy_betrag
DOUBLE PRECISION,DIMENSION(:,:,:),ALLOCATABLE:: kvec
DOUBLE COMPLEX  ,Dimension(:,:)  ,ALLOCATABLE:: in,in2,psi,out,ableitung,psiin,k1,k2,k3,k4
DOUBLE COMPLEX  ,Dimension(:,:)  ,ALLOCATABLE:: dummy_c_feld,dummy_feld
!_______________________________________________________________________________
ALLOCATE(in(N,N),in2(N,N),psi(N,N),psiin(N,N),out(N,N),ableitung(N,N))
ALLOCATE(dummy_c_feld(N,N),dummy_feld(N,N),k1(N,N),k2(N,N),k3(N,N),k4(N,N))
ALLOCATE(kvec(N,N,3),dummy2_feld(N,N),dummy_betrag(N,N))
!_______________________________________________________________________________
!Programm zur Berechnung der Swift-Hohenberg-Gleichung
!_______________________________________________________________________________
!Initialisierung des Feldes
in=0
Call initial(in)
!_______________________________________________________________________________
! Wellenzahlen initialisieren
CALL kvec_init(kvec)
!_______________________________________________________________________________
! Dislin initialisieren
IF(screen) THEN
! CALL METAFL ('XWIN')    ! Dislin Initialisierung
! CALL PAGE(5000,2800)
! CALL SCRMOD('REVERS')
! CALL SCLMOD('FULL')
! CALL DISINI()

CALL METAFL ('XWIN')    ! Dislin Initialisierung
CALL PAGE(2*ps,ps)
CALL SCRMOD('REVERS')
CALL SCLMOD('FULL')
CALL WINSIZ(ps/2,ps/4)
CALL DISINI()
CALL AX3LEN((2*ps)/3,(2*ps)/3,(2*ps)/3)

END IF
!_______________________________________________________________________________
! Plaene initialisieren
CALL dfftw_plan_dft_2d(plan_hin,N,N,in,out,FFTW_FORWARD,FFTW_ESTIMATE)
CALL dfftw_plan_dft_2d(plan_her,N,N,out,in,FFTW_BACKWARD,FFTW_ESTIMATE)
CALL dfftw_plan_dft_2d(plan_dummy,N,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))
!_______________________________________________________________________________
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))

     dummy2_feld =REAL(dummy_feld)
     dummy_betrag=ABS(dummy_feld)

     CALL display(dummy2_feld,dummy_betrag)
   END IF

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

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,i
DOUBLE PRECISION, DIMENSION(N,N,3):: kvec
!_______________________________________________________________________________

DO k=1,N
   DO i=1,N
      if(k.le.N/2+1.and.i.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.and.i.le.N/2+1) kvec(i,k,1:2)=(/ 2.0*pi/length*(i-1)   ,2.0*pi/length*(-N+k-1) /)
      if(k.le.N/2+1.and.i.gt.N/2+1) kvec(i,k,1:2)=(/ 2.0*pi/length*(-N+i-1),2.0*pi/length*(k-1)    /)
      if(k.gt.N/2+1.and.i.gt.N/2+1) kvec(i,k,1:2)=(/ 2.0*pi/length*(-N+i-1),2.0*pi/length*(-N+k-1) /)
   END DO 
END DO 

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


!_______________________________________________________________________________
END SUBROUTINE
!***********************************************************************************
!
!
!
!***********************************************************************************
SUBROUTINE display(real_feld,abst_feld)
USE para
IMPLICIT NONE
!_______________________________________________________________________________
REAL :: mini, maxi
DOUBLE PRECISION, DIMENSION(N,N) :: real_feld,abst_feld
CHARACTER(150)  ::  filename
!_______________________________________________________________________________
! Fuer z-Skalierung Minimum und Maximum bestimmen

!_______________________________________________________________________________
! Graph formatieren, ausgeben und beenden
IF(.NOT.screen) THEN
   CALL METAFL ('PNG')    ! Dislin Initialisierung
   CALL PAGE(2*ps,ps)
   CALL SCRMOD('REVERS')
   CALL SCLMOD('FULL')
   CALL WINSIZ(ps/2,ps/4)
   filename=TRIM("B")//".png"
   CALL SETFIL(TRIM(filename))
   CALL DISINI()
   CALL AX3LEN((2*ps)/3,(2*ps)/3,(2*ps)/3)
END IF

maxi=MAXVAL(real_feld)
mini=MINVAL(real_feld)



! CALL ERASE()
! CALL SETVLT ('TEMP')
! CALL AXSPOS(ps+ps/10,(5*ps)/6)
! CALL AXSSCL('LIN','XY')
! CALL NAME('','X')
! CALL NAME('','Y')
! !CALL GRAF3(0.0, REAL(dmx/2), 0.0, REAL(dmx/4) , 0.0, REAL(dmy/2),0.0, REAL(dmy/4), mini, maxi, mini, (maxi-mini)/10)
! CALL GRAF3(0.0, REAL(dmx), 0.0, REAL(dmx/4) , 0.0, REAL(dmy) ,0.0,
! REAL(dmy/4), mini, maxi, mini, (maxi-mini)/10)
! CALL TITLE()
! CALL CRVMAT(real_feld(dmx/4+1:3*dmx/4,dmy/4+1:3*dmy/4),dmx/2,dmy/2,8,8)
! !CALL CRVMAT(real_feld,dmx,dmy,8,8)
! CALL COLOR('WHITE')
! CALL ENDGRF()
! 
! CALL SETVLT ('RAIN')
! CALL AXSPOS(ps/10,(5*ps)/6)
! CALL AXSSCL('LIN','XY')
! CALL NAME('','X')
! CALL NAME('','Y')
! CALL GRAF3(0.0, REAL(dmx), 0.0, REAL(dmx/4) , 0.0, REAL(dmy) ,0.0,REAL(dmy/4), mini, maxi, mini, (maxi-mini)/10)
! CALL TITLE()
! CALL CRVMAT(real_feld,dmx,dmy,8,8)
! CALL COLOR('WHITE')
! CALL ENDGRF()






CALL ERASE()
CALL AXSPOS(ps/10,(5*ps)/6)
CALL SETVLT ('RAIN')
CALL AXSSCL('LIN','XY')
CALL NAME('','X')
CALL NAME('','Y')
CALL GRAF3(0.0, REAL(N), 0.0, REAL(N/4), 0.0, REAL(N) ,0.0, REAL(N/4), mini, maxi, mini, (maxi-mini)/10)
CALL TITLE()
CALL CRVMAT(REAL(real_feld),N,N,8,8)
CALL COLOR('WHITE')
CALL ENDGRF()

maxi=MAXVAL(abst_feld)
mini=MINVAL(abst_feld)



CALL AXSPOS(ps+ps/10,(5*ps)/6)
CALL SETVLT ('RAIN')
CALL AXSSCL('LIN','XY')
CALL NAME('','X')
CALL NAME('','Y')
CALL GRAF3(0.0, REAL(N), 0.0, REAL(N/4), 0.0, REAL(N) ,0.0, REAL(N/4), mini, maxi, mini, (maxi-mini)/10)
CALL TITLE()
CALL CRVMAT(REAL(abst_feld),N,N,8,8)
CALL COLOR('WHITE')
CALL ENDGRF()

IF(.NOT.screen)  CALL DISFIN()

!_______________________________________________________________________________
END SUBROUTINE
!***********************************************************************************
!
!
!
!***********************************************************************************
SUBROUTINE EulerImplizit(psi,out,kvec,plan_hin,plan_her)
USE para
IMPLICIT NONE
!_______________________________________________________________________________
INTEGER                          :: i,k
INTEGER*8                           plan_hin,plan_her
DOUBLE PRECISION,DIMENSION(N,N,3):: kvec
DOUBLE COMPLEX  ,Dimension(N,N)  :: psi,out,psinichtlin
!_______________________________________________________________________________
CALL Nichtlinear(psi,psinichtlin,plan_hin,plan_her)

DO k=1,N
DO i=1,N
   out(i,k)=(psi(i,k)/deltat-psinichtlin(i,k))/(1.0d0/deltat-1.0d0+(1.0d0+imi*alpha)*(kvec(i,k,3)))
END DO
END DO		
!_______________________________________________________________________________
END SUBROUTINE
!***********************************************************************************
!
!
!
!***********************************************************************************
SUBROUTINE Nichtlinear(in,psiquadratpsi,plan_hin,plan_her)
USE para
IMPLICIT NONE
!_______________________________________________________________________________
#include'fftw3.f'
!_______________________________________________________________________________
INTEGER                       :: i,k
INTEGER*8                        plan_hin,plan_her
DOUBLE COMPLEX  ,Dimension(N,N) :: psiquadratpsi,in,psiquadratpsirealraum,reserve,psirealraum
!_______________________________________________________________________________
reserve=in
!___________________________________________________________
!Übergang in Realraum zur Berechnung der Nichtlinearität
CALL dfftw_execute_dft(plan_her,reserve,psirealraum)
psirealraum=psirealraum/(DBLE(N))
!_______________________________________________________________________________
DO k=1,N
DO i=1,N
     psiquadratpsirealraum(i,k)=(1.0d0+imi*beta)*psirealraum(i,k)*(DBLE(psirealraum(i,k))**2+DIMAG(psirealraum(i,k))**2)
END DO
END DO
!_______________________________________________________________________________
!Übergang zurück in den Fourier-Raum
CALL dfftw_execute_dft(plan_hin,psiquadratpsirealraum,psiquadratpsi)! 
psiquadratpsi=psiquadratpsi/(DBLE(N))

!_______________________________________________________________________________
END SUBROUTINE Nichtlinear
!***********************************************************************************
!
!
!
!***********************************************************************************
SUBROUTINE Initial(in)
USE para
IMPLICIT NONE
!_______________________________________________________________________________
INTEGER                    :: k,i
DOUBLE COMPLEX  ,DIMENSION(N,N):: in
DOUBLE PRECISION,DIMENSION(N,N):: realfeld,imagfeld
!_______________________________________________________________________________

DO k=1,N
DO i=1,N
      CALL RANDOM_NUMBER(realfeld(i,k))
      CALL RANDOM_NUMBER(imagfeld(i,k))
      in(i,k)=amp*((realfeld(i,k)-0.5)+imi*(imagfeld(i,k)-0.5))
END DO
END DO
!_______________________________________________________________________________
END SUBROUTINE Initial
!***********************************************************************************
!