subroutine 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,mx,dadt)

implicit none

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

real,parameter::rd = 6.02214e26_r8*1.38065e-23_r8/28.966_r8
real,parameter::rl = 2.501e6_r8
real,parameter::cp = 1004.64_r8
real,parameter::grav = 9.80616_r8
real,parameter::eps1 = 18.016_r8/28.966_r8

real(r8),intent(in)::p(nt,nz)

real(r8),intent(in)::q(nt,nz)
real(r8),intent(in)::qd(nt,nz)
real(r8),intent(in)::ql(nt,nz)
real(r8),intent(in)::qu(nt,nz)
real(r8),intent(in)::qhat(nt,nz)
real(r8),intent(in)::qstp(nt,nz)

real(r8),intent(in)::t(nt,nz)
real(r8),intent(in)::tp(nt,nz)
real(r8),intent(in)::tl(nt)

real(r8),intent(in)::s(nt,nz)
real(r8),intent(in)::sd(nt,nz)
real(r8),intent(in)::su(nt,nz)
real(r8),intent(in)::shat(nt,nz)

real(r8),intent(in)::mc(nt,nz)
real(r8),intent(in)::mu(nt,nz)
real(r8),intent(in)::md(nt,nz)

real(r8),intent(in)::dp(nt,nz)
real(r8),intent(in)::du(nt,nz)
real(r8),intent(in)::zf(nt,nz+1)

integer,intent(in)::lcl(nt)
integer,intent(in)::jt(nt)
integer,intent(in)::jt2(nt)
integer,intent(in)::mx(nt)

real(r8),intent(out)::dadt(nt)

real(r8) dtpdt(nt,nz)
real(r8) dqsdtp(nt,nz)
real(r8) dtmdt(nt,nz)
real(r8) dqmdt(nt,nz)
real(r8) dboydt(nt,nz)
real(r8) thetavp(nt,nz)
real(r8) thetavm(nt,nz)
real(r8) dtbdt(nt),dqbdt(nt),dtldt(nt)
real(r8) beta
real(r8) eb
real(r8) debdt

integer i,k,kmin,kmax

!114-128
do i = 1,nt
     eb = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i)))
     dtbdt(i) = 0._r8
     dqbdt(i) = 0._r8
     debdt = eps1*p(i,mx(i))/ (eps1+q(i,mx(i)))**2*dqbdt(i)
     dtldt(i) = -2840._r8* (3.5_r8/t(i,mx(i))*dtbdt(i)-debdt/eb)/ &
                (3.5_r8*log(t(i,mx(i)))-log(eb)-4.805_r8)**2
end do                
!132-138
do k = 1,nz
    do i = 1,nt
        dtmdt(i,k) = 0._r8
        dqmdt(i,k) = 0._r8
    end do
end do    
!140-149
do k = 1,nz-1
    do i = 1,nt
        if(k==jt2(i)) then
            dtmdt(i,k) = (1._r8/dp(i,k))*(mu(i,k+1)* (su(i,k+1)-shat(i,k+1)- &
                          rl/cp*ql(i,k+1))+md(i,k+1)* (sd(i,k+1)-shat(i,k+1)))
            dqmdt(i,k) = (1._r8/dp(i,k))*(mu(i,k+1)* (qu(i,k+1)- &
                          qhat(i,k+1)+ql(i,k+1))+md(i,k+1)*(qd(i,k+1)-qhat(i,k+1)))
        end if
    end do
end do    
!151-169
beta = 0._r8
do k = 1,nz-1
    do i = 1,nt
        if(k>jt2(i).and.k<mx(i)) then
            dtmdt(i,k) = (mc(i,k)* (shat(i,k)-s(i,k))+mc(i,k+1)* (s(i,k)-shat(i,k+1)))/ &
                          dp(i,k) - rl/cp*du(i,k)*(beta*ql(i,k)+ (1-beta)*ql(i,k+1))
            dqmdt(i,k) = (mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)+cp/rl* (su(i,k+1)-s(i,k)))- &
                          mu(i,k)* (qu(i,k)-qhat(i,k)+cp/rl*(su(i,k)-s(i,k)))+md(i,k+1)* &
                         (qd(i,k+1)-qhat(i,k+1)+cp/rl*(sd(i,k+1)-s(i,k)))-md(i,k)* &
                         (qd(i,k)-qhat(i,k)+cp/rl*(sd(i,k)-s(i,k))))/dp(i,k) + &
                          du(i,k)* (beta*ql(i,k)+(1-beta)*ql(i,k+1))
        end if
    end do
end do    
!171-192
do k = 1,nz
    do i = 1,nt
        if(k>=jt(i).and.k<=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))
            dqsdtp(i,k) = qstp(i,k)* (1._r8+qstp(i,k)/eps1)*eps1*rl/(rd*tp(i,k)**2)

            dtpdt(i,k) = tp(i,k)/ (1._r8+rl/cp* (dqsdtp(i,k)-qstp(i,k)/tp(i,k)))* &
                         (dtbdt(i)/t(i,mx(i))+rl/cp* (dqbdt(i)/tl(i)-q(i,mx(i))/ &
                          tl(i)**2*dtldt(i)))

            dboydt(i,k) = ((dtpdt(i,k)/tp(i,k)+1._r8/(1._r8+1.608_r8*qstp(i,k)-q(i,mx(i)))* &
                          (1.608_r8 * dqsdtp(i,k) * dtpdt(i,k) -dqbdt(i))) - (dtmdt(i,k)/t(i,k)+0.608_r8/ &
                          (1._r8+0.608_r8*q(i,k))*dqmdt(i,k)))*grav*thetavp(i,k)/thetavm(i,k)
        end if
    end do
end do    
!194-207
do k = 1,nz
    do i = 1,nt
        if(k>lcl(i).and.k<mx(i)) then
            thetavp(i,k) = tp(i,k)*(1000._r8/p(i,k)**(rd/cp)*(1._r8+0.608_r8*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)))

            dboydt(i,k) = (dtbdt(i)/t(i,mx(i))+0.608_r8/(1._r8+0.608_r8*q(i,mx(i)))*dqbdt(i)- &
                           dtmdt(i,k)/t(i,k)-0.608_r8/(1._r8+0.608_r8*q(i,k))*dqmdt(i,k))* &
                           grav*thetavp(i,k)/thetavm(i,k)
        end if
    end do
end do    
!212-221
dadt = 0._r8
kmin = minval(jt)
kmax = maxval(mx)-1
do k = kmin,kmax
    do i = 1,nt
        if(k>=jt(i).and.k<=mx(i)-1) then
            dadt(i) = dadt(i)+dboydt(i,k)*(zf(i,k)-zf(i,k+1))
        end if
    end do
end do    

end subroutine dadt_mth2
