calthv Subroutine

  subroutine calthv
    use modglobal, only : lmoist,ib,ie,jb,je,kb,ke,kh,zf,zh,dzh,rlv,rd,rv,cp,eps1
    use modfields, only : thl0,thl0h,ql0,ql0h,qt0,qt0h,sv0,exnf,exnh,thv0h,dthvdz
    use modsurfdata,only : dthldz,dqtdz

    implicit none

    integer i, j, k
    real    qs
    real    c1,c2,dq,dth,dthv,temp
    real    a_dry, b_dry, a_moist, b_moist, c_liquid, epsilon, eps_I, chi_sat, chi
    real    del_thv_sat, del_thv_dry
    dthvdz = 0.0
    if (lmoist) then

       do  k=kb,ke+kh
          do  j=jb,je
             do  i=ib,ie
                thv0h(i,j,k) = (thl0h(i,j,k)+rlv*ql0h(i,j,k)/(cp*exnh(k))) &
                     *(1+(rv/rd-1)*qt0h(i,j,k)-rv/rd*ql0h(i,j,k))
             end do
          end do
       end do

       do k=kb+1,ke
          do j=jb,je
             do i=ib,ie
                !
                !         default thv jump computed unsaturated
                !
                epsilon = rd/rv
                eps_I = 1/epsilon - 1.  !cstep approx 0.608

                a_dry = 1. + eps_I * qt0(i,j,k)
                b_dry = eps_I * thl0(i,j,k)

                dth = thl0(i,j,k+1)-thl0(i,j,k-1)
                dq  = qt0(i,j,k+1)-qt0(i,j,k-1)

                del_thv_dry = a_dry   * dth + b_dry * dq

                dthv = del_thv_dry

                if  (ql0(i,j,k)> 0) then  !include moist thermodynamics
                   temp = thl0(i,j,k)*exnf(k)+(rlv/cp)*ql0(i,j,k)
                   qs   = qt0(i,j,k) - ql0(i,j,k)

                   a_moist = (1.-qt0(i,j,k)+qs/epsilon*(1.+rlv/(rv*temp))) &
                        /(1.+rlv**2*qs/(cp*rv*temp**2))
                   b_moist = a_moist*rlv/cp-temp
                   c_liquid = a_dry * rlv / cp - thl0(i,j,k) / epsilon

                   del_thv_sat = a_moist * dth + b_moist * dq

                   chi     = 2*chi_half*(zf(k) - zf(k-1))/(dzh(k)+dzh(k+1))
                   chi_sat = c_liquid * ql0(i,j,k) / (del_thv_dry - del_thv_sat)

                   if (chi < chi_sat) then  !mixed parcel is saturated
                      dthv = del_thv_sat
                   end if
                end if

                dthvdz(i,j,k) = dthv/(dzh(k+1)+dzh(k))
             end do
          end do
       end do

       do j=jb,je
          do i=ib,ie
            dthvdz(i,j,kb)=0.
          end do
       end do

    else
       !      thv0h = thl0h
       thv0h = thl0h(:,:,kb:ke+kh)
       do k=kb+1,ke
          do j=jb,je
             do i=ib,ie
                dthvdz(i,j,k) = (thl0(i,j,k+1)-thl0(i,j,k-1))/(dzh(k+1)+dzh(k))
             end do
          end do
       end do
       do  j=jb,je
          do  i=ib,ie
            dthvdz(i,j,kb)=0.
          end do
       end do
    end if

    !CvH remove WHERE
    !where (abs(dthvdz)<eps1)
    !  dthvdz = sign(eps1,dthvdz)
    !end where
    do k=kb,ke
       do j=jb,je
          do i=ib,ie
             if(abs(dthvdz(i,j,k)) < eps1) then
                dthvdz(i,j,k) = sign(eps1, dthvdz(i,j,k))
             end if
          end do
       end do
    end do



  end subroutine calthv