****************************************************************************************************
* 02/10/2017
* WE STUDY THE CONVERGENCE ORDER OF THE DOUGLAS, MODIFIED CRAIG-SNEYD AND HUNDSDORFER-VERWER METHODS 
* ON 3d LINEAR PDES OF PARABOLIC TYPE WHERE TIME DEPENDENT BCs AND MIXED DERIVATIVES ARE ALLOWED
****************************************************************************************************
      PROGRAM DOMCSHV3D  

      IMPLICIT REAL*8 (A-H,O-Z)
	CHARACTER LEY*70,FILENAME*64

	PARAMETER (ND=128,NN=20)  
               DIMENSION U(ND,ND,ND),U0(ND,ND,ND),SOL(ND,ND,ND)
      DIMENSION LEY(NN),EGH2(NN),EGHU(NN),PH2(NN),PHU(NN)
	INTEGER IVMET(NN),LINSYST(NN)
	 REAL*4 TIMENS, FINISH, START, TIME(NN) !COMPUTING CPU TIMES WITH CPU_TIME

	COMMON /BLOCK1/D11,D22,D33,D12,D13,D23,A1,A2,A3,CK 
	COMMON /BLOCK2/NX,NY,NZ
	COMMON /BLOCK3/NLINSYST

* SOLVING: 
* U_t=D11*U_xx+D22*U_yy+D33*U_zz
*     +2*D12*U_xy+2*D13*U_xz+2*D23*U_yz
*     +A1*U_x+A2*U_y+A3*U_z+G(x,y,z,t) 
*  -> U'(t)= F(t,U)+G(t) (ODEs)
* (x,y,z) \in (0,1)^3, t\in (0,TEND).
* DIRICHLET BCs AND INITIAL CONDITION 
* SO THAT THE EXACT SOLUTION IS
*     U(x,y,z,t)=(4^3x(1-x)y(1-y)z(1-z) 
*      +k((x+1/3)^2+(y+1/4)^2)+(z+1/5)^2))*exp(t) 
* k=0 Homogeneous BCs, k\ne 0 NonHomog. BCs

* MOL APPROACH
* SPATIAL DISCRETIZATION USING STANDARD 
* CENTRAL DIFFERENCES OVER A UNIFORM GRID: 
* DX=1/(NX+1), DY=1/(NY+1), DZ=1/(NZ+1)

      LEY(1)=' U_t=D11*U_xx+D22*U_yy+D33*U_zz+'
      LEY(2)='     +2*D12*U_xy+2*D13*U_xz+2*D23*U_yz'
      LEY(3)='     +A1*U_x+A2*U_y+A3*U_z+G(x,y,z,t)'
      LEY(4)=' (x,y,z) \in (0,1)^3, t \in (0,TEND]'
      LEY(5)=' EXACT SOLUTION u=(4^3x(1-x)y(1-y)z(1-z)+k*h(x,y,z))*exp(t)'
      LEY(6)=' k\in{0,1}, h(x,y,z)=(x+1/3)^2+(y+1/4)^2+(z+1/5)^2'
      LEY(7)=' Dirichlet BCs and Init.Cond. '
      LEY(8)=' MOL approach. Central Differences in space '
      LEY(9)=' spatial mesh-with:  DX=1/(NX+1),DY=1/(NY+1),DZ=1/(NZ+1)'
      LEY(10)=' Dimensional Splitting'
      LEY(11)=' TAU=TEND/NSTEPS, Time integrations are performed for'
      LEY(12)=' NSTEP(J)=NSTEP0*2^(I-1), I=1,2,...,NJ'

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

* CHOOSE COEFFICIENTS OF THE PDE PROBLEM
	D11=1.D0
	D22=1.D0
	D33=1.D0
	D12=0.9D0
	D13=D12
	D23=D12
	A1=0D0
	A2=0D0
	A3=0D0

* TIME INTERVAL
	T0=0.D0
	TEND=1.D0

* CHOOSE WHICH METHODS ARE TO BE USED
* METHOD=1: DOUGLAS
* METHOD=2: MODIFIED CRAIG-SNEYD
* METHOD=3: HUNDSDORFER-VERWER 	

      NMET=1
      IVMET(1)=3 ! 1:DOUGLAS, 2: MCS, 3: HV
	           
     
