Shallow Water Equations

Ref: http://en.wikipedia.org/wiki/Shallow_water_equations

The equations are

({{(Du/Dt) - f v, =, - g ∂η/∂x - b u}, {(Dv/Dt) + f u, =, - g ∂η/∂y - b v}, {((∂η)/(∂t)), =, -((∂[u(H+η)])/(∂x))-((∂[v(H+η)])/(∂y))}})

The Coriolis force term is negligible so f=0.

({{(Du/Dt) , =, - g ∂η/∂x - b u}, {(Dv/Dt) , =, - g ∂η/∂y - b v}, {((∂η)/(∂t)), =, -((∂[u(H+η)])/(∂x))-((∂[v(H+η)])/(∂y))}})

If everything is independent of y then

({{(Du/Dt) , =, - g ∂η/∂x - b u}, {(Dv/Dt) , =, - b v}, {((∂η)/(∂t)), =, -((∂[u(H+η)])/(∂x))}})

and we can have a consistent solution with v=0.

({{(Du/Dt) , =, - g ∂η/∂x - b u}, {((∂η)/(∂t)), =, -((∂[u(H+η)])/(∂x))}})

Expanding the material time derivative (D/Dt) we get the coupled system.

({{((∂u)/(∂t))+u ((∂u)/(∂x)) , =, - g ∂η/∂x - b u}, {((∂η)/(∂t)), =, -((∂[u(H+η)])/(∂x))}})

An explicit  approximation for this (with (Subscript[𝒲, j]^n)=((Subscript[ℋ, j])+(Subscript[η, j]^n))) is

({{(((Subscript[𝒰, j]^(n+1))-(Subscript[𝒰, j]^n))/(Δt))+(Subscript[𝒰, j]^n) (((Subscript[∇, x])(Subscript[𝒰, j]^n))/(Δx)) , =, - g (((Subscript[∇, x])(Subscript[η, j]^n))/(Δx)) - b (Subscript[𝒰, j]^n)}, {(((Subscript[η, j]^(n+1))-(Subscript[η, j]^n))/(Δt)), =, -((Subscript[𝒲, j]^n)(((Subscript[∇, x])(Subscript[𝒰, j]^n))/(Δx))+(Subscript[𝒰, j]^n)(((Subscript[∇, x])(Subscript[𝒲, j]^n))/(Δx)))}})

I need initial conditions giving (Subscript[𝒰, j]^1) and (Subscript[η, j]^1).  With the differencing I also need boundary conditions giving (Subscript[𝒰, 1]^n) and (Subscript[η, 1]^n):  I guess the simplest thing would be to use zero conditions for these.  Finally I need the value for g=9.8 m/(s^2) and a profile (Subscript[ℋ, j]) for the ocean bottom. This is relatively easy to implement if we rewrite it as

({{(Subscript[𝒰, j]^(n+1)), =, (1-b Δt)(Subscript[𝒰, j]^n)-((Δt)/(Δx))(Subscript[𝒰, j]^n) (Subscript[∇, x])(Subscript[𝒰, j]^n)- g ((Δt)/(Δx)) (Subscript[∇, x])(Subscript[η, j]^n)}, {(Subscript[η, j]^(n+1)), =, (Subscript[η, j]^n)-((Δt)/(Δx))((Subscript[𝒲, j]^n)(Subscript[∇, x])(Subscript[𝒰, j]^n)+(Subscript[𝒰, j]^n)(Subscript[∇, x])(Subscript[𝒲, j]^n))}})

DMinus[VecIn_]:=Module[
{VecTemp},
VecTemp=VecIn-RotateRight[VecIn];
(* VecTemp[[1]]*=0.0; *)
VecTemp
]
DPlus[VecIn_]:=Module[
{VecTemp},
VecTemp=RotateLeft[VecIn]-VecIn;
(* VecTemp[[-1]]*=0.0; *)
VecTemp
]
DUpWind[VecIn_, WindIn_]:= Module[
{Indicator, D1, D2, DOut},
Indicator=Positive[WindIn]/.{
True->1,
False->0
};
D1=DMinus[VecIn];
D2=DPlus[VecIn];
DOut=D1*Indicator+D2*(1-Indicator);
DOut[[-1]]*=0.0; 
DOut
]

Checking

{xMin,xMax,Δx}={0.0,10.0,0.1};
mid=5.0;
{a, Δt}={1.0,0.05};
u0=Table[ Sech[x-mid],{x, xMin,xMax,Δx}]-0.5;
u=u0;
Do[
u=u-a ((Δt)/(Δx))*u*DUpWind[u,u],
{23}]
ListPlot[{u0,u}]

Explicit Code

All that remains is to add the boundary conditions.  We want to have a wave arrive at the "x=0" end of the array.  In some sense the simplest way to do this is to just add this as a specific value at the end of the update stage. I am going to add an elevation (Subscript[η, 1]) and a zero velocity (Subscript[u, 1])=0.

