! ! Common block and statement functions for saturation vapor pressure ! look-up procedure, J. J. Hack, February 1990 ! ! $Id: wv_saturation.F90,v 1.1.12.2.10.1 2014-04-14 16:04:56 dbarahon Exp $ ! module wv_saturation #ifdef GEOS5 use MAPL_ConstantsMod, r8 => MAPL_R8 #endif #ifdef NEMS_GSM use funcphys, only : fpvsl, fpvsi, fpvs ! saturation vapor pressure for water & ice #endif use machine, only : r8 => kind_phys !++jtb (comm out) !--jtb implicit none private save ! ! Public interfaces ! public gestbl public estblf public aqsat public aqsatd public vqsatd public fqsatd public qsat_water public vqsat_water public qsat_ice public vqsat_ice public vqsatd_water public aqsat_water public vqsatd2_water public vqsatd2_water_single public vqsatd2_ice_single public vqsatd2 public vqsatd2_single public polysvp ! ! Data used by cldwat ! public hlatv, tmin, hlatf, rgasv, pcf, cp, epsqs, ttrice ! ! Data ! integer plenest parameter (plenest=250) ! ! Table of saturation vapor pressure values es from tmin degrees ! to tmax+1 degrees k in one degree increments. ttrice defines the ! transition region where es is a combination of ice & water values ! 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 integer, parameter :: iulog=6 contains real(r8) function estblf( td ) ! ! Saturation vapor pressure table lookup ! 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 ) !----------------------------------------------------------------------- ! ! Purpose: ! Builds saturation vapor pressure table for later lookup procedure. ! ! Method: ! Uses Goff & Gratch (1946) relationships to generate the table ! according to a set of free parameters defined below. Auxiliary ! routines are also included for making rapid estimates (well with 1%) ! of both es and d(es)/dt for the particular table configuration. ! ! Author: J. Hack ! !----------------------------------------------------------------------- !------------------------------Arguments-------------------------------- ! ! Input arguments ! 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 ! !---------------------------Local variables----------------------------- ! real(r8) t integer n integer lentbl integer itype ! 1 -> ice phase, no transition ! -x -> ice phase, x degree transition logical ip ! !----------------------------------------------------------------------- ! ! Set es table parameters ! tmin = tmn tmax = tmx ttrice = trice icephs = ip ! ! Set physical constants required for es calculation ! 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(*,*) "AHHH wv_sat" STOP end if ! ! Begin building es table. ! Check whether ice phase requested. ! If so, set appropriate transition range for temperature ! 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),tmelt,itype) end do ! do n=lentbl+1,plenest estbl(n) = -99999.0_r8 end do ! ! Table complete -- Set coefficients for polynomial approximation of ! difference between saturation vapor press over water and saturation ! pressure over ice for -ttrice < t < 0 (degrees C). NOTE: polynomial ! is valid in the range -40 < t < 0 (degrees C). ! ! --- Degree 5 approximation --- ! 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 ! ! --- Degree 6 approximation --- ! !-----pcf(1) = 7.63285250063e-02 !-----pcf(2) = -5.86048427932e+00 !-----pcf(3) = -4.38660831780e-01 !-----pcf(4) = -1.37898276415e-02 !-----pcf(5) = -2.14444472424e-04 !-----pcf(6) = -1.36639103771e-06 ! !++jtb (comment out) ! if (masterproc) then ! !!write(iulog,*)' *** SATURATION VAPOR PRESSURE TABLE COMPLETED ***' ! end if !--jtb 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 ) !----------------------------------------------------------------------- ! ! Purpose: ! Utility procedure to look up and return saturation vapor pressure from ! precomputed table, calculate and return saturation specific humidity ! (g/g),for input arrays of temperature and pressure (dimensioned ii,kk) ! This routine is useful for evaluating only a selected region in the ! vertical. ! ! Method: ! <Describe the algorithm(s) used in the routine.> ! <Also include any applicable external references.> ! ! Author: J. Hack ! !------------------------------Arguments-------------------------------- ! ! Input arguments ! 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) ! ! Output arguments ! real(r8), intent(out) :: es(ii,kk) real(r8), intent(out) :: qs(ii,kk) ! !---------------------------Local workspace----------------------------- ! real(r8) omeps integer i, k ! !----------------------------------------------------------------------- ! omeps = 1.0_r8 - epsqs do k=kstart,kend do i=1,ilen es(i,k) = min(estblf(t(i,k)),p(i,k)) ! ! Saturation specific humidity ! qs(i,k) = min(1.0_r8, epsqs*es(i,k)/(p(i,k)-omeps*es(i,k))) ! ! The following check is to avoid the generation of negative values ! that can occur in the upper stratosphere and mesosphere ! ! 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 !++xl subroutine aqsat_water(t, p, es, qs, ii, ilen, kk, kstart,kend) !----------------------------------------------------------------------- ! ! Purpose: ! Utility procedure to look up and return saturation vapor pressure from ! precomputed table, calculate and return saturation specific humidity ! (g/g),for input arrays of temperature and pressure (dimensioned ii,kk) ! This routine is useful for evaluating only a selected region in the ! vertical. ! ! Method: ! <Describe the algorithm(s) used in the routine.> ! <Also include any applicable external references.> ! ! Author: J. Hack ! !------------------------------Arguments-------------------------------- ! ! Input arguments ! 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) ! ! Output arguments ! real(r8), intent(out) :: es(ii,kk) real(r8), intent(out) :: qs(ii,kk) ! !---------------------------Local workspace----------------------------- ! 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)) #ifdef GEOS5 es(i,k) = min(polysvp(t(i,k),0), p(i,k)) #endif #ifdef NEMS_GSM es(i,k) = min(fpvsl(t(i,k)), p(i,k)) #endif ! ! Saturation specific humidity ! qs(i,k) = min(1.0_r8, epsqs*es(i,k)/(p(i,k)-omeps*es(i,k))) ! ! The following check is to avoid the generation of negative values ! that can occur in the upper stratosphere and mesosphere ! ! 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_water !--xl subroutine aqsatd(t, p, es, qs, gam, ii, ilen, kk, kstart, kend) !----------------------------------------------------------------------- ! ! Purpose: ! Utility procedure to look up and return saturation vapor pressure from ! precomputed table, calculate and return saturation specific humidity ! (g/g). ! ! Method: ! Differs from aqsat by also calculating and returning ! gamma (l/cp)*(d(qsat)/dT) ! Input arrays temperature and pressure (dimensioned ii,kk). ! ! Author: J. Hack ! !------------------------------Arguments-------------------------------- ! ! Input arguments ! integer, intent(in) :: ii integer, intent(in) :: ilen integer, intent(in) :: kk integer, intent(in) :: kstart integer, intent(in) :: kend real(r8), intent(in) :: t(ii,kk) real(r8), intent(in) :: p(ii,kk) ! ! Output arguments ! real(r8), intent(out) :: es(ii,kk) real(r8), intent(out) :: qs(ii,kk) real(r8), intent(out) :: gam(ii,kk) ! !---------------------------Local workspace----------------------------- ! logical lflg integer i integer k 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 k=kstart,kend do i=1,ilen es(i,k) = min(p(i,k), estblf(t(i,k))) ! ! Saturation specific humidity ! qs(i,k) = min(1.0_r8, epsqs*es(i,k)/(p(i,k)-omeps*es(i,k))) ! ! The following check is to avoid the generation of negative qs ! values which can occur in the upper stratosphere and mesosphere ! ! ! 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 ! ! "generalized" analytic expression for t derivative of es ! accurate to within 1 percent for 173.16 < t < 373.16 ! trinv = 0.0_r8 if ((.not. icephs) .or. (ttrice == 0.0_r8)) go to 10 trinv = 1.0_r8/ttrice ! do k=kstart,kend do i=1,ilen ! ! Weighting of hlat accounts for transition from water to ice ! polynomial expression approximates difference between es over ! water and es over ice from 0 to -ttrice (C) (min of ttrice is ! -40): required for accurate estimate of es derivative in transition ! range from ice to water also accounting for change of hlatv with t ! above freezing where constant slope is given by -2369 j/(kg c) =cpv - cw ! tc = t(i,k) - 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,k) < 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,k)/(rgasv*t(i,k)*t(i,k)) + tterm*trinv gam(i,k) = hltalt*qs(i,k)*p(i,k)*desdt/(cp*es(i,k)*(p(i,k) & & - omeps*es(i,k))) if (qs(i,k) == 1.0_r8) gam(i,k) = 0.0_r8 end do end do ! go to 20 ! ! No icephs or water to ice transition ! 10 do k=kstart,kend do i=1,ilen ! ! Account for change of hlatv with t above freezing where ! constant slope is given by -2369 j/(kg c) = cpv - cw ! hlatvp = hlatv - 2369.0_r8*(t(i,k)-tmelt) if (icephs) then hlatsb = hlatv + hlatf else hlatsb = hlatv end if if (t(i,k) < tmelt) then hltalt = hlatsb else hltalt = hlatvp end if desdt = hltalt*es(i,k)/(rgasv*t(i,k)*t(i,k)) gam(i,k) = hltalt*qs(i,k)*p(i,k)*desdt/(cp*es(i,k)*(p(i,k) & & - omeps*es(i,k))) if (qs(i,k) == 1.0_r8) gam(i,k) = 0.0_r8 end do end do ! 20 return end subroutine aqsatd subroutine vqsatd(t ,p ,es ,qs ,gam , len ) !----------------------------------------------------------------------- ! ! Purpose: ! Utility procedure to look up and return saturation vapor pressure from ! precomputed table, calculate and return saturation specific humidity ! (g/g), and calculate and return gamma (l/cp)*(d(qsat)/dT). The same ! function as qsatd, but operates on vectors of temperature and pressure ! ! Method: ! ! Author: J. Hack ! !------------------------------Arguments-------------------------------- ! ! Input arguments ! integer, intent(in) :: len real(r8), intent(in) :: t(len) real(r8), intent(in) :: p(len) ! ! Output arguments ! real(r8), intent(out) :: es(len) real(r8), intent(out) :: qs(len) real(r8), intent(out) :: gam(len) ! !--------------------------Local Variables------------------------------ ! 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) = min(estblf(t(i)), p(i)) ! ! Saturation specific humidity ! qs(i) = epsqs*es(i)/(p(i) - omeps*es(i)) ! ! The following check is to avoid the generation of negative ! values that can occur in the upper stratosphere and mesosphere ! 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 ! ! "generalized" analytic expression for t derivative of es ! accurate to within 1 percent for 173.16 < t < 373.16 ! 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 ! ! Weighting of hlat accounts for transition from water to ice ! polynomial expression approximates difference between es over ! water and es over ice from 0 to -ttrice (C) (min of ttrice is ! -40): required for accurate estimate of es derivative in transition ! range from ice to water also accounting for change of hlatv with t ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw ! 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 ! ! No icephs or water to ice transition ! 10 do i=1,len ! ! Account for change of hlatv with t above freezing where ! constant slope is given by -2369 j/(kg c) = cpv - cw ! 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 !++xl subroutine vqsatd_water(t, p, es, qs, gam, len) !------------------------------Arguments-------------------------------- ! ! Input arguments ! integer, intent(in) :: len real(r8), intent(in) :: t(len) real(r8), intent(in) :: p(len) ! ! Output arguments ! real(r8), intent(out) :: es(len) real(r8), intent(out) :: qs(len) real(r8), intent(out) :: gam(len) ! !--------------------------Local Variables------------------------------ ! ! 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 #ifdef NEMS_GSM es(i) = min(fpvsl(t(i)), p(i)) #else es(i) = min(polysvp(t(i),0), p(i)) #endif ! ! Saturation specific humidity ! qs(i) = min(1.0_r8, epsqs*es(i) / (p(i)-omeps*es(i))) ! ! The following check is to avoid the generation of negative ! values that can occur in the upper stratosphere and mesosphere ! ! 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 ! ! No icephs or water to ice transition ! do i=1,len ! ! Account for change of hlatv with t above freezing where ! constant slope is given by -2369 j/(kg c) = cpv - cw ! 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,typ) ! Compute saturation vapor pressure by using ! function from Goff and Gatch (1946) ! Polysvp returned in units of pa. ! T is input in units of K. ! type refers to saturation with respect to liquid (0) or ice (1) !!DONIFF Changed to Murphy and Koop (2005) (03/04/14) real(r8) dum real(r8) t,polysvp integer typ if (.true.) then !ice if (typ == 1) then polysvp = MurphyKoop_svp_ice(t) end if if (typ == 0) then polysvp = MurphyKoop_svp_water(t) end if else ! ice if (typ.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 (typ.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 if end function polysvp !--xl 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 vqsat_water(t,p,qsat_water,len) integer, intent(in) :: len real(r8) t(len) real(r8) p(len) real(r8) qsat_water(len) real(r8) es real(r8), parameter :: t0inv = 1._r8/273._r8 real(r8) coef integer :: i coef = hlatv/rgasv do i=1,len es = 611._r8*exp(coef*(t0inv-1./t(i))) qsat_water(i) = epsqs*es/(p(i)-(1.-epsqs)*es) if(qsat_water(i) < 0.) qsat_water(i) = 1. enddo return end subroutine vqsat_water real(r8) function qsat_ice(t,p) real(r8) t real(r8) p real(r8) es real(r8), parameter :: t0inv = 1._r8/273._r8 es = 611.*exp((hlatv+hlatf)/rgasv*(t0inv-1./t)) qsat_ice = epsqs*es/(p-(1.-epsqs)*es) if(qsat_ice < 0.) qsat_ice = 1. end function qsat_ice subroutine vqsat_ice(t,p,qsat_ice,len) integer,intent(in) :: len real(r8) t(len) real(r8) p(len) real(r8) qsat_ice(len) real(r8) es real(r8), parameter :: t0inv = 1._r8/273._r8 real(r8) coef integer :: i coef = (hlatv+hlatf)/rgasv do i=1,len es = 611.*exp(coef*(t0inv-1./t(i))) qsat_ice(i) = epsqs*es/(p(i)-(1.-epsqs)*es) if(qsat_ice(i) < 0.) qsat_ice(i) = 1. enddo return end subroutine vqsat_ice ! Sungsu ! Below two subroutines (vqsatd2_water,vqsatd2_water_single) are by Sungsu ! Replace 'gam -> dqsdt' ! Sungsu subroutine vqsatd2_water(t ,p ,es ,qs ,dqsdt , len ) !------------------------------Arguments-------------------------------- ! ! Input arguments ! integer, intent(in) :: len real(r8), intent(in) :: t(len) real(r8), intent(in) :: p(len) ! ! Output arguments ! real(r8), intent(out) :: es(len) real(r8), intent(out) :: qs(len) real(r8), intent(out) :: dqsdt(len) ! !--------------------------Local Variables------------------------------ ! ! integer i ! real(r8) omeps real(r8) hltalt ! real(r8) hlatsb real(r8) hlatvp real(r8) desdt real(r8) gam(len) ! !----------------------------------------------------------------------- ! omeps = 1.0_r8 - epsqs do i=1,len #ifdef GEOS5 es(i) = min(polysvp(t(i),0), p(i)) #endif #ifdef NEMS_GSM es(i) = min(fpvsl(t(i)), p(i)) #endif ! ! Saturation specific humidity ! qs(i) = epsqs*es(i)/(p(i) - omeps*es(i)) ! ! The following check is to avoid the generation of negative ! values that can occur in the upper stratosphere and mesosphere ! 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 ! ! No icephs or water to ice transition ! do i=1,len ! ! Account for change of hlatv with t above freezing where ! constant slope is given by -2369 j/(kg c) = cpv - cw ! 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 dqsdt(i) = (cp/hltalt)*gam(i) end do ! return ! end subroutine vqsatd2_water subroutine vqsatd2_water_single(t ,p ,es ,qs ,dqsdt) !------------------------------Arguments-------------------------------- ! ! Input arguments ! real(r8), intent(in) :: t, p ! ! Output arguments ! real(r8), intent(out) :: es, qs, dqsdt ! !--------------------------Local Variables------------------------------ ! ! integer i ! real(r8) omeps, hltalt, hlatsb, hlatvp, desdt, gam ! !----------------------------------------------------------------------- ! omeps = 1.0_r8 - epsqs ! do i=1,len #ifdef GEOS5 es = min(p, polysvp(t,0)) #endif #ifdef NEMS_GSM es = min(p, fpvsl(t)) #endif ! ! Saturation specific humidity ! qs = min(1.0_r8, epsqs*es/(p-omeps*es)) ! ! The following check is to avoid the generation of negative ! values that can occur in the upper stratosphere and mesosphere ! ! if (qs < 0.0_r8) then ! qs = 1.0_r8 ! es = p ! end if ! end do ! ! No icephs or water to ice transition ! ! do i=1,len ! ! Account for change of hlatv with t above freezing where ! constant slope is given by -2369 j/(kg c) = cpv - cw ! hlatvp = hlatv - 2369.0_r8*(t-tmelt) hlatsb = hlatv 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 ! end do ! return ! end subroutine vqsatd2_water_single subroutine vqsatd2(t ,p ,es ,qs ,dqsdt , len ) !----------------------------------------------------------------------- ! Sungsu : This is directly copied from 'vqsatd' but 'dqsdt' is output ! instead of gam for use in Sungsu's equilibrium stratiform ! macrophysics scheme. ! ! Purpose: ! Utility procedure to look up and return saturation vapor pressure from ! precomputed table, calculate and return saturation specific humidity ! (g/g), and calculate and return gamma (l/cp)*(d(qsat)/dT). The same ! function as qsatd, but operates on vectors of temperature and pressure ! ! Method: ! ! Author: J. Hack ! !------------------------------Arguments-------------------------------- ! ! Input arguments ! integer, intent(in) :: len real(r8), intent(in) :: t(len) real(r8), intent(in) :: p(len) ! ! Output arguments ! real(r8), intent(out) :: es(len) real(r8), intent(out) :: qs(len) real(r8), intent(out) :: dqsdt(len) ! !--------------------------Local Variables------------------------------ ! 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 #ifdef GEOS5 es(i) = min(p(i), estblf(t(i))) #endif #ifdef NEMS_GSM es(i) = min(p(i), fpvsi(t(i))) #endif ! ! Saturation specific humidity ! qs(i) = epsqs*es(i)/(p(i) - omeps*es(i)) ! ! The following check is to avoid the generation of negative ! values that can occur in the upper stratosphere and mesosphere ! 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 ! ! "generalized" analytic expression for t derivative of es ! accurate to within 1 percent for 173.16 < t < 373.16 ! trinv = 0.0_r8 if ((.not. icephs) .or. (ttrice == 0.0_r8)) go to 10 trinv = 1.0_r8/ttrice do i=1,len ! ! Weighting of hlat accounts for transition from water to ice ! polynomial expression approximates difference between es over ! water and es over ice from 0 to -ttrice (C) (min of ttrice is ! -40): required for accurate estimate of es derivative in transition ! range from ice to water also accounting for change of hlatv with t ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw ! 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 ! ! No icephs or water to ice transition ! 10 do i=1,len ! ! Account for change of hlatv with t above freezing where ! constant slope is given by -2369 j/(kg c) = cpv - cw ! 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 ! Below routine is by Sungsu subroutine vqsatd2_single(t ,p ,es ,qs ,dqsdt) !----------------------------------------------------------------------- ! Sungsu : This is directly copied from 'vqsatd' but 'dqsdt' is output ! instead of gam for use in Sungsu's equilibrium stratiform ! macrophysics scheme. ! ! Purpose: ! Utility procedure to look up and return saturation vapor pressure from ! precomputed table, calculate and return saturation specific humidity ! (g/g), and calculate and return gamma (l/cp)*(d(qsat)/dT). The same ! function as qsatd, but operates on vectors of temperature and pressure ! ! Method: ! ! Author: J. Hack ! !------------------------------Arguments-------------------------------- ! ! Input arguments ! real(r8), intent(in) :: t, p ! ! Output arguments ! real(r8), intent(out) :: es, qs, dqsdt ! !--------------------------Local Variables------------------------------ ! logical lflg ! ! integer i ! index for vector calculations ! 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 ! do i=1,len #ifdef GEOS5 es = estblf(t) #endif #ifdef NEMS_GSM es = min(fpvs(t), p) #endif ! ! Saturation specific humidity ! qs = epsqs*es/(p - omeps*es) ! ! The following check is to avoid the generation of negative ! values that can occur in the upper stratosphere and mesosphere ! qs = min(1.0_r8,qs) ! ! if (qs < 0.0_r8) then ! qs = 1.0_r8 ! es = p ! end if ! end do ! ! "generalized" analytic expression for t derivative of es ! accurate to within 1 percent for 173.16 < t < 373.16 ! trinv = 0.0_r8 if ((.not. icephs) .or. (ttrice == 0.0_r8)) go to 10 trinv = 1.0_r8/ttrice ! do i=1,len ! ! Weighting of hlat accounts for transition from water to ice ! polynomial expression approximates difference between es over ! water and es over ice from 0 to -ttrice (C) (min of ttrice is ! -40): required for accurate estimate of es derivative in transition ! range from ice to water also accounting for change of hlatv with t ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw ! 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 ! end do return ! ! No icephs or water to ice transition ! 10 continue !10 do i=1,len ! ! Account for change of hlatv with t above freezing where ! constant slope is given by -2369 j/(kg c) = cpv - cw ! 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 ! end do ! return ! end subroutine vqsatd2_single !---------------------------------------------------------------------- !---------------------------------------------------------------------- subroutine gffgch(t ,es ,tmelt ,itype ) !----------------------------------------------------------------------- ! ! Purpose: ! Computes saturation vapor pressure over water and/or over ice using ! Goff & Gratch (1946) relationships. ! <Say what the routine does> ! ! Method: ! T (temperature), and itype are input parameters, while es (saturation ! vapor pressure) is an output parameter. The input parameter itype ! serves two purposes: a value of zero indicates that saturation vapor ! pressures over water are to be returned (regardless of temperature), ! while a value of one indicates that saturation vapor pressures over ! ice should be returned when t is less than freezing degrees. If itype ! is negative, its absolute value is interpreted to define a temperature ! transition region below freezing in which the returned ! saturation vapor pressure is a weighted average of the respective ice ! and water value. That is, in the temperature range 0 => -itype ! degrees c, the saturation vapor pressures are assumed to be a weighted ! average of the vapor pressure over supercooled water and ice (all ! water at 0 c; all ice at -itype c). Maximum transition range => 40 c ! ! Author: J. Hack ! !----------------------------------------------------------------------- implicit none !------------------------------Arguments-------------------------------- ! ! Input arguments ! real(r8), intent(in) :: t ,tmelt ! ! Output arguments ! integer, intent(inout) :: itype real(r8), intent(out) :: es ! !---------------------------Local variables----------------------------- ! real(r8) e1 real(r8) e2 real(r8) eswtr real(r8) f real(r8) f1 real(r8) f2 real(r8) f3 real(r8) f4 real(r8) f5 real(r8) ps real(r8) t0 real(r8) term1 real(r8) term2 real(r8) term3 real(r8) tr real(r8) ts real(r8) weight integer itypo ! !----------------------------------------------------------------------- ! ! Check on whether there is to be a transition region for es ! if (itype < 0) then tr = abs(real(itype,r8)) itypo = itype itype = 1 else tr = 0.0_r8 itypo = itype end if if (tr > 40.0_r8) then write(iulog,900) tr end if ! if(t < (tmelt - tr) .and. itype == 1) go to 10 ! ! Water ! 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 eswtr = es ! if(t >= tmelt .or. itype == 0) go to 20 ! ! Ice ! 10 continue t0 = tmelt term1 = 2.01889049_r8/(t0/t) term2 = 3.56654_r8*log(t0/t) term3 = 20.947031_r8*(t0/t) es = 575.185606e10_r8*exp(-(term1 + term2 + term3)) ! if (t < (tmelt - tr)) go to 20 ! ! Weighted transition between water and ice ! weight = min((tmelt - t)/tr,1.0_r8) es = weight*es + (1.0_r8 - weight)*eswtr ! 20 continue itype = itypo return ! 900 format('GFFGCH: FATAL ERROR ******************************',/, & & 'TRANSITION RANGE FOR WATER TO ICE SATURATION VAPOR', ' PRESSURE, & & TR, EXCEEDS MAXIMUM ALLOWABLE VALUE OF', ' 40.0 DEGREES C',/, & & ' TR = ',f7.2) ! end subroutine gffgch !!DONIF USe Murphy and Koop (2005) (Written by Andrew Gettelman) function MurphyKoop_svp_water(tx) result(es) real(r8), intent(in) :: tx real(r8) :: es real(r8):: t t=min(tx, 332.0_r8) t=max(123.0_r8, tx) es = exp(54.842763_r8 - (6763.22_r8 / t) - (4.210_r8 * log(t)) + & & (0.000367_r8 * t) + (tanh(0.0415_r8 * (t - 218.8_r8)) * & & (53.878_r8 - (1331.22_r8 / t) - (9.44523_r8 * log(t)) + & & 0.014025_r8 * t))) end function MurphyKoop_svp_water function MurphyKoop_svp_ice(tx) result(es) real(r8), intent(in) :: tx real(r8) :: t real(r8) :: es t=max(100.0_r8, tx) t=min(274.0_r8, tx) es = exp(9.550426_r8 - (5723.265_r8 / t) + (3.53068_r8 * & & log(t)) - (0.00728332_r8 * t)) end function MurphyKoop_svp_ice !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine vqsatd2_ice_single(t ,p ,es ,qs ,dqsdt) !------------------------------Arguments-------------------------------- ! ! Input arguments ! real(r8), intent(in) :: t, p ! ! Output arguments ! real(r8), intent(out) :: es, qs, dqsdt ! !--------------------------Local Variables------------------------------ ! ! integer i ! real(r8) omeps, hltalt, hlatsb, hlatvp, desdt, gam ! !----------------------------------------------------------------------- ! omeps = 1.0_r8 - epsqs ! do i=1,len #ifdef GEOS5 es = min(polysvp(t,1),p) #endif #ifdef NEMS_GSM es = min(fpvsi(t),p) #endif ! ! Saturation specific humidity ! qs = min(1.0_r8, epsqs*es/(p-omeps*es)) ! ! The following check is to avoid the generation of negative ! values that can occur in the upper stratosphere and mesosphere ! ! if (qs < 0.0_r8) then ! qs = 1.0_r8 ! es = p ! end if ! end do ! ! No icephs or water to ice transition ! ! do i=1,len ! ! Account for change of hlatv with t above freezing where ! constant slope is given by -2369 j/(kg c) = cpv - cw ! hltalt = hlatv + hlatf desdt = hltalt*es/(rgasv*t*t) if (qs < 1.0_r8) then gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es)) else gam = 0.0_r8 endif dqsdt = (cp/hltalt)*gam ! end do ! return ! end subroutine vqsatd2_ice_single !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end module wv_saturation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!