C *************************************************************************
C .........................................................................
C            1D HYDRODYNAMIC SIMULATION OF THE FREE-FALL COLLAPSE
C            OF A HOMOGENEOUS SPHERE
C
C            (c) Jordi Jose, 2016
C            Stellar Explosions: Hydrodynamics and Nucleosynthesis
C            CRC/Taylor and Francis
C .........................................................................
C *************************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER (N=100)
      COMMON/F0/F1(N+1),F2(N+1),F5(N+1)
      COMMON/F1/DF1W(N+1),DF1RI(N+1),DF1RI0(N+1)
      COMMON/F21/DF2UI(N+1),DF2RI(N+1)
      COMMON/F5/DF5UI(N+1),DF5RI(N+1)
      COMMON/CONST/PI,G
      COMMON/PHYSICS/U(N+1),UU(N+1),R(N+1),RR(N+1)
      DIMENSION UUU(N+1),RRR(N+1)
C
C -------------------------------------------------------------------------
C .............N      = Number of shells of the computational domain.
C                       N is defined in Main Program and in subroutines 
C                       HENYEY3 and ENERGYC 
C ............ F1     = Conservation of mass in finite differences
C ............ F2     = Conservation of momentum in finite differences
C ............ F5     = Lagrangian velocity in finite differences
C ............ DF1W   = Partial derivative of function F1 with respect
C                       to Wi-1/2
C ............ DF1RI  = Partial derivative of function F1 with respect
C                       to Ri
C ............ DF1RI0 = Partial derivative of function F1 with respect
C                       to Ri-1
C ............ DF2UI  = Partial derivative of function F2 with respect
C                       to Ui
C ............ DF2RI  = Partial derivative of function F2 with respect
C                       to Ri
C ............ DF5RI  = Partial derivative of function F5 with respect
C                       to Ri
C ............ DF5UI  = Partial derivative of function F5 with respect
C                       to Ui
C ............ G      = Gravitational constant           (dyn cm^2 g^-2)
C ............ U      = Velocity                         (cm s^-1)
C                       [while U represents velocity at time t_n,
C                       UU is the velocity evaluated at time t_n-1 and
C                       UUU is the velocity at time t_n-2]
C ............ R      = Radius                           (cm)
C                       [while R represents the radius at time t_n,
C                       RR is the radius at time t_n-1 and
C                       RRR is the radius at time t_n-2]
C -------------------------------------------------------------------------
C
      DIMENSION RM(N+1),RHO(N),V(N),W(N)
      DIMENSION RRM(N+1),VV(N),VVV(N),WW(N)
      DIMENSION DM(N+1),Q(N+1),S(N+1),AMM(N+1),QM(N+1)
      DIMENSION FF(N+1),F(N+1),DFRI(N+1)
      DIMENSION AX1(N+1),AX2(N+1)
      DIMENSION HW(N),HU(N+1),HR(N+1)
C
C -------------------------------------------------------------------------
C ............ RM     = Ln(R)
C                       [while RM is evaluated at time t_n,
C                       RRM corresponds to time t_n-1]
C ............ RHO    = Density (1/V)                       (g cm^-3)
C ............ V      = Specific volume (1/RHO)             (cm^3 g^-1)
C                       [while V is evaluated at time t_n,
C                       VV corresponds to time t_n-1, and VVV 
C                       is evaluated at time t_n-2]
C ............ W      = Ln(V)
C                       [while W is evaluated at time t_n,
C                       WW corresponds to time t_n-1]
C ............ DM     = Mass interior to the i-th interface (g)
C ............ Q      = 1 - DM/Total Mass
C ............ S      = Ln(Q)
C ............ AMM    = DSQRT[DMi*DMi-1]                    (g)
C ............ QM     = 1 - AMM/Total Mass
C ............ F      = Auxiliar function F used in the equation of
C                       momentum conservation at time t_n
C ............ FF     = Auxiliar function F at time t_n-1
C ............ DFRI   = Partial derivative of function F with
C                       respect to Ri
C ............ AX1    = Auxiliar function for F1
C ............ AX2    = Auxiliar function for F2
C ............ HW     = Corrections obtained for W
C ............ HU     = Corrections obtained for U
C ............ HR     = Corrections obtained for R
C -------------------------------------------------------------------------
C
      OPEN(UNIT=5,FILE='FREE.INI',STATUS='OLD')
        READ(5,*) ICOI,ICOF,IWRT,DACUR
        READ(5,*) TIMEI,TIMEF,DTIME,EPS
        READ(5,*) BETA,RINI,RHOI
      OPEN(UNIT=7,FILE='OUTFIL7.OUT1',STATUS='NEW')
      OPEN(UNIT=10,FILE='PLOT.OUT1',STATUS='NEW')
      OPEN(UNIT=11,FILE='RADTIM.OUT1',STATUS='NEW')
      OPEN(UNIT=13,FILE='ANALYTIC.OUT1',STATUS='NEW')
