program QG2LW c c this program is used to solve quasi-geostrophic 2 - level omega c equation incorporating subgrid scale stable condensation heating. c parameter (l=21,m=13,np=7,l1=l-1,l2=l-2,m1=m-1,m2=m-2) c c declare variables. c real bb(l,m),q(l,m),es(l,m),qs(l,m),phi(m),dx(l),cor(m) real omega(l,m),omegah(l,m),sigmah(l,m),a(l,m),b(l,m) real force(l,m),fzeta(l,m),ftheta(l,m),c(l,m),d(l,m) real datag (l,m,np),datat(l,m,np),datar(l,m,np) real z500(l,m),z1000(l,m),t750(l,m),rh750(l,m) c c define constants. c data f0,g,alfa,eps/1.e-04,9.81,-1.6,.1e-08/ data sigma,dy,dphi,dp/.02,2.5e05,2.5,500./ pi = 4.0*atan(1.0) phi(1) = -15. c1 = pi/180. c c open input files . c open (20,file='geop21.dat',status='old') open (30,file='temp21.dat',status='old') open (40,file='rel21.dat ',status='old') c c read multilevel input data c 15 format(10f8.2) c do 3200 ip = 1, np read (20,15) ((datag(i,j,ip),i=1,l),j=1,m) read (30,15) ((datat(i,j,ip),i=1,l),j=1,m) read (40,15) ((datar(i,j,ip),i=1,l),j=1,m) 3200 continue do 3202 i = 1, l do 3202 j = 1, m c c extract geopotential(meters) at 1000 and 500 mb . c z1000 (i,j) = datag (i,j,1) z500 (i,j) = datag (i,j,4) c c extract temperature(kelvin) at 750 mb assumed to be =(temp. at 700)+2. c t750 (i,j) = datat(i,j,3) + 2. c c extract relative humidity(%) at 750 mb assumed to be the same as 700 mb. c rh750 (i,j) = datar(i,j,3) 3202 continue c c compute the grid distance and Coriolis parameter . c do 3204 j = 2, m 3204 phi(j) = phi(j-1) + dphi do 3206 j = 1, m dx(j) = dy*cos(phi(j)*c1) 3206 cor(j) = 2.*7.29/1.e05*sin(phi(j)*c1) c c compute streamfunctions at level 1000 and 500 mb.for midlatitude c synoptic scales motions the streamfunction is approximated by : c psi = g*z/fo . c do 3208 j = 1, m do 3208 i = 1, l z1000 (i,j) = z1000(i,j)*g/f0 z500 (i,j) = z500(i,j)*g/f0 3208 continue c c differential vorticity advection. c perform the laplacian of the streamfunction c and compute the jacobian of the streamfunction c and the absolute vorticity at both levels c call LAP (a,z1000(1,1),dx,dy,l,m,l1,m1,l2,m2) call LAP (b,z500 (1,1),dx,dy,l,m,l1,m1,l2,m2) c do 3210 j = 1, m do 3210 i = 1, l a(i,j) = a(i,j) + cor(j) 3210 b(i,j) = b(i,j) + cor(j) c call JAC (c,z1000(1,1),a,dx,dy,l,m,l1,m1,l2,m2) call JAC (d, z500(1,1),b,dx,dy,l,m,l1,m1,l2,m2) c c forcing due to differential vorticity advection c do 3212 j = 1, m do 3212 i = 1, l sigmah (i,j) = sigma 3212 fzeta (i,j) = (d(i,j)-c(i,j))/dp*f0 c c laplacian of thermal advection. c performs the jacobian of streamfunction and vertical gradient c of streamfunction then form the laplacian c do 3214 j = 1, m do 3214 i = 1, l c(i,j) = (z500 (i,j) - z1000(i,j))/dp 3214 d(i,j) = (z1000(i,j) + z500 (i,j))*.5 c call JAC (a,d,c,dx,dy,l,m,l1,m1,l2,m2) call LAP (b,a,dx,dy,l,m,l1,m1,l2,m2) c c forcing due to laplacian of thermal advection c do 3216 j = 1, m do 3216 i = 1, l 3216 ftheta(i,j) = -f0*b(i,j) c c compute omega due to each forcing component c through relaxation method .scale omega. c call RELAXW (omega,fzeta,sigmah,l,m,dx,dy,l2,m2,eps,alfa) c do 3218 i = 1, l do 3218 j = 1, m 3218 omega (i,j) = omega(i,j)*1e+05 print *,((omega(i,j),j=7,7),i=1,l,6) do 3220 j = 1,m do 3220 i = 1,l a(i,j) = omega(i,j) 3220 continue c call RELAXW (omega,ftheta,sigmah,l,m,dx,dy,l2,m2,eps,alfa) c do 3222 i = 1, l do 3222 j = 1, m 3222 omega(i,j) = omega(i,j)*1e+05 print *,((omega(i,j),j=7,7),i=1,l,6) c c compute total forcing and omega (dry) due to both c forcing contributions c do 3224 j = 1, m do 3224 i = 1, l omega(i,j) = omega(i,j) + a(i,j) 3224 force(i,j) = fzeta(i,j) + ftheta(i,j) do 3226 j = 1, m do 3226 i = 1, l if ( omega(i,j).ge.0. ) go to 70 c c stable heating is invoked where rising motion is present, c the modified stability parameter is computed by sigmal. c call SIGMAL (sigmah,t750,rh750,sigma,l,m,bb,q,es,qs) go to 3226 70 sigmah(i,j) = sigma 3226 continue c c compute omega for the moist case. c omega moist is the difference between the total c omega with heating and the dry omega . c call RELAXW (omegah,force,sigmah,l,m,dx,dy,l2,m2,eps,alfa) c c scale and display omega c do 3228 i = 1, l do 3228 j = 1, m 3228 omegah(i,j) = omegah(i,j)*1e+05 print *,((omegah(i,j),j=7,7),i=1,l,6) c do 3230 j = 1, m do 3230 i = 1, l omegah(i,j) = omegah(i,j) - omega(i,j) 3230 continue stop end subroutine RELAXW (u,y,sigmah,l,m,dx,dy,l2,m2,eps,alfa) c c This subroutine computes the vertical velocity using the relaxa c -tion technique. a parabola fit of the form w = a*p**2+b*p+c is c used . c real u(l,m),y(l,m),sigmah(l,m),dx(m) data ia,f0,dp/500,1.e-04,500./ ct = 8./3. do 3210 i = 1, l do 3210 j = 1, m 3210 u(i,j) = 0. npts = (l-4)*(m-4) lsc = -1 nsc = 0 24 nrel = 0 do 3212 i = 3, l2 do 3212 j = 3, m2 c bb = (u(i+1,j)*sigmah(i+1,j)+u(i-1,j)*sigmah(i-1,j)-2. & *u(i,j)*sigmah(i,j))/dx(j)**2 +(u(i,j+1)*sigmah(i,j+1) & +u(i,j-1)*sigmah(i,j-1)-2.*u(i,j)*sigmah(i,j))/dy**2 c cc = -ct*u(i,j)/(dp*dp)*f0*f0 dd = y(i,j) pmu = 2.*sigmah(i,j)*(1.+dx(j)*dx(j)/(dy*dy)) & +(ct*f0*f0/(dp*dp)*dx(j)*dx(j)) c r = (-bb-cc+dd)*dx(j)*dx(j)/pmu if (lsc - nsc) 29,29,30 29 u(i,j) = u(i,j)+alfa*r 30 d = abs (r) if (d-eps) 32,32,33 32 nrel = nrel + 1 33 continue 3212 continue nsc = nsc + 1 if (nrel - npts) 45,43,43 43 continue if (lsc-nsc) 44,46,46 44 lsc = nsc + 1 45 if (nsc - ia) 24,46,46 46 continue write(6,1000) write(6,1001)eps,alfa write(6,1002)nsc,nrel 1000 format (1h1/14x,22hprogress of relaxation) 1001 format (1h0/2x,8hepsilon=e11.3,10h, alpha=f7.3) 1002 format (1h-/2x,11hscan count=i5,19h, points relaxed=i4) 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 SIGMAL (sigmah,t,rh,sigmad,l,m,bb,q,es,qs) c c This subroutine computes the stability change c due to stable condensation heat release. c real sigmah(l,m),t(l,m),rh(l,m) real bb (l,m),q(l,m),es(l,m),qs(l,m) data ca,cb,to,eps,eso /17.27,35.86,273.16,0.622,6.11 / data p,cp,rho,sigmao/750.,1004.,0.65,0.002/ rg = 2008./7. hl = 4.18*597.3e03 do 3220 j = 1, m do 3220 i = 1, l es(i,j)= eso*exp(ca*(t(i,j)-to)/(t(i,j)-cb)) qs(i,j)= eps*es(i,j)/(p-es(i,j)) if (rh(i,j).le.rho) go to 50 c c stable heating is invoked only where the grid c relative humidity is larger than rho. c al = (rh(i,j)-rho)/(1.-rho) bn = (cp*t(i,j)-qs(i,j)*hl)/(eps*hl) bd = cp*rg*t(i,j)*t(i,j)/(qs(i,j)*hl*(eps+qs(i,j))) bb(i,j) = al*rg*t(i,j)*(1.-bn)/(p*(hl+bd)) sigmah(i,j) = sigmad-rg*hl*bb(i,j)/(cp*p) if (sigmah(i,j).lt.sigmao) sigmah(i,j) = sigmao go to 3220 50 sigmah(i,j) = sigmad 3220 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