***************************************************************************************
* 02/10/2017
***************************************************************************************

	SUBROUTINE AMFW_ONESTEP(THETA,T,TAU,U,XJAC,YJAC,ZJAC)
* 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(NX,NY,NZ): NUMERICAL SOLUTION AT T
*        XJAC, YJAC,...: ARE TRIDIAGONAL MATRICES TO SOLVE LINEAR SYSTEMS
*        
* OUTPUT: U((NX,NY,NZ): 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,NX),YJAC(3,NY),ZJAC(3,NZ)
	DIMENSION U(NX,NY,NZ),UA(NX,NY,NZ)
	DIMENSION G(NX,NY,NZ),GA(NX,NY,NZ)
      DIMENSION B1(NX,NY,NZ),B2(NX,NY,NZ),B3(NX,NY,NZ) !B0(NX,NY,NZ),
      DIMENSION FXX(NX,NY,NZ),FYY(NX,NY,NZ),FZZ(NX,NY,NZ)
      DIMENSION FXYZ(NX,NY,NZ)!,FSAXYZ(NX,NY,NZ)
      DIMENSION S1(NX,NY,NZ),S2(NX,NY,NZ),S3(NX,NY,NZ)
      DIMENSION S4(NX,NY,NZ),S5(NX,NY,NZ),SA(NX,NY,NZ)

            	COMMON /COEFFICIENTS/A,Q,C,B,HB,RO,IQ,ISTAGE
	COMMON /BLOCK2/NX,NY,NZ
	COMMON /BLOCK3/NLINSYST

      FAC=THETA*TAU
      
      CALL GDER(T,G) ! TIME DERIVATIVE OF F(t,U) AT (t_n,U_n)         

* BOUNDARY VALUES ASSOCIATED TO D11*U_xx+A1*U_X AT THE TIME T.
                CALL FDERXX(T,0.D0*U,B1)	
* BOUNDARY VALUES ASSOCIATED TO D22*U_yy+A2*U_y AT THE TIME T.                
                CALL FDERYY(T,0.D0*U,B2)
* BOUNDARY VALUES ASSOCIATED TO D33*U_zz+A3*U_z AT THE TIME T.                
                CALL FDERZZ(T,0.D0*U,B3)
      
* FIRST STAGE
      UA=U
      TA=T
	GA=G

      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)
      CALL FDERXYZ(TA,UA,FXYZ)
                      
      S1=TAU*(FXX+FYY+FZZ+FXYZ+GA)!+FAC*TAU*RO(1)*(G)
      SA=S1+FAC*TAU*RO(1)*B1 
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(1)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(1)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA) 
      S1=SA  

      NLINSYST=NLINSYST+3	

	IF (ISTAGE.GE.2) THEN
* SECOND STAGE
      UA=U+A(2,1)*S1
      TA=T+C(2)*TAU
	CALL GDER(TA,GA)
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)   
      CALL FDERXYZ(TA,UA,FXYZ)      

      S2=TAU*(FXX+FYY+FZZ+FXYZ+GA)+Q(2,1)*S1!+FAC*TAU*RO(2)*(G)
      SA=S2+FAC*TAU*RO(2)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(2)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(2)*B3  
      CALL SOLVE_DIREC_Z(ZJAC,SA) 
      S2=SA         
      
      NLINSYST=NLINSYST+3
	ENDIF

	IF (ISTAGE.GE.3) THEN
* THIRD STAGE
      UA=U+A(3,1)*S1+A(3,2)*S2
      TA=T+C(3)*TAU
	CALL GDER(TA,GA)
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)
      CALL FDERXYZ(TA,UA,FXYZ)

      S3=TAU*(FXX+FYY+FZZ+FXYZ+GA)+Q(3,1)*S1+Q(3,2)*S2!+FAC*TAU*RO(3)*(G)
      SA=S3+FAC*TAU*RO(3)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(3)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(3)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)       
      S3=SA 
      
      NLINSYST=NLINSYST+3
	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
	CALL GDER(TA,GA)
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)
      CALL FDERXYZ(TA,UA,FXYZ)

      S4=TAU*(FXX+FYY+FZZ+FXYZ+GA)+Q(4,1)*S1+Q(4,2)*S2+Q(4,3)*S3
C     &               +FAC*TAU*RO(4)*(G) 
      SA=S4+FAC*TAU*RO(4)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(4)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(4)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)
      S4=SA     

      NLINSYST=NLINSYST+3
	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
	CALL GDER(TA,GA)
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)
      CALL FDERXYZ(TA,UA,FXYZ)

      S5=TAU*(FXX+FYY+FZZ+FXYZ+GA)+Q(5,1)*S1
     &   +Q(5,2)*S2+Q(5,3)*S3+Q(5,4)*S4!+FAC*TAU*RO(5)*(G)
      SA=S5+FAC*TAU*RO(5)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(5)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(5)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA) 
      S5=SA  
      
      NLINSYST=NLINSYST+3
	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,ZJAC)
* 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(NX,NY,NZ): NUMERICAL SOLUTION AT T
*        XJAC, YJAC,...: TRIDIAGONAL MATRICES TO SOLVE LINEAR SYSTEMS
*        
* OUTPUT: U(NX,NY,NZ): 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,NX),YJAC(3,NY),ZJAC(3,NZ)
	DIMENSION U(NX,NY,NZ),UA(NX,NY,NZ)
	DIMENSION G(NX,NY,NZ),GA(NX,NY,NZ)
      DIMENSION B0(NX,NY,NZ),B1(NX,NY,NZ),B2(NX,NY,NZ),B3(NX,NY,NZ)
      DIMENSION FXX(NX,NY,NZ),FYY(NX,NY,NZ),FZZ(NX,NY,NZ)
      DIMENSION FXYZ(NX,NY,NZ),FSAXYZ(NX,NY,NZ)
      DIMENSION S1(NX,NY,NZ),S2(NX,NY,NZ),S3(NX,NY,NZ)
      DIMENSION S4(NX,NY,NZ),S5(NX,NY,NZ),SA(NX,NY,NZ)

            	COMMON /COEFFICIENTS/A,Q,C,B,HB,RO,IQ,ISTAGE
	COMMON /BLOCK2/NX,NY,NZ	
	COMMON /BLOCK3/NLINSYST

      FAC=THETA*TAU
      
      CALL GDER(T,G) ! TIME DERIVATIVE OF F(t,U) AT (t_n,U_n)         
	
* BOUNDARY VALUES ASSOCIATED TO D11*U_xx+A1*U_X AT THE TIME T.
               CALL FDERXX(T,0.D0*U,B1)
* BOUNDARY VALUES ASSOCIATED TO D22*U_yy+A2*U_y AT THE TIME T.                
                CALL FDERYY(T,0.D0*U,B2)
* BOUNDARY VALUES ASSOCIATED TO D33*U_zz+A3*U_z AT THE TIME T.                
                CALL FDERZZ(T,0.D0*U,B3)
* BOUNDARY VALUES ASSOCIATED TO 2*D12*U_xy+2*D13*U_xz+2*D23*U_yz AT THE TIME T.
               CALL FDERXYZ(T,0.D0*U,B0)

* FIRST STAGE
      UA=U
      TA=T
	GA=G
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)
      CALL FDERXYZ(TA,UA,FXYZ)

      S1=TAU*(FXX+FYY+FZZ+FXYZ+GA)
      SA=S1+FAC*TAU*RO(1)*B1 
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(1)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(1)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)     
      CALL FDERXYZ(T,SA,FSAXYZ) 
      SA=S1+FAC*(FSAXYZ-B0)+FAC*TAU*RO(1)*(G+B0)+FAC*TAU*RO(1)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(1)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(1)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)
      S1=SA  

      NLINSYST=NLINSYST+6	

	IF (ISTAGE.GE.2) THEN
* SECOND STAGE
      UA=U+A(2,1)*S1
      TA=T+C(2)*TAU
	CALL GDER(TA,GA)
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)   
      CALL FDERXYZ(TA,UA,FXYZ)      

      S2=TAU*(FXX+FYY+FZZ+FXYZ+GA)+Q(2,1)*S1
      SA=S2+FAC*TAU*RO(2)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(2)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(2)*B3  
      CALL SOLVE_DIREC_Z(ZJAC,SA)     
      CALL FDERXYZ(T,SA,FSAXYZ)       
      SA=S2+FAC*(FSAXYZ-B0)+FAC*TAU*RO(2)*(G+B0)+FAC*TAU*RO(2)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(2)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(2)*B3  
      CALL SOLVE_DIREC_Z(ZJAC,SA)    
      S2=SA         
      
      NLINSYST=NLINSYST+6
	ENDIF

	IF (ISTAGE.GE.3) THEN
* THIRD STAGE
                UA=U+A(3,1)*S1+A(3,2)*S2
      TA=T+C(3)*TAU
	CALL GDER(TA,GA)
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)
      CALL FDERXYZ(TA,UA,FXYZ)

      S3=TAU*(FXX+FYY+FZZ+FXYZ+GA)+Q(3,1)*S1+Q(3,2)*S2
      SA=S3+FAC*TAU*RO(3)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(3)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(3)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)      
      CALL FDERXYZ(T,SA,FSAXYZ)
      SA=S3+FAC*(FSAXYZ-B0)+FAC*TAU*RO(3)*(G+B0)+FAC*TAU*RO(3)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(3)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(3)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)       
      S3=SA 
      
      NLINSYST=NLINSYST+6
	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
	CALL GDER(TA,GA)
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)
      CALL FDERXYZ(TA,UA,FXYZ)

      S4=TAU*(FXX+FYY+FZZ+FXYZ+GA)+Q(4,1)*S1+Q(4,2)*S2+Q(4,3)*S3
      SA=S4+FAC*TAU*RO(4)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(4)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(4)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)      
      CALL FDERXYZ(T,SA,FSAXYZ)
      SA=S4+FAC*(FSAXYZ-B0)+FAC*TAU*RO(4)*(G+B0)+FAC*TAU*RO(4)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(4)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(4)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)        
      S4=SA     

      NLINSYST=NLINSYST+6	
	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
	CALL GDER(TA,GA)
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)
      CALL FDERXYZ(TA,UA,FXYZ)

      S5=TAU*(FXX+FYY+FZZ+FXYZ+GA)+Q(5,1)*S1
     &   +Q(5,2)*S2+Q(5,3)*S3+Q(5,4)*S4
      SA=S5+FAC*TAU*RO(5)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(5)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(5)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)      
      CALL FDERXYZ(T,SA,FSAXYZ)
      SA=S5+FAC*(FSAXYZ-B0)+FAC*TAU*RO(5)*(G+B0)+FAC*TAU*RO(5)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(5)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(5)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)        
      S5=SA  
      
      NLINSYST=NLINSYST+6
	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,ZJAC)
* 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(NX,NY,NZ): NUMERICAL SOLUTION AT T
*        XJAC, YJAC,...: TRIDIAGONAL MATRICES TO SOLVE LINEAR SYSTEMS
*        
* OUTPUT: U(NX,NY,NZ): 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,NX),YJAC(3,NY),ZJAC(3,NZ)
	DIMENSION U(NX,NY,NZ),UA(NX,NY,NZ)

	DIMENSION G(NX,NY,NZ),GA(NX,NY,NZ)
      DIMENSION B0(NX,NY,NZ),B1(NX,NY,NZ),B2(NX,NY,NZ),B3(NX,NY,NZ)
      DIMENSION FXX(NX,NY,NZ),FYY(NX,NY,NZ),FZZ(NX,NY,NZ)
      DIMENSION FXYZ(NX,NY,NZ)
      DIMENSION S1(NX,NY,NZ),S2(NX,NY,NZ),S3(NX,NY,NZ)
      DIMENSION S4(NX,NY,NZ),S5(NX,NY,NZ),SA(NX,NY,NZ)
      DIMENSION FSAX(NX,NY,NZ),FSAY(NX,NY,NZ),FSAZ(NX,NY,NZ)
      DIMENSION FSA0(NX,NY,NZ)
      REAL*8 MU

	COMMON /COEFFICIENTS/A,Q,C,B,HB,RO,IQ,ISTAGE
	COMMON /BLOCK2/NX,NY,NZ
	COMMON /BLOCK3/NLINSYST

      FAC=THETA*TAU
      FACMU=MU*TAU
      
      CALL GDER(T,G) ! TIME DERIVATIVE OF F(t,U) AT (t_n,U_n)         

* BOUNDARY VALUES ASSOCIATED TO D11*U_xx+A1*U_X AT THE TIME T.	
                CALL FDERXX(T,0.D0*U,B1)
* BOUNDARY VALUES ASSOCIATED TO D22*U_yy+A2*U_y AT THE TIME T.                
                CALL FDERYY(T,0.D0*U,B2)
* BOUNDARY VALUES ASSOCIATED TO D33*U_zz+A3*U_z AT THE TIME T.                
                CALL FDERZZ(T,0.D0*U,B3)
* BOUNDARY VALUES ASSOCIATED TO 2*D12*U_xy+2*D13*U_xz+2*D23*U_yz AT THE TIME T.
               CALL FDERXYZ(T,0.D0*U,B0)
   
* FIRST STAGE
      UA=U
      TA=T
	GA=G
 
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)
      CALL FDERXYZ(TA,UA,FXYZ)
                      
      S1=TAU*(FXX+FYY+FZZ+FXYZ+GA)
      SA=S1+FAC*TAU*RO(1)*B1 
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(1)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(1)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)  
      CALL FDERXX(T,SA,FSAX)
      CALL FDERYY(T,SA,FSAY)
      CALL FDERZZ(T,SA,FSAZ)
      CALL FDERXYZ(T,SA,FSA0) 
      S1=(2.D0)*S1+FACMU*TAU*RO(1)*(G+B0+B1+B2+B3)-SA
     &    +FACMU*(FSAX-B1+FSAY-B2+FSAZ-B3+FSA0-B0)
      SA=S1+FAC*TAU*RO(1)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(1)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(1)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)
      S1=SA 

      NLINSYST=NLINSYST+6		

	IF (ISTAGE.GE.2) THEN
* SECOND STAGE
      UA=U+A(2,1)*S1
      TA=T+C(2)*TAU
	CALL GDER(TA,GA)
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)   
      CALL FDERXYZ(TA,UA,FXYZ)      

      S2=TAU*(FXX+FYY+FZZ+FXYZ+GA)+Q(2,1)*S1
      SA=S2+FAC*TAU*RO(2)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(2)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(2)*B3  
      CALL SOLVE_DIREC_Z(ZJAC,SA) 
      CALL FDERXX(T,SA,FSAX)
      CALL FDERYY(T,SA,FSAY)
      CALL FDERZZ(T,SA,FSAZ)
      CALL FDERXYZ(T,SA,FSA0)
      S2=(2.D0)*S2+FACMU*TAU*RO(2)*(G+B0+B1+B2+B3)-SA
     &    +FACMU*(FSAX-B1+FSAY-B2+FSAZ-B3+FSA0-B0)
      SA=S2+FAC*TAU*RO(2)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(2)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(2)*B3  
      CALL SOLVE_DIREC_Z(ZJAC,SA)    
      S2=SA   
      
      NLINSYST=NLINSYST+6	
	ENDIF

	IF (ISTAGE.GE.3) THEN
* THIRD STAGE
      UA=U+A(3,1)*S1+A(3,2)*S2
      TA=T+C(3)*TAU
	CALL GDER(TA,GA)
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)
      CALL FDERXYZ(TA,UA,FXYZ)

      S3=TAU*(FXX+FYY+FZZ+FXYZ+GA)+Q(3,1)*S1+Q(3,2)*S2
      SA=S3+FAC*TAU*RO(3)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(3)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(3)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA) 
      CALL FDERXX(T,SA,FSAX)
      CALL FDERYY(T,SA,FSAY)
      CALL FDERZZ(T,SA,FSAZ)
      CALL FDERXYZ(T,SA,FSA0)
      S3=(2.D0)*S3+FACMU*TAU*RO(3)*(G+B0+B1+B2+B3)-SA
     &    +FACMU*(FSAX-B1+FSAY-B2+FSAZ-B3+FSA0-B0)
      SA=S3+FAC*TAU*RO(3)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(3)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(3)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)       
      S3=SA 

      NLINSYST=NLINSYST+6	
	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
	CALL GDER(TA,GA)
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)
      CALL FDERXYZ(TA,UA,FXYZ)

      S4=TAU*(FXX+FYY+FZZ+FXYZ+GA)+Q(4,1)*S1+Q(4,2)*S2+Q(4,3)*S3
      SA=S4+FAC*TAU*RO(4)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(4)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(4)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)
      CALL FDERXX(T,SA,FSAX)
      CALL FDERYY(T,SA,FSAY)
      CALL FDERZZ(T,SA,FSAZ)
      CALL FDERXYZ(T,SA,FSA0)      
      S4=(2.D0)*S4+FACMU*TAU*RO(4)*(G+B0+B1+B2+B3)-SA
     &    +FACMU*(FSAX-B1+FSAY-B2+FSAZ-B3+FSA0-B0)
      SA=S4+FAC*TAU*RO(4)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(4)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(4)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)        
      S4=SA     

      NLINSYST=NLINSYST+6
	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
	CALL GDER(TA,GA)
      CALL FDERXX(TA,UA,FXX)
      CALL FDERYY(TA,UA,FYY)
      CALL FDERZZ(TA,UA,FZZ)
      CALL FDERXYZ(TA,UA,FXYZ)

      S5=TAU*(FXX+FYY+FZZ+FXYZ+GA)+Q(5,1)*S1
     &   +Q(5,2)*S2+Q(5,3)*S3+Q(5,4)*S4
      SA=S5+FAC*TAU*RO(5)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(5)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(5)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA) 
      CALL FDERXX(T,SA,FSAX)
      CALL FDERYY(T,SA,FSAY)
      CALL FDERZZ(T,SA,FSAZ)
      CALL FDERXYZ(T,SA,FSA0)
      S5=(2.D0)*S5+FACMU*TAU*RO(5)*(G+B0+B1+B2+B3)-SA
     &    +FACMU*(FSAX-B1+FSAY-B2+FSAZ-B3+FSA0-B0)
      SA=S5+FAC*TAU*RO(5)*B1
      CALL SOLVE_DIREC_X(XJAC,SA)
      SA=SA+FAC*TAU*RO(5)*B2
      CALL SOLVE_DIREC_Y(YJAC,SA)
      SA=SA+FAC*TAU*RO(5)*B3
      CALL SOLVE_DIREC_Z(ZJAC,SA)        
      S5=SA
      
      NLINSYST=NLINSYST+6
	ENDIF 

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

      RETURN
	END
