anchor.diff:          logical,intent(in)::debug
           
           real,dimension(nz)::fall_rate
  -        real::cfl,snowfall
  -        integer::i
  +        real::cfl,snowfall, qvsat
  +        integer::i,cfl_step
  +        real :: L_evap,L_subl,L_melt
  +        
  +        qvsat=1
  +        L_melt=-1*LH_liquid  ! J/kg (should change with temperature)
           
           if ((debug).and.(dt<1)) print*, "internal dt=",dt
           do i=1,nz
  
positive.code:         real :: mass2temp,delta
          mass2temp=Lheat/heat_capacity!*(p/(R*t)*dV))
          
          delta=q1-(q1*change_rate)
          ! hopefully we don't over shoot saturation (use a 10% buffer)
          if (delta>((qmax-q2)*0.90)) then
              delta=(qmax-q2)*0.90
          endif
          
          q1=q1-delta
  		if (q1<0) then 
  			print*, "phase_change"
  			print*, q1,q2,delta,qmax,change_rate
  			stop
  		endif
          q2=q2+delta
          t=t+delta*mass2temp
      
      end subroutine
      
      subroutine mp_conversions(p,t,qv,qc,qr,qs,dt)
  
negative.code:             cfl=ceiling(maxval(dt/dz*fall_rate))
              fall_rate=dt*fall_rate/cfl
              ! substepping to satisfy CFL criteria
              do cfl_step=1,nint(cfl)
                  rain=rain+sediment(qr,fall_rate,rho,dz,nz)
  				! allow any rain that reached a unsaturated layer to evaporate
  				do i=1,nz
  			        L_evap=-1*(LH_vapor+(373.15-t(i))*dLHvdt)   ! J/kg
  					qvsat=sat_mr(t(i),p(i))
  		            if (qv(i)<qvsat) then
  		                if (qr(i)>SMALL_VALUE) then
  		                    ! evaporate rain
  		                    call phase_change(p(i),t(i),qr(i),qvsat,qv(i),L_evap,1.0) !cloud2rain)
  		                endif
  					endif
  				enddo
              enddo
          endif
          
  