C SHELLIG.FOR, Version 2.0, January 1992
C
C 11/1/91 SHELLG: lowest starting point for B0 search is 2  
C 1/27/92 Adopted to IGRF-91 coeffcients model
C 2/5/92  Reduce variable-names: INTER(P)SHC,EXTRA(P)SHC,INITI(ALI)ZE
C
C*********************************************************************
C  SUBROUTINES FINDB0, SHELLG, STOER, FELDG, FELDCOF, GETSHC,        *
C	INTERSHC, EXTRASHC, INITIZE                                  *
C*********************************************************************
C*********************************************************************
C
C
      SUBROUTINE FINDB0(STPS,BDEL,VALUE,BEQU,RR0)
C--------------------------------------------------------------------
C FINDS SMALLEST MAGNETIC FIELD STRENGTH ON FIELD LINE
C
C INPUT:   STPS   STEP SIZE FOR FIELD LINE TRACING
C  	COMMON/FIDB0/
C	   SP     DIPOLE ORIENTED COORDINATES FORM SHELLG; P(1,*),
C 		  P(2,*), P(3,*) CLOSEST TO MAGNETIC EQUATOR 
C	   BDEL   REQUIRED ACCURACY  = [ B(LAST) - BEQU ] / BEQU  
C		  B(LAST)  IS FIELD STRENGTH BEFORE BEQU
C
C OUTPUT:  VALUE  =.FALSE., IF BEQU IS NOT MINIMAL VALUE ON FIELD LINE
C	   BEQU	  MAGNETIC FIELD STRENGTH AT MAGNETIC EQUATOR
C	   RR0	  EQUATORIAL RADIUS NORMALIZED TO EARTH RADIUS
C	   BDEL	  FINAL ACHIEVED ACCURACY
C--------------------------------------------------------------------
      IMPLICIT NONE
c      DIMENSION 	P(8,4),SP(3)
c      LOGICAL		VALUE
c      COMMON/FIDB0/	SP
c
c --- added by RY ---
c
      real SP,STPS,BDEL,BEQU,RR0,STEP,P,BQ2,R2,BQ3,R3
      real BQ1,R1,ZZ,BMIN,BOLD,B,SQRT,ROLD,BDELTA,STEP12
      integer I,J,N,IRUN

      DIMENSION        P(8,4),SP(3)
      LOGICAL          VALUE
      COMMON/FIDB0/    SP

C
      BOLD = 0.0
      ROLD = 0.0
      BDELTA = 0.0
C
      STEP=STPS
      IRUN=0
7777  IRUN=IRUN+1
      IF(IRUN.GT.5) THEN
	VALUE=.FALSE.
	GOTO 8888
	ENDIF
C*********************FIRST THREE POINTS 
      P(1,2)=SP(1)
      P(2,2)=SP(2)
      P(3,2)=SP(3)
      STEP=-SIGN(STEP,P(3,2))
      CALL STOER(P(1,2),BQ2,R2)
      P(1,3)=P(1,2)+0.5*STEP*P(4,2)
      P(2,3)=P(2,2)+0.5*STEP*P(5,2)
      P(3,3)=P(3,2)+0.5*STEP
      CALL STOER(P(1,3),BQ3,R3)
      P(1,1)=P(1,2)-STEP*(2.*P(4,2)-P(4,3))
      P(2,1)=P(2,2)-STEP*(2.*P(5,2)-P(5,3))
      P(3,1)=P(3,2)-STEP
      CALL STOER(P(1,1),BQ1,R1)
      P(1,3)=P(1,2)+STEP*(20.*P(4,3)-3.*P(4,2)+P(4,1))/18.
      P(2,3)=P(2,2)+STEP*(20.*P(5,3)-3.*P(5,2)+P(5,1))/18.
      P(3,3)=P(3,2)+STEP
      CALL STOER(P(1,3),BQ3,R3)
C******************INVERT SENSE IF REQUIRED
      IF(BQ3.LE.BQ1) GOTO 2
      	STEP=-STEP
      	R3=R1
      	BQ3=BQ1
      	DO 1 I=1,5
      		ZZ=P(I,1)
      		P(I,1)=P(I,3)
1     		P(I,3)=ZZ
C******************INITIALIZATION 
2	STEP12=STEP/12.
	VALUE=.TRUE.
	BMIN=1.E4
	BOLD=1.E4
C******************CORRECTOR (FIELD LINE TRACING)
	N=0
5555  P(1,3)=P(1,2)+STEP12*(5.*P(4,3)+8.*P(4,2)-P(4,1))
	N=N+1
      P(2,3)=P(2,2)+STEP12*(5.*P(5,3)+8.*P(5,2)-P(5,1))
C******************PREDICTOR (FIELD LINE TRACING)
      P(1,4)=P(1,3)+STEP12*(23.*P(4,3)-16.*P(4,2)+5.*P(4,1))
      P(2,4)=P(2,3)+STEP12*(23.*P(5,3)-16.*P(5,2)+5.*P(5,1))
      P(3,4)=P(3,3)+STEP
      CALL STOER(P(1,4),BQ3,R3)
	DO 1111 J=1,3
	DO 1111 I=1,8
1111	P(I,J)=P(I,J+1)
	B=SQRT(BQ3)
	IF(B.LT.BMIN) BMIN=B
	IF(B.LE.BOLD) THEN
		BOLD=B
		ROLD=1./R3
		SP(1)=P(1,4)
		SP(2)=P(2,4)
		SP(3)=P(3,4)
		GOTO 5555
		ENDIF
	IF(BOLD.NE.BMIN) THEN
		VALUE=.FALSE.
		ENDIF
	BDELTA=(B-BOLD)/BOLD
	IF(BDELTA.GT.BDEL) THEN
		STEP=STEP/10.
		GOTO 7777
		ENDIF
8888	RR0=ROLD
	BEQU=BOLD
	BDEL=BDELTA
	RETURN
	END
C
C
      SUBROUTINE SHELLG(GLAT,GLON,ALT,DIMO,FL,ICODE,B0)
C--------------------------------------------------------------------
C CALCULATES L-VALUE FOR SPECIFIED GEODAETIC COORDINATES, ALTITUDE
C AND GEMAGNETIC FIELD MODEL.
C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTER, INTERNAL NOTE 
C      NO. 67, 1970.
C      G. KLUGE, COMPUTER PHYSICS COMMUNICATIONS 3, 31-35, 1972
C--------------------------------------------------------------------
C CHANGES (D. BILITZA, NOV 87):
C   - USING CORRECT DIPOL MOMENT I.E.,DIFFERENT COMMON/MODEL/
C   - USING IGRF EARTH MAGNETIC FIELD MODELS FROM 1945 TO 1990
C--------------------------------------------------------------------
C  INPUT:  ENTRY POINT SHELLG
C	 	GLAT  GEODETIC LATITUDE IN DEGREES (NORTH)
C         	GLON  GEODETIC LONGITUDE IN DEGREES (EAST)
C         	ALT   ALTITUDE IN KM ABOVE SEA LEVEL
C
C	   ENTRY POINT SHELLC
C		V(3)  CARTESIAN COORDINATES IN EARTH RADII (6371.2 KM)
C			X-AXIS POINTING TO EQUATOR AT 0 LONGITUDE
C			Y-AXIS POINTING TO EQUATOR AT 90 LONG.
C			Z-AXIS POINTING TO NORTH POLE
C
C	   DIMO	      DIPOL MOMENT IN GAUSS (NORMALIZED TO EARTH RADIUS) 
C
C	   COMMON 
C		X(3)	NOT USED
C		H(144)	FIELD MODEL COEFFICIENTS ADJUSTED FOR SHELLG
C-----------------------------------------------------------------------
C  OUTPUT: FL   	L-VALUE
C	   ICODE  	=1 NORMAL COMPLETION
C			=2 UNPHYSICAL CONJUGATE POINT (FL MEANINGLESS)
C			=3 SHELL PARAMETER GREATER THAN LIMIT UP TO
C			   WHICH ACCURATE CALCULATION IS REQUIRED;
C			   APPROXIMATION IS USED.
C 	   B0   	MAGNETIC FIELD STRENGTH IN GAUSS
C-----------------------------------------------------------------------
      IMPLICIT NONE
c      DIMENSION 	V(3),U(3,3),P(8,100),SP(3)
c      COMMON 		X(3),H(144)
c      COMMON/FIDB0/	SP
c      COMMON/GENER/	UMR,ERA,AQUAD,BQUAD

c
c --- added by RY ---
c
      real SP,X,H,UMR,ERA,AQUAD,BQUAD,GLAT,GLON,ALT,DIMO
      real FL,B0,RMIN,RMAX,STEP,BEQU,CT,RLAT,ST,D,V,U,R3
      real RQ,P,BQ3,BQ1,R1,ZZ,FI,ORADIK,OTERM
      real STP,C1,D0,E0,T,ZQ,Z,R,HLI,FF,RADIK,TERM
      real STEQ,RLON,R3H,BQ2,R2,STEP12,STEP2,C0,C2,C3
      real D1,E1,E2,XX,D2,DIMOB0,GG
      integer ICODE,I,N,IEQU

      DIMENSION         V(3),U(3,3),P(8,100),SP(3)
      COMMON            X(3),H(144)
      COMMON/FIDB0/     SP
      COMMON/GENER/     UMR,ERA,AQUAD,BQUAD