* LOOPS FOR SEVERAL VALUES OF CK(=k) AND SEVERAL METHODS

      DO IK=0,1
      DO IMET=1,NMET
	CK=1.D0*IK
	METHOD=IVMET(IMET)

* CHOOSE THETA FOR THE CORRESPONDING METHOD	
	IF (METHOD.EQ.1) THEN
                              THETA=1.D0/2.D0
                ELSEIF (METHOD.EQ.2) THEN                               
                               THETA=1.D0/3.D0
                ELSEIF (METHOD.EQ.3) THEN
                               THETA=0.4020D0
                ENDIF
                                         
* NUMBER OF SPATIAL AND TEMPORAL INTEGRATIONS
      
      NUMINT=1 ! NUMBER OF SPATIAL INTEGRATIONS
	NJ=10    ! NUMBER OF TIME INTEGRATIONS
      NSTEP0=4 ! INITIAL NUMBER OF TIME STEPS

* FILE FOR RESULTS
                ITAG1=FLOOR(10*D12)
	ITAG2=FLOOR(100*THETA)
	
	IF (METHOD.EQ.1) THEN
               WRITE(filename,'(A,I1,A,I1,A,I2,A)')
     &"DOU-K",FLOOR(CK),"-D0",ITAG1,"-TH0",ITAG2,".txt"
               ELSEIF (METHOD.EQ.2) THEN               
                     WRITE(filename,'(A,I1,A,I1,A,I2,A)')
     &"MCS-K",FLOOR(CK),"-D0",ITAG1,"-TH0",ITAG2,".txt"
               ELSEIF (METHOD.EQ.3) THEN               
                     WRITE(filename,'(A,I1,A,I1,A,I2,A,I2,A)')
     &"HV-K",FLOOR(CK),"-D0",ITAG1,"-TH0",ITAG2,".txt"
                    ENDIF               
      OPEN(1,FILE=filename)
      REWIND 1

* WRITING THE HEADING OF THE FILE FOR THE INTEGRATION STATISTICS
	DO L=1,12
            WRITE(1,'(A)') LEY(L)
      ENDDO
      WRITE(1,'(A)') '--------------------------------------------------
     &---'
     	IF (METHOD.EQ.1) THEN
     	WRITE(1,20) THETA
20	FORMAT('    DOUGLAS     THETA=',F8.4)                              
                ELSEIF (METHOD.EQ.2) THEN                               
               WRITE(1,21) THETA
21	FORMAT('    MCS     THETA=',F8.4)   
                ELSEIF (METHOD.EQ.3) THEN
               WRITE(1,22) THETA
22	FORMAT('    HV     THETA=',F8.4)   
                ENDIF
            WRITE(1,'(A)') '--------------------------------------------------
     &---------'

	WRITE(1,'(A),/') ''

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

		WRITE(1,75) D11,D22,D33,D12,D13,D23,A1,A2,A3
75	FORMAT('D11=',F8.3,'	D22=',F8.3,'	D33=',F8.3,
     &'	D12=',F8.3,'	D13=',F8.3,'	D23=',F8.3,
     &'	 A1=',F8.3,'	 A2=',F8.3,'	 A3=',F8.3)

	WRITE(1,'(A),/') ''
	WRITE(1,'(A)') '---------------------------'
	WRITE(1,'(1X,A,D10.3)') ' K=',CK
	WRITE(1,'(A)') '---------------------------'
	WRITE(1,'(A),/') ''

* LOOP FOR INCREASING SPATIAL DIMENSIONS
C	DO INN=1,NUMINT
      DO INN=6,6
	
	U0=0.D0
	
      PRINT*,'SPATIAL-INTEGRATION=',INN
		NX=2**(INN+1)
		NY=2**(INN+1)
		NZ=2**(INN+1)

	WRITE(1,'(A),/') ''
      WRITE(1,'(A)') '--------------------------------------------------
     &---'
    	WRITE(1,100) NX,NY,NZ
100	FORMAT('NX=',I4,'    NY=',I4,'    NZ=',I4)