C
C -------------------------------------------------------------------------
C ............ ICOI   = Initial model
C ............ ICOF   = Final model
C ............ IWRT   = Label that controls the frequency
C                       of data output
C ............ DACUR  = Precision requirement for the 
C                       iterative procedure
C ............ TIMEI  = Initial time                        (s)
C ............ TIMEF  = Final time                          (s)
C ............ DTIME  = Initial time-step                   (s)
C ............ EPS    = Maximum variation allowance for
C                       physical variables R and V
C ............ BETA   = Parameter of the scheme
C              ...... 0 : Explicit scheme
C              ...... 1 : Implicit scheme
C ............ RINI   = Initial radius of the sphere        (cm)
C ............ RHOI   = Initial density of the sphere       (g cm^-3)
C      
C -------------------------------------------------------------------------
C
      PI=4.D0*DATAN(1.D0)
      DMENV=4.D0*PI*RHOI*(RINI**3.D0)/3.D0
      G=6.67D-8
C
C -------------------------------------------------------------------------
C Physical constants
C ............ DMENV  = Total mass of the sphere           (g)
C ............ G      = Gravitational constant             (dyne cm^2 g^-2)
C -------------------------------------------------------------------------
C
C Discretization of the sphere in N concentric shells. Subroutine ZONING
C assigns a mass DMi+1 - DMi to each individual shell i
C
      CALL ZONING(N,DMENV,DM)
C
C Boundary conditions
C
      R(1)=0.D0
      U(1)=0.D0
C
C Assumption of hydrostatic equilibrium (Ui = 0) and homogeneity of the
C sphere (Vi = 1/RHOI).
C Determination of the radius of each intershell, Ri+1, from mass 
C conservation, using the initial density RHOI and the mass distribution
C obtained in subroutine ZONING
C
      DO 90 I=1,N
      U(I+1)=0.D0
      V(I)=1.D0/RHOI
      R(I+1)=(R(I)**3.D0+(3.D0*(DM(I+1)-DM(I))/4.D0/PI/RHOI))
     &**(1.D0/3.D0)
  90  CONTINUE
C
C -------------------------------------------------------------------------
C Initialization of counters and variables
C ............ ICO   =  Model number
C ............ INDY  =  Number of iterations needed for the current model 
C                       (step)
C ............ IWW   =  Output counter
C ............ TIME  =  Total computed time                 (s)
C -------------------------------------------------------------------------
C
      ICO=ICOI
      INDY=0
      IWW=1
      TIME=TIMEI
C
C Storage of input values and initial model in Unit 7
C
      IF(ICO.EQ.1) THEN
       WRITE(7,1800) ICOI,ICOF,IWRT,DACUR
       WRITE(7,1802) TIMEI,TIMEF,DTIME,EPS
       WRITE(7,1801) BETA,DMENV,RINI,RHOI
       WRITE(7,1900)
       WRITE(7,1902)
       DO 94 J=1,N
       WRITE(7,905) J,1.D0/V(J),R(J+1),U(J+1)
   94  CONTINUE
       WRITE(7,777)
      ENDIF
C
C -------------------------------------------------------------------------
C Analytical solution for the free-fall collapse problem 
C [S.A. Colgate and R.H. White, The Astrophysical Journal (1966), 
C vol. 143, p. 632]. 
C ............ RX    =  Ratio of current total radius of the sphere R(N+1)
C                       over the initial radius RINI
C ............ FUNC  =  TIME*SQRT(8*PI*G*RHOI/3) = SQRT(1-RX)*SQRT(RX) +
C                       ASIN[SQRT(1-RX)]                      
C ...........  FUNX  =  PI/2 - FUNC
C ............ DRX   =  Step used to generate the output in RX
C -------------------------------------------------------------------------
C
      RX=1.0005D0
      DRX=0.0025D0
C 
C Storage of RX and FUNX in Unit 13 
C
      DO 756 I=1,400
      RX=RX-DRX
      FUNC=DASIN((1.D0-RX)**0.5D0)+((1.D0-RX)**0.5D0)*(RX**0.5D0)
      FUNX=PI/2.D0-FUNC
      WRITE(13,790) FUNX,RX
 756  CONTINUE
C
C For each variable Y, we identify Y=YY (i.e., previous value of Y at
C time t_n-1 is equal to the current value of Y at time t_n), since no 
C earlier values for Y are available. 
C
C
      RR(1)=R(1)
      UU(1)=U(1)
C
      DO 8 I=1,N
      RR(I+1)=R(I+1)
      RM(I+1)=DLOG(R(I+1))
      RRM(I+1)=DLOG(RR(I+1))
      VV(I)=V(I)
      UU(I+1)=U(I+1)
      W(I)=DLOG(V(I))
      WW(I)=DLOG(VV(I))
    8 CONTINUE
C     
C Definition of several mass variables used in the calculations
C
      DO 4 I=1,N+1
      Q(I)=1.D0-DM(I)/DMENV
      S(I)=DLOG(Q(I))
    4 CONTINUE
C
      DO 101 I=2,N+1
      AMM(I)=DSQRT(DM(I)*DM(I-1))
      QM(I)=1.D0-AMM(I)/DMENV
  101 CONTINUE
