subroutine dcape_mth2(q, t, p, tp, qstp, zf, mx, lcl, lel, ptendq, ptendt, &
                      rd, cp, grav, &
                      dcape)

implicit none

integer, parameter :: r8    = selected_real_kind(12)
!integer, parameter :: pcols = 1440
integer, parameter :: pcols = 1
integer, parameter :: pver  = 40

real(r8), intent(inout) :: q(pcols,pver) ! spec humidity
real(r8), intent(inout) :: t(pcols,pver) ! temperature
real(r8), intent(inout) :: p(pcols,pver) ! pressure (mb)

real(r8), intent(in) :: tp(pcols,pver)   ! parcel temp
real(r8), intent(in) :: zf(pcols,pver+1) ! height of interface levels
real(r8), intent(in) :: qstp(pcols,pver) ! spec humidity of parcel

real(r8), intent(in) :: ptendq(pcols,pver) ! q tendency
real(r8), intent(in) :: ptendt(pcols,pver) ! t tendency

integer, intent(in) :: mx(pcols)    ! base of updraft
integer, intent(in) :: lcl(pcols)   ! index of lcl
integer, intent(in) :: lel(pcols)   ! index of launch level

real(r8) rd
real(r8) cp
real(r8) grav

integer i
integer k

real(r8) thetavp(pcols,pver)
real(r8) thetavm(pcols,pver)
real(r8) dbdtls(pcols,pver)

real(r8), intent(out) :: dcape(pcols)

!
do k = 1,pver
    do i = 1,pcols
        dbdtls(i,k) = 0._r8
    end do
end do
!
do k = 1,pver
    do i = 1,pcols
        if(k.ge.lel(i) .and. k.le.lcl(i)) then
            thetavp(i,k) = tp(i,k)*(1000._r8/p(i,k))**(rd/cp)* &
                           (1._r8+1.608_r8*qstp(i,k)-q(i,mx(i)))
            thetavm(i,k) = t(i,k)*(1000._r8/p(i,k))**(rd/cp)* &
                           (1._r8+0.608_r8*q(i,k))
            dbdtls(i,k)  = -(1._r8/t(i,k)*ptendt(i,k)+ &
                           0.608_r8/(1._r8+0.608_r8*q(i,k))*ptendq(i,k)) &
                           *grav*thetavp(i,k)/thetavm(i,k)
        end if
    end do
end do    
!
do i = 1,pcols
    dcape(i) = 0._r8
    do k = lel(i),mx(i)-1
        dcape(i) = dcape(i)+dbdtls(i,k)*(zf(i,k)-zf(i,k+1))
    end do
end do    

return
end subroutine dcape_mth2
