      SUBROUTINE TABLEX                                                 TABX001 

C        CALCULATES ELECTRON DENSITY AND GRADIENT FROM PROFILES HAVING  TABX002 

C        THE SAME FORM AS THOSE USED BY CROFTS RAY TRACING PROGRAM      TABX003 

C MAKES AN EXPONENTIAL EXTRAPOLATION DOWN USING THE BOTTOM TWO POINTS   TABX004 

C        NEEDS SUBROUTINE GAUSEL                                        TABX005 

      DIMENSION HPC(250),FN2C(250),ALPHA(250),BETA(250),GAMMA(250),     TABX006 

     1 DELTA(250),SLOPE(250),MAT(4,5)                                   TABX007 

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RAD,K,DUM(2)                     TABX008 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  TABX009 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 TABX010 

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(READFN,W(100)),(PERT,W(150))  TABX011 

      REAL MAT,K                                                        TABX012 

      DATA (MODX(1)=6HTABLEX)                                           TABX013 

      ENTRY ELECTX                                                      TABX014 

      IF (READFN.EQ.0.) GO TO 10                                        TABX015 

      READFN=0.                                                         TABX016 

      READ 1000, NOC,(HPC(I),FN2C(I),I=1,NOC)                           TABX017 

 1000 FORMAT (I4/(F8.2,E12.4))                                          TABX018 

      PRINT 1200, (HPC(I),FN2C(I),I=1,NOC)                              TABX019 

 1200 FORMAT(1H1,14X,6HHEIGHT,4X,16HELECTRON DENSITY/(1X,F20.10,E20.10))TABX020 

      A=0.                                                              TABX021 

      IF(FN2C(1).NE.0.) A=ALOG(FN2C(2)/FN2C(1))/(HPC(2)-HPC(1))         TABX022 

      FN2C(1)=K*FN2C(1)                                                 TABX023 

      FN2C(2)=K*FN2C(2)                                                 TABX024 

      SLOPE(1)=A*FN2C(1)                                                TABX025 

      SLOPE(NOC)=0.                                                     TABX026 

      NMAX=1                                                            TABX027 

      DO 6 I=2,NOC                                                      TABX028 

      IF (FN2C(I).GT.FN2C(NMAX)) NMAX=I                                 TABX029 

      IF (I.EQ.NOC) GO TO 4                                             TABX030 

      FN2C(I+1)=K*FN2C(I+1)                                             TABX031 

      DO 3 J=1,3                                                        TABX032 

      M=I+J-2                                                           TABX033 

      MAT(J,1)=1.                                                       TABX034 

      MAT(J,2)=HPC(M)                                                   TABX035 

      MAT(J,3)=HPC(M)**2                                                TABX036 

  3   MAT(J,4)=FN2C(M)                                                  TABX037 

      CALL GAUSEL (MAT,4,3,4,NRANK)                                     TABX038 

      IF (NRANK.LT.3) GO TO 60                                          TABX039 

      SLOPE(I)=MAT(2,4)+2.*MAT(3,4)*HPC(I)                              TABX040 

    4 DO 5 J=1,2                                                        TABX041 

      M=I+J-2                                                           TABX042 

      MAT(J,1)=1.                                                       TABX043 

      MAT(J,2)=HPC(M)                                                   TABX044 

      MAT(J,3)=HPC(M)**2                                                TABX045 

      MAT(J,4)=HPC(M)**3                                                TABX046 

      MAT(J,5)=FN2C(M)                                                  TABX047 

      L=J+2                                                             TABX048 

      MAT(L,1)=0.                                                       TABX049 

      MAT(L,2)=1.                                                       TABX050 

      MAT(L,3)=2.*HPC(M)                                                TABX051 

      MAT(L,4)=3.*HPC(M)**2                                             TABX052 

    5 MAT(L,5)=SLOPE(M)                                                 TABX053 

      CALL GAUSEL (MAT,4,4,5,NRANK)                                     TABX054 

      IF (NRANK.LT.4) GO TO 60                                          TABX055 

      ALPHA(I)=MAT(1,5)                                                 TABX056 

      BETA(I)=MAT(2,5)                                                  TABX057 

      GAMMA(I)=MAT(3,5)                                                 TABX058 

    6 DELTA(I)=MAT(4,5)                                                 TABX059 

      HMAX=HPC(NMAX)                                                    TABX060 

      NH=2                                                              TABX061 

   10 H=R(1)-EARTHR                                                     TABX062 

      F2=F*F                                                            TABX063 

      PXPR=PXPTH=PXPPH=0.                                               TABX064 

      IF (H.GE.HPC(1)) GO TO 12                                         TABX065 

   11 NH=2                                                              TABX066 

      X=0.                                                              TABX067 

      IF(FN2C(1).EQ.0.) GO TO 50                                        TABX068 

      X=FN2C(1)*EXP(A*(H-HPC(1)))/F2                                    TABX069 

      PXPR=A*X                                                          TABX070 

      GO TO 50                                                          TABX071 

   12 IF (H.GE.HPC(NOC)) GO TO 18                                       TABX072 

      NSTEP=1                                                           TABX073 

      IF (H.LT.HPC(NH-1)) NSTEP=-1                                      TABX074 

   15 IF (HPC(NH-1).LE.H.AND.H.LT.HPC(NH)) GO TO 16                     TABX075 

      NH=NH+NSTEP                                                       TABX076 

      GO TO 15                                                          TABX077 

   16 X=(ALPHA(NH)+H*(BETA(NH)+H*(GAMMA(NH)+H*DELTA(NH))))/F2           TABX078 

      PXPR=(BETA(NH)+H*(2.*GAMMA(NH)+H*3.*DELTA(NH)))/F2                TABX079 

      GO TO 50                                                          TABX080 

   18 X=FN2C(NOC)/F2                                                    TABX081 

   50 IF (PERT.NE.0.) CALL ELECT1                                       TABX082 

      RETURN                                                            TABX083 

   60 PRINT 6000, I,HPC(I)                                              TABX084 

 6000 FORMAT(4H THE,I4,55HTH POINT IN THE ELECTRON DENSITY PROFILE HAS TTABX085 

     1HE HEIGHT,F8.2,40H KM, WHICH IS THE SAME AS ANOTHER POINT.)       TABX086 

      CALL EXIT                                                         TABX087 

         END                                                            TABX088 

      SUBROUTINE GAUSEL (C,NRD,NRR,NCC,NSF)                             GAUS001 

