new file mode 100644
--- /dev/null
+++ b/waves.f95
@@ -0,0 +1,250 @@
+! ******************** PROGRAM EXPLANATION ***************************
+!
+! Purpose: to model a water surface using the shallow water equations under
+!          the effect of several splashes.
+!
+! This file contains the following (in order):
+!   -The module "wavesmod" which contains:
+!      - Constants and parameters which control the behaviour of the program
+!         (eg number of timesteps and length of timestep)
+!      - The subroutine calcu which calculates the horizontal water velocity u
+!         from the gradient of the water height
+!      - The subroutine calcv which calculates the horizontal water velocity v
+!         from the gradient of the water heighte
+!      - The subroutine calch which calculates the water height h
+!         from the convergence and divergence of u and v 
+! 
+!  -The main program which contains:
+!      - Main program variable initialisation and memory allocation
+!      - Specification of the location and timing of splashes
+!      - Opening output file
+!      - The main model loop including model output to file
+
+
+! *********************** WAVESMOD ***********************************
+
+MODULE wavesmod
+!
+! Purpose:
+!   To share data and subroutines within the waves.f95 program
+!
+IMPLICIT NONE
+SAVE
+
+! Declare constants - Change any of these if you want the program to behave differently
+REAL, PARAMETER :: g = 9.81				! Gravity (m s-2)
+REAL, PARAMETER :: PI = 3.14159			! PI
+INTEGER, PARAMETER :: numboxes = 200 	! Number of grid boxes wide the domain is
+REAL, PARAMETER :: poolwidth = 0.2		! Pool width (m)
+REAL, PARAMETER :: drag = 2.0			! Viscous drag (s-1)
+REAL, PARAMETER :: dt = 0.0001			! Timestep (s)
+REAL, PARAMETER :: depth = 0.01			! Pool depth (m)
+INTEGER, PARAMETER :: ntimesteps = 30000! The total number of timesteps that the program will run for
+INTEGER, PARAMETER :: numsplashes = 6	! The number of splashes - The splash locations are defined in the main program
+REAL, PARAMETER :: sprad = 0.03			! Splash radius
+REAL, PARAMETER :: spheight = 0.01		! Splash height
+INTEGER, PARAMETER :: numtomiss = 100	! How timesteps do we skip before writing to file
+
+! Initialise derived constants
+REAL :: dx								! Distance between grid boxes in x dimension (m)
+REAL :: dy								! Distance between grid boxes in y dimension (m)
+
+! Initialise arrays
+REAL, ALLOCATABLE, DIMENSION(:,:) :: u0 ! column water speed in x direction at previous timestep (m s-1)
+REAL, ALLOCATABLE, DIMENSION(:,:) :: u1 ! column water speed in x direction at next timestep (m s-1)
+REAL, ALLOCATABLE, DIMENSION(:,:) :: v0 ! column water speed in y direction at previous timestep (m s-1)
+REAL, ALLOCATABLE, DIMENSION(:,:) :: v1 ! column water speed in y direction at next timestep (m s-1)
+REAL, ALLOCATABLE, DIMENSION(:,:) :: h0 ! column water height at previous timestep (m)
+REAL, ALLOCATABLE, DIMENSION(:,:) :: h1 ! column water height at next timestep (m)
+
+! Include subroutines
+CONTAINS
+
+! ************************* CALCU ******************************************
+SUBROUTINE calcu
+!
+! Purpose:
+!   To calculate u at the next timestep using the u shallow water equation
+!
+IMPLICIT NONE
+
+! Initialise arrays
+REAL, DIMENSION(numboxes-1,numboxes) :: hxp1   ! column water height at grid boxes at x plus 1 (m)
+REAL, DIMENSION(numboxes-1,numboxes) :: dh	   ! change in water height between grid boxes in x direction (m)
+
+u1=0.0
+
+! Calculate the change in p in the x dimension
+hxp1 = h0(2:,:)
+dh = hxp1-h0(:numboxes-1,:)
+
+! Run shallow water equation
+u1(2:numboxes,:) = (-g*dh/dx - u0(2:numboxes,:)*drag) * dt + u0(2:numboxes,:)
+
+END SUBROUTINE calcu
+
+! ************************* CALCV ******************************************
+SUBROUTINE calcv
+!
+! Purpose:
+!   To calculate v at the next timestep using the v shallow water equation
+!
+IMPLICIT NONE
+
+! Initialise arrays
+REAL, DIMENSION(numboxes,numboxes-1) :: hyp1   ! column water height at grid boxes at y plus 1 (m)
+REAL, DIMENSION(numboxes,numboxes-1) :: dh	   ! change in water height between grid boxes in y direction (m)
+
+v1=0.0
+
+! Calculate the change in h in the y dimension
+hyp1 = h0(:,2:)
+dh = hyp1-h0(:,:numboxes-1)
+
+! Run shallow water equation
+v1(:,2:numboxes) = (-g*dh/dy - v0(:,2:numboxes)*drag) * dt + v0(:,2:numboxes)
+
+END SUBROUTINE calcv
+
+! ************************ CALCH ******************************************
+SUBROUTINE calch
+!
+! Purpose:
+!   To calculate h at the next timestep as an average of convergence between the last and this time point
+!
+IMPLICIT NONE
+
+! Initialise arrays
+REAL, DIMENSION(numboxes,numboxes) :: uxp1   ! time average u at grid boxes x plus 1 (m2 s-2)
+REAL, DIMENSION(numboxes,numboxes) :: vyp1   ! time average v at grid boxes y plus 1 (m2 s-2)
+REAL, DIMENSION(numboxes,numboxes) :: du	 ! change in u between grid boxes in the x dimension (m2 s-2)
+REAL, DIMENSION(numboxes,numboxes) :: dv	 ! change in v between grid boxes in the y dimension (m2 s-2)
+
+! Calcualate the change in ubar in the x dimension and the change in vbar in the y dimension
+uxp1 = u1(2:numboxes+1,:)
+vyp1 = v1(:,2:numboxes+1)
+du = uxp1-u1(:numboxes,:)
+dv = vyp1-v1(:,:numboxes)
+
+h1 = -depth*(du/dx + dv/dy)*dt + h0
+
+END SUBROUTINE calch
+
+END MODULE wavesmod
+
+! ************************* MAIN PROGRAM *****************************
+
+PROGRAM waves
+! Purpose:
+!  To model waves in a pool
+USE wavesmod  ! Use wavesmod - including array declarations and procedures
+IMPLICIT NONE
+
+! Declare variables
+INTEGER :: t		! The timestep number
+REAL :: time		! The time since the program started
+INTEGER :: status	! The status of memory allocation procedures
+INTEGER :: loopcount	! How many times through the loop between
+                        !  time counts are you 
+REAL, DIMENSION(3,numsplashes) :: splashes ! Where shall we make the
+                                           !  splashes (m,m,s) 
+REAL, DIMENSION(numboxes,numboxes) :: sp  ! The splash profile
+REAL, DIMENSION(numboxes) :: x,y	  ! The x and y dimensions
+REAL :: spdist      ! The distance from the centre of the splash
+INTEGER :: i,j	    ! Indices used in loops
+INTEGER :: splashcount=0	! How many splashes have we had
+INTEGER :: Id=1		!  File ID number
+CHARACTER(len=3) :: numboxes_str  ! The number of boxes converted to a string
+CHARACTER(len=10) :: fmt,dim_fmt ! The format to output to file
+INTEGER :: noutputcounts  ! How many time outputs are there
+
+! Allocate memory to the arrays
+ALLOCATE(u0(numboxes+1,numboxes), u1(numboxes+1,numboxes), v0(numboxes,numboxes+1), &
+         v1(numboxes,numboxes+1), h0(numboxes,numboxes), h1(numboxes,numboxes), STAT=status)
+IF (status /= 0) THEN
+  WRITE(*,*) 'Error allocating memory for arrays!'
+END IF
+
+! Calculate derived constants
+dx = poolwidth/real(numboxes) 	
+dy = poolwidth/real(numboxes)
+noutputcounts = ntimesteps/numtomiss
+
+! Calculate the x and y dimensions
+DO i=1,numboxes
+	x(i)=dx*(i-1)
+    y(i)=dy*(i-1)
+END DO
+
+! Declare where and when the splashes occur
+! They are in the format: (/ x location (m), y location (m), splash time (s) /)
+! with each row of the array containing a different splash.
+! The last splash is a dummy splash that should never be executed but should remain
+! there with a large splash time for the program to work.
+splashes(1:3,1) = (/ 0.04,0.04,0.04 /)  
+splashes(1:3,2) = (/ 0.04,0.04,0.54 /) 
+splashes(1:3,3) = (/ 0.04,0.04,1.04 /) 
+splashes(1:3,4) = (/ 0.04,0.04,1.54 /) 
+splashes(1:3,5) = (/ 0.04,0.04,2.04 /) 
+splashes(1:3,6) = (/ 0.04,0.04,400.04 /) 
+
+! Initialise arrays
+u0=0.0
+v0=0.0
+h0=0.0
+loopcount=0
+
+! Make a format variable for describing data input to files
+numboxes_str = ' '
+write (numboxes_str,'(I3)') numboxes
+dim_fmt= '(' // numboxes_str // 'F10.3)'
+fmt = '(' // numboxes_str // 'F15.8)'
+
+! Write some constants needed for file processing
+OPEN (UNIT=Id, FILE='waves.dat', STATUS='REPLACE')
+
+! The main program loop
+DO t = 1,ntimesteps
+  time = real(t)*dt
+
+  ! Add a splash if required
+  if (time >= splashes(3,splashcount+1)) then
+    sp = 0.0
+    splashcount=splashcount+1
+  	DO i=1,numboxes
+      DO j=1,numboxes
+  		spdist = ((x(i)-splashes(1,splashcount))**2 + (y(j)-splashes(2,splashcount))**2)**0.5
+        if (spdist < sprad/3.0) then	! In the ridge
+          sp(i,j) = spheight * cos(3.0*PI*spdist/(2.0*sprad))
+        endif
+      END DO
+    END DO
+    h0 = h0 + sp
+  end if
+
+  ! Make the timestep calculations
+  CALL calcu
+  CALL calcv
+  CALL calch
+  
+  ! Write to an output file every 100th timestep
+  loopcount = loopcount + 1
+  IF (loopcount >= numtomiss) THEN 
+    WRITE(*,*) 'Time=', time
+    do j=1,numboxes
+       WRITE (Id,fmt) h1(1:numboxes,j)
+    enddo
+    loopcount=0
+  END IF
+    
+  ! Shift the timestep so that h0 contains h1, u0 contains u1 and v0 contains v1
+  h0=h1
+  u0=u1
+  v0=v1
+
+END DO
+
+! Finish up.
+END PROGRAM waves
+
+