trisubs.f

linux环境下利用fortran语言开发

                zs2 = real(z(4))
        else
                zp2 = real(z(4))
        endif
      else
        error = 2
        return
      endif
      
*   finds the correct root according to mode and scattering type

*   QP slowness is real: test for scattering direction

        do444 i=1,3
                p1(i) = vtan(i)+real(zp1)*normal(i)
                p2(i) = vtan(i)+real(zp2)*normal(i)
444        continue
        call groupvel(p1,axis,2,elast,w1)
        call groupvel(p2,axis,2,elast,w2)
        dot1 = 0.e0
        dot2 = 0.e0
        do555 i=1,3
                dot1 = dot1 + w1(i)*normal(i)
                dot2 = dot2 + w2(i)*normal(i)
555        continue

        if(iscat .EQ. 0) then
                if(dot1 .GT. 0.e0) then
                        pnorm = zp1
                else if(dot2 .GT. 0.e0) then
                        pnorm = zp2
                else
*               no reflected wave
                        error = 3
                        return
                endif
        else
                if(dot1 .GT. 0.e0) then
                        pnorm = zp2
                else if(dot2 .GT. 0.e0) then
                        pnorm = zp1
                else
*               no transmitted wave
                        error = 3
                        return
                endif
        endif

      do200 i=1,3
        pscat2(i) = vtan(i)+pnorm*normal(i)
200   continue

*   QS slowness is real: test for scattering direction

        do666 i=1,3
                p1(i) = vtan(i)+real(zs1)*normal(i)
                p2(i) = vtan(i)+real(zs2)*normal(i)
666        continue

        call groupvel(p1,axis,3,elast,w1)
        call groupvel(p2,axis,3,elast,w2)

        dot1 = 0.e0
        dot2 = 0.e0

        do777 i=1,3
                dot1 = dot1 + w1(i)*normal(i)
                dot2 = dot2 + w2(i)*normal(i)
777        continue

        if(iscat .EQ. 0) then
                if(dot1 .GT. 0.e0) then
                        pnorm = zs1
                else if(dot2 .GT. 0.e0) then
                        pnorm = zs2
                else
*               no reflected wave
                        error = 3
                        return
                endif
        else
                if(dot1 .GT. 0.e0) then
                        pnorm = zs2
                else if(dot2 .GT. 0.e0) then
                        pnorm = zs1
                else
*               no transmitted wave
                        error = 3
                        return
                endif
        endif

        do300 i=1,3
          pscat3(i) = vtan(i)+pnorm*normal(i)
300     continue

        return

      endif

      if(ireal .EQ. 2) then

*   check that the real roots are associated with the QS wave

        ip = 0
        is = 0

        pa2 = 2.e0*ptan*na*ta*real(z(1))
        pr2 = pa2
        pa2 = pa2+ptan2*ta2+real(z(1))*real(z(1))*na2
        pr2 = ptan2*ta1-pr2+real(z(1))*real(z(1))*na1

        alpha = (A+L)*pr2+(C+L)*pa2
        beta = (A-L)*pr2-(C-L)*pa2
        beta = beta*beta + 4.e0*pr2*pa2*(L+F)*(L+F)
        beta = csqrt(beta)

        del1 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha+beta)
        del2 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha-beta)

        if(cabs(del1) .LE. cabs(del2)) then
          zp1 = real(z(1))
          ip = ip+1               
        else
          zs1 = real(z(1))
          is = is+1
        endif

        pa2 = 2.e0*ptan*na*ta*real(z(2))
        pr2 = pa2
        pa2 = pa2+ptan2*ta2+real(z(2))*real(z(2))*na2
        pr2 = ptan2*ta1-pr2+real(z(2))*real(z(2))*na1

        alpha = (A+L)*pr2+(C+L)*pa2
        beta = (A-L)*pr2-(C-L)*pa2
        beta = beta*beta + 4.e0*pr2*pa2*(L+F)*(L+F)
        beta = csqrt(beta)

          del1 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha+beta)
          del2 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha-beta)

        if(cabs(del1) .LE. cabs(del2)) then
          if(ip .EQ. 1) then
            zp2 = real(z(2))
          else
            zp1 = real(z(2))
          endif
          ip = ip+1
        else
          if(is .EQ. 1) then
            zs2 = real(z(2))
          else
            zs1 = real(z(2))
          endif
          is = is+1
        endif

        if(is .NE. 2) then

          error=2
          return

        else

*   QS slowness is real: test for scattering direction

          do1666 i=1,3
            p1(i) = vtan(i)+real(zs1)*normal(i)
            p2(i) = vtan(i)+real(zs2)*normal(i)
