SUBROUTINE KWAVE(im10)
c
c   FILE KWAVE_BASE8QJ.F  FILE  TO  COMPUTE  momentum deposition from thermally
c      damped Kelvin wave for generation of SAO residual circulation/tracer transport
c      (ie, double peak structure) in 2D model following Gray and Pyle (1987)  (EF 6/13/96)
c
C   ALL CALCULTIONS DONE IN MKS UNITS  ==>  AND DOUBLE PRECISION, for 59 LEVELS
C

#include "com2d.h"


       REAL*8 ab1, ab2, aa3, kk3, cc3, sum1, xa, yf
       REAL*8 xlatr(37), thd(59), rr0(59), pp0(59), nnn(37,59)

       REAL*8 rad, om1, delz, hh, rr1, kap, xpi, qp1, xlat5(37)  
       REAL*8 pres59(59), zz59(59), rho59(59), gr59(59)
       REAL*8 XN2(37,59), UBAR1(37,59), DUDZ(37,59), XKZZT(37,59)
       REAL*8 FXTT(37,59), FK(37,59)

       COMMON/NEWSF1/XN2, UBAR1, DUDZ, XKZZT, FXTT, FK

       COMMON/NEWSFC/rad, om1, delz, hh, rr1, kap, xlat5, pres59,
     c               zz59, rho59, gr59


c                                      convert degrees latitude to radians
       do 101 ij=1,37 
 101      xlatr(ij) = xlat5(ij)*qp1


c  Slow damping rates as defined in Dunkerton, 1979
c          ab1 = 7.5257d-4
c          ab2 = .69897d0
c
c  Fast damping rates as defined in Dunkerton, 1979
c       ab1 = 1.7474d-3
c       ab2 = .30103d0

c  Modified Fast damping rates (Not quite as fast as the fast rates above)
       ab1 = 1.5d-3
       ab2 = .43103d0


c define thermal damping profile  (in 1/sec), set to zero in troposphere               ;; zz59 in km

       do 120 ik=1,5
 120      thd(ik) = 0.0d0

       do 121 ik=6,15
 121      thd(ik) = 0.4d-6 + 0.8d-6*((zz59(ik)-17.)/13.)

       do 122 ik=16,59
 122      thd(ik) = 1./(86400.*DEXP(2.3*(ab1*(zz59(ik)-50.)**2 + ab2)))

c
c  Seems like better results in HF FOR ALL SEASONS are achieved with a modified FAST damping rate,
c     and 1.5 times the Gray and Pyle tropopause mom fluxes (theirs was 7.e-3), as listed below
c     this seems a bit better overall compared with the Base8qk and Base8ql runs
c
c  Define vertical momentum flux at tropopause, zonal wave# (k=1), and phase speed

       aa3 = 1.5*7.0d-3
       kk3 = 1./rad
       cc3 = 50.d0


c    COMPUTE R (rr0)  and  P (pp0) AT EQUATOR only, damp out in latitude

       do 150 ik=1,59
          rr0(ik) = 0.0
          pp0(ik) = 0.0
       do 151 ij=1,37
          fk(ij,ik) = 0.0
 151      nnn(ij,ik) = DSQRT(xn2(ij,ik))
 150   continue


       do 201 ik=1,59
 201      rr0(ik) = thd(ik)*nnn(19,ik)/(kk3*((ubar1(19,ik)-cc3)**2))

       do 300 ik=2,59 
           sum1 = 0.0d0
              do 305 ii = 1,ik-1
                 xa = (rr0(ii) + rr0(ii+1))/2.
                 sum1 = sum1 + xa*delz
 305          continue
          pp0(ik) = sum1
 300  continue 

c                                                          ;;  FK(37,59) is du/dt (in m/sec2), eq (6) in GP87
       do 400 ik=8,59
 400      fk(19,ik) = 
     c      aa3*DEXP((zz59(ik)-zz59(8))*1000./hh)*rr0(ik)*DEXP(-pp0(ik))

       do 450 ik=8,59
       do 450 ij=1,37
          if (ij .ne. 19) then
             yf = DEXP(-2.*om1*rad*(xlatr(ij)**2)/(cc3 - ubar1(19,ik)))
             if (yf .lt. 1.d-20) yf = 0.0d0
             fk(ij,ik) = fk(19,ik)*yf
          endif
 450   continue


 	SAVE
	RETURN
	END