C***********************************************************************
       PROGRAM csdcalc2
C***********************************************************************
C Program to calculate a CSD.
C Written by......: Ronald G. Resmini
C Date............: May 28, 1992
C Modified........: December 30, 1994
C Place...........: Department of Earth and Planetary Sciences
C                   Johns Hopkins University
C                   Baltimore, Maryland 21218
C                   (410) 516-7034
C
C Modified........: December 30, 2000
C Place...........: 8604 Larkhaven Terrace, Fairfax Station, VA 22039
C                   (703) 690-4120
C***********************************************************************

       IMPLICIT NONE

       REAL         pctmelt,
     +              grthrate,
     +              crystsize(2590),
     +              crystpop(2590),
     +              currtemp,
     +              priortemp,
     +              liquidtemp,
     +              solidtemp,
     +              coolrate,
     +              spacing,
     +              timeinc,
     +              secshour,
     +              pi,
     +              nucleii,
     +              starttime,
     +              stoptime,
     +              timestep,
     +              z


       INTEGER      numiter,
     +              opiter,
     +              position,
     +              k


C*****************************************************************************
C                        Main processing statements:
C*****************************************************************************
C secshour...............seconds per hour, conversion factor
C crystsize()............array of crystal sizes
C crystpop().............array of crystal populations of each crystal size
C currtemp...............temperature at current time step (C)
C priortemp..............temperature at prior time step (C)
C liquidtemp.............liquidus temperature (C)
C solidtemp..............solidus temperature (C)
C coolrate...............cooling rate (C/hour)
C timeinc................hours
C grthrate...............crystal growth rate in meters/hour
C spacing................node edge length (meters)
C pctmelt................percent melt in node at current step
C numiter................number of program iterations - not currently used
C opiter.................number of iterations between outputs
C position...............distance node is from sill contact (meters)
C z......................main DO loop index
C k......................number of crystsize() and crystpop() elements
C starttime..............
C stoptime...............
C timestep...............time step in the temperature file

       secshour   =  3600.0
       pi         =  3.141592654
       grthrate   =  1.8e-08
       liquidtemp =  1140.0
       solidtemp  =   900.0
       priortemp  =  liquidtemp
       timeinc    =  2.00
       position   =  3.00
       k          =  2590
       spacing    =  0.25
       starttime  =  90.0
       stoptime   =  5210.0
       opiter     =  5210.0

       OPEN(7, FILE='c:\sill\node412\COOL412.T1',
     +         ACCESS='SEQUENTIAL',STATUS='OLD')

       OPEN(8, FILE='c:\sill\node412\GRTH412.TM1',
     +         ACCESS='SEQUENTIAL',STATUS='NEW')

       DO 1000 z = starttime, stoptime, timeinc

           READ(7,*) timestep, currtemp

           WRITE(*,*) timestep, currtemp

           coolrate = (priortemp - currtemp) / timeinc

         pctmelt = ((currtemp - solidtemp) / (liquidtemp - solidtemp))

           IF (coolrate .GT. 0.0) THEN
               crystpop(1) = nucleii(position, coolrate, timeinc)
           ENDIF
   
           crystsize(1) = 0.5e-06

           IF (coolrate .GT. 0.0) THEN
              CALL growsolids(crystsize, crystpop, pctmelt,
     +             timeinc, grthrate, spacing, k, z, pi)
              WRITE(8,*) z, char(9), coolrate, char(9), grthrate
           ENDIF

            priortemp = currtemp

             IF (MOD(z, 50) .EQ. 0) THEN
                 write(*,*) z
             ENDIF

            IF (MOD(z, opiter) .EQ. 0) THEN
                CALL outcryst(crystsize, crystpop, k)
            ENDIF

 1000  CONTINUE

       CLOSE(7)
       CLOSE(8)

       STOP
       END

C*****************************************************************************
       SUBROUTINE growsolids(crystsize,
     +                       crystpop,
     +                       pctofmelt,
     +                       timeinc,
     +                       grthrate,
     +                       spacing,
     +                       k,
     +                       z,
     +                       pi)

       IMPLICIT NONE

       INTEGER k,
     +         j

       REAL   epsilon,
     +        masspresent,
     +        massinnode,
     +        masstotallo,
     +        masstotalmid,
     +        masstotalhi,
     +        crystsize(k),
     +        crystpop(k),
     +        grthrate,
     +        pctofmelt,
     +        grthratelo,
     +        grthratemid,
     +        grthratehi,
     +        grthamthi,
     +        grthamtmid,
     +        grthamtlo,
     +        timeinc,
     +        density,
     +        pi,
     +        volfactor,
     +        spacing,
     +        z


       CHARACTER*1 shiftflag

C k...........corresponding crystpop() and crystsize() array addresses
C j...........an in-general loop index variable

       epsilon = 1.0e-01
       density = 2640000.0
       volfactor = (4.0/3.0) * pi * density
       shiftflag = 'N'

       IF (pctofmelt .LT. 0.00) THEN
           pctofmelt = 0.00
       END IF
       masspresent  = ((1.0 - pctofmelt) * (spacing**3) * density)
       grthratemid  = grthrate
       grthratelo   = grthratemid / 2.00
       grthratehi   = grthratemid * 2.00
       masstotalmid = 0.00

       j = 1
       massinnode = 0.00
       DO WHILE ((crystpop(j) .NE. 0.0) .AND. (j .LE. z+1))
          massinnode = massinnode + (crystsize(j)**3
     +                               * volfactor * crystpop(j))
          j = j + 1
       END DO
       
       DO WHILE ((ABS(masstotalmid - masspresent) .GE. epsilon) .AND.
     +          (massinnode .LT. masspresent) .AND.
     +          (ABS(massinnode - masspresent) .GE. epsilon))
       shiftflag = 'Y'

