c#######################################################################
c This is subroutine MARS2.FOR  written by Harold Geller
c
c This subroutine displays the Earth temperatures over latitude
c based an old homework assignment, 
c where the albedo changes to the ice albedo (or rather the
c 1-albedo value) when sin(latitude) is 0.95 or greater.
c This is used as a comparison to the values found for Mars
c under calculation in the third part of the MARSMAIN menu.
c#######################################################################
        SUBROUTINE MARS2

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(181), ydata(181)

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 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 = ((Q*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 = ((Q*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 = -100.0
        ymax =  60.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 at given Latitude'
        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 SUBROUTINE MARS2
c######################################################################