C***********************************************************************
       PROGRAM melt1d3
C***********************************************************************
C Program to calculate 1D cooling and solidification with true
C solidification fronts; not a Stefan model.
C***********************************************************************

       IMPLICIT NONE

       REAL         temold(960),
     +              tempnew(960),
     +              thcond(960),
     +              heat(960),
     +              thparm(960),
     +              pctmelt(960),
     +              heatmelt(960),
     +              top,
     +              bottom,
     +              spacing,
     +              timeinc,
     +              totltime,
     +              secshour,
     +              micrsecs

       CHARACTER*1  litholo(960),
     +              state(960),
     +              prtscrn

       INTEGER      row,
     +              numiter,
     +              opiter,
     +              z

       CHARACTER*45 temfile,
     +              heatfile,
     +              lithfile,
     +              statefile,
     +              pctmltfile


C*****************************************************************************
C                        Main processing statements:
C*****************************************************************************
C secshour...................seconds per hour
C micrsecs...................microseconds per second
C row........................INTEGER
C spacing....................meters
C timeinc....................hours
C totltime...................hours
C top........................temperature in degrees Celsius
C bottom.....................temperature in degrees Celsius

       secshour   =  3600.0
       micrsecs   =  1.0e-06
       row        =  960.0
       spacing    =  0.25
       timeinc    =  2.00
       totltime   =  5300.0
       top        =  100.0
       bottom     =  100.0
       temfile    =  'c:\eve\t1.dat'
       heatfile   =  'c:\eve\h1.dat'
       lithfile   =  'c:\eve\l1.dat'
       statefile  =  'c:\eve\state1.dat'
       pctmltfile =  'c:\eve\pctmlt1.dat'
       opiter     =   5300
       prtscrn    =  'n'

       CALL setparms(numiter, totltime, timeinc)
       CALL inlith(litholo, row, lithfile)
       CALL intemp(temold, row, temfile)
       CALL inheat(heat, row, heatfile, micrsecs, secshour)
       CALL instate(state, row, statefile)
       CALL inpctmlt(pctmelt, row, pctmltfile)

       OPEN(7, FILE='c:\scratch\COOL412.T1',
     +         ACCESS='SEQUENTIAL',STATUS='NEW')

       DO 1000 z = 1, numiter
           CALL updteparms(temold, thparm, thcond, litholo, state,
     +                     heatmelt, pctmelt, spacing, row)

           CALL finite(temold, tempnew, thcond,
     +                 thparm, heat, heatmelt, pctmelt, state,
     +                 litholo, top, bottom, timeinc,
     +                 spacing, row)

            CALL extemps(temold, tempnew, row)

            IF ((prtscrn .EQ. 'y') .OR. (prtscrn .EQ. 'Y')) THEN
                CALL prntemps(temold, row)
                PRINT *, 'End of iteration'
                PRINT *
            ENDIF

            IF (MOD(z, opiter) .EQ. 0) THEN
C               CALL outtemps(temold, row)
C               CALL outthcond(thcond, row)
C               CALL outpctmlt(pctmelt, row)
C               CALL outstate(state, row)
            ENDIF

            IF (MOD(z,1) .EQ. 0) THEN
                CALL printtemps(temold(412), z)
            ENDIF

             IF (MOD(z,10) .EQ. 0) THEN
                 WRITE(*,*) z
             ENDIF

C            CALL scrnstate(state, row, z)

 1000  CONTINUE

       CLOSE(7)

       STOP
       END

C*****************************************************************************
C             Establish the number of program iterations:
C*****************************************************************************
       SUBROUTINE setparms (numiter, totltime, timeinc)

            IMPLICIT NONE

            INTEGER numiter

            REAL   totltime, timeinc

            numiter = totltime/timeinc

       RETURN
       END

