      SUBROUTINE CONSTY                                                 CONY001 

C         CONSTANT DIP AND GYROFREQUENCY                                CONY002 

      COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRCONY003 

     1,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP                                CONY004 

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

      EQUIVALENCE (F,W(6)),(FH,W(201)),(DIP,W(202))                     CONY006 

      DATA (MODY=6HCONSTY)                                              CONY007 

      ENTRY MAGY                                                        CONY008 

      Y=FH/F                                                            CONY009 

      YR=Y*SIN(DIP)                                                     CONY010 

      YTH=Y*COS(DIP)                                                    CONY011 

      RETURN                                                            CONY012 

         END                                                            CONY013-

      SUBROUTINE DIPOLY                                                 DIPO001 

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

      COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRDIPO003 

     1,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP                                DIPO004 

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

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(FH,W(201))                    DIPO006 

      DATA (MODY=6HDIPOLY)                                              DIPO007 

      ENTRY MAGY                                                        DIPO008 

      SINTH=SIN(R(2))                                                   DIPO009 

      COSTH=SIN(PID2-R(2))                                              DIPO010 

      TERM9=SQRT(1.+3.*COSTH**2)                                        DIPO011 

      T1=FH*(EARTHR/R(1))**3/F                                          DIPO012 

      Y=T1*TERM9                                                        DIPO013 

      YR= 2.*T1*COSTH                                                   DIPO014 

      YTH= T1*SINTH                                                     DIPO015 

      PYRPR=-3.*YR/R(1)                                                 DIPO016 

      PYRPT=-2.*YTH                                                     DIPO017 

      PYTPR=-3.*YTH/R(1)                                                DIPO018 

      PYTPT=.5*YR                                                       DIPO019 

      PYPR=-3.*Y/R(1)                                                   DIPO020 

      PYPTH=-3.*Y*SINTH*COSTH/TERM9**2                                  DIPO021 

      RETURN                                                            DIPO022 

         END                                                            DIPO023-

       SUBROUTINE CUBEY                                                 CUBE001 

C        CONSTANT DIP.                                                  CUBE002 

C        GYROFREQ DECREASES AS CUBE OF DISTANCE FROM CENTER OF EARTH.   CUBE003 

C        THIS MODEL HAS SAME HEIGHT VARIATION AS A DIPOLE FIELD.        CUBE004 

      COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRCUBE005 

     1,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP                                CUBE006 

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

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(FH,W(201)),(DIP,W(202))       CUBE008 

      DATA(MODY=5HCUBEY)                                                CUBE009 

      ENTRY MAGY                                                        CUBE010 

       Y=(EARTHR/R)**3 *FH/F                                            CUBE011 

      YR= Y*SIN(DIP)                                                    CUBE012 

      YTH= Y*COS(DIP)                                                   CUBE013 

      PYPR=-3.*Y/R                                                      CUBE014 

      PYRPR=-3.*YR/R                                                    CUBE015 

      PYTPR=-3.*YTH/R                                                   CUBE016 

      RETURN                                                            CUBE017 

         END                                                            CUBE018-

      SUBROUTINE HARMONY                                                HARM001 

C   MODEL OF THE EARTH S MAGNETIC FIELD BASED ON A HARMONIC ANALYSIS    HARM002 

      DIMENSION PHPTH(7,7),PGPTH(7,7),A1(7,7),B1(7,7)                   HARM003 

      DIMENSION H(7,7),G(7,7),GG(7,7),HH(7,7),SINP(7),COSP(7)           HARM004 

      COMMON /YY/ MODY,Y,PYPR,PYPTH,PYPPH,YR,PYRPR,PYRPT,PYRPP,YTH,PYTPRHARM005 

     1,PYTPT,PYTPP,YPH,PYPPR,PYPPT,PYPPP                                HARM006 

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

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

      EQUIVALENCE (THETA,R(2)),(PHI,R(3))                               HARM009 

      EQUIVALENCE (EARTHR,W(2)),(F,W(6)),(READFH,W(200))                HARM010 

