**********************************************************************************
* CORRECTION 15/11/2017: USE OF THE SPLITTING IN THE NEW PAPER OF AMFR-W-METHODS
*     F(T,U)= G(U) + FXY(T,U) + FXX(T,U) + FYY(T,U)
*     G(U)=-RD*U, 
*  FXY(T,U)=A0*U+G0*EXP(-RF*T), FXX(T,U)=AS*U+G1*EXP(-RF*T), FYY(T,U)=AV*U+G2*EXP(-RF*T)
**********************************************************************************
**********************************************************************************
* 19/10/2017
* WE STUDY THE CONVERGENCE ORDER OF AMFW-, PDEW- AND AMFRW-METHODS 
* ON THE 2D- HESTON PROBLEM 
* IT HAS MIXED DERIVATIVES AND TIME DEPENDENT BCs 
**********************************************************************************

      PROGRAM AMFWMETHODS2D

      IMPLICIT REAL*8 (A-H,O-Z)
	CHARACTER LEY*99,FILENAME*99,FILINP*140

	PARAMETER (MD=200,ND=2*MD,NN=20,NS=5 ) 
        PARAMETER (ND3=3*ND,MD5=5*MD,NMX3=3*ND*MD)

	DIMENSION U(ND,0:MD-1),U0(ND,0:MD-1),SOL(ND,0:MD-1)
       DIMENSION LEY(NN),EGH2(NN),EGHU(NN),PH2(NN),PHU(NN)
       DIMENSION A(NS,NS),Q(NS,NS),C(NS),B(NS),HB(NS),RO(NS)
      
        DIMENSION SMESH(0:ND),VMESH(0:MD)
	DIMENSION AS(3,ND,0:MD-1),AV(5,0:MD-1)

      
      INTEGER IVMET(NN),LINSYST(NN)
      REAL*4 TIMENS, FINISH, START, TIME(NN) !COMPUTING CPU TIMES WITH CPU_TIME
            REAL*8 MU
            
        COMMON /COEFFICIENTS/A,Q,C,B,HB,RO,IQ,ISTAGE
        COMMON /PARAMS/XKAPPA,ETA,SIGMA,RHO,RD,RF,XK
	COMMON /BLOCK2/N,M
	COMMON /BLOCK3/NLINSYST
	COMMON /MESHX/SMESH
	COMMON /MESHY/VMESH

* SOLVING: HESTON PROBLEM -> U(t)= F(t,U) (ODEs)
* (s,v) \in (0,Send)x(0,Vend), t\in (0,TEND).
* 
* MOL APPROACH
* SPATIAL DISCRETIZATION USING STANDARD CENTRAL DIFFERENCES OVER A UNIFORM GRID: DX=1/(N+1), DY=1/(M+1)
* TIME INTEGRATION BY MEANS OF W-PDE-METHODS (WITH DIMENSIONAL SPLITTING)

      LEY(1)=' Heston Problem:'
      LEY(2)=' (s,y) \in (0,Send)x(0,Vend), t \in (0,TEND]'
      LEY(3)=' MOL approach. Central Differences in space '
      LEY(4)=' non-uniform spatial mesh'
      LEY(5)=' Dimensional Splitting'
      LEY(6)=' TAU=TEND/NSTEPS, Time integrations are performed for'
      LEY(7)=' NSTEP(J)=NSTEP0*2^(I-1), I=1,2,...,NJ'

	DO I=1,7
	WRITE(*,'(A)') LEY(I)
	ENDDO

* TIME INTERVAL	
	T0=0.D0
	!TEND depends on the case
	
* CHOOSE WHICH METHODS ARE TO BE USED
      NMET=2
      IVMET(1)=0    !1-stage- w-method 
      IVMET(2)=61 	

*DIFFERENT SPATIAL DIMENSIONS
       ! IEMES(1)=25
       ! IEMES(2)=50
       ! IEMES(3)=100
        

* CHOOSE WHICH AMF METHOD IS TO BE CONSIDERED (METOPT)
* METOPT=1: AMF-W METHODS
* METOPT=2: PDE-W METHODS
* METOPT=3: AMFR-W METHODS
* LOOP FOR SEVERAL AMF OPTIONS
               DO METOPT=1,3
               
	               IF (METOPT.EQ.1) THEN
                                             LEY(10)=' AMF-W METHOD'
                               ELSEIF (METOPT.EQ.2) THEN
                                             LEY(10)=' PDE-W METHOD'
                               ELSEIF (METOPT.EQ.3) THEN
                                             LEY(10)=' AMFR-W METHOD'
                               ENDIF   
               PRINT*, ' '        
	WRITE(*,'(A)') LEY(10)     
     
* LOOP FOR SEVERAL VALUES OF THE PARAMETERS OF THE PROBLEM 
      DO IK=0,2 !1,1 !2
      
*Parameters of the problem	
	CALL PROBPAR(IK,XKAPPA,ETA,SIGMA,RHO,RD,RF,XK,TEND)
	SEND=30*XK
	VEND=15.D0   
  
*LOOP FOR SEVERAL METHODS  

*********************    
      DO IMET=1,NMET
      
	METHOD=IVMET(IMET)
	
* OBTAIN COEFFICIENTS OF THE METHOD
  	CALL W_AMF_METHODS(METHOD,THETA,MU)
	

* NUMBER OF SPATIAL AND TEMPORAL INTEGRATIONS
       NUMINT=1 !2 !2 ! NUMBER OF SPATIAL INTEGRATIONS
	NJ=10    ! NUMBER OF TIME INTEGRATIONS
      NSTEP0=4 ! INITIAL NUMBER OF TIME STEPS

* FILE FOR RESULTS
        
	ITAG2=FLOOR(100*THETA)
	ITAG3=FLOOR(100*MU)
	
	IF (METOPT.EQ.1) THEN
               WRITE(filename,'(A,I1,A,I2,A,I2,A)')
     &"AMFW-C",IK,"-M",METHOD,"-TH0",ITAG2,"hes400.txt"
               ELSEIF (METOPT.EQ.2) THEN               
                     WRITE(filename,'(A,I1,A,I2,A,I2,A)')
     &"PDEW-C",IK,"-M",METHOD,"-TH0",ITAG2,"hes400.txt"
               ELSEIF (METOPT.EQ.3) THEN               
                     WRITE(filename,'(A,I1,A,I2,A,I2,A,I2,A)')
     &"AMFRW-C",IK,"-M",METHOD,"-TH0",ITAG2,"-MU0",ITAG3,"hes400.txt"
                    ENDIF               
      OPEN(1,FILE=filename)
      REWIND 1

* WRITING THE HEADING OF THE FILE FOR THE INTEGRATION STATISTICS
	DO L=1,7
            WRITE(1,'(A)') LEY(L)
      ENDDO
      WRITE(1,'(A)') '------------------------------------------------'
    
                    WRITE(1,'(A)') LEY(10)
      WRITE(1,'(A)') '------------------------------------------------'
      
*	WRITE(1,'(A),/') 

	WRITE(1,50) NSTEP0,NJ,TEND
