c
      parameter (nx=600,nz=600, dx=2.5)
      parameter (nmax=750, dt=0.5e-3)
      parameter (nsrc=1, ixs= 300, izs = 199)
      parameter (nsnap=2, ifsnap=650, idsnap=650)
      parameter (isotype=1)
      parameter (iabmax = 80)
      parameter (topbc='abbc')
c
      dimension vx(0:nx,0:nz), vz(0:nx,0:nz)
      dimension txx(0:nx,0:nz), tzz(0:nx,0:nz)
      dimension txz(0:nx,0:nz) 
c
      dimension bx(0:nx,0:nz), bz(0:nx,0:nz)
      dimension rlamu(0:nx,0:nz), rla(0:nx,0:nz)
      dimension rmuxz(0:nx,0:nz), rho(0:nx,0:nz) 
c
      dimension source(nmax), eponge(0:nx)
c
      g1 = -1.0/24.0
      g2 = 9.0/8.0
c
      dtdx=dt/dx
      soufac = dt/(dx*dx)
c
      isnap = 0
c
c ... source signature
c
      call wavelet(nmax,dt,source)
      print *,'wavelet'
c
c ... coefficients for finite-difference
c
      call model(rlamu,rla,rmuxz,rho,bx,bz,dtdx,nx,nz,
     +           bx0,bz0,bz1,rmuxz1,rla0,rla1,rlamu0,rlamu1)
      print *,'model'
c
c ... the damping coefficients for bc
c
      a0=0.3/float(iabmax)
      do ia=0,iabmax
         eponge(ia)=exp(-(a0*float(iabmax-ia))**2)
      enddo
c
c ...     initializes displacements, pressures, boundary conditions
c
      do i=0,nx
         do k=0,nz
            vx(i,k)=0.0
            vz(i,k)=0.0
            txx(i,k)=0.0
            tzz(i,k)=0.0
            txz(i,k)=0.0
         enddo
      enddo
c
c ===============================================================
c     Equations of momentum conservation 
c ===============================================================
c
      do n=1,nmax
c
      do k=2,nz-2
         do i=2,nx-2
            vx(i,k)=vx(i,k)+bx(i,k)*
     *        (g1*(txx(i+1,k)-txx(i-2,k))+
     *        g2*(txx(i,k)-txx(i-1,k))+
     *        g1*(txz(i,k+2)-txz(i,k-1))+
     *        g2*(txz(i,k+1)-txz(i,k)))
            vz(i,k)=vz(i,k)+bz(i,k)*
     *        (g1*(txz(i+2,k)-txz(i-1,k))+
     *        g2*(txz(i+1,k)-txz(i,k))+
     *        g1*(tzz(i,k+1)-tzz(i,k-2))+
     *        g2*(tzz(i,k)-tzz(i,k-1)))
         enddo
      enddo
c
c loose ends
c
      do k=1,1
         do i=2,nx-2
            vx(i,k)=vx(i,k)+bx(i,k)*
     *        (g1*(txx(i+1,k)-txx(i-2,k))+
     *        g2*(txx(i,k)-txx(i-1,k))+
     *        g1*(txz(i,k+2)-txz(i,k-1))+
     *        g2*(txz(i,k+1)-txz(i,k)))
         enddo
      enddo
c
      do k=2,nz-2
         do i=1,1
            vz(i,k)=vz(i,k)+bz(i,k)*
     *        (g1*(txz(i+2,k)-txz(i-1,k))+
     *        g2*(txz(i+1,k)-txz(i,k))+
     *        g1*(tzz(i,k+1)-tzz(i,k-2))+
     *        g2*(tzz(i,k)-tzz(i,k-1)))
         enddo
      enddo
c
      if(mod(n,100).eq.0) print *, 'Timestep=',n
