c This code contains several subroutines for
c manipulating arrays in spherical coordinates. Most of the
c routines require the variables in common blocks /cons/
c and /ncepll/.
c
c Included routines:
c ** ddx
c ** ddy
c ** del2
c ** pson
c
subroutine ddx(h,scale,dhdx)
c This routine calculates the x-derivative of h, multiplied by
c scale to give dhdx, where x = a*cos(lat)*lon
c
parameter (nx=61,ny=61)
dimension h(nx,ny), dhdx(nx,ny)
c
dimension rlond(nx),rlatd(ny)
dimension rlonr(nx),rlatr(ny)
dimension sinlat(ny),coslat(ny),tanlat(ny)
c
common /cons/ pi,g,rd,dtr,erad,erot
common /ncepll/ rlond,rlonr,rlatd,rlatr,dlon,dlat,
+ dlonr,dlatr,sinlat,coslat,tanlat
c
c Interior points
cf = scale/(erad*2.0*dlonr)
c
do 10 j=1,ny
do 10 i=2,nx-1
dhdx(i,j) = cf*(h(i+1,j)-h(i-1,j))/coslat(j)
10 continue
c
c x-boundary points
cf = 2.0*cf
do 20 j=1,ny
dhdx( 1,j) = cf*(h( 2,j) - h( 1,j))/coslat(j)
dhdx(nx,j) = cf*(h(nx,j) - h(nx-1,j))/coslat(j)
20 continue
c
return
end
subroutine ddy(h,scale,dhdy)
c This routine calculates the y-derivative of h, multiplied by
c scale to give dhdy, where y = a*lat
c
parameter (nx=61,ny=61)
dimension h(nx,ny), dhdy(nx,ny)
c
dimension rlond(nx),rlatd(ny)
dimension rlonr(nx),rlatr(ny)
dimension sinlat(ny),coslat(ny),tanlat(ny)
c
common /cons/ pi,g,rd,dtr,erad,erot
common /ncepll/ rlond,rlonr,rlatd,rlatr,dlon,dlat,
+ dlonr,dlatr,sinlat,coslat,tanlat
c
c Interior points
cf = scale/(erad*2.0*dlatr)
c
do 10 i=1,nx
do 10 j=2,ny-1
dhdy(i,j) = cf*(h(i,j+1)-h(i,j-1))
10 continue
c
c y-boundary points
cf = 2.0*cf
do 20 i=1,nx
dhdy(i, 1) = cf*(h(i, 2) - h(i, 1))
dhdy(i,ny) = cf*(h(i,ny) - h(i,ny-1))
20 continue
c
return
end
subroutine del2(h,scale,del2h)
c This routine calculates the Laplacian of h, multiplied by
c scale, to give del2h
c
parameter (nx=61,ny=61)
dimension h(nx,ny), del2h(nx,ny)
c
dimension rlond(nx),rlatd(ny)
dimension rlonr(nx),rlatr(ny)
dimension sinlat(ny),coslat(ny),tanlat(ny)
c
common /cons/ pi,g,rd,dtr,erad,erot
common /ncepll/ rlond,rlonr,rlatd,rlatr,dlon,dlat,
+ dlonr,dlatr,sinlat,coslat,tanlat
c
cfx2 = scale/(erad*erad*dlonr*dlonr)
cfy2 = scale/(erad*erad*dlatr*dlatr)
cfy1 = -scale/(erad*erad*2.0*dlatr)
c
c Interior points
do 10 j=2,ny-1
cli2 = 1.0/(coslat(j)**2)
do 10 i=2,nx-1
del2h(i,j) = cfx2*(h(i+1,j)+h(i-1,j) - 2.0*h(i,j))/cli2 +
+ cfy2*(h(i,j+1)+h(i,j-1) - 2.0*h(i,j)) +
+ cfy1*(h(i,j+1)-h(i,j-1))*tanlat(j)
10 continue
c
c
c Bottom and top boundary points
do 15 i=2,nx-1
j = 1
cli2 = 1.0/(coslat(j)**2)
del2h(i,j) = cfx2*(h(i+1,j)+h(i-1,j) - 2.0*h(i,j ))/cli2 +
+ cfy2*(h(i,j+2)+h(i ,j) - 2.0*h(i,j+1)) +
+ cfy1*(h(i,j+1)-h(i,j))*2.0*tanlat(j)
c
j = ny
cli2 = 1.0/(coslat(j)**2)
del2h(i,j) = cfx2*(h(i+1,j)+h(i-1,j ) - 2.0*h(i,j ))/cli2 +
+ cfy2*(h(i ,j)+h(i ,j-2) - 2.0*h(i,j-1)) +
+ cfy1*(h(i,j)-h(i,j-1))*2.0*tanlat(j)
15 continue
c
c left and right boundary points
do 20 j=2,ny-1
i = 1
cli2 = 1.0/(coslat(j)**2)
del2h(i,j) = cfx2*(h(i+2,j )+h(i,j ) - 2.0*h(i+1,j))/cli2 +
+ cfy2*(h(i ,j+1)+h(i,j-1) - 2.0*h(i,j)) +
+ cfy1*(h(i ,j+1)-h(i,j-1))*tanlat(j)
c
i = nx
cli2 = 1.0/(coslat(j)**2)
del2h(i,j) = cfx2*(h(i,j )+h(i-2,j ) - 2.0*h(i-1,j))/cli2 +
+ cfy2*(h(i,j+1)+h(i ,j-1) - 2.0*h(i ,j)) +
+ cfy1*(h(i,j+1)-h(i ,j-1))*tanlat(j)
c
20 continue
c
c Upper-left corner
i = 1
j = ny
cli2 = 1.0/(coslat(j)**2)
del2h(i,j) = cfx2*(h(i+2,j)+h(i,j ) - 2.0*h(i+1,j ))/cli2 +
+ cfy2*(h(i ,j)+h(i,j-2) - 2.0*h(i ,j-1)) +
+ cfy1*(h(i ,j)-h(i,j-1))*tanlat(j)
c
c Lower-left corner
i = 1
j = 1
1 cli2 = 1.0/(coslat(j)**2)
del2h(i,j) = cfx2*(h(i+2,j )+h(i,j) - 2.0*h(i+1,j ))/cli2 +
+ cfy2*(h(i ,j+2)+h(i,j) - 2.0*h(i ,j+1)) +
+ cfy1*(h(i ,j+1)-h(i,j))*tanlat(j)
c
c Upper-right corner
i = nx
j = ny
cli2 = 1.0/(coslat(j)**2)
del2h(i,j) = cfx2*(h(i-2,j )+h(i,j) - 2.0*h(i-1,j ))/cli2 +
+ cfy2*(h(i ,j)+h(i,j-2) - 2.0*h(i ,j-1)) +
+ cfy1*(h(i ,j)-h(i,j-1))*tanlat(j)
c
c Lower-right corner
i = nx
j = 1
cli2 = 1.0/(coslat(j)**2)
del2h(i,j) = cfx2*(h(i-2,j )+h(i,j) - 2.0*h(i-1,j ))/cli2 +
+ cfy2*(h(i ,j+2)+h(i,j) - 2.0*h(i ,j+1)) +
+ cfy1*(h(i ,j+1)-h(i,j))*tanlat(j)
c
return
end
subroutine pson(f,hfg,h,ierr)
c This routine solves (del**2)h = f in spherical geometry using
c over-relaxation. The boundary values for h and the first guess
c for h at the interior points are assumed to be in the array hfg.
c If the relaxation converges, ierr=0, otherwise, ierr=1.
c
parameter (nx=61,ny=61)
dimension f(nx,ny),hfg(nx,ny),h(nx,ny)
c
dimension rlond(nx),rlatd(ny)
dimension rlonr(nx),rlatr(ny)
dimension sinlat(ny),coslat(ny),tanlat(ny)
c
common /cons/ pi,g,rd,dtr,erad,erot
common /log/ lulog
common /ncepll/ rlond,rlonr,rlatd,rlatr,dlon,dlat,
+ dlonr,dlatr,sinlat,coslat,tanlat
c
c Local variables for storing coefficients
dimension c1(ny),c2(ny),c3(ny),c4(ny)
c
c Specify max number of iterations and error check increment
nit = 200
iecal = 10
emax = 2.0e-5
c
c Calculate constants for over-relaxation
def = 1.0 - 2.0*((sin(pi/(2.0*float(ny-1))))**2)
omega = 2.0/(1.0 + sqrt(1.0-def*def))
oomega = 1.0-omega
c
c Move first guess h to h array
do 10 j=1,ny
do 10 i=1,nx
h(i,j) = hfg(i,j)
10 continue
c
c Find the maximum magnitude of f for normalizing the error
enorm = 0.0
do 15 j=2,ny-1
do 15 i=2,nx-1
if (abs(f(i,j)) .gt. enorm) enorm = abs(f(i,j))
15 continue
c
if (enorm .le. 0.0) then
c The forcing term is zero. Scale error from boundary values.
bsum = 0.0
count = 0.0
do 16 j=1,ny
do 16 i=1,nx
if (j .ne. 1 .and. j .ne. ny .and.
+ i .ne. 1 .and. i .ne. nx ) go to 16
count = count + 1.0
bsum = bsum + abs(h(i,j))
16 continue
bsum = bsum/count
c
enorm = bsum/(dlatr*dlatr*erad*erad)
endif
c
c Calculate common factors for iteration
do 20 j=1,ny
sdlonr = coslat(j)*dlonr
t1 = ( (dlatr*sdlonr)**2 )/
+ ( dlatr**2 + sdlonr**2 )
c
c1(j) = 0.5*t1/(sdlonr*sdlonr)
c2(j) = 0.5*t1/(dlatr*dlatr)
c3(j) = -0.25*t1/(dlatr)
c4(j) = -0.5*t1*erad*erad
20 continue
c
c Perform iteration
errtp = 1.0e+10
c
do 25 k=1,nit
do 30 j=2,ny-1
do 30 i=2,nx-1
h(i,j) = oomega*h(i,j) + omega*(
+ c1(j)*(h(i+1,j)+h(i-1,j)) +
+ c2(j)*(h(i,j+1)+h(i,j-1)) +
+ c3(j)*(h(i,j+1)-h(i,j-1)) +
+ c4(j)*(f(i,j)) )
30 continue
c
c Check for convergence
if (mod(k,iecal) .eq. 0) then
call rchk(f,h,enorm,errt)
c
write(lulog,888) k,errt,emax,enrom
888 format(1x,'k=',i3,' errt,emax,enorm= ',3(e11.4))
c
if (errt .le. emax .or. errt .gt. errtp) go to 1000
errtp = errt
endif
c
25 continue
c
ierr = 1
return
c
1000 continue
ierr = 0
c
return
end
subroutine rchk(f,h,enorm,errt)
c This routine calculates the relative error between (del**2)h and f for
c determining if the relaxation in routine pson has converged.
c
parameter (nx=61,ny=61)
dimension f(nx,ny),h(nx,ny)
c
c Local array
dimension d2h(nx,ny)
c
rpts = float( (nx-1)*(ny-1) )
errt = 0.0
scale = 1.0
c
call del2(h,scale,d2h)
do 10 j=2,ny
do 10 i=2,nx
errt = errt + (f(i,j)-d2h(i,j))**2
10 continue
c
errt = ( sqrt(errt/rpts) )/enorm
c
return
end