anchor.diff:                  ! W_m = (domain%w_cu + domain%w)                     * (dt/dx)
                   ! call rebalance_cu_winds(U_m,V_m,W_m)
               else
  -                U_m = u(2:nx,:,:) * (dt/dx)
  -                V_m = v(:,:,2:ny) * (dt/dx)
  -                ! note, even though dz!=dx, W is computed from the divergence in U/V so it is scaled by dx/dz already
  -                W_m = w   * (dt/dx)
  +                U_m = u(2:nx,:,:) * ((dz(1:nx-1,:,:)+dz(2:nx,:,:))/2 * dt/dx)
  +                V_m = v(:,:,2:ny) * ((dz(:,:,1:ny-1)+dz(:,:,2:ny))/2 * dt/dx)
  +                W_m = w * dt
               endif
           ! endif
   
  
positive.code: !
  !     end subroutine flux2
  
      subroutine advect3d(q,u,v,w,rho,dz,nx,nz,ny,options)
          implicit none
          real, dimension(1:nx,  1:nz,1:ny),  intent(inout)   :: q
          real, dimension(1:nx,  1:nz,1:ny),  intent(in)      :: w
          real, dimension(1:nx-1,1:nz,1:ny),  intent(in)      :: u
          real, dimension(1:nx,  1:nz,1:ny-1),intent(in)      :: v
          real, dimension(:, :, :),  intent(in)      :: rho
          real, dimension(:, :, :),  intent(in)      :: dz
          integer,                            intent(in)      :: ny, nz, nx
          type(options_t), intent(in)::options
  
  
negative.code:     !
      ! end subroutine advect_cu_winds
  
      subroutine setup_module_winds(u,v,w, dz, dx, options, dt)
          implicit none
          real,               intent(in)  :: u(:,:,:)
          real,               intent(in)  :: v(:,:,:)
          real,               intent(in)  :: w(:,:,:)
          real,               intent(in)  :: dz(:,:,:)
          real,               intent(in)  :: dx
          type(options_t),    intent(in)  :: options
          real,               intent(in)  :: dt
  