	PROGRAM DISPLACEMENT
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!double precision is used to possible intensive iterations
	EXTERNAL FUNC
!fUNC defines the Laplace transformed solution
	COMMON/DA/Z
!	OPEN(5,FILE='DISPLACEMENTOUT.TXT',STATUS='OLD')
    OPEN(5,FILE='stdin',STATUS='OLD')
	S=2.0
!S is the real time.Its value can be the physical time with a dimension(t)
!or the dimensionless time without a dimension
!This sample program computes the distribution of lagging temperature in
!space(z) from(z0) to (zf) at S 
	NTERMS=40000
!Maximum number of terms used in the Riemann-sum approximation
!default value is set to 20000
	GMMA=0.0
!GMMA is the value for the product of gmma*beta, A value of 0.0 intrigues the 
!default value of gmma*beta=4.7
	z0=0.0
	zf=10.01
!Define the initial(z0) and final(zf) position for the temperature distribution
	DO 10 I=1,21
	Z=z0+(zf-z0)*(I-1)/20
!Z is the space variable. This sample program discretizes the physical 
!domain from 0 to 2.2 into 100 intervals
	CALL LAPINV(FUNC,GMMA,S,RESULTT,NTERMS)
!LAPINV is the subroutine for the Riemann-sum approximation of the Laplace inversion.
!RESULTT=temperature at (z,s)
	WRITE(5,666)Z,RESULTT
666	FORMAT(1X,6E12.4)
10    CONTINUE
 	CLOSE(5)
 	END
!The function subroutine DISPLACEMENT(P) defines the solutions of displacement in
!the Laplace transform domain. This block needs to be modified for
!defferent problems with different solutions.
 	FUNCTION FUNC(P)
 	IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 	COMMON/DA/Z
 	COMPLEX P,C1,C2,C3,C4,C5,RK1,RK2,RM1,RM2,RM3
!The Laplace transform variable P must be complex.Any other functions 
!of P used in defining the transformed solution of temperature(FACZ in
!this Program for example) mist be declared to be complex accordingly.
	TAU0=0.03
      TAU1=0.05
	ALPHA0=0.
	ALPHA1=0.
	T0=293.0
	V=2.0
	ROU=8954.0
	CE=383.1
	Q0=10.0
      RL=10.0
	ALPHA=1
      K=386
      DLMIDAE=7.76E10
	GGE=3.86E10
	ALPHAS=1.78E-5
      ALPHA11=ALPHA!
     
      W=Q0/V
      DLMIDA=DLMIDAE*(1+ALPHA0*P)
	GG=GGE*(1+ALPHA1*P)
      BETTA1E=(3.*DLMIDAE+2.*GGE)*ALPHAS
      BETTA1=(3.*DLMIDAE*ALPHA0+2.*GGE*ALPHA1)*ALPHAS/BETTA1E
      BETTA=BETTA1E*(1+BETTA1*P)

	g1=(DLMIDAE*(1+ALPHA0*P)+2.*GGE*(1+ALPHA1*P))/(DLMIDAE+2.*GGE)
	g2=((1+BETTA1*P)*(1+TAU0**ALPHA*P**ALPHA/ALPHA11))/(DLMIDAE+2.*GGE)
	g3=P*(1+TAU1**ALPHA*P**ALPHA/ALPHA11)
	g4=(BETTA1E**2*T0*P*(1+BETTA1*P)*(1+TAU1**ALPHA*P**ALPHA/ALPHA11))/(ROU*CE*(DLMIDAE+2.*GGE))
	g5=(BETTA1E*T0*W*(1+TAU1**ALPHA*P**ALPHA/ALPHA11))/(DLMIDAE+2.*GGE)


    RM1=(g1*g3+g2*g4+P**2)/g1
	RM2=g3*P**2/g1
	RM3=g2*g5*P/(g1*V)
	RK1=SQRT((RM1+SQRT(RM1**2-4.*RM2))/2.)
	RK2=SQRT((RM1-SQRT(RM1**2-4.*RM2))/2.)
	C5=RM3/((P/V)**4-RM1*(P/V)**2+RM2)

	C1=((RK2**2-(P/V)**2)*(EXP(RK1*RL)-EXP(-P*RL/V))*C5)/((RK1**2-RK2**2)*(EXP(RK1*RL)-EXP(-RK1*RL)))
    C2=-((RK2**2-(P/V)**2)*(EXP(-RK1*RL)-EXP(-P*RL/V))*C5)/((RK1**2-RK2**2)*(EXP(RK1*RL)-EXP(-RK1*RL)))
	C3=-((RK1**2-(P/V)**2)*(EXP(RK2*RL)-EXP(-P*RL/V))*C5)/((RK1**2-RK2**2)*(EXP(RK2*RL)-EXP(-RK2*RL)))
	C4=((RK1**2-(P/V)**2)*(EXP(-RK2*RL)-EXP(-P*RL/V))*C5)/((RK1**2-RK2**2)*(EXP(RK2*RL)-EXP(-RK2*RL)))

	FUNC=C1*EXP(-RK1*Z)+C2*EXP(RK1*Z)+C3*EXP(-RK2*Z)+C4*EXP(RK2*Z)+C5*EXP(-P*Z/V)

 	RETURN
 	END

