program JACOBIAN c c this program computes the arakawa jacobian over a domain of grid c -ed data.It uses a nine point fouth order scheme. c parameter (l=10,m=20,l1=l-1,m1=m-1,l2=l-2,m2=m-2) real psi(l,m),zta(l,m),a(l,m),x(l),y(m),dx(m) pi = 4.*atan(1.0) h = 200. dy = h do 2300 j = 1, m 2300 dx(j) = dy yk = 2.*pi/1000. yl = pi / 1000. x(1) = 0. y(1) = 0. do 2302 i = 2, l im1 = i-1 do 2302 j = 2, m jm1 = j-1 x(i) = x(im1) + h y(j) = y(jm1) + h 2302 continue sum = 0. do 2304 i = 1, l do 2304 j = 1, m c c define the analytical functions psi and zta. c psi(i,j) = sin(yk*x(i)) * sin(yl*y(j)) + cos(yl*y(j)) zta(i,j) = -(yk**2+yl**2)*sin(yk*x(i))*sin(yl*y(j))- & yl**2 * cos(yl*y(j)) a(1,j) = zta(1,j) a(l,j) = zta(l,j) a(i,1) = zta(i,1) a(i,m) = zta(i,m) sum = sum + (zta(i,j) / (l*m))**2 2304 continue c c compute the jacobian c call JAC (a,psi,zta,dx,dy,l,m,l1,m1,l2,m2) c c display outputs for 1 column . c write (6,1000) write (6,1001) write (6,1002)((i,j,a(i,j),i=1,l),j=4,4) 1000 format(//,20x,'arakawa jacobian scheme.',//) 1001 format(2x,'i j',10x, 'estimated jacobian',//) 1002 format( (2i3,8x,e15.8) ) stop end subroutine JAC (a,b,c,dx,dy,l,m,l1,m1,l2,m2) c c This subroutine performs the so-called Arakawa c Jacobian .This scheme satisfies the integral c relations of invariance for the total kinetic c energy and the mean square vorticity. c The following code uses a nine-point c second order Arakawa Jacobian. c c variable definitions c c input : b and c c output : a c dx : east - west grid distance. c dy : north-south grid distance. c l,m : first and second dimension of a,b,c c real a(l,m,1), b(l,m,1), dx(m), c(l,m,1) k = 1 do 2310 j = 2, m1 dm = 12.*dy*dx(j) do 2310 i = 1, l if (i-1) 80,80,81 80 im1 = l1 ip1 = 2 go to 83 81 if (i-l) 82,80,80 82 im1 = i-1 ip1 = i+1 83 continue c a(i,j,k) = (b(i,j-1,k)+b(ip1,j-1,k)-b(i,j+1,k)-b(ip1,j+1,k)) & *(c(ip1,j,k)-c(i,j,k))+(b(im1,j-1,k)+b(i,j-1,k)-b(im1,j+1,k & )-b(i,j+1,k))*(c(i,j,k)-c(im1,j,k))+(b(ip1,j,k)+b(ip1,j+1,k & )-b(im1,j,k)-b(im1,j+1,k))*(c(i,j+1,k)-c(i,j,k))+(b(ip1,j-1 & ,k)+b(ip1,j,k)-b(im1,j-1,k)-b(im1,j,k))*(c(i,j,k)-c(i,j-1,k & ))+(b(ip1,j,k)-b(i,j+1,k))*(c(ip1,j+1,k)-c(i,j,k))+(b(i,j-1 & ,k)-b(im1,j,k))*(c(i,j,k)-c(im1,j-1,k))+(b(i,j+1,k)-b(im1,j & ,k))*(c(im1,j+1,k)-c(i,j,k))+(b(ip1,j,k)-b(i,j-1,k))*(c(i,j & ,k)-c(ip1,j-1,k)) c a(i,j,k) = a(i,j,k)/dm 2310 continue do 2311 i = 1, l if (i-1) 70,70,71 70 im1 = l1 ip1 = 2 go to 73 71 if (i-l) 72,70,70 72 im1 = i-1 ip1 = i+1 73 continue dm = 12.*dy*dx(1) c a(i,1,k) = (b(i,1,k)+b(ip1,1,k)-b(i,2,k)-b(ip1,2,k))*(c(i,1,k)+ & c(ip1,1,k))-(b(im1,1,k)+b(i,1,k)-b(im1,2,k)-b(i,2,k))*(c(im1 & ,1,k)+c(i,1,k))+(b(ip1,1,k)+b(ip1,2,k)-b(im1,1,k)-b(im1,2,k) & )*(c(i,1,k)+c(i,2,k))+(b(ip1,1,k)-b(i,2,k))*(c(i,1,k)+c(ip1, & 2,k))+(b(i,2,k)-b(im1,1,k))*(c(im1,2,k)+c(i,1,k)) c a(i,1,k) = a(i,1,k)/dm dm = 12.*dy*dx(m) c a(i,m,k) = (b(i,m-1,k)+b(ip1,m-1,k)-b(i,m,k)-b(ip1,m,k))*(c(i,m, & k)+c(ip1,m,k))-(b(im1,m-1,k)+b(i,m-1,k)-b(im1,m,k)-b(i,m,k))* & (c(im1,m,k)+c(i,m,k))-(b(ip1,m-1,k)+b(ip1,m,k)-b(im1,m-1,k)-b & (im1,m,k))*(c(i,m-1,k)+c(i,m,k))-(b(i,m-1,k)-b(im1,m,k))*(c(i & m1,m-1,k)+c(i,m,k))-(b(ip1,m,k)-b(i,m-1,k))*(c(i,m,k)+c(ip1,m & -1,k)) c a(i,m,k) = a(i,m,k)/dm 2311 continue return end