PROGRAM ADVEKTION

USE Dislin

Implicit none
!_______________________________________________________________________________
#include'fftw3.f'

!Befehlzum Kompilieren: ifort Advektion.f90 -fpp -lfftw3 -WARN -ldislin  
!_______________________________________________________________________________
INTEGER k,i
PARAMETER N=2048, n_iter=80000, screen_every=100
DOUBLE PRECISION:: pi,deltat,c
DOUBLE COMPLEX:: imi
INTEGER*8 plan_hin,plan_her, plan_dummy
DOUBLE PRECISION, DIMENSION(N):: in, dummy_feld
DOUBLE PRECISION, DIMENSION(N/2+1):: kvec
DOUBLE COMPLEX, Dimension(N/2+1):: out,K1,K2,K3,K4,us, dummy_c_feld

!_______________________________________________________________________________
!Programm zur Berechnung der Advektionsgleichung
 c=0.01
imi=(0.0d0,1.0d0)
pi=DACOS(-1.0d0)
Deltat=0.01
!_______________________________________________________________________________

!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)
	  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, out)
out=out/DBLE(N)**0.5d0

!_______________________________________________________________________________
! Runge-Kutta-Verfahren für Zeitableitung

do k=1,N/2+1
us(k)=out(k)
end do


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,i,c,deltat)

!Pause
  END IF

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

  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 ADVEKTION

!*******************************************************************************
!*******************************************************************************


SUBROUTINE RechteHand(in,out,n,c,kvec)
!_______________________________________________________________________________
!Rechte-Hand-Seite

INTEGER :: n,k
DOUBLE COMPLEX,DIMENSION(N/2+1):: in,out
DOUBLE PRECISION,DIMENSION(N/2+1):: kvec
DOUBLE PRECISION:: c
DOUBLE COMPLEX:: imi

 
imi=DCMPLX(0.0d0,1.0d0)


! Ableitung berechnen

	DO k=1,N/2+1
	  out(k)=-c*imi*kvec(k)*in(k)
	END DO
!_______________________________________________________________________________

END SUBROUTINE RechteHand


!***********************************************************************************
!***********************************************************************************

SUBROUTINE display(feld_in,n,zeit,c,deltat)

USE dislin

IMPLICIT NONE

INTEGER :: i,n,zeit,IC
DOUBLE PRECISION, DIMENSION(n)  :: feld_in,xwert
REAL, DIMENSION(n) :: x_achse, y_achse,y2_achse,y3_achse
REAL :: pi, mini, maxi
DOUBLE PRECISION:: c,deltat


pi=ACOS(-1.0)



DO i=1,n
  x_achse(i)=2.0*pi/REAL(n)*(i-1)
  xwert(i)=2.0*pi*(i-1)/REAL(N)-(c*zeit*deltat)
	DO WHILE (xwert(i)<0)
	xwert(i)=xwert(i)+2*pi
	END DO

  y_achse(i)=REAL(feld_in(i))
  y2_achse(i)=REAL(EXP(-2.0*pi*(xwert(i)-pi/2.0)**2))
  y3_achse(i)=REAL(EXP(-2.0*pi*(2.0*pi*(i-1)/REAL(N)-pi/2.0)**2))
END DO


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('Advektionsgleichung',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()

IC=INTRGB(0.95,0.95,0.95)
CALL AXSBGD(IC)

!CALL SETRGB(0.7,0.7,0.7)
CALL GRID(1,1)


CALL COLOR('Green')
CALL CURVE(x_achse,y_achse,n)

CALL COLOR('Blue')
CALL CURVE(x_achse,y2_achse,n)

CALL COLOR('Red')
CALL CURVE(x_achse,y3_achse,n)

CALL ENDGRF()


END SUBROUTINE


!***********************************************************************************
!***********************************************************************************