!
      module mgw_module
!
! Gravity wave parameterization module.
! This code is modified from tgcm24, which initially used code from CCM,
!   to call the waccm2/cam gw module.
! 2/1/05 btf: merge Ray's mods from yr92ncep/yrrun (timegcm1.1) to make
!   timegcm1.2.
!
      use params_module,only: nlev,nlat,nlon,nlonp4,nlevp1
      use gw_share,only: pgwv,cpair,pdel,rpdel,pmid,pint,
     |  lnpint,swap2d,invert2d,pcnst,pnats,
     |  usorla,ulsorla,faca,facnh,facsh,facoron,facoros,uth,faca1,
     |  wspnh, wspsh
      use gw_drag,only: gw_intr
      use cons_module,only: calc_gw
      use addfld_module,only: addfld
      implicit none
!
      real :: 
     |  rayk(nlev),          ! rayleigh friction
     |  difk(nlevp1),
     |  dift(nlevp1),
     |  difv(nlevp1)         ! formerly xmue
!
! Arrays on subdomains (allocated by alloc_gw)
      real,allocatable,dimension(:,:) ::
     |  raykk,               ! rayleigh friction (nlevp1,lat0:lat1)
     |  gwsclh,              ! scale height      (nlevp1,lon0:lon1)
     |  difk_gw              ! output by gw_intr (nlevp1,lon0:lon1)
                             ! (not defined yet as of 10/03, see d() in gw_drag.F)
      real :: facnhe,facshe
!
! Diffusivity terms difkk, difkt, and difdv are 3d fields (see fields.F)
! are returned to sub dynamics by sub mgw.
!
      contains
!-----------------------------------------------------------------------
      subroutine backgrnd
!
! Called once per timestep from advance.
!
      use params_module,only: dz,zibot
      use init_module,only: iter,iday
      use cons_module,only: dt,pi
!
      real,parameter ::
     |  eddylb  = 4.0e-8,
     |  zplb    = -17.,
     |  eddymb1 = 6.0e-8,
     |  zpmb1   = -13.,
     |  eddymb2 = 8.0e-8,
     |  zpmb2   = -9.,
     |  eddytp  = 1.0e-7,
     |  zptp    = -6.5,
     |  eddtlb  = 4.0e-8,
     |  eddtmb1 = 6.0e-8,
     |  eddtmb2 = 8.0e-8,
     |  eddttp  = 1.0e-7
!
!    |  eddylb  = 5.0e-8,
!    |  zplb    = -17.,
!    |  eddymb1 = 1.0e-7,
!    |  zpmb1   = -13.,
!    |  eddymb2 = 3.0e-7,
!    |  zpmb2   = -9.,
!    |  eddytp  = 5.0e-7,
!    |  zptp    = -6.5,
!    |  eddtlb  = eddylb,
!    |  eddtmb1 = eddymb1,
!    |  eddtmb2 = eddymb2,
!    |  eddttp  = eddytp
!
! Rayleigh friction background damping time 1/116 days constant
! with altitude (1.E-7 S-1)
!
      real,parameter ::
     |  usorle  = -5.0,
     |  ulsorle = +2.0,
     |  usorls  = +3.0,
     |  ulsorls = -3.0,
!    |  face    = 0.004,
     |  face    = 0.015,
     |  facs    = 0.040,
     |  fac1e   = 0.005,
     |  fac1s   = 0.002
      real :: 
     |  raykle ,
     |  raykme1,
     |  raykme2,
     |  raykte ,
     |  raykls ,
     |  raykms1,
     |  raykms2,
     |  raykts
!
      real :: dday,day1,day2,rayklb,raykmb1,raykmb2,rayktb,
     |  ds,s,expds,col1(nlevp1)
      integer :: k,klb,kmb1,kmb2,kturb
!
	raykle  = 1.5e-7
	raykme1 = 1.0e-7
	raykme2 = 0.75e-7
	raykte  = 0.5e-7
	raykls  = 5.0e-7
	raykms1 = 5.0e-7
	raykms2 = 3.0e-7
	raykts  = 2.0e-7
!
! Calculate eddy diffusion and rayleigh friction profiles
!
      dday = float(iday)+amod(float(iter)*dt,86400.)/86400.
      uth = amod(float(iter)*dt,86400.)/3600.
      day1 = abs(cos(2.*pi*(dday-355.)/365.))
      day2 = cos(2.*pi*(dday-355.)/365.)
!
      rayklb = raykle+(raykls-raykle)*day1
      raykmb1 = raykme1+(raykms1-raykme1)*day1
      raykmb2 = raykme2+(raykms2-raykme2)*day1
      rayktb = raykte+(raykts-raykte)*day1
