anchor.diff:        prod = (qv(k)-qvs)/(1.d0+qvs*f5/(pk(k)*theta(k)-36.d0)**2)
   
         ! Evaporation rate following KW eq. 2.14a,b
  -      ern = amin1(dt0*(((1.6d0+124.9d0*(r(k)*qr(k))**.2046)  &
  +      ern = min(dt0*(((1.6d0+124.9d0*(r(k)*qr(k))**.2046)  &
               *(r(k)*qr(k))**.525)/(2550000d0*pc(k)            &
               /(3.8d0 *qvs)+540000d0))*(dim(qvs,qv(k))         &
  -            /(r(k)*qvs)),amax1(-prod-qc(k),0.),qr(k))
  +            /(r(k)*qvs)),max(-prod-qc(k),0.D0),qr(k))
   
         ! Saturation adjustment following KW eq. 3.10
  -      theta(k)= theta(k) + 2500000d0/(1003.d0*pk(k))*(amax1( prod,-qc(k))-ern)
  -      qv(k) = amax1(qv(k)-max(prod,-qc(k))+ern,0.)
  +      theta(k)= theta(k) + 2500000d0/(1003.d0*pk(k))*(max( prod,-qc(k))-ern)
  +      qv(k) = max(qv(k)-max(prod,-qc(k))+ern,0.D0)
         qc(k) = qc(k)+max(prod,-qc(k))
         qr(k) = qr(k)-ern
       end do
  
positive.code: !
  !    Klemp, J. B., W. C. Skamarock, W. C., and S.-H. Park, 2015:
  !    Idealized Global Nonhydrostatic Atmospheric Test Cases on a Reduced
  !    Radius Sphere. Journal of Advances in Modeling Earth Systems. 
  !    doi:10.1002/2015MS000435
  !
  !=======================================================================
  
negative.code:   !   Input / output parameters
    !------------------------------------------------
  
    INTEGER, INTENT(IN) :: nz ! Number of thermodynamic levels in the column
    REAL(8), DIMENSION(nz), INTENT(INOUT) :: &
              theta   ,     & ! Potential temperature (K)
              qv      ,     & ! Water vapor mixing ratio (gm/gm)
  