C********* SAME AS SUBROUTINE GAUSSEL WRITTEN BY L. DAVID LEWIS ********GAUS002 

      DIMENSION C(NRD,NCC),L(128,2)                                     GAUS003 

C          BITS = 2.**-18                                               GAUS004 

      DATA (BITS=3.8146972656E-6)                                       GAUS005 

      NR=NRR                                                            GAUS006 

      NC=NCC                                                            GAUS007 

      IF(NC.LT.NR.OR.NR.GT.128.OR.NR.LE.0) CALL EXIT                    GAUS008 

C                                                                       GAUS009 

C          INITIALIZE.                                                  GAUS010 

      NSF=0                                                             GAUS011 

      NRM=NR-1                                                          GAUS012 

      NRP=NR+1                                                          GAUS013 

      D=1.                                                              GAUS014 

      LSD=1                                                             GAUS015 

      DO 1 KR=1,NR                                                      GAUS016 

      L(KR,1)=KR                                                        GAUS017 

1     L(KR,2)=0                                                         GAUS018 

      IF(NR.EQ.1) GO TO 42                                              GAUS019 

C                                                                       GAUS020 

C          ELIMINATION PHASE.                                           GAUS021 

      DO 41 KP=1,NRM                                                    GAUS022 

      KPP=KP+1                                                          GAUS023 

      PM=0.                                                             GAUS024 

      MPN=0                                                             GAUS025 

C                                                                       GAUS026 

C          SEARCH COLUMN KP FROM DIAGONAL DOWN FOR MAX PIVOT.           GAUS027 

      DO 2 KR=KP,NR                                                     GAUS028 

      LKR=L(KR,1)                                                       GAUS029 

      PT=ABS(C(LKR,KP))                                                 GAUS030 

      IF(PT.LE.PM) GO TO 2                                              GAUS031 

      PM=PT                                                             GAUS032 

      MPN=KR                                                            GAUS033 

      LMP=LKR                                                           GAUS034 

2     CONTINUE                                                          GAUS035 

C                                                                       GAUS036 

C          IF MAX PIVOT IS ZERO, MATRIX IS SINGULAR.                    GAUS037 

      IF(MPN.EQ.0) GO TO 9                                              GAUS038 

      NSF=NSF+1                                                         GAUS039 

      IF(MPN.EQ.KP) GO TO 3                                             GAUS040 

C                                                                       GAUS041 

C          NEW ROW NUMBER KP HAS MAX PIVOT.                             GAUS042 

      LSD=-LSD                                                          GAUS043 

      L(KP,2)=L(MPN,1)=L(KP,1)                                          GAUS044 

      L(KP,1)=LMP                                                       GAUS045 

C                                                                       GAUS046 

C          ROW OPERATIONS TO ZERO COLUMN KP BELOW DIAGONAL.             GAUS047 

3     MKP=L(KP,1)                                                       GAUS048 

      P=C(MKP,KP)                                                       GAUS049 

      D=D*P                                                             GAUS050 

      DO 41 KR=KPP,NR                                                   GAUS051 

      MKR=L(KR,1)                                                       GAUS052 

      Q=C(MKR,KP)/P                                                     GAUS053 

      IF(Q.EQ.0.) GO TO 41                                              GAUS054 

C                                                                       GAUS055 

C          SUBTRACT Q * PIVOT ROW FROM ROW KR.                          GAUS056 

      DO 4 LC=KPP,NC                                                    GAUS057 

      R=Q*C(MKP,LC)                                                     GAUS058 

      C(MKR,LC)=C(MKR,LC)-R                                             GAUS059 

4     IF(ABS(C(MKR,LC)).LT.ABS(R)*BITS) C(MKR,LC)=0.                    GAUS060 

41    CONTINUE                                                          GAUS061 

C                                                                       GAUS062 

C          LOWER RIGHT HAND CORNER.                                     GAUS063 

42    LNR=L(NR,1)                                                       GAUS064 

      P=C(LNR,NR)                                                       GAUS065 

      IF(P.EQ.0.) GO TO 9                                               GAUS066 

      NSF=NSF+1                                                         GAUS067 

      D=D*P*LSD                                                         GAUS068 

      IF(NR.EQ.NC) GO TO 8                                              GAUS069 

C                                                                       GAUS070 

C          BACK SOLUTION PHASE.                                         GAUS071 

      DO 61 MC=NRP,NC                                                   GAUS072 

      C(LNR,MC)=C(LNR,MC)/P                                             GAUS073 

      IF(NR.EQ.1) GO TO 61                                              GAUS074 

      DO 6 LL=1,NRM                                                     GAUS075 

      KR=NR-LL                                                          GAUS076 

      MR=L(KR,1)                                                        GAUS077 

      KRP=KR+1                                                          GAUS078 

      DO 5 MS=KRP,NR                                                    GAUS079 

      LMS=L(MS,1)                                                       GAUS080 

      R=C(MR,MS)*C(LMS,MC)                                              GAUS081 

      C(MR,MC)=C(MR,MC)-R                                               GAUS082 

5     IF(ABS(C(MR,MC)).LT.ABS(R)*BITS) C(MR,MC)=0.                      GAUS083 

6     C(MR,MC)=C(MR,MC)/C(MR,KR)                                        GAUS084 

61    CONTINUE                                                          GAUS085 

C                                                                       GAUS086 

