anchor.diff:          ny = domain%grid%ny
   
           ! calculate U,V,W normalized for dt/dx (dx**2 for density advection so we can skip a /dx in the actual advection code)
  -        call setup_module_winds(domain%u%data_3d, domain%v%data_3d, domain%w%data_3d, domain%dx, options, dt)
  +        call setup_module_winds(domain%u%data_3d, domain%v%data_3d, domain%w%data_3d, domain%dz_interface%data_3d, domain%dx, options, dt)
   
           ! lastqv_m=domain%qv
   
  -        if (options%vars_to_advect(kVARS%water_vapor)>0)                  call advect3d(domain%water_vapor%data_3d,             U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, nx,nz,ny, options)
  -        if (options%vars_to_advect(kVARS%cloud_water)>0)                  call advect3d(domain%cloud_water_mass%data_3d,        U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, nx,nz,ny, options)
  -        if (options%vars_to_advect(kVARS%rain_in_air)>0)                  call advect3d(domain%rain_mass%data_3d,               U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, nx,nz,ny, options)
  -        if (options%vars_to_advect(kVARS%snow_in_air)>0)                  call advect3d(domain%snow_mass%data_3d,               U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, nx,nz,ny, options)
  -        if (options%vars_to_advect(kVARS%potential_temperature)>0)        call advect3d(domain%potential_temperature%data_3d,   U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, nx,nz,ny, options)
  -        if (options%vars_to_advect(kVARS%cloud_ice)>0)                    call advect3d(domain%cloud_ice_mass%data_3d,          U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, nx,nz,ny, options)
  -        if (options%vars_to_advect(kVARS%graupel_in_air)>0)               call advect3d(domain%graupel_mass%data_3d,            U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, nx,nz,ny, options)
  -        if (options%vars_to_advect(kVARS%ice_number_concentration)>0)     call advect3d(domain%cloud_ice_number%data_3d,        U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, nx,nz,ny, options)
  -        if (options%vars_to_advect(kVARS%rain_number_concentration)>0)    call advect3d(domain%rain_number%data_3d,             U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, nx,nz,ny, options)
  -        if (options%vars_to_advect(kVARS%snow_number_concentration)>0)    call advect3d(domain%snow_number%data_3d,             U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, nx,nz,ny, options)
  -        if (options%vars_to_advect(kVARS%graupel_number_concentration)>0) call advect3d(domain%graupel_number%data_3d,          U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, nx,nz,ny, options)
  +        if (options%vars_to_advect(kVARS%water_vapor)>0)                  call advect3d(domain%water_vapor%data_3d,             U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, domain%dx, nx,nz,ny, options)
  +        if (options%vars_to_advect(kVARS%cloud_water)>0)                  call advect3d(domain%cloud_water_mass%data_3d,        U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, domain%dx, nx,nz,ny, options)
  +        if (options%vars_to_advect(kVARS%rain_in_air)>0)                  call advect3d(domain%rain_mass%data_3d,               U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, domain%dx, nx,nz,ny, options)
  +        if (options%vars_to_advect(kVARS%snow_in_air)>0)                  call advect3d(domain%snow_mass%data_3d,               U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, domain%dx, nx,nz,ny, options)
  +        if (options%vars_to_advect(kVARS%potential_temperature)>0)        call advect3d(domain%potential_temperature%data_3d,   U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, domain%dx, nx,nz,ny, options)
  +        if (options%vars_to_advect(kVARS%cloud_ice)>0)                    call advect3d(domain%cloud_ice_mass%data_3d,          U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, domain%dx, nx,nz,ny, options)
  +        if (options%vars_to_advect(kVARS%graupel_in_air)>0)               call advect3d(domain%graupel_mass%data_3d,            U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, domain%dx, nx,nz,ny, options)
  +        if (options%vars_to_advect(kVARS%ice_number_concentration)>0)     call advect3d(domain%cloud_ice_number%data_3d,        U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, domain%dx, nx,nz,ny, options)
  +        if (options%vars_to_advect(kVARS%rain_number_concentration)>0)    call advect3d(domain%rain_number%data_3d,             U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, domain%dx, nx,nz,ny, options)
  +        if (options%vars_to_advect(kVARS%snow_number_concentration)>0)    call advect3d(domain%snow_number%data_3d,             U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, domain%dx, nx,nz,ny, options)
  +        if (options%vars_to_advect(kVARS%graupel_number_concentration)>0) call advect3d(domain%graupel_number%data_3d,          U_m,V_m,W_m, domain%density%data_3d, domain%dz_interface%data_3d, domain%dx, nx,nz,ny, options)
   
           ! if (options%physics%convection > 0) then
           !     call advect_cu_winds(domain, options, dt)
  
positive.code:     !
      ! end subroutine advect_cu_winds
  
      subroutine setup_module_winds(u,v,w, dx, options, dt)
          implicit none
          real,               intent(in)  :: u(:,:,:)
          real,               intent(in)  :: v(:,:,:)
          real,               intent(in)  :: w(:,:,:)
          real,               intent(in)  :: dx
          type(options_t),    intent(in)  :: options
          real,               intent(in)  :: dt
  
negative.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,:,:) * ((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
  
  