({{(Subscript[𝒰, j]^(n+1)), =, (1-b Δt)(Subscript[𝒰, j]^n)-((Δt)/(Δx))(Subscript[𝒰, j]^(n+1)) ((Subscript[𝒰, j]^n)-(Subscript[𝒰, j-1]^n))- g ((Δt)/(Δx)) ((Subscript[η, j]^n)-(Subscript[η, j-1]^n))}, {(Subscript[η, j]^(n+1)), =, (Subscript[η, j]^n)-((Δt)/(Δx))((Subscript[𝒲, j]^n)(((Subscript[∇, x])(Subscript[𝒰, j]^n))/(Δx))+(Subscript[𝒰, j]^n)(((Subscript[∇, x])(Subscript[𝒲, j]^n))/(Δx)))}})

Clear[ExplicitStepSWEq]
ExplicitStepSWEq[ {UIn_,ηIn_},Δt_,Δx_,{b_,HIn_},ηBC_]:= Module[
{g,rat,WTemp,wDiff,uDiff1,uDiff2,ηDiff, uOut, ηOut},
g=9.8;
rat=Δt/Δx;
WTemp= (HIn +ηIn) ;
wDiff= DUpWind[WTemp,UIn];
uDiff1=DUpWind[UIn,UIn];
uDiff2=DUpWind[UIn,WTemp];
ηDiff=DMinus[ηIn];
uOut=(1-b*Δt)UIn-rat*( UIn*uDiff1+ g*ηDiff);
ηOut=ηIn - rat*(WTemp*uDiff2+UIn*wDiff);
uOut[[1]]=0.0;
ηOut[[1]]=ηBC;

{uOut,ηOut}
]

{Δt,Δx}={0.05,0.2};
xMax=20;
b=0.1;
g=9.8;
H=Table[ 10-0.1x,{x, 0.0, xMax, Δx}];
H=1+ 0.0 H;
U0=0.0*H;
η0=Table[ Sech[x-10],{x, 0.0, xMax,Δx}];
η=0.0η0;
U=U0;
Nt=200;
n=1;
While[
And[n<Nt,Max[Abs[U]+Abs[η]]<100], 
(*
Print[
StringForm["n=``: \nStability Criteria are: \nu*Δt/Δx=``, g*Δt/Δx=``, w*Δt/Δx=``",
n,Max[Abs[U]]*Δt/Δx, g*Δt/Δx, Max[U+H]*Δt/Δx]
]; 
*)
n++;
{U,η}=ExplicitStepSWEq[{U,η},Δt, Δx,{b,H},0.2 Sin[n Δt]]
]
Print[{n*Δt, Nt*Δt}]

Table[ {n Δt,0.2 Sin[n Δt]},{n, 1,200}]//ListPlot

ListPlot[η,
PlotRange->All,
PlotLabel->"η",
PlotJoined->True
]
ListPlot[U,
PlotRange->All,
PlotLabel->"U",
PlotJoined->True]

ListPlot[{H,H+η0},
PlotRange->{0,12}
]

Allan's New Stuff

Matrices and Sparse Matrices and Speed (Oh My!)
Differencing and Sparse Matrices

I am going to implement a few difference operators and give them semi-meaningful names.  First I am going to assume that I have a table of values that are naturally arranged in a grid of x,y values.  I am going to call the grid of values U[[i,j]] and we are going to assume that U[[i,j]]≈u[(Subscript[x, i]),(Subscript[y, j])] for some function u.

AδPlusX[InMat_]:= Module[
{Blah,i},
Table[
InMat[[i+1]]-InMat[[i]],
{i,1, Length[InMat]-1}
]
]
AδMinusX[InMat_]:= Module[
{Blah,i},
Table[
InMat[[i]]-InMat[[i-1]],
{i,2, Length[InMat]}
]
]

or ...

AδPlusX[InMat_]:= Module[
{Blah,i},
Table[
InMat[[i+1]]-InMat[[i]],
{i,1, Length[InMat]-1}
]
]
AδMinusX[InMat_]:= Module[
{Blah,i},
Table[
InMat[[i]]-InMat[[i-1]],
{i,2, Length[InMat]}
]
]

Or (probably better and measurably faster)

δPlusX[InMat_]:= Module[
{Blah},
RotateLeft[InMat]-InMat
]
δMinusX[InMat_]:= Module[
{Blah,i},
InMat-RotateRight[InMat]
]

What about the y-derivatives

δPlusY[InMat_]:= Module[
{Temp},
Temp=Transpose[InMat];
Transpose[δPlusX[Temp]]
]
δMinusY[InMat_]:= Module[
{Temp},
Temp=Transpose[InMat];
Transpose[δMinusX[Temp]]
]

U=Table[u[x,y],{x,0, 10,0.01},{y,2,3,0.02}];
Dimensions[U]
Timing[AδPlusX[U];]
Timing[δPlusX[U];]
Timing[δPlusY[U];]

Of course we can use matrices. (first the version that does not wrap around!)

δPlusXMat=SparseArray[
Flatten[
Table[ {
{i,i}->1,
{i,i+1}->-1
},
{i, Length[U]-1}]
]
]
Dimensions[U]
Timing[Dimensions[δPlusXMat.U]]
NonSparseδ=Normal[δPlusXMat];
Timing[Dimensions[NonSparseδ.U]]

The conclusion is that you should not use full matrices for sparse things!    Now a version that does wrap around but without the full matrix computation!