1666      continue

          call groupvel(p1,axis,3,elast,w1)
          call groupvel(p2,axis,3,elast,w2)

          dot1 = 0.e0
          dot2 = 0.e0

          do1777 i=1,3
            dot1 = dot1 + w1(i)*normal(i)
            dot2 = dot2 + w2(i)*normal(i)
1777      continue

        if(iscat .EQ. 0) then
          if(dot1 .GT. 0.e0) then
            pnorm = zs1
          else if(dot2 .GT. 0.e0) then
            pnorm = zs2
          else
*           no reflected wave
            error = 3
            return
          endif
        else
          if(dot1 .GT. 0.e0) then
            pnorm = zs2
          else if(dot2 .GT. 0.e0) then
            pnorm = zs1
          else
*           no transmitted wave
            error = 3
            return
          endif
        endif
        
      endif

      do1300 i=1,3
        pscat3(i) = vtan(i)+pnorm*normal(i)
1300  continue

*   QP scattering is complex: test for decay
*   bearing in mind that we are considering positive
*   frequencies(!).

        if(iscat .EQ. 0) then
                if(aimag(z(3)) .GT. 0.) then
                        pnorm = z(3)
                else
                        pnorm = z(4)
                endif
        else
                if(aimag(z(3)) .GT. 0.) then
                        pnorm = z(4)
                else
                        pnorm = z(3)
                endif
        endif

        do1200 i=1,3
          pscat2(i) = vtan(i)+pnorm*normal(i)
1200    continue

        return
 
      endif

      if(ireal .EQ. 0) then

*   QP and QS modes are both complex: reshuffle the roots 
*   so that the two with positive imaginary parts are in front ranks

        do1111 i=1,4
          if(aimag(z(i)) .GT. 0.e0) then
            zz = z(1)
            z(1) = z(i)
            z(i) = zz
            goto 1222
          endif
1111    continue

1222    do1333 i=2,4
           if(aimag(z(i)) .GT. 0.e0) then
             zz = z(2)
             z(2) = z(i)
             z(i) = zz
           endif
1333    continue          

        if(iscat .EQ. 0.e0) then

*   The evanescent slownesses have positive imaginary part.
*   Now discriminate between QP and QS:

        
          pa2 = 2.e0*ptan*na*ta*z(1)
          pr2 = pa2
          pa2 = pa2+ptan2*ta2+z(1)*z(1)*na2
          pr2 = ptan2*ta1-pr2+z(1)*z(1)*na1
    
          alpha = (A+L)*pr2+(C+L)*pa2
          beta = (A-L)*pr2-(C-L)*pa2
          beta = beta*beta + 4.e0*pr2*pa2*(L+F)*(L+F)
          beta = csqrt(beta)

          del1 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha+beta)
          del2 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha-beta)

          if(cabs(del1) .LE. cabs(del2)) then
            zp1 = z(1)
            zs1 = z(2)
          else
            zp1 = z(2)
            zs1 = z(1)
          endif

          do2300 i=1,3
            pscat3(i) = vtan(i)+zs1*normal(i)
2300      continue

          do2200 i=1,3
              pscat2(i) = vtan(i)+zp1*normal(i)
2200      continue

          return

        else

*   The evanescent slownesses have negative imaginary part.
*   Now discriminate between QP and QS:

        
          pa2 = 2.e0*ptan*na*ta*z(3)
          pr2 = pa2
          pa2 = pa2+ptan2*ta2+z(3)*z(3)*na2
          pr2 = ptan2*ta1-pr2+z(3)*z(3)*na1
    
          alpha = (A+L)*pr2+(C+L)*pa2
          beta = (A-L)*pr2-(C-L)*pa2
          beta = beta*beta + 4.e0*pr2*pa2*(L+F)*(L+F)
          beta = csqrt(beta)

          del1 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha+beta)
          del2 = cmplx(1.e0,0.e0) - 2.e0*rho/(alpha-beta)

          if(cabs(del1) .LE. cabs(del2)) then
            zp1 = z(3)
            zs1 = z(4)
          else
            zp1 = z(4)
            zs1 = z(3)
          endif

          do3300 i=1,3
            pscat3(i) = vtan(i)+zs1*normal(i)
3300      continue

          do3200 i=1,3
              pscat2(i) = vtan(i)+zp1*normal(i)
3200      continue

          return

        endif

      endif

      return
 
      end



      subroutine cshotprof()

