anchor.diff:  !   cloud2rain=exp(-1.0*rain_formation_time_const*dt)
   !   
   !   do i=1,10
  -!       call mp_simple(p,t,qv,qc,qr,qs,rain,snow,dt,dx2,dz,nz,.True.)
  -! !         write(*,*) "p=",p
  -!       write(*,*) t(2),qv(2),qc(2),qr(2),qs(2),rain,snow
  -! !         write(*,*) "t=",t
  +!       call mp_simple(pressure,temperature,qv,qc,qr,qs,rain,snow,dt,dx2,dz,nz,.True.)
  +! !         write(*,*) "pressure=",pressure
  +!       write(*,*) temperature(2),qv(2),qc(2),qr(2),qs(2),rain,snow
  +! !         write(*,*) "temperature=",temperature
   ! !         write(*,*) "qv=",qv
   ! !         write(*,*) "qc=",qc
   ! !         write(*,*) "qr=",qr
  
positive.code:     private
      public::mp_simple_driver
      
      real, parameter :: LH_vapor=2.26E6 ! J/kg
      real, parameter :: dLHvdt = 2400! APPROXIMATE increase in latent heat with a decrease in temperature (below 373.15K)
                                      ! derived from various curves plotted online
                                      ! not sure why this isn't 4186.0 (J/kg/K) (i.e. specific heat of water)
                                      ! NOTE, ice=2110 (J/kg/K), steam=2080 (J/kg/K)
      real, parameter :: LH_liquid=3.34E5 ! J/kg
      real, parameter :: heat_capacity = 1006.0 ! air heat capacity J/kg/K
      real, parameter :: SMALL_VALUE = 1E-15
      real, parameter :: SMALL_PRESSURE = 1000. ! in Pa this is actually pretty small... 
  !     real, parameter :: mp_R=287.058 ! J/(kg K) specific gas constant for air
  !     real, parameter :: mp_g=9.81 ! gravity m/s^2
  
  ! arbitrary calibratable timescales default values as used in the linear model
  ! these should be pulled out to a parameter file for calibration purposes 
  ! but this is approximately how they are implemented in SB04
      real,parameter :: snow_formation_time_const=1/2000.0 ! [1/s]
      real,parameter :: rain_formation_time_const=1/500.0  ! [1/s]
      real,parameter :: freezing_threshold=273.15          ! [K]
      real,parameter :: snow_fall_rate=1.5                 ! [m/s]   for a water vapor scale height of 3750m corresponds to tau_f = 2500
      real,parameter :: rain_fall_rate=10.0                ! [m/s]   for a water vapor scale height of 3750m corresponds to tau_f = 375
      real,parameter :: snow_cloud_init=0.00001            ! [kg/kg] cloud ice content before snow will start to form
      real,parameter :: rain_cloud_init=0.00001            ! [kg/kg] cloud water content before rain will start to form
      
      
  !   these are recalculated every call because they are a function of dt
  !   conversion "time" = exp( -1 * time_constant * dt)
      real :: rain_evap=0.999
      real :: snow_evap=0.999
      real :: snow_melt=0.999
      real :: cloud2rain=0.999
      real :: cloud2snow=0.999
      !$omp threadprivate(cloud2rain,cloud2snow,snow_melt,snow_evap,rain_evap)
      
      contains
      real function sat_mr(t,p)
      ! Calculate the saturated mixing ratio at a temperature (K), pressure (Pa)
          implicit none
          real,intent(in) :: t,p
          real :: e_s,mr_s,a,b
  
  !       from http://www.dtic.mil/dtic/tr/fulltext/u2/778316.pdf
  !           Lowe, P.R. and J.M. Ficke., 1974: THE COMPUTATION OF SATURATION VAPOR PRESSURE 
  !               Environmental Prediction Research Facility, Technical Paper No. 4-74
  !       which references:
  !           Murray, F. W., 1967: On the computation of saturation vapor pressure. 
  !               Journal of Applied Meteorology, Vol. 6, pp. 203-204.
  !       Also notes a 6th order polynomial and look up table as viable options. 
          if (t<freezing_threshold) then
              a=21.8745584
              b=7.66
          else
              a=17.2693882
              b=35.86
          endif
          e_s = 610.78* exp(a*(t-273.16)/(t-b)) !(Pa)
  
  !       alternate formulations
  !       Polynomial:
  !       e_s = ao + t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+a6*t))))) a0-6 defined separately for water and ice
  !       e_s = 611.2*exp(17.67*(t-273.15)/(t-29.65)) ! (Pa)
          !from : http://www.srh.noaa.gov/images/epz/wxcalc/vaporPressure.pdf
  !       e_s = 611.0*10.0**(7.5*(t-273.15)/(t-35.45))
          
          
  !       if ((p-e_s)<0) then
  !           print*, "e_s more than p!"
  !           print*, p,e_s,t,p-e_s
  !       endif
  !       e_s=min(e_s,p-SMALL_PRESSURE)
          if ((p-e_s)<=0) then
              e_s=p*0.99999
          endif
          !from : http://www.srh.noaa.gov/images/epz/wxcalc/mixingRatio.pdf
          sat_mr=0.6219907*e_s/(p-e_s) !(kg/kg)
      end function sat_mr
      
  
      subroutine cloud_conversion(p,t,qv,qc,qvsat,dt)
          implicit none
          real,intent(inout)::t,qv,qc,qvsat
          real,intent(in)::dt,p
          real :: vapor2temp,excess,deltat,pre_qc,pre_qv,pre_t,lastqv
          integer :: iteration
          real :: maxerr
          
          maxerr=1e-6
          iteration=0
          lastqv=qv+maxerr*2
          vapor2temp=(LH_vapor+(373.15-t)*dLHvdt)/heat_capacity
          pre_qc=qc !DEBUG
          pre_qv=qv !DEBUG
          pre_t=t !DEBUG
          excess=0
          
          do while ((abs(lastqv-qv)>maxerr).and.(iteration<5))
              iteration=iteration+1
              lastqv=qv
              ! calculate the saturating mixing ratio
              qvsat=sat_mr(t,p)
              ! if saturated create clouds
              if (qv>qvsat) then
                  excess=(qv-qvsat)*0.95
                  ! temperature change if all vapor is converted
                  t=t+(excess*vapor2temp)
                  qv=qv-excess
                  qc=qc+excess
              ! if unsaturated and clouds exist, evaporate clouds
              else if (qc>0) then
                  excess=(qvsat-qv)*0.95
                  if (excess<qc) then
                      t=t-(excess*vapor2temp)
                      qv=qv+excess
                      qc=qc-excess
                  else
                      qv=qv+qc
                      t=t-(qc*vapor2temp)
                      excess=qc
                      qc=0.
                  endif
                  excess=excess*(-1) !DEBUG
              endif
          enddo
          if (iteration==5) then
  !           print*, iteration
  !           print*, pre_qc,pre_qv,pre_t,p
  !           print*, qc, qv,t, qvsat
  !           print*, sat_mr(pre_t,p),sat_mr(t,p), sat_mr(pre_t,p)-sat_mr(t,p)
              qv=sat_mr(pre_t,p)
              t=pre_t
              qc=pre_qc
          endif
  !       if (t<170) then
  !           t=170
  !       else if (t>350) then
  !           t=350
  !       endif
          
          qc=max(qc,0.)
          if ((t>350).or.(qc>0.01).or.(qvsat>1)) then
              deltat=excess*vapor2temp
              print*, "mp_simple: data out of bounds"
              print*, "iter=",iteration
              print*, "preqc=",pre_qc,"preqv=",pre_qv,"pret=",pre_t,"p=",p
              print*, "qc=",qc, "qv=",qv,"t=",t, "qvsat=",qvsat
              print*, "mrs_pret=",sat_mr(pre_t,p),"mrs_t=",sat_mr(t,p), "mrs_delta=",sat_mr(pre_t,p)-sat_mr(t,p)
  
              print*, "qv=",qv, "qvs=",qvsat, "qc=",qc, "preqc=",pre_qc, "t=",t, "excess=",excess, "vapor2temp=",vapor2temp
              print*, "t-deltat",t-deltat, "p=",p, "mrs_t-dT=",sat_mr(t-deltat,p)
          endif
          
      end subroutine
  
  
      subroutine cloud2hydrometeor(qc,q,conversion,qcmin)
          implicit none
          real,intent(inout) :: qc,q
  
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
  