anchor.diff:                  maxwind3d = maxwind3d * three_d_cfl
                   
                   ! to ensure we are stable for 1D advection:
  -                maxwind1d = max( maxval(abs(u)), maxval(abs(v)))
  -                maxwind1d = max( maxwind1d, maxval(abs(w)))
  +                if (use_density) then
  +                    maxwind1d = max( maxval(abs(u(2:,:,:) / (rho*dz*dx) )), maxval(abs(v(:,:,2:) / (rho*dz*dx))) )
  +                    maxwind1d = max( maxwind1d, maxval(abs(w/(rho*dz*dx))) )
  +                else
  +                    maxwind1d = max( maxval(abs(u)), maxval(abs(v)))
  +                    maxwind1d = max( maxwind1d, maxval(abs(w)))
  +                endif
                   ! also insure stability for 1D advection
                   maxwind3d = max(maxwind1d,maxwind3d)
   
  
positive.code:             ! compute internal timestep dt to maintain stability
              ! courant condition for 3D advection. Note that w is normalized by dx/dz
              ! pick the minimum dt from the begining or the end of the current timestep
              dt = compute_dt(domain%dx, domain%u, domain%v, domain%w, CFL, cfl_strictness=options%cfl_strictness)
              ! set an upper bound on dt to keep microphysics and convection stable (?) not sure what time is required here. 
              dt = min(dt,120.0) !better min=180?
              if (options%debug) write(*,"(A,f5.1,A,f5.1,A$)") char(13), 100-max(0.0,(end_time-model_time-dt)/options%in_dt*100)," %  dt=",dt,"s  "
  
negative.code:     !! @return dt [ scalar ]        Maximum stable time step    [s]
      !!
      !!------------------------------------------------------------
      function compute_dt(dx, u, v, w, rho, dz, CFL, cfl_strictness, use_density) result(dt)
          real,       intent(in)                   :: dx
          real,       intent(in), dimension(:,:,:) :: u, v, w, rho, dz
          real,       intent(in)                   :: CFL
          integer,    intent(in)                   :: cfl_strictness
          logical,    intent(in)                   :: use_density
          ! output value
          real :: dt
          ! locals
  