wv_saturation.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

! Author: J. Hack
! 
!------------------------------Arguments--------------------------------
!
! Input arguments
!
   integer, intent(in) :: ii            ! I dimension of arrays t, p, es, qs
   integer, intent(in) :: ilen          ! Vector length in I direction
   integer, intent(in) :: kk            ! K dimension of arrays t, p, es, qs
   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
   real(r8), intent(out) :: gam(ii,kk)       ! (l/cp)*(d(qs)/dt)
!
!---------------------------Local workspace-----------------------------
!
   logical lflg          ! True if in temperature transition region
   integer i             ! i index for vector calculations
   integer k             ! k index
   real(r8) omeps            ! 1. - 0.622
   real(r8) trinv            ! Reciprocal of ttrice (transition range)
   real(r8) tc               ! Temperature (in degrees C)
   real(r8) weight           ! Weight for es transition from water to ice
   real(r8) hltalt           ! Appropriately modified hlat for T derivatives
   real(r8) hlatsb           ! hlat weighted in transition region
   real(r8) hlatvp           ! hlat modified for t changes above freezing
   real(r8) tterm            ! Account for d(es)/dT in transition region
   real(r8) desdt            ! d(es)/dT
!
!-----------------------------------------------------------------------
!
   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 qs
! values which 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
!
! "generalized" analytic expression for t derivative of es
! accurate to within 1 percent for 173.16 < t < 373.16
!
   trinv = 0.0
   if ((.not. icephs) .or. (ttrice.eq.0.0)) go to 10
   trinv = 1.0/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)
         weight = min(-tc*trinv,1.0_r8)
         hlatsb = hlatv + weight*hlatf
         hlatvp = hlatv - 2369.0*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
         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) gam(i,k) = 0.0
      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*(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) gam(i,k) = 0.0
      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       ! vector length
   real(r8), intent(in) :: t(len)       ! temperature
   real(r8), intent(in) :: p(len)       ! pressure
!
! Output arguments
!
   real(r8), intent(out) :: es(len)   ! saturation vapor pressure
   real(r8), intent(out) :: qs(len)   ! saturation specific humidity
   real(r8), intent(out) :: gam(len)  ! (l/cp)*(d(qs)/dt)
!
!--------------------------Local Variables------------------------------
!
   logical lflg   ! true if in temperature transition region
!
   integer i      ! index for vector calculations
!
   real(r8) omeps     ! 1. - 0.622
   real(r8) trinv     ! reciprocal of ttrice (transition range)
   real(r8) tc        ! temperature (in degrees C)
   real(r8) weight    ! weight for es transition from water to ice
   real(r8) hltalt    ! appropriately modified hlat for T derivatives
!
   real(r8) hlatsb    ! hlat weighted in transition region
   real(r8) hlatvp    ! hlat modified for t changes above freezing
   real(r8) tterm     ! account for d(es)/dT in transition region
   real(r8) desdt     ! d(es)/dT
!
!-----------------------------------------------------------------------
!
   omeps = 1.0 - epsqs
   do i=1,len
      es(i) = estblf(t(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) then
         qs(i) = 1.0
         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
   if ((.not. icephs) .or. (ttrice.eq.0.0)) go to 10
   trinv = 1.0/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)
      weight = min(-tc*trinv,1.0_r8)
      hlatsb = hlatv + weight*hlatf
      hlatvp = hlatv - 2369.0*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
      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) gam(i) = 0.0
   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*(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) gam(i) = 0.0
   end do
!
   return
!
end subroutine vqsatd








end module wv_saturation