Ask Question

Name:
Title:
Your Question:

Answer Question

Name:
Your Answer:
User Submitted Source Code!


Description:
  Solid liquid equilibria using UNIFAC models ([email protected])
Language: FORTRAN
Code:

      PROGRAM UNIFAC
c
      implicit double precision (A-H,O-Z)
      character*80 ntext,NTEXTE,nom
      character hefic*60,fichs*60,RFIC*60,RDON*60,GFIC*60,GDON*60
      character param*60,rep4*10
      dimension R(20),KG(20),NTEXT(40),NUM(20),IPAR(8),DEV(300,2)
      dimension NUME(20),NTEXTE(40),KGE(20),ECART(300,2),nom(40)
      dimension JENS(20),RS(20),QS(20),F(10),Vw(20)
      dimension RSE(20),QSE(20),ER(300,2),delta(20),Vm(20)
      dimension X(10,10),XB(10),XS(10),XM(10),XE(10),XK(10),XX(10)
     dimension XR(10),Vet(20),Cet(20)!,GMC1(300,2)
          
c
      common/points/ NG,NGR,NOBS,NPTS
      common/param/ A(20,20),B(20,20),C(20,20),rl,gs,ghe,Q(20)
      common/ABC/ T(300),NM(300,2),XXX(300,2),NNY(20,20),GMAX,NME(300,2)
     s          ,XXV(300,2)
      common/HE/ TE(300),XXXE(300,2),HEEXP(300,2),HECAL(300,2),HMAX
     s          ,NNYE(20,20),NU(300,2),NUE(300,2)
      common/CGM/ GME(300,2),GMC(300,2),GMR(300,2),IP(300),rep2,rep0
c     common/CGM/ GME(300,2),GMC(300,2),IP(300),rep2,rep0
c
       NP   = 6
       NRES = 7
c
      call system ("cls")
c
      write(*,*)
      write(*,*)' Voulez-vous utiliser le modèle de : '
      write(*,*)'                  -  Rasmussen              ( 1 )  '
      write(*,*)'                  -  Gmehling               ( 2 )  '
      write(*,*)
      read *,rep2
c
      write(*,'(A\)')' Nom de fichier de solubilitè         :   '
      read (*,212)fichs
      write(*,*)
      write(*,'(A\)')' Nom de fichier des enthalpies d excès :   '
      read (*,212)hefic
      write(*,*)
      write(*,'(A\)')' Matrice des paramètres                :   '
      read (*,212)param

 212  format(8A)
      if(rep2.EQ.1.or.rep2.EQ.3) then
           RFIC  = fichs
           RDON  = hefic
c          PARAM = 'param_r'
c
           open(UNIT= 5,FILE=RFIC,STATUS ='OLD')
           open(UNIT= 3,FILE=RDON,STATUS ='OLD')
c
      elseif(rep2.EQ.2) then
           GFIC  = fichs
           GDON  = hefic
c          PARAM = 'param_g'
c
           open(UNIT= 1,FILE=GFIC,STATUS ='OLD')
           open(UNIT= 2,FILE=GDON,STATUS ='OLD')
c
      elseif(rep2.EQ.3) then
           RFIC  = fichs
           RDON  = hefic
c          PARAM = 'param_s'
c
           open(UNIT= 5,FILE=RFIC,STATUS ='OLD')
           open(UNIT= 3,FILE=RDON,STATUS ='OLD')

      endif
c
      open(UNIT= 55,FILE=PARAM,STATUS ='OLD')
      open(UNIT= NP,FILE='sol.r')
      open(UNIT= NRES,FILE='he.r')
      open(UNIT= 33,FILE='enthalpi.vap')
      open(UNIT= 4,FILE='volume.mol')
      open(UNIT= 22,FILE='d_volume.bis')
      open(UNIT= 44,FILE='c_volume.txt')
c
      if(rep2.EQ.1.or.rep2.EQ.3) then
           NC = 5
           ND = 3
      elseif(rep2.EQ.2) then
           NC = 1
           ND = 2
      endif
