anchor.diff:          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 this only gets called once for each forcing input
  +        !------------------------------------------------------------
  +        ! 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))
  
positive.code:         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
          do i=1,nz
              ! if we are incorporating density into the advection equation (as we should?)
              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
                  ! 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
                  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
                  ! 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
                  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)
                  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)
                  endif
              else
                  ! 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
                  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
                  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
                  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. 
      subroutine make_winds_grid_relative(domain)
          type(domain_type), intent(inout) :: domain
          real,dimension(:),allocatable :: u,v
  
negative.code:         
      end subroutine allocate_winds
  
      !>------------------------------------------------------------
      !! Provides a routine to deallocate memory allocated in allocate_winds
      !!
      !!------------------------------------------------------------
      subroutine finalize_winds(domain)
          type(domain_type), intent(inout) :: domain
          
  