anchor.diff:       1                                                    +v(3)*dsin(la)
         vn = vn*111.19d0
         ve = ve*111.19d0
  -      print *,stlo,stla,ve,vn
  +      write(*,*) stlo,stla,ve,vn
         goto 102
     101 continue
         END
   
   C----------------------------------------------------------------------C
   
  -      SUBROUTINE morvel(pole,polec)
  -      IMPLICIT none
  -      CHARACTER*80 polec,plate(2)
  -      REAL*8 pole(3),x(2,3)
  -      INTEGER i
  -      read(polec,*) plate(1),plate(2)
  -      do 101 i = 1,2
  -          if (plate(i).eq.'AM') then
  -              x(i,1) = 0.04829
  -              x(i,2) = -0.37592
  -              x(i,3) = 0.84817
  -          elseif (plate(i).eq.'AN') then
  -              x(i,1) = 0.07204
  -              x(i,2) = -0.35489
  -              x(i,3) = 0.81026
  -          elseif (plate(i).eq.'AR') then
  -              x(i,1) = 0.48451
  -              x(i,2) = -0.31742
  -              x(i,3) = 1.00388
  -          elseif (plate(i).eq.'AU') then
  -              x(i,1) = 0.53463
  -              x(i,2) = 0.05931
  -              x(i,3) = 0.93478
  -          elseif (plate(i).eq.'CA') then
  -              x(i,1) = 0.11032
  -              x(i,2) = -0.49683
  -              x(i,3) = 0.74786
  -          elseif (plate(i).eq.'CO') then
  -              x(i,1) = -0.48084
  -              x(i,2) = -1.14521
  -              x(i,3) = 1.12543
  -          elseif (plate(i).eq.'CP') then
  -              x(i,1) = 0.52040
  -              x(i,2) = -0.09291
  -              x(i,3) = 1.00856
  -          elseif (plate(i).eq.'EU') then
  -              x(i,1) = 0.07941
  -              x(i,2) = -0.40365
  -              x(i,3) = 0.75055
  -          elseif (plate(i).eq.'IN') then
  -              x(i,1) = 0.46728
  -              x(i,2) = -0.28306
  -              x(i,3) = 1.00170
  -          elseif (plate(i).eq.'JF') then
  -              x(i,1) = 0.49382
  -              x(i,2) = 0.38353
  -              x(i,3) = -0.00681
  -          elseif (plate(i).eq.'LW') then
  -              x(i,1) = 0.18270
  -              x(i,2) = -0.42829
  -              x(i,3) = 0.80758
  -          elseif (plate(i).eq.'NA') then
  -              x(i,1) = 0.15483
  -              x(i,2) = -0.46836
  -              x(i,3) = 0.56521
  -          elseif (plate(i).eq.'NB') then
  -              x(i,1) = 0.19329
  -              x(i,2) = -0.44604
  -              x(i,3) = 0.79881
  -          elseif (plate(i).eq.'MQ') then
  -              x(i,1) = 0.85484
  -              x(i,2) = -0.11982
  -              x(i,3) = 1.44874
  -          elseif (plate(i).eq.'NZ') then
  -              x(i,1) = 0.02874
  -              x(i,2) = -0.73542
  -              x(i,3) = 1.08528
  -          elseif (plate(i).eq.'PS') then
  -              x(i,1) = 0.66046
  -              x(i,2) = -0.59193
  -              x(i,3) = -0.07185
  -          elseif (plate(i).eq.'RI') then
  -              x(i,1) = -1.14352
  -              x(i,2) = -4.32660
  -              x(i,3) = 2.15245
  -          elseif (plate(i).eq.'SA') then
  -              x(i,1) = 0.08198
  -              x(i,2) = -0.35590
  -              x(i,3) = 0.54095
  -          elseif (plate(i).eq.'SC') then
  -              x(i,1) = 0.08341
  -              x(i,2) = -0.39306
  -              x(i,3) = 0.63896
  -          elseif (plate(i).eq.'SM') then
  -              x(i,1) = 0.14189
  -              x(i,2) = -0.48048
  -              x(i,3) = 0.84261
  -          elseif (plate(i).eq.'SR') then
  -              x(i,1) = 0.08813
  -              x(i,2) = -0.34742
  -              x(i,3) = 0.52526
  -          elseif (plate(i).eq.'SU') then
  -              x(i,1) = 0.10210
  -              x(i,2) = -0.47858
  -              x(i,3) = 0.84111
  -          elseif (plate(i).eq.'SW') then
  -              x(i,1) = 1.06493
  -              x(i,2) = -0.97101
  -              x(i,3) = -0.09671
  -          elseif (plate(i).eq.'YZ') then
  -              x(i,1) = 0.05326
  -              x(i,2) = -0.39837
  -              x(i,3) = 0.88007
  -          elseif (plate(i).eq.'ITRF') then
  -              x(i,1) = 0.11257
  -              x(i,2) = -0.28145
  -              x(i,3) = 0.60534
  -          elseif (plate(i).eq.'PA') then
  -              x(i,1) = 0.0
  -              x(i,2) = 0.0
  -              x(i,3) = 0.0
  -          else
  -              write(0,*) '!! Error: no plate named '//trim(plate(i))
  -              write(0,*) 'AM: Amur'
  -              write(0,*) 'AN: Antarctica'
  -              write(0,*) 'AR: Arabia'
  -              write(0,*) 'AU: Australia'
  -              write(0,*) 'CA: Caribbean'
  -              write(0,*) 'CO: Cocos'
  -              write(0,*) 'CP: Capricorn'
  -              write(0,*) 'EU: Eurasia'
  -              write(0,*) 'IN: India'
  -              write(0,*) 'ITRF: International reference frame'
  -              write(0,*) 'JF: Juan de Fuca'
  -              write(0,*) 'LW: Lwandle'
  -              write(0,*) 'MQ: Macquarie'
  -              write(0,*) 'NA: North America'
  -              write(0,*) 'NB: Nubia'
  -              write(0,*) 'NZ: Nazca'
  -              write(0,*) 'PA: Pacific'
  -              write(0,*) 'PS: Philippine Sea'
  -              write(0,*) 'RI: Rivera'
  -              write(0,*) 'SA: South America'
  -              write(0,*) 'SC: Scotia'
  -              write(0,*) 'SM: Somalia'
  -              write(0,*) 'SR: Sur'
  -              write(0,*) 'SU: Sunda'
  -              write(0,*) 'SW: Sandwich'
  -              write(0,*) 'YZ: Yangtze'
  -              call usage('')
  +      SUBROUTINE getpole(pole,polec,model,verbos)
  +      IMPLICIT none
  +      CHARACTER*80 polec
  +      REAL*8 pole(3),x(3)
  +      INTEGER i,isave,model,verbos
  +
  +      if (verbos.eq.1) then
  +          write(0,*)
  +          write(0,*) 'getpole says: Getting Euler pole'
  +      endif
  +
  +      ! Print available plates if requested
  +      if (model.eq.-1) then
  +          call morvelplates()
  +      elseif (model.eq.-2) then
  +          call nuvel1aplates()
  +      elseif (model.eq.-3) then
  +          call morvel56plates()
  +      elseif (model.lt.0) then
  +          write(0,*) '!! Error: no plate model number ',abs(model)
  +          call usage('!! Exiting')
  +      endif
  +
  +      ! Parse the pole from the command line
  +      if (verbos.eq.1) then
  +          write(0,*) 'getpole says: Parsing plates/pole from input: ',
  +     1               trim(polec)
  +      endif
  +      if (polec.eq.'none') call usage('!! Error: no plate/pole defined')
  +
  +      ! Input format is either "P1/P2" or "lon/lat/vel"
  +      ! Check for first "/"
  +      i = index(polec,'/')
  +      if (i.gt.0) then
  +          polec(i:i) = ' '
  +          isave = i
  +      else
  +          call usage('!! Error: input must contain at least one "/"')
  +      endif
  +
  +      ! Check for second "/" to determine operation mode
  +      i = index(polec,'/')
  +      if (i.gt.0) then
  +          polec(i:i) = ' '
  +          i = -1
  +      else
  +          i = 1
  +      endif
  +
  +      ! Compute Cartesian Euler pole
  +      if (i.gt.0) then
  +          if (verbos.eq.1) then
  +              write(0,*) 'getpole says: Using plate circuit model'
             endif
  -  101 continue
  -      x(1,1) = x(2,1) - x(1,1)
  -      x(1,2) = x(2,2) - x(1,2)
  -      x(1,3) = x(2,3) - x(1,3)
  -      pole(3) = dsqrt(x(1,1)*x(1,1)+x(1,2)*x(1,2)+x(1,3)*x(1,3))
  -      pole(2) = dacos(x(1,3)/pole(3))
  -      pole(2) = 90.0d0 - pole(2)*1.8d2/3.14159265d0
  -      pole(1) = datan2(x(1,2),x(1,1))*1.8d2/3.14159265d0
  +          call model2pole(pole,polec,model,verbos)
  +      elseif (i.lt.0) then
  +          if (verbos.eq.1) then
  +              write(0,*) 'getpole says: Reading pole from command line'
  +          endif
  +          read(polec,*) pole(1),pole(2),pole(3)
  +      else
  +          call usage('!! Error in parsing pole inputs')
  +      endif
  +      call lalomag2xyz(x(1),x(2),x(3),pole(2),pole(1),pole(3))
  +      pole(1) = x(1)
  +      pole(2) = x(2)
  +      pole(3) = x(3)
  +      if (verbos.eq.1) then
  +          write(0,1003) pole(1),pole(2),pole(3)
  +      endif
  + 1003 format(' getpole says: Euler pole is: ',3(1PE14.6))
  +
         RETURN
         END
   
  
positive.code: 
        SUBROUTINE usage(str)
        IMPLICIT none
        INTEGER d
        CHARACTER str*(*)
        d = 0
        if (str.eq.'d') d = 1
        write(*,*)
       1 'Usage: platemotion -plate P1/P2 | -pole lon/lat/vel',
       2                       ' [-f IFILE] [-o OFILE] [-h]'
        write(*,*)
       1 '  -plate P1/P2      Define plates (P2 moving w.r.t. fixed P1)'
        write(*,*)
       1 '  -pole LON/LAT/VEL Define Euler pole (velocity in deg/Ma)'
        write(*,*)
       1 '  -f IFILE          Input file (default: run with stdin)'
        write(*,*)
       1 '                        lon lat'
  
negative.code:       write(*,*)
       1 '                        v_east v_north (mm/yr)'
        write(*,*)
       1 '  -verbose          Turn on verbose mode (for debugging)'
        write(*,*)
       1 '  -h                Online help'
        write(*,*)
        STOP
        END
  