bordeaux-model/brd.f

      PROGRAM BORDEAUX
C----------------------------------------------------------- 

      IMPLICIT NONE
      REAL*8 AA(21), P(21), T(21), T0(21), Q(21), S(21), Y(21), VAR(21)
      REAL*8 COU(21),DIS(21),VIT(21)
      REAL*8 CES2, PCES2,ALPHA,A,BPR
      REAL*8 VAL,VOL,DELTAVOL,DELTAP,CAP(8)
      REAL*8 POL(8),CONG(8), OCC(8),TOTPOL,TOTCONG,TOTPOL0,TOTCONG0
      REAL*8 BPRK(8),BPRA(8),BPRB(8),BPRT0(8),BPRT(8),BPRCOEF(8)
      REAL*8 ZERO,ONE,TWO, EPS
      REAL*8 VOLTC0,VOLTC1,VARTC,VOLVP0,VOLVP1,VARVP
      INTEGER CALIB,I, IARGC, ITER, ITERMAX
      CHARACTER(LEN=20) LINE
      CHARACTER(LEN=50) FICHIER
      CHARACTER(LEN=36) SCENARIO
      CHARACTER(LEN=12)  MODES (8)

      ZERO=0.0D0
      ONE=1.D0 
      TWO=2.D0
      CALIB=1 
      ITERMAX=50
      EPS=0.1D-9

      P = ZERO
      T = ZERO

C =============== LES DONNEES

      A=ONE
      DATA  MODES/ 'VP(CBD->CBD)', 'TC(CBD->CBD)',
     1             'VP(CBD->SUB)', 'TC(CBD->SUB)',
     1             'VP(SUB->CBD)', 'TC(SUB->CBD)',
     1             'VP(SUB->SUB)', 'TC(SUB->SUB)'/

C LES DONNEES SONT DANS UN FICHIER TEXTE EXTERNE
      IF(IARGC().NE.1) THEN 
         PRINT*, 'Un fichier de données doit être fournie'
         STOP
      ENDIF

      CALL GETARG(1, FICHIER)
C     PRINT*, FICHIER
      OPEN(UNIT=1,FILE=FICHIER)
      PRINT*, "Reading file: ", FICHIER
 6    READ(1,*,END=7) LINE
      IF(LINE(1:4).EQ.':vol') THEN
        READ(1,*) T(1:8)
      ELSEIF(LINE(1:4).EQ.':hou') THEN
        READ(1,*) T(13), T(15)
      ELSEIF(LINE(1:4).EQ.':ren') THEN
        READ(1,*) P(13), P(15)
      ELSEIF(LINE(1:4).EQ.':alp') THEN
        READ(1,*) ALPHA
      ELSEIF(LINE(1:4).EQ.':cou') THEN
        READ(1,*) COU(1:8)
      ELSEIF(LINE(1:4).EQ.':dis') THEN
        READ(1,*) DIS(1:8)
      ELSEIF(LINE(1:4).EQ.':vit') THEN
        READ(1,*) VIT(1:8)
      ELSEIF(LINE(1:4).EQ.':pri') THEN
        READ(1,*) P(20)
      ELSEIF(LINE(1:4).EQ.':val') THEN
        READ(1,*) VAL
      ELSEIF(LINE(1:4).EQ.':ela') THEN
        READ(1,*) S(09:12),S(14),S(16),S(17:19),S(21)
      ELSEIF(LINE(1:4).EQ.':pol') THEN
        READ(1,*) POL
      ELSEIF(LINE(1:4).EQ.':con') THEN
        READ(1,*) CONG
      ELSEIF(LINE(1:4).EQ.':occ') THEN
        READ(1,*) OCC
      ELSEIF(LINE(1:4).EQ.':abp') THEN
        READ(1,*) BPRA
      ELSEIF(LINE(1:4).EQ.':bbp') THEN
        READ(1,*) BPRB
      ELSEIF(LINE(1:4).EQ.':tbp') THEN
        READ(1,*) BPRCOEF
      ELSEIF(LINE(1:4).EQ.':sce') THEN
        READ(1,*) SCENARIO
      ENDIF
      GOTO 6 
 7    CLOSE(1)              

