      PROGRAM ELEPOT(INPUT,OUTPUT,TAPE3,TAPE4)                          
C     THIS PROGRAM SOLVES FOR ELECTRIC POTENTIAL ON THE SURFACE OF A    
C     SPHERE.                                                           
      PARAMETER (NLMB=73,NTHET=37,NTHM=36,NTHETA=19,LLENGT=209364)      
      COMMON/COFPM/ISUM,SCTA(NTHET),FPARM(4)                            
      COMMON/CONDU/SOO(NLMB,NTHM),SOT(NLMB,NTHM),STO(NLMB,NTHM),        
     1  STT(NLMB,NTHM),DCO(NLMB,NTHM),DCT(NLMB,NTHM),DXO(NTHM),         
     2  DAY,UT,ALAMP,PHIP,ZPO,RREF,DLMBDA,DTHETA                        
      COMMON/SPACSV/ETHETA(NLMB,NTHET),ELMBDA(NLMB,NTHET),              
     1  DFAC(NLMB,NTHET),SGYY(NLMB,NTHET),PLTB(NLMB,NTHET)              
      DIMENSION PFRHS(NLMB,NTHET),PHI(NLMB,NTHET),FRHS(NLMB,NTHETA)     
      DIMENSION WORK(LLENGT),IPARM(11)                                  
      COMMON/XYCV/XX(100),YY(100)                                       
      CHARACTER*5 LABX,LABY                                             
      CHARACTER*56 LABL                                                 
      COMMON/CVLX/LABX                                                  
      COMMON/CVLY/LABY                                                  
      COMMON/CVLBL/LABL                                                 
      COMMON/STDCV/MXR,MX,MY,LTYPE, XC, XD, MAGX, MINX,xLABX, YC,       
     1YD, MAGY, MINY,xLABY,xLABL                                        
      EXTERNAL COFLT,BNDYC                                              
C***********************************************************************
C   NLMB=NUMBER OF LONGITUDE POINTS AROUND A LATITUDE CIRCLE,360 DEGREES
C   NTHET=NUMBER OF LATITUDE POINTS USED BY ELEPOT SOLVER (POLE TO POLE)
C   NTHETA=NUMBER OF LATITUDE POINTS USED BY ELEPOT SOLVER (POLE TO EQUA
C***********************************************************************
      CALL OPNGKS                                                       
      CALL DYNCOND                                                      
      PII=4.0*ATAN(1.0)                                                 
      RAD = 180./PII                                                    
      FPARM(1)= -UT*15./RAD -ALAMP                                      
      FPARM(2)=2.0*PII + FPARM(1)                                       
      FPARM(3)=0.                                                       
      FPARM(4)=90.*PII/180.                                             
      DTHETA=(FPARM(4)-FPARM(3))/(NTHETA-1)                             
      DLMBDA=(FPARM(2)-FPARM(1))/(NLMB-1)                               
      IPARM(1)=0                                                        
      IPARM(2)=NLMB                                                     
      IPARM(3)=NLMB                                                     
      IPARM(4)=NTHETA                                                   
      IPARM(5)=0                                                        
      IPARM(6)=0                                                        
      IPARM(7)=1                                                        
      IPARM(8)=2                                                        
      IPARM(9)=4                                                        
      IPARM(9)=2                                                        
      IPARM(10)=33                                                      
      IPARM(11)=LLENGT                                                  
      RE=6471.E+3                                                       
      A = RREF                                                          
      PHIINF = 0.                                                       
      IWHOLE = 0                                                        
      IF (NTHETA.NE.NTHET) IWHOLE = 1                                   
C     SET RHS                                                           
      DO 107 J=1,NTHET                                                  
      COLAT=FPARM(3)+(J-1)*DTHETA                                       
      YY(J) = 90. - COLAT*RAD                                           
      SCT=SIN(COLAT)                                                    
      CCT=COS(COLAT)                                                    
      SCTA(J)=SCT                                                       
  107 CONTINUE                                                          
      YY(NTHET+1) = YY(NTHET)                                           
      DO 117 L=1,NLMB                                                   
      XX(L) =  (L - 1)*DLMBDA*RAD/15.                                   
  117 CONTINUE                                                          
      XX(NLMB+1) = 25.                                                  
      DO 127 J=2,NTHETA                                                 
      JCON = NTHET - J + 1                                              
      SCM = .5*(SCTA(J-1) + SCTA(J))                                    
      SCP = SCM                                                         
      IF (J.LT.NTHETA) SCP = .5*(SCTA(J) + SCTA(J+1))                   
      DO 126 L=1,NLMB                                                   
      SGYY(L,J) =  .5*(SCM*STT(L,J-1) + SCP*STT(L,J))                   
      IF (IWHOLE.GT.0) SGYY(L,JCON) =  .5*(SCM*STT(L,JCON) +            
     1  SCP*STT(L,JCON-1))                                              
  126 CONTINUE                                                          
  127 CONTINUE                                                          
      XAN = 0.                                                          
      XAS = 0.                                                          
      NLMM = NLMB - 1                                                   
      DO 137 L=1,NLMM                                                   
      XAN = XAN + STT(L,1)                                              
      XAS = XAS + STT(L,NTHM)                                           
  137 CONTINUE                                                          
      X = .5*SCTA(2)/NLMM                                               
      XAN = XAN*X                                                       
      XAS = XAS*X                                                       
      DO 147 L=1,NLMB                                                   
      SGYY(L,1) = XAN                                                   
      IF (IWHOLE.GT.0) SGYY(L,NTHET) = XAS                              
  147 CONTINUE                                                          
      DO 157 J=1,NTHETA                                                 
      JCON = NTHET - J + 1                                              
      DO 156 L=1,NLMB                                                   
      PHI(L,J)=PHIINF                                                   
      FRHS(L,J) = SGYY(L,J)*SCTA(J)                                     
      IF (IWHOLE.GT.0) FRHS(L,J) = FRHS(L,J) + SCTA(J)*SGYY(L,JCON)     
  156 CONTINUE                                                          
  157 CONTINUE                                                          
C     CALL ASOUT(FRHS,4HFRHS,NLMB,NLMB,NTHETA,18,1)                     
      CALL ASOUT(FRHS,'FRHS',NLMB,NLMB,NTHETA,18,1)                     
      ISUM = IWHOLE                                                     
      CALL CROSEL(IPARM,FPARM,COFLT,BNDYC,WORK,FRHS,PHI,IER)            
      WRITE(6,4491) IER                                                 
 4491 FORMAT(1X,4HIER=,I4)                                              
C     CALL ASOUT(PHI,3HPHI,NLMB,NLMB,NTHETA,18,1)                       
      CALL ASOUT(PHI,'PHI',NLMB,NLMB,NTHETA,18,1)                       
      IPARM(1) = 1                                                      
      CALL JPARPOL(1,DMFAC,0,PHI)                                       
      CALL JPARPOL(NTHM,X,0,PHI)                                        
      IF (IWHOLE.GT.0) DMFAC = DMFAC + X                                
      KM = 1                                                            
      IF (IWHOLE.GT.0) KM = 2                                           
      NTM = NTHETA - 1                                                  
      NTMM = NTHETA - 2                                                 
      DO 503 L=1,NLMB                                                   
      LM1 = L - 1                                                       
      IF (L.EQ.1) LM1 = NLMB - 1                                        
      LP1 = L + 1                                                       
      IF (L.EQ.NLMB) LP1 = 2                                            
      FRHS(L,1) = 0.                                                    
      DO 503 J=2,NTHETA                                                 
      JJ = J                                                            
      FRHS(L,J) = 0.                                                    
      DO 390 K=1,KM                                                     
      DKODXO = (DCO(L,JJ) - DCO(L,JJ-1))/DTHETA                         
      DKTDXT = (DCT(LP1,JJ) + DCT(LP1,JJ-1) - DCT(LM1,JJ) - DCT(LM1,JJ -
     1  1))/(4.*DLMBDA)                                                 
      PFRHS(L,JJ) = DKODXO + DKTDXT                                     
      FRHS(L,J) = FRHS(L,J) + PFRHS(L,JJ)                               
  390 JJ = NTHET - J + 1                                                
  503 CONTINUE                                                          
C     CALL ASOUT(FRHS,4HFRHS,NLMB,NLMB,NTHETA,18,1)                     
      CALL ASOUT(FRHS,'FRHS',NLMB,NLMB,NTHETA,18,1)                     
      MXR = NLMB                                                        
      MX = NLMB                                                         
      MY = NTHETA                                                       
      LTYPE=1                                                           
      XC=XX(1)                                                          
      XD=XX(NLMB)                                                       
      YC=0.                                                             
      YD=90.                                                            
      MAGX=8                                                            
      MINX=3                                                            
      MAGY=3                                                            
      MINY=3                                                            
C     LABX= '   MLT'                                                    
C     LABY= '  GMLAT'                                                   
      LABX= ' MLT'                                                      
      LABY= 'GMLAT'                                                     
      LABL = ' '                                                        
C     LABL(2)=8H                                                        
C     LABL(3)=8H                                                        
C     LABL(4)=8H                                                        
C     LABL(5)=8H                                                        
C     LABL(6)=8H                                                        
C     LABL(7)=8H                                                        
C     LABL(1)=8H FRHS                                                   
      LABL(1:5)= ' FRHS'                                                
      CALL STDCON(FRHS,0.,0)                                            
C     COMPUTE SOLUTION                                                  
      CALL CROSEL(IPARM,FPARM,COFLT,BNDYC,WORK,FRHS,PHI,IER)            
      WRITE(6,4491) IER                                                 
C     CALL ASOUT(PHI,3HPHI,NLMB,NLMB,NTHETA,18,1)                       
      CALL ASOUT(PHI,'PHI',NLMB,NLMB,NTHETA,18,1)                       
      CALL JPARPOL(1,FAC,1,PHI)                                         
      CALL JPARPOL(NTHM,X,1,PHI)                                        
      IF (IWHOLE.GT.0) FAC = FAC + X                                    
      SYY = SGYY(1,1)                                                   
      IF (IWHOLE.GT.0) SYY = SYY + SGYY(1,NTHET)                        
      ALPAD = FAC/(DMFAC + SYY)                                         
      SM1 = 0.                                                          
      SM2 = 0.                                                          
      DO 507 L=1,NLMB                                                   
      FRHS(L,1) = 0.                                                    
  507 CONTINUE                                                          
      DO 660 J=2,NTHETA                                                 
      X = 1.                                                            
      IF (J.EQ.NTHETA) X = .5                                           
      DO 660 L=1,NLMB                                                   
      FRHS(L,J) = 0.                                                    
      JJ = J                                                            
      DO 660 K=1,KM                                                     
      XDM = ALPAD*SCTA(JJ)*SGYY(L,JJ)                                   
      IF (L.EQ.1) GO TO 659                                             
      SM1 = SM1 + X*(PFRHS(L,JJ))**2/SCTA(JJ)                           
      SM2 = SM2 + X*XDM*XDM/SCTA(J)                                     
  659 CONTINUE                                                          
      IF (J.EQ.NTHETA.AND.K.EQ.2) GO TO 658                             
      PFRHS(L,JJ) = PFRHS(L,JJ) - XDM                                   
  658 FRHS(L,J) = FRHS(L,J) + PFRHS(L,JJ)                               
  660 JJ = NTHET - J + 1                                                
C     CALL ASOUT(FRHS,4HFRHS,NLMB,NLMB,NTHETA,18,1)                     
      CALL ASOUT(FRHS,'FRHS',NLMB,NLMB,NTHETA,18,1)                     
      XDM = SQRT(SM2/SM1)                                               
      WRITE(6,661) XDM,ALPAD                                            
C 661 FORMAT (1H ,20HRMS CHANGE TO RHS = ,E10.3,10H    ALPAD=,E10.3)    
  661 FORMAT (1X,'RMS CHANGE TO RHS = ',E10.3,'    ALPAD=',E10.3)       
      CALL CROSEL(IPARM,FPARM,COFLT,BNDYC,WORK,FRHS,PHI,IER)            
      WRITE(6,4491) IER                                                 
C     CALL ASOUT(PHI,3HPHI,NLMB,NLMB,NTHETA,18,1)                       
      CALL ASOUT(PHI,'PHI',NLMB,NLMB,NTHETA,18,1)                       
      IF (NTHETA.EQ.NTHET) GO TO 714                                    
      DO 713 J=1,NTHETA                                                 
      JCON = NTHET - J + 1                                              
      DO 713 I=1,NLMB                                                   
  713 PHI(I,JCON) = PHI(I,J)                                            
  714 CONTINUE                                                          
C     CALL ASOUT(PHI,3HPHI,NLMB,NLMB,NTHET ,18,1)                       
      CALL ASOUT(PHI,'PHI',NLMB,NLMB,NTHET ,18,1)                       
C     LABL(1)=8HPOTENTIA                                                
C     LABL(2)=8HL (VOLTS                                                
C     LABL(3)=8H)                                                       
      LABL(1:24) = ' '                                                  
C89/6/26 At this point, PHI holds the electric potential, in Volts,     
C over a uniformly spaced                                               
C grid NLMB X NTHET (= MXR X MY, I think) at approximately              
C 140 km altitude, EXCEPT at the equator, where the altitude is 90 km.  
C At two later point in the program PHI is again plotted out, once at   
C 300 km and another time at 90 km.  For both of these later plots,     
C the values of PHI in the array are the same, but the latitudinal      
C locations of the grid points are changed to reflect the mapping along 
C dipole field lines.                                                   
      CALL STDCON(PHI,500.,0)                                           
      DIST = 2.*A*DTHETA                                                
C     DERIVE AND PROCESS VECTORS FROM PHI                               
      DO 717 J=1,NTHETA                                                 
      JCON = NTHET - J + 1                                              
      DSTY = 2.*A*DLMBDA*SCTA(J)                                        
      JM = J - 1                                                        
      LA = 0                                                            
      JL = J                                                            
      LS = 1                                                            
      IF (J.EQ.1) THEN                                                  
      DSTY = DIST                                                       
      JL = 2                                                            
      JM = 2                                                            
      LA = NLMB/2                                                       
      LS = NLMB/4                                                       
      END IF                                                            
      JP = J + 1                                                        
      IF (J.EQ.NTHETA) JP = NTHETA - 1                                  
      DO 716 L=1,NLMB                                                   
      LJ = L + LA                                                       
      IF (LJ.GT.NLMB) LJ = LJ - NLMB + 1                                
      LM = L - LS                                                       
      IF (LM.LT.1) LM = LM + NLMB - 1                                   
      LP = L + LS                                                       
      IF (LP.GT.NLMB) LP = LP - NLMB + 1                                
      ETHETA(L,J) = (PHI(LJ,JM) - PHI(L,JP))/DIST                       
      ELMBDA(L,J) = (PHI(LM,JL) - PHI(LP,JL))/DSTY                      
      IF (NTHETA.NE.NTHET) THEN                                         
      ETHETA(L,JCON) = - ETHETA(L,J)                                    
      ELMBDA(L,JCON) = ELMBDA(L,J)                                      
      END IF                                                            
  716 CONTINUE                                                          
  717 CONTINUE                                                          
C     LABL(1)=8HELECTRIC                                                
C     LABL(2)=8H FIELD A                                                
C     LABL(3)=8HT RREF                                                  
      LABL(1:24) = 'ELECTRIC FIELD AT RREF'                             
      DO 751 J=1,NTHM                                                   
  751 YY(J) = 90. - (FLOAT(J) - .5)*DTHETA*RAD                          
      YD = YY(1) + .5*DTHETA*RAD                                        
      YC = YY(NTHM) - .5*DTHETA*RAD                                     
      DIST = A*DTHETA                                                   
      DO 750 L=1,NLMB                                                   
      DO 750 J=1,NTHM                                                   
      ELMBDA(L,J) = .5*(ELMBDA(L,J) + ELMBDA(L,J+1))                    
  750 ETHETA(L,J) = (PHI(L,J) - PHI(L,J+1))/DIST                        
      DO 753 J=1,NLMB,2                                                 
      DO 753 K=1,NTHM                                                   
      PLTB(J,K) = 1.E30                                                 
      IF (J.NE.NLMB) PLTB(J+1,K) = -ETHETA(J+1,K)                       
  753 CONTINUE                                                          
      MAGY = 6                                                          
      MY = NTHM                                                         
      CALL STDVCT(ELMBDA,PLTB,.048)                                     
      DO  767 J=1,NTHM                                                  
      SCT = .5*(SCTA(J) + SCTA(J+1))                                    
      DO 766 L=1,NLMB                                                   
      ET=ETHETA(L,J)                                                    
      EL = ELMBDA(L,J)*SCT                                              
      ETHETA(L,J) = (SOO(L,J)*ET + SOT(L,J)*EL)/SCT                     
      ELMBDA(L,J) =  STO(L,J)*ET + STT(L,J)*EL                          
  766 CONTINUE                                                          
  767 CONTINUE                                                          
C     LABL(1) = 8HELECTRIC                                              
C     LABL(2) = 8H-FIELD-D                                              
C     LABL(3) = 8HRIVEN CU                                              
C     LABL(4) = 8HRRENT                                                 
      LABL(1:32) = 'ELECTRIC-FIELD-DRIVEN CURRENT'                      
      DO 854 K=1,NTHM                                                   
      DO 853 J=2,NLMB,2                                                 
  853 PLTB(J,K) = -ETHETA(J,K)                                          
  854 CONTINUE                                                          
      CALL STDVCT(ELMBDA,PLTB,2.4)                                      
C  CALCULATE DOWNWARD CURRENT DENSITY                                   
      ISUM = 0                                                          
      DO 867 L=1,NLMB                                                   
      XLMBDA = (L-1)*DLMBDA + FPARM(1)                                  
      LM = L - 1                                                        
      IF (LM.LT.1) LM = NLMB - 1                                        
      LP = L + 1                                                        
      IF (LP.GT.NLMB) LP = 2                                            
      DO 866 J=2,NTHM                                                   
      X = RREF*RREF*SCTA(J)                                             
      COLAT = (J - 1)*DTHETA                                            
      CALL COFLT(XLMBDA,COLAT,CLL,CLT,CTT,CL,CT,CF)                     
      PLL = (PHI(LM,J) - 2.*PHI(L,J) + PHI(LP,J))/(DLMBDA*DLMBDA)       
      PLT = (PHI(LP,J+1) - PHI(LP,J-1) - PHI(LM,J+1) + PHI(LM,J-1))/    
     1  (4.*DLMBDA*DTHETA)                                              
      PTT = (PHI(L,J-1) - 2.*PHI(L,J) + PHI(L,J+1))/(DTHETA*DTHETA)     
      PT = (PHI(L,J+1) - PHI(L,J-1))/(2.*DTHETA)                        
      PL = (PHI(LP,J) - PHI(LM,J))/(2.*DLMBDA)                          
      DFAC(L,J) = (PFRHS(L,J) - CLL*PLL - CLT*PLT - PTT*CTT - CL*PL -   
     1  CT*PT)/X                                                        
  866 CONTINUE                                                          
  867 CONTINUE                                                          
      CALL JPARPOL(1,FAC,1,PHI)                                         
      FAC = (FAC - ALPAD*SGYY(1,1))/(RREF*RREF)                         
      CALL JPARPOL(NTHM,DMFAC,1,PHI)                                    
      DMFAC = (DMFAC - ALPAD*SGYY(1,NTHET))/(RREF*RREF)                 
      DO 870 L=1,NLMB                                                   
      DFAC(L,1) = FAC                                                   
      IF (IWHOLE.GT.0) DFAC(L,NTHET) = DMFAC                            
      DFAC(L,NTHETA)=0.                                                 
  870 CONTINUE                                                          
C COMPUTE TOTAL IONOSPHERIC CURRENT                                     
      SQALTR = SQRT(RREF/6.671E6)                                       
      DO 880 J=1,NTM                                                    
      JCON = NTHET - J                                                  
      TH2 = (FLOAT(J) - .5)*DTHETA                                      
      SCT = SIN(TH2)                                                    
      TH1 = ASIN(SQALTR*SCT)                                            
      DELTH = TH2 - TH1                                                 
      IF (DELTH.GT.DTHETA) GO TO 875                                    
      JSV = J                                                           
      JCSV = NTHET - J                                                  
      DXM = DELTH*A/(2.*DTHETA)                                         
      FXP = DXM*(DTHETA - DELTH)                                        
      FXM = DXM*(DTHETA + DELTH)                                        
C CONTRIBUTION TO NORTHWARD CURRENT FROM FAC BETWEEN RREF AND 300 KM    
      DO 874 L=1,NLMB                                                   
      SVJN = FXM*DFAC(L,J) + FXP*DFAC(L,J+1)                            
      ETHETA(L,J) = ETHETA(L,J) + DCO(L,J)/(A*SCT) - SVJN               
      IF (NTHET.NE.NTHETA) ETHETA(L,JCON)=ETHETA(L,JCON) + DCO(L,JCON)/ 
     1  (A*SCT) - SVJN                                                  
  874 CONTINUE                                                          
      GO TO 879                                                         
  875 DXM = JCSV - JSV                                                  
      FXM = (JCSV - J)/DXM                                              
      FXP = (J - JSV)/DXM                                               
      DO 878 L=1,NLMB                                                   
      X = -ETHETA(L,JSV)                                                
      IF (NTHET.NE.NTHETA) X = ETHETA(L,JCSV)                           
      ETHETA(L,J) = FXM*ETHETA(L,JSV) + FXP*X                           
      IF (NTHET.NE.NTHETA) ETHETA(L,JCON) = FXP*ETHETA(L,JSV) +         
     1  FXM*ETHETA(L,JCSV)                                              
  878 CONTINUE                                                          
  879 CONTINUE                                                          
      DO 880 L=1,NLMB                                                   
      ELMBDA(L,J) = ELMBDA(L,J) + DCT(L,J)/A                            
      IF (NTHET.NE.NTHETA) ELMBDA(L,JCON)=ELMBDA(L,JCON) + DCT(L,JCON)/A
  880 CONTINUE                                                          
      CALL GSCLIP(0)                                                    
      DO 881 J=1,NTHM                                                   
      DO 881 L=2,NLMB,2                                                 
  881 PLTB(L,J) = -ETHETA(L,J)                                          
C     LABL(1) = 8HTOTAL HO                                              
C     LABL(2) = 8HRIZONTAL                                              
C     LABL(3) = 8H CURRENT                                              
C     LABL(4) = 8H BETWEEN                                              
C     LABL(5) = 8H 90 KM A                                              
C     LABL(6) = 8HND 300 K                                              
C     LABL(7) = 8HM                                                     
      LABL = 'TOTAL HORIZONTAL CURRENT BETWEEN 0 KM AND 300 KM'         
      CALL STDVCT(ELMBDA,PLTB,2.4)                                      
C     LABL(1)=8HDOWNWARD                                                
C     LABL(2)=8H CURRENT                                                
C     LABL(3)=8H AT RREF                                                
C     LABL(4)=8H (A/M2)                                                 
C     LABL(5)=8H                                                        
C     LABL(6)=8H                                                        
C     LABL(7)=8H                                                        
      LABL = 'DOWNWARD CURRENT AT RREF (A/M2)'                          
      DO 882 J=1,NTHET                                                  
  882 YY(J) = 90. - (J - 1)*DTHETA*RAD                                  
      YD = YY(1)                                                        
      YC = YY(NTHET)                                                    
      MY = NTHET                                                        
      CALL STDCON(DFAC,2.E-9,0)                                         
C  CONVERT LATITUDES AND DFAC ALONG FIELD LINES TO 300 KM               
      DO 890 J=1,NTM                                                    
      DXM = SCTA(J)/SQALTR                                              
      IF (DXM.GE.1.) GO TO 891                                          
      REFLAT = RAD*ACOS(SCTA(J))                                        
      YY(J) = RAD*ACOS(DXM)                                             
      X = YY(J)/REFLAT                                                  
      DO 890 L=1,NLMB                                                   
      PLTB(L,J) = PHI(L,J)                                              
  890 DFAC(L,J) = X*DFAC(L,J)                                           
  891 J = J - 1                                                         
      MY = 2*J                                                          
      PHIEQ = 0.                                                        
      DO 892 L=2,NLMB                                                   
  892 PHIEQ = PHIEQ + PHI(L,J)                                          
      PHIEQ = PHIEQ/FLOAT(NLMB - 1)                                     
      DO 893 K=1,J                                                      
      JCON = MY - K + 1                                                 
      YY(JCON) = -YY(K)                                                 
      DO 893 L=1,NLMB                                                   
      PLTB(L,K) = PLTB(L,K) - PHIEQ                                     
      IF (NTHET.EQ.NTHETA) GO TO 893                                    
      PLTB(L,JCON) = PLTB(L,K)                                          
      DFAC(L,JCON) = -DFAC(L,K)                                         
  893 CONTINUE                                                          
      WRITE(6,8932) (YY(I),I=1,J)                                       
C8932 FORMAT (1H0,18F5.1)                                               
 8932 FORMAT (1X/18F5.1)                                                
C     LABL(1) = 8HDOWNWARD                                              
C     LABL(2) = 8H CURRENT                                              
C     LABL(3) = 8H AT 300                                               
C     LABL(4) = 8HKM                                                    
      LABL(1:32) = 'DOWNWARD CURRENT AT 300 KM'                         
      CALL STDCON(DFAC,2.E-9,0)                                         
C     LABL(1) = 8HPOTENTIA                                              
C     LABL(2) = 8HL AT 300                                              
C     LABL(3) = 8H KM (VOL                                              
C     LABL(4) = 8HTS)                                                   
      LABL(1:32) = 'POTENTIAL AT 300 KM (VOLTS)'                        
      CALL STDCON(PLTB,500.,0)                                          
C  PLOT ELECTRODYNAMIC DRIFT CURVES                                     
      CALL VIPLT(PLTB,270.)                                             
C  PLOT POTENTIAL FROM EMPIRICAL MODEL                                  
      I = MY/2                                                          
      DO 8935 J=1,I                                                     
      DO 8935 L=1,NLMB                                                  
      XMLON = (FPARM(1) + (L-1)*DLMBDA)*RAD                             
      UTS = UT + (ALAMP+69.)/15.                                        
 8935 CALL EFIELD(YY(J),XMLON,DAY,UTS,0,0,PLTB(L,J),X,DXM)              
      X = 0.                                                            
      DO 8936 L=2,NLMB                                                  
 8936 X = X + PLTB(L,I)                                                 
      X = X/FLOAT(NLMB - 1)                                             
      DO 8937 J=1,I                                                     
      JCON = MY - J + 1                                                 
      DO 8937 L=1,NLMB                                                  
      PLTB(L,J) = PLTB(L,J) - X                                         
      IF (NTHET.NE.NTHETA) PLTB(L,JCON) = PLTB(L,J)                     
 8937 CONTINUE                                                          
C     LABL(1) = 8HPOTENTIA                                              
C     LABL(2) = 8HL AT 300                                              
C     LABL(3) = 8H KM (VOL                                              
C     LABL(4) = 8HTS) (EMP                                              
C     LABL(5) = 8HIRICAL)                                               
C     LABL(6) = 8H                                                      
C     LABL(7) = 8H                                                      
      LABL = 'POTENTIAL AT 300 KM (VOLTS) (EMPIRICAL)'                  
      CALL STDCON(PLTB,500.,0)                                          
C  CONVERT LATITUDES ALONG FIELD LINES TO RBOT                          
      RBOT = 6.371E6 + ZPO*1.E3                                         
      X = SQRT(RBOT/RREF)                                               
      YY(NTHETA) = 0.                                                   
      DO 894 J=2,NTM                                                    
      JCON = NTHET -J + 1                                               
      YY(J) = RAD*ACOS(X*SCTA(J))                                       
      IF (NTHET.NE.NTHETA) YY(JCON) = -YY(J)                            
  894 CONTINUE                                                          
      DO 895 J=1,NTHET                                                  
      DO 895 L=1,NLMB                                                   
  895 PHI(L,J) = PHI(L,J) - PHIEQ                                       
C     LABL(2) = 8HL AT  90                                              
C     LABL(4) = 8HTS)                                                   
C     LABL(5) = 8H                                                      
      LABL(15:40) = '90 KM (VOLTS)'                                     
      MY = NTHET                                                        
      CALL STDCON(PHI,500.,0)                                           
      CALL EQVCUR(ETHETA,ELMBDA)                                        
      CALL CLSGKS                                                       
      END                                                               
      SUBROUTINE DYNCOND                                                
      PARAMETER(NLAT=36,NTHM=36,NLONG=73,NLMB=73,NTHET=37,BEQ=2.82E-5)  
