anchor.diff:                          ! 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
  -                        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)))
  +                        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
  
positive.code:         maxwind1d = 0
          maxwind3d = 0
          
          nx = size(w,1)
          nz = size(w,2)
          ny = size(w,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)))
              
              maxwind3d = maxwind1d * sqrt3
          else if (cfl_strictness==5) then
              maxwind3d = maxval(abs(u)) + maxval(abs(v)) + maxval(abs(w))
          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
  
negative.code:                 maxwind3d = maxwind3d * three_d_cfl
                  
                  ! to ensure we are stable for 1D advection:
                  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)
  
  