c
c ...  absorbing bc (left and right)
c
      do ka=0,nz
         do ia=0,iabmax
            vx(ia,ka)=vx(ia,ka)*eponge(ia)
            vz(ia,ka)=vz(ia,ka)*eponge(ia)
            vx(nx-ia,ka)=vx(nx-ia,ka)*eponge(ia)
            vz(nx-1-ia,ka)=vz(nx-1-ia,ka)*eponge(ia)
         enddo
      enddo
c
c ... absorbing bc (top and bottom) 
c
      do ka=0,iabmax
         do ia=0,nx
            if(topbc.eq.'abbc') then 
               vx(ia,ka)=vx(ia,ka)*eponge(ka)
               vz(ia,ka)=vz(ia,ka)*eponge(ka)
            endif 
            vx(ia,nz-1-ka)=vx(ia,nz-1-ka)*eponge(ka)
            vz(ia,nz-ka)=vz(ia,nz-ka)*eponge(ka)
         enddo
      enddo
c
c ... free surface boundary condition 
c
      if(topbc.eq.'free') then 
         do i=2,nx-2
            vx(i,0)=vx(i,0)+bx0*(g1*(txx(i+1,0)-txx(i-2,0))+
     *         g2*(txx(i,0)-txx(i-1,0))+
     *         g1*(txz(i,2)+txz(i,1))+
     *         g2*(txz(i,1)-txz(i,0)))
         enddo
         do i=1,nx-2
            vz(i,0)=vz(i,0)+2*bz0*
     *        (g1*tzz(i,1)+g2*tzz(i,0))
            vz(i,1)=vz(i,1)+bz1*
     *        (g1*(tzz(i,2)+tzz(i,0))+
     *         g2*(tzz(i,1)-tzz(i,0))+
     *         g1*(txz(i+2,1)-txz(i-1,1))+
     *         g2*(txz(i+1,1)-txz(i,1)))
         enddo
      endif 
c 
c .... vertical force source
c
      if(isotype.eq.4) then
         addsou=0.5*soufac*source(n)/rho(ixs,izs)
         vx(ixs,izs)=vx(ixs,izs)+addsou
         vx(ixs,izs+1)=vx(ixs,izs+1)+addsou
      endif
c
c .... add vertical force source
c
      if(isotype.eq.3) then
         addsou=0.5*soufac*source(n)/rho(ixs-1,izs+1)
         vz(ixs-1,izs+1)=vz(ixs-1,izs+1)+addsou
         vz(ixs,izs+1)=vz(ixs,izs+1)+addsou
      endif
c      
c ===============================================================
c     Stress-strain relations for an isotropic elastic medium 
c ===============================================================
c
c      
      do k=2,nz-2
         do i=1,nx-2
c
c...  xx strain
c
            exx=g1*(vx(i+2,k)-vx(i-1,k))+g2*(vx(i+1,k)-vx(i,k))
c
c...  yy strains
c
            ezz=g1*(vz(i,k+2)-vz(i,k-1))+g2*(vz(i,k+1)-vz(i,k))
c
c...  xy strain
c
            if(k.eq.(nz-2))then
               exz=vx(i,k)-vx(i,k-1)+
     *              g1*(vz(i+1,k)-vz(i-2,k))+g2*(vz(i,k)-vz(i-1,k))
               if(i.eq.1) exz=vx(1,k)-vx(1,k-1)+vz(1,k)-vz(0,k)
            endif
            if(k.ge.2.and.k.le.(nz-3))then
               exz=g1*(vx(i,k+1)-vx(i,k-2))+g2*(vx(i,k)-vx(i,k-1))+
     *                 g1*(vz(i+1,k)-vz(i-2,k))+g2*(vz(i,k)-vz(i-1,k))
               if(i.eq.1) then
                  exz=g1*(vx(1,k+1)-vx(1,k-2))+g2*(vx(1,k)-vx(1,k-1))+
     *                vz(1,k)-vz(0,k)
               endif
            endif
c
c...  stresses
c
            txx(i,k)=txx(i,k) + rlamu(i,k)*exx+rla(i,k)*ezz
            tzz(i,k)=tzz(i,k) + rla(i,k)*exx+rlamu(i,k)*ezz
            txz(i,k)=txz(i,k) + rmuxz(i,k)*exz
         enddo
      enddo
