! 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   returnend 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,iter10    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)   returnend subroutine findspend module cldwat