50	FORMAT('NSTEP0=',I2,'	NJ=',I2,'	TEND=',F5.3)

	WRITE(1,75) XKAPPA,ETA,SIGMA,RHO,RD,RF,XK 
75	FORMAT('kappa=',F8.3,'	eta=',F8.3,' sigma=',F8.3,
     &'	 rho=',F8.3,' r_d=',F8.3,' r_f=',F8.3,' K=',F8.3)
     
     	IF (METOPT.EQ.1) THEN
     	WRITE(1,80) THETA
80	FORMAT('THETA=',F8.4)
               ELSEIF (METOPT.EQ.2) THEN  
     	WRITE(1,81) THETA
81	FORMAT('THETA=',F8.4)
               ELSEIF  (METOPT.EQ.3) THEN  
	WRITE(1,82) THETA, MU
82	FORMAT('THETA=',F8.4,'	MU=',F8.4)
               ENDIF

*	WRITE(1,'(A),/') 
	WRITE(1,'(A)') '---------------------------'
	WRITE(1,'(1X,A,I3,4X,A,I3)') 'METHOD=',METHOD,' CASE=',IK
	WRITE(1,'(A)') '---------------------------'
*	WRITE(1,'(A),/') 

* LOOP FOR INCREASING SPATIAL DIMENSIONS
	DO INN=1,NUMINT

      PRINT*,'SPATIAL-INTEGRATION=',INN
		M=200  !IEMES(INN)    !2**(INN+1)-1
		N=2*M

*	WRITE(1,'(A),/') 
      WRITE(1,'(A)') '-------------------------------------------------'
      
    	WRITE(1,100) N,M
100	FORMAT('N=',I4,'    M=',I4)


*****NON-UNIFORM SPATIAL MESH
*Computation of the non-uniform mesh and the corresponding increments 
	CALL SNET(TEND,SEND,SMESH)
	CALL VNET(VEND,VMESH)
*******

*In this problem, these two matrices are constant in the whole program. Then, they are
*only computed once and share with the COMMON/MATRICES/ 
*They are used in the subroutine RHSHES that compute the derivative function
       CALL ASMAT(AS)
       CALL AVMAT(AV)

* COMPUTING THE INITIAL VALUE 
	CALL INIHES(U0)

* COMPUTING  THE EXACT SOLUTION AT END-POINT
        CALL HESTEXACT(M,IK,FILINP)

	PRINT*,'NOMBRE DEL FICHERO=',FILINP

	OPEN(2,FILE=FILINP)
	DO JJ=0,M-1
	DO II=1,N
	 READ(2,'(D24.17)') SOL(II,JJ)
*	PRINT*,'LEE BIEN I=',II,', J=',JJ,' SOL= ',SOL(II,JJ)
	ENDDO  
	ENDDO  
	CLOSE(2)
*	PRINT*,'SALIO DEL FICHERO BIEN'

*********
	WRITE(1,150) 'NSTEPS','EG2','EGU','ORD2','ORDU','TIME','#LINSYST'
150   FORMAT (A6,T13,A3,T26,A3,T38,A4,T46,A4,T55,A4,T65,A7)

* PERFORMING ALL THE TIME INTEGRATIONS
      DO I=1,NJ

	NSTEPS=NSTEP0*2**(I-1)
	U=U0
               NLINSYST=0
	CALL CPU_TIME(START)
        CALL AMFWMETH(METOPT,THETA,MU,METHOD,NSTEPS,TEND,U,
     &      AS,AV)   
        CALL CPU_TIME(FINISH)
      TIMENS=FINISH-START

* COMPUTING THE GLOBAL ERRORS
	U=SOL-U
	CALL GLOBAL_ERRORS(U,GE2,GEU)

* COMPUTING THE ORDER OF CONVERGENCE P AS POWER OF TAU
	ORD2=0.D0
	ORDU=0.D0

	IF (I.GT.1) THEN
	      ORD2=DLOG(AE2/(1.D-30+GE2))/DLOG(2.D0)
            ORDU=DLOG(AEU/(1.D-30+GEU))/DLOG(2.D0)		  
      ENDIF
      AE2=GE2
      AEU=GEU

* COMPUTING THE ORDER OF CONVERGENCE P IN PDE SENSE (TIME STEPSIZE=GRIDSIZE)
	IF (I.EQ.INN) THEN
		EGH2(I)=GE2
		PH2(I)=DLOG(GE2)/DLOG(2.D0)
		EGHU(I)=GEU
		PHU(I)=DLOG(GEU)/DLOG(2.D0)
		TIME(I)=FINISH-START
		LINSYST(I)=NLINSYST
	ENDIF

*	WRITING STATISTICS IN A FILE
	WRITE(1,200) NSTEP0*2**(I-1),GE2,GEU,ORD2,ORDU,TIMENS,NLINSYST
200   FORMAT(1X,I5,3X,2(D10.4,3X),2(F6.3,2X),D10.4,3X,I6)

	PRINT*,'NSTEPS=',NSTEPS,'DONE!'

      END DO ! END I LOOP
	END DO ! END INN LOOP

*	WRITE(1,'(A),/') 
*	WRITE(1,'(A),/') 
	WRITE(1,'(A)') '----------------------------------'
	WRITE(1,'(1X,A)') 'ORDER PDE (TAU=H, NSTEPS=N+1=M+1)'
	WRITE(1,'(A)') '----------------------------------'
*	WRITE(1,'(A),/') 

	WRITE(1,250) 'NSTEPS','EG2','EGU','PH2','PHU','TIME','#LINSYST'
250   FORMAT (A6,T13,A3,T26,A3,T38,A4,T46,A4,T55,A4,T65,A7)
	WRITE(1,300) NSTEP0,EGH2(1),EGHU(1),0.D0,0.D0,TIME(1),LINSYST(1)
	DO I=2,NUMINT
	WRITE(1,300) NSTEP0*2**(I-1),EGH2(I),EGHU(I),
     &            -PH2(I)+PH2(I-1),-PHU(I)+PHU(I-1),TIME(I),LINSYST(I)
300   FORMAT(1X,I5,3X,2(D10.4,3X),2(F6.3,2X),D10.4,3X,I6)
	ENDDO

	CLOSE(1)

      END DO ! END IMET LOOP
      END DO ! END IK LOOP
      END DO ! END METOPT LOOP

	STOP
	END

***************************************************************************************

      SUBROUTINE AMFWMETH(METOPT,THETA,MU,METHOD,NSTEPS,TEND,U,
     &      AS,AV)
      
* IN THIS ROUTINE AMF-W, PDE-W, AND AMFR-W METHODS WITH FIXED TIME STEPSIZE ARE APPLIED FOR THE TIME INTEGRATION OF
* LINEAR PARABOLIC PDES WITH POSSIBLY MIXED DERIVATIVES. NON-HOMOGENEOUS DIRICHLET BCs ARE ALLOWED.

      IMPLICIT REAL*8 (A-H,O-Z)

      PARAMETER (NS=5)

      DIMENSION A(NS,NS),Q(NS,NS),C(NS),B(NS),HB(NS),RO(NS)
	DIMENSION XJAC(3,N,0:M-1),YJAC(5,0:M-1)
      DIMENSION U(N,0:M-1),WORK(3,N)
      DIMENSION AS(3,N,0:M-1),AV(5,0:M-1)
       REAL*8 MU

        COMMON /COEFFICIENTS/A,Q,C,B,HB,RO,IQ,ISTAGE
	COMMON /PARAMS/XKAPPA,ETA,SIGMA,RHO,RD,RF,XK
	COMMON /BLOCK2/N,M
	
	
      TAU=TEND/NSTEPS