C LE VOLUME DE TRANSPORT AVEC LES FLUX INITIAUX
      VOL=ZERO
      DO 230, I=1,8
        VOL=VOL+T(I)/OCC(I)*CONG(I)
230   CONTINUE
C     PRINT*, VOL 

C CALIBRER LA VALEUR DE K POUR CHAQUE MODE
C
C The BPR function is
C 
C    T = T_0 ( 1 + a (V/K)^b )
C
C So:
C
C    K = V ( ( T/T_0 - 1 ) * (1/a) ) ^ (-1/b)
C
      DO 22, I=1,8
        BPRT(I)=DIS(I)/VIT(I)
        BPRT0(I)=BPRT(I)/BPRCOEF(I)
        BPRK(I)=VOL*((ONE/BPRA(I)*(BPRT(I)/BPRT0(I)-ONE)))
     2           **(-ONE/BPRB(I))
C       PRINT '(I5,7F7.2)',I,BPRA(I),BPRB(I),BPRT0(I),BPRT(I),
C    2      BPRCOEF(I), BPRK(I)
22    CONTINUE
C        DO 28, I=1,9
C           P(I)=COU(I)*DIS(I)+
C     2          BPR(VOL,BPRT0(I),BPRA(I),BPRB(I),BPRK(I))*VAL
C           PRINT '(I5,F7.3)', I, P(I)
C28      CONTINUE

C CALCUL DES COUTS GENERALISES
      DO 19, I=1,8
        P(I)=COU(I)*DIS(I)+DIS(I)/VIT(I)*VAL
C     PRINT '(I5,6F9.3)', I,VIT(I),DIS(I),COU(I),VAL,P(I),
C    +                    COU(I)*DIS(I)/P(I)
 19   CONTINUE


      ITER=0
C ============= LA CALIBRATION DES PARAMETRES (SI CALIB = 1)
C ============= CALCUL DES PRIX, DES VOLUMES

100   CONTINUE 

C DEBUT DE LA BOUCLE PRINCIPALE. CETTE BOUCLE CONTIENT QUATRE NIVEAUX
C (CORRESPONDANT A L ARBORESCENCE DU MODELE). ELLE EST PASSEE DEUX FOIS,
C LA PREMIERE POUR LA CALIBRATION (CALIB=1), ET LA DEUXIEME POUR LA
C SIMULATION (CALIB=0) DES CHANGEMENTS DANS LES PARAMETRES PRIX OU
C ELASTICITES

C == NIVEAU 1

      IF (CALIB.EQ.1) THEN 

         CALL CAL2(P(01),P(02),T(01),T(02),S(09),AA(01),AA(02))
         CALL CAL2(P(03),P(04),T(03),T(04),S(10),AA(03),AA(04))
         CALL CAL2(P(05),P(06),T(05),T(06),S(11),AA(05),AA(06))
         CALL CAL2(P(07),P(08),T(07),T(08),S(12),AA(07),AA(08))

         Y(09)=P(1)*T(1)+P(2)*T(2)
         Y(10)=P(3)*T(3)+P(4)*T(4)
         Y(11)=P(5)*T(5)+P(6)*T(6)
         Y(12)=P(7)*T(7)+P(8)*T(8) 

      ENDIF

      P(09)=PCES2(P(01),P(02),A,AA(01),AA(02),S(09))
      P(10)=PCES2(P(03),P(04),A,AA(03),AA(04),S(10))
      P(11)=PCES2(P(05),P(06),A,AA(05),AA(06),S(11))
      P(12)=PCES2(P(07),P(08),A,AA(07),AA(08),S(12))
      T(09)= CES2(T(01),T(02),A,AA(01),AA(02),S(09))
      T(10)= CES2(T(03),T(04),A,AA(03),AA(04),S(10))
      T(11)= CES2(T(05),T(06),A,AA(05),AA(06),S(11))
      T(12)= CES2(T(07),T(08),A,AA(07),AA(08),S(12)) 

