program gwdrive c Driver for the gwforce subroutine developed for SKYHI. parameter(nn=40,iz0=32,pi=3.14159,nk=1) real zp(nn+1),p(nn+1),uj(nn+1),rj(nn+1),bj(nn+1) real gwf(nn+1),gwd(nn),ked(nn+1),Ded(nn) real lat,f0,Bt,Bw,cw,ki,cmax,umax,dc,uin,Bi,dzt integer nc,flag flag=1 Bt=.01 Bw=.1 cw=20. lat=40. f0=2.*sin(pi*lat/180.)*7.292e-5 print *,'f0 =',f0 umax=0. open(unit=4,status='unknown',file='skyhi_in.dat') read(4,80) (zp(i+1),i=1,nn) read(4,80) (p(i+1),i=1,nn) read(4,80) (uj(i+1),i=1,nn) read(4,80) (rj(i+1),i=1,nn) read(4,80) (bj(i+1),i=1,nn) 80 format (1P6E13.5) dzt=zp(2)-zp(3) zp(1)=zp(2)+dzt uj(1)=2.*uj(2)-uj(3) rj(1)=rj(2)*rj(2)/rj(3) bj(1)=bj(2) print *,'Source Level Pressure (mb)=',p(iz0) open(unit=10,status='new',file='force_test.d', & access='sequential',form='formatted') 90 format (6G13.5) do i=1,nk Bi=Bt/nk ki=2.*pi/(30*(10.**i))/1.e3 print *,'i=',i print *,'k=',ki call gwfc(flag,Bi,Bw,cw,ki,rj,uj,bj,zp,f0,iz0,gwf,ked) do j=2,nn+1 gwd(j-1)=gwd(j-1)+gwf(j) Ded(j-1)=Ded(j-1)+ked(j) enddo enddo c Write out the results to a file. write(10,80) (gwd(i),i=1,nn) write(10,80) (zp(i),i=2,nn+1) write(10,80) (uj(i),i=2,nn+1) write(10,80) (rj(i),i=2,nn+1) write(10,80) (bj(i),i=2,nn+1) write(10,80) (p(i),i=2,nn+1) write(10,80) (Ded(i),i=1,nn) end c ------------------------------------------------------------------------- subroutine gwfc(flag,Bt,Bw,cw,k,rho,u,bf,z,f0,iz0,gwf,ked) c Version for SKYHI -- 10 March 1999 -- M.J. Alexander c Computes the gravity wave-driven-forcing given vertical profiles c of wind, density, and buoyancy frequency. The parameterization c is described in Alexander and Dunkerton [1998]. All quantities c in mks units unless otherwise specified. c c The gravity wave source is placed at the iz0 altitude index. c The forcing is computed at each level higher than this. c------------------------------------------------------------------ c INPUT: c A. Gravity wave source parameters (assumed Gaussian shape in c0) c flag = 1 for peak flux at c0=0 c flag = 0 for peak flux at ci=0 c Bt = sum|momentum flux| for all +/-c (kg/m/s^2) c Bw = amplitude for the spectrum (m^2/s^2) ~ u'w' c cw = half-width for the spectrum in c (m/s) c k = horizontal wavenumber (/m) = 2Pi/(wavelength) c dc = spectral resolution (m/s) c nc = number of spectral points (symmetric about c0=0) c B. Atmospheric variables c z = altitude array (km) c rho(z) = density (kg/m^3) c u(z) = horizontal wind (m/s) c bf(z) = buoyancy frequency (/s) c C. Other c f0 = Coriolis parameter (/s) = 1.458e-4 * sin(latitude) c OUTPUT: c gwf(z) = gravity-wave-driven force on the mean flow (m/s/day) c ked(z) = eddy diffusion coefficient (m^2/s) c (Note: ked should be reduced by some fraction < 1 before using.) c------------------------------------------------------------------ c INTERNAL VARIABLES: c nz = number of vertical grid points c iz0 = source level index c A. Variables defined at the lower boundary: c B0(c0) = amplitude spectrum as a function of c0 (m/s)^2 c c0(nc) = ground-relative phase speed c Bsum = total mag of flux at source level. c eps = intermittency factor c u0 = wind at the source level c r0 = density at the source level c cmax = maximum phase speed in the source spectrum c cc = phase speed where the source function peaks c ci = intrinsic phase speed c cmsk = 0 if ci=0 c flg0 = 0 if ci=0 occurs in the phase speed array c cj0 = index of the last negative phase speed c B. Variables defined at each grid level: c omc = critical frequency that marks total internal reflection c alp2 = scale height factor 1/(4H)^2 c icm = max and min phase speed indices not yet reflected c rmsk = 0 if wave has been reflected, 1 if not reflected c msk(c0) = 0 if wave breaks or reflects below, 1 if still propagating c C. Variables used in the loop through phase speed j c c0j=c0(j), B0j=B0(j) c c D. Variables used in the computation of the vertical gradient c dfac = conversion factor to m/s/day c rhof = density factor c dz = vertical increment c------------------------------------------------------------------ parameter(nz=41,nc=333,dc=.6) integer iz0 real z(nz),rho(nz),u(nz),bf(nz),gwf(nz),ked(nz) real fm,fe,Bsum real Nb,Hb,rb,ub,alp2,k2,f2,fac,Foc,u0,N0,r0,H0,omc real B0(nc),c0(nc),Bt,Bw,cw,k,f0 real c,ci,cc,cmax,c0j,B0j,cj0,L2 real dz,dBf,eps,dfac,rhof integer flag,rflg,flg0,rmsk,cmsk,icm(2),ict(2),sgn integer msk(nc) c Lower boundary set up u0=u(iz0) N0=bf(iz0) r0=rho(iz0) H0=-(z(iz0-1)-z(iz0))/alog(rho(iz0-1)/r0) flg0=1 do i=iz0,nz gwf(i)=0. enddo cc=u0*(1-flag) cmax=(nc-1)*dc/2. fac=k/2./N0 k2=k*k f2=f0*f0 alp2=1./(4.e6*(H0*H0)) omc=sqrt((N0*N0*k2)/(k2+alp2)) rflg=1 Bsum=0. rmsk=0 L2=alog(2.) c Loop through phase speeds do j=1,nc msk(j)=1 cmsk=1 c0j=(j-1)*dc-cmax c0(j)=c0j ci=c0j-u0 c=c0j*flag+ci*(1-flag) if (ci.eq.0) cmsk=0 sgn=cmsk*ci/(abs(ci)+(1-cmsk)) B0j=sgn*(Bw*exp(-L2*(c/cw)**2)) B0(j)=B0j if (sgn.le.0) j0=j if ((c0j.gt.u0).and.(rflg.gt.0)) rflg=0 if (c0j.eq.u0) msk(j)=0 if ((abs(c0j-u0)*k).lt.omc) then rmsk=1 icm(rflg+1)=j rflg=0 else rmsk=0 endif msk(j)=msk(j)*rmsk if (msk(j).eq.1) then Foc=B0j/(c0j-u0)**3 if ((Foc-fac).ge.0.0) msk(j)=0 endif Bsum=Bsum+abs(B0j) enddo if (Bsum.eq.0.0) then print *,'Zero flux input' return endif cj0=(j0-1)*dc-cmax if (cj0.eq.u0) flg0=0 c Intermittency factor eps=Bt/(r0*Bsum) c Conversion factor for /km/s to /m/day dfac=24.*3600./1.e3 c print *,'Intermittency=',eps,' dfac=',dfac c Altitude loop do i=iz0-1,1,-1 fm=0. fe=0. dz=z(i)-z(i+1) Hb=-dz/alog(rho(i)/rho(i+1)) rb=rho(i) rbh=sqrt(rb*rho(i+1)) Nb=bf(i) ub=u(i) fac=rb/r0*k/2./Nb alp2=1./(4e6*Hb*Hb) omc=sqrt((Nb*Nb*k2)/(k2+alp2)) rflg=1 c Phase speed loop do j=icm(2),icm(1) if (msk(j).eq.1) then B0j=B0(j) c0j=c0(j) if ((c0j.gt.ub).and.(rflg.gt.0)) then ict(rflg+1)=j0+flg0 rflg=0 endif if ((abs(c0j-ub)*k).lt.omc) then rmsk=1 ict(rflg+1)=j rflg=0 else rmsk=0 endif endif msk(j)=msk(j)*rmsk if (c0j.eq.ub) msk(j)=0 if (msk(j).eq.1) then Foc=B0j/(c0j-ub)**3 if((Foc-fac.ge.0).or.((c0j-u0)*(c0j-ub).le.0))msk(j)=0 fm=fm+(1-msk(j))*B0j mske=(1-msk(j)) c if the breaking level and critical level are not resolved c then no eddy diffusion if (((c0j-u0)*(c0j-ub)).lt.0) mske=0 fe=fe+mske*(c0j-ub)*B0j endif enddo icm(1)=min(ict(1),icm(1)) icm(2)=max(ict(2),icm(2)) c Compute the force and eddy diffusion coeff. gwf(i)=(dfac*r0/rbh)*fm*eps/dz gwf(i+1)=(gwf(i+1)+gwf(i))/2. ked(i)=(r0/rbh)*fe*eps/(dz*1e3*Nb*Nb) ked(i+1)=(ked(i+1)+ked(i))/2. enddo return end