c layer one and two are flipped now
      implicit real*8(a-h,o-z)

      parameter(ntyp=139)
      dimension amapt(72),amapth(72),amaps(72),amapi(72)
      dimension amapl(6,72),flev(ntyp,7),flev2(7)
      dimension elev(72,36),ilon(72),il(36)

      character*2 ctype(ntyp),line*506,dum*1,dum0*5
      character*2 types(72)


c... read in elevation map
c......................... 
      open(3,file='CNelevatio.txt')
      open(19,file='map.topo')
      read(3,899)line
      do 200 j=1,36
c        read(3,*)ilon
         read(3,899)line
         read(line,*)il(j),(elev(i,j),i=1,72)
         write(19,902)(elev(i,j)/1000.,i=37,72),
     +   (elev(i,j)/1000.,i=1,36)
 200  continue
      close(3)
      close(19)

      open(2,file='CNtype_key.txt')
      open(7,file='CNtype.txt')
      open(8,file='map.t7')
      open(16,file='map.t0')
      open(9,file='map.t1')
      open(10,file='map.t2')
      open(11,file='map.t3')
      open(12,file='map.t4')
      open(13,file='map.t5')
      open(14,file='map.t6')
      open(15,file='map.thick')
      open(17,file='map.sed')
      open(18,file='map.ice')

c... read in key for crust types
c...............................
      read(2,890)dum
      do 101 i=1,ntyp
         read(2,899)ctype(i)
         print 899,ctype(i)
         read(2,891)dum
         read(2,899)line
         read(line,*)(flev(i,l),l=1,7),dum0,fthick
c. switch layer 1 and 2
c. new layer 1 is water
c. new layer 2 is ice
         aux=flev(i,2)
         flev(i,2)=flev(i,1)
         flev(i,1)=aux
         felev=0.
         do 102 l=1,7
 102     felev=felev+flev(i,l)
         if(felev.ne.fthick)then
           print*,' read error in thickness of crust'
           print 905,' type: ',i,ctype(i)
           print*,felev,fthick 
           stop
         endif
 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,ntyp
            if(types(i).eq.ctype(l))then
              do 104 m=1,7
 104          flev2(m)=flev(l,m)
              amaps(i)=flev2(3)+flev2(4)
              amapi(i)=flev2(2)
              if(elev(i,j).lt.0.)then
                flev2(1)=-elev(i,j)/1000.
                elev(i,j)=0.      
              endif
              do 103 k=2,7
 103          flev2(k)=flev2(k)+flev2(k-1)
              amapt(i)=elev(i,j)/1000.-flev2(7) 
              amapth(i)=flev2(7) 
c uncomment next line if thickness without water is desired
              if(elev(i,j).eq.0.)amapth(i)=amapth(i)-flev2(1)
              do 30 k=1,6
 30           amapl(k,i)=elev(i,j)/1000.-flev2(k)
              goto 10
            endif
 20         continue
            print*,' crust type code not found: ',types(i)
            print*,' latitude: ',ilat,' long index: ',i
 10      continue
         write(8,902)(amapt(l),l=37,72),(amapt(l),l=1,36)
         write(15,902)
     +        (amapth(l),l=37,72),(amapth(l),l=1,36)
         write(16,902)(elev(l,j)/1000.,l=37,72)
     +        ,(elev(l,j)/1000.,l=1,36)
         write(17,902)(amaps(l),l=37,72),(amaps(l),l=1,36)
         write(18,902)(amapi(l),l=37,72),(amapi(l),l=1,36)
         do 50 k=1,6
 50      write(8+k,902)(amapl(k,l),l=37,72),(amapl(k,l),l=1,36)
c     print 900, types
 40   continue

 890  format(////a)
 891  format(//a)
 899  format(a)
 900  format(72(1x,a2))
 901  format(72f7.0)
 902  format(30(8e15.5/))
 903  format(73i7)
 904  format(i4,72(2x,a2))
 905  format(a,i6,1x,a2)
 99   continue
      close(7)
      end