C***********************************************************************
C   NLAT=NUMBER OF GEOGRAPHIC LATITUDE POINTS IN IONOSPHERE CALCS       
C   NTHM=NTHET-1,WHERE NTHET=NUMBER OF GEOMAGNETIC LATITUDE POINTS      
C    IN DYNAMO CALCULATIONS.  NTHM MUST BE EVEN.                        
C   NLONG=NUMBER OF GEOGRAPHIC LONGITUDE POINTS IN IONOSPHERE COMPS     
C   NLMB=NUMBER OF GEOMAGNETIC LONGITUDE POINTS IN DYNAMO COMPUTATIONS  
C   BEQ=EQUATORIAL MAGNETIC FIELD STRENGTH IN TESLA                     
C***********************************************************************
      COMMON/CONDU/SOO(NLMB,NTHM),SOT(NLMB,NTHM),STO(NLMB,NTHM),        
     1  STT(NLMB,NTHM),DCO(NLMB,NTHM),DCT(NLMB,NTHM),DXO(NTHM),         
     2  DAY,UT,ALAMP,PHIP,ZPO,RREF,DLMBDA,DTHETA                        
      DIMENSION FLIS1(NLONG,NLAT),FLIS2(NLONG,NLAT),FLIUS1(NLONG,NLAT), 
     1  FLIUS2(NLONG,NLAT),ESOO(NLONG,2),ESOT(NLONG,2),ESTO(NLONG,2),   
     2  ESTT(NLONG,2),EDCO(NLONG,2),EDCT(NLONG,2),LATSV(NLONG,2)        
     3  ,FLIAS1(NLONG,NLAT),FLIAS2(NLONG,NLAT),GMLON(NLONG,NLAT)        
      COMMON/SPACSV/UH(NLMB,NTHET),VH(NLMB,NTHET),UP(NLMB,NTHET),       
     1  VP(NLMB,NTHET),RAT(NLMB,NTHET)                                  
      COMMON/XYCV/XX(100),YY(100)                                       
      CHARACTER*5 LABX,LABY                                             
      CHARACTER*56 LABL                                                 
      COMMON/CVLX/LABX                                                  
      COMMON/CVLY/LABY                                                  
      COMMON/CVLBL/LABL                                                 
      COMMON/STDCV/MXR,MX,MY,LTYPE, XC, XD, MAGX, MINX,xLABX, YC,       
     1YD, MAGY, MINY,xLABY,xLABL                                        
      DTR=ATAN(1.)/45.                                                  
      PI=3.1415926535898                                                
      READ(3) DAY,UT,ALAMP,PHIP,ZPO,RREF,FLIS1,FLIS2,FLIUS1,FLIUS2      
     1  ,FLIAS1,FLIAS2,ESOO,ESOT,ESTO,ESTT,EDCO,EDCT,LATSV,GMLON        
      WRITE (6,12344) DAY,UT,ALAMP,PHIP,ZPO,RREF                        
