C
      SUBROUTINE CNSTNT
      implicit none
C     ****
C     ****     DEFINES AND CALCULATES CONSTANTS NEEDED BY PROGRAM
C     ****
C     ****     PARAMETERS DEFINE GEOGRAPHIC AND GEOMAGNETIC GRIDS
C     ****
      include "params.h"
      include "consts.h"
C     ****
C     ****    NOW SET CONSTANTS
C     ****
      real :: tabl2,tanth0,tanths,theta0,hamh0
      COMMON/BLNK/TABL2(91,3:5),TANTH0(JMAXM),TANTHS(JMAXM),
     1  THETA0(JMAXM),HAMH0(JMAXM)
!
! Local:
      integer :: i,j,n
      real :: e,alfa,dtheta,tanths2
!
      DATA E/1.E-6/
      ALFA = 1.668
      R0 = RE+H0
      R00 = RE+H00
      RS = RE+HS
      PI = 4.*ATAN(1.)
      DTR = PI/180.
      RTD = 180./PI
      DLATG = PI/FLOAT(JMAXG)
      DLONG = 2.*PI/FLOAT(IMAXG)
      DLATM = PI/FLOAT(JMAXM-1)
      DLONM = 2.*PI/FLOAT(IMAXM)
C     ****
C     ****     FILL ARRAY YLATG
C     ****
      DO 1 J = 1,JMAXG
        YLATG(J)=-.5*(PI-DLATG)+FLOAT(J-1)*DLATG
    1 CONTINUE
      YLATG(0) = -PI/2.+E
      YLATG(JMAXGP) = PI/2.-E
C     ****
C     ****     FILL YLONG
C     ****
      DO 2 I = 1,IMAXGP
        YLONG(I) = -PI+FLOAT(I-1)*DLONG
    2 CONTINUE
C     ****
C     ****     FILL ARRAY YLATM (EQUALLY SPACED IN THETAO BUT HOLDS
C     ****       CORRESPONDING VALUE OF THETAS)
C     ****
      DO 3 J = 1,JMAXM
C       ****    THETA0 = EQUALLY SPACED GRID VALUES
        THETA0(J) = -PI/2.+FLOAT(J-1)*DLATM
    3 CONTINUE
      DO 4 J = 2,JMAXM-1
C       ****     TANTH0 = ABS(TAN(THETA0))
        TANTH0(J) = ABS(TAN(THETA0(J)))
C       ****     HAMH0 = HA-H0
        HAMH0(J) = R1*TANTH0(J)+R0*TANTH0(J)**(2.+2.*ALFA)/
     1    (1.+TANTH0(J)**2)**ALFA
C       ****     TANTHS = ABS(TAN(THETAS))
        TANTHS(J) = SQRT(HAMH0(J)/R0)
C       ****     YLATM = TANTHS
        YLATM(J) = SIGN(ATAN(TANTHS(J)),THETA0(J))
C       ****     RCOS0S = COS(THETA0)/COS(THETAS)
        RCOS0S(J) = SQRT((1.+TANTHS(J)**2)/(1.+TANTH0(J)**2))
C       ****     DT0DTS = D(THETA0)/D(THETAS)
c99/2/12b
C old   DT0DTS(J) = (2.*R0*TANTHS(J)*(1.+TANTHS(J)**2))/
      tanths2 = TANTHS(J)**2
C DT1DTS = DT0DTS divided by abs[sin(Im)].  Remains finite
C   and nonzero at magnetic equator.
      DT1DTS(J) = (R0*sqrt(1.+4.*tanths2)*(1.+tanths2))/
     1    (R1*(1.+TANTH0(J)**2)+2.*R0*TANTH0(J)**(2.*ALFA+1.)*
     2    (1.+ALFA+TANTH0(J)**2)/(1.+TANTH0(J)**2)**ALFA)
      DT0DTS(J) = DT1DTS(J)*2.*TANTHS(J)/sqrt(1.+4.*tanths2)
c99/2/12e
    4 CONTINUE
C     ****
C     ****     NOW DO POLES
C     ****
      YLATM(1) = THETA0(1)
      YLATM(JMAXM) = THETA0(JMAXM)
      RCOS0S(1) = 1.
      RCOS0S(JMAXM) = 1.
      DT0DTS(1) = 1.
      DT0DTS(JMAXM) = 1.
c99/2/12b
      DT1DTS(1) = 1.
      DT1DTS(JMAXM) = 1.
c99/2/12e
C     ****
C     ****     FILL YLONM
C     ****
      DO 5 I = 1,IMAXMP
        YLONM(I) = -PI+FLOAT(I-1)*DLONM
    5 CONTINUE
      DTHETA = PI/(2.*90.)
      TABLE(1,1) = 0.
      TABLE(1,2) = 0.
      DO 6 I = 2,91
        TABLE(I,1) = TABLE(I-1,1)+DTHETA
    6 CONTINUE
      DO 7 I = 2,90
        TABL2(I,4) = TAN(TABLE(I,1))
        TABLE(I,2) = TABLE(I,1)
    7 CONTINUE
      DO 8 N = 1,7
        DO 9 I = 2,90
          TABL2(I,3) = TABLE(I,2)
          TABLE(I,2) = TAN(TABL2(I,3))
          TABL2(I,5) = SQRT(R1/R0*TABLE(I,2)+TABLE(I,2)**
     1      (2.*(1.+ALFA))/(1.+TABLE(I,2)**2)**ALFA)
          TABLE(I,2) = TABL2(I,3)-(TABL2(I,5)-TABL2(I,4))*2.*TABL2(I,5)/
     2      (R1/R0*(1.+TABLE(I,2)**2)+2.*TABLE(I,2)**(2.*ALFA+1.)*
     3      (1.+ALFA+TABLE(I,2)**2)/(1.+TABLE(I,2)**2)**ALFA)
    9   CONTINUE
    8 CONTINUE
      RETURN
      END
C