!
! This is coupled to mgw.f
      usorla  = usorle +(usorls -usorle )*day1
      ulsorla = ulsorle+(ulsorls-ulsorle)*day1
!     write(6,"('backgrnd: day1=',e12.4,' usorla=',e12.4,
!    |  ' ulsorla=',e12.4)") day1,usorla,ulsorla ! -1.468, 0.5
!
! This is coupled to mgwbgnd.f
      faca = face+(facs-face)*day1
      faca1 = fac1e+(fac1s-fac1e)*day1
!      facnh =1.0+0.5*day2
!      facsh =1.0-0.5*day2
       facnh =1.0
       facsh =1.0
!
       facnhe =1.0+0.5*day2
       facshe =1.0-0.6*day2
!
      wspnh = 25.+ 5.*DAY2
      wspsh = 25.- 5.*DAY2
!     facoron = (1.+day2)**2.
!     facoros = (1.-day2)**2.
      facoron = (1.+0.5*day2)
      facoros = (1.-0.5*day2)
!
      klb = 1
      kmb1 = (zpmb1-zibot)/dz+1.0001
      kmb2 = (zpmb2-zibot)/dz+1.0001
      kturb = (zptp-zibot)/dz+1.0001
      difk(1) = eddylb
      dift(1) = eddtlb
      difv(1) = eddylb
      rayk(1) = rayklb
!
      do k=2,kmb1
        difk(k) = eddylb*exp(alog(eddymb1/eddylb)*(k-klb)/(kmb1-klb))
        dift(k) = eddtlb*exp(alog(eddtmb1/eddtlb)*(k-klb)/(kmb1-klb))
        difv(k) = difk(k)
        rayk(k) = rayklb*exp(alog(raykmb1/rayklb)*(k-klb)/(kmb1-klb))
      enddo ! k=2,kmb1
!
      do k=kmb1,kmb2
        difk(k)=eddymb1*exp(alog(eddymb2/eddymb1)*(k-kmb1)/(kmb2-kmb1))
        dift(k)=eddtmb1*exp(alog(eddtmb2/eddtmb1)*(k-kmb1)/(kmb2-kmb1))
        difv(k) = difk(k)
        rayk(k)=raykmb1*exp(alog(raykmb2/raykmb1)*(k-kmb1)/(kmb2-kmb1))
      enddo ! k=kmb1,kmb2
!
      do k=kmb2+1,kturb
        difk(k)=eddymb2*exp(alog(eddytp/eddymb2)*(k-kmb2)/(kturb-kmb2))
        dift(k)=eddtmb2*exp(alog(eddttp/eddtmb2)*(k-kmb2)/(kturb-kmb2))
        difv(k) = difk(k)
        rayk(k)=raykmb2*exp(alog(rayktb/raykmb2)*(k-kmb2)/(kturb-kmb2))
      enddo ! k=kmb2+1,kturb
!
      do k=kturb+1,nlevp1
        col1(k) = zibot+(k-1)*dz
        difk(k) = eddytp*exp(zptp-col1(k))
        dift(k) = eddttp*exp(zptp-col1(k))
        difv(k) = difk(k)
      enddo ! k=kturb+1,nlevp1
!
      do k=kturb+1,nlev
        rayk(k) = rayktb*exp(zptp-col1(k))
!       rayk(k) = rayktb
      enddo ! k=kturb+1,nlev
!
!     write(6,"('  difk=',/,(6e12.4))") difk
!     write(6,"('  dift=',/,(6e12.4))") dift
!     write(6,"('  difv=',/,(6e12.4))") difv
      end subroutine backgrnd
!-----------------------------------------------------------------------
      subroutine gwsource(z,lev0,lev1,lon0,lon1,lat)
!
! Called from latitude loop in dynamics. 
!
      use cons_module,only: dtr
      use init_module,only: glat,zpint
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  z  ! geopotential
!
! Local:
      integer :: k,i,lonz
      real :: pi_g,al0,dlat0,rlat,zp80
!
! Garcia type analytic gwsource latitude distribution:
!
      real,parameter :: rayksh=1.0, rayknh=1.0
!
! Calculate gravity wave parameter range in pressure coord.
! Rayleigh friction from ccm2 specification:
!
!     real,parameter :: raykred = 0.
      real,parameter :: raykred = 0.10
!
      real :: 
     |  raykb(nlev),
     |  raykccm(nlev),
     |  ablat(nlat)
      real :: rayfrc(lev0:lev1,lon0:lon1) ! for diag
