thermodynamics Subroutine

  subroutine thermodynamics
    use modglobal, only : lmoist, timee, kb, ke, kh, ib, ih, ie, jb, jh, je,rlv, cp, rslabs, rd, rv, libm, eps1
    use modfields, only : thl0,thl0h,qt0,qt0h,ql0,ql0h,presf,presh,exnf,exnh,thvh,thv0h,qt0av,ql0av,thvf,rhof,IIc,IIw,IIcs,IIws
    use modmpi,    only : slabsum,avexy_ibm,myid
!ILS13 added variables behind "exnh"
    implicit none
    integer :: k
    if (timee==0) call diagfld
    if (lmoist) then
       call thermo(thl0,qt0,ql0,presf,exnf)
    end if

    call diagfld
    call calc_halflev !calculate halflevel values of qt0 and thl0
    if (lmoist) then
       call thermo(thl0h,qt0h,ql0h,presh,exnh)
    end if
    call calthv

!ILS13 introduced from DALES4.0   13.05.2015
    thvh=0.
!    call slabsum(thvh,kb,ke+kh,thv0h(:,:,kb:ke+kh),ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh) !redefine halflevel thv using calculated thv
!    thvh = thvh/rslabs
    call avexy_ibm(thvh(kb:ke+kh),thv0h(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIw(ib:ie,jb:je,kb:ke+kh),IIws(kb:ke+kh),.false.)

!    if (libm) then
!      call avexy_ibm(thvh(kb:ke),thv0h(ib:ie,jb:je,kb:ke),ib,ie,jb,je,kb,ke,IIw(ib:ie,jb:je,kb:ke),IIws(kb:ke))    
!    else
!      call slabsum(thvh,kb,ke+kh,thv0h(:,:,kb:ke+kh),ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
!     !redefine halflevel thv using calculated thv
!     thvh = thvh/rslabs
!    end if

    thvh(kb) = th0av(kb)*(1+(rv/rd-1)*qt0av(kb)-rv/rd*ql0av(kb)) ! override first level
    if (abs(thvh(kb+1))<eps1) then
      thvh(kb+1) = th0av(kb+1)*(1+(rv/rd-1)*qt0av(kb+1)-rv/rd*ql0av(kb+1)) ! override second level if all blocks at kb
    end if
!    where (thvh==0) !override slabs completely covered by blocks
!      thvh = th0av(kb)*(1+(rv/rd-1)*qt0av(kb)-rv/rd*ql0av(kb))
!    endwhere

    do k=kb,ke+kh
!    thv0(ib+ih:ie,jb+jh:je,k) = (thl0(ib+ih:ie,jb+ih:je,k)+rlv*ql0(ib+ih:ie,jb+ih:je,k)/(cp*exnf(k)))*(1+(rv/rd-1)*qt0(ib+ih:ie,jb+ih:je,k)-rv/rd*ql0(ib+ih:ie,jb+ih:je,k))
    thv0(ib:ie,jb:je,k) = (thl0(ib:ie,jb:je,k)+rlv*ql0(ib:ie,jb:je,k)/(cp*exnf(k)))*(1+(rv/rd-1)*qt0(ib:ie,jb:je,k)-rv/rd*ql0(ib:ie,jb:je,k))
    enddo
    thvf = 0.0

    !write(*,*) "thv0",thv0
!    call slabsum(thvf,kb,ke+kh,thv0,ib,ie+ih,jb,je+jh,kb,ke+kh,ib+ih,ie,jb+ih,je,kb,ke+kh)
!    call slabsum(thvf,kb,ke+kh,thv0,ib,ie,jb,je,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
    call avexy_ibm(thvf(kb:ke+kh),thv0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc(ib:ie,jb:je,kb:ke+kh),IIcs(kb:ke+kh),.false.)
!    write(*,*) 'IIc(2,2,:), myid' , IIc(12,2,:), myid

!    where (thvf==0) !override slabs completely covered by blocks
!      thvf = th0av(kb)*(1+(rv/rd-1)*qt0av(kb)-rv/rd*ql0av(kb))
!    endwhere

!    thvf = thvf/rslabs
    !write(*,*) "thvf",thvf
    !write(*,*) "exnf",exnf
   
!    do k=1,k1
!      rhof(k) = presf(k)/(rd*thvf(k)*exnf(k))
!    end do

  end subroutine thermodynamics