C == NIVEAU 2

      IF (CALIB.EQ.1) THEN

         CALL CAL2(P(09),P(10),T(09),T(10),S(14),AA(09),AA(10))
         CALL CAL2(P(11),P(12),T(11),T(12),S(16),AA(11),AA(12))

         Y(14)=Y(09)+Y(10)
         Y(16)=Y(11)+Y(12) 

      ENDIF

      P(14)=PCES2(P(09),P(10),A,AA(09),AA(10),S(14))
      P(16)=PCES2(P(11),P(12),A,AA(11),AA(12),S(16))
      T(14)= CES2(T(09),T(10),A,AA(09),AA(10),S(14))
      T(16)= CES2(T(11),T(12),A,AA(11),AA(12),S(16))

C == NIVEAU 3

      IF (CALIB.EQ.1) THEN

         CALL CAL2(P(13),P(14),T(13),T(14),S(17),AA(13),AA(14))
         CALL CAL2(P(15),P(16),T(15),T(16),S(18),AA(15),AA(16))
         
         Y(17)=P(13)*T(13)+Y(14)
         Y(18)=P(15)*T(15)+Y(16) 

      ENDIF

      P(17)=PCES2(P(13),P(14),A,AA(13),AA(14),S(17))
      P(18)=PCES2(P(15),P(16),A,AA(15),AA(16),S(18))
      T(17)= CES2(T(13),T(14),A,AA(13),AA(14),S(17))
      T(18)= CES2(T(15),T(16),A,AA(15),AA(16),S(18))

C == NIVEAU 4

      IF (CALIB.EQ.1) THEN

         CALL CAL2(P(17),P(18),T(17),T(18),S(19),AA(17),AA(18))
         
         Y(19)=Y(17)+Y(18)
         T(20)=(ONE-ALPHA)/ALPHA*Y(19)

      ENDIF

      P(19)=PCES2(P(17),P(18),A,AA(17),AA(18),S(19))
      T(19)= CES2(T(17),T(18),A,AA(17),AA(18),S(19)) 

C == NIVEAU 5


      IF (CALIB.EQ.1) THEN

         CALL CAL2(P(19),P(20),T(19),T(20),S(21),AA(19),AA(20))
         
         Y(21)=Y(19)+P(20)*T(20)

C        SAVEGARDER LES FLUX DE TRANSPORT INITIAUX DANS UN TABLEAU
         DO I=1, 8
           T0(I)=T(I)
         ENDDO
         T0(13)=T(13)
         T0(15)=T(15)
         T0(20)=T(20) 
         T0(21)=CES2(T(19),T(20),A,AA(19),AA(20),S(21)) 

      ENDIF

      P(21)=PCES2(P(19),P(20),A,AA(19),AA(20),S(21))
      T(21)= CES2(T(19),T(20),A,AA(19),AA(20),S(21))

C  -----------------------------------------------------------
C  ============ VERIFICATION DU CALCUL ET VOLUMES DE TRANSPORT
C  -----------------------------------------------------------

C CALCUL DES UTILITES (INDICES QUANTITES) A TOUS LES NIVEAUX
      CALL QUANTITIES(AA,P,S,Y,A,POL,CONG,OCC,Q,TOTPOL,TOTCONG) 

      IF(CALIB.EQ.1) THEN
        CALIB=0 
C  ON REPREND LE FICHIER DES DONNEES POUR LIRE LES CHANGEMENTS DE PRIX
C  OU D ELASTICITE ET ON EXECUTE UNE DEUXIEME FOIS LA BOUCLE
        OPEN(UNIT=1,FILE=FICHIER)
 16     READ(1,*,END=17) LINE
        IF(LINE(1:4).EQ.':2di') THEN
          READ(1,*) DIS(1:8)
        ELSEIF(LINE(1:4).EQ.':2el') THEN
          READ(1,*) S(09:12),S(14),S(16),S(17:19),S(21)
        ELSEIF(LINE(1:4).EQ.':2vi') THEN
          READ(1,*) VIT(1:8)
        ELSEIF(LINE(1:4).EQ.':2re') THEN
          READ(1,*) P(13), P(15)
        ELSEIF(LINE(1:4).EQ.':2pr') THEN
          READ(1,*) P(20)
        ELSEIF(LINE(1:4).EQ.':2co') THEN
          READ(1,*) COU(1:8)
        ELSEIF(LINE(1:4).EQ.':2va') THEN
          READ(1,*) VAL
        ELSEIF(LINE(1:4).EQ.':2oc') THEN
          READ(1,*) OCC(1:8)
        ENDIF
        GOTO 16 
 17     CLOSE(1)              

        DO 18, I=1,8
           P(I)=COU(I)*DIS(I)+DIS(I)/VIT(I)*VAL