C2344 FORMAT(1H ,6E11.2)                                                
12344 FORMAT(1X,6E11.2)                                                 
      RAD = 1./DTR                                                      
      DLMB = 2.*PI/(NLMB-1)                                             
      DLON = 2.*PI/(NLONG-1)                                            
      DLAT = PI/FLOAT(NLAT)                                             
      DTHT = PI/FLOAT(NTHM)                                             
      DO 5 K=1,NTHM                                                     
    5 DXO(K) = DTHT                                                     
      XX(NLMB+1) = 25.                                                  
      YY(NTHM+1) = -90.                                                 
      DO 4000 K=1,NTHM                                                  
      RMLAT = (.5 + NTHM/2.- K)*DTHT                                    
      YY(K) = RMLAT/DTR                                                 
      COLAT= (K-.5)*DTHT                                                
      CTHT = COS(COLAT)                                                 
      CTHS = CTHT*CTHT                                                  
      BOBE = SQRT(1. + 3.*CTHS)                                         
      BWPM = BEQ*BOBE                                                   
      SI = ABS(2.*CTHT/BOBE)                                            
      KODE = 0                                                          
      IF (ABS(FLOAT(NTHM+1)/2.E0 - K).LT.1.) KODE = 1                   
      DO 4000 J=1,NLMB                                                  
      RMLON = -UT*15.*DTR + (J-1)*DLMB - ALAMP                          
      CALL MTG(RLATG,RLONG,PHIP,ALAMP,RMLAT,RMLON)                      
      XLON = (RLONG + PI)/DLON + 1.                                     
      LON = XLON                                                        
      LON = MAX0(1,LON)                                                 
      LON = MIN0(NLONG-1,LON)                                           
      LONAL = LON                                                       
      LONAR = LON + 1                                                   
      LONBL = LON                                                       
      LONBR = LON + 1                                                   
      XLON = XLON - LON                                                 
      XLAT =  (RLATG + .5*PI)/DLAT + .5                                 
      LAT = XLAT                                                        
      LATB = LAT                                                        
      LATA = LAT + 1                                                    
      FAL = XLAT - LAT                                                  
      IF (LAT.GT.0) GO TO 450                                           
      LATB = 1                                                          
      LONBL = LONBL + (NLONG-1)/2                                       
      IF (LONBL.GT.NLONG) LONBL = LONBL - NLONG + 1                     
      LONBR = LONBL - 1                                                 
  450 CONTINUE                                                          
      IF (LAT.LT.NLAT) GO TO 460                                        
      LATA = NLAT                                                       
      LONAL = LONAL + (NLONG-1)/2                                       
      IF (LONAL.GT.NLONG) LONAL = LONAL - NLONG + 1                     
      LONAR = LONAL - 1                                                 
  460 CONTINUE                                                          
      XAL = FAL*(1.-XLON)                                               
      XAR = FAL*XLON                                                    
      XBL = (1.-FAL)*(1.-XLON)                                          
      XBR = (1.-FAL)*XLON                                               
      FLISP = XAL*FLIS1(LONAL,LATA) + XAR*FLIS1(LONAR,LATA) +           
     1        XBL*FLIS1(LONBL,LATB) + XBR*FLIS1(LONBR,LATB)             
      FLISH = XAL*FLIS2(LONAL,LATA) + XAR*FLIS2(LONAR,LATA) +           
     1        XBL*FLIS2(LONBL,LATB) + XBR*FLIS2(LONBR,LATB)             
      IF (ABS(RMLAT).GT.1.4) GO TO 462                                  
      FLIUSP = XAL*FLIUS1(LONAL,LATA) + XAR*FLIUS1(LONAR,LATA) +        
     1        XBL*FLIUS1(LONBL,LATB) + XBR*FLIUS1(LONBR,LATB)           
      FLIUSH = XAL*FLIUS2(LONAL,LATA) + XAR*FLIUS2(LONAR,LATA) +        
     1        XBL*FLIUS2(LONBL,LATB) + XBR*FLIUS2(LONBR,LATB)           
      FLIASP = XAL*FLIAS1(LONAL,LATA) + XAR*FLIAS1(LONAR,LATA) +        
     1        XBL*FLIAS1(LONBL,LATB) + XBR*FLIAS1(LONBR,LATB)           
      FLIASH = XAL*FLIAS2(LONAL,LATA) + XAR*FLIAS2(LONAR,LATA) +        
     1        XBL*FLIAS2(LONBL,LATB) + XBR*FLIAS2(LONBR,LATB)           
      GO TO 463                                                         
  462 ROT = GMLON(LONAL,LATA)/RAD - RMLON                               
      CAL = COS(ROT)                                                    
      SAL = SIN(ROT)                                                    
      ROT = GMLON(LONAR,LATA)/RAD - RMLON                               
      CAR = COS(ROT)                                                    
      SAR = SIN(ROT)                                                    
      ROT = GMLON(LONBL,LATB)/RAD - RMLON                               
      CBL = COS(ROT)                                                    
      SBL = SIN(ROT)                                                    
      ROT = GMLON(LONBR,LATB)/RAD - RMLON                               
      CBR = COS(ROT)                                                    
      SBR = SIN(ROT)                                                    
      DAL = FLIUS1(LONAL,LATA)*CAL - FLIAS1(LONAL,LATA)*SAL             
      DAR = FLIUS1(LONAR,LATA)*CAR - FLIAS1(LONAR,LATA)*SAR             
      DBL = FLIUS1(LONBL,LATB)*CBL - FLIAS1(LONBL,LATB)*SBL             
      DBR = FLIUS1(LONBR,LATB)*CBR - FLIAS1(LONBR,LATB)*SBR             
      FLIUSP = XAL*DAL + XAR*DAR + XBL*DBL + XBR*DBR                    
      DAL = FLIAS1(LONAL,LATA)*CAL + FLIUS1(LONAL,LATA)*SAL             
      DAR = FLIAS1(LONAR,LATA)*CAR + FLIUS1(LONAR,LATA)*SAR             
      DBL = FLIAS1(LONBL,LATB)*CBL + FLIUS1(LONBL,LATB)*SBL             
      DBR = FLIAS1(LONBR,LATB)*CBR + FLIUS1(LONBR,LATB)*SBR             
      FLIASP = XAL*DAL + XAR*DAR + XBL*DBL + XBR*DBR                    
      DAL = FLIUS2(LONAL,LATA)*CAL - FLIAS2(LONAL,LATA)*SAL             
      DAR = FLIUS2(LONAR,LATA)*CAR - FLIAS2(LONAR,LATA)*SAR             
      DBL = FLIUS2(LONBL,LATB)*CBL - FLIAS2(LONBL,LATB)*SBL             
      DBR = FLIUS2(LONBR,LATB)*CBR - FLIAS2(LONBR,LATB)*SBR             
      FLIUSH = XAL*DAL + XAR*DAR + XBL*DBL + XBR*DBR                    
      DAL = FLIAS2(LONAL,LATA)*CAL + FLIUS2(LONAL,LATA)*SAL             
      DAR = FLIAS2(LONAR,LATA)*CAR + FLIUS2(LONAR,LATA)*SAR             
      DBL = FLIAS2(LONBL,LATB)*CBL + FLIUS2(LONBL,LATB)*SBL             
      DBR = FLIAS2(LONBR,LATB)*CBR + FLIUS2(LONBR,LATB)*SBR             
      FLIASH = XAL*DAL + XAR*DAR + XBL*DBL + XBR*DBR                    
  463 FLIDCA = (FLIUSP + FLIASH)*BWPM                                   
      FLIDCP = (FLIUSH - FLIASP)*BWPM                                   
      RAT(J,K) = ALOG10(FLISH/FLISP)                                    
      IF (MOD(J,2).NE.0) GO TO 465                                      
      UP(J,K) = FLIUSP/FLISP                                            
      UH(J,K) = FLIUSH/FLISH                                            
      VP(J,K) = FLIASP/FLISP/SI                                         
      VH(J,K) = FLIASH/FLISH/SI                                         
      GO TO 467                                                         
  465 UP(J,K) = 0.                                                      
      UH(J,K) = 0.                                                      
      VP(J,K) = 0.                                                      
      VH(J,K) = 0.                                                      
  467 CONTINUE                                                          
      IF (KODE.EQ.1) GO TO 480                                          
      SOO(J,K) = FLISP/SI                                               
      SOT(J,K) = FLISH                                                  
      STO(J,K) = -FLISH                                                 
      STT(J,K) = FLISP*SI                                               
      DCO(J,K) = -(XAL*FLIDCA+XAR*FLIDCA+                               
     1  XBL*FLIDCA+XBR*FLIDCA)                                          
      DCT(J,K) =(XAL*FLIDCP+XAR*FLIDCP+                                 
     1  XBL*FLIDCP+XBR*FLIDCP)*SI                                       
      GO TO 3950                                                        
  480 CONTINUE                                                          
      FAL = (XLAT - LATSV(LON,1))/FLOAT(LATSV(LON,2) - LATSV(LON,1))    
      FAR = (XLAT - LATSV(LONAR,1))/FLOAT(LATSV(LONAR,2)-LATSV(LONBR,1))
      XAL = FAL*(1. - XLON)                                             
      XAR = FAR*XLON                                                    
      XBL = (1. - FAL)*(1. - XLON)                                      
      XBR = (1. - FAR)*XLON                                             
      SOO(J,K) = XAL*ESOO(LON,2) + XAR*ESOO(LONAR,2)                    
     1         + XBL*ESOO(LON,1) + XBR*ESOO(LONBR,1)                    
      SOT(J,K) = XAL*ESOT(LON,2) + XAR*ESOT(LONAR,2)                    
     1         + XBL*ESOT(LON,1) + XBR*ESOT(LONBR,1)                    
      STO(J,K) = XAL*ESTO(LON,2) + XAR*ESTO(LONAR,2)                    
     1         + XBL*ESTO(LON,1) + XBR*ESTO(LONBR,1)                    
      STT(J,K) = XAL*ESTT(LON,2) + XAR*ESTT(LONAR,2)                    
     1         + XBL*ESTT(LON,1) + XBR*ESTT(LONBR,1)                    
      DCO(J,K) = XAL*EDCO(LON,2) + XAR*EDCO(LONAR,2)                    
     1         + XBL*EDCO(LON,1) + XBR*EDCO(LONBR,1)                    
      DCT(J,K) = XAL*EDCT(LON,2) + XAR*EDCT(LONAR,2)                    
     1         + XBL*EDCT(LON,1) + XBR*EDCT(LONBR,1)                    
      XX(J) = 24.*(J-1)/FLOAT(NLMB-1)                                   
 3950 CONTINUE                                                          
      IF (CTHT.GT.0.) GO TO 4000                                        
      SOT(J,K) = -SOT(J,K)                                              
      STO(J,K) = -STO(J,K)                                              
      DCO(J,K) = -DCO(J,K)                                              
      VH(J,K) = -VH(J,K)                                                
      VP(J,K) = -VP(J,K)                                                
 4000 CONTINUE                                                          
      NTHMO2 = NTHM/2                                                   
      NTHETA = NTHMO2 + 1                                               
      K1 = NTHMO2 - 1                                                   
      K2 = NTHETA + 1                                                   
      DO 4003 J=2,NLMB,2                                                
      VP(J,NTHMO2) = (VP(J,K2) + 2.*VP(J,K1))/3.E0                      
      VP(J,NTHETA) = (2.*VP(J,K2) + VP(J,K1))/3.E0                      
      VH(J,NTHMO2) = (VH(J,K2) + 2.*VH(J,K1))/3.E0                      
 4003 VH(J,NTHETA) = (2.*VH(J,K2) + VH(J,K1))/3.E0                      
      KM = NTHM                                                         
C     CALL ASOUT(SOO,3HSOO,NLMB,NLMB,KM,18,1)                           
      CALL ASOUT(SOO,'SOO',NLMB,NLMB,KM,18,1)                           
      CALL ASOUT(SOT,'SOT',NLMB,NLMB,KM,18,1)                           
      CALL ASOUT(STO,'STO',NLMB,NLMB,KM,18,1)                           
      CALL ASOUT(STT,'STT',NLMB,NLMB,KM,18,1)                           
      CALL ASOUT(DCO,'DCO',NLMB,NLMB,KM,18,1)                           
      CALL ASOUT(DCT,'DCT',NLMB,NLMB,KM,18,1)                           
      MXR = NLMB                                                        
      MX = NLMB                                                         
      MY = NTHM                                                         
      LTYPE = 1                                                         
      MAGX= 8                                                           
      MINX = 3                                                          
      MAGY = 6                                                          
      MINY = 3                                                          
      XC = XX(1)                                                        
      XD = XX(NLMB)                                                     
      YC = YY(NTHM) - 90./FLOAT(NTHM)                                   
      YD = YY(1) + 90./FLOAT(NTHM)                                      
C     LABX = 8H   MLT                                                   
C     LABY = 8H  GMLAT                                                  
C     LABL(1) = 8HLOG10(HA                                              
C     LABL(2) = 8HLL TO PE                                              
C     LABL(3) = 8HDERSEN C                                              
C     LABL(4) = 8HONDUCTAN                                              
C     LABL(5) = 8HCE RATIO                                              
C     LABL(6) = 8H)                                                     
C     LABL(7) = 8H                                                      
C     LABX = '   MLT'                                                   
C     LABY = '  GMLAT'                                                  
      LABX = ' MLT'                                                     
      LABY = 'GMLAT'                                                    
      LABL = 'LOG10(HALL TO PEDERSEN CONDUCTANCE RATIO)'                
      CALL STDCON(RAT,.1,0)                                             
C     LABL(1) = 8HSTT  (SI                                              
C     LABL(2) = 8HEMENS)                                                
C     LABL(3) = 8H                                                      
C     LABL(4) = 8H                                                      
C     LABL(5) = 8H                                                      
C     LABL(6) = 8H                                                      
      LABL = 'STT  (SIEMENS)'                                           
      CALL STDCON(STT,5.,0)                                             
C     LABL(1) = 8HPEDERSEN                                              
C     LABL(2) = 8H-WEIGHTE                                              
C     LABL(3) = 8HD WIND                                                
      LABL = 'PEDERSEN-WEIGHTED WIND'                                   
      CALL STDVCT(UP,VP,2400.)                                          
C     LABL(1) = 8H    HALL                                              
      LABL(1:8) = '    HALL'                                            
      CALL STDVCT(UH,VH,2400.)                                          
C  REVERSE DCO SIGN FOR VECTOR PLOT                                     
      DO 4355 K=1,KM                                                    
      DO 4355 J=1,NLMB,2                                                
      IF (J.NE.NLMB) DCO(J+1,K) = -DCO(J+1,K)                           
C  MAKE EVERY OTHER POINT SMALL ENOUGH NOT TO BE PLOTTED                
      DCT(J,K) = 1.E-4*DCT(J,K)                                         
 4355 DCO(J,K) = -1.E-4*DCO(J,K)                                        
C     LABL(1) = 8HWIND-DRI                                              
C     LABL(2) = 8HVEN CURR                                              
C     LABL(3) = 8HENT                                                   
      LABL = 'WIND-DRIVEN CURRENT'                                      
      CALL STDVCT(DCT,DCO,2.4)                                          
      THM = 0.                                                          
      DO 4010 K=1,KM                                                    
      THP = THM + DXO(K)                                                
      X = .5*(THM + THP)                                                
      X = SIN(X)                                                        
      DO 4005 J=1,NLMB                                                  
      SOO(J,K) = SOO(J,K)*X                                             
      STT(J,K) = STT(J,K)/X                                             
      DCT(J,K) = DCT(J,K)*RREF                                          
 4005 DCO(J,K) =-DCO(J,K)*X*RREF                                        
      DO 4008 J=1,NLMB,2                                                
      DCO(J,K) = 1.E4*DCO(J,K)                                          
 4008 DCT(J,K) = 1.E4*DCT(J,K)                                          
 4010 THM = THP                                                         
      RETURN                                                            
      END                                                               
      SUBROUTINE GTM(RLATG,RLONG,RLATMP,RLONMP,RLATM,RLONM,DIP,DEC)     
