anchor.diff:          real::delta
           
           if (qc > qcmin) then
  -            delta=qc-(qc*conversion)
  +            delta = qc-(qc*conversion)
           else
  -            delta=0
  +            delta = 0
           endif
           
           if (delta<qc) then
  -            qc=qc-delta
  -            q=q+delta
  +            qc = qc-delta
  +            q = q+delta
           else
  -            q=q+qc
  -            qc=0.
  -        endif
  -        qc=max(qc,0.)
  +            q = q+qc
  +            qc = 0.
  +        endif
  +        qc = max(qc,0.)
       end subroutine
       
  -    subroutine phase_change(p,t,q1,qmax,q2,Lheat,change_rate)
  -        ! convert from q1 (e.g. rain) to q2 (e.g. vapor)
  -        ! qmax = the maximum value for q2 (e.g. qv_sat)
  -        ! p = pressure (not used)
  -        ! t = temperature
  -        ! Lheat = latent heat associated with the phase change
  -        ! change_rate = fraction of q1 that can change in a timestep
  -        implicit none
  -        real, intent(inout)::t,q1,q2
  -        real,intent(in) :: p,qmax,change_rate,Lheat
  +    !>----------------------------------------------------------
  +    !!  Change the "phase" of a hydrometeor (e.g. evaporate rain)
  +    !! 
  +    !!  Written to apply generically to any conversion. 
  +    !!  Convert from q1 (e.g. rain) to q2 (e.g. vapor) using Lheat to affect temperature
  +    !!  and change_rate to control the conversion time scale
  +    !!
  +    !! @param pressure    air pressure (not used)                       [Pa]
  +    !! @param temperature air temperature                               [K]
  +    !! @param q1          hydrometeor to convert from                   [kg/kg]
  +    !! @param qmax        the maximum value for q2 (e.g. qv_sat)        [kg/kg]
  +    !! @param q2          hydrometeor to convert to                     [kg/kg]
  +    !! @param Lheat       latent heat associated with the phase change  [J/kg]
  +    !! @param change_rate fraction of q1 that can change in a timestep  []
  +    !!
  +    !!----------------------------------------------------------
  +    subroutine phase_change(pressure,temperature,q1,qmax,q2,Lheat,change_rate)
  +        implicit none
  +        real, intent(inout)::temperature,q1,q2
  +        real,intent(in) :: pressure,qmax,change_rate,Lheat
           real :: mass2temp,delta
  -        mass2temp=Lheat/heat_capacity!*(p/(R*t)*dV))
  -        
  -        delta=(qmax-q2)*change_rate
  -        if (delta>q1) delta=q1
  +        mass2temp = Lheat/heat_capacity!*(pressure/(R*temperature)*dV))
  +        
  +        delta = (qmax-q2)*change_rate
  +        if (delta>q1) delta = q1
           ! hopefully we don't over shoot saturation (use a 1% buffer)
           if (delta>((qmax-q2)*0.99)) then
  -            delta=(qmax-q2)*0.99
  -        endif
  -        
  -        q1=q1-delta
  +            delta = (qmax-q2)*0.99
  +        endif
  +        
  +        q1 = q1-delta
  +        ! bounds checking
           if (q1<0) then 
               if ((q1+SMALL_VALUE)<0) then
  -                q1=0
  +                q1 = 0
               else
                   print*, "phase_change"
                   print*, q1,q2,delta,qmax,change_rate
                   stop
               endif
           endif
  -        q2=q2+delta
  -        t=t+delta*mass2temp
  +        q2 = q2+delta
  +        temperature = temperature+delta*mass2temp
           
       end subroutine
       
  -    subroutine mp_conversions(p,t,qv,qc,qr,qs,dt)
  -        implicit none
  -        real, intent(inout) :: p,t,qv,qc,qr,qs
  +    !>----------------------------------------------------------
  +    !!  Compute microphysical conversions
  +    !! 
  +    !!  Convert cloud water to and from vapor 
  +    !!  Convert cloud water to rain, and cloud ice to snow
  +    !!  Cloud water is assumed to be ice if temperature is less than 0C
  +    !!  Conversions also include latent heat feedbacks to temperature
  +    !! 
  +    !!  @param pressure     air pressure             [Pa]
  +    !!  @param temperature  air temperature          [K]
  +    !!  @param qv           water vapor mixing ratio [kg/kg]
  +    !!  @param qc           cloud water mixing ratio [kg/kg]
  +    !!  @param qr           rain mixing ratio        [kg/kg]
  +    !!  @param qs           snow mixing ratio        [kg/kg]
  +    !!  @param dt           time step                [sec.]
  +    !!
  +    !!----------------------------------------------------------
  +    subroutine mp_conversions(pressure,temperature,qv,qc,qr,qs,dt)
  +        implicit none
  +        real, intent(inout) :: pressure,temperature,qv,qc,qr,qs
           real,intent(in)::dt
           real :: qvsat,L_evap,L_subl,L_melt
           
  -        qvsat=1
  -        L_melt=-1*LH_liquid  ! J/kg (should change with temperature)
  -        L_evap=-1*(LH_vapor+(373.15-t)*dLHvdt)   ! J/kg
  -        L_subl=L_melt+L_evap ! J/kg
  -        !convert cloud water to and from water vapor
  -        call cloud_conversion(p,t,qv,qc,qvsat,dt)
  +        qvsat = 1
  +        L_melt = -1*LH_liquid  ! J/kg (should change with temperature)
  +        L_evap = -1*(LH_vapor+(373.15-temperature)*dLHvdt)   ! J/kg
  +        L_subl = L_melt+L_evap ! J/kg
  +        ! convert cloud water to and from water vapor
  +        call cloud_conversion(pressure,temperature,qv,qc,qvsat,dt)
           ! if there are no species to process we will just return
           if ((qc+qr+qs) >SMALL_VALUE) then
               if (qc>SMALL_VALUE) then
  -                if (t>freezing_threshold) then
  +                if (temperature>freezing_threshold) then
                       ! convert cloud water to rain drops
                       call cloud2hydrometeor(qc,qr,cloud2rain,rain_cloud_init)
   !                   if (qs>SMALL_VALUE) then
   !                       ! it is above freezing, so start melting any snow if present
  -!                       call phase_change(p,t,qs,100.,qr,L_melt,dt/snow_melt_time)
  -! !                         write(*,*) "Snow Melt",t
  +!                       call phase_change(pressure,temperature,qs,100.,qr,L_melt,dt/snow_melt_time)
  +! !                         write(*,*) "Snow Melt",temperature
   !                   endif
                   else
                       ! convert cloud water to snow flakes
  