18      CONTINUE

        TOTCONG0=TOTCONG
        TOTPOL0=TOTPOL
        GOTO 100 
      ENDIF 

      DO 14, I=1,8
      VAR(I) =(Q( I)-T0( I))/T0( I)*100.D0
14    CONTINUE
      VAR(13)=(Q(13)-T0(13))/T0(13)*100.D0
      VAR(15)=(Q(15)-T0(15))/T0(15)*100.D0
      VAR(20)=(Q(20)-T0(20))/T0(20)*100.D0
      VAR(21)=(Q(21)-T0(21))/T0(21)*100.D0
C FIN DE LA BOUCLE PRINCIPALE 

C VOLUME TOTAL DE TRANSPORT AVEC LES NOUVEAUX FLUX
      VOL=ZERO
      DO 23, I=1,8
        VOL=VOL+Q(I)/OCC(I)*CONG(I)
23    CONTINUE 

      DO 28, I=1,8
         P(I)=COU(I)*DIS(I)+
     2        BPR(VOL,BPRT0(I),BPRA(I),BPRB(I),BPRK(I))*VAL
C        PRINT '(I5,3F9.3)', I, BPRT0(I), BPRK(I), P(I)
28    CONTINUE 

      VOLVP0=T0(1)+T0(3)+T0(5)+T0(7)
      VOLVP1=Q(1)+Q(3)+Q(5)+Q(7)
      VOLTC0=T0(2)+T0(4)+T0(6)+T0(8)
      VOLTC1=Q(2)+Q(4)+Q(6)+Q(8)
      VARTC = ( VOLTC1 - VOLTC0 ) / VOLTC0
      VARVP = ( VOLVP1 - VOLVP0 ) / VOLVP0

      ITER=ITER+1
      DELTAP=ZERO
      DO 60, I=1,8
         DELTAP=DELTAP+(Q(I)-T(I))*(Q(I)-T(I))
60    CONTINUE
      PRINT '(9X,A11,I4,A6,E10.2)', ' Itération', ITER,'---> ', DELTAP
      IF((DELTAP.GT.EPS).AND.(ITER.LT.ITERMAX)) THEN
        DO 61, I=1,8
          T(I)=Q(I)
61      CONTINUE
        T(20)=Q(20)
        GOTO 100
      ENDIF 

      IF(ITER.EQ.1) GOTO 100

      VOLVP0 = T0(1) + T0(3) + T0(5) + T0(7)
      VOLTC0 = T0(2) + T0(4) + T0(6) + T0(8)
      VOLVP1 = Q(1)  +  Q(3) +  Q(5) +  Q(7)
      VOLTC1 = Q(2)  +  Q(4) +  Q(6) +  Q(8)
      VARTC  =  ( VOLTC1 - VOLTC0 ) / VOLTC0 * 100.0D0
      VARVP  =  ( VOLVP1 - VOLVP0 ) / VOLVP0 * 100.0D0

