module wv_saturation
use shr_kind_mod, only: r8 => shr_kind_r8
use module_cam_support, only: endrun, iulog
use module_wrf_error
implicit none
private
save
public gestbl
public estblf
public aqsat
public vqsatd
public fqsatd
public qsat_water
public vqsatd_water
public vqsatd2
public vqsatd2_single
public polysvp
public hlatv, tmin, hlatf, rgasv, pcf, cp, epsqs, ttrice
integer plenest
parameter (plenest=250)
real(r8) estbl(plenest)
real(r8) tmin
real(r8) tmax
real(r8) ttrice
real(r8) pcf(6)
real(r8) epsqs
real(r8) rgasv
real(r8) hlatf
real(r8) hlatv
real(r8) cp
real(r8) tmelt
logical icephs
contains
real(r8) function estblf( td )
real(r8), intent(in) :: td
real(r8) :: e
real(r8) :: ai
integer :: i
e = max(min(td,tmax),tmin)
i = int(e-tmin)+1
ai = aint(e-tmin)
estblf = (tmin+ai-e+1._r8)* &
estbl(i)-(tmin+ai-e)* &
estbl(i+1)
end function estblf
subroutine gestbl(tmn ,tmx ,trice ,ip ,epsil , &
latvap ,latice ,rh2o ,cpair ,tmeltx )
use module_cam_support, only: masterproc
use module_cam_gffgch
real(r8), intent(in) :: tmn
real(r8), intent(in) :: tmx
real(r8), intent(in) :: epsil
real(r8), intent(in) :: trice
real(r8), intent(in) :: latvap
real(r8), intent(in) :: latice
real(r8), intent(in) :: rh2o
real(r8), intent(in) :: cpair
real(r8), intent(in) :: tmeltx
real(r8) t
integer n
integer lentbl
integer itype
logical ip
tmin = tmn
tmax = tmx
ttrice = trice
icephs = ip
epsqs = epsil
hlatv = latvap
hlatf = latice
rgasv = rh2o
cp = cpair
tmelt = tmeltx
lentbl = INT(tmax-tmin+2.000001_r8)
if (lentbl .gt. plenest) then
write(iulog,9000) tmax, tmin, plenest
call wrf_message(iulog)
call endrun ('GESTBL')
end if
if (icephs) then
if (ttrice /= 0.0_r8) then
itype = -ttrice
else
itype = 1
end if
else
itype = 0
end if
t = tmin - 1.0_r8
do n=1,lentbl
t = t + 1.0_r8
call gffgch(t,estbl(n),itype)
end do
do n=lentbl+1,plenest
estbl(n) = -99999.0_r8
end do
pcf(1) = 5.04469588506e-01_r8
pcf(2) = -5.47288442819e+00_r8
pcf(3) = -3.67471858735e-01_r8
pcf(4) = -8.95963532403e-03_r8
pcf(5) = -7.78053686625e-05_r8
if (masterproc) then
write(iulog,*)' *** SATURATION VAPOR PRESSURE TABLE COMPLETED ***'
call wrf_message(iulog)
end if
return
9000 format('GESTBL: FATAL ERROR *********************************',/, &
' TMAX AND TMIN REQUIRE A LARGER DIMENSION ON THE LENGTH', &
' OF THE SATURATION VAPOR PRESSURE TABLE ESTBL(PLENEST)',/, &
' TMAX, TMIN, AND PLENEST => ', 2f7.2, i3)
end subroutine gestbl
subroutine aqsat(t ,p ,es ,qs ,ii , &
ilen ,kk ,kstart ,kend )
integer, intent(in) :: ii
integer, intent(in) :: kk
integer, intent(in) :: ilen
integer, intent(in) :: kstart
integer, intent(in) :: kend
real(r8), intent(in) :: t(ii,kk)
real(r8), intent(in) :: p(ii,kk)
real(r8), intent(out) :: es(ii,kk)
real(r8), intent(out) :: qs(ii,kk)
real(r8) omeps
integer i, k
omeps = 1.0_r8 - epsqs
do k=kstart,kend
do i=1,ilen
es(i,k) = estblf(t(i,k))
qs(i,k) = epsqs*es(i,k)/(p(i,k) - omeps*es(i,k))
qs(i,k) = min(1.0_r8,qs(i,k))
if (qs(i,k) < 0.0_r8) then
qs(i,k) = 1.0_r8
es(i,k) = p(i,k)
end if
end do
end do
return
end subroutine aqsat
subroutine vqsatd(t ,p ,es ,qs ,gam , &
len )
integer, intent(in) :: len
real(r8), intent(in) :: t(len)
real(r8), intent(in) :: p(len)
real(r8), intent(out) :: es(len)
real(r8), intent(out) :: qs(len)
real(r8), intent(out) :: gam(len)
logical lflg
integer i
real(r8) omeps
real(r8) trinv
real(r8) tc
real(r8) weight
real(r8) hltalt
real(r8) hlatsb
real(r8) hlatvp
real(r8) tterm
real(r8) desdt
omeps = 1.0_r8 - epsqs
do i=1,len
es(i) = estblf(t(i))
qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
qs(i) = min(1.0_r8,qs(i))
if (qs(i) < 0.0_r8) then
qs(i) = 1.0_r8
es(i) = p(i)
end if
end do
trinv = 0.0_r8
if ((.not. icephs) .or. (ttrice.eq.0.0_r8)) go to 10
trinv = 1.0_r8/ttrice
do i=1,len
tc = t(i) - tmelt
lflg = (tc >= -ttrice .and. tc < 0.0_r8)
weight = min(-tc*trinv,1.0_r8)
hlatsb = hlatv + weight*hlatf
hlatvp = hlatv - 2369.0_r8*tc
if (t(i) < tmelt) then
hltalt = hlatsb
else
hltalt = hlatvp
end if
if (lflg) then
tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) + tc*pcf(5))))
else
tterm = 0.0_r8
end if
desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) + tterm*trinv
gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
if (qs(i) == 1.0_r8) gam(i) = 0.0_r8
end do
return
10 do i=1,len
hlatvp = hlatv - 2369.0_r8*(t(i)-tmelt)
if (icephs) then
hlatsb = hlatv + hlatf
else
hlatsb = hlatv
end if
if (t(i) < tmelt) then
hltalt = hlatsb
else
hltalt = hlatvp
end if
desdt = hltalt*es(i)/(rgasv*t(i)*t(i))
gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
if (qs(i) == 1.0_r8) gam(i) = 0.0_r8
end do
return
end subroutine vqsatd
subroutine vqsatd_water(t ,p ,es ,qs ,gam , &
len )
integer, intent(in) :: len
real(r8), intent(in) :: t(len)
real(r8), intent(in) :: p(len)
real(r8), intent(out) :: es(len)
real(r8), intent(out) :: qs(len)
real(r8), intent(out) :: gam(len)
integer i
real(r8) omeps
real(r8) hltalt
real(r8) hlatsb
real(r8) hlatvp
real(r8) desdt
omeps = 1.0_r8 - epsqs
do i=1,len
es(i) = polysvp(t(i),0)
qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
qs(i) = min(1.0_r8,qs(i))
if (qs(i) < 0.0_r8) then
qs(i) = 1.0_r8
es(i) = p(i)
end if
end do
do i=1,len
hlatvp = hlatv - 2369.0_r8*(t(i)-tmelt)
hlatsb = hlatv
if (t(i) < tmelt) then
hltalt = hlatsb
else
hltalt = hlatvp
end if
desdt = hltalt*es(i)/(rgasv*t(i)*t(i))
gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
if (qs(i) == 1.0_r8) gam(i) = 0.0_r8
end do
return
end subroutine vqsatd_water
function polysvp (T,type)
real(r8) dum
real(r8) T,polysvp
integer type
if (type.eq.1) then
polysvp = 10._r8**(-9.09718_r8*(273.16_r8/t-1._r8)-3.56654_r8* &
log10(273.16_r8/t)+0.876793_r8*(1._r8-t/273.16_r8)+ &
log10(6.1071_r8))*100._r8
end if
if (type.eq.0) then
polysvp = 10._r8**(-7.90298_r8*(373.16_r8/t-1._r8)+ &
5.02808_r8*log10(373.16_r8/t)- &
1.3816e-7_r8*(10._r8**(11.344_r8*(1._r8-t/373.16_r8))-1._r8)+ &
8.1328e-3_r8*(10._r8**(-3.49149_r8*(373.16_r8/t-1._r8))-1._r8)+ &
log10(1013.246_r8))*100._r8
end if
end function polysvp
integer function fqsatd(t ,p ,es ,qs ,gam , len )
integer, intent(in) :: len
real(r8), intent(in) :: t(len)
real(r8), intent(in) :: p(len)
real(r8), intent(out) :: es(len)
real(r8), intent(out) :: qs(len)
real(r8), intent(out) :: gam(len)
call vqsatd(t ,p ,es ,qs ,gam , len )
fqsatd = 1
return
end function fqsatd
real(r8) function qsat_water(t,p)
real(r8) t
real(r8) p
real(r8) es
real(r8) ps, ts, e1, e2, f1, f2, f3, f4, f5, f
ps = 1013.246_r8
ts = 373.16_r8
e1 = 11.344_r8*(1.0_r8 - t/ts)
e2 = -3.49149_r8*(ts/t - 1.0_r8)
f1 = -7.90298_r8*(ts/t - 1.0_r8)
f2 = 5.02808_r8*log10(ts/t)
f3 = -1.3816_r8*(10.0_r8**e1 - 1.0_r8)/10000000.0_r8
f4 = 8.1328_r8*(10.0_r8**e2 - 1.0_r8)/1000.0_r8
f5 = log10(ps)
f = f1 + f2 + f3 + f4 + f5
es = (10.0_r8**f)*100.0_r8
qsat_water = epsqs*es/(p-(1.-epsqs)*es)
if(qsat_water < 0.) qsat_water = 1.
end function qsat_water
subroutine vqsatd2(t ,p ,es ,qs ,dqsdt , &
len )
integer, intent(in) :: len
real(r8), intent(in) :: t(len)
real(r8), intent(in) :: p(len)
real(r8), intent(out) :: es(len)
real(r8), intent(out) :: qs(len)
real(r8), intent(out) :: dqsdt(len)
logical lflg
integer i
real(r8) omeps
real(r8) trinv
real(r8) tc
real(r8) weight
real(r8) hltalt
real(r8) hlatsb
real(r8) hlatvp
real(r8) tterm
real(r8) desdt
real(r8) gam(len)
omeps = 1.0_r8 - epsqs
do i=1,len
es(i) = estblf(t(i))
qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
qs(i) = min(1.0_r8,qs(i))
if (qs(i) < 0.0_r8) then
qs(i) = 1.0_r8
es(i) = p(i)
end if
end do
trinv = 0.0_r8
if ((.not. icephs) .or. (ttrice.eq.0.0_r8)) go to 10
trinv = 1.0_r8/ttrice
do i=1,len
tc = t(i) - tmelt
lflg = (tc >= -ttrice .and. tc < 0.0_r8)
weight = min(-tc*trinv,1.0_r8)
hlatsb = hlatv + weight*hlatf
hlatvp = hlatv - 2369.0_r8*tc
if (t(i) < tmelt) then
hltalt = hlatsb
else
hltalt = hlatvp
end if
if (lflg) then
tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) + tc*pcf(5))))
else
tterm = 0.0_r8
end if
desdt = hltalt*es(i)/(rgasv*t(i)*t(i)) + tterm*trinv
gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
if (qs(i) == 1.0_r8) gam(i) = 0.0_r8
dqsdt(i) = (cp/hltalt)*gam(i)
end do
return
10 do i=1,len
hlatvp = hlatv - 2369.0_r8*(t(i)-tmelt)
if (icephs) then
hlatsb = hlatv + hlatf
else
hlatsb = hlatv
end if
if (t(i) < tmelt) then
hltalt = hlatsb
else
hltalt = hlatvp
end if
desdt = hltalt*es(i)/(rgasv*t(i)*t(i))
gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
if (qs(i) == 1.0_r8) gam(i) = 0.0_r8
dqsdt(i) = (cp/hltalt)*gam(i)
end do
return
end subroutine vqsatd2
subroutine vqsatd2_single(t ,p ,es ,qs ,dqsdt)
real(r8), intent(in) :: t
real(r8), intent(in) :: p
real(r8), intent(out) :: es
real(r8), intent(out) :: qs
real(r8), intent(out) :: dqsdt
logical lflg
real(r8) omeps
real(r8) trinv
real(r8) tc
real(r8) weight
real(r8) hltalt
real(r8) hlatsb
real(r8) hlatvp
real(r8) tterm
real(r8) desdt
real(r8) gam
omeps = 1.0_r8 - epsqs
es = estblf(t)
qs = epsqs*es/(p - omeps*es)
qs = min(1.0_r8,qs)
if (qs < 0.0_r8) then
qs = 1.0_r8
es = p
end if
trinv = 0.0_r8
if ((.not. icephs) .or. (ttrice.eq.0.0_r8)) go to 10
trinv = 1.0_r8/ttrice
tc = t - tmelt
lflg = (tc >= -ttrice .and. tc < 0.0_r8)
weight = min(-tc*trinv,1.0_r8)
hlatsb = hlatv + weight*hlatf
hlatvp = hlatv - 2369.0_r8*tc
if (t < tmelt) then
hltalt = hlatsb
else
hltalt = hlatvp
end if
if (lflg) then
tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) + tc*pcf(5))))
else
tterm = 0.0_r8
end if
desdt = hltalt*es/(rgasv*t*t) + tterm*trinv
gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
if (qs == 1.0_r8) gam = 0.0_r8
dqsdt = (cp/hltalt)*gam
return
10 continue
hlatvp = hlatv - 2369.0_r8*(t-tmelt)
if (icephs) then
hlatsb = hlatv + hlatf
else
hlatsb = hlatv
end if
if (t < tmelt) then
hltalt = hlatsb
else
hltalt = hlatvp
end if
desdt = hltalt*es/(rgasv*t*t)
gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
if (qs == 1.0_r8) gam = 0.0_r8
dqsdt = (cp/hltalt)*gam
return
end subroutine vqsatd2_single
end module wv_saturation