C*****************************************************************************
C          The bulk solid properties of the ten lithologies:
C*****************************************************************************
C NOTE: All lithology bulk and thermal properties are to produce values with
C       the following units: Density:              Kg/m3,
C                            Cp:                   J/Kg.C, and
C                            Thermal conductivity: J/m.hr.C
C                            Melt temperature:     C
C                            Solidus temperature:  C
C                            Melt heat:            J/Kg
C*****************************************************************************
       SUBROUTINE lithprops (temperature,
     +                       density,
     +                       Cp,
     +                       thcond,
     +                       litholo,
     +                       melttemp,
     +                       tsolid,
     +                       meltheat)

       IMPLICIT NONE

       REAL temperature,
     +        density,
     +        Cp,
     +        thcond,
     +        melttemp,
     +        tsolid,
     +        meltheat,
     +        t

       CHARACTER*1 litholo

            t = temperature + 273.15
            IF (litholo .EQ. 'A') THEN
                density = 2640.00
                Cp = 710.50
                thcond = 10260.00
                melttemp = 1140.0
                tsolid   =  900.0
                meltheat = 250000.0
            ELSEIF (litholo .EQ. 'B') THEN
                density = 2800.0
                Cp = 7646.127 - (1.075987*t) + (3711871.0/t) -
     +               (284465.7/(t**.5)) - (2.332948e+07/(t**2))
                thcond = 9.5125855e+07 - (170825.71 * t) +
     +               (7.29367210e+08/(t**.5)) + (100.528494 * (t**2))
                melttemp = 850.00
                tsolid   = 800.0
                meltheat = 205000.00
            ELSEIF (litholo .EQ. 'C') THEN
                density = 3030.0
                Cp = 2168.439 - (0.143958*t) + (583292.0/t) -
     +               (50970.02/(t**.5)) - (4.35851e+07/(t**2))
                thcond = 9.3257633e+07 - (154671.0 * t) +
     +               (6.50567e+08/(t**.5)) + (84.3217 * (t**2))
                melttemp = 850.00
                tsolid   = 800.0
                meltheat = 205000.00
            ELSEIF (litholo .EQ. 'D') THEN
                density = 2640.0
                Cp = 3595.382 - (0.3759056*t) + (1244499.0/t)-
     +               (104494.2/(t**.5)) - (8.036352e+07/(t**2))
                thcond = 8.647805e+07 - (155296.1 * t) +
     +               (6.630611e+08/(t**.5)) + (91.38945 * (t**2))
                melttemp = 850.00
                tsolid   = 800.00
                meltheat = 205000.00
            ELSEIF (litholo .EQ. 'E') THEN
                density = 2640
                 Cp = 710.50
                 thcond = 6.9587e+07
                melttemp = 850.00
                tsolid   = 800.00
                meltheat = 205000.00
            ELSEIF (litholo .EQ. 'F') THEN
                density = 2640
                Cp = 710.50
                thcond = 6.9587e+07
                melttemp = 850.00
                tsolid   = 800.00
                meltheat = 205000.00
            ELSEIF (litholo .EQ. 'G') THEN
                density = 2640
                Cp = 710.50
                thcond = 6.9587e+07
                melttemp = 850.00
                tsolid   = 800.00
                meltheat = 205000.00
            ELSEIF (litholo .EQ. 'H') THEN
                density = 2640
                Cp = 710.50
                thcond = 6.9587e+07
                melttemp = 850.00
                tsolid   = 800.00
                meltheat = 205000.00
            ELSEIF (litholo .EQ. 'I') THEN
                density = 2640
                Cp = 710.50
                thcond = 6.9587e+07
                melttemp = 850.00
                tsolid   = 800.00
                meltheat = 205000.00
            ELSE
                density = 2640
                Cp = 710.50
                thcond = 6.9587e+07
                melttemp = 850.00
                tsolid   = 800.00
                meltheat = 205000.00
            ENDIF

       RETURN
       END

