anchor.diff:          v = cos(direction) * magnitude
       end function calc_v
   
  -    ! return the moist / saturated adiabatic lapse rate for a given
  -    ! Temperature and mixing ratio (really MR could be calculated as f(T))
  -    ! from http://glossary.ametsoc.org/wiki/Saturation-adiabatic_lapse_rate
  +    !>----------------------------------------------------------
  +    !! Calculate the saturated adiabatic lapse rate from a T/Moisture input
  +    !!
  +    !! return the moist / saturated adiabatic lapse rate for a given
  +    !! Temperature and mixing ratio (really MR could be calculated as f(T))
  +    !! from http://glossary.ametsoc.org/wiki/Saturation-adiabatic_lapse_rate
  +    !!
  +    !!----------------------------------------------------------
       pure function calc_sat_lapse_rate(T,mr) result(sat_lapse)
           implicit none
  -        real, intent(in) :: T,mr  ! inputs T in K (?) and mr in kg/kg
  +        real, intent(in) :: T,mr  ! inputs T in K and mr in kg/kg
           real :: L
           real :: sat_lapse
   
  
positive.code:                     bottom = max(1, j - (vsmooth - (top-j)) )
                      
                      domain%nsquared(i,j,k) = calc_stability(domain%th(i,bottom,k), domain%th(i,top,k),  &
                                              domain%pii(i,bottom,k),domain%pii(i,top,k), &
                                              domain%z(i,bottom,k),  domain%z(i,top,k),   &
                                              domain%qv(i,bottom,k), domain%qv(i,top,k),  &
                                              domain%cloud(i,j,k)+domain%ice(i,j,k)+domain%qrain(i,j,k)+domain%qsnow(i,j,k))
                      domain%nsquared(i,j,k) = max(min(domain%nsquared(i,j,k) * nsq_calibration(i,k), max_stability), min_stability)
                  end do
              end do
  
negative.code:         u = sin(direction) * magnitude
      end function calc_u
  
      !>----------------------------------------------------------
      !! Calculate the strength of the v wind field given a direction [0-2*pi] and magnitude
      !!
      !!----------------------------------------------------------
      pure function calc_v(direction, magnitude) result(v)
          implicit none
          real, intent(in) :: direction, magnitude
  