c
c Visualizations of linear QG vorticity waves
c
c Apr.-3, 2024 Andreas Muenchow, University of Delaware
c		compiles with gfortran (GNU fortran 13.2.0)
c		on Apple's macos 12.7.4 (Monterey)
c
	real arg,x,y,z
	real m,dm,h,N,zero
	real m0,m1,m2,arg0,arg1,arg2
	real k,l,mz0,mz1,mz2,Ld
c
c Parameters in SI units
c
	H = 1000.
        N = 0.01
	g = 9.81
	fact = N*N*H/g
	write(6,*)'N^2*H/g: ',fact
c
c 43.45 N latitude (17.45 hours inertial period)
c
	f0 = 0.0001
	beta = 1e-11
	pi = atan2(1.,1.)*4.
c	write(6,*)'Inertial Period (17.449 hrs): ',2.*pi/f0/3600.
	zero = 0.
c
c Vertical Structure Functions from Boundary Conditions
c
	m = 0.
	dm = pi/100./2.
	do 1 i=1,600*2
          m = m+dm/H
	  arg = m*H
	  x = cos(arg)
	  y = sin(arg)
          if (x*x .gt. 0.000001) then
             z = y/x
	     y = fact/arg
	     write(16,101)arg/pi,z,y/fact
          end if
1	continue
101	format(f9.3,2f12.5)
c
c Stream Function Amplitude a(z) for modes 0,1,2:
c	determined from solving transzendental equation
c
	m0 = 0.275/H * pi
	m1 = 1.090/H * pi
	m2 = 2.050/H * pi
	z = 0
	dz = H/1000. 	
	do 2 i=1,1000
	   z = z+dz
	   arg0 = m0*z
	   arg1 = m1*z
	   arg2 = m2*z
	   a0 = cos(arg0)
	   a1 = cos(arg1)
	   a2 = cos(arg2)
	   write(17,102)z,a0,a1,a2,zero
2	continue
102	format(f9.3,3f12.5,f3.0)
c
c Dispersion for l=0
c
	Ld = N*H/f0
	write(6,*)'Rossby Radius N*H/f0 (km): ',Ld/1000.
	l = 0.
	k = 0.
	   con = 24.*3600.
	   yr = con*365.25
	mz0 = m0*(f0/N)
	mz1 = m1*(f0/N)
	mz2 = m2*(f0/N)
	dk = 1./Ld/50.
	do 3 i=1,1000
	   k = k-dk
	   w0 = -beta*k/(k*k+l*l+mz0*mz0)
	   w1 = -beta*k/(k*k+l*l+mz1*mz1)
	   w2 = -beta*k/(k*k+l*l+mz2*mz2)
	   write(18,103)k*Ld,(w0/f0),(w1/f0),(w2/f0)
	   write(20,103)k*Ld,(w0*yr),(w1*yr),(w2*yr)
	   write(19,103)k*Ld,1./(w0*con),1./(w1*con),1./(w2*con)
3	continue
103	format(f10.5,3f15.8)
	stop
	end