!Subroutine LAPINV performs the Riemann-sum to approximate the Laplace inversion
!of the function specified in FUNC(P)
 	SUBROUTINE LAPINV(FUNC,GMMA,S,RESULTT,NTERMS)
 	IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 	EXTERNAL FUNC
 	COMMON/DA/Z
 	COMPLEX GAM,B,CPR,PARTB,CHKCON
 	EPS=1.0D-10
!EPS defines the convergence threshold for the ration test of partial sums,
!EPS=(Temperature(N+1)-Temperatur(N)/Temperature(N),
!with N denoting the partial sum of the first N terms
 	GAM=(0.0,0.0)
!Avoid the use of initial condition at S=0 in this program.
!If needed,select a very small value of S,such as 0.001, to evalue the 
!initial temperature for validating purposes.
 	IF(S.EQ.0.0) THEN
 	WRITE(*,*)'LAPLACE VARIABLE CANNOT BE ZERO!'
 	RETURN
 	ENDIF
!The default value of GMMA*S=4.7
 	IF(GMMA.EQ.0.0)GMMA=4.7/S
 	GAM=GMMA
!Default number of terms used in the Riemann sum is 20000 terms
 	IF(NTERMS.EQ.0) NTERMS=20000
	PI=ACOS(-1.)
	B=(0.0,1.0)
	FIRST=(1./S)*EXP(GAM*S)
	PARTA=0.5*FUNC(GAM)
      PARTB=(0.,0.)
	I=0
!Check convergence for the first NTERMS in the Riemann-sum
!Raise warning flags if EPS is larger than the specified value
5     IF(I.EQ.NTERMS)THEN
 	WRITE(*,*)'NO CONVERGENCE FOR Z=',Z
	GO TO 15
	ENDIF
	I=I+1
	CPR=GAM+B*(I*PI/S)
	CHKODD=MOD(I,2)
	CHKCON=PARTB
	IF(CHKODD.EQ.0)THEN
	PARTB=PARTB+FUNC(CPR)
	ELSE
	PARTB=PARTB+FUNC(CPR)*(-1)
	ENDIF
	RESULT2=FIRST*(PARTA+REAL(PARTB))
	RESULT1=FIRST*(PARTA+REAL(CHKCON))
!If summation is zero,apply different convergence check
 	IF(RESULT1.EQ.0.0)THEN
	IF(ABS(RESULT2).LT.EPS)GO TO 15
	GO TO 5
	ENDIF
!The convergence check is then abs(abs(f(n)-f(n-1))/f(n)).
!first avoid divide by zero error
 	IF(RESULT2.EQ.0.0)THEN
	IF(ABS(RESULT1).LT.EPS)GO TO 15
	GO TO 5
	ENDIF
	CCON=ABS(ABS((RESULT2-RESULT1))/RESULT2)
	IF(CCON.LE.EPS)GO TO 15
	GO TO 5
!If the Laplace variable is changed in successive calculations, such as
!the multiple-time calculations,make sure to reset to zero for the next
!case.
15    GMMA=0.0
 	RESULTT=RESULT1
	RETURN
	END