subroutine thermo (thl,qt,ql,pressure,exner)
! use modglobal, only : ih,jh,i1,j1,k1,es0,at,bt,rd,rv,rlv,cp,tmelt
use modglobal, only : ih,jh,ib,ie,jb,je,kb,ke,kh,es0,at,bt,rd,rv,rlv,cp,tmelt
use modsurfdata, only : thls
implicit none
integer i, j, k
real tl, es, qs, qsl, b1
! real, intent(in) :: qt(ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh),thl(ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh),exner(kb:ke+kh),pressure(kb:ke+kh)
real, intent(in) :: qt(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh),thl(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh),exner(kb:ke+kh),pressure(kb:ke+kh)
real, intent(out) :: ql(ib-ih:ie+ih,jb-jh:je+jh,kb:ke+kh)
real :: Tnr,qsatur,Tnr_old
integer :: niter,nitert
if (lqlnr) then
!mc calculation of T with Newton-Raphson method
!mc first guess is Tnr=tl
!mc
nitert = 0
do k=kb,ke+kh
do j=jb,je
do i=ib,ie
tl = thl(i,j,k)*exner(k)
Tnr=tl
Tnr_old=0.
do while (abs(Tnr-Tnr_old)/Tnr>1e-5)
niter = niter+1
Tnr_old = Tnr
es = es0*exp(at*(Tnr-tmelt)/(Tnr-bt))
qsatur= rd/rv*es/(pressure(k)-(1-rd/rv)*es)
Tnr = Tnr - (Tnr+(rlv/cp)*qsatur-tl- &
(rlv/cp)*qt(i,j,k))/(1+(rlv**2*qsatur)/ &
(rv*cp*Tnr**2))
end do
nitert =max(nitert,niter)
niter = 0
ql(i,j,k) = dim(qt(i,j,k)-qsatur,0.)
end do
end do
end do
else
do k=kb,ke+kh
do j=jb,je
do i=ib,ie
tl = thl(i,j,k)*exner(k)
es = es0*exp(at*(tl-tmelt)/(tl-bt))
qsl = rd/rv*es/(pressure(k)-(1-rd/rv)*es)
b1 = rlv**2/(tl**2*cp*rv)
qs = qsl*(1.+b1*qt(i,j,k))/(1.+b1*qsl)
ql(i,j,k) = dim(qt(i,j,k)-qs,0.)
end do
end do
end do
end if
return
end subroutine thermo