C     ****     TRANSFORMS GEOGRAPHIC COORDINATES TO GEOMAGNETIC,        
C     ****     CALCULATES DIP AND DECLINATION.                          
C     ****                                                              
C     ****     RLATG, RLONG --  GEOGRAPHIC LATITUDE AND LONGITUDE (RADS)
C     ****                                                              
C     ****     RLATMP, RLONMP --  GEOGRAPHIC COORDINATES OF GEOMAGNETIC 
C     ****     NORTH POLE. (RADIANS)                                    
C     ****                                                              
C     ****     RLATM, RLONM --  OUTPUT GEOMAGNETIC COORDINATES (RADS)   
C     ****                                                              
C     ****     DIP, DEC -- OUTPUT DIP AND DECLINATION.                  
C     ****                                                              
      DATA E/1.E-10/                                                    
      SINLAM=SIN(RLATMP)                                                
      COSLAM=COS(RLATMP)                                                
      W1=SIN(RLATG)                                                     
      W2=COS(RLATG)                                                     
      W3=SIN(RLONG+E-RLONMP)                                            
      W4=COS(RLONG+E-RLONMP)                                            
      W4=SINLAM*W1+COSLAM*W2*W4                                         
      RLATM=ASIN(W4)                                                    
      DEC=ATAN2(COSLAM*W2*W3,SINLAM-W1*W4)                              
      DEC=-DEC                                                          
      RLONM=ATAN2(W3*W2,(W4*SINLAM-W1)/COSLAM)                          
C  ASSUME DIP NOT NEEDED                                                
      GO TO 1                                                           
      DIP=ATAN(2.*TAN(RLATM))                                           
    1 CONTINUE                                                          
      RETURN                                                            
      ENTRY MTG(RLATG,RLONG,RLATMP,RLONMP,RLATM,RLONM)                  
C     ****     INVERSE TRANSFORMATION                                   
C     ****                                                              
C     ****     PARAMETERS AS ABOVE EXCEPT INPUT --  RLATM,RLONM         
C     ****                                OUTPUT --  RLATG,RLONG        
C     ****                                                              
      SINLAM=SIN(RLATMP)                                                
      COSLAM=COS(RLATMP)                                                
      SINLOM=SIN(RLONMP)                                                
      COSLOM=COS(RLONMP)                                                
      W1=SIN(RLATM)                                                     
      W2=COS(RLATM)                                                     
      W3=SIN(RLONM)                                                     
      W4=COS(RLONM)                                                     
      RLATG=SINLAM*W1-COSLAM*W2*W4                                      
      W3=W3*W2                                                          
      W4=-(RLATG*SINLAM-W1)/COSLAM                                      
      RLATG=ASIN(RLATG)                                                 
      RLONG=ATAN2(SINLOM*W4+COSLOM*W3,COSLOM*W4-SINLOM*W3)              
    2 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE COFLT(XLMBDA,COLAT,CLL,CLT,CTT,CL,CT,CF)               
      PARAMETER (NLMB=73,NTHET=37,NTHM=36,IWHOLE=1)                     
      COMMON/COFPM/ISUM,SCTA(NTHET),FPARM(4)                            
      COMMON/CONDU/SOO(NLMB,NTHM),SOT(NLMB,NTHM),STO(NLMB,NTHM),        
     1  STT(NLMB,NTHM),DCO(NLMB,NTHM),DCT(NLMB,NTHM),DXO(NTHM),         
     2  DAY,UT,ALAMP,PHIP,ZPO,RREF,DLMBDA,DTHETA                        
      KM = 1                                                            
      IF (ISUM.GT.0) KM = 2                                             
      L=INT((XLMBDA-FPARM(1))/DLMBDA+0.5)+1                             
      J=INT((COLAT-FPARM(3))/DTHETA+0.5)+1                              
      IF (J.EQ.1) GO TO 900                                             
      IF (IWHOLE.EQ.1.AND.J.EQ.NTHET) GO TO 900                         
      LP = L+1                                                          
      IF (LP.GT.NLMB) LP = LP - NLMB + 1                                
      LM = L - 1                                                        
      IF (LM.LT.1) LM = LM + NLMB - 1                                   
      DO 10 K=1,KM                                                      
      PA = SCTA(J-1)                                                    
      PB = SCTA(J)                                                      
      PC = SCTA(J+1)                                                    
      FO = (3.*PA + PC)*(PA + 3.*PC)/(8.*PB*(PA + PC))                  
      DEN = 4.*(PA + PC)*(DXO(J) + DXO(J-1))                            
      F1P = (PA + 3.*PC)* DXO(J)/DEN                                    
      F1M = (3.*PA + PC)*DXO(J-1)/DEN                                   
      FT = 2.*PA*PC/(PB*(PA + PC))                                      
      DEN = PA*DXO(J-1) + PC*DXO(J)                                     
      F3P = PC*DXO(J)/DEN                                               
      F3M = PA*DXO(J-1)/DEN                                             
      SXX = .5*(SOO(L,J) + SOO(L,J-1))                                  
      DSXXDXO = (SOO(L,J) - SOO(L,J-1))/DTHETA                          
      DSYXDXT = (STO(LP,J) + STO(LP,J-1) - STO(LM,J) - STO(LM,J-1))/    
     1  (4.*DLMBDA)                                                     
      DSXYDXO = (FO*(SOT(L,J) - SOT(L,J-1)) - .5* (1. - FT)*            
     1  (STO(L,J) - STO(L,J-1)))/DTHETA                                 
      DSYYDXT = FO*(STT(LP,J) + STT(LP,J-1) - STT(LM,J) - STT(LM,J-1))/ 
     1  (4.*DLMBDA)                                                     
      TSC = 2.*(F1P*SOT(L,J) + F1M*SOT(L,J-1)) + F3P*STO(L,J) +         
     1  F3M*STO(L,J-1)                                                  
      SYY = .5*FO*(STT(L,J) + STT(L,J-1))                               
      IF (K.EQ.2) GO TO 5                                               
      CLL = SYY                                                         
      CLT = TSC                                                         
      CTT = SXX                                                         
      CL = DSXYDXO + DSYYDXT                                            
      CT = DSXXDXO + DSYXDXT                                            
      GO TO 10                                                          
    5 CLL = CLL + SYY                                                   
      CLT = CLT - TSC                                                   
      CTT = CTT + SXX                                                   
      CL = CL + DSXYDXO + DSYYDXT                                       
      CT = CT - DSXXDXO - DSYXDXT                                       
   10 J = NTHET - J + 1                                                 
      GO TO 999                                                         
  900 CLL = 1.                                                          
      CLT = 0.                                                          
      CTT = 1.                                                          
      CL = 0.                                                           
      CT = 0.                                                           
  999 CONTINUE                                                          
      CF=0.                                                             
      RETURN                                                            
      END                                                               
      SUBROUTINE JPARPOL(K,FAC,INCLVB,PHI)                              
C  FAC = FIELD ALIGNED CURRENT TIMES (RADIUS OF EARTH)**2 IN AMP AT POLE
C  K = 1 IS FOR NORTH POLE, K = NTHET - 1 IS FOR SOUTH POLE.            
C  INCLVB = 0 LEAVES OUT VXB PART.  INCLVB = 1 INCLUDES VXB PART.       
      PARAMETER (NLMB=73,NTHM=36,NTHET=37)                              
      COMMON/CONDU/SOO(NLMB,NTHM),SOT(NLMB,NTHM),STO(NLMB,NTHM),        
     1  STT(NLMB,NTHM),DCO(NLMB,NTHM),DCT(NLMB,NTHM),DXO(NTHM),         
     2  DAY,UT,ALAMP,PHIP,ZPO,RREF,DLMBDA,DTHETA                        
      DIMENSION PHI(NLMB,NTHET)                                         
      NLM = NLMB - 1                                                    
      FAC = 0.                                                          
      DO 10 L=1,NLM                                                     
      LM = L - 1                                                        
      IF (L.EQ.1) LM = NLM                                              
      DYNCR = 0.                                                        
      IF (INCLVB.GT.0) DYNCR = DCO(L,K)                                 
      HDSOTDX = .25*(SOT(L+1,K) - SOT(LM,K))/DLMBDA                     
      IF (K.EQ.1) GO TO 10                                              
      DYNCR = -DYNCR                                                    
      HDSOTDX = -HDSOTDX                                                
   10 FAC = FAC + DYNCR -(SOO(L,K)/DTHETA-HDSOTDX)*(PHI(L,2) - PHI(L,1))
      FAC = 8.*FAC/(NLM*DTHETA*DTHETA)                                  
      RETURN                                                            
      END                                                               
      SUBROUTINE BNDYC(KBDY,XORY,ALFA,BETA,GAM,GBDY)                    
      ALFA=0.                                                           
      BETA=1.                                                           
      GAM=0.                                                            
      GBDY=0.                                                           
      RETURN                                                            
      END                                                               
      SUBROUTINE EQVCUR(JS,JE)                                          
      PARAMETER (NTHET=37,NLMB=73,IWHOLE=1,NTHPLT=39)                   
      REAL JS,JE,JSC,JEC,JSS,JSA,JES,JEA,JTH,JPH                        
      INTEGER TMPO                                                      
C                                                                       
      COMMON/EQCR/R(49,25),PMOPMMO(25),P(49,25),CONTINT,PSI(NLMB,NTHPLT)
C                                                                       
      COMMON/FAR/F(49),SQ2                                              
      CHARACTER*5 LABX,LABY                                             
      CHARACTER*56 LABL                                                 
      COMMON/CVLX/LABX                                                  
      COMMON/CVLY/LABY                                                  
      COMMON/CVLBL/LABL                                                 
      COMMON/STDCV/MXR,MX,MY,LTYPE, XC, XD, MAGX, MINX,XLABX, YC,       
     1YD, MAGY, MINY,XLABY,XLABL                                        
      COMMON/XYCV/XX(100),YY(100)                                       
      DIMENSION JS(NLMB,NTHET),JE(NLMB,NTHET),DP(49,25),QOST(7,7),C1(49,
     1  49),C2(7,13)                                                    
C                                                                       
      DIMENSION FAK(73)                                                 
      CONTINT = 2.E4                                                    
C                                                                       
      NTH = NTHET - 1                                                   
      QCNST1 = 3.*SQRT(10.)/16.E0                                       
      QCNST2 = 3.*SQRT(70.)/32.E0                                       
      NLM=NLMB                                                          
      RI = 6.481E6                                                      
      SQ2 = SQRT(2.E0)                                                  
      NM = 48                                                           
      MM = 24                                                           
      TMPO = 2*MM + 1                                                   
      NMP = NM + 1                                                      
      MMP = MM + 1                                                      
      NY = NLM                                                          
      DO 50 MP = 1, MMP                                                 
      M = MP - 1                                                        
      XMS = M*M                                                         
      IF (MP.NE.1) PMOPMMO(MP) = SQRT(1. + .5/FLOAT(MP-1))              
      DO 50 NP=MP,NMP                                                   
      XNS = (NP - 1)**2                                                 
      DEN = 4.*XNS - 1                                                  
      IF (DEN.LT.1) DEN = 1.                                            
      R(NP,MP) = SQRT((XNS  - XMS)/DEN)                                 
      C1(NP,MMP-M) = 0.                                                 
   50 C1(NP,MMP+M) = 0.                                                 
      DO 60 MQ=1,13                                                     
      DO 60 NP=1,7                                                      
   60 C2(NP,MQ) = 0.                                                    
      NLMM = NLM - 1                                                    
      JMX = NTH                                                         
      IF (IWHOLE.GT.0) JMX = NTH/2                                      
      DTH = 3.1415926535898/(2.*JMX)                                    
      RDPHOFP = .5*RI/NLMM                                              
      DPH = 2.*3.1415926535898/NLMM                                     
      DO 70 NP=1,7                                                      
      DO 70 MP=1,7                                                      
   70 QOST(NP,MP) = 0.                                                  
      CTM=1.                                                            
      DO 200,J=1,JMX                                                    
      JCON = NTH - J + 1                                                
      TH = (J - .5)*DTH                                                 
      CTP = COS(TH + .5*DTH)                                            
      IF (J.EQ.JMX) CTP = 0.                                            
      DAOFPR = RDPHOFP*(CTM - CTP)                                      
      CTM = CTP                                                         
      CT = COS(TH)                                                      
      ST = SIN(TH)                                                      
      STS = ST*ST                                                       
C  FIRST LET P = PNM/ST                                                 
      P(1,1) = 1./ST                                                    
      DO 80 MP=1,MMP                                                    
      M = MP - 1                                                        
      IF (MP.NE.1) P(MP,MP) = PMOPMMO(MP)*P(M,M)*ST                     
      DP(MP,MP) = M*CT*P(MP,MP)                                         
      PM2 = 0.                                                          
      DO 80 N=MP,NM                                                     
      NP = N + 1                                                        
      P(NP,MP) = (CT*P(N,MP) - R(N,MP)*PM2)/R(NP,MP)                    
      PM2 = P(N,MP)                                                     
   80 DP(NP,MP) = N*CT*P(NP,MP) - (2*N + 1)*R(NP,MP)*PM2                
      XL = 1. + .5*STS*ALOG((1. + CT)/(1. - CT))/CT                     
      QOST(2,2) = 1.2247449*CT                                          
      QOST(3,3) = 1.3693064*CT*ST*XL                                    
      QOST(4,2) =  .2291288*CT*(4. - STS*(3. + 6.*STS))                 
      QOST(4,4) = 1.4790199*CT*STS*(1. + 2.*STS)                        
      X = STS*(2. + 3.*STS*XL)                                          
      QOST(5,3) = QCNST1*CT*ST*(4. - X)                                 
      QOST(5,5) = QCNST2*CT*ST*X                                        
      X = STS*(3. + STS*(4. + 8.*STS))                                  
      QOST(6,2) = .018021728*CT*(56. - STS*(188. - 13.*X))              
      QOST(6,4) =  .5255737*CT*STS*(8. - X)                             
      QOST(6,6) =  .5484353*CT*STS*X                                    
      X = STS*(8. + STS*(10. + 15.*STS*XL))                             
      QOST(7,3) = .057727979*CT*ST*(64. - STS*(168. - 3.*X))            
      QOST(7,5) = .079047290*CT*ST*STS*(80.  - 3.*X)                    
      QOST(7,7) =  .2140611*CT*ST*STS*X                                 
      DO 150 L=1,NLMM                                                   
      PH = (L - 1)*DPH                                                  
      CALL FCMP(PH,MMP)                                                 
      JSC = - JS(L,J)                                                   
      JEC = JE(L,J)                                                     
      IF (IWHOLE.EQ.0) GO TO 90                                         
      JSC = JS(L,JCON)                                                  
      JEC = JE(L,JCON)                                                  
   90 JSS = JS(L,J) - JSC                                               
      JSA = JS(L,J) + JSC                                               
      JES = JE(L,J) + JEC                                               
      JEA = JE(L,J) - JEC                                               
      DO 120 MQ=1,TMPO                                                  
      MQF = MQ - MM + 6                                                 
      M = MQ - MMP                                                      
      MP = 1 + IABS(M)                                                  
      MC = MMP - M                                                      
      IAS = -1                                                          
      DO 120 NP=MP,NMP                                                  
      N = NP - 1                                                        
      IAS = - IAS                                                       
      IF (IAS.EQ.1) GO TO 105                                           
      JPH = JES                                                         
      JTH = JSS                                                         
      GO TO 110                                                         
  105 JPH = JEA                                                         
      JTH = JSA                                                         
      IF (NP.GT.7) GO TO 110                                            
      S = 1.                                                            
      IF (M.LT.0) S= -1.                                                
      C2(NP,MQF) = C2(NP,MQF) + DAOFPR*((M*P(NP,MP) - QOST(NP,MP)*2.*   
     1   CT*S*N)*F(MC)*JTH + S*M*QOST(NP,MP)*F(MQ)*JPH)                 
  110 C1(NP,MQ) = C1(NP,MQ)+ DAOFPR*(DP(NP,MP)*F(MQ)*JPH - M*P(NP,MP)*  
     1  F(MC)*JTH)                                                      
  120 CONTINUE                                                          
  150 CONTINUE                                                          
  200 CONTINUE                                                          
      DO 250 MQ = 1,TMPO                                                
      MQF = MQ - MM + 6                                                 
      M = IABS(MQ - MMP)                                                
      MP = M + 1                                                        
      DO 250 NP=MP,NMP                                                  
      N = NP - 1                                                        
      DEN = N*(N+1)*(1 + ((N*(N+1) - M*M)*DTH*DTH + M*M*DPH*DPH)/12.E0) 
      IF (DEN.EQ.0.) DEN = 1.                                           
      C1(NP,MQ ) = C1(NP,MQ)/DEN                                        
      IF (NP.LE.7) C2(NP,MQF) = C2(NP,MQF)/DEN                          
  250 CONTINUE                                                          
      MY = NTHPLT                                                       
      NHFPLT = (NTHPLT - 1)/2                                           
      PLINC = 180./FLOAT(NTHPLT-3)                                      
      DO 251 J=1,NHFPLT                                                 
      YY(J) = -90. + (J-1)*PLINC                                        
      JCON = NTHPLT +1 - J                                              
  251 YY(JCON) = -YY(J)                                                 
      YY(NHFPLT) = -PLINC/2.                                            
      YY(NHFPLT+1) = 0.                                                 
      YY(NHFPLT+2) = PLINC/2.                                           
      YY(NTHPLT+1) = 91.                                                