* JACOBIAN SPLITTING FOR X-DIRECTION
        DO J=0,M-1
	
        DO I=1,N
              WORK(1,I)=-THETA*TAU*AS(1,I,J) ! UPPER DIAGONAL
	      WORK(2,I)=1.D0-THETA*TAU*AS(2,I,J) ! MAIN DIAGONAL
	      WORK(3,I)=-THETA*TAU*AS(3,I,J) ! LOWER DIAGONAL
	ENDDO
		CALL TRIDIAGONAL(N,WORK)
        DO I=1,N
              XJAC(1,I,J)=WORK(1,I) ! UPPER DIAGONAL
	      XJAC(2,I,J)=WORK(2,I) ! MAIN DIAGONAL
	      XJAC(3,I,J)=WORK(3,I) ! LOWER DIAGONAL
	ENDDO
	
	ENDDO

* JACOBIAN SPLITTING FOR Y-DIRECTION
	DO J=0,M-1
	      YJAC(1,J)= -THETA*TAU*AV(1,J) ! SECOND-UPPER DIAGONAL
	      YJAC(2,J)= -THETA*TAU*AV(2,J)! UPPER DIAGONAL
	      YJAC(3,J)=1.D0-THETA*TAU*AV(3,J) ! MAIN DIAGONAL
	      YJAC(4,J)= -THETA*TAU*AV(4,J)! LOWER DIAGONAL
	      YJAC(5,J)= -THETA*TAU*AV(5,J)! SECOND-LOWER DIAGONAL
	ENDDO

* CROUT-LU DECOMPOSITION OF THE TRIDIAGONAL MATRICES
	CALL PENTADIAG(M,YJAC)

* THE TIME INTEGRATION WITH THE CHOSEN AMF VERSION:

            IF (METOPT.EQ.1) THEN
                     DO I=1,NSTEPS
		T=(I-1)*TAU
             CALL AMFW_ONESTEP(THETA,T,TAU,U,XJAC,YJAC,
     &        AS,AV)
                     ENDDO         
                              
            ELSEIF (METOPT.EQ.2) THEN
                     DO I=1,NSTEPS
		T=(I-1)*TAU	
             CALL PDEW_ONESTEP(THETA,T,TAU,U,XJAC,YJAC,
     &        AS,AV)
	      ENDDO
	      
               ELSEIF (METOPT.EQ.3) THEN
                      DO I=1,NSTEPS
		T=(I-1)*TAU	
            CALL AMFRW_ONESTEP(THETA,MU,T,TAU,U,XJAC,YJAC,
     &        AS,AV)
	      ENDDO               
               ENDIF

	RETURN
	END

***************************************************************************************

      SUBROUTINE SOLVE_DIREC_X(XJAC,W0)
* THIS ROUTINE SOLVE M TRIDIAGONAL LINEAR SYSTEMS OF DIMENSION N (X-DIRECTION) 
* XJAC(J)*X= W0(*,I), I=1,2,...,N, FOR J=0,1,...,M-1
* THE SOLUTION IS STORED IN THE SAME VECTOR W0.
* INPUT: XJAC(3,N,0:M-1),N,M,W0(N,M)
* OUTPUT: W0(N,M)
* WX IS AN AUXILIARY VECTOR

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION XJAC(3,N,0:M-1),W0(N,0:M-1),WX(N),WORK(3,N)
	COMMON /BLOCK2/N,M

	DO J=0,M-1

		DO K=1,N
			WX(K)=W0(K,J) !AUXILIAR VECTOR
			
			DO II=1,3   !TRIDIAGONAL MATRIX
                        WORK(II,K)=XJAC(II,K,J) 
	      		ENDDO 
		ENDDO
        
		CALL TRIDSOLVE(N,WORK,WX)

		DO K=1,N
			W0(K,J)=WX(K)
		ENDDO
	ENDDO

	RETURN
	END

***************************************************************************************

      SUBROUTINE SOLVE_DIREC_Y(YJAC,W0)
* THIS ROUTINE SOLVE N PENTADIAGONAL LINEAR SYSTEMS OF DIMENSION M (Y-DIRECTION) 
* YJAC*X= W0(K,*), K=1,2,...,N
* THE SOLUTION IS STORED IN THE SAME VECTOR W0.
* INPUT: YJAC(5,0:M-1),N,M,W0(N,0:M-1)
* OUTPUT: W0(N,M)
* WX IS AN AUXILIARY VECTOR

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION YJAC(5,0:M-1),W0(N,0:M-1),WX(0:M-1)
	COMMON /BLOCK2/N,M

	DO K=1,N
		DO J=0,M-1
			WX(J)=W0(K,J) !AUXILIAR VECTOR
		ENDDO

		CALL PENTASOLVE(M,YJAC,WX)

		DO J=0,M-1
			W0(K,J)=WX(J)
		ENDDO
	ENDDO

	RETURN
	END


***************************************************************************************

      SUBROUTINE FDERXX(T,AS,V,F)
*	THIS ROUTINE CALCULATES THE SPLIT-XX FUNCTION OF THE SEMI-DISCRETIZED ODE SYSTEM AT ALL MESH-POINTS
*	F=FX(T,V)= AS*V  WITHOUT THE TERM G1*EXP(-RF*T)
*	W IS AN ARRAY INCLUDING THE BOUNDARY VALUES FOR V AT TIME T.

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION V(N,0:M-1),F(N,0:M-1)
	DIMENSION AS(3,N,0:M-1)

	COMMON /BLOCK2/N,M
	COMMON /PARAMS/XKAPPA,ETA,SIGMA,RHO,RD,RF,XK
 

	DO J=0,M-1

                F(1,J)=AS(2,1,J)*V(1,J)+AS(1,2,J)*V(2,J)
	
		DO I=2,N-1
		F(I,J)=AS(3,I-1,J)*V(I-1,J)+AS(2,I,J)*V(I,J)+
     &            AS(1,I+1,J)*V(I+1,J)
		ENDDO

		F(N,J)=AS(3,N-1,J)*V(N-1,J)+AS(2,N,J)*V(N,J)
	
	ENDDO
        
        RETURN
	END

***************************************************************************************

      SUBROUTINE FDERYY(T,AV,V,F)
*	THIS ROUTINE CALCULATES THE SPLIT-YY FUNCTION OF THE SEMI-DISCRETIZED ODE SYSTEM AT ALL MESH-POINTS
*	F=F2(T,V)=AV*V WITHOUT THE TERM G2*EXP(-RF*T)
*	W IS AN ARRAY INCLUDING THE BOUNDARY VALUES FOR V AT TIME T.

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION V(N,0:M-1),F(N,0:M-1)
	DIMENSION AV(5,0:M-1)

	COMMON /BLOCK2/N,M
	COMMON /PARAMS/XKAPPA,ETA,SIGMA,RHO,RD,RF,XK
 
	DO I=1,N

