      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

