;c ------------------------------------------------------------------------- pro gwfc, flag,Bt,Bw,cw,k,rho,u,bf,z,f0,iz0,gwf,ked ;c ------------------------------------------------------------------------- ;c 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------------------------------------------------------------------ nz=n_elements(z) nc=333 dc=.6 gwf=fltarr(nz) & ked=fltarr(nz) B0=fltarr(nc) c0=fltarr(nc) icm=intarr(2) ict=intarr(2) msk=intarr(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 for i=iz0,nz-1 do gwf(i)=0. cc=u0*(1-flag) cmax=nc*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 for j=0,nc-1 do begin msk(j)=1 cmsk=1 c0j=j*dc-cmax c0(j)=c0j ci=c0j-u0 c=c0j*flag+ci*(1-flag) if (ci eq 0) then cmsk=0 sgn=cmsk*ci/(abs(ci)+(1-cmsk)) B0j=sgn*(Bw*exp(-L2*(c/cw)^2)) B0(j)=B0j if (sgn le 0) then j0=j if ((c0j gt u0) and (rflg gt 0)) then rflg=0 if (c0j eq u0) then msk(j)=0 if ((abs(ci)*k) lt omc) then begin rmsk=1 icm(rflg)=j rflg=0 endif else begin rmsk=0 endelse msk(j)=msk(j)*rmsk if (msk(j) eq 1) then begin Foc=B0j/ci^3 if ((Foc-fac) ge 0.0) then msk(j)=0 endif Bsum=Bsum+abs(B0j) endfor if (Bsum eq 0.0) then begin print,'Zero flux input' return endif cj0=j0*dc-cmax if (cj0 eq u0) then flg0=0 ;c Intermittency factor eps=Bt/(r0*Bsum) ;c Conversion factor for /km/s to /m/day dfac=24.*3600./1.e3 ;c Altitude loop for i=iz0+1,nz-1 do begin 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 for j=icm(1),icm(0) do begin if (msk(j) eq 1) then begin B0j=B0(j) c0j=c0(j) if ((c0j gt ub) and (rflg gt 0)) then begin ict(rflg)=j0+flg0 rflg=0 endif if ((abs(c0j-ub)*k) lt omc) then begin rmsk=1 ict(rflg)=j rflg=0 endif else begin rmsk=0 endelse endif msk(j)=msk(j)*rmsk if (c0j eq ub) then msk(j)=0 if (msk(j) eq 1) then begin Foc=B0j/(c0j-ub)^3 if((Foc-fac ge 0) or ((c0j-u0)*(c0j-ub) le 0)) then 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) then mske=0 fe=fe+mske*(c0j-ub)*B0j endif endfor icm(0)=ict(0)icm(1) ;c Compute the force and eddy diffusion coeff. ;c (Eddy diffusion here is simply energy dissipation rate over ;c buoyancy frequency squared.) 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. endfor return end ;