*J=0
        	F(I,0)=AV(3,0)*V(I,0)+AV(2,1)*V(I,1)+AV(1,2)*V(I,2) 		
*J=1
		F(I,1)=AV(4,0)*V(I,0)+AV(3,1)*V(I,1)+
     &            AV(2,2)*V(I,2)+AV(1,3)*V(I,3) 		
	
		DO J=2,M-3
		F(I,J)=AV(5,J-2)*V(I,J-2)+AV(4,J-1)*V(I,J-1)+
     &            AV(3,J)*V(I,J)+AV(2,J+1)*V(I,J+1)+AV(1,J+2)*V(I,J+2) 				
		ENDDO
*J=M-2		
		F(I,M-2)=AV(5,M-4)*V(I,M-4)+AV(4,M-3)*V(I,M-3)+
     &            AV(3,M-2)*V(I,M-2)+AV(2,M-1)*V(I,M-1) 		
*J=M-1
		F(I,M-1)=AV(5,M-3)*V(I,M-3)+AV(4,M-2)*V(I,M-2)+
     &            AV(3,M-1)*V(I,M-1)				
	ENDDO

	RETURN
	END

***************************************************************************************

      SUBROUTINE FDERXY(T,V,F)
*	THIS ROUTINE CALCULATES THE SPLIT-XY FUNCTION OF THE SEMI-DISCRETIZED ODE SYSTEM AT ALL MESH-POINTS
*	F=G(V)+FXY(T,V)= (A0-RD*I)*V -WITHOUT THE TERM G0*EXP(-RF*T)
*	W IS AN ARRAY INCLUDING THE BOUNDARY VALUES FOR V AT TIME T.
        
        PARAMETER(ND=200,MD=100)
       IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION V(N,0:M-1),F(N,0:M-1)
	DIMENSION SMESH(0:ND),VMESH(0:MD),SDIF(N),VDIF(M)
	DIMENSION BES(-1:1,N-1), BEV(-1:1,M-1)

        COMMON /BLOCK2/N,M
	COMMON /PARAMS/XKAPPA,ETA,SIGMA,RHO,RD,RF,XK
	COMMON /MESHX/SMESH
	COMMON /MESHY/VMESH

		
	CALL DIFFER(N,SMESH,SDIF)
	CALL DIFFER(M,VMESH,VDIF)

*MIXED DERIVATIVE
*Here the components F0(I,0) are not necessary since they are multiplied by v_0=0.
*Besides, when I=N, that is, the point (S,v_j) the mixed derivative is 0 due to the 
*Neumann condition (du/ds)(S,v,t)=e^{-rf t}
*Then, we only computes the components for i=1:N-1, j=1:M-1

*Computation of \beta_{i,-1,0,1}, i=1,N-1
	DO I=1,N-1

	CTE1=SDIF(I)+SDIF(I+1)

	!these are \beta_{i,+1}, i=1:N-1
	BES(1,I)= SDIF(I)/(SDIF(I+1)*CTE1)

	!these are \beta_{i,0}, i=1:N-1
	BES(0,I)=(SDIF(I+1) -SDIF(I))/(SDIF(I)*SDIF(I+1))

	!these are \beta_{i,-1}, i=1:N-1
	IF(I.EQ.1) THEN
	BES(-1,I)=0.D0
	ELSE
	BES(-1,I)= -BES(0,I)-BES(1,I)
	ENDIF

	ENDDO

	DO J=1,M-1

	CTE1=VDIF(J)+VDIF(J+1)

	!these are \hat{\beta}_{j,-1}, j=1:M-1
	BEV(-1,J)=-VDIF(J+1)/(VDIF(J)*CTE1)

	!these are \hat{\beta}_{j,0}, j=1:M-1
	BEV(0,J) =(VDIF(J+1)-VDIF(J))/(VDIF(J)*VDIF(J+1))

	!these are \hat{\beta}_{j,1}, j=1:M-1
	BEV(1,J)= -BEV(0,J)- BEV(-1,J)

	ENDDO

	F=0.D0

*J=0
*	DO I=1,N
*	F(I,0)=0.D0
*	ENDDO

*I=N
*	DO J=1,M-1
*	F(N,J)=0.D0
*	ENDDO

*Otherwise
	CTE1=RHO*SIGMA

*J=1,...,M-2
	DO J=1,M-2

*	 I=1
		I=1
		F(I,J)=0.D0

		DO KI=0,1
		DO KJ=-1,1
		F(I,J)= F(I,J) + BES(KI,I)*BEV(KJ,J)*V(I+KI,J+KJ)
		ENDDO
		ENDDO
		F(I,J)= CTE1*SMESH(I)*VMESH(J)*F(I,J)

*	 I=2,...,N-1
		DO I=2,N-1
		F(I,J)=0.D0

		DO KI=-1,1
		DO KJ=-1,1
		F(I,J)= F(I,J) + BES(KI,I)*BEV(KJ,J)*V(I+KI,J+KJ)
		ENDDO
		ENDDO
		F(I,J)= CTE1*SMESH(I)*VMESH(J)*F(I,J)

		ENDDO

	ENDDO  !END OF J=1,M-2

*J=M-1
	J=M-1

*	 I=1
		I=1
		F(I,J)=0.D0

		DO KI=0,1
		DO KJ=-1,0
		F(I,J)= F(I,J) + BES(KI,I)*BEV(KJ,J)*V(I+KI,J+KJ)
		ENDDO
		ENDDO
		F(I,J)= CTE1*SMESH(I)*VMESH(J)*F(I,J)

*	 I=2,...,N-1
		DO I=2,N-1
		F(I,J)=0.D0

		DO KI=-1,1
		DO KJ=-1,0
		F(I,J)= F(I,J) + BES(KI,I)*BEV(KJ,J)*V(I+KI,J+KJ)
		ENDDO
		ENDDO
		F(I,J)= CTE1*SMESH(I)*VMESH(J)*F(I,J)

		ENDDO
		
        !ADDING THE G(V)=-RD*V TERM	
		F=F-RD*V
		
	RETURN

	END	

***************************************************************************************

*COMPUTATION OF THE AS-MATRIX FOR THE LINEAR SYSTEMS WITH RESPECTO TO S
	SUBROUTINE ASMAT(AS)

        PARAMETER(MD=200,ND=100)
	IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION AS(3,M,0:N-1),WORK1(3,M),WORK2(3,M)
	DIMENSION SMESH(0:MD),SDIF(M),VMESH(0:ND)
	
	DIMENSION DELTA(-1:1,M-1)
	DIMENSION BETA(-1:1,M-1)

        COMMON /PARAMS/XKAPPA,ETA,SIGMA,RHO,RD,RF,XK
	COMMON /BLOCK2/M,N
	COMMON /MESHX/SMESH
	COMMON /MESHY/VMESH

	
	CALL DIFFER(M,SMESH,SDIF)

