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