C
C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3
C-- STEP IS STEP SIZE FOR FIELD LINE TRACING
C-- STEQ IS STEP SIZE FOR INTEGRATION
C 
      DATA RMIN,RMAX	/0.05,1.01/
      DATA STEP,STEQ	/0.20,0.03/
	BEQU=1.E10
C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES        SHEL0080
      RLAT=GLAT*UMR
      CT=SIN(RLAT)                                                      
      ST=COS(RLAT)                                                      
      D=SQRT(AQUAD-(AQUAD-BQUAD)*CT*CT)
      X(1)=(ALT+AQUAD/D)*ST/ERA
      X(3)=(ALT+BQUAD/D)*CT/ERA
      RLON=GLON*UMR
      X(2)=X(1)*SIN(RLON)                                               
      X(1)=X(1)*COS(RLON)                                               
      GOTO9                                                             
      ENTRY SHELLC(V,FL,B0)                                             
C*****ENTRY POINT  SHELLC  TO BE USED WITH CARTESIAN CO-ORDINATES       SHEL0200
      X(1)=V(1)                                                         
      X(2)=V(2)                                                         
      X(3)=V(3)                                                         
C*****CONVERT TO DIPOL-ORIENTED CO-ORDINATES                            SHEL0240
      DATA U/                 +0.3511737,-0.9148385,-0.1993679,         
     A                        +0.9335804,+0.3583680,+0.0000000,         
     B                        +0.0714471,-0.1861260,+0.9799247/         
9     RQ=1./(X(1)*X(1)+X(2)*X(2)+X(3)*X(3))
      R3H=SQRT(RQ*SQRT(RQ))                                             
      P(1,2)=(X(1)*U(1,1)+X(2)*U(2,1)+X(3)*U(3,1))*R3H                  
      P(2,2)=(X(1)*U(1,2)+X(2)*U(2,2)            )*R3H                  
      P(3,2)=(X(1)*U(1,3)+X(2)*U(2,3)+X(3)*U(3,3))*RQ                   
C*****FIRST THREE POINTS OF FIELD LINE                                  SHEL0330
      STEP=-SIGN(STEP,P(3,2))                                           
      CALL STOER(P(1,2),BQ2,R2)                                         
      B0=SQRT(BQ2)                                                      
      P(1,3)=P(1,2)+0.5*STEP*P(4,2)                                     
      P(2,3)=P(2,2)+0.5*STEP*P(5,2)                                     
      P(3,3)=P(3,2)+0.5*STEP                                            
      CALL STOER(P(1,3),BQ3,R3)                                         
      P(1,1)=P(1,2)-STEP*(2.*P(4,2)-P(4,3))                             
      P(2,1)=P(2,2)-STEP*(2.*P(5,2)-P(5,3))                             
      P(3,1)=P(3,2)-STEP                                                
      CALL STOER(P(1,1),BQ1,R1)                                         
      P(1,3)=P(1,2)+STEP*(20.*P(4,3)-3.*P(4,2)+P(4,1))/18.              
      P(2,3)=P(2,2)+STEP*(20.*P(5,3)-3.*P(5,2)+P(5,1))/18.              
      P(3,3)=P(3,2)+STEP                                                
      CALL STOER(P(1,3),BQ3,R3)                                         
C*****INVERT SENSE IF REQUIRED                                          SHEL0490
      IF(BQ3.LE.BQ1)GOTO2                                               
      STEP=-STEP                                                        
      R3=R1                                                             
      BQ3=BQ1                                                           
      DO 1 I=1,7                                                        
      ZZ=P(I,1)                                                         
      P(I,1)=P(I,3)                                                     
1     P(I,3)=ZZ                                                         
C*****SEARCH FOR LOWEST MAGNETIC FIELD STRENGTH
2     IF(BQ1.LT.BEQU) THEN
	BEQU=BQ1
	IEQU=1
	ENDIF
      IF(BQ2.LT.BEQU) THEN
	BEQU=BQ2
	IEQU=2
	ENDIF
      IF(BQ3.LT.BEQU) THEN
	BEQU=BQ3
	IEQU=3
	ENDIF
C*****INITIALIZATION OF INTEGRATION LOOPS                               SHEL0580
      STEP12=STEP/12.
      STEP2=STEP+STEP                                                   
      STEQ=SIGN(STEQ,STEP)                                              
      FI=0.                                                             
      ICODE=1                                                           
      ORADIK=0.                                                         
      OTERM=0.                                                          
      STP=R2*STEQ                                                       
      Z=P(3,2)+STP                                                      
      STP=STP/0.75
      P(8,1)=STEP2*(P(1,1)*P(4,1)+P(2,1)*P(5,1))                        
      P(8,2)=STEP2*(P(1,2)*P(4,2)+P(2,2)*P(5,2))                        
C*****MAIN LOOP (FIELD LINE TRACING)                                    SHEL0710
      DO 3 N=3,3333                                                     
C*****CORRECTOR (FIELD LINE TRACING)                                    SHEL0730
      P(1,N)=P(1,N-1)+STEP12*(5.*P(4,N)+8.*P(4,N-1)-P(4,N-2))           
      P(2,N)=P(2,N-1)+STEP12*(5.*P(5,N)+8.*P(5,N-1)-P(5,N-2))           
C*****PREPARE EXPANSION COEFFICIENTS FOR INTERPOLATION                  SHEL0760
C*****OF SLOWLY VARYING QUANTITIES                                      SHEL0770
      P(8,N)=STEP2*(P(1,N)*P(4,N)+P(2,N)*P(5,N))                        
      C0=P(1,N-1)**2+P(2,N-1)**2                                        
      C1=P(8,N-1)                                                       
      C2=(P(8,N)-P(8,N-2))*0.25                                         
      C3=(P(8,N)+P(8,N-2)-C1-C1)/6.0
      D0=P(6,N-1)                                                       
      D1=(P(6,N)-P(6,N-2))*0.5                                          
      D2=(P(6,N)+P(6,N-2)-D0-D0)*0.5                                    
      E0=P(7,N-1)
      E1=(P(7,N)-P(7,N-2))*0.5                                          
      E2=(P(7,N)+P(7,N-2)-E0-E0)*0.5                                    

C*****INNER LOOP (FOR QUADRATURE)                                       SHEL0890
4     T=(Z-P(3,N-1))/STEP                                               
      IF(T.GT.1.)GOTO5                                                  
      HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                                   
      ZQ=Z*Z
      R=HLI+SQRT(HLI*HLI+ZQ)
      IF(R.LE.RMIN)GOTO30                                               
      RQ=R*R
      FF=SQRT(1.+3.*ZQ/RQ)                                              
      RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF                                 
      IF(R-RMAX)44,44,45                                                
45    ICODE=2                                                           
      RADIK=RADIK-12.*(R-RMAX)**2                                       
44    IF(RADIK+RADIK.LE.ORADIK) GOTO 10
      TERM=SQRT(RADIK)*FF*((E2*T+E1)*T+E0)/(RQ+ZQ)                      
      FI=FI+STP*(OTERM+TERM)                                            
      ORADIK=RADIK                                                      
      OTERM=TERM                                                        
      STP=R*STEQ                                                        
      Z=Z+STP                                                           
      GOTO4                                                             
C*****PREDICTOR (FIELD LINE TRACING)                                    SHEL1100
5     P(1,N+1)=P(1,N)+STEP12*(23.*P(4,N)-16.*P(4,N-1)+5.*P(4,N-2))      
      P(2,N+1)=P(2,N)+STEP12*(23.*P(5,N)-16.*P(5,N-1)+5.*P(5,N-2))      
      P(3,N+1)=P(3,N)+STEP                                              
      CALL STOER(P(1,N+1),BQ3,R3)                                       
C*****SEARCH FOR LOWEST MAGNETIC FIELD STRENGTH
      IF(BQ3.LT.BEQU) THEN
	IEQU=N+1
	BEQU=BQ3
	ENDIF
3     CONTINUE
10    IF(IEQU.lt.2) IEQU=2 
      SP(1)=P(1,IEQU-1)
      SP(2)=P(2,IEQU-1)
      SP(3)=P(3,IEQU-1)
      IF(ORADIK.LT.1E-15)GOTO11                                         
      FI=FI+STP/0.75*OTERM*ORADIK/(ORADIK-RADIK)              