C*****************************************************************************
C         Get the bulk solid properties of the ten lithologies:
C*****************************************************************************
       SUBROUTINE getsldprops (density,
     +                         Cp,
     +                         thcond,
     +                         melttemp,
     +                         tsolid,
     +                         meltheat,
     +                         temp,
     +                         lithtype)

            IMPLICIT NONE

            REAL   density,
     +               Cp,
     +               thcond,
     +               melttemp,
     +               tsolid,
     +               meltheat,
     +               temp

            CHARACTER*1 lithtype

            CALL lithprops(temp, density, Cp, thcond, lithtype,
     +                     melttemp, tsolid, meltheat)

       RETURN
       END

C*****************************************************************************
C          The bulk liquid properties of the ten lithologies:
C*****************************************************************************
C NOTE: All lithology bulk and thermal properties are to produce values with
C       the following units: Density:              Kg/m3,
C                            Cp:                   J/Kg.C, and
C                            Thermal conductivity: J/m.hr.C
C                            Melt temperature:     C
C                            Solidus temperature:  C
C*****************************************************************************
       SUBROUTINE liqdprops (temperature,
     +                       density,
     +                       Cp,
     +                       thcond,
     +                       litholo,
     +                       melttemp,
     +                       tsolid)

       IMPLICIT NONE

       REAL   temperature,
     +          density,
     +          Cp,
     +          thcond,
     +          melttemp,
     +          tsolid,
     +          t

       CHARACTER*1 litholo

            t = temperature + 273.15
            IF (litholo .EQ. 'A') THEN
                density = 2640.00
                Cp = 1752.1667
                thcond = 10260.00
                melttemp = 1140.0
                tsolid   =  900.0
            ELSEIF (litholo .EQ. 'B') THEN
                density = 2800.0
                Cp = 7646.127 - (1.075987*t) + (3711871.0/t) -
     +               (284465.7/(t**.5)) - (2.332948e+07/(t**2))
                thcond = 9.5125855e+07 - (170825.71 * t) +
     +               (7.29367210e+08/(t**.5)) + (100.528494 * (t**2))
                melttemp = 850.00
                tsolid   = 800.00
            ELSEIF (litholo .EQ. 'C') THEN
                density = 3030.0
                Cp = 2168.439 - (0.143958*t) + (583292.0/t) -
     +               (50970.02/(t**.5)) - (4.35851e+07/(t**2))
                thcond = 9.3257633e+07 - (154671.0 * t) +
     +               (6.50567e+08/(t**.5)) + (84.3217 * (t**2))
                melttemp = 850.00
                tsolid   = 800.00
            ELSEIF (litholo .EQ. 'D') THEN
                density = 2640.0
                Cp = 3595.382 - (0.3759056*t) + (1244499.0/t)-
     +               (104494.2/(t**.5)) - (8.036352e+07/(t**2))
                thcond = 8.647805e+07 - (155296.1 * t) +
     +               (6.630611e+08/(t**.5)) + (91.38945 * (t**2))
                melttemp = 850.00
                tsolid   = 800.00
            ELSEIF (litholo .EQ. 'E') THEN
                density = 2640
                Cp = 710.50
                thcond = 6.9587e+07
                melttemp = 850.00
                tsolid   = 800.00
            ELSEIF (litholo .EQ. 'F') THEN
                density = 2640
                Cp = 710.50
                thcond = 6.9587e+07
                melttemp = 850.00
                tsolid   = 800.00
            ELSEIF (litholo .EQ. 'G') THEN
                density = 2640
                Cp = 710.50
                thcond = 6.9587e+07
                melttemp = 850.00
                tsolid   = 800.00
            ELSEIF (litholo .EQ. 'H') THEN
                density = 2640
                Cp = 710.50
                thcond = 6.9587e+07
                melttemp = 850.00
                tsolid   = 800.00
            ELSEIF (litholo .EQ. 'I') THEN
                density = 2640
                Cp = 710.50
                thcond = 6.9587e+07
                melttemp = 850.00
                tsolid   = 800.00
            ELSE
                density = 2640
                Cp = 710.50
                thcond = 6.9587e+07
                melttemp = 850.00
                tsolid   = 800.00
            ENDIF

       RETURN
       END