C          SHUFFLE SOLUTION ROWS BACK TO NATURAL ORDER.                 GAUS087 

      DO 71 LL=1,NRM                                                    GAUS088 

      KR=NR-LL                                                          GAUS089 

      MKR=L(KR,2)                                                       GAUS090 

      IF(MKR.EQ.0) GO TO 71                                             GAUS091 

      MKP=L(KR,1)                                                       GAUS092 

      DO 7 LC=NRP,NC                                                    GAUS093 

      Q=C(MKR,LC)                                                       GAUS094 

      C(MKR,LC)=C(MKP,LC)                                               GAUS095 

7     C(MKP,LC)=Q                                                       GAUS096 

71    CONTINUE                                                          GAUS097 

C                                                                       GAUS098 

C          NORMAL AND SINGULAR RETURNS.  GOOD SOLUTION COULD HAVE D=0.  GAUS099 

8     C(1,1)=D                                                          GAUS100 

      GO TO 91                                                          GAUS101 

9     C(1,1)=0.                                                         GAUS102 

   91 RETURN                                                            GAUS103 

         END                                                            GAUS104-

      SUBROUTINE CHAPX                                                  CHAP001 

C          CHAPMAN LAYER WITH TILTS, RIPPLES, AND GRADIENTS             CHAP002 

C          W(104) = 0.5 FOR AN ALPHA-CHAPMAN LAYER                      CHAP003 

C                 = 1.0 FOR A BETA-CHAPMAN LAYER                        CHAP004 

      COMMON /CONST/ PI,PIT2,PID2,DUM(5)                                CHAP005 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  CHAP006 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 CHAP007 

      EQUIVALENCE (THETA,R(2))                                          CHAP008 

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(FC,W(101)),(HM,W(102)),       CHAP009 

     1 (SH,W(103)),(ALPHA,W(104)),(A,W(105)),(B,W(106)),(C,W(107)),     CHAP010 

     2 (E,W(108)),(PERT,W(150))                                         CHAP011 

      DATA (MODX(1)=6H CHAPX)                                           CHAP012 

      ENTRY ELECTX                                                      CHAP013 

      THETA2=THETA-PID2                                                 CHAP014 

      HMAX=HM+EARTHR*E*THETA2                                           CHAP015 

      H=R(1)-EARTHR                                                     CHAP016 

      Z=(H-HMAX)/SH                                                     CHAP017 

      D=0.                                                              CHAP018 

      IF (B.NE.0.) D=PIT2/B                                             CHAP019 

      TEMP=1.+A*SIN(D*THETA2)+C*THETA2                                  CHAP020 

      EXZ=1.-EXP(-Z)                                                    CHAP021 

      X=(FC/F)**2*TEMP*EXP(ALPHA*(EXZ-Z))                               CHAP022 

      PXPR=-ALPHA*X*EXZ/SH                                              CHAP023 

      PXPTH=X*(D*A*SIN(PID2-D*THETA2)+C)/TEMP-PXPR*EARTHR*E             CHAP024 

      PXPPH=0.                                                          CHAP0245

      IF (PERT.NE.0.) CALL ELECT1                                       CHAP025 

      RETURN                                                            CHAP026 

         END                                                            CHAP 27-

      SUBROUTINE VCHAPX                                                 VCHA001 

C          CHAPMAN LAYER WITH VARIABLE SCALE HEIGHT                     VCHA002 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  VCHA003 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 VCHA004 

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(FC,W(101)),(HM,W(102)),       VCHA005 

     1 (CHI,W(103)),(PERT,W(150))                                       VCHA006 

      DATA (MODX(1)=6HVCHAPX)                                           VCHA007 

      ENTRY ELECTX                                                      VCHA008 

      HMAX=HM                                                           VCHA009 

      X=PXPR=PXPTH=PXPPH=0.                                             VCHA010 

      H=R(1)-EARTHR                                                     VCHA011 

      IF (H.LE.0.) GO TO 50                                             VCHA012 

      TAU=(HM/H)**CHI                                                   VCHA013 

      X=(FC/F)**2*SQRT(TAU)*EXP(0.5*(1.-TAU))                           VCHA014 

      PXPR=.5*X*(TAU-1.)*CHI/H                                          VCHA015 

   50 IF (PERT.NE.0.) CALL ELECT1                                       VCHA016 

      RETURN                                                            VCHA017 

         END                                                            VCHA018-

      SUBROUTINE DCHAPT                                                 DCHA001 

C          TWO CHAPMAN LAYERS WITH TILTS                                DCHA002 

      COMMON /CONST/ PI,PIT2,PID2,DUM(5)                                DCHA003 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  DCHA004 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 DCHA005 

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(FC1,W(101)),(HM1,W(102)),     DCHA006 

     1 (SH1,W(103)),(FC2,W(104)),(HM2,W(105)),(SH2,W(106)),(C,W(107)),  DCHA007 

     2 (E,W(108)),(PERT,W(150))                                         DCHA008 

      DATA (MODX(1)=6HDCHAPT)                                           DCHA009 

      ENTRY ELECTX                                                      DCHA010 

      EARTHE=EARTHR*E                                                   DCHA011 

      THETA2=R(2)-PID2                                                  DCHA012 

      HMAX=HM1+EARTHE*THETA2                                            DCHA013 

      X=PXPR=PXPTH=PXPPH=0.                                             DCHA014 

      H=R(1)-EARTHR                                                     DCHA015 

      IF (H.LT.0.) GO TO 50                                             DCHA016 

      Z1=(H-HMAX)/SH1                                                   DCHA017 

      EXPZ1=1.-EXP(-Z1)                                                 DCHA018 

      TEMP=1.+C*THETA2                                                  DCHA019 

      X=(FC1/F)**2*TEMP*EXP(.5*(EXPZ1-Z1))                              DCHA020 

      PXPR=-0.5*X*EXPZ1/SH1                                             DCHA021 

      PXPTH=X*C/TEMP-PXPR*EARTHE                                        DCHA022 

      IF (FC2.EQ.0.) GO TO 50                                           DCHA023 

      Z2=(H-HM2-EARTHE*THETA2)/SH2                                      DCHA024 

      EXPZ2=1.-EXP(-Z2)                                                 DCHA025 

      X2=(FC2/F)**2*TEMP*EXP(.5*(EXPZ2-Z2))                             DCHA026 

      X=X+X2                                                            DCHA027 

      PXPR2=-0.5*X2*EXPZ2/SH2                                           DCHA028 

      PXPR=PXPR+PXPR2                                                   DCHA029 

      PXPTH=PXPTH+X2*C/TEMP-PXPR2*EARTHE                                DCHA030 

   50 IF (PERT.NE.0.) CALL ELECT1                                       DCHA031 

      RETURN                                                            DCHA032 

         END                                                            DCHA033-

      SUBROUTINE LINEAR                                                 LINE001 

