C          Transfer coordinates:
C          ART=0   -   Geographic  -> Geomagnetic
C          ART=1   -   Geomagnetic -> Geographic
C
C             DOLG,DOLM  - GG. and GM. longitude
C             TET ,TETM  - GG. and GM. colatitude
C
      SUBROUTINE GGMRAW(ART,DOLG,TET,DOLM,TETM)
      INTEGER ART
      DOUBLE PRECISION ZPI,FAKTOR,CBG,CI,SI,XLM,BM,CBM,SBM,
     *      CLM,SLM,SBG,BG,SLG,CLG,XLG,YLG
        real DOLG,TET,DOLM,TETM
      ZPI=6.28318530718
      FAKTOR=0.01745329252
      CBG=11.4*FAKTOR
      CI=DCOS(CBG)
      SI=DSIN(CBG)
      IF(ART.EQ.0) GO TO 10
        XLM=DOLM
        BM=90.-TETM
        CBM=DCOS(BM*FAKTOR)
        SBM=DSIN(BM*FAKTOR)
        CLM=DCOS(XLM*FAKTOR)
        SLM=DSIN(XLM*FAKTOR)
        SBG=SBM*CI-CBM*CLM*SI
        BG=DASIN(SBG)
        CBG=DCOS(BG)
        SLG=(CBM*SLM)/CBG
        CLG=(SBM*SI+CBM*CLM*CI)/CBG
        IF(CLG.GT.1..OR.(1.-CLG).LT.1.E-10)CLG=1.
        IF(CLG.LT.-1..OR.(1.+CLG).LT.1.E-10)CLG=-1.
        XLG=DACOS(CLG)
        IF(SLG.LT.0.0) XLG=ZPI-DACOS(CLG)
        BG=BG/FAKTOR
        XLG=XLG/FAKTOR
        XLG=XLG-69.8
        IF(XLG.LT.0.0) XLG=XLG+360.0
        TET=90.-BG
        DOLG=XLG
        GO TO 20
   10   BG=90.-TET
        XLG=DOLG
        YLG=XLG+69.8
        CBG=DCOS(BG*FAKTOR)
        SBG=DSIN(BG*FAKTOR)
        CLG=DCOS(YLG*FAKTOR)
        SLG=DSIN(YLG*FAKTOR)
        SBM=SBG*CI+CBG*CLG*SI
        BM=DASIN(SBM)
        CBM=DCOS(BM)
        SLM=(CBG*SLG)/CBM
        CLM=(-SBG*SI+CBG*CLG*CI)/CBM
        if(abs(clm).gt.1.d0)then
          CLM=sign(1.d0,CLM)
        endif
        XLM=DACOS(CLM)
        IF(SLM.LT.0.0) XLM=ZPI-DACOS(CLM)
        BM=BM/FAKTOR
        XLM=XLM/FAKTOR
        DOLM=XLM
        IF(ABS(TET-180.).LT.1.E-3)DOLM=0.
        TETM=90.-BM
   20 CONTINUE
      RETURN
      END