Add elastic demand to scbd (fixed frequencies)

src/main.f

@@ -1,5 +1,8 @@
       PROGRAM SCBD
       IMPLICIT REAL*8 (A-H,O-Z)
+      INTEGER IT,ITMAX
+      DOUBLE PRECISION DEMAND, QUADERR
+      DOUBLE PRECISION MUSC, MUCE, MUSE, POPSC0, POPCE0, POPSE0, EPS
       COMMON /PB/ NPROB
 
       INCLUDE "param.inc"
@@ -10,7 +13,10 @@
     
       ZERO=0.0D0
       ONE=1.0D0
-      TWO=2.0D0
+      TWO=2.0D0 
+
+      EPS = 1.0D-5
+      ITMAX = 100
 
       CALL  PARAM("scbd.dat",PSC,PCE,PSE,BB,TAUSC,TAUCE,TAUSE,
      1  FSC,FCE,FSE,AA) 
@@ -22,6 +28,14 @@
       ASER = VAL * ASER
 
 
+      POPSC0 = POPSC
+      POPCE0 = POPCE
+      POPSE0 = POPSE
+
+      MUSC = 0.01D0
+      MUCE = 0.01D0
+      MUSE = 0.01D0 
+
       DO 1, I=1,90
       LINE0(I)='-'
       LINE1(I)='='
@@ -84,16 +98,13 @@ C  1 -- COMPUTE THE OPTIMUM
 C  ------------------------
 
       ID=1
-      CALL OPTIMUM (TC,USCT,UCET,USET,FSC,FCE,FSE,AA,
-     1              PSC,PCE,BB,TAUSC,TAUCE,TAUSE,ID)
 
-      WRITE (*, '(A12,35X,F8.4)' ) 'Optimum', TC
-      POPTRAIN = USCT+UCET+USET
-      POPCAR = POP-POPTRAIN
-      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
-     &   USCT,UCET,USET,POPTRAIN,
-     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
+      IT = 0
+
+ 305  CONTINUE
 
+      CALL OPTIMUM (TC,USCT,UCET,USET,FSC,FCE,FSE,AA,
+     1              PSC,PCE,BB,TAUSC,TAUCE,TAUSE,ID) 
       C1T=FSCT+CW/(TWO*FSC)+ASCT*(USCT+USET)/FSC
       C2T=FCET+CW/(TWO*FCE)+ACET*(UCET+USET)/FCE 
       C3T=C1T+C2T-AA*CW/(TWO*FSC)+(ONE-AA)*SC
@@ -103,14 +114,40 @@ C  ------------------------
       C3R=FSER+ASER*(POPSE-USET)
       CR =((POPSC-USCT)*C1R+(POPCE-UCET)*C2R+(POPSE-USET)*C3R)/
      &    (POPSC+POPCE+POPSE-USCT-UCET-USET) 
+      POPSCOLD = POPSC
+      POPCEOLD = POPCE
+      POPSEOLD = POPSE
+      POPSC = DEMAND(POPSC0, MUSC , FSCT, C1T)
+      POPCE = DEMAND(POPCE0, MUCE , FCET, C2T)
+      POPSE = DEMAND(POPSE0, MUSE , FSCT+FCET, C3T)
+      ERR = ZERO
+      ERR = ERR + QUADERR(POPSC,POPSCOLD)
+      ERR = ERR + QUADERR(POPCE,POPCEOLD)
+      ERR = ERR + QUADERR(POPSE,POPSEOLD)
+      IT = IT + 1
+      IF(IT.LE.ITMAX.AND.ERR.GT.EPS) GOTO 305
+
+      WRITE (*, '(A12,35X,F8.4)' ) 'Optimum', TC
+      POPTRAIN = USCT+UCET+USET
+      POPCAR = POP-POPTRAIN
+      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
+     &   USCT,UCET,USET,POPTRAIN,
+     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
+
       WRITE (*, '(11X,A5,4F8.2,5X,4F8.2,5X)')  'Costs', 
      &   C1T,C2T,C3T,CT,C1R,C2R,C3R,CR
       WRITE (*, '(11X,A5,4F8.2)') 'Fares',PSC,PCE,BB*(PSC+PCE),
      &     PSC*USCT+PCE*UCET+PSE*USET
       WRITE (*, '(11X,A5,37X,4F8.2)') 'Tolls',TAUSC,TAUCE,TAUSE,
      &  (POPSC-USCT)*TAUSC+(POPCE-UCET)*TAUCE+(POPSE-USET)*TAUSE