C
C-- The minimal allowable value of FI was changed from 1E-15 to 1E-12,
C-- because 1E-38 is the minimal allowable arg. for ALOG in our envir.
C-- D. Bilitza, Nov 87.
C
11    FI=0.5*ABS(FI)/SQRT(B0)+1E-12                       
C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.  
C
C-- Correct dipole moment is used here. D. Bilitza, Nov 87.
C
      DIMOB0=DIMO/B0
      XX=ALOG(FI*FI*FI/DIMOB0)
      IF(XX.GT.23.0) GOTO 776   
      IF(XX.GT.11.7) GOTO 775  
      IF(XX.GT.+3.0) GOTO 774    
      IF(XX.GT.-3.0) GOTO 773   
      IF(XX.GT.-22.) GOTO 772  
  771 GG=3.33338E-1*XX+3.0062102E-1                                     
      GOTO777                                                           
  772 GG=((((((((-8.1537735E-14*XX+8.3232531E-13)*XX+1.0066362E-9)*XX+  
     18.1048663E-8)*XX+3.2916354E-6)*XX+8.2711096E-5)*XX+1.3714667E-3)* 
     2XX+1.5017245E-2)*XX+4.3432642E-1)*XX+6.2337691E-1                 
      GOTO777                                                           
  773 GG=((((((((2.6047023E-10*XX+2.3028767E-9)*XX-2.1997983E-8)*XX-    
     15.3977642E-7)*XX-3.3408822E-6)*XX+3.8379917E-5)*XX+1.1784234E-3)* 
     2XX+1.4492441E-2)*XX+4.3352788E-1)*XX+6.228644E-1                  
      GOTO777                                                           
  774 GG=((((((((6.3271665E-10*XX-3.958306E-8)*XX+9.9766148E-07)*XX-    
     11.2531932E-5)*XX+7.9451313E-5)*XX-3.2077032E-4)*XX+2.1680398E-3)* 
     2XX+1.2817956E-2)*XX+4.3510529E-1)*XX+6.222355E-1                  
      GOTO777                                                           
  775 GG=(((((2.8212095E-8*XX-3.8049276E-6)*XX+2.170224E-4)*XX-6.7310339
     1E-3)*XX+1.2038224E-1)*XX-1.8461796E-1)*XX+2.0007187E0             
      GOTO777                                                           
  776 GG=XX-3.0460681E0                                                 
  777 FL=EXP(ALOG((1.+EXP(GG))*DIMOB0)/3.0)
      RETURN                                                            
C*****APPROXIMATION FOR HIGH VALUES OF L.                               SHEL1450
30    ICODE=3                                                           
      T=-P(3,N-1)/STEP                                                  
      FL=1./(ABS(((C3*T+C2)*T+C1)*T+C0)+1E-15)                          
      RETURN                                                            
      END                                                               
C
C
      SUBROUTINE STOER(P,BQ,R)                                          
C*******************************************************************
C* SUBROUTINE USED FOR FIELD LINE TRACING IN SHELLG                *
C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *
C*******************************************************************
      IMPLICIT NONE
c      DIMENSION 	P(7),U(3,3)
c      COMMON 		XI(3),H(144)
c
c --- added by RY - define all variables ---
c
      real XI,H,P,BQ,R,ZM,FLI,U,Q,DX,DY,DZ,RQ,WR,XM,YM
      real DXM,DYM,DZM,DR,DSQ

      DIMENSION         P(7),U(3,3)
      COMMON            XI(3),H(144)

C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES          SHEL1540
      ZM=P(3)                                                           
      FLI=P(1)*P(1)+P(2)*P(2)+1E-15
      R=0.5*(FLI+SQRT(FLI*FLI+(ZM+ZM)**2))
      RQ=R*R
      WR=SQRT(R)                                                        
      XM=P(1)*WR                                                        
      YM=P(2)*WR                                                        
C*****TRANSFORM TO GEOGRAPHIC CO-ORDINATE SYSTEM                        SHEL1620
      DATA U/                 +0.3511737,-0.9148385,-0.1993679,         
     A                        +0.9335804,+0.3583680,+0.0000000,         
     B                        +0.0714471,-0.1861260,+0.9799247/         
      XI(1)=XM*U(1,1)+YM*U(1,2)+ZM*U(1,3)                               
      XI(2)=XM*U(2,1)+YM*U(2,2)+ZM*U(2,3)                               
      XI(3)=XM*U(3,1)          +ZM*U(3,3)                               
C*****COMPUTE DERIVATIVES                                               SHEL1690
      CALL FELDI(XI,H)                                                  
      Q=H(1)/RQ                                                         
      DX=H(3)+H(3)+Q*XI(1)                                              
      DY=H(4)+H(4)+Q*XI(2)                                              
      DZ=H(2)+H(2)+Q*XI(3)                                              
C*****TRANSFORM BACK TO GEOMAGNETIC CO-ORDINATE SYSTEM                  SHEL1750
      DXM=U(1,1)*DX+U(2,1)*DY+U(3,1)*DZ                                 
      DYM=U(1,2)*DX+U(2,2)*DY                                           
      DZM=U(1,3)*DX+U(2,3)*DY+U(3,3)*DZ                                 
      DR=(XM*DXM+YM*DYM+ZM*DZM)/R                                       
C*****FORM SLOWLY VARYING EXPRESSIONS                                   SHEL1800
      P(4)=(WR*DXM-0.5*P(1)*DR)/(R*DZM)                                 
      P(5)=(WR*DYM-0.5*P(2)*DR)/(R*DZM)                                 
      DSQ=RQ*(DXM*DXM+DYM*DYM+DZM*DZM)
      BQ=DSQ*RQ*RQ
      P(6)=SQRT(DSQ/(RQ+3.*ZM*ZM))                                      
      P(7)=P(6)*(RQ+ZM*ZM)/(RQ*DZM)                                     
      RETURN                                                            
      END                                                               
C
C
      SUBROUTINE FELDG(GLAT,GLON,ALT,BNORTH,BEAST,BDOWN,BABS)           
C-------------------------------------------------------------------
C CALCULATES EARTH MAGNETIC FIELD FROM SPHERICAL HARMONICS MODEL
C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61, 
C      1970.
C--------------------------------------------------------------------
C CHANGES (D. BILITZA, NOV 87):
C   - FIELD COEFFICIENTS IN BINARY DATA FILES INSTEAD OF BLOCK DATA
C   - CALCULATES DIPOL MOMENT
C--------------------------------------------------------------------
C  INPUT:  ENTRY POINT FELDG
C	 	GLAT  GEODETIC LATITUDE IN DEGREES (NORTH)
C         	GLON  GEODETIC LONGITUDE IN DEGREES (EAST)
C         	ALT   ALTITUDE IN KM ABOVE SEA LEVEL
C
C	   ENTRY POINT FELDC
C		V(3)  CARTESIAN COORDINATES IN EARTH RADII (6371.2 KM)
C			X-AXIS POINTING TO EQUATOR AT 0 LONGITUDE
C			Y-AXIS POINTING TO EQUATOR AT 90 LONG.
C			Z-AXIS POINTING TO NORTH POLE
C
C	   COMMON BLANK AND ENTRY POINT FELDI ARE NEEDED WHEN USED
C	     IN CONNECTION WITH L-CALCULATION PROGRAM SHELLG.
C	
C	   COMMON /MODEL/ AND /GENER/
C		UMR     = ATAN(1.0)*4./180.   <DEGREE>*UMR=<RADIANT>
C		ERA	EARTH RADIUS FOR NORMALIZATION OF CARTESIAN 
C			COORDINATES (6371.2 KM)
C		AQUAD, BQUAD   SQUARE OF MAJOR AND MINOR HALF AXIS FOR 
C			EARTH ELLIPSOID AS RECOMMENDED BY INTERNATIONAL 
C			ASTRONOMICAL UNION (6378.160, 6356.775 KM).
C		NMAX    MAXIMUM ORDER OF SPHERICAL HARMONICS
C		TIME	YEAR (DECIMAL: 1973.5) FOR WHICH MAGNETIC 
C			FIELD IS TO BE CALCULATED
C		G(M)	NORMALIZED FIELD COEFFICIENTS (SEE FELDCOF)
C			M=NMAX*(NMAX+2)
C------------------------------------------------------------------------
C  OUTPUT: BABS   MAGNETIC FIELD STRENGTH IN GAUSS
C	   BNORTH, BEAST, BDOWN   COMPONENTS OF THE FIELD WITH RESPECT
C		  TO THE LOCAL GEODETIC COORDINATE SYSTEM, WITH AXIS
C		  POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
C		  AND DOWNWARD.   
C-----------------------------------------------------------------------
      IMPLICIT NONE
c      DIMENSION 	V(3),B(3)   
c      CHARACTER*12 	NAME
c      COMMON 		XI(3),H(144)
c      COMMON/MODEL/	NAME,NMAX,TIME,G(144)  
c      COMMON/GENER/	UMR,ERA,AQUAD,BQUAD
c
c --- added by RY - define variables ---
c
      real XI,H,UMR,ERA,AQUAD,BQUAD,TIME,G,GLAT,GLON
      real ALT,BNORTH,BEAST,BDOWN,BABS,CT,RLAT,ST
      real D,CP,SP,V,B,XXX,YYY,F,T,RQ,RLON,ZZZ,RHO,X,Y,Z
      real BXXX,BZZZ,BYYY,S,BRHO
      integer I,IHMAX,LAST,K,IMAX,M,IL,IS,IH,NMAX

      DIMENSION         V(3),B(3)
      CHARACTER*12      NAME
      COMMON            XI(3),H(144)
      COMMON/SHEL_MODEL/     NAME,NMAX,TIME,G(144)
      COMMON/GENER/     UMR,ERA,AQUAD,BQUAD