C
C Special treatment for the innermost shell to avoid numerical problems
C when DM1=0 (i.e., simulation extending all the way down to the center 
C of the sphere)
C
      IF (DM(1).EQ.0.D0) THEN
       AMM(2)=DM(2)/2.D0
       QM(2)=1.D0-AMM(2)/DMENV
      ENDIF
C     
C Definition of auxiliar functions AX1 and AX2 that appear in the 
C system of finite difference equations
C
      DO 907 I=2,N+1
      AX1(I)=4.D0*PI/(3.D0*DMENV*(Q(I)-Q(I-1)))
  907 CONTINUE
C
      DO 5 I=2,N
      AX2(I)=4.D0*PI/(DMENV*(QM(I+1)-QM(I)))
      FF(I)=-G*DM(I)/RR(I)/RR(I)
    5 CONTINUE
      AX2(N+1)=4.D0*PI/(DMENV*(Q(N+1)-QM(N+1)))
      FF(N+1)=-G*DMENV/RR(N+1)/RR(N+1)
C
C -------------------------------------------------------------------------
C Energy budget 
C ............ EKIN  =  Kinetic energy                      (erg)
C ............ EGRAV =  Gravitational potential energy      (erg)
C ............ EOLD  =  Total initial energy                (erg)
C -------------------------------------------------------------------------
C
      EKIN=0.D0
      EGRAV=0.D0
C      
      DO 813 I=2,N+1
      EKIN=EKIN+DABS(UU(I)*UU(I-1))*(DM(I)-DM(I-1))/2.D0
      IF (RR(I-1).NE.0.D0) THEN
       EGRAV=EGRAV-G*(DM(I)-DM(I-1))*AMM(I)/DSQRT(RR(I)*RR(I-1))
      ELSE
       EGRAV=0.D0
      ENDIF
  813 CONTINUE
      EOLD=EKIN+EGRAV
      WRITE(7,793) EOLD
C
C -------------------------------------------------------------------------
C                        M A I N           L O O P
C -------------------------------------------------------------------------
C For each physical magnitud Y, the program determines the set of
C corrections HY that should be applied to the first guess value,
C by means of an iterative procedure
C
  309 CONTINUE
C
      INDY=INDY+1
C
C -------------------------------------------------------------------------
C Innermost and Intermediate Shells
C -------------------------------------------------------------------------
C All physical magnitudes, such as V (or W), that are evaluated at the 
C middle of a shell, i+1/2, are written as:
C     V(I-1/2) --> V(I-1)
C     V(I+1/2) --> V(I),             I=1, N
C      
      DO 160 I=2,N
C
C Function F1 and its derivatives
C
      F1(I)=AX1(I)*(R(I)*R(I)*R(I)-R(I-1)*R(I-1)*R(I-1))+V(I-1)
      DF1W(I)=V(I-1)
      DF1RI(I)=AX1(I)*3.D0*R(I)*R(I)*R(I)
      DF1RI0(I)=-AX1(I)*3.D0*R(I-1)*R(I-1)*R(I-1)
C
C Function F2 and its derivatives
C
      F(I)=-G*DM(I)/R(I)/R(I)
      F2(I)=(U(I)-UU(I))/DTIME-(1.D0-BETA)*FF(I)-BETA*F(I)
      DFRI(I)=2.D0*G*DM(I)/R(I)/R(I)/R(I)
      DF2UI(I)=1.D0/DTIME
      DF2RI(I)=-R(I)*BETA*DFRI(I)
C
C Function F5 and its derivatives
C
      F5(I)=(RM(I)-RRM(I))/DTIME-BETA*U(I)/R(I)-(1.D0-BETA)*UU(I)/
     &RR(I)
      DF5UI(I)=-BETA/R(I)
      DF5RI(I)=1.D0/DTIME+BETA*U(I)/R(I)
C
 160  CONTINUE
C
C -------------------------------------------------------------------------
C Surface Shell
C -------------------------------------------------------------------------
C
C Function F1 at the surface and its derivatives
C
      F1(N+1)=AX1(N+1)*(R(N+1)*R(N+1)*R(N+1)-R(N)*R(N)*R(N))+V(N)
      DF1W(N+1)=V(N)
      DF1RI(N+1)=AX1(N+1)*3.D0*R(N+1)*R(N+1)*R(N+1)
      DF1RI0(N+1)=-AX1(N+1)*3.D0*R(N)*R(N)*R(N)
C
C Function F2 at the surface and its derivatives
C
      F(N+1)=-G*DMENV/R(N+1)/R(N+1)
      F2(N+1)=(U(N+1)-UU(N+1))/DTIME-(1.D0-BETA)*FF(N+1)-BETA*F(N+1)
      DFRI(N+1)=2.D0*G*DMENV/R(N+1)/R(N+1)/R(N+1)
      DF2UI(N+1)=1.D0/DTIME
      DF2RI(N+1)=-R(N+1)*BETA*DFRI(N+1)
C
C Function F5 at the surface and its derivatives
C
      F5(N+1)=(RM(N+1)-RRM(N+1))/DTIME-BETA*U(N+1)/R(N+1)-
     &(1.D0-BETA)*UU(N+1)/RR(N+1)
      DF5UI(N+1)=-BETA/R(N+1)
      DF5RI(N+1)=1.D0/DTIME+BETA*U(N+1)/R(N+1)