+      PRINT*, 'Populations --> ', POPSC, POPCE, POPSE, IT
       PRINT '(90A)', LINE0
 
+C Reset population
+      POPSC = POPSC0
+      POPCE = POPCE0
+      POPSE = POPSE0 
+
 C  ------------------------
 C  2 -- BETA LIMITED: BETA <= BETA_MAX
 C  ------------------------
@@ -153,14 +190,14 @@ CC  2 -- UNPRICED EQUILIBRIUM
 CC  -------------------------
 
       ID = 2
+
+      IT = 0
+ 300  CONTINUE 
+C     PRINT*, "-->", IT, POPSC, POPCE, POPSE, ERR
       CALL EQUILIBRIUM (TC,USCT,UCET,USET,FSC,FCE,FSE,AA,
      1     ZERO,ZERO,ZERO,ZERO,ZERO,ZERO,ID) 
       POPTRAIN = USCT+UCET+USET
       POPCAR = POP-POPTRAIN
-      WRITE (*, '(A12,35X,F8.4)' ) 'Unpriced', TC
-      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
-     &   USCT,UCET,USET,POPTRAIN,
-     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
 
       C1T=FSCT+CW/(TWO*FSC)+ASCT*(USCT+USET)/FSC
       C2T=FCET+CW/(TWO*FCE)+ACET*(UCET+USET)/FCE 
@@ -171,10 +208,33 @@ CC  -------------------------
       C3R=FSER+ASER*(POPSE-USET)
       CR =((POPSC-USCT)*C1R+(POPCE-UCET)*C2R+(POPSE-USET)*C3R)/
      &    (POPSC+POPCE+POPSE-USCT-UCET-USET) 
+      POPSCOLD = POPSC
+      POPCEOLD = POPCE
+      POPSEOLD = POPSE
+      POPSC = DEMAND(POPSC0, MUSC , FSCT, C1T)
+      POPCE = DEMAND(POPCE0, MUCE , FCET, C2T)
+      POPSE = DEMAND(POPSE0, MUSE , FSCT+FCET, C3T)
+      ERR = ZERO
+      ERR = ERR + QUADERR(POPSC,POPSCOLD)
+      ERR = ERR + QUADERR(POPCE,POPCEOLD)
+      ERR = ERR + QUADERR(POPSE,POPSEOLD)
+      IT = IT + 1
+      IF(IT.LE.ITMAX.AND.ERR.GT.EPS) GOTO 300
+
+      WRITE (*, '(A12,35X,F8.4)' ) 'Unpriced', TC
+      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
+     &   USCT,UCET,USET,POPTRAIN,
+     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
       WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Costs', 
      &   C1T,C2T,C3T,CT,C1R,C2R,C3R,CR
+      PRINT*, 'Populations --> ', POPSC, POPCE, POPSE, IT
       PRINT '(90A)', LINE0
 
+C Reset population
+      POPSC = POPSC0
+      POPCE = POPCE0
+      POPSE = POPSE0 
+
 CC  -----------------------
 CC  2 -- PRICED EQUILIBRIUM
 CC  -----------------------
@@ -194,15 +254,13 @@ CC  -----------------------
 CC  2 -- ROAD PROFIT
 CC  -----------------------
 
+      IT = 0
+
+ 310  CONTINUE 
+
       CALL PRFROAD (PRFR,TC,USCT,UCET,USET,FSC,FCE,AA,
      1     ZERO,ZERO,ZERO,TAUSC,TAUCE,TAUSE) 
 C     WRITE (*,'(8F9.4)') TC, USCT,UCET,USET,PRFR,TAUSC,TAUCE,TAUSE 
-      WRITE (*, '(A12,35X,F8.4)' ) 'Road prf', TC
-      POPTRAIN = USCT+UCET+USET
-      POPCAR = POP-POPTRAIN
-      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
-     &   USCT,UCET,USET,POPTRAIN,
-     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
 
       C1T=FSCT+CW/(TWO*FSC)+ASCT*(USCT+USET)/FSC
       C2T=FCET+CW/(TWO*FCE)+ACET*(UCET+USET)/FCE 