C        LINEAR ELECTRON DENSITY MODEL                                  LINE002 

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RAD,K,DUM(2)                     LINE003 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  LINE004 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 LINE005 

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(FACT,W(101)),(HM,W(102)),     LINE006 

     1 (HMIN,W(103)),(PERT,W(150))                                      LINE007 

      REAL K                                                            LINE008 

      DATA (MODX(1)=6HLINEAR)                                           LINE009 

      ENTRY ELECTX                                                      LINE010 

      H=R(1)-EARTHR                                                     LINE011 

      HMAX=HM                                                           LINE012 

      X=PXPR=PXPTH=PXPPH=0.                                             LINE013 

      IF (H.LE.HMIN) GO TO 50                                           LINE014 

      PXPR=K*FACT/F**2                                                  LINE015 

      X=PXPR*(H-HMIN)                                                   LINE016 

   50 IF (PERT.NE.0.) CALL ELECT1                                       LINE017 

      RETURN                                                            LINE018 

         END                                                            LINE019-

      SUBROUTINE QPARAB                                                 PARA001 

C        PLAIN PARABOLIC OR QUASI-PARABOLIC PROFILE                     PARA002 

C        W(104) = 0. FOR A PLAIN PARABOLIC PROFILE                      PARA003 

C               = 1. FOR A QUASI-PARABOLIC PROFILE                      PARA004 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  PARA005 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 PARA006 

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(FC,W(101)),(HM,W(102)),       PARA007 

     1 (YM,W(103)),(QUASI,W(104)),(PERT,W(150))                         PARA008 

      DATA (MODX(1)=6HQPARAB)                                           PARA009 

      ENTRY ELECTX                                                      PARA010 

      HMAX=HM                                                           PARA011 

      X=PXPR=PXPTH=PXPPH=0.                                             PARA012 

      H=R(1)-EARTHR                                                     PARA013 

      FCF2=(FC/F)**2                                                    PARA014 

      CONST=1.                                                          PARA015 

      IF (QUASI.EQ.1.) CONST=(EARTHR+HM-YM)/R(1)                        PARA016 

      Z=(H-HM)/YM*CONST                                                 PARA017 

      X=MAX1F(0.,FCF2*(1.-Z*Z))                                         PARA018 

      IF (X.EQ.0.) GO TO 50                                             PARA019 

      IF (QUASI.EQ.1.) CONST=(EARTHR+HM)*(EARTHR+HM-YM)/R(1)**2         PARA020 

      PXPR=-2.*Z*FCF2/YM*CONST                                          PARA021 

   50 IF (PERT.NE.0.) CALL ELECT1                                       PARA022 

      RETURN                                                            PARA023 

         END                                                            PARA024-

      SUBROUTINE BULGE                                                  BULG001 

C        ANALYTICAL MODEL OF THE VARIATION OF THE EQUATORIAL F2 LAYER   BULG002 

C        IN GEOMAGNETIC LATITUDE  (EQUATORIAL BULGE AND ANOMALY)        BULG003 

C        SEE FIGURE 3.18B, PAGE 133 IN DAVIES (1965).                   BULG004 

C        THIS MODEL HAS NO VARIATION IN GEOMAGNETIC LONGITUDE.          BULG005 

      COMMON /CONST/ PI,PIT2,PID2,DUM(5)                                BULG006 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  BULG007 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 BULG008 

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(PERT,W(150))                  BULG009 

      DATA (MODX(1)=6H BULGE)                                           BULG010 

      ENTRY ELECTX                                                      BULG011 

      H=R(1)-EARTHR                                                     BULG012 

      PHMPTH=PFC2PTH=0.                                                 BULG013 

      HMAX=350.                                                         BULG014 

      FC2=225.                                                          BULG015 

      IF(H.LT.100.) GO TO 2                                             BULG016 

C        EQUATORIAL BULGE                                               BULG017 

      BULLAT=7.5*(PID2-R(2))                                            BULG018 

      IF(ABS(BUL LAT).GE.PI)  GO TO 1                                   BULG019 

      HMAX=430.+80.*COS(BULLAT)                                         BULG020 

      PHMPTH=600.*SIN(BUL LAT)                                          BULG021 

C        EQUATORIAL ANOMALY                                             BULG022 

    1 ANMLAT=22.5*(PID2-R(2))/PI                                        BULG023 

      POW=2.-ABS(ANM LAT)                                               BULG024 

      FC2=50.*ANM LAT**2*EXP(  POW     ) + 40.                          BULG025 

      PFC2PTH=-1125./PI*POW*ANMLAT*EXP(POW)                             BULG026 

