! -----------------------------------------------------------------------------
! Milk as a mixture of conventional components:
! mixture H2O (Aspen Plus) + FAT + PROTEIN + LACTOSE (Aspen Plus) + ASH
!   - Liquid surface tension
! -----------------------------------------------------------------------------
C
      SUBROUTINE SIG2U (T, P, Z, N, IDX, IRW, IIW, KCALC, KOP,
     &NDS, KDIAG, QMX, DQMX, KER)
C
        IMPLICIT NONE
        DIMENSION Z(N), IDX(N), KOP(10)
        INTEGER N, DMS_KCCIDC, ICW, IMF,IMP,IMC,IMA
        INTEGER IDX, IOFF
        INTEGER LMWW, LMWF, LMWP,LMWC, LMWA
        INTEGER IRW, IIW, KCALC, KOP, NDS, KDIAG, KER     
        INTEGER I, II, IIMW,IIMF, IIMP, IIMC, IIMA, J
        REAL*8 Z, WC, P, T, TN, THESIGMX, QMX, DQMX, B(1)
        REAL*8 SIGWATER, SIGDIP, TRI
        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
        !XMW(I) = DMS_IFCMNC('MW') + I ! Retrieve molecular weight
        SIGDIP(I,J) = DMS_IFCMNC('SIGDIP') + 11*(J -1) + I
        ! Retrieving components index:
        ICW = DMS_KCCIDC('H2O')
        !IMF = DMS_KCCIDC('MILKFAT')
        !IMP = DMS_KCCIDC('MILKPROT')
        !IMC = DMS_KCCIDC('MILKCARB')
        !IMA = DMS_KCCIDC('MILKASH')
C
        ! Looking for components index in the IDX vector
        DO 10 I=1,N
            II = IDX(I)
            IF (II.EQ.ICW) IIMW = I
            IF (II.EQ.IMF) IIMF = I
            IF (II.EQ.IMP) IIMP = I
            IF (II.EQ.IMC) IIMC = I
            IF (II.EQ.IMA) IIMA = I
10      CONTINUE
        ! Accessing molar (weight) fractions
        WC = Z(IIMW)
        ! If WC = 1 then do not use other components
        IF (WC.EQ.1) THEN
          ! Retrieve surface tension of water
          ! DO NOT call mixture surface tension monitor !!!
          ! --> It would call this subroutine !
          TRI = T/647.096
          SIGWATER = B(SIGDIP(1,IIMW))*(1 - TRI) **( B(SIGDIP(2,IIMW))
     &    + B(SIGDIP(3,IIMW))*TRI + B(SIGDIP(4,IIMW))*TRI**2
     &    + B(SIGDIP(5,IIMW))*TRI**3 ) 
          THESIGMX = SIGWATER
        ELSE
!1015    FORMAT('MassFracWater: ', E11.4)
!        WRITE (MAXWRT_MAXBUF, 1015) FRACW
!        CALL DMS_WRTTRM(1)
!1016    FORMAT('MassFracFat: ', E11.4)
!        WRITE (MAXWRT_MAXBUF, 1016) FRACF
!        CALL DMS_WRTTRM(1)
!1010    FORMAT('MUWATER: ', E11.4)
!        WRITE (MAXWRT_MAXBUF, 1010) MUWATER
!        CALL DMS_WRTTRM(1)
          ! Temperature in °C
          TN = T - 273.15
          IF (KOP(1).EQ.0) THEN
              ! Calculate mixture surface tension as a function
              ! of temperature only : Bertsch, 1983
              THESIGMX = 1e-3 * (1.8e-4 * TN**2 -0.163 * TN + 55.6)
          ELSE IF (KOP(1).EQ.2) THEN
		  ! For sensitivity studies: use surface tension model from Bertsch (1983)
		  ! with +10% increase of the calculated surface tension value
              THESIGMX = (1e-3 * (1.8e-4 * TN**2 -0.163 * TN + 55.6))
     &                   *1.1
          ENDIF
!1001    FORMAT('THESIGMX: ', E11.4)
!        WRITE (MAXWRT_MAXBUF, 1001) THESIGMX
!        CALL DMS_WRTTRM(1)
1101      FORMAT('SIGMAMX: ', E11.4)
          write (USER_NHSTRY, 1101) THESIGMX
!******************************************************
        ! Final property result
        ENDIF
        QMX = THESIGMX
C
      END