c --------------------------------------------------------------------- c version by Alex Pogoreltsev, 2005 - see comments in the text c---------------------------------------------------------------------- subroutine grwaves_RSHU C ================================================================= C GRAVITY WAVE DRAG parameterization based on an analytical solution C (WKB approximation) of the vertical structure equation for the GWs C in the atmosphere with arbitrary background wind U0(z), V0(z) and C realistic radiative damping (Newtonian cooling). C---------------------------------------------------------------------- C Eddy diffusion coefficient is estimated using the idea of GW breaking C due to instability proposed by Lindzen (1981). C Heating/cooling rates are calculated following C Medvedev and Klaassen (2002), however, the efficiency of mechanical C energy conversion into heat is taken into account, and turbulent C Prandl number = 3 (Huang und Smith, 1991). C Heat flux and its divergence are calculated explicitly, C using the polarization relations for GWs and analytical solution. C ================================================================ include 'com_main.fc' real u0(kgit),v0(kgit),t0(kgit),ac0(kgit),ac1(kgit), $ Fheat1(kgit), heatc1(kgit), Fheatc(kgit), heatco(kgit), $ u0z(kgit),v0z(kgit),t0z(kgit), u0zz(kgit), v0zz(kgit), $ cphint(kgit), omb2(kgit),derkz2(kgit), $ debg0(kgit),debgm(kgit),geddym(kgit,igit,nb), $ degw(kgit),alfaNC(kgit),cphase(6),teta(8),w0(48), $ weight(6),ww(6,36) real forceU(kgit),forceV(kgit),forceT(kgit), $ work1(igit),work2(igit),work3(igit),work4(igit) c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c cphase = phase speed [m/s] c xl = horizontal wave length [m] c teta = azimuth of GW [deg] c ww = vertical velocity perturbations at lower bondary [m/s] c pr = 3! Prandl number (Huang and Smith, 1991) c eff = efficiency of impuls fluxes c Rg = gas constant c gam = Cp/Cv c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ data nwaves /48/, eff0/1./ data cphase/5.,10.,15.,20.,25.,30./, $ teta/0.,45.,90.,135.,180.,225.,270.,315./ c----------------------------------------------------- the amplitudes data Rg/287.04/,gam/1.4/,xl/300000./,pr/3./ data omega0/1.e-4/,alpha/0.8333/,gamma/1.5/,c1/10./,c2/4./ c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rkx=2.*pi/xl c------------------------------ Gavrilov do i=1,6 omegai=cphase(i)*rkx omrel = omegai/omega0 c1rel = cphase(i)/c1 c2rel = c2/cphase(i) weight(i)=cphase(i)/(1.+omrel**alpha)/(1.+c1rel+c2rel**gamma) end do do j=1,nb do i=1,6 c--------------------------- for January c wampl0 = 0.012*(1.+1./3.*tanh((float(j)-18.5)/6.)) c--------------------------- for July c wampl0 = 0.015*(1.-1./3.*tanh((float(j)-18.5)/6.)) c-------------------------- OLD parameters c wampl0 = 0.01-0.003*tanh((float(j)-18.5)/9.) c wampl0 = 0.015-0.005*tanh((float(j)-18.5)/9.) c--------------------------- for April c wampl0 = 0.015 wampl0 = 0.0125 ww(i,j)=wampl0*weight(i) end do end do c-------------------------------------------------------------------- c h = 7000. c Cp = Rg*gam/(gam-1.) c------------------- Prandtl numbers for molecular duffusion Pr_mol=4.*gam/(9.*gam-5.) eff = eff0/float(nwaves)*(1.-exp(-float(nsec)/5./86400.)) efheat=0.5 fdec = EXP(-dz/h) kgitm1=kgit-1 kgitm2=kgit-2 c------------------------------------------ COMMA to GW levels do k=1,kgitm1 kp1=k+1 alfaNC(k)=(alfa (k)+alfa (kp1))/2. debg0 (k)=(geddy0(k)+geddy0(kp1))/2. end do c-------------------------------------------------- LATITUDES do 1000 j=1,nb do nt=1,8 w0(nt) = ww(1,j) w0(nt+ 8) = ww(2,j) w0(nt+16) = ww(3,j) w0(nt+24) = ww(4,j) w0(nt+32) = ww(5,j) w0(nt+40) = ww(6,j) end do c-------------------------------------------------- LONGITUDES do 1000 i= 1,igit do k= 1,kgit forceU(k)= 0. forceV(k)= 0. forceT(k)= 0. geddy(k,i,j)=geddy0(k) rheat=0.015/rm(k)*tn1(j,k,i)**0.667 rvisk=Pr_mol/C_p*rheat geddym(k,i,j)=rvisk/rhos/rou(k) end do c------------------------------- background profiles do k=1,kgitm1 kp1=k+1 wzonk =un1(j,k,i) wmerk =vn1(j,k,i) tempk =tn1(j,k,i) wzonk1=un1(j,kp1,i) wmerk1=vn1(j,kp1,i) tempk1=tn1(j,kp1,i) u0 (k)=(wzonk1+wzonk)/2. u0z(k)=(wzonk1-wzonk)/dz v0 (k)=(wmerk1+wmerk)/2. v0z(k)=(wmerk1-wmerk)/dz t0 (k)=(tempk1+tempk)/2. t0z(k)=(tempk1-tempk)/dz omb2(k) = Rg/h*((gam-1.)/gam*t0(k)/h+t0z(k)) debgm(k)= (geddym(k,i,j)+geddym(kp1,i,j))/2. enddo do k=2,kgit-2 km1=k-1 kp1=k+1 u0zz(k)=(u0(km1)-2.*u0(k)+u0(kp1))/dz/dz v0zz(k)=(v0(km1)-2.*v0(k)+v0(kp1))/dz/dz end do c------------------------------- over all waves considered do 2000 nw=1,nwaves do k= 1,kgit degw(k)=0. Fheat1(k)=0. heatc1(k)=0. Fheatc(k)=0. heatco(k)=0. end do mc= int((nw-1)/7.9) + 1 mpha= nw - (mc-1)*8 cost= cos(pi*teta(mpha)/180.) sint= sin(pi*teta(mpha)/180.) c----------------------------- critical level (kcrit) do k=4,kgitm1 if (omb2(k).le.0.) go to 10 c----------------------------- intrinsic phase speed if(k.eq.4)cphkm1 =cphase(mc)-u0(3)*cost-v0(3)*sint cphint(k)=cphase(mc)-u0(k)*cost-v0(k)*sint if(abs(cphint(k)).le.1.) go to 10 sq= cphint(k)*cphkm1 C-------------------------------------------------------------- C the perturbation approach is applicable if only the radiative C damping is weak, so "alfaNC*SQRT(omb2)*h/rkx/cphint/cphint" C has to be small (<< 1) in other cases we assume that we are C close to the critical level C-------------------------------------------------------------- geddef=debg0(k)*(1.+1./pr)+debgm(k)*(1.+1./Pr_mol) alfaEF=alfaNC(k)+geddef*omb2(k)/cphint(k)/cphint(k) dcrit = alfaEF*SQRT(omb2(k))*h/rkx if (sq.le.dcrit) go to 10 cphkm1=cphint(k) end do 10 kcrit = k-1 if (kcrit.le.4) go to 2000 kcrtp1=kcrit+1 kcrtp2=kcrit+2 wamkm1 = w0(nw) C---------------- the vertical wave number Kz at the lower boundary C (altitude about of 10 km) rkz0 = SQRT(omb2(4)/cphint(4)/cphint(4)) rkzkm1 = rkz0 ompkm1 = rkx*cphint(4) realL1 = -1./h do 3000 k=5,kcrit omb = SQRT(omb2(k)) omega = rkx*cphase(mc) ompls = rkx*cphint(k ) C-------------------------------------------------------------------- C Eqn. [d^2/dz^2+L(z)*d/dz+M(z)]W=0 C-------------------------------------------------------------------- realL = -1./h omplsz= -rkx*(u0z (k)*cost+v0z (k)*sint) omplzz= -rkx*(u0zz(k)*cost+v0zz(k)*sint) realM = rkx*rkx*omb2(k)/ompls/ompls-omplsz/ompls/h-omplzz/ompls C---------------------- W.K.B. solution ----------------------------- C rkz2 = realM-realL*realL/4.-d(realL)/dz/2. C-------------------------------------------------------------------- rkz2 = realM-realL*realL/4. if(rkz2.gt.0.) go to 50 wampl = wamkm1 rkz=rkzkm1 go to 40 50 continue rkz = sqrt(rkz2) C------------------- correction to the L(z) due to Newtonian cooling C and background diffusion (has to be < 1/h) geddef=debg0(k)*(1.+1./pr)+debgm(k)*(1.+1./Pr_mol) realL = realL+ $ rkz*(alfaNC(k)+geddef*rkz*rkz)/abs(ompls) C------------------------- W.K.B. solution for amplitude !W'(z)! wampl = wamkm1*SQRT(rkzkm1/rkz)*EXP(-(realL+realL1)*dz/4.) wbreak = abs(ompls)/rkz if(wampl.le.wbreak) go to 40 wampl = wbreak c--------------- Eddy diffusion coefficient due to breaking of GW rkzave = (rkz +rkzkm1)/2. omplav = (ompls+ompkm1)/2. C-------------------------------------------------- COMMA levels geddef=geddy0(k)*(1.+1./pr)+geddym(k,i,j)*(1.+1./Pr_mol) degw(k) = abs(omplav)/(1.+1./pr)/rkzave/rkzave/rkzave*(1./h- $ rkzave*(alfa(k)+geddef*rkzave*rkzave)/abs(omplav)- $ 2.*ALOG(wampl/wamkm1*SQRT(rkz/rkzkm1))/dz) if(degw(k).gt.1000.)degw(k)=1000. 40 continue C--------------------------------------------------------------- C Force per unit mass due to divergency of GW momentum flux C The simplified Eqn. U'=i/Kx*(d/dz-1/h)W' has been used, C 0.5 is due to averaging over period ave.(W'W')= 0.5*!W'!^2 C--------------------------------------------------------------- C-------------------------------------------------- COMMA levels ac0(k) = -0.5/rkx* $ ((rkz*wampl*wampl-rkzkm1*wamkm1*wamkm1)/dz- $ (rkz*wampl*wampl+rkzkm1*wamkm1*wamkm1)/2./h) if(cphint(3).lt.0.) ac0(k)=-ac0(k) C--------------------------------------------------------------- C Heating/Cooling terms and flux of heat due to GW C--------------------------------------------------------------- wampla = (wampl+wamkm1)/2. geddef = degw(k)+geddy0(k)+geddym(k,i,j) if(geddef.gt.1000.)geddef=1000. heatc1(k) = efheat*wampla*wampla*rkzave*rkzave/2./C_p/rkx/rkx* $ geddef*rkzave*rkzave derkz2(k) = geddef*rkzave*rkzave C---------------------------------------------------- GWs levels Fheat1(k) = - wampl*wampl*h*omb2(k)/2./Rg/ompls/ompls wamkm1 = wampl rkzkm1 = rkz ompkm1 = ompls realL1 = realL 3000 continue do k=1,kgit-1 C--------------------------- COMMA to GWs levels Fheat1(k) = Fheat1(k)*(alfaNC(k)+(derkz2(k+1)+derkz2(k))/2./pr) end do do k=2,kgit C--------------------------- GWs to COMMA levels heatc1(k) = heatc1(k)+(Fheat1(k)+Fheat1(k-1))/2./gam/h end do C------------------- Smoothing of heating/cooling and flux terms if(kcrit.le.kgitm1) then ac0 (kcrtp1)=ac0 (kcrit)*fdec Fheat1(kcrtp1)=Fheat1(kcrit)*fdec heatc1(kcrtp1)=heatc1(kcrit)*fdec degw (kcrtp1)=degw (kcrit)*fdec end if if(kcrit.le.kgitm2) then ac0 (kcrtp2)=ac0 (kcrtp1)*fdec Fheat1(kcrtp2)=Fheat1(kcrtp1)*fdec heatc1(kcrtp2)=heatc1(kcrtp1)*fdec degw (kcrtp2)=degw (kcrtp1)*fdec end if c--------------------------------- smoothing in altitude do k=8,kgitm2 km1=k-1 km2=k-2 kp1=k+1 kp2=k+2 c--------------------------------------- strong smoothing version ac1(k)=(ac0(km2)+ac0(km1)+ac0(k)+ac0(kp2)+ac0(kp1))/5. Fheatc(k)=(Fheat1(km2)+Fheat1(km1)+Fheat1(k)+ $ Fheat1(kp2)+Fheat1(kp1))/5. heatco(k)=(heatc1(km2)+heatc1(km1)+heatc1(k)+ $ heatc1(kp2)+heatc1(kp1))/5. geddy(k,i,j)=geddy(k,i,j)+ $ eff*(degw(km2)+degw(km1)+degw(k)+ $ degw(kp2)+degw(kp1))/5. c--------------------------------------- weak smoothing version c ac1(k)=(ac0(km1)+ac0(k)+ac0(kp1))/3. c Fheatc(k)=(Fheat1(km1)+Fheat1(k)+Fheat1(kp1))/3. c heatco(k)=(heatc1(km1)+heatc1(k)+heatc1(kp1))/3. c geddy(k,i,j)=geddy(k,i,j)+ c $ eff*(degw(km1)+degw(k)+degw(kp1))/3. c------------------------------------------------------------- if(geddy(k,i,j).gt.1000.) geddy(k,i,j)=1000. end do c------------------------------------------------------------------ ac1(kgitm1)=(ac0(kgitm2)+ac0(kgitm1)+ac0(kgit))/3. Fheatc(kgitm1)=(Fheat1(kgitm2)+Fheat1(kgitm1)+Fheat1(kgit))/3. heatco(kgitm1)=(heatc1(kgitm2)+heatc1(kgitm1)+heatc1(kgit))/3. do k=5,kcrit km1=k-1 heatGW = heatco(k)-(Fheatc(k)-Fheatc(km1))/dz forceU(k) = forceU(k) + eff*ac1(k)*cost forceV(k) = forceV(k) + eff*ac1(k)*sint forceT(k) = forceT(k) + eff*heatGW end do forceU(kcrtp1)= forceU(kcrit)*fdec forceV(kcrtp1)= forceV(kcrit)*fdec forceT(kcrtp1)= forceT(kcrit)*fdec 2000 continue do k=1,kgit fgru(k,i,j)= forceU(k) fgrv(k,i,j)= forceV(k) fgrt(k,i,j)= forceT(k) end do 1000 continue c------------------------------------ smoothing in latitude c Alex Pogoreltsev 13.03.2007 c---------------------------------------------------------- do i=1,igit do k=8,kgit do j=1,nb work1(j)=geddy(k,i,j) work2(j)=fgru (k,i,j) work3(j)=fgrv (k,i,j) work4(j)=fgrt (k,i,j) end do do j=2,nb-1 jp1=j+1 jm1=j-1 geddy(k,i,j)= (work1(jm1)+work1(j)+work1(jp1))/3. fgru (k,i,j)= (work2(jm1)+work2(j)+work2(jp1))/3. fgrv (k,i,j)= (work3(jm1)+work3(j)+work3(jp1))/3. fgrt (k,i,j)= (work4(jm1)+work4(j)+work4(jp1))/3. end do end do end do c------------------------------------ smoothing in longitude do j=1,nb c do k=4,kgit do k=8,kgit do i=1,igit work1(i)=geddy(k,i,j) work2(i)=fgru (k,i,j) work3(i)=fgrv (k,i,j) work4(i)=fgrt (k,i,j) end do do i=1,igit ip1=i+1 ip2=i+2 im1=i-1 im2=i-2 c--------------------------- if(i.eq.1) then im1=igit im2=igit-1 end if if(i.eq.2)im2=igit c--------------------------- if(i.eq.igit) then ip1=1 ip2=2 end if if(i.eq.igit-1)ip2=1 geddy(k,i,j)= (work1(im2)+work1(im1)+work1(i)+ $ work1(ip2)+work1(ip1))/5. fgru (k,i,j)= (work2(im2)+work2(im1)+work2(i)+ $ work2(ip2)+work2(ip1))/5. fgrv (k,i,j)= (work3(im2)+work3(im1)+work3(i)+ $ work3(ip2)+work3(ip1))/5. fgrt (k,i,j)= (work4(im2)+work4(im1)+work4(i)+ $ work4(ip2)+work4(ip1))/5. end do end do end do return end