C    FORCING PLASMA FREQ AT 100 KM TO BE 2 MHZ IN ORDER TO CALCULATE SH BULG027 

    2 ALPHA=2.*ALOG(FC2/4.)+1.                                          BULG028 

      Z100=-ALOG(ALPHA)                                                 BULG029 

      DO 3 I=1,5                                                        BULG030 

    3 Z100=-ALOG(ALPHA-Z100)                                            BULG031 

      SH=(100.-HMAX)/Z100                                               BULG032 

      Z=(H-HMAX)/SH                                                     BULG033 

      EXZ=1.-EXP(-Z)                                                    BULG034 

      X=FC2*EXP(.5*(EXZ-Z))/F**2                                        BULG035 

      PXPR=-0.5*X*EXZ/SH                                                BULG036 

      PXPTH=-PXPR*(1.-Z/Z100)*PHMPTH+(1.-Z*EXZ/(Z100*(1.-EXP(-Z100))))  BULG037 

     1      *X/FC2*PFC2PTH                                              BULG038 

      PXPPH=0.                                                          BULG0385

      IF (PERT.NE.0.) CALL ELECT1                                       BULG039 

      RETURN                                                            BULG040 

         END                                                            BULG041-

      SUBROUTINE EXPX                                                   EXPX001 

C        EXPONENTIAL ELECTRON DENSITY MODEL                             EXPX002 

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RAD,K,DUM(2)                     EXPX003 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  EXPX004 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 EXPX005 

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),                               EXPX006 

     1 (N0,W(101)),(H0,W(102)),(A,W(103)),(PERT,W(150))                 EXPX007 

      REAL N, N0 ,K                                                     EXPX008 

      DATA (MODX(1)=4HEXPX),(HMAX=350.)                                 EXPX009 

      ENTRY ELECTX                                                      EXPX010 

      H=R(1)-EARTHR                                                     EXPX011 

      N =N0 * EXP(A*(H-H0))                                             EXPX012 

      X=K*N/F**2                                                        EXPX013 

      PXPR=A*X                                                          EXPX014 

      PXPTH=PXPPH=0.                                                    EXPX0145

      IF (PERT.NE.0.) CALL ELECT1                                       EXPX015 

      RETURN                                                            EXPX016 

         END                                                            EXPX017-

      SUBROUTINE ELECT1                                                 ELEC001 

C          USE WHEN AN ELECTRON DENSITY PERTURBATION IS NOT WANTED      ELEC002 

      COMMON /XX/ MODX(2),X(6)                                          ELEC003 

      COMMON /WW/ ID(10),W0,W(400)                                      ELEC004 

      EQUIVALENCE (PERT,W(150))                                         ELEC005 

      DATA (MODX(2)=6H NONE )                                           ELEC006 

      PERT=0.                                                           ELEC007 

      RETURN                                                            ELEC008 

         END                                                            ELEC009-

      SUBROUTINE TORUS                                                  TOR 001 

      COMMON /CONST/ PI,PIT2,PID2,DUM(5)                                TOR 002 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  TOR 003 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 TOR 004 

      EQUIVALENCE (EARTHR,W(2)),(C0,W(151)),(A,W(152)),(B,W(153)),      TOR 005 

     1 (BETA,W(154)),(H0,W(155))                                        TOR 006 

      REAL LAMBDA                                                       TOR 007 

      DATA (PDPP=0.),(MODX(2)=6H TORUS)                                 TOR 008 

      ENTRY ELECT1                                                      TOR 009 

      IF (X.EQ.0..AND.PXPR.EQ.0..AND.PXPTH.EQ.0..AND.PXPPH.EQ.0.) RETURNTOR 010 

      IF (C0.EQ.0.) RETURN                                              TOR 011 

      R0=EARTHR+H0                                                      TOR 012 

      Z=R(1)-R0                                                         TOR 013 

      LAMBDA=R0*(R(2)-PID2)                                             TOR 014 

      SINBET=SIN(BETA)                                                  TOR 015 

      COSBET=COS(BETA)                                                  TOR 016 

      P=LAMBDA*COSBET+Z*SINBET                                          TOR 017 

      Y=Z*COSBET-LAMBDA*SINBET                                          TOR 018 

      DELTA=C0*EXP(-(P/A)**2-(Y/B)**2)                                  TOR 019 

      DEL1=DELTA+1.                                                     TOR 020 

      PDPR=-2.*DELTA*(P*SINBET/A**2+Y*COSBET/B**2)                      TOR 021 

      PDPT=-2.*DELTA*(P*R0*COSBET/A**2-Y*R0*SINBET/B**2)                TOR 022 

      PXPR=PXPR*DEL1+X*PDPR                                             TOR 023 

      PXPTH=PXPTH*DEL1+X*PDPT                                           TOR 024 

      PXPPH=PXPPH*DEL1+X*PDPP                                           TOR 025 

      X=X*DEL1                                                          TOR 026 

      RETURN                                                            TOR 027 

         END                                                            TOR 028-

      SUBROUTINE DTORUS                                                 DTOR001 

      COMMON /CONST/ PI,PIT2,PID2,DUM(5)                                DTOR002 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT                       DTOR003 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 DTOR004 

      EQUIVALENCE (EARTHR,W(2)),(C1,W(151)),(A1,W(152)),(B1,W(153)),    DTOR005 

     1 (BETA,W(154)),(H1,W(155)),(C2,W(156)),(A2,W(157)),(B2,W(158)),   DTOR006 

     2 (H2,W(159))                                                      DTOR007 

      REAL LAMBDA1,LAMBDA2                                              DTOR008 

      DATA (MODX(2)=6HDTORUS),(PDPP=0.)                                 DTOR009 

      ENTRY ELECT1                                                      DTOR010 

      IF (X.EQ.0..AND.PXPR.EQ.0..AND.PXPTH.EQ.0..AND.PXPPH.EQ.0.) RETURNDTOR011 

      IF (C1.EQ.0.) RETURN                                              DTOR012 

      R1=EARTHR+H1                                                      DTOR013 

      R2=EARTHR+H2                                                      DTOR014 

      Z1=R(1)-R1                                                        DTOR015 

      Z2=R(1)-R2                                                        DTOR016 

      LAMBDA1=R1*(R(2)-PID2)                                            DTOR017 

      LAMBDA2=R2*(R(2)-PID2+(H1-H2)/R2/TANF(BETA))                      DTOR018 

      SINBET=SIN(BETA)                                                  DTOR019 

      COSBET=COS(BETA)                                                  DTOR020 

      P1=LAMBDA1*COSBET+Z1*SINBET                                       DTOR021 

      P2=LAMBDA2*COSBET+Z2*SINBET                                       DTOR022 

      Y1=Z1*COSBET-LAMBDA1*SINBET                                       DTOR023 

      Y2=Z2*COSBET-LAMBDA2*SINBET                                       DTOR024 

      DELTA1=C1*EXP(-(P1/A1)**2-(Y1/B1)**2)                             DTOR025 

      DELTA2=C2*EXP(-(P2/A2)**2-(Y2/B2)**2)                             DTOR026 

      DEL1=1.+DELTA1+DELTA2                                             DTOR027 

      PDPR1=-2.*DELTA1*(P1*SINBET/A1**2+Y1*COSBET/B1**2)                DTOR028 

      PDPR2=-2.*DELTA2*(P2*SINBET/A2**2+Y2*COSBET/B2**2)                DTOR029 

      PDPT1=-2.*DELTA1*(P1*R1*COSBET/A1**2-Y1*R1*SINBET/B1**2)          DTOR030 

      PDPT2=-2.*DELTA2*(P2*R2*COSBET/A2**2-Y2*R2*SINBET/B2**2)          DTOR031 

      PXPR=PXPR*DEL1+X*(PDPR1+PDPR2)                                    DTOR032 

      PXPTH=PXPTH*DEL1+X*(PDPT1+PDPT2)                                  DTOR033 

      PXPPH=PXPPH*DEL1*PDPP                                             DTOR034 

      X=X*DEL1                                                          DTOR035 

      RETURN                                                            DTOR036 

         END                                                            DTOR037-

      SUBROUTINE TROUGH                                                 TROU001 

