```
c#######################################################################
c       Subroutine MARS6.FOR
c                               by Harold Geller
c  The case where the solar flux was some 70% greater at
c    an epoch some 5 billion years ago
c The goal here is to determine the effect of the solar flux constant
c on the placement of the snow/ice line.
c In order to do this, I assumed that I must not fix the latitude
c of the occurence of the snow/ice line and allow it to drift
c with the temperature.
c Therefore, I allow the temperature to determine the latitude of
c the snow/ice line for all cases and here will plot the temperature
c based on a 70% increase in the solar constant over a period
c of 5 billion years.
c#######################################################################
SUBROUTINE MARS6

c       the following represents the variables
c       albedocomp => 1-albedo (0.7 except for snow/ice then its 0.4)
c       A => the intercept of the linearized infrared cooling (211.1)
c       B => the slope of the linearized infrared cooling (1.55)
c       Q =>= the solar flux given as 334.4 W/m^2
c       xlat => the latitude (north or south) on the earth's surface
c       Tx => the temperature at a given latitude

c       for plotting purposes I use a routine written by
c       Liam Gumley formerly of Curtin University of Technology
c       now with Research Data Systems Corp. of Greenbelt, MD

c       include definition file for logical*1 'init' variables
c       and integer*2 'color' variables

include    'plot.fd'

real*8     albedocomp, A, B, Q, xlat, Tx, sinxlat, sofx

c       make subscripted variables to hold various values for plotting

real*4     xdata(1086), ydata(1086)

c       initialize things for the first plot

init graphics = .true.
init border   = .true.
init window   = .true.
init axes     = .true.
init fonts    = .true.
plot method   = .true.
title only    = .false.
font height   = 10
font width    = 5

c       set up colors for just black and white
c        so that printing can be performed more clearly
c         on my printer

border color = 7
window color = 0
axis color   = 7
label color  = 7
plot color   = 7
title color  = 7

c       calculate the temperatures for the current senario

albedocomp = 0.7
A = 211.1
B = 1.55
Q = 334.4

c       set counter for number of points for plot.for

number = 0

c   in this part I set 1-albedo to 0.7 and it
c     remains at that value until the temperature of a latitude
c       is less than the -10 degrees Celsius of the snow/ice line

DO 10 xlat= 0,90,1

sinxlat = sin((ABS(xlat)*3.14159)/180)
sofx = 1-((0.241)*(3*sinxlat**2-1))
Tx = ((Q*sofx*albedocomp)-A)/B
ydata(91+xlat) = Tx
xdata(91+xlat) = xlat

IF (Tx .LT. -10.0) THEN
albedocomp = 0.4
END IF

number = number + 1

10      continue

c       Here I RESET 1-albedo to 0.7 and it
c        remains at that value until the temperature of the latitude
c         is less than the snow/ice line

albedocomp = 0.7

DO 20 xlat= -1,-90,-1

sinxlat = sin((ABS(xlat)*3.14159)/180)
sofx = 1-((0.241)*(3*sinxlat**2-1))
Tx = ((Q*sofx*albedocomp)-A)/B
ydata(91+xlat) = Tx
xdata(91+xlat) = xlat

IF (Tx .LT. -10.0) THEN
albedocomp = 0.4
END IF

number = number + 1

20      continue

c   NOW the plot for the case as it was 5 billion years ago
c    assuming the solar constant to be 70% that which it is now

c   FIRST I adjust the solar constant Q

Q = 1.7*Q

c   in this part I set 1-albedo to 0.7 and it
c     remains at that value until the temperature of a latitude
c       is less than the -10 degrees Celsius of the snow/ice line

albedocomp = 0.7

DO 30 xlat= 0,90,1

sinxlat = sin((ABS(xlat)*3.14159)/180)
sofx = 1-((0.241)*(3*sinxlat**2-1))
Tx = ((Q*sofx*albedocomp)-A)/B
ydata(272+xlat) = Tx
xdata(272+xlat) = xlat

IF (Tx .LT. -10.0) THEN
albedocomp = 0.4
END IF

number = number + 1

30      continue

c       Here I RESET 1-albedo to 0.7 and it
c        remains at that value until the temperature of the latitude
c         is less than the snow/ice line

albedocomp = 0.7

DO 40 xlat= -1,-90,-1

sinxlat = sin((ABS(xlat)*3.14159)/180)
sofx = 1-((0.241)*(3*sinxlat**2-1))
Tx = ((Q*sofx*albedocomp)-A)/B
ydata(272+xlat) = Tx
xdata(272+xlat) = xlat

IF (Tx .LT. -10.0) THEN
albedocomp = 0.4
END IF

number = number + 1

40      continue

c  Here's the case for 4 billion years ago

c   FIRST I adjust the solar constant Q

Q = 1.56*334.4

c   in this part I set 1-albedo to 0.7 and it
c     remains at that value until the temperature of a latitude
c       is less than the -10 degrees Celsius of the snow/ice line

albedocomp = 0.7

DO 50 xlat= 0,90,1

sinxlat = sin((ABS(xlat)*3.14159)/180)
sofx = 1-((0.241)*(3*sinxlat**2-1))
Tx = ((Q*sofx*albedocomp)-A)/B
ydata(453+xlat) = Tx
xdata(453+xlat) = xlat

IF (Tx .LT. -10.0) THEN
albedocomp = 0.4
END IF

number = number + 1

50      continue

c       Here I RESET 1-albedo to 0.7 and it
c        remains at that value until the temperature of the latitude
c         is less than the snow/ice line

albedocomp = 0.7

DO 60 xlat= -1,-90,-1

sinxlat = sin((ABS(xlat)*3.14159)/180)
sofx = 1-((0.241)*(3*sinxlat**2-1))
Tx = ((Q*sofx*albedocomp)-A)/B
ydata(453+xlat) = Tx
xdata(453+xlat) = xlat

IF (Tx .LT. -10.0) THEN
albedocomp = 0.4
END IF

number = number + 1

60      continue

c  Here's the case for 3 billion years ago

c   FIRST I adjust the solar constant Q

Q = 1.42*334.4

c   in this part I set 1-albedo to 0.7 and it
c     remains at that value until the temperature of a latitude
c       is less than the -10 degrees Celsius of the snow/ice line

albedocomp = 0.7

DO 70 xlat= 0,90,1

sinxlat = sin((ABS(xlat)*3.14159)/180)
sofx = 1-((0.241)*(3*sinxlat**2-1))
Tx = ((Q*sofx*albedocomp)-A)/B
ydata(634+xlat) = Tx
xdata(634+xlat) = xlat

IF (Tx .LT. -10.0) THEN
albedocomp = 0.4
END IF

number = number + 1

70      continue

c       Here I RESET 1-albedo to 0.7 and it
c        remains at that value until the temperature of the latitude
c         is less than the snow/ice line

albedocomp = 0.7

DO 80 xlat= -1,-90,-1

sinxlat = sin((ABS(xlat)*3.14159)/180)
sofx = 1-((0.241)*(3*sinxlat**2-1))
Tx = ((Q*sofx*albedocomp)-A)/B
ydata(634+xlat) = Tx
xdata(634+xlat) = xlat

IF (Tx .LT. -10.0) THEN
albedocomp = 0.4
END IF

number = number + 1

80      continue

c  Here's the case for 2 billion years ago

c   FIRST I adjust the solar constant Q

Q = 1.28*334.4

c   in this part I set 1-albedo to 0.7 and it
c     remains at that value until the temperature of a latitude
c       is less than the -10 degrees Celsius of the snow/ice line

albedocomp = 0.7

DO 90 xlat= 0,90,1

sinxlat = sin((ABS(xlat)*3.14159)/180)
sofx = 1-((0.241)*(3*sinxlat**2-1))
Tx = ((Q*sofx*albedocomp)-A)/B
ydata(815+xlat) = Tx
xdata(815+xlat) = xlat

IF (Tx .LT. -10.0) THEN
albedocomp = 0.4
END IF

number = number + 1

90      continue

c       Here I RESET 1-albedo to 0.7 and it
c        remains at that value until the temperature of the latitude
c         is less than the snow/ice line

albedocomp = 0.7

DO 100 xlat= -1,-90,-1

sinxlat = sin((ABS(xlat)*3.14159)/180)
sofx = 1-((0.241)*(3*sinxlat**2-1))
Tx = ((Q*sofx*albedocomp)-A)/B
ydata(815+xlat) = Tx
xdata(815+xlat) = xlat

IF (Tx .LT. -10.0) THEN
albedocomp = 0.4
END IF

number = number + 1

100      continue

c  Here's the case for 1 billion years ago

c   FIRST I adjust the solar constant Q

Q = 1.14*334.4

c   in this part I set 1-albedo to 0.7 and it
c     remains at that value until the temperature of a latitude
c       is less than the -10 degrees Celsius of the snow/ice line

albedocomp = 0.7

DO 110 xlat= 0,90,1

sinxlat = sin((ABS(xlat)*3.14159)/180)
sofx = 1-((0.241)*(3*sinxlat**2-1))
Tx = ((Q*sofx*albedocomp)-A)/B
ydata(996+xlat) = Tx
xdata(996+xlat) = xlat

IF (Tx .LT. -10.0) THEN
albedocomp = 0.4
END IF

number = number + 1

110      continue

c       Here I RESET 1-albedo to 0.7 and it
c        remains at that value until the temperature of the latitude
c         is less than the snow/ice line

albedocomp = 0.7

DO 120 xlat= -1,-90,-1

sinxlat = sin((ABS(xlat)*3.14159)/180)
sofx = 1-((0.241)*(3*sinxlat**2-1))
Tx = ((Q*sofx*albedocomp)-A)/B
ydata(996+xlat) = Tx
xdata(996+xlat) = xlat

IF (Tx .LT. -10.0) THEN
albedocomp = 0.4
END IF

number = number + 1

120      continue

c       set up the plot limits

xmin = - 90.0
xmax =   90.0
ymin = -130.0
ymax = 200.0
x tic space = 20.0
y tic space = 20.0

c       position this plot at the top left of the screen

x bottom left = 0.0
y bottom left = 0.0
x top right   = 1.0
y top right   = 1.0

c       give it a title

title = 'Temperature of Earth Over 5 Billion Years (70% Case)'
x axis title = 'Latitude (degrees)'
y axis title = 'Temperature (degree C)'

call            plot( init graphics, init border, init window,
&                  init axes, init fonts,
&                  xdata, ydata, number, plot method,
&                  invert x, invert y, xmin, xmax, ymin, ymax,
&                  x tic space, y tic space,
&                  x bottom left, y bottom left,
&                  x top right, y top right,
&                  border color, window color, axis color,
&                  label color, plot color, title color,
&                  x axis title, y axis title, title, title only,
&                  font height, font width )

c       wait for user to hit return

call endgraphics()

RETURN
END

c#######################################################################
c       END OF MARS6.FOR
c#######################################################################
```