PROGRAM REGIONTES
      IMPLICIT NONE
      INTEGER N,NOBS,NARGS,NLQ
      PARAMETER(N=1000)
      REAL*8 EMPL(N),EMPG(N),PROD(N),PRODL(N)
      REAL*8 TES(N,N),COEF(N,N),LQ(N)
      REAL*8 MATLQ(N,N),LEONTIEF(N,N),MULTIPLIERS(N)
      REAL*8 DELTA,BMULT,BEMPL
      INTEGER ORDM(N),DEBUG,LMAT
      CHARACTER*25 FILENAME
      CHARACTER*20 REGIONAME
      CHARACTER*40 LABELS(N)
      CHARACTER*5 CODES(N) 
      CHARACTER*5 METHOD
      CHARACTER*6 SORT
      CHARACTER*10 ARGS
      CHARACTER*10 FMT
      REAL*8 EMPL09(N),EMPL30(9),EMPL1(N),EMPL2(N)
      REAL*8 PRODL09(N),PRODL30(N),PRODL1(N),PRODL2(N)
      REAL*8 F09(N),F30(N),F1(N),F2(N) 
      INTEGER I,J 
********************************
      REAL*8 F(N),Y(N),WORK(N),Z(N),RCOND,LL(N,N),L(N,N)
      REAL*8 DET, NEWPRODL(N)
      INTEGER IPVT(N)

      NARGS=IARGC()
      IF(NARGS.LT.1) THEN
        PRINT*,'YOU SHOULD PROVIDE THE INPUT FILE NAME'
        STOP
      ENDIF 
      CALL GETARG(1,FILENAME)

      NOBS=N
      CALL INPUTDATA(FILENAME,NOBS,TES,EMPL,EMPG,PROD,
     &       REGIONAME,DEBUG,METHOD,DELTA,LABELS,CODES,SORT) 
C     PRINT*,'step 1'
C     WRITE(*,'(36F9.2)') TES(25,1:NOBS)

C The command line options overwrite those 
C given in the input file

      IF(NARGS.GT.1) THEN
      DO 10 I=1,NARGS
         CALL GETARG(I,ARGS)
         IF(ARGS.EQ.'-L=YES') THEN
            LMAT=1
         ELSEIF(ARGS.EQ.'-L=NO') THEN
            LMAT=0
         ELSEIF(ARGS.EQ.'-d=YES') THEN
            DEBUG=1
         ELSEIF(ARGS.EQ.'-d=NO') THEN
            DEBUG=0
         ELSEIF(ARGS.EQ.'-s=LQ') THEN
            SORT='LQ'
         ELSEIF(ARGS.EQ.'-s=MULT') THEN
            SORT='METHOD'
         ELSEIF(ARGS.EQ.'-s=NONE') THEN
            SORT='NONE'
         ELSEIF(ARGS.EQ.'-m=SLQ') THEN
            METHOD='SLQ'
         ELSEIF(ARGS.EQ.'-m=CILQ') THEN
            METHOD='CILQ'
         ELSEIF(ARGS.EQ.'-m=RLQ') THEN
            METHOD='RLQ'
         ELSEIF(ARGS.EQ.'-m=CILQ') THEN
            METHOD='CILQ'
         ELSEIF(ARGS.EQ.'-m=FLQ') THEN
            METHOD='FLQ'
         ELSEIF(ARGS.EQ.'-m=NAT') THEN
            METHOD='NAT'
         ENDIF
 10   CONTINUE
      ENDIF

      CALL MATCOEF(TES,PROD,COEF,NOBS,N,DEBUG)
      CALL LOCQ(METHOD,EMPL,EMPG,NOBS,N,DELTA,MATLQ,LQ,DEBUG)
C     PRINT*,'step 3'
      CALL REGIONAL(N,NOBS,COEF,MATLQ,LEONTIEF,MULTIPLIERS,DEBUG)