@@ -213,23 +271,47 @@ C     WRITE (*,'(8F9.4)') TC, USCT,UCET,USET,PRFR,TAUSC,TAUCE,TAUSE
       C3R=FSER+ASER*(POPSE-USET)
       CR =((POPSC-USCT)*C1R+(POPCE-UCET)*C2R+(POPSE-USET)*C3R)/
      &    (POPSC+POPCE+POPSE-USCT-UCET-USET) 
+
+      POPSCOLD = POPSC
+      POPCEOLD = POPCE
+      POPSEOLD = POPSE
+      POPSC = DEMAND(POPSC0, MUSC , FSCT, C1T)
+      POPCE = DEMAND(POPCE0, MUCE , FCET, C2T)
+      POPSE = DEMAND(POPSE0, MUSE , FSCT+FCET, C3T)
+      ERR = ZERO
+      ERR = ERR + QUADERR(POPSC,POPSCOLD)
+      ERR = ERR + QUADERR(POPCE,POPCEOLD)
+      ERR = ERR + QUADERR(POPSE,POPSEOLD)
+      IT = IT + 1
+      IF(IT.LE.ITMAX.AND.ERR.GT.EPS) GOTO 310
+
+      WRITE (*, '(A12,35X,F8.4)' ) 'Road prf', TC
+      POPTRAIN = USCT+UCET+USET
+      POPCAR = POP-POPTRAIN
+      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
+     &   USCT,UCET,USET,POPTRAIN,
+     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
       WRITE (*, '(11X,A5,4F8.2,5X,4F8.2,5X)')  'Costs', 
      &   C1T,C2T,C3T,CT,C1R,C2R,C3R,CR
       WRITE (*, '(11X,A5,37X,4F8.2)') 'Tolls',TAUSC,TAUCE,TAUSE,PRFR
+      PRINT*, 'Populations --> ', POPSC, POPCE, POPSE, IT
       PRINT '(90A)', LINE0
 
+C Reset population
+      POPSC = POPSC0
+      POPCE = POPCE0
+      POPSE = POPSE0 
+
 CC  --------------------
 CC  2 -- TRAIN PROFIT
 CC  --------------------
 
+      IT = 0
+
+ 320  CONTINUE 
+
       CALL PRFTRAIN (PRFT,TC,USCT,UCET,USET,FSC,FCE,AA,
      1     PSC,PCE,BB,ZERO,ZERO,ZERO) 
-      WRITE (*, '(A12,35X,F8.4)' ) 'Train prf', TC
-      POPTRAIN = USCT+UCET+USET
-      POPCAR = POP-POPTRAIN
-      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
-     &   USCT,UCET,USET,POPTRAIN,
-     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
 
       C1T=FSCT+CW/(TWO*FSC)+ASCT*(USCT+USET)/FSC
       C2T=FCET+CW/(TWO*FCE)+ACET*(UCET+USET)/FCE 
@@ -240,25 +322,50 @@ CC  --------------------
       C3R=FSER+ASER*(POPSE-USET)
       CR =((POPSC-USCT)*C1R+(POPCE-UCET)*C2R+(POPSE-USET)*C3R)/
      &    (POPSC+POPCE+POPSE-USCT-UCET-USET) 
+
+      POPSCOLD = POPSC
+      POPCEOLD = POPCE
+      POPSEOLD = POPSE
+      POPSC = DEMAND(POPSC0, MUSC , FSCT, C1T + PSC)
+      POPCE = DEMAND(POPCE0, MUCE , FCET, C2T + PCE)
+      POPSE = DEMAND(POPSE0, MUSE , FSCT+FCET, C3T + PSC + PCE )
+      ERR = ZERO
+      ERR = ERR + QUADERR(POPSC,POPSCOLD)
+      ERR = ERR + QUADERR(POPCE,POPCEOLD)
+      ERR = ERR + QUADERR(POPSE,POPSEOLD)
+      IT = IT + 1
+      IF(IT.LE.ITMAX.AND.ERR.GT.EPS) GOTO 320 
+
+      WRITE (*, '(A12,35X,F8.4)' ) 'Train prf', TC
+      POPTRAIN = USCT+UCET+USET
+      POPCAR = POP-POPTRAIN
+      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
+     &   USCT,UCET,USET,POPTRAIN,
+     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
       WRITE (*, '(11X,A5,4F8.2,5X,4F8.2,5X)')  'Costs', 
      &   C1T,C2T,C3T,CT,C1R,C2R,C3R,CR
       WRITE (*, '(11X,A5,4F8.2)') 'Fares',PSC,PCE,BB*(PSC+PCE),PRFT