C
C Solution of the system based on Henyey method. Corrections HY 
C are applied to the first guess values of Y
C
      CALL HENYEY3(HW,HU,HR)
C
      DO 35 I=1,N
      W(I)=W(I)+HW(I)
      U(I+1)=U(I+1)+HU(I+1)
      RM(I+1)=RM(I+1)+HR(I+1)
   35 CONTINUE
C
      DO 41 I=1,N
      V(I)=EXP(W(I))
      R(I+1)=EXP(RM(I+1))
   41 CONTINUE
C
C -------------------------------------------------------------------------
C Checkpoint for the corrections. If HY > DACUR, convergence
C has not yet been reached and the procedure is iterated with
C the corrected guess values of Y, Y' = Y + HY.
C ............ ICONT  = Accuracy counter
C              ...... 0 : Accuracy criterion is satisfied
C              ...... 2 : Accuracy criterion in not satisfied and
C                         model is forced to another iteration
C -------------------------------------------------------------------------
C
      ICONT=0
      DO 260 K=1,N
      IF(DABS(HW(K)*2.D0/(W(K)+WW(K))).GT.DACUR) THEN
       ICONT=2
      ENDIF
      IF(DABS(HU(K+1)*2.D0/(U(K+1)+UU(K+1))).GT.DACUR) THEN
       ICONT=2
      ENDIF
      IF(DABS(HR(K+1)*2.D0/(RM(K+1)+RRM(K+1))).GT.DACUR) THEN
       ICONT=2
      ENDIF
  260 CONTINUE
C
C If convergence is not reached after 10 iterations, the adopted
C time-step is reduced by a factor of 10 and the overall iteration
C process is reinitialized for this model
C
      IF(INDY.GE.10) THEN
      DTIME=DTIME/10.D0
      ICONT=2
      ENDIF
C
      IF(ICONT.EQ.2) GOTO 309
C
C -------------------------------------------------------------------------
C Storage of results in Units 7, 10, and 11, whenever IWW = IWRT
C UNIT  7: output file containing density (1/V), radius (R), and velocity
C          (V) at t=TIME.
C UNIT 10: output file aimed at comparing the numerical and the  
C          analytical solutions of the free-fall collapse problem 
C          [see UNIT 13]. It also includes the density at the innermost 
C          and outermost shells to check whether the sphere remains 
C          homogeneous during collapse, as proved analytically.
C
C ............ DA    =  Ratio of current total radius of the sphere 
C                       over the initial radius [DA = RX]
C ............ DB    =  PI/2 - TIME*SQRT(8*PI*G*RHOI/3) = 
C                       PI/2 - SQRT(1-RX)*SQRT(RX)-ASIN[SQRT(1-RX)] 
C                       [DB = FUNX]
C UNIT 11: output data listing the overall size of the sphere,
C          R(N+1), as a function of time (TIME-TIMEI). 
C -------------------------------------------------------------------------
C
      WRITE(7,*) '********* ICO=',ICO
      WRITE(7,*) '........ INDY=',INDY
C      
      IF(IWW.EQ.IWRT) THEN
       WRITE(7,1806) ICO,INDY,TIME,DTIME
       WRITE(7,1902)
       DO 92 J=1,N
       WRITE(7,905) J,1.D0/V(J),R(J+1),U(J+1)
   92  CONTINUE
C
       DA=R(N+1)/RINI
       DB=PI/2.D0-(TIME-TIMEI)*((8.D0*PI*G*RHOI/3.D0)**(1./2.))
       WRITE(10,990) DB,DA,1.D0/V(1),1.D0/V(N)
       WRITE(11,790) TIME-TIMEI,R(N+1)
C
C -------------------------------------------------------------------------
C Check for energy conservation
C ............ ENEW  =  Total energy at t=TIME                  (erg)
C -------------------------------------------------------------------------
C
      CALL ENERGYC(DM,AMM,EOLD,ENEW)
C
C Initialization of counter IWW
C
       IWW=0
      ENDIF
C
C Writing last configuration after completion of the simulation
C
      IF((ICO.GT.ICOF).OR.(TIME.GT.TIMEF)) GOTO 310
C
C -------------------------------------------------------------------------
C Time-step evaluation
C ............ DDTIME  = Current time-step                      (s) 
C ............ DTIME   = First guess for the new time-step      (s) 
C -------------------------------------------------------------------------
C
      DDTIME=DTIME
      CALL TSTEP(DDTIME,EPS,V,VV,R,RR,DTIME,N)
C
C Redefinition of variables at previous time and first guess;
C linear extrapolation of the new values at t(new)=t(old)+DTIME
C
      DO 78 I=1,N
      RRR(I+1)=RR(I+1)
      UUU(I+1)=UU(I+1)
      VVV(I)=VV(I)
   78 CONTINUE
C
      DO 76 I=1,N
      RR(I+1)=R(I+1)
      RRM(I+1)=RM(I+1)
      UU(I+1)=U(I+1)
      VV(I)=V(I)
      WW(I)=W(I)
      FF(I+1)=F(I+1)
   76 CONTINUE
