c  ADDPROJ - written by Jeff Phillips
      character*50 ifile, image, ofile
      character*80 line
      character*1 ptype
      dimension id(14),pgm(2)
      print*
      print*,'Append projection information to a PDS image label file.'
      print*
      print '(a\)',' Enter image name:'
      read(5,100) image
  100 format(a)
      l=len_trim(image)
      ofile=image(1:l)//'.lbl'
c      open(10,file=ofile,form='formatted',status='old',access='append')
      open(10,file=ofile,form='formatted',status='old')
c      ofile=image(1:l)//'.dat'
c      print *,'enter color image name:'
c      read(5,100) image
c      l=len_trim(image)
c      ofile=image(1:l)//'.lbl'
c      open(11,file=ofile,form='formatted',status='unknown')
   10 read(10,100,end=11) line
      if(line(1:10).eq.'Projection') then
        backspace(10)
        go to 12
      else if (line(1:3 ).eq.'END') then
        backspace(10)
        go to 12
      else
        go to 10
      endif
   11 continue
      close(10)
      open(10,file=ofile,form='formatted',status='old',access='append')
   12 continue
      print*,'Available projections:'
      print*,'   Albers equal area elliptical'
      print*,'   Ellipsoidal transverse mercator'
      print*,'   None'
      print*
   13 print '(a\)',' Enter projection type (A,E,N):'
      read(5,100) ptype
      if(ptype.eq.'a'.or.ptype.eq.'A') then
        ptype='a'
        write(10,101)
      else if(ptype.eq.'e'.or.ptype.eq.'E') then
        ptype='e'
        write(10,102)
      else if(ptype.eq.'n'.or.ptype.eq.'N') then
        ptype='n'
        write(10,103)
      else
        print*,'invalid choice'
        go to 13
      endif
  101 format('Projection = Albers equal area elliptical')
  102 format('Projection = Ellipsoidal transverse mercator')
  103 format('Projection = None')
      if(ptype.ne.'n') then
        print '(a\)',' Enter base latitude: '
        read(5,*) baslat
        write(10,104) baslat
  104 format(f10.3,'  Lat0')
        if(ptype.eq.'a') then
          print*,'Standard latitudes for USGS Albers projections:'
          print*,'29.5 45.5 US'
          print*,'55.0 65.0 Alaska'
          print*,' 8.0 18.0 Hawaii'
          print*
          print '(a\)',' Enter 1st standard latitude: '
          read(5,*) baslat
          write(10,105) baslat
  105 format(f10.3,'  Lat1')
          print '(a\)',' Enter 2nd standard latitude: '
          read(5,*) baslat
          write(10,106) baslat
  106 format(f10.3,'  Lat2')
        endif
        print '(a\)',' Enter central longitude: '
        read(5,*) baslat
        write(10,107) baslat
  107 format(f10.3,'  Lon0')
        if(ptype.eq.'e') then
          print '(a\)',' Enter scale on the central meridian (0.9996 for
     1 the UTM projection): '
          read(5,*) baslat
          write(10,108) baslat
  108 format(f10.4,'  k0')
        endif
      endif
   14 print '(a\)',' Enter the grid file used to make the image: '
      read(5,100) ifile
      open(11,file=ifile,form='unformatted',status='old',err=15)
      go to 16
   15 print*,'file not found'
      go to 14
   16 read(11)id,pgm,nc,nr,nz,xo,dx,yo,dy
      baslat=(xo-dx/2.)*1000.
      write(10,109) baslat
  109 format(f10.1,'  X0')
      baslat=(yo+(nr-1)*dy+dy/2.)*1000.
      write(10,110) baslat
  110 format(f10.1,'  Y0')
      baslat=(dx+dy)*500.
      write(10,111) baslat
  111 format(f10.1,'  PixelSize (meters)')
      if(ptype.ne.'n') then
   20   print*,'Available ellipsoids:'
        print*,'0  GRS 80        1980 6378137.0 6356752.3  Newly adopted
     1'
        print*,'1  WGS 72        1972 6378135.0 6356750.5  NASA, DoD, oi
     1l Co.'
        print*,'2  Australian    1965 6378160.0 6356774.7  Australia'
        print*,'3  Krasovsky     1940 6378245.0 6356863.0  Soviet Union'
        print*,'4  International 1924 6378388.0 6356911.9  Remainder of
     1the World'
        print*,'5  Hayford       1924 6378388.0 6356911.9  Remainder of
     1the World'
        print*,'6  Clarke        1880 6378249.1 6356514.9  Most of Afric
     1a; France'
        print*,'7  Clarke        1866 6378206.4 6356583.8  North America
     1;'
        print*,'                                           Philippines'
        print*,'8  Airy          1830 6377563.4 6356256.9  Great Britain
     1'
        print*,'9  Everest       1830 6377276.3 6356075.4  India; Burma;
     1 Pak.;'
        print*,'                                           Afgan.; Thail
     1and;etc.'
        print '(a\)',' Enter the ellipsoid number: '
        read(5,*) num
        if(num.eq.0) then
          write(10,112) num
  112 format(8x,i2,'  ellipsoid (GRS 80 1980)')
        else if(num.eq.1) then
          write(10,113) num
  113 format(8x,i2,'  ellipsoid (WGS 72 1972)')
        else if(num.eq.2) then
          write(10,114) num
  114 format(8x,i2,'  ellipsoid (Australian 1965)')
        else if(num.eq.3) then
          write(10,115) num
  115 format(8x,i2,'  ellipsoid (Krasovsky 1940)')
        else if(num.eq.4) then
          write(10,116) num
  116 format(8x,i2,'  ellipsoid (International 1924)')
        else if(num.eq.5) then
          write(10,117) num
  117 format(8x,i2,'  ellipsoid (Hayford 1924)')
        else if(num.eq.6) then
          write(10,118) num
  118 format(8x,i2,'  ellipsoid (Clarke 1880)')
        else if(num.eq.7) then
          write(10,119) num
  119 format(8x,i2,'  ellipsoid (Clarke 1866)')
        else if(num.eq.8) then
          write(10,120) num
  120 format(8x,i2,'  ellipsoid (Airy 1830)')
        else if(num.eq.9) then
          write(10,121) num
  121 format(8x,i2,'  ellipsoid (Everest 1830)')
        else
          print*,'invalid response'
          go to 20
        endif
      endif
   30 print '(a\)',' How many vector sets are to be specified?: '
      read(5,*) num
      if(num.le.0) go to 60
      write(10,122) num
  122 format('VECTOR_SETS = ',i3)
      do 50 i=1,num
      print '(a,i3,a\)',' Enter label to be used for vector set',i,': '
      read(5,100) ifile
      l=len_trim(ifile)
      print '(a\)',' Enter number of files in the set: '
      read(5,*) num1
      write(10,123) ifile(1:l),num1
  123 format(2x,a,' = ',i3)
      print '(a\)',' Enter the color palette number for the set: '
      read(5,*) num2
      do 40 j=1,num1
      print '(a,i3,a\)',' Enter file prefix ',j,': '
      read(5,100) ifile
      write(10,124) ifile,num2
  124 format(5x,a8,' 120.0 0.0 ',i3)
   40 continue
   50 continue
   60 print '(a\)',' How many caption lines are to be specified?: '
      read(5,*) num
      if(num.le.0) go to 90
      write(10,125) num
  125 format('CAPTION = ',i3,' LINES')
      print '(a\)',' Enter the color palette number for the caption: '
      read(5,*) num2
      write(10,126) num2
  126 format('TEXT_COLOR = ',i3)
      print '(a\)',' Enter the text size for the caption (1(max) to 4(mi
     1n)): '
      read(5,*) num2
      write(10,127) num2
  127 format('TEXT_SIZE = ',i1)
      print*,'Enter lines for caption:'
      do 70 i=1,num
      read(5,100) line
      write(10,100) line
   70 continue
   90 write(10,128)
  128 format('END')
      stop
      end
