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'
       read(*,*) theta

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'

       read(*,*) menu_choice

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

        read(*,*)

        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

        read(*,*)

        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

        read(*,*)

        call endgraphics()

       END IF

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

       IF (inquire .EQ. 1) THEN
          GOTO 5
       END IF

       END