phys/module_cam_wv_saturation.F

   1 #define WRF_PORT
   2 #define MODAL_AERO
   3 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
   4 !
   5 ! Common block and statement functions for saturation vapor pressure
   6 ! look-up procedure, J. J. Hack, February 1990
   7 !
   8 ! Ported to WRF by William.Gustafson@pnl.gov, Nov. 2009
   9 ! Updated to version from CESM1_0_1, Nov. 2010
  10 ! Updated to version CESM1.0.3 (CAM5.1.01)- Balwinder.Singh@pnnl.gov
  11 !
  12 ! $Id$
  13 !
  14 module wv_saturation
  15   use shr_kind_mod, only: r8 => shr_kind_r8
  16 #ifndef WRF_PORT
  17   use abortutils,   only: endrun
  18   use cam_logfile,  only: iulog
  19 #else
  20   use module_cam_support, only: endrun, iulog
  21   use module_wrf_error
  22 #endif
  23
  24   implicit none
  25   private
  26   save
  27 !
  28 ! Public interfaces
  29 !
  30   public gestbl   ! Initialization subroutine
  31   public estblf   ! saturation pressure table lookup
  32   public aqsat    ! Returns saturation vapor pressure
  33 #ifndef WRF_PORT
  34   public aqsatd   ! Same as aqsat, but also returns a temperature derivitive
  35 #endif
  36   public vqsatd   ! Vector version of aqsatd
  37   public fqsatd   ! Function version of vqsatd
  38   public qsat_water  ! saturation mixing ration with respect to liquid water
  39 #ifndef WRF_PORT
  40   public vqsat_water ! vector version of qsat_water
  41   public qsat_ice    ! saturation mixing ration with respect to ice
  42   public vqsat_ice   ! vector version of qsat_ice
  43 #endif
  44   public vqsatd_water
  45 #ifndef WRF_PORT
  46   public aqsat_water
  47   public vqsatd2_water         ! Variant of vqsatd_water to print out dqsdT
  48   public vqsatd2_water_single  ! Single value version of vqsatd2_water
  49 #endif
  50   public vqsatd2
  51   public vqsatd2_single
  52   public polysvp
  53 !
  54 ! Data used by cldwat
  55 !
  56   public hlatv, tmin, hlatf, rgasv, pcf, cp, epsqs, ttrice
  57 !
  58 ! Data
  59 !
  60   integer plenest  ! length of saturation vapor pressure table
  61   parameter (plenest=250)
  62 !
  63 ! Table of saturation vapor pressure values es from tmin degrees
  64 ! to tmax+1 degrees k in one degree increments.  ttrice defines the
  65 ! transition region where es is a combination of ice & water values
  66 !
  67   real(r8) estbl(plenest)      ! table values of saturation vapor pressure
  68   real(r8) tmin       ! min temperature (K) for table
  69   real(r8) tmax       ! max temperature (K) for table
  70   real(r8) ttrice     ! transition range from es over H2O to es over ice
  71   real(r8) pcf(6)     ! polynomial coeffs -> es transition water to ice
  72   real(r8) epsqs      ! Ratio of h2o to dry air molecular weights
  73   real(r8) rgasv      ! Gas constant for water vapor
  74   real(r8) hlatf      ! Latent heat of vaporization
  75   real(r8) hlatv      ! Latent heat of fusion
  76   real(r8) cp         ! specific heat of dry air
  77   real(r8) tmelt      ! Melting point of water (K)
  78   logical icephs  ! false => saturation vapor press over water only
  79
  80 contains
  81
  82    real(r8) function estblf( td )
  83 !
  84 ! Saturation vapor pressure table lookup
  85 !
  86    real(r8), intent(in) :: td         ! Temperature for saturation lookup
  87 !
  88    real(r8) :: e       ! intermediate variable for es look-up
  89    real(r8) :: ai
  90    integer  :: i
  91 !
  92    e = max(min(td,tmax),tmin)   ! partial pressure
  93    i = int(e-tmin)+1
  94    ai = aint(e-tmin)
  95    estblf = (tmin+ai-e+1._r8)* &
  96             estbl(i)-(tmin+ai-e)* &
  97             estbl(i+1)
  98    end function estblf
  99
 100 subroutine gestbl(tmn     ,tmx     ,trice   ,ip      ,epsil   , &
 101                   latvap  ,latice  ,rh2o    ,cpair   ,tmeltx   )
 102 !-----------------------------------------------------------------------
 103 !
 104 ! Purpose:
 105 ! Builds saturation vapor pressure table for later lookup procedure.
 106 !
 107 ! Method:
 108 ! Uses Goff & Gratch (1946) relationships to generate the table
 109 ! according to a set of free parameters defined below.  Auxiliary
 110 ! routines are also included for making rapid estimates (well with 1%)
 111 ! of both es and d(es)/dt for the particular table configuration.
 112 !
 113 ! Author: J. Hack
 114 !
 115 !-----------------------------------------------------------------------
 116 #ifndef WRF_PORT
 117    use spmd_utils, only: masterproc
 118 #else
 119    use module_cam_support, only: masterproc
 120    use module_cam_gffgch
 121 #endif
 122 !------------------------------Arguments--------------------------------
 123 !
 124 ! Input arguments
 125 !
 126    real(r8), intent(in) :: tmn           ! Minimum temperature entry in es lookup table
 127    real(r8), intent(in) :: tmx           ! Maximum temperature entry in es lookup table
 128    real(r8), intent(in) :: epsil         ! Ratio of h2o to dry air molecular weights
 129    real(r8), intent(in) :: trice         ! Transition range from es over range to es over ice
 130    real(r8), intent(in) :: latvap        ! Latent heat of vaporization
 131    real(r8), intent(in) :: latice        ! Latent heat of fusion
 132    real(r8), intent(in) :: rh2o          ! Gas constant for water vapor
 133    real(r8), intent(in) :: cpair         ! Specific heat of dry air
 134    real(r8), intent(in) :: tmeltx        ! Melting point of water (K)
 135 !
 136 !---------------------------Local variables-----------------------------
 137 !
 138    real(r8) t             ! Temperature
 139    integer n          ! Increment counter
 140    integer lentbl     ! Calculated length of lookup table
 141    integer itype      ! Ice phase: 0 -> no ice phase
 142 !            1 -> ice phase, no transition
 143 !           -x -> ice phase, x degree transition
 144    logical ip         ! Ice phase logical flag
 145 !
 146 !-----------------------------------------------------------------------
 147 !
 148 ! Set es table parameters
 149 !
 150    tmin   = tmn       ! Minimum temperature entry in table
 151    tmax   = tmx       ! Maximum temperature entry in table
 152    ttrice = trice     ! Trans. range from es over h2o to es over ice
 153    icephs = ip        ! Ice phase (true or false)
 154 !
 155 ! Set physical constants required for es calculation
 156 !
 157    epsqs  = epsil
 158    hlatv  = latvap
 159    hlatf  = latice
 160    rgasv  = rh2o
 161    cp     = cpair
 162    tmelt  = tmeltx
 163 !
 164    lentbl = INT(tmax-tmin+2.000001_r8)
 165    if (lentbl .gt. plenest) then
 166       write(iulog,9000) tmax, tmin, plenest
 167 #ifdef WRF_PORT
 168       call wrf_message(iulog)
 169 #endif
 170       call endrun ('GESTBL')    ! Abnormal termination
 171    end if
 172 !
 173 ! Begin building es table.
 174 ! Check whether ice phase requested.
 175 ! If so, set appropriate transition range for temperature
 176 !
 177    if (icephs) then
 178       if (ttrice /= 0.0_r8) then
 179          itype = -ttrice
 180       else
 181          itype = 1
 182       end if
 183    else
 184       itype = 0
 185    end if
 186 !
 187    t = tmin - 1.0_r8
 188    do n=1,lentbl
 189       t = t + 1.0_r8
 190       call gffgch(t,estbl(n),itype)
 191    end do
 192 !
 193    do n=lentbl+1,plenest
 194       estbl(n) = -99999.0_r8
 195    end do
 196 !
 197 ! Table complete -- Set coefficients for polynomial approximation of
 198 ! difference between saturation vapor press over water and saturation
 199 ! pressure over ice for -ttrice < t < 0 (degrees C). NOTE: polynomial
 200 ! is valid in the range -40 < t < 0 (degrees C).
 201 !
 202 !                  --- Degree 5 approximation ---
 203 !
 204    pcf(1) =  5.04469588506e-01_r8
 205    pcf(2) = -5.47288442819e+00_r8
 206    pcf(3) = -3.67471858735e-01_r8
 207    pcf(4) = -8.95963532403e-03_r8
 208    pcf(5) = -7.78053686625e-05_r8
 209 !
 210 !                  --- Degree 6 approximation ---
 211 !
 212 !-----pcf(1) =  7.63285250063e-02
 213 !-----pcf(2) = -5.86048427932e+00
 214 !-----pcf(3) = -4.38660831780e-01
 215 !-----pcf(4) = -1.37898276415e-02
 216 !-----pcf(5) = -2.14444472424e-04
 217 !-----pcf(6) = -1.36639103771e-06
 218 !
 219    if (masterproc) then
 220       write(iulog,*)' *** SATURATION VAPOR PRESSURE TABLE COMPLETED ***'
 221 #ifdef WRF_PORT
 222       call wrf_message(iulog)
 223 #endif
 224    end if
 225
 226    return
 227 !
 228 9000 format('GESTBL: FATAL ERROR *********************************',/, &
 229             ' TMAX AND TMIN REQUIRE A LARGER DIMENSION ON THE LENGTH', &
 230             ' OF THE SATURATION VAPOR PRESSURE TABLE ESTBL(PLENEST)',/, &
 231             ' TMAX, TMIN, AND PLENEST => ', 2f7.2, i3)
 232 !
 233 end subroutine gestbl
 234
 235 subroutine aqsat(t       ,p       ,es      ,qs        ,ii      , &
 236                  ilen    ,kk      ,kstart  ,kend      )
 237 !-----------------------------------------------------------------------
 238 !
 239 ! Purpose:
 240 ! Utility procedure to look up and return saturation vapor pressure from
 241 ! precomputed table, calculate and return saturation specific humidity
 242 ! (g/g),for input arrays of temperature and pressure (dimensioned ii,kk)
 243 ! This routine is useful for evaluating only a selected region in the
 244 ! vertical.
 245 !
 246 ! Method:
 247 ! <Describe the algorithm(s) used in the routine.>
 248 ! <Also include any applicable external references.>
 249 !
 250 ! Author: J. Hack
 251 !
 252 !------------------------------Arguments--------------------------------
 253 !
 254 ! Input arguments
 255 !
 256    integer, intent(in) :: ii             ! I dimension of arrays t, p, es, qs
 257    integer, intent(in) :: kk             ! K dimension of arrays t, p, es, qs
 258    integer, intent(in) :: ilen           ! Length of vectors in I direction which
 259    integer, intent(in) :: kstart         ! Starting location in K direction
 260    integer, intent(in) :: kend           ! Ending location in K direction
 261    real(r8), intent(in) :: t(ii,kk)          ! Temperature
 262    real(r8), intent(in) :: p(ii,kk)          ! Pressure
 263 !
 264 ! Output arguments
 265 !
 266    real(r8), intent(out) :: es(ii,kk)         ! Saturation vapor pressure
 267    real(r8), intent(out) :: qs(ii,kk)         ! Saturation specific humidity
 268 !
 269 !---------------------------Local workspace-----------------------------
 270 !
 271    real(r8) omeps             ! 1 - 0.622
 272    integer i, k           ! Indices
 273 !
 274 !-----------------------------------------------------------------------
 275 !
 276    omeps = 1.0_r8 - epsqs
 277    do k=kstart,kend
 278       do i=1,ilen
 279          es(i,k) = estblf(t(i,k))
 280 !
 281 ! Saturation specific humidity
 282 !
 283          qs(i,k) = epsqs*es(i,k)/(p(i,k) - omeps*es(i,k))
 284 !
 285 ! The following check is to avoid the generation of negative values
 286 ! that can occur in the upper stratosphere and mesosphere
 287 !
 288          qs(i,k) = min(1.0_r8,qs(i,k))
 289 !
 290          if (qs(i,k) < 0.0_r8) then
 291             qs(i,k) = 1.0_r8
 292             es(i,k) = p(i,k)
 293          end if
 294       end do
 295    end do
 296 !
 297    return
 298 end subroutine aqsat
 299
 300 !++xl
 301 #ifndef WRF_PORT
 302 subroutine aqsat_water(t       ,p       ,es      ,qs        ,ii      , &
 303                        ilen    ,kk      ,kstart  ,kend      )
 304 !-----------------------------------------------------------------------
 305 !
 306 ! Purpose:
 307 ! Utility procedure to look up and return saturation vapor pressure from
 308 ! precomputed table, calculate and return saturation specific humidity
 309 ! (g/g),for input arrays of temperature and pressure (dimensioned ii,kk)
 310 ! This routine is useful for evaluating only a selected region in the
 311 ! vertical.
 312 !
 313 ! Method:
 314 ! <Describe the algorithm(s) used in the routine.>
 315 ! <Also include any applicable external references.>
 316 !
 317 ! Author: J. Hack
 318 !
 319 !------------------------------Arguments--------------------------------
 320 !
 321 ! Input arguments
 322 !
 323    integer, intent(in) :: ii             ! I dimension of arrays t, p, es, qs
 324    integer, intent(in) :: kk             ! K dimension of arrays t, p, es, qs
 325    integer, intent(in) :: ilen           ! Length of vectors in I direction which
 326    integer, intent(in) :: kstart         ! Starting location in K direction
 327    integer, intent(in) :: kend           ! Ending location in K direction
 328    real(r8), intent(in) :: t(ii,kk)          ! Temperature
 329    real(r8), intent(in) :: p(ii,kk)          ! Pressure
 330 !
 331 ! Output arguments
 332 !
 333    real(r8), intent(out) :: es(ii,kk)         ! Saturation vapor pressure
 334    real(r8), intent(out) :: qs(ii,kk)         ! Saturation specific humidity
 335 !
 336 !---------------------------Local workspace-----------------------------
 337 !
 338    real(r8) omeps             ! 1 - 0.622
 339    integer i, k           ! Indices
 340 !
 341 !-----------------------------------------------------------------------
 342 !
 343    omeps = 1.0_r8 - epsqs
 344    do k=kstart,kend
 345       do i=1,ilen
 346 !        es(i,k) = estblf(t(i,k))
 347          es(i,k) = polysvp(t(i,k),0)
 348 !
 349 ! Saturation specific humidity
 350 !
 351          qs(i,k) = epsqs*es(i,k)/(p(i,k) - omeps*es(i,k))
 352 !
 353 ! The following check is to avoid the generation of negative values
 354 ! that can occur in the upper stratosphere and mesosphere
 355 !
 356          qs(i,k) = min(1.0_r8,qs(i,k))
 357 !
 358          if (qs(i,k) < 0.0_r8) then
 359             qs(i,k) = 1.0_r8
 360             es(i,k) = p(i,k)
 361          end if
 362       end do
 363    end do
 364 !
 365    return
 366 end subroutine aqsat_water
 367 !--xl
 368
 369
 370 subroutine aqsatd(t       ,p       ,es      ,qs      ,gam     , &
 371                   ii      ,ilen    ,kk      ,kstart  ,kend    )
 372 !-----------------------------------------------------------------------
 373 !
 374 ! Purpose:
 375 ! Utility procedure to look up and return saturation vapor pressure from
 376 ! precomputed table, calculate and return saturation specific humidity
 377 ! (g/g).
 378 !
 379 ! Method:
 380 ! Differs from aqsat by also calculating and returning
 381 ! gamma (l/cp)*(d(qsat)/dT)
 382 ! Input arrays temperature and pressure (dimensioned ii,kk).
 383 !
 384 ! Author: J. Hack
 385 !
 386 !------------------------------Arguments--------------------------------
 387 !
 388 ! Input arguments
 389 !
 390    integer, intent(in) :: ii            ! I dimension of arrays t, p, es, qs
 391    integer, intent(in) :: ilen          ! Vector length in I direction
 392    integer, intent(in) :: kk            ! K dimension of arrays t, p, es, qs
 393    integer, intent(in) :: kstart        ! Starting location in K direction
 394    integer, intent(in) :: kend          ! Ending location in K direction
 395
 396    real(r8), intent(in) :: t(ii,kk)         ! Temperature
 397    real(r8), intent(in) :: p(ii,kk)         ! Pressure
 398
 399 !
 400 ! Output arguments
 401 !
 402    real(r8), intent(out) :: es(ii,kk)        ! Saturation vapor pressure
 403    real(r8), intent(out) :: qs(ii,kk)        ! Saturation specific humidity
 404    real(r8), intent(out) :: gam(ii,kk)       ! (l/cp)*(d(qs)/dt)
 405 !
 406 !---------------------------Local workspace-----------------------------
 407 !
 408    logical lflg          ! True if in temperature transition region
 409    integer i             ! i index for vector calculations
 410    integer k             ! k index
 411    real(r8) omeps            ! 1. - 0.622
 412    real(r8) trinv            ! Reciprocal of ttrice (transition range)
 413    real(r8) tc               ! Temperature (in degrees C)
 414    real(r8) weight           ! Weight for es transition from water to ice
 415    real(r8) hltalt           ! Appropriately modified hlat for T derivatives
 416    real(r8) hlatsb           ! hlat weighted in transition region
 417    real(r8) hlatvp           ! hlat modified for t changes above freezing
 418    real(r8) tterm            ! Account for d(es)/dT in transition region
 419    real(r8) desdt            ! d(es)/dT
 420 !
 421 !-----------------------------------------------------------------------
 422 !
 423    omeps = 1.0_r8 - epsqs
 424    do k=kstart,kend
 425       do i=1,ilen
 426          es(i,k) = estblf(t(i,k))
 427 !
 428 ! Saturation specific humidity
 429 !
 430          qs(i,k) = epsqs*es(i,k)/(p(i,k) - omeps*es(i,k))
 431 !
 432 ! The following check is to avoid the generation of negative qs
 433 ! values which can occur in the upper stratosphere and mesosphere
 434 !
 435          qs(i,k) = min(1.0_r8,qs(i,k))
 436 !
 437          if (qs(i,k) < 0.0_r8) then
 438             qs(i,k) = 1.0_r8
 439             es(i,k) = p(i,k)
 440          end if
 441       end do
 442    end do
 443 !
 444 ! "generalized" analytic expression for t derivative of es
 445 ! accurate to within 1 percent for 173.16 < t < 373.16
 446 !
 447    trinv = 0.0_r8
 448    if ((.not. icephs) .or. (ttrice.eq.0.0_r8)) go to 10
 449    trinv = 1.0_r8/ttrice
 450 !
 451    do k=kstart,kend
 452       do i=1,ilen
 453 !
 454 ! Weighting of hlat accounts for transition from water to ice
 455 ! polynomial expression approximates difference between es over
 456 ! water and es over ice from 0 to -ttrice (C) (min of ttrice is
 457 ! -40): required for accurate estimate of es derivative in transition
 458 ! range from ice to water also accounting for change of hlatv with t
 459 ! above freezing where constant slope is given by -2369 j/(kg c) =cpv - cw
 460 !
 461          tc     = t(i,k) - tmelt
 462          lflg   = (tc >= -ttrice .and. tc < 0.0_r8)
 463          weight = min(-tc*trinv,1.0_r8)
 464          hlatsb = hlatv + weight*hlatf
 465          hlatvp = hlatv - 2369.0_r8*tc
 466          if (t(i,k) < tmelt) then
 467             hltalt = hlatsb
 468          else
 469             hltalt = hlatvp
 470          end if
 471          if (lflg) then
 472             tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) + tc*pcf(5))))
 473          else
 474             tterm = 0.0_r8
 475          end if
 476          desdt    = hltalt*es(i,k)/(rgasv*t(i,k)*t(i,k)) + tterm*trinv
 477          gam(i,k) = hltalt*qs(i,k)*p(i,k)*desdt/(cp*es(i,k)*(p(i,k) - omeps*es(i,k)))
 478          if (qs(i,k) == 1.0_r8) gam(i,k) = 0.0_r8
 479       end do
 480    end do
 481 !
 482    go to 20
 483 !
 484 ! No icephs or water to ice transition
 485 !
 486 10 do k=kstart,kend
 487       do i=1,ilen
 488 !
 489 ! Account for change of hlatv with t above freezing where
 490 ! constant slope is given by -2369 j/(kg c) = cpv - cw
 491 !
 492          hlatvp = hlatv - 2369.0_r8*(t(i,k)-tmelt)
 493          if (icephs) then
 494             hlatsb = hlatv + hlatf
 495          else
 496             hlatsb = hlatv
 497          end if
 498          if (t(i,k) < tmelt) then
 499             hltalt = hlatsb
 500          else
 501             hltalt = hlatvp
 502          end if
 503          desdt    = hltalt*es(i,k)/(rgasv*t(i,k)*t(i,k))
 504          gam(i,k) = hltalt*qs(i,k)*p(i,k)*desdt/(cp*es(i,k)*(p(i,k) - omeps*es(i,k)))
 505          if (qs(i,k) == 1.0_r8) gam(i,k) = 0.0_r8
 506       end do
 507    end do
 508 !
 509 20 return
 510 end subroutine aqsatd
 511 #endif
 512 subroutine vqsatd(t       ,p       ,es      ,qs      ,gam      , &
 513                   len     )
 514 !-----------------------------------------------------------------------
 515 !
 516 ! Purpose:
 517 ! Utility procedure to look up and return saturation vapor pressure from
 518 ! precomputed table, calculate and return saturation specific humidity
 519 ! (g/g), and calculate and return gamma (l/cp)*(d(qsat)/dT).  The same
 520 ! function as qsatd, but operates on vectors of temperature and pressure
 521 !
 522 ! Method:
 523 !
 524 ! Author: J. Hack
 525 !
 526 !------------------------------Arguments--------------------------------
 527 !
 528 ! Input arguments
 529 !
 530    integer, intent(in) :: len       ! vector length
 531    real(r8), intent(in) :: t(len)       ! temperature
 532    real(r8), intent(in) :: p(len)       ! pressure
 533 !
 534 ! Output arguments
 535 !
 536    real(r8), intent(out) :: es(len)   ! saturation vapor pressure
 537    real(r8), intent(out) :: qs(len)   ! saturation specific humidity
 538    real(r8), intent(out) :: gam(len)  ! (l/cp)*(d(qs)/dt)
 539 !
 540 !--------------------------Local Variables------------------------------
 541 !
 542    logical lflg   ! true if in temperature transition region
 543 !
 544    integer i      ! index for vector calculations
 545 !
 546    real(r8) omeps     ! 1. - 0.622
 547    real(r8) trinv     ! reciprocal of ttrice (transition range)
 548    real(r8) tc        ! temperature (in degrees C)
 549    real(r8) weight    ! weight for es transition from water to ice
 550    real(r8) hltalt    ! appropriately modified hlat for T derivatives
 551 !
 552    real(r8) hlatsb    ! hlat weighted in transition region
 553    real(r8) hlatvp    ! hlat modified for t changes above freezing
 554    real(r8) tterm     ! account for d(es)/dT in transition region
 555    real(r8) desdt     ! d(es)/dT
 556 !
 557 !-----------------------------------------------------------------------
 558 !
 559    omeps = 1.0_r8 - epsqs
 560    do i=1,len
 561       es(i) = estblf(t(i))
 562 !
 563 ! Saturation specific humidity
 564 !
 565       qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
 566 !
 567 ! The following check is to avoid the generation of negative
 568 ! values that can occur in the upper stratosphere and mesosphere
 569 !
 570       qs(i) = min(1.0_r8,qs(i))
 571 !
 572       if (qs(i) < 0.0_r8) then
 573          qs(i) = 1.0_r8
 574          es(i) = p(i)
 575       end if
 576    end do
 577 !
 578 ! "generalized" analytic expression for t derivative of es
 579 ! accurate to within 1 percent for 173.16 < t < 373.16
 580 !
 581    trinv = 0.0_r8
 582    if ((.not. icephs) .or. (ttrice.eq.0.0_r8)) go to 10
 583    trinv = 1.0_r8/ttrice
 584    do i=1,len
 585 !
 586 ! Weighting of hlat accounts for transition from water to ice
 587 ! polynomial expression approximates difference between es over
 588 ! water and es over ice from 0 to -ttrice (C) (min of ttrice is
 589 ! -40): required for accurate estimate of es derivative in transition
 590 ! range from ice to water also accounting for change of hlatv with t
 591 ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
 592 !
 593       tc     = t(i) - tmelt
 594       lflg   = (tc >= -ttrice .and. tc < 0.0_r8)
 595       weight = min(-tc*trinv,1.0_r8)
 596       hlatsb = hlatv + weight*hlatf
 597       hlatvp = hlatv - 2369.0_r8*tc
 598       if (t(i) < tmelt) then
 599          hltalt = hlatsb
 600       else
 601          hltalt = hlatvp
 602       end if
 603       if (lflg) then
 604          tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) + tc*pcf(5))))
 605       else
 606          tterm = 0.0_r8
 607       end if
 608       desdt  = hltalt*es(i)/(rgasv*t(i)*t(i)) + tterm*trinv
 609       gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
 610       if (qs(i) == 1.0_r8) gam(i) = 0.0_r8
 611    end do
 612    return
 613 !
 614 ! No icephs or water to ice transition
 615 !
 616 10 do i=1,len
 617 !
 618 ! Account for change of hlatv with t above freezing where
 619 ! constant slope is given by -2369 j/(kg c) = cpv - cw
 620 !
 621       hlatvp = hlatv - 2369.0_r8*(t(i)-tmelt)
 622       if (icephs) then
 623          hlatsb = hlatv + hlatf
 624       else
 625          hlatsb = hlatv
 626       end if
 627       if (t(i) < tmelt) then
 628          hltalt = hlatsb
 629       else
 630          hltalt = hlatvp
 631       end if
 632       desdt  = hltalt*es(i)/(rgasv*t(i)*t(i))
 633       gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
 634       if (qs(i) == 1.0_r8) gam(i) = 0.0_r8
 635    end do
 636 !
 637    return
 638 !
 639 end subroutine vqsatd
 640
 641 !++xl
 642 subroutine vqsatd_water(t       ,p       ,es      ,qs      ,gam      , &
 643                         len     )
 644
 645 !------------------------------Arguments--------------------------------
 646 !
 647 ! Input arguments
 648 !
 649    integer, intent(in) :: len       ! vector length
 650    real(r8), intent(in) :: t(len)       ! temperature
 651    real(r8), intent(in) :: p(len)       ! pressure
 652
 653 !
 654 ! Output arguments
 655 !
 656    real(r8), intent(out) :: es(len)   ! saturation vapor pressure
 657    real(r8), intent(out) :: qs(len)   ! saturation specific humidity
 658    real(r8), intent(out) :: gam(len)  ! (l/cp)*(d(qs)/dt)
 659
 660 !
 661 !--------------------------Local Variables------------------------------
 662 !
 663 !
 664    integer i      ! index for vector calculations
 665 !
 666    real(r8) omeps     ! 1. - 0.622
 667    real(r8) hltalt    ! appropriately modified hlat for T derivatives
 668 !
 669    real(r8) hlatsb    ! hlat weighted in transition region
 670    real(r8) hlatvp    ! hlat modified for t changes above freezing
 671    real(r8) desdt     ! d(es)/dT
 672 !
 673 !-----------------------------------------------------------------------
 674 !
 675    omeps = 1.0_r8 - epsqs
 676    do i=1,len
 677       es(i) = polysvp(t(i),0)
 678 !
 679 ! Saturation specific humidity
 680 !
 681       qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
 682 !
 683 ! The following check is to avoid the generation of negative
 684 ! values that can occur in the upper stratosphere and mesosphere
 685 !
 686       qs(i) = min(1.0_r8,qs(i))
 687 !
 688       if (qs(i) < 0.0_r8) then
 689          qs(i) = 1.0_r8
 690          es(i) = p(i)
 691       end if
 692    end do
 693 !
 694 ! No icephs or water to ice transition
 695 !
 696    do i=1,len
 697 !
 698 ! Account for change of hlatv with t above freezing where
 699 ! constant slope is given by -2369 j/(kg c) = cpv - cw
 700 !
 701       hlatvp = hlatv - 2369.0_r8*(t(i)-tmelt)
 702       hlatsb = hlatv
 703       if (t(i) < tmelt) then
 704          hltalt = hlatsb
 705       else
 706          hltalt = hlatvp
 707       end if
 708       desdt  = hltalt*es(i)/(rgasv*t(i)*t(i))
 709       gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
 710       if (qs(i) == 1.0_r8) gam(i) = 0.0_r8
 711    end do
 712 !
 713    return
 714 !
 715 end subroutine vqsatd_water
 716
 717       function polysvp (T,type)
 718 !  Compute saturation vapor pressure by using
 719 ! function from Goff and Gatch (1946)
 720
 721 !  Polysvp returned in units of pa.
 722 !  T is input in units of K.
 723 !  type refers to saturation with respect to liquid (0) or ice (1)
 724
 725       real(r8) dum
 726
 727       real(r8) T,polysvp
 728
 729       integer type
 730
 731 ! ice
 732
 733       if (type.eq.1) then
 734
 735 ! Goff Gatch equation (good down to -100 C)
 736
 737          polysvp = 10._r8**(-9.09718_r8*(273.16_r8/t-1._r8)-3.56654_r8* &
 738           log10(273.16_r8/t)+0.876793_r8*(1._r8-t/273.16_r8)+ &
 739           log10(6.1071_r8))*100._r8
 740
 741       end if
 742
 743 ! Goff Gatch equation, uncertain below -70 C
 744
 745       if (type.eq.0) then
 746          polysvp = 10._r8**(-7.90298_r8*(373.16_r8/t-1._r8)+ &
 747              5.02808_r8*log10(373.16_r8/t)- &
 748              1.3816e-7_r8*(10._r8**(11.344_r8*(1._r8-t/373.16_r8))-1._r8)+ &
 749              8.1328e-3_r8*(10._r8**(-3.49149_r8*(373.16_r8/t-1._r8))-1._r8)+ &
 750              log10(1013.246_r8))*100._r8
 751          end if
 752
 753
 754       end function polysvp
 755 !--xl
 756
 757 integer function fqsatd(t    ,p    ,es    ,qs   ,gam   , len     )
 758   !-----------------------------------------------------------------------
 759   ! Purpose:
 760   ! This is merely a function interface vqsatd.
 761   !------------------------------Arguments--------------------------------
 762   ! Input arguments
 763   integer, intent(in) :: len       ! vector length
 764   real(r8), intent(in) :: t(len)       ! temperature
 765   real(r8), intent(in) :: p(len)       ! pressure
 766   ! Output arguments
 767   real(r8), intent(out) :: es(len)   ! saturation vapor pressure
 768   real(r8), intent(out) :: qs(len)   ! saturation specific humidity
 769   real(r8), intent(out) :: gam(len)  ! (l/cp)*(d(qs)/dt)
 770   ! Call vqsatd
 771   call vqsatd(t       ,p       ,es      ,qs      ,gam  , len     )
 772   fqsatd = 1
 773   return
 774 end function fqsatd
 775
 776 real(r8) function qsat_water(t,p)
 777   !  saturation mixing ratio w/respect to liquid water
 778   real(r8) t ! temperature
 779   real(r8) p ! pressure (Pa)
 780   real(r8) es ! saturation vapor pressure (Pa)
 781   real(r8) ps, ts, e1, e2, f1, f2, f3, f4, f5, f
 782   !  real(r8) t0inv ! 1/273.
 783   !  data t0inv/0.003663/
 784   !  save t0inv
 785   !  es = 611.*exp(hlatv/rgasv*(t0inv-1./t))
 786
 787   ps = 1013.246_r8
 788   ts = 373.16_r8
 789   e1 = 11.344_r8*(1.0_r8 - t/ts)
 790   e2 = -3.49149_r8*(ts/t - 1.0_r8)
 791   f1 = -7.90298_r8*(ts/t - 1.0_r8)
 792   f2 = 5.02808_r8*log10(ts/t)
 793   f3 = -1.3816_r8*(10.0_r8**e1 - 1.0_r8)/10000000.0_r8
 794   f4 = 8.1328_r8*(10.0_r8**e2 - 1.0_r8)/1000.0_r8
 795   f5 = log10(ps)
 796   f  = f1 + f2 + f3 + f4 + f5
 797   es = (10.0_r8**f)*100.0_r8
 798
 799   qsat_water = epsqs*es/(p-(1.-epsqs)*es) ! saturation w/respect to liquid only
 800   if(qsat_water < 0.) qsat_water = 1.
 801
 802 end function qsat_water
 803 #ifndef WRF_PORT
 804 subroutine vqsat_water(t,p,qsat_water,len)
 805   !  saturation mixing ratio w/respect to liquid water
 806   integer, intent(in)  :: len
 807   real(r8) t(len) ! temperature
 808   real(r8) p(len) ! pressure (Pa)
 809   real(r8) qsat_water(len)
 810   real(r8) es ! saturation vapor pressure (Pa)
 811   real(r8), parameter :: t0inv = 1._r8/273._r8
 812   real(r8) coef
 813   integer :: i
 814
 815   coef = hlatv/rgasv
 816   do i=1,len
 817      es = 611._r8*exp(coef*(t0inv-1./t(i)))
 818      qsat_water(i) = epsqs*es/(p(i)-(1.-epsqs)*es) ! saturation w/respect to liquid only
 819      if(qsat_water(i) < 0.) qsat_water(i) = 1.
 820   enddo
 821
 822   return
 823
 824 end subroutine vqsat_water
 825
 826 real(r8) function qsat_ice(t,p)
 827   !  saturation mixing ratio w/respect to ice
 828   real(r8) t ! temperature
 829   real(r8) p ! pressure (Pa)
 830   real(r8) es ! saturation vapor pressure (Pa)
 831   real(r8), parameter :: t0inv = 1._r8/273._r8
 832   es = 611.*exp((hlatv+hlatf)/rgasv*(t0inv-1./t))
 833   qsat_ice = epsqs*es/(p-(1.-epsqs)*es) ! saturation w/respect to liquid only
 834   if(qsat_ice < 0.) qsat_ice = 1.
 835
 836 end function qsat_ice
 837
 838 subroutine vqsat_ice(t,p,qsat_ice,len)
 839   !  saturation mixing ratio w/respect to liquid water
 840   integer,intent(in) :: len
 841   real(r8) t(len) ! temperature
 842   real(r8) p(len) ! pressure (Pa)
 843   real(r8) qsat_ice(len)
 844   real(r8) es ! saturation vapor pressure (Pa)
 845   real(r8), parameter :: t0inv = 1._r8/273._r8
 846   real(r8) coef
 847   integer :: i
 848
 849   coef = (hlatv+hlatf)/rgasv
 850   do i=1,len
 851      es = 611.*exp(coef*(t0inv-1./t(i)))
 852      qsat_ice(i) = epsqs*es/(p(i)-(1.-epsqs)*es) ! saturation w/respect to liquid only
 853      if(qsat_ice(i) < 0.) qsat_ice(i) = 1.
 854   enddo
 855
 856   return
 857
 858 end subroutine vqsat_ice
 859
 860 ! Sungsu
 861 ! Below two subroutines (vqsatd2_water,vqsatd2_water_single) are by Sungsu
 862 ! Replace 'gam -> dqsdt'
 863 ! Sungsu
 864
 865 subroutine vqsatd2_water(t       ,p       ,es      ,qs      ,dqsdt      , &
 866                         len     )
 867
 868 !------------------------------Arguments--------------------------------
 869 !
 870 ! Input arguments
 871 !
 872    integer, intent(in) :: len       ! vector length
 873    real(r8), intent(in) :: t(len)       ! temperature
 874    real(r8), intent(in) :: p(len)       ! pressure
 875
 876 !
 877 ! Output arguments
 878 !
 879    real(r8), intent(out) :: es(len)   ! saturation vapor pressure
 880    real(r8), intent(out) :: qs(len)   ! saturation specific humidity
 881  ! real(r8), intent(out) :: gam(len)  ! (l/cp)*(d(qs)/dt)
 882  ! Sungsu
 883    real(r8), intent(out) :: dqsdt(len)  ! (d(qs)/dt)
 884  ! End by Sungsu
 885
 886 !
 887 !--------------------------Local Variables------------------------------
 888 !
 889 !
 890    integer i      ! index for vector calculations
 891 !
 892    real(r8) omeps     ! 1. - 0.622
 893    real(r8) hltalt    ! appropriately modified hlat for T derivatives
 894 !
 895    real(r8) hlatsb    ! hlat weighted in transition region
 896    real(r8) hlatvp    ! hlat modified for t changes above freezing
 897    real(r8) desdt     ! d(es)/dT
 898
 899  ! Sungsu
 900    real(r8) gam(len)  ! (l/cp)*(d(qs)/dt)
 901  ! End by Sungsu
 902
 903 !
 904 !-----------------------------------------------------------------------
 905 !
 906    omeps = 1.0_r8 - epsqs
 907    do i=1,len
 908       es(i) = polysvp(t(i),0)
 909 !
 910 ! Saturation specific humidity
 911 !
 912       qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
 913 !
 914 ! The following check is to avoid the generation of negative
 915 ! values that can occur in the upper stratosphere and mesosphere
 916 !
 917       qs(i) = min(1.0_r8,qs(i))
 918 !
 919       if (qs(i) < 0.0_r8) then
 920          qs(i) = 1.0_r8
 921          es(i) = p(i)
 922       end if
 923    end do
 924 !
 925 ! No icephs or water to ice transition
 926 !
 927    do i=1,len
 928 !
 929 ! Account for change of hlatv with t above freezing where
 930 ! constant slope is given by -2369 j/(kg c) = cpv - cw
 931 !
 932       hlatvp = hlatv - 2369.0_r8*(t(i)-tmelt)
 933       hlatsb = hlatv
 934       if (t(i) < tmelt) then
 935          hltalt = hlatsb
 936       else
 937          hltalt = hlatvp
 938       end if
 939       desdt  = hltalt*es(i)/(rgasv*t(i)*t(i))
 940       gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
 941       if (qs(i) == 1.0_r8) gam(i) = 0.0_r8
 942     ! Sungsu
 943       dqsdt(i) = (cp/hltalt)*gam(i)
 944     ! End by Sungsu
 945    end do
 946 !
 947    return
 948 !
 949 end subroutine vqsatd2_water
 950
 951 subroutine vqsatd2_water_single(t       ,p       ,es      ,qs      ,dqsdt)
 952
 953 !------------------------------Arguments--------------------------------
 954 !
 955 ! Input arguments
 956 !
 957
 958    real(r8), intent(in) :: t       ! temperature
 959    real(r8), intent(in) :: p       ! pressure
 960
 961 !
 962 ! Output arguments
 963 !
 964    real(r8), intent(out) :: es   ! saturation vapor pressure
 965    real(r8), intent(out) :: qs   ! saturation specific humidity
 966  ! real(r8), intent(out) :: gam  ! (l/cp)*(d(qs)/dt)
 967  ! Sungsu
 968    real(r8), intent(out) :: dqsdt  ! (d(qs)/dt)
 969  ! End by Sungsu
 970
 971 !
 972 !--------------------------Local Variables------------------------------
 973 !
 974 !
 975    integer i      ! index for vector calculations
 976 !
 977    real(r8) omeps     ! 1. - 0.622
 978    real(r8) hltalt    ! appropriately modified hlat for T derivatives
 979 !
 980    real(r8) hlatsb    ! hlat weighted in transition region
 981    real(r8) hlatvp    ! hlat modified for t changes above freezing
 982    real(r8) desdt     ! d(es)/dT
 983
 984  ! Sungsu
 985    real(r8) gam  ! (l/cp)*(d(qs)/dt)
 986  ! End by Sungsu
 987
 988 !
 989 !-----------------------------------------------------------------------
 990 !
 991    omeps = 1.0_r8 - epsqs
 992 !  do i=1,len
 993       es = polysvp(t,0)
 994 !
 995 ! Saturation specific humidity
 996 !
 997       qs = epsqs*es/(p - omeps*es)
 998 !
 999 ! The following check is to avoid the generation of negative