C     LABL(1) = 8HEQUIVALE                                              
C     LABL(2) = 8HNT CURRE                                              
C     LABL(3) = 8HNT FUNCT                                              
C     LABL(4) = 8HION (TOR                                              
C     LABL(5) = 8HOIDAL IO                                              
C     LABL(6) = 8HNOSPHERI                                              
C     LABL(7) = 8HC PART)                                               
      LABL = 'EQUIVALENT CURRENT FUNCTION (TOROIDAL IONOSPHERIC PART'   
      CALL PLT(C1,NM,MM,0)                                              
      IF (IWHOLE.EQ.0) GO TO 400                                        
C     LABL(4) = 8HION (ANT                                              
C     LABL(5) = 8HISYMMETR                                              
C     LABL(6) = 8HIC PART)                                              
C     LABL(7) = 8H                                                      
      LABL(29:56) = '(ANTISYMMETRIC PART)'                              
C      CALL PLT(C1,NM,MM,1)                                             
C     LABL(4) = 8HION (SYM                                              
C     LABL(5) = 8HMETRIC T                                              
C     LABL(6) = 8HOROIDAL                                               
C     LABL(7) = 8HPART)                                                 
      LABL(29:56) = '(SYMMETRIC TOROIDAL PART)'                         
C      CALL PLT(C1,NM,MM,-1)                                            
      DO 300 MQF = 1,13                                                 
      MQ = MQF - 6 + MM                                                 
      MP = IABS(MQF - 7) + 1                                            
      DO 300 NP=MP,7                                                    
  300 C1(NP,MQ) = C1(NP,MQ) + C2(NP,MQF)                                
C     LABL(4) = 8HION (TOT                                              
C     LABL(5) = 8HAL)                                                   
C     LABL(6) = 8H                                                      
C     LABL(7) = 8H                                                      
      LABL(29:56) = '(TOTAL)'                                           
      CALL PLT(C1,NM,MM,0)                                              
C     LABL(4) = 8HION (FIE                                              
C     LABL(5) = 8HLD-ALIGN                                              
C     LABL(6) = 8HED PART)                                              
      LABL(29:56) = '(FIELD-ALIGNED PART)'                              
      CALL PLT(C2,6,6,-1)                                               
C     LABL(4) = 8HION (TOT                                              
C     LABL(5) = 8HAL SYMME                                              
C     LABL(6) = 8HTRIC)                                                 
      LABL(29:56) = '(TOTAL SYMMETRIC)'                                 
C      CALL PLT(C1,NM,MM,-1)                                            
C                                                                       
  400 CONTINUE                                                          
      AOB = 6371./6481.                                                 
      AOBN = 1.                                                         
      DO 410 NP=1,NMP                                                   
      AOBN = AOBN*AOB                                                   
      N = NP - 1                                                        
  410 FAK(NP) = AOBN                                                    
      DO 450 MQ = 1,TMPO                                                
      M = IABS(MQ - MMP)                                                
      MP = M + 1                                                        
      DO 450 NP=MP,NMP                                                  
      C1(NP,MQ ) = C1(NP,MQ)*FAK(NP)                                    
  450 CONTINUE                                                          
C     LABL(4) = 8HION (GRO                                              
C     LABL(5) = 8HUND)                                                  
C     LABL(6) = 8H                                                      
      LABL(29:56) = '(GROUND)'                                          
      CALL PLT(C1,NM,MM,0)                                              
      UM = 3.1415926535898*4.E-7                                        
      DO 510 NP=1,NMP                                                   
      N = NP - 1                                                        
  510 FAK(NP) = UM*NP/FLOAT(2*N + 1)                                    
      DO 550 MQ = 1,TMPO                                                
      M = IABS(MQ - MMP)                                                
      MP = M + 1                                                        
      DO 550 NP=MP,NMP                                                  
      C1(NP,MQ ) = C1(NP,MQ)*FAK(NP)                                    
  550 CONTINUE                                                          
      CONTINT = .02                                                     
C     LABL(1) = 8HMAGNETIC                                              
C     LABL(2) = 8H POTENTI                                              
C     LABL(3) = 8HAL (EXTE                                              
C     LABL(4) = 8HRNAL COM                                              
C     LABL(5) = 8HPONENT)                                               
      LABL = 'MAGNETIC POTENTIAL (EXTERNAL COMPONENT)'                  
      CALL PLT(C1,NM,MM,0)                                              
      COA = 5771./6371.                                                 
      COATNPO = 1./COA                                                  
      COAS = COA*COA                                                    
      DO 610 NP=1,NMP                                                   
      COATNPO = COATNPO*COAS                                            
      N = NP - 1                                                        
  610 FAK(NP) = 1. + N*COATNPO/FLOAT(NP)                                
      DO 650 MQ = 1,TMPO                                                
      M = IABS(MQ - MMP)                                                
      MP = M + 1                                                        
      DO 650 NP=MP,NMP                                                  
      C1(NP,MQ ) = C1(NP,MQ)*FAK(NP)                                    
  650 CONTINUE                                                          
C     LABL(3) = 8HAL (TOTA                                              
C     LABL(4) = 8HL)                                                    
C     LABL(5) = 8H                                                      
      LABL(20:56) = '(TOTAL)'                                           
      CALL PLT(C1,NM,MM,0)                                              
      CALL HDPLT(PSI,360.E-9)                                           
C                                                                       
      RETURN                                                            
      END                                                               
      SUBROUTINE FCMP(PH,MMP)                                           
      COMMON/FAR/F(49),SQ2                                              
      MM = MMP - 1                                                      
      F(MMP ) = 1.                                                      
      CP = COS(PH)                                                      
      SP = SIN(PH)                                                      
      F(MMP - 1) = SQ2*CP                                               
      F(MMP + 1) = SQ2*SP                                               
      DO 85 M=2,MM                                                      
      MMO = M - 1                                                       
      F(MMP - M) = CP*F(MMP - MMO) - SP*F(MMP + MMO)                    
   85 F(MMP + M) = CP*F(MMP + MMO) + SP*F(MMP - MMO)                    
      RETURN                                                            
      END                                                               
      SUBROUTINE PLT(C,NM,MM,ISM)                                       
C  ISM=-1 GIVES SYMMETRIC COMPONENT                                     
C  ISM=0  GIVES TOTAL                                                   
C  ISM=1  GIVES ANTISYMMETRIC COMPONENT                                 
      PARAMETER (NLMB=73,NTHPLT=39)                                     
      INTEGER TMPO                                                      
C                                                                       
      COMMON/EQCR/R(49,25),PMOPMMO(25),P(49,25),CONTINT,PSI(NLMB,NTHPLT)
C                                                                       
      COMMON/FAR/F(49),SQ2                                              
      COMMON/XYCV/XX(100),YY(100)                                       
C                                                                       
      DIMENSION C(2401)                                                 
C                                                                       
      RAD = 180./3.1415926535898                                        
      NTH = NTHPLT                                                      
      NLM = NLMB                                                        
      JMX = (NTH + 1)/2                                                 
      TMPO = 2.*MM + 1                                                  
      MMP = MM + 1                                                      
      NMP = NM + 1                                                      
      DO 100 J=1,JMX                                                    
      JCON = 2*JMX - J                                                  
C  COMPUTE IN N. HEMISPHERE (OPPOSITE TO HEMISPHERE OF YY VALUES)       
      CT = COS((90.+YY(J))/RAD)                                         
      ST = SQRT(1.- CT*CT)                                              
      P(1,1) = 1.                                                       
      DO 80 MP=1,MMP                                                    
      M = MP - 1                                                        
      IF (MP.NE.1) P(MP,MP) = PMOPMMO(MP)*P(M,M)*ST                     
      PM2 = 0.                                                          
      DO 80 N=MP,NM                                                     
      NP = N + 1                                                        
      P(NP,MP) = (CT*P(N,MP) - R(N,MP)*PM2)/R(NP,MP)                    
   80 PM2 = P(N,MP)                                                     
      DO 100 L=1,NLM                                                    
      PH = XX(L)*15./RAD                                                
      CALL FCMP(PH,MMP)                                                 
      PSI(L,J ) = 0.                                                    
      PSI(L,JCON) = 0.                                                  
      DO 90 MQ = 1,TMPO                                                 
      MP = 1 + IABS(MQ - MMP)                                           
      IAS = -1                                                          
      DO 90 NP=MP,NMP                                                   
      INDX = (MQ-1)*NMP + NP                                            
      N = NP - 1                                                        
      IAS = - IAS                                                       
      IF (ISM.LT.0) GO TO 85                                            
      IF (IAS.LT.0) PSI(L,JCON) = PSI(L,JCON) + C(INDX)*P(NP,MP)*F(MQ)  
   85 CONTINUE                                                          
      IF (ISM.GT.0) GO TO 90                                            
      IF (IAS.GT.0) PSI(L,J) = PSI(L,J) + C(INDX)*P(NP,MP)*F(MQ)        
   90 CONTINUE                                                          
      IF (ISM) 91,92,93                                                 
   91 PSI(L,JCON) = PSI(L,J)                                            
      GO TO 100                                                         
   92 X = PSI(L,J)                                                      
      PSI(L,J) = PSI(L,J) - PSI(L,JCON)                                 
      PSI(L,JCON) = PSI(L,JCON) +X                                      
      GO TO 100                                                         
   93 PSI(L,J) = -PSI(L,JCON)                                           
  100 CONTINUE                                                          
C                                                                       
      CALL STDCON(PSI,CONTINT,0)                                        
C                                                                       
      RETURN                                                            
      END                                                               
      SUBROUTINE VIPLT(VAR,FULSCL)                                      
      PARAMETER (BEQ=2.82E-5)                                           
C   VAR = ELECTRIC POTENTIAL AT 300 KM.                                 
C   FULSCL = DRIFT VALUE (M/S) CORRESPONDING TO 90 DEGREES LATITUDE.    
C   XX,YY= X AND Y AXIS VALUES THAT CORRESPOND WITH MX AND MY POINTS    
C     IN VAR.                                                           
C                                                                       
      CHARACTER*5 LABX,LABY,LINFMT*6                                    
      CHARACTER*56 LABL                                                 
      COMMON/CVLX/LABX                                                  
      COMMON/CVLY/LABY                                                  
      COMMON/CVLBL/LABL                                                 
      COMMON/STDCV/MXR,MX,MY,LTYPE, XC, XD, MAGX, MINX,XLABX, YC,       
     1YD, MAGY, MINY,XLABY,XLABL                                        
      COMMON/XYCV/XX(100), YY(100)                                      
      DIMENSION VAR(MXR, MY),ZZ(100),YMID(100)                          
C***********************************************************************
      DIMENSION TSHIFT(3),VN(3,48),VE(3,48),PLTLAT(3)                   
     1  ,NN(3,48),NE(3,48)                                              
      CHARACTER*1 ISYMB(3),JSYMB                                        
C***********************************************************************
      RAD = 180./3.1415926535898                                        
      LINFMT= '(F5.0)'                                                  
      SCL = YD/FULSCL                                                   
      RI = 6.671E6                                                      
      DO 33 M=2,MY/2                                                    
   33 YMID(M) = (YY(M-1) + YY(M+1))/2.E0                                
      PRINT 35,(YMID(M),M=2,MY/2)                                       
C  35 FORMAT (1H ,16F5.1)                                               
   35 FORMAT (1X,16F5.1)                                                
      CALL SET (0.1, 0.95, 0.2, 0.8, XC, XD, 0, YD, LTYPE)              
      CXA = (XD-XC)/(.85*1024.)                                         
      CXB = (.95*XC-.1*XD)/.85                                          
      CYA = YD/(.60*1024.)                                              
      CYB = -.2*YD/.60                                                  
      LABL(1:36) = 'UPWARD/POLEWARD DRIFT   (1 DEGREE = '               
      X = 1./SCL                                                        
C     ENCODE (8,46,LABL(5)) X                                           
C  46 FORMAT (4HE = ,F3.1,1H )                                          
      WRITE (LABL(37:39),46)                                            
   46 FORMAT (F3.1)                                                     
      LABL(40:56) = ' M/S)            '                                 
C     CALL PWRITY (105,860, LABL, 56, 2, 0, -1)                         
      CALL PWRITY (105.*CXA+CXB,860.*CYA+CYB, LABL, 56, 2, 0, -1)       
      CALL LABMOD (LINFMT, LINFMT, 5, 5, 1, 1, 0, 0, 0)                 
      CALL PERIML (MAGX, MINX, 9, 2)                                    
C     CALL PWRY (15, 500,LABY, 8, 2.5, 1,1)                             
C     CALL PWRY (525, 140, LABX, 8, 2.5, 0, 1)                          
      CALL PWRITY (15.*CXA+CXB,.5*YD,LABY,5,2.5,90,0)                   
      CALL PWRITY (.5*(XC+XD),140.*CYA+CYB,LABX,5,2.5,0,0)              
      BEQ3 = BEQ*(6.481E6/RI)**3                                        
      DO 200 M=2,MY/2                                                   
      ST = COS(YMID(M)/RAD)                                             
      B = BEQ3*SQRT(4. - 3.*ST*ST)                                      
      RST = RI*ST                                                       
      DO 180 N=2,MX-1                                                   
      TDLON = (XX(N+1) - XX(N-1))*15./RAD                               
  180 ZZ(N) = (VAR(N-1,M) - VAR(N+1,M))*SCL/(TDLON*RST*B) + YMID(M)     
  200 CALL CURVE(XX(2),ZZ(2),MX-2)                                      
      CALL FRAME                                                        
      LABL(1:24) = 'EASTWARD DRIFT'                                     
C     LABL(2) = 8H DRIFT                                                
C     LABL(3) = 8H                                                      
      CALL PWRITY (105,860, LABL, 56, 2, 0, -1)                         
      CALL LABMOD (LINFMT, LINFMT, 5, 5, 1, 1, 0, 0, 0)                 
      CALL PERIML (MAGX, MINX, 9, 2)                                    
      CALL PWRY (15, 500,LABY, 8, 2.5, 1,1)                             
      CALL PWRY (525, 140, LABX, 8, 2.5, 0, 1)                          
      DO 300 M=2,MY/2                                                   
      TDLAT = (YY(M+1) - YY(M-1))/RAD                                   
      CT = SIN(YMID(M)/RAD)                                             
C     BDOWN = BEQ3*2.*CT/SQRT(1. + 3.*CT*CT)                            
C ABOVE IN ERROR. CORRECTED HERE 86 JUNE 19.                            
      BDOWN = BEQ3*4.*CT/SQRT(1. + 3.*CT*CT)                            
      DO 280 N=2,MX-1                                                   
