anchor.diff:             !                          / rho(2:nx-1,nz,i) / dz(2:nx-1,nz,i)
              ! else
                  ! perform horizontal advection
  -               q(2:nx-1,:,i)      = q(2:nx-1,:,i)       - ((f1(2:nx-1,:) - f1(1:nx-2,:)) + (f3 - f4))
  +               q(2:nx-1,:,i)      = q(2:nx-1,:,i)       - ((f1(2:nx-1,:) - f1(1:nx-2,:)) + (f3 - f4)) / (dx * dz(2:nx-1,:,i))
                  ! then vertical (order doesn't matter because fluxes f1-6 are calculated before applying them)
                  ! add fluxes to middle layers
  -               q(2:nx-1,2:nz-1,i) = q(2:nx-1,2:nz-1,i)  - (f5(:,2:nz-1) - f5(:,1:nz-2))
  +               q(2:nx-1,2:nz-1,i) = q(2:nx-1,2:nz-1,i)  - (f5(:,2:nz-1) - f5(:,1:nz-2)) / dz(2:nx-1,2:nz-1,i)
                  ! add fluxes to bottom layer
  -               q(2:nx-1,1,i)      = q(2:nx-1,1,i)       - f5(:,1)
  +               q(2:nx-1,1,i)      = q(2:nx-1,1,i)       - f5(:,1) / dz(2:nx-1,1,i)
                  ! add fluxes to top layer
  -               q(2:nx-1,nz,i)     = q(2:nx-1,nz,i)      - (qin(2:nx-1,nz,i) * w(2:nx-1,nz,i) - f5(:,nz-1))
  +               q(2:nx-1,nz,i)     = q(2:nx-1,nz,i)      - (qin(2:nx-1,nz,i) * w(2:nx-1,nz,i) - f5(:,nz-1)) / dz(2:nx-1,nz,i)
              ! endif
           enddo
           !$omp end do
  
positive.code:                 ! 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)
              endif
          ! endif
  
  
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
  