program BAROUT c c This program reads in the initial and forecast c streamfunction fields and computes the components c of the wind field.The option of computing the height c field is also built in .the height field is deduced c through the non linear reverse balance law. c geostrophic balance is assumed for the southern and c northern boundaries . The height field is not written c out to the output field . Wind and height field are c computed within subroutine ZFIELD. c c Originator : florida state university.last revision c L. Bounoua ,1994 c c input : c c (1) initial and forecast streamfunction fields from c logical unit 17. c c output : c c (1) logical unit 23 contains the initial and c sequences of forecast u,v fields for each forecast hour. c c subroutines : bound , const , jac , lap, relax, zfield. c called c c parameter nmap defines the number of sets of c streamfunction fields to be processed. c c tend : is forecast time of wbaro.for c nint : is output times from wbaro.for c nmap : is derived from the above 2 variables c parameter(l = 51, m = 21) real psi(l,m),u(l,m),v(l,m) real z(l,m),wrka(l,m),cor(m),dx(m) m1 = m-1 m2 = m-2 tend = 24 itend = int(tend) nint = 12 nmap = tend/nint + 1 print *,' number of output psi ',nmap c call CONST (m,dx,dy,cor,beta) call CONST (m,dx,dy,cor,beta,slat,dphi) c open(17,file='baro_out.dat',status='old') open(23,file='ufield_out',status='unknown') c do 10300 n = 1, nmap c c Read the streamfunction field . c read(17,878) ((psi(i,j),i=1,l),j=1,m) 878 format(6e13.6) c c subroutine 'zfield' is called to compute the c height field from the streamfunction field using c non-linear balance. the arrays u(i,j) and v(i,j) c contains the component u and v of the wind field . c zbar = 3000. zbar2 = 0. call ZFIELD (psi,u,v,l,m,dx,dy,cor,wrka,z,1.,1.,zbar,beta) do 10302 j = 1, m do 10302 i = 1, l-1 zbar2 = zbar2 + z(i,j) 10302 continue zbar2 = zbar2 / (m*(l-1)) do 10304 j = 1, m do 10304 i = 1, l z(i,j) = z(i,j)-zbar2 + zbar 10304 continue c c The southern and northern boundary values of the c u and v fields are set by linear interpolation. c do 10306 i = 1, l u(i,1) = 2.*u(i,2) - u(i,3) v(i,1) = 2.*v(i,2) - v(i,3) u(i,m) = 2.*u(i,m1) - u(i,m2) v(i,m) = 2.*v(i,m1) - v(i,m2) 10306 continue c c Write u,v components to unit 23 . c write (23,1000) ((u(i,j),i=1,l-6),j=1,m) write (23,1000) ((v(i,j),i=1,l-6),j=1,m) 10300 continue 1000 format(6e13.6) stop end subroutine CONST (m,dx,dy,cor,beta,slat,dphi) c c This subroutine defines the constants and c parameters required by program INFIELD . users of c this program need only to modify the parameter c statement in program INFIELD and this subroutine c to their requirements. c c Definitions : c c m : number of grid points in the meridional c direction. c dx(m) : grid size in meters for each row of c grid points in the zonal direction. c dy : grid size in meters for the meridional c direction. c cor(m) : Coriolis parameter for each row of grid c points. c beta : beta plane parameter (=df/dy). c c omega : earth rotation (rad/s) c er : earth radius (m ) c real cor(m),dx(m) dphi = 2.5 slat = 0.0 anlat = slat+(m-1)*dphi omega = 7.292e-05 er = 6.37122e06 pi = 4.0*atan(1.0) rad = pi/180. dy = dphi*111.1*1000.0 phi = slat do 10120 j= 1, m dx(j) = dy*cos(phi*rad) cor(j) = 2.0*omega*sin(phi*rad) phi = phi+dphi 10120 continue the = (slat+anlat)/2.0 beta = 2.0*omega*cos(the*rad)/er 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