+      PRINT*, 'Populations --> ', POPSC, POPCE, POPSE, IT
       PRINT '(90A)', LINE0 
 
+C Reset population
+      POPSC = POPSC0
+      POPCE = POPCE0
+      POPSE = POPSE0 
+
+
 CC  --------------------
 CC  2 -- SEMI-PUBLIC
 CC  --------------------
       PSC = ZERO
       PCE = ZERO
       PSE = ZERO
+
+      IT = 0
+ 330  CONTINUE
+
       CALL SEMIPUBLIC (PRFT,PRFR,TC,USCT,UCET,USET,FSC,FCE,AA,
      1                 PSC,PCE,BB,TAUSC,TAUCE,TAUSE)
-      WRITE (*, '(A12,35X,F8.4)' ) 'Semi-pub', TC
-      POPTRAIN = USCT+UCET+USET
-      POPCAR = POP-POPTRAIN
-      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
-     &   USCT,UCET,USET,POPTRAIN,
-     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
 
       C1T=FSCT+CW/(TWO*FSC)+ASCT*(USCT+USET)/FSC
       C2T=FCET+CW/(TWO*FCE)+ACET*(UCET+USET)/FCE 
@@ -269,24 +376,48 @@ CC  --------------------
       C3R=FSER+ASER*(POPSE-USET)
       CR =((POPSC-USCT)*C1R+(POPCE-UCET)*C2R+(POPSE-USET)*C3R)/
      &    (POPSC+POPCE+POPSE-USCT-UCET-USET) 
+
+      POPSCOLD = POPSC
+      POPCEOLD = POPCE
+      POPSEOLD = POPSE
+      POPSC = DEMAND(POPSC0, MUSC , FSCT, C1T+PSC)
+      POPCE = DEMAND(POPCE0, MUCE , FCET, C2T+PCE)
+      POPSE = DEMAND(POPSE0, MUSE , FSCT+FCET, C3T+PSC+PCE)
+      ERR = ZERO
+      ERR = ERR + QUADERR(POPSC,POPSCOLD)
+      ERR = ERR + QUADERR(POPCE,POPCEOLD)
+      ERR = ERR + QUADERR(POPSE,POPSEOLD)
+      IT = IT + 1
+      IF(IT.LE.ITMAX.AND.ERR.GT.EPS) GOTO 330
+
+      WRITE (*, '(A12,35X,F8.4)' ) 'Semi-pub', TC
+      POPTRAIN = USCT+UCET+USET
+      POPCAR = POP-POPTRAIN
+      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
+     &   USCT,UCET,USET,POPTRAIN,
+     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
       WRITE (*, '(11X,A5,4F8.2,5X,4F8.2,5X)')  'Costs', 
      &   C1T,C2T,C3T,CT,C1R,C2R,C3R,CR
       WRITE (*, '(11X,A5,4F8.2)') 'Fares',PSC,PCE,BB*(PSC+PCE),PRFT
       WRITE (*, '(11X,A5,37X,4F8.2)') 'Tolls',TAUSC,TAUCE,TAUSE,PRFR
+      PRINT*, 'Populations --> ', POPSC, POPCE, POPSE, IT
       PRINT '(90A)', LINE0
 
+C Reset population
+      POPSC = POPSC0
+      POPCE = POPCE0
+      POPSE = POPSE0 
+
 CC  --------------------
 CC  2 -- DUOPOLY
 CC  --------------------
 
+      IT = 0
+
+ 340  CONTINUE
+
       CALL DUOPOLY (PRFT,PRFR,TC,USCT,UCET,USET,FSC,FCE,AA,
      1              PSC,PCE,BB,TAUSC,TAUCE,TAUSE)
-      WRITE (*, '(A12,35X,F8.4)' ) 'Duopoly', TC
-      POPTRAIN = USCT+UCET+USET
-      POPCAR = POP-POPTRAIN
-      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
-     &   USCT,UCET,USET,POPTRAIN,
-     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
 
       C1T=FSCT+CW/(TWO*FSC)+ASCT*(USCT+USET)/FSC
       C2T=FCET+CW/(TWO*FCE)+ACET*(UCET+USET)/FCE 
@@ -297,30 +428,55 @@ CC  --------------------
       C3R=FSER+ASER*(POPSE-USET)
       CR =((POPSC-USCT)*C1R+(POPCE-UCET)*C2R+(POPSE-USET)*C3R)/
      &    (POPSC+POPCE+POPSE-USCT-UCET-USET) 
