!***********************************************************************************
MODULE para
IMPLICIT NONE
!_______________________________________________________________________________
INTEGER          :: N=256!Anzahl
DOUBLE COMPLEX   :: imi=(0.0d0,1.0d0) !Imaginäre Einheit
DOUBLE PRECISION :: pi=DACOS(-1.0d0)
DOUBLE PRECISION :: deltat=0.001d0!Zeitschritt
INTEGER          :: screen_every=20!für dislin-Ausgabe
INTEGER          :: n_iter=1000!Anzahl der Zeitschritte
DOUBLE PRECISION :: nu=0.01d0
DOUBLE PRECISION :: amp=100.0d0
Logical          :: screen=.FALSE.!.TRUE.!.FALSE.
Integer          :: ps=2500
INTEGER          :: Startwert1x=128
INTEGER          :: Startwert1y=128-32
INTEGER          :: Startwert2x=128
INTEGER          :: Startwert2y=128+32
DOUBLE PRECISION :: mu=0.07d0
END MODULE
!***********************************************************************************
!
!***********************************************************************************
PROGRAM wirbel
USE para
Implicit none
!_______________________________________________________________________________
#include'fftw3.f'
!_______________________________________________________________________________
INTEGER                                         i
INTEGER*8                                       plan_hin,plan_her, plan_dummy
DOUBLE PRECISION,DIMENSION(:,:)  ,ALLOCATABLE:: omega_initial,omega_feld
DOUBLE PRECISION,DIMENSION(:,:,:),ALLOCATABLE:: kvec
DOUBLE COMPLEX  ,Dimension(:,:)  ,ALLOCATABLE:: omega,out
DOUBLE COMPLEX  ,Dimension(:,:)  ,ALLOCATABLE:: omega_c_feld,k1,k2,k3,k4,omega_dummy_feld
!_______________________________________________________________________________
ALLOCATE(omega_initial(N,N),omega(N/2+1,N),out(N/2+1,N),kvec(N/2+1,N,3),omega_dummy_feld(N/2+1,N))
ALLOCATE(omega_c_feld(N/2+1,N),omega_feld(N,N),k1(N/2+1,N),k2(N/2+1,N),k3(N/2+1,N),k4(N/2+1,N))
!_______________________________________________________________________________
! Wellenzahlen initialisieren
CALL kvec_init(kvec)
!_______________________________________________________________________________
!Initialisierung des Feldes der Wirbelstärke mit Zufallszahlen
Call initial(omega_initial)
!_______________________________________________________________________________
! Dislin initialisieren
IF(screen) THEN

CALL METAFL ('XWIN')
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_r2c_2d(plan_hin,N,N,omega_initial,out,FFTW_ESTIMATE)
CALL dfftw_plan_dft_c2r_2d(plan_her,N,N,out,omega_initial,FFTW_ESTIMATE)
CALL dfftw_plan_dft_c2r_2d(plan_dummy,N,N,omega_c_feld,omega_feld,FFTW_ESTIMATE)
!_______________________________________________________________________________
! omega_initial den Fourierraum
CALL dfftw_execute_dft_r2c(plan_hin,omega_initial,omega)
omega=omega/DBLE(N)
!omega ist Wirbelstärke im Fourier-Raum

CALL Dealising(omega,kvec)

!_______________________________________________________________________________
DO i=0,n_iter

   IF(MOD(i,screen_every).EQ.0) THEN

     write(*,*) i*deltat

     omega_c_feld=omega
     CALL dfftw_execute_dft_c2r(plan_dummy,omega_c_feld,omega_feld)
     omega_feld=omega_feld/DBLE(N)

     omega_dummy_feld=omega

     CALL display(omega_feld,omega_dummy_feld)

   END IF

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

   omega=omega+(k1/6.0+k2/3.0+k3/3.0+k4/6.0)*Deltat

!   CALL Dealising(omega,kvec)

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                               :: 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)=(/ (i-1),(k-1)    /)
      if(k.gt.N/2+1) kvec(i,k,1:2)=(/ (i-1),(-N+k-1) /)
   END DO 
END DO 
!_______________________________________________________________________________
!Betragsquadrat
DO i=1,N/2+1
   DO k=1,N
      kvec(i,k,3)= kvec(i,k,1)**2+kvec(i,k,2)**2 
   END DO
END DO
!_______________________________________________________________________________
END SUBROUTINE
!*******************************************************************************
!
!
!
!***********************************************************************************
SUBROUTINE display(omega_feld,feld_in)
USE para
IMPLICIT NONE
!_______________________________________________________________________________
INTEGER                               :: i,k
REAL                                  :: mini,maxi
DOUBLE PRECISION, DIMENSION(N,N)      :: omega_feld
DOUBLE COMPLEX  , DIMENSION (N/2+1,N) :: feld_in
REAL            , DIMENSION(N,N)      :: fourierfeld
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(omega_feld)
mini=MINVAL(omega_feld)

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(omega_feld),N,N,8,8)
CALL COLOR('WHITE')
CALL ENDGRF()


! Fouriertransformierte plotten >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
DO k=1,N
  DO i=1,N/2+1
        IF(k.LE.N/2+1)  THEN
          fourierfeld(i+N/2-1,k+N/2-1)=REAL(feld_in(i,k))**2+AIMAG(feld_in(i,k))**2
        ELSE
          fourierfeld(i+N/2-1,k-N/2-1)=REAL(feld_in(i,k))**2+AIMAG(feld_in(i,k))**2
        END IF
  END DO
END DO