!
! lonz = arbitrary longitude for Z (hardwired at 37 in tgcm24)
! (Can set rakccm = 0. for comparison with tgcm24)
!
! 2/22/07 btf: Using Z at lonz==lon0 creates discontinuities in
!   raykk between mpi tasks at their lon0 boundaries. To prevent 
!   this, use zp instead of z to calculate raykccm. Assume 4 scale 
!   heights across the 80 km transition.
!
!     lonz = lon0
      zp80 = -10. ! zp at approx 80 km
      do k=lev0,lev1-1
        raykccm(k)=1.0/1.5*(1.-tanh((zpint(k)-zp80)/4.))/86400.*raykred
        raykb(k) = rayk(k)+raykccm(k)
      enddo
!
      pi_g  =4.*atan(1.)
      al0   =30.*pi_g/180.
      dlat0 =20.*pi_g/180.
      rlat = glat(lat)*dtr
!     ablat(lat) =
!    |  (0.5*(1.+tanh( (rlat-al0)/dlat0))
!    |  +0.5*(1.+tanh(-(rlat+al0)/dlat0)))
!     ablat(lat)=max(ablat(lat),0.1)
      ablat(lat) = 1.
      do k=lev0,lev1-1
        if (lat <= nlat/2) then ! south hem
          raykk(k,lat) = raykb(k)*rayksh*ablat(lat)
        else                ! north hem
          raykk(k,lat) = raykb(k)*rayknh*ablat(lat)
        endif
! Save raykk to secondary history in local rayfrc:
        rayfrc(k,:) = raykk(k,lat)
      enddo ! k=1,nlev-1
!     write(6,"('gwsource: lat=',i3,' raykk(:,lat)=',/,(6e12.4))")
!    |  lat,raykk(:,lat)
!
!     write(6,"('Completed gwsource: lat=',i3)") lat

!     write(6,"('gwsource: lat=',i3,' rayfrc(:,1)=',/,(6e12.4))")
!    |  lat,rayfrc(:,1) 
      call addfld('RAYK',' ',' ',rayfrc,
     |  'lev',lev0,lev1,'lon',lon0,lon1,lat)

      end subroutine gwsource
!-----------------------------------------------------------------------
      subroutine mgw(tn,un,vn,barm,z,cp,gwu,gwv,gwt,
     |  difkk,difkt,difkv,lev0,lev1,lon0,lon1,lat)
!
! Driver for gravity wave parameterization. This calls cam driver gw_intr.
! Called from lat loop in dynamics.
!
      use cons_module,only: dtr,gask,dt,grav
      use gw_share,only: prndtl,prndtlt,prndtlm
      use init_module,only: glat
      use init_module,only: iter,iday

!+Secondary GW forcing in thermosphere
      use sgwf_module,only: secondary_gwf
!-Secondary GW forcing in thermosphere

#ifdef MPI
      use mpi_module,only: mp_gatherlons_f3d,mp_scatterlons_f3d,mytidi
      implicit none
#else
      implicit none
      integer :: mytidi=0
#endif
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  tn,    ! neutral temperature (deg K)
     |  un,    ! zonal wind velocity (cm/s)
     |  vn,    ! meridional wind velocity (cm/s)
     |  barm,  ! mean molecular mass
     |  z,     ! geopotential
     |  cp     ! specific heat
!
! Outputs at current latitude:
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
     |  gwu,   ! gravity wave zonal tendency (output)
     |  gwv,   ! gravity wave meridional tendency (output)
     |  gwt,   ! gravity wave temperature tendency (output)
     |  difkk, ! eddy viscosity
     |  difkt, ! eddy thermal diffusion
     |  difkv  ! eddy diffusion for constituents
!
! Local:
      integer :: k,i,i0,i1,nk,nkm1,nlevs,nlons
      real :: rlat,zhi,zlow,hlat,dzhl
      real :: pi_g,al0,dlat0,al1,dlat1,flat_gw
      real :: flat_gw1,hlateq
      real :: usorl,ulsorl,uhsorh
      real,dimension(lev0:lev1,lon0:lon1) :: usaf
      real,dimension(lev0:lev1) :: dmpf
      real,dimension(lev0:lev1,lon0:lon1) :: dkk
      real,dimension(lev0:lev1,lon0:lon1) :: dkkt
      real,dimension(lev0:lev1,lon0:lon1) :: dkkv
      real,dimension(lev0:lev1,nlonp4) :: utglb
      real,dimension(lev0:lev1) :: utzm
      real,dimension(nlonp4) :: utpt
