123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523 |
- 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
|