C             A PERTURBATION TO AN ELECTRON DENSITY MODEL               TROU002 

      COMMON /CONST/ PI,PIT2,PID2,DUM(5)                                TROU003 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  TROU004 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 TROU005 

      EQUIVALENCE (A,W(151)),(B,W(152)),(ALAT,W(153)),(FACTOR,W(154))   TROU006 

      DATA (MODX(2)=6HTROUGH)                                           TROU007 

      ENTRY ELECT1                                                      TROU008 

      IF (X.EQ.0..AND.PXPR.EQ.0..AND.PXPTH.EQ.0..AND.PXPPH.EQ.0.) RETURNTROU009 

      IF (A.EQ.0.) RETURN                                               TROU010 

      ANGLE=R(2)+ALAT-PID2                                              TROU011 

      WIDTH=B                                                           TROU012 

      IF (ANGLE.GT.0.)   WIDTH=FACTOR*B                                 TROU013 

      ANGLE=ANGLE/WIDTH                                                 TROU014 

      DELTA=A*EXP(-ANGLE**2)                                            TROU015 

      DEL1=DELTA+1.                                                     TROU016 

      PXPR=PXPR*DEL1                                                    TROU017 

      PXPTH=PXPTH*DEL1-2.*X*ANGLE*DELTA/WIDTH                           TROU018 

      PXPPH=PXPPH*DEL1                                                  TROU019 

      X=X*DEL1                                                          TROU020 

      RETURN                                                            TROU021 

         END                                                            TROU022-

      SUBROUTINE SHOCK                                                  SHOC001 

C   A PERTURBATION TO AN ELECTRON DENSITY MODEL SIMULATING A SHOCK WAVE SHOC002 

      COMMON /CONST/ PI,PIT2,PID2,DUM(5)                                SHOC003 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  SHOC004 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 SHOC005 

      EQUIVALENCE (EARTHR,W(2)),(P,W(151)),(WW,W(152)),(ALAT,W(153)),   SHOC006 

     1 (ALON,W(154)),(S,W(155)),(H0,W(156))                             SHOC007 

      REAL LAT,LON                                                      SHOC008 

      DATA (MODX(2)=6H SHOCK)                                           SHOC009 

      ENTRY ELECT1                                                      SHOC010 

      IF (X.EQ.0..AND.PXPR.EQ.0..AND.PXPTH.EQ.0..AND.PXPPH.EQ.0.) RETURNSHOC011 

      IF (P.EQ.0..OR.WW.EQ.0.)  RETURN                                  SHOC012 

      H=R(1)-EARTHR                                                     SHOC013 

      RHOC=S*(H-H0)-WW                                                  SHOC014 

      LON=R(3)-ALON                                                     SHOC015 

      LAT=PID2-R(2)-ALAT                                                SHOC016 

      COSLON=COS(LON)                                                   SHOC017 

      COSLAT=COS(LAT)                                                   SHOC018 

      U=COSLON*COSLAT                                                   SHOC019 

      RHO=R(1)*ACOS(U)                                                  SHOC020 

      DIF=RHOC-RHO                                                      SHOC021 

      CON=-9./WW**2                                                     SHOC022 

      CONS=P*EXP(CON*DIF**2)                                            SHOC023 

      CONST=1.+CONS                                                     SHOC024 

      CON=2.*CON*CONS*DIF                                               SHOC025 

      PXPR=PXPR*CONST+X*CON*(S-RHO/R(1))                                SHOC026 

      CONS=R(1)*(1./SQRT(1.-U**2))                                      SHOC027 

      PXPTH=PXPTH*CONST+X*CON*CONS*COSLON*SIN(LAT)                      SHOC028 

      PXPPH=PXPPH*CONST-X*CON*CONS*COSLAT*SIN(LON)                      SHOC029 

      X=X*CONST                                                         SHOC030 

      RETURN                                                            SHOC031 

         END                                                            SHOC032-

      SUBROUTINE WAVE                                                   WAVE001 