*work(1,2:M)  =DIAG+1
*work(2,1:M)  =MAIN DIAG 
*work(3,1:M-1)=DIAG-1

	DELTA(-1,1)=0.D0
	DELTA(0,1) =-2.D0/(SDIF(1)*SDIF(2))
	DELTA(1,1) = 2.D0/(SDIF(2)*(SDIF(1)+SDIF(2)))

	BETA(-1,1)=0.D0
	BETA(0,1) =(SDIF(2)-SDIF(1))/(SDIF(1)*SDIF(2))
	BETA(1,1) = SDIF(1)/(SDIF(2)*(SDIF(1)+SDIF(2)))

	DO I=2,M-1

	FAC1=SDIF(I)+SDIF(I+1)
	
	!these are \delta_{i,-1}, i=2:M-1
	DELTA(-1,I)= 2.D0/(SDIF(I)*FAC1)

	!these are \delta_{i,+1}, i=2:M-1
	DELTA(1,I)= 2.D0/(SDIF(I+1)*FAC1)

	!these are \delta_{i,0}, i=2:M-1
*	DELTA(I,0)= -2.D0/(SDIF(I)*SDIF(I+1))
	DELTA(0,I)= -DELTA(-1,I)-DELTA(1,I) !As Karels Matlab program 
	
	!these are \beta_{i,-1}, i=2:M-1
	BETA(-1,I)= -SDIF(I+1)/(SDIF(I)*FAC1)

	!these are \beta_{i,+1}, i=2:M-1
	BETA(1,I)= SDIF(I)/(SDIF(I+1)*FAC1)

	!these are \beta_{i,0}, i=2:M-1
*	BETA(I,0)=(SDIF(I+1) -SDIF(I))/(SDIF(I)*SDIF(I+1))
	BETA(0,I)=-BETA(-1,I)-BETA(1,I)  !As Karels Matlab program

	ENDDO


*DIFFUSION TERM	
	WORK1=0.D0
		
	WORK1(2,1)=SMESH(1)**2 *DELTA(0,1)
	WORK1(3,1)=SMESH(2)**2 *DELTA(-1,2)
	
	DO I=2,M-2
	WORK1(1,I)=SMESH(I-1)**2 *DELTA(1,I-1)
	WORK1(2,I)=SMESH(I)**2   *DELTA(0,I)
	WORK1(3,I)=SMESH(I+1)**2 *DELTA(-1,I+1)
	ENDDO

	FAC1=2.D0*(SMESH(M)/SDIF(M))**2

	WORK1(1,M-1)=SMESH(M-2)**2 *DELTA(1,M-2)
	WORK1(2,M-1)=SMESH(M-1)**2 *DELTA(0,M-1)
	WORK1(3,M-1)=FAC1
	
	WORK1(1,M)=SMESH(M-1)**2 *DELTA(1,M-1)
	WORK1(2,M)= -FAC1


*ADVECTION TERM
	WORK2=0.D0

	WORK2(2,1)=SMESH(1)*BETA(0,1)
	WORK2(3,1)=SMESH(2)*BETA(-1,2)

	DO I=2,M-2
	WORK2(1,I)=SMESH(I-1)*BETA(1,I-1)
	WORK2(2,I)=SMESH(I)*BETA(0,I)
	WORK2(3,I)=SMESH(I+1)*BETA(-1,I+1)
	ENDDO
	
	WORK2(1,M-1)=SMESH(M-2)*BETA(1,M-2)
	WORK2(2,M-1)=SMESH(M-1)*BETA(0,M-1)
	
	WORK2(1,M)=SMESH(M-1)*BETA(1,M-1)

	
*THE WHOLE MATRIX


*	IF(J.GE.1) THEN !FOR J=0 V0=0
        AS=0.D0

	DO I=1,M
	AS(1,I,0)= (RD-RF)*WORK2(1,I)
	AS(2,I,0)= (RD-RF)*WORK2(2,I) !!!-(RD/2.D0)  !MAIN DIAGONAL
	AS(3,I,0)= (RD-RF)*WORK2(3,I) 
	ENDDO
	
	DO J=1,N-1
	DO I=1,M
	AS(1,I,J)= VMESH(J)*WORK1(1,I)/2.D0+(RD-RF)*WORK2(1,I)
	AS(2,I,J)= VMESH(J)*WORK1(2,I)/2.D0+(RD-RF)*WORK2(2,I) !!!- 
****     &           (RD/2.D0)  !MAIN DIAGONAL
	AS(3,I,J)= VMESH(J)*WORK1(3,I)/2.D0+(RD-RF)*WORK2(3,I) 
	ENDDO
	ENDDO 

	RETURN

	END

*******************************************************************************

*COMPUTATION OF THE AV-MATRIX FOR THE LINEAR SYSTEMS WITH RESPECTO TO V

	SUBROUTINE AVMAT(AV)

        PARAMETER(ND=100)
        
	IMPLICIT REAL*8 (A-H,O-Z) 
	INTEGER JIND

	DIMENSION AV(5,0:N-1),WORK1(5,0:N-1),WORK2(5,0:N-1)
	DIMENSION VMESH(0:ND),VDIF(N)
	DIMENSION DELTA(-1:1,N-1)
	DIMENSION ABETA(-1:1,N-1)

	COMMON /PARAMS/XKAPPA,ETA,SIGMA,RHO,RD,RF,XK
	COMMON /BLOCK2/M,N
	COMMON /MESHY/VMESH

	CALL DIFFER(N,VMESH,VDIF)

*work(1,2:N-1)= DIAG+2 
*work(2,1:N-1)= DIAG+1 
*work(3,0:N-1)= MAIN DIAG 
*work(4,0:N-2)= DIAG-1
*work(5,0:N-3)= DIAG-2
 

*DIFFUSION TERM (TRIDIAGONAL)

	DO J=1,N-2
	
	FAC1=VDIF(J)+VDIF(J+1)
	!these are \hat{\delta}_{j,1}, j=1:N-2
	DELTA(1,J)=2.D0/(VDIF(J+1)*FAC1)
	!these are \hat{\delta}_{j,-1}, j=1:N-2
	DELTA(-1,J)=2.D0/(VDIF(J)*FAC1)
	!these are \hat{\delta}_{j,0}, j=1:N-2
	DELTA(0,J)=-DELTA(1,J)-DELTA(-1,J)

	ENDDO

	DELTA(-1,N-1)= 2.D0/(VDIF(N-1)*(VDIF(N-1)+VDIF(N)))
	DELTA(0,N-1) =-2.D0/(VDIF(N-1)*VDIF(N))
	DELTA(1,N-1) = 0.D0


        WORK1=0.D0
*J=0        
	WORK1(4,0)= VMESH(1)*DELTA(-1,1)
*J=1
	WORK1(3,1)= VMESH(1)*DELTA(0,1)
	WORK1(4,1)= VMESH(2)*DELTA(-1,2)

	DO J=2,N-2
	WORK1(2,J)=VMESH(J-1)*DELTA(1,J-1)
	WORK1(3,J)=VMESH(J)*DELTA(0,J)
	WORK1(4,J)=VMESH(J+1)*DELTA(-1,J+1)
	ENDDO