C     PRINT*,COEF(25,1:NOBS)
C     PRINT*,'step 4'
      NLQ=0
      DO I=1,NOBS
         IF(LQ(I).GT.1.D0) THEN
                 NLQ=NLQ+1
                 BEMPL=BEMPL+EMPL(I)
         ENDIF
      ENDDO
      BMULT=SUM(EMPL(1:NOBS))/BEMPL
      PRINT*, '======================================'
      PRINT*, 'REGION/CITY NAME: ', REGIONAME
      PRINT*, '--------------------------------------'
      PRINT '(X,A31,I3)',   'NUMBER OF SECTORS WITH LQ > 1 :', NLQ
      PRINT*, '--------------------------------------'
      PRINT '(X,A31,F7.3)', 'BASE MULTIPLIER               :', BMULT
      PRINT*, '--------------------------------------'
      IF(SORT.EQ.'LQ') THEN
         CALL SORTRX(NOBS,LQ,ORDM)
      ELSEIF(SORT.EQ.'METHOD')THEN
         CALL SORTRX(NOBS,MULTIPLIERS,ORDM)
      ELSE
         DO I=1,NOBS
            ORDM(I)=NOBS-I+1
         ENDDO
      ENDIF
      PRINT'(18X,A6)',METHOD                 
      PRINT*,' CODES   LQ       MULTIP.   DESCR.'
      PRINT*,'======================================'
      DO I=1,NOBS
         J=ORDM(NOBS-I+1)
         PRINT '((5X)(A4)(f6.3)(3X)(F6.3)(5X)(A35))', 
     &           CODES(J), LQ(J),MULTIPLIERS(J),LABELS(J)
      ENDDO 
      PRINT*,'======================================'
      PRINT* 
      IF(LMAT.EQ.1) THEN
        WRITE(FMT,'("("I0,"F9.3)")') NOBS
        PRINT*,'LEONTIEF MATRIX (I-A)^(-1)'
        DO I=1,NOBS
          WRITE(*,FMT),(LEONTIEF(I,J),J=1,NOBS)
        ENDDO
      ENDIF
      PRINT* 

************************************************
************************************************
* This part is added for computations to the
* case of IDF.
************************************************
************************************************
C     DO I=1,NOBS
C        PRODL(I)=PROD(I)/EMPG(I)*EMPL(I)
*        PRINT '(A2,3F15.3)', CODES(I),PRODL(I),PROD(I),PROD(I)/EMPG(I)
C     ENDDO
* LL = (I-A) and
* LEONTIEF = (I-A)^-1
* PRODL : the initial production
* EMPL : the initial employment
* NEWPRODL : the new production
* NEWEMPL : the new employment
* 
C     LL=LEONTIEF
C     CALL DGECO(LL,N,NOBS,IPVT,RCOND,Z)
C     CALL DGEDI(LL,N,NOBS,IPVT,DET,WORK,01)

* Now LL is the inverse of LEONTIEF
      CALL DGEMV('N',NOBS,NOBS,1.0D0,LL,N,PRODL,1,0.0D0,Y,1)
* Y contains the final consumption
 
*     PRINT*, 'Emploi final à l échelle locale (M€)'
*     DO I=1,NOBS
*         PRINT '(F12.3)', Y(I)
*     ENDDO

* here a test of whether the computation of 
* the inverse matrix are correct 
*     CALL DGEMM('N','N',NOBS,NOBS,NOBS,1.D0,LL,N,LEONTIEF,N,0.D0,L,N)

* the transport sector has an increase of 10% in final demand
*     PRINT '(F12.3)',Y(20)
C     Y(20)=1.1D0*Y(20)
*     PRINT '(F12.3)',Y(20)

* the induced change in production is

C     CALL DGEMV('N',NOBS,NOBS,1.D0,LEONTIEF,N,Y,1,0.0D0,NEWPRODL,1)

C     DO I=1,NOBS
C        PRINT '(A2,",",F12.3,",",F12.3)', CODES(I),LQ(I),
C    &   (NEWPRODL(I)-PRODL(I))/(PROD(I)/EMPG(I))
C     ENDDO 

*     PRINT*, LEONTIEF(1,1:NOBS) 
*     DO I=1,5
*        PRINT '(5F12.3)', LL(I,1:5)
*     ENDDO
*     PRINT*
*     DO I=1,5
*        PRINT '(5F12.3)', LEONTIEF(I,1:5)
*     ENDDO
*     PRINT*
*     DO I=1,5
*        PRINT '(5F12.3)', L(I,1:5)
*     ENDDO
*     Y(20)=Y(20)*1.1D0

      END