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

        read(*,*)

        call endgraphics()

        RETURN
        END

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