+
+      POPSCOLD = POPSC
+      POPCEOLD = POPCE
+      POPSEOLD = POPSE
+      POPSC = DEMAND(POPSC0, MUSC , FSCT, C1T+PSC)
+      POPCE = DEMAND(POPCE0, MUCE , FCET, C2T+PCE)
+      POPSE = DEMAND(POPSE0, MUSE , FSCT+FCET, C3T+PSC+PCE)
+      ERR = ZERO
+      ERR = ERR + QUADERR(POPSC,POPSCOLD)
+      ERR = ERR + QUADERR(POPCE,POPCEOLD)
+      ERR = ERR + QUADERR(POPSE,POPSEOLD)
+      IT = IT + 1
+      IF(IT.LE.ITMAX.AND.ERR.GT.EPS) GOTO 340 
+
+      WRITE (*, '(A12,35X,F8.4)' ) 'Duopoly', TC
+      POPTRAIN = USCT+UCET+USET
+      POPCAR = POP-POPTRAIN
+      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
+     &   USCT,UCET,USET,POPTRAIN,
+     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
       WRITE (*, '(11X,A5,4F8.2,5X,4F8.2,5X)')  'Costs', 
      &   C1T,C2T,C3T,CT,C1R,C2R,C3R,CR
       WRITE (*, '(11X,A5,4F8.2)') 'Fares',PSC,PCE,BB*(PSC+PCE),PRFT
       WRITE (*, '(11X,A5,37X,4F8.2)') 'Tolls',TAUSC,TAUCE,TAUSE,PRFR
-
+      PRINT*, 'Populations --> ', POPSC, POPCE, POPSE, IT 
       PRINT '(90A)', LINE3 
+
       WRITE (*, '(16X,4A8,5X,4A8)') 'SC','CE','SE','AGR', 
      &                              'SC','CE','SE','AGR' 
       PRINT '(90A)', LINE0
+
+C Reset population
+      POPSC = POPSC0
+      POPCE = POPCE0
+      POPSE = POPSE0 
+
 C  ---------------------------------------------
 C  1 -- COMPUTE THE OPTIMUM WITH DIRECT TRAIN SE
 C  ---------------------------------------------
 
       ID=10
+
+      IT = 0
+
+ 400  CONTINUE
+
       CALL OPTIMUM (TC,USCT,UCET,USET,FSC,FCE,FSE,AA,
      1              PSC,PCE,BB,TAUSC,TAUCE,TAUSE,ID)
 
-      WRITE (*, '(A12,35X,F8.4)' ) 'Optimum', TC
-      POPTRAIN = USCT+UCET+USET
-      POPCAR = POP-POPTRAIN
-      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
-     &   USCT,UCET,USET,POPTRAIN,
-     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
-
       C1T=FSCT+CW/(TWO*FSC)+ASCT*USCT/FSC
       C2T=FCET+CW/(TWO*FCE)+ACET*UCET/FCE 
       C3T=FSET+CW/(TWO*FSE)+ASET*USET/FCE 
@@ -330,27 +486,52 @@ C  ---------------------------------------------
       C3R=FSER+ASER*(POPSE-USET)
       CR =((POPSC-USCT)*C1R+(POPCE-UCET)*C2R+(POPSE-USET)*C3R)/
      &    (POPSC+POPCE+POPSE-USCT-UCET-USET) 
+
+      POPSCOLD = POPSC
+      POPCEOLD = POPCE
+      POPSEOLD = POPSE
+      POPSC = DEMAND(POPSC0, MUSC , FSCT, C1T)
+      POPCE = DEMAND(POPCE0, MUCE , FCET, C2T)
+      POPSE = DEMAND(POPSE0, MUSE , FSCT+FCET, C3T)
+      ERR = ZERO
+      ERR = ERR + QUADERR(POPSC,POPSCOLD)
+      ERR = ERR + QUADERR(POPCE,POPCEOLD)
+      ERR = ERR + QUADERR(POPSE,POPSEOLD)
+      IT = IT + 1
+      IF(IT.LE.ITMAX.AND.ERR.GT.EPS) GOTO 400
+
+      WRITE (*, '(A12,35X,F8.4)' ) 'Optimum', TC
+      POPTRAIN = USCT+UCET+USET
+      POPCAR = POP-POPTRAIN
+      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
+     &   USCT,UCET,USET,POPTRAIN,
+     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
+
       WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Costs', 
      &   C1T,C2T,C3T,CT,C1R,C2R,C3R,CR
       WRITE (*, '(11X,A5,4F8.2)') 'Fares',PSC,PCE,BB,
      &     PSC*USCT+PCE*UCET+PSE*USET
       WRITE (*, '(11X,A5,37X,4F8.2)') 'Tolls',TAUSC,TAUCE,TAUSE,
      &  (POPSC-USCT)*TAUSC+(POPCE-UCET)*TAUCE+(POPSE-USET)*TAUSE
