cldwat.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

! the following lines calculate the slope parameter and snow mixing ratio
! from the precip rate using the equations found in lin et al 83
! in the most natural form, but it is expensive, so after some tedious
! algebraic manipulation you can use the cheaper form found below
!            vfalls = c*gam4pd/(6*lamdas**d)*sqrt(rhonot/rhocgs)
!     $               *0.01   ! convert from cm/s to m/s
!            snowmr(i) = snowfr*precab(i)/(rho(i)*vfalls*cldpr(i))
!            snowmr(i) = ( prscgs(i)*mcon02 * (rhocgs**mcon03) )**mcon04
!            lamdas = (prhonos/max(rhocgs*snowmr(i),small))**0.25
!            csacw = mcon01*sqrt(rhonot/rhocgs)/(lamdas**thrpd)
!
! coefficient for collection by snow independent of phase
!
      csacx = mcon07*rhocgs**mcon08*prscgs(i)**mcon05
!
! collection of liquid by snow (lin et al 1983)
!
      psacw = csacx*liqmr(i)*esw
#ifdef PERGRO
! this is necessary for pergro
      psacw = 0.
#endif
!
! collection of ice by snow (lin et al 1983)
!
      psaci = csacx*icemr(i)*esi
!
! total conversion of condensate to precipitate
!
      ptot = pwaut + psaut + pracw + psacw + psaci
!
! the recipricol of cloud water amnt (or zero if no cloud water)
!
!         rcwm =  totmr(i)/(max(totmr(i),small)**2)
!
! turn the tendency back into a loss rate (1/seconds)
!
      if (totmr(i) > 0.) then
         coef(i) = ptot/totmr(i)
      else
         coef(i) = 0.
      endif

      fwaut(i) = pwaut/max(ptot,small)
      fsaut(i) = psaut/max(ptot,small)
      fracw(i) = pracw/max(ptot,small)
      fsacw(i) = psacw/max(ptot,small)
      fsaci(i) = psaci/max(ptot,small)

   end do
#ifdef DEBUG
   i = icollook(nlook)
   if (lchnk == lchnklook(nlook) ) then
      write (6,*)
      write (6,*) '------', k
      write (6,*) ' liqmr, rainmr,precab ', liqmr(i), rainmr(i), precab(i)*8.64e4
      write (6,*) ' frac: waut,saut,racw,sacw,saci ', &
           fwaut(i), fsaut(i), fracw(i), fsacw(i), fsaci(i)
   endif
#endif

   return
end subroutine findmcnew

!##############################################################################

subroutine findsp (lchnk, ncol, q, t, p, tsp, qsp)
!----------------------------------------------------------------------- 
! 
! 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
!
   integer, intent(in) :: lchnk                 ! chunk identifier
   integer, intent(in) :: ncol                  ! number of atmospheric columns

   real(r8), intent(in) :: q(pcols,pver)        ! water vapor (kg/kg)
   real(r8), intent(in) :: t(pcols,pver)        ! temperature (K)
   real(r8), intent(in) :: p(pcols,pver)        ! pressure    (Pa)
!
! output arguments
!
   real(r8), intent(out) :: tsp(pcols,pver)      ! saturation temp (K)
   real(r8), intent(out) :: qsp(pcols,pver)      ! saturation mixing ratio (kg/kg)
!
! local variables
!
   integer i                 ! work variable
   integer k                 ! work variable
   logical lflg              ! work variable
   integer iter              ! work variable
   integer l                 ! work variable

   real(r8) omeps                ! 1 minus epsilon
   real(r8) trinv                ! work variable
   real(r8) es                   ! sat. vapor pressure
   real(r8) desdt                ! change in sat vap pressure wrt temperature
!     real(r8) desdp                ! change in sat vap pressure wrt pressure
   real(r8) dqsdt                ! change in sat spec. hum. wrt temperature
   real(r8) dgdt                 ! work variable
   real(r8) g                    ! work variable
   real(r8) weight(pcols)        ! work variable
   real(r8) hlatsb               ! (sublimation)
   real(r8) hlatvp               ! (vaporization)
   real(r8) hltalt(pcols,pver)   ! lat. heat. of vap.
   real(r8) tterm                ! work var.
   real(r8) qs                   ! spec. hum. of water vapor
   real(r8) tc                   ! crit temp of transition to ice

! work variables
   real(r8) t1, q1, dt, dq
   real(r8) dtm, dqm
   real(r8) qvd, a1, tmp
   real(r8) rair
   real(r8) r1b, c1, c2, c3
   real(r8) denom
   real(r8) dttol
   real(r8) dqtol
   integer doit(pcols) 
   real(r8) enin(pcols), enout(pcols)
   real(r8) tlim(pcols)

   omeps = 1.0 - epsqs
   trinv = 1.0/ttrice
   a1 = 7.5*log(10._r8)
   rair =  287.04
   c3 = rair*a1/cp
   dtm = 0.    ! needed for iter=0 blowup with f90 -ei
   dqm = 0.    ! needed for iter=0 blowup with f90 -ei
   dttol = 1.e-4 ! the relative temp error tolerance required to quit the iteration
   dqtol = 1.e-4 ! the relative moisture error tolerance required to quit the iteration
!  tmin = 173.16 ! the coldest temperature we can deal with
!
! max number of times to iterate the calculation
   iter = 8
