Drying Rate Calculator

program drying_rate
    implicit none
    integer :: mode, i, n_points, iostat_val
    double precision :: Ms, A, X0, Xc, Xe, Xf, Nc, kfall, Tair, RH
    double precision :: tConst, tFall, tTotal, waterRemoved, finalRate, XstartFall
    double precision :: t, X, rate, removed, frac
    character(len=80) :: regime, period

    read(*,*,iostat=iostat_val) mode
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid drying mode input.'
        stop
    end if
    read(*,*,iostat=iostat_val) Ms
    read(*,*,iostat=iostat_val) A
    read(*,*,iostat=iostat_val) X0
    read(*,*,iostat=iostat_val) Xc
    read(*,*,iostat=iostat_val) Xe
    read(*,*,iostat=iostat_val) Xf
    read(*,*,iostat=iostat_val) Nc
    read(*,*,iostat=iostat_val) kfall
    read(*,*,iostat=iostat_val) Tair
    read(*,*,iostat=iostat_val) RH
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all drying inputs.'
        stop
    end if
    if (Ms <= 0.0d0 .or. A <= 0.0d0 .or. Nc <= 0.0d0 .or. kfall <= 0.0d0) then
        write(*,*) 'ERROR: Dry mass, area, constant rate, and falling coefficient must be positive.'
        stop
    end if
    if (X0 < 0.0d0 .or. Xc < 0.0d0 .or. Xe < 0.0d0 .or. Xf < 0.0d0) then
        write(*,*) 'ERROR: Moisture contents cannot be negative.'
        stop
    end if
    if (Xf >= X0) then
        write(*,*) 'ERROR: Final moisture must be below initial moisture.'
        stop
    end if
    if (Xf <= Xe) then
        write(*,*) 'ERROR: Final moisture must be greater than equilibrium moisture for finite falling-rate time.'
        stop
    end if

    tConst = 0.0d0
    if (mode == 1 .and. X0 > Xc .and. Xf < Xc) then
        tConst = Ms*(X0-Xc)/(A*Nc)
        XstartFall = Xc
        regime = 'Constant-rate followed by falling-rate drying'
    else if (mode == 1 .and. X0 > Xc .and. Xf >= Xc) then
        tConst = Ms*(X0-Xf)/(A*Nc)
        XstartFall = Xf
        regime = 'Constant-rate drying only to requested final moisture'
    else
        XstartFall = X0
        regime = 'Falling-rate drying only'
    end if

    if (Xf < XstartFall) then
        tFall = log((XstartFall-Xe)/(Xf-Xe))/kfall
    else
        tFall = 0.0d0
    end if
    tTotal = tConst + tFall
    waterRemoved = Ms*(X0-Xf)
    if (Xf >= Xc .and. mode == 1) then
        finalRate = Nc
    else
        finalRate = kfall*Ms/A*(Xf-Xe)
    end if

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   DRYING RATE CALCULATOR ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,I8)')       '  Drying Mode              = ', mode
    write(*,'(A,ES12.4,A)') '  Dry Solid Mass           = ', Ms, ' kg'
    write(*,'(A,ES12.4,A)') '  Drying Area              = ', A, ' m2'
    write(*,'(A,ES12.4)')   '  Initial Moisture X0      = ', X0
    write(*,'(A,ES12.4)')   '  Critical Moisture (X_c)  = ', Xc
    write(*,'(A,ES12.4)')   '  Equilibrium Moisture Xe  = ', Xe
    write(*,'(A,ES12.4)')   '  Final Moisture Xf        = ', Xf
    write(*,'(A,ES12.4,A)') '  Constant Drying Rate (N_c) = ', Nc, ' kg/m2.h'
    write(*,'(A,ES12.4,A)') '  Falling Rate Coeff       = ', kfall, ' 1/h'
    write(*,'(A,F12.3,A)')  '  Air Temperature          = ', Tair, ' deg-C'
    write(*,'(A,F12.3,A)')  '  Relative Humidity        = ', RH, ' percent'
    write(*,*)
    write(*,'(A)') '--- DRYING RESULTS ------------------------------------------'
    write(*,'(A,ES12.4,A)') '  Total Drying Time        = ', tTotal, ' h'
    write(*,'(A,ES12.4,A)') '  Constant Rate Time       = ', tConst, ' h'
    write(*,'(A,ES12.4,A)') '  Falling Rate Time        = ', tFall, ' h'
    write(*,'(A,ES12.4,A)') '  Water Removed            = ', waterRemoved, ' kg'
    write(*,'(A,ES12.4,A)') '  Final Drying Rate        = ', finalRate, ' kg/m2.h'
    write(*,'(A,A)')        '  Drying Regime            = ', trim(regime)
    write(*,*)
    write(*,'(A)') '--- DRYING CURVE PROFILE ------------------------------------'
    write(*,'(A)') '  time[h]      X[kg/kg]      rate[kg/m2.h] removed[kg]    period'
    write(*,'(A)') '  ----------------------------------------------------------------'
    n_points = 80
    do i=0,n_points-1
        frac = dble(i)/dble(n_points-1)
        t = frac*tTotal
        if (t <= tConst .and. tConst > 0.0d0) then
            X = X0 - A*Nc/Ms*t
            rate = Nc
            period = 'CONST'
        else
            if (tFall > 0.0d0) then
                X = Xe + (XstartFall-Xe)*exp(-kfall*(t-tConst))
                rate = kfall*Ms/A*(X-Xe)
            else
                X = Xf
                rate = finalRate
            end if
            period = 'FALL'
        end if
        if (X < Xf) X = Xf
        removed = Ms*(X0-X)
        write(*,'(F10.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,A)') t, X, rate, removed, trim(period)
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Constant-rate period: t = Ms*(X_initial-X_critical)/(A*Nc).'
    write(*,'(A)') '  Falling-rate model: dX/dt = -kfall*(X-Xe).'
    write(*,'(A)') '  Falling-rate time: ln((Xstart-Xe)/(Xf-Xe))/kfall.'

end program drying_rate