anchor.diff:                        + qr(k) * CVW(I_QR)
                endif
   
  -             fq(ij,kk,I_QV) = ( qv(k) - q(ij,kk,I_QV) ) / dt
  -             fe(ij,kk)      = ( cv(k) * t(ij,k) - ein(ij,kk) ) / dt
  +             fq(ij,kk,I_QV) = fq(ij,kk,I_QV) + ( qv(k) - q(ij,kk,I_QV) ) / dt
  +             fe(ij,kk)      = fe(ij,kk)      + ( cv(k) * t(ij,k) - ein(ij,kk) ) / dt
             enddo
          enddo
   
  
positive.code:        do k = 1, vlayer
            kk = kmax - k + 1 ! reverse order
  
            fvx(:,kk) = ( u(:,k) * ix(:) + v(:,k) * jx(:) - vx(:,kk) ) / dt
            fvy(:,kk) = ( u(:,k) * iy(:) + v(:,k) * jy(:) - vy(:,kk) ) / dt
            fvz(:,kk) = ( u(:,k) * iz(:) + v(:,k) * jz(:) - vz(:,kk) ) / dt
         enddo
  
         do k = 1, vlayer
  
negative.code:                 + qc(:) * CVW(I_QC) &
                  + qr(:) * CVW(I_QR)
  
            fq(ij,kmin:kmax,I_QV) = fq(ij,kmin:kmax,I_QV) + ( qv(:) - q(ij,kmin:kmax,I_QV) ) / dt
            fq(ij,kmin:kmax,I_QC) = fq(ij,kmin:kmax,I_QC) + ( qc(:) - q(ij,kmin:kmax,I_QC) ) / dt
            fq(ij,kmin:kmax,I_QR) = fq(ij,kmin:kmax,I_QR) + ( qr(:) - q(ij,kmin:kmax,I_QR) ) / dt
            fe(ij,kmin:kmax)      = fe(ij,kmin:kmax)      + ( cv(:) * theta(:) * pk(:) - ein(ij,kmin:kmax) ) / dt
         enddo
      endif
  
  