C AFFICHE LE RESULTAT
      PRINT*
      PRINT '(10X,"===========================================")'
      PRINT '(10X,A36)', SCENARIO
      PRINT '(10X,"*******************************************")'
      PRINT '(12X,A9,3X,3A9)', 'Modes', 'Base', 'New','%'
      PRINT '(10X,"-------------------------------------------")'
      DO I=1,8,2
        PRINT '(12X,A12,3F9.2)', MODES(I), T0(I), Q(I), VAR(I)
      ENDDO
      PRINT '(10X,"-----")'
      PRINT '(12X,A12,3F9.2)','Total VP   ',VOLVP0,VOLVP1,VARVP
      PRINT '(12X,A12,3F9.2)','Congestion ',
     &        TOTCONG0,TOTCONG,(TOTCONG-TOTCONG0)/TOTCONG0*100
      PRINT '(12X,A12,3F9.2)','Pollution ',
     &        TOTPOL0,TOTPOL,(TOTPOL-TOTPOL0)/TOTPOL0*100
      PRINT '(10X,"...........................................")'
      DO I=2,8,2
        PRINT '(12X,A12,3F9.2)', MODES(I), T0(I), Q(I), VAR(I)
      ENDDO
      PRINT '(10X,"-----")'
      PRINT '(12X,A12,3F9.2)','Total TC   ',VOLTC0,VOLTC1,VARTC
      PRINT '(10X,"-------------------------------------------")'
      PRINT '(12X,A12,3A9)', 'Housing    ', 'Base', 'New','% '
      PRINT '(10X,"-------------------------------------------")'
      I=13
      PRINT '(12X,A12,3F9.2)', ' H_CBD    ', T0(I), Q(I), VAR(I)
      I=15
      PRINT '(12X,A12,3F9.2)', ' H_SUB    ', T0(I), Q(I), VAR(I)
      PRINT '(10X,"===========================================")'
      PRINT* 

C FIN DU PROGRAMME
      END 


