PROGRAM Burgers

USE Dislin

Implicit none
!_______________________________________________________________________________
#include'fftw3.f'
!Befehlzum Kompilieren: ifort Burgers.f90 -fpp -lfftw3 -WARN -ldislin  
!_______________________________________________________________________________
INTEGER                              k,i,g
PARAMETER                            N=2048, n_iter=80000, screen_every=200
INTEGER*8                            plan_hin,plan_her, plan_dummy
DOUBLE PRECISION                  :: pi,deltat,c,amplitude
DOUBLE PRECISION, DIMENSION(N)    :: in, dummy_feld,sinusse
DOUBLE PRECISION, DIMENSION(N/2+1):: kvec
DOUBLE COMPLEX                    :: imi
DOUBLE COMPLEX, Dimension(N/2+1)  :: out,K1,K2,K3,K4,us, dummy_c_feld
!_______________________________________________________________________________
!Programm zur Berechnung der Advektionsgleichung
 c=0.01 !Viskosität
imi=(0.0d0,1.0d0)
pi=DACOS(-1.0d0)
Deltat=0.0001
!_______________________________________________________________________________
!Initialisierung des Feldes

open(10,File='Exponentialfunktion.dat')
	DO k=1,N
	!  in(k)=exp(-2.0*pi*(2.0*pi*(k-1)/real(N)-pi/2.0)**2)
         DO g=1,12
           amplitude=(-1.0)**g!*(real(g))**(0.1)
           sinusse(k)=amplitude*sin(g*2.0*pi*(k-1)/real(N)+3.0*real(g))
           in(k)=in(k)+sinusse(k)
         END DO
       in(k)=in(k)+0.2
	  write(10,*) 2.0*pi/REAL(n)*(k-1),in(k)
	END DO
close(10)
!_______________________________________________________________________________
! Wellenzahlen initialisieren
DO k=1,N/2+1
   kvec(k)=k-1
END DO
!_______________________________________________________________________________
! 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_r2c_1d(plan_hin,N,in,out,FFTW_ESTIMATE)
CALL dfftw_plan_dft_c2r_1d(plan_her,N,out,in,FFTW_ESTIMATE)
CALL dfftw_plan_dft_c2r_1d(plan_dummy,N,dummy_c_feld,dummy_feld,FFTW_ESTIMATE)
!_______________________________________________________________________________
! in den Fourierraum
CALL dfftw_execute_dft_r2c(plan_hin,in,us)

us=us/DBLE(N)**0.5d0



!_______________________________________________________________________________
! Runge-Kutta-Verfahren für Zeitableitung
DO i=1,n_iter
  ! dislin, eine dummy Kopie vom Feld
  IF(MOD(i,screen_every).EQ.0) THEN
    dummy_c_feld=us
    CALL dfftw_execute_dft_c2r(plan_dummy,dummy_c_feld,dummy_feld)
    dummy_feld=dummy_feld/DBLE(N)**0.5d0
    CALL display(dummy_feld,n)
  END IF


  CALL RechteHand(us,k1,N,c,kvec,plan_hin,plan_her)
  CALL RechteHand(us+deltat/2.0*k1,k2,N,c,kvec,plan_hin,plan_her)
  CALL RechteHand(us+deltat/2.0*k2,k3,N,c,kvec,plan_hin,plan_her)
  CALL RechteHand(us+deltat*k3,k4,N,c,kvec,plan_hin,plan_her)

  us=us+(k1/6.0+k2/3.0+k3/3.0+k4/6.0)*Deltat
END DO
!_______________________________________________________________________________
! in den Ortsraum
CALL dfftw_execute_dft_c2r(plan_her,us,in)
in=in/DBLE(N)**0.5d0

open(14,File='Ergebnis.dat')
	DO k=1,N
	  write(14,*) 2.0*pi/REAL(n)*(k-1),in(k)
	END DO
close(14)
!_______________________________________________________________________________
!Zerstörung der Plaene
CALL dfftw_destroy_plan(plan_hin)
CALL dfftw_destroy_plan(plan_her)
!_______________________________________________________________________________
CALL DISFIN()

END PROGRAM Burgers
!*******************************************************************************
!*******************************************************************************
SUBROUTINE RechteHand(in,out,n,c,kvec,plan_hin,plan_her)

#include'fftw3.f'
!_______________________________________________________________________________
!Rechte-Hand-Seite

INTEGER :: N,k,i
INTEGER*8 plan_hin,plan_her
DOUBLE COMPLEX,DIMENSION(N/2+1):: in,out,dusus,duscomp,reserve
DOUBLE PRECISION,DIMENSION(N/2+1):: kvec
DOUBLE PRECISION,DIMENSION(N):: dususreal,usreal,dusreal
DOUBLE PRECISION:: c
DOUBLE COMPLEX:: imi

imi=DCMPLX(0.0d0,1.0d0)

!___________________________________________________________
! 1. Ableitung berechnen
 
  	DO k=1,N/2+1
            duscomp(k)=imi*kvec(k)*in(k)
            reserve(k)=in(k)
          END DO
!___________________________________________________________
!Übergang in Realraum zur Berechnung der Nichtlinearität
  CALL dfftw_execute_dft_c2r(plan_her,reserve,usreal)
  usreal=usreal/DBLE(N)**0.5d0
  CALL dfftw_execute_dft_c2r(plan_her,duscomp,dusreal)
  dusreal=dusreal/DBLE(N)**0.5d0
! 
! 
  DO i=1,N
    dususreal(i)=usreal(i)*dusreal(i)
  END DO
!  
!  !____________________________________________________________
!  !Übergang zurück in den Fourier-Raum
  CALL dfftw_execute_dft_r2c(plan_hin,dususreal,dusus)
  dusus=dusus/DBLE(N)**0.5d0

  DO k=1,N/2+1
    out(k)=-c*(kvec(k))**2*in(k)-dusus(k)
  End DO	

END SUBROUTINE RechteHand
!***********************************************************************************
!***********************************************************************************
SUBROUTINE display(feld_in,N)
USE dislin
IMPLICIT NONE

INTEGER :: i,N
DOUBLE PRECISION, DIMENSION(N)  :: feld_in
REAL, DIMENSION(n) :: x_achse, y_achse
REAL :: pi, mini, maxi
pi=ACOS(-1.0)

DO i=1,N
  x_achse(i)=2.0*pi/REAL(n)*(i-1)
  y_achse(i)=REAL(feld_in(i))
END DO

CALL ERASE()
CALL AXSPOS(450,1800)
CALL AXSLEN(1200,1200)
CALL AXSSCL('LIN','XY')
CALL NAME('x-axis','X')
CALL NAME('y-axis','Y')
CALL LABDIG(-1,'X')
CALL TICKS(10,'XY')
CALL TITLIN('Demonstration',1)
CALL TITLIN('Waermeleitungsgleichung',2)


mini=MINVAL(y_achse)
maxi=MAXVAL(y_achse)+0.1
CALL GRAF(-0.1, 2.0*pi, 0.0, 1.0 , mini, maxi ,mini, (maxi-mini)/10.0)
CALL SETRGB(0.7,0.7,0.7)
CALL GRID(1,1)
CALL COLOR('FORE')
CALL TITLE()
CALL GRID(1,1)

CALL COLOR('Red')
CALL CURVE(x_achse,y_achse,n)
CALL ENDGRF()

END SUBROUTINE
!***********************************************************************************
!***********************************************************************************