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)
!! X. Long: This is a fix to tackle incorrect thls input. The reason why tl is going less than 100K (unphysical)
!! is that it is calculated from thl which dose not change over time here.Probably this is happening at
!'calc_halflev and call thermo(thl0h,qt0h,ql0h,presh,exnh)'. This problem should be fixed later.
if (tl<100.0) then
tl=100.0
end if
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