IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
C      PARAMETER (INMAX=75,JNMIN=5,JNMAX=75,KNMIN=2,KNMAX=35)
      PARAMETER (INMAX=75,JNMIN=5,JNMAX=75,KNMIN=2,KNMAX=34)
      PARAMETER (NFMAX=1258)
      COMMON /QPARAM/RAD,DEG,RE,RREF,RBTM,RREFER,RBTMER,STEP0,DTH0,JMIN,
     &JMAX,DTHMSH,KMIN,KMAX,IHEM,IYEAR,IMNTH,IDATE
      COMMON /QMESH /PH0(INMAX),TH0(JNMIN:JNMAX),THMSH(KNMIN:KNMAX),JJJ(
     &NFMAX),KKK(NFMAX),ID(NFMAX*2+1),L(NFMAX*2+1,3),THREF(NFMAX),PHREF(
     &NFMAX),BBBTM(NFMAX),indmsh
      COMMON /QIGRF /RR(NFMAX*2+1),TH(NFMAX*2+1),PH(NFMAX*2+1),YR(NFMAX*
     &2+1),YT(NFMAX*2+1),YP(NFMAX*2+1),BB(NFMAX*2+1),STEP(NFMAX*2+1)

      CALL INIT0(IN,NF)
      
      open(indmsh,file='MESH',form='unformatted')
      DO 10 II=1,IN
         write(*,*) II,'/',IN,' done'
         CALL INIT1(II,NF)
         IF(II.EQ.1) CALL OUT0(IN,NF)
         JF2=NF*2
         CALL INIT2(NF)
         CALL NEXT(JF2)

 1       CONTINUE

         DO 1010 JF=1,JF2
            IF(L(JF,3).EQ.0) THEN
               IF(L(JF,1).EQ.0) THEN
                  IF(ABS(RR(JF)-RREF).LE.RREFER) THEN
                     IiF=(JF+1)/2
                     THREF(IiF)=TH(JF)*DEG
                     PHREF(IiF)=PH(JF)*DEG
                     write(30,*) jf,rr(jf),TH0(JJJ(iiF)),thmsh(kkk(iif))
     &                    ,th(jf)*deg,ph(jf)*deg
                     L(JF,1)=1
                     ID(JF)=MOD(JF,2)*2-1
                     STEP(JF)=STEP0*ID(JF)
                  ELSE IF(((RR(JF).LT.RREF).AND.(RR(JF)-YR(JF).GT.RREF))
     &                    .OR.((RR(JF).GT.RREF).AND.
     &                    (RR(JF)-YR(JF).LT.RREF))) THEN
                     ID(JF)=-ID(JF)
                     STEP(JF)=STEP(JF)*(-0.5)
                  END IF
               END IF

               IF(L(JF,2).EQ.0) THEN
                  IF(ABS(RR(JF)-RBTM).LE.RBTMER) THEN
                     IiF=(JF+1)/2
                     BBBTM(IiF)=BB(JF)
                     L(JF,2)=1
                     ID(JF)=MOD(JF,2)*2-1
                     STEP(JF)=STEP0*ID(JF)
                  ELSE IF(((RR(JF).LT.RBTM).AND.(RR(JF)-YR(JF).GT.RBTM))
     &                    .OR.((RR(JF).GT.RBTM).AND.
     &                    (RR(JF)-YR(JF).LT.RBTM))) THEN
                     ID(JF)=-ID(JF)
                     STEP(JF)=STEP(JF)*(-0.5)
                  END IF
               END IF
            END IF

            IF((RR(JF).LT.RE).AND.(YR(JF).LT.0.0D0)) THEN
               STEP(JF)=0.0D0
               L(JF,3)=1
            END IF
 1010    CONTINUE

         DO 1020 JF=JF2,1,-1
            IF(L(JF,3).EQ.0) THEN
               JF2=JF
               CALL NEXT(JF2)
               GOTO 1
            END IF
 1020    CONTINUE

         CALL OUT1(II,NF)
 10   CONTINUE
      close(indmsh)

      STOP
      END



      SUBROUTINE IGRF90(IYEAR,IMNTH,IDATE)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      COMMON /QIGCOE/GG(1:8,0:8),HH(1:8,0:8),SS(0:8,0:8)
      DIMENSION G(7,1:8,0:8),H(7,1:8,0:8)

      IF((IYEAR.LT.1965).OR.(IYEAR.GT.1995)) GOTO 999
      DATE=IYEAR+(30.4176D0*(IMNTH-1)+(IDATE-1))/365.0D0

      IF(IYEAR.GE.1990) THEN
        L=6
        P=1.0D0
        Q=DATE-1990.0D0
      ELSE IF(IYEAR.GE.1985) THEN
        L=5
        Q=(DATE-1985.0D0)/5.0D0
        P=1.0D0-Q
      ELSE IF(IYEAR.GE.1980) THEN
        L=4
        Q=(DATE-1980.0D0)/5.0D0
        P=1.0D0-Q
      ELSE IF(IYEAR.GE.1975) THEN
        L=3
        Q=(DATE-1975.0D0)/5.0D0
        P=1.0D0-Q
      ELSE IF(IYEAR.GE.1970) THEN
        L=2
        Q=(DATE-1970.0D0)/5.0D0
        P=1.0D0-Q
      ELSE IF(IYEAR.GE.1965) THEN
        L=1
        Q=(DATE-1965.0D0)/5.0D0
        P=1.0D0-Q
      END IF

      DO 20 N=1,8
        DO 30 M=0,N
          GG(N,M)=(G(L,N,M)*P+G(L+1,N,M)*Q)*1.0D-9
          HH(N,M)=(H(L,N,M)*P+H(L+1,N,M)*Q)*1.0D-9
   30 CONTINUE
   20 CONTINUE

      SS(0,0)=1.0D0
      DO 40 N=1,8
        SS(N,0)=SS(N-1,0)*(2*N-1)/DBLE(N)
        SS(N,1)=SS(N,0)*SQRT(2*N/DBLE(N+1))
        DO 50 M=2,N
          SS(N,M)=SS(N,M-1)*SQRT((N-M+1)/DBLE(N+M))
   50 CONTINUE
   40 CONTINUE

  999 RETURN

      DATA (G(1,1,M),M=0,1)/-30334.,-2119./(G(1,2,M),M=0,2)/-1662.,2997.
     &,1594./(G(1,3,M),M=0,3)/1297.,-2038.,1292.,856./(G(1,4,M),M=0,4)/9
     &57.,804.,479.,-390.,252./(G(1,5,M),M=0,5)/-219.,358.,254.,-31.,-15
     &7.,-62./(G(1,6,M),M=0,6)/45.,61.,8.,-228.,4.,1.,-111./(G(1,7,M),M=
     &0,7)/75.,-57.,4.,13.,-26.,-6.,13.,1./(G(1,8,M),M=0,8)/13.,5.,-4.,-
     &14.,0.,8.,-1.,11.,4./
      DATA (H(1,1,M),M=0,1)/0.,5776./(H(1,2,M),M=0,2)/0.,-2016.,114./(H(
     &1,3,M),M=0,3)/0.,-404.,240.,-165./(H(1,4,M),M=0,4)/0.,148.,-269.,1
     &3.,-269./(H(1,5,M),M=0,5)/0.,19.,128.,-126.,-97.,81./(H(1,6,M),M=0
     &,6)/0.,-11.,100.,68.,-32.,-8.,-7./(H(1,7,M),M=0,7)/0.,-61.,-27.,-2
     &.,6.,26.,-23.,-12./(H(1,8,M),M=0,8)/0.,7.,-12.,9.,-16.,4.,24.,-3.,
     &-17./
      DATA (G(2,1,M),M=0,1)/-30220.,-2068./(G(2,2,M),M=0,2)/-1781.,3000.
     &,1611./(G(2,3,M),M=0,3)/1287.,-2091.,1278.,838./(G(2,4,M),M=0,4)/9
     &52.,800.,461.,-395.,234./(G(2,5,M),M=0,5)/-216.,359.,262.,-42.,-16
     &0.,-56./(G(2,6,M),M=0,6)/43.,64.,15.,-212.,2.,3.,-112./(G(2,7,M),M
     &=0,7)/72.,-57.,1.,14.,-22.,-2.,13.,-2./(G(2,8,M),M=0,8)/14.,6.,-2.
     &,-13.,-3.,5.,0.,11.,3./
      DATA (H(2,1,M),M=0,1)/0.,5737./(H(2,2,M),M=0,2)/0.,-2047.,25./(H(2
     &,3,M),M=0,3)/0.,-366.,251.,-196./(H(2,4,M),M=0,4)/0.,167.,-266.,26
     &.,-279./(H(2,5,M),M=0,5)/0.,26.,139.,-139.,-91.,83./(H(2,6,M),M=0,
     &6)/0.,-12.,100.,72.,-37.,-6.,1./(H(2,7,M),M=0,7)/0.,-70.,-27.,-4.,
     &8.,23.,-23.,-11./(H(2,8,M),M=0,8)/0.,7.,-15.,6.,-17.,6.,21.,-6.,-1
     &6./
      DATA (G(3,1,M),M=0,1)/-30100.,-2013./(G(3,2,M),M=0,2)/-1902.,3010.
     &,1632./(G(3,3,M),M=0,3)/1276.,-2144.,1260.,830./(G(3,4,M),M=0,4)/9
     &46.,791.,438.,-405.,216./(G(3,5,M),M=0,5)/-218.,356.,264.,-59.,-15
     &9.,-49./(G(3,6,M),M=0,6)/45.,66.,28.,-198.,1.,6.,-111./(G(3,7,M),M
     &=0,7)/71.,-56.,1.,16.,-14.,0.,12.,-5./(G(3,8,M),M=0,8)/14.,6.,-1.,
     &-12.,-8.,4.,0.,10.,1./
      DATA (H(3,1,M),M=0,1)/0.,5675./(H(3,2,M),M=0,2)/0.,-2067.,68./(H(3
     &,3,M),M=0,3)/0.,-333.,262.,-223./(H(3,4,M),M=0,4)/0.,191.,-265.,39
     &.,-288./(H(3,5,M),M=0,5)/0.,31.,148.,-152.,-83.,88./(H(3,6,M),M=0,
     &6)/0.,-13.,99.,75.,-41.,-4.,11./(H(3,7,M),M=0,7)/0.,-77.,-26.,-5.,
     &10.,22.,-23.,-12./(H(3,8,M),M=0,8)/0.,6.,-16.,4.,-19.,6.,18.,-10.,
     &-17./
      DATA (G(4,1,M),M=0,1)/-29992.,-1956./(G(4,2,M),M=0,2)/-1997.,3027.
     &,1663./(G(4,3,M),M=0,3)/1281.,-2180.,1251.,833./(G(4,4,M),M=0,4)/9
     &38.,782.,398.,-419.,199./(G(4,5,M),M=0,5)/-218.,357.,261.,-74.,-16
     &2.,-48./(G(4,6,M),M=0,6)/48.,66.,42.,-192.,4.,14.,-108./(G(4,7,M),
     &M=0,7)/72.,-59.,2.,21.,-12.,1.,11.,-2./(G(4,8,M),M=0,8)/18.,6.,0.,
     &-11.,-7.,4.,3.,6.,-1./
      DATA (H(4,1,M),M=0,1)/0.,5604./(H(4,2,M),M=0,2)/0.,-2129.,-200./(H
     &(4,3,M),M=0,3)/0.,-336.,271.,-252./(H(4,4,M),M=0,4)/0.,212.,-257.,
     &53.,-297./(H(4,5,M),M=0,5)/0.,46.,150.,-151.,-78.,92./(H(4,6,M),M=
     &0,6)/0.,-15.,93.,71.,-43.,-2.,17./(H(4,7,M),M=0,7)/0.,-82.,-27.,-5
     &.,16.,18.,-23.,-10./(H(4,8,M),M=0,8)/0.,7.,-18.,4.,-22.,9.,16.,-13
     &.,-15./
      DATA (G(5,1,M),M=0,1)/-29873.,-1905./(G(5,2,M),M=0,2)/-2072.,3044.
     &,1687./(G(5,3,M),M=0,3)/1296.,-2208.,1247.,829./(G(5,4,M),M=0,4)/9
     &36.,780.,361.,-424.,170./(G(5,5,M),M=0,5)/-214.,355.,253.,-93.,-16
     &4.,-46./(G(5,6,M),M=0,6)/53.,65.,51.,-185.,4.,16.,-102./(G(5,7,M),
     &M=0,7)/74.,-62.,3.,24.,-6.,4.,10.,0./(G(5,8,M),M=0,8)/21.,6.,0.,-1
     &1.,-9.,4.,4.,4.,-4./
      DATA (H(5,1,M),M=0,1)/0.,5500./(H(5,2,M),M=0,2)/0.,-2197.,-306./(H
     &(5,3,M),M=0,3)/0.,-310.,284.,-297./(H(5,4,M),M=0,4)/0.,232.,-249.,
     &69.,-297./(H(5,5,M),M=0,5)/0.,47.,150.,-154.,-75.,95./(H(5,6,M),M=
     &0,6)/0.,-16.,88.,69.,-48.,-1.,21./(H(5,7,M),M=0,7)/0.,-83.,-27.,-2
     &.,20.,17.,-23.,-7./(H(5,8,M),M=0,8)/0.,8.,-19.,5.,-23.,11.,14.,-15
     &.,-11./
      DATA (G(6,1,M),M=0,1)/-29775.,-1851./(G(6,2,M),M=0,2)/-2136.,3058.
     &,1693./(G(6,3,M),M=0,3)/1315.,-2240.,1246.,807./(G(6,4,M),M=0,4)/9
     &39.,782.,324.,-423.,142./(G(6,5,M),M=0,5)/-211.,353.,244.,-111.,-1
     &66.,-37./(G(6,6,M),M=0,6)/61.,64.,60.,-178.,2.,17.,-96./(G(6,7,M),
     &M=0,7)/77.,-64.,4.,28.,1.,6.,10.,0./(G(6,8,M),M=0,8)/22.,5.,-1.,-1
     &1.,-12.,4.,4.,3.,-6./
      DATA (H(6,1,M),M=0,1)/0.,5411./(H(6,2,M),M=0,2)/0.,-2278.,-380./(H
     &(6,3,M),M=0,3)/0.,-287.,293.,-348./(H(6,4,M),M=0,4)/0.,248.,-240.,
     &87.,-299./(H(6,5,M),M=0,5)/0.,47.,153.,-154.,-69.,98./(H(6,6,M),M=
     &0,6)/0.,-16.,83.,68.,-52.,2.,27./(H(6,7,M),M=0,7)/0.,-81.,-27.,1.,
     &20.,16.,-23.,-5./(H(6,8,M),M=0,8)/0.,10.,-20.,7.,-22.,12.,11.,-16.
     &,-11./
      DATA (G(7,1,M),M=0,1)/18.0,10.6/(G(7,2,M),M=0,2)/-12.9,2.4,0./(G(7
     &,3,M),M=0,3)/3.3,-6.7,0.1,-5.9/(G(7,4,M),M=0,4)/0.5,0.6,-7.,0.5,-5
     &.5/(G(7,5,M),M=0,5)/0.6,-0.1,-1.6,-3.1,-0.1,2.3/(G(7,6,M),M=0,6)/1
     &.3,-0.2,1.8,1.3,-0.2,0.1,1.2/(G(7,7,M),M=0,7)/0.6,-0.5,-0.3,0.6,1.
     &6,0.2,0.2,0.3/(G(7,8,M),M=0,8)/0.2,-0.7,-0.2,0.1,-1.1,0.,-0.1,-0.5
     &,-0.6/
      DATA (H(7,1,M),M=0,1)/0.,-16.1/(H(7,2,M),M=0,2)/0.,-15.8,-13.8/(H(
     &7,3,M),M=0,3)/0.,4.4,1.6,-10.6/(H(7,4,M),M=0,4)/0.,2.6,1.8,3.1,-1.
     &4/(H(7,5,M),M=0,5)/0.,-0.1,0.5,0.4,1.7,0.4/(H(7,6,M),M=0,6)/0.,0.2
     &,-1.3,0.,-0.9,0.5,1.2/(H(7,7,M),M=0,7)/0.,0.6,0.2,0.8,-0.5,-0.2,0.
     &,0./(H(7,8,M),M=0,8)/0.,0.5,-0.2,0.3,0.3,0.4,-0.5,-0.3,0.6/

      END

      SUBROUTINE INIT0(IN,NF)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER (INMAX=75,JNMIN=5,JNMAX=75,KNMIN=2,KNMAX=34)
      PARAMETER (NFMAX=1258)
      COMMON /QPARAM/RAD,DEG,RE,RREF,RBTM,RREFER,RBTMER,STEP0,DTH0,JMIN,
     &JMAX,DTHMSH,KMIN,KMAX,IHEM,IYEAR,IMNTH,IDATE
      COMMON /QMESH /PH0(INMAX),TH0(JNMIN:JNMAX),THMSH(KNMIN:KNMAX),JJJ(
     &NFMAX),KKK(NFMAX),ID(NFMAX*2+1),L(NFMAX*2+1,3),THREF(NFMAX),PHREF(
     &NFMAX),BBBTM(NFMAX),indmsh

      RAD = ASIN(1.0D0)/90.0D0
      DEG = 1.0D0/RAD
      RE = 6370.0D0
      RREF = RE+1000.0D0
      RBTM = RE+3000.0D0
      RREFER = 0.5D0
      RBTMER = 0.5D0
      STEP0 = 20.0D0

      IN=75
      DO 99 II=1,IN
        PH0(II)=-185.0+5.0*(II-1)
   99 CONTINUE

      DTH0 = 1.0D0
      JMIN = 5.0D0/DTH0
      JMAX = 75.0D0/DTH0
      DO 10 JJ=JMIN,JMAX
        TH0(JJ)=DTH0*JJ
   10 CONTINUE

      DTHMSH = 5.0D0
      KMIN = 10.0D0/DTHMSH
      KMAX = 170.0D0/DTHMSH
      DO 20 KK=KMIN,KMAX
        if (DTHMSH*(KK-1).ge.90.0d0) then
          THMSH(KK)=DTHMSH*(KK-1)
        else
          THMSH(KK)=DTHMSH*KK
        endif
   20 CONTINUE

      NF=0
      DO 30 JJ=JMIN,JMAX
        DO 40 KK=KMIN,KMAX
          IF((THMSH(KK).GE. TH0(JJ)).AND.(THMSH(KK).LE.180.0D0-TH0(JJ)))
     & THEN
            IF((SIN(THMSH(KK)*RAD)/SIN(TH0(JJ)*RAD))**2.LT.12.0d0) THEN
              NF = NF+1
              JJJ(NF) = JJ
              KKK(NF) = KK
            END IF
          END IF
   40 CONTINUE
   30 CONTINUE

      IHEM =-1

      IYEAR = 1990
      IMNTH = 1
      IDATE = 1

      CALL DGTOM0
      CALL IGRF90(IYEAR,IMNTH,IDATE)

      RETURN
      END

      SUBROUTINE INIT1(II,NF)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER (INMAX=75,JNMIN=5,JNMAX=75,KNMIN=2,KNMAX=34)
      PARAMETER (NFMAX=1258)
      COMMON /QPARAM/RAD,DEG,RE,RREF,RBTM,RREFER,RBTMER,STEP0,DTH0,JMIN,
     &JMAX,DTHMSH,KMIN,KMAX,IHEM,IYEAR,IMNTH,IDATE
      COMMON /QMESH /PH0(INMAX),TH0(JNMIN:JNMAX),THMSH(KNMIN:KNMAX),JJJ(
     &NFMAX),KKK(NFMAX),ID(NFMAX*2+1),L(NFMAX*2+1,3),THREF(NFMAX),PHREF(
     &NFMAX),BBBTM(NFMAX),indmsh
      COMMON /QRMAX /RMAX(NFMAX*2+1)
      COMMON /QIGRF /RR(NFMAX*2+1),TH(NFMAX*2+1),PH(NFMAX*2+1),YR(NFMAX*
     &2+1),YT(NFMAX*2+1),YP(NFMAX*2+1),BB(NFMAX*2+1),STEP(NFMAX*2+1)

      DO 1010 IFF=1,NF
       JF1=IFF*2-1
       JF2=IFF*2

       RR(JF1) = RE*(SIN(THMSH(KKK(IFF))*RAD)/SIN(TH0(JJJ(IFF))*RAD))**2
     &
       RR(JF2) = RE*(SIN(THMSH(KKK(IFF))*RAD)/SIN(TH0(JJJ(IFF))*RAD))**2
     &
       TH(JF1) = THMSH(KKK(IFF))*RAD
       TH(JF2) = THMSH(KKK(IFF))*RAD
       PH(JF1) = PH0(II) *RAD
       PH(JF2) = PH0(II) *RAD
       STEP(JF1) = STEP0
       STEP(JF2) =-STEP0

       RMAX(JF1) = RR(JF1)
       RMAX(JF2) = RR(JF2)
 1010 CONTINUE

      RETURN
      END

      SUBROUTINE INIT2(NF)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER (INMAX=75,JNMIN=5,JNMAX=75,KNMIN=2,KNMAX=34)
      PARAMETER (NFMAX=1258)
      COMMON /QPARAM/RAD,DEG,RE,RREF,RBTM,RREFER,RBTMER,STEP0,DTH0,JMIN,
     &JMAX,DTHMSH,KMIN,KMAX,IHEM,IYEAR,IMNTH,IDATE
      COMMON /QMESH /PH0(INMAX),TH0(JNMIN:JNMAX),THMSH(KNMIN:KNMAX),JJJ(
     &NFMAX),KKK(NFMAX),ID(NFMAX*2+1),L(NFMAX*2+1,3),THREF(NFMAX),PHREF(
     &NFMAX),BBBTM(NFMAX),indmsh
      COMMON /QIGRF /RR(NFMAX*2+1),TH(NFMAX*2+1),PH(NFMAX*2+1),YR(NFMAX*
     &2+1),YT(NFMAX*2+1),YP(NFMAX*2+1),BB(NFMAX*2+1),STEP(NFMAX*2+1)

      DO 1010 IFF=1,NF
       JF1=IFF*2-1
       JF2=IFF*2

       ID(JF1)= 1
       IF(TH0(JJJ(IFF)).GT.ASIN(SQRT(RE/RREF))*deg) then
          L(JF1,1)=1
       ELSEIF((RR(JF1).GE.RREF).AND.(THMSH(KKK(IFF)).lt.90.0d0)) then
          L(JF1,1)=0
       ELSEIF((RR(JF1).LT.RREF).AND.(THMSH(KKK(IFF)).gt.90.0d0))then
          L(JF1,1)=0
       ELSEIF(THMSH(KKK(IFF)).eq.90.0d0.AND.THMSH(KKK(IFF-1)).ne.90.0d0)
     & then
          L(JF1,1)=0
       ELSE
          L(JF1,1)=1
       END IF

       IF(TH0(JJJ(IFF)).GT.ASIN(SQRT(RE/RBTM))*deg) then
          L(JF1,2)=1
       ELSEIF((RR(JF1).GE.RBTM).AND.(THMSH(KKK(IFF)).lt.90.0d0)) then
          L(JF1,2)=0
       ELSEIF((RR(JF1).LT.RBTM).AND.(THMSH(KKK(IFF)).gt.90.0d0)) then
          L(JF1,2)=0
       ELSEIF(THMSH(KKK(IFF)).eq.90.0d0.AND.THMSH(KKK(IFF-1)).ne.90.0d0)
     & then
          L(JF1,2)=0
       ELSE
          L(JF1,2)=1
       END IF
       L(JF1,3)=0

       ID(JF2)=-1
       IF(TH0(JJJ(IFF)).GT.ASIN(SQRT(RE/RREF))*deg) then
          L(JF2,1)=1
       ELSEIF((RR(JF2).GE.RREF).AND.(THMSH(KKK(IFF)).gt.90.0d0)) then
          L(JF2,1)=0
       ELSEIF((RR(JF2).LT.RREF).AND.(THMSH(KKK(IFF)).lt.90.0d0)) then
          L(JF2,1)=0
       ELSEIF(THMSH(KKK(IFF)).eq.90.0d0.AND.THMSH(KKK(IFF-1)).eq.90.0d0)
     & then
          L(JF2,1)=0
       ELSE
          L(JF2,1)=1
       END IF

       IF(TH0(JJJ(IFF)).GT.ASIN(SQRT(RE/RBTM))*deg) then
          L(JF2,2)=1
       ELSEIF((RR(JF2).GE.RBTM).AND.(THMSH(KKK(IFF)).gt.90.0d0)) then
          L(JF2,2)=0
       ELSEIF((RR(JF2).LT.RBTM).AND.(THMSH(KKK(IFF)).lt.90.0d0)) then
          L(JF2,2)=0
       ELSEIF(THMSH(KKK(IFF)).eq.90.0d0.AND.THMSH(KKK(IFF-1)).eq.90.0d0)
     & then
          L(JF2,2)=0
       ELSE
          L(JF2,2)=1
       END IF
       L(JF2,3)=0

       THREF(IFF)=999.0D0
       PHREF(IFF)=999.0D0
       BBBTM(IFF)= 0.0D0

       write(40,*) jf1,rr(jf1),th0(JJJ(iff)),thmsh(KKK(iff)),(L(jf1,i),i
     &=1,3)
       write(40,*) jf2,rr(jf2),th0(JJJ(iff)),thmsh(KKK(iff)),(L(jf2,i),i
     &=1,3)
 1010 CONTINUE

      RETURN
      END

      SUBROUTINE OUT0(IN,NF)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER (INMAX=75,JNMIN=5,JNMAX=75,KNMIN=2,KNMAX=34)
      PARAMETER (NFMAX=1258)
      COMMON /QPARAM/RAD,DEG,RE,RREF,RBTM,RREFER,RBTMER,STEP0,DTH0,JMIN,
     &JMAX,DTHMSH,KMIN,KMAX,IHEM,IYEAR,IMNTH,IDATE
      COMMON /QMESH /PH0(INMAX),TH0(JNMIN:JNMAX),THMSH(KNMIN:KNMAX),JJJ(
     &NFMAX),KKK(NFMAX),ID(NFMAX*2+1),L(NFMAX*2+1,3),THREF(NFMAX),PHREF(
     &NFMAX),BBBTM(NFMAX),indmsh
      COMMON /QIGRF /RR(NFMAX*2+1),TH(NFMAX*2+1),PH(NFMAX*2+1),YR(NFMAX*
     &2+1),YT(NFMAX*2+1),YP(NFMAX*2+1),BB(NFMAX*2+1),STEP(NFMAX*2+1)

      DIMENSION RMSH(NFMAX)

      DO 10 IFF=1,NF
        RMSH(IFF)=RR(IFF*2)
   10 CONTINUE

      WRITE(indmsh) RREF,RBTM,RREFER,RBTMER,STEP0,DTH0,JMIN,JMAX,DTHMSH,
     &              KMIN,KMAX,IHEM,IYEAR,IMNTH,IDATE
      WRITE(indmsh) IN,NFMAX,NF
      WRITE(indmsh) JJJ,KKK
      WRITE(indmsh) RMSH

      RETURN
      END

      SUBROUTINE OUT1(II,NF)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER (INMAX=75,JNMIN=5,JNMAX=75,KNMIN=2,KNMAX=34)
      PARAMETER (NFMAX=1258)
      COMMON /QPARAM/RAD,DEG,RE,RREF,RBTM,RREFER,RBTMER,STEP0,DTH0,JMIN,
     &JMAX,DTHMSH,KMIN,KMAX,IHEM,IYEAR,IMNTH,IDATE
      COMMON /QMESH /PH0(INMAX),TH0(JNMIN:JNMAX),THMSH(KNMIN:KNMAX),JJJ(
     &NFMAX),KKK(NFMAX),ID(NFMAX*2+1),L(NFMAX*2+1,3),THREF(NFMAX),PHREF(
     &NFMAX),BBBTM(NFMAX),indmsh
      COMMON /QRMAX /RMAX(NFMAX*2+1)

      DIMENSION XLVAL(NFMAX)

      DO 10 IFF=1,NF
         XLVAL(IFF)=MAX(RMAX(IFF*2-1),RMAX(IFF*2))/RE
   10 CONTINUE

      IND=11
      WRITE(indmsh) THREF
      WRITE(indmsh) PHREF
      WRITE(indmsh) XLVAL
      WRITE(indmsh) BBBTM
      WRITE(indmsh) PH0(II)

      RETURN
      END

      SUBROUTINE NEXT(JF2)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER (NFMAX=1258)
      COMMON /QRMAX /RMAX(NFMAX*2+1)
      COMMON /QIGRF /RR(NFMAX*2+1),TH(NFMAX*2+1),PH(NFMAX*2+1),YR(NFMAX*
     &2+1),YT(NFMAX*2+1),YP(NFMAX*2+1),BB(NFMAX*2+1),STEP(NFMAX*2+1)
      COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0

      DO 1010 JF=1,JF2
       CALL DMTOG1(TH(JF),PH(JF))
 1010 CONTINUE

      CALL IGRF(JF2)

      DO 1020 JF=1,JF2
       YR(JF) = YR(JF) *STEP(JF)
       YT(JF) = YT(JF)/RR(JF) *STEP(JF)
       YP(JF) = YP(JF)/RR(JF)/SIN(TH(JF))*STEP(JF)
       RR(JF) = RR(JF)+YR(JF)
       TH(JF) = TH(JF)+YT(JF)
       PH(JF) = PH(JF)+YP(JF)
       CALL DGTOM1(TH(JF),PH(JF))
       RMAX(JF)=MAX(RR(JF),RMAX(JF))
 1020 CONTINUE

      RETURN
      END

      SUBROUTINE IGRF(JF2)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)
      PARAMETER (NFMAX=1258)
      PARAMETER (NIGRF=8)
      COMMON /QIGRF /RR(NFMAX*2+1),TH(NFMAX*2+1),PH(NFMAX*2+1),YR(NFMAX*
     &2+1),YT(NFMAX*2+1),YP(NFMAX*2+1),BB(NFMAX*2+1),STEP(NFMAX*2+1)
      COMMON /QIGCOE/GG(1:8,0:8),HH(1:8,0:8),SS(0:8,0:8)

      DIMENSION AARR (NFMAX*2+1,3:NIGRF+2),COSTH (NFMAX*2+1),SINTH(NFMAX
     &*2+1),COTTH(NFMAX*2+1),COSMPH(NFMAX*2+1,0:NIGRF),SINMPH(NFMAX*2+1,
     &0:NIGRF),P (NFMAX*2+1,0:NIGRF,0:NIGRF),DP(NFMAX*2+1,0:NIGRF,0:NIGR
     &F),BR(NFMAX*2+1),BT(NFMAX*2+1),BP(NFMAX*2+1),XX(NFMAX*2+1),YY(NFMA
     &X*2+1),ZZ(NFMAX*2+1)

      DATA AA/6371.2D0/

      DO 1010 JF=1,JF2
         P (JF,0,0) = 1.0D0
         DP(JF,0,0) = 0.0D0
         BR(JF) = 0.0D0
         BT(JF) = 0.0D0
         BP(JF) = 0.0D0

         AARR (JF,3) = (AA/RR(JF))**3
         COSTH(JF) = COS(TH(JF))
         SINTH(JF) = SIN(TH(JF))
         COTTH(JF) = COSTH(JF)/SINTH(JF)
 1010 CONTINUE

      DO 10 N=4,NIGRF+2
       DO 1020 JF=1,JF2
         AARR(JF,N) = (AA/RR(JF))*AARR(JF,N-1)
 1020 CONTINUE
   10 CONTINUE

      DO 20 M=0,NIGRF
       DO 1030 JF=1,JF2
         COSMPH(JF,M) = COS(M*PH(JF))
         SINMPH(JF,M) = SIN(M*PH(JF))
 1030 CONTINUE
   20 CONTINUE

      DO 30 N=1,NIGRF
        DO 1040 JF=1,JF2
         XX(JF) = 0.0D0
         YY(JF) = 0.0D0
         ZZ(JF) = 0.0D0

         P (JF,N,N ) = SINTH(JF)*P (JF,N-1,N-1)
         DP(JF,N,N ) = SINTH(JF)*DP(JF,N-1,N-1)+COSTH(JF)*P (JF,N-1,N-1)
     &
         P (JF,N,N-1) = COSTH(JF)*P (JF,N-1,N-1)
         DP(JF,N,N-1) = COSTH(JF)*DP(JF,N-1,N-1)-SINTH(JF)*P (JF,N-1,N-1
     &)
 1040 CONTINUE

        DO 40 M=N-2,0,-1
          REALK = DBLE((N-1)**2-M**2)/DBLE((2*N-3)*(2*N-1))
          DO 1050 JF=1,JF2
           P (JF,N,M) = COSTH(JF)*P (JF,N-1,M)-REALK *P (JF,N-2,M)
           DP(JF,N,M) = COSTH(JF)*DP(JF,N-1,M)-SINTH(JF)*P (JF,N-1,M)-RE
     &ALK *DP(JF,N-2,M)
 1050 CONTINUE
   40 CONTINUE

        DO 50 M=0,N
          DO 1060 JF=1,JF2
           PP = SS(N,M)*P (JF,N,M)
           DPP = SS(N,M)*DP(JF,N,M)

           GCHS = GG(N,M)*COSMPH(JF,M) + HH(N,M)*SINMPH(JF,M)
           GSHC = GG(N,M)*SINMPH(JF,M) - HH(N,M)*COSMPH(JF,M)

           XX(JF) = XX(JF) + GCHS*PP
           YY(JF) = YY(JF) + GCHS*DPP
           ZZ(JF) = ZZ(JF) + M*GSHC*PP
 1060 CONTINUE
   50 CONTINUE

        DO 1070 JF=1,JF2
         BR(JF) = BR(JF) + (N+1)*AARR(JF,N+2)*XX(JF)
         BT(JF) = BT(JF) + AARR(JF,N+2)*YY(JF)
         BP(JF) = BP(JF) + AARR(JF,N+2)*ZZ(JF)
 1070 CONTINUE
   30 CONTINUE

      DO 1080 JF=1,JF2
       BT(JF) =-BT(JF)
       BP(JF) = BP(JF)/SINTH(JF)
       BB(JF) = SQRT(BR(JF)**2+BT(JF)**2+BP(JF)**2)

       YR(JF) = BR(JF)/BB(JF)
       YT(JF) = BT(JF)/BB(JF)
       YP(JF) = BP(JF)/BB(JF)
 1080 CONTINUE

      RETURN
      END

      SUBROUTINE DGTOM1(TH,PH)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0

      COSTH = COS(TH)
      SINTH = SIN(TH)
      COSPH = COS(PH)
      SINPH = SIN(PH)
      QUI = SINTH*(COSPH*COSP0-SINPH*SINP0)
      TH = ACOS(COSTH*COST0+SINT0*QUI)
      AAA = SINTH*(SINPH*COSP0+COSPH*SINP0)
      BBB = COST0*QUI-SINT0*COSTH
      PH = ATAN2(AAA,BBB)

      RETURN
      END

      SUBROUTINE DMTOG1(TH,PH)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0

      COSTH = COS(TH)
      SINTH = SIN(TH)
      COSPH = COS(PH)
      SINPH = SIN(PH)
      TH = ACOS(COST0*COSTH-SINT0*SINTH*COSPH)
      AAA = -COST0*SINP0*SINTH*COSPH+COSP0*SINTH*SINPH-SINT0*SINP0*COSTH
     &
      BBB = COST0*COSP0*SINTH*COSPH+SINP0*SINTH*SINPH+SINT0*COSP0*COSTH
      PH = ATAN2(AAA,BBB)

      RETURN
      END

      SUBROUTINE DGTOM0
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /QTHPH0/COST0,SINT0,COSP0,SINP0

      RAD = ASIN(1.0D0)/90.0D0
      TH0 = 11.2 D0*RAD
      PH0 = 70.75D0*RAD
      COST0 = COS(TH0)
      SINT0 = SIN(TH0)
      COSP0 = COS(PH0)
      SINP0 = SIN(PH0)

      RETURN
      END