! Compile with
! f95 -O4 -r8 -o d1 d1.f90
!
! You can see various useful plots by using gnuplot with the commands:
!
! set data style line
! plot "d1Kt.dat"
! plot "d1E.dat"
! plot "d1u.dat"
! plot "d1.dat"
!

program d1
  !
  ! simple 1-D fluid equations with Lax-Wendroff
  ! GWH, Oct. 23, 2001

  implicit none

  real, dimension(:), allocatable ::  x
  real, dimension(:), allocatable ::  den,  u , p , E , phi
  real, dimension(:), allocatable ::  den2, u2, p2, E2, phi2

  real :: dt,dx, time, Temp, nu
  real, parameter :: pi = 3.1415926535897931

  integer, parameter :: Mx=66
  integer :: j, n, nplot

  ! initializations

  allocate(x(Mx))
  allocate(den(Mx) , u(Mx) , p(Mx) , E(Mx) , phi(Mx))
  allocate(den2(Mx), u2(Mx), p2(Mx), E2(Mx), phi2(Mx))
  den=0; u=0; p=0; E=0; phi=0 ; time=0.0 ; Temp=1; nu=0.0

  open(unit=20, file="d1.dat", status="unknown")
  write(20,*) "# x den(x)  (with blank lines separating different times)"

  open(unit=23, file="d1Kt.dat", status="unknown")
  write(20,*) "# t <u**2>(t)"

  do j=1,Mx
     x(j)=(j-1)
  enddo

!  den=-(x-Mx/2)*exp(-(x-Mx/2)**2/(Mx/2/8)**2)/abs(Mx/2/8)
!  den=den-volavg(den)

  u=-(x-Mx/2)*exp(-(x-Mx/2)**2/(Mx/2/8)**2)/abs(Mx/2/8)

  call Esolve

  dx=1; dt=0.001
  nplot=1.0/dt  ! make plots every t=1.0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Time loop

  do n=1,10000000
     call step
     time=time+dt
     if(mod(n,nplot) .eq. 0) then
        do j=1,Mx-1
           write(20,*) x(j), den(j)
        enddo
        write(20,*)
        write(23,*) volavg(u*u)
     endif
     if(time .gt. 2*pi*10) exit
  enddo

  print *, "Completed to time ", time

  open(unit=21,file="d1u.dat",status="unknown")
  write(21,*) "# x    u(x)"
  do j=1,Mx
     write(21,*) x(j), u(j)
  enddo

  open(unit=22,file="d1E.dat",status="unknown")
  write(22,*) "# x   E(j)"
  do j=1,Mx
     write(22,*) x(j), E(j)
  enddo

!
! End of Main program block
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine step

  ! First step of Lax-Wendroff
    do j=2,Mx-1
       
!      den2(j)= (den(j)+den(j+1))/2 -dt/2/dx*(u(j+1)-u(j))
!      u2(j)=   (u(j)+u(j+1))/2     -dt/2/dx*Temp*(den(j+1)-den(j)) -dt/2*E(j)

       den2(j)= (den(j)+den(j+1))/2 -dt/2/dx*(u(j+1)-u(j))
       u2(j)=   (u(j)+u(j+1))/2     -dt/2/dx*(p(j+1)-p(j)) -dt/2*E(j)
       p2(j)=   (p(j)+p(j+1))/2     -dt/2/dx*3*(u(j+1)-u(j)) &
            -dt/2*nu*(p(j)-den(j))
       
    enddo
    
    ! Periodic B.C.'s:
    den2(Mx)= den2(2)
    u2(Mx)=   u2(2)
    p2(Mx)=   p2(2)
    den2(1)= den2(Mx-1)
    u2(1)=   u2(Mx-1)
    p2(1)=   p2(Mx-1)

    call Esolve2
    
    ! Second step of Lax-Wendroff
    do j=2,Mx-1

!       den(j)= den(j) -dt/dx*(u2(j)-u2(j-1))
!       u(j)= u(j) -dt/dx*Temp*(den2(j)-den2(j-1)) -dt*E(j)

       den(j)= den(j) -dt/dx*(u2(j)-u2(j-1))
       u(j)=     u(j) -dt/dx*(p2(j)-p2(j-1)) -dt*E(j)
       p(j)=     p(j) -dt/dx*3*(u2(j)-u2(j-1)) -dt*nu*(p2(j)-den2(j))

    enddo

    ! Periodic B.C.'s:
    den(Mx)=den(2)
    u(Mx)=    u(2)
    p(Mx)=    p(2)
    den(1)= den(Mx-1)
    u(1)=     u(Mx-1)
    p(1)=     p(Mx-1)

    call Esolve
    
  end subroutine step

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine Esolve
    ! Solve for E, the electric field
    ! This is solving for E(x(j+1/1)), half way 
    ! in between grid points, which we denote
    ! as E(j) here.
    E(1)=0.0
    do j=2,Mx-1
       E(j)=E(j-1)-dx*den(j)
    enddo
    E(Mx)=E(2)
    E=E-volavg(E)
    ! The requirement that volavg(E)=0 comes from
    ! restriction that phi be periodic, 
  end subroutine Esolve

  subroutine Esolve2
    ! Solve for E, the electric field
    ! This is solving for E(x(j))=E(j), evaluated
    ! on the grid points.
    E(1)=0.0
    do j=2,Mx-1
       E(j)=E(j-1)-dx*den2(j-1)
    enddo
    E=E-volavg(E)
  end subroutine Esolve2

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  real function volavg(a)
    ! Calculate the volume average of the field "a".
    real :: a(:)
    volavg=0
    do j=2,Mx-1
       volavg=volavg+a(j)
    enddo
    volavg=volavg/(Mx-2)
  end function volavg

end program d1