c
c Rossby Adjustment Problem (Transients)
c	Gill (1982) chapter 7.3
c
c	compiles with fortran g77 om MacBook Pro OS X 10.5.8
c
c Oct.-4, 2005 Andreas Muenchow, University of Delaware
c Mar.-16, 2012
c
	double precision u(50,400),eta(50,400),v(50,400)
	real eta0,dx,arg,arg1,arg2,x,x1,x2
	real a,c,f,h
	open (unit=7,file='adjust.dat')
c
c All units are SI
c
	eta0 = 1.
	f = 0.0001
	g = 9.81
	H = 100.
	c = sqrt(g*H)
	a = c/f
c
c "it" represents the time stepping
c "ix" represents the space stepping
c
	do 1 it=1,20
c	do 2 ix = 1,4
	do 2 ix = 1,400
	   t = 0.5/f*float(it)*4
	   x = 0.5*a*float(ix-1)*0.1
c
c Keep track of where the "front" is
c
	   if (x.lt.c*t) then
c
c Eq. (7.3.14)
c
	   	arg = f*sqrt(t*t-x*x/c/c)
	   	u(it,ix) = g*eta0/c*BesJ0(arg)
c
c calculate v(x,t) as the integral of u(x,t) as d(v)/dt=-f*u from y-momentum
c calculate eta(x,t) as the integral of u(x,t) as d(eta)/dt=-h*d(u)/dx from continuity
c
		v(it,ix) = 0.
		eta(it,ix) = 0.
		dx = a*0.05
		do 3 itt=1,it*100
	           t = 0.5/f*float(itt)*0.04
		   t2 = t*t
		   c2 = c*c
		   x1 = x-dx
		   x2 = x+dx
		   xx = sqrt(max(x1*x1,x2*x2))
		   if (max(x,xx) .lt. c*t) then
	   	      arg  = f*sqrt(t2-x*x/c2)
	   	      arg1 = f*sqrt(t2-x1*x1/c2)
	   	      arg2 = f*sqrt(t2-x2*x2/c2)
	   	      v(it,ix) = v(it,ix)-BesJ0(arg)
c		      test1 = BesJ0(arg2)
c		      test2 = -BesJ0(arg1)
c	write(6,*)test1,test2,test1+test2
		      eta(it,ix) = eta(it,ix)-eta0*(BesJ0(arg2)-BesJ0(arg1))
		   else
		      eta(it,ix) = 0.
		   end if
3		continue
c                            scaling              time-step   space-step
		v(it,ix) =   g*eta0/c* v(it,ix)   *0.5*0.04
		eta(it,ix) = c*eta0*   eta(it,ix) *0.5*0.04/f /dx/2.
	   else
		u(it,ix) = 0.
		eta(it,ix) = -eta0
		v(it,ix) = 0.
	   end if
2	continue
101	format(f6.2,10f9.3)
102	format(f6.2,10f9.5)
1	continue
c
c For some unknown reason eta(it,ix) becomes a "nan" on occasion
c	for some even "ix" (akm Mar.-15, 2012)
c
	do 4 ix=1,400,2
	   write(7,101)(ix-1)*0.5*0.1,(u(it,ix),it=1,20,2)
	   write(8,102)(ix-1)*0.5*0.1,(eta(it,ix),it=1,20,2)
	   write(9,101)(ix-1)*0.5*0.1,(v(it,ix),it=1,20,2)
4	continue
	write(10,103)(float(it)/4,u(it,4),v(it,4),it=1,50)
103	format(f5.1,2f9.3)
	stop
	end