*J=N-1
	WORK1(2,N-1)= VMESH(N-2)*DELTA(1,N-2)
	WORK1(3,N-1)= VMESH(N-1)*DELTA(0,N-1)


*ADVECTION TERM (PENTADIAGONAL)

*Detection of the index JIND sucha that V_{JIND)<1 AND V_{JIND+1}>=1
 	JIND=0
	DO J=1,N
	IF(VMESH(J).GE.1.D0) THEN
	JIND=J-1
	EXIT  !exit of the do loop
	ENDIF
	ENDDO

*Coefficients of the finite difference formulas

*j=0
	FAC1=VDIF(1)+VDIF(2)
	GAM00=(-2.D0*VDIF(1)-VDIF(2))/(VDIF(1)*FAC1)
	GAM01=FAC1/(VDIF(1)*VDIF(2))
	GAM02=-VDIF(1)/(VDIF(2)*FAC1)


*j=1,JIND  (v_{jind}<1)
	DO J=1,JIND
	
	FAC1=VDIF(J)+VDIF(J+1)
	!these are \hat{\beta}_{j,-1}, j=1:JIND
	ABETA(-1,J)=-VDIF(J+1)/(VDIF(J)*FAC1)

	!these are \hat{\beta}_{j,1}, j=1:JIND
	ABETA(1,J)=VDIF(J)/(VDIF(J+1)*FAC1)

	!these are \hat{\beta}_{j,0}, j=1:JIND
	ABETA(0,J)= -ABETA(-1,J)-ABETA(1,J)
	
	ENDDO


*J=JIND-1,N-1

	DO J=JIND+1,N-1

	FAC1=VDIF(J-1)+VDIF(J)

	!these are \hat{\alpha}_{j,0}, j=JIND+1:N-1
	ABETA(0,J)=(VDIF(J-1)+2.D0*VDIF(J))/(VDIF(J)*FAC1)		

	!these are \hat{\alpha}_{j,-1}, j=JIND+1:N-1
	ABETA(-1,J)=(-VDIF(J-1)-VDIF(J))/(VDIF(J-1)*VDIF(J))		

	!these are \hat{\alpha}_{j,-2}, j=JIND+1:N-1
	ABETA(1,J)= -ABETA(0,J)- ABETA(-1,J)	
	ENDDO


        WORK2=0.D0

	WORK2(3,0)=ETA*GAM00
	WORK2(4,0)=(ETA-VMESH(1))*ABETA(-1,1)

	WORK2(2,1)= ETA*GAM01
	WORK2(3,1)= (ETA-VMESH(1))*ABETA(0,1)
	WORK2(4,1)= (ETA-VMESH(2))*ABETA(-1,2)

	WORK2(1,2)= ETA*GAM02   !only value for the first row

	DO J=2,JIND-2
	WORK2(2,J)= (ETA-VMESH(J-1))*ABETA(1,J-1)
	WORK2(3,J)= (ETA-VMESH(J))*ABETA(0,J)
	WORK2(4,J)= (ETA-VMESH(J+1))*ABETA(-1,J+1)
	ENDDO

	DO J=JIND-1,JIND+1
	WORK2(2,J)= (ETA-VMESH(J-1))*ABETA(1,J-1)
	WORK2(3,J)= (ETA-VMESH(J))*ABETA(0,J)
	WORK2(4,J)= (ETA-VMESH(J+1))*ABETA(-1,J+1)
	WORK2(5,J)= (ETA-VMESH(J+2))*ABETA(1,J+2)  !ABETA(1,J+2)=Alpha(-2,J+2)
        ENDDO 

	DO J=JIND+2,N-3
	WORK2(3,J)= (ETA-VMESH(J))*ABETA(0,J)
	WORK2(4,J)= (ETA-VMESH(J+1))*ABETA(-1,J+1)
	WORK2(5,J)= (ETA-VMESH(J+2))*ABETA(1,J+2) !ABETA(1,J+2)=Alpha(-2,J+2)
	ENDDO

	WORK2(3,N-2)= (ETA-VMESH(N-2))*ABETA(0,N-2)
	WORK2(4,N-2)= (ETA-VMESH(N-1))*ABETA(-1,N-1)

	WORK2(3,N-1)= (ETA-VMESH(N-1))*ABETA(0,N-1)

*THE WHOLE MATRIX
	FAC1=SIGMA**2/2.D0

	DO J=0,N-1
	AV(1,J)= FAC1*WORK1(1,J)+XKAPPA*WORK2(1,J)
	AV(2,J)= FAC1*WORK1(2,J)+XKAPPA*WORK2(2,J)
	AV(3,J)= FAC1*WORK1(3,J)+XKAPPA*WORK2(3,J) !!!!-(RD/2.D0)  !MAIN DIAGONAL
	AV(4,J)= FAC1*WORK1(4,J)+XKAPPA*WORK2(4,J)
	AV(5,J)= FAC1*WORK1(5,J)+XKAPPA*WORK2(5,J)
	ENDDO

	RETURN

	END

***************************************************************************************
      SUBROUTINE GES(G0,G1,G2)  
 
       PARAMETER(ND=200,MD=100)
         
      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION SMESH(0:ND),VMESH(0:MD)
        DIMENSION G0(N,0:M-1),G1(N,0:M-1),G2(N,0:M-1)

        COMMON /PARAMS/XKAPPA,ETA,SIGMA,RHO,RD,RF,XK
	COMMON /BLOCK2/N,M
	COMMON /MESHX/SMESH
	COMMON /MESHY/VMESH

        G0=0.D0
        G1=0.D0
        G2=0.D0
     
*Non-homogeneous term of the mixed derivative
       CTE0=SIGMA*RHO*VMESH(M-1)
       CTE0=CTE0*(VMESH(M-1)-VMESH(M-2))
       CTE0=CTE0/((VMESH(M)-VMESH(M-1))*(VMESH(M)-VMESH(M-2)))

	DO I=1,N-1
       G0(I,M-1)=CTE0*SMESH(I)
       ENDDO

*Non-homogeneous term on the s-direction
	DO J=0,M-1
	G1(N,J)=VMESH(J)*(SMESH(N)**2)/(SMESH(N)-SMESH(N-1)) +
     &        (RD-RF)*SMESH(N)
	ENDDO

*Non-homogeneous term on the v-direction
       CTE2=(SIGMA**2)*VMESH(M-1)
       CTE2=CTE2/((VMESH(M)-VMESH(M-1))*(VMESH(M)-VMESH(M-2))) 

        DO I=1,N
        G2(I,M-1)=CTE2*SMESH(I)
        ENDDO
		
        RETURN

	END


