**********************************************************************************
* CORRECTION 15/11/2017: USE OF THE SPLITTING IN THE NEW PAPER OF AMFR-W-METHODS
**********************************************************************************
***************************************************************************************
* 19/10/2017
***************************************************************************************

	SUBROUTINE AMFW_ONESTEP(THETA,T,TAU,U,XJAC,YJAC,
     &        AS,AV)
* ONE STEP WITH AN AMF-W-METHOD WITH COEFFICIENTS GIVEN IN COMMON/COEFFICIENTS IS PERFORMED
*
* INPUT:      THETA: DIAGONAL COEFFICIENT OF THE  W-METHOD
*	   T,TAU: CURRENT T-POINT AND STEP-SIZE
*        U(N,M): NUMERICAL SOLUTION AT T
*        XJAC, YJAC,...: ARE TRIDIAGONAL MATRICES TO SOLVE LINEAR SYSTEMS
*
* OUTPUT: U(N,M): NUMERICAL SOLUTION AT POINT T+TAU

        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),UA(N,0:M-1)
      DIMENSION FXX(N,0:M-1),FYY(N,0:M-1),FXY(N,0:M-1)
      DIMENSION SA(N,0:M-1),S1(N,0:M-1),S2(N,0:M-1)
      DIMENSION S3(N,0:M-1),S4(N,0:M-1),S5(N,0:M-1)

      DIMENSION AS(3,N,0:M-1),AV(5,0:M-1)
      DIMENSION G0(N,0:M-1),G1(N,0:M-1),G2(N,0:M-1)

       	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


        SA=0.D0
        S1=0.D0
        S2=0.D0
        S3=0.D0
        S4=0.D0
        S5=0.D0
      
*Constant vectors for the boundary conditions
        CALL GES(G0,G1,G2)

        DEXPO=-RF*DEXP(-RF*T)

      FAC=THETA*TAU
      COEFD=1.D0+FAC*RD
      
* FIRST STAGE
      UA=U
      TA=T
      
      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      SA=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)
      SA=SA+FAC*TAU*RO(1)*G1*DEXPO  
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(1)*G2*DEXPO 
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S1=SA

      NLINSYST=NLINSYST+2

	IF (ISTAGE.GE.2) THEN
* SECOND STAGE
      UA=U+A(2,1)*S1
      TA=T+C(2)*TAU

      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      SA=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)+Q(2,1)*S1
      SA=SA+FAC*TAU*RO(2)*G1*DEXPO 
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(2)*G2*DEXPO 
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S2=SA

      NLINSYST=NLINSYST+2
	ENDIF

	IF (ISTAGE.GE.3) THEN
* THIRD STAGE
      UA=U+A(3,1)*S1+A(3,2)*S2
      TA=T+C(3)*TAU

      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      SA=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)+
     &        Q(3,1)*S1+Q(3,2)*S2
      SA=SA+FAC*TAU*RO(3)*G1*DEXPO       
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(3)*G2*DEXPO      
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S3=SA
      
      NLINSYST=NLINSYST+2
	ENDIF

	IF (ISTAGE.GE.4) THEN
* FOURTH STAGE
      UA=U+A(4,1)*S1+A(4,2)*S2+A(4,3)*S3
      TA=T+C(4)*TAU

      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      SA=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)+
     &        Q(4,1)*S1+Q(4,2)*S2+Q(4,3)*S3
      SA=SA+FAC*TAU*RO(4)*G1*DEXPO     
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(4)*G2*DEXPO      
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S4=SA


      NLINSYST=NLINSYST+2
	ENDIF

	IF (ISTAGE.GE.5) THEN
* FIFTH STAGE
      UA=U+A(5,1)*S1+A(5,2)*S2+A(5,3)*S3+A(5,4)*S4
      TA=T+C(5)*TAU

      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      SA=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)+
     &    Q(5,1)*S1+Q(5,2)*S2+Q(5,3)*S3+Q(5,4)*S4
      SA=SA+FAC*TAU*RO(5)*G1*DEXPO     
         CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(5)*G2*DEXPO     
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S5=SA

      NLINSYST=NLINSYST+2
	ENDIF

* ADVANCING SOLUTION
      U=U+B(1)*S1+B(2)*S2+B(3)*S3+B(4)*S4+B(5)*S5

      RETURN
	END

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

	SUBROUTINE PDEW_ONESTEP(THETA,T,TAU,U,XJAC,YJAC,
     &        AS,AV)
