   ====================================================================
c              READ MODEL 3d
c   ====================================================================
c
      subroutine rdmd3d(vpe,vse,den,nxt,nyt,nzt,ifre)

ccc   routine to read the model 
c     nxt     # of node points in x-direction (integer)(sent)
c     nyt     # of node points in y-direction (integer)(sent)
c     nzt     # of node points in z-direction (integer)(sent)
c     vpe(1)  lowest p-velocity               (real)   (returned)
c     vpe(2)  highest p-velocity              (real)   (returned)
c     vse(1)  lowest s-velocity               (real)   (returned)
c     vse(2)  highest s-velocity              (real)   (returned)
c     den(1)  lowest density                  (real)   (returned)
c     den(2)  highest density                 (real)   (returned)
c     ifre    flag; =0 flat fs, <>0  topo     (integer)(returned)

      include 'parstat'
      dimension vpe(*),vse(*),den(*)
      vpe(2)=0.
      vpe(1)=1.0E10 
      vse(2)=0.
      den(2)=0.
      den(1)=1.0E10 
c
c   HOMOGENOUS MEDIUM
c
c    open(10,file='PVELOCITY',form='unformatted')
c      open(11,file='SVELOCITY',form='unformatted')
c      open(12,file='DENSITY',form='unformatted')


     


c     do 61 k=nz,1,-1
c        do 61 j=1,ny
c           do 61 i=1,nx
c               vp=1.73
c               vs=1.
c               dd=1.
c              vp=6000.
c              vs=3464.
c              dd=2670.
c$$$               read(10) vp
c$$$               read(11) vs
c$$$               read(12) dd
 
c LAYERED MEDIUM 

open(10,file='model.dat')

do 61 k=nz,1,-1
read(10,*) vp,vs,dd
do 61 j=1,ny
do 61 i=1,nx

               
               if (vp.gt.vpe(2)) vpe(2)=vp
               if (vp.lt.vpe(1)) vpe(1)=vp
               if (vs.gt.vse(2)) vse(2)=vs
               if (vs.lt.vse(1)) vse(1)=vs
               if (dd.gt.den(2)) den(2)=dd
               if (dd.lt.den(1)) den(1)=dd
               
               lam1(i,j,k)=dd*(vp**2-2.*vs**2)
               mu1(i,j,k)=dd*vs**2
               d1(i,j,k)=dd
 61         continue
 
 close(10)