!+Secondary GW forcing in thermosphere
      real,dimension(lev0:lev1,lon0:lon1) :: gwsfx,gwsfy
!-Secondary GW forcing in thermosphere
!
! waccm/cam module inputs:
      real,dimension(lon0:lon1) :: 
     |  sgh,           ! standard deviation of orography
     |  pblh,          ! planetary boundary layer height
     |  landfrac       ! land fraction
      real,dimension(lon0:lon1,lev0:lev1) :: ! inputs for cam in cam dim order
     |  zm,            ! height at midpoints
     |  dse,           ! dry static energy at midpoints (state%s)
     |  pm,            ! pressure at midpoints (pmid)
     |  pd,            ! delta pressure (pdel)
     |  rpd            ! 1/pd
      real,dimension(lon0:lon1,lev0:lev1) :: ! inputs for cam in cam dim order
     |  pi,            ! pressure at interfaces (pint)
     |  lnpi           ! log(pi) (lnpint)
      real,dimension(1:lon1-lon0+1,lev0:lev1,pcnst+pnats) :: q_cam ! constituent(s)
      real,dimension(1:lon1-lon0+1,lev0:lev1) ::
     |  t_cam,u_cam,v_cam
      real,dimension(1:lon1-lon0+1,lev0:lev1) ::
     |  u_cam_0,v_cam_0               ! these are 0 background wind in zonal and meridional directions
!
! waccm/cam module output at cam dim order:
      real,dimension(1:lon1-lon0+1,lev0:lev1) ::
     |  ttend_cam, utend_cam, vtend_cam, difk_cam
      real,dimension(1:lon1-lon0+1,lev0:lev1,pcnst+pnats) :: qtend_cam
!
! waccm/cam module output at tgcm dim order:
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  ttend, utend, vtend, ttend_degday,utend_msday,vtend_msday,
     |  ttend1,utend1,vtend1,difk_gw1
!
      i0 = lon0 
      i1 = lon1
      nk = lev1-lev0+1
      nkm1 = nk-1
!
!+time-gcm+
      rlat = glat(lat)*dtr
      zhi  = 1.20e5-0.05e5*sin(abs(rlat))
      zlow = 1.00e5-0.05e5*sin(abs(rlat))
!     zhi  = 1.40e5-0.20e5*sin(abs(rlat))
!     zlow = 1.15e5-0.20e5*sin(abs(rlat))
!     zhi  = 1.60e5
!     zlow = 1.30e5
      hlat = 30./57.295
      dzhl = zhi-zlow
      dmpf = 0.

      u_cam_0 = 0
      v_cam_0 = 0
!
! Initialize tendencies.
!
      utend_cam = 0
      vtend_cam = 0
      ttend_cam = 0
      difk_cam = 0
      qtend_cam = 0
!-time-gcm-

!
! Include dependence on latitude:
! The lat function was obtained by RR Garcia as currently 
!   used in his 2D model
! [Added by F. Sassi on May 30, 1997, 
!  modified for timegcm1 by B. Foster, 10/3/03]
!
      pi_g  =4.*atan(1.)
      al0   =20.*pi_g/180.
      dlat0 =10.*pi_g/180.
      al1   =65.*pi_g/180.
      dlat1 =45.*pi_g/180.
	flat_gw=
     |  (0.5*(1.+tanh( (rlat-al0)/dlat0))*facnhe
     |  +0.5*(1.+tanh(-(rlat+al0)/dlat0))*facshe)
!    |       *abs(sin(2.*rlat))
       hlateq = 30.*dtr
       flat_gw1 = 0.04*exp(-(rlat/hlateq)**2)
!      flat_gw1 = 0.
       flat_gw = flat_gw+flat_gw1
!      flat_gw=1.
!
      usorl  = usorla*100./86400.*exp(-(rlat/hlat)**2)
      ulsorl =ulsorla*100./86400.*exp(-(rlat/hlat)**2)
      uhsorh = 0.0
