c $Id: igrf.f 3304 2011-01-17 15:25:59Z brideout $        
C
C Subroutine IGRF_SUB to compute IGRF parameters for IRI. 
C Rewritten by Bill Rideout to simply call other geolib methods.  
C
        subroutine igrf_sub(xlat,xlong,year,height,
     &          xl,icode,dipl,babs)
c-----------------------------------------------------------------------        
c INPUT:
c    xlat      geodatic latitude in degrees
c    xlong     geodatic longitude in degrees
c    year      decimal year (year+month/12.0-0.5 or 
c                  year+day-of-year/365 or ../366 if leap year) 
c    height    height in km
c OUTPUT:
c    xl        L value
c    icode      =1  L is correct; =2  L is not correct;
c               =3  an approximation is used
c    dipl      dip latitude in degrees
c    babs     magnetic field strength in gauss
c-----------------------------------------------------------------------        

      REAL              xlat,xlong,xl,dipl,babs,year,height
      INTEGER           icode
      
      DOUBLE PRECISION  TM,GDLAT,GLON,GDALT,GCLAT,RKM,T,CT,ST,P,CP,SP,
     *                  BR,BT,BP,B,XB,YB,ZB,GCALT,BB,RL,
     *                  BDOWN,BNORTH,BEAST
C     .. External Subroutines ..
      EXTERNAL CONVRT,MILMAG,INVAR
      
C     .. Data statements ..
      DATA DTR/0.0174532925199D0/
      
        
C     convert to doubles 
      TM = DBLE(year)
      GDLAT = DBLE(xlat)
      GLON = DBLE(xlong)
      GDALT = DBLE(height)
      CALL CONVRT(1,GDLAT,GDALT,GCLAT,RKM)
      
C     .....calculate magnetic field at observation point.....
      T = DTR*(90.D0-GCLAT)
      CT = DCOS(T)
      ST = DSIN(T)
      P = DTR*GLON
      CP = DCOS(P)
      SP = DSIN(P)
      CALL MILMAG(TM,RKM,ST,CT,SP,CP,BR,BT,BP,B)
      CALL GDV(GDLAT,GCLAT,BR,BT,BP,XB,YB,ZB)
      BNORTH = XB
      BEAST = YB
      BDOWN = ZB
C
C     .....calculate dip latitude.....
      DIPL=REAL(ATAN(BDOWN/2.0/sqrt(BNORTH*BNORTH+BEAST*BEAST))/DTR)
      
C     .....calculate l-shell parameter.....
      GCALT = RKM - 6378.16D0/DSQRT(1.D0+.0067397D0*ST*ST)
      CALL INVAR(TM,GCLAT,GLON,GCALT,0.01D0,BB,RL)    
      RL = MAX(RL,1.0D0)      
  
C     convert results
      xl=REAL(RL)
      icode = 1
      babs=REAL(B)
      RETURN
      END
C
C

C
      SUBROUTINE GEODIP(IYR,SLA,SLO,DLA,DLO,J)

C  Calculates dipole geomagnetic coordinates from geocentric coordinates
C  or vice versa.

C                     J=0           J=1
C		INPUT:     J,SLA,SLO     J,DLA,DLO
C		OUTPUT:     DLA,DLO       SLA,SLO

C  Last revision: November 2005 (Vladimir Papitashvili)
C  The code is modifed from GEOCOR written by V.Popov and V.Papitashvili
C  in mid-1980s. 
      REAL SLA,SLO,DLA,DLO
      INTEGER IYR,J
      DOUBLE PRECISION R,COL,RLO,X,Y,Z,RM,TH,PF,XM,YM,ZM
      DOUBLE PRECISION SZM,DCO

C     .. Data statements ..
      DATA UMR /0.0174532925199D0/ 

C  Earth's radius (km) RE = 6371.2

C  The radius of the sphere to compute the coordinates (in Re)
C        RH = (RE + HI)/RE
         R = 1.

         if(j.gt.0) goto 1234
        
         COL = DBLE((90.- SLA)*UMR)
         RLO = DBLE(SLO*UMR)
      CALL SPHCAR(R,COL,RLO,X,Y,Z,1)
      CALL GEOMAG(X,Y,Z,XM,YM,ZM,1,IYR)
      CALL SPHCAR(RM,TH,PF,XM,YM,ZM,-1)
         SZM = ZM
         DLO = REAL(PF/UMR)
         DCO = TH/UMR
         DLA = REAL(90.- DCO)
         RETURN
          
1234     continue
      COL = DBLE((90.- DLA)*UMR)
      RLO = DBLE(DLO*UMR)
      CALL SPHCAR(R,COL,RLO,XM,YM,ZM,1)
      CALL GEOMAG(X,Y,Z,XM,YM,ZM,-1,IYR)
      CALL SPHCAR(RM,TH,PF,X,Y,Z,-1)
        SZM = ZM
        SLO = REAL(PF/UMR)
        SCO = TH/UMR
        SLA = REAL(90.- SCO)

      RETURN
      END
c
c
C
      subroutine igrf_dip(xlat,xlong,year,height,dip,dipl,ymodip)
c-----------------------------------------------------------------------        
c INPUT:
c    xlat      geodatic latitude in degrees
c    xlong     geodatic longitude in degrees
c    year      decimal year (year+month/12.0-0.5 or 
c                  year+day-of-year/365 or ../366 if leap year) 
c    height    height in km
c OUTPUT:
c    dip       magnetic inclination (dip) in degrees
c    dipl      dip latitude in degrees
c    ymodip    modified dip latitude = asin{dip/sqrt[dip^2+cos(LATI)]} 
c-----------------------------------------------------------------------        
      REAL              xlat,xlong,dipl,year,height
      
      DOUBLE PRECISION  TM,GDLAT,GLON,GDALT,GCLAT,RKM,T,CT,ST,P,CP,SP,
     *                  BR,BT,BP,B,XB,YB,ZB,HB,AINC,
     *                  dipdiv,SMODIP,BDOWN,BNORTH,BEAST
C     .. External Subroutines ..
      EXTERNAL CONVRT,MILMAG,INVAR
      
C     .. Data statements ..
      DATA DTR/0.0174532925199D0/
      
        
C     convert to doubles to use coord 
      TM = DBLE(year)
      GDLAT = DBLE(xlat)
      GLON = DBLE(xlong)
      GDALT = DBLE(height)
      CALL CONVRT(1,GDLAT,GDALT,GCLAT,RKM)
      
C     .....calculate magnetic field at observation point.....
      T = DTR*(90.D0-GCLAT)
      CT = DCOS(T)
      ST = DSIN(T)
      P = DTR*GLON
      CP = DCOS(P)
      SP = DSIN(P)
      CALL MILMAG(TM,RKM,ST,CT,SP,CP,BR,BT,BP,B)
      CALL GDV(GDLAT,GCLAT,BR,BT,BP,XB,YB,ZB)
      BNORTH = XB
      BEAST = YB
      BDOWN = ZB      
C
C     .....calculate inclination.....
      HB = DSQRT(XB*XB+YB*YB)
      AINC = DATAN2(ZB,HB) 
      
      dipdiv=AINC/SQRT(AINC*AINC+cos(XLAT*DTR))
      IF(ABS(dipdiv).GT.1.D0) dipdiv=SIGN(1.D0,dipdiv)
      SMODIP=ASIN(dipdiv)

      DIPL=REAL(ATAN(BDOWN/2.0/sqrt(BNORTH*BNORTH+BEAST*BEAST))/DTR)
      YMODIP=REAL(SMODIP/DTR)                            
      DIP=REAL(AINC/DTR)
      RETURN
      END