This simulation is for present Greenland data.

 INITIAL CONDITIONS:
 TEMP AT 0    METERS =  -31.0000000000000       DEG C
 TEMP AT 10   METERS =  -30.7800000000000       DEG C
 TEMP AT 1000 METERS =  -9.00000000000000       DEG C

 ENTER NEW SURFACE TEMPERATURE
 -21

  Year  Day  Depth(m)  degC/10m
   249  358    880     0.00023     This is when this level
   250   13    870     0.00021     warms up to -3 deg C.
   250   75    890     0.00026
   250  179    860     0.00018
   250  231    900     0.00028
   251   79    910     0.00030
   251  205    850     0.00014
   251  336    920     0.00032
   252  262    930     0.00034
   253  215    940     0.00035
   253  259    840     0.00011
   254  191    950     0.00037
   255  185    960     0.00038
   256  192    970     0.00039
   257  211    980     0.00040
   258   83    830     0.00006
   258  238    990     0.00041
   259  272   1000     0.00041


* THIS PROGRAM MODELS A 1 METER SQUARE COLUMN OF ICE 1000 METERS HIGH
* THE COLUMN IS BROKEN UP INTO 10 METER LONG SEGMENTS.  THE DENSITY
* OF ICE IS ASSUMED TO BE 919 KGM PER CUBIC METER.  THE CONDUCTIVITY
* OF ICE IS ASSUMED TO BE 2.3 WATTS PER METER PER DEGREES CELSIUS.
* THE SPECIFIC HEAT OF ICE IS ASSUMED TO BE 2 KILOJOULES PER KILOGRAM
* PER DEGREES CELSIUS.  THE PURPOSE OF THIS PROGRAM IS TO DETERMINE
* THE EFFECT OF SURFACE TEMPERATURES AND THEIR EFFECT ON THE MELTING
* OF GREENLAND'S ICE.  THE CURRENT TEMPERATURE IS -31 DEG C AND HAS BEEN
* AT THAT TEMPERATURE FOR 10,000 YEARS.  BEFORE 10K YEARS AGO THE AVERAGE
* SURFACE TEMPERATURE WAS ABOUT -51 DEG C OVER GREENLAND.  THE INITIAL
* CONDITIONS FOR THIS PROGRAM ASSUME THE SURFACE IS -31 DEGREES C AND THE
* TEMPERATURE AT 1000 METERS IS -9 DEGREES C.  THIS IS CONSISTENT WITH
* THE THERMAL CONDUCTIVITY STATED ABOVE FOR A STEADY STATE CONDITION,
* WHICH SEEMS REASONABLE SINCE THE TIME CONSTANT FOR 1000 METERS OF ICE
* IS ONLY A FEW HUNDRED YEARS VERSUS THE 10,000 YEARS THE TEMPERATURE
* HAS BEEN AT -31 DEGREES C ON THE SURFACE.  NOW WE CAN MODEL THE EFFECT
* OF A TEMPERATURE RISE ON THE MELTING OF ICE AT THE BASE OF THE ICE.
*
      REAL*8 TEMP(100),T,Q,DT,SLOPE,B,TS
* SET TEMPERATURES TO SO -31 C AT 0 M AND -9 DEG C AT 1000 METERS
      SLOPE=(31.-9.)/1000.D0
      B=-31.D0
      DO 1 I=1,100
   1  TEMP(I) = SLOPE*I*10 + B
      WRITE(*,*) 'INITIAL CONDITIONS:'
      WRITE(*,*) 'TEMP AT 0    METERS =',B,        ' DEG C'
      WRITE(*,*) 'TEMP AT 10   METERS =',TEMP(1),  ' DEG C'
      WRITE(*,*) 'TEMP AT 1000 METERS =',TEMP(100),' DEG C'
      WRITE(*,*)
      WRITE(*,*) 'ENTER NEW SURFACE TEMPERATURE'
      READ(*,*) TS
* PRINT HEADING
      WRITE(*,*) ' Year  Day  Depth(m)  degC/10m'
* DEFINE THE YEARS LOOP AND RUN IT UNTIL THE BOTTOM IS AT 0 DEGREES C
      DO 2 IYEAR=1,1000
* DO THE CALCULATION EACH DAY SO THE TEMPERATURE CHANGES ARE VERY SMALL
      DO 3 IDAY=1,365
* USE LADDER CALCULATION METHOD FROM TOP NODE I TO BOTTOM NODE I
      T=TS
      DO 4 I=1,100
* CALCULATE THE 24 HOUR ENERGY THAT FLOWS THROUGH K TO HEAT UP THE BLOCK
* DELT TEMP IN DEGREES CELSIUS      2.3 WATTS/M/DEG C   1SQM/10M
* WATTS ARE JOULES PER SECOND AND THERE ARE 86400 SECONDS IN A DAY
      Q=(T-TEMP(I))*2.3/10.*24.*60.*60.  ! JOULES PER DAY INTO BLOCK I
* ACCOUNT FOR EARTH FLOW OF HEAT UPWARD
      Q=Q+4371.84
* RAISES THE TEMPERATURE OF NODE I BY DT DEGREES (2000 J / KGM / øC)
* THERE ARE 919 KGM/CUBIC METER AND 9190 KGM/10 CUBIC METERS
      DT=Q/(2000.*9190.)
      T=TEMP(I)+DT
      IF(TEMP(I).LT.-3..AND.T.GE.-3.) THEN
        WRITE(*,'(I6,I5,I7,F12.5)') IYEAR,IDAY,I*10,DT
        IF(I.EQ.100) STOP
      ENDIF
      TEMP(I)=T
   4  CONTINUE
   3  CONTINUE
      IDAY=365
   2  CONTINUE
      WRITE(*,*) 'YEAR 1000 TEMP AT 1000 M=',TEMP(100)
      STOP
      END