Solid-Liquid Leaching & Washing

program solid_liquid_leaching
    implicit none
    integer :: mode, N, i, reqStages, iostat_val
    double precision :: S, X0, targetRecovery, V, Uinput, U, m, washRatio, E, residualFrac, recovery
    double precision :: xSolid, recovered, loss, extract, totalSolute, soluteRecovered, raffLoss

    read(*,*,iostat=iostat_val) mode
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid mode input.'
        stop
    end if
    read(*,*,iostat=iostat_val) S
    read(*,*,iostat=iostat_val) X0
    read(*,*,iostat=iostat_val) targetRecovery
    read(*,*,iostat=iostat_val) N
    read(*,*,iostat=iostat_val) V
    read(*,*,iostat=iostat_val) Uinput
    read(*,*,iostat=iostat_val) m
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all leaching inputs.'
        stop
    end if
    if (S <= 0.0d0 .or. X0 < 0.0d0 .or. V < 0.0d0 .or. Uinput <= 0.0d0 .or. m <= 0.0d0) then
        write(*,*) 'ERROR: Solid mass, solvent, underflow and equilibrium factor must be valid.'
        stop
    end if
    if (N < 1) N = 1
    if (targetRecovery < 0.0d0 .or. targetRecovery >= 100.0d0) then
        write(*,*) 'ERROR: Target recovery must be between 0 and less than 100 percent.'
        stop
    end if

    if (Uinput > 10.0d0) then
        U = Uinput
    else
        U = Uinput*S
    end if
    totalSolute = S*X0
    washRatio = V/max(1.0d-30,U)
    E = washRatio/m

    if (mode == 1) then
        residualFrac = 1.0d0/(1.0d0 + E)
    else
        residualFrac = kremser_residual(E,N)
    end if
    recovery = (1.0d0 - residualFrac)*100.0d0
    xSolid = X0*residualFrac
    raffLoss = totalSolute*residualFrac
    soluteRecovered = totalSolute - raffLoss

    reqStages = 1
    do i=1,200
        if (mode == 1) then
            residualFrac = 1.0d0/(1.0d0 + E)
        else
            residualFrac = kremser_residual(E,i)
        end if
        if ((1.0d0-residualFrac)*100.0d0 >= targetRecovery) then
            reqStages = i
            exit
        end if
        reqStages = i
    end do

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   SOLID-LIQUID LEACHING ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,I8)')       '  Mode                     = ', mode
    write(*,'(A,ES12.4,A)') '  Dry Insoluble Solids     = ', S, ' kg'
    write(*,'(A,ES12.4,A)') '  Initial Loading X0       = ', X0, ' kg/kg dry solid'
    write(*,'(A,ES12.4,A)') '  Fresh Solvent V          = ', V, ' kg'
    write(*,'(A,ES12.4,A)') '  Retained Underflow U     = ', U, ' kg solution'
    write(*,'(A,ES12.4)')   '  Equilibrium Factor m     = ', m
    write(*,'(A,I8)')       '  Specified Stages         = ', N
    write(*,*)
    write(*,'(A)') '--- DESIGN RESULTS ------------------------------------------'
    write(*,'(A,ES12.4)')   '  Wash Ratio               = ', washRatio
    write(*,'(A,ES12.4)')   '  Kremser Factor           = ', E
    write(*,'(A,ES12.4,A)') '  Solute Recovery          = ', recovery, ' percent'
    write(*,'(A,I8)')       '  Required Stage Count     = ', reqStages
    write(*,'(A,ES12.4,A)') '  Final Solid Loading      = ', xSolid, ' kg/kg dry solid'
    write(*,'(A,ES12.4,A)') '  Solute Recovered         = ', soluteRecovered, ' kg'
    write(*,'(A,ES12.4,A)') '  Raffinate Loss           = ', raffLoss, ' kg'
    write(*,*)
    write(*,'(A)') '--- STAGE PROFILE -------------------------------------------'
    write(*,'(A)') '  stage        X_solid       recovery[%]   raff_loss     extract_solute'
    write(*,'(A)') '  ---------------------------------------------------------------------'
    do i=1,max(N,reqStages)
        if (mode == 1) then
            residualFrac = 1.0d0/(1.0d0 + E)
        else
            residualFrac = kremser_residual(E,i)
        end if
        xSolid = X0*residualFrac
        recovered = (1.0d0-residualFrac)*100.0d0
        loss = totalSolute*residualFrac
        extract = totalSolute-loss
        write(*,'(I8,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') i, xSolid, recovered, loss, extract
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Wash ratio W = V/U. Kremser extraction factor E = W/m.'
    write(*,'(A)') '  Single stage residual fraction = 1/(1+E).'
    write(*,'(A)') '  Countercurrent residual fraction = (E-1)/(E^(N+1)-1), with E=1 limit = 1/(N+1).'

contains
    double precision function kremser_residual(Ein, Nin)
        implicit none
        double precision, intent(in) :: Ein
        integer, intent(in) :: Nin
        if (abs(Ein-1.0d0) < 1.0d-10) then
            kremser_residual = 1.0d0/dble(Nin+1)
        else if (Ein > 0.0d0) then
            kremser_residual = (Ein-1.0d0)/(Ein**dble(Nin+1)-1.0d0)
            if (kremser_residual < 0.0d0) kremser_residual = abs(kremser_residual)
            if (kremser_residual > 1.0d0) kremser_residual = 1.0d0
        else
            kremser_residual = 1.0d0
        end if
    end function kremser_residual
end program solid_liquid_leaching