C
C-- IS RECORDS ENTRY POINT
C
C*****ENTRY POINT  FELDG  TO BE USED WITH GEODETIC CO-ORDINATES         SHEL1920
      IS=1                                                              
      RLAT=GLAT*UMR
      CT=SIN(RLAT)                                                      
      ST=COS(RLAT)                                                      
      D=SQRT(AQUAD-(AQUAD-BQUAD)*CT*CT)                                 
      RLON=GLON*UMR
      CP=COS(RLON)                                                      
      SP=SIN(RLON)                                                      
       ZZZ=(ALT+BQUAD/D)*CT/ERA
       RHO=(ALT+AQUAD/D)*ST/ERA
       XXX=RHO*CP                                                       
       YYY=RHO*SP                                                       
      GOTO10                                                            
      ENTRY FELDC(V,B)                                                  
C*****ENTRY POINT  FELDC  TO BE USED WITH CARTESIAN CO-ORDINATES        SHEL2070
      IS=2                                                              
      XXX=V(1)                                                          
      YYY=V(2)                                                          
      ZZZ=V(3)                                                          
10    RQ=1./(XXX*XXX+YYY*YYY+ZZZ*ZZZ) 
      XI(1)=XXX*RQ                                                      
      XI(2)=YYY*RQ                                                      
      XI(3)=ZZZ*RQ                                                      
      GOTO20                                                            
      ENTRY FELDI                                                       
C*****ENTRY POINT  FELDI  USED FOR L COMPUTATION                        SHEL2190
      IS=3                                                              
20    IHMAX=NMAX*NMAX+1                                                 
      LAST=IHMAX+NMAX+NMAX                                              
      IMAX=NMAX+NMAX-1                                                  
      DO 8 I=IHMAX,LAST                                                 
8     H(I)=G(I)                                                         
      DO 6 K=1,3,2                                                      
      I=IMAX                                                            
      IH=IHMAX                                                          
1     IL=IH-I                                                           
      F=2./FLOAT(I-K+2)                                                 
      X=XI(1)*F                                                         
      Y=XI(2)*F                                                         
      Z=XI(3)*(F+F)                                                     
      I=I-2                                                             
      IF(I-1)5,4,2                                                      
2     DO 3 M=3,I,2                                                      
      H(IL+M+1)=G(IL+M+1)+Z*H(IH+M+1)+X*(H(IH+M+3)-H(IH+M-1))           
     A                               -Y*(H(IH+M+2)+H(IH+M-2))           
3     H(IL+M)=G(IL+M)+Z*H(IH+M)+X*(H(IH+M+2)-H(IH+M-2))                 
     A                         +Y*(H(IH+M+3)+H(IH+M-1))                 
4     H(IL+2)=G(IL+2)+Z*H(IH+2)+X*H(IH+4)-Y*(H(IH+3)+H(IH))             
      H(IL+1)=G(IL+1)+Z*H(IH+1)+Y*H(IH+4)+X*(H(IH+3)-H(IH))             
5     H(IL)=G(IL)+Z*H(IH)+2.*(X*H(IH+1)+Y*H(IH+2))                      
      IH=IL                                                             
      IF(I.GE.K)GOTO1                                                   
6     CONTINUE                                                          
      IF(IS.EQ.3)RETURN                                                 
      S=.5*H(1)+2.*(H(2)*XI(3)+H(3)*XI(1)+H(4)*XI(2))                   
      T=(RQ+RQ)*SQRT(RQ)                                                
      BXXX=T*(H(3)-S*XXX)                                               
      BYYY=T*(H(4)-S*YYY)                                               
      BZZZ=T*(H(2)-S*ZZZ)                                               
      IF(IS.EQ.2)GOTO7                                                  
      BABS=SQRT(BXXX*BXXX+BYYY*BYYY+BZZZ*BZZZ)
      BEAST=BYYY*CP-BXXX*SP                                             
      BRHO=BYYY*SP+BXXX*CP                                              
      BNORTH=BZZZ*ST-BRHO*CT                                            
      BDOWN=-BZZZ*CT-BRHO*ST                                            
      RETURN                                                            
7     B(1)=BXXX                                                         
      B(2)=BYYY                                                         
      B(3)=BZZZ                                                         
      RETURN                                                            
      END                                                               
C
C
	SUBROUTINE FELDCOF(YEAR,DIMO,IER)
C------------------------------------------------------------------------
C  DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS
C
C	INPUT:  YEAR	DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO
C			BE CALCULATED
C	OUTPUT:	DIMO	GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED 
C			TO EARTH'S RADIUS) AT THE TIME (YEAR)
C  D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771, 
C	(301)286-9536   NOV 1987.
C-----------------------------------------------------------------------
        IMPLICIT NONE
	CHARACTER*12	FILMOD, FIL1, FIL2           
c	DIMENSION	GH1(144),GH2(120),GHA(144),FILMOD(11),DTEMOD(11)
c	DOUBLE PRECISION X,F0,F 
c	COMMON/MODEL/	FIL1,NMAX,TIME,GH1
c	COMMON/GENER/	UMR,ERAD,AQUAD,BQUAD
c	DATA		FILMOD /'dgrf45.dat', 'dgrf50.dat',            
c     1			'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',      
c     2			'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat',      
c     3			'dgrf85.dat', 'igrf90.dat','igrf90s.dat'/
c	DATA		DTEMOD / 1945., 1950., 1955., 1960., 1965.,           
c     1			1970., 1975., 1980., 1985., 1990., 1995./      
c
c --- added by RY ---
c
        real UMR,ERAD,AQUAD,BQUAD,TIME,GH1,YEAR,DIMO
        real DTEMOD,DTE1,DTE2,GH2,GHA,SQRT2
        integer IU,IS,IYEA,L,IER,J,I,N,M,NMAX,NUMYE
        integer NMAX1,NMAX2

        DIMENSION       GH1(144),GH2(144),GHA(144),FILMOD(11),DTEMOD(11)
        DOUBLE PRECISION X,F0,F
        COMMON/SHEL_MODEL/   FIL1,NMAX,TIME,GH1
        COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD
        DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',
     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',
     2                  'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat',
     3                  'dgrf85.dat', 'igrf90.dat','igrf90s.dat'/
        DATA            DTEMOD / 1945., 1950., 1955., 1960., 1965.,
     1                  1970., 1975., 1980., 1985., 1990., 1995./
        character*30 errstr
        character*80 subinfo
         
c
c  numye is number of years represented by IGRF models
c
	NUMYE=10
        errstr = ' ERROR : FELDCOF'
C
C  IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION
C  IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS
C
	IU = 10
	IS = 0
        do L=1,144
          GH1(L) = 0.0
          GH2(L) = 0.0
          GHA(L) = 0.0
        enddo
C-- DETERMINE IGRF-YEARS FOR INPUT-YEAR
	TIME = YEAR
	IYEA = INT(YEAR/5.)*5
	L = (IYEA - 1945)/5 + 1
	IF(L.LT.1) L=1
	IF(L.GT.NUMYE) L=NUMYE         
	DTE1 = DTEMOD(L)   
	FIL1 = FILMOD(L)   
	DTE2 = DTEMOD(L+1) 
	FIL2 = FILMOD(L+1) 
C-- GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
	CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
	    IF (IER .NE. 0) THEN
              subinfo = errstr//' getting IGRF coefficients'
              call fcecho(subinfo)
              return
              ENDIF 
	CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
	    IF (IER .NE. 0) THEN
              subinfo = errstr//' getting IGRF coefficients'
              call fcecho(subinfo) 
              return 
            ENDIF
                    
C-- DETERMINE IGRF COEFFICIENTS FOR YEAR
	IF (L .LE. NUMYE-1) THEN                        
	  CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2, 
     1		NMAX2, GH2, NMAX, GHA)                        
	ELSE               
	  CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,     
     1		GH2, NMAX, GHA)                                    
	ENDIF 
C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
	F0=0.D0
	DO 1234 J=1,3
	   F = GHA(J) * 1.D-5
	   F0 = F0 + F * F
1234	CONTINUE
	DIMO = DSQRT(F0)

	GH1(1) =  0.0
	I=2          
      	F0=1.D-5                
      	IF(IS.EQ.0) F0=-F0 
      	SQRT2=SQRT(2.)      

      DO 9 N=1,NMAX           
	X = N
      	F0 = F0 * X * X / (4.D0 * X - 2.D0)               
      	IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
	F = F0 * 0.5D0                                    
      	IF(IS.EQ.0) F = F * SQRT2
	GH1(I) = GHA(I-1) * F0
	I = I+1                                         
      DO 9 M=1,N                                    
      	F = F * (X + M) / (X - M + 1.D0)                 
	IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M))      	
	GH1(I) = GHA(I-1) * F
	GH1(I+1) = GHA(I) * F
        I=I+2