C*****************************************************************************
C         Get the bulk liquid properties of the ten lithologies:
C*****************************************************************************
       SUBROUTINE getliqprops (density,
     +                         Cp,
     +                         thcond,
     +                         melttemp,
     +                         tsolid,
     +                         temp,
     +                         lithtype)

            IMPLICIT NONE

            REAL   density,
     +               Cp,
     +               thcond,
     +               melttemp,
     +               tsolid,
     +               temp

            CHARACTER*1 lithtype

            CALL liqdprops(temp, density, Cp, thcond,
     +                     lithtype, melttemp, tsolid)

       RETURN
       END

C*****************************************************************************
C Procedure to calculate the temperature dependent nature of the grid
C lithologies:
C*****************************************************************************
       SUBROUTINE updteparms (temold,
     +                        thparm,
     +                        thcond,
     +                        litholo,
     +                        state,
     +                        heatmelt,
     +                        pctmelt,
     +                        spacing,
     +                        row)

       IMPLICIT NONE

       INTEGER row,
     +         i,
     +         j

       REAL  Cp,
     +         density,
     +         thermcond,
     +         spacing,
     +         thparm(row),
     +         thcond(row),
     +         temold(row),
     +         heatmelt(row),
     +         pctmelt(row),
     +         melttemp,
     +         tsolid,
     +         meltheat,
     +         amount

       CHARACTER*1 litholo(row),
     +             state(row)

       DO 1030 i = 1, row
          IF (state(i) .EQ. 'S') THEN
              CALL getsldprops(density, Cp, thermcond,
     +              melttemp, tsolid, meltheat, temold(i), litholo(i))
              amount = density * spacing * spacing * spacing
              thparm(i) = (1.0 / (Cp * amount))
              thcond(i) = thermcond
              heatmelt(i) = (tsolid - temold(i)) * Cp * amount
          ELSEIF (state(i) .EQ. 'M') THEN
              CALL getliqprops(density, Cp, thermcond,
     +              melttemp, tsolid, temold(i), litholo(i))
              amount = density * spacing * spacing * spacing
              thparm(i) = (1.0 / (Cp * amount))
              thcond(i) = thermcond
              CALL getsldprops(density, Cp, thermcond,
     +              melttemp, tsolid, meltheat, temold(i), litholo(i))
              heatmelt(i) = (((1.0 - pctmelt(i)) * amount) * meltheat)
          ENDIF
 1030  CONTINUE

       RETURN
       END

C*****************************************************************************
C              Calculation of volumetric heat production:
C*****************************************************************************
       FUNCTION heatprod (heat,
     +                    timeinc,
     +                    spacing)

       IMPLICIT NONE

       REAL heat,
     +        timeinc,
     +        spacing,
     +        heatprod

       heatprod = heat * timeinc * spacing * spacing * spacing

       RETURN
       END

C****************************************************************************
       SUBROUTINE extemps(temold,
     +                    tempnew,
     +                    row)

       IMPLICIT NONE

       INTEGER row,
     +         i,
     +         j

       REAL temold(row),
     +        tempnew(row)

       DO 1050 i = 1, row
          temold(i) = tempnew(i)
 1050  CONTINUE

       RETURN
       END

C*****************************************************************************
C                   Grid print-to-screen routine:
C*****************************************************************************
       SUBROUTINE prntemps(temold,
     +                     row)

       IMPLICIT NONE

       INTEGER  i,
     +          j,
     +          row

       REAL     temold(row)

 1100  FORMAT()
       DO 1080 i = 1, row
           WRITE(*, *) temold(i)
 1080  CONTINUE

       RETURN
       END