C          WRITE(*,*) 'z = ', z, ' k = ', k, ' growsolids '

          masstotalhi  = 0.00
          masstotallo  = 0.00
          masstotalmid = 0.00

          grthamthi = timeinc * grthratehi
          j = 1
          DO WHILE ((crystpop(j) .NE. 0.0) .AND. (j .LE. z+1))
            masstotalhi = masstotalhi + ((crystsize(j) + grthamthi)**3
     +       * volfactor * crystpop(j))
             j = j + 1
          END DO

          grthamtlo = timeinc * grthratelo
          j = 1
          DO WHILE ((crystpop(j) .NE. 0.0) .AND. (j .LE. z+1))
            masstotallo = masstotallo + ((crystsize(j) + grthamtlo)**3
     +       * volfactor * crystpop(j))
             j = j + 1
          END DO

          grthamtmid = timeinc * grthratemid
          j = 1
          DO WHILE ((crystpop(j) .NE. 0.0) .AND. (j .LE. z+1))
         masstotalmid = masstotalmid + ((crystsize(j) + grthamtmid)**3
     +       * volfactor * crystpop(j))
             j = j + 1
          END DO

          IF (masstotalhi .LT. masspresent) THEN
             grthratelo = grthratehi
             grthratehi = (2.00 * grthratehi)
             grthratemid  = ((grthratehi - grthratelo) / 2.0)
     +                                          + grthratelo
          ELSE IF (masstotallo .GT. masspresent) THEN
             grthratehi = grthratelo
             grthratelo = (grthratelo / 2.00)
             grthratemid  = ((grthratehi - grthratelo) / 2.0)
     +                                          + grthratelo
          ELSE IF ((masstotalhi .GE. masspresent) .AND.
     +            (masstotalmid .LT. masspresent)) THEN
             grthratelo = grthratemid
             grthratemid  = ((grthratehi - grthratelo) / 2.0)
     +                                          + grthratelo
          ELSE IF ((masstotallo .LE. masspresent) .AND.
     +            (masstotalmid .GT. masspresent)) THEN
             grthratehi = grthratemid
             grthratemid  = ((grthratehi - grthratelo) / 2.0)
     +                                          + grthratelo
          END IF

       END DO     

       IF (shiftflag .EQ. 'Y') THEN
           grthrate = grthratemid
           CALL shift(crystsize, crystpop, grthamtmid, k, z)
           shiftflag = 'N'
       ENDIF

       RETURN
       END

C*****************************************************************************
       SUBROUTINE shift(crystsize, 
     +                  crystpop,
     +                  grthamt,
     +                  k,
     +                  z)

       IMPLICIT NONE

       INTEGER j,
     +         k

       REAL   crystsize(k),
     +        crystpop(k),
     +        grthamt,
     +        z
 

C       WRITE(*,*) ' z = ', z, ' k = ', k, ' shift  ', ' grthamt = ', grthamt

       DO 1033 j = k, 2, -1
          IF (crystsize(j-1) .NE. 0.0) THEN
              crystpop(j) = crystpop(j-1)
              crystpop(j-1) = 0.00
              crystsize(j) = (crystsize(j-1) + grthamt)
              crystsize(j-1) = 0.00
          ENDIF
 1033  CONTINUE

       RETURN
       END

C*****************************************************************************
       REAL FUNCTION nucleii(position, coolrate, timeinc)

       IMPLICIT NONE

       INTEGER position

       REAL coolrate, timeinc, nucslope

C      nucslope = ((((position)) * (-2.20506e+07)) + 4.60385e+08)
C      nucleii = ((nucslope * coolrate) * timeinc)
C      nucleii = (4.2857e+09 * timeinc)

C      nucleii = (((1.0804e+08 * coolrate)
C    +     + (5.8078e+07 * coolrate * coolrate)) * timeinc)

C      nucleii = (4.2857e+07 * coolrate) * timeinc

C      nucleii = DNINT((4.2857e+13) * timeinc)
C      nucleii = DNINT(((4.170837d+15*coolrate*coolrate) - 
C    +           (6.280511d+11*coolrate) + (1.055224d+09))*timeinc)

       nucleii = (((3600.0d+0*((10.0d+0)**((1.37d+0*dlog10(coolrate))
     +                 + 0.47d+0)))*15625.0d+0)*timeinc)

       RETURN
       END

C*****************************************************************************
C                       File creation routines:
C*****************************************************************************
       SUBROUTINE outcryst(crystsize,
     +                     crystpop,
     +                     k)

       IMPLICIT NONE

       INTEGER k,
     +         i,
     +         j

       REAL   crystsize(k),
     +        crystpop(k)

 2061  FORMAT(1X, E17.9, 1A, E17.9)

        OPEN(9, FILE='c:\sill\node412\csd05p.TM1',
     +         ACCESS='SEQUENTIAL',STATUS='NEW',RECL=36,
     +         FORM='FORMATTED')

       DO 2042 j = 1, k, 1
        IF ((crystsize(j) .NE. 0.0) .AND. (crystpop(j) .NE. 0.0)) THEN
             WRITE(9,2061) crystsize(j), ',', crystpop(j)
        ENDIF
 2042  CONTINUE
       CLOSE(9)

       RETURN
       END

C****************************************************************************