* COMPUTING THE INITIAL VALUE AND THE EXACT SOLUTION AT END-POINT
	CALL SOLUT(T0,U0)
	CALL SOLUT(TEND,SOL)

	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

	U=0.D0

	NSTEPS=NSTEP0*2**(I-1)
	U=U0	
               
               NLINSYST=0
               CALL CPU_TIME(START)
               CALL METH(THETA,METHOD,NSTEPS,TEND,U)
                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-60+GE2))/DLOG(2.D0)
            ORDU=DLOG(AEU/(1.D-60+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=NX+1=NY+1=NZ+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

	STOP
	END

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

      SUBROUTINE METH(THETA,METHOD,NSTEPS,TEND,U)
* IN THIS ROUTINE DOUGLAS, MCS AND HV 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 XJAC(3,NX),YJAC(3,NY),ZJAC(3,NZ)
      DIMENSION U(NX,NY,NZ)

	COMMON /BLOCK1/D11,D22,D33,D12,D13,D23,A1,A2,A3,CK   
	COMMON /BLOCK2/NX,NY,NZ

	DX=1.D0/(NX+1)
	DY=1.D0/(NY+1) 
	DZ=1.D0/(NZ+1) 
	
	TAU=TEND/NSTEPS 

* PREPARING THE CROUT-LU DECOMPOSITION FOR THE TRIDIAGONAL MATRICES
      DD11=(2*D11/DX**2)*TAU
      UD11= -DD11/2
	UA11=-(A1/(2*DX))*TAU

	DD22=(2*D22/DY**2)*TAU
      UD22= -DD22/2
      UA22=-(A2/(2*DY))*TAU

	DD33=(2*D33/DZ**2)*TAU
      UD33= -DD33/2
      UA33=-(A3/(2*DZ))*TAU	

* JACOBIAN SPLITTING FOR X-DIRECTION
	DO J=1,NX
		XJAC(2,J)=1.D0+DD11*THETA ! MAIN DIAGONAL
		XJAC(1,J)= (UD11+UA11)*THETA !UPPER DIAGONAL
		XJAC(3,J)= (UD11-UA11)*THETA !LOWER DIAGONAL
	ENDDO 
* JACOBIAN SPLITTING FOR Y-DIRECTION
	DO J=1,NY
		YJAC(2,J)=1.D0+DD22*THETA ! MAIN DIAGONAL
		YJAC(1,J)= (UD22+UA22)*THETA !UPPER DIAGONAL
		YJAC(3,J)= (UD22-UA22)*THETA !LOWER DIAGONAL
	ENDDO 
* JACOBIAN SPLITTING FOR Z-DIRECTION	
	DO J=1,NZ
	      ZJAC(2,J)=1.D0+DD33*THETA ! MAIN DIAGONAL
            ZJAC(1,J)= (UD33+UA33)*THETA !UPPER DIAGONAL
            ZJAC(3,J)= (UD33-UA33)*THETA !LOWER DIAGONAL
      ENDDO

* CROUT-LU DECOMPOSITION OF THE TRIDIAGONAL MATRICES
	CALL TRIDIAGONAL(NX,XJAC)
	CALL TRIDIAGONAL(NY,YJAC)
	CALL TRIDIAGONAL(NZ,ZJAC)

* THE TIME INTEGRATION WITH THE CHOSEN AMF VERSION:

               IF (METHOD.EQ.1) THEN
                     DO I=1,NSTEPS
		T=(I-1)*TAU
                          CALL DOUGLAS(THETA,T,TAU,U,XJAC,YJAC,ZJAC) 
                     ENDDO         
                              
               ELSEIF (METHOD.EQ.2) THEN
                     DO I=1,NSTEPS
		T=(I-1)*TAU	
		CALL MCS(THETA,T,TAU,U,XJAC,YJAC,ZJAC) 
	      ENDDO
	      
               ELSEIF (METHOD.EQ.3) THEN
                      DO I=1,NSTEPS
		T=(I-1)*TAU	
		CALL HV(THETA,T,TAU,U,XJAC,YJAC,ZJAC) 
	      ENDDO               
               ENDIF

	RETURN
	END

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

      SUBROUTINE SOLVE_DIREC_X(XJAC,W0)
* THIS ROUTINE SOLVE NY*NZ TRIDIAGONAL 
* LINEAR SYSTEMS OF DIMENSION NX (X-DIRECTION):
* XJAC*X= W0(*,K,L), K=1,...,NY, L=1,...,NZ.
* THE SOLUTION IS STORED IN THE SAME VECTOR W0.
* INPUT: XJAC(3,NX),W0(NX,NY,NZ)
* OUTPUT: W0(NX,NY,NZ)
* W IS AN AUXILIARY VECTOR

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION XJAC(3,NX),W0(NX,NY,NZ),W(NX)
	COMMON /BLOCK2/NX,NY,NZ

	DO K=1,NY
		DO L=1,NZ
		      DO J=1,NX
      			W(J)=W0(J,K,L) !AUXILIARY VECTOR
		      ENDDO

		      CALL TRIDSOLVE(NX,XJAC,W)

		      DO J=1,NX
			      W0(J,K,L)=W(J)
		      ENDDO
      	ENDDO
     	ENDDO
     	
	RETURN
	END

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

      SUBROUTINE SOLVE_DIREC_Y(YJAC,W0)
* THIS ROUTINE SOLVE NX*NZ TRIDIAGONAL 
* LINEAR SYSTEMS OF DIMENSION NY (Y-DIRECTION):
* YJAC*X= W0(J,*,L), J=1,...,NX, L=1,...,NZ.
* THE SOLUTION IS STORED IN THE SAME VECTOR W0.
* INPUT: YJAC(3,NY),W0(NX,NY,NZ)
* OUTPUT: W0(NX,NY,NZ)
* W IS AN AUXILIARY VECTOR

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION YJAC(3,NY),W0(NX,NY,NZ),W(NY)
	COMMON /BLOCK2/NX,NY,NZ

	DO J=1,NX
		DO L=1,NZ
		      DO K=1,NY
      			W(K)=W0(J,K,L) !AUXILIARY VECTOR
		      ENDDO

		      CALL TRIDSOLVE(NY,YJAC,W)

		      DO K=1,NY
			      W0(J,K,L)=W(K)
		      ENDDO
      	ENDDO
	ENDDO

	RETURN
	END
	
***************************************************************************************

      SUBROUTINE SOLVE_DIREC_Z(ZJAC,W0)
* THIS ROUTINE SOLVE NX*NY TRIDIAGONAL 
* LINEAR SYSTEMS OF DIMENSION NZ (Z-DIRECTION):
* ZJAC*X= W0(J,K,*), J=1,...,NX, K=1,...,NY.
* THE SOLUTION IS STORED IN THE SAME VECTOR W0.
* INPUT: ZJAC(3,NZ),NX,NY,NZ,W0(NX,NY,NZ)
* OUTPUT: W0(NX,NY,NZ)
* W IS AN AUXILIARY VECTOR

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION ZJAC(3,NZ),W0(NX,NY,NZ),W(NZ)
	COMMON /BLOCK2/NX,NY,NZ

	DO J=1,NX
		DO K=1,NY
		      DO L=1,NZ
      			W(L)=W0(J,K,L) !AUXILIARY VECTOR
		      ENDDO

		      CALL TRIDSOLVE(NZ,ZJAC,W)

		      DO L=1,NZ
			      W0(J,K,L)=W(L)
		      ENDDO
      	ENDDO
	ENDDO	

	RETURN
	END	

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

	SUBROUTINE GDER(T,G)
* THIS ROUTINE CALCULATES THE FUNCTION G(x,y,z,t) 
* AT THE INTERIOR POINTS (WITHOUT BOUNDARY CONDITIONS)

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION G(NX,NY,NZ)

	COMMON /BLOCK1/D11,D22,D33,D12,D13,D23,A1,A2,A3,CK
	COMMON /BLOCK2/NX,NY,NZ
	
		IFAC=4**3

      DX=1.D0/(NX+1)
	DY=1.D0/(NY+1)
	DZ=1.D0/(NZ+1)

      DD11=2*D11
	DD22=2*D22
	DD33=2*D33
	DD12=2*D12
	DD13=2*D13
	DD23=2*D23

	ET=DEXP(T)

	DD=2*(D11+D22+D33)

   	DO L=1,NZ
      Z=L*DZ
     	PZ=Z*(1.D0-Z)
     	Z2=(Z+0.20D0)**2
     	QZ=1.D0-2*Z
     	RZ=2.D0*(Z+1.D0/5)
         	DO K=1,NY
      	Y=K*DY
	      PY=Y*(1.D0-Y)
     		Y2=(Y+0.25D0)**2
     		QY=1.D0-2*Y
           	RY=2.D0*(Y+1.D0/4)
	      	DO J=1,NX
		      X=J*DX
		      X2=(X+1.D0/3)**2
		      PX=X*(1.D0-X)
		     	QX=1.D0-2*X
     	            RX=2.D0*(X+1.D0/3)

             G(J,K,L)=IFAC*(PX*PY*PZ+DD11*PY*PZ+DD22*PX*PZ+DD33*PX*PY)
		      G(J,K,L)=G(J,K,L)+CK*(-DD+X2+Y2+Z2)
        G(J,K,L)=G(J,K,L)
     & -IFAC*(DD12*QX*QY*PZ+DD13*QX*PY*QZ+DD23*PX*QY*QZ)
        G(J,K,L)=G(J,K,L)-IFAC*(A1*QX*PY*PZ+A2*PX*QY*PZ+A3*PX*PY*QZ)
                  G(J,K,L)=G(J,K,L)-CK*(A1*RX+A2*RY+A3*RZ)
		      G(J,K,L)=G(J,K,L)*ET			      
		      ENDDO
      	ENDDO
	ENDDO

	RETURN
	END

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

      SUBROUTINE FDERXX(T,V,F)
* THIS ROUTINE CALCULATES THE SPLIT-XX FUNCTION (PURE DIFF. AND ADV.)
* OF THE SEMI-DISCRETIZED ODE SYSTEM AT ALL MESH-POINTS
* F=FX(T,V)
* W IS AN ARRAY INCLUDING THE BOUNDARY VALUES FOR V AT TIME T.

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION V(NX,NY,NZ),F(NX,NY,NZ),W(0:NX+1,0:NY+1,0:NZ+1)

	COMMON /BLOCK1/D11,D22,D33,D12,D13,D23,A1,A2,A3,CK
	COMMON /BLOCK2/NX,NY,NZ

      DX=1.D0/(NX+1)
	DX2=1.D0/DX**2

      COEFX2=D11*DX2
      COEFX1=A1/(2.D0*DX)

* EXTENSION OF V TO THE BOUNDARY AT THE TIME T
	CALL BOUND3(T,V,W) 

      DO L=1,NZ
      	DO K=1,NY
	      	DO J=1,NX
		      F(J,K,L)=COEFX2*(W(J+1,K,L)-2*W(J,K,L)+W(J-1,K,L))
     &                  +COEFX1*(W(J+1,K,L)-W(J-1,K,L))
                  ENDDO
		ENDDO
	ENDDO

	RETURN
	END

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

      SUBROUTINE FDERYY(T,V,F)
* THIS ROUTINE CALCULATES THE SPLIT-YY FUNCTION (PURE DIFF. AND ADV.) 
* OF THE SEMI-DISCRETIZED ODE SYSTEM AT ALL MESH-POINTS
* F=FY(T,V)
* W IS AN ARRAY INCLUDING 
* THE BOUNDARY VALUES FOR V AT TIME T.

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION V(NX,NY,NZ),F(NX,NY,NZ),W(0:NX+1,0:NY+1,0:NZ+1)

	COMMON /BLOCK1/D11,D22,D33,D12,D13,D23,A1,A2,A3,CK
	COMMON /BLOCK2/NX,NY,NZ

	DY=1.D0/(NY+1)
	DY2=1.D0/DY**2

	COEFY2=D22*DY2
	COEFY1=A2/(2.D0*DY)

* EXTENSION OF V TO THE BOUNDARY AT THE TIME T
	CALL BOUND3(T,V,W) 

      DO L=1,NZ
      	DO K=1,NY
      		DO J=1,NX
	      	F(J,K,L)=COEFY2*(W(J,K+1,L)-2*W(J,K,L)+W(J,K-1,L))
     &                  +COEFY1*(W(J,K+1,L)-W(J,K-1,L))
                  ENDDO
		ENDDO
	ENDDO

	RETURN
	END
	
***************************************************************************************

      SUBROUTINE FDERZZ(T,V,F)
* THIS ROUTINE CALCULATES THE SPLIT-ZZ FUNCTION (PURE DIFF. AND ADV.)
* OF THE SEMI-DISCRETIZED ODE SYSTEM AT ALL MESH-POINTS
* F=FZ(T,V)
* W IS AN ARRAY INCLUDING 
* THE BOUNDARY VALUES FOR V AT TIME T.

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION V(NX,NY,NZ),F(NX,NY,NZ),W(0:NX+1,0:NY+1,0:NZ+1)

	COMMON /BLOCK1/D11,D22,D33,D12,D13,D23,A1,A2,A3,CK
	COMMON /BLOCK2/NX,NY,NZ

	DZ=1.D0/(NZ+1)
	DZ2=1.D0/DZ**2

	COEFZ2=D33*DZ2
	COEFZ1=A3/(2.D0*DZ)

* EXTENSION OF V TO THE BOUNDARY AT THE TIME T
	CALL BOUND3(T,V,W) 

      DO L=1,NZ
      	DO K=1,NY
      		DO J=1,NX
	      	F(J,K,L)=COEFZ2*(W(J,K,L+1)-2*W(J,K,L)+W(J,K,L-1))
     &                  +COEFZ1*(W(J,K,L+1)-W(J,K,L-1))
                  ENDDO
		ENDDO
	ENDDO

	RETURN
	END	

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

      SUBROUTINE FDERXYZ(T,V,F)
* THIS ROUTINE CALCULATES THE SPLIT-XYZ FUNCTION (MIXED DIFF.)
* OF THE SEMI-DISCRETIZED ODE SYSTEM AT ALL MESH-POINTS
* F=FXY(T,V)
* W IS AN ARRAY INCLUDING 
* THE BOUNDARY VALUES FOR V AT TIME T.

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION V(NX,NY,NZ),F(NX,NY,NZ),W(0:NX+1,0:NY+1,0:NZ+1)
	DIMENSION FXY(NX,NY,NZ),FXZ(NX,NY,NZ),FYZ(NX,NY,NZ)

	COMMON /BLOCK1/D11,D22,D33,D12,D13,D23,A1,A2,A3,CK
	COMMON /BLOCK2/NX,NY,NZ

      DX=1.D0/(NX+1)
	DY=1.D0/(NY+1)
	DZ=1.D0/(NZ+1)

* EXTENSION OF V TO THE BOUNDARY AT THE TIME T
	CALL BOUND3(T,V,W) 

* DISCRETIZATION FOR U_xy
	
	DXY=1.D0/(4.D0*DX*DY)
	COEFXY=2.D0*D12*DXY
      
      DO L=1,NZ
	DO K=1,NY
	DO J=1,NX
      FXY(J,K,L)=W(J+1,K+1,L)-W(J+1,K-1,L)-W(J-1,K+1,L)+W(J-1,K-1,L)
      FXY(J,K,L)=COEFXY*FXY(J,K,L)
	ENDDO
	ENDDO
	ENDDO

* DISCRETIZATION FOR U_xz
	
	DXZ=1.D0/(4.D0*DX*DZ)
	COEFXZ=2.D0*D13*DXZ
	
	DO L=1,NZ
	DO K=1,NY
	DO J=1,NX
      FXZ(J,K,L)=W(J+1,K,L+1)-W(J+1,K,L-1)-W(J-1,K,L+1)+W(J-1,K,L-1)
      FXZ(J,K,L)=COEFXZ*FXZ(J,K,L)
	ENDDO
	ENDDO
	ENDDO
	
* DISCRETIZATION FOR U_yz

	DYZ=1.D0/(4.D0*DY*DZ)
	COEFYZ=2.D0*D23*DYZ
      
      DO L=1,NZ
	DO K=1,NY
	DO J=1,NX
      FYZ(J,K,L)=W(J,K+1,L+1)-W(J,K+1,L-1)-W(J,K-1,L+1)+W(J,K-1,L-1)
      FYZ(J,K,L)=COEFYZ*FYZ(J,K,L)
	ENDDO
	ENDDO
	ENDDO

	F=FXY+FXZ+FYZ
	
	RETURN
	END

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

	SUBROUTINE BOUND3(T,V,W)
* THIS ROUTINE CALCULATES THE EXTENSION W TO THE BOUNDARY OF THE "INTERIOR" ARRAY V AT TIME T.

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION V(NX,NY,NZ), W(0:NX+1,0:NY+1,0:NZ+1)

	COMMON /BLOCK1/D11,D22,D33,D12,D13,D23,A1,A2,A3,CK
	COMMON /BLOCK2/NX,NY,NZ

      DX=1.D0/(NX+1)
	DY=1.D0/(NY+1)
	DZ=1.D0/(NZ+1)

	DO L=1,NZ
		DO K=1,NY
			DO J=1,NX
				W(J,K,L)=V(J,K,L)
			ENDDO
		ENDDO
	ENDDO

	DO K=0,NY+1
		DO J=0,NX+1
			X=J*DX
			Y=K*DY
			W(J,K,0)=EXACT(X,Y,0.D0,CK,T)
			W(J,K,NZ+1)=EXACT(X,Y,1.D0,CK,T)
		ENDDO
	ENDDO

	DO L=1,NZ
		DO J=0,NX+1
			X=J*DX
			Z=L*DZ
			W(J,0,L)=EXACT(X,0.D0,Z,CK,T)
			W(J,NY+1,L)=EXACT(X,1.D0,Z,CK,T)
		ENDDO
	ENDDO

	DO L=1,NZ
		DO K=1,NY
			Y=K*DY
			Z=L*DZ
			W(0,K,L)=EXACT(0.D0,Y,Z,CK,T)
			W(NX+1,K,L)=EXACT(1.D0,Y,Z,CK,T)
		ENDDO
	ENDDO

	RETURN
	END

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

      REAL*8 FUNCTION EXACT(X,Y,Z,C,T)
      IMPLICIT REAL*8 (A-H,O-Z)
* THIS ROUTINE COMPUTES THE EXACT SOLUTION OF THE PDE AT A GIVEN POINT (X,Y,Z,T)

                               IFAC=4**3

	EXACT=IFAC*X*(1.D0-X)*Y*(1.D0-Y)*Z*(1.D0-Z)
      EXACT=EXACT+C*(X+1.D0/3)**2+C*(Y+0.25D0)**2+C*(Z+0.20D0)**2
	EXACT=EXACT*DEXP(T)

	RETURN
	END

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

	SUBROUTINE SOLUT(T,SOL)
* THIS ROUTINE COMPUTES THE EXACT SOLUTION OF THE PDE AT ALL MESH-POINTS

      IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION SOL(NX,NY,NZ)

	COMMON /BLOCK1/D11,D22,D33,D12,D13,D23,A1,A2,A3,CK
	COMMON /BLOCK2/NX,NY,NZ

                IFAC=4**3
                              
      DX=1.D0/(NX+1)
	DY=1.D0/(NY+1)
	DZ=1.D0/(NZ+1)

	ET=DEXP(T)

      DO L=1,NZ
            Z=L*DZ
		PZ=Z*(1.D0-Z)
		Z2=(Z+0.20D0)**2
   	      DO K=1,NY
		      Y=K*DY
      		PY=Y*(1.D0-Y)
      		Y2=(Y+0.25D0)**2
	      	DO J=1,NX
		      	X=J*DX
			      PX=X*(1.D0-X)
			      X2=(X+1.D0/3)**2
            		SOL(J,K,L)=(IFAC*PX*PY*PZ+CK*(X2+Y2+Z2))*ET
	      	ENDDO
      	ENDDO
	ENDDO

	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(NX,NY,NZ)
	COMMON /BLOCK2/NX,NY,NZ

      GEU=0.D0
	GE2=0.D0

	NXYZ=NX*NY*NZ

	DO L=1,NZ
	      DO K=1,NY
      	      DO J=1,NX
	      	      GEU=MAX(GEU,DABS(U(J,K,L)))
            		GE2=GE2+U(J,K,L)**2/NXYZ
            	ENDDO
	      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

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