c layer one and two flipped, after the read statement!
c layer 1: water
c layer 2: ice

      parameter(ityp=139)
      dimension fvel(ityp,8),fvels(ityp,8),frho(ityp,8)
      dimension amapvp(8,72,36),amaprho(8,72,36),
     +          amapvs(8,72,36)
      character*2 ctype(ityp),line*506,dum*1,dum0*5
      character*2 types(72)
      open(2,file='CNtype_key.txt')
      open(7,file='CNtype.txt')

c... read in key for crust types
c...............................
      read(2,890)dum
      do 101 i=1,ityp
         read(2,899)ctype(i)
         print 899,ctype(i)
         read(2,899)line
         read(line,*)(fvel(i,l),l=1,8)
         read(2,899)line
         read(line,*)(fvels(i,l),l=1,8)
         read(2,899)line
         read(line,*)(frho(i,l),l=1,8)
c flip layers
         aux=fvel(i,1)
         fvel(i,1)=fvel(i,2)
         fvel(i,2)=aux
         aux=fvels(i,1)
         fvels(i,1)=fvels(i,2)
         fvels(i,2)=aux
         aux=frho(i,1)
         frho(i,1)=frho(i,2)
         frho(i,2)=aux
         read(2,899)dum
 101  continue

c... read CNtype file
c...............................
      read(7,899)line
      read(line,*)flons
      do 40 j=1,36
         read(7,899)line       
         read(line,*)ilat,types
         print*,ilat
         do 10 i=1,72
            do 20 l=1,ityp
            if(types(i).eq.ctype(l))then
              do 30 k=1,8
              amapvp(k,i,j)=fvel(l,k)
              amapvs(k,i,j)=fvels(l,k)
              amaprho(k,i,j)=frho(l,k)
c             amapvp(1,i,j)=fvel(l,1)
c             amapvp(2,i,j)=fvel(l,2)
c             do 30 k=3,8
c             amapvp(k,i,j)=fvel(l,k)
c take poisson's ratio 0.27 for average crust
c             amapvs(k,i,j)=fvel(l,k)/1.78
c take average value of birch's law from table 6
c             amaprho(k,i,j)=0.391+0.391*fvel(l,k)
c30           continue
c put in values of water layer explicitly
c             amapvs(1,i,j)=0.
c             amaprho(1,i,j)=1.03
c put in values of ice layer explicitly
c             amapvs(2,i,j)=1.9
c             amaprho(2,i,j)=0.98
 30           continue
              goto 10
            endif
 20         continue
            print*,' crust type code not found: ',types(i)
            print*,' latitude: ',ilat,' long index: ',i
 10      continue
 40   continue
      do 50 k=1,8
         write(dum,'(i1)')k
         open(8+k,file='map.vp'//dum)
         open(18+k,file='map.vs'//dum)
         open(28+k,file='map.rho'//dum)
         do 60 i=1,36
         write(8+k,902)(amapvp(k,l,i),l=37,72),
     +                 (amapvp(k,l,i),l=1,36)
         write(18+k,902)(amapvs(k,l,i),l=37,72),
     +                 (amapvs(k,l,i),l=1,36)
         write(28+k,902)(amaprho(k,l,i),l=37,72),
     +                 (amaprho(k,l,i),l=1,36)
 60      continue
         close(8+k)
         close(18+k)
         close(28+k)
 50   continue

 890  format(////a)
 891  format(//a)
 899  format(a)
 900  format(72(1x,a2))
 901  format(72f7.0)
 902  format(30(8e15.5/))
 99   continue
      close(7)
      end