Freeze Drying Sublimation (Pikal Lyophilization)

program lyophilization_freeze_drying
    implicit none
    integer :: iostat_val, step, num_steps
    double precision :: V_fill_mL, Av_cm2, cs_pct, L0_mm, Pch_mTorr, Tshelf_C, Kv_Wm2K, Rp0_Torr
    double precision :: Pch_Pa, dH_sub_Jg, rho_ice_gcm3, Tcollapse_C
    double precision :: Ldry_mm, Tice_C, Tice_K, Psat_Pa, Psat_Torr, Rp_Torr, dmdt_gcm2h
    double precision :: dt_step_h, total_time_h, q_flux_Wm2, Ldry_step
    character(len=32) :: collapse_status

    ! Read inputs
    read(*,*,iostat=iostat_val) V_fill_mL      ! Vial Fill Volume [mL] (e.g. 3.0)
    read(*,*,iostat=iostat_val) Av_cm2         ! Vial Base Area [cm2] (e.g. 3.8)
    read(*,*,iostat=iostat_val) cs_pct         ! Solid Formulation Content [%] (e.g. 5.0)
    read(*,*,iostat=iostat_val) Pch_mTorr      ! Chamber Vacuum [mTorr] (e.g. 100.0)
    read(*,*,iostat=iostat_val) Tshelf_C       ! Shelf Temperature [deg C] (e.g. -10.0)
    read(*,*,iostat=iostat_val) Kv_Wm2K        ! Vial Heat Transfer Coeff [W/(m2.K)] (e.g. 25.0)
    read(*,*,iostat=iostat_val) Rp0_Torr       ! Base Dry Cake Resistance [Torr.hr.cm2/g] (e.g. 1.8)
    read(*,*,iostat=iostat_val) Tcollapse_C    ! Critical Collapse Temp [deg C] (e.g. -28.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for lyophilization calculation.'
        stop
    end if

    if (V_fill_mL <= 0.0d0 .or. Av_cm2 <= 0.0d0 .or. Pch_mTorr <= 0.0d0) then
        write(*,*) 'ERROR: Fill volume, area, and chamber pressure must be positive.'
        stop
    end if

    Pch_Pa = Pch_mTorr * 0.133322d0
    dH_sub_Jg = 2835.0d0 ! J/g sublimation latent heat
    rho_ice_gcm3 = 0.917d0
    L0_mm = (V_fill_mL / Av_cm2) * 10.0d0 ! mm cake height

    ! Numerical integration of sublimation front progression
    num_steps = 100
    Ldry_step = L0_mm / dble(num_steps)
    total_time_h = 0.0d0

    ! Estimate steady sublimation front temperature Tice [deg C]
    ! Balance: Kv * (Tshelf - Tice) = dmdt * dH_sub
    ! Sublimation front typically stays around -30 to -20 C during primary drying
    Tice_C = Tshelf_C - 18.0d0
    if (Tice_C > -5.0d0) Tice_C = -5.0d0
    if (Tice_C < -45.0d0) Tice_C = -45.0d0
    Tice_K = Tice_C + 273.15d0

    Psat_Pa = exp(28.89d0 - 6140.0d0 / Tice_K)
    Psat_Torr = Psat_Pa / 133.322d0

    do step = 1, num_steps
        Ldry_mm = (dble(step) - 0.5d0) * Ldry_step
        ! Dynamic cake resistance increases with dry layer thickness
        Rp_Torr = Rp0_Torr + 1.2d0 * (Ldry_mm / 10.0d0)
        dmdt_gcm2h = max(0.001d0, (Psat_Torr - (Pch_mTorr * 1.0d-3)) / Rp_Torr)
        
        ! Step time dt = (rho_ice * (1 - cs) * dL) / dmdt
        dt_step_h = (rho_ice_gcm3 * (1.0d0 - cs_pct/100.0d0) * (Ldry_step * 0.1d0)) / dmdt_gcm2h
        total_time_h = total_time_h + dt_step_h
    end do

    q_flux_Wm2 = Kv_Wm2K * (Tshelf_C - Tice_C)

    if (Tice_C < Tcollapse_C) then
        collapse_status = 'SAFE (No Product Collapse)'
    else
        collapse_status = 'WARNING: Micro-Collapse Risk!'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — LYOPHILIZATION SUBLIMATION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Vial Fill / Cake Depth    = ', V_fill_mL, ' mL / ', L0_mm, ' mm'
    write(*,'(A,F10.1,A,F10.1,A)') 'Chamber P / Shelf Temp    = ', Pch_mTorr, ' mTorr / ', Tshelf_C, ' C'
    write(*,'(A,F10.2,A,F10.1,A)') 'Vial Heat Transfer Kv     = ', Kv_Wm2K, ' W/(m2.K) (q = ', q_flux_Wm2, ' W/m2)'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'PRIMARY DRYING DURATION   = ', total_time_h, ' hours'
    write(*,'(A,F10.2,A)')  'Sublimation Interface T   = ', Tice_C, ' deg C (Collapse Limit = ', Tcollapse_C, ' C)'
    write(*,'(A,A)')        'Cake Structural Integrity = ', trim(collapse_status)
    write(*,'(A,F10.2,A)')  'Ice Saturation Vapor P    = ', Psat_Torr * 1000.0d0, ' mTorr'
    write(*,'(A,F10.3,A)')  'Average Sublimation Flux  = ', (V_fill_mL * rho_ice_gcm3) / (Av_cm2 * total_time_h), ' g/(cm2.h)'
    write(*,'(A)') '============================================================'

end program lyophilization_freeze_drying