C 280 ZZ(N) = (VAR(N,M-1) - VAR(N,M+1))*SCL/(RI*BDOWN*TDLAT) + YMID(M)  
C ABOVE IN ERROR. CORRECTED HERE 86 JUNE 19.                            
  280 ZZ(N) = (VAR(N,M+1) - VAR(N,M-1))*SCL/(RI*BDOWN*TDLAT) + YMID(M)  
  300 CALL CURVE(XX(2),ZZ(2),MX-2)                                      
      CALL FRAME                                                        
C READ OBSERVED VALUES HERE                                             
      NSTA = 3                                                          
      DO 1120 M=1,NSTA                                                  
      READ (5,1110) CLAT,XLON,ISYMB(M),PLTLAT(M)                        
 1110 FORMAT(2F10.0,7X,A1,F5.2)                                         
      WRITE (6,1111) CLAT,XLON,ISYMB(M),PLTLAT(M)                       
C1111 FORMAT(1H ,2F10.0,4X,A1,F5.2)                                     
 1111 FORMAT(1X,2F10.0,4X,A1,F5.2)                                      
      TSHIFT(M) = (XLON - 69.)/15.                                      
      READ (5,1112) ((VN(M,I),NN(M,I),I=1,48),(VE(M,I),NE(M,I),I=1,48)) 
 1112 FORMAT (8(F4.0,3X,I3))                                            
      WRITE (6,1113) ((VN(M,I),I=1,48),(VE(M,I),I=1,48))                
C1113 FORMAT (1H ,24F5.1)                                               
 1113 FORMAT (1X,24F5.1)                                                
 1120 CONTINUE                                                          
      CALL SET (0.1, 0.95, 0.2, 0.8, XC, XD, 0, YD, LTYPE)              
      CXA = (XD-XC)/(.85*1024.)                                         
      CXB = (.95*XC-.1*XD)/.85                                          
      CYA = YD/(.60*1024.)                                              
      CYB = -.2*YD/.60                                                  
C     LABL(1) = 8HUPWARD/P                                              
C     LABL(2) = 8HOLEWARD                                               
C     LABL(3) = 8HDRIFT                                                 
C     CALL PWRITY (105,860, LABL, 56, 2, 0, -1)                         
C     CALL LABMOD (LINFMT, LINFMT, 5, 5, 1, 1, 0, 0, 0)                 
C     CALL PERIML (MAGX, MINX, 9, 2)                                    
C     CALL PWRY (15, 500,LABY, 8, 2.5, 1,1)                             
C     CALL PWRY (525, 140, LABX, 8, 2.5, 0, 1)                          
      DO 1200 M=1,NSTA                                                  
      NPTS = 0                                                          
      DO 1180 N=1,48                                                    
      IF (NN(M,N).EQ.0) GO TO 1180                                      
      NPTS = NPTS + 1                                                   
      X = .5*N + TSHIFT(M)                                              
      IF (X.GT.24.) X = X - 24.                                         
      IF (X.LT.0.) X = X + 24.                                          
      YMID(NPTS) = X                                                    
      ZZ(NPTS) = VN(M,N)*SCL + PLTLAT(M)                                
 1180 CONTINUE                                                          
      JSYMB = ISYMB(M)                                                  
      CALL LINE(0.,PLTLAT(M),24.,PLTLAT(M))                             
C1200 CALL POINTS(YMID,ZZ,NPTS,JSYMB,0)                                 
 1200 CALL POINTS(YMID,ZZ,NPTS,ICHAR(JSYMB),0)                          
      CALL FRAME                                                        
C     LABL(1) = 8HEASTWARD                                              
C     LABL(2) = 8H DRIFT                                                
C     LABL(3) = 8H                                                      
C     CALL PWRITY (105,860, LABL, 56, 2, 0, -1)                         
C     CALL LABMOD (LINFMT, LINFMT, 5, 5, 1, 1, 0, 0, 0)                 
C     CALL PERIML (MAGX, MINX, 9, 2)                                    
C     CALL PWRY (15, 500,LABY, 8, 2.5, 1,1)                             
C     CALL PWRY (525, 140, LABX, 8, 2.5, 0, 1)                          
      DO 1400 M=1,NSTA                                                  
      NPTS = 0                                                          
      DO 1380 N=1,48                                                    
      IF (NE(M,N).EQ.0) GO TO 1380                                      
      NPTS = NPTS + 1                                                   
      X = .5*N + TSHIFT(M)                                              
      IF (X.GT.24.) X = X - 24.                                         
      IF (X.LT.0.) X = X + 24.                                          
      YMID(NPTS) = X                                                    
      ZZ(NPTS) = VE(M,N)*SCL + PLTLAT(M)                                
 1380 CONTINUE                                                          
      CALL LINE(0.,PLTLAT(M),24.,PLTLAT(M))                             
      JSYMB = ISYMB(M)                                                  
C1400 CALL POINTS(YMID,ZZ,NPTS,JSYMB,0)                                 
 1400 CALL POINTS(YMID,ZZ,NPTS,ICHAR(JSYMB),0)                          
      CALL FRAME                                                        
      RETURN                                                            
      END                                                               
      SUBROUTINE HDPLT(VAR,FULSCL)                                      
C   VAR(MXR, MY) = MAGNETIC POTENTIAL                                   
C   FULSCL = VAR VALUE CORRESPONDING TO DISTANCE BETWEEN TOP AND BOTTOM.
      CHARACTER*8 LABX,LABY,LINFMT*6                                    
      CHARACTER*56 LABL                                                 
      COMMON/CVLX/LABX                                                  
      COMMON/CVLY/LABY                                                  
      COMMON/CVLBL/LABL                                                 
      COMMON/STDCV/MXR,MX,MY,LTYPE, XC, XD, MAGX, MINX,XLABX, YC,       
     1YD, MAGY, MINY,XLABY,XLABL                                        
      COMMON/XYCV/XX(100), YY(100)                                      
      DIMENSION VAR(MXR, MY),ZZ(100),YMID(100)                          
C***********************************************************************
      DIMENSION TSHIFT(5),HH(5,24),DD(5,24),PLTLAT(5)                   
      CHARACTER*1 ISYMB(5),JSYMB                                        
C***********************************************************************
      LINFMT= '(F5.0)'                                                  
      RAD = 180./3.1415926535898                                        
      SCL = ABS(YC - YD)/FULSCL                                         
      RE = 6.371E6                                                      
      DO 33 M=2,MY-1                                                    
   33 YMID(M) = (YY(M-1) + YY(M+1))/2.E0                                
      PRINT 35,(YMID(M),M=2,MY-1)                                       
C  35 FORMAT (1H ,18F5.1)                                               
   35 FORMAT (1X,18F5.1)                                                
      CALL SET (0.1, 0.95, 0.2, 0.8, XC, XD, YC, YD, LTYPE)             
      CXA = (XD-XC)/(.85*1024.)                                         
      CXB = (.95*XC-.1*XD)/.85                                          
      CYA = (YD-YC)/(.60*1024.)                                         
      CYB = (.8*YC-.2*YD)/.60                                           
      LABL(1:28) = 'H VARIATION     1 DEGREE = '                        
C     LABL(2) = 8HION                                                   
C     LABL(3) = 8H(1 DEGRE                                              
      X = 1.E9/SCL                                                      
C     ENCODE (8,46,LABL(4)) X                                           
C  46 FORMAT (4HE = ,F3.1,1H )                                          
      WRITE (LABL(29:31),46)                                            
   46 FORMAT (F3.1)                                                     
      LABL(32:56) = ' NT'                                               
C     LABL(5) = 8HNT)                                                   
C     LABL(6) = 8H                                                      
C     LABL(7) = 8H                                                      
C     CALL PWRITY (105,860, LABL, 56, 2, 0, -1)                         
      CALL PWRITY (105.*CXA+CXB,860.*CYA+CYB, LABL, 56, 2, 0, -1)       
      CALL LABMOD (LINFMT, LINFMT, 5, 5, 1, 1, 0, 0, 0)                 
      CALL PERIML (MAGX, MINX, MAGY, MINY)                              
C     CALL PWRY (15, 500,LABY, 8, 2.5, 1,1)                             
C     CALL PWRY (525, 140, LABX, 8, 2.5, 0, 1)                          
      CALL PWRITY (15.*CXA+CXB,.5*YD,LABY,5,2.5,90,0)                   
      CALL PWRITY (.5*(XC+XD),140.*CYA+CYB,LABX,5,2.5,0,0)              
      MID = MY/2 + 1                                                    
      DO 200 M=2,MY-1                                                   
      IF (M.EQ.MID-1 .OR. M.EQ.MID+1) GO TO 200                         
      TDLAT = (YY(M+1) - YY(M-1))/RAD                                   
      DO 180 N=2,MX-1                                                   
  180 ZZ(N) = (VAR(N,M-1) - VAR(N,M+1))*SCL/(RE*TDLAT) + YMID(M)        
      CALL CURVE(XX(2),ZZ(2),MX-2)                                      
  200 CONTINUE                                                          
      CALL FRAME                                                        
C     LABL(1) = 8HD VARIAT                                              
      LABL(1:1) = 'D'                                                   
      CALL PWRITY (105,860, LABL, 56, 2, 0, -1)                         
      CALL LABMOD (LINFMT, LINFMT, 5, 5, 1, 1, 0, 0, 0)                 
      CALL PERIML (MAGX, MINX, MAGY, MINY)                              
      CALL PWRY (15, 500,LABY, 8, 2.5, 1,1)                             
      CALL PWRY (525, 140, LABX, 8, 2.5, 0, 1)                          
      DO 300 M=2,MY-1                                                   
      RST = RE*COS(YMID(M)/RAD)                                         
      DO 280 N=2,MX-1                                                   
      TDLON = (XX(N+1) - XX(N-1))*15./RAD                               
  280 ZZ(N) = (VAR(N-1,M) - VAR(N+1,M))*SCL/(RST*TDLON) + YMID(M)       
  300 CALL CURVE(XX(2),ZZ(2),MX-2)                                      
      CALL FRAME                                                        
C READ OBSERVED VALUES HERE                                             
      NSTA = 5                                                          
      DO 1120 M=1,NSTA                                                  
      READ (5,1110) XLAT,XLON,ISYMB(M),PLTLAT(M),NZEROH,SLOPE           
 1110 FORMAT(2F10.0,7X,A1,F5.2,I7,F10.0)                                
      WRITE (6,1111) XLAT,XLON,ISYMB(M),PLTLAT(M)                       
C1111 FORMAT(1H ,2F10.0,4X,A1,F5.2)                                     
 1111 FORMAT(1X,2F10.0,4X,A1,F5.2)                                      
      TSHIFT(M) = (XLON - 69.)/15.                                      
      READ (5,1112) ((HH(M,I),I=1,24),(DD(M,I),I=1,24))                 
 1112 FORMAT (12F5.0)                                                   
      DO 1115 I=1,24                                                    
 1115 HH(M,I) = HH(M,I) + (SLOPE/24.)*(I-NZEROH)                        
      WRITE (6,1113) ((HH(M,I),I=1,24),(DD(M,I),I=1,24))                
C1113 FORMAT (1H ,24F5.1)                                               
 1113 FORMAT (1X,24F5.1)                                                
 1120 CONTINUE                                                          
      CALL SET (0.1, 0.95, 0.2, 0.8, XC, XD,YC, YD, LTYPE)              
      CXA = (XD-XC)/(.85*1024.)                                         
      CXB = (.95*XC-.1*XD)/.85                                          
      CYA = (YD-YC)/(.60*1024.)                                         
      CYB = (.8*YC-.2*YD)/.60                                           
C     LABL(1) = 8HUPWARD/P                                              
C     LABL(2) = 8HOLEWARD                                               
C     LABL(3) = 8HDRIFT                                                 
C     CALL PWRITY (105,860, LABL, 56, 2, 0, -1)                         
C     CALL LABMOD (LINFMT, LINFMT, 5, 5, 1, 1, 0, 0, 0)                 
C     CALL PERIML (MAGX, MINX, 9, 2)                                    
C     CALL PWRY (15, 500,LABY, 8, 2.5, 1,1)                             
C     CALL PWRY (525, 140, LABX, 8, 2.5, 0, 1)                          
      DO 1200 M=1,NSTA                                                  
      NPTS = 0                                                          
      DO 1180 N=1,24                                                    
      NPTS = NPTS + 1                                                   
      X = N-.5 + TSHIFT(M)                                              
      IF (X.GT.24.) X = X - 24.                                         
      IF (X.LT.0.) X = X + 24.                                          
      YMID(NPTS) = X                                                    
      ZZ(NPTS) = HH(M,N)*SCL*1.E-9 + PLTLAT(M)                          
 1180 CONTINUE                                                          
      JSYMB = ISYMB(M)                                                  
      CALL LINE(0.,PLTLAT(M),24.,PLTLAT(M))                             
C1200 CALL POINTS(YMID,ZZ,NPTS,JSYMB,0)                                 
 1200 CALL POINTS(YMID,ZZ,NPTS,ICHAR(JSYMB),0)                          
      CALL FRAME                                                        
C     LABL(1) = 8HEASTWARD                                              
C     LABL(2) = 8H DRIFT                                                
C     LABL(3) = 8H                                                      
C     CALL PWRITY (105,860, LABL, 56, 2, 0, -1)                         
C     CALL LABMOD (LINFMT, LINFMT, 5, 5, 1, 1, 0, 0, 0)                 
C     CALL PERIML (MAGX, MINX, 9, 2)                                    
C     CALL PWRY (15, 500,LABY, 8, 2.5, 1,1)                             
C     CALL PWRY (525, 140, LABX, 8, 2.5, 0, 1)                          
      DO 1400 M=1,NSTA                                                  
      NPTS = 0                                                          
      DO 1380 N=1,24                                                    
      NPTS = NPTS + 1                                                   
      X = N-.5 + TSHIFT(M)                                              
      IF (X.GT.24.) X = X - 24.                                         
      IF (X.LT.0.) X = X + 24.                                          
      YMID(NPTS) = X                                                    
      ZZ(NPTS) = DD(M,N)*SCL*1.E-9 + PLTLAT(M)                          
 1380 CONTINUE                                                          
      CALL LINE(0.,PLTLAT(M),24.,PLTLAT(M))                             
      JSYMB = ISYMB(M)                                                  
C1400 CALL POINTS(YMID,ZZ,NPTS,JSYMB,0)                                 
 1400 CALL POINTS(YMID,ZZ,NPTS,ICHAR(JSYMB),0)                          
      CALL FRAME                                                        
      RETURN                                                            
      END                                                               
      SUBROUTINE EFIELD(XMLAT,XMLON,DAYNO,UT,ISEASAV,IUTAV,POT,VU,VE)   