δMinusXMat=SparseArray[
Flatten[
{
{{1,1}->-1,{1,Length[U]}->1},
Table[ {
{i,i-1}->1,
{i,i}->-1
},
{i, 2,Length[U]}]
}
]
]
Dimensions[U]
Timing[Dimensions[δMinusXMat.U]]

You can do the same thing with the y derivatives but there is not much point since it is slower than the "rotate" version.

Sparse vs Full Matrices

The question came up why not just use sparse matrices.

n=1000;
A= Table[RandomReal[],{n},{n}];
AA=SparseArray[A]
B=Table[RandomReal[],{n},{n}];
BB=SparseArray[B]
Timing[A.B;]
Timing[AA.BB;]

The memory storage is also worse.

ByteCount[A]
ByteCount[AA]

But for genuinely sparse matrices the sparse matrix stuff is a bit faster.

n=(10^4);
i=1;
SparseSample=SparseArray[
Flatten[
{
Table[{k,k}->RandomReal[{10, 100}],
{k,1,n}],
Table[
i=RandomInteger[{i,n-1}];
{i, RandomInteger[{i+1,n}]}->RandomReal[{-1,1}],
{100}
]
}
],
{n,n}
]
b= Table[RandomReal[{0,1}],{n}];
Timing[LinearSolve[SparseSample,b];]
Timing[LinearSolve[Normal[SparseSample],b];]

Old Stuff

Discretization

Note that:
u is the velocity in the x dimension,
v is the velocity in the y dimension,
H is the mean height of the horizontal pressure surface,
ψ is the deviation of the horizontal pressure surface from its mean,
g is the acceleration of gravity,
f is the term corresponding to the Coriolis force,
b is the viscous drag.
We will be considering x-dimensions where 0≤x≤1mile and y-dimensions where 0≤y≤(100/5280)miles.

First Equation

We will deal with the first equation from wikipedia.org.
(Du/Dt)-f v=-g ((∂ψ)/(∂x))-b u     =>
((∂u)/(∂t))-f v=-g ((∂ψ)/(∂x))-b u .
For now, let f=0  and  (* b=1*)  and α=-(( Δt)/(Δx)).
((∂u)/(∂t))-f v=-g ((∂ψ)/(∂x))-b u     =>
((∂u)/(∂t))=-g ((∂ψ)/(∂x))-u       =>
(((Subscript[U, J,K]^(n+1))-(Subscript[U, J,K]^n))/(Δt))=((-g)/(Δx))((Subscript[Ψ, J+1,K]^n)-(Subscript[Ψ, J,K]^n))-b (Subscript[U, J,K]^n)     =>
(Subscript[U, J,K]^(n+1))-(Subscript[U, J,K]^n)=α g((Subscript[Ψ, J+1,K]^n)-(Subscript[Ψ, J,K]^n))-Δt b (Subscript[U, J,K]^n)    =>
(Subscript[U, J,K]^(n+1))=α g((Subscript[Ψ, J+1,K]^n)-(Subscript[Ψ, J,K]^n))-(Subscript[U, J,K]^n)(Δt b-1)

Second Equation

Now we will consider a new function:
(Dv/Dt)-f u=-g ((∂ψ)/(∂y))-b v     =>
((∂v)/(∂t))-f u=-g ((∂ψ)/(∂y))-b v .
For now, let f=0  and  b=1  and β=-(( Δt)/(Δy)).
((∂v)/(∂t))-f u=-g ((∂ψ)/(∂y))-b v     =>
((∂v)/(∂t))=-g ((∂ψ)/(∂y))-b v       =>
(((Subscript[V, J,K]^(n+1))-(Subscript[V, J,K]^n))/(Δt))=((-g)/(Δy))((Subscript[Ψ, J,K+1]^n)-(Subscript[Ψ, J,K]^n))-b (Subscript[V, J,K]^n)     =>
(Subscript[V, J,K]^(n+1))-(Subscript[V, J,K]^n)=β g((Subscript[Ψ, J,K+1]^n)-(Subscript[Ψ, J,K]^n))-Δt b (Subscript[V, J,K]^n)    =>
(Subscript[V, J,K]^(n+1))=β g((Subscript[Ψ, J,K+1]^n)-(Subscript[Ψ, J,K]^n))-(Subscript[V, J,K]^n)(Δt b-1).

{nx,ny,nt}={85,85,500};
Hmm=Table[
RandomReal[], 
{nx+1},
{ny+1},
{234}];
ByteCount[Hmm]/(10^6.)

Final Equation

Making One Equation

We will now combine the first two equations with the final one:

(Overscript[(Subscript[Ψ, J]^K), n+1])=α (Subscript[H, J]^K) (Overscript[(Subscript[U, J+1]^K), n])+(Overscript[(Subscript[U, J]^K), n])(α (Overscript[(Subscript[Ψ, J+1]^K), n])-α (Subscript[H, J]^K) -Δt)+β (Subscript[H, J]^K)(Overscript[(Subscript[V, J]^(K+1)), n])+β(Overscript[(Subscript[V, J]^K), n])((Overscript[(Subscript[Ψ, J]^(K+1)), n])-(Subscript[H, J]^K))+(Overscript[(Subscript[Ψ, J]^K), n])(1+β((Overscript[(Subscript[V, J]^(K+1)), n])-2 (Overscript[(Subscript[V, J]^K), n]))+α((Overscript[(Subscript[U, J+1]^K), n])-2 (Overscript[(Subscript[U, J]^K), n])))      =>
(Overscript[(Subscript[Ψ, J]^K), n+1])=α (Subscript[H, J]^K) (Overscript[(Subscript[U, J+1]^K), n])+(Overscript[(Subscript[U, J]^K), n])(α ((((Overscript[(Subscript[U, J]^K), n+1])+(Overscript[(Subscript[U, J]^K), n])(Δt-1))/(α g))+(Overscript[(Subscript[Ψ, J]^K), n]))-α (Subscript[H, J]^K) -Δt)+β (Subscript[H, J]^K)(Overscript[(Subscript[V, J]^(K+1)), n])+β(Overscript[(Subscript[V, J]^K), n])((((Overscript[(Subscript[V, J]^K), n+1])+(Overscript[(Subscript[V, J]^K), n])(Δt-1))/(β g))+(Overscript[(Subscript[Ψ, J]^K), n])-(Subscript[H, J]^K))+(Overscript[(Subscript[Ψ, J]^K), n])(1+β((Overscript[(Subscript[V, J]^(K+1)), n])-2 (Overscript[(Subscript[V, J]^K), n]))+α((Overscript[(Subscript[U, J+1]^K), n])-2 (Overscript[(Subscript[U, J]^K), n])))     =>
(Overscript[(Subscript[Ψ, J]^K), n+1])=α (Subscript[H, J]^K) (Overscript[(Subscript[U, J+1]^K), n])+(Overscript[(Subscript[U, J]^K), n])((((Overscript[(Subscript[U, J]^K), n+1])+(Overscript[(Subscript[U, J]^K), n])(Δt-1))/g)+α (Overscript[(Subscript[Ψ, J]^K), n])-α (Subscript[H, J]^K) -Δt)+β (Subscript[H, J]^K)(Overscript[(Subscript[V, J]^(K+1)), n])+(Overscript[(Subscript[V, J]^K), n])((((Overscript[(Subscript[V, J]^K), n+1])+(Overscript[(Subscript[V, J]^K), n])(Δt-1))/g)+(Overscript[β (Subscript[Ψ, J]^K), n])-β (Subscript[H, J]^K))+(Overscript[(Subscript[Ψ, J]^K), n])(1+β((Overscript[(Subscript[V, J]^(K+1)), n])-2 (Overscript[(Subscript[V, J]^K), n]))+α((Overscript[(Subscript[U, J+1]^K), n])-2 (Overscript[(Subscript[U, J]^K), n])))     =>
(Overscript[(Subscript[Ψ, J]^K), n+1])=α (Subscript[H, J]^K) (Overscript[(Subscript[U, J+1]^K), n])+(Overscript[(Subscript[U, J]^K), n])((((Overscript[(Subscript[U, J]^K), n+1])+(Overscript[(Subscript[U, J]^K), n])(Δt-1))/g)-α (Subscript[H, J]^K) -Δt)+α (Overscript[(Subscript[U, J]^K), n]) (Overscript[(Subscript[Ψ, J]^K), n])+β (Subscript[H, J]^K)(Overscript[(Subscript[V, J]^(K+1)), n])+(Overscript[(Subscript[V, J]^K), n])((((Overscript[(Subscript[V, J]^K), n+1])+(Overscript[(Subscript[V, J]^K), n])(Δt-1))/g)-β (Subscript[H, J]^K))+β(Overscript[(Subscript[V, J]^K), n])  (Overscript[(Subscript[Ψ, J]^K), n])+(Overscript[(Subscript[Ψ, J]^K), n])(1+β((Overscript[(Subscript[V, J]^(K+1)), n])-2 (Overscript[(Subscript[V, J]^K), n]))+α((Overscript[(Subscript[U, J+1]^K), n])-2 (Overscript[(Subscript[U, J]^K), n])))     =>
(Overscript[(Subscript[Ψ, J]^K), n+1])=α (Subscript[H, J]^K) (Overscript[(Subscript[U, J+1]^K), n])+(Overscript[(Subscript[U, J]^K), n])((((Overscript[(Subscript[U, J]^K), n+1])+(Overscript[(Subscript[U, J]^K), n])(Δt-1))/g)-α (Subscript[H, J]^K) -Δt)+β (Subscript[H, J]^K)(Overscript[(Subscript[V, J]^(K+1)), n])+(Overscript[(Subscript[V, J]^K), n])((((Overscript[(Subscript[V, J]^K), n+1])+(Overscript[(Subscript[V, J]^K), n])(Δt-1))/g)-β (Subscript[H, J]^K))+(Overscript[(Subscript[Ψ, J]^K), n])(1+β((Overscript[(Subscript[V, J]^(K+1)), n])- (Overscript[(Subscript[V, J]^K), n]))+α((Overscript[(Subscript[U, J+1]^K), n])- (Overscript[(Subscript[U, J]^K), n])))     =>
(Overscript[(Subscript[Ψ, J]^K), n+1])=α (Subscript[H, J]^K) (Overscript[(Subscript[U, J+1]^K), n])-α (Subscript[H, J]^K)(Overscript[(Subscript[U, J]^K), n]) +(Overscript[(Subscript[U, J]^K), n])((((Overscript[(Subscript[U, J]^K), n+1])+(Overscript[(Subscript[U, J]^K), n])(Δt-1))/g)-Δt)+β (Subscript[H, J]^K)(Overscript[(Subscript[V, J]^(K+1)), n])+(Overscript[(Subscript[V, J]^K), n])((((Overscript[(Subscript[V, J]^K), n+1])+(Overscript[(Subscript[V, J]^K), n])(Δt-1))/g)-β (Subscript[H, J]^K))+(Overscript[(Subscript[Ψ, J]^K), n])(1+β((Overscript[(Subscript[V, J]^(K+1)), n])- (Overscript[(Subscript[V, J]^K), n]))+α((Overscript[(Subscript[U, J+1]^K), n])- (Overscript[(Subscript[U, J]^K), n])))     =>
(Overscript[(Subscript[Ψ, J]^K), n+1])=α (Subscript[H, J]^K) ((Overscript[(Subscript[U, J+1]^K), n])-(Overscript[(Subscript[U, J]^K), n])) +(Overscript[(Subscript[U, J]^K), n])((((Overscript[(Subscript[U, J]^K), n+1])+(Overscript[(Subscript[U, J]^K), n])(Δt-1))/g)-Δt)+β (Subscript[H, J]^K)(Overscript[(Subscript[V, J]^(K+1)), n])+(Overscript[(Subscript[V, J]^K), n])((((Overscript[(Subscript[V, J]^K), n+1])+(Overscript[(Subscript[V, J]^K), n])(Δt-1))/g)-β (Subscript[H, J]^K))+(Overscript[(Subscript[Ψ, J]^K), n])(1+β((Overscript[(Subscript[V, J]^(K+1)), n])- (Overscript[(Subscript[V, J]^K), n]))+α((Overscript[(Subscript[U, J+1]^K), n])- (Overscript[(Subscript[U, J]^K), n])))     =>
(Overscript[(Subscript[Ψ, J]^K), n+1])=α (Subscript[H, J]^K) ((Overscript[(Subscript[U, J+1]^K), n])-(Overscript[(Subscript[U, J]^K), n])) +(Overscript[(Subscript[U, J]^K), n])((((Overscript[(Subscript[U, J]^K), n+1])+(Overscript[(Subscript[U, J]^K), n])(Δt-1))/g)-Δt)+β (Subscript[H, J]^K)((Overscript[(Subscript[V, J]^(K+1)), n])-(Overscript[(Subscript[V, J]^K), n]))+(Overscript[(Subscript[V, J]^K), n])((((Overscript[(Subscript[V, J]^K), n+1])+(Overscript[(Subscript[V, J]^K), n])(Δt-1))/g))+(Overscript[(Subscript[Ψ, J]^K), n])(1+β((Overscript[(Subscript[V, J]^(K+1)), n])- (Overscript[(Subscript[V, J]^K), n]))+α((Overscript[(Subscript[U, J+1]^K), n])- (Overscript[(Subscript[U, J]^K), n])))     =>

