anchor.diff:          maxwind1d = 0
           maxwind3d = 0
           
  -        nx = size(w,1)
  -        nz = size(w,2)
  -        ny = size(w,3)
  +        nx = size(rho,1)
  +        nz = size(rho,2)
  +        ny = size(rho,3)
   
           if (cfl_strictness==1) then
               ! 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
  +                        
               maxwind3d = maxwind1d * sqrt3
           else if (cfl_strictness==5) then
  -            maxwind3d = maxval(abs(u)) + maxval(abs(v)) + maxval(abs(w))
  +            
  +            if (use_density) then
  +                maxwind1d = maxval(abs(u(2:,:,:) / (rho*dz*dx) )) &
  +                          + maxval(abs(v(:,:,2:) / (rho*dz*dx) )) &
  +                          + maxval(abs(w(:,:, :) / (rho*dz*dx) ))
  +            else
  +                maxwind3d = maxval(abs(u)) + maxval(abs(v)) + maxval(abs(w))
  +            endif
  +
           else
               ! to ensure we are stable for 3D advection we'll use the average "max" wind speed
               ! but that average has to be divided by sqrt(3) for stability in 3 dimensional advection
  
positive.code:     !! @return dt [ scalar ]        Maximum stable time step    [s]
      !!
      !!------------------------------------------------------------
      function compute_dt(dx, u, v, w, CFL, cfl_strictness) result(dt)
          real,       intent(in)                   :: dx
          real,       intent(in), dimension(:,:,:) :: u, v, w
          real,       intent(in)                   :: CFL
          integer,    intent(in)                   :: cfl_strictness
          ! output value
          real :: dt
          ! locals
  
negative.code:                         ! just compute the sum of the wind speeds, but take the max over the two
                          ! faces of the grid cell (e.g. east and west sides)
                          ! this will be divided by 3 later by three_d_cfl
                          if (use_density) then
                              current_wind = (max(abs(u(i,k,j)), abs(u(i+1,k,j))) &
                                            + max(abs(v(i,k,j)), abs(v(i,k,j+1))) & 
                                            + max(abs(w(i,k,j)), abs(w(i,k+zoffset,j))) ) &
                                            / (rho(i,k,j) * dz(i,k,j) * dx)
                          else
                              current_wind = max(abs(u(i,k,j)), abs(u(i+1,k,j))) &
                                            +max(abs(v(i,k,j)), abs(v(i,k,j+1))) & 
                                            +max(abs(w(i,k,j)), abs(w(i,k+zoffset,j)))
                          endif
                          maxwind3d = max(maxwind3d, current_wind)
                      ENDDO
                  ENDDO
  