9     CONTINUE                                          
	RETURN
	END
C
C
	SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)           
                                                                        
C ===============================================================               
C                                                                               
C	Version 1.01                                                 
C                                                                               
C	Reads spherical harmonic coefficients from the specified     
C	file into an array.                                          
C                                                                               
C	Input:                                                       
C	    IU    - Logical unit number                              
C	    FSPEC - File specification                               
C                                                                               
C	Output:                                                      
C	    NMAX  - Maximum degree and order of model                
C	    ERAD  - Earth's radius associated with the spherical     
C		    harmonic coefficients, in the same units as      
C		    elevation                                        
C	    GH    - Schmidt quasi-normal internal spherical          
C		    harmonic coefficients                            
C	    IER   - Error number: =  0, no error                     
C				  = -2, records out of order         
C			     	  = FORTRAN run-time error number    
C                                                                    
C	A. Zunde                                                     
C	USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225    
C                                                                               
C ===============================================================               
                                                                        
      IMPLICIT NONE 
	CHARACTER	FSPEC*(*)                                    
	DIMENSION	GH(*)                                        
                                                                        
C ---------------------------------------------------------------               
C	Open coefficient file. Read past first header record.        
C	Read degree and order of model and Earth's radius.           
C ---------------------------------------------------------------               
c --- added by RY ---
c

      real ERAD,GH
      integer IU,NMAX,IER
c     integer NN,MM,N,M,G,H,I


c      OPEN (IU, FILE=comfile(1:len), STATUS='OLD', 
c     &     IOSTAT=IER, ERR=999)     
c	READ (IU, *, IOSTAT=IER, ERR=999)                            
c	READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD                 
                                                                        
C ---------------------------------------------------------------               
C	Read the coefficient file, arranged as follows:              
C                                                                               
C					N     M     G     H          
C					----------------------       
C				    /   1     0    GH(1)  -          
C				   /	1     1    GH(2) GH(3)       
C				  /	2     0    GH(4)  -          
C				 /	2     1    GH(5) GH(6)       
C	    NMAX*(NMAX+3)/2 	/	2     2    GH(7) GH(8)       
C	       records		\	3     0    GH(9)  -          
C				 \      .     .     .     .          
C				  \	.     .     .     .          
C	    NMAX*(NMAX+2)	   \	.     .     .     .          
C	    elements in GH	    \  NMAX  NMAX   .     .          
C                                                                               
C	N and M are, respectively, the degree and order of the       
C	coefficient.                                                 
C ---------------------------------------------------------------               
                                                                        
c	I = 0                                                        
c	DO 2211 NN = 1, NMAX                                              
c	    DO 2233 MM = 0, NN                                            
c		READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H         
c		IF (NN .NE. N .OR. MM .NE. M) THEN                   
c		    IER = -2                                         
c		    GOTO 999                                         
c		ENDIF                                                
c		I = I + 1                                            
c		GH(I) = G                                            
c		IF (M .NE. 0) THEN                                   
c		    I = I + 1                                        
c		    GH(I) = H                                        
c		ENDIF                                                
c2233 	    CONTINUE                                                    
c2211	CONTINUE                                                        
                                                                        
c999	CLOSE (IU)                                                   

c The above has been deleted as FTOOLS do NOT support data in ASCII
c format. The data is now stored as DATA statements, later it is
c anticipated that the data will be stored in FITS format files,
c accesible via the CALDB access software.

         IF (FSPEC.EQ.'dgrf45.dat') THEN
           call dgrf45_dat_dat(NMAX,ERAD,GH,IER)
         ELSEIF (FSPEC.EQ.'dgrf55.dat') THEN
           call dgrf55_dat_dat(NMAX,ERAD,GH,IER)
         ELSEIF (FSPEC.EQ.'dgrf60.dat') THEN
           call dgrf60_dat_dat(NMAX,ERAD,GH,IER)
         ELSEIF (FSPEC.EQ.'dgrf50.dat') THEN
           call dgrf50_dat_dat(NMAX,ERAD,GH,IER)
         ELSEIF (FSPEC.EQ.'dgrf70.dat') THEN
           call dgrf70_dat_dat(NMAX,ERAD,GH,IER)
         ELSEIF (FSPEC.EQ.'dgrf75.dat') THEN
           call dgrf75_dat_dat(NMAX,ERAD,GH,IER)
         ELSEIF (FSPEC.EQ.'dgrf65.dat') THEN
           call dgrf65_dat_dat(NMAX,ERAD,GH,IER)
         ELSEIF (FSPEC.EQ.'dgrf85.dat') THEN
           call dgrf85_dat_dat(NMAX,ERAD,GH,IER)
         ELSEIF (FSPEC.EQ.'dgrf80.dat') THEN
           call dgrf80_dat_dat(NMAX,ERAD,GH,IER)
         ELSEIF (FSPEC.EQ.'igrf90s.dat') THEN
           call igrf90s_dat_dat(NMAX,ERAD,GH,IER)
         ELSEIF (FSPEC.EQ.'igrf90.dat') THEN
           call igrf90_dat_dat(NMAX,ERAD,GH,IER)
         ENDIF
	RETURN                                                       
	END                                                          
C
C
	SUBROUTINE INTERSHC (DATE, DTE1, NMAX1, GH1, DTE2,          
     1			      NMAX2, GH2, NMAX, GH)                  
                                                                        
C ===============================================================               
C                                                                               
C	Version 1.01                                                 
C                                                                               
C	Interpolates linearly, in time, between two spherical        
C	harmonic models.                                             
C                                                                               
C	Input:                                                       
C	    DATE  - Date of resulting model (in decimal year)        
C	    DTE1  - Date of earlier model                            
C	    NMAX1 - Maximum degree and order of earlier model        
C	    GH1   - Schmidt quasi-normal internal spherical          
C		    harmonic coefficients of earlier model           
C	    DTE2  - Date of later model                              
C	    NMAX2 - Maximum degree and order of later model          
C	    GH2   - Schmidt quasi-normal internal spherical          
C		    harmonic coefficients of later model             
C                                                                               
C	Output:                                                      
C	    GH    - Coefficients of resulting model                  
C	    NMAX  - Maximum degree and order of resulting model      
C                                                                               
C	A. Zunde                                                     
C	USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225    
C                                                                               
C ===============================================================               
c
c --- added by RY ---
c
        real DATE,DTE2,GH1,DTE1,GH2
        real GH,FACTOR
        integer K,L,NMAX,I,NMAX1,NMAX2
        
                                                                        
	DIMENSION	GH1(*), GH2(*), GH(*)                        
                                                                        
C ---------------------------------------------------------------               
C	The coefficients (GH) of the resulting model, at date        
C	DATE, are computed by linearly interpolating between the     
C	coefficients of the earlier model (GH1), at date DTE1,       
C	and those of the later model (GH2), at date DTE2. If one     
C	model is smaller than the other, the interpolation is        
C	performed with the missing coefficients assumed to be 0.     
C ---------------------------------------------------------------               
                                                                        
	FACTOR = (DATE - DTE1) / (DTE2 - DTE1)                       
                                                                        
	IF (NMAX1 .EQ. NMAX2) THEN                                   
	    K = NMAX1 * (NMAX1 + 2)                                  
	    NMAX = NMAX1                                             
	ELSE IF (NMAX1 .GT. NMAX2) THEN                              
	    K = NMAX2 * (NMAX2 + 2)                                  
	    L = NMAX1 * (NMAX1 + 2)                                  
	    DO 1122 I = K + 1, L                                          
1122		GH(I) = GH1(I) + FACTOR * (-GH1(I))                  
	    NMAX = NMAX1                                             
	ELSE                                                         
	    K = NMAX1 * (NMAX1 + 2)                                  
	    L = NMAX2 * (NMAX2 + 2)                                  
	    DO 1133 I = K + 1, L                                          
1133		GH(I) = FACTOR * GH2(I)                              
	    NMAX = NMAX2                                             
	ENDIF                                                        
                                                                        
	DO 1144 I = 1, K                                                  
1144	    GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))              
                                                                        
	RETURN                                                       
	END                                                          
C
C
	SUBROUTINE EXTRASHC (DATE, DTE1, NMAX1, GH1, NMAX2,           
     1			      GH2, NMAX, GH)                           
                                                                        
C ===============================================================               
C                                                                               
C	Version 1.01                                                   
C                                                                               
C	Extrapolates linearly a spherical harmonic model with a        
C	rate-of-change model.                                          
C                                                                               
C	Input:                                                         
C	    DATE  - Date of resulting model (in decimal year)          
C	    DTE1  - Date of base model                                 
C	    NMAX1 - Maximum degree and order of base model             
C	    GH1   - Schmidt quasi-normal internal spherical            
C		    harmonic coefficients of base model                
C	    NMAX2 - Maximum degree and order of rate-of-change         
C		    model                                              
C	    GH2   - Schmidt quasi-normal internal spherical            
C		    harmonic coefficients of rate-of-change model      
C                                                                               
C	Output:                                                        
C	    GH    - Coefficients of resulting model                    
C	    NMAX  - Maximum degree and order of resulting model        
C                                                                               
C	A. Zunde                                                       
C	USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225      
C                                                                               
C ===============================================================               
c
c --- added by RY ---
c
        IMPLICIT NONE
        real DATE,DTE1,GH1,GH2,FACTOR,GH
        integer I,K,L,NMAX1,NMAX2,NMAX
                                                                        
	DIMENSION	GH1(*), GH2(*), GH(*)                          
                                                                        
