anchor.diff:  !!   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 T for Latent Heating)
  -!!   cloud2hydrometeor   - converts cloud water to rain or snow (T dependant)
  +!!                          (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: p,th,pii,rho,qv,qc,qr,qs,rain,snow,dt,dz,nx,ny,nz
  -!!   p   = 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  - seconds    - 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
  +!! 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
  
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:         real::delta
          
          if (qc > qcmin) then
              delta = qc-(qc*conversion)
          else
              delta = 0
          endif
          
          if (delta<qc) then
              qc = qc-delta
              q = q+delta
          else
              q = q+qc
              qc = 0.
          endif
          qc = max(qc,0.)
      end subroutine
      
      !>----------------------------------------------------------
      !!  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!*(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
          ! bounds checking
          if (q1<0) then 
              if ((q1+SMALL_VALUE)<0) then
                  q1 = 0
              else
                  print*, "phase_change"
                  print*, q1,q2,delta,qmax,change_rate
                  stop
              endif
          endif
          q2 = q2+delta
          temperature = temperature+delta*mass2temp
          
      end subroutine
      
      !>----------------------------------------------------------
      !!  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-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 (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(pressure,temperature,qs,100.,qr,L_melt,dt/snow_melt_time)
  ! !                         write(*,*) "Snow Melt",temperature
  !                   endif
                  else
                      ! convert cloud water to snow flakes
  