subroutine zm_convr(q0,t0,p,z,pf,zf,q1,t1, &
                    cape0,cape,dcape,jb,jt0,jt,mb,mb0,cape00,negcape)
implicit none

! Constants
integer,parameter::r8 = selected_real_kind(12)
!integer,parameter::nt = 1440
integer,parameter::nt = 1
integer,parameter::nz = 40

real(r8)::cp = 1004.64_r8
!real(r8)::dt = 1._r8
real(r8)::dt = 360._r8
real(r8)::g  = 9.80616_r8
real(r8)::rgas = 6.02214e26_r8*1.38065e-23_r8/28.966_r8
real(r8)::rl = 2.501e6_r8
!real(r8)::tpert(nt) = 0.5_r8
real(r8)::tpert(nt) = 0._r8 ! wx: 2019-05-12
! Variables
real(r8),intent(in)::q0(nt,nz)
real(r8),intent(in)::t0(nt,nz)
real(r8),intent(in)::q1(nt,nz)
real(r8),intent(in)::t1(nt,nz)
real(r8),intent(in)::p(nt,nz)
real(r8),intent(in)::z(nt,nz)
real(r8),intent(in)::pf(nt,nz+1)
real(r8),intent(in)::zf(nt,nz+1)


real(r8),intent(out)::cape0(nt)
real(r8),intent(out)::cape00(nt)
real(r8),intent(out)::cape(nt)
real(r8),intent(out)::dcape(nt)
real(r8),intent(out)::negcape(nt)
real(r8),intent(out)::mb(nt)
real(r8),intent(out)::mb0(nt)

integer,intent(out)::jb(nt)
integer,intent(out)::jt0(nt)
integer,intent(out)::jt(nt)

integer::jtj(nt)

real(r8) tp0(nt,nz)
real(r8) qstp0(nt,nz)
real(r8) tl0(nt)
integer jb0(nt)
integer mx0(nt)
integer lcl0(nt)
integer lon0(nt)

real(r8) q(nt,nz)
real(r8) t(nt,nz)
real(r8) tp(nt,nz)
real(r8) qstp(nt,nz)
real(r8) tl(nt)
real(r8) negcape2(nt)

integer lcl(nt)
integer lon(nt)

integer i
integer k

!!
real(r8) s(nt,nz)
real(r8) shat(nt,nz)
real(r8) qhat(nt,nz)
real(r8) hsat(nt,nz)
real(r8) dp(nt,nz)
real(r8) sdifr
real(r8) qdifr
real(r8) mumax(nt)
!real(r8) est(nt)
real(r8) qst(nt,nz)
real(r8) mu(nt,nz)
real(r8) eu(nt,nz)
real(r8) du(nt,nz)
real(r8) md(nt,nz)
real(r8) ed(nt,nz)
real(r8) sd(nt,nz)
real(r8) qd(nt,nz)
real(r8) mc(nt,nz)
real(r8) qu(nt,nz)
real(r8) su(nt,nz)
real(r8) qs(nt,nz)
real(r8) hmn(nt,nz)
real(r8) ql(nt,nz)
real(r8) cmeg(nt,nz)
integer jlcl(nt)
integer j0(nt)
integer jd(nt)
integer jt2(nt)
real(r8) pflx(nt,nz+1)
real(r8) evp(nt,nz)
real(r8) cu(nt,nz)
real(r8) rprd(nt,nz)
!!

! Routine
! Initialization
!424-431
do i = 1,nt
    dcape(i) = 0._r8
    cape(i)  = 0._r8
    cape0(i) = 0._r8
    mb(i)    = 0._r8
    mb0(i)    = 0._r8
    jb(i)    = 1
    jb0(i)   = 1
    !mx0(i)   = 1
    lcl0(i)  = 1
    lcl(i)   = 1
    jt0(i)   = 40
    jt(i)    = 40
end do    
!405-414
do k = 1,nz
    do i = 1,nt
        !s(i,k) = t(i,k)+(g/cp)*z(i,k)
        s(i,k) = t0(i,k)+(g/cp)*z(i,k)
        tp(i,k) = 0.0_r8
        shat(i,k) = s(i,k)
        !qhat(i,k) = q(i,k)
        qhat(i,k) = q0(i,k)
    end do
end do    
! cape0
call buoyan_dilute(1,nt, &
                   q0,t0,p,z,pf, &
                   tp0,qstp0,tl0,rl,cape0, &
                   lcl0,jt0,lon0,jb0, &
                   rgas,g,cp,0, &
                   tpert,negcape,mx0,cape00)
do i = 1,nt
    mx0(i) = jb0(i)
    do k = 1,nz
        if(k<mx0(i)) then
            q(i,k) = q1(i,k)
            if(q(i,k).le.0._r8) q(i,k) = 1.e-13_r8
            t(i,k) = t1(i,k)
        end if
        if(k>=mx0(i)) then
            q(i,k) = q0(i,k)
            t(i,k) = t0(i,k)
        end if
    end do
end do
!!
! cape
call buoyan_dilute(2,nt, &
                   q,t,p,z,pf, &
                   tp,qstp,tl,rl,cape, &
                   lcl,jt,lon,jb, &
                   rgas,g,cp,0, &
                   tpert,negcape2,mx0)
! dcape
do i = 1,nt
    dcape(i) = (cape(i)-cape0(i))/dt
end do

return
end subroutine zm_convr
