anchor.diff:              if (qv<qvsat) then
                   if (qr>SMALL_VALUE) then
                       ! evaporate rain
  -                    call phase_change(p,t,qr,qvsat,qv,L_evap,1.0) !cloud2rain)
  +                    call phase_change(pressure,temperature,qr,qvsat,qv,L_evap,1.0) !cloud2rain)
                   endif
                   if (qs>SMALL_VALUE) then
                       ! sublimate snow
  -                    call phase_change(p,t,qs,qvsat,qv,L_subl,1.0) !cloud2snow)
  +                    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
  
positive.code:         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
              ! convert specific humidity to mixing ratio
  !           qv(i)=qv(i)/(1-qv(i))
              call mp_conversions(p(i),t(i),qv(i),qc(i),qr(i),qs(i),dt)
          enddo
  
          ! SEDIMENTATION for rain
          if (maxval(qr)>SMALL_VALUE) then
              fall_rate=rain_fall_rate
              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
  
negative.code: !!   mp_simple           - calls mp_conversions for all z, then calls sediment
  !!   mp_conversions      - handles all microphysics conversions (e.g. vapor->cloud->rain[->snow->vapor],...)
  !!   cloud_conversion    - uses sat_mr to calculate sub or supersaturation and changes vapor and cloud water to balance 
  !!                          (adjusts temperature for Latent Heating)
  !!   cloud2hydrometeor   - converts cloud water to rain or snow (temperature dependant)
  !!   phase_change        - handles rain and snow evaporation (can take a time constant)
  !!   sediment            - advects falling rain and snow due to gravity (not wind)
  !!   sat_mr              - calculate saturated water vapor mixing ratio
  !!
  !! Driver inputs: pressure,th,pii,rho,qv,qc,qr,qs,rain,snow,dt,dz,nx,ny,nz
  !!   pressure   = pressure                      - 3D - input  - Pa     - (nx,nz,ny)
  !!   th         = potential temperature         - 3D - in/out - K      - (nx,nz,ny)
  !!   pii        = exner function                - 3D - input  - []     - (nx,nz,ny)
  !!   rho        = air density                   - 3D - input  - kg/m^3 - (nx,nz,ny)
  !!   qv         = specific humidity             - 3D - in/out - kg/kg  - (nx,nz,ny)
  !!   qc         = cloud water content           - 3D - in/out - kg/kg  - (nx,nz,ny)
  !!   qr         = rain water content            - 3D - in/out - kg/kg  - (nx,nz,ny)
  !!   qs         = snow water content            - 3D - in/out - kg/kg  - (nx,nz,ny)
  !!   rain       = accumulated rain              - 2D - output - mm     - (nx,ny)
  !!   snow       = accumulated snow              - 2D - output - mm     - (nx,ny)
  !!   dt         = time step                     - 0D - input  - sec.   - scalar
  !!   nx         = number of ew grid cells       - 0D - input  - n      - scalar
  !!   ny         = number of ns grid cells       - 0D - input  - n      - scalar
  !!   nz         = number of vertical grid cells - 0D - input  - n      - scalar
  !! </pre>
  !!
  !!  @author
  