c
c
c... free surface boundary conditions 
c
      if(topbc.eq.'free')then 
        do i=2,nx-2
           txz(i,0)=0.0
           txz(i,1)=txz(i,1)+
     +           rmuxz1*(g1*(vx(i,2)+vz(i-1,0)-vz(i,0)-vx(i,0))+
     +           g2*(vx(i,1)-vx(i,0))+
     +           g1*(vz(i+1,1)-vz(i-2,1))+g2*(vz(i,1)-vz(i-1,1)))
        enddo
        do i=0,nx
           tzz(i,0)=tzz(i,0)+
     +      rla0*(g1*(vx(i+2,0)-vx(i-1,0))+g2*(vx(i+1,0)-vx(i,0)))+
     +      rlamu0*(vz(i,1)-vz(i,0))
           tzz(i,1)=tzz(i,1)+
     +      rla1*(g1*(vx(i+2,1)-vx(i-1,1))+g2*(vx(i+1,1)-vx(i,1)))+
     +      rlamu1*(g1*(vz(i,3)-vz(i,0))+g2*(vz(i,2)-vz(i,1)))
        enddo
        do i=0,nx-1
           txx(i,0)=txx(i,0) +
     +      rlamu0*(g1*(vx(i+2,0)-vx(i-1,0))+g2*(vx(i+1,0)-vx(i,0)))+
     +      rla0*(vz(i,1)-vz(i,0))
           txx(i,1)=txx(i,1)+
     +      rlamu1*(g1*(vx(i+2,1)-vx(i-1,1))+g2*(vx(i+1,1)-vx(i,1)))+
     +      rla1*(g1*(vz(i,3)-vz(i,0))+g2*(vz(i,2)-vz(i,1)))
        enddo
      endif
c
c ... absorbing bc (left and right)
c
      do ka=0,nz
         do ia=0,iabmax
            txz(ia,ka)=txz(ia,ka)*eponge(ia)
            txx(ia,ka)=txx(ia,ka)*eponge(ia)
            tzz(ia,ka)=tzz(ia,ka)*eponge(ia)
            txz(nx-ia,ka)=txz(nx-ia,ka)*eponge(ia)
            txx(nx-1-ia,ka)=txx(nx-1-ia,ka)*eponge(ia)
            tzz(nx-1-ia,ka)=tzz(nx-1-ia,ka)*eponge(ia)
         enddo
      enddo
c
c ... absorbing bc (top and bottom)
c
      do ka=0,iabmax
         do ia=0,nx
            if(topbc .eq. 'abbc') then
               txz(ia,ka)=txz(ia,ka)*eponge(ka)
               txx(ia,ka)=txx(ia,ka)*eponge(ka)
               tzz(ia,ka)=tzz(ia,ka)*eponge(ka)
            endif
            txz(ia,nz-ka)=txz(ia,nz-ka)*eponge(ka)
            txx(ia,nz-1-ka)=txx(ia,nz-1-ka)*eponge(ka)
            tzz(ia,nz-1-ka)=tzz(ia,nz-1-ka)*eponge(ka)
         enddo
      enddo
c
c...  pressure sources
c
      if(isotype.eq.1) then
         addsou=soufac*source(n)*0.25
         txx(ixs-1,izs)=txx(ixs-1,izs)+addsou
         txx(ixs,izs)=txx(ixs,izs)+addsou
         txx(ixs-1,izs+1)=txx(ixs-1,izs+1)+addsou
         txx(ixs,izs+1)=txx(ixs,izs+1)+addsou
c
         tzz(ixs-1,izs)=tzz(ixs-1,izs)+addsou
         tzz(ixs,izs)=tzz(ixs,izs)+addsou
         tzz(ixs-1,izs+1)=tzz(ixs-1,izs+1)+addsou
         tzz(ixs,izs+1)=tzz(ixs,izs+1)+addsou
       endif
