subroutine fromztop
use modglobal, only : kmax,kb,ke,kh,dzf,dzh,rv,rd,cp,tmelt,zf,grav,pref0,lEB
use modfields, only : qt0av,ql0av,presf,presh,thvh,thvf,IIcs
use modsurfdata,only : ps,thvs
implicit none
integer k
real rdocp
real,allocatable,dimension (:) :: thetah, qth, qlh
allocate(thetah(kb:ke+kh), qth(kb:ke+kh), qlh(kb:ke+kh))
rdocp = rd/cp
!**************************************************
! 1.0 Determine theta and qt at half levels *
!**************************************************
do k=kb+1,ke+kh
thetah(k) = (th0av(k)*dzf(k-1) + th0av(k-1)*dzf(k))/(2*dzh(k))
qth (k) = (qt0av(k)*dzf(k-1) + qt0av(k-1)*dzf(k))/(2*dzh(k))
qlh (k) = (ql0av(k)*dzf(k-1) + ql0av(k-1)*dzf(k))/(2*dzh(k))
end do
!**************************************************
! 2.1 calculate pressures at full levels *
! assuming hydrostatic equilibrium *
!**************************************************
! 1: lowest level: use thvs
thvh(kb) = thvs
presf(kb) = ps**rdocp - &
grav*(pref0**rdocp)*zf(kb) /(cp*thvh(kb))
presf(kb) = presf(kb)**(1./rdocp)
! 2: higher levels
do k=kb+1,ke+kh
thvh(k) = thetah(k)*(1+(rv/rd-1)*qth(k)-rv/rd*qlh(k))
presf(k) = presf(k-1)**rdocp - grav*(pref0**rdocp)*dzh(k) /(cp*thvh(k))
presf(k) = presf(k)**(1./rdocp)
end do
!**************************************************
! 2.2 calculate pressures at half levels *
! assuming hydrostatic equilibrium *
!**************************************************
presh(kb) = ps
thvf(kb) = th0av(kb)*(1+(rv/rd-1)*qt0av(kb)-rv/rd*ql0av(kb))
do k=kb+1,ke+kh
thvf(k) = th0av(k)*(1+(rv/rd-1)*qt0av(k)-rv/rd*ql0av(k))
presh(k) = presh(k-1)**rdocp - grav*(pref0**rdocp)*dzf(k-1) / (cp*thvf(k-1))
presh(k) = presh(k)**(1./rdocp)
end do
deallocate(thetah, qth, qlh)
return
end subroutine fromztop