module wv_saturation !--------------------------------------------------------------------! ! Module Overview: ! ! ! ! This module provides an interface to wv_sat_methods, providing ! ! saturation vapor pressure and related calculations to CAM. ! ! ! ! The original wv_saturation codes were introduced by J. J. Hack, ! ! February 1990. The code has been extensively rewritten since then, ! ! including a total refactoring in Summer 2012. ! ! ! !--------------------------------------------------------------------! ! Methods: ! ! ! ! Pure water/ice saturation vapor pressures are calculated on the ! ! fly, with the specific method determined by a runtime option. ! ! Mixed phase SVP is interpolated from the internal table, estbl, ! ! which is created during initialization. ! ! ! ! The default method for calculating SVP is determined by a namelist ! ! option, and used whenever svp_water/ice or qsat are called. ! ! ! !--------------------------------------------------------------------! use shr_kind_mod, only: r8 => shr_kind_r8 use physconst, only: epsilo, & latvap, & latice, & rh2o, & cpair, & tmelt, & h2otrip use wv_sat_methods, only: & svp_to_qsat => wv_sat_svp_to_qsat implicit none private save ! Public interfaces ! Namelist, initialization, finalization public wv_sat_readnl public wv_sat_init public wv_sat_final ! Saturation vapor pressure calculations public svp_water public svp_ice ! Mixed phase (water + ice) saturation vapor pressure table lookup public estblf public svp_to_qsat ! Subroutines that return both SVP and humidity ! Optional arguments do temperature derivatives public qsat ! Mixed phase public qsat_water ! SVP over water only public qsat_ice ! SVP over ice only ! Wet bulb temperature solver public :: findsp_vc, findsp ! Data ! This value is slightly high, but it seems to be the value for the ! steam point of water originally (and most frequently) used in the ! Goff & Gratch scheme. real(r8), parameter :: tboil = 373.16_r8 ! Table of saturation vapor pressure values (estbl) from tmin to ! tmax+1 Kelvin, in one degree increments. ttrice defines the ! transition region, estbl contains a combination of ice & water ! values. ! Make these public parameters in case another module wants to see the ! extent of the table. real(r8), public, parameter :: tmin = 127.16_r8 real(r8), public, parameter :: tmax = 375.16_r8 real(r8), parameter :: ttrice = 20.00_r8 ! transition range from es over H2O to es over ice integer :: plenest ! length of estbl real(r8), allocatable :: estbl(:) ! table values of saturation vapor pressure real(r8) :: omeps ! 1.0_r8 - epsilo real(r8) :: c3 ! parameter used by findsp ! 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). real(r8) :: pcf(5) = (/ & 5.04469588506e-01_r8, & -5.47288442819e+00_r8, & -3.67471858735e-01_r8, & -8.95963532403e-03_r8, & -7.78053686625e-05_r8 /) ! --- Degree 6 approximation --- ! real(r8) :: pcf(6) = (/ & ! 7.63285250063e-02, & ! 5.86048427932e+00, & ! 4.38660831780e-01, & ! 1.37898276415e-02, & ! 2.14444472424e-04, & ! 1.36639103771e-06 /) contains !--------------------------------------------------------------------- ! ADMINISTRATIVE FUNCTIONS !--------------------------------------------------------------------- subroutine wv_sat_readnl(nlfile) !------------------------------------------------------------------! ! Purpose: ! ! Get runtime options for wv_saturation. ! !------------------------------------------------------------------! use wv_sat_methods, only: wv_sat_get_scheme_idx, & wv_sat_valid_idx, & wv_sat_set_default use spmd_utils, only: masterproc use namelist_utils, only: find_group_name use units, only: getunit, freeunit use mpishorthand use cam_abortutils, only: endrun character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input ! Local variables integer :: unitn, ierr character(len=32) :: wv_sat_scheme = "GoffGratch" character(len=*), parameter :: subname = 'wv_sat_readnl' namelist /wv_sat_nl/ wv_sat_scheme !----------------------------------------------------------------------------- if (masterproc) then unitn = getunit() open( unitn, file=trim(nlfile), status='old' ) call find_group_name(unitn, 'wv_sat_nl', status=ierr) if (ierr == 0) then read(unitn, wv_sat_nl, iostat=ierr) if (ierr /= 0) then call endrun(subname // ':: ERROR reading namelist') return end if end if close(unitn) call freeunit(unitn) end if #ifdef SPMD call mpibcast(wv_sat_scheme, len(wv_sat_scheme) , mpichar, 0, mpicom) #endif if (.not. wv_sat_set_default(wv_sat_scheme)) then call endrun('wv_sat_readnl :: Invalid wv_sat_scheme.') return end if end subroutine wv_sat_readnl subroutine wv_sat_init !------------------------------------------------------------------! ! Purpose: ! ! Initialize module (e.g. setting parameters, initializing the ! ! SVP lookup table). ! !------------------------------------------------------------------! use wv_sat_methods, only: wv_sat_methods_init, & wv_sat_get_scheme_idx, & wv_sat_valid_idx use spmd_utils, only: masterproc use cam_logfile, only: iulog use cam_abortutils, only: endrun use shr_assert_mod, only: shr_assert_in_domain use error_messages, only: handle_errmsg integer :: status ! For wv_sat_methods error reporting. character(len=256) :: errstring ! For generating internal SVP table. real(r8) :: t ! Temperature integer :: i ! Increment counter ! Precalculated because so frequently used. omeps = 1.0_r8 - epsilo ! Transition range method is only valid for transition temperatures at: ! -40 deg C < T < 0 deg C call shr_assert_in_domain(ttrice, ge=0._r8, le=40._r8, varname="ttrice",& msg="wv_sat_init: Invalid transition temperature range.") ! This parameter uses a hardcoded 287.04_r8? c3 = 287.04_r8*(7.5_r8*log(10._r8))/cpair ! Init "methods" module containing actual SVP formulae. call wv_sat_methods_init(r8, tmelt, h2otrip, tboil, ttrice, & epsilo, errstring) call handle_errmsg(errstring, subname="wv_sat_methods_init") ! Add two to make the table slightly too big, just in case. plenest = ceiling(tmax-tmin) + 2 ! Allocate SVP table. allocate(estbl(plenest), stat=status) if (status /= 0) then call endrun('wv_sat_init :: ERROR allocating saturation vapor pressure table') return end if do i = 1, plenest estbl(i) = svp_trans(tmin + real(i-1,r8)) end do if (masterproc) then write(iulog,*)' *** SATURATION VAPOR PRESSURE TABLE COMPLETED ***' end if end subroutine wv_sat_init subroutine wv_sat_final !------------------------------------------------------------------! ! Purpose: ! ! Deallocate global variables in module. ! !------------------------------------------------------------------! use cam_abortutils, only: endrun integer :: status if (allocated(estbl)) then deallocate(estbl, stat=status) if (status /= 0) then call endrun('wv_sat_final :: ERROR deallocating table') return end if end if end subroutine wv_sat_final !--------------------------------------------------------------------- ! DEFAULT SVP FUNCTIONS !--------------------------------------------------------------------- ! Compute saturation vapor pressure over water elemental function svp_water(t) result(es) use wv_sat_methods, only: & wv_sat_svp_water real(r8), intent(in) :: t ! Temperature (K) real(r8) :: es ! SVP (Pa) es = wv_sat_svp_water(T) end function svp_water ! Compute saturation vapor pressure over ice elemental function svp_ice(t) result(es) use wv_sat_methods, only: & wv_sat_svp_ice real(r8), intent(in) :: t ! Temperature (K) real(r8) :: es ! SVP (Pa) es = wv_sat_svp_ice(T) end function svp_ice ! Compute saturation vapor pressure with an ice-water transition elemental function svp_trans(t) result(es) use wv_sat_methods, only: & wv_sat_svp_trans real(r8), intent(in) :: t ! Temperature (K) real(r8) :: es ! SVP (Pa) es = wv_sat_svp_trans(T) end function svp_trans !--------------------------------------------------------------------- ! UTILITIES !--------------------------------------------------------------------- ! Does linear interpolation from nearest values found ! in the table (estbl). elemental function estblf(t) result(es) real(r8), intent(in) :: t ! Temperature real(r8) :: es ! SVP (Pa) integer :: i ! Index for t in the table real(r8) :: t_tmp ! intermediate temperature for es look-up real(r8) :: weight ! Weight for interpolation t_tmp = max(min(t,tmax)-tmin, 0._r8) ! Number of table entries above tmin i = int(t_tmp) + 1 ! Corresponding index. weight = t_tmp - aint(t_tmp, r8) ! Fractional part of t_tmp (for interpolation). es = (1._r8 - weight)*estbl(i) + weight*estbl(i+1) end function estblf ! Get enthalpy based only on temperature ! and specific humidity. elemental function tq_enthalpy(t, q, hltalt) result(enthalpy) real(r8), intent(in) :: t ! Temperature real(r8), intent(in) :: q ! Specific humidity real(r8), intent(in) :: hltalt ! Modified hlat for T derivatives real(r8) :: enthalpy enthalpy = cpair * t + hltalt * q end function tq_enthalpy !--------------------------------------------------------------------- ! LATENT HEAT OF VAPORIZATION CORRECTIONS !--------------------------------------------------------------------- elemental subroutine no_ip_hltalt(t, hltalt) !------------------------------------------------------------------! ! Purpose: ! ! Calculate latent heat of vaporization of pure liquid water at ! ! a given temperature. ! !------------------------------------------------------------------! ! Inputs real(r8), intent(in) :: t ! Temperature ! Outputs real(r8), intent(out) :: hltalt ! Appropriately modified hlat hltalt = latvap ! Account for change of latvap with t above freezing where ! constant slope is given by -2369 j/(kg c) = cpv - cw if (t >= tmelt) then hltalt = hltalt - 2369.0_r8*(t-tmelt) end if end subroutine no_ip_hltalt elemental subroutine calc_hltalt(t, hltalt, tterm) !------------------------------------------------------------------! ! Purpose: ! ! Calculate latent heat of vaporization of water at a given ! ! temperature, taking into account the ice phase if temperature ! ! is below freezing. ! ! Optional argument also calculates a term used to calculate ! ! d(es)/dT within the water-ice transition range. ! !------------------------------------------------------------------! ! Inputs real(r8), intent(in) :: t ! Temperature ! Outputs real(r8), intent(out) :: hltalt ! Appropriately modified hlat ! Term to account for d(es)/dT in transition region. real(r8), intent(out), optional :: tterm ! Local variables real(r8) :: tc ! Temperature in degrees C real(r8) :: weight ! Weight for es transition from water to ice ! Loop iterator integer :: i if (present(tterm)) tterm = 0.0_r8 call no_ip_hltalt(t,hltalt) if (t < tmelt) then ! Weighting of hlat accounts for transition from water to ice. tc = t - tmelt if (tc >= -ttrice) then weight = -tc/ttrice ! polynomial expression approximates difference between es ! over water and es over ice from 0 to -ttrice (C) (max of ! ttrice is 40): required for accurate estimate of es ! derivative in transition range from ice to water if (present(tterm)) then do i = size(pcf), 1, -1 tterm = pcf(i) + tc*tterm end do tterm = tterm/ttrice end if else weight = 1.0_r8 end if hltalt = hltalt + weight*latice end if end subroutine calc_hltalt !--------------------------------------------------------------------- ! OPTIONAL OUTPUTS !--------------------------------------------------------------------- ! Temperature derivative outputs, for qsat_* elemental subroutine deriv_outputs(t, p, es, qs, hltalt, tterm, & gam, dqsdt) ! Inputs real(r8), intent(in) :: t ! Temperature real(r8), intent(in) :: p ! Pressure real(r8), intent(in) :: es ! Saturation vapor pressure real(r8), intent(in) :: qs ! Saturation specific humidity real(r8), intent(in) :: hltalt ! Modified latent heat real(r8), intent(in) :: tterm ! Extra term for d(es)/dT in ! transition region. ! Outputs real(r8), intent(out), optional :: gam ! (hltalt/cpair)*(d(qs)/dt) real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt) ! Local variables real(r8) :: desdt ! d(es)/dt real(r8) :: dqsdt_loc ! local copy of dqsdt if (qs == 1.0_r8) then dqsdt_loc = 0._r8 else desdt = hltalt*es/(rh2o*t*t) + tterm dqsdt_loc = qs*p*desdt/(es*(p-omeps*es)) end if if (present(dqsdt)) dqsdt = dqsdt_loc if (present(gam)) gam = dqsdt_loc * (hltalt/cpair) end subroutine deriv_outputs !--------------------------------------------------------------------- ! QSAT (SPECIFIC HUMIDITY) PROCEDURES !--------------------------------------------------------------------- elemental subroutine qsat(t, p, es, qs, gam, dqsdt, enthalpy) !------------------------------------------------------------------! ! Purpose: ! ! Look up and return saturation vapor pressure from precomputed ! ! table, then calculate and return saturation specific humidity. ! ! Optionally return various temperature derivatives or enthalpy ! ! at saturation. ! !------------------------------------------------------------------! ! Inputs real(r8), intent(in) :: t ! Temperature real(r8), intent(in) :: p ! Pressure ! Outputs real(r8), intent(out) :: es ! Saturation vapor pressure real(r8), intent(out) :: qs ! Saturation specific humidity real(r8), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt) real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q ! Local variables real(r8) :: hltalt ! Modified latent heat for T derivatives real(r8) :: tterm ! Account for d(es)/dT in transition region es = estblf(t) qs = svp_to_qsat(es, p) ! Ensures returned es is consistent with limiters on qs. es = min(es, p) ! Calculate optional arguments. if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then ! "generalized" analytic expression for t derivative of es ! accurate to within 1 percent for 173.16 < t < 373.16 call calc_hltalt(t, hltalt, tterm) if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt) call deriv_outputs(t, p, es, qs, hltalt, tterm, & gam=gam, dqsdt=dqsdt) end if end subroutine qsat elemental subroutine qsat_water(t, p, es, qs, gam, dqsdt, enthalpy) !------------------------------------------------------------------! ! Purpose: ! ! Calculate SVP over water at a given temperature, and then ! ! calculate and return saturation specific humidity. ! ! Optionally return various temperature derivatives or enthalpy ! ! at saturation. ! !------------------------------------------------------------------! use wv_sat_methods, only: wv_sat_qsat_water ! Inputs real(r8), intent(in) :: t ! Temperature real(r8), intent(in) :: p ! Pressure ! Outputs real(r8), intent(out) :: es ! Saturation vapor pressure real(r8), intent(out) :: qs ! Saturation specific humidity real(r8), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt) real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q ! Local variables real(r8) :: hltalt ! Modified latent heat for T derivatives call wv_sat_qsat_water(t, p, es, qs) if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then ! "generalized" analytic expression for t derivative of es ! accurate to within 1 percent for 173.16 < t < 373.16 call no_ip_hltalt(t, hltalt) if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt) ! For pure water/ice transition term is 0. call deriv_outputs(t, p, es, qs, hltalt, 0._r8, & gam=gam, dqsdt=dqsdt) end if end subroutine qsat_water elemental subroutine qsat_ice(t, p, es, qs, gam, dqsdt, enthalpy) !------------------------------------------------------------------! ! Purpose: ! ! Calculate SVP over ice at a given temperature, and then ! ! calculate and return saturation specific humidity. ! ! Optionally return various temperature derivatives or enthalpy ! ! at saturation. ! !------------------------------------------------------------------! use wv_sat_methods, only: wv_sat_qsat_ice ! Inputs real(r8), intent(in) :: t ! Temperature real(r8), intent(in) :: p ! Pressure ! Outputs real(r8), intent(out) :: es ! Saturation vapor pressure real(r8), intent(out) :: qs ! Saturation specific humidity real(r8), intent(out), optional :: gam ! (l/cpair)*(d(qs)/dt) real(r8), intent(out), optional :: dqsdt ! (d(qs)/dt) real(r8), intent(out), optional :: enthalpy ! cpair*t + hltalt*q ! Local variables real(r8) :: hltalt ! Modified latent heat for T derivatives call wv_sat_qsat_ice(t, p, es, qs) if (present(gam) .or. present(dqsdt) .or. present(enthalpy)) then ! For pure ice, just add latent heats. hltalt = latvap + latice if (present(enthalpy)) enthalpy = tq_enthalpy(t, qs, hltalt) ! For pure water/ice transition term is 0. call deriv_outputs(t, p, es, qs, hltalt, 0._r8, & gam=gam, dqsdt=dqsdt) end if end subroutine qsat_ice !--------------------------------------------------------------------- ! FINDSP (WET BULB TEMPERATURE) PROCEDURES !--------------------------------------------------------------------- subroutine findsp_vc(q, t, p, use_ice, tsp, qsp) use cam_logfile, only: iulog use cam_abortutils, only: endrun ! Wrapper for findsp which is 1D and handles the output status. ! Changing findsp to elemental restricted debugging output. ! If that output is needed again, it's preferable *not* to copy findsp, ! but to change the existing version. ! input arguments real(r8), intent(in) :: q(:) ! water vapor (kg/kg) real(r8), intent(in) :: t(:) ! temperature (K) real(r8), intent(in) :: p(:) ! pressure (Pa) logical, intent(in) :: use_ice ! flag to include ice phase in calculations ! output arguments real(r8), intent(out) :: tsp(:) ! saturation temp (K) real(r8), intent(out) :: qsp(:) ! saturation mixing ratio (kg/kg) integer :: status(size(q)) ! flag representing state of output ! 0 => Successful convergence ! 1 => No calculation done: pressure or specific ! humidity not within usable range ! 2 => Run failed to converge ! 4 => Temperature fell below minimum ! 8 => Enthalpy not conserved integer :: n, i n = size(q) call findsp(q, t, p, use_ice, tsp, qsp, status) ! Currently, only 2 and 8 seem to be treated as fatal errors. do i = 1,n if (status(i) == 2) then write(iulog,*) ' findsp not converging at i = ', i write(iulog,*) ' t, q, p ', t(i), q(i), p(i) write(iulog,*) ' tsp, qsp ', tsp(i), qsp(i) call endrun ('wv_saturation::FINDSP -- not converging') else if (status(i) == 8) then write(iulog,*) ' the enthalpy is not conserved at i = ', i write(iulog,*) ' t, q, p ', t(i), q(i), p(i) write(iulog,*) ' tsp, qsp ', tsp(i), qsp(i) call endrun ('wv_saturation::FINDSP -- enthalpy is not conserved') endif end do end subroutine findsp_vc elemental subroutine findsp (q, t, p, use_ice, tsp, qsp, status) !----------------------------------------------------------------------- ! ! Purpose: ! find the wet bulb temperature for a given t and q ! in a longitude height section ! wet bulb temp is the temperature and spec humidity that is ! just saturated and has the same enthalpy ! if q > qs(t) then tsp > t and qsp = qs(tsp) < q ! if q < qs(t) then tsp < t and qsp = qs(tsp) > q ! ! Method: ! a Newton method is used ! first guess uses an algorithm provided by John Petch from the UKMO ! we exclude points where the physical situation is unrealistic ! e.g. where the temperature is outside the range of validity for the ! saturation vapor pressure, or where the water vapor pressure ! exceeds the ambient pressure, or the saturation specific humidity is ! unrealistic ! ! Author: P. Rasch ! !----------------------------------------------------------------------- ! ! input arguments ! real(r8), intent(in) :: q ! water vapor (kg/kg) real(r8), intent(in) :: t ! temperature (K) real(r8), intent(in) :: p ! pressure (Pa) logical, intent(in) :: use_ice ! flag to include ice phase in calculations ! ! output arguments ! real(r8), intent(out) :: tsp ! saturation temp (K) real(r8), intent(out) :: qsp ! saturation mixing ratio (kg/kg) integer, intent(out) :: status ! flag representing state of output ! 0 => Successful convergence ! 1 => No calculation done: pressure or specific ! humidity not within usable range ! 2 => Run failed to converge ! 4 => Temperature fell below minimum ! 8 => Enthalpy not conserved ! ! local variables ! integer, parameter :: iter = 8 ! max number of times to iterate the calculation integer :: l ! iterator real(r8) es ! sat. vapor pressure real(r8) gam ! change in sat spec. hum. wrt temperature (times hltalt/cpair) real(r8) dgdt ! work variable real(r8) g ! work variable real(r8) hltalt ! lat. heat. of vap. real(r8) qs ! spec. hum. of water vapor ! work variables real(r8) t1, q1, dt, dq real(r8) qvd real(r8) r1b, c1, c2 real(r8), parameter :: dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration real(r8), parameter :: dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration real(r8) enin, enout ! Saturation specific humidity at this temperature if (use_ice) then call qsat(t, p, es, qs) else call qsat_water(t, p, es, qs) end if ! make sure a meaningful calculation is possible if (p <= 5._r8*es .or. qs <= 0._r8 .or. qs >= 0.5_r8 & .or. t < tmin .or. t > tmax) then status = 1 ! Keep initial parameters when conditions aren't suitable tsp = t qsp = q enin = 1._r8 enout = 1._r8 return end if ! Prepare to iterate status = 2 ! Get initial enthalpy if (use_ice) then call calc_hltalt(t,hltalt) else call no_ip_hltalt(t,hltalt) end if enin = tq_enthalpy(t, q, hltalt) ! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch) c1 = hltalt*c3 c2 = (t + 36._r8)**2 r1b = c2/(c2 + c1*qs) qvd = r1b * (q - qs) tsp = t + ((hltalt/cpair)*qvd) ! Generate qsp, gam, and enout from tsp. if (use_ice) then call qsat(tsp, p, es, qsp, gam=gam, enthalpy=enout) else call qsat_water(tsp, p, es, qsp, gam=gam, enthalpy=enout) end if ! iterate on first guess do l = 1, iter g = enin - enout dgdt = -cpair * (1 + gam) ! New tsp t1 = tsp - g/dgdt dt = abs(t1 - tsp)/t1 tsp = t1 ! bail out if past end of temperature range if ( tsp < tmin ) then tsp = tmin ! Get latent heat and set qsp to a value ! that preserves enthalpy. if (use_ice) then call calc_hltalt(tsp,hltalt) else call no_ip_hltalt(tsp,hltalt) end if qsp = (enin - cpair*tsp)/hltalt enout = tq_enthalpy(tsp, qsp, hltalt) status = 4 exit end if ! Re-generate qsp, gam, and enout from new tsp. if (use_ice) then call qsat(tsp, p, es, q1, gam=gam, enthalpy=enout) else call qsat_water(tsp, p, es, q1, gam=gam, enthalpy=enout) end if dq = abs(q1 - qsp)/max(q1,1.e-12_r8) qsp = q1 ! if converged at this point, exclude it from more iterations if (dt < dttol .and. dq < dqtol) then status = 0 exit endif end do ! Test for enthalpy conservation if (abs((enin-enout)/(enin+enout)) > 1.e-4_r8) status = 8 end subroutine findsp end module wv_saturation