C*********************************************************************** SUBROUTINE INITIAL(DT,DX,M,Nb,Nf,Q,Delta,Ndir) IMPLICIT NONE INTEGER j, Nb, Nf, Ndir REAL*8 DT, DX, Q(-1:1,0:Nf), Delta, M, expres, r0 C Set the initial conditions of the Solution IF ( Ndir .EQ. -1 ) THEN DO j = 0, Nf Q(1,j) = 0.04D0*EXP( -( DX*j - 2.0D0 )**2 / Delta**2 ) Q(0,j) = 0.04D0*EXP( -( DX*j - 2.0D0 - DT )**2 / Delta**2 ) ENDDO ELSEIF ( M .LT. 0.5D0 ) THEN DO j = 0, Nf expres = DX*j - 2.0D0 / + 4.0D0*M*DLOG( ( DX*j + 1.0D0 - 2.0D0*M ) / / ( 3.0D0 - 2.0D0*M ) ) IF ( DABS(expres/Delta) .GT. 15.0D0 ) THEN Q(1,j) = 0.0D0 ELSE Q(1,j) = 0.04D0*EXP( -1.0D0*expres**2 / Delta**2 ) ENDIF IF ( DABS((expres+DT)/Delta) .GT. 15.0D0 ) THEN Q(0,j) = 0.0D0 ELSE Q(0,j) = 0.04D0*EXP( -1.0D0*( expres + DT )**2 / / Delta**2 ) ENDIF ENDDO ELSEIF ( M .GT. 2.5D0 ) THEN DO j = 0, Nf expres = DX*j - 2.0D0 / + 4.0D0*M*DLOG( ( 2.0D0*M - DX*j - 1.0D0 ) / / ( 2.0D0*M - 3.0D0 ) ) IF ( DABS(expres/Delta) .GT. 15.0D0 ) THEN Q(1,j) = 0.0D0 ELSE Q(1,j) = 0.04D0*EXP( -1.0D0*expres**2 / Delta**2 ) ENDIF IF ( DABS((expres+DT)/Delta) .GT. 15.0D0 ) THEN Q(0,j) = 0.0D0 ELSE Q(0,j) = 0.04D0*EXP( -1.0D0*( expres + DT )**2 / / Delta**2 ) ENDIF ENDDO ELSE Q(1,Nb) = 0.0D0 Q(0,Nb) = 0.0D0 IF ( Nb .NE. Nf ) THEN r0 = Nb*DX + 2.0D0 DO j = Nb+1, Nf expres = DX*j - Nb*DX - 1.0D0 / + 4.0D0*M*DLOG( ( DX*j + 1.0D0 - 2.0D0*M ) / / ( r0 - 2.0D0*M ) ) IF ( DABS(expres/Delta) .GT. 15.0D0 ) THEN Q(1,j) = 0.0D0 ELSE Q(1,j) = 0.04D0*EXP( -1.0D0*expres**2 / Delta**2 ) ENDIF IF ( DABS((expres+DT)/Delta) .GT. 15.0D0 ) THEN Q(0,j) = 0.0D0 ELSE Q(0,j) = 0.04D0*EXP( -1.0D0*( expres + DT )**2 / / Delta**2 ) ENDIF ENDDO ENDIF IF ( Nb .NE. 0 ) THEN r0 = Nb*DX + 0.5D0 DO j = 0, Nb - 1 expres = DX*j - Nb*DX + 0.5D0 / + 4.0D0*M*DLOG( ( 2.0D0*M - DX*j - 1.0D0 ) / / ( 2.0D0*M - r0 ) ) IF ( DABS(expres/Delta) .GT. 15.0D0 ) THEN Q(1,j) = 0.0D0 ELSE Q(1,j) = 0.04D0*EXP( -1.0D0*expres**2 / Delta**2 ) ENDIF IF ( DABS((expres+DT)/Delta) .GT. 15.0D0 ) THEN Q(0,j) = 0.0D0 ELSE Q(0,j) = 0.04D0*EXP( -1.0D0*( expres + DT )**2 / / Delta**2 ) ENDIF ENDDO ENDIF ENDIF RETURN END C*********************************************************************** SUBROUTINE MATRIX(DT,DX,M,thetarr,thetatt,thetar,thetat,Nb,Nuse, / Mat,Work,Ord) IMPLICIT NONE INTEGER j, k, Nb, Nuse, info REAL*8 DT, DX, rho, M, thetarr, thetatt, thetar, thetat, beta, / alpinv, gammar, gammat, r, Mat(Nuse,Nuse), Work(Nuse), / rcond, thxa, thxb, thxc INTEGER Ord(Nuse) C Set the initial conditions of the matrix. DO j = 1, Nuse Work(j) = 0.0D0 DO k = 1, Nuse Mat(j,k) = 0.0D0 ENDDO ENDDO rho = DT / DX IF ( M .GE. 2.5D0 ) THEN Mat(1,Nuse-2) = rho Mat(1,Nuse-1) = -4.0D0*rho Mat(1,Nuse) = 3.0D0*( rho - 1.0D0 ) ELSEIF ( M .GT. 0.5D0 ) THEN Mat(1,Nb-1) = rho Mat(1,Nb) = -4.0D0*rho Mat(1,Nb+1) = 3.0D0*(rho - 1.0D0 ) ELSE Mat(1,1) = 3.0D0*( rho + 1.0D0 ) Mat(1,2) = -4.0D0*rho Mat(1,3) = rho ENDIF DO j = 2, Nuse-1 r = ( j - 1.0D0 )*DX + 1.0D0 beta = 2.0D0*M / r alpinv = 1.0D0 + beta gammar = -1.0D0*beta / r gammat = -1.0D0*gammar thxa = alpinv*thetatt + 0.5D0*gammat*DT*thetat / - thetarr*( 1.0D0 - beta )*rho**2 thxb = beta - 0.5D0*gammar*DT*thetar thxc = alpinv*( 1.0D0 - thetatt ) / + 0.5D0*gammat*DT*( 1.0D0 - thetat ) / + thetarr*( 1.0D0 - beta )*rho**2 Mat(j,j-1) = 0.5D0*( thxa + rho*thxb ) Mat(j,j) = thxc Mat(j,j+1) = 0.5D0*( thxa - rho*thxb ) ENDDO beta = 2.0D0*M / 5.0D0 Mat(Nuse,Nuse-2) = rho*( 1.0D0 - beta ) Mat(Nuse,Nuse-1) = -4.0D0*rho*( 1.0D0 - beta ) Mat(Nuse,Nuse) = 3.0D0*( rho*( 1.0D0 - beta ) + 1.0D0 + beta ) C CALL DGEFA(Mat,Nuse,Nuse,Ord,INFO) C IF ( INFO .NE. 0 ) THEN C PRINT *, 'The Inversion has failed.' C STOP C ENDIF CALL DGECO(Mat,Nuse,Nuse,Ord,rcond,work) PRINT 15, 1.0D0 / rcond 15 FORMAT('The condition number = ',G15.7) RETURN END C*********************************************************************** SUBROUTINE EVOL(DT,DX,M,thetarr,thetatt,thetar,thetat,Nb,Nf,Q, / Nuse,Mat,Vec,Ord) IMPLICIT NONE INTEGER j, Nb, Nf, Nuse REAL*8 DT, DX, rho, M, r, thetarr, thetatt, thetar, thetat, / Q(-1:1,0:Nf), Mat(Nuse,Nuse), Vec(Nuse), alpinv, beta, / gammar, gammat, thxa, thxb, thxc, thxd, thxe, thxf INTEGER Ord(Nuse) C Evolve to the next timestep. DO j = 0, Nf Q(-1,j) = Q(0,j) Q(0,j) = Q(1,j) ENDDO rho = DT / DX IF ( M .GE. 2.5D0 ) THEN Vec(1) = Q(-1,Nf) - 4.0D0*Q(0,Nf) ELSEIF ( M .GT. 0.5D0 ) THEN Vec(1) = Q(-1,Nb) - 4.0D0*Q(0,Nb) ELSE Vec(1) = 4.0D0*Q(0,0) - Q(-1,0) ENDIF DO j = 1, Nf - 1 r = DX*j + 1.0D0 beta = 2.0D0*M / r alpinv = 1.0D0 + beta gammar = -1.0D0*beta / r gammat = -1.0D0*gammar thxa = alpinv*thetatt / + ( 1.0D0 - beta )*( 1.0D0 - thetarr )*rho**2 thxb = 0.5D0*gammar*rho*DT*( 1.0D0 - thetar ) thxc = alpinv*( 1.0D0 - thetatt ) / - ( 1.0D0 - beta )*( 1.0D0 - thetarr )*rho**2 thxd = thetarr*( 1.0D0 - beta )*rho**2 - alpinv*thetatt / + 0.5D0*gammat*thetat*DT thxe = 0.5D0*gammar*thetar*DT + beta thxf = 0.5D0*gammat*DT*( 1.0D0 - thetat ) / - alpinv*( 1.0D0 - thetatt ) / - thetarr*( 1.0D0 - beta )*rho**2 Vec(j+1) = thxa*( Q(0,j+1) + Q(0,j-1) ) + 2.0D0*thxc*Q(0,j) / - thxb*( Q(0,j+1) - Q(0,j-1) ) + thxf*Q(-1,j) / + 0.5D0*thxd*( Q(-1,j+1) + Q(-1,j-1) ) / - 0.5D0*rho*thxe*( Q(-1,j+1) - Q(-1,j-1) ) ENDDO beta = 2.0D0*M / 5.0D0 Vec(Nuse) = ( 1.0D0 + beta )*( 4.0D0*Q(0,Nf) - Q(-1,Nf) ) CALL DGESL(Mat,Nuse,Nuse,Ord,Vec,0) DO j = 0, Nf Q(1,j) = Vec(j+1) ENDDO RETURN END