C  GIVES QUIET-DAY IONOSPHERIC ELECTROSTATIC PSEUDO-POTENTIAL AND E X B 
C  DRIFTS AT 300 KM FOR SOLAR MINIMUM CONDITIONS.  SEE RICHMOND ET AL. (
C  1980, P.     ) FOR DEFINITIONS OF MAGNETIC COORDINATES, PSEUDO-POTENT
C  AND DRIFT COMPONENTS.                                                
C                                                                       
C  INPUT PARAMETERS -                                                   
C    XMLAT, XMLON ARE MAGNETIC LATITUDE AND EAST LONGITUDE IN DEGREES.  
C    DAYNO IS DAY NUMBER OF THE YEAR FROM 1. TO 365.24, WITH 1. BEING JA
C    UT IS UNIVERSAL TIME IN HOURS.                                     
C    ISEASAV IS 0 IF NO SEASONAL AVERAGING IS DESIRED.                  
C            IS 1 FOR AVERAGE OVER NOV. - FEB.                          
C            IS 2 FOR AVERAGE OVER MAY -AUG.                            
C            IS 3 FOR AVERAGE OVER MAR., APR., SEPT., OCT.              
C            IS 4 FOR AVERAGE OVER ENTIRE YEAR.                         
C      IF ISEASAV.NE.0, DAYNO IS IGNORED.                               
C    IUTAV IS 0 IF NO UT AVERAGING IS DESIRED.                          
C          IS 1 FOR AVERAGE OVER ALL UT AT THE FIXED LOCAL TIME GIVEN BY
C            UT + (XMLON - 69.)/15.                                     
C                                                                       
C  OUTPUT PARAMETERS -                                                  
C    POT IS THE ELECTROSTATIC PSEUDO-POTENTIAL IN VOLTS.                
C    VU IS THE DRIFT VELOCITY COMPONENT PERPENDICULAR TO THE GEOMAGNETIC
C      IN THE UPWARD/POLEWARD DIRECTION IN THE MAGNETIC MERIDIAN PLANE, 
C    VE IS THE DRIFT VELOCITY COMPONENT IN THE MAGNETIC EASTWARD DIRECTI
C      M/S.                                                             
C  THE OUTPUT VALUES ARE GEOPHYSICALLY MEANINGFUL ONLY FOR LATITUDES BET
C    ABOUT -65 AND +65 DEGREES.  IF ISEASAV OR IUTAV IS OUT OF RANGE POT
C    AND VE ARE SET TO -1/0.                                            
      DIMENSION KF(128),LF(128),MF(128),NF(128),JF(128),Q(5),RS(16),    
     1   FUT(5),RR(16),RSRS(16),P(16),VP(16),PA(5,9),VA(5,9),FS(3),     
     2  FT(3,3),SML(4),CML(4),              A(128),PB(9),VB(9)          
      DATA IFIRST/1/,ISVP/5/,A/                                         
     * -70.,-183.,  31.,-112.,  19., -39.,  -2.,   2., -33.,   2.,   2.,
     *-111.,  46.,  -4.,  -5.,   7.,   9., -17.,   2.,   9., -10.,   2.,
     *  -9.,  22., 145., -57., -42.,  -6.,   6.,  -5.,  -2.,  20.,  16.,
     *  16., -77., -18.,  13.,  -8.,  16., -52., -10.,   7.,   2.,  11.,
     * -28.,   2., -85., -82.,   3.,-281., -71., -25., -57., -50.,  21.,
     * -10.,  10., -81.,  24.,   7.,   5.,  30.,  32.,   5.,  -5.,  11.,
     * -31.,   8.,  10.,  20., -15., -42.,  32.,   7., -19.,   7.,  34.,
     * -11., -15.,  26.,  21.,   1.,  22.,  12.,  -2., 275., 777.,-318.,
     *-320.,-208.,  47., 429.,-523.,   8., -35.,-224.,-450., -66.,  -7.,
     *  -8.,-231.,  55.,   6., -28., -51., -81.,  48.,   9.,   2., -10.,
     *  54.,  16., 112.,  69., -33., 120., -47.,   5., -19., -17., -23.,
     * -40., -22., -21.,  -7., -30.,  15.,   3./                        
      IF (IFIRST.EQ.0) GO TO 510                                        
      IFIRST = 0                                                        
C  (THROUGH STATEMENT 510) SET UP CONSTANT PARAMETERS IN FIRST CALL TO  
C    SUBROUTINE.                                                        
C  SET UP VALUES OF XMLATP, TUP, TLP, AND DAYNOP WHICH ARE NOT EQUAL TO 
C    UT, MAGNETIC LOCAL TIME, AND DAYNO, RESPECTIVELY.                  
      XMLATP = -361.                                                    
      IF (XMLATP.EQ.XMLAT) XMLATP = 0.                                  
      TUP = -25.                                                        
      IF (TUP.EQ.UT) TUP = 0.                                           
      TL = UT + (XMLON - 69.)/15.                                       
      TLP = -25.                                                        
      IF (TLP.EQ.TL) TLP = 0.                                           
      DAYNOP = -366.                                                    
      IF (DAYNOP.EQ.DAYNO) DAYNOP = 0.                                  
C  (THROUGH STATEMENT 100) SELECT ONLY THOSE TERMS IN SERIES FOR WHICH  
C    COEFFICIENTS A ARE DEFINED TO BE NON-ZERO.                         
      I = 0                                                             
      DO 100 KP=1,3                                                     
      LST = 4 - KP                                                      
      LEND = 2 + KP                                                     
      DO 100 LP = LST,LEND                                              
      L = LP - 3                                                        
      IF (MOD(KP+LP,2) .NE.0) GO TO 100                                 
      DO 90 NP=2,8                                                      
      N = NP - 1                                                        
      IF (KP.NE.3.AND.N.GT.6) GO TO 90                                  
      IF (IABS(L).EQ.2.AND.N.GT.5) GO TO 90                             
      DO 80 MP = 1,9                                                    
      M = MP - 5                                                        
      IF (IABS(M).GT.N) GO TO 80                                        
      IF (IABS(M).GT.3.AND.IABS(L).EQ.2) GO TO 80                       
      IF (MOD(N-M,2).NE.0) GO TO 80                                     
      I = I + 1                                                         
      KF(I) = KP                                                        
      LF (I)= LP                                                        
      NF(I) = NP                                                        
      MF(I) = MP                                                        
   80 CONTINUE                                                          
   90 CONTINUE                                                          
  100 CONTINUE                                                          
      IMAX = I                                                          
      FT(1,1) = .75*SQRT(6.E0)/3.1415926535898                          
      FT(1,2) = 2.E0*FT(1,1)                                            
      FT(1,3) = 1.E0                                                    
      FT(2,1) = FT(1,1)                                                 
      FT(2,2) = -FT(1,2)                                                
      FT(2,3) = 1.E0                                                    
      FT(3,1) = FT(2,2)                                                 
      FT(3,2) = 0.                                                      
      FT (3,3) = 1.E0                                                   
      HRANG= 3.1415926535898/12.                                        
      DANG = 3.1415926535898/182.62                                     
      SQ2 = SQRT(2.E0)                                                  
      RAD = 180./3.1415926535898                                        
      RB = -6.671E6*5.2E-5                                              
C RB IS -(EARTH RADIUS + 3.E5 M) TIMES DIPOLE MAGNETIC FIELD AT POLE AT 
      DO 400 I=1,IMAX                                                   
      MM = IABS(MF(I) - 5)                                              
C  JF GIVES APPROPRIATE INDEX OF LEGENDRE POLYNOMIALS AS ORDERED BETWEEN
C    STATEMENTS 530 AND 595.                                            
  400 JF(I) = (2*(NF(I) + 7*MM) - (MM - 1)**2 + 4)/4                    
C  (THROUGH STATEMENT 500) COMPUTE RR (DEFINED AS R(N,M)*R(N-1,M)) AND R
C    (DEFINED AS R(N-1,M)**2 + R(N-2,M)**2) NEEDED FOR LEGENDRE POLYNOMI
C    GENERATING RECURSION RELATIONS, WHERE R(N,M) IS DEFINED AS         
C    SQRT(N**2 - M**2)/SQRT(4*N**2 - 1).  ORDERING IS SAME AS FOR P(N,M)
      J = 0                                                             
      DO 500 MP=1,5                                                     
      M = MP - 1                                                        
      XM = M                                                            
      IF (M.NE.0) Q(MP) = SQRT((2.*XM + 1.)/(2.*XM))                    
      DO 500 NP = MP,8,2                                                
      N = NP - 1                                                        
      XNS = N*N                                                         
      XNMS = (N-1)**2                                                   
      J = J + 1                                                         
      RS(J) = (XNS - XM*XM)/(4.*XNS - 1.)                               
      RSM    = (XNMS - XM*XM)/(4.*XNMS - 1.)                            
      RR(J) = SQRT(RS(J)*RSM)                                           
      IF (NP.NE.MP) RSRS(J) = RSM    + RS(J - 1)                        
  500 CONTINUE                                                          
  510 CONTINUE                                                          
      IF (IUTAV.LT.0.OR.IUTAV.GT.1) GO TO 988                           
      IF (IABS(ISEASAV-2).GT.2) GO TO 988                               
C  (THROUGH STATEMENT 530) IF SEASONAL HARMONIC AMPLITUDES ARE DIFFERENT
C    PREVIOUS CALL TO THIS SUBROUTINE, RECOMPUTE THEM.                  
      ICPT = 0                                                          
      IF (ISEASAV.EQ.ISVP) GO TO 524                                    
      ICPT = 1                                                          
      ISVP = ISEASAV                                                    
C  (THROUGH STATEMENT 530) FS(1), FS(2), FS(3) ARE FACTORS FOR AMPLITUDE
C    SEMIANNUAL, ANNUAL, YEARLY AVERAGE COMPONENTS, RESPECTIVELY.       
C  IF ISEASAV = 0, COMPUTE FS FOR GIVEN DAY OF THE YEAR.                
      IF (ISEASAV.EQ.0)GO TO 525                                        
C  IF ISEASAV = 4, USE ONLY YEARLY AVERAGE COMPONENT.                   
      IF (ISEASAV.EQ.4) GO TO 527                                       
C  IF ISEASAV = 1 - 3, COMPUTE FS FOR APPROPRIATE SEASONAL AVERAGE.     
      DO 521 K=1,3                                                      
  521 FS(K) = FT(ISEASAV,K)                                             
      GO TO 530                                                         
  524 CONTINUE                                                          
      IF (ISEASAV.NE.0) GO TO 530                                       
      IF (DAYNO.EQ.DAYNOP) GO TO 530                                    
  525 DAYNOP = DAYNO                                                    
      ICPT = 1                                                          
      ANG = (DAYNO + 9.)*DANG                                           
      FS(1) = SQ2*COS(2.*ANG)                                           
      FS(2) = SQ2*COS(ANG)                                              
      FS(3) = 1.                                                        
      GO TO 530                                                         
  527 FS(1) = 0.                                                        
      FS(2) = 0.                                                        
      FS(3) = 1.                                                        
  530 CONTINUE                                                          
C  IF MAGNETIC LATITUDE IS SAME AS IN PREVIOUS CALL TO THIS SUBROUTINE, 
C    STATEMENT 596.                                                     
      IF (XMLAT.EQ.XMLATP) GO TO 596                                    
      ICPT = 1                                                          
      XMLATP = XMLAT                                                    
      CT = SIN(XMLAT/RAD)                                               
      CTS = CT*CT                                                       
      ST = SQRT(1. - CTS)                                               
      RBT = RB*SQRT(.25 + .75*CTS)                                      
C  (THROUGH STATEMENT 595) CALCULATE LEGENDRE POLYNOMIALS P(N,M) AS WELL
C    DEFINED AS (DP(N,M)/D(COLATITUDE))/(RB*CT).                        
      J = 0                                                             
      DO 595 MP = 1,5                                                   
C  FOR MP=1, P=P(N,M).  FOR MP.GT.1, P=P(N,M)/ST.  THIS DIFFERENCE IS SO
C    PROGRAM NEVER DIVIDES BY ST.                                       
      J = J + 1                                                         
      XM = MP - 1                                                       
      MPP = MP + 2                                                      
      IF (MP.GT.1) GO TO 550                                            
      X = 1.                                                            
      P(1) = 1.E0                                                       
      XZ = 2.*ST/RB                                                     
      VP(1) = 0.                                                        
      GO TO 560                                                         
  550 P(J) = Q(MP)*X                                                    
      X = P(J)*ST                                                       
      VP(J) = XM*P(J)/RB                                                
      XZ = 2.*ST*ST/RB                                                  
  560 DO 590 NP = MPP,8,2                                               
      J = J + 1                                                         
      Z = 0.                                                            
      ZP = 0.                                                           
      IF (NP.EQ.MPP) GO TO 570                                          
      Z = RR(J-1)*P(J-2)                                                
      ZP = RR(J-1) *VP(J-2)                                             
  570 P(J) = ((CTS - RSRS(J))*P(J-1) - Z)/RR(J)                         
  590 VP(J) = ((CTS - RSRS(J))*VP(J-1) - ZP - XZ*P(J-1))/RR(J)          
  595 CONTINUE                                                          
  596 CONTINUE                                                          
C  (THROUGH STATEMENT 600) CALCULATE ARRAYS OF FOURIER COEFFICIENTS, WIT
C    INDICATING HARMONIC OF UT AND MPF INDICATING HARMONIC OF MAGNETIC L
C    TIME.  IF ARRAYS ARE SAME AS IN PREVIOUS CALL TO THIS SUBROUTINE, S
C    STATEMENT 601.                                                     
      IF (ICPT.EQ.0) GO TO 601                                          
      DO 597 MPF = 1,9                                                  
      DO 597 LP=1,5                                                     
      PA(LP,MPF) = 0.                                                   
  597 VA(LP,MPF) = 0.                                                   
      DO 600 I=1,IMAX                                                   
      LP = LF(I)                                                        
      MPF = MF(I)                                                       
      J = JF(I)                                                         
      X = A(I)*FS(KF(I))                                                
      PA(LP,MPF) = PA(LP,MPF) + X*P(J)                                  
      VA(LP,MPF) = VA(LP,MPF) + X*VP(J)                                 
  600 CONTINUE                                                          
  601 CONTINUE                                                          
C  (THROUGH STATEMENT 607) CALCULATE FOURIER COEFFICIENTS PB AND VB AT G
C    FOR HARMONICS OF TL.                                               
C  IF UT IS SAME AS IN PREVIOUS CALL TO THIS SUBROUTINE, SKIP TO STATEME
      IF (UT.EQ.TUP) GO TO 603                                          
      TUP = UT                                                          
      ICPT = 1                                                          
      TUA = UT*HRANG                                                    
      SL = SIN(TUA)                                                     
      CL = COS(TUA)                                                     
      FUT(3) = 1.                                                       
      FUT(2) = SQ2*CL                                                   
      FUT(4) = SQ2*SL                                                   
      FUT(1) = CL*FUT(2) - SL*FUT(4)                                    
      FUT(5) = CL*FUT(4) + SL*FUT(2)                                    
  603 CONTINUE                                                          
      IF (ICPT.EQ.0) GO TO 608                                          
      DO 607 MPF = 1,9                                                  
      PB(MPF) = 0.                                                      
      VB(MPF) = 0.                                                      
      DO 605 LP=1,5                                                     
      IF (IUTAV.NE.0.AND.LP.NE.3) GO TO 605                             
      PB(MPF) = PB(MPF) + FUT(LP)*PA(LP,MPF)                            
      VB(MPF) = VB(MPF) + FUT(LP)*VA(LP,MPF)                            
  605 CONTINUE                                                          
  607 CONTINUE                                                          
C  TL IS MAGNETIC LOCAL TIME.                                           
  608 TL = UT + (XMLON - 69.)/15.                                       
C  IF TL IS SAME AS IN PREVIOUS CALL TO THIS SUBROUTINE, SKIP TO STATEME
      IF (TL.EQ.TLP) GO TO 630                                          
      TLP = TL                                                          
C  (THROUGH STATEMENT 610) CALCULATE SINES AND COSINES, TIMES SQ2, OF HA
C    OF MAGNETIC LOCAL TIME.                                            
      TLA = TL*HRANG                                                    
      SL = SIN(TLA)                                                     
      CL = COS(TLA)                                                     
      SML(1) = SQ2*SL                                                   
      CML(1) = SQ2*CL                                                   
      DO 610 M=2,4                                                      
      SML(M) = CL*SML(M-1) + SL*CML(M-1)                                
  610 CML(M) = CL*CML(M-1) - SL*SML(M-1)                                
