anchor.diff:      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 :: 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 
  +! 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]
  
positive.code:         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 (temperature < freezing_threshold) then
              a = 21.8745584
              b = 7.66
  
negative.code:         real,intent(inout),dimension(nx,ny)    :: rain,snow
          real,intent(in)    :: dt
          integer,intent(in) :: nx,ny,nz
  
          ! local variables
          real,allocatable,dimension(:) :: temperature
          integer :: i,j
  
  !       calculate these once for every call because they are only a function of dt
          cloud2snow = exp(-1.0*snow_formation_time_const*dt)
          cloud2rain = exp(-1.0*rain_formation_time_const*dt)
  