C     User Activity Coefficient Subroutine
C     23/09/2014 Martial Madoumier
C     07/11/2014 Version with "true" molar fractions works with
C                MW <> 1 for dry components.
C
      SUBROUTINE GMU (T,       P,      Z,      N,      IDX,
     2                IRW,     IIW,    KCALC,  KOP,    NDS,
     3                KDIAG,   GAMMA,  DGAMMA, KER,    NSUB,
     4                NSUP,    IDXSUB, IDXSUP, GAMSC,  DGAMSC)
C
      IMPLICIT NONE
C
C     DECLARE VARIABLES USED IN DIMENSIONING
C
      INTEGER N, NSUB, NSUP
C
C     DECLARE ARGUMENTS
C
      INTEGER IDX(N),KOP(10),      IDXSUB(NSUB),
     +        IDXSUP(NSUP), IRW,   IIW,   KCALC, NDS,
     +        KDIAG, KER,   I
      REAL*8 Z(N),  GAMMA(N),     DGAMMA(N),
     +       GAMSC(NSUB,NSUP),    DGAMSC(NSUB,NSUP),
     +       T,     P
C
C     DECLARE LOCAL VARIABLES
C
      INTEGER DMS_IFCMNC,
     +        DMS_KCCIDC,M,
     +        IMW, IMF, IMP, IMC, IMA,
     +        IIMW, IIMF, IIMP, IIMC, IIMA

      REAL*8 B(1), AntA,AntB,AntC,
     +       WC,FC,PC,CC,AC,BPE,
     +       MWW, MWF, MWP, MWC, MWA,
     +       FPSATWATER,FPSATSOLID,THET,THETW,
     +       THEGAMMA(100),PANT,TSATWATER,AW,
     +       Rw, SUMMOLFRACDC,Xsol,HVAPWATER,TR,
     +       TSATGOAL,PSATGOAL,PSATOLD,PSATSOLIDOLD,
     +       PSat(100),PSATMIX,PSATWNEW,PSATDRY,PSATSOLID,
     +       MASSTOT,MASSW,MASSF,MASSP,MASSC,MASSA,
     +       FRACW,FRACF,FRACP,FRACC,FRACA,
     +       MOLTOT,MOLW,MOLF,MOLP,MOLC,MOLA
      INTEGER IER

      EQUIVALENCE (B(1), IB(1))
      DATA AntA,AntB,AntC /5.11564,1687.537,-42.98/
C
#include "dms_ncomp.cmn"
#include "ppexec_user.cmn"
#include "pputl_ppglob.cmn"
#include "dms_ppwork.cmn"
#include "dms_ipwork.cmn"
#include "dms_plex.cmn"
C     To display messages in control panel
#include "dms_maxwrt.cmn"
C
C      FUNCTION STATEMENTS
C
      FPSATWATER(THETW) = exp(73.649
     &     - 7258.2/(THETW  + 0)
     &     + 0*THETW -7.3037 *log(THETW)
     &     + 4.1653e-6 *THETW**2)
C     
      FPSATSOLID(THET) = exp(1182.2 
     &     - 84682/(THET  + 0)
     &     + 0.1564*THET -175.85 *log(THET)
     &     - 2.38e-5 *THET**2)
C     
C
C     BEGIN EXECUTABLE CODE
C
      IER = 0
C
!     Retrieving components index:
      CALL COMPO(N, IDX, Z, IMW, IMF, IMP, IMC, IMA,
     *     IMW, IMF, IIMP, IIMC, IIMA,
     *     WC, FC, PC, CC, AC,
     *     MWW, MWF, MWP, MWC, MWA, IER)
      IF(IER.NE.0) RETURN
C
C     COMPUTE THE NATURAL LOGARITHM OF ACTIVITY COEFFICIENTS
C
C     Calculate Gamma <> 1 only if N > 1
      if (WC.EQ.1) then
C       If water only is present
 13     THEGAMMA(1) = 1d0
        if (N.GT.1d0) then
C         If other components involved
          DO I = 2,N
            THEGAMMA(I) = 1d0
          ENDDO
        endif     
      else
