C         A PROGRAM TO CONVERT BETWEEN SOME COORDINATE SYSTEMS
$declare
$large
      INTERFACE TO SUBROUTINE
     1             clearscreen[FAR,C,ALIAS:"__clearscreen"] (area)
      INTEGER*2 area
      END

      INTEGER        KIND,incode,outcode,istat,i,datnum
      REAL*8         X(7000),Y(7000),Z(7000),pi
      CHARACTER*10   COMM1(7000),COMM2(7000),COMM3(7000)
      character*40   inname,outname,tempname
      character*5    junk
      character*40   infmt,outfmt,outfmt1,outfmt2
      logical        Is84in,Is84out,change,FileIn,jump,FileOpen
      PI=3.141592653589793D0
      open(10,file='transfer.cfg',status='old',iostat=istat)
      if(istat .eq. 0)then
        read(10,'(a40)')infmt
        read(10,'(a40)')outfmt1
        read(10,'(a40)')outfmt2
        close(10)
      else
        infmt='(3(f13.0,1x),3a10)'
        outfmt1='(3(f13.8,1x),3a10)'
        outfmt2='(3(f13.3,1x),3a10)'
        write(*,*)' Config File transfer.cfg not found, use default'
        write(*,*)' Input format (3(f13.0,1x),3a10) '
        write(*,*)' Output format for PLH output (3(f13.8,1x),3a10)'
        write(*,*)' Output format for X,Y,N,E output (3(f13.3,1x),3a10)'
        read(*,'(a)')junk
      endif
      FileOpen=.false.
1     CALL clearscreen(0)
      if( FileOpen)then
         write(*,'(//20x,a,a)')'Process File : ',inname
      else
         write(*,'(//20x,a)')'Process File :    None'
      endif
      WRITE(*,1001)
