!! Common block and statement functions for saturation vapor pressure! look-up procedure, J. J. Hack, February 1990!! $Id: wv_saturation.F90,v 1.1.4.1 2001/12/18 19:40:05 eaton Exp $!module wv_saturation  use precision, only: r8  private  save!! Public interfaces!  public gestbl   ! Initialization subroutine  public estblf   ! saturation pressure table lookup  public aqsat    ! Returns saturation vapor pressure  public aqsatd   ! Same as aqsat, but also returns a temperature derivitive  public vqsatd   ! Vector version of aqsatd!! Data used by cldwat!  public hlatv, tmin, hlatf, rgasv, pcf, cp, epsqs, ttrice!! Data!  integer plenest  ! length of saturation vapor pressure table  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)      ! table values of saturation vapor pressure  real(r8) tmin       ! min temperature (K) for table  real(r8) tmax       ! max temperature (K) for table  real(r8) ttrice     ! transition range from es over H2O to es over ice  real(r8) pcf(6)     ! polynomial coeffs -> es transition water to ice  real(r8) epsqs      ! Ratio of h2o to dry air molecular weights   real(r8) rgasv      ! Gas constant for water vapor  real(r8) hlatf      ! Latent heat of vaporization  real(r8) hlatv      ! Latent heat of fusion  real(r8) cp         ! specific heat of dry air  real(r8) tmelt      ! Melting point of water (K)  logical icephs  ! false => saturation vapor press over water onlycontains   real(r8) function estblf( td )!! Saturation vapor pressure table lookup!   real(r8), intent(in) :: td         ! Temperature for saturation lookup!   real(r8) e       ! intermediate variable for es look-up! partial pressure statement function   e(td) = max(min(td,tmax),tmin)!   estblf = (tmin+int(e(td)-tmin)-e(td)+1.)* &            estbl(int(e(td)-tmin)+1)-(tmin+int(e(td)-tmin)-e(td))* &            estbl(int(e(td)-tmin)+2)   end function estblfsubroutine 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! !-----------------------------------------------------------------------   use pmgrid, only: masterproc!------------------------------Arguments--------------------------------!! Input arguments!   real(r8), intent(in) :: tmn           ! Minimum temperature entry in es lookup table   real(r8), intent(in) :: tmx           ! Maximum temperature entry in es lookup table   real(r8), intent(in) :: epsil         ! Ratio of h2o to dry air molecular weights   real(r8), intent(in) :: trice         ! Transition range from es over range to es over ice   real(r8), intent(in) :: latvap        ! Latent heat of vaporization   real(r8), intent(in) :: latice        ! Latent heat of fusion   real(r8), intent(in) :: rh2o          ! Gas constant for water vapor   real(r8), intent(in) :: cpair         ! Specific heat of dry air   real(r8), intent(in) :: tmeltx        ! Melting point of water (K)!!---------------------------Local variables-----------------------------!   real(r8) t             ! Temperature   integer n          ! Increment counter   integer lentbl     ! Calculated length of lookup table   integer itype      ! Ice phase: 0 -> no ice phase!            1 -> ice phase, no transition!           -x -> ice phase, x degree transition   logical ip         ! Ice phase logical flag!!-----------------------------------------------------------------------!! Set es table parameters!   tmin   = tmn       ! Minimum temperature entry in table   tmax   = tmx       ! Maximum temperature entry in table   ttrice = trice     ! Trans. range from es over h2o to es over ice   icephs = ip        ! Ice phase (true or false)!! 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)   if (lentbl .gt. plenest) then      write(6,9000) tmax, tmin, plenest      call endrun    ! Abnormal termination   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) then         itype = -ttrice      else         itype = 1      end if   else      itype = 0   end if!   t = tmin - 1.0   do n=1,lentbl      t = t + 1.0      call gffgch(t,estbl(n),itype)   end do!   do n=lentbl+1,plenest      estbl(n) = -99999.0   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   pcf(2) = -5.47288442819e+00   pcf(3) = -3.67471858735e-01   pcf(4) = -8.95963532403e-03   pcf(5) = -7.78053686625e-05!!                  --- 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!   if (masterproc) then      write(6,*)' *** SATURATION VAPOR PRESSURE TABLE COMPLETED ***'   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 gestblsubroutine 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             ! I dimension of arrays t, p, es, qs   integer, intent(in) :: kk             ! K dimension of arrays t, p, es, qs   integer, intent(in) :: ilen           ! Length of vectors in I direction which   integer, intent(in) :: kstart         ! Starting location in K direction   integer, intent(in) :: kend           ! Ending location in K direction   real(r8), intent(in) :: t(ii,kk)          ! Temperature   real(r8), intent(in) :: p(ii,kk)          ! Pressure!! Output arguments!   real(r8), intent(out) :: es(ii,kk)         ! Saturation vapor pressure   real(r8), intent(out) :: qs(ii,kk)         ! Saturation specific humidity!!---------------------------Local workspace-----------------------------!   real(r8) omeps             ! 1 - 0.622   integer i, k           ! Indices!!-----------------------------------------------------------------------!   omeps = 1.0 - epsqs   do k=kstart,kend      do i=1,ilen         es(i,k) = estblf(t(i,k))!! Saturation specific humidity!         qs(i,k) = 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!         qs(i,k) = min(1.0_r8,qs(i,k))!         if (qs(i,k) < 0.0) then            qs(i,k) = 1.0            es(i,k) = p(i,k)         end if      end do   end do!   returnend subroutine aqsatsubroutine 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).!