* ONE STEP WITH A PDE-W-METHOD WITH COEFFICIENTS GIVEN IN COMMON/COEFFICIENTS IS PERFORMED
*
* INPUT:      THETA: DIAGONAL COEFFICIENT OF THE  W-METHOD
*	   T,TAU: CURRENT T-POINT AND STEP-SIZE
*        U(N,M): NUMERICAL SOLUTION AT T
*        XJAC, YJAC,...: ARE TRIDIAGONAL MATRICES TO SOLVE LINEAR SYSTEMS
*
* OUTPUT: U(N,M): NUMERICAL SOLUTION AT POINT T+TAU

      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),UA(N,0:M-1)
      DIMENSION FXX(N,0:M-1),FYY(N,0:M-1),FXY(N,0:M-1)
      DIMENSION SA(N,0:M-1),FSA(N,0:M-1),S1(N,0:M-1),S2(N,0:M-1)
      DIMENSION S3(N,0:M-1),S4(N,0:M-1),S5(N,0:M-1)

      DIMENSION AS(3,N,0:M-1),AV(5,0:M-1)
      DIMENSION G0(N,0:M-1),G1(N,0:M-1),G2(N,0:M-1)

        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

        SA=0.D0
        S1=0.D0
        S2=0.D0
        S3=0.D0
        S4=0.D0
        S5=0.D0
        
*Constant vectors for the boundary conditions
        CALL GES(G0,G1,G2)
       
        DEXPO=-RF*DEXP(-RF*T)

      FAC=THETA*TAU

* FIRST STAGE
      UA=U
      TA=T
     
      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      S1=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)
      SA=S1+FAC*TAU*RO(1)*G1*DEXPO   
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(1)*G2*DEXPO   
      CALL SOLVE_DIREC_Y(YJAC,SA)

      CALL FDERXY(TA,SA,FSA)
      SA=S1+FAC*FSA+FAC*TAU*RO(1)*G0*DEXPO 
      SA=SA+FAC*TAU*RO(1)*G1*DEXPO   
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(1)*G2*DEXPO   
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S1=SA

      NLINSYST=NLINSYST+4

	IF (ISTAGE.GE.2) THEN
* SECOND STAGE
      UA=U+A(2,1)*S1
      TA=T+C(2)*TAU
      
      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      S2=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)+Q(2,1)*S1
      SA=S2+FAC*TAU*RO(2)*G1*DEXPO  
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(2)*G2*DEXPO  
      CALL SOLVE_DIREC_Y(YJAC,SA)

      CALL FDERXY(T,SA,FSA)
      SA=S2+FAC*FSA+FAC*TAU*RO(2)*G0*DEXPO  
      SA=SA+FAC*TAU*RO(2)*G1*DEXPO  
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(2)*G2*DEXPO  
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S2=SA

      NLINSYST=NLINSYST+4
	ENDIF

	IF (ISTAGE.GE.3) THEN
* THIRD STAGE
      UA=U+A(3,1)*S1+A(3,2)*S2
      TA=T+C(3)*TAU

      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      S3=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)+
     &        Q(3,1)*S1+Q(3,2)*S2
      SA=S3+FAC*TAU*RO(3)*G1*DEXPO 
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(3)*G2*DEXPO 
      CALL SOLVE_DIREC_Y(YJAC,SA)

      CALL FDERXY(T,SA,FSA)
      SA=S3+FAC*FSA+FAC*TAU*RO(3)*G0*DEXPO
      SA=SA+FAC*TAU*RO(3)*G1*DEXPO  
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(3)*G2*DEXPO  
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S3=SA

      NLINSYST=NLINSYST+4
	ENDIF

	IF (ISTAGE.GE.4) THEN
* FOURTH STAGE
      UA=U+A(4,1)*S1+A(4,2)*S2+A(4,3)*S3
      TA=T+C(4)*TAU

      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      S4=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)+
     &        Q(4,1)*S1+Q(4,2)*S2+Q(4,3)*S3
      SA=S4+FAC*TAU*RO(4)*G1*DEXPO
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(4)*G2*DEXPO     
      CALL SOLVE_DIREC_Y(YJAC,SA)

      CALL FDERXY(T,SA,FSA)
      SA=S4+FAC*FSA+FAC*TAU*RO(4)*G0*DEXPO
      SA=SA+FAC*TAU*RO(4)*G1*DEXPO   
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(4)*G2*DEXPO       
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S4=SA

      NLINSYST=NLINSYST+4
	ENDIF

	IF (ISTAGE.GE.5) THEN
* FIFTH STAGE
      UA=U+A(5,1)*S1+A(5,2)*S2+A(5,3)*S3+A(5,4)*S4
      TA=T+C(5)*TAU
      
      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      S5=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)
     &  +Q(5,1)*S1+Q(5,2)*S2+Q(5,3)*S3+Q(5,4)*S4
      SA=S5+FAC*TAU*RO(5)*G1*DEXPO
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(5)*G2*DEXPO   
      CALL SOLVE_DIREC_Y(YJAC,SA)

      CALL FDERXY(T,SA,FSA)
        SA=S5+FAC*FSA+FAC*TAU*RO(5)*G0*DEXPO
        SA=SA+FAC*TAU*RO(5)*G1*DEXPO   
      CALL SOLVE_DIREC_X(XJAC,SA)
       SA=SA+FAC*TAU*RO(5)*G2*DEXPO   
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S5=SA

      NLINSYST=NLINSYST+4
	ENDIF

* ADVANCING SOLUTION
      U=U+B(1)*S1+B(2)*S2+B(3)*S3+B(4)*S4+B(5)*S5

      RETURN
	END

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

	SUBROUTINE AMFRW_ONESTEP(THETA,MU,T,TAU,U,XJAC,YJAC,
     &        AS,AV)
* ONE STEP WITH A AMFR-W-METHOD WITH COEFFICIENTS GIVEN IN COMMON/COEFFICIENTS IS PERFORMED
*
* INPUT:     THETA: DIAGONAL COEFFICIENT OF THE  W-METHOD
*                MU: PARAMETER FOR LINEAR REFINEMENT (TYPICALLY MU=THETA)
*	   T,TAU: CURRENT T-POINT AND STEP-SIZE
*        U(N,M): NUMERICAL SOLUTION AT T
*        XJAC, YJAC,...: ARE TRIDIAGONAL MATRICES TO SOLVE LINEAR SYSTEMS
*
* OUTPUT: U(N,M): NUMERICAL SOLUTION AT POINT T+TAU

      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),UA(N,0:M-1)
      DIMENSION FXX(N,0:M-1),FYY(N,0:M-1),FXY(N,0:M-1)
      DIMENSION SA(N,0:M-1),S1(N,0:M-1),S2(N,0:M-1)
      DIMENSION S3(N,0:M-1),S4(N,0:M-1),S5(N,0:M-1)
      DIMENSION FSAX(N,0:M-1),FSAY(N,0:M-1),FSA0(N,0:M-1)

      DIMENSION AS(3,N,0:M-1),AV(5,0:M-1)
      DIMENSION G0(N,0:M-1),G1(N,0:M-1),G2(N,0:M-1)

      REAL*8 MU

        COMMON /PARAMS/XKAPPA,ETA,SIGMA,RHO,RD,RF,XK
	COMMON /COEFFICIENTS/A,Q,C,B,HB,RO,IQ,ISTAGE
	COMMON /BLOCK2/N,M
	COMMON /BLOCK3/NLINSYST

*Constant vectors for the boundary conditions
        CALL GES(G0,G1,G2)

        DEXPO=-RF*DEXP(-RF*T)

      FAC=THETA*TAU
      FACMU=MU*TAU

* FIRST STAGE
      UA=U
      TA=T
      
      EXPO=DEXP(-RF*TA)      
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      S1=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)
      SA=S1+FAC*TAU*RO(1)*G1*DEXPO
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(1)*G2*DEXPO
      CALL SOLVE_DIREC_Y(YJAC,SA)

      CALL FDERXX(TA,AS,SA,FSAX)
      CALL FDERYY(TA,AV,SA,FSAY)
      CALL FDERXY(TA,SA,FSA0)

      SA=2.D0*S1-SA+FACMU*(FSAX+FSAY+FSA0)
      SA=SA+ FACMU*TAU*RO(1)*(G0+G1+G2)*DEXPO

      SA=SA+FAC*TAU*RO(1)*G1*DEXPO
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(1)*G2*DEXPO
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S1=SA

      NLINSYST=NLINSYST+4

	IF (ISTAGE.GE.2) THEN
* SECOND STAGE
      UA=U+A(2,1)*S1
      TA=T+C(2)*TAU
      
      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      S2=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)+Q(2,1)*S1
      SA=S2+FAC*TAU*RO(2)*G1*DEXPO
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(2)*G2*DEXPO
      CALL SOLVE_DIREC_Y(YJAC,SA)

      CALL FDERXX(TA,AS,SA,FSAX)
      CALL FDERYY(TA,AV,SA,FSAY)
      CALL FDERXY(TA,SA,FSA0)

      SA=2.D0*S2-SA+FACMU*(FSAX+FSAY+FSA0)
      SA=SA+ FACMU*TAU*RO(2)*(G0+G1+G2)*DEXPO

      SA=SA+FAC*TAU*RO(2)*G1*DEXPO
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(2)*G2*DEXPO
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S2=SA

      NLINSYST=NLINSYST+4
	ENDIF

	IF (ISTAGE.GE.3) THEN
* THIRD STAGE
      UA=U+A(3,1)*S1+A(3,2)*S2
      TA=T+C(3)*TAU
      
      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      S3=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)+
     &        Q(3,1)*S1+Q(3,2)*S2
      SA=S3+FAC*TAU*RO(3)*G1*DEXPO
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(3)*G2*DEXPO
      CALL SOLVE_DIREC_Y(YJAC,SA)

      CALL FDERXX(TA,AS,SA,FSAX)
      CALL FDERYY(TA,AV,SA,FSAY)
      CALL FDERXY(TA,SA,FSA0)

      SA=2.D0*S3-SA+FACMU*(FSAX+FSAY+FSA0)
      SA=SA+ FACMU*TAU*RO(3)*(G0+G1+G2)*DEXPO   
      
      SA=SA+FAC*TAU*RO(3)*G1*DEXPO
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(3)*G2*DEXPO
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S3=SA

      NLINSYST=NLINSYST+4
	ENDIF

	IF (ISTAGE.GE.4) THEN
* FOURTH STAGE
      UA=U+A(4,1)*S1+A(4,2)*S2+A(4,3)*S3
      TA=T+C(4)*TAU
     
      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      S4=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)+
     &   Q(4,1)*S1+Q(4,2)*S2+Q(4,3)*S3
      SA=S4+FAC*TAU*RO(4)*G1*DEXPO
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(4)*G2*DEXPO
      CALL SOLVE_DIREC_Y(YJAC,SA)

      CALL FDERXX(TA,AS,SA,FSAX)
      CALL FDERYY(TA,AV,SA,FSAY)
      CALL FDERXY(TA,SA,FSA0)

      SA=2.D0*S4-SA+FACMU*(FSAX+FSAY+FSA0)
      SA=SA+ FACMU*TAU*RO(4)*(G0+G1+G2)*DEXPO   
      
      SA=SA+FAC*TAU*RO(4)*G1*DEXPO
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(4)*G2*DEXPO
        CALL SOLVE_DIREC_Y(YJAC,SA)
      S4=SA

      NLINSYST=NLINSYST+4
	ENDIF

	IF (ISTAGE.GE.5) THEN
* FIFTH STAGE
      UA=U+A(5,1)*S1+A(5,2)*S2+A(5,3)*S3+A(5,4)*S4
      TA=T+C(5)*TAU
      
      EXPO=DEXP(-RF*TA)
      CALL FDERXX(TA,AS,UA,FXX)
      CALL FDERYY(TA,AV,UA,FYY)
      CALL FDERXY(TA,UA,FXY)

      S5=TAU*(FXX+FYY+FXY+(G0+G1+G2)*EXPO)+
     &   Q(5,1)*S1+Q(5,2)*S2+Q(5,3)*S3+Q(5,4)*S4
      SA=S5+FAC*TAU*RO(5)*G1*DEXPO
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(5)*G2*DEXPO
      CALL SOLVE_DIREC_Y(YJAC,SA)

      CALL FDERXX(TA,AS,SA,FSAX)
      CALL FDERYY(TA,AV,SA,FSAY)
      CALL FDERXY(TA,SA,FSA0)

      SA=2.D0*S5-SA+FACMU*(FSAX+FSAY+FSA0)
      SA=SA+ FACMU*TAU*RO(5)*(G0+G1+G2)*DEXPO   

      SA=SA+FAC*TAU*RO(5)*G1*DEXPO    
        CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(5)*G2*DEXPO    
      CALL SOLVE_DIREC_Y(YJAC,SA)
      S5=SA

      NLINSYST=NLINSYST+4
	ENDIF

* ADVANCING SOLUTION
      U=U+B(1)*S1+B(2)*S2+B(3)*S3+B(4)*S4+B(5)*S5

      RETURN
	END
