subroutine zm_convr(q0,t0,p,z,pf,zf,ptendq,ptendt, &
                    cape0,cape,dcape,jb,jt0,jt,mb,dadt,jt3,mb0,cape00)
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)::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(in)::ptendq(nt,nz)
real(r8),intent(in)::ptendt(nt,nz)


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)::mb(nt)
real(r8),intent(out)::mb0(nt)

integer,intent(out)::jb(nt)
integer,intent(out)::jt0(nt)
integer,intent(out)::jt(nt)
integer,intent(out)::jt3(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) negcape(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)
real(r8),intent(out)::dadt(nt)
!!

! 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)
!!
!wx: 20191227
do i = 1,nt
    mx0(i) = jb0(i)
    jb(i)  = mx0(i)
    !do k = 1,nz-1
    !    if(p(i,k+1).ge.100.and.p(i,k).lt.100) then
    !        jt(i) = k+1
    !    end if    
    !end do    
    jt(i)  = jt0(i)
    lcl(i) = lcl0(i)
    tl(i)  = tl0(i)
    do k = 1,nz
        q(i,k) = q0(i,k)
        if(q(i,k).le.0_r8) q(i,k) = 1.e-13_r8
        t(i,k) = t0(i,k)
        tp(i,k) = tp0(i,k)
        qstp(i,k) = qstp0(i,k)
    end do 
end do    

!do i = 1,nt
!    mx0(i) = jb0(i)
!    do k = 1,nz
!        if(k<mx0(i)) then
!            q(i,k) = q0(i,k)+ptendq(i,k)*dt
!            if(q(i,k).le.0._r8) q(i,k) = 1.e-13_r8
!            t(i,k) = t0(i,k)+ptendt(i,k)*dt
!        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,negcape,mx0)
! dcape
!do i = 1,nt
!    dcape(i) = (cape(i)-cape0(i))/dt
!end do

!call dcape_mth2(q,t,p,tp,qstp,zf,mx0,lcl,jt,ptendq,ptendt, &
!                rgas,cp,g,dcape)

!405-414
!do k = 1,nz
!    do i = 1,nt
!        s(i,k) = t(i,k)+(g/cp)*z(i,k)
!        shat(i,k) = s(i,k)
!        qhat(i,k) = q(i,k)
!    end do
!end do    
!506
do k = 1,nz
    do i = 1,nt
        dp(i,k) = pf(i,k+1)-pf(i,k)
    end do
end do    
!555-575
do k = 2,nz
    do i = 1,nt
        sdifr = 0._r8
        qdifr = 0._r8
        if(s(i,k)>0._r8.or.s(i,k-1)>0._r8) &
            sdifr = abs((s(i,k)-s(i,k-1))/max(s(i,k-1),s(i,k)))
        if(q(i,k)>0._r8.or.q(i,k-1)>0._r8) &
            qdifr = abs((q(i,k)-q(i,k-1))/max(q(i,k-1),q(i,k)))
        if(sdifr>1.E-6_r8) then
            shat(i,k) = log(s(i,k-1)/s(i,k))*s(i,k-1)*s(i,k)/(s(i,k-1)-s(i,k))
        else
            shat(i,k) = 0.5_r8*(s(i,k)+s(i,k-1))
        end if
        if(qdifr>1.E-6_r8) then
            qhat(i,k) = log(q(i,k-1)/q(i,k))*q(i,k-1)*q(i,k)/(q(i,k-1)-q(i,k))
        else
            qhat(i,k) = 0.5_r8*(q(i,k)+q(i,k-1))
        end if
    end do
end do    

!est,qst
!do k = 1,nz
!    do i = 1,nt
!        call qmmr_hPa(t(i,k),p(i,k),est(i),qst(i))
!    end do
!end do    

!593-603
call cldprp_wx(1,&
               q,t,p, &
               z,s,mu,eu,du, &
               md,ed,sd,qd,mc, &
               qu,su,zf,qst,hmn, &
               hsat,shat,ql, &
               cmeg,mx0,jt,jt2,jlcl, &
               mx0,j0,jd,rl,nt, &
               rgas,g,cp,0, &
               pflx,evp,cu,rprd,1,1)
!!wx: 20200115
do i = 1,nt
    jt3(i) = jt2(i)
    jtj(i) = 1 ! Deep convection
    if((p(i,jb(i))-p(i,jt2(i))).le.200._r8.and.(p(i,jb(i))-p(i,jt2(i))).gt.0._r8) then
        jtj(i) = 2 ! Shallow convection
    end if
end do    
!!
!607-617
do k = 1,nz
    do i = 1,nt
         du(i,k) = du(i,k)*(zf(i,k)-zf(i,k+1))/dp(i,k)
    end do
end do    
!
do i = 1,nt
    do k = 1,nz-1
        if(p(i,k+1).ge.150.and.p(i,k).lt.150) then
            jt(i) = k+1
            jt2(i) = jt(i)
        end if    
    end do
end do

call dcape_mth2(q,t,p,tp,qstp,zf,mx0,lcl,jt,ptendq,ptendt, &
                rgas,cp,g,dcape)
!
call dadt_mth2(q,t,p,s, &
               tp,qu,su,mc, &
               du,mu,md,qd,sd, &
               qhat,shat,dp,qstp,zf, &
               ql,tl,lcl,jt,jt2,mx0,dadt)

do i = 1,nt
    if((dcape(i))>0._r8.and.dadt(i)<0._r8.and.p(i,mx0(i)).ge.600._r8) then
    !if((dcape(i)+dcape0(i))>0._r8.and.dadt(i)<0._r8.and.p(i,mx0(i)).ge.600._r8.and.jtj(i).eq.1) then
        mb(i) = max(-(dcape(i))/dadt(i),0._r8)*100._r8/g
        mb0(i) = max(-(dcape(i))/dadt(i),0._r8)*100._r8/g
    else
        mb(i) = 0._r8
        mb0(i) = 0._r8
    end if
    if(jtj(i).eq.2) then
        mb(i) = -1._r8
    end if    
end do
!!
!do i = 1,nt
!    dcape(i) = dcape(i)+dcape0(i)
!end do
!!
!
!do i = 1,nt
!    mumax(i) = 0
!end do
!do k = 2,nz
!    do i = 1,nt
!        mumax(i) = max(mumax(i),mu(i,k)/dp(i,k))
!    end do
!end do
!
!do i = 1,nt
!    if(mumax(i)>0._r8) then
!        mb(i) = min(mb(i),0.5_r8/(dt*mumax(i)))
!    else
!        mb(i) = 0._r8
!    end if
!end do
! mb/s => kg/m^2/s
!do k = 1,nz
!    do i = 1,nt
!        mu(i,k) = mu(i,k)*mb(i)*100._r8/g
!    end do
!end do    
!
!do i = 1,nt
!    mb(i) = mu(i,mx0(i))
!end do

!print*,"dcape =",dcape
!print*,"dadt=",dadt
!print*,"mb=",mb
return
end subroutine zm_convr
