anchor.diff:                      i=ny
                       u_hat(1:nx-1,i,z) = (u_hat(1:nx-1,i,z) + u_hat(2:nx,i,z)) /2
                   endif
  +                
  +                ! If we are using density in the advection calculations, modify the linear perturbation
  +                ! to get the vertical velocities closer to what they would be without density (boussinesq)
  +                ! need to check if this makes the most sense when close to the surface
                   if (useDensity) then
  -                    ! if we are using density in the advection calculations, modify the linear perturbation
  -                    ! to get the vertical velocities closer to what they would be without density (boussinesq)
  -                    ! need to check if this makes the most sense when close to the surface
                       if (debug) then
                           write(*,*) "Using a density correction in linear winds"
                       endif
  
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: 
      end function calc_stability
  
      !>----------------------------------------------------------
      !! Calculate a single brunt vaisala frequency effectively averaged across the entire domain
      !!
      !!----------------------------------------------------------
      pure function calc_domain_stability(domain) result(BV_freq)
          implicit none
          class(linearizable_type),intent(in)::domain
  