c
c...  shear sources
c
      if(isotype.eq.2) then
         addsou=soufac*source(n)*0.25
         txx(ixs-1,izs)=txx(ixs-1,izs)+addsou
         txx(ixs,izs)=txx(ixs,izs)+addsou
         txx(ixs-1,izs+1)=txx(ixs-1,izs+1)+addsou
         txx(ixs,izs+1)=txx(ixs,izs+1)+addsou
c
         tzz(ixs-1,izs)=tzz(ixs-1,izs)-addsou
         tzz(ixs,izs)=tzz(ixs,izs)-addsou
         tzz(ixs-1,izs+1)=tzz(ixs-1,izs+1)-addsou
         tzz(ixs,izs+1)=tzz(ixs,izs+1)-addsou
       endif
c
       if(mod(n,idsnap).eq.0.and.isnap.lt.nsnap) then
          do iz=0,nz-1,2
             do ix=0,nx-1,2
                aaa = tzz(ix,iz)+txx(ix,iz)
                write(38+isnap,*) aaa
             enddo
          enddo
          isnap = isnap+1
       endif

       enddo
c
       stop
       end
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+     model calculates coefficients for finite-differences eqns   +
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
       subroutine model(rlamu,rmu,rmuxz,rho,bx,bz,
     *                 dtdx,nx,nz,bx0,bz0,bz1,rmuxz1,
     *                 rmu0,rmu1,rlamu0,rlamu1)
c

c
       parameter (nmat=3, eps=1.e-6)
       dimension rlamu(0:nx,0:nz), rmu(0:nx,0:nz) 
       dimension rho(0:nx,0:nz), rmuxz(0:nx,0:nz)
       dimension bx(0:nx,0:nz), bz(0:nx,0:nz)
       dimension alpha(nmat), beta(nmat), dens(nmat)
c
       data alpha(1),beta(1),dens(1)/1500.0,   0.0,1000.0/
       data alpha(2),beta(2),dens(2)/2200.0,1150.0,2050.0/
       data alpha(3),beta(3),dens(3)/2700.0,1550.0,2750.0/
c
c
c ...  begin construction of the geological
       do ix=0,nx
          do iz=0,nz
             rlamu(ix,iz) = alpha(1)*alpha(1)*dens(1)
             rmu(ix,iz)   = beta(1)*beta(1)*dens(1)
             rho(ix,iz)   = dens(1)
          enddo
       enddo
c
       do ix=0,nx
          izxa = 100
          do iz=izxa+1,nz
             rlamu(ix,iz) = alpha(2)*alpha(2)*dens(2)
             rmu(ix,iz)   = beta(2)*beta(2)*dens(2)
             rho(ix,iz)   = dens(2)
          enddo
          if(ix.le.(nx/2)) then
             izxc = 300
             do iz=izxc+1,nz
                rlamu(ix,iz) = alpha(3)*alpha(3)*dens(3)
                rmu(ix,iz)   = beta(3)*beta(3)*dens(3)
                rho(ix,iz)   = dens(3)
             enddo
          endif
          if(ix.gt.(nx/2)) then
             izxc = 400
             do iz=izxc+1,nz
                rlamu(ix,iz) = alpha(3)*alpha(3)*dens(3)
                rmu(ix,iz)   = beta(3)*beta(3)*dens(3)
                rho(ix,iz)   = dens(3)
             enddo
          endif
       enddo
c ...  end construction of the geology 
c
       do i=0,nx
          do k=0,nz
             im=max0(0,i-1)
             km=max0(0,k-1)
             ii=min0(i,nx-1)
             kk=min0(k,nz-1)
c
c ... b_x
             bx(i,k) = 0.5*dtdx*((1.0/rho(im,kk))+(1.0/rho(ii,kk)))
c
c ... b_z 
             rhobz  =.5*(rho(ii,km)+rho(ii,kk))
             bz(i,k) = 0.5*dtdx*((1.0/rho(ii,km))+(1.0/rho(ii,kk)))