C
      DO 77 J=1,N
      R(J+1)=RR(J+1)+((RR(J+1)-RRR(J+1))*DTIME/DDTIME)
      U(J+1)=UU(J+1)+((UU(J+1)-UUU(J+1))*DTIME/DDTIME)
      V(J)=VV(J)+((VV(J)-VVV(J))*DTIME/DDTIME)
   77 CONTINUE
C
C Update on total time and model counters 
C
      TIME=TIME+DTIME
      ICO=ICO+1
      IWW=IWW+1
      INDY=0
C
C -------------------------------------------------------------------------
C Density control  
C ............ INEG  = check counter
C              ....... 0 : All shells have positive densities
C              ....... 5 : At least one shell has a negative density,
C                          and execution is aborted
C -------------------------------------------------------------------------
C
      INEG=0
      DO 3934 K=1,N
      IF(V(K).LT.0.D0) INEG=5
 3934 CONTINUE
C
      IF(INEG.EQ.5) THEN
      WRITE(7,*) 'Numerical problem (RHO<0). Execution aborted'
      GOTO 310
      ENDIF
C
      GOTO 309
C
  310 CONTINUE
C
C     Formats
C
  777 FORMAT(70('*'))
  790 FORMAT(2(2X,1PD12.5))
  793 FORMAT(2X,'INITIAL ENERGY=',1PD11.4)
  905 FORMAT(2X,I8,3(2X,1PD12.5))
  990 FORMAT(2(2X,1PD12.5),2(2X,1PD14.7))
 1800 FORMAT(2X,'INITIAL MODEL =',I8,2X,'FINAL MODEL =',I8,/,2X,'PRINT
     &STEP =',I8,2X,'PRECISION =',1PD10.3)
 1801 FORMAT(2X,'BETA =',1PD10.3,2X,'MASS =',1PD10.3,/,2X,'RADIUS =',
     &1PD10.3,2X,'DENSITY =',1PD10.3)
 1802 FORMAT(2X,'INITIAL TIME =',1PD10.3,2X,'FINAL TIME =',1PD10.3,/,
     &2X,'TIME-STEP =',1PD10.3,2X,'TIME-STEP VARIATION =',1PD10.3)
 1806 FORMAT(2X,'ICO = ',I8,2X,'ITERATIONS/MODEL =',I8,2X,
     &/,'TIME =',1PD12.5,2X,'DTIME =',1PD12.5)
 1900 FORMAT(2X,'INITIAL MODEL:')
 1902 FORMAT(2X,'     SHELL   DENSITY       RADIUS        VELOCITY')
 4780 FORMAT(2X,I5,3(2X,1PD10.3))
C
      STOP
      END
C -------------------------------------------------------------------------
C -------------------------------------------------------------------------
C                         SUBROUTINES
C -------------------------------------------------------------------------
C -------------------------------------------------------------------------
C ******************************************************************
      SUBROUTINE ZONING(N,DMV,DM)
C 
C Discretization of the sphere in N concentric shells, with N+1
C intershells. For simplicity, all shells are assumed to contain
C the same mass. Note that for other applications, it may be 
C convenient to use an unequal mass grid (i.e., lower-mass shells 
C close to the center or at the outer edge, where more resolution 
C is often needed).
C ...........N    = Number of shells of the computational
C                   domain
C .......... DMV  = Total mass of the sphere                (g)
C .......... DMS  = Mass per shell                          (g)
C .......... DM   = Mass interior to the i-th intershell    (g)
C ******************************************************************
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DM(N+1)
C
      DMS=DMV/DFLOAT(N)
C
      DM(1)=0.D0
      DO 401 L=1,N
      DM(L+1)=DFLOAT(L)*DMS
  401 CONTINUE
C  
      RETURN
      END
C
C ******************************************************************
      SUBROUTINE SCALING(XC,DFAC,ND)
C
C Re-scaling of matrix X to improve accuracy 
C .......... XC   = Input matrix
C .......... ND   = Dimension of matrix XC [NDxND]
C .......... DFAC = Re-scaling factor per row
C ******************************************************************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION DFAC(3),XC(3,3)
C
      DO 13 I=1,ND
      DFAC(I)=DABS(XC(I,1))
      DO 14 K=1,ND
      IF(DABS(XC(I,K)).GE.DFAC(I)) DFAC(I)=DABS(XC(I,K))
  14  CONTINUE
  13  CONTINUE
C
      SUP=DFAC(1)
C
      DO 15 I=1,ND
      IF (DFAC(I).GE.SUP) SUP=DFAC(I)
  15  CONTINUE
C
      DO 16 I=1,ND
      IF (DFAC(I).NE.0.D0) THEN
      DFAC(I)=SUP/DFAC(I)
       ELSE
      DFAC(I)=1.D0
      ENDIF
  16  CONTINUE
C
      RETURN
      END
C
C ***************************************************************************
      SUBROUTINE ENERGYC(DM,AMM,EOLD,ENEW)
C
C This subroutine checks whether energy is conserved
C .......... N      = Number of shells of the computational
C                     domain
C .......... ENEW   = Total energy at t=TIME                  (erg)
C .......... EKIN   = Kinetic energy                          (erg)
C .......... EGRAV  = Gravitational potential energy          (erg)
C .......... EOLD   = Total initial energy                    (erg)
C .......... ERATIO = Relative error (in %) in energy 
C                     conservation
C .......... G      = Gravitational constant                  (dyn cm^2 g-2)
C .......... U      = Velocity                                (cm s^-1)
C                     [while U represents velocity at time 
C                     t_n, UU is the velocity at time t_n-1] 
C .......... R      = Radius                                  (cm)
C                     [R represents radius at time t_n;
C                     RR is the radius at time t_n-1] 
C .......... DM     = Mass interior to the i-th intershell    (g)
C .......... AMM    = DSQRT[DMi*DMi-1]                        (g)
C ***************************************************************************
C
      PARAMETER (N=100)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON/CONST/PI,G
      COMMON/PHYSICS/U(N+1),UU(N+1),R(N+1),RR(N+1)
      DIMENSION AMM(N+1),DM(N+1)
C
C Initialization of kinetic and gravitational potential
C energy terms
C
      EKIN=0.D0
      EGRAV=0.D0
C
      DO 10 I=2,N+1
      EKIN=EKIN+DABS(U(I)*U(I-1))*(DM(I)-DM(I-1))/2.D0
      IF (RR(I-1).NE.0.D0) THEN
       EGRAV=EGRAV-G*(DM(I)-DM(I-1))*AMM(I)/DSQRT(R(I)*R(I-1))
      ELSE
       EGRAV=0.D0
      ENDIF
 10   CONTINUE
C
      ENEW=EKIN+EGRAV
C
      ERATIO=(ENEW-EOLD)*100.D0/EOLD
      WRITE(7,*) '*********** ERATIO (%)=',ERATIO
C  
      RETURN
      END
C
C *****************************************************************************
      SUBROUTINE HENYEY3(DW,DU,DR)
C
C Subroutine that solves the matrix equation A*B = C based on Henyey
C method [L.G. Henyey, J.E. Forbes and N.L. Gould, The Astrophysical 
C Journal (1964), vol. 139, p. 306; L.G. Henyey, L. Wilets, K.H. Bohm, 
C R. LeLevier and R.D. Levee, The Astrophysical Journal (1959), vol. 129, 
C p. 628]. 
C
C .......... N      = Number of shells of the computational domain
C .......... DW     = Corrections obtained for W
C .......... DU     = Corrections obtained for U
C .......... DR     = Corrections obtained for R
C .......... F1     = Conservation of mass in finite differences
C .......... F2     = Conservation of momentum in finite differences
C .......... F5     = Lagrangian velocity in finite differences
C .......... DF1W   = Partial derivative of function F1 with respect to Wi-1/2
C .......... DF1RI  = Partial derivative of function F1 with respect to Ri
C .......... DF1RI0 = Partial derivative of function F1 with respect to Ri-1
C .......... DF2UI  = Partial derivative of function F2 with respect to Ui
C .......... DF2RI  = Partial derivative of function F2 with respect to Ri
C .......... DF5RI  = Partial derivative of function F5 with respect to Ri
C .......... DF5UI  = Partial derivative of function F5 with respect to Ui
C .......... X      = Matrix that contains the partial derivatives of F1, F2,
C                     and F5, for the innermost shell
C .......... Z      = Matrix that contains the values of F1, F2, and F5, for 
C                     the innermost shell
C .......... Y      = Matrix of unknowns DW, DU, and DR, for the innermost 
C                     shell [ X * Y = Z ]
C .......... XBUR   = Matrix X after rescaling
C .......... ZBUR   = Matrix Z after rescaling
C .......... W      = Matrix that contains the partial derivatives of F1, F2,
C                     and F5, for the N-2 intermediate shells
C .......... S      = Matrix that contains the values of F1, F2, and F5, for 
C                     the intermediate shells
C .......... T      = Matrix of unknowns DW, DU, and DR, for the intermediate 
C                     shells [ W * T = S ]
C .......... WBUR   = Matrix W after rescaling
C .......... SBUR   = Matrix S after rescaling
C .......... F      = Matrix that contains the partial derivatives of F1, F2,
C                     and F5, for the outermost shell
C .......... G      = Matrix that contains the values of F1, F2, and F5, for 
C                     the outermost shell
C .......... V      = Matrix of unknowns DW, DU, and DR, for the outermost 
C                     shell [ F * V = G ]
C .......... FBUR    = Matrix F after rescaling
C .......... GBUR    = Matrix G after rescaling
C .......... DFC     = Rescaling factor for matrices XBUR and ZBUR
C .......... DFC2    = Rescaling factor for matrices WBUR and SBUR
C .......... DFC3    = Rescaling factor for matrices FBUR and GBUR
C *****************************************************************************
C
      PARAMETER (N=100)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON/F0/F1(N+1),F2(N+1),F5(N+1)
      COMMON/F1/DF1W(N+1),DF1RI(N+1),DF1RI0(N+1)
      COMMON/F21/DF2UI(N+1),DF2RI(N+1)
      COMMON/F5/DF5UI(N+1),DF5RI(N+1)
      DIMENSION X(3,3),Z(3),Y(3)
      DIMENSION W(3,3),S(3),T(3)
      DIMENSION F(3,3),G(3),V(3)
      DIMENSION DW(N),DR(N+1),DU(N+1)
      DIMENSION XBUR(3,3),ZBUR(3),DFC(3)
      DIMENSION WBUR(3,3),SBUR(3),DFC2(3)
      DIMENSION FBUR(3,3),GBUR(3),DFC3(3)
C
C Matrix form of the linearized system of equations for the innermost 
C shell of the sphere [ X * Y = Z ]
C
      X(1,1)=DF1W(2)
      X(1,2)=0.D0
      X(1,3)=DF1RI(2)
      X(2,1)=0.D0
      X(2,2)=DF2UI(2)
      X(2,3)=DF2RI(2)
      X(3,1)=0.D0
      X(3,2)=DF5UI(2)
      X(3,3)=DF5RI(2)
C
      Z(1)=-F1(2)
      Z(2)=-F2(2)
      Z(3)=-F5(2)
C
C Rescaling of matrices X and Z to improve accuracy 
C
        CALL SCALING(X,DFC,3)
        DO 55 J=1,3
        DO 56 I=1,3
        XBUR(J,I)=X(J,I)*DFC(J)
  56    CONTINUE
  55    CONTINUE
        DO 57 K=1,3    
        ZBUR(K)=Z(K)*DFC(K)
  57    CONTINUE
C
C Solution of the matrix system XBUR * Y = ZBUR 
C through Gaussian elimination
C
      CALL GAUSS(XBUR,ZBUR,Y,3)
C
      DW(1)=Y(1)
      DU(2)=Y(2)
      DR(2)=Y(3)
C
C Matrix form of the linearized system of equations for the N-2 
C intermediate shells of the sphere [ W * T = S ]. Correcting terms
C obtained for the innermost shell are moved to matrix S
C
      DO 20 I=3,N
      W(1,1)=DF1W(I)
      W(1,2)=0.D0
      W(1,3)=DF1RI(I)
      W(2,1)=0.D0
      W(2,2)=DF2UI(I)
      W(2,3)=DF2RI(I)
      W(3,1)=0.D0
      W(3,2)=DF5UI(I)
      W(3,3)=DF5RI(I)
C     
      S(1)=-F1(I)-DF1RI0(I)*DR(I-1)
      S(2)=-F2(I)
      S(3)=-F5(I)
C
C Rescaling of matrices W and S to improve accuracy 
C
      CALL SCALING(W,DFC2,3)
      DO 855 J=1,3
      DO 856 IL=1,3
      WBUR(J,IL)=W(J,IL)*DFC2(J)
  856 CONTINUE
  855 CONTINUE
      DO 857 K=1,3    
      SBUR(K)=S(K)*DFC2(K)
  857 CONTINUE
C
C Solution of the matrix system WBUR * T = SBUR
C through Gaussian elimination
C
      CALL GAUSS(WBUR,SBUR,T,3)
C
      DW(I-1)=T(1)
      DU(I)=T(2)
      DR(I)=T(3)
  20  CONTINUE
C
C Matrix form of the linearized system of equations for the outermost 
C shell of the sphere [ F * V = G ]
C
      F(1,1)=DF1W(N+1)
      F(1,2)=0.D0
      F(1,3)=DF1RI(N+1)
      F(2,1)=0.D0
      F(2,2)=DF2UI(N+1)
      F(2,3)=DF2RI(N+1)
      F(3,1)=0.D0
      F(3,2)=DF5UI(N+1)
      F(3,3)=DF5RI(N+1)
C     
      G(1)=-F1(N+1)-DF1RI0(N+1)*DR(N)
      G(2)=-F2(N+1)
      G(3)=-F5(N+1)      
C
C Rescaling of matrices F and G to improve accuracy 
C
      CALL SCALING(F,DFC3,3)
      DO 955 J=1,3
      DO 136 I=1,3
      FBUR(J,I)=F(J,I)*DFC3(J)
  136 CONTINUE
  955 CONTINUE
      DO 255 J=1,3
      GBUR(J)=G(J)*DFC3(J)
  255 CONTINUE
C
C Solution of the matrix system FBUR * V = GBUR
C through Gaussian elimination
C
      CALL GAUSS(FBUR,GBUR,V,3)
C
      DW(N)=V(1)
      DU(N+1)=V(2)
      DR(N+1)=V(3)
C
      RETURN
      END
C
C *************************************************************************
      SUBROUTINE TSTEP(DTOLD,EPS,V,VV,R,RR,DTIME,N)
C      
C .......... N      = Number of shells of the computational 
C                     domain
C .......... DTOLD  = Current time-step                        (s) 
C .......... DTIME  = First guess for the new time-step        (s) 
C .......... V      = Specific volume (1/RHO)                  (cm^3 g^-1)
C                     [while V is evaluated at time tn,
C                     VV corresponds to time tn-1]
C .......... R      = Radius                                   (cm)
C                     [while R represents radius at time t_n,
C                     RR is the radius at time t_n-1]
C .......... EPS    = Maximum variation allowance for variables 
C                     R and V
C .......... AM2    = Auxiliar function used in the 
C                     determination of the new time-step
C .......... AM4    = Auxiliar function used in the 
C                     determination of the new time-step
C .......... DTMIN  = Auxiliar function used in the 
C                     determination of the new time-step
C *************************************************************************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION AM2(N),AM4(N),DTF(N)
      DIMENSION V(N),VV(N),R(N+1),RR(N+1)
C
      DTIME=DTOLD
C
C First guess of the new time-step based on the maximum allowed 
C variation of the physical variables V and R (U is not considered 
C since it dramatically slows down the calculation. For other 
C applications however, it may be wise to add an additional
C control based on the variation of the Lagrangian fluid velocity).
C From the current variation of a variable Y, DY/Dt, the maximum
C allowed variation of Y, EPS*Y, yields an estimate of the new 
C time-step
C
      DO 290 JK=1,N
      IF ((V(JK)-VV(JK)).NE.0.D0) THEN
       AM2(JK)=DABS(EPS*V(JK)*DTIME/(V(JK)-VV(JK)))
      ELSE
       AM2(JK)=1.D40
      ENDIF
  290 CONTINUE
C
      DO 490 JK=1,N
      IF ((R(JK+1)-RR(JK+1)).NE.0.D0) THEN
       AM4(JK)=DABS(EPS*R(JK+1)*DTIME/(R(JK+1)-RR(JK+1)))
      ELSE
       AM4(JK)=1.D40
      ENDIF
  490 CONTINUE
C
C Determination of the minimum time-step 
C
      DO 300 I=1,N
      DTF(I)=DMIN1(AM2(I),AM4(I))
  300 CONTINUE
C      
      DTMIN=DTF(1)
      DO 350 I=2,N
      IF (DTF(I).LT.DTMIN) DTMIN=DTF(I)
 350  CONTINUE
C
C DTMIN is the new time-step, provided that it does not
C exceed the previous time-step by more than 50% 
C
      DTM=1.5D0*DTIME
      IF (DTMIN.GT.DTM) THEN
      DTIME=DTM
      ELSE
      DTIME=DTMIN
      ENDIF
C
      RETURN
      END
C
C ***********************************************************************
      SUBROUTINE GAUSS(X,Z,Y,N)
C
C Solution of the matrix system of equations X(N,N) * Y(N) = Z(N) 
C through Gaussian elimination.
C
C .......... X      = Input matrix containing the partial derivatives 
C                     of F1, F2, and F5
C .......... Z      = Input matrix containing the values of F1, F2, and 
C                     F5
C .......... Y      = Output matrix of unknowns DW, DU, and DR
C .......... DFAC   = Rescaling factor for matrix X, per row
C ***********************************************************************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C      
      DIMENSION X(N,N),Z(N),DFAC(N),Y(N),IJJ(N)
C
C Determination of rescaling factors for matrix X
C
      DO 20 I=1,N
      SUP=0.D0
      DO 30 J=1,N
      SUP=DMAX1(SUP,DABS(X(I,J)))
   30 CONTINUE
      DFAC(I)=SUP
   20 CONTINUE
C
C Identification of the largest pivoting element. For each column,
C the largest rescaled value of X(i,j)/DFAC(i) is stored as DLPE and its
C location is registered by counter L
C
      DO 10 K=1,N
      IJJ(K)=K
   10 CONTINUE

      DO 40 J=1,N-1
      XPI=0.D0
C
      DO 50 I= J,N
      DLPE=DABS(X(IJJ(I),J))/DFAC(IJJ(I))
      IF (DLPE.GT.XPI) THEN
       XPI=DLPE
       L=I
      ELSE
      ENDIF
   50 CONTINUE
C
C Row exchange through IJJ counter
C
      IAUX=IJJ(J)
      IJJ(J)=IJJ(L)
      IJJ(L)=IAUX
C      
      DO 60 I=J+1,N
      XPO=X(IJJ(I),J)/X(IJJ(J),J)
      X(IJJ(I),J)=XPO
C
      DO 70 M=J+1,N
      X(IJJ(I),M)=X(IJJ(I),M)-XPO*X(IJJ(J),M)
   70 CONTINUE
   60 CONTINUE
C
   40 CONTINUE
C
      DO 80 I=1,N-1
      DO 90 J= I+1,N
      Z(IJJ(J))=Z(IJJ(J))-X(IJJ(J),I)*Z(IJJ(I))
   90 CONTINUE
   80 CONTINUE
C
      Y(N)=Z(IJJ(N))/X(IJJ(N),N)
      DO 100 I=N-1,1,-1
      Y(I)= Z(IJJ(I))
      DO 110 J=I+1,N
      Y(I)=Y(I)-X(IJJ(I),J)*Y(J)
  110 CONTINUE
      Y(I)=Y(I)/X(IJJ(I),I)
  100 CONTINUE
C
      RETURN
      END
C      