C   RATIO OF CHARGE TO MASS FOR ELECTRON                                HARM011 

      DATA(EOM=1.7589E7)                                                HARM012 

      DATA (SET=0.),(H=1.,48(0.)),(G=49(0.)),(PHPTH=49(0.))             HARM013 

     1 ,(PGPTH=49(0.)),(MODY=7HHARMONY)                                 HARM014 

      ENTRY MAGY                                                        HARM015 

      IF(SET) GO TO 2                                                   HARM016 

      DO 1 M=1,7                                                        HARM017 

      DO 1 N=1,7                                                        HARM018 

      B1(M,N)=(N+M-1)*(N-M+1)/(2*N-1.)                                  HARM019 

    1 A1(M,N)=B1(M,N)/(2*N+1)                                           HARM020 

      SET=1.                                                            HARM021 

    2 IF(READFH.EQ.0.) GO TO 3                                          HARM022 

      READ 2000,GG,HH                                                   HARM023 

 2000 FORMAT(1X,F9.4,6F10.4)                                            HARM024 

      PRINT 2100,GG                                                     HARM025 

 2100 FORMAT(1H1,10X,1H0,14X,1H1,14X,1H2,14X,1H3,14X,1H4,14X,1H5,14X,1H6HARM026 

     1 /9X,7(1HG,14X)/10X,7(1HN,14X)//(1X,7F15.6))                      HARM027 

      PRINT 2200,HH                                                     HARM028 

 2200 FORMAT(//  11X,1H0,14X,1H1,14X,1H2,14X,1H3,14X,1H4,14X,1H5,14X,1H6HARM029 

     1 /9X,7(1HH,14X)/10X,7(1HN,14X)//(1X,7F15.6))                      HARM030 

      READFH=0.                                                         HARM031 

  3   COSTHE=COS(THETA)                                                 HARM032 

      SINTHE=SIN(THETA)                                                 HARM033 

      AOR=EARTHR/R(1)                                                   HARM034 

      PAORPR=-AOR/R(1)                                                  HARM035 

      CNST2=AOR                                                         HARM036 

      PCNSPR=PAORPR                                                     HARM037 

      FIN1=PFIN1R=PFIN1T=PFIN1P=0.                                      HARM038 

      FIN2=PFIN2R=PFIN2T=PFIN2P=0.                                      HARM039 

      FIN3=PFIN3R=PFIN3T=PFIN3P=0.                                      HARM040 

      DO 4 M=1,7                                                        HARM041 

      SINP(M)=SIN((M-1)*PHI)                                            HARM042 

  4   COSP(M)=COS((M-1)*PHI)                                            HARM043 

      H(1,2)=COSTHE                                                     HARM044 

      H(2,2)=SINTHE                                                     HARM045 

      DO 5 M=1,5                                                        HARM046 

      H(M+1,M+2)=COSTHE*H(M+1,M+1)                                      HARM047 

      H(M+2,M+2)=SINTHE*H(M+1,M+1)                                      HARM048 

      DO 5 N=M,5                                                        HARM049 

  5   H(M,N+2)=COSTHE*H(M,N+1)-A1(M,N)*H(M,N)                           HARM050 

      DO 6 M=1,6                                                        HARM051 

      G(M+1,M+1)=-M*COSTHE*H(M+1,M+1)                                   HARM052 

      PHPTH(M+1,M+1)=-G(M+1,M+1)/SINTHE                                 HARM053 

      PGPTH(M+1,M+1)=M*SINTHE*H(M+1,M+1)-M*COSTHE*PHPTH(M+1,M+1)        HARM054 

      DO 6 N=M,6                                                        HARM055 

      G(M,N+1)=-N*COSTHE*H(M,N+1)+B1(M,N)*H(M,N)                        HARM056 

      PHPTH(M,N+1)=-G(M,N+1)/SINTHE                                     HARM057 

  6   PGPTH(M,N+1)=N*SINTHE*H(M,N+1)-N*COSTHE*PHPTH(M,N+1)+B1(M,N)*PHPTHHARM058 

     1 (M,N)                                                            HARM059 

      DO 8 N=1,7                                                        HARM060 

      CR=PCRPTH=PCRPPH=0.                                               HARM061 

      CTH=PCTHPT=PCTHPP=0.                                              HARM062 

      CPH=PCPHPT=PCPHPP=0.                                              HARM063 

      DO 7 M=1,N                                                        HARM064 

      TEMP1=GG(M,N)*COSP(M)+HH(M,N)*SINP(M)                             HARM065 

      TEMP2=(M-1)*(HH(M,N)*COSP(M)-GG(M,N)*SINP(M))                     HARM066 

      CR   =CR   +H(M,N)*TEMP1                                          HARM067 

      PCRPTH=PCRPTH+PHPTH(M,N)*TEMP1                                    HARM068 

      PCRPPH=PCRPPH+H(M,N)*TEMP2                                        HARM069 

      CTH   =CTH   +G(M,N)*TEMP1                                        HARM070 

      PCTHPT=PCTHPT+PGPTH(M,N)*TEMP1                                    HARM071 

      PCTHPP=PCTHPP+G(M,N)*TEMP2                                        HARM072 

      CPH   =CPH   +H(M,N)*TEMP2                                        HARM073 

      PCPHPT=PCPHPT+PHPTH(M,N)*TEMP2                                    HARM074 

  7   PCPHPP=PCPHPP-H(M,N)*(M-1)**2*TEMP1                               HARM075 

      CNST2=CNST2*AOR                                                   HARM076 

      PCNSPR=CNST2*PAORPR+AOR*PCNSPR                                    HARM077 

      FIN1=FIN1+N*CNST2*CR                                              HARM078 

      PFIN1R=PFIN1R+N*PCNSPR*CR                                         HARM079 

      PFIN1T=PFIN1T+N*CNST2*PCRPTH                                      HARM080 

      PFIN1P=PFIN1P+N*CNST2*PCRPPH                                      HARM081 

      FIN2=FIN2+CNST2*CTH                                               HARM082 

      PFIN2R=PFIN2R+PCNSPR*CTH                                          HARM083 

      PFIN2T=PFIN2T+CNST2*PCTHPT                                        HARM084 

      PFIN2P=PFIN2P+CNST2*PCTHPP                                        HARM085 

      FIN3=FIN3+CNST2*CPH                                               HARM086 

      PFIN3R=PFIN3R+PCNSPR*CPH                                          HARM087 

      PFIN3T=PFIN3T+CNST2*PCPHPT                                        HARM088 

  8   PFIN3P=PFIN3P+CNST2*PCPHPP                                        HARM089 

      HTHETA=-FIN2/SINTHE                                               HARM090 

      HPHI=FIN3/SINTHE                                                  HARM091 

C********* CONVERT FROM MAG FIELD IN GAUSS TO GYROFREQ IN MHZ           HARM092 

      CONST=-EOM/PIT2*1.E-6/F                                           HARM093 

      YR=-CONST*FIN1                                                    HARM094 

      YTH=CONST*HTHETA                                                  HARM095 

      YPH=CONST*HPHI                                                    HARM096 

      Y=SQRT(YR**2+YTH**2+YPH**2)                                       HARM097 

      PYRPR=-CONST*PFIN1R                                               HARM098 

      PYTPR=-CONST*PFIN2R/SINTHE                                        HARM099 

      PYPPR=CONST*PFIN3R/SINTHE                                         HARM100 

      PYPR=(YR*PYRPR+YTH*PYTPR+YPH*PYPPR)/Y                             HARM101 

      PYRPT=-CONST*PFIN1T                                               HARM102 

      PYTPT=-CONST*(PFIN2T/SINTHE+HTHETA*COSTHE/SINTHE)                 HARM103 

      PYPPT=CONST*(PFIN3T/SINTHE-HPHI*COSTHE/SINTHE)                    HARM104 

      PYPTH=(YR*PYRPT+YTH*PYTPT+YPH*PYPPT)/Y                            HARM105 

      PYRPP=-CONST*PFIN1P                                               HARM106 

      PYTPP=-CONST*PFIN2P/SINTHE                                        HARM107 

      PYPPP=CONST*PFIN3P/SINTHE                                         HARM108 

      PYPPH=(YR*PYRPP+YTH*PYTPP+YPH*PYPPP)/Y                            HARM109 

      RETURN                                                            HARM110 

C   COEFFICIENTS IN GAUSSIAN UNITS FROM JONES AND MELOTTE (1953).       HARM111 

C   THE FOLLOWING 14 CARDS CAN BE USED AS DATA CARDS FOR THIS SUBROUTINEHARM112 

C   0.                                                                  HARM113 

C.3039    .0218                                                         HARM114 

C.0176    -.0509    -.0135                                              HARM115 

C-.0255    .0515    -.0236    -.0074                                    HARM116 

C-.0398   -.0397    -.0238     .0087    -.0018                          HARM117 

C.0293    -.0329    -.0130     .0031     .0034     .0005                HARM118 

C-.0211   -.0073    -.0007     .0210     .0017    -.0004     .0006      HARM119 

C 0.                                                                    HARM120 

C         -.0555                                                        HARM121 

C          .0260    -.0044                                              HARM122 

C          .0190    -.0033    -.0001                                    HARM123 

C         -.0139     .0076     .0019     .0010                          HARM124 

C          .0057    -.0018     .0009     .0032    -.0004                HARM125 

C         -.0026    -.0204     .0018     .0009     .0004     .0002      HARM126 

C   THE FOLLOWING SET OF GAUSS NORMALIZED COEFFICIENTS WERE CONVERTED   HARM127 

C   FROM THE SCHMIDT NORMALIZED COEFFICIENTS CALCULATED BY LINEARLY     HARM128 

C   EXTRAPOLATING TO EPOCH 1974 THE COEFFICIENTS PUBLISHED FOR EPOCH    HARM129 

C  1960 BY CAIN AND SWEENEY (1970).  (USES EARTH RADIUS = 6371.2)       HARM130 

C .000000                                                               HARM131 

C+.300953 +.020298                                                      HARM132 

C+.028106 -.05214   -.014435                                            HARM133 

C-.0308   +.06560   -.025252  -.006952                                  HARM134 

C-.041243 -.043956  -.016897  +.008021  -.002525                        HARM135 

C+.014742 -.037078  -.018906  +.002819  +.003656  +.000036              HARM136 

C-.006713 -.012234  -.004364  +.02137   +.001593  -.000072  +.00068     HARM137 

C .000000                                                               HARM138 

C .000000 -.057886                                                      HARM139 

C .000000 +.035942  +.001129                                            HARM140 

C .000000 +.011084  -.004421  +.001180                                  HARM141 

C .000000 -.010299  +.008794  -.000086  +.002256                        HARM142 

C .000000 -.003849  -.012615  +.007845  +.002207  -.000328              HARM143 

C .000000 +.003157  -.012670  -.009281  +.002286  -.000135  +.000243    HARM144 

         END                                                            HARM145-