If  the velocity in the y direction is 0 with explicit H

(Overscript[(Subscript[Ψ, J]^K), n+1])=α(100-x) (Overscript[(Subscript[U, J+1]^K), n])+(Overscript[(Subscript[U, J]^K), n])((((Overscript[(Subscript[U, J]^K), n+1])+(Overscript[(Subscript[U, J]^K), n])(Δt-1))/g)-α (100-x) -Δt)+β (100-x)(Overscript[(Subscript[V, J]^(K+1)), n])+(Overscript[(Subscript[V, J]^K), n])((((Overscript[(Subscript[V, J]^K), n+1])+(Overscript[(Subscript[V, J]^K), n])(Δt-1))/g)-β (100-x))+(Overscript[(Subscript[Ψ, J]^K), n])(1+β((Overscript[(Subscript[V, J]^(K+1)), n])- (Overscript[(Subscript[V, J]^K), n]))+α((Overscript[(Subscript[U, J+1]^K), n])-(Overscript[(Subscript[U, J]^K), n])))      =>
(Overscript[(Subscript[Ψ, J]^K), n+1])=α(100-x) (Overscript[(Subscript[U, J+1]^K), n])+(Overscript[(Subscript[U, J]^K), n])((((Overscript[(Subscript[U, J]^K), n+1])+(Overscript[(Subscript[U, J]^K), n])(Δt-1))/g)-α (100-x) -Δt)+β (100-x)×0+0×(((0+0×(Δt-1))/g)-β (100-x))+(Overscript[(Subscript[Ψ, J]^K), n])(1+β(0-0)+α((Overscript[(Subscript[U, J+1]^K), n])-(Overscript[(Subscript[U, J]^K), n])))      =>
(Overscript[(Subscript[Ψ, J]^K), n+1])=α(100-x) (Overscript[(Subscript[U, J+1]^K), n])+(Overscript[(Subscript[U, J]^K), n])((((Overscript[(Subscript[U, J]^K), n+1])+(Overscript[(Subscript[U, J]^K), n])(Δt-1))/g))-α (Overscript[(Subscript[U, J]^K), n]) (100-x) -Δt (Overscript[(Subscript[U, J]^K), n])+(Overscript[(Subscript[Ψ, J]^K), n])(1+α((Overscript[(Subscript[U, J+1]^K), n])-(Overscript[(Subscript[U, J]^K), n])))      =>
(Overscript[(Subscript[Ψ, J]^K), n+1])=α(100-x) ((Overscript[(Subscript[U, J+1]^K), n])-(Overscript[(Subscript[U, J]^K), n]))+(Overscript[(Subscript[U, J]^K), n])((((Overscript[(Subscript[U, J]^K), n+1])+(Overscript[(Subscript[U, J]^K), n])(Δt-1))/g)) -Δt (Overscript[(Subscript[U, J]^K), n])+(Overscript[(Subscript[Ψ, J]^K), n])(1+α((Overscript[(Subscript[U, J+1]^K), n])-(Overscript[(Subscript[U, J]^K), n])))      =>