C ----------------------------------------
C LES FONCTIONS ET SOUS-ROUTINES UTILISEES
C DANS LE PROGRAMME
C ----------------------------------------

      FUNCTION CES2(X,Y,A,A1,A2,S)
      IMPLICIT NONE
      REAL*8 CES2
      REAL*8 X,Y,A,A1,A2,S
      REAL*8 R,ONE

      ONE=1.D0

      R=(S-ONE)/S

      CES2= A * ( A1**(ONE-R) * X**R +  
     1            A2**(ONE-R) * Y**R )**(ONE/R)

      RETURN
      END 

      FUNCTION PCES2(P1,P2,A,A1,A2,S)
      IMPLICIT NONE
      REAL*8 PCES2
      REAL*8 P1,P2,A,A1,A2,S
      REAL*8 R,RP
      REAL*8 ONE

      ONE=1.D0 
      
      R=ONE-S 
      PCES2= A * ( A1 * P1**R +  
     1             A2 * P2**R )**(ONE/R)

      RETURN
      END 

      SUBROUTINE CAL2(P1,P2,T1,T2,S,A1,A2)
      IMPLICIT NONE
      REAL*8 P1,P2,T1,T2,S,A1,A2
      REAL*8 R

      R=1.D0/S

      A1=(P1*T1**R/(P1*T1**R+P2*T2**R))**S 
      A2=(P2*T2**R/(P1*T1**R+P2*T2**R))**S 

      RETURN
      END


      SUBROUTINE QUANTITIES(AA,P,S,Y,A,POL,CONG,OCC,Q,TOTPOL,TOTCONG)

      IMPLICIT NONE
      REAL*8 AA(21), P(21), S(21), Y(21), A, POL(8), CONG(8), OCC(8)
      REAL*8 Q(21), TOTPOL, TOTCONG
      REAL*8 CES2
      REAL*8 ZERO
      INTEGER I
    
      ZERO=0.D0

      Q(1)=Y(21)/P(21)* 
     1    AA(19)*(P(21)/P(19))**S(21)*
     1    AA(17)*(P(19)/P(17))**S(19)*
     2    AA(14)*(P(17)/P(14))**S(17)*
     3    AA(09)*(P(14)/P(09))**S(14)*
     4    AA(01)*(P(09)/P(01))**S(09)

      Q(2)=Y(21)/P(21)* 
     1    AA(19)*(P(21)/P(19))**S(21)*
     1    AA(17)*(P(19)/P(17))**S(19)*
     2    AA(14)*(P(17)/P(14))**S(17)*
     3    AA(09)*(P(14)/P(09))**S(14)*
     4    AA(02)*(P(09)/P(02))**S(09)

      Q(3)=Y(21)/P(21)* 
     1    AA(19)*(P(21)/P(19))**S(21)*
     1    AA(17)*(P(19)/P(17))**S(19)*
     2    AA(14)*(P(17)/P(14))**S(17)*
     3    AA(10)*(P(14)/P(10))**S(14)*
     4    AA(03)*(P(10)/P(03))**S(10)

      Q(4)=Y(21)/P(21)* 
     1    AA(19)*(P(21)/P(19))**S(21)*
     1    AA(17)*(P(19)/P(17))**S(19)*
     2    AA(14)*(P(17)/P(14))**S(17)*
     3    AA(10)*(P(14)/P(10))**S(14)*
     4    AA(04)*(P(10)/P(04))**S(10)

      Q(5)=Y(21)/P(21)* 
     1    AA(19)*(P(21)/P(19))**S(21)*
     1    AA(18)*(P(19)/P(18))**S(19)*
     2    AA(16)*(P(18)/P(16))**S(18)*
     3    AA(11)*(P(16)/P(11))**S(16)*
     4    AA(05)*(P(11)/P(05))**S(11)

      Q(6)=Y(21)/P(21)* 
     1    AA(19)*(P(21)/P(19))**S(21)*
     1    AA(18)*(P(19)/P(18))**S(19)*
     2    AA(16)*(P(18)/P(16))**S(18)*
     3    AA(11)*(P(16)/P(11))**S(16)*
     4    AA(06)*(P(11)/P(06))**S(11)

      Q(7)=Y(21)/P(21)* 
     1    AA(19)*(P(21)/P(19))**S(21)*
     1    AA(18)*(P(19)/P(18))**S(19)*
     2    AA(16)*(P(18)/P(16))**S(18)*
     3    AA(12)*(P(16)/P(12))**S(16)*
     4    AA(07)*(P(12)/P(07))**S(12)

      Q(8)=Y(21)/P(21)* 
     1    AA(19)*(P(21)/P(19))**S(21)*
     1    AA(18)*(P(19)/P(18))**S(19)*
     2    AA(16)*(P(18)/P(16))**S(18)*
     3    AA(12)*(P(16)/P(12))**S(16)*
     4    AA(08)*(P(12)/P(08))**S(12)

      Q(13)=Y(21)/P(21)* 
     1    AA(19)*(P(21)/P(19))**S(21)*
     1    AA(17)*(P(19)/P(17))**S(19)*
     2    AA(13)*(P(17)/P(13))**S(17) 

      Q(15)=Y(21)/P(21)* 
     1    AA(19)*(P(21)/P(19))**S(21)*
     1    AA(18)*(P(19)/P(18))**S(19)*
     2    AA(15)*(P(18)/P(15))**S(18) 
      
      Q(20)=Y(21)/P(21)*AA(20)*(P(21)/P(20))**S(21) 

      Q(09)=CES2(Q(01),Q(02),A,AA(01),AA(02),S(09))
      Q(10)=CES2(Q(03),Q(04),A,AA(03),AA(04),S(10))
      Q(11)=CES2(Q(05),Q(06),A,AA(05),AA(06),S(11))
      Q(12)=CES2(Q(07),Q(08),A,AA(07),AA(08),S(12))
      Q(14)=CES2(Q(09),Q(10),A,AA(09),AA(10),S(14))
      Q(16)=CES2(Q(11),Q(12),A,AA(11),AA(12),S(16))
      Q(17)=CES2(Q(13),Q(14),A,AA(13),AA(14),S(17))
      Q(18)=CES2(Q(15),Q(16),A,AA(15),AA(16),S(18))
      Q(19)=CES2(Q(17),Q(18),A,AA(17),AA(18),S(19))
      Q(21)=CES2(Q(19),Q(20),A,AA(19),AA(20),S(21))

C CALCUL DES EMISSIONS ET DU VOLUME EQUIVALENT VOITURE (CONGESTION)
      TOTPOL=ZERO
      TOTCONG=ZERO
      DO 10, I=1, 8
        TOTPOL=TOTPOL+POL(I)*Q(I)/OCC(I)
        TOTCONG=TOTCONG+CONG(I)*Q(I)/OCC(I)
10    CONTINUE 

      RETURN
      END


      FUNCTION BPR(V,T0,A,B,K)
      IMPLICIT NONE
      REAL*8 BPR, V, T0, A, B, K
      
      BPR=T0*(1.D0+A*(V/K)**B)

      RETURN
      END