c
      READ(NC,'(A)')
      READ(ND,'(A)')
      READ(NC,*)NCOMP,NG,NOBS,LAAF,NL
      READ(ND,*)NCONS,NGR,NPTS,NLIN
      write(*,*)
      write(*,*)' Voulez-vous :'
      write(*,*)'             - Une estimation des param�tres : (1) '
      write(*,*)'             - Un calcul simple              : (2) '
      read(*,*)NOIT
      write(*,*)

      IF(NOIT.EQ.2) GOTO 1301
      WRITE(NRES,15)
      WRITE(NP,15)
   15 FORMAT(1H ,'                  PARAMETER ESTIMATION'//)
      GOTO 1302
 1301 WRITE(NP,1303)
      WRITE(NRES,1303)
 1303 FORMAT(/1H ,20X,'CALCUL DES COEFFICIENTS D ACTIVITE'//)
 1302 CONTINUE
      DO 14 I = 1,NL
         READ(NC,13)NTEXT(I)
         WRITE(NP,13)NTEXT(I)
c        WRITE(*,13)NTEXT(I)
   13    FORMAT(25A80)
   14 CONTINUE
      write(NP,*)
      READ(NC,'(A)')
      READ(NC,*)(NUM(I),I=1,NCOMP)
c
      read *
      DO 145 I = 1,NLIN
         READ(ND,13)NTEXTE(I)
         WRITE(NRES,13)NTEXTE(I)
c        WRITE(*,13)NTEXTE(I)
  145 CONTINUE
      write(NRES,*)
      READ(ND,'(A)')
      READ(ND,*)(NUME(I),I=1,NCONS)
c
c  ---------------------------------------------------------------------
c  |     NUM   gives the numbers attached to                           |
c  |           the different components                                |
c  ---------------------------------------------------------------------
c
      READ(NC,'(A)')
      READ(ND,'(A)')
      DO 1304 I = 1,NCOMP
 1304  READ(NC,*)(NNY(I,K),K = 1,NG)
c
      DO 1314 I = 1,NCONS
 1314 READ(ND,*)(NNYE(I,K),K = 1,NGR)

c
c  ---------------------------------------------------------------------
c  |      NNY(I,K)  is the matrix giving the number of groups of       |
c  |                kind K in molecule I                               |
c  ---------------------------------------------------------------------
c
      READ(NC,'(A)')
      READ(NC,*)(KG(K),K = 1,NG)
c
      READ(ND,'(A)')
      READ(ND,*)(KGE(K),K = 1,NGR)
c
c  -----------------------------------------------------------------------
c  |       KG(K)   is the number attached to group K                     |
c  -----------------------------------------------------------------------
c
      READ(NC,'(A)')
      READ(ND,'(A)')
      DO 100 N = 1,NOBS
         READ(NC,*)IP(N),T(N),(NU(N,I),XXX(N,I),GME(N,I),I = 1,2)
  100 CONTINUE
c
      read(22,*)
      read(22,*)
      read(44,*)
      do i =1,10
        read(22,*)nom(i),Vw(i),Vm(i),delta(i)
        read(44,*)nom(i),Vm(i),Vet(i),Cet(i)
      enddo
c
       do N = 1,NOBS

        XXV(N,1) = XXX(N,1) * Vm(NU(N,1)) / (  XXX(N,1) * Vm(NU(N,1))
     s                      + ( XXX(N,2) ) * Vm( NU(N,2) ) )
        XXV(N,2) = XXX(N,2) * Vm(NU(N,2)) / (  XXX(N,1) * Vm(NU(N,1))
     s                      + ( XXX(N,2) ) * Vm( NU(N,2) ) )
       enddo
c
      DO 114 N = 1,NPTS
         READ(ND,*)TE(N),(NUE(N,I),XXXE(N,I),HEEXP(N,I),I = 1,2)
  114 CONTINUE
      READ(ND,'(A)')
      read(ND,*)HMAX
c
      IF(LAAF-1)820,820,821
  821 CONTINUE
      DO 822 N = 1,NOBS
        DO 822 I = 1,2
           GGG      = GME(N,I)
           GME(N,I) = DEXP(GGG)
  822 CONTINUE
  820 CONTINUE
c      write(*,*) ' NG =',NG
c      pause
      
      READ(55,'(A)')
      DO 705 I = 1,NG
         READ(55,*)R(I),Q(I)
  705 CONTINUE
c
      READ(55,'(A)')
      DO 706 I = 1,NG
         READ(55,*)(A(I,J),J=1,NG)
         READ(55,*)(B(I,J),J=1,NG)
         READ(55,*)(C(I,J),J=1,NG)
         READ(55,*)
  706 CONTINUE
  726 continue
c
      DO 8 I  = 1,NCOMP
          RS(I) = 0.D0
          QS(I) = 0.D0
      DO 8 J = 1,NG
            RS(I) = RS(I) + NNY(I,J) * R(J)
            QS(I) = QS(I) + NNY(I,J) * Q(J)
    8 CONTINUE
      WRITE(NP,5)
    5 FORMAT(1H ,'    COMPONENT /   GROUPS',/)
c
      DO 181 I  = 1,NCONS
          RSE(I) = 0.D0
          QSE(I) = 0.D0
        DO 181 J = 1,NGR
            RSE(I) = RSE(I) + NNYE(I,J) * R(J)
            QSE(I) = QSE(I) + NNYE(I,J) * Q(J)
  181 CONTINUE
      WRITE(NRES,5)
c
      WRITE(NP,12)(KG(I),I = 1,NG)
      WRITE(NRES,12)(KGE(I),I = 1,NGR)
   12 FORMAT(17X,20I4)
      DO 10 I = 1,NCOMP
   10 WRITE(NP,6)NUM(I),(NNY(I,J),J = 1,NG)
c
      DO 110 I = 1,NCONS
  110 WRITE(NRES,6)NUME(I),(NNYE(I,J),J = 1,NGR)
c
    6 FORMAT(10X,I3,4X,20I4)
      WRITE(NP,3)
      WRITE(NRES,3)
    3 FORMAT(//1H ,'  GROUP NO    GROUP R      GROUP Q',/)
      DO 11 I = 1,NG
   11 WRITE(NP,4)KG(I),R(I),Q(I)
c
      DO 211 I = 1,NGR
  211 WRITE(NRES,4)KGE(I),R(I),Q(I)
c
    4 FORMAT(4X,I3,5X,2F12.4)
      WRITE(NP,1400)
      WRITE(NRES,1400)
 1400 FORMAT(//1H ,'  GROUP INTERACTION PARAMETERS',/)
      write (NP,'(A\)') ' Param�tres  A(i,j) :'
      DO 1401 I = 1,NG
 1401 WRITE(NP,1402)(A(I,J),J=1,NG)
      write (NP,'(A\)') ' Param�tres  B(i,j) :'
      do 1501 I = 1,NG
 1501 write(NP,1402)(B(I,J),J=1,NG)
      write (NP,'(A\)') ' Param�tres  C(i,j) :'
      do 1502 I = 1,NG
 1502 write(NP,804)(C(I,J),J=1,NG)
 1402 format(20F15.5)
c
      write (NRES,'(A\)') ' Param�tres  A(i,j) :'
      DO 2401 I = 1,NGR
 2401 WRITE(NRES,1402)(A(I,J),J=1,NGR)
      write (NRES,'(A\)') ' Param�tres  B(i,j) :'
      do 2501 I = 1,NGR
 2501 write(NRES,1402)(B(I,J),J=1,NGR)
      write (NRES,'(A\)') ' Param�tres  C(i,j) :'
      do 2502 I = 1,NGR
 2502 write(NRES,804)(C(I,J),J=1,NGR)
c
      DO 102 N = 1,NOBS
      DO 102 I = 1,2
      DO 104 J = 1,NCOMP
        IF(NU(N,I)-NUM(J))104,103,104
  103 NM(N,I) = J
      GOTO 102
  104 CONTINUE
  102 CONTINUE
c
      DO 302 N = 1,NPTS
      DO 302 I = 1,2
      DO 304 J = 1,NCONS
        IF(NUE(N,I)-NUME(J))304,303,304
  303 NME(N,I) = J
       GOTO 302
  304 CONTINUE
  302 CONTINUE
c
      rl = 2.D0/3.D0
      print *,' rl = 2/3'
      write(*,*)
      write(*,'(A\)')' Modification de rl  ( o \ n ) ---> '
      read (*,'(A\)')rep4
      write(*,*)
         if(rep4.EQ.'o') then
             write(*,'(A\)')' Nouvelle valeur de rl --------> '
             read *,rl
         endif
      CALL PFAC3(RS,QS,delta,Vw,Vm)
c
c 138 continue
c
      IF( NOIT - 1 ) 830,830,831
  831 CONTINUE
c
         CALL PFAC4
         CALL HE
c
      DO 832 NR = 1,NOBS
        DO 832 I = 1,2
          GMR(NR,I) = GMC(NR,I) + GMR(NR,I)
c          GMR(NR,I) =  GMR(NR,I)  
          GMR(NR,I) = DEXP( GMR(NR,I) )
  832 CONTINUE
      GOTO 833
  830 CONTINUE
      KRIT = 1
      write(*,*)
      write(*,'(A\)')' Nombre de param�tres � estimer  NPAR =  '
      READ(*,*)NPAR
      write(*,*)
      write(*,*)'  Estimation du param�tre A -----> ( 0 )   '
      write(*,*)'  Estimation du param�tre B -----> ( 1 )   '
      write(*,'(A\)')'   Estimation de A et B ----------> ( 2 )   '
      READ(*,*)rep0
      write(*,*)
      write(*,'(A\)')' Nombre de param�tres identiques    IDEN =  '
      READ(*,*)IDEN
      write(*,*)
c      write(*,'(A\)')' Facteur de pond�ration (solubilit�) gs :  '
c      READ(*,*)gs
      gs = 1.D0
c      write(*,'(A\)')' Facteur de pond�ration (hE)         ghe:  '
c      READ(*,*)ghe
      ghe = 1.D0

c
      write(*,*)
      write(*,*)' ------------------------------------------------ '
      write(*,*)' |           CH3 : 1          ACH   :  5        | '
      write(*,*)' |           CH2 : 2          AC    :  6        | '
      write(*,*)' |           CH  : 3                            | '
      write(*,*)' |           CH0 : 4                            | '
      write(*,*)' ------------------------------------------------ '
      write(*,*)
      write(*,*)' Vous voulez estimer les param�tres de quels groupes ?'
      write(*,*)
c
            do I = 1,2
               write(*,*)' Groupe ',I,' :     '
               READ(*,*)IPAR(I)
            enddo
c
      IF(IDEN)950,950,951
  951 IDEN = 2 * IDEN
      write(*,*)
      write(*,*)' Quels sont les param�tres identiques ? '
      write(*,*)
          READ(*,*)(JENS(J),J = 1,IDEN)
  950 CONTINUE
c
c  ---------------------------------------------------------------------
c  |     JENS is the vector indicating the identical                   |
c  |          pairs of parameters                                      |
c  ---------------------------------------------------------------------
c
      NN  = NPAR + 1
      N   = NPAR
      write(*,*)
      write(*,*)' Valeurs initiales des param�tres :'
      write(*,*)
         do i = 1,NPAR
               read(*,*)X(1,i)
         enddo
      SA = 1.D-6
c

      DO 201 J = 2,NN
           DO 201 I = 1,N
                IF(J-I-1)202,203,202
  203                X(J,I) = 1.1D0 * X(1,I)
                     GOTO 201
  202                X(J,I) = X(1,I)
  201 CONTINUE
      WRITE(NP,300)(IPAR(I),I=1,N)
  300 FORMAT(/1H ,'   INITIAL PARAMETERS',8I3/)
      DO 204 J = 1,NN
          WRITE(NP,400)(X(J,I),I=1,N)
  204 continue
      DO 1 J = 1,NN
         DO 21 I = 1,N
   21       XX(I) = X(J,I)
      CALL FMIN(NPAR,IPAR,XX,FF,KRIT,JENS,IDEN)
    1 F(J)  = FF
      NF    = NN
c
c  ---------------------------------------------------------------------
c  |     NF is the number of calculations of F                         |
c  ---------------------------------------------------------------------
c
      ALFA  = 1.0D0
      BETA  = 0.5D0
      GAMMA = 2.0D0
      ITER  = 0
      JPR   = 0
  400 FORMAT(24F10.3)
c
c ----------------------------------------------------------------------
c |      ESTIMATION OF LOWEST VALUE OF F=FB                            |
c ----------------------------------------------------------------------
c
   25 FB = F(1)
      DO 98 I = 1,N
   98      XB(I) = X(1,I)
      DO 31 J = 2,NN
         IF(FB-F(J))31,31,108
  108    FB = F(J)
      DO 41 I = 1,N
   41      XB(I) = X(J,I)
   31 CONTINUE
c
c  ---------------------------------------------------------------------
c  |     ESTIMATION OF THE HIGHEST VALUE OF F=FS                       |
c  ---------------------------------------------------------------------
c
      FS = F(1)
      DO 51 I = 1,N
            XS(I) = X(1,I)
   51 continue
      JS = 1
      DO 61 J = 2,NN
       IF(FS-F(J))111,61,61
  111  FS = F(J)
       JS = J
       DO 71 I = 1,N
   71        XS(I) = X(J,I)
   61 CONTINUE
c
c  ---------------------------------------------------------------------
c  |     CALCULATION OF THE CENTROID XM(I) OF POINTS                   |
c  |     EXCLUDING XS(I)                                               |
c  ---------------------------------------------------------------------
c
      DO 81 I = 1,N
   81       XM(I) = - XS(I)
      DO  9 J = 1,NN
        DO 122 I = 1,N
  122          XM(I) = XM(I) + X(J,I)
    9 CONTINUE
      DO 121 I = 1,N
  121     XM(I) = XM(I) / FLOAT(N)
c
c  ---------------------------------------------------------------------
c  |     REFLECTION                                                    |
c  ---------------------------------------------------------------------
c
      DO 131 I = 1,N
  131       XR(I) = XM(I) + ALFA * ( XM(I) - XS(I) )
      CALL FMIN(NPAR,IPAR,XR,FR,KRIT,JENS,IDEN)
      NF = NF + 1
c
c  ---------------------------------------------------------------------
c  |     EXPANSION                                                     |
c  ---------------------------------------------------------------------
c
      IF(FR-FB)141,151,151
  141 DO 161 I = 1,N
  161     XE(I) = XM(I) + GAMMA * ( XR(I) - XM(I) )
      CALL FMIN(NPAR,IPAR,XE,FE,KRIT,JENS,IDEN)
      NF = NF +1
      IF(FE-FB)17,18,18
   17 DO 19 I = 1,N
           X(JS,I)   = XE(I)
   19      XS(I)     = XE(I)
      F(JS) = FE
c
c  ---------------------------------------------------------------------
c  |     CALCULATION OF HALTING CRITERION                              |
c  ---------------------------------------------------------------------
c
   27 FM   = 0.D0
      DO 20 J = 1,NN
   20 FM   = FM + F(J)
      FM   = FM / FLOAT(NN)
      FRMS = 0.D0
      DO 22 J = 1,NN
   22 FRMS = ( F(J) - FM )**2 + FRMS
      RMS  =  SQRT( FRMS / FLOAT(N) )
      ITER = ITER + 1
      JPR  = JPR  + 1
      IF(ITER-500)500,500,501
  501 continue
      write(*,*)
      write(*,*)' On d�passe 500 it�rations '
      write(*,*)' Arr�t de l''estimation   '
      write(*,*)
      goto 23
  500 CONTINUE
      IF(JPR-1)902,902,903
  903 CONTINUE
      IF(JPR-6)901,904,904
  904 JPR = 1
  902 CONTINUE
c      WRITE(NP,107)ITER,NF
      WRITE(*,118)
  118 format(/,'--------------------------------------------------------
     s----------------')
      WRITE(*,107)ITER,NF
  107 FORMAT(/,1H ,'   ITERATION',I4,'   NUMBER OF CALLS FOR THE SUBROUT
     sINE',I5)
      WRITE(*,109)
  109 FORMAT(/,'   PARAMETRES')
c
      WRITE(*,444)(X(JS,I),I=1,N)
  444 format(5X,20F14.5)
c
      WRITE(*,106)F(JS),RMS
  106 FORMAT(/,1H ,'  FMIN =',F14.5,'             SD =',E14.5)
  901 CONTINUE
      IF(RMS-SA)23,23,25
c
c  ---------------------------------------------------------------------
c  |     NEW SIMPLEX                                                   |
c  |     FE GREATER THAN FB                                            |
c  ---------------------------------------------------------------------
c
   18 DO 26 I = 1,N
           X(JS,I)   =   XR(I)
   26      XS(I)     =   XR(I)
      F(JS) = FR
      FS    = FR
      GO TO 27
c
c  ---------------------------------------------------------------------
c  |     NEW SIMPLEX                                                   |
c  |     FR GREATER THAN FB                                            |
c  ---------------------------------------------------------------------
c
  151 DO 30 J = 1,NN
        IF(J-JS)28,30,28
   28   IF(FR-F(J))18,18,30
   30 CONTINUE
        IF(FR-FS)91,91,32
   91 DO 33 I = 1,N
         X(JS,I)   = XR(I)
   33    XS(I)     = XR(I)
      F(JS) = FR
      FS    = FR
   32 DO 34 I = 1,N
   34    XK(I) = XM(I) + BETA * ( XS(I) - XM(I) )
      CALL FMIN(NPAR,IPAR,XK,FK,KRIT,JENS,IDEN)
      NF = NF + 1
      IF(FK-FS)35,35,36
   35 DO 37 I = 1,N
         X(JS,I)   = XK(I)
   37    XS(I)     = XK(I)
      F(JS) = FK
      FS    = FK
      GO TO 27
   36 DO 38 J = 1,NN
         DO 39 I = 1,N
   39         X(J,I) = ( X(J,I) + XB(I) ) / 2.D0
   38 CONTINUE
      GO TO 27
   23 WRITE(NP,905)
      WRITE(NRES,905)
  905 FORMAT(//'     FINAL PARAMETERS : ')
      WRITE(NP,444)(X(JS,I),I=1,N)
      WRITE(NP,*)
      WRITE(NP,106)F(JS),RMS
c
      WRITE(NRES,444)(X(JS,I),I=1,N)
      WRITE(NRES,*)
      WRITE(NRES,106)F(JS),RMS
  833 CONTINUE
      np1 = 0
      np2 = 0
      DO 906 N = 1,NOBS
         if(IP(N).EQ.1) then
              np1 = np1 + 1
              DEV(N,1) = ( GMR(N,1) - GME(N,1) ) / GME(N,1)
              DEV(N,1) =  DEV(N,1) * 100.D0
         elseif(IP(N).EQ.2) then
              np2 = np2 + 1
              DEV(N,2) = ( GMR(N,2) - GME(N,2) ) / GME(N,2)
              DEV(N,2) =  DEV(N,2) * 100.D0
         endif
  906 CONTINUE
c
      DO 916 N = 1,NPTS
          ECART(N,1) = HECAL(N,1)   - HEEXP(N,1)
          ER(N,1)    = ( HECAL(N,1) - HEEXP(N,1) ) / HMAX
c          ER(N,1)    = ( HECAL(N,1) - HEEXP(N,1) ) / HEEXP(N,1)
          ER(N,1)    = ER(N,1)    * 100.D0
          ER(N,1)    = DABS( ER(N,1) )
  916 CONTINUE
c
        write(NP,*)
                      write(NP,*)'--------------------------------------
     s-------------------------'
        write(NP,207)
  207   format(/,'Temp /K',5X,'N�',4X,'x(1)',3X,'Ga1cal',5X
     s           ,'DEV1%',5X,'Ga2cal',5X,'DEV2%'/)
                      write(NP,*)'--------------------------------------
     s-------------------------'
        write(NP,*)
c
        write(NRES,*)
                      write(NRES,*)'--------------------------------------
     s------------------------------------------'
        write(NRES,307)
  307   format(/,'N�',5X,'Temp /K',6X,'x(1)',4X,'hE (exp)',7X,
     s           'hE (cal)',7X,'ECART',10X,'%'/)
                      write(NRES,*)'--------------------------------------
     s------------------------------------------'
        write(NRES,*)
        No        = 0
        DO 333 N = 1,NPTS
        WRITE(NRES,397)NUE(N,1),TE(N),XXXE(N,1),HEEXP(N,1),HECAL(N,1)
     s                   ,ECART(N,1),ER(N,1)
                  if( NUE(N,1).NE.NUE(N+1,1) ) then
               RMSD31  = 0.D0
               NPOINT = N-No
               RMSD31= 100.D0 * ( (RMSD31/(NPOINT))**0.5D0 )
               NPOINT = N-No
               No  =  N
                   write(NRES,*)
                   write(NRES,'(A19,F10.4)')'RMSD (%)         = ',RMSD31
                   write(NRES,'(A19,I5)')'Nombre de points = ',NPOINT
                      write(NRES,*)
                      write(NRES,*)'--------------------------------------
     s------------------------------------------'
                       write(NRES,*)
                  elseif( NUE(N,2).NE.NUE(N+1,2) ) then
               RMSD31    = 0.D0
             do i =No+1,N
c               RMSD31 = RMSD31 +((HEEXP(i,1) - HECAL(i,1))/HEEXP(i,1) )
     s          **2.D0
               RMSD31 = RMSD31 +((HEEXP(i,1) - HECAL(i,1))/HMAX )**2.D0
             enddo
               NPOINT = N-No
               RMSD31= 100.D0 * ( (RMSD31/(NPOINT))**0.5D0 )
               No  =  N
                   write(NRES,*)
                   write(NRES,'(A19,F10.4)')'RMSD (%)         = ',RMSD31
                   write(NRES,'(A19,I5)')'Nombre de points = ',NPOINT
                      write(NRES,*)
                      write(NRES,*)'--------------------------------------
     s------------------------------------------'
                       write(NP,*)
                  endif
  397    FORMAT(I2,6X,F6.2,F10.4,2(F12.4,3X),F10.4,5X,F10.2)
  333   CONTINUE
c
         DO 222 N = 1,NOBS
c            WRITE(NP,97)T(N),XXX(N,1),GMR(N,1),XXX(N,2),GMR(N,2)
c            goto 987
         if(IP(N).EQ.1) then
                 WRITE(NP,96)T(N),NU(N,1),XXX(N,1),GMR(N,1),DEV(N,1)
c                WRITE(NP,96)T(N),NU(N,1),XXX(N,1),GMR(N,1),GME(N,1)
                  if( NU(N,1).NE.NU(N+1,1) ) then
                      write(NP,*)
                      write(NP,*)'--------------------------------------
     s-------------------------'
c
                       write(NP,*)
                  elseif( NU(N,2).NE.NU(N+1,2) ) then
                      write(NP,*)
                      write(NP,*)'--------------------------------------
     s-------------------------'
c
                       write(NP,*)
                  endif
         elseif(IP(N).EQ.2) then
               WRITE(NP,198)T(N),NU(N,2),XXX(N,1),GMR(N,2),DEV(N,2)
c              WRITE(NP,198)T(N),NU(N,2),XXX(N,1),GMR(N,2),GME(N,2)
                  if( NU(N,1).NE.NU(N+1,1) ) then
                      write(NP,*)
                      write(NP,*)'--------------------------------------
     s-------------------------'
                       write(NP,*)
                  elseif( NU(N,2).NE.NU(N+1,2) ) then
                      write(NP,*)
                      write(NP,*)'--------------------------------------
     s-------------------------'
                       write(NP,*)
                  endif
         endif
  987    continue
   97    FORMAT(F6.2,F10.4,F10.4,F10.4,F10.4)
   96    FORMAT(F6.2,I6,F10.4,F10.4,F10.4)
  198    FORMAT(F6.2,I6,F10.4,20X,F10.4,F10.4)
  222   CONTINUE
        WRITE(NP,77)NOBS
        WRITE(NRES,77)NPTS
   77   format(/,5X,'Nombre de points total  = ',I3,/)
        RMSRD1  = 0.D0
        RMSRD2  = 0.D0
        RMSRD3  = 0.D0
c
        DO 848 N = 1,NPTS
         RMSRD3 = RMSRD3 +((HEEXP(N,1) - HECAL(N,1))/HMAX )**2
c         RMSRD3 = RMSRD3 + ((HEEXP(N,1) - HECAL(N,1))/ HEEXP(N,1))
     s       **2
    
  848   CONTINUE
        DO 838 N = 1,NOBS
            if(IP(N).EQ.1) then
               RMSRD1 = RMSRD1 + ( (GMR(N,1) - GME(N,1))/GME(N,1) )**2
            elseif(IP(N).EQ.2) then
               RMSRD2 = RMSRD2 + ( (GMR(N,2) - GME(N,2))/GME(N,2) )**2
            endif
  838   CONTINUE
c
      call system ("cls")
      write(*,*)
      write(*,*)'-------------------------------------------------------
     s---------------------'
      write(*,*)
      if(rep2.EQ.1) then
           write(*,*)' Mod�le de RASMUSSEN '
      elseif(rep2.EQ.2) then
           write(*,*)' Mod�le de GMEHLING  '
      endif
c
      write(*,*)
      write(*,909)fichs,hefic
  909 format('  Fichiers de donn�es  :  ',A12,5X,A12)
      write(*,*)
c
      WRITE(*,209)F(JS),RMS
  209 FORMAT(/,1H ,'  FMIN =',F14.5,'             SD =',E14.5)
c
        RMSRD1 = 100 * ( (RMSRD1/np1)**0.5D0 )
        RMSRD2 = 100 * ( (RMSRD2/np2)**0.5D0 )
        RMSRD3 = 100 * ( (RMSRD3/NPTS)**0.5D0 )
        write(NRES,849)RMSRD3
  849   format(//,' Ecart quadratique moyen (hE) = ',F8.4,'  % '/)
        write(NP,839)RMSRD1,RMSRD2
  839   format(//,' Ecart quadratique moyen (sur gamma 1) = ',F8.4, 
     s '  % '/,' Ecart quadratique moyen (sur gamma 2) = ',F8.4,'  % '/)
        write(*,739)RMSRD1,RMSRD2,RMSRD3
  739   format(/,'  Ecart quadratique moyen (sur gamma 1) = ',F8.4, 
     s '  % '/,'  Ecart quadratique moyen (sur gamma 2) = ',F8.4,
     s '  % '/,'  Ecart quadratique moyen (sur   hE   ) = ',F8.4,'  % ')
c
      write(*,*)
c     write(*,777)Vm(1),Vm(2)
      write(*,777)A(IPAR(1),IPAR(2)),A(IPAR(2),IPAR(1)),B(IPAR(1),IPAR(2
     s)),B(IPAR(2),IPAR(1))
c  777 format('  Nouveau param�tre : '/,5X,' Vm(1) : ',F15.5,5X,' Vm(2) :
c     s ',F15.5)
  777 format('  Nouveaux param�tres : '/,5X,' A(i,j) : ',2F15.5,/
     s,5X,' B(i,j) : ',2F15.5)
      write(*,*)
      write(*,*)'-------------------------------------------------------
     s---------------------'
        IF(NOIT-1)835,835,836
  835 CONTINUE
      WRITE(NP,1400)
      write(NP,'(A\)')'Nouveaux param�tres A(i,j) :'
      DO 2001 I = 1,NG
 2001 WRITE(NP,1402)(A(I,J),J=1,NG)
      write(NP,'(A\)')'Nouveaux param�tres B(i,j) :'
      do 2002 I = 1,NG
 2002 write(NP,1402)(B(I,J),J=1,NG)
      write(NP,'(A\)')'Nouveaux param�tres C(i,j) :'
      do 2003 I = 1,NG
 2003 write(NP,804)(C(I,J),J=1,NG)
c
      WRITE(NP,1400)
      write(NRES,'(A\)')'Nouveaux param�tres A(i,j) :'
      DO 3001 I = 1,NGR
 3001 WRITE(NRES,1402)(A(I,J),J=1,NGR)
      write(NRES,'(A\)')'Nouveaux param�tres B(i,j) :'
      do 3002 I = 1,NG
 3002 write(NRES,1402)(B(I,J),J=1,NGR)
      write(NRES,'(A\)')'Nouveaux param�tres C(i,j) :'
      do 3003 I = 1,NG
 3003 write(NRES,804)(C(I,J),J=1,NG)
c
      WRITE(NP,996)
      WRITE(NRES,996)
  804 FORMAT(20F15.5)
  996 FORMAT(//'    THE SUM OF THE SQUARED DIFFERENCES BETWEEN THE')
      IF(KRIT-1)997,997,998
  998 WRITE(NP,995)
      WRITE(NRES,995)
  995 FORMAT('      LOGARITMS TO THE')
  997 WRITE(NP,999)
      WRITE(NRES,999)
  999 FORMAT('EXPERIMENTAL AND CALCULATED ACTIVITY COEFFICIENTS')
                                                                                          
  836 CONTINUE
       !Stop
     
      END Program 
c
c  **********************************************************************
c  *                                                                    *
c  *                          SUBROUTINE PFAC3                          *
c  *                                                                    *
c  **********************************************************************
      SUBROUTINE PFAC3(RS,QS,delta,Vw,Vm)
c
      implicit double precision (A-H,O-Z)
      DIMENSION RS(20),QS(20),VPRIM(20),V(20),F(20),Vw(20)
      DIMENSION VM(20),delta(20),phi(20),GMC2(300,2)
c
      common/points/ NG,NGR,NOBS,NPTS
      common/param/ A(20,20),B(20,20),C(20,20),rl,gs,ghe,Q(20)
      common/ABC/ T(300),NM(300,2),XXX(300,2),NNY(20,20),GMAX,
     s            NME(300,2),XXV(300,2)
      common/HE/ TE(300),XXXE(300,2),HEEXP(300,2),HECAL(300,2),HMAX
     s          ,NNYE(20,20),NU(300,2),NUE(300,2)
c     common/CGM/ GME(300,2),GMC(300,2),IP(300),rep2,rep0
      common/CGM/ GME(300,2),GMC(300,2),GMR(300,2),IP(300),rep2,rep0
c
      if(rep2.EQ.1.or.rep2.EQ.3) then
         goto 10
      else
          goto 20
      endif
   10 continue
c
      DO 13 N = 1,NOBS
                  S3       =  0.D0
                  S4       =  0.D0
      DO 12 I = 1,2
             J        =  NM(N,I)
c              RS(J)    =  Vm(J) -    Vw(J)
                 S3     =  S3    +  ( RS(J)**rl ) * XXX(N,I)
                 S4     =  S4    +    Vm(J)       * XXX(N,I)
   12 continue
      DO 13 I = 1,2
              J        =  NM(N,I)
c              RS(J)    =  Vm(J) - Vw(J)
                 VPRIM(I)  = ( RS(J)**rl )   /  S3
                 GMC(N,I)  =  1  -  VPRIM(I)  +  DLOG( VPRIM(I) )
c
                 phi(I)   =  Vm(J) * XXX(N,I) / S4
                GMC2(N,I) = Vm(J) * ((1-phi(I))**2.D0)
     s                             * ((delta(1)-delta(2))**2.D0)
c                 GMC2(N,I) = GMC2(N,I) / ( 8.314D0 * T(J) )
                  GMC(N,I)  = GMC(N,I)
c
   13 CONTINUE
      RETURN
   20 continue
c
      DO 3 N = 1,NOBS
              S1 = 0.D0
              S2 = 0.D0
              S3 = 0.D0
              S4 = 0.D0
      DO 2 I = 1,2
             J   = NM(N,I)
             S4  = S4     +       Vm(J) * XXX(N,I)
             S3  = S3     + (RS(J)**rl) * XXX(N,I)
             S2  = S2     +       QS(J) * XXX(N,I)
             S1  = S1     +       RS(J) * XXX(N,I)
    2  Continue
       DO 3 I = 1,2
            J    = NM(N,I)
            F(I) = QS(J)  / S2
            V(I) = RS(J)  / S1
        phi(I)   =  Vm(J) * XXX(N,I) / S4
        VPRIM(I) = ( RS(J)**rl) / S3
        GMC(N,I) =  1 - VPRIM(I) + DLOG( VPRIM(I) ) - 5.D0 * QS(J) *
     s                   ( 1 - ( V(I) / F(I) ) + DLOG( V(I) / F(I) ) )
c
                 GMC2(N,I) = Vm(J) * ((1-phi(I))**2.D0)
     s                             * ((delta(1)-delta(2))**2.D0)
                 GMC2(N,I) = GMC2(N,I) / ( 8.314D0 * T(J) )
        GMC(N,I)  = GMC(N,I)
    3 CONTINUE
      RETURN
      END

c
c  *************************************************************************
c  *                                                                       *
c  *                        SUBROUTINE   PFAC4                             *
c  *                                                                       *
c  *************************************************************************
      SUBROUTINE PFAC4
c
      implicit double precision (A-H,O-Z)
      dimension GMOL(20),ATET(20),P(20,20),HRT(20)
      dimension ANYK(20),BNYK(20),GK(20,3)
c
      common/points/ NG,NGR,NOBS,NPTS
      common/param/ A(20,20),B(20,20),C(20,20),rl,gs,ghe,Q(20)
      common/ABC/ T(300),NM(300,2),XXX(300,2),NNY(20,20),GMAX,NME(300,2)
     s            ,XXV(300,2)
      common/HE/ TE(300),XXXE(300,2),HEEXP(300,2),HECAL(300,2),HMAX
     s          ,NNYE(20,20),NU(300,2),NUE(300,2)
      common/CGM/ GME(300,2),GMC(300,2),GMR(300,2),IP(300),rep2,rep0

      T0     =  298.15D0
      
      do ij = 1,4 
        HRT(ij) = 88430.09563
      enddo
      
      do jk = 5,6
        HRT(jk) = 13798.17848
      enddo

      do ik = 7,10
        HRT(ik) = 0.D0
      enddo

      DO 250 NR = 1,NOBS
       DO 7 I = 1,NG
        DO 7 J = 1,NG
c

      if(rep2.EQ.1) then
         
         P(I,J) = DEXP( -( A(I,J) + B(I,J)  *   ( T(NR) - T0 )
     s                  + C(I,J) * ( T(NR) * DLOG( T0/T(NR) )
     s                  +            T(NR) - T0 ) ) / T(NR) )
      elseif(rep2.EQ.2) then
        P(I,J) = DEXP( -( A(I,J) + B(I,J) * T(NR)
     s                           + C(I,J) * T(NR)**2.D0 ) / T(NR) )
      elseif(rep2.EQ.3) then
        P(I,J) = DEXP( -( A(I,J) + B(I,J) / T(NR) ) / T(NR) )
      endif

   7  Continue
      DO 105 II = 1,3
       IF(II-2)100,100,101
  100 SNYK = 0.D0
c
c-----------------------------------------------------------------------
c            PURE COMPONENT                                            !
c-----------------------------------------------------------------------
c
          J   = NM(NR,II)
      DO 12 K = 1,NG
   12  SNYK   = SNYK + FLOAT(NNY(J,K))
      DO 13 K = 1,NG
   13  GMOL(K)= FLOAT(NNY(J,K)) / SNYK
       GOTO 102
  101 SNYK = 0.D0
c
c-----------------------------------------------------------------------
c           MIXTURE                                                    !
c-----------------------------------------------------------------------
c
      DO 2 I = 1,2
         J   = NM(NR,I)
      DO 2 K = 1,NG
    2  SNYK  = SNYK + FLOAT(NNY(J,K)) * XXX(NR,I)
      DO 3 K = 1,NG
        GNYK = 0.D0
        DO 4 I = 1,2
           J    = NM(NR,I)
    4      GNYK = GNYK + FLOAT(NNY(J,K)) * XXX(NR,I)
    3 GMOL(K)= GNYK / SNYK
c
c-----------------------------------------------------------------------
c          CALCULATION OF GROUP AREA FRACTIONS                         !
c-----------------------------------------------------------------------
c
  102 SNYK = 0.D0
      DO 5 K = 1,NG
    5   SNYK = SNYK + Q(K) * GMOL(K)
      DO 6 K = 1,NG
    6 ATET(K)= Q(K) * GMOL(K) / SNYK
c
c-----------------------------------------------------------------------
c         CALCULATION OF ACTIVITY COEFFICIENTS                         !
c-----------------------------------------------------------------------
c
      DO 9  K   = 1,NG
                 ANYK(K)   = 0.D0
      DO 10 M = 1,NG
                 SNYK      = 0.D0
      DO 8 N = 1,NG
                 SNYK      = SNYK    + ATET(N)*P(N,M)
    8 CONTINUE
                 ANYK(K)   = ANYK(K) + ATET(M)*P(K,M)  / SNYK
   10 continue
                 BNYK(K)   = 0.D0
      DO 11 M = 1,NG
   11            BNYK(K)   = BNYK(K) + ATET(M)*P(M,K)
                 BNYK(K)   = DLOG( BNYK(K) )
                 GK(K,II)  = Q(K)    * ( 1.D0 - BNYK(K) - ANYK(K) )
    9 continue
c
  105 CONTINUE
      DO 250 I = 1,2
                 J   = NM(NR,I)
                 S8  = 0.D0
      DO 248 K = 1,NG
                 S8  = S8  +  FLOAT( NNY(J,K) ) * ( GK(K,3) - GK(K,I) )
  248  continue
                 GMR(NR,I)  = S8
  250 CONTINUE
      RETURN
      END
c
c  *************************************************************************
c  *                                                                       *
c  *                        SUBROUTINE   HE                                *
c  *                                                                       *
c  *************************************************************************
      SUBROUTINE HE
c
      implicit double precision (A-H,O-Z)
      dimension GMOL(20),ATET(20),P(20,20),TDGAMAT(20,3)
c
      common/points/ NG,NGR,NOBS,NPTS
      common/param/ A(20,20),B(20,20),C(20,20),rl,gs,ghe,Q(20)
      common/ABC/ T(300),NM(300,2),XXX(300,2),NNY(20,20),GMAX,NME(300,2)
     s            ,XXV(300,2)
      common/HE/ TE(300),XXXE(300,2),HEEXP(300,2),HECAL(300,2),HMAX
     s          ,NNYE(20,20),NU(300,2),NUE(300,2)
c     common/CGM/ GME(300,2),GMC(300,2),IP(300),rep2,rep0
      common/CGM/ GME(300,2),GMC(300,2),GMR(300,2),IP(300),rep2,rep0

      T0  =  298.15D0
c
      DO 250 NR = 1,NPTS
       DO 7 I = 1,NGR
        DO 7 J = 1,NGR
c
      if(rep2.EQ.1) then
        P(I,J) = DEXP( -( A(I,J) + B(I,J)  *   ( TE(NR) - T0 )
     s                  + C(I,J) * ( TE(NR) * DLOG( T0/TE(NR) )
     s                  +             TE(NR) - T0 ) ) / TE(NR) )
      elseif(rep2.EQ.2) then
        P(I,J) = DEXP( -( A(I,J) + B(I,J) * TE(NR)
     s                            + C(I,J) * TE(NR)**2.D0 ) / TE(NR) )
      elseif(rep2.EQ.3) then
        P(I,J) = DEXP( -( A(I,J) + B(I,J) / TE(NR) ) / TE(NR) )
      endif
c
   7  Continue
c
c-----------------------------------------------------------------------
c             CALCULATION OF GROUP MOLE FRACTIONS                      !
c-----------------------------------------------------------------------
c
      DO 105 II = 1,3
       IF(II-2)100,100,101
  100 SNYK = 0.D0
c
c-----------------------------------------------------------------------
c            PURE COMPONENT                                            !
c-----------------------------------------------------------------------
c
      J = NME(NR,II)
      DO 12 K = 1,NGR
   12  SNYK   = SNYK + FLOAT(NNYE(J,K))
      DO 13 K = 1,NGR
   13  GMOL(K)= FLOAT(NNYE(J,K)) / SNYK
       GOTO 102
  101 SNYK = 0.D0
c
c-----------------------------------------------------------------------
c           MIXTURE                                                    !
c-----------------------------------------------------------------------
c
      DO 2 I = 1,2
           J = NME(NR,I)
      DO 2 K = 1,NGR
    2  SNYK  = SNYK + FLOAT(NNYE(J,K)) * XXXE(NR,I)
      DO 3 K = 1,NGR
        GNYK = 0.D0
        DO 4 I = 1,2
           J = NME(NR,I)
    4      GNYK = GNYK + FLOAT(NNYE(J,K)) * XXXE(NR,I)
    3 GMOL(K)= GNYK / SNYK
c
c-----------------------------------------------------------------------
c          CALCULATION OF GROUP AREA FRACTIONS                         !
c-----------------------------------------------------------------------
c
  102 SNYK = 0.D0
      DO 5 K = 1,NGR
    5   SNYK = SNYK + Q(K) * GMOL(K)
      DO 6 K = 1,NGR
    6 ATET(K)= Q(K) * GMOL(K) / SNYK
c
c  ---------------------------------------------------------------------
c  |              CALCULATION OF EXCESS ENTHALPIES                     |
c  ---------------------------------------------------------------------
c
      DO 9  K   = 1,NGR
        S  6      = 0.D0
c
      DO 10 M = 1,NGR
        S1      = 0.D0
        S2      = 0.D0
        S3      = 0.D0
c
      DO 8 N = 1,NGR
        S1    = S1   + ATET(N) * P(N,K)
        S3    = S3   + ATET(N) * P(N,M)
   
    8 CONTINUE
        if(rep2.EQ.1) then
             S6      = S6 + ATET(M) * ( ( P(M,K) * ( B(M,K) +  C(M,K)
     s                    * DLOG( T0 / TE(NR) ) + DLOG( P(M,K) ) ) / S1)
     s                    + ( P(K,M) * S2 / ( S3**2.D0) ) )
        elseif(rep2.EQ.2) then
            S6      = S6 + ATET(M) * ( ( P(M,K) * ( B(M,K) +  C(M,K)
     s                   * 2.D0 * TE(NR) + DLOG( P(M,K) ) ) / S1 )
     s                   + ( P(K,M) * S2 / ( S3**2.D0) ) )
        elseif(rep2.EQ.3) then
            S6      = S6 + ATET(M) * (  P(M,K) * (  DLOG( P(M,K) )
     s                   - B(M,K) / TE(NR)**2.D0 )  / S1
     s                   + P(K,M) * S2 /  S3**2.D0 )
        endif
   10 continue
    9 TDGAMAT(K,II) = Q(K) * S6
c
  105 CONTINUE
      S7   = 0.D0
      do 249 I = 1,2
         SNYK = 0.D0
           J  = NME(NR,I)
          DO 248 K = 1,NGR
  248 SNYK = SNYK + FLOAT( NNYE(J,K) ) * XXXE(NR,I) * ( TDGAMAT(K,3)
     s                                            -   TDGAMAT(K,I) )
      S7          = S7      + SNYK
  249 HECAL(NR,1) = (-1.D0) * (8.314D0) * TE(NR) * S7
  250 CONTINUE
      RETURN
      END
c
c  *************************************************************************
c  *                                                                       *
c  *                          SUBROUTINE   FMIN                            *
c  *                                                                       *
c  *************************************************************************
c
      SUBROUTINE FMIN(NPAR,IPAR,XX,FF,KRIT,JENS,IDEN)
c
      implicit double precision (A-H,O-Z)
      dimension IPAR(10),XX(10),JENS(20)
c
      common/points/ NG,NGR,NOBS,NPTS
      common/param/ A(20,20),B(20,20),C(20,20),rl,gs,ghe,Q(20)
      common/ABC/ T(300),NM(300,2),XXX(300,2),NNY(20,20),GMAX,NME(300,2)
     s            ,XXV(300,2)
      common/HE/ TE(300),XXXE(300,2),HEEXP(300,2),HECAL(300,2),HMAX
     s          ,NNYE(20,20),NU(300,2),NUE(300,2)
c     common/CGM/ GME(300,2),GMC(300,2),IP(300),rep2,rep0
      common/CGM/ GME(300,2),GMC(300,2),GMR(300,2),IP(300),rep2,rep0

c
      DO 2 I = 1,NPAR
         KI       = IPAR(I)
         KJ       = IPAR(I+1)
        if(rep0.EQ.0) then
               A(KI,KJ) = XX(I)
            if(NPAR.GT.1) then
               A(KJ,KI) = XX(I+1)
            endif
        elseif(rep0.EQ.1) then
             B(KI,KJ) = XX(I)
             B(KJ,KI) = XX(I+1)
        elseif(rep0.EQ.2) then
             A(KI,KJ) = XX(I)
             A(KJ,KI) = XX(I+1)
             B(KI,KJ) = XX(I+2)
             B(KJ,KI) = XX(I+3)
        endif
    2 continue
      IF(IDEN)9,9,8
c
    8      KKI    = JENS(1)
           KKJ    = JENS(2)
             DO 7 J  = 3,IDEN,2
               IKI        = JENS(J)
               IKJ        = JENS(J+1)
                  if(rep0.EQ.0) then
                     A(IKI,IKJ) = A(KKI,KKJ)
                     A(IKJ,IKI) = A(KKJ,KKI)
                  elseif(rep0.EQ.1) then
                     B(IKI,IKJ) = B(KKI,KKJ)
                     B(IKJ,IKI) = B(KKJ,KKI)
                  elseif(rep0.EQ.2) then
                     A(IKI,IKJ) = A(KKI,KKJ)
                     A(IKJ,IKI) = A(KKJ,KKI)
                     B(IKI,IKJ) = B(KKI,KKJ)
                     B(IKJ,IKI) = B(KKJ,KKI)
                  endif
    7        CONTINUE
    9 CONTINUE
c
        CALL PFAC4
        CALL HE
c
      do 200 NR = 1,NOBS
          do 200 I =1,2
              GMR(NR,I) = GMC(NR,I) + GMR(NR,I)
              GMR(NR,I) = DEXP(GMR(NR,I))
  200 continue
c
      FF  = 0.D0
      FF1 = 0.D0
      FF2 = 0.D0
      FF3 = 0.D0
      DO 3 N = 1,NOBS
      DO 3 I =1,2
   10  IF(KRIT-1)10,10,20
   
       if(IP(N).EQ.1)then
             FF1 = FF1 + ( ( GMR(N,1) - GME(N,1) ) / GME(N,1) )**2.D0
             FF1 = FF1 * gs
       elseif(IP(N).EQ.2) then
             FF2 = FF2 + ( ( GMR(N,2) - GME(N,2) ) / GME(N,2) )**2.D0
             FF2 = FF2 * gs
       endif
      GOTO 3
c
   20  if(IP(N).EQ.1)then
             FF1 = FF1 + ( ( DLOG( GMR(N,1) ) - DLOG(GME(N,1) ) ) /
     s                       DLOG(GME(N,1) ) )**2.D0
       
       elseif(IP(N).EQ.2) then
             FF2 = FF2 + ( ( DLOG( GMR(N,2) ) - DLOG(GME(N,2) ) ) /
     s                       DLOG(GME(N,2) ) )**2.D0
       endif
    3 CONTINUE
c
      DO 13 N = 1,NPTS
      DO 13 I = 1,2
          FF3 = FF3 + ( ( HEEXP(N,1) - HECAL(N,1) ) / HEEXP(N,1))**2.D0
          FF3 = FF3 * ghe
  13  continue
       if(rep0.EQ.0) then
         FF3 = 0.D0
       elseif (rep0.EQ.1) then
         FF1 = 0.D0
         FF2 = 0.D0
       endif
c        FF  =0.0D0*( FF1 + FF2 ) / NOBS  +  FF3 / NPTS
         FF  =( FF1 + FF2 ) / NOBS  +  FF3 / NPTS
      RETURN
      END
     
Comments: