****************************************************************************************************
* 02/10/2017
* DOUGLAS, MODIFIED CRAIG-SNEYD AND HUNDSDOFER-VERWER SCHEMES 
****************************************************************************************************

      SUBROUTINE DOUGLAS(THETA,T,TAU,U,XJAC,YJAC,ZJAC) 
* TIME INTEGRATION WITH THE DOUGLAS METHOD
*
* INPUT:  THETA: PARAMETER OF THE METHOD
*		T,TAU:   CURRENT POINT AND TIME STEPSIZE
*         U(NX,NY,NZ): NUMERICAL SOLUTION AT T
*               XJAC,YJAC,ZJAC: JACOBIAN AT THE CURRENT POINT
*
* OUTPUT: U(NX,NY,NZ): NUMERICAL SOLUTION AT POINT T+TAU
	
      IMPLICIT REAL*8 (A-H,O-Z)
	COMMON /BLOCK1/D11,D22,D33,D12,D13,D23,A1,A2,A3,CK   
	COMMON /BLOCK2/NX,NY,NZ
	COMMON /BLOCK3/NLINSYST

	DIMENSION XJAC(3,NX),YJAC(3,NY),ZJAC(3,NZ)
	DIMENSION U(NX,NY,NZ),V(NX,NY,NZ)
	DIMENSION FX1(NX,NY,NZ),FY1(NX,NY,NZ),FZ1(NX,NY,NZ),F01(NX,NY,NZ)
	DIMENSION BX1(NX,NY,NZ),BY1(NX,NY,NZ),BZ1(NX,NY,NZ),FG1(NX,NY,NZ)

* G-FUNCTIONS NECESSARY AT POINT (T,U)
	T1=T
	CALL GDER(T1,FG1)
	T2=T+TAU 
	
	CALL FDERXX(T1,U,FX1) 
	CALL FDERYY(T1,U,FY1) 
	CALL FDERZZ(T1,U,FZ1)
	CALL FDERXYZ(T1,U,F01) 

	V=U+TAU*(FX1+FY1+FZ1+F01+FG1)

	CALL FDERXX(T2,0.D0*U,BX1)
	V=V+(TAU*THETA)*(BX1-FX1)
	CALL SOLVE_DIREC_X(XJAC,V) 

	CALL FDERYY(T2,0.D0*U,BY1)
	V=V+(TAU*THETA)*(BY1-FY1)
	CALL SOLVE_DIREC_Y(YJAC,V)
	
	CALL FDERZZ(T2,0.D0*U,BZ1)
	V=V+(TAU*THETA)*(BZ1-FZ1)
	CALL SOLVE_DIREC_Z(ZJAC,V)
	U=V

	NLINSYST=NLINSYST+3

      RETURN
	END

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

      SUBROUTINE MCS(THETA,T,TAU,U,XJAC,YJAC,ZJAC)
* TIME INTEGRATION WITH THE MCS METHOD
*
* INPUT:  THETA: PARAMETER OF THE METHOD
*		T,TAU:   CURRENT POINT AND TIME STEPSIZE
*         U(NX,NY,NZ): NUMERICAL SOLUTION AT T
*               XJAC,YJAC,ZJAC: JACOBIAN AT THE CURRENT POINT
*
* OUTPUT: U(NX,NY,NZ): NUMERICAL SOLUTION AT POINT T+TAU
	
      IMPLICIT REAL*8 (A-H,O-Z)
	COMMON /BLOCK1/D11,D22,D33,D12,D13,D23,A1,A2,A3,CK
	COMMON /BLOCK2/NX,NY,NZ
	COMMON /BLOCK3/NLINSYST

	DIMENSION XJAC(3,NX),YJAC(3,NY),ZJAC(3,NZ)
	DIMENSION U(NX,NY,NZ),V(NX,NY,NZ),V0(NX,NY,NZ)
	DIMENSION FX1(NX,NY,NZ),FY1(NX,NY,NZ),FZ1(NX,NY,NZ),F01(NX,NY,NZ)
	DIMENSION FX2(NX,NY,NZ),FY2(NX,NY,NZ),FZ2(NX,NY,NZ),F02(NX,NY,NZ)
	DIMENSION FG1(NX,NY,NZ),FG2(NX,NY,NZ),F1(NX,NY,NZ),F2(NX,NY,NZ)
	DIMENSION BX1(NX,NY,NZ),BY1(NX,NY,NZ),BZ1(NX,NY,NZ)	

* G-FUNCTIONS NECESSARY AT POINT (T,U)	
	T1=T
	CALL GDER(T1,FG1)
	T2=T+TAU  
	CALL GDER(T2,FG2)
	
	CALL FDERXX(T1,U,FX1) 
	CALL FDERYY(T1,U,FY1) 
	CALL FDERZZ(T1,U,FZ1) 
	CALL FDERXYZ(T1,U,F01)
	F1=FX1+FY1+FZ1+F01+FG1

	V0=U+TAU*F1

                CALL FDERXX(T2,0.D0*U,BX1)
	V=V0+(TAU*THETA)*(BX1-FX1)
	CALL SOLVE_DIREC_X(XJAC,V) 

                CALL FDERYY(T2,0.D0*U,BY1)	
	V=V+(TAU*THETA)*(BY1-FY1)
	CALL SOLVE_DIREC_Y(YJAC,V)
	
	CALL FDERZZ(T2,0.D0*U,BZ1)
	V=V+(TAU*THETA)*(BZ1-FZ1)
	CALL SOLVE_DIREC_Z(ZJAC,V)
	
	CALL FDERXX(T2,V,FX2) 
	CALL FDERYY(T2,V,FY2) 
	CALL FDERZZ(T2,V,FZ2)
	CALL FDERXYZ(T2,V,F02)
	F2=FX2+FY2+FZ2+F02+FG2

	V=V0+(TAU*THETA)*(F02+FG2-F01-FG1)
	V=V+(1.D0/2.D0-THETA)*TAU*(F2-F1)
	
	V=V+(TAU*THETA)*(BX1-FX1)
	CALL SOLVE_DIREC_X(XJAC,V) 

	V=V+(TAU*THETA)*(BY1-FY1)
	CALL SOLVE_DIREC_Y(YJAC,V)
	
	V=V+(TAU*THETA)*(BZ1-FZ1)
	CALL SOLVE_DIREC_Z(ZJAC,V)
	U=V
	
	NLINSYST=NLINSYST+6
	
      RETURN
	END

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

      SUBROUTINE HV(THETA,T,TAU,U,XJAC,YJAC,ZJAC) 
* TIME INTEGRATION WITH THE HV METHOD
*
* INPUT:  THETA: PARAMETER OF THE METHOD
*		T,TAU:   CURRENT POINT AND TIME STEPSIZE
*         U(NX,NY,NZ): NUMERICAL SOLUTION AT T
*               XJAC,YJAC,ZJAC: JACOBIAN AT THE CURRENT POINT
*
* OUTPUT: U(NX,NY,NZ): NUMERICAL SOLUTION AT POINT T+TAU
	
      IMPLICIT REAL*8 (A-H,O-Z)
	COMMON /BLOCK1/D11,D22,D33,D12,D13,D23,A1,A2,A3,CK 
	COMMON /BLOCK2/NX,NY,NZ
	COMMON /BLOCK3/NLINSYST

	DIMENSION XJAC(3,NX),YJAC(3,NY),ZJAC(3,NZ)
	DIMENSION U(NX,NY,NZ),V(NX,NY,NZ),V0(NX,NY,NZ)
	DIMENSION FX1(NX,NY,NZ),FY1(NX,NY,NZ),FZ1(NX,NY,NZ),F01(NX,NY,NZ)
	DIMENSION FX2(NX,NY,NZ),FY2(NX,NY,NZ),FZ2(NX,NY,NZ),F02(NX,NY,NZ)
	DIMENSION FG1(NX,NY,NZ),FG2(NX,NY,NZ),F1(NX,NY,NZ),F2(NX,NY,NZ)
	DIMENSION BX1(NX,NY,NZ),BY1(NX,NY,NZ),BZ1(NX,NY,NZ)

* G-FUNCTIONS NECESSARY AT POINT (T,U)
	T1=T
	CALL GDER(T1,FG1)
	T2=T+TAU  
	CALL GDER(T2,FG2)
	
	CALL FDERXX(T1,U,FX1) 
	CALL FDERYY(T1,U,FY1) 
	CALL FDERZZ(T1,U,FZ1) 
	CALL FDERXYZ(T1,U,F01)
	F1=FX1+FY1+FZ1+F01+FG1

	V0=U+TAU*F1
	
	CALL FDERXX(T2,0.D0*U,BX1)
	V=V0+(TAU*THETA)*(BX1-FX1)
	CALL SOLVE_DIREC_X(XJAC,V) 

	CALL FDERYY(T2,0.D0*U,BY1)
	V=V+(TAU*THETA)*(BY1-FY1)
	CALL SOLVE_DIREC_Y(YJAC,V)
	
	CALL FDERZZ(T2,0.D0*U,BZ1)
	V=V+(TAU*THETA)*(BZ1-FZ1)
	CALL SOLVE_DIREC_Z(ZJAC,V)
	
	CALL FDERXX(T2,V,FX2) 
	CALL FDERYY(T2,V,FY2) 
	CALL FDERZZ(T2,V,FZ2) 
	CALL FDERXYZ(T2,V,F02)
	F2=FX2+FY2+FZ2+F02+FG2

	V=V0+(1.D0/2.D0)*TAU*(F2-F1)
	
	V=V+(TAU*THETA)*(BX1-FX2)
	CALL SOLVE_DIREC_X(XJAC,V) 

	V=V+(TAU*THETA)*(BY1-FY2)
	CALL SOLVE_DIREC_Y(YJAC,V)
	
	V=V+(TAU*THETA)*(BZ1-FZ2)
	CALL SOLVE_DIREC_Z(ZJAC,V)
	U=V
	
	NLINSYST=NLINSYST+6

      RETURN
	END