DO k=1,N
  DO i=1,N/2
     fourierfeld(i,k)=fourierfeld(N-i,k)
  END DO
END DO

DO k=1,N
  DO i=1,N
    IF(fourierfeld(i,k).NE.0.0) THEN
      fourierfeld(i,k)=LOG(fourierfeld(i,k))
    ELSE
      fourierfeld(i,k)=0
    END IF
  END DO
END DO

!fourierfeld=(fourierfeld)**0.5

! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

! Fuer z-Skalierung Minimum und Maximum bestimmen >>>>>>>>>>>>>>>>>>>>>>>>>>>>>
maxi=MAXVAL(fourierfeld)
mini=MINVAL(fourierfeld)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

! Graph formatieren, ausgeben und beenden >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
!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(N/2), 0.0, REAL(N/4) , 0.0, REAL(N/2) ,0.0, REAL(N/4), mini, maxi, mini, (maxi-mini)/10)
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(fourierfeld(N/4+1:3*N/4,N/4+1:3*N/4),N/2,N/2,8,8)
CALL CRVMAT(fourierfeld,N,N,8,8)
CALL COLOR('WHITE')
CALL ENDGRF()
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

IF(.NOT.screen)  CALL DISFIN()
!_______________________________________________________________________________
END SUBROUTINE
!***********************************************************************************
!
!
!
!***********************************************************************************
SUBROUTINE Dealising(omega,kvec)
USE para
INTEGER                              :: i,k
DOUBLE PRECISION,DIMENSION(N/2+1,N,3):: kvec
DOUBLE COMPLEX  ,Dimension(N/2+1,N)  :: omega

DO k=1,N
   DO i=1,N/2+1
      IF(kvec(i,k,3).gt.REAL(N**2)/9.0) then
         omega(i,k)=(0.0,0.0)
      ELSE
      END IF
   End DO	
END DO

END SUBROUTINE
!***********************************************************************************
!
!
!
!***********************************************************************************
SUBROUTINE RechteHand(omega,out,plan_hin,plan_her,kvec)
USE para
#include'fftw3.f'
!_______________________________________________________________________________
!Rechte-Hand-Seite
INTEGER*8 plan_hin,plan_her
INTEGER :: k,i
DOUBLE COMPLEX,DIMENSION(N/2+1,N)    :: omega,omega_copy ,out  ,Nichtlinx,Nichtliny ,ux,uy
DOUBLE PRECISION,DIMENSION(N/2+1,N,3):: kvec
DOUBLE PRECISION,DIMENSION(N,N)      :: omegareal ,uyreal,uxreal ,Nichtlinrealx,Nichtlinrealy

omega_copy=omega

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

      IF(kvec(i,k,3).eq.0) THEN
        ux(i,k)=0.0d0
      ELSE
        ux(i,k)= imi*kvec(i,k,2)*omega(i,k)/kvec(i,k,3)
      END IF

      IF(kvec(i,k,3).eq.0) THEN
        uy(i,k)=0.0d0
      ELSE
        uy(i,k)=-imi*kvec(i,k,1)*omega(i,k)/kvec(i,k,3)
      END IF

   END DO
END DO
!___________________________________________________________
!Übergang omega_initial Realraum zur Berechnung der Nichtlinearität
CALL dfftw_execute_dft_c2r(plan_her,omega_copy,omegareal)
CALL dfftw_execute_dft_c2r(plan_her,ux,uxreal)
CALL dfftw_execute_dft_c2r(plan_her,uy,uyreal)
omegareal=omegareal/DBLE(N)
uxreal=uxreal/DBLE(N)
uyreal=uyreal/DBLE(N)
!_______________________________________________________________________________
DO k=1,N
   DO i=1,N
        Nichtlinrealx(i,k)= uyreal(i,k)*omegareal(i,k)
        Nichtlinrealy(i,k)=-uxreal(i,k)*omegareal(i,k)
   END DO
END DO
!_______________________________________________________________________________
!Übergang zurück omega_initial den Fourier-Raum
CALL dfftw_execute_dft_r2c(plan_hin,Nichtlinrealx,Nichtlinx)
Nichtlinx=Nichtlinx/DBLE(N)
CALL dfftw_execute_dft_r2c(plan_hin,Nichtlinrealy,Nichtliny)
Nichtliny=Nichtliny/DBLE(N)		

!_______________________________________________________________________________
CALL Dealising(Nichtlinx,kvec)
CALL Dealising(Nichtliny,kvec)

DO k=1,N
   DO i=1,N/2+1
      out(i,k)=-nu*kvec(i,k,3)*omega(i,k) + imi* (kvec(i,k,1)*Nichtliny(i,k)-kvec(i,k,2)*Nichtlinx(i,k))
   End DO	
END DO

END SUBROUTINE RechteHand
!***********************************************************************************
!
!!
!***********************************************************************************
SUBROUTINE Initial(omega_initial)
USE para
!_______________________________________________________________________________
INTEGER                        :: i,k
DOUBLE PRECISION,DIMENSION(N,N):: omega_initial
!_______________________________________________________________________________
DO k=1,N
   DO i=1,N
  !    CALL RANDOM_NUMBER(omega_initial(i,k))

      omega_initial(i,k)= &
                         amp*exp(REAL(-((i-Startwert1x)/REAL(mu*N))**2-((k-Startwert1y)/REAL(mu*N))**2))&
                        +amp*exp(REAL(-((i-Startwert2x)/REAL(mu*N))**2-((k-Startwert2y)/REAL(mu*N))**2))
!(omega_initial(i,k)-0.5d0)


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