! -----------------------------------------------------------------------------
! Milk as a mixture of conventional components:
! mixture H2O (Aspen Plus) + FAT + PROTEIN + LACTOSE (Aspen Plus) + ASH
!   - Liquid dynamic viscosity
! -----------------------------------------------------------------------------
C
      SUBROUTINE MUL2U (T, P, Z, N, IDX, IRW, IIW, KCALC, KOP,
     &     NDS, KDIAG, QMX, DQMX, KER)
C     
      IMPLICIT NONE
      INTEGER N, DMS_KCCIDC, ICW, IMF,IMP,IMC,IMA
      INTEGER IDX(N), IOFF
      INTEGER LMWW, LMWF, LMWP,LMWC, LMWA
      INTEGER IRW, IIW, KCALC, KOP(10), NDS, KDIAG, KER     
      INTEGER I, J,II, IIMW,IIMF, IIMP, IIMC, IIMA,
     &     IMSF, IIMSF
      REAL*8 Z(n), P, FC, PC, AC, CC, SFC
      REAL*8 X, WC, T, TN, MULMX, 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 ACARB, AFAT, APROT, MUWATER, MULDIP,RESALC
      EQUIVALENCE (B(1), IB(1))
      INTEGER DMS_IFCMNC
      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"
!     Statement functions
      MULDIP(I,J) = DMS_IFCMNC('MULDIP') + 11*(J -1) + I
!     Retrieving components index:
      CALL COMPO(N, IDX, Z, ICW, 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
!     Retrieve viscosity of water
      MUWATER = exp(B(MULDIP(1,IIMW)) + B(MULDIP(2,IIMW)) / T
     &     + B(MULDIP(3,IIMW))*log(T) +B(MULDIP(4,IIMW)) * T**2
     &     + B(MULDIP(5,IIMW)) / T**2)
!     If WC = 1 then do not use other components
      IF (WC.EQ.1) THEN
 13     MULMX = MUWATER
      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) Convert molar fraction to mass fraction
        FRACW = MASSW / MASSTOT
        FRACF = MASSF / MASSTOT
        FRACP = MASSP / MASSTOT
        FRACC = MASSC / MASSTOT
!       Temperature in °C
        TN = T - 273.15
!       Viscosity of the mixture (Pa.s)
        IF (KOP(1).EQ.0) THEN
!         Morison et al. (2013)
          if (TN > 100) then    ! Protection in case of out-of-range conditions
            TN = 100
          endif
          ACARB = 3.35 -2.38e-2 *TN + 1.25e-4 * TN**2
          AFAT = 3.46 -0.025 *TN +1.6e-4 * TN**2
          APROT = 15.367 -0.178 *TN +0.0017 * TN**2
!         Protection in case of out-of-range conditions
          IF (FRACW < 0.5) THEN
            FRACW = 0.5
          ENDIF
          MULMX = MUWATER*exp(((ACARB*FRACC + 
     &         AFAT*FRACF + APROT*FRACP)/FRACW))	
        ELSE IF (KOP(1).EQ.1) THEN
!         Alcantara et al. (2012) (whole milk)
          RESALC = 174300e-9 * exp((-8.9589 * FRACW)
     &         + (25.8339e3 /(8.3145 *T)))
          MULMX = RESALC*1000
        ELSE IF (KOP(1).EQ.2) THEN
          if (TN > 100) then    ! Protection in case of out-of-range conditions
            TN = 100
          endif
!         For sensitivity studies: use viscosity model from Morison et al. (2013)
!         with +10% increase of the calculated viscosity value
          ACARB = 3.35 -2.38e-2 *TN + 1.25e-4 * TN**2
          AFAT = 3.46 -0.025 *TN +1.6e-4 * TN**2
          APROT = 15.367 -0.178 *TN +0.0017 * TN**2
!         Protection in case of out-of-range conditions
          IF (FRACW < 0.5) THEN
            FRACW = 0.5
          ENDIF
          MULMX = (MUWATER*exp(((ACARB*FRACC + 
     &         AFAT*FRACF + APROT*FRACP)/FRACW)))*1.1
        ELSE
!         Consider pure water
          GOTO 13
        ENDIF
      ENDIF
!     Final property result
      QMX = MULMX
C
      RETURN
      END