C*****************************************************************************
       SUBROUTINE prntmelt(pctmelt,
     +                     row)

       IMPLICIT NONE

       INTEGER  i,
     +          j,
     +          row

       REAL     pctmelt(row)

       DO 1087 i = 1, row
          WRITE(*, *) pctmelt(i)
 1087  CONTINUE

       RETURN
       END

C*****************************************************************************
       SUBROUTINE prntstate(state,
     +                      row)

       IMPLICIT NONE

       INTEGER  i,
     +          j,
     +          row

       CHARACTER*1  state(row)

       DO 1089 i = 1, row
           WRITE(*, *) state(i)
 1089  CONTINUE

       RETURN
       END

C*****************************************************************************
C Subroutine to calculate the new temperatures of the nodes:
C*****************************************************************************
       SUBROUTINE newtemp (jle, jletomelt, spacing, percentmelt,
     +                     tempnew, tempold, thermparm,
     +                     phase, lithtype)

         IMPLICIT NONE

         REAL   jle,
     +          jletomelt,
     +          density,
     +          Cp,
     +          thermcond,
     +          melttemp,
     +          tsolid,
     +          spacing,
     +          percentmelt,
     +          meltjle,
     +          tempnew,
     +          tempold,
     +          amount,
     +          amtmelted,
     +          thermparm,
     +          heatdiff

         CHARACTER*1 phase,
     +               lithtype

         IF (phase .EQ. 'S') THEN
             CALL getsldprops(density, Cp, thermcond, melttemp,
     +                        tsolid, meltjle, tempold, lithtype)
             amount = (density * spacing * spacing * spacing)
             thermparm = (1.0 / (Cp * amount))
             tempnew = (jle * thermparm) + tempold
         ELSEIF (phase .EQ. 'M') THEN
             CALL getliqprops(density, Cp, thermcond, melttemp,
     +                        tsolid, tempold, lithtype)
             amount = density * spacing * spacing * spacing
             thermparm = (1.0 / (Cp * amount))
             IF (((tsolid - tempold) * Cp * amount) .GT. jle) THEN
                heatdiff = jle - ((tsolid - tempold) * Cp * amount)
                CALL getsldprops(density, Cp, thermcond, melttemp,
     +                           tsolid, meltjle, tempold, lithtype)
                amount = (density * spacing * spacing * spacing)
                thermparm = (1.0 / (Cp * amount))
                tempnew = (heatdiff * thermparm) + tsolid
                percentmelt = 0.0
                phase = 'S'
             ELSE
                IF (((tsolid - tempold) * Cp * amount) .EQ. jle) THEN
                   tempnew = tsolid
                   percentmelt = 0.0
                   phase = 'S'
                ELSE
                   CALL getliqprops(density, Cp, thermcond, melttemp,
     +                              tsolid, tempold, lithtype)
                   amount = density * spacing * spacing * spacing
                   thermparm = (1.0 / (Cp * amount))
                   tempnew = (jle * thermparm) + tempold
         percentmelt = ((tempnew - melttemp)/(melttemp - tsolid)) + 1.0
                   phase = 'M'
                ENDIF
             ENDIF
         ENDIF

       RETURN
       END

C*****************************************************************************
C           The finite difference calculation routines:
C*****************************************************************************
       SUBROUTINE finite (temold,
     +                    tempnew,
     +                    thcond,
     +                    thparm,
     +                    heat,
     +                    heatmelt,
     +                    pctmelt,
     +                    state,
     +                    litholo,
     +                    lefttemp,
     +                    rghttemp,
     +                    timeinc,
     +                    spacing,
     +                    row)

       IMPLICIT NONE

       INTEGER i,
     +         row

       REAL  joules,
     +         temold(row),
     +         tempnew(row),
     +         thcond(row),
     +         thparm(row),
     +         heat(row),
     +         heatmelt(row),
     +         pctmelt(row),
     +         lefttemp,
     +         rghttemp,
     +         timeinc,
     +         spacing,
     +         heatprod

       CHARACTER*1  litholo(row),
     +              state(row)

            joules = ((lefttemp - temold(1)) *
     +            (thcond(1) * spacing * timeinc)) +
     +            ((temold(2) - temold(1)) * spacing * timeinc *
     +            ((thcond(2) + thcond(1)) / 2.0)) 
