anchor.diff:                  ! 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 (debug) then
  -                        !$omp critical (print_lock)
  -                        write(*,*) "Using a density correction in linear winds"
  -                        !$omp end critical (print_lock)
  -                    endif
  -                    u_hat(buffer:realnx_u+buffer,buffer:realny+buffer,z) = &
  -                        2*real(u_hat(buffer:realnx_u+buffer,buffer:realny+buffer,z))! / domain%rho(1:realnx,z,1:realny)
  -                    v_hat(buffer:realnx+buffer,buffer:realny_v+buffer,z) = &
  -                        2*real(v_hat(buffer:realnx+buffer,buffer:realny_v+buffer,z))! / domain%rho(1:realnx,z,1:realny)
  -                endif
  +                ! if (useDensity) then
  +                !     if (debug) then
  +                !         !$omp critical (print_lock)
  +                !         write(*,*) "Using a density correction in linear winds"
  +                !         !$omp end critical (print_lock)
  +                !     endif
  +                !     u_hat(buffer:realnx_u+buffer,buffer:realny+buffer,z) = &
  +                !         2*real(u_hat(buffer:realnx_u+buffer,buffer:realny+buffer,z))! / domain%rho(1:realnx,z,1:realny)
  +                !     v_hat(buffer:realnx+buffer,buffer:realny_v+buffer,z) = &
  +                !         2*real(v_hat(buffer:realnx+buffer,buffer:realny_v+buffer,z))! / domain%rho(1:realnx,z,1:realny)
  +                ! endif
                   
                   ! If we are removing linear winds from a low res field, subtract u_hat v_hat instead
                   ! real(real()) extracts real component of complex, then converts to a real data type (may not be necessary except for IO?)
  
positive.code:                         v = domain%v(vi,j,k)
                      endif
                      n = n * ((top-bottom)+1)
  
                      ! Calculate the direction of the current grid cell wind
                      dpos = 1
                      curdir = calc_direction( u, v )
  
negative.code:                                     + linear_update_fraction * (sweight*wind_first + (1-sweight)*wind_second)
  
                          if (reverse) then
                              domain%v(i,j,k) = domain%v(i,j,k) - v_perturbation(i,j,k) * linear_contribution
                          else
                              ! for the high res domain, linear_mask should incorporate linear_contribution
                              domain%v(i,j,k) = domain%v(i,j,k) + v_perturbation(i,j,k) * linear_mask(min(nx,i),min(ny,k))
                          endif
                      endif
                  end do
  