+      PRINT*, 'Populations --> ', POPSC, POPCE, POPSE, IT
       PRINT '(90A)', LINE0
 
+C Reset population
+      POPSC = POPSC0
+      POPCE = POPCE0
+      POPSE = POPSE0 
+
 CC  ----------------------------------------------
 CC  2 -- UNPRICED EQUILIBRIUM WITH DIRECT TRAIN SE
 CC  ----------------------------------------------
 
       ID = 11
+
+      IT = 0
+
+ 410  CONTINUE 
       CALL EQUILIBRIUM (TC,USCT,UCET,USET,FSC,FCE,FSE,AA,
      1     ZERO,ZERO,ZERO,ZERO,ZERO,ZERO,ID) 
-      WRITE (*, '(A12,35X,F8.4)' ) 'Unpriced', TC
-      POPTRAIN = USCT+UCET+USET
-      POPCAR = POP-POPTRAIN
-      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
-     &   USCT,UCET,USET,POPTRAIN,
-     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
 
       C1T=FSCT+CW/(TWO*FSC)+ASCT*USCT/FSC
       C2T=FCET+CW/(TWO*FCE)+ACET*UCET/FCE 
@@ -361,16 +542,56 @@ CC  ----------------------------------------------
       C3R=FSER+ASER*(POPSE-USET)
       CR =((POPSC-USCT)*C1R+(POPCE-UCET)*C2R+(POPSE-USET)*C3R)/
      &    (POPSC+POPCE+POPSE-USCT-UCET-USET) 
+      POPSCOLD = POPSC
+      POPCEOLD = POPCE
+      POPSEOLD = POPSE
+      POPSC = DEMAND(POPSC0, MUSC , FSCT, C1T)
+      POPCE = DEMAND(POPCE0, MUCE , FCET, C2T)
+      POPSE = DEMAND(POPSE0, MUSE , FSCT+FCET, C3T)
+      ERR = ZERO
+      ERR = ERR + QUADERR(POPSC,POPSCOLD)
+      ERR = ERR + QUADERR(POPCE,POPCEOLD)
+      ERR = ERR + QUADERR(POPSE,POPSEOLD)
+      IT = IT + 1
+      IF(IT.LE.ITMAX.AND.ERR.GT.EPS) GOTO 410
+
+
+      WRITE (*, '(A12,35X,F8.4)' ) 'Unpriced', TC
+      POPTRAIN = USCT+UCET+USET
+      POPCAR = POP-POPTRAIN
+      WRITE (*, '(11X,A5,4F8.2,5X,4F8.2)')  'Users', 
+     &   USCT,UCET,USET,POPTRAIN,
+     &   POPSC-USCT,POPCE-UCET,POPSE-USET,POPCAR 
+
       WRITE (*, '(11X,A5,4F8.2,5X,4F8.2,5X)')  'Costs', 
      &   C1T,C2T,C3T,CT,C1R,C2R,C3R,CR
+      PRINT*, 'Populations --> ', POPSC, POPCE, POPSE, IT
       PRINT '(90A)', LINE1
 
+C Reset population (No need at the end of the program)
+C     POPSC = POPSC0
+C     POPCE = POPCE0
+C     POPSE = POPSE0 
 
       END
 
 
+      FUNCTION DEMAND(POP, MU, C0, C)
+      IMPLICIT NONE
+      DOUBLE PRECISION DEMAND, POP, MU, C0, C 
 
+      DEMAND = POP ** ( 1.0D0 - MU * ( C/C0 - 1.0D0) )
+
+      RETURN
+      END
 
 
+      FUNCTION QUADERR(X1, X2)
+      IMPLICIT NONE
+      DOUBLE PRECISION QUADERR, X1, X2
 
+      QUADERR = DSQRT( (X1-X2) ** 2 )
+
+      RETURN
+      END