!
   do k = k1mb,pver

!
! first guess on the wet bulb temperature
!
      do i = 1,ncol

#ifdef DEBUG
         if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
            write (6,*) ' '
            write (6,*) ' level, t, q, p', k, t(i,k), q(i,k), p(i,k)
         endif
#endif
! limit the temperature range to that relevant to the sat vap pres tables
         tlim(i) = min(max(t(i,k),173._r8),373._r8)
         es = estblf(tlim(i))
         denom = p(i,k) - omeps*es
         qs = epsqs*es/denom
         doit(i) = 0
         enout(i) = 1.
! make sure a meaningful calculation is possible
         if (p(i,k) > 5.*es .and. qs > 0. .and. qs < 0.5) then
!
! Saturation specific humidity
!
             qs = min(epsqs*es/denom,1._r8)
!
! "generalized" analytic expression for t derivative of es
!  accurate to within 1 percent for 173.16 < t < 373.16
!
! 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     = tlim(i) - t0
             lflg   = (tc >= -ttrice .and. tc < 0.0)
             weight(i) = min(-tc*trinv,1.0_r8)
             hlatsb = hlatv + weight(i)*hlatf
             hlatvp = hlatv - 2369.0*tc
             if (tlim(i) < t0) then
                hltalt(i,k) = hlatsb
             else
                hltalt(i,k) = hlatvp
             end if
             enin(i) = cp*tlim(i) + hltalt(i,k)*q(i,k)

! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
             tmp =  q(i,k) - qs
             c1 = hltalt(i,k)*c3
             c2 = (tlim(i) + 36.)**2
             r1b    = c2/(c2 + c1*qs)
             qvd   = r1b*tmp
             tsp(i,k) = tlim(i) + ((hltalt(i,k)/cp)*qvd)
#ifdef DEBUG
             if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
                write (6,*) ' relative humidity ', q(i,k)/qs
                write (6,*) ' first guess ', tsp(i,k)
             endif
#endif
             es = estblf(tsp(i,k))
             qsp(i,k) = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
          else
             doit(i) = 1
             tsp(i,k) = tlim(i)
             qsp(i,k) = q(i,k)
             enin(i) = 1.
          endif
       end do   ! end do i
!
! now iterate on first guess
!
      do l = 1, iter
         dtm = 0
         dqm = 0
         do i = 1,ncol
            if (doit(i) == 0) then
               es = estblf(tsp(i,k))
!
! Saturation specific humidity
!
               qs = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
!
! "generalized" analytic expression for t derivative of es
! accurate to within 1 percent for 173.16 < t < 373.16
!
! 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     = tsp(i,k) - t0
               lflg   = (tc >= -ttrice .and. tc < 0.0)
               weight(i) = min(-tc*trinv,1.0_r8)
               hlatsb = hlatv + weight(i)*hlatf
               hlatvp = hlatv - 2369.0*tc
               if (tsp(i,k) < t0) then
                  hltalt(i,k) = hlatsb
               else
                  hltalt(i,k) = 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(i,k)*es/(rgasv*tsp(i,k)*tsp(i,k)) + tterm*trinv
               dqsdt = (epsqs + omeps*qs)/(p(i,k) - omeps*es)*desdt
!              g = cp*(tlim(i)-tsp(i,k)) + hltalt(i,k)*q(i,k)- hltalt(i,k)*qsp(i,k)
               g = enin(i) - (cp*tsp(i,k) + hltalt(i,k)*qsp(i,k))
               dgdt = -(cp + hltalt(i,k)*dqsdt)
               t1 = tsp(i,k) - g/dgdt
               dt = abs(t1 - tsp(i,k))/t1
               tsp(i,k) = max(t1,tmin)
               es = estblf(tsp(i,k))
               q1 = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
               dq = abs(q1 - qsp(i,k))/max(q1,1.e-12_r8)
               qsp(i,k) = q1
#ifdef DEBUG
               if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
                  write (6,*) ' rel chg lev, iter, t, q ', k, l, dt, dq, g
               endif
#endif
               dtm = max(dtm,dt)
               dqm = max(dqm,dq)
! if converged at this point, exclude it from more iterations
               if (dt < dttol .and. dq < dqtol) then
                  doit(i) = 2
               endif
               enout(i) = cp*tsp(i,k) + hltalt(i,k)*qsp(i,k)
! bail out if we are too near the end of temp range
               if (tsp(i,k) < 174.16) then
                  doit(i) = 4
               endif
            else
            endif
         end do              ! do i = 1,ncol

         if (dtm < dttol .and. dqm < dqtol) then
            go to 10
         endif

      end do                 ! do l = 1,iter
10    continue

      if (dtm > dttol .or. dqm > dqtol) then
         do i = 1,ncol
            if (doit(i) == 0) then
               write (6,*) ' findsp not converging at point i, k ', i, k
               write (6,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
               write (6,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
               call endrun ()
            endif
         end do
      endif
      do i = 1,ncol
         if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4) then
            write (6,*) ' the enthalpy is not conserved for point ', &
                 i, k, enin(i), enout(i)
            write (6,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
            write (6,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
            call endrun ()
         endif
      end do
      
   end do                    ! level loop (k=1,pver)

   return
end subroutine findsp

end module cldwat