***************************************************************************************
*COMPUTATION OF THE NON-UNIFORM SPATIAL MESH-IN THE S-DIRECTION 
	SUBROUTINE SNET(TEND,SEND,SMESH)

        PARAMETER(MD=200)     	
        IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION SMESH(0:MD)	
	
	COMMON /PARAMS/XKAPPA,ETA,SIGMA,RHO,RD,RF,XK
	COMMON /BLOCK2/M,N
	
	C=XK/1.D1

	PL=MAX(0.5D0,DEXP(-0.1D0*TEND))*XK
	PR=MIN(1.5D0,DEXP(+0.1D0*TEND))*XK

	XMIN= DASINH(-PL/C)
	XINT= (PR-PL)/C
	XMAX= XINT+ DASINH((SEND-PR)/C)

	XINCRE=(XMAX-XMIN)/M

	DO I=0,M
		XI=XMIN+I*XINCRE

		IF (XI.LE.0.D0) THEN
			SMESH(I)= PL+ C* DSINH(XI)
		ELSEIF ((XI.GT.0.D0).AND.(XI.LT.XINT)) THEN
			SMESH(I)= PL+ C* XI
		ELSE
			SMESH(I)= PR+ C* DSINH(XI -XINT)
		ENDIF

		SMESH(I)=MAX(SMESH(I),0.D0)

	ENDDO

	RETURN
	END

***************************************************************************************

*COMPUTATION OF THE NON-UNIFORM SPATIAL MESH-IN THE V-DIRECTION 
	SUBROUTINE VNET(VEND,VMESH)
	
        PARAMETER(ND=100)         	
        IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION VMESH(0:ND)
	
	COMMON /PARAMS/XKAPPA,ETA,SIGMA,RHO,RD,RF,XK
	COMMON /BLOCK2/M,N
	
	D=VEND/5.D2

	XINCRE=DASINH(VEND/D)/N
	
	DO J=0,N
		VMESH(J)= D*DSINH(J*XINCRE)
		VMESH(J)=MAX(VMESH(J),0.D0)
	ENDDO

	RETURN
	END


***************************************************************************************

*COMPUTATION OF THE INITIAL CONDITIONS
	SUBROUTINE INIHES(U0)
	
        PARAMETER(MD=200)
        IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION U0(M,0:N-1),SMESH(0:MD)

	COMMON /PARAMS/XKAPPA,ETA,SIGMA,RHO,RD,RF,XK
	COMMON /BLOCK2/M,N
	COMMON /MESHX/SMESH
	
	IND=0
	DO I=1,M
	IF(SMESH(I) .GT. XK) THEN
	IND=I
	EXIT
	ENDIF
	ENDDO

	IF(DABS(SMESH(IND)-XK).GT.DABS(SMESH(IND-1)-XK)) IND=IND-1
	

	SL=(SMESH(IND-1)+SMESH(IND))/2.D0
	SR=(SMESH(IND+1)+SMESH(IND))/2.D0
		

	DO J=0,N-1
	DO I=1,M

	IF(I.LT.IND) THEN
	U0(I,J)=0.D0
	ELSEIF(I.EQ.IND) THEN
	U0(I,J)= 0.5D0*((SR-XK)**2)/(SR-SL)
	ELSE
	U0(I,J)=SMESH(I)-XK
	ENDIF

	ENDDO
	ENDDO

	RETURN

	END

***************************************************************************************
*THE INVERSE OF THE HYPERBOLIC SINE IS NOT AVAILABLE IN FORTRAN 90

	DOUBLE PRECISION FUNCTION DASINH(X)

	DOUBLE PRECISION X

	DASINH= DLOG(X+DSQRT(X*X+1.D0))

	RETURN

	END

***************************************************************************************
        SUBROUTINE DIFFER(N,V,DIFV)

        DOUBLE PRECISION V(0:N),DIFV(N)
        
        DO K=1,N
        DIFV(K)=V(K)-V(K-1)
        ENDDO
        
        RETURN
        END
        
***************************************************************************************

	SUBROUTINE PROBPAR(ICASE,XKAPPA,ETA,SIGMA,RHO,RD,RF,XK,TEND)

	REAL*8 XKAPPA,ETA,SIGMA,RHO,RD,RF,XK,TEND
                        
	IF (ICASE.EQ.0) THEN
	XKAPPA=1.5D0
	ETA=0.2D-1
	SIGMA=0.62D0
	RHO=-0.67D0
	RD=1.D-2
	RF=2.D-2
	XK=1.D2
	TEND=1.D0

	ELSEIF (ICASE.EQ.1) THEN
	XKAPPA=1.5D0
	ETA=0.4D-1
	SIGMA=0.3D0
	RHO=-0.9D0
	RD=0.25D-1
	RF=0.D0
	XK=1.D2
	TEND=1.D0


	ELSEIF (ICASE.EQ.2) THEN
	XKAPPA=3.D0
	ETA=0.12D0
	SIGMA=0.4D-1
	RHO= 0.6D0
	RD=0.1D-1
	RF=0.4D-1
	XK=1.D2
	TEND=1.D0

	ELSEIF (ICASE.EQ.3) THEN
	XKAPPA=0.6067D0
	ETA=0.0707D0
	SIGMA=0.2928D-1
	RHO= -0.7571D0
	RD=0.3D-1
	RF=0.0D0
	XK=1.D2
	TEND=3.D0


	ELSEIF (ICASE.EQ.4) THEN
	XKAPPA=2.5D0
	ETA=0.06D0
	SIGMA=0.5D0
	RHO= -0.1D0
	RD=0.0507D0
	RF=0.0469D0
	XK=1.D2
	TEND=0.25D0

	ENDIF

	RETURN
	END

***************************************************************************************
	
	SUBROUTINE HESTEXACT(N,ICASE,FILINP)

	CHARACTER*140 FILINP
*	CHARACTER DIREC*18 !THE DIRECTORY WHERE THE REFERENCE SOLUTION IS PLACED	

*	DIREC='../Ref_sol_buenas/'
*        DIREC='./Ref_sol_buenas/'
		IF(ICASE.EQ.1) THEN
              IF(N.EQ.25) THEN
                     FILINP='REFNEW-25-C1.TXT' !'karel_25_c1.txt' !'HES_N25C1.TXT'              
                ELSEIF(N.EQ.50) THEN
                     FILINP='REFNEW-50-C1.TXT'
              ELSEIF(N.EQ.100) THEN
                     FILINP='REFNEW-100-C1.TXT'
               ELSEIF(N.EQ.200) THEN
                     FILINP='REFSOL-200-C1-15.txt'
              ENDIF

       
       ELSEIF(ICASE.EQ.2) THEN
              IF(N.EQ.25) THEN
                    FILINP='REFNEW-25-C2.TXT' !'karel_25_c2.txt' !'HES_N25C1.TXT'                     
              ELSEIF(N.EQ.50) THEN
                     FILINP='REFNEW-50-C2.TXT'
              ELSEIF(N.EQ.100) THEN
                     FILINP='REFNEW-100-C2.TXT'
              ELSEIF(N.EQ.200) THEN
                     FILINP='REFSOL-200-C2-15.txt'
             ENDIF
              
*       ELSEIF(ICASE.EQ.3) THEN
*              IF(N.EQ.25) THEN
*                     FILINP=DIREC//'REFNEW-25-C3-15.TXT' 
*              ELSEIF(N.EQ.50) THEN
*                     FILINP=DIREC//'REFNEW-50-C3-15.TXT'
*              ELSEIF(N.EQ.100) THEN
*                     FILINP=DIREC//'REFNEW-100-C3-15.TXT'
*              ENDIF
       
