anchor.diff:      do k=1,nz
   
         ! Autoconversion and accretion rates following KW eq. 2.13a,b
  -      qrprod = qc(k) - (qc(k)-dt0*amax1(.001*(qc(k)-.001d0),0.))/(1.d0+dt0*2.2d0*qr(k)**.875)
  -      qc(k) = amax1(qc(k)-qrprod,0.)
  -      qr(k) = amax1(qr(k)+qrprod+sed(k),0.)
  +      qrprod = qc(k) - (qc(k)-dt0*max(.001*(qc(k)-.001d0),0.D0))/(1.d0+dt0*2.2d0*qr(k)**.875)
  +      qc(k) = max(qc(k)-qrprod,0.D0)
  +      qr(k) = max(qr(k)+qrprod+sed(k),0.D0)
   
         ! Saturation vapor mixing ratio (gm/gm) following KW eq. 2.11
         qvs = pc(k)*exp(f2x*(pk(k)*theta(k)-273.d0)   &
  
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)
  