anchor.diff:          allocate(dv(nx-2,ny-2))
           allocate(divergence(nx-2,ny-2))
           
  -! If this becomes a bottle neck in the code it could be parallelized over y  
  -!       loop over domain levels
  +        ! If this becomes a bottle neck in the code it could be parallelized over y  
  +        ! loop over domain levels
           do i=1,nz
  -            ! if we are incorporating density into the advection equation (as we should?)
  +            !------------------------------------------------------------
  +            ! If we are incorporating density into the advection equation
  +            ! then it needs to be incorporated when balancing the wind field
  +            ! 
  +            ! Note that the "else" case below does the same thing without density
  +            ! and is much easier to understand
  +            !------------------------------------------------------------
               if (options%advect_density) then
                   ! first calculate density on the u and v grid points
  -                rhou=(domain%rho(1:nx-1,i,2:ny-1) + domain%rho(2:nx,i,2:ny-1))/2
  -                rhov=(domain%rho(2:nx-1,i,1:ny-1) + domain%rho(2:nx-1,i,2:ny))/2
  +                rhou = (domain%rho(1:nx-1,i,2:ny-1) + domain%rho(2:nx,i,2:ny-1))/2
  +                rhov = (domain%rho(2:nx-1,i,1:ny-1) + domain%rho(2:nx-1,i,2:ny))/2
  +                
                   ! for most grid points interpolate between the current grid and the one above it
                   if (i<nz) then
  -                    rhow=(domain%rho(2:nx-1,i,2:ny-1) + domain%rho(2:nx-1,i+1,2:ny-1))/2
  +                    rhow = (domain%rho(2:nx-1,i,2:ny-1) + domain%rho(2:nx-1,i+1,2:ny-1))/2
                   else
                       ! for the top grid cell extrapolate upwards based on the current grid and the one below it
  -                    rhow=2*domain%rho(2:nx-1,i,2:ny-1)-domain%rho(2:nx-1,i-1,2:ny-1)
  -                endif
  +                    rhow = 2*domain%rho(2:nx-1,i,2:ny-1) - domain%rho(2:nx-1,i-1,2:ny-1)
  +                endif
  +                
                   ! calculate horizontal divergence based on the wind * density field
  -                domain%vr(2:nx-1,i,2:ny)=rhov*domain%v(2:nx-1,i,2:ny) * domain%dz_inter(1,i,1)*domain%dx
  -                dv=domain%vr(2:nx-1,i,3:ny) - domain%vr(2:nx-1,i,2:ny-1)
  -                
  -                domain%ur(2:nx,i,2:ny-1)=rhou*domain%u(2:nx,i,2:ny-1) * domain%dz_inter(1,i,1)*domain%dx
  -                du=domain%ur(3:nx,i,2:ny-1) - domain%ur(2:nx-1,i,2:ny-1)
  -                
  -                divergence=du+dv
  +                domain%vr(2:nx-1,i,2:ny) = rhov * domain%v(2:nx-1,i,2:ny) * domain%dz_inter(1,i,1) * domain%dx
  +                dv = domain%vr(2:nx-1,i,3:ny) - domain%vr(2:nx-1,i,2:ny-1)
  +                
  +                domain%ur(2:nx,i,2:ny-1) = rhou * domain%u(2:nx,i,2:ny-1) * domain%dz_inter(1,i,1) * domain%dx
  +                du = domain%ur(3:nx,i,2:ny-1) - domain%ur(2:nx-1,i,2:ny-1)
  +                
  +                divergence = du + dv
                   if (i==1) then
                       ! if this is the first model level start from 0 at the ground
  -                    domain%wr(2:nx-1,i,2:ny-1)=0-divergence
  -                    domain%w(2:nx-1,i,2:ny-1)=0-divergence/rhow/(domain%dx**2)
  +                    domain%wr(2:nx-1,i,2:ny-1) = 0 - divergence
  +                    domain%w (2:nx-1,i,2:ny-1) = 0 - divergence / rhow / (domain%dx**2)
                   else
                       ! else calculate w as a change from w at the level below
  -                    domain%wr(2:nx-1,i,2:ny-1)=domain%wr(2:nx-1,i-1,2:ny-1)-divergence
  -                    domain%w(2:nx-1,i,2:ny-1)=domain%w(2:nx-1,i-1,2:ny-1)-divergence/rhow/(domain%dx**2)
  +                    domain%wr(2:nx-1,i,2:ny-1) = domain%wr(2:nx-1,i-1,2:ny-1) - divergence
  +                    domain%w (2:nx-1,i,2:ny-1) = domain%w (2:nx-1,i-1,2:ny-1) - divergence / rhow / (domain%dx**2)
                   endif
               else
  +                !------------------------------------------------------------
  +                ! If we are not incorporating density, this is much easier to understand
  +                !------------------------------------------------------------
                   ! calculate horizontal divergence
  -                dv=domain%v(2:nx-1,i,3:ny) - domain%v(2:nx-1,i,2:ny-1)
  -                du=domain%u(3:nx,i,2:ny-1) - domain%u(2:nx-1,i,2:ny-1)
  -                divergence=du+dv
  +                dv = domain%v(2:nx-1,i,3:ny) - domain%v(2:nx-1,i,2:ny-1)
  +                du = domain%u(3:nx,i,2:ny-1) - domain%u(2:nx-1,i,2:ny-1)
  +                divergence = du + dv
  +                
  +                ! Then calculate w to balance
                   if (i==1) then
                       ! if this is the first model level start from 0 at the ground
  -                    domain%w(2:nx-1,i,2:ny-1)=0-divergence
  +                    domain%w(2:nx-1,i,2:ny-1) = 0 - divergence
                   else
                       ! else calculate w as a change from w at the level below
  -                    domain%w(2:nx-1,i,2:ny-1)=domain%w(2:nx-1,i-1,2:ny-1)-divergence
  +                    domain%w(2:nx-1,i,2:ny-1) = domain%w(2:nx-1,i-1,2:ny-1) - divergence
                   endif
               endif
           enddo
  -        ! NOTE w is scaled by dx/dz because it is balancing divergence through a dx*dz cell face with a flow through a dx*dx face
  -        ! could rescale below but for now this is left in the advection code (the only place it matters for now)
  -        ! this makes it easier to play with varying options for including varying dz and rho in advection. 
  -        ! domain%w=domain%w/domain%dx * domain%dz_inter
  -        
  -        deallocate(du,dv,divergence)
  -        if (options%advect_density) then
  -            deallocate(rhou,rhov,rhow)
  -        endif
  +        
       end subroutine balance_uvw
       
  -!   Correct for a grid that is locally rotated with respect to EW,NS (e.g. at the edges of the domain)
  -!   Assumes forcing winds are EW,NS relative, not grid relative. 
  +    !>------------------------------------------------------------
  +    !! Correct for a grid that is locally rotated with respect to EW,NS
  +    !!
  +    !! Assumes forcing winds are EW, NS relative, not grid relative.
  +    !! 
  +    !!------------------------------------------------------------
       subroutine make_winds_grid_relative(domain)
           type(domain_type), intent(inout) :: domain
           real,dimension(:),allocatable :: u,v
  
positive.code:         
      end subroutine allocate_winds
  
  ! deallocate memory
      subroutine finalize_winds(domain)
          type(domain_type), intent(inout) :: domain
          
  
negative.code:         nz=size(domain%w,2)
          ny=size(domain%w,3)
          
          !------------------------------------------------------------
          ! These could be module level to prevent lots of allocation/deallocation/reallocations
          ! but these are relatively small, and should be allocated efficiently on the stack
          !------------------------------------------------------------
          if (options%advect_density) then
              allocate(rhou(nx-1,ny-2))
              allocate(rhov(nx-2,ny-1))
  