C          PERTURBATION TO AN ALPHA-CHAPMAN ELECTRON DENSITY MODEL      WAVE002 

      COMMON /CONST/ PI,PIT2,PID2,DUM(5)                                WAVE003 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT,HMAX                  WAVE004 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 WAVE005 

      EQUIVALENCE (EARTHR,W(2)),(Z0,W(151)),(SH,W(152)),(DELTA,W(153)), WAVE006 

     1 (VSUBX,W(154)),(LAMBDAX,W(155)),(LAMBDAZ,W(156)),(TP,W(157))     WAVE007 

      REAL LAMBDAX,LAMBDAZ                                              WAVE008 

      DATA (MODX(2)=6H WAVE )                                           WAVE009 

      ENTRY ELECT1                                                      WAVE010 

      IF (X.EQ.0..AND.PXPR.EQ.0..AND.PXPTH.EQ.0..AND.PXPPH.EQ.0.) RETURNWAVE011 

      IF (DELTA.EQ.0..OR.SH.EQ.0.)   RETURN                             WAVE012 

      H=R(1)-EARTHR                                                     WAVE013 

      EXPQ=EXP(-((H-Z0)/SH)**2)                                         WAVE014 

      TMP=PIT2*(TP+(PID2-R(2))*EARTHR/LAMBDAX+H/LAMBDAZ)                WAVE015 

      SINW=SIN(TMP)                                                     WAVE016 

      COSW=SIN(PID2-TMP)                                                WAVE017 

      CONS=1.0+DELTA*EXPQ*COSW                                          WAVE018 

      IF (H.NE.0.) PXPR=PXPR*CONS-X*DELTA*EXPQ*(2.0/SH**2*(H-Z0)*COSW   WAVE019 

     1 +PIT2/LAMBDAZ*SINW)                                              WAVE020 

      PXPTH=PXPTH*CONS+X*DELTA*PIT2*EARTHR/LAMBDAX*SINW*EXPQ            WAVE021 

      PXPPH=PXPPH*CONS                                                  WAVE022 

      PXPT=0.                                                           WAVE023 

      IF (VSUBX.NE.0.) PXPT=-PIT2*VSUBX/LAMBDAX*X*DELTA*EXPQ*SINW       WAVE024 

      X=X*CONS                                                          WAVE025 

      RETURN                                                            WAVE026 

         END                                                            WAVE027-

      SUBROUTINE WAVE2                                                  WAV2001 

C          PERTURBATION TO AN ANY ELECTRON DENSITY MODEL                WAV2002 

      COMMON /CONST/ PI,PIT2,PID2,DUM(5)                                WAV2003 

      COMMON /XX/ MODX(2),X,PXPR,PXPTH,PXPPH,PXPT                       WAV2004 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 WAV2005 

      EQUIVALENCE (EARTHR,W(2)),(Z0,W(151)),(SH,W(152)),(DELTA,W(153)), WAV2006 

     1 (VSUBX,W(154)),(LAMBDAX,W(155)),(LAMBDAZ,W(156)),(TP,W(157)),    WAV2007 

     2 (TH00,W(158)),(THC,W(159)),(PHIC,W(160)),(VGX,W(161))            WAV2008 

      REAL LAMBDAX,LAMBDAZ                                              WAV2009 

      DATA (MODX(2)=6H WAVE2)                                           WAV2010 

      ENTRY ELECT1                                                      WAV2011 

      IF (X.EQ.0..AND.PXPR.EQ.0..AND.PXPTH.EQ.0..AND.PXPPH.EQ.0.) RETURNWAV2012 

      IF (DELTA.EQ.0..OR.SH.EQ.0.)   RETURN                             WAV2013 

      H=R(1)-EARTHR                                                     WAV2014 

      TH0=TH00+LAMBDAX*TP*VGX/VSUBX/EARTHR                              WAV2015 

      EXPR=EXP(-((H-Z0)/SH)**2)                                         WAV2016 

      EXPTH=EXP(-((R(2)-TH0)/THC)**2)                                   WAV2017 

      EXPPHI=EXP(-(R(3)/PHIC)**2)                                       WAV2018 

      WW=PIT2*(TP+(PID2-R(2))*EARTHR/LAMBDAX+H/LAMBDAZ)                 WAV2019 

      SINW=SIN(WW)                                                      WAV2020 

      COSW=COS(WW)                                                      WAV2021 

      E=DELTA*EXPR*EXPTH*EXPPHI                                         WAV2022 

      CONS=1.0+E*COSW                                                   WAV2023 

      PXPR=PXPR*CONS-X*E*2.*(COSW*(H-Z0)/SH**2+PI/LAMBDAZ*SINW)         WAV2024 

      PXPTH=PXPTH*CONS+2.*E*(X*PI*EARTHR*SINW/LAMBDAX-(R(2)-TH0)/       WAV2025 

     1 THC**2*COSW)                                                     WAV2026 

      PXPPH=PXPPH*CONS-X*2.*E*R(3)/PHIC**2*COSW                         WAV2027 

      PXPT=-PIT2*VSUBX*E/LAMBDAX*SINW+2.0*E*VGX/EARTHR*COSW*(R(2)-TH0   WAV2028 

     1 -LAMBDAX*TP/EARTHR)/THC                                          WAV2029 

      X=X*CONS                                                          WAV2030 

      RETURN                                                            WAV2031 

         END                                                            WAV2032-

      SUBROUTINE DOPPLER                                                DOPP001 

C      COMPUTES DN/DT                 FROM PROFILES HAVING THE SAME FORMDOPP002 

C      AS THOSE USED BY SUBROUTINE TABLE X                              DOPP003 

C MAKES AN EXPONENTIAL EXTRAPOLATION DOWN USING THE BOTTOM TWO POINTS   DOPP004 