C    +         +   heatprod(heat(1), timeinc, spacing)
            CALL  newtemp(joules, heatmelt(1), spacing, pctmelt(1),
     +                    tempnew(1), temold(1), thparm(1),
     +                    state(1), litholo(1))

            DO 1500 i = 2, (row - 1)
                   joules =  ((temold(i-1) - temold(i)) *
     +                      spacing * timeinc *
     +                    ((thcond(i-1) + thcond(i)) / 2.0)) +
     +                    ((temold(i+1) - temold(i)) *
     +                      spacing * timeinc *
     +                    ((thcond(i+1) + thcond(i)) / 2.0)) 
C    +         +   heatprod(heat(i), timeinc, spacing)
                  CALL newtemp(joules, heatmelt(i), spacing, pctmelt(i),
     +                         tempnew(i), temold(i), thparm(i),
     +                         state(i), litholo(i))
 1500       CONTINUE

            joules = ((rghttemp - temold(row)) *
     +             (thcond(row) * spacing * timeinc)) +
     +            ((temold(row-1) - temold(row)) *
     +              spacing * timeinc *
     +            ((thcond(row-1) + thcond(row)) / 2.0))
C    +    +  heatprod(heat(row), timeinc, spacing)
            CALL newtemp(joules, heatmelt(row), spacing, pctmelt(row),
     +                   tempnew(row), temold(row), thparm(row),
     +                   state(row), litholo(row))

       RETURN
       END

C*****************************************************************************
C               File input of program parameters:
C*****************************************************************************
       SUBROUTINE intemp (temold,
     +                    row,
     +                    temfile)

       IMPLICIT NONE

       INTEGER i,
     +         j,
     +         row

       REAL    temold(row)

       CHARACTER*45 temfile

       OPEN(7, FILE=temfile, ACCESS='SEQUENTIAL',STATUS='OLD')

 1150  FORMAT(F7.2)
       DO 1120 i = 1, row
          READ(7,*) temold(i)
 1120  CONTINUE

       CLOSE(7)

       RETURN
       END

C*****************************************************************************
       SUBROUTINE inpctmlt (pctmelt,
     +                      row,
     +                      pctmltfile)

       IMPLICIT NONE

       INTEGER i,
     +         j,
     +         row

       REAL    pctmelt(row)

       CHARACTER*45 pctmltfile

       OPEN(7, FILE=pctmltfile, ACCESS='SEQUENTIAL',STATUS='OLD')

 1152  FORMAT(F7.2)
       DO 1122 i = 1, row
          READ(7,*) pctmelt(i)
 1122  CONTINUE

       CLOSE(7)

       RETURN
       END

C*****************************************************************************
       SUBROUTINE inheat (heat,
     +                    row,
     +                    heatfile,
     +                    micrsecs,
     +                    secshour)

       IMPLICIT NONE

       INTEGER i,
     +         j,
     +         row

       REAL    heat(row),
     +           micrsecs,
     +           secshour

       CHARACTER*45 heatfile

       OPEN(7, FILE=heatfile, ACCESS='SEQUENTIAL',STATUS='OLD')

 1250  FORMAT(F7.2)
       DO 1220 i = 1, row
          READ(7,*) heat(i)
          heat(i) = (heat(i) * micrsecs * secshour)
 1220  CONTINUE

       CLOSE(7)

       RETURN
       END