*       ELSEIF(ICASE.EQ.4) THEN
*              IF(N.EQ.25) THEN
*                     FILINP=DIREC//'REFNEW-25-C4-15.TXT' !'karel_25_c1.txt' !'HES_N25C1.TXT'
*              ELSEIF(N.EQ.50) THEN
*                     FILINP=DIREC//'REFNEW-50-C4-15.TXT'
*              ELSEIF(N.EQ.100) THEN
*                     FILINP=DIREC//'REFNEW-100-C4-15.TXT'
*              ENDIF
              
        ELSEIF(ICASE.EQ.0) THEN
              IF(N.EQ.25) THEN
                     FILINP='REFSOL-25-C0-15.txt'
              ELSEIF(N.EQ.50) THEN
                     FILINP='REFSOL-50-C0-15.txt'
              ELSEIF(N.EQ.100) THEN
                     FILINP='REFSOL-100-C0-15.txt'
               ELSEIF(N.EQ.200) THEN
                     FILINP='REFSOL-200-C0-15.txt'
              ENDIF
       
	ENDIF 
	
	RETURN
	END

***************************************************************************************

	SUBROUTINE GLOBAL_ERRORS(U,GE2,GEU)
* THIS ROUTINE COMPUTES THE L2-NORM AND THE L_{\infty}-NORM OF U

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION U(N,0:M-1)
	COMMON /BLOCK2/N,M

      GEU=0.D0
	GE2=0.D0

	NM=N*M

	DO K=0,M-1
		DO J=1,N
			GEU=MAX(GEU,DABS(U(J,K)))
			GE2=GE2+U(J,K)**2/NM
		ENDDO
	ENDDO

      GE2=DSQRT(GE2)

	RETURN
	END

***************************************************************************************
  
      SUBROUTINE TRIDIAGONAL(N,A)
* THIS ROUTINE MAKES THE CROUT-LU DECOMPOSITION OF THE TRIDIAGONAL MATRIX A
* INPUT: A(3,N) MATRIX ENTERED IN BANDED FORM
*              A(1,K), K=2,3,...,N STORES THE UPPER DIAGONAL
*              A(2,K), K=1,2,...,N STORES THE MAIN DIAGONAL
*              A(3,K), K=1,2,...,N-1 STORES THE LOWER DIAGONAL
* OUTPUT: A(3,N) MATRIX  IN BANDED FORM
*              A(1,K), K=2,3,...,N STORES THE UPPER DIAGONAL OF CROUT LU
*              A(2,K), K=1,2,...,N STORES THE MAIN DIAGONAL
*              A(3,K), K=1,2,...,N-1 STORES THE LOWER DIAGONAL (IT IS UNCHANGED)

	IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION A(3,N)

	DO K=2,N
		A(1,K)=A(1,K)/A(2,K-1)
		A(2,K)=A(2,K)-A(3,K-1)*A(1,K)
	ENDDO

	RETURN
	END

***************************************************************************************

      SUBROUTINE TRIDSOLVE(N,A,B)
* THIS ROUTINE SOLVE LINEAR SYSTEMS A*X=B, WITH A BEING A TRIDIAGONAL MATRIX FACTORED IN TRIDIAGONAL()
* THE SOLUTION X IS STORED IN B AS OUTPUT
* INPUT: A(3,N) FROM ROUTINE TRIDIAGONAL()
*        B=B(N)
* OUTPUT: B=B(N) IS THE SOLUTION X OF A*X=B

	IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION A(3,N),B(N)

      B(1)=B(1)/A(2,1)

	DO K=2,N
		B(K)=(B(K)-A(3,K-1)*B(K-1))/A(2,K)
	ENDDO

	DO K=N-1,1,-1
		B(K)=B(K)-A(1,K+1)*B(K+1)
	ENDDO

	RETURN
	END

***************************************************************************************
	
	SUBROUTINE PENTADIAG(N,A)
* THIS ROUTINE MAKES THE CROUT-LU DECOMPOSITION OF THE TRIDIAGONAL MATRIX A
* INPUT: A(5,N) MATRIX ENTERED IN BANDED FORM
*              A(1,K), K=3,4,...,N STORES THE SECOND-UPPER DIAGONAL
*              A(2,K), K=2,3,...,N STORES THE UPPER DIAGONAL
*              A(3,K), K=1,2,...,N STORES THE MAIN DIAGONAL
*              A(4,K), K=1,2,...,N-1 STORES THE LOWER DIAGONAL
*              A(5,K), K=1,2,...,N-2 STORES THE SECOND-LOWER DIAGONAL
* OUTPUT: A(5,N) MATRIX  OF CROUT LU IN BANDED FORM
*              A(1,K), K=3,4,...,N STORES THE SECOND-UPPER DIAGONAL
*              A(2,K), K=2,3,...,N STORES THE UPPER DIAGONAL
*              A(3,K), K=1,2,...,N STORES THE MAIN DIAGONAL
*              A(4,K), K=1,2,...,N-1 STORES THE LOWER DIAGONAL
*              A(5,K), K=1,2,...,N-2 STORES THE SECOND-LOWER DIAGONAL
*                                    (IT IS UNCHANGED)

	IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION A(5,N)

*Crouts algorithm FOR FIVE DIAGONALS- WITHOUT PIVOTATION

        DO K=2,N
        
        IF (K.GT.2)  A(1,K)=A(1,K)/A(3,K-2)
   
        IF (K.EQ.2) THEN   
           A(2,K)=A(2,K)/A(3,K-1)
           A(3,K)=A(3,K)-A(4,K-1)*A(2,K)
        ELSE  
           A(2,K)=(A(2,K)-A(4,K-2)*A(1,K))/A(3,K-1)    
           A(3,K)=A(3,K)-A(5,K-2)*A(1,K)-A(4,K-1)*A(2,K)
        ENDIF 
    
        IF (K.LT.N) A(4,K)=A(4,K)-A(5,K-1)*A(2,K)
            
        ENDDO

	RETURN
	
	END



***************************************************************************************

        SUBROUTINE PENTASOLVE(N,A,B)
* THIS ROUTINE SOLVE LINEAR SYSTEMS A*X=B, WITH A BEING A PENTADIAGONAL MATRIX FACTORED IN TRIDIAGONAL()
* THE SOLUTION X IS STORED IN B AS OUTPUT
* INPUT: A(5,N) FROM ROUTINE PENTADIAGONAL()
*        B=B(N)
* OUTPUT: B=B(N) IS THE SOLUTION X OF A*X=B

	IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION A(5,N),B(N)
	
        B(1)=B(1)/A(3,1)
        B(2)=(B(2)-A(4,1)*B(1))/A(3,2)
        
        DO K=3,N
        B(K)=(B(K)-A(4,K-1)*B(K-1)-A(5,K-2)*B(K-2))/A(3,K)
        ENDDO

        B(N-1)=B(N-1)-A(2,N)*B(N)

        DO K=N-2,1,-1
		B(K)=B(K)-A(2,K+1)*B(K+1)-A(1,K+2)*B(K+2)
	ENDDO

	RETURN

	END