C        NEEDS SUBROUTINE GAUSEL                                        DOPP005 

      DIMENSION HPC(250),FN2C(250),ALPHA(250),BETA(250),GAMMA(250),     DOPP006 

     1 DELTA(250),SLOPE(250),MAT(4,5)                                   DOPP007 

      COMMON /CONST/ PI,PIT2,PID2,DEGS,RAD,K,DUM(2)                     DOPP008 

      COMMON /XX/ MODX(2),XDUM,PXPR,PXPTH,PXPPH,X,HMAX                  DOPP009 

      COMMON R(6) /WW/ ID(10),W0,W(400)                                 DOPP010 

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(READFN,W(151))                DOPP011 

      REAL MAT,K                                                        DOPP012 

      DATA (MODX(2)=7HDOPPLER)                                          DOPP013 

      ENTRY ELECT1                                                      DOPP014 

      IF (READFN.EQ.0.) GO TO 10                                        DOPP015 

      READFN=0.                                                         DOPP016 

      READ 1000, NOC,(HPC(I),FN2C(I),I=1,NOC)                           DOPP017 

 1000 FORMAT (I4/(F8.2,E12.4))                                          DOPP018 

      PRINT 1200, (HPC(I),FN2C(I),I=1,NOC)                              DOPP019 

 1200 FORMAT(1H1,14X,6HHEIGHT,4X,16H     DN/DT      /(1X,F20.10,E20.10))DOPP020 

      A=0.                                                              DOPP021 

      IF(FN2C(1).NE.0.) A=ALOG(FN2C(2)/FN2C(1))/(HPC(2)-HPC(1))         DOPP022 

      FN2C(1)=K*FN2C(1)                                                 DOPP023 

      FN2C(2)=K*FN2C(2)                                                 DOPP024 

      SLOPE(1)=A*FN2C(1)                                                DOPP025 

      SLOPE(NOC)=0.                                                     DOPP026 

      NMAX=1                                                            DOPP027 

      DO 6 I=2,NOC                                                      DOPP028 

      IF (FN2C(I).GT.FN2C(NMAX)) NMAX=I                                 DOPP029 

      IF (I.EQ.NOC) GO TO 4                                             DOPP030 

      FN2C(I+1)=K*FN2C(I+1)                                             DOPP031 

      DO 3 J=1,3                                                        DOPP032 

      M=I+J-2                                                           DOPP033 

      MAT(J,1)=1.                                                       DOPP034 

      MAT(J,2)=HPC(M)                                                   DOPP035 

      MAT(J,3)=HPC(M)**2                                                DOPP036 

  3   MAT(J,4)=FN2C(M)                                                  DOPP037 

      CALL GAUSEL (MAT,4,3,4,NRANK)                                     DOPP038 

      IF (NRANK.LT.3) GO TO 60                                          DOPP039 

      SLOPE(I)=MAT(2,4)+2.*MAT(3,4)*HPC(I)                              DOPP040 

    4 DO 5 J=1,2                                                        DOPP041 

      M=I+J-2                                                           DOPP042 

      MAT(J,1)=1.                                                       DOPP043 

      MAT(J,2)=HPC(M)                                                   DOPP044 

      MAT(J,3)=HPC(M)**2                                                DOPP045 

      MAT(J,4)=HPC(M)**3                                                DOPP046 

      MAT(J,5)=FN2C(M)                                                  DOPP047 

      L=J+2                                                             DOPP048 

      MAT(L,1)=0.                                                       DOPP049 

      MAT(L,2)=1.                                                       DOPP050 

      MAT(L,3)=2.*HPC(M)                                                DOPP051 

      MAT(L,4)=3.*HPC(M)**2                                             DOPP052 

    5 MAT(L,5)=SLOPE(M)                                                 DOPP053 

      CALL GAUSEL (MAT,4,4,5,NRANK)                                     DOPP054 

      IF (NRANK.LT.4) GO TO 60                                          DOPP055 

      ALPHA(I)=MAT(1,5)                                                 DOPP056 

      BETA(I)=MAT(2,5)                                                  DOPP057 

      GAMMA(I)=MAT(3,5)                                                 DOPP058 

    6 DELTA(I)=MAT(4,5)                                                 DOPP059 

      HMAX=AMAX1(HMAX,HPC(NMAX))                                        DOPP060 

      NH=2                                                              DOPP061 

   10 H=R(1)-EARTHR                                                     DOPP062 

      F2=F*F                                                            DOPP063 

      IF (H.GE.HPC(1)) GO TO 12                                         DOPP064 

   11 NH=2                                                              DOPP065 

      X=0.                                                              DOPP066 

      IF(FN2C(1).EQ.0.) GO TO 50                                        DOPP067 

      X=FN2C(1)*EXP(A*(H-HPC(1)))/F2                                    DOPP068 

      GO TO 50                                                          DOPP069 

   12 IF (H.GE.HPC(NOC)) GO TO 18                                       DOPP070 

      NSTEP=1                                                           DOPP071 

      IF (H.LT.HPC(NH-1)) NSTEP=-1                                      DOPP072 

   15 IF (HPC(NH-1).LE.H.AND.H.LT.HPC(NH)) GO TO 16                     DOPP073 

      NH=NH+NSTEP                                                       DOPP074 

      GO TO 15                                                          DOPP075 

   16 X=(ALPHA(NH)+H*(BETA(NH)+H*(GAMMA(NH)+H*DELTA(NH))))/F2           DOPP076 

      GO TO 50                                                          DOPP077 

   18 X=FN2C(NOC)/F2                                                    DOPP078 

   50 CONTINUE                                                          DOPP079 

      RETURN                                                            DOPP080 

   60 PRINT 6000, I,HPC(I)                                              DOPP081 

 6000 FORMAT(4H THE,I4,55HTH POINT IN THE   DN/DT          PROFILE HAS TDOPP082 

     1HE HEIGHT,F8.2,40H KM, WHICH IS THE SAME AS ANOTHER POINT.)       DOPP083 

      CALL EXIT                                                         DOPP084 

         END                                                            DOPP085 

