! -----------------------------------------------------------------------------
! Milk as a mixture of conventional components:
! mixture H2O (Aspen Plus) + FAT + PROTEIN + LACTOSE (Aspen Plus) + ASH
!   - Liquid molar enthalpy
! -----------------------------------------------------------------------------
C
      SUBROUTINE HL2U (T, P, Z, N, IDX, IRW, IIW, KCALC, KOP,
     &NDS, KDIAG, QMX, DQMX, KER) 
C
      IMPLICIT NONE
      DIMENSION Z(N), IDX(N), KOP(10), NBOPST(6), NAME(2)
      DIMENSION PHI(N), H(N), S(N), G(N), V(N), B(1)
      DIMENSION DPHI(N), DH(N), DS(N), DG(N), DV(N)
      INTEGER N, DMS_KCCIDC, IMW, IMF,IMP,IMC,IMA
      INTEGER IDX, I, II, IIMW, IIMF, IIMP, IIMC, IIMA
      INTEGER KDIAG, KER, IRW, IIW, KCALC, KOP, NDS
      INTEGER NBOPST, NAME, KV, KH, KS, KG, KBASE 
      INTEGER IOFF, LMWW, LMWF, LMWP, LMWC, LMWA
      INTEGER KPHI, DMS_IFCMNC
      REAL*8 B, Z, T, TN, TREF, P
      REAL*8 MWW, MWF, MWP, MWC, MWA, AVG_MW,
     &       MASSW, MASSF, MASSP, MASSC, MASSA, MASSTOT
      REAL*8 FRACW, FRACF, FRACP, FRACC, FRACA
      REAL*8 WC, FC, PC, CC, AC
      REAL*8 PHI, H, S, G, V, DPHI, DH, DS, DG, DV
      REAL*8 THEHLMXM, THEHLMX, HREF, HVAR
      REAL*8 CPMILK, CPN,  QMX, DQMX
      REAL*8 DELT1, DELT2, DELT3, DELT4, DELT5,
     +     INCPW, INCPF, INCPP, INCPC, INCPA,
     +     FA1,FA2,FA3,PA1,PA2,PA3,CA1,CA2,CA3,AA1,AA2,AA3
      EQUIVALENCE (B(1), IB(1))
      INTEGER IER
C     
#include "dms_ncomp.cmn"
#include "ppexec_user.cmn"
#include "dms_ppwork.cmn"
#include "dms_ipwork.cmn"
#include "dms_plex.cmn"
#include "dms_maxwrt.cmn"
C
      IER = 0
C     Retrieving components index:
      CALL COMPO(N, IDX, Z, IMW, IMF, IMP, IMC, IMA,
     *     IIMW, IIMF, IIMP, IIMC, IIMA,
     *     WC, FC, PC, CC, AC, 
     *     MWW, MWF, MWP, MWC, MWA, IER)
      IF(IER.NE.0) RETURN
C     If WC = 1 then do not use other components
      IF (WC.EQ.1) THEN
C       Reference temperature and enthalpy
 13     TREF = 25d0
        HREF = 0d0
C       Temperature in °C
        TN = T -273.15d0
C       Temperature differences
        DELT1 = TN - TREF
        DELT2 = ((TN**2) - (TREF**2))/2.D0
        DELT3 = ((TN**3) - (TREF**3))/3.D0
        DELT4 = ((TN**4) - (TREF**4))/4.D0
        DELT5 = ((TN**5) - (TREF**5))/5.D0
C       Compute water Cp
        INCPW = 4.23 *DELT1
     &       -2.62E-3 *DELT2
     &       +4.18E-5 *DELT3 -2.15E-7 *DELT4 +5.20E-10 *DELT5
C       Mixture enthalpy (J.kg^-1)
        THEHLMXM = (HREF + INCPW)*1000d0
!       Mixture enthalpy (J.kmol^-1) 
!       THEHLMX = (HREF + INCPW)*1000
        THEHLMX = THEHLMXM*MWW
      ELSE
!       Conversion to mass fraction:
!       4) Calculate mass of components in 1 mole
        MASSW = MWW*WC
        MASSF = MWF*FC
        MASSP = MWP*PC
        MASSC = MWC*CC
        MASSA = MWA*AC
!       MASSTOT = Total mass of 1 mole of milk mixture
        MASSTOT = MASSW+MASSF+MASSP+MASSC+MASSA
!       5) Calculate mass fractions
        FRACW = MASSW / MASSTOT
        FRACF = MASSF / MASSTOT
        FRACP = MASSP / MASSTOT
        FRACC = MASSC / MASSTOT
        FRACA = MASSA / MASSTOT
        IF (KOP(1).EQ.0) THEN
!         Calculate with Choi & Okos (1986) model
!         Reference temperature and enthalpy
          TREF = 25d0
          HREF = 0d0
!         Temperature in °C
          TN = T -273.15d0
!         Temperature differences
          DELT1 = TN - TREF
          DELT2 = (TN**2 - TREF**2)/2.D0
          DELT3 = (TN**3 - TREF**3)/3.D0
          DELT4 = (TN**4 - TREF**4)/4.D0
          DELT5 = (TN**5 - TREF**5)/5.D0