1000 ! values that can occur in the upper stratosphere and mesosphere
1001 !
1002       qs = min(1.0_r8,qs)
1003 !
1004       if (qs < 0.0_r8) then
1005          qs = 1.0_r8
1006          es = p
1007       end if
1008 !  end do
1009 !
1010 ! No icephs or water to ice transition
1011 !
1012 !  do i=1,len
1013 !
1014 ! Account for change of hlatv with t above freezing where
1015 ! constant slope is given by -2369 j/(kg c) = cpv - cw
1016 !
1017       hlatvp = hlatv - 2369.0_r8*(t-tmelt)
1018       hlatsb = hlatv
1019       if (t < tmelt) then
1020          hltalt = hlatsb
1021       else
1022          hltalt = hlatvp
1023       end if
1024       desdt  = hltalt*es/(rgasv*t*t)
1025       gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
1026       if (qs == 1.0_r8) gam = 0.0_r8
1027     ! Sungsu
1028       dqsdt = (cp/hltalt)*gam
1029     ! End by Sungsu
1030 !  end do
1031 !
1032    return
1033 !
1034 end subroutine vqsatd2_water_single
1035
1036 #endif
1037 subroutine vqsatd2(t       ,p       ,es      ,qs      ,dqsdt      , &
1038                    len     )
1039 !-----------------------------------------------------------------------
1040 ! Sungsu : This is directly copied from 'vqsatd' but 'dqsdt' is output
1041 !          instead of gam for use in Sungsu's equilibrium stratiform
1042 !          macrophysics scheme.
1043 !
1044 ! Purpose:
1045 ! Utility procedure to look up and return saturation vapor pressure from
1046 ! precomputed table, calculate and return saturation specific humidity
1047 ! (g/g), and calculate and return gamma (l/cp)*(d(qsat)/dT).  The same
1048 ! function as qsatd, but operates on vectors of temperature and pressure
1049 !
1050 ! Method:
1051 !
1052 ! Author: J. Hack
1053 !
1054 !------------------------------Arguments--------------------------------
1055 !
1056 ! Input arguments
1057 !
1058    integer, intent(in) :: len       ! vector length
1059    real(r8), intent(in) :: t(len)       ! temperature
1060    real(r8), intent(in) :: p(len)       ! pressure
1061 !
1062 ! Output arguments
1063 !
1064    real(r8), intent(out) :: es(len)   ! saturation vapor pressure
1065    real(r8), intent(out) :: qs(len)   ! saturation specific humidity
1066  ! real(r8), intent(out) :: gam(len)  ! (l/cp)*(d(qs)/dt)
1067  ! Sungsu
1068    real(r8), intent(out) :: dqsdt(len)  ! (d(qs)/dt)
1069  ! End by Sungsu
1070
1071 !
1072 !--------------------------Local Variables------------------------------
1073 !
1074    logical lflg   ! true if in temperature transition region
1075 !
1076    integer i      ! index for vector calculations
1077 !
1078    real(r8) omeps     ! 1. - 0.622
1079    real(r8) trinv     ! reciprocal of ttrice (transition range)
1080    real(r8) tc        ! temperature (in degrees C)
1081    real(r8) weight    ! weight for es transition from water to ice
1082    real(r8) hltalt    ! appropriately modified hlat for T derivatives
1083 !
1084    real(r8) hlatsb    ! hlat weighted in transition region
1085    real(r8) hlatvp    ! hlat modified for t changes above freezing
1086    real(r8) tterm     ! account for d(es)/dT in transition region
1087    real(r8) desdt     ! d(es)/dT
1088
1089  ! Sungsu
1090    real(r8) gam(len)  ! (l/cp)*(d(qs)/dt)
1091  ! End by Sungsu
1092 !
1093 !-----------------------------------------------------------------------
1094 !
1095    omeps = 1.0_r8 - epsqs
1096    do i=1,len
1097       es(i) = estblf(t(i))
1098 !
1099 ! Saturation specific humidity
1100 !
1101       qs(i) = epsqs*es(i)/(p(i) - omeps*es(i))
1102 !
1103 ! The following check is to avoid the generation of negative
1104 ! values that can occur in the upper stratosphere and mesosphere
1105 !
1106       qs(i) = min(1.0_r8,qs(i))
1107 !
1108       if (qs(i) < 0.0_r8) then
1109          qs(i) = 1.0_r8
1110          es(i) = p(i)
1111       end if
1112    end do
1113 !
1114 ! "generalized" analytic expression for t derivative of es
1115 ! accurate to within 1 percent for 173.16 < t < 373.16
1116 !
1117    trinv = 0.0_r8
1118    if ((.not. icephs) .or. (ttrice.eq.0.0_r8)) go to 10
1119    trinv = 1.0_r8/ttrice
1120    do i=1,len
1121 !
1122 ! Weighting of hlat accounts for transition from water to ice
1123 ! polynomial expression approximates difference between es over
1124 ! water and es over ice from 0 to -ttrice (C) (min of ttrice is
1125 ! -40): required for accurate estimate of es derivative in transition
1126 ! range from ice to water also accounting for change of hlatv with t
1127 ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
1128 !
1129       tc     = t(i) - tmelt
1130       lflg   = (tc >= -ttrice .and. tc < 0.0_r8)
1131       weight = min(-tc*trinv,1.0_r8)
1132       hlatsb = hlatv + weight*hlatf
1133       hlatvp = hlatv - 2369.0_r8*tc
1134       if (t(i) < tmelt) then
1135          hltalt = hlatsb
1136       else
1137          hltalt = hlatvp
1138       end if
1139       if (lflg) then
1140          tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) + tc*pcf(5))))
1141       else
1142          tterm = 0.0_r8
1143       end if
1144       desdt  = hltalt*es(i)/(rgasv*t(i)*t(i)) + tterm*trinv
1145       gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
1146       if (qs(i) == 1.0_r8) gam(i) = 0.0_r8
1147     ! Sungsu
1148       dqsdt(i) = (cp/hltalt)*gam(i)
1149     ! End by Sungsu
1150    end do
1151    return
1152 !
1153 ! No icephs or water to ice transition
1154 !
1155 10 do i=1,len
1156 !
1157 ! Account for change of hlatv with t above freezing where
1158 ! constant slope is given by -2369 j/(kg c) = cpv - cw
1159 !
1160       hlatvp = hlatv - 2369.0_r8*(t(i)-tmelt)
1161       if (icephs) then
1162          hlatsb = hlatv + hlatf
1163       else
1164          hlatsb = hlatv
1165       end if
1166       if (t(i) < tmelt) then
1167          hltalt = hlatsb
1168       else
1169          hltalt = hlatvp
1170       end if
1171       desdt  = hltalt*es(i)/(rgasv*t(i)*t(i))
1172       gam(i) = hltalt*qs(i)*p(i)*desdt/(cp*es(i)*(p(i) - omeps*es(i)))
1173       if (qs(i) == 1.0_r8) gam(i) = 0.0_r8
1174     ! Sungsu
1175       dqsdt(i) = (cp/hltalt)*gam(i)
1176     ! End by Sungsu
1177    end do
1178 !
1179    return
1180 !
1181 end subroutine vqsatd2
1182
1183
1184 ! Below routine is by Sungsu
1185
1186 subroutine vqsatd2_single(t       ,p       ,es      ,qs      ,dqsdt)
1187 !-----------------------------------------------------------------------
1188 ! Sungsu : This is directly copied from 'vqsatd' but 'dqsdt' is output
1189 !          instead of gam for use in Sungsu's equilibrium stratiform
1190 !          macrophysics scheme.
1191 !
1192 ! Purpose:
1193 ! Utility procedure to look up and return saturation vapor pressure from
1194 ! precomputed table, calculate and return saturation specific humidity
1195 ! (g/g), and calculate and return gamma (l/cp)*(d(qsat)/dT).  The same
1196 ! function as qsatd, but operates on vectors of temperature and pressure
1197 !
1198 ! Method:
1199 !
1200 ! Author: J. Hack
1201 !
1202 !------------------------------Arguments--------------------------------
1203 !
1204 ! Input arguments
1205 !
1206    real(r8), intent(in) :: t       ! temperature
1207    real(r8), intent(in) :: p       ! pressure
1208 !
1209 ! Output arguments
1210 !
1211    real(r8), intent(out) :: es     ! saturation vapor pressure
1212    real(r8), intent(out) :: qs     ! saturation specific humidity
1213  ! real(r8), intent(out) :: gam    ! (l/cp)*(d(qs)/dt)
1214  ! Sungsu
1215    real(r8), intent(out) :: dqsdt  ! (d(qs)/dt)
1216  ! End by Sungsu
1217
1218 !
1219 !--------------------------Local Variables------------------------------
1220 !
1221    logical lflg   ! true if in temperature transition region
1222 !
1223 !  integer i      ! index for vector calculations
1224 !
1225    real(r8) omeps     ! 1. - 0.622
1226    real(r8) trinv     ! reciprocal of ttrice (transition range)
1227    real(r8) tc        ! temperature (in degrees C)
1228    real(r8) weight    ! weight for es transition from water to ice
1229    real(r8) hltalt    ! appropriately modified hlat for T derivatives
1230 !
1231    real(r8) hlatsb    ! hlat weighted in transition region
1232    real(r8) hlatvp    ! hlat modified for t changes above freezing
1233    real(r8) tterm     ! account for d(es)/dT in transition region
1234    real(r8) desdt     ! d(es)/dT
1235
1236  ! Sungsu
1237    real(r8) gam       ! (l/cp)*(d(qs)/dt)
1238  ! End by Sungsu
1239 !
1240 !-----------------------------------------------------------------------
1241 !
1242    omeps = 1.0_r8 - epsqs
1243
1244 !  do i=1,len
1245
1246       es = estblf(t)
1247 !
1248 ! Saturation specific humidity
1249 !
1250       qs = epsqs*es/(p - omeps*es)
1251 !
1252 ! The following check is to avoid the generation of negative
1253 ! values that can occur in the upper stratosphere and mesosphere
1254 !
1255       qs = min(1.0_r8,qs)
1256 !
1257       if (qs < 0.0_r8) then
1258          qs = 1.0_r8
1259          es = p
1260       end if
1261
1262 !  end do
1263 !
1264 ! "generalized" analytic expression for t derivative of es
1265 ! accurate to within 1 percent for 173.16 < t < 373.16
1266 !
1267    trinv = 0.0_r8
1268    if ((.not. icephs) .or. (ttrice.eq.0.0_r8)) go to 10
1269    trinv = 1.0_r8/ttrice
1270
1271 !  do i=1,len
1272 !
1273 ! Weighting of hlat accounts for transition from water to ice
1274 ! polynomial expression approximates difference between es over
1275 ! water and es over ice from 0 to -ttrice (C) (min of ttrice is
1276 ! -40): required for accurate estimate of es derivative in transition
1277 ! range from ice to water also accounting for change of hlatv with t
1278 ! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
1279 !
1280       tc     = t - tmelt
1281       lflg   = (tc >= -ttrice .and. tc < 0.0_r8)
1282       weight = min(-tc*trinv,1.0_r8)
1283       hlatsb = hlatv + weight*hlatf
1284       hlatvp = hlatv - 2369.0_r8*tc
1285       if (t < tmelt) then
1286          hltalt = hlatsb
1287       else
1288          hltalt = hlatvp
1289       end if
1290       if (lflg) then
1291          tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3) + tc*(pcf(4) + tc*pcf(5))))
1292       else
1293          tterm = 0.0_r8
1294       end if
1295       desdt  = hltalt*es/(rgasv*t*t) + tterm*trinv
1296       gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
1297       if (qs == 1.0_r8) gam = 0.0_r8
1298     ! Sungsu
1299       dqsdt = (cp/hltalt)*gam
1300     ! End by Sungsu
1301 !  end do
1302    return
1303 !
1304 ! No icephs or water to ice transition
1305 !
1306
1307 10 continue
1308
1309 !10 do i=1,len
1310 !
1311 ! Account for change of hlatv with t above freezing where
1312 ! constant slope is given by -2369 j/(kg c) = cpv - cw
1313 !
1314       hlatvp = hlatv - 2369.0_r8*(t-tmelt)
1315       if (icephs) then
1316          hlatsb = hlatv + hlatf
1317       else
1318          hlatsb = hlatv
1319       end if
1320       if (t < tmelt) then
1321          hltalt = hlatsb
1322       else
1323          hltalt = hlatvp
1324       end if
1325       desdt  = hltalt*es/(rgasv*t*t)
1326       gam = hltalt*qs*p*desdt/(cp*es*(p - omeps*es))
1327       if (qs == 1.0_r8) gam = 0.0_r8
1328     ! Sungsu
1329       dqsdt = (cp/hltalt)*gam
1330     ! End by Sungsu
1331
1332 !  end do
1333 !
1334    return
1335 !
1336 end subroutine vqsatd2_single
1337
1338
1339 end module wv_saturation