! -----------------------------------------------------------------------------
! Milk as a mixture of conventional components:
! mixture H2O (Aspen Plus) + FAT + PROTEIN + LACTOSE (Aspen Plus) + ASH
!   - Liquid molar volume
! -----------------------------------------------------------------------------
C
      SUBROUTINE VL2U (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), X(N)
      DIMENSION V(N)
      INTEGER N, DMS_KCCIDC, IMW, IMF,IMP,IMC,IMA
      INTEGER IDX
      INTEGER IRW, IIW, KCALC, KOP, NDS, KDIAG, KER     
      INTEGER I, II, IIMW,IIMF, IIMP, IIMC, IIMA, J
      INTEGER NBOPST, NAME, KV, KH, KS, KG, KBASE
      INTEGER KPHI
      REAL*8 Z, P, X, T, TN, RHO, RHON, QMX, DQMX, B(1)
      REAL*8 MWW, MWF, MWP, MWC, MWA
      REAL*8 MASSW, MASSF, MASSP, MASSC, MASSA, MASSTOT
      REAL*8 FRACW, FRACF, FRACP, FRACC, FRACA
      REAL*8 RHOW,RHOF,RHOP,RHOC,RHOA,RHOTOT
      REAL*8 TAUW, RHOWATER, THEVLMX, VMX, DVMX
      REAL*8 PHI, H, S, G, V, DPHI, DH, DS, DG, DV
      INTEGER IER
      REAL*8 WC, FC, PC, CC, AC
!     REAL*8 DNLDIP
      EQUIVALENCE (B(1), IB(1))
      INTEGER DMS_IFCMNC
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
! Statement functions
!     DNLDIP(I,J) = DMS_IFCMNC('DNLDIP') + 11*(J -1) + I
!     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
!     If WC = 1 then do not use other components
      IF (WC.EQ.1) THEN
!       Retrieve molar volume of water
!       TAUW = 1 - (T/647.096)
!       RHOWATER = B(DNLDIP(1,IMW)) + B(DNLDIP(2,IMW))*(TAUW**0.35)
!       &    + B(DNLDIP(3,IMW))*(TAUW**(2/3)) + B(DNLDIP(4,IMW))*TAUW
!     &    + B(DNLDIP(5,IMW))*(TAUW**(4/3))
          !THEVLMX =  1/RHOWATER
 13     CALL PPUTL_GOPSET (NBOPST, NAME)
        KV = 1
        KBASE = 0
        CALL PPMON_LTHRM (T, P, 1, IDX, NBOPST, KDIAG, KBASE,
     &       KPHI, KH, KS, KG, KV, PHI, H, S, G, V,
     &       DPHI, DH, DS, DG, DV, KER)
        THEVLMX = V(1)
      ELSE
        FC = Z(IIMF)
!       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) Convert molar fraction to mass fraction
        FRACW = MASSW / MASSTOT
        FRACF = MASSF / MASSTOT
        FRACP = MASSP / MASSTOT
        FRACC = MASSC / MASSTOT
        FRACA = MASSA / MASSTOT
!       Temperature in °C
        TN = T -273.15d0
!       Density of the mixture (kg/m3)
        IF (KOP(1).EQ.0) THEN
!         Choi & Okos (1986)
          RHOW = 997.18 + 3.1439E-3 *TN -3.7574e-3 *TN**2
          RHOF = 925.59 -4.1757E-1 *TN
          RHOP = 1329.9 -5.184E-1 *TN
          RHOC = 1599.1 -3.1046E-1 *TN
          RHOA = 2423.8 -2.8063E-1 *TN
          RHOTOT = 1/(FRACW/RHOW + FRACF/RHOF + FRACP/RHOP
     &         + FRACC/RHOC + FRACA/RHOA)
          RHO = RHOTOT
        ELSE IF (KOP(1).EQ.1) THEN
!         Minim et al. (2002) model - rho (kg/m3)
!         FC = Fat Content ; WC = Water Content
          RHO = 1185.64-0.341*TN+FRACF*(-58.107)+FRACW*(-58.239)
        ELSE IF (KOP(1).EQ.2) THEN
!         For sensitivity studies: use density model from Choi & Okos (1986)
!         with +10% increase of the calculated density value
          RHOW = 997.18 + 3.1439E-3 *TN -3.7574e-3 *TN**2
          RHOF = 925.59 -4.1757E-1 *TN
          RHOP = 1329.9 -5.184E-1 *TN
          RHOC = 1599.1 -3.1046E-1 *TN
          RHOA = 2423.8 -2.8063E-1 *TN
          RHOTOT = 1/(FRACW/RHOW + FRACF/RHOF + FRACP/RHOP
     &         + FRACC/RHOC + FRACA/RHOA)
          RHO = RHOTOT*1.1
        ELSE
!         Consider pure water
          GOTO 13
        ENDIF
!       Molar density (reasoning for 1 mole)
        RHON = (1*RHO)/MASSTOT
!       RHON = 1000
!       Conversion to liquid molar volume
        THEVLMX = 1/RHON
      ENDIF
      QMX = THEVLMX
      RETURN
      END