positive.code:         integer :: i
          
  !       calculate the mass of material falling out of the bottom model level
          sediment=v(1)*q(1)*rho(1) ![m] * [kg/kg] * [kg/m^3] = [kg/m^2]
  !       remove that from the bottom model layer. 
          q(1)=q(1)-(sediment/dz(1)/rho(1)) ! [kg/m^2] / [m] / [kg/m^3] = [kg/kg]
          do i=1,n-1
              flux(i)=v(i+1)*q(i+1)*rho(i+1)
          enddo
          do i=1,n-1
              q(i)=q(i)+flux(i)/(rho(i)*dz(i))
              q(i+1)=q(i+1)-flux(i)/(rho(i+1)*dz(i+1))
          enddo
      
      end function
  
      subroutine mp_simple(p,t,rho,qv,qc,qr,qs,rain,snow,dt,dz,nz,debug)
          implicit none
          real,intent(inout),dimension(nz)::p,t,rho,qv,qc,qr,qs
          real,intent(inout)::rain,snow
          real,intent(in),dimension(nz)::dz
          real,intent(in)::dt
  
negative.code:             if (qv<qvsat) then
                  if (qr>SMALL_VALUE) then
                      ! evaporate rain
                      call phase_change(pressure,temperature,qr,qvsat,qv,L_evap,1.0) !cloud2rain)
                  endif
                  if (qs>SMALL_VALUE) then
                      ! sublimate snow
                      call phase_change(pressure,temperature,qs,qvsat,qv,L_subl,1.0) !cloud2snow)
                  endif
              endif
          endif
      end subroutine
      
      !>----------------------------------------------------------
      !!  Compute sedimentation in a column
      !! 
      !!  Takes a mixing ratio of some species along with the fall velocities,
      !!  air densities and layer thicknesses.  
      !! 
      !!  @param  q   1D microphysical species to sediment (snow, rain, etc) [kg/kg]
      !!  @param  v   1D vertical velocity [m/s]
      !!  @param  rho 1D air density [kg/m^3]
      !!  @param  dz  1D layer thickness [m]
      !!  @param  n   0D number of layers 
      !!
      !!----------------------------------------------------------
      real function sediment(q,v,rho,dz,n)
          implicit none
          real,intent(inout),dimension(n)::q
  