C*****************************************************************************
       SUBROUTINE inlith (litholo,
     +                    row,
     +                    lithfile)

       IMPLICIT NONE

       INTEGER i,
     +         j,
     +         row

       CHARACTER*1  litholo(row)

       CHARACTER*45 lithfile

       OPEN(7, FILE=lithfile, ACCESS='SEQUENTIAL',STATUS='OLD')

 1350  FORMAT(A)
       DO 1320 i = 1, row
          READ(7,1350) litholo(i)
 1320  CONTINUE

       CLOSE(7)

       RETURN
       END

C*****************************************************************************
       SUBROUTINE instate (state,
     +                     row,
     +                     statefile)

       IMPLICIT NONE

       INTEGER i,
     +         j,
     +         row

       CHARACTER*1  state(row)

       CHARACTER*45 statefile

       OPEN(7, FILE=statefile, ACCESS='SEQUENTIAL',STATUS='OLD')

 1352  FORMAT(A)
       DO 1322 i = 1, row
          READ(7,1352) state(i)
 1322  CONTINUE

       CLOSE(7)

       RETURN
       END

C*****************************************************************************
C                       File creation routines:
C*****************************************************************************
       SUBROUTINE outtemps (temold,
     +                      row)

       IMPLICIT NONE

       INTEGER i,
     +         j,
     +         row

       REAL    temold(row)

       OPEN(7, FILE='c:\eve\TEMOUT.TMP',
     +      ACCESS='SEQUENTIAL',STATUS='NEW')

 2020  FORMAT(F8.2)
       DO 1800 i = 1, row
          WRITE(7, 2020) temold(i)
 1800  CONTINUE

       CLOSE(7)

       RETURN
       END

C*****************************************************************************
       SUBROUTINE printtemps(temperature,
     +                       z)

       IMPLICIT NONE

       INTEGER z

       REAL    temperature

C 2023  FORMAT(I6, 2A, F7.2) Old format stmt.

 2023  FORMAT(I6, 2A, F10.5)
          WRITE(7, 2023) z, char(9), ',', temperature

       RETURN
       END

C*****************************************************************************
       SUBROUTINE outthcond(thcond,
     +                      row)

       IMPLICIT NONE

       INTEGER i,
     +         j,
     +         row

       REAL    thcond(row)

       OPEN(7, FILE='c:\eve\CONOUT.TMP',
     +      ACCESS='SEQUENTIAL',STATUS='NEW')

 2200  FORMAT(F6.2)
       DO 2800 i = 1, row
          WRITE(7, 2200) (thcond(i) / 3600)
 2800  CONTINUE

       CLOSE(7)

       RETURN
       END

C*****************************************************************************
       SUBROUTINE outpctmlt (pctmelt,
     +                       row)

       IMPLICIT NONE

       INTEGER i,
     +         j,
     +         row

       REAL    pctmelt(row)

       OPEN(7, FILE='c:\eve\PCTMELT.OUT',
     +      ACCESS='SEQUENTIAL',STATUS='NEW')

 2010  FORMAT(F6.2)
       DO 1803 i = 1, row
          WRITE(7, 2010) pctmelt(i)
 1803  CONTINUE

       CLOSE(7)

       RETURN
       END

C*****************************************************************************
       SUBROUTINE outstate (state,
     +                      row)

       IMPLICIT NONE

       INTEGER i,
     +         j,
     +         row

       CHARACTER*1    state(row)

       OPEN(7, FILE='c:\eve\STATE.TMP',
     +      ACCESS='SEQUENTIAL',STATUS='NEW')

 2000  FORMAT(A)
       DO 1805 i = 1, row
          WRITE(7, 2000) state(i)
 1805  CONTINUE

       CLOSE(7)

       RETURN
       END

C*****************************************************************************
       SUBROUTINE scrnstate(state,
     +                     row,
     +                     z)

       IMPLICIT NONE

       INTEGER  i,
     +          z,
     +          row

       CHARACTER*1  state(row)

       IF (state(35) .EQ. 'S') THEN 
           WRITE(*,*) z
       ENDIF

       RETURN
       END
C*****************************************************************************