!     write(6,"('mgw: glat=',f7.2,' rlat=',f6.2,' hlat=',f6.2,
!    |  ' usorl=',1pe12.4,' ulsorl=',e12.4)") glat(lat),rlat,hlat,
!    |  usorl,ulsorl
!
      do i=lon0,lon1
        do k=lev0+1,lev1-1
          gwsclh(k,i) = (gask*0.5*(tn(k-1,i)+tn(k,i))/
     |                  (grav*barm(k,i)))*1.e-2      ! meters 
        enddo ! k=lev0+1,lev1-1
        gwsclh(1,:) = 2.*gwsclh(2,:)-gwsclh(3,:)
        gwsclh(lev1,:) = 1.5*gwsclh(lev1-1,:)-0.5*gwsclh(lev1-2,:)
      enddo ! i=lon0,lon1

      sgh(:) = 10.5     ! standard deviation of orography
      pblh(:) = 0.      ! planetary boundary height
      landfrac(:) = 0.  ! land fraction
      q_cam(:,:,:) = 0. ! constituent(s) (# const = pcnst+pnats)
!
! zm = Heights at midpoints.
! dse = Dry static energy at midpoints.
!
      do i=lon0,lon1
        do k=lev0+1,lev1 ! 2,45 -> lev1-k+1 = 44,1,-1
          zm (i,lev1-k+1) = 0.5*(z(k,i)+z(k-1,i))*1.e-2
          dse(i,lev1-k+1) = cpair*tn(k,i)+grav*1.e-2*zm(i,lev1-k+1)
        enddo
      enddo
      zm (:,lev1) = zm (:,lev1-1)
      dse(:,lev1) = dse(:,lev1-1)
!
! Set subdomain pressures for gw_intr (pint,pmid,etc are defined at the
! full domain by set_gw_params (they are redundant in longitude)).
!
! 1/5/05 btf: changed pint,pmid, etc 2nd dimensions from pver to pver+1
!     to avoid out-of-bounds here.
! 5/15/06 btf: Fixed dimensions of pint and lnpi on right-hand side
!              (in gw_share they are dimensioned 0:pver+1
!
      pm(:,:)   = pmid  (lon0:lon1,lev0:lev1)
      pd(:,:)   = pdel  (lon0:lon1,lev0:lev1)
      rpd(:,:)  = rpdel (lon0:lon1,lev0:lev1)
      pi(:,:)   = pint  (lon0:lon1,0:lev1-1)
      lnpi(:,:) = lnpint(lon0:lon1,0:lev1-1)
!
! Set t,u,v inputs for gw_intr:
! real,dimension(1:lon1-lon0+1,lev0:lev1) :: t_cam,u_cam,v_cam
!
      do i=lon0,lon1
         do k=lev0,lev1
            t_cam(i-lon0+1,lev1-k+1) = tn(k,i)
            u_cam(i-lon0+1,lev1-k+1) = un(k,i)*1.e-2
            v_cam(i-lon0+1,lev1-k+1) = vn(k,i)*1.e-2   ! zonal index shift so that it starts from 1
         enddo
!
! Overwrite top level (cam k=1, gcm k=lev1) because historically in the 
! gcm, the bottom boundary of t,u,v (k=1) is stored in top slot (k=lev1).
! 5/4/06 btf: bottom boundary is no longer in top slot, but now top slot
!   contains spval, so let this overwrite stand:
        t_cam(:,1) = t_cam(:,2)
        u_cam(:,1) = u_cam(:,2)
        v_cam(:,1) = v_cam(:,2)
      enddo
!
      nlevs = lev1-lev0+1 ! # subdomain levels (pver  in waccm/cam module)
      nlons = lon1-lon0+1 ! # subdomain longs  (pcols in waccm/cam module)
!
!     call addfld('T_GCM',' ',' ',
!    |  swap2d(invert2d(t_cam,nlons,nlevs,2),nlons,nlevs),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('U_GCM',' ',' ',
!    |  swap2d(invert2d(u_cam,nlons,nlevs,2),nlons,nlevs),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('V_GCM',' ',' ',
!    |  swap2d(invert2d(v_cam,nlons,nlevs,2),nlons,nlevs),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Call gravity wave driver of waccm/cam1 parameterization. Pass 2-d tgcm 
!   arrays, with dimensions swapped (from (nlevs,nlons) to (nlons,nlevs)), 
!   and with the vertical dimension nlevs inverted (from bot->top to top->bot).
! (u,v are converted from cm to m before going into gw_intr)
! Output: utend_cam,vtend_cam,ttend_cam,difk_cam,qtend_cam
!
      if (calc_gw) then ! calc_gw is in cons.F
!
! First call gw_intr in meridional direction, zonal wind set to 0.
!
        call gw_intr(
     |    u_cam_0,v_cam,t_cam,q_cam,
     |    dse,zm,sgh,pblh,dt,landfrac,pm,pi,lnpi,pd,rpd,rlat,
     |    utend_cam,vtend_cam,ttend_cam,difk_cam,qtend_cam,
     |    nlons,nlev,lon0,lon1,lat)
!
! then call gw_intr in zonal direction, meridional wind set to 0.
!
        call gw_intr(
     |    u_cam,v_cam_0,t_cam,q_cam,
     |    dse,zm,sgh,pblh,dt,landfrac,pm,pi,
     |    lnpi,pd,rpd,rlat,
     |    utend_cam,vtend_cam,ttend_cam,difk_cam,qtend_cam,
     |    nlons,nlev,lon0,lon1,lat)
!
! Transform waccm output to swapped and inverted tgcm dims, 
!   store in ttend, utend, vtend.
!
        ttend(lev0:lev1-1,:)=swap2d(invert2d(ttend_cam,nlons,nlev
     |                       ,2),nlons,nlev)*flat_gw
        utend(lev0:lev1-1,:)=swap2d(invert2d(utend_cam,nlons,nlev
     |                       ,2),nlons,nlev)*flat_gw
        vtend(lev0:lev1-1,:)=swap2d(invert2d(vtend_cam,nlons,nlev
     |                       ,2),nlons,nlev)*flat_gw
        difk_gw=swap2d(invert2d(difk_cam,nlons,nlevs,2),nlons,
     |                        nlevs)*flat_gw


      do k = lev0,lev1-1
         if(zm(lon0,lev1-k+1).le.zlow) then
             dmpf(k) = 1.
          else if(zm(lon0,lev1-k+1).gt.zhi) then
             dmpf(k) = 0.
          else
             dmpf(k) = (cos(3.14159/2.*(zm(lon0,lev1-k+1)
     |                -zlow)/dzhl))**2
         endif
      enddo

      do i=lon0,lon1
         do k=lev0,lev1-1
            utglb(k,i) = utend(k,i)
         enddo
      enddo

#ifdef MPI
      call mp_gatherlons_f3d(utglb,lev0,lev1,lon0,lon1,1,1,1)
#endif

      if(mytidi == 0) then
        do k=lev0,lev1-1
           utzm(k) = sum(utglb(k,3:nlonp4-2))/nlon
           utpt(:) = utglb(k,:)-utzm(k)
!          utglb(k,:) = utzm(k) + utpt(:)*0.0075
	   utglb(k,:) = utzm(k) + utpt(:)*0.03
        enddo
      endif ! mytidi==0
#ifdef MPI
      call mp_scatterlons_f3d(utglb,lev0,lev1,lon0,lon1,1,1,1)
#endif

      do i=lon0,lon1
       do k = lev0,lev1-1
         utend(k,i) = utglb(k,i)
	 utend(k,i) = utend(k,i)*dmpf(k)
	 vtend(k,i) = vtend(k,i)*dmpf(k)
	 ttend(k,i) = ttend(k,i)*dmpf(k)
	 difk_gw(k,i) = difk_gw(k,i)*dmpf(k)
	 if(utend(k,i)*8.64e4.gt.150.) utend(k,i) = 150./8.64e4
	 if(vtend(k,i)*8.64e4.gt.150.) vtend(k,i) = 150./8.64e4
	 if(ttend(k,i)/1004.*8.64e4.gt. 2.0) ttend(k,i) =  2.0/8.64e4*1004.
	 if(utend(k,i)*8.64e4.lt.-150.) utend(k,i) = -150./8.64e4
	 if(vtend(k,i)*8.64e4.lt.-150.) vtend(k,i) = -150./8.64e4
	 if(ttend(k,i)/1004.*8.64e4.le. -2.0) ttend(k,i) = -2.0/8.64e4*1004.
	 if(difk_gw(k,i).gt.150.) difk_gw(k,i) = 150.
       enddo
      enddo
      difk_gw(lev1,:) = 0.
      utend(lev1,:) = 0.
      vtend(lev1,:) = 0.
      ttend(lev1,:) = 0.

!     call addfld('TTEND',' ',' ',ttend,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('UTEND',' ',' ',utend,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('VTEND',' ',' ',vtend,'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('DIFKGW',' ',' ',difk_gw,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Transfer to output arrays:

!+Apply secondary forcing in the thermosphere

      call secondary_gwf(z,lev0,lev1,lon0,lon1,lat,gwsfx,gwsfy)

!-apply secondary forcing in the thermosphere

        do i=lon0,lon1
          do k=lev0,lev1-1 
!
! U tendency:
            gwu(k,i) = utend(k,i)*1.e2+                      ! Total U-tend (cm/s/s)
     |        usorl*exp(1.-(z(k,i)-85.0e5)/17.e5-
     |              exp(  -(z(k,i)-85.0e5)/17.e5))+
     |       ulsorl*exp( -((z(k,i)-45.e5)/15.e5)**2)
     |       + gwsfx(k,i)*1.e2                               ! secondary forcing in thermosphere
            utend_msday(k,i) = utend(k,i)*86400.            ! m/s/day (for diag)
!
! Artificial momentum source for SAO (diagnostic)
	    usaf(k,i) =                                      ! momentum (m/s/day)
     |        (usorl*exp(1.-(z(k,i)-85.0e5)/17.e5-
     |               exp(  -(z(k,i)-85.0e5)/17.e5))+
     |        ulsorl*exp( -((z(k,i)-45.e5)/15.e5)**2))*8.64e2
!
! V tendency:
            gwv(k,i) = vtend(k,i)*1.e2    ! cm/s/s
     |                +gwsfy(k,i)*1.e2                      ! secondary forcing in thermosphere
            vtend_msday(k,i) = vtend(k,i)*86400.            ! m/s/day (for diag)
!
! Convert ttend from J/kg/s, as output by gw_drag, to K/day for diagnostic:
! cp(tgcm) = ergs/deg/gm (specific heat at constant pressure)
! (1 Joule = 1.e7 ergs, and 1kg = 1000 g)
!           gwt(k,i) = ttend(k,i)/cp(k,i)*1.e7/1000.*86400.    ! deg/day
!
! Leave gwt in ergs/gm/s for later use in dt.
! (10/03: in roble's tgcm24co2d/modsrc.jsolr, this conversion is in dt)
!
            gwt(k,i) = ttend(k,i)*1.e4                         ! ergs/gm/sec (to dt)
            ttend_degday(k,i) = ttend(k,i)/cp(k,i)*1.e4*86400. ! deg/day (for diag)
          enddo ! k=lev0,lev1 
        enddo ! i=lon0,lon1

      call addfld('CPGW',' ',' ',cp (lev0:lev1-1,i0:i1),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('GWU' ,' ',' ',gwu(lev0:lev1-1,i0:i1),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('GWV' ,' ',' ',gwv(lev0:lev1-1,i0:i1),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('GWT' ,' ',' ',ttend_degday(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!
      call addfld('UMOM'  ,' ',' ',usaf(lev0:lev1-1,i0:i1),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('UTNDCY',' ',' ',utend_msday(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('VTNDCY',' ',' ',vtend_msday(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
      call addfld('TTNDCY',' ',' ',ttend_degday(lev0:lev1-1,:),
     |  'lev',lev0,lev1-1,'lon',i0,i1,lat)

      else !  not calc_gw: waccm/ccm gw_intr not called
!       if (lat==1)
!    |    write(6,"('mgw: lat=',i3,' calc_gw=',l1,
!    |      ' setting gwu,v,t = 0.')") lat,calc_gw
        gwu = 0. ! whole-array op
        gwv = 0.
        gwt = 0.
      endif
!
      if (calc_gw) then
        do i=lon0,lon1
          do k=lev0,lev1
            difkk(k,i) =  difk(k)+difk_gw(k,i)/(gwsclh(k,i)**2)/
     |        prndtl

            difkt(k,i) = difk(k)+difk_gw(k,i)/
     |                   (gwsclh(k,i)**2)/prndtlt

            difkv(k,i) = difk(k)+difk_gw(k,i)/
     |                   (gwsclh(k,i)**2)/prndtlm
            dkk(k,i) = difkk(k,i)*gwsclh(k,i)**2
            dkkt(k,i) = difkt(k,i)*gwsclh(k,i)**2
            dkkv(k,i) = difkv(k,i)*gwsclh(k,i)**2

          enddo ! k=lev0,lev1

 !        write(6,"('mgw: lat=',i2,' difk=',/,(6e12.4))") lat,difk
 !        write(6,"('mgw: lat=',i2,' dift=',/,(6e12.4))") lat,dift

        enddo ! i=lon0,lon1
      else ! not calculating gw
        do k=lev0,lev1
          difkk(k,:) = difk(k) 
          difkt(k,:) = difk(k) 
          difkv(k,:) = difk(k) 
          dkk  (k,:) = difk(k) 
          dkkt (k,:) = difk(k) 
          dkkv (k,:) = difk(k) 
        enddo
      endif

!     write(6,"('mgw: lat=',i2,' flat_gw=',e12.4,' prndtl,t,m=',
!    |  3e12.4)") lat,flat_gw,prndtl,prndtlt,prndtlm

!     call addfld('GWDIFKK1',' ',' ',difkk(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('GWDIFKK2',' ',' ',dkk  (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('GWDIFKT',' ',' ',dkkt (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('GWDIFKV',' ',' ',dkkv (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!
      call addfld('DIFK_GW',' ',' ',difk_gw(:,i0:i1),
     |  'lev',lev0,lev1,'lon',i0,i1,lat)
      call addfld('DIFKK',' ',' ',dkk(:,i0:i1),
     |  'lev',lev0,lev1,'lon',i0,i1,lat)
      call addfld('DIFKT',' ',' ',dkkt(:,i0:i1),
     |  'lev',lev0,lev1,'lon',i0,i1,lat)
      call addfld('DIFKV',' ',' ',dkkv(:,i0:i1),
     |  'lev',lev0,lev1,'lon',i0,i1,lat)
! 
      end subroutine mgw
!-----------------------------------------------------------------------
      subroutine set_gw_params
      use gw_drag,only: gw_inti
      use gw_share,only: cpair,cpwv,grav,rair,pmid,pint,pdel,hypi,
     |  lnpint,rpdel
      use params_module,only: dz,zibot
      use params_module,only: pcols=>nlonp4, pver=>nlev
!
! Local
      integer :: k,kk
      real :: zpc

! Set pressures for gw_drag module (cam style, i.e., top->bot)
!     real :: pmid(pcols,pver)        ! midpoint pressures
!     real :: pint(pcols,0:pver)      ! interface pressures
!     real :: hypi(0:pver)            ! interface pressures
!
! real :: pmid(pcols,pver) ! midpoint pressures (gw_share.F)
      do k=1,nlev     ! 1,44
        kk = nlev-k+1 ! nlev,1,-1
        zpc = zibot+0.5*dz+float(k-1)*dz 
        pmid(:,kk)   = 5.e-5*exp(-zpc)
      enddo

!     write(6,"('set_gw_params: pmid(1,:)=',/,(6e12.4))")
!    |  pmid(1,:)

! real :: pint(pcols,0:pver) ! interface pressures (gw_share.F)
      do k=0,nlev
        kk = nlev-k   ! nlev,0,-1
        zpc = zibot+float(k)*dz 
        pint(:,kk) = 5.e-5*exp(-zpc)
      enddo

!     write(6,"('set_gw_params: pint(1,:)=',/,(6e12.4))")
!    |  pint(1,:)

      do k=1,nlev
        pdel(:,k) = pint(:,k)-pint(:,k-1)
        rpdel(:,k) = 1./pdel(:,k)
      enddo
!
! 5/4/06 btf: Check pint before taking log to avoid fpe during debug
!
!     real :: pint(pcols,0:pver+1)      ! interface pressures
!     write(6,"('set_gw_params before log(pint): pint=',/,(6e12.4))")
!    |  pint

!     lnpint = log(pint)
      do k=0,pver+1
        do kk=1,pcols
          if (pint(kk,k) > 1.e-20) then
            lnpint(kk,k) = log(pint(kk,k))
          else
            lnpint(kk,k) = 1.e-20
          endif
        enddo
      enddo

      hypi(:) = pint(1,:)
      call gw_inti(cpair,cpwv,grav,rair,hypi,nlev)
! 
      end subroutine set_gw_params
!-----------------------------------------------------------------------
      subroutine alloc_gw(lon0,lon1,lat0,lat1)
!
! Args:
      integer,intent(in) :: lon0,lon1,lat0,lat1
!
! Local:
      integer :: istat
!
! Allocate arrays needed in gravity wave parameterization:
!
      allocate(raykk(nlevp1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_gw: error allocating',
     |  ' raykk: stat=',i3)") istat
!
      allocate(gwsclh(nlevp1,lon0:lon1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_gw: error allocating',
     |  ' gwsclh: stat=',i3)") istat
!
      allocate(difk_gw(nlevp1,lon0:lon1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_gw: error allocating',
     |  ' difk_gw: stat=',i3)") istat
!
!     allocate(difkk(nlevp1,lon0:lon1),stat=istat)
!     if (istat /= 0) write(6,"('>>> alloc_gw: error allocating',
!    |  ' difkk: stat=',i3)") istat
!
!     allocate(difkt(nlevp1,lon0:lon1),stat=istat)
!     if (istat /= 0) write(6,"('>>> alloc_gw: error allocating',
!    |  ' difkk: stat=',i3)") istat
!
!     allocate(difkv(nlevp1,lon0:lon1),stat=istat)
!     if (istat /= 0) write(6,"('>>> alloc_gw: error allocating',
!    |  ' difkv: stat=',i3)") istat
!
      end subroutine alloc_gw
!-----------------------------------------------------------------------
      end module mgw_module