c--------------------------------------------------------------------
      subroutine dslcalc (tmk,ght,ter,t2,arr
     $     ,miy,mjx,mkzh,rheight,maxdtdz)
c
c     David Ovens addition of Deep Stable Layer
c     Defined in June 2014, after Wolyn and McKee MWR March 1989
c     values are percentage of 1.5-km depth that are stable
c     based on their criterion of lapse rate <= 2.5C/km
c
c     tmk  = temperature 3-d array on mass/half-eta levels
c     ght  = geopotential height 3-d array in m on mass/half-eta levels
c     t2   = 2m temperature 2-d array
c     rheight = depth of layer in meters (1.5km is the default)
c     maxdtdz = maximum lapse rate 2.5C/1000m = .0025 
c     NOTE: will use tbot-ttop and hgttop-hgtbot to keep positive dtdz

      real rheight, maxdtdz, dtdz, deltat, deltaz
      integer id,jd
      dimension tmk(miy,mjx,mkzh), ght(miy,mjx,mkzh)
     $     , ter(miy,mjx), t2(miy,mjx)
     $     , arr(miy,mjx)
c
      include 'comconst'
c
c     comment out next line
c      print *,'ok to here in dslcalc'
      id = 120
      jd = 300
c      id = 172
c      jd = 100
      do j = 1, mjx-1
      do i = 1, miy-1

         do k = 1, mkzh         ! descending from top of model atmos to the sfc
            km1=max(1,k-1)
            rhtk = ght(i,j,k)-ter(i,j)
            rhtkm1 = ght(i,j,km1)-ter(i,j)
            if (rhtk.gt.rheight) then
               continue
            else
               if( rhtk.le.rheight .and. rhtkm1.gt.rheight )
     $              then
c     first height that is below the top level (1.5km)
c     do a linear in height vertical interpolation of the temperature
c     y' = y(x1) + (y(x2)-y(x1)) * (x'-x1)/(x2-x1)
                  tmkt = tmk(i,j,km1) + (tmk(i,j,k)-tmk(i,j,km1)) *
     &                 (rheight-rhtkm1) /
     &                 (rhtk-rhtkm1)
                  deltat = tmk(i,j,k) - tmkt
                  deltaz = rheight - (ght(i,j,k) - ter(i,j))
               elseif (rhtk.le.rheight) then
                  if (k+1.gt.mkzh) then
                     deltat = t2(i,j) - tmk(i,j,mkzh)
                     deltaz = ght(i,j,mkzh) - ter(i,j) - 2.
                  else
                     deltat = tmk(i,j,k+1) - tmk(i,j,k)
                     deltaz = ght(i,j,k) - ght(i,j,k+1)
                  endif
               endif
               dtdz = deltat / deltaz
               if (dtdz.le.maxdtdz) then
                  arr(i,j) = arr(i,j) + 100. * deltaz / rheight
               endif
c     comment out next line
c               if (i.eq.id .and. j.eq.jd) then
c                  print '("i,j,k,deltat,deltaz,dtdz,dslcalc="
c     $                 ,3(i3,1x),5(f10.3,1x))'
c     $                 ,id,jd,k,deltat,dtdz,deltaz,arr(i,j)
c               endif
           endif
         enddo
c         if (i.eq.id .and. j.eq.jd) then
c            print '("out  ,deltat,dtdz,deltaz,dslcalc="
c     $           ,3(i3,1x),5(f10.3,1x))'
c     $           ,id,jd,k,deltat,dtdz,deltaz,arr(i,j)
c         endif
c
      enddo
      enddo
c
      return
      end