*
*   generates a shotprofile over a tabular transversely 
*   isotropic medium.
*
*   GLOBAL VARIABLES:
*
*       block "model"   elastic and geometric model parameters
*       block "array"   source-receiver configuration
*       block "path"    raypath parameters
*       block "output"  graphics/verbose parameters
*       block "io"      input/output logical unit numbers
*
*   OUTPUT FILES:
*
*       uh1.dat         first horizontal displacement
*       uh2.dat         second horizontal displacement
*       uv.dat          vertical displacement
*

*     initializing a few parameters:
*        maxtrace, maxsample: max numbers of traces and time samples
*                              on a shot record
*        maxrad, maxang: maxrad*(maxang-2) is the maximum number of 
*                        points in the interpolation map. maxang must 
*                        be more than 3 and less than 100000.
*        difdip: initial increment for take-off dip angle


      integer maxtrace, maxsample
      integer maxrad, maxang
      real difdip

      parameter(maxtrace=200, maxsample=2048)
      parameter(maxang=43, maxrad=1000)
      parameter(difdip=1.e0)

*      layer boundary color
      parameter(lcolor=4)

*     variable declaration

      real param(900), geom(0:99), dt, scale, xmin, ymin
      real xs, ys, cazi, cdip, cg1, dx, cr, cdg, ail, acl, xr, yr
      real wave(255), hilwav(255), fl, fh
      integer nlayer, nsegment(200), mode(200,100), tsamp, trace
      integer interface(200,0:100), triso(100), ng, npath, verbose
      integer nsweep, ntaper, leng
      integer stdin, stderr, stdout, temp1, temp2, temp3
      common /model/ param, geom, nlayer, triso
      common /path/ mode, interface, nsegment, npath
      common /array/ xs,ys,cazi,cdip,cr,cdg,xr,yr,cg1,
     &               ail,acl,dx,ng,tsamp,dt
      common /output/ verbose, trace, scale, xmin, ymin
      common /wavelet/ nsweep, ntaper, leng, fl, fh, wave, hilwav
      common /io/ stdin, stdout, stderr, temp1, temp2, temp3

      complex disp(3), tempil, tempcl
      real uil(maxtrace,maxsample), ucl(maxtrace,maxsample)
      real uv(maxtrace,maxsample), sg2(maxrad), dg2
      real deg2rad, dxx, dyy, cail, sail, cacl, sacl, xp, yp, zp
      real x, y, g1, g2, t, xx(maxang,maxrad), yy(maxang,maxrad)
      real xmax, ymax, sg1(maxang), ddip
      integer ipath, ig, error, time, j, imax, imin, ii, i, ng2
      integer maxa, maxr

      call getoutput()
      call getmodel()
      call getarray()
      call getpath()

*     checks record parameters against max array bounds

      write(stderr,"(/,' *** RUNNING ...',/)")

      if(ng .GT. maxtrace) then
	    stop ' chshotprof: ng too large'
      endif

      if(tsamp .GT. maxsample) then
	    stop ' chshotprof: tsamp too large'
      endif

      deg2rad = 1.7453292e-2

      dxx = dx*cos(cg1*deg2rad)
      dyy = dx*sin(cg1*deg2rad)
      cail = cos(ail*deg2rad)
      sail = sin(ail*deg2rad)
      cacl = cos(acl*deg2rad)
      sacl = sin(acl*deg2rad)

      if(trace .EQ. 1) then

*        compute some graphics parameters
         xmin = amin1(xr,xr+(ng-1)*dxx,xs)
         ymin = amin1(yr,yr+(ng-1)*dyy,ys)
         xmax = amax1(xr,xr+(ng-1)*dxx,xs)
         ymax = amax1(yr,yr+(ng-1)*dyy,ys)

*        trace interfaces

         do1000 j=1,nlayer

            write(stdout,*) 5, lcolor

            xp = 0.
            yp = 0.
            zp = geom(j-1)
            write(stdout,*) xp,yp,zp

            xp = xmax-xmin
            write(stdout,*) xp,yp,zp

            yp = ymax-ymin
            write(stdout,*) xp,yp,zp

            xp = 0.
            write(stdout,*) xp,yp,zp

            yp = 0.
            write(stdout,*) xp,yp,zp

1000     continue

      endif

      do10 ipath=1,npath

         if(verbose .EQ. 1) 
     &      write(stderr,"(/,' *** RAYPATH #',i3,':',/)") ipath

         if(verbose .EQ. 1) write(stderr,"('   map started...')")

         maxr = maxrad
         maxa = maxang
         ddip = difdip

         call shootmap(xx,yy,sg1,sg2,ipath,ng2,maxr,maxa,ddip)

         if(verbose .EQ. 1) 
     &      write(stderr,"('              ...map completed',/)")