PROGRAM uebung
!*********************************************************************
!
! Description:
!  Sheet 5, Exercise 1
!
!**********************************************************************

! mo_output provides a simplified interface for writing netCDF data
!     see supplement.pdf for details
  USE mo_output,   ONLY:  open_out_file,   &
                          init_out,        &
                          write_out,       &
                          close_out_file

  IMPLICIT NONE

! 1. variable declaration-----------------------------------

   ! number of grid-points in x-direction
  INTEGER, PARAMETER :: jmax = 201
   ! number of time steps
  INTEGER, PARAMETER ::  nmax = 120
   ! grid spacing in x-direction
  DOUBLE PRECISION, PARAMETER :: dx   = 5.D3
   ! length of time step
  DOUBLE PRECISION, PARAMETER :: dt   = 10.0D0
   ! velocity
 ! DOUBLE PRECISION, PARAMETER :: c    = 0.5D0
   ! mean depth of fluid
  DOUBLE PRECISION, PARAMETER :: H    = 1.D3
   ! gravitaional acceleration
  DOUBLE PRECISION, PARAMETER :: g    = 9.81D0

   ! netcdf output unit
  INTEGER :: ncid
    ! "time step" index
  INTEGER :: nstep
    ! the function (3 time level scheme)
  DOUBLE PRECISION, POINTER, DIMENSION(:) :: ufct
  DOUBLE PRECISION, POINTER, DIMENSION(:) :: ufct_old
  DOUBLE PRECISION, POINTER, DIMENSION(:) :: ufct_new
  DOUBLE PRECISION, POINTER, DIMENSION(:) :: hfct
  DOUBLE PRECISION, POINTER, DIMENSION(:) :: hfct_old
  DOUBLE PRECISION, POINTER, DIMENSION(:) :: hfct_new

  INTEGER :: j
  DOUBLE PRECISION :: fact

!-end variable declaration --------------------------------

    ! write to screen:
  WRITE(6,*) 'This is uebung.exe'

  ! 1. allocate (add ghost points at the boundaries)
   ALLOCATE( ufct(0:jmax+3))
   ALLOCATE( ufct_old(0:jmax+1))
   ALLOCATE( ufct_new(0:jmax+1))
   ALLOCATE( hfct(0:jmax+3))
   ALLOCATE( hfct_old(0:jmax+1))
   ALLOCATE( hfct_new(0:jmax+1))


  ! 2. initialize the function
    DO j=0, jmax+1
     IF ( DBLE(j-1)*dx .GE. 450.D3 .AND. DBLE(j-1)*dx .LE. 550.D3 ) THEN
       hfct(j) = 1.D0
     ELSE
       hfct(j) = 0.D0
     END IF
    END DO

print*, hfct(91)

   ! 4. open a netCDF output file
  CALL open_out_file(ncid, "out.nc")

   ! 5. add the variable "ufct" to a list of output variables
  CALL init_out(ncid, ufct, "ufct", "none","x", timedep=.TRUE.)
  CALL init_out(ncid, hfct, "hfct", "none","x", timedep=.TRUE.)

   ! 6. write all the variables in the list of output variables to the output file
   !   here: initial values
  CALL write_out(ncid)

   ! 7. start time loop
  DO nstep=1,nmax

    ! periodic boundary cond.

     ufct(0)=ufct(jmax)
     ufct(jmax+1)=ufct(1)
     hfct(0)=hfct(jmax)
     hfct(jmax+1)=hfct(1)

   ! first step should be Euler forward
       fact=1.D0
     IF ( nstep .EQ. 1 ) THEN
       ufct_old = ufct(1:jmax)
       hfct_old = hfct(1:jmax)
       fact=0.5D0
     END IF

   ! calculate values at  step n+1 ACHTUNG: Halb-Schritte j beginnt also bei 1/2
    DO j = 1, jmax
     ufct_new(j) = ufct_old(j) - fact * dt * g * ( hfct(j+3) - hfct(j-1) ) / dx
     hfct_new(j) = hfct_old(j) - fact * dt * H * ( ufct(j+1) - ufct(j-1) ) / dx
    END DO

    ! apply the Robert Asselin filter
    DO j = 1, jmax
     ufct_old(j) = ufct(j) + 0.1d0 * (ufct_old(j) - 2.d0* ufct(j) + ufct_new(j))
     hfct_old(j) = hfct(j) + 0.1d0 * (hfct_old(j) - 2.d0* hfct(j) + hfct_new(j))
    END DO

    ! update ufct
      ufct(1:jmax-100)=ufct_new(2:jmax+1:2)
      hfct(1:jmax-100)=hfct_new(2:jmax+1:2)
 print*, nstep, hfct(45)
    ! write output
    CALL write_out(ncid)

       ! update ufct
      ufct(1:jmax)=ufct_new
      hfct(1:jmax)=hfct_new

 END DO


   ! 8 close the output file
  CALL close_out_file(ncid)

   ! 9. free memory
  DEALLOCATE(ufct)
  DEALLOCATE(ufct_new)
  DEALLOCATE(ufct_old)
  DEALLOCATE(hfct)
  DEALLOCATE(hfct_new)
  DEALLOCATE(hfct_old)

  WRITE(6,*) 'uebung.exe: DONE.'

END PROGRAM uebung