1001   FORMAT(/15X,'*****************************************'/
     *        15X,'*   Coordinate Transformation Porgram   *'/
     *        15X,'*****************************************'//
     *        15X,'******   Coordinate System Table   ******'/
     *        15X,'[1]  WGS(84)(Lon,Lat,H) , [2]  WGS(84)(X,Y,Z)'/
     *        15X,'[3]  TM2(84)(E,N,H) , [4]  UTM(84)(E,N,H)'//
     *        15X,'[5]  GRS(67)(Lon,Lat,H) , [6]  GRS(67)(X,Y,Z)'/
     *        15X,'[7]  TM2(67)(E,N,H) , [8]  UTM(67)(E,N,H)'//)
      Write(*,'(10x,a\)')'Specify input file name([K] for KeyBoard) : '
        read(*,'(a)')tempname
      if( tempname .eq. 'k' .or. tempname .eq. 'K')then
         FileIn=.false.
      else if( tempname .eq. ' ' .and. (.not. FileOpen))then
         write(*,'(/15x,a)')'**** Error input Specified ****'
         read(*,'(a)')junk
      else if( tempname .eq. ' ' .and. FileOpen)then
         goto 12
      else
         inname=tempname
         open(1,file=inname,status='old',iostat=istat)
         if( istat .ne. 0)then
           write(*,'(//15x,a)')'** File Not Found **, Press Enter'
           read(*,'(a)')junk
           goto 1
         endif
         datnum=1
         do i=1,7000
           read(1,fmt=infmt,end=11)x(i),y(i),z(i),comm1(i),
     *           comm2(i),comm3(i)
           datnum=i
         enddo
11       if (datnum .gt. 6999)then
            write(*,'(//15x,a)')'*** Too Much data ***, Reduce to',
     *       ' 6999 lines'
            read(*,'(a)')junk
            stop
         else
            write(*,'(/15x,a,i8)')'Total data point : ',datnum
         endif
         FileIn=.true.
      endif
12    continue
      write(*,'(10x,a\)')'Enter Input Coordinate System kind : '
        read(*,*)incode
        if( incode .le. 4 .and. incode .gt. 0)then
            Is84in=.true.
        else if( incode .ge. 5 .and. incode .le. 8)then
            Is84in=.false.
        else
            write(*,'(/15x,a)')'**** Error **** , Press Enter'
            read(*,'(a)')junk
            goto 1
        endif
      write(*,'(10x,a\)')'Enter Output Coordiante System kind : '
        read(*,*)outcode
        if( outcode .le. 4 .and. outcode .gt. 0)then
            Is84out=.true.
        else if( outcode .le. 8 .and. outcode .gt. 4)then
            Is84out=.false.
        else
            write(*,'(/15x,a)')'**** Error **** , Press Enter'
            read(*,'(a)')junk
            goto 1
        endif
        if( outcode .eq. 1 .or. outcode .eq. 5)then
           outfmt=outfmt1
        else
           outfmt=outfmt2
        endif
      if( FileIn )then
        write(*,'(/10x,a\)')'Enter Output file name : '
        read(*,'(a)')outname
        open(2,file=outname,status='unknown')
        write(*,*)'            Process your data, Please stand by....'
      else
        write(*,'(1X,a\)')'Enter Coordinate : '
        read(*,*)x(1),y(1),z(1)
        datnum=1
      endif
      if(Is84in .and. Is84out .or. (.not.(Is84in .or. Is84out)))then
        change=.false.
      else
        change=.true.
      endif
      do i=1,datnum
        jump=.false.
        if( incode .eq. 1)then
           x(i)=pi*x(i)/180.d0
           y(i)=pi*y(i)/180.d0
           kind=1
        else if(incode .eq. 5)then
           x(i)=pi*x(i)/180.d0
           y(i)=pi*y(i)/180.d0
           kind=0
        else if( incode .eq. 2 )then
           kind=1
           if( .not. change )then
              call xyz2ph(x(i),y(i),z(i),kind)
           else
              jump=.true.
           endif
        else if( incode .eq. 6)then
           kind=0
           if( .not. change )then
              call xyz2ph(x(i),y(i),z(i),kind)
           else
              jump=.true.
           endif
        else if( incode .eq. 3)then
           kind=1
           call tm2ph(x(i),y(i),kind,0)
        else if( incode .eq. 7)then
           kind=0
           call tm2ph(x(i),y(i),kind,0)
        else if( incode .eq. 4)then
           kind=1
           call tm2ph(x(i),y(i),kind,1)
        else if( incode .eq. 8)then
           kind=0
           call tm2ph(x(i),y(i),kind,1)
        endif
        if( change )then
          if(.not. jump) call ph2xyz(x(i),y(i),z(i),kind)
          call xyz8467(x(i),y(i),z(i),kind)
          if( kind .eq. 0)then
                kind=1
          else
                kind=0
          endif
          if( outcode .ne. 2 .and. outcode .ne. 6)
     *       call xyz2ph(x(i),y(i),z(i),kind)
        endif
        if( outcode .eq. 1 .or. outcode .eq. 5)then
          x(i)=x(i)*180.d0/pi
          y(i)=y(i)*180.d0/pi
        else if( outcode .eq. 2)then
           kind=1
           if( .not. change)call ph2xyz(x(i),y(i),z(i),kind)
        else if( outcode .eq. 6)then
           kind=0
           if( .not. change)call ph2xyz(x(i),y(i),z(i),kind)
        else if( outcode .eq. 3)then
           kind=1
           call ph2tm(x(i),y(i),z(i),kind,0)
        else if( outcode .eq. 7)then
           kind=0
           call ph2tm(x(i),y(i),z(i),kind,0)
        else if( outcode .eq. 4)then
           kind=1
           call ph2tm(x(i),y(i),z(i),kind,1)
        else if( outcode .eq. 8)then
           kind=0
           call ph2tm(x(i),y(i),z(i),kind,1)
        endif
      enddo
        if(FileIn)then
           do i=1,datnum
             write(2,fmt=outfmt)x(i),y(i),z(i),comm1(i),comm2(i)
     *             ,comm3(i)
           enddo
           close(2)
        else
             write(*,fmt=outfmt)x(1),y(1),z(1)
        endif
        write(*,'(15x,a\)')'Do it again (Y/N) : '
        read(*,'(a)')junk
        if(junk .eq. 'y' .or. junk .eq. 'Y')goto 1
        stop
        end



C     KIND 1 CONVERT  XYZ84 --> XYZ67
C     KIND 0 CONVERT  XYZ67 --> XYZ84

      SUBROUTINE XYZ8467(X,Y,Z,KIND)
      REAL*8  XT,YT,ZT,XO,YO,ZO,EX,EY,EZ,S,XK,YK,ZK,x,y,z
      integer kind

        XO=764.558D0
        YO=361.229D0
        ZO=178.374D0
        EX=0.0000011698d0
        EY=-.0000018398d0
        EZ=-.0000009822d0
        S =-.00002329
        XK=-3000808.169d0
        YK= 5011437.601d0
        ZK= 2550970.447d0

       IF (KIND .EQ. 1)THEN
          X=X-XK
          Y=Y-YK
          Z=Z-ZK
          XT=XO+(1.d0+S)*(X+EZ*Y-EY*Z)+XK
          YT=YO+(1.d0+S)*(-EZ*X+Y+EX*Z)+YK
          ZT=ZO+(1.d0+S)*(EY*X-EX*Y+Z)+ZK
       ELSE
          XT=(X-XO-XK)/(1.D0+S)-EZ*(Y-YO-YK)/(1.D0+S)+EY*
     *       (Z-ZO-ZK)/(1.D0+S)+XK
          YT=EZ*(X-XO-XK)/(1.D0+S)+(Y-YO-YK)/(1.D0+S)-EX*
     *       (Z-ZO-ZK)/(1.D0+S)+YK
          ZT=-EY*(X-XO-XK)/(1.D0+S)+EX*(Y-YO-YK)/(1.D0+S)+
     *       (Z-ZO-ZK)/(1.D0+S)+ZK
       ENDIF
          X=XT
          Y=YT
          Z=ZT
       RETURN
       END



C   KIND = 0 -> 67
C   KIND = 1 -> 84
C   XYZ -> PH(RADIAN)

       SUBROUTINE XYZ2PH(X,Y,Z,KIND)
       INTEGER      KIND,I
       REAL*8   X,Y,Z,PI,A,F,EE,EE2,P,L,HY(6),N(6),DHY,KEY,NT,H
       PI=3.141592653589793D0
       IF (KIND .EQ. 0) THEN
          A=6378160.D0
          F=298.25D0
       ELSE
          A=6378137.D0
          F=298.257223563D0
       ENDIF
       EE=1.D0-(1.D0-1.D0/F)**2.D0
       EE2=1.D0/(1.D0-EE)-1.D0
       P=(X**2.D0+Y**2.D0)**0.5D0
       L=DATAN(Y/X)
       HY(1)=DATAN(Z/((1.d0-EE)*P))
       DO 66 I=1,5
          N(I)=A/((1.D0-EE*DSIN(HY(I))**2.D0)**0.5D0)
          HY(I+1)=DATAN((Z+EE*N(I)*DSIN(HY(I)))/P)
          DHY=DABS(HY(I+1)-HY(I))
          IF(DHY .LT. 1.D-11)THEN
             KEY=HY(I+1)
             GO TO 100
          END IF
 66    CONTINUE
C
 100   NT=A/((1.D0-EE*DSIN(KEY)**2.D0)**0.5D0)
       H=P/DCOS(KEY)-NT
       IF(X .LT. 0)THEN
          L=L+PI
       ELSE IF(X .GT. 0 .AND. Y .LT. 0)THEN
          L=L+2.D0*PI
       END IF
       X = L
       Y = KEY
       Z = H
       RETURN
       END

C     A SUBROUTINE FOR ****  PH(RADIAN) -> XYZ  ****
C     KIND = 0 ->67
C     KIND = 1 -> 84

       SUBROUTINE PH2XYZ(LD,PHI,H,KIND)
       INTEGER    KIND
       REAL*8 PHI,LD,H,AS,F,BS,EE,DM,RN,X,Y,Z
       IF (KIND.EQ.0) THEN
          AS=6378160.D0
          F=298.25D0
       ELSE
          AS=6378137.D0
          F=298.257223563D0
       END IF
       F=1.d0/F
       BS=AS*(1.D0-F)
       EE=1.D0-(BS/AS)**2.D0
       DM=DSQRT(1.D0-EE*(DSIN(PHI))**2.D0)
       RN=AS/DM
       X=(RN+H)*DCOS(PHI)*DCOS(LD)
       Y=(RN+H)*DCOS(PHI)*DSIN(LD)
       Z=(RN*(BS/AS)**2.d0+H)*DSIN(PHI)
       PHI= Y
       LD = X
       H  = Z
       RETURN
       END


C      KIND = 0 -> 67
C      KIND = 1 -> 84
C      TorU = 0 -> TM2
C      TorU = 1 -> UTM
C       CONVERT PH(RADIAN) -> TM2 OR UTM

      SUBROUTINE PH2TM(LM2,PH2,HH,KIND,TorU)
      INTEGER      TorU,KIND
      REAL*8   PH2,HH,A,F,EE,EE2,SHIFT,LM1,SCALE,PH1,RAP,LM2,N
      REAL*8   UU,T,PH,PM
      IF (KIND .EQ. 0) THEN
         A=6378160.D0
         F=298.25D0
      ELSE
         A=6378137.D0
         F=298.257223563D0
      END IF
      EE=1.D0-(1.D0-1.D0/F)**2.D0
      EE2=1.D0/(1.D0-EE)-1.D0
      IF( TORU .EQ. 0) THEN
         SHIFT=250000.D0
         LM1=121.D0
         SCALE=0.9999D0
      ELSE
         SHIFT=500000.D0
         LM1=123.D0
         SCALE=0.9996D0
      END IF
      LM1=LM1/57.29577951308D0
      PH1=0.d0
      RAP=LM2-LM1
      UU=EE2*DCOS(PH2)**2.D0
      N=A/DSQRT(1.d0-EE*DSIN(PH2)**2.D0)
      T=DSIN(PH2)/DCOS(PH2)
      LM2=N*DCOS(PH2)*RAP+N*(1.d0-T**2.D0+UU)*DCOS(PH2)**3.D0*
     *    RAP**3.D0/6.D0+N*(5.D0-18.D0*T**2.d0+T**4.d0+14.D0*UU
     *    -58.D0*UU*T*T)*DCOS(PH2)**5.d0*RAP**5.d0/120.D0
      LM2=SCALE*LM2+SHIFT
      PH=PM(PH2,PH1,EE,A)+N*T*RAP**2.d0*DCOS(PH2)**2.d0/2.d0+
     *   N*T*(5.D0-T**2.d0+9.D0*UU+4.d0*UU**2.d0)*DCOS(PH2)**4.d0*RAP
     *   **4.d0/24.D0+N*T*DCOS(PH2)**6.d0*RAP**6.d0*(61.D0-58.D0*T*T+
     *   T**4.d0+270.D0*UU-330.D0*UU*T*T)/720.D0
      HH=(DSIN(PH2)*RAP+DSIN(PH2)*DCOS(PH2)**2.d0*(1.d0+3.D0*UU+
     *   2.d0*UU**2.d0)*RAP**3.d0/3.D0+T*DCOS(PH2)**5.d0*(2.d0-
     *   T**2.d0)*RAP**5.d0/15.D0)*206264.80625D0
      PH2=SCALE*PH
      RETURN
      END


C      CONVERT    TM2 OR UTM -> PHL
C      KIND = 0 -> 67
C      KIND = 1 -> 84
C      TorU = 0 -> TM2
C      TorU = 1 -> UTM

      SUBROUTINE TM2PH(CE,CN,KIND,TorU)
      INTEGER      TorU,KIND
      REAL*8   CN,CE,A,F,E,SHIFT,LM1,SCALE,A1,AA,RM,RN,YM,S,T,U,V
      REAL*8   BB,CC,DD,C1,C2,PH,DYM,R1,R2,T1,Z,CPH,T2,CLA
      RM(PH,A1,E)=A1*(1.d0-E)/(1.d0-E*DSIN(PH)**2.d0)**1.5D0
      RN(PH,A1,E)=A1/(1.d0-E*DSIN(PH)**2.d0)**0.5D0
      YM(PH,S,T,U,V)=S*PH+T*DSIN(2.d0*PH)+U*DSIN(4.d0*PH)+
     *               V*DSIN(6.d0*PH)
      IF (KIND .EQ. 0) THEN
         A=6378160.D0
         F=298.25D0
      ELSE
         A=6378137.D0
         F=298.257223563D0
      END IF
      IF( TORU .EQ. 0) THEN
         SHIFT=250000.D0
         LM1=121.D0
         SCALE=0.9999D0
      ELSE
         SHIFT=500000.D0
         LM1=123.D0
         SCALE=0.9996D0
      END IF
      E=2.D0/F-1.D0/F**2.d0
      A1=A*(1.d0-E)
      AA=A1*(1.D0+0.75D0*E+45.D0*E**2.D0/64.D0+175.D0*E**3.D0/256.D0)
      BB=A1*(0.75D0*E+15.D0*E**2.D0/16.D0+525.D0*E**3.D0/512.D0)*
     *   (-0.5D0)
      CC=A1*(15.D0*E**2.D0/64.D0+106.D0*E**3.D0/256.D0)*0.25D0
      DD=A1*(35.D0*E**3.D0/512.D0)/(-6.D0)
      C1=CN/SCALE
      C2=(CE-SHIFT)/SCALE
      PH=C1/AA
   10 DYM=C1-YM(PH,AA,BB,CC,DD)
      PH=PH+DYM/(AA+2.D0*BB*DCOS(2.D0*PH)+4.D0*CC*DCOS(4.D0*PH)
     *   +6.D0*DD*DCOS(6.D0*PH))
      IF(ABS(DYM)-0.0001D0)15,15,10
   15 R1=RM(PH,A,E)
      R2=RN(PH,A,E)
      T1=DTAN(PH)
      Z=C2/R2
      CPH=PH-Z*C2/R1*0.5D0*T1+Z**3.D0*C2/R1/24.D0*T1*
     *    (5.D0+3.D0*T1**2.D0)
      T2=T1**2.D0
      CLA=(Z-Z**3.D0/6.D0*(R2/R1+2.D0*T2)+Z**5.D0/120.D0*
     *    (5.D0+28.D0*T2+24.D0*T2**2.D0))/DCOS(PH)
      CLA=CLA+LM1*3.141592653589793D0/180.D0
      CN=CPH
      CE=CLA
      RETURN
      END


C
C     A FUNCTION USED BY PH2TM
C
      FUNCTION PM(P2,P1,E1,RA)
      REAL*8 P2,P1,E1,RA,A,B,C,D,E,PM
      A=1.D0+0.75D0*E1+45.D0*E1**2.D0/64.D0+175.D0*E1**3.D0/256.D0+
     *  11025.D0*E1**4.D0/16384.D0
      B=0.75D0*E1+15.D0*E1**2.D0/16.D0+525.D0*E1**3.D0/512.D0+
     *  2205.D0*E1**4.D0/2048.D0
      C=15.D0*E1**2.D0/64.D0+105.D0*E1**3.D0/256.D0+
     *  2205.D0*E1**4.D0/4096.D0
      D=35.D0*E1**3.D0/512.D0+315.D0*E1**4.D0/2048.D0
      E=315.D0*E1**4.D0/16384.D0
      PM=RA*(1.d0-E1)*(A*(P2-P1)-B*(DSIN(2.d0*P2)-DSIN(2.d0*P1))/2.d0+C*
     *   (DSIN(4.d0*P2)-DSIN(4.d0*P1))/4.d0-D*(DSIN(6.D0*P2)-
     *   DSIN(6.D0*P1))/6.D0+E*(DSIN(8.d0*P2)-DSIN(8.d0*P1))/8.d0)
      RETURN
      END
***