Unaltered Final Equation

We will now consider the final equation given:
((∂ψ)/(∂t))=-((∂(u (H+ψ)))/(∂x))-((∂(v (H+ψ)))/(∂y))       =>
((∂ψ)/(∂t))=-(((∂u)/(∂x))(H+ψ))-(((∂(H+ψ))/(∂y))v)      =>
((∂ψ)/(∂t))=-(((∂u)/(∂x))(H+ψ)+(((∂H)/(∂x))+((∂ψ)/(∂x)))u)-(((∂v)/(∂y))(H+ψ)+(((∂H)/(∂y))+((∂ψ)/(∂y)))v)   =>
(((Subscript[Ψ, J,K]^(n+1))-(Subscript[Ψ, J,K]^n))/(Δt))=-((((Subscript[U, J+1,K]^n)-(Subscript[U, J,K]^n))/(Δx))((Subscript[H, J,K])-(Subscript[Ψ, J,K]^n))+((((Subscript[H, J+1,K])-(Subscript[H, J,K]))/(Δx))+(((Subscript[Ψ, J+1,K]^n)-(Subscript[Ψ, J,K]^n))/(Δx)))(Subscript[U, J,K]^n))-((((Subscript[V, J,K+1]^n)-(Subscript[V, J,K]^n))/(Δy))((Subscript[H, J,K])-(Subscript[Ψ, J,K]^n))+((((Subscript[H, J,K+1])-(Subscript[H, J,K]))/(Δy))+(((Subscript[Ψ, J,K+1]^n)-(Subscript[Ψ, J,K]^n))/(Δy)))(Subscript[V, J,K]^n))   =>
(Subscript[Ψ, J,K]^(n+1))-(Subscript[Ψ, J,K]^n)=-Δt ((((Subscript[U, J+1,K]^n)-(Subscript[U, J,K]^n))/(Δx))((Subscript[H, J,K])-(Subscript[Ψ, J,K]^n))+((((Subscript[H, J+1,K])-(Subscript[H, J,K]))/(Δx))+(((Subscript[Ψ, J+1,K]^n)-(Subscript[Ψ, J,K]^n))/(Δx)))(Subscript[U, J,K]^n))-Δt ((((Subscript[V, J,K+1]^n)-(Subscript[V, J,K]^n))/(Δy))((Subscript[H, J,K])-(Subscript[Ψ, J,K]^n))+((((Subscript[H, J,K+1])-(Subscript[H, J,K]))/(Δy))+(((Subscript[Ψ, J,K+1]^n)-(Subscript[Ψ, J,K]^n))/(Δy)))(Subscript[V, J,K]^n))   =>
(Subscript[Ψ, J,K]^(n+1))=-Δt ((((Subscript[U, J+1,K]^n)-(Subscript[U, J,K]^n))/(Δx))((Subscript[H, J,K])-(Subscript[Ψ, J,K]^n))+((((Subscript[H, J+1,K])-(Subscript[H, J,K]))/(Δx))+(((Subscript[Ψ, J+1,K]^n)-(Subscript[Ψ, J,K]^n))/(Δx)))(Subscript[U, J,K]^n))-Δt ((((Subscript[V, J,K+1]^n)-(Subscript[V, J,K]^n))/(Δy))((Subscript[H, J,K])-(Subscript[Ψ, J,K]^n))+((((Subscript[H, J,K+1])-(Subscript[H, J,K]))/(Δy))+(((Subscript[Ψ, J,K+1]^n)-(Subscript[Ψ, J,K]^n))/(Δy)))(Subscript[V, J,K]^n))+(Subscript[Ψ, J,K]^n) =>
(Subscript[Ψ, J,K]^(n+1))=α((Subscript[U, J+1,K]^n)-(Subscript[U, J,K]^n))((Subscript[H, J,K])-(Subscript[Ψ, J,K]^n))+(Subscript[U, J,K]^n)(α((Subscript[H, J+1,K])-(Subscript[H, J,K]))+α((Subscript[Ψ, J+1,K]^n)-(Subscript[Ψ, J,K]^n)))+ β((Subscript[V, J,K+1]^n)-(Subscript[V, J,K]^n))((Subscript[H, J,K])-(Subscript[Ψ, J,K]^n))+(Subscript[V, J,K]^n)(β((Subscript[H, J,K+1])-(Subscript[H, J,K]))+β((Subscript[Ψ, J,K+1]^n)-(Subscript[Ψ, J,K]^n)))+(Subscript[Ψ, J,K]^n) .

