module mgw
!
! 1/27/06 btf: Introduce some terms from timegcm mgw.F for the new
!              zp-10 lbc.
      use params_module,only: nlev,nlevp1
      implicit none
!
      real ::
     |  rayk(nlev),          ! rayleigh friction
     |  difk(nlevp1),
     |  dift(nlevp1)
!
! Arrays on subdomains (allocated by alloc_gw)
      real,allocatable,dimension(:,:,:) ::
     |  raykk                ! rayleigh friction (lev,lon,lat) subdomains
      real,allocatable,dimension(:,:,:) ::
     |  difkk                ! eddy viscosity (lev,lon,lat) subdomains
!
      contains
!-----------------------------------------------------------------------
      subroutine backgrnd
      use params_module,only: zibot,dz
      implicit none

      real,parameter ::
     |  eddylb  = 5.0e-7,
     |  zplb    = -10.,
     |  eddymb1 = 6.0e-7,
     |  zpmb1   = -9.,
     |  eddymb2 = 9.e-7,
     |  zpmb2   = -8.,
     |  eddytp  = 3.0e-6,
     |  zptp    = -6.0,
     |  eddtlb  = 5.0e-7,
     |  eddtmb1 = 6.0e-7,
     |  eddtmb2 = 9.0e-7,
     |  eddttp  = 5.0e-7  ! timegcm1.2 value is 5.0e-7
!
      real,parameter ::
     |  rayklb  = 2.5e-5,
     |  raykmb1 = 1.0e-5,
     |  raykmb2 = 0.7e-5,
     |  rayktb  = 0.2e-5
!
!    |  rayklb  = 2.5e-5,
!    |  raykmb1 = 2.5e-5,
!    |  raykmb2 = 1.5e-5,
!    |  rayktb  = 0.8e-5
!
      integer :: k,klb,kmb1,kmb2,kturb
!
      real ::
     |  col1(nlevp1)
!
      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
      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))
        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))
        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))
        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))
      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,"('mgw backgrnd: difk=',/,(6e12.4))") difk
!     write(6,"('mgw backgrnd: dift=',/,(6e12.4))") dift
!     write(6,"('mgw backgrnd: rayk=',/,(6e12.4))") rayk

      end subroutine backgrnd
!-----------------------------------------------------------------------
      subroutine gwsource(tn,barm,rfricu,rfricv,eddy,lev0,lev1,
     |  lon0,lon1,lat)
      use addfld_module,only: addfld
      use input_module,only: gswm_rayfric,sgar_eddy
      use cons_module,only: grav,gask
!
! Args:
      real,dimension(lev0:lev1,lon0-2:lon1+2),
     |  intent(in) :: 
     |  tn,    ! neutral temperature on gcm subdomain (K)
     |  barm   ! mean mass
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0:lon1),intent(in) :: 
     |  rfricu,rfricv   ! rayleigh friction (s-1) from gswm_mlt.F
      real,dimension(lev0:lev1,lon0:lon1),intent(inout) :: 
     |  eddy            ! eddy diffusion (m^2/s) from sol-gar (gswm_mlt.F)
!
! Local:
      integer :: i,k
      real :: scht(lev0:lev1)
!
! Rayleigh friction (note RAYKK is saved to secondary hist in duv.F)
      do i=lon0,lon1
        raykk(:,i,lat) = rayk(:)
!
! 2/28/06 btf: add rayleigh friction from gswm_mlt.F:
! (rfricu and rfricv are identical -- just add the U component for now)
! See duv.F for raykk diagnostic in m/s/day.
!
        if (gswm_rayfric > 0) 
     |    raykk(:,i,lat) = raykk(:,i,lat)+rfricu(:,i)
      enddo
!
! Eddy viscosity:
      do i=lon0,lon1
        difkk(:,i,lat) = difk(:)
!
! Add sol-gar eddy diffusion from gswm_mlt.F:
! Use scale height to convert sgar eddy from m^2/s to 1/s
!
        if (sgar_eddy > 0) then ! eddy is from gswm_mlt
          do k=lev0+1,lev1-1
            scht(k) = (gask*0.5*(tn(k-1,i)+tn(k,i))/
     |        (grav*barm(k,i)))*1.e-2  ! meters
          enddo
          scht(lev0) = 2.*scht(lev0+1)-scht(lev0+2)
          scht(lev1) = 1.5*scht(lev1-1)-0.5*scht(lev1-2)
          difkk(:,i,lat) = difkk(:,i,lat)+eddy(:,i)/scht(:)**2
        endif ! sgar_eddy > 0
      enddo ! i=lon0,lon1
!
! Save sgar eddy (s-1) to secondary history:
      call addfld('DIFKK','Eddy viscosity',
     |  's-1',difkk(:,:,lat),'lev',lev0,lev1,'lon',lon0,lon1,lat)
      end subroutine gwsource
!-----------------------------------------------------------------------
      subroutine alloc_mgw(nlev,lon0,lon1,lat0,lat1)
!
! Args:
      integer,intent(in) :: nlev,lon0,lon1,lat0,lat1
!
! Local:
      integer :: istat
!
! Allocate arrays needed in gravity wave parameterization:
!
      allocate(raykk(nlev,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_mgw: error allocating',
     |  ' raykk: stat=',i3)") istat
!
      allocate(difkk(nlev,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_mgw: error allocating',
     |  ' difkk: stat=',i3)") istat
!
!     allocate(difkt(nlev,lon0:lon1),stat=istat)
!     if (istat /= 0) write(6,"('>>> alloc_mgw: error allocating',
!    |  ' difkk: stat=',i3)") istat
!
!     allocate(difkv(nlev,lon0:lon1),stat=istat)
!     if (istat /= 0) write(6,"('>>> alloc_mgw: error allocating',
!    |  ' difkv: stat=',i3)") istat
!
      end subroutine alloc_mgw
!-----------------------------------------------------------------------
      end module mgw