program GEOPOTENTIAL c c this program computes the geopotential height field c from the streamfunction using 3 different methods : c c i : the geostrophic balance. c ii : the non-linear balance . c iii : the linear balance . c parameter (l=21,m=13,np=7) real dx (m),cor(m),a2(l,m),z500(l,m) real psi(l,m),upsi(l,m),vpsi(l,m) c c set the mean height at 500 mband define some constants. c data zbar,beta,slat,grid /5600.,2.0e-11,-15.,2.5 / c c Note that the streamfunction calculated by program c STREAM isused as input . c open (30,file='psirf.dat',status='old') open (40,file='geop.dat',status='unknown') c c read streamfunction field from tape30. c 1000 format(10e13.6) read(30,1000)psi c c define the meridional grid spacing and c the Coriolis parameter. c pi = 4.0*atan(1.0) rad = pi/180. dy = 111.1 * 1000. * grid do 4200 j = 1, m alat = (slat + (j-1)*grid)*rad dx(j) = dy * cos(alat) cor(j) = 2.*7.29/1e05*sin(alat) 4200 continue c c compute geopotential via the geostrophic balance. c the flags fk1 and fk2 are set before the call c fk1 = 0. fk2 = 0. call ZFIELD ( psi,upsi,vpsi,l,m,dx,dy,cor, & a2, z500, fk1, fk2,zbar,beta ) c c write output . c write (40,1000)z500 c c reinitialise the work arrays to zero. c call ZERO (upsi,vpsi,l,m) call ZERO (a2 ,z500,l,m) c c compute geopotential via the linear balance. c fk1 = 1. fk2 = 0. call ZFIELD ( psi,upsi,vpsi,l,m,dx,dy,cor, & a2, z500, fk1, fk2,zbar,beta ) c c write output . c write (40,1000)z500 c c reinitialise the work arrays to zero. c call ZERO (upsi,vpsi,l,m) call ZERO (a2 ,z500,l,m) c c compute geopotential via the non-linear balance. c fk1 = 1. fk2 = 1. call ZFIELD ( psi,upsi,vpsi,l,m,dx,dy,cor, & a2,z500, fk1, fk2,zbar,beta) c c write output . c write (40,1000)z500 stop end subroutine ZERO (a,b,l,m) c c This subroutine initializes c fields to zero c real a(l,m),b(l,m) do 4220 j = 1, m do 4220 i = 1, l a(i,j) = 0.0 4220 b(i,j) = 0.0 return end subroutine ZFIELD (psi,upsi,vpsi,l,m,dx,dy, & cor,a2,z,fk1,fk2,zbar,beta) c c This subroutine computes the geopotential height c from the streamfunction using different balance models c controlled by the following switches. c c : fk1 = 0 and fk2=0 for geostrophic balance c : fk1 = 1 and fk2=0 for linear balance c : fk1 = 1 and fk2=1 for non linear balance c real z(l,m),psi(l,m),dx(m),zinv(100),z1(100) real vpsi(l,m),upsi(l,m),a2(l,m),cor(m) l1 = l-1 l2 = l-2 m1 = m-1 m2 = m-2 ginv = 1.0/9.81 do 4210 j = 1, m z1(j) = dx(j)**2 4210 zinv(j)= 1./z1(j) zz = dy**2 zzinv = 1./zz sum = 0. do 4211 i = 1, l do 4211 j = 1 ,m sum = sum + psi(i,j) 4211 continue psibar = sum/float(l*m) print *,'psibar=',psibar do 4212 j = 2, m1 do 4212 i = 1, l im1 = i-1 ip1 = i+1 if (im1.lt.1) im1 = l1 if (ip1.gt.l) ip1 = 2 vpsi(i,j) = (psi(ip1,j)-psi(im1,j))/(2.*dx(j)) upsi(i,j) = -(psi(i,j+1)*dx(j)/dx(j+1)- & psi(i,j-1)*dx(j)/dx(j-1))/(2.*dy) 4212 continue call LAP (z ,psi,dx,dy,l,m,l1,m1,l2,m2) call JAC (a2,upsi,vpsi,dx,dy,l,m,l1,m1,l2,m2) do 4213 i = 1, l do 4213 j = 2, m1 betau = upsi(i,j)*beta fvort = cor(j)*z (i,j) a2(i,j)= fvort*ginv-betau*ginv*fk1+2.*a2(i,j)*ginv*fk2 4213 continue do 4214 i = 1, l a2(i,1)= 0. a2(i,m)= 0. 4214 continue 10 continue do 4215 j = 1, m do 4215 i = 1, l1 z(i,j) = zbar+cor(j)*ginv*(psi(i,j) -psibar) 4215 continue do 4216 j = 1, m 4216 z(l,j) = z(1,j) c c solve Poisson equation using relaxation . c call RELAX (z,zzinv,z1,a2 ,zinv,l,l1,m1,m) return end subroutine RELAX (x,zzinv,z,y,zinv,l,l1,m1,m) c c This subroutine solve the Poisson equation using c the overrelaxation method . c c definition of variables : c c x contains first guess and eventually c the output of the field to be relaxed c y forcing function of l by m matrix c z latitudinal increment square c zinv inverse of dx**2 c zzinv inverse of dy**2 c l east-west dimension c m north-south dimension c c definition of constants : c npts number of points to be relaxed c nrel number of points relaxed c alfa relaxation factor c ia maximum number of iterations allowed c eps tolerance error c nsc count of number of scan for convergence c lsc last scan after convergence c real x(l,m),y(l,m),z(m),zinv(m) npts = l*(m-2) nlax = 1 mm = 2 mmm = m1 alfa = .46 ia = 1000 eps = 1. nsc = 0 lsc = -1 15 nrel = 0 do 4110 j = mm, mmm do 4110 i = 1, l im1 = i-1 ip1 = i+1 if (im1.lt.1) im1 = l1 if (ip1.gt.l) ip1 = 2 r = (x(ip1,j)+x(im1,j)-2.*x(i,j))*zinv(j)+ & (x(i,j+1)+x(i,j-1)-2.*x(i,j))*zzinv r = (r-y(i,j))*z(j) if (lsc-nsc) 29,29,30 29 x(i,j) = x(i,j) + alfa*r 30 if (abs(r).le.eps) nrel = nrel+1 4110 continue c c nrel gives the number of points that have converged . c nsc gives the number of scan made over the domain and c if it is less then the maximum number of iteration c allowed for complete convergence it keeps going to c statement number 15. lsc is the final check once c convergence has been achieved it does one more loop c before finally jumping out to statement 300 . c nsc = nsc+1 if (nrel-npts) 13,14,14 14 if (lsc .ge. nsc) go to 300 18 lsc = nsc+1 13 if (nsc.lt.ia) go to 15 201 format(50h progress of relaxation npts,nrel,nsc,ia ) 300 continue write(6,201) 200 format(6x,4i9) write (6,200)npts,nrel,nsc,ia return 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 subroutine LAP (a,b,dx,dy,l,m,l1,m1,l2,m2) c c This subroutine calculates the Laplacian of any function. c the standard second-order finite differencing scheme is c used . c This subroutine is different from those used in chapt2 . c It written in a more genral way for a 3D variable . It c also uses different grid spacing ( dx ,dy). c real a(l,m,1), b(l,m,1), dx(m) k = 1 do 3230 j = 2, m1 do 3230 i = 2, l1 3230 a(i,j,k) = (((b(i+1,j,k)+b(i-1,j,k)-2.*b(i,j,k)) & /dx(j)**2+(b(i,j+1,k)+b(i,j-1,k)-2.*b(i,j,k))/ & dy**2)) do 3231 j = 2, m1 a(1,j,k) = (((b(2,j,k)+b(l1,j,k)-2.*b(1,j,k)) & /dx(j)**2+(b(1,j+1,k)+b(1,j-1,k)-2.*b(1,j,k)) & /dy**2)) 3231 a(l,j,k) = a(1,j,k) 1000 continue return end