C ---------------------------------------------------------------               
C	The coefficients (GH) of the resulting model, at date          
C	DATE, are computed by linearly extrapolating the coef-         
C	ficients of the base model (GH1), at date DTE1, using          
C	those of the rate-of-change model (GH2), at date DTE2. If      
C	one model is smaller than the other, the extrapolation is      
C	performed with the missing coefficients assumed to be 0.       
C ---------------------------------------------------------------               
                                                                        
	FACTOR = (DATE - DTE1)                                         
                                                                        
	IF (NMAX1 .EQ. NMAX2) THEN                                     
	    K = NMAX1 * (NMAX1 + 2)                                    
	    NMAX = NMAX1                                               
	ELSE IF (NMAX1 .GT. NMAX2) THEN                                
	    K = NMAX2 * (NMAX2 + 2)                                    
	    L = NMAX1 * (NMAX1 + 2)                                    
	    DO 1155 I = K + 1, L                                            
1155		GH(I) = GH1(I)                                         
	    NMAX = NMAX1                                               
	ELSE                                                           
	    K = NMAX1 * (NMAX1 + 2)                                    
	    L = NMAX2 * (NMAX2 + 2)                                    
	    DO 1166 I = K + 1, L                                            
1166		GH(I) = FACTOR * GH2(I)                                
	    NMAX = NMAX2                                               
	ENDIF                                                          
                                                                        
	DO 1177 I = 1, K                                                    
1177	    GH(I) = GH1(I) + FACTOR * GH2(I)                           
                                                                        
	RETURN                                                         
	END                                                            
C
C
	SUBROUTINE INITIZE
C----------------------------------------------------------------
C Initializes the parameters in COMMON/GENER/
C
C	UMR     = ATAN(1.0)*4./180.   <DEGREE>*UMR=<RADIANT>
C	ERA	EARTH RADIUS FOR NORMALIZATION OF CARTESIAN 
C			COORDINATES (6371.2 KM) 
C	EREQU	MAJOR HALF AXIS FOR EARTH ELLIPSOID (6378.160 KM)
C	ERPOL	MINOR HALF AXIS FOR EARTH ELLIPSOID (6356.775 KM)
C	AQUAD	SQUARE OF MAJOR HALF AXIS FOR EARTH ELLIPSOID
C	BQUAD   SQUARE OF MINOR HALF AXIS FOR EARTH ELLIPSOID
C
C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL 
C ASTRONOMICAL UNION .
C-----------------------------------------------------------------
c
c --- added by RY ---
c
        IMPLICIT NONE
        real UMR,ERA,AQUAD,BQUAD,EREQU,ERPOL

      	COMMON/GENER/	UMR,ERA,AQUAD,BQUAD
	ERA=6371.2
	EREQU=6378.16
	ERPOL=6356.775
	AQUAD=EREQU*EREQU
	BQUAD=ERPOL*ERPOL
	UMR=ATAN(1.0)*4./180.
	RETURN
	END


       subroutine dgrf45_dat_dat(NMAX,ERAD,GH,IER)
 
c DESCRIPTION 
c This subroutine contains DATA in the form of a data
c statement. For the file dgrf45.dat  .
c 
c added by RY 11 Nov 1994;
c later may be stored as FITS files
c 
        integer NMAX,IER,i,max
        real ERAD,GH(*),GHD(120)

        NMAX =   10
        ERAD =     6371.20

      DATA GHD /-30594.00,-2285.00,5810.00,-1244.00,  2990.00,
     & -1702.00,1578.00,477.00,1282.00,-1834.00,-499.00,1255.00,
     & 186.00,913.00,-11.00,944.00,776.00,144.00,544.00,-276.00,
     & -421.00, -55.00,304.00,-178.00,-253.00,346.00,-12.00,194.00,
     & 95.00,-20.00,-67.00,-142.00,-119.00,-82.00,82.00,59.00,
     & 57.00,     6.00,     6.00,   100.00,  -246.00,    16.00,
     & -25.00,    -9.00,    21.00,   -16.00,  -104.00,   -39.00,
     &  70.00,   -40.00,   -45.00,     0.00,   -18.00,     0.00,
     &  2.00,   -29.00,     6.00,   -10.00,    28.00,    15.00,
     & -17.00,    29.00,   -22.00,    13.00,     7.00,    12.00,
     & -8.00,   -21.00,    -5.00,   -12.00,     9.00,    -7.00,
     &  7.00,     2.00,   -10.00,    18.00,     7.00,     3.00,
     &  2.00,   -11.00,     5.00,   -21.00,   -27.00,     1.00,
     &  17.00,   -11.00,    29.00,     3.00,    -9.00,    16.00,
     &  4.00,    -3.00,     9.00,    -4.00,     6.00,    -3.00,
     &  1.00,    -4.00,     8.00,    -3.00,    11.00,     5.00,
     &  1.00,     1.00,     2.00,   -20.00,    -5.00,    -1.00,
     & -1.00,    -6.00,     8.00,     6.00,    -1.00,    -4.00,
     & -3.00,    -2.00,     5.00,     0.00,    -2.00,    -2.00/

        max = 120
        do i=1,max
         GH(i) = GHD(i)
        enddo
 
        return
        end


       subroutine dgrf55_dat_dat(NMAX,ERAD,GH,IER)
 
c DESCRIPTION 
c This subroutine contains DATA in the form of a data
c statement. For the file dgrf55.dat  .
c 
c added by RY 11 Nov 1994;
c later may be stored as FITS files
c 
        integer nmax,max,ier,i
        real gh(*),GHD(120),erad

        nmax =   10
        erad =     6371.20

      DATA GHD /-30500.00,-2215.00,5820.00,-1440.00,3003.00,-1898.00,
     & 1581.00,291.00,1302.00,-1944.00,-462.00,1288.00,216.00,882.00,
     & -83.00,958.00,796.00,133.00,510.00,-274.00,-397.00,-23.00,
     & 290.00, -230.00,-229.00,360.00,15.00,230.00,110.00,-23.00,
     & -98.00, -152.00, -121.00, -69.00,78.00,47.00,57.00,-9.00,3.00,
     & 96.00,-247.00,48.00, -8.00,-16.00,7.00,-12.00,-107.00,-24.00,
     & 65.00, -56.00, -50.00,2.00, -24.00,10.00,-4.00,-32.00,8.00,
     & -11.00, 28.00, 9.00, -20.00,18.00,-18.00,11.00,9.00,10.00,
     & -6.00,   -15.00,   -14.00,     5.00,     6.00,   -23.00,
     & 10.00,     3.00,    -7.00,    23.00,     6.00,    -4.00,
     & 9.00,   -13.00,     4.00,     9.00,   -11.00,    -4.00,
     & 12.00,    -5.00,     7.00,     2.00,     6.00,     4.00,
     & -2.00,     1.00,    10.00,     2.00,     7.00,     2.00,
     & -6.00,     5.00,     5.00,    -3.00,    -5.00,    -4.00,
     & -1.00,     0.00,     2.00,    -8.00,    -3.00,    -2.00,
     &  7.00,    -4.00,     4.00,     1.00,    -2.00,    -3.00,
     &  6.00,     7.00,    -2.00,    -1.00,     0.00,    -3.00/

        max = 120
        do i=1,max
         GH(i) = GHD(i)
        enddo
 
        return
        end
       subroutine dgrf60_dat_dat(NMAX,ERAD,GH,IER)
 
c DESCRIPTION 
c This subroutine contains DATA in the form of a data
c statement. For the file dgrf60.dat  .
c 
c added by RY 11 Nov 1994;
c later may be stored as FITS files
c 
        integer NMAX,IER,i,max
        real GH(*),ERAD,GHD(120)

        nmax =   10
        erad =     6371.20

      DATA GHD /-30421.00,-2169.00,5791.00,-1555.00,3002.00,
     & -1967.00,1590.00,206.00,1302.00,-1992.00,-414.00,1289.00,
     &  224.00,878.00,-130.00,957.00,800.00,135.00,504.00,-278.00,
     &  -394.00,3.00,269.00,-255.00,-222.00,362.00,16.00,242.00,125.00,
     & -26.00,-117.00,-156.00, -114.00, -63.00, 81.00, 46.00,
     &    58.00,   -10.00,     1.00,    99.00,  -237.00,    60.00,
     &    -1.00,   -20.00,    -2.00,   -11.00,  -113.00,   -17.00,
     &    67.00,   -56.00,   -55.00,     5.00,   -28.00,    15.00,
     &    -6.00,   -32.00,     7.00,    -7.00,    23.00,    17.00,
     &   -18.00,     8.00,   -17.00,    15.00,     6.00,    11.00,
     &    -4.00,   -14.00,   -11.00,     7.00,     2.00,   -18.00,
     &    10.00,     4.00,    -5.00,    23.00,    10.00,     1.00,
     &     8.00,   -20.00,     4.00,     6.00,   -18.00,     0.00,
     &    12.00,    -9.00,     2.00,     1.00,     0.00,     4.00,
     &    -3.00,    -1.00,     9.00,    -2.00,     8.00,     3.00,
     &     0.00,    -1.00,     5.00,     1.00,    -3.00,     4.00,
     &     4.00,     1.00,     0.00,     0.00,    -1.00,     2.00,
     &     4.00,    -5.00,     6.00,     1.00,     1.00,    -1.00,
     &    -1.00,     6.00,     2.00,     0.00,     0.00,    -7.00/

        max = 120
        do i=1,max
          GH(i) = GHD(i)
        enddo 
        return
        end


       subroutine dgrf50_dat_dat(NMAX,ERAD,GH,IER)
 
