module wv_saturation use shr_kind_mod, only: r8 => shr_kind_r8 use module_cam_support, only: & masterproc, & endrun, & iulog use module_wrf_error implicit none private save public gestbl public estblf public aqsat public fqsatd 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_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_debug(100,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 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 end module wv_saturation