c 
c .... mu_xz
             akan = 1.0/(rmu(ii,kk)+eps)
             akan = akan + (1.0/(rmu(im,kk)+eps))
             akan = akan + (1.0/(rmu(ii,km)+eps))
             akan = akan + (1.0/(rmu(im,km)+eps))
             rmuxz(i,k) = (dtdx*4.0)/akan
             if(rmuxz(i,k).lt.1.0) rmuxz(i,k) = 0.0
          enddo
       enddo
c
       do i=0,nx-1
          do k=0,nz-1
c
c.... lambda
             rmu(i,k)   = dtdx*(rlamu(i,k)-(2.0*rmu(i,k)))
c
c ... lambda + 2 mu
             rlamu(i,k) = dtdx*rlamu(i,k)

          enddo
       enddo
c
c ...  fd coefficients for free surface
       bx0    = bx(0,0)
       bz0    = bz(0,0)
       bz1    = bz(1,0)
       rlamu0 = rlamu(0,0)
       rlamu1 = rlamu(1,0)
       rmu0   = rmu(0,0)
       rmu1   = rmu(1,0)
       rmuxz1 = rmuxz(1,0) 
c
c .... some adition conditions to reduce the ``if'' statements in fdm code
c .... vx 
       do i=0,nx
          bx(i,0) = 0.0
          bx(i,nz-1) = 0.0
          bx(i,nz) = 0.0
       enddo
       do k=0,nz
          bx(0,k) = 0.0
          bx(1,k) = 0.0
          bx(nx-1,k) = 0.0
          bx(nx,k) = 0.0
       enddo
c .... vz 
       do i=0,nx
          bz(i,0) = 0.0
          bz(i,1) = 0.0
          bz(i,nz-1) = 0.0
          bz(i,nz) = 0.0
       enddo
c
       do k=0,nz
          bz(0,k) = 0.0
          bz(nx-1,k) = 0.0
          bz(nx,k) = 0.0
       enddo
c ... txx and tzz
       do i=0,nx
          rlamu(i,0) = 0.0
          rlamu(i,1) = 0.0
          rlamu(i,nz-2) = 0.0
          rlamu(i,nz-1) = 0.0
          rlamu(i,nz) = 0.0
          rmu(i,0) = 0.0
          rmu(i,1) = 0.0
          rmu(i,nz-2) = 0.0
          rmu(i,nz-1) = 0.0
          rmu(i,nz) = 0.0
       enddo
       do k=0,nz
          rlamu(0,k)    = 0.0
          rlamu(nx-1,k) = 0.0
          rlamu(nx,k)   = 0.0
          rmu(0,k)    = 0.0
          rmu(nx-1,k) = 0.0
          rmu(nx,k)   = 0.0
       enddo
c ... txz
       do i=0,nx
          rmuxz(i,0) = 0.0
          rmuxz(i,1) = 0.0
          rmuxz(i,nz-1) = 0.0
          rmuxz(i,nz) = 0.0
       enddo
       do k=0,nz
          rmuxz(0,k)    = 0.0
          rmuxz(nx,k)   = 0.0
       enddo
c
       return
       end
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+    wavelet calculates the source signature                      +
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
       subroutine wavelet(nt,dt,wav)
c
       dimension wav(*)
C
       do i=1,nt
          wav(i) = 0.0
       enddo
C
       sfreq=25.
       lwav   = 1 + int(1.55/(sfreq*dt))
c
       do i=1,lwav
          wav(i) = ricker(float(i-1)*dt,sfreq)
       enddo
c
       return 
       end
C
       function ricker(time,sfreq)
       pix = 4.05367
       g0  = 0.5
       g1  = 0.565
       g2  = 0.105
       t0  = time*sfreq*pix 
       ricker = 0.0
       ricker = g0*cos(t0)-g1*cos(2*t0)+g2*cos(3*t0)
       return
       end
c