Implementations

For[k=1,k≤(Nx+1)*(Ny+1),k++,
M[[k]]=RotateRight[M[[k]],1]
]

Clear[a,b,c,d,Ny,Nx,Δx,Δy,g,r,h,t,f,M,A,B,x,y,α,n,e,Δt,gra,U1,U,Ψ1,Ψ2,H,H1,V,V1,K,q,SM]
{a,b}={0.,1.}; (*length to the beach:  1 mile*)
{c,d}={0.,1000/5280.}; (*width of the beach:  1000 feet*)
{e,f}={0.,123.};(*time: 1 hour*)
(*NOTE:  Nx*Ny must be less than Nt*)
Ny=3.; (*number of steps in the y direction*)
Nx=12.; (*number of steps in the x direction*)
Nt=40.; (*number of steps in the t direction*)
Δx=((b-a)/Nx);
Δy=((d-c)/Ny);
Δt=((e-f)/Nt);
α=-(( Δt)/(Δx));
β=-(( Δt)/(Δy));
gra=1.;
b=1.;
g[x_,y_,t_]=30/5280. (E^(-100.(((x-t))^2)));
φ[x_,y_]=(100./5280.)-(100./5280.)x;
(*(Subscript[V, J,K]^(n+1))=β g((Subscript[Ψ, J,K+1]^n)-(Subscript[Ψ, J,K]^n))-(Subscript[V, J,K]^n)(Δt b-1)*)
(*(Subscript[Ψ, J,K]^(n+1))=α((Subscript[U, J+1,K]^n)-(Subscript[U, J,K]^n))((Subscript[H, J,K])-(Subscript[Ψ, J,K]^n))+(Subscript[U, J,K]^n)(α((Subscript[H, J+1,K])-(Subscript[H, J,K]))+α((Subscript[Ψ, J+1,K]^n)-(Subscript[Ψ, J,K]^n)))+ β((Subscript[V, J,K+1]^n)-(Subscript[V, J,K]^n))((Subscript[H, J,K])-(Subscript[Ψ, J,K]^n))+(Subscript[V, J,K]^n)(β((Subscript[H, J,K+1])-(Subscript[H, J,K]))+β((Subscript[Ψ, J,K+1]^n)-(Subscript[Ψ, J,K]^n)))+(Subscript[Ψ, J,K]^n) *)
(*y derivative is 0:
(Subscript[Ψ, J,K]^(n+1))=α((Subscript[U, J+1,K]^n)-(Subscript[U, J,K]^n))((Subscript[H, J,K])-(Subscript[Ψ, J,K]^n))+(Subscript[U, J,K]^n)(α((Subscript[H, J+1,K])-(Subscript[H, J,K]))+α((Subscript[Ψ, J+1,K]^n)-(Subscript[Ψ, J,K]^n)))+(Subscript[Ψ, J,K]^n) *)
(*(Subscript[U, J,K]^(n+1))=α g((Subscript[Ψ, J+1,K]^n)-(Subscript[Ψ, J,K]^n))-(Subscript[U, J,K]^n)(Δt b-1) *)
SM=Table[g[x,y,0],{y,c,d,Δy},{x,a,b,Δx}];
MatrixForm[SM];
(*initial conditions on the height of the waves at time 0*)
M=SparseArray[
{{i_,j_,q_}/;i==1->g[(q/(((Nx+1)*(Ny+1)))),(j/(((Nx+1)*(Ny+1)))),((i-1)/(Nt+1))]},{(Nt+1),(Nx+1)*(Ny+1),(Nx+1)*(Ny+1)}];
(*The (Subscript[Ψ, J+1,K]^n) sparse array matrix*)
Ψ1=SparseArray[{{i_,i_,i_}->0},{Nt+1,(Nx+1)*(Ny+1),(Nx+1)*(Ny+1)}];
(*The (Subscript[Ψ, J,K+1]^n) sparse array matrix*)
Ψ2=SparseArray[{{i_,i_,i_}->0},{Nt+1,(Nx+1)*(Ny+1),(Nx+1)*(Ny+1)}];
U=SparseArray[{{i_,j_,q_}/;i==1->((q/((Nx+1)*(Ny+1))))},{Nt+1,(Nx+1)*(Ny+1),(Nx+1)*(Ny+1)}];
U1=SparseArray[{{i_,i_,i_}->0},{Nt+1,(Nx+1)*(Ny+1),(Nx+1)*(Ny+1)}];
H=SparseArray[{{i_,j_,q_}->φ[(q/(((Nx+1)*(Ny+1)))),(j/(((Nx+1)*(Ny+1))))]},{Nt+1,(Nx+1)*(Ny+1),(Nx+1)*(Ny+1)}];
H1=SparseArray[{{i_,i_,i_}->0},{Nt+1,(Nx+1)*(Ny+1),(Nx+1)*(Ny+1)}];
(*H1=SparseArray[{{i_,i_,i_}->0},{(Nx+1)*(Ny+1),(Nx+1)*(Ny+1),(Nt+1)}];*)
(*M[[1]]=Table[g[x,y,0],{y,c,d,Δy},{x,a,b,Δx}]*)
(*U[[1]]=Table[1,{x,a,b,Δx},{y,c,d,Δy}];*)
(*guide of placements in 3D matrix: K[[t,y,x]]*)
Do[Ψ2[[n-1]]=RotateLeft[M[[n-1]],1];
For[k=1,k≤(Nx+1)*(Ny+1),k++,
Ψ1[[n-1,k]]=RotateLeft[Ψ1[[n-1,k]],1];
Ψ1[[n-1,k,(Nx+1)*(Ny+1)]]=0];
For[r=1,r≤(Nx+1)*(Ny+1),r++,
U1[[n-1,r]]=RotateLeft[U[[n-1,r]],1];
U1[[n-1,r,(Nx+1)*(Ny+1)]]=0];
For[h=1,h≤(Nx+1)*(Ny+1),h++,
H1[[n-1,h]]=RotateLeft[H[[n-1,h]],1];
H1[[n-1,h,(Nx+1)*(Ny+1)]]=0];
U[[n]]=α*gra*(Ψ1[[n-1]]-M[[n-1]])-U[[n-1]]*(Δt*b-1);
M[[n]]=α*(U1[[n-1]]-U[[n-1]])*(H[[n-1]]-M[[n-1]])+U[[n-1]]*(α(H1[[n-1]]-H[[n-1]])+α*(Ψ1[[n-1]]-M[[n-1]]))+M[[n-1]],
{n,2,Nt+1}]
MatrixForm[α*gra*(Ψ1[[2]]-M[[2]])-U[[2]]*(Δt*b-1)];
MatrixForm[M[[3]]];
Manipulate[ListPlot3D[M[[κ]],
PlotRange->{((-1)/5280),(60/5280)}
],{κ,1,Nt+1,1}]
(*Animate[ListPlot3D[M[[ρ]],
PlotRange->{0,3}],{ρ,1,Nt+2,1}]*)
Plot3D[g[x,y,0],{x,0,1},{y,0,1},
PlotRange->{0,(60/5280)}]

