program cfl_criterion
    implicit none
    double precision :: u,v,w,c_s,dx,dy,dz,dt_in,alph,nu_v,gam
    integer :: scheme, ios, ndim, i, npts
    double precision :: Vmag, Mach
    double precision :: cfl_conv, cfl_acou, cfl_tot
    double precision :: d_th, d_vi
    double precision :: inv_dx2, cfl_max_scheme
    double precision :: dt_cfl, dt_diff, dt_visc, dt_comb, dt_use
    double precision :: Re_cell, Pe_cell, combined
    character(len=30) :: regime, sch_name, stab_type
    double precision :: tc, dc, vc, cc, st_val

    read(*,*,iostat=ios) u;      if(ios/=0)then;write(*,*)'ERROR: Invalid u.';stop;end if
    read(*,*,iostat=ios) v;      if(ios/=0)then;write(*,*)'ERROR: Invalid v.';stop;end if
    read(*,*,iostat=ios) w;      if(ios/=0)then;write(*,*)'ERROR: Invalid w.';stop;end if
    read(*,*,iostat=ios) c_s;    if(ios/=0)then;write(*,*)'ERROR: Invalid c_sound.';stop;end if
    read(*,*,iostat=ios) dx;     if(ios/=0)then;write(*,*)'ERROR: Invalid dx.';stop;end if
    read(*,*,iostat=ios) dy;     if(ios/=0)then;write(*,*)'ERROR: Invalid dy.';stop;end if
    read(*,*,iostat=ios) dz;     if(ios/=0)then;write(*,*)'ERROR: Invalid dz.';stop;end if
    read(*,*,iostat=ios) dt_in;  if(ios/=0)then;write(*,*)'ERROR: Invalid dt.';stop;end if
    read(*,*,iostat=ios) alph;   if(ios/=0)then;write(*,*)'ERROR: Invalid alpha.';stop;end if
    read(*,*,iostat=ios) nu_v;   if(ios/=0)then;write(*,*)'ERROR: Invalid nu.';stop;end if
    read(*,*,iostat=ios) gam;    if(ios/=0)then;write(*,*)'ERROR: Invalid gamma.';stop;end if
    read(*,*,iostat=ios) scheme; if(ios/=0)then;write(*,*)'ERROR: Invalid scheme.';stop;end if

    if(dx<=0d0)then;write(*,*)'ERROR: dx must be > 0.';stop;end if

    ! Determine dimensionality
    ndim = 1
    if (dy > 0d0) ndim = 2
    if (dy > 0d0 .and. dz > 0d0) ndim = 3
    if (dy <= 0d0) dy = 1d30
    if (dz <= 0d0) dz = 1d30

    Vmag = sqrt(u**2 + v**2 + w**2)
    if (c_s > 0d0) then
        Mach = Vmag / c_s
    else
        Mach = 0d0
        c_s = 0d0
    end if

    ! Scheme info
    select case(scheme)
    case(1)
        sch_name = 'Euler Explicit'
        cfl_max_scheme = 1.0d0
        stab_type = 'Conditionally Stable'
    case(2)
        sch_name = 'Leapfrog'
        cfl_max_scheme = 1.0d0
        stab_type = 'Conditionally Stable'
    case(3)
        sch_name = 'Runge-Kutta 4'
        cfl_max_scheme = 2.83d0
        stab_type = 'Conditionally Stable'
    case(4)
        sch_name = 'Implicit (1st order)'
        cfl_max_scheme = 1d10
        stab_type = 'Unconditionally Stable'
    case default
        sch_name = 'Euler Explicit'
        cfl_max_scheme = 1.0d0
        stab_type = 'Conditionally Stable'
    end select

    ! Inverse grid spacing squared sum
    inv_dx2 = 1d0/dx**2
    if (ndim >= 2) inv_dx2 = inv_dx2 + 1d0/dy**2
    if (ndim >= 3) inv_dx2 = inv_dx2 + 1d0/dz**2

    ! Max stable time steps
    ! CFL-based
    if (Vmag + c_s > 0d0) then
        dt_cfl = cfl_max_scheme * dx / (Vmag + c_s)
    else
        dt_cfl = 1d30
    end if

    ! Diffusion-based
    if (alph > 0d0) then
        dt_diff = 0.5d0 / (alph * inv_dx2)
    else
        dt_diff = 1d30
    end if

    ! Viscous-based
    if (nu_v > 0d0) then
        dt_visc = 0.5d0 / (nu_v * inv_dx2)
    else
        dt_visc = 1d30
    end if

    dt_comb = min(dt_cfl, dt_diff, dt_visc)

    ! Determine dt to use
    if (dt_in > 0d0) then
        dt_use = dt_in
    else
        dt_use = dt_comb * 0.9d0
    end if

    ! CFL numbers
    cfl_conv = abs(u)*dt_use/dx
    if (ndim >= 2) cfl_conv = cfl_conv + abs(v)*dt_use/dy
    if (ndim >= 3) cfl_conv = cfl_conv + abs(w)*dt_use/dz

    if (c_s > 0d0) then
        cfl_acou = (abs(u)+c_s)*dt_use/dx
        if (ndim >= 2) cfl_acou = cfl_acou + (abs(v)+c_s)*dt_use/dy
        if (ndim >= 3) cfl_acou = cfl_acou + (abs(w)+c_s)*dt_use/dz
    else
        cfl_acou = cfl_conv
    end if

    cfl_tot = cfl_acou

    ! Diffusion numbers
    d_th = alph * dt_use * inv_dx2
    d_vi = nu_v * dt_use * inv_dx2

    ! Combined
    combined = cfl_conv + 2d0*max(d_th, d_vi)

    ! Cell numbers
    if (nu_v > 0d0) then
        Re_cell = abs(u)*dx/nu_v
    else
        Re_cell = 0d0
    end if
    if (alph > 0d0) then
        Pe_cell = abs(u)*dx/alph
    else
        Pe_cell = 0d0
    end if

    ! Stability status
    if (scheme == 4) then
        regime = 'STABLE (Implicit)'
    else if (cfl_tot <= cfl_max_scheme .and. d_th <= 0.5d0 .and. d_vi <= 0.5d0) then
        regime = 'STABLE'
    else
        regime = 'UNSTABLE'
    end if

    ! ---- Output ------------------------------------------------------
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   CFL (COURANT-FRIEDRICHS-LEWY) CRITERION CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)

    write(*,'(A)') '--- INPUT CONDITIONS ----------------------------------------'
    write(*,'(A,F14.6,A)') '  u (x-velocity)          = ', u, ' m/s'
    write(*,'(A,F14.6,A)') '  v (y-velocity)          = ', v, ' m/s'
    write(*,'(A,F14.6,A)') '  w (z-velocity)          = ', w, ' m/s'
    write(*,'(A,F14.4,A)') '  c (speed of sound)      = ', c_s, ' m/s'
    write(*,'(A,F14.6,A)') '  |V| (magnitude)         = ', Vmag, ' m/s'
    write(*,'(A,F14.6)')   '  Mach number             = ', Mach
    write(*,*)

    write(*,'(A)') '--- GRID INFORMATION ----------------------------------------'
    write(*,'(A,ES14.6,A)') '  dx                      = ', dx, ' m'
    if (ndim >= 2) write(*,'(A,ES14.6,A)') '  dy                      = ', dy, ' m'
    if (ndim >= 3) write(*,'(A,ES14.6,A)') '  dz                      = ', dz, ' m'
    write(*,'(A,I2,A)')     '  Dimensions              = ', ndim, 'D'
    if (ndim >= 2) write(*,'(A,F14.4)')   '  Aspect ratio dy/dx      = ', dy/dx
    write(*,*)

    write(*,'(A)') '--- SCHEME INFORMATION --------------------------------------'
    write(*,'(A,A)')        '  Scheme                  = ', trim(sch_name)
    write(*,'(A,F14.4)')    '  CFL_max (scheme)        = ', min(cfl_max_scheme, 999d0)
    write(*,'(A,A)')        '  Stability type          = ', trim(stab_type)
    write(*,*)

    write(*,'(A)') '--- TIME STEP -----------------------------------------------'
    write(*,'(A,ES14.6,A)') '  dt (used)               = ', dt_use, ' s'
    write(*,'(A,ES14.6,A)') '  dt_max (CFL limit)      = ', dt_cfl, ' s'
    write(*,'(A,ES14.6,A)') '  dt_max (diffusion)      = ', dt_diff, ' s'
    write(*,'(A,ES14.6,A)') '  dt_max (viscous)        = ', dt_visc, ' s'
    write(*,'(A,ES14.6,A)') '  dt_max (combined)       = ', dt_comb, ' s'
    write(*,*)

    write(*,'(A)') '--- CFL NUMBERS ---------------------------------------------'
    write(*,'(A,F14.6)')    '  CFL (convective)        = ', cfl_conv
    write(*,'(A,F14.6)')    '  CFL (acoustic)          = ', cfl_acou
    write(*,'(A,F14.6)')    '  CFL_max (scheme limit)  = ', min(cfl_max_scheme, 999d0)
    write(*,'(A,F14.2,A)')  '  CFL margin              = ', (cfl_max_scheme - cfl_tot)/cfl_max_scheme*100d0, ' %'
    write(*,*)

    write(*,'(A)') '--- DIFFUSION NUMBERS (Fourier) -----------------------------'
    write(*,'(A,F14.6)')    '  d_thermal = alpha*dt/dx2 = ', d_th
    write(*,'(A,F14.6)')    '  d_viscous = nu*dt/dx2    = ', d_vi
    write(*,'(A,F14.6)')    '  Limit (explicit)         = 0.500000'
    write(*,*)

    write(*,'(A)') '--- COMBINED STABILITY --------------------------------------'
    write(*,'(A,F14.6)')    '  CFL + 2*d_max           = ', combined
    write(*,'(A,F14.6)')    '  Limit                    = 1.000000'
    write(*,*)

    write(*,'(A)') '--- CELL NUMBERS --------------------------------------------'
    write(*,'(A,F14.4)')    '  Re_cell = |u|*dx/nu     = ', Re_cell
    write(*,'(A,F14.4)')    '  Pe_cell = |u|*dx/alpha   = ', Pe_cell
    if (Re_cell > 2d0) then
        write(*,'(A)')      '  WARNING: Re_cell > 2, may need upwinding'
    end if
    write(*,*)

    write(*,'(A)') '--- STABILITY STATUS ----------------------------------------'
    write(*,'(A,A)')        '  STATUS                   = ', trim(regime)
    write(*,*)

    ! ---- Profile sweep -----------------------------------------------
    write(*,'(A)') '--- PROFILE vs dt ------------------------------------------'
    write(*,'(A)') '  dt          CFL_conv    CFL_acou    d_thermal   d_viscous   CFL+2d      Status'
    write(*,'(A)') '  --------------------------------------------------------------------------------------'

    npts = 40
    do i = 1, npts
        tc = dt_comb/10d0 + (dt_comb*3d0 - dt_comb/10d0)*dble(i-1)/dble(npts-1)
        if (tc <= 0d0) tc = 1d-15

        cc = abs(u)*tc/dx
        if (ndim >= 2) cc = cc + abs(v)*tc/dy
        if (ndim >= 3) cc = cc + abs(w)*tc/dz

        if (c_s > 0d0) then
            st_val = (abs(u)+c_s)*tc/dx
            if (ndim >= 2) st_val = st_val + (abs(v)+c_s)*tc/dy
            if (ndim >= 3) st_val = st_val + (abs(w)+c_s)*tc/dz
        else
            st_val = cc
        end if

        dc = alph * tc * inv_dx2
        vc = nu_v * tc * inv_dx2

        write(*,'(ES12.4,2X,F10.6,2X,F10.6,2X,F10.6,2X,F10.6,2X,F10.6,2X,A)') &
            tc, cc, st_val, dc, vc, cc+2d0*max(dc,vc), &
            merge('OK    ','UNSTAB', (st_val<=cfl_max_scheme .and. dc<=0.5d0 .and. vc<=0.5d0))
    end do

    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') '  CFL = |u|*dt/dx + |v|*dt/dy + |w|*dt/dz  (Eq. 2.86)'
    write(*,'(A)') '  CFL_acoustic = (|u|+c)*dt/dx'
    write(*,'(A)') '  d = alpha*dt/dx^2  (diffusion/Fourier number)'
    write(*,'(A)') '  Combined: CFL + 2d <= 1'
    write(*,'(A)') '  dt_cfl = CFL_max * dx / (|u|+c)'
    write(*,'(A)') '  dt_diff = 0.5 / (alpha * sum(1/dxi^2))'
    write(*,'(A)') '  Re_cell = |u|*dx/nu, Pe_cell = |u|*dx/alpha'
    write(*,'(A)') '============================================================'

end program cfl_criterion