C         
C         Here the MW has an impact
C         
!         Coefficients for heat capacities of food components
          FA1 = 1.9842*MWF
          FA2 = 0.0014733*MWF
          FA3 = -0.0000048008*MWF
          PA1 = 2.0082*MWP
          PA2 = 0.0012089*MWP
          PA3 = -0.0000013129*MWP
          CA1 = 1.5488*MWC
          CA2 = 0.0019625*MWC
          CA3 = -0.0000059399*MWC
          AA1 = 1.0926*MWA
          AA2 = 0.0018896*MWA
          AA3 = -0.0000036817*MWA
!         Calculate integral of individual heat capacities
          INCPW = WC*(76.15060135*DELT1 -0.0471862879*DELT2
     &         +0.00075231463*DELT3 
     *         -3.8782287E-06*DELT4 +9.3701E-09*DELT5)
          INCPF = FC*(FA1*DELT1 +FA2*DELT2
     &         +FA3*DELT3)
          INCPP = PC*(PA1*DELT1 +PA2*DELT2
     &         +PA3*DELT3)
          INCPC = CC*(CA1*DELT1 +CA2*DELT2
     &         +CA3*DELT3)
          INCPA = AC*(AA1*DELT1 +AA2*DELT2
     &         +AA3*DELT3)
C         
C         Variable term of enthalpy
          HVAR = INCPW + INCPF + INCPP + INCPC + INCPA
!         Mixture enthalpy (J.kmol^-1) 
          THEHLMX = (HREF + HVAR)*1000d0
          write(user_nhstry,*) 'THEHLMX : ', THEHLMX
        ELSE IF (KOP(1).EQ.1) THEN
!         Calculate with Minim et al. (2002) model
!         FC = Fat Content ; WC = Water Content
!         Heat capacity of the mixture (kJ.kg^-1.K^-1)
          CPMILK = 1.4017 +0.0021*T -1.7430*FRACF +2.1816*FRACW
!         Molar heat capacity
          CPN = CPMILK*MASSTOT
!         Reference temperature (K) and enthalpy
          TREF = 298.15d0
          HREF = 0d0
!         Temperature differences
          DELT1 = T - TREF
          DELT2 = (T**2 - TREF**2)/2.D0
!         Integral of heat capacity from TREF to T
          HVAR = 1.4017*DELT1 +0.0021*DELT2 -1.7430*FRACF*DELT1
     &         +2.1816*FRACW*DELT1
!         Mixture enthalpy (J.kg^-1)
          THEHLMXM = (HREF + HVAR)*1000d0
!         Mixture enthalpy (J.kmol^-1)
          AVG_MW = MWW*WC + 1*FC + 1*PC + 1*CC + 1*AC
          THEHLMX = THEHLMXM*AVG_MW
        ELSE IF (KOP(1).EQ.2) THEN
!         For sensitivity studies: use Cp model from Choi & Okos (1986)
!         with +10% increase of the calculated Cp value
!         Reference temperature and enthalpy
          TREF = 25d0
          HREF = 0d0
C         Temperature in °C
          TN = T -273.15d0
C         Temperature differences
          DELT1 = TN - TREF
          DELT2 = (TN**2 - TREF**2)/2.D0
          DELT3 = (TN**3 - TREF**3)/3.D0
          DELT4 = (TN**4 - TREF**4)/4.D0
          DELT5 = (TN**5 - TREF**5)/5.D0
!         Coefficients for heat capacities of food components
          FA1 = 1.9842*MWF
          FA2 = 0.0014733*MWF
          FA3 = -0.0000048008*MWF
          PA1 = 2.0082*MWP
          PA2 = 0.0012089*MWP
          PA3 = -0.0000013129*MWP
          CA1 = 1.5488*MWC
          CA2 = 0.0019625*MWC
          CA3 = -0.0000059399*MWC
          AA1 = 1.0926*MWA
          AA2 = 0.0018896*MWA
          AA3 = -0.0000036817*MWA
C         Calculate integral of individual heat capacities
          INCPW = WC*(76.15060135*DELT1 -0.0471862879*DELT2
     &         +0.00075231463*DELT3 -3.8782287E-06*DELT4 
     *         +9.3701E-09*DELT5)
          INCPF = FC*(FA1*DELT1 +FA2*DELT2
     &         +FA3*DELT3)
          INCPP = PC*(PA1*DELT1 +PA2*DELT2
     &         +PA3*DELT3)
          INCPC = CC*(CA1*DELT1 +CA2*DELT2
     &         +CA3*DELT3)
          INCPA = AC*(AA1*DELT1 +AA2*DELT2
     &         +AA3*DELT3)
C         Variable term of enthalpy
          HVAR = INCPW + INCPF + INCPP + INCPC + INCPA
C         Mixture enthalpy (J.kmol^-1) 
          THEHLMX = ((HREF + HVAR)*1000d0)*1.1
        ELSE
C         Consider pure water
          GOTO 13
        ENDIF
      ENDIF
C     Final result
      QMX = THEHLMX
 999  CONTINUE
!     check for IER to write to control panel
      IF(IER.NE.0)THEN
        WRITE (MAXWRT_MAXBUF, 1000)
        CALL DMS_WRTTRM(1)
      ENDIF
C
      RETURN
 1000 FORMAT(' Error in HL2USR - check history file')
      END