c DESCRIPTION 
c This subroutine contains DATA in the form of a data
c statement. For the file dgrf50.dat  .
c 
c added by RY 11 Nov 1994;
c later may be stored as FITS files
c 
        integer NMAX,IER,i,max
        real ERAD, GH(*),GHD(120)

        nmax =   10
        erad =     6371.20

      DATA GHD /-30554.00,-2250.00,5815.00,-1341.00,2998.00,-1810.00,
     & 1576.00,381.00,1297.00,-1889.00,-476.00,1274.00,206.00,896.00,
     & -46.00,954.00,792.00,136.00,528.00,-278.00,-408.00, -37.00,
     & 303.00,-210.00,-240.00,349.00,3.00,211.00,103.00,-20.00,-87.00,
     & -147.00,-122.00,-76.00,80.00,54.00,57.00,-1.00,4.00,99.00,
     & -247.00,33.00,-16.00,-12.00,12.00,-12.00,-105.00,-30.00,65.00,
     & -55.00,-35.00,2.00,-17.00,1.00,0.00,-40.00,10.00,-7.00,36.00,
     & 5.00,  -18.00, 19.00,   -16.00,    22.00,    15.00,     5.00,
     &    -4.00,   -22.00,    -1.00,     0.00,    11.00,   -21.00,
     &    15.00,    -8.00,   -13.00,    17.00,     5.00,    -4.00,
     &    -1.00,   -17.00,     3.00,    -7.00,   -24.00,    -1.00,
     &    19.00,   -25.00,    12.00,    10.00,     2.00,     5.00,
     &     2.00,    -5.00,     8.00,    -2.00,     8.00,     3.00,
     &   -11.00,     8.00,    -7.00,    -8.00,     4.00,    13.00,
     &    -1.00,    -2.00,    13.00,   -10.00,    -4.00,     2.00,
     &     4.00,    -3.00,    12.00,     6.00,     3.00,    -3.00,
     &     2.00,     6.00,    10.00,    11.00,     3.00,     8.00/

        max = 120
        do i=1,max
          GH(i) = GHD(i)
        enddo 
        return
        end


       subroutine dgrf70_dat_dat(NMAX,ERAD,GH,IER)
 
c DESCRIPTION 
c This subroutine contains DATA in the form of a data
c statement. For the file dgrf70.dat  .
c 
c added by RY 11 Nov 1994;
c later may be stored as FITS files
c
        integer NMAX,i,max,ier
        real ERAD,GH(*),GHD(120)
 
        NMAX =   10
        ERAD =     6371.20

      DATA GHD /-30220.00,-2068.00,5737.00,-1781.00,3000.00,-2047.00,
     & 1611.00, 25.00,1287.00,-2091.00, -366.00,1278.00,251.00,838.00,
     & -196.00, 952.00, 800.00,167.00,461.00,-266.00,-395.00, 26.00,
     & 234.00,-279.00,-216.00, 359.00, 26.00, 262.00, 139.00,-42.00,
     &  -139.00,  -160.00,   -91.00,   -56.00,    83.00,    43.00,
     &    64.00,   -12.00,    15.00,   100.00,  -212.00,    72.00,
     &     2.00,   -37.00,     3.00,    -6.00,  -112.00,     1.00,
     &    72.00,   -57.00,   -70.00,     1.00,   -27.00,    14.00,
     &    -4.00,   -22.00,     8.00,    -2.00,    23.00,    13.00,
     &   -23.00,    -2.00,   -11.00,    14.00,     6.00,     7.00,
     &    -2.00,   -15.00,   -13.00,     6.00,    -3.00,   -17.00,
     &     5.00,     6.00,     0.00,    21.00,    11.00,    -6.00,
     &     3.00,   -16.00,     8.00,    10.00,   -21.00,     2.00,
     &    16.00,   -12.00,     6.00,    10.00,    -4.00,    -1.00,
     &    -5.00,     0.00,    10.00,     3.00,    11.00,     1.00,
     &    -2.00,    -1.00,     1.00,    -3.00,    -3.00,     1.00,
     &     2.00,     1.00,    -5.00,     3.00,    -1.00,     4.00,
     &     6.00,    -4.00,     4.00,     0.00,     1.00,    -1.00,
     &     0.00,     3.00,     3.00,     1.00,    -1.00,    -4.00/

        max = 120
        do i=1,max
          GH(i) = GHD(i)
        enddo 
        return
        end


       subroutine dgrf75_dat_dat(NMAX,ERAD,GH,IER)
 
c DESCRIPTION 
c This subroutine contains DATA in the form of a data
c statement. For the file dgrf75.dat  .
c 
c added by RY 11 Nov 1994;
c later may be stored as FITS files
c 
        integer NMAX,IER,max,i
        real ERAD,GH(*),GHD(120)

        NMAX =   10
        ERAD =     6371.20

      DATA GHD /-30100.00,-2013.00,5675.00,-1902.00,3010.00,-2067.00,
     & 1632.00,-68.00,1276.00,-2144.00,-333.00,1260.00,262.00,830.00,
     & -223.00,946.00,791.00,191.00,438.00,-265.00,-405.00,39.00,216.00,
     & -288.00,-218.00,356.00,31.00,264.00,148.00,-59.00,
     &  -152.00,  -159.00,   -83.00,   -49.00,    88.00,    45.00,
     &    66.00,   -13.00,    28.00,    99.00,  -198.00,    75.00,
     &     1.00,   -41.00,     6.00,    -4.00,  -111.00,    11.00,
     &    71.00,   -56.00,   -77.00,     1.00,   -26.00,    16.00,
     &    -5.00,   -14.00,    10.00,     0.00,    22.00,    12.00,
     &   -23.00,    -5.00,   -12.00,    14.00,     6.00,     6.00,
     &    -1.00,   -16.00,   -12.00,     4.00,    -8.00,   -19.00,
     &     4.00,     6.00,     0.00,    18.00,    10.00,   -10.00,
     &     1.00,   -17.00,     7.00,    10.00,   -21.00,     2.00,
     &    16.00,   -12.00,     7.00,    10.00,    -4.00,    -1.00,
     &    -5.00,    -1.00,    10.00,     4.00,    11.00,     1.00,
     &    -3.00,    -2.00,     1.00,    -3.00,    -3.00,     1.00,
     &     2.00,     1.00,    -5.00,     3.00,    -2.00,     4.00,
     &     5.00,    -4.00,     4.00,    -1.00,     1.00,    -1.00,
     &     0.00,     3.00,     3.00,     1.00,    -1.00,    -5.00/

        max = 120
        do i=1,max
         GH(i) = GHD(i)
        enddo 
        return
        end


       subroutine dgrf65_dat_dat(NMAX,ERAD,GH,IER)
 
c DESCRIPTION 
c This subroutine contains DATA in the form of a data
c statement. For the file dgrf65.dat  .
c 
c added by RY 11 Nov 1994;
c later may be stored as FITS files
c 
        integer NMAX,IER,i,max
        real ERAD,GH(*),GHD(120)

        NMAX =   10
        ERAD =     6371.20

      DATA GHD /-30334.00,-2119.00, 5776.00,-1662.00,2997.00,-2016.00,
     & 1594.00,114.00,1297.00,-2038.00, -404.00,1292.00,240.00,856.00,
     & -165.00,   957.00,   804.00,   148.00,
     &   479.00,  -269.00,  -390.00,    13.00,   252.00,  -269.00,
     &  -219.00,   358.00,    19.00,   254.00,   128.00,   -31.00,
     &  -126.00,  -157.00,   -97.00,   -62.00,    81.00,    45.00,
     &    61.00,   -11.00,     8.00,   100.00,  -228.00,    68.00,
     &     4.00,   -32.00,     1.00,    -8.00,  -111.00,    -7.00,
     &    75.00,   -57.00,   -61.00,     4.00,   -27.00,    13.00,
     &    -2.00,   -26.00,     6.00,    -6.00,    26.00,    13.00,
     &   -23.00,     1.00,   -12.00,    13.00,     5.00,     7.00,
     &    -4.00,   -12.00,   -14.00,     9.00,     0.00,   -16.00,
     &     8.00,     4.00,    -1.00,    24.00,    11.00,    -3.00,
     &     4.00,   -17.00,     8.00,    10.00,   -22.00,     2.00,
     &    15.00,   -13.00,     7.00,    10.00,    -4.00,    -1.00,
     &    -5.00,    -1.00,    10.00,     5.00,    10.00,     1.00,
     &    -4.00,    -2.00,     1.00,    -2.00,    -3.00,     2.00,
     &     2.00,     1.00,    -5.00,     2.00,    -2.00,     6.00,
     &     4.00,    -4.00,     4.00,     0.00,     0.00,    -2.00,
     &     2.00,     3.00,     2.00,     0.00,     0.00,    -6.00/

        max = 120
        do i=1,max
         GH(i) = GHD(i)
        enddo 
        return
        end


       subroutine dgrf85_dat_dat(NMAX,ERAD,GH,IER)
 