Max[M[[-1]]]
Min[M[[-1]]]

ListPlot3D[M[[-1]],
PlotRange->All
]

Clear[a,b,c,d,Ny,Nx,Δx,Δy,g,r,h,t,f,M,A,B,x,y,α,n,e,Δt,gra,U1,U,Ψ1,Ψ2,H,H1,V,V1,K,q,SM]
{a,b}={0.,1.}; (*length to the beach:  1 mile*)
{c,d}={0.,1000/5280.}; (*width of the beach:  1000 feet*)
{e,f}={0.,1.};(*time: 1 hour*)
Ny=12; (*number of steps in the y direction*)
Nx=13; (*number of steps in the y direction*)
Nt=6; (*number of steps in the y direction*)
Δx=((b-a)/Nx);
Δy=((d-c)/Ny);
Δt=((e-f)/Nt);
α=-(( Δt)/(Δx));
β=-(( Δt)/(Δy));
gra=9.8;
b=1.;
g[x_,y_,t_]=(E^(-100.(((x-t))^2)));
K=SparseArray[{{i_,j_,q_}/;i==1->g[(((b-a)*j)/((Nx+1)*(Ny+1))),(((d-c)*q)/((Nx+1)*(Ny+1))),0]},{(Nt+1),(Nx+1)*(Ny+1),(Nx+1)*(Ny+1)}]

T=SparseArray[{{i_,i_,i_}->0},{Nt+1,(Nx+1)*(Ny+1),(Nx+1)*(Ny+1)}];
Manipulate[ListPlot3D[K[[κ]],
PlotRange->{0,3}],{κ,1,Nt+2,1}]
Plot3D[g[x,y,0],{x,0,1},{y,0,1},
PlotRange->{0,3}]

Plot3D[(E^(-100(((x-.2))^2))),{x,0,1},{y,0,100/5280},
PlotRange->{0,3}]

Things

(Subscript[U, i,j]^(n+1))=α g((Subscript[Ψ, i,j+1]^n)-(Subscript[Ψ, i,j]^n))-(Subscript[U, i,j]^n)(Δt b-1)

U0=RandomReal[ {-3,3},{4,3}];
V0= RandomReal[ {-3,3},{4,3}];

U=U0;
V=V0;
Do[
W=RotateLeft[V]-V;
U=U -Sin[ U *W];
V=V+W,
{123}]
U

Clear[U,V]
U = Table[(Subscript[u, i,j]),{i,4},{j,5}];
MatrixForm[
RotateLeft[U]-U
]

W= Transpose[U];
MatrixForm[
Transpose[RotateLeft[W]-W]
]

DiffX[MatIn_] :=Module[
{i,j},
{n,m}=Dimensions[MatIn];
Table[
MatIn[i+1,j]-Mat[i,j],
{i,n},
{j,m}
]
]