C  CALCULATE POT, VE, AND VU BY SUMMING FOURIER COEFFICIENTS MULTIPLIED 
C    APPROPRIATE SINES AND COSINES.                                     
  630 POT = 0.                                                          
      VE = 0.                                                           
      VU = 0.                                                           
      DO 640 M=1,4                                                      
      MPF = M+5                                                         
      MMF = 5 - M                                                       
      XM = M                                                            
      POT = POT + PB(MPF)*SML(M) + PB(MMF)*CML(M)                       
      VE = VE + VB(MPF)*SML(M) +   VB(MMF)*CML(M)                       
  640 VU = VU + XM*(PB(MPF)*CML(M) - PB(MMF)*SML(M))                    
      POT = POT*ST + PB(5)                                              
      VE = VE + VB(5)                                                   
      VU = VU/RBT                                                       
      GO TO 999                                                         
  988 POT = -1/IFIRST                                                   
      VU = -1/IFIRST                                                    
      VE = -1/IFIRST                                                    
  999 RETURN                                                            
      END                                                               
      SUBROUTINE ASOUT(A,NAME,IA,NROW,NCOL,NR,NC)                       
      DIMENSION A(IA,NCOL)                                              
C     SPARSE OUTPUT OF A 2-D ARRAY                                      
C     REAL NAME                                                         
      CHARACTER*(*) NAME                                                
C  REMOVE FOLLOWING STATEMENT TO ACTIVATE THIS SUBROUTINE.              
      RETURN                                                            
      PRINT 100,NAME                                                    
C 100 FORMAT(1H A8)                                                     
  100 FORMAT(1X,A8)                                                     
      DO 1 I=1,NROW,NR                                                  
      PRINT 200,I,(A(I,J),J=1,NCOL,NC)                                  
C 200 FORMAT(1H *ROW *I2,/(1H 10E10.3))                                 
  200 FORMAT(1X,'ROW ',I2,/(1X,10E10.3))                                
    1 CONTINUE                                                          
      RETURN                                                            
      END                                                               
      SUBROUTINE STDCON(VAR, CINT, IPRNT)                               
C                                                                       
C THIS SUBROUTINE SETS UP AXES, LABELS THEM,                            
C   LABELS THE PLOT, AND CALLS CONREC.                                  
C PARAMETERS ARE PASSED IN THE CALLING LIST AND IN THE                  
C   COMMON BLOCKS CONCAL AND XYCV. (XYCV IS ALSO IN CONREC)             
C THE INPUT PARAMETERS ARE ---                                          
C   VAR(MXR, MY) = THE ARRAY TO BE CONTOURED                            
C   CINT = CONTOUR INTERVAL. IF ZERO, WILL BE SET BY CONREC.            
C     (IS BEST TO SET CINT ONESELF TO INSURE DESIRABLE SPACING).        
C   IPRNT= FLAG TO PRINT OUT VAR ARRAY (ALONG WITH XX AND YY).          
C     IF= 0, WILL NOT PRINT OUT.                                        
C   MXR = X DIMENSION OF VAR                                            
C   MX, MY = THE NUMBER OF X AND Y POINTS OF VAR                        
C     THAT ARE ACTUALLY PLOTTED. MX. LE. MXR.                           
C   LTYPE= TYPE OF PLOT. LABELS ARE IN F5.0 FORMAT FOR                  
C   LINEAR AND IN E8.1 FORMAT FOR LOG.                                  
C        = 1, LINEAR X, LINEAR Y                                        
C        = 2, LINEAR X, LOG Y                                           
C        = 3, LOG X, LINEAR Y                                           
C        = 4, LOG X, LOG Y                                              
C   XC, YC= MIN VALUE ON X(Y) AXIS                                      
C   XD, YD= MAX VALUE ON X(Y) AXIS                                      
C   MAGX, MAGY = NO. OF LABELLED SPACES ON X(Y) AXIS                    
C     (MUST BE 1 FOR LOG PLOTS)                                         
C   MINX, MINY= NO. OF SPACES BETWEEN LABELLED                          
C     TICK MARKS ON X(Y)- AXIS                                          
C     (SHOULD BE 9 FOR LOG PLOTS)                                       
C   LABX, LABY= LABELS FOR X(Y) AXIS. 8 CHARACTERS                      
C   LABL(7)= TITLE ON TOP OF PLOT. FOR CRAY- 1,                         
C     IS 7*8= 56 CHARACTERS LONG                                        
C   XX,YY= X AND Y AXIS VALUES THAT CORRESPOND WITH MX AND MY POINTS    
C     IN VAR.                                                           
C                                                                       
      CHARACTER*6 LINFMT,LOGFMT                                         
      CHARACTER*8 LABX,LABY                                             
      CHARACTER*56 LABL                                                 
      COMMON/CVLX/LABX                                                  
      COMMON/CVLY/LABY                                                  
      COMMON/CVLBL/LABL                                                 
      COMMON/STDCV/MXR,MX,MY,LTYPE, XC, XD, MAGX, MINX,XLABX, YC,       
     1YD, MAGY, MINY,XLABY,XLABL                                        
      COMMON/XYCV/XX(100), YY(100)                                      
      DIMENSION VAR(MXR, MY)                                            
C                                                                       
      IF (MX .LT. 100 .AND. MY .LT. 100) GO TO 20                       
      WRITE (6, 100) MX, MY                                             
  100 FORMAT (1X,'TRANSFORMATION ARRAYS ARE UNDERDIMENSIONED            
     1MX AND MY ', 2I5)                                                 
      STOP                                                              
  20  IF(IPRNT.EQ.0) GO TO 10                                           
C     WRITE (6, 101) (LABL(I), I=1,7)                                   
C 101 FORMAT (1X// 1X, 7A8/)                                            
C     WRITE (6, 101) (LABL)                                             
C 101 FORMAT (1X// 1X, A56/)                                            
      WRITE (6,102) (XX(I),I=1,MX)                                      
  102 FORMAT (1X, 7X, 10E12.4)                                          
      WRITE (6,103)                                                     
  103 FORMAT (1X)                                                       
      DO 5 J=1,MY                                                       
    5 WRITE (6,104) YY(J), (VAR(I,J),I=1,MX)                            
  104 FORMAT (1X, F7.2, 10E12.4/ (8X, 10E12.4))                         
   10 CONTINUE                                                          
      CALL SET (0.1, 0.95, 0.2, 0.8, XC, XD, YC, YD, LTYPE)             
      CXA = (XD-XC)/(.85*1024.)                                         
C     CXB = (.95*XC-.1*XD)/(.85*1024.)                                  
      CXB = (.95*XC-.1*XD)/.85                                          
      CYA = (YD-YC)/(.60*1024.)                                         
C     CYB = (.8*YC-.2*YD)/(.60*1024.)                                   
      CYB = (.8*YC-.2*YD)/.60                                           
C     CALL PWRITY (105,860, LABL, 56, 2, 0, -1)                         
      CALL PWRITY (105.*CXA+CXB,860.*CYA+CYB, LABL, 56, 2, 0, -1)       
C     LINFMT= 6H(F5.0)                                                  
C     LOGFMT= 6H(E8.1)                                                  
      LINFMT= '(F5.0)'                                                  
      LOGFMT= '(E8.1)'                                                  
      GO TO (1, 2, 3, 4), LTYPE                                         
    1 CALL LABMOD (LINFMT, LINFMT, 5, 5, 1, 1, 0, 0, 0)                 
      GO TO 6                                                           
    2 CALL LABMOD (LINFMT, LOGFMT, 5, 8, 1, 1, 0, 0, 0)                 
      GO TO 6                                                           
    3 CALL LABMOD (LOGFMT, LINFMT, 8, 5, 1, 1, 0, 0, 0)                 
      GO TO 6                                                           
    4 CALL LABMOD (LOGFMT, LOGFMT,8, 8, 1, 1, 0, 0, 0)                  
    6 CALL PERIML (MAGX, MINX, MAGY, MINY)                              
C     CALL PWRY (15, 500,LABY, 8, 2.5, 1,1)                             
C     CALL PWRY (525, 140, LABX, 8, 2.5, 0, 1)                          
      CALL PWRITY (15.*CXA+CXB,.5*YD,LABY,5,2.5,90,0)                   
      CALL PWRITY (.5*(XC+XD),140.*CYA+CYB,LABX,5,2.5,0,0)              
C     CALL CONREC (VAR,MXR, MX,MY,0.,0.,CINT,1,0,-1635B)                
C SUPPRESS WRITING OF H AND L                                           
      CALL CPSETI('MAP',1)                                              
      CALL CPCNRC (VAR,MXR, MX,MY,0.,0.,CINT,1,-1,-1635B)               
      CALL FRAME                                                        
      RETURN                                                            
      END                                                               
      SUBROUTINE CPMPXY(IMAP,X,Y,FX,FY)                                 
      COMMON/XYCV/XX(100),YY(100)                                       
      FX     =XX(IFIX(X))+(XX(IFIX(X)+1)-XX(IFIX(X)))*(X-AINT(X))       
      FY     =YY(IFIX(Y))+(YY(IFIX(Y)+1)-YY(IFIX(Y)))*(Y-AINT(Y))       
      RETURN                                                            
      END                                                               
      SUBROUTINE STDVCT(U,V,VNORM)                                      
C                                                                       
C THIS SUBROUTINE SETS UP AXES, LABELS THEM,                            
C   LABELS THE PLOT, AND CALLS VELVCT.                                  
C PARAMETERS ARE PASSED IN THE CALLING LIST AND IN THE                  
C   COMMON BLOCKS STDCV AND XYCV. (XYCV IS IN VELVCT)                   
C THE INPUT PARAMETERS ARE ---                                          
C   U,V(MXR, MY) = THE ARRAYS TO BE PLOTTED IN THE X AND Y DIRECTIONS   
C   VNORM = VECTOR INTENSITY CORRESPONDING TO X-AXIS LENGTH             
C     IF ZERO, WILL BE SET BY VELVCT.                                   
C     (IS BEST TO SET VNORM ONESELF TO INSURE DESIRABLE LENGTHS).       
C   MXR = X DIMENSION OF VAR                                            
C   MX, MY = THE NUMBER OF X AND Y POINTS OF VAR                        
C     THAT ARE ACTUALLY PLOTTED. MX. LE. MXR.                           
C   NOCHNGE = UNUSED COMMON VARIABLE.                                   
C   XC, YC= MIN VALUE ON X(Y) AXIS                                      
C   XD, YD= MAX VALUE ON X(Y) AXIS                                      
C   MAGX, MAGY = NO. OF LABELLED SPACES ON X(Y) AXIS                    
C   MINX, MINY= NO. OF SPACES BETWEEN LABELLED                          
C     TICK MARKS ON X(Y)- AXIS                                          
C   LABX, LABY= LABELS FOR X(Y) AXIS. 8 CHARACTERS                      
C   LABL(7)= TITLE ON TOP OF PLOT. FOR CRAY- 1,                         
C     IS 7*8= 56 CHARACTERS LONG                                        
C   XX,YY: X AND Y AXIS VALUES THAT CORRESPOND WITH MX AND MY POINTS    
C     IN VAR.                                                           
C   FLO: MINIMUM VECTOR LENGTH TO BE PLOTTED.  SET HERE TO .002*VNORM.  
C   SPVL : IF EITHER U(I,J) = SPVL(1) OR V(I,J) = SPVL(2) THE VECTOR    
C     WILL NOT BE PLOTTED .                                             
C                                                                       
      CHARACTER*6 LINFMT,LOGFMT                                         
      CHARACTER*8 LABX,LABY                                             
      CHARACTER*56 LABL                                                 
      COMMON/CVLX/LABX                                                  
      COMMON/CVLY/LABY                                                  
      COMMON/CVLBL/LABL                                                 
      COMMON/STDCV/MXR,MX,MY,NOCHNGE, XC, XD, MAGX, MINX,XLABX, YC,     
     1YD, MAGY, MINY,XLABY,XLABL                                        
      DIMENSION U(MXR, MY),V(MXR,MY),SPV(2)                             
      LTYPE = 1                                                         
      SPVL = 1.E30                                                      
      SPV(1) = SPVL                                                     
      SPV(2) = SPVL                                                     
      HI = .125* VNORM                                                  
      LENGTH = 109                                                      
      IF (HI.EQ.0.) LENGTH = 0                                          
      FLO = .002*VNORM                                                  
C                                                                       
      IF (MX .LT. 100 .AND. MY .LT. 100) GO TO 10                       
      WRITE (6, 100) MX, MY                                             
  100 FORMAT (1X,'TRANSFORMATION ARRAYS ARE UNDERDIMENSIONED            
     1MX AND MY ', 2I5)                                                 
      STOP                                                              
   10 CONTINUE                                                          
      CALL SET (0.1, 0.95, 0.2, 0.8, XC, XD, YC, YD, LTYPE)             
      CXA = (XD-XC)/(.85*1024.)                                         
C     CXB = (.95*XC-.1*XD)/(.85*1024.)                                  
      CXB = (.95*XC-.1*XD)/.85                                          
      CYA = (YD-YC)/(.60*1024.)                                         
C     CYB = (.8*YC-.2*YD)/(.60*1024.)                                   
      CYB = (.8*YC-.2*YD)/.60                                           
C     CALL PWRITY (105,860, LABL, 56, 2, 0, -1)                         
      CALL PWRITY (105.*CXA+CXB,860.*CYA+CYB, LABL, 56, 2, 0, -1)       
      LINFMT= '(F5.0)'                                                  
      CALL LABMOD (LINFMT, LINFMT, 5, 5, 1, 1, 0, 0, 0)                 
      CALL PERIML (MAGX, MINX, MAGY, MINY)                              
C     CALL PWRY (15, 500,LABY, 8, 2.5, 1,1)                             
C     CALL PWRY (525, 140, LABX, 8, 2.5, 0, 1)                          
      CALL PWRITY (15.*CXA+CXB,.5*YD,LABY,5,2.5,90,0)                   
      CALL PWRITY (.5*(XC+XD),140.*CYA+CYB,LABX,5,2.5,0,0)              
      CALL SET (.1,.95,.2,.8,1.,REAL(MX),.5,.5+REAL(MY), LTYPE)         
      CALL VELVCT(U,MXR,V,MXR, MX,MY,FLO,HI,1,LENGTH,3,SPV)             
      CALL FRAME                                                        
      MID = MY/2                                                        
      X = 0.                                                            
      N = 0                                                             
      DO 40 I=1,MX                                                      
      IF (V(I,MID).EQ.SPVL) GO TO 40                                    
      N = N + 1                                                         
      X = X + V(I,MID)                                                  
   40 CONTINUE                                                          
      X = X/FLOAT(N)                                                    
      PRINT 50, X                                                       
C  50 FORMAT (1H ,34HAVERAGE UPWARD COMPONENT AT MY/2 =,E10.3)          
   50 FORMAT (1X,'AVERAGE UPWARD COMPONENT AT MY/2 =',E10.3)            
      RETURN                                                            
      END                                                               
      BLOCK DATA VELDAT                                                 
C THIS 'ROUTINE' DEFINES THE DEFAULT VALUES OF THE VELVCT PARAMETERS.   
      COMMON /VEC1/   ASH        ,EXT        ,ICTRFG     ,ILAB       ,  
     +                IOFFD      ,IOFFM      ,ISX        ,ISY        ,  
     +                RMN        ,RMX        ,SIDE       ,SIZE       ,  
     +                XLT        ,YBT        ,ZMN        ,ZMX           
      COMMON /VEC2/   BIG        ,INCX       ,INCY                      
      DATA     EXT /    0.25 /                                          
      DATA  ICTRFG /    0    /                                          
      DATA    ILAB /    0    /                                          
      DATA   IOFFD /    0    /                                          
      DATA   IOFFM /    0    /                                          
      DATA     RMN /  160.00 /                                          
      DATA     RMX / 6400.00 /                                          
      DATA    SIDE /    0.90 /                                          
      DATA    SIZE /  256.00 /                                          
      DATA     XLT /    0.05 /                                          
      DATA     YBT /    0.05 /                                          
      DATA     ZMX /    0.00 /                                          
      DATA    INCX /    1     /                                         
      DATA    INCY /    1     /                                         
      END                                                               
C EQUINOX DRIFTS (6 LINES FOR POLEWARD, THEN 6 FOR EASTWARD) FOR MH,    
C  ARECIBO, JICAMARCA.  FIRST VALUE IS 0.5 UT.                          