c DESCRIPTION 
c This subroutine contains DATA in the form of a data
c statement. For the file dgrf85.dat  .
c 
c added by RY 11 Nov 1994;
c later may be stored as FITS files
c 
        integer NMAX,IER,max,i
        real ERAD,GH(*),GHD(120)

        nmax =   10
        erad =     6371.20

      DATA GHD /-29873.00,-1905.00,5500.00,-2072.00,3044.00,-2197.00,
     & 1687.00,-306.00,1296.00,-2208.00,-310.00,1247.00,284.00,829.00,
     & -297.00,936.00,780.00,232.00,361.00,-249.00,-424.00,69.00,
     & 170.00,-297.00,-214.00,355.00,47.00,253.00,150.00,-93.00,
     & -154.00,-164.00,-75.00,-46.00,95.00,53.00,
     & 65.00,-16.00,51.00,88.00,-185.00,69.00,
     & 4.00,-48.00,16.00,-1.00,-102.00,21.00,
     &    74.00,   -62.00,   -83.00,     3.00,   -27.00,    24.00,
     &    -2.00,    -6.00,    20.00,     4.00,    17.00,    10.00,
     &   -23.00,     0.00,    -7.00,    21.00,     6.00,     8.00,
     &     0.00,   -19.00,   -11.00,     5.00,    -9.00,   -23.00,
     &     4.00,    11.00,     4.00,    14.00,     4.00,   -15.00,
     &    -4.00,   -11.00,     5.00,    10.00,   -21.00,     1.00,
     &    15.00,   -12.00,     9.00,     9.00,    -6.00,    -3.00,
     &    -6.00,    -1.00,     9.00,     7.00,     9.00,     1.00,
     &    -7.00,    -5.00,     2.00,    -4.00,    -4.00,     1.00,
     &     3.00,     0.00,    -5.00,     3.00,    -2.00,     6.00,
     &     5.00,    -4.00,     3.00,     0.00,     1.00,    -1.00,
     &     2.00,     4.00,     3.00,     0.00,     0.00,    -6.00/

        max = 120
        do i=1,max
          GH(i) = GHD(i)
        enddo 
        return
        end


       subroutine dgrf80_dat_dat(NMAX,ERAD,GH,IER)
 
c DESCRIPTION 
c This subroutine contains DATA in the form of a data
c statement. For the file dgrf80.dat  .
c 
c added by RY 11 Nov 1994;
c later may be stored as FITS files
c
         integer NMAX,IER,max,i
         real ERAD,GH(*),GHD(120)
 
        NMAX =   10
        ERAD =     6371.20

      DATA GHD /-29992.00,-1956.00,5604.00,-1997.00,3027.00,-2129.00,
     & 1663.00,-200.00,1281.00,-2180.00,-336.00,1251.00,271.00,833.00,
     & -252.00,938.00,782.00,212.00,398.00,-257.00,-419.00,53.00,
     & 199.00,-297.00,-218.00,357.00,46.00,261.00,150.00,-74.00,
     & -151.00,-162.00,-78.00,-48.00,92.00,48.00,66.00,-15.00,42.00,
     & 93.00,-192.00,71.00,4.00, -43.00, 14.00, -2.00,-108.00, 17.00,
     &    72.00,   -59.00,   -82.00,     2.00,   -27.00,    21.00,
     &    -5.00,   -12.00,    16.00,     1.00,    18.00,    11.00,
     &   -23.00,    -2.00,   -10.00,    18.00,     6.00,     7.00,
     &     0.00,   -18.00,   -11.00,     4.00,    -7.00,   -22.00,
     &     4.00,     9.00,     3.00,    16.00,     6.00,   -13.00,
     &    -1.00,   -15.00,     5.00,    10.00,   -21.00,     1.00,
     &    16.00,   -12.00,     9.00,     9.00,    -5.00,    -3.00,
     &    -6.00,    -1.00,     9.00,     7.00,    10.00,     2.00,
     &    -6.00,    -5.00,     2.00,    -4.00,    -4.00,     1.00,
     &     2.00,     0.00,    -5.00,     3.00,    -2.00,     6.00,
     &     5.00,    -4.00,     3.00,     0.00,     1.00,    -1.00,
     &     2.00,     4.00,     3.00,     0.00,     0.00,    -6.00/

        max = 120
        do i=1,max
         GH(i) = GHD(i)
        enddo 
        return
        end


       subroutine igrf90s_dat_dat(NMAX,ERAD,GH,IER)
 
c DESCRIPTION 
c This subroutine contains DATA in the form of a data
c statement. For the file igrf90s.dat .
c 
c added by RY 11 Nov 1994;
c later may be stored as FITS files
c 
        integer NMAX,IER,max,i
        real ERAD,GH(*),GHD(80)

        NMAX =   8
        ERAD =     6371.20

      DATA GHD /18.00,11.00,-16.00,-13.00,2.00,-16.00,0.00,
     & -14.00,3.00,-7.00,4.00,0.00,2.00,-6.00,-11.00,1.00,1.00,
     & 3.00,-7.00,2.00,1.00,3.00,-6.00,-1.00,1.00,0.00,0.00,-2.00,
     & 1.00,-3.00,0.00,0.00, 2.00, 2.00, 0.00,     1.00,
     &     0.00,     0.00,     2.00,    -1.00,     1.00,     0.00,
     &     0.00,    -1.00,     0.00,     1.00,     1.00,     1.00,
     &     1.00,    -1.00,     1.00,     0.00,     0.00,     1.00,
     &     1.00,     2.00,    -1.00,     0.00,     0.00,     0.00,
     &     0.00,     0.00,     0.00,     0.00,    -1.00,     1.00,
     &     0.00,     0.00,     0.00,     0.00,    -1.00,     0.00,
     &     0.00,     0.00,     0.00,    -1.00,    -1.00,     0.00,
     &    -1.00,     1.00/

        max = 80
        do i=1,max
         GH(i) = GHD(i)
        enddo 
        return
        end


       subroutine igrf90_dat_dat(NMAX,ERAD,GH,IER)
 
c DESCRIPTION 
c This subroutine contains DATA in the form of a data
c statement. For the file igrf90.dat  .
c 
c added by RY 11 Nov 1994;
c later may be stored as FITS files
c
         integer NMAX,IER,max,i
         real GH(*),ERAD,GHD(120)
 
        NMAX =   10
        ERAD =     6371.20

      DATA GHD/-29775.00,-1851.00,5411.00,-2136.00,3058.00,-2278.00,
     & 1693.00,-380.00,1315.00,-2240.00, -287.00,1246.00,293.00,
     & 807.00, -348.00,939.00,782.00,248.00,324.00,-240.00,-423.00,
     & 87.00,142.00,-299.00,-211.00,353.00,47.00,244.00,153.00,
     & -111.00,-154.00,-166.00, -69.00, -37.00, 98.00, 61.00,
     &    64.00,   -16.00,    60.00,    83.00,  -178.00,    68.00,
     &     2.00,   -52.00,    17.00,     2.00,   -96.00,    27.00,
     &    77.00,   -64.00,   -81.00,     4.00,   -27.00,    28.00,
     &     1.00,     1.00,    20.00,     6.00,    16.00,    10.00,
     &   -23.00,     0.00,    -5.00,    22.00,     5.00,    10.00,
     &    -1.00,   -20.00,   -11.00,     7.00,   -12.00,   -22.00,
     &     4.00,    12.00,     4.00,    11.00,     3.00,   -16.00,
     &    -6.00,   -11.00,     4.00,    10.00,   -21.00,     1.00,
     &    15.00,   -12.00,    10.00,     9.00,    -6.00,    -4.00,
     &    -6.00,    -1.00,     9.00,     7.00,     9.00,     2.00,
     &    -7.00,    -6.00,     2.00,    -4.00,    -4.00,     1.00,
     &     2.00,     0.00,    -5.00,     3.00,    -2.00,     6.00,
     &     4.00,    -4.00,     3.00,     0.00,     1.00,    -1.00,
     &     2.00,     4.00,     3.00,     0.00,     0.00,    -6.00/
        max = 120
        do i=1,max
          GH(i) = GHD(i)
        enddo 
        return
        end