C       Protection in case WC considered = 0 (problem with log(1 - Xsol)
        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
        FRACW = MASSW / MASSTOT
        FRACF = MASSF / MASSTOT
        FRACP = MASSP / MASSTOT
        FRACC = MASSC / MASSTOT
        FRACA = MASSA / MASSTOT
        if (FRACW.LT.0.2) then
          FRACW = 0.2
          MOLW = FRACW/MWW
          MOLF = FRACF/MWF
          MOLP = FRACP/MWP
          MOLC = FRACC/MWC
          MOLA = FRACA/MWA
          MOLTOT = MOLW+MOLF+MOLP+MOLC+MOLA
          WC = MOLW/MOLTOT
        endif
C       Saturation temperature of water
        PANT = P*1e-5
        TSATWATER = -(AntB/(log10(PANT) - AntA))-(AntC)
C
C       BPE from literature model
C
!       Compute variables for BPE model
        Rw = 8.3145 /(MWW *1e-3)
C        SUMMOLFRACDC = TmolFracF+TmolFracP+TmolFracC+TmolFracA
C        Xsol = SUMMOLFRACDC/(TmolFracW + SUMMOLFRACDC)
        SUMMOLFRACDC = FC+PC+CC+AC
        Xsol = SUMMOLFRACDC/(WC + SUMMOLFRACDC)
        TR = T/647.096
        if (TR.GE.1) then
          IER = 1
          write(user_nhstry,*) 'Error - TR >= 1'
          GOTO 999
        endif
        write(user_nhstry,*) 'TR : ', TR
        HVAPWATER = (56600000*(1 - TR)**(0.61204 -0.6257 *TR
     &       +0.3988 *TR**2 + 0*TR**3))/MWW
        write(user_nhstry,*) 'HVAPWATER : ', HVAPWATER
        CALL BPELMX(TSATWATER, Xsol, HVAPWATER,
     &              KOP(1), 1, BPE)
        write(user_nhstry,*) 'BPE : ', BPE
C         IF (KOP(1).EQ.0) THEN
C !         Winchester (2000)
C           BPE = (1 /(1/TSATWATER + Rw*log(1 - Xsol) /HVAPWATER))
C      &         - TSATWATER
C !         Attribute boundaries to BPE
C           if (BPE.GE.10) then
C             BPE = 10
C           elseif (BPE.LT.0) then
C             BPE = 0
C             goto 13
C           endif
C         ELSE IF (KOP(1).EQ.2) THEN
C !         For sensitivity studies: use BPE model from Winchester (2000)
C !         with +10% increase of the calculated BPE value
C           BPE = ((1 /(1/TSATWATER + Rw*log(1 - Xsol) /HVAPWATER))
C      &         - TSATWATER)*1.1
C !         Attribute boundaries to BPE
C           if (BPE.GE.10) then
C             BPE = 10
C           elseif (BPE.LT.0) then
C             BPE = 0
C             goto 13
C           endif
C         ELSE
C !         No difference with pure water
C           BPE = 0
C           goto 13
C         ENDIF
C
C       End of BPE models
C
!       Desired true water boiling temperature in milk
        TSATGOAL = TSATWATER + BPE
!       Calculate current vapour pressure (no gamma)
        PSATSOLIDOLD = FPSATSOLID(TSATWATER)
        PSATOLD = WC*FPSATWATER(TSATWATER)
     &       +FC*PSATSOLIDOLD
     &       +PC*PSATSOLIDOLD+CC*PSATSOLIDOLD
     &       +AC*PSATSOLIDOLD		
!       Calculate corresponding bubble pressure
        PSATSOLID = FPSATSOLID(TSATGOAL)
        PSATGOAL = WC*FPSATWATER(TSATGOAL)
     &       +FC*PSATSOLID
     &       +PC*PSATSOLID+CC*PSATSOLID
     &       +AC*PSATSOLID	
!       Calculate vapour pressure from Antoine model
        PSat(1) = FPSATWATER(T)
        DO M = 2,N
          PSat(M) = FPSATSOLID(T)
        ENDDO
!       Bubble pressure at stream temperature
        PSATMIX = WC*PSat(1)+FC*PSat(2)
     &       +PC*PSat(3)+CC*PSat(4)
     &       +AC*PSat(5)		
C       New partial pressure of water at saturation
        PSATWNEW = PSATGOAL - (PSATMIX - WC*PSat(1))
C       Saturation pressure of dry matter components
        PSATDRY = FC*PSATSOLID
     &       +PC*PSATSOLID+CC*PSATSOLID
     &       +AC*PSATSOLID
        if (WC.EQ.0) then
C         No water ERROR
          THEGAMMA(1) = 1
          DO I = 2,N
            THEGAMMA(I) = 1
          ENDDO
        else
C         Water + dry matter components
C          THEGAMMA(1) =  (P - (PSATDRY))/
C     &         (WC*FPSATWATER(TSATGOAL))

          ! Compute aw from Berry et al. (1980) formula
          AW = exp((BPE * HVAPWATER) / (-Rw * (TSATWATER**2)))
          write(user_nhstry,*) 'AW : ', AW
          THEGAMMA(1) = AW / WC
          write(user_nhstry,*) 'THEGAMMA(1) : ', THEGAMMA(1)
          DO I = 2,N
            THEGAMMA(I) = 1
          ENDDO
        endif
      endif
      
C     Final result
      DO I = 1, N
        GAMMA(I) = log(THEGAMMA(I))
      ENDDO
 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 GMU - check history file')
      END
