```
c *********************************************************************
c This is --- MARSMAIN.FOR ---	    written by	   Harold A. Geller
c
c This main portion of the program calls other subroutines after
c displaying an initial menu, from which the user chooses what he
c desires in way of simulation or display of data
c This program has successfully been compiled on an IBM PC compatible
c using Microsoft FORTRAN
c
c *********************************************************************

program MARSMAIN

c Include PLOT.FD for the purposes of plotting

include    'plot.fd'

c      First, declare the variables and their values.  These will
c	be transferred to the various subroutines as required.

c      I will use the following parameters:
c	diameter of Mars,   d = 6794 kilometers
c	radius of Mars,     r = 3397 kilometers
c	mass of Mars,	    m = 6.578 x 10^23 kilograms
c	density of Mars,    rho = 3.9 grams per cubic centimeter
c	length of day,	    day = 8.864 x 10^4 seconds
c	obliquity,	    obl = 25.0 degrees
c	eccentricity,	    ecc = 0.093
c	solar orbit period, sop = 687 Earth days
c	scale height,	    H = 10
c	albedo, 	    albedo = 0.3
c       semi-major axis,    axis = 227.9 x 10^6 kilometers
c          albedocomp => 1-albedo
c    intercept of the linearized infrared cooling curve fit, A=211.1
c    slope of the linearized infrared cooling curve fit, B=1.55
c    the Martian solar constant must be calculated based on the
c     distance from the sun given the orbital angle and ellipse data
c      the solar constant for earth, Qe=334.4 Watts per square meter
c    Martian latitude,   xlat ==> expressed in degrees
c    temperature at a given latititude, Tx ==> expressed in degrees C
c    theta, the Mars/Sun angle of the orbit,   theta ==> in degrees
c    dist, the Mars/Sun calculated distance,   dist  ==> in kilometers

real*8 d,r,m,rho,day,obl,ecc,sop,H,albedo,albedocomp,Qe,Qm
real*8 axis,xlat,Tx,theta,dist,A,B,sofx,sinxlat,tQm

c       make subscripted variables to hold various values for plotting

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

d = 6.794e3
r = 3.397e3
m = 6.578e23
rho = 3.9
day = 8.864e4
obl = 25.0
ecc = 0.093
sop = 687.0
H = 10.0
albedo = 0.3
albedocomp = 1-albedo
Qe = 334.4
axis = 227.9e6

c       initialize things for the plotting

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   Before I display the menu, I must find the distance of
c    Mars from the Sun for the specific position of the planet
c     along the ellipse
c   To do this I call the MARSDIST subroutine, which calculates
c    the distance using Kepler's first law of planetary motion
c   This distance will then be used in the ensuing calculations
c   The Earth calculations assume a circular orbit.

5      write(*,*) 'Enter Sun/Mars ellipse angle'

c   Now calculate the distance
c    and the distance will be printed as well as the
c     solar radiation to be used in the other subroutines

c Now I calculate the Sun/Mars distance using Kepler's
c  first law where distance is semi-major axis times
c   1 minus the eccentricity squared divided by 1 plus
c    the eccentricity times the cosine of the Sun/Mars angle

dist = (axis*(1-ecc**2)) / (1+ecc*cos(theta*3.14159/180))

write(*,*) 'axis=',axis
write(*,*) 'ecc=',ecc
write(*,*) 'theta=',theta
write(*,*) 'cos(theta)=',cos(theta*3.14159/180)

c Now I calculate the solar constant for this distance

Qm = (Qe) / ((dist / 1.5e8)**2)

c Now print the values that will be used for the Qm and dist

write(*,*) 'The calculated Sun/Mars distance is:',dist
write(*,*) 'The derived solar constant at this distance is:',Qm

c   Now I set up the main menu

write(*,*) 'Enter the number corresponding to your choice'
write(*,*) ' '
write(*,*) '1) Display the physical characteristics'
write(*,*) '2) Display the current Earth temperature profile'
write(*,*) '3) Display the current Mars temperature profile'
write(*,*) '4) Display the Earth temps over 5 b.y. for 30% case'
write(*,*) '5) Display the Mars temps over 5 b.y. for 30% case'
write(*,*) '6) Display the Earth temps over 5 b.y. for 70% case'
write(*,*) '7) Display the Mars temps over 5 b.y. for 70% case'

c   Now I check to see what was chosen from the menu

IF (menu_choice .EQ. 1) THEN

c   Here I display all parameters

write(*,*) 'The diameter of Mars (kilometers):        ',d
write(*,*) 'The radius of Mars (kilometers):          ',r
write(*,*) 'The mass of Mars (kilograms):             ',m
write(*,*) 'The density of Mars (grams/cm^3):         ',rho
write(*,*) 'The length of a Mars day (seconds):       ',day
write(*,*) 'The obliquity of Mars (degrees):          ',obl
write(*,*) 'The eccentricity of Mars:                 ',ecc
write(*,*) 'The Mars solar orbit period (Earth days): ',sop
write(*,*) 'The scale height of Mars:                 ',H
write(*,*) 'The albedo of Mars:                       ',albedo
write(*,*) 'The semi-major axis of Mars (kilometers): ',axis

END IF

IF (menu_choice .EQ. 2) THEN
CALL MARS2
END IF

IF (menu_choice .EQ. 3) THEN

c       calculate the temperatures for the current senario

albedocomp = 0.7
A = 211.1
B = 1.55

c       set counter for number of points for plot.for

number = 0

c       in this program I set 1-albedo to 0.7 and it
c        remains at that value until the sin of the latitude
c         is greater than 0.95

DO 10 xlat= 0,90,1

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

IF (sinxlat .GT. 0.95) 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 sin of the latitude
c         is greater than 0.95

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 = ((Qm*sofx*albedocomp)-A)/B
ydata(91+xlat) = Tx
xdata(91+xlat) = xlat

IF (sinxlat .GT. 0.95) THEN
albedocomp = 0.4
END IF

number = number + 1

20      continue

c       set up the plot limits

xmin = - 90.0
xmax =   90.0
ymin = -130.0
ymax =  -30.0
x tic space = 20.0
y tic space = 10.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 = 'Mars Temp vs. Latitude for Specified M-S'
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()

END IF

IF (menu_choice .EQ. 4) THEN
CALL MARS4
END IF

IF (menu_choice .EQ. 5) THEN

c  reset equation parameters if not set already
albedocomp = 0.7
A = 211.1
B = 1.55

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 -78.5 degrees Celsius of the snow/ice line

DO 105 xlat= 0,90,1

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

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

number = number + 1

105      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 200 xlat= -1,-90,-1

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

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

number = number + 1

200      continue

c   NOW the plot for the case as it was 5 billion years ago
c    assuming the solar constant was 30% greater than now

c   FIRST I adjust the solar constant Q

tQm = 1.3*Qm

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 300 xlat= 0,90,1

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

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

number = number + 1

300      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 400 xlat= -1,-90,-1

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

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

number = number + 1

400      continue

c  Here's the case for 4 billion years ago

c   FIRST I adjust the solar constant Q

tQm = 1.24*Qm

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 500 xlat= 0,90,1

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

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

number = number + 1

500      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 600 xlat= -1,-90,-1

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

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

number = number + 1

600      continue

c  Here's the case for 3 billion years ago

c   FIRST I adjust the solar constant Q

tQm = 1.18*Qm

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 700 xlat= 0,90,1

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

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

number = number + 1

700      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 800 xlat= -1,-90,-1

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

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

number = number + 1

800      continue

c  Here's the case for 2 billion years ago

c   FIRST I adjust the solar constant Q

tQm = 1.12*Qm

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 900 xlat= 0,90,1

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

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

number = number + 1

900      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 1000 xlat= -1,-90,-1

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

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

number = number + 1

1000      continue

c  Here's the case for 1 billion years ago

c   FIRST I adjust the solar constant Q

tQm = 1.06*Qm

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 1100 xlat= 0,90,1

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

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

number = number + 1

1100      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 1200 xlat= -1,-90,-1

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

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

number = number + 1

1200      continue

c       set up the plot limits

xmin = - 90.0
xmax =   90.0
ymin = -130.0
ymax = 0.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 = 'Mars Temps Over 5 b.yr. (30% Case per M-S)'
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()
END IF

IF (menu_choice .EQ. 6) THEN
CALL MARS6
END IF

IF (menu_choice .EQ. 7) THEN
c  This is the case where the solar flux was 70% greater 5 b.yr. ago
c  reset equation parameters if not set already
albedocomp = 0.7
A = 211.1
B = 1.55

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 -78.5 degrees Celsius of the snow/ice line

DO 110 xlat= 0,90,1

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

IF (Tx .LT. -78.5) 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 210 xlat= -1,-90,-1

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

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

number = number + 1

210      continue

c   NOW the plot for the case as it was 5 billion years ago
c    assuming the solar constant was 70% greater than now

c   FIRST I adjust the solar constant Q

tQm = 1.7*Qm

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 310 xlat= 0,90,1

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

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

number = number + 1

310      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 410 xlat= -1,-90,-1

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

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

number = number + 1

410      continue

c  Here's the case for 4 billion years ago

c   FIRST I adjust the solar constant Q

tQm = 1.56*Qm

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 510 xlat= 0,90,1

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

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

number = number + 1

510      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 610 xlat= -1,-90,-1

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

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

number = number + 1

610      continue

c  Here's the case for 3 billion years ago

c   FIRST I adjust the solar constant Q

tQm = 1.42*Qm

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 710 xlat= 0,90,1

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

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

number = number + 1

710      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 810 xlat= -1,-90,-1

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

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

number = number + 1

810      continue

c  Here's the case for 2 billion years ago

c   FIRST I adjust the solar constant Q

tQm = 1.28*Qm

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 910 xlat= 0,90,1

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

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

number = number + 1

910      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 1010 xlat= -1,-90,-1

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

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

number = number + 1

1010      continue

c  Here's the case for 1 billion years ago

c   FIRST I adjust the solar constant Q

tQm = 1.14*Qm

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 1110 xlat= 0,90,1

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

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

number = number + 1

1110      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 1210 xlat= -1,-90,-1

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

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

number = number + 1

1210      continue

c       set up the plot limits

xmin = - 90.0
xmax =   90.0
ymin = -130.0
ymax = 40.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 = 'Mars Temps Over 5 b.yr. (70% Case per M-S)'
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()

END IF

write(*,*) 'Care to see more? (If yes enter 1)'