开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 并行计算 > 查看源码
higdon2d.f
资源名称：wave_forward.rar　[点击查看]
上传用户：ls4004p
上传日期：2007-08-05
资源大小：2314k
文件大小：12k
源码类别：

并行计算

开发平台：

Matlab

higdon2d.f：源码内容
  1. C==========================================================================
  2. C    Routinen InitBound , higx1, higxfs1, higx2, higxfs2,
  3. C                       , higz1, higzfs1, higz2, higzfs2
  4. C    07.04.1995  by Joachim Falk
  5. C    implementation of  2D-Fields
  6. C==========================================================================
  7.       subroutine initbound(lamu,mu,rho,dx,dz,
  8.      &                     tvp,tvs,bvp,bvs,lvp,lvs,rvp,rvs,
  9.      &                     nx,nz,dt,fs)
  10.       implicit none
  12.       include 'param.f'
  14.       integer fs, nx , nz, r, i
  15.       real mu(MAXCOLS,MAXROWS),lamu(MAXCOLS,MAXROWS),
  16.      &     rho(MAXCOLS,MAXROWS),
  17.      &     dx, dz,
  18.      &     tvp(MAXCOLS),tvs(MAXCOLS),bvp(MAXCOLS),bvs(MAXCOLS),
  19.      &     lvp(MAXROWS),lvs(MAXROWS),rvp(MAXROWS),rvs(MAXROWS),
  20.      &     dt,hlp, hlp2
  21.       
  22.       hlp2 = sqrt(0.5)
  24.       do r=1, nz
  25.           IF (mu(1,r).gt.1.e5) THEN 
  26.              hlp = sqrt(mu(1,r) / rho(1,r)) * dt / dx
  27.              lvp(r) = (1.-hlp) / (1.+hlp)
  28.              lvs(r) = (hlp2-hlp) / (hlp2+hlp)
  29. C             hlp = sqrt(lamu(1,r) / rho(1,r)) * dt / dx
  30. C             lvp(r) = (1.-hlp) / (1.+hlp)
  31. C             hlp = sqrt(mu(1,r) / rho(1,r)) * dt / dx
  32. C             lvs(r) = (1.-hlp) / (1.+hlp)
  33.           ELSE
  34.              hlp = sqrt(lamu(1,r) / rho(1,r)) * dt / dx
  35.              lvp(r) = (1.-hlp) / (1.+hlp)
  36.              lvs(r) = (hlp2-hlp) / (hlp2+hlp)
  37.           ENDIF
  38.           hlp = lvp(r)
  39.           lvp(r) = hlp + lvs(r)
  40.           lvs(r) = hlp * lvs(r)
  42.           IF (mu(nx,r).gt.1.e5) THEN 
  43.              hlp = sqrt(mu(nx,r) / rho(nx,r)) * dt / dx
  44.              rvp(r) = (1.-hlp) / (1.+hlp)
  45.              rvs(r) = (hlp2-hlp) / (hlp2+hlp)
  46. C             hlp = sqrt(lamu(nx,r) / rho(nx,r)) * dt / dx
  47. C             rvp(r) = (1.-hlp) / (1.+hlp)
  48. C             hlp = sqrt(mu(nx,r) / rho(nx,r)) * dt / dx
  49. C             rvs(r) = (1.-hlp) / (1.+hlp)
  50.           ELSE
  51.              hlp = sqrt(lamu(nx,r) / rho(nx,r)) * dt / dx
  52.              rvp(r) = (1.-hlp) / (1.+hlp)
  53.              rvs(r) = (hlp2-hlp) / (hlp2+hlp)
  54.           ENDIF
  55.           hlp = rvp(r)
  56.           rvp(r) = hlp + rvs(r)
  57.           rvs(r) = hlp * rvs(r)
  58.         enddo
  60.         do i = 1 , nx
  61.            IF (mu(i,nz).gt.1.e5) THEN 
  62.               hlp = sqrt(mu(i,nz) / rho(i,nz)) * dt / dz
  63.               bvp(i) = (1.-hlp) / (1.+hlp)
  64.               hlp = sqrt(mu(i,nz) / rho(i,nz)) * dt / dz
  65.               bvs(i) = (hlp2-hlp) / (hlp2+hlp)
  66. C              hlp = sqrt(lamu(i,nz) / rho(i,nz)) * dt / dz
  67. C              bvp(i) = (1.-hlp) / (1.+hlp)
  68. C              hlp = sqrt(mu(i,nz) / rho(i,nz)) * dt / dz
  69. C              bvs(i) = (1.-hlp) / (1.+hlp)
  70.            ELSE
  71.               hlp = sqrt(lamu(i,nz) / rho(i,nz)) * dt / dz
  72.               bvp(i) = (1.-hlp) / (1.+hlp)
  73.               bvs(i) = (hlp2-hlp) / (hlp2+hlp)
  74.            ENDIF
  75.            hlp = bvp(i)
  76.            bvp(i) = hlp + bvs(i)
  77.            bvs(i) = hlp * bvs(i)
  78.         enddo
  80.         if (fs.ne.1) then   
  81.           do i = 1 , nx
  82.              IF (mu(i,1).gt.1.e5) THEN 
  83.                 hlp = sqrt(mu(i,1) / rho(i,1)) * dt / dz
  84.                 tvp(i) = (1.-hlp) / (1.+hlp)
  85.                 tvs(i) = (hlp2-hlp) / (hlp2+hlp)
  86. C                hlp = sqrt(lamu(i,1) / rho(i,1)) * dt / dz
  87. C                tvp(i) = (1.-hlp) / (1.+hlp)
  88. C                hlp = sqrt(mu(i,1) / rho(i,1)) * dt / dz
  89. C                tvs(i) = (1.-hlp) / (1.+hlp)
  90.              ELSE
  91.                 hlp = sqrt(lamu(i,1) / rho(i,1)) * dt / dz
  92.                 tvp(i) = (1.-hlp) / (1.+hlp)
  93.                 tvs(i) = (hlp2-hlp) / (hlp2+hlp)
  94.              ENDIF
  95.             hlp = tvp(i)
  96.             tvp(i) = hlp + tvs(i)
  97.             tvs(i) = hlp * tvs(i)
  98.           enddo
  99.         endif
  100.       return
  101.       end
  103. C==========================================================================
  104.       subroutine HigX2(Ux,Uxo,tvp,tvs,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
  105.       implicit none
  107.       include 'param.f'
  109.       integer in, nx, nz, r, c
  110.       real Ux(MAXCOLS,MAXROWS),Uxo(MAXCOLS,MAXROWS),
  111.      &     tvp(MAXCOLS),tvs(MAXCOLS),bvp(MAXCOLS),bvs(MAXCOLS),
  112.      &     lvp(MAXROWS),lvs(MAXROWS),rvp(MAXROWS),rvs(MAXROWS)
  113.           
  114.       do c = 2 , nx-2
  115.          Ux(c,1) = Ux(c,1) - tvp(c) * Ux(c,2) - tvs(c) * Ux(c,3)
  116.       enddo
  118.       in = nx-1                 ! Spalte nx-1 Rand rechts (2.Durchlauf)
  119.       do r = 2, nz-1            ! Schleife ueber re und li Rand (2.Durchlauf)
  120.          Ux(1,r) = Ux(1,r) - lvp(r) * Ux(2,r) - lvs(r) * Ux(3,r)
  121.          Ux(in,r) = Ux(in,r) - rvp(r)*Ux(in-1,r)-rvs(r)*Ux(in-2,r)
  122.       enddo
  123.       
  124.       do c=2, nx-2
  125.          Ux(c,nz) = Ux(c,nz) - bvp(c)*Ux(c,nz-1) - bvs(c)*Ux(c,nz-2)
  126.       enddo
  127.       return
  128.       end
  129. C==========================================================================
  130.       subroutine HigZ1(Uz,Uzo,tvp,tvs,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
  131.       implicit none
  133.       include 'param.f'
  135.       integer in, nx, nz, r, c
  136.       real Uz(MAXCOLS,MAXROWS),Uzo(MAXCOLS,MAXROWS),
  137.      &     tvp(MAXCOLS),tvs(MAXCOLS),bvp(MAXCOLS),bvs(MAXCOLS),
  138.      &     lvp(MAXROWS),lvs(MAXROWS),rvp(MAXROWS),rvs(MAXROWS),hlp
  139.         
  140.       do c=2, nx-1   
  141.          hlp = Uz(c,1)
  142.          Uz(c,1) = 2.* Uz(c,2) - Uzo(c,3)
  143.      &        - tvp(c) * (Uzo(c,2) - Uz(c,3) - Uz(c,1))
  144.      &        - tvs(c) * (Uzo(c,1) - 2.* Uz(c,2))
  145.          Uzo(c,1) = hlp
  146.       enddo
  147.       in = nx
  148.       do r=2, nz-2              ! Rand links (erster Durchlauf)
  149.          hlp = Uz(1,r)
  150.          Uz(1,r) = 2. * Uz(2,r) - Uzo(3,r)
  151.      &        - lvp(r) * (Uzo(2,r) - Uz(3,r) - Uz(1,r))
  152.      &        - lvs(r) * (Uzo(1,r) - 2.*Uz(2,r))
  153.          Uzo(1,r) = hlp
  154.          
  155.          hlp = Uz(in,r)
  156.          Uz(in,r) = 2. * Uz(in-1,r) - Uzo(in-2,r)
  157.      &        - rvp(r) * (Uzo(in-1,r) - Uz(in-2,r) - Uz(in,r))
  158.      &        - rvs(r) * (Uzo(in,r) - 2.*Uz(in-1,r))
  159.          Uzo(in,r) = hlp
  160.       enddo
  161.       
  162.       in = nz-1
  163.       do c=2, nx-1              ! Zeile nz-1
  164.          hlp = Uz(c,in)
  165.          Uz(c,in) = 2.* Uz(c,in-1) - Uzo(c,in-2)
  166.      &        - bvp(c) * (Uzo(c,in-1) - Uz(c,in-2) - Uz(c,in))
  167.      &        - bvs(c) * (Uzo(c,in) - 2.* Uz(c,in-1))
  168.          Uzo(c,in) = hlp
  169.       enddo
  170.       return
  171.       end
  172. C==========================================================================
  173.       subroutine HigZ2(Uz,Uzo,tvp,tvs,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
  174.       implicit none
  176.       include 'param.f'
  178.       integer in, nx, nz, r, c
  179.       real Uz(MAXCOLS,MAXROWS),Uzo(MAXCOLS,MAXROWS),
  180.      &     tvp(MAXCOLS),tvs(MAXCOLS),bvp(MAXCOLS),bvs(MAXCOLS),
  181.      &     lvp(MAXROWS),lvs(MAXROWS),rvp(MAXROWS),rvs(MAXROWS)
  182.       do c=2, nx-1              ! Zeile 1
  183.          Uz(c,1) = Uz(c,1) - tvp(c) * Uz(c,2) - tvs(c) * Uz(c,3)
  184.       enddo
  185.       in = nx
  186.       do r= 2 , nz-2   
  187.          Uz(1,r) = Uz(1,r) - lvp(r) * Uz(2,r) - lvs(r) * Uz(3,r)
  188.          Uz(in,r) = Uz(in,r)-rvp(r) * Uz(in-1,r)-rvs(r)*Uz(in-2,r)
  189.       enddo
  190.       in = nz-1
  191.       do c=2, nx-1              ! Zeile nz-1
  192.          Uz(c,in) = Uz(c,in)-bvp(c)*Uz(c,in-1)-bvs(c)*Uz(c,in-2)
  193.       enddo
  194.       return
  195.       end
  196. C============================================================================
  197.       subroutine HigX1(Ux,Uxo,tvp,tvs,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
  198.       implicit none
  200.       include 'param.f'
  202.       integer in, nx, nz, r, c
  203.       real Ux(MAXCOLS,MAXROWS),Uxo(MAXCOLS,MAXROWS),
  204.      &     tvp(MAXCOLS),tvs(MAXCOLS),bvp(MAXCOLS),bvs(MAXCOLS),
  205.      &     lvp(MAXROWS),lvs(MAXROWS),rvp(MAXROWS),rvs(MAXROWS),hlp
  206.           
  207.           do c=2, nx-2
  208.              hlp = Ux(c,1)
  209.              Ux(c,1) = 2.* Ux(c,2) - Uxo(c,3)
  210.      &             - tvp(c) * (Uxo(c,2) - Ux(c,3) - Ux(c,1))
  211.      &             - tvs(c) * (Uxo(c,1) - 2.*Ux(c,2))
  212.              Uxo(c,1) = hlp
  213.           enddo
  214.                                                                    
  215.           in = nx-1
  216.           do r=2 , nz-1
  217.             hlp = Ux(1,r)
  218.             Ux(1,r) = 2.* Ux(2,r) - Uxo(3,r)
  219.      &             - lvp(r) * (Uxo(2,r) - Ux(3,r) - Ux(1,r))
  220.      &             - lvs(r) * (Uxo(1,r) - 2.* Ux(2,r))
  221.             Uxo(1,r) = hlp
  223.             hlp = Ux(in,r)
  224.             Ux(in,r) = 2.* Ux(in-1,r) - Uxo(in-2,r)
  225.      &             - rvp(r)  * (Uxo(in-1,r) - Ux(in-2,r) - Ux(in,r))
  226.      &             - rvs(r)  * (Uxo(in,r) - 2.* Ux(in-1,r))
  227.             Uxo(in,r) = hlp
  228.           enddo
  229.           
  230.           do c=2, nx-2          !  Rand unten (erster Durchlauf)
  231.              hlp = Ux(c,nz)
  232.              Ux(c,nz) = 2.* Ux(c,nz-1) - Uxo(c,nz-2)
  233.      &             - bvp(c) * (Uxo(c,nz-1) - Ux(c,nz-2) - Ux(c,nz))
  234.      &             - bvs(c) * (Uxo(c,nz) - 2.*Ux(c,nz-1))
  235.              Uxo(c,nz) = hlp
  236.           enddo
  237.         return  
  238.         end
  239. C=========================================================================
  240.       subroutine HigfsX1(Ux,Uxo,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
  241.       implicit none
  243.       include 'param.f'
  245.       integer in, nx, nz, r, c
  246.       real Ux(MAXCOLS,MAXROWS),Uxo(MAXCOLS,MAXROWS),
  247.      &     bvp(MAXCOLS),bvs(MAXCOLS),
  248.      &     lvp(MAXROWS),lvs(MAXROWS),rvp(MAXROWS),rvs(MAXROWS),hlp
  249.          
  250.          in = nx-1
  251.          do r=1 , nz-1
  252.             hlp = Ux(1,r)
  253.             Ux(1,r) = 2.* Ux(2,r) - Uxo(3,r)
  254.      &             - lvp(r) * (Uxo(2,r) - Ux(3,r) - Ux(1,r))
  255.      &             - lvs(r) * (Uxo(1,r) - 2.* Ux(2,r))
  256.             Uxo(1,r) = hlp
  257.             
  258.             hlp = Ux(in,r)
  259.             Ux(in,r) = 2.* Ux(in-1,r) - Uxo(in-2,r)
  260.      &             - rvp(r)  * (Uxo(in-1,r) - Ux(in-2,r) - Ux(in,r))
  261.      &             - rvs(r)  * (Uxo(in,r) - 2.* Ux(in-1,r))
  262.             Uxo(in,r) = hlp
  263.          enddo
  264.           
  265.           in = nz
  266.           do c=2, nx-2 !  Rand unten (erster Durchlauf)
  267.              hlp = Ux(c,in)
  268.              Ux(c,in) = 2.* Ux(c,in-1) - Uxo(c,in-2)
  269.      &             - bvp(c) * (Uxo(c,in-1) - Ux(c,in-2) - Ux(c,in))
  270.      &             - bvs(c) * (Uxo(c,in) - 2.*Ux(c,in-1))
  271.              Uxo(c,in) = hlp
  272.           enddo
  273.         return  
  274.         end
  275. C============================================================================
  276.       subroutine HigfsX2(Ux,Uxo,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
  277.       implicit none
  279.       include 'param.f'
  281.       integer in, nx, nz, r, c
  282.       real Ux(MAXCOLS,MAXROWS),Uxo(MAXCOLS,MAXROWS),
  283.      &     bvp(MAXCOLS),bvs(MAXCOLS),
  284.      &     lvp(MAXROWS),lvs(MAXROWS),rvp(MAXROWS),rvs(MAXROWS)
  285.          
  286.       in = nx-1
  287.       do r= 1 , nz-1            ! Schleife ueber re und li Rand (2.Durchlauf)
  288.          Ux(1,r)  = Ux(1,r) - lvp(r) * Ux(2,r) - lvs(r) * Ux(3,r)
  289.          Ux(in,r) = Ux(in,r)- rvp(r)*Ux(in-1,r) - rvs(r)*Ux(in-2,r)
  290.       enddo
  291.       do c=2, nx-2
  292.          Ux(c,nz) = Ux(c,nz)-bvp(c)*Ux(c,nz-1) - bvs(c)*Ux(c,nz-2)
  293.       enddo
  294.       return
  295.       end
  296. C============================================================================
  297.       subroutine HigfsZ1(Uz,Uzo,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
  298.       implicit none
  300.       include 'param.f'
  302.       integer in, nx, nz, r, c
  303.       real Uz(MAXCOLS,MAXROWS),Uzo(MAXCOLS,MAXROWS),
  304.      &     bvp(MAXCOLS),bvs(MAXCOLS),
  305.      &     lvp(MAXROWS),lvs(MAXROWS),rvp(MAXROWS),rvs(MAXROWS),hlp
  307.       do r=1 , nz-2  
  308.          hlp = Uz(1,r)
  309.          Uz(1,r) = 2.0 * Uz(2,r) - Uzo(3,r)
  310.      &        - lvp(r) * (Uzo(2,r) - Uz(3,r) - Uz(1,r))
  311.      &        - lvs(r) * (Uzo(1,r) - 2.0 * Uz(2,r) )
  312.          Uzo(1,r) = hlp
  313.          hlp = Uz(nx,r)
  314.          Uz(nx,r) = 2.0 * Uz(nx-1,r) - Uzo(nx-2,r)
  315.      &        - rvp(r) * (Uzo(nx-1,r) - Uz(nx-2,r) - Uz(nx,r))
  316.      &        - rvs(r) * (Uzo(nx,r) - 2.0 * Uz(nx-1,r) )
  317.          Uzo(nx,r) = hlp
  318.       enddo
  319.       in = nz-1
  320.       do c=2, nx-1 
  321.          hlp = Uz(c,in)
  322.          Uz(c,in) = 2.0 * Uz(c,in-1) - Uzo(c,in-2)
  323.      &        - bvp(c) * (Uzo(c,in-1) - Uz(c,in-2) - Uz(c,in))
  324.      &        - bvs(c) * (Uzo(c,in) - 2.0 * Uz(c,in-1) )
  325.          Uzo(c,in) = hlp
  326.       enddo
  327.       return
  328.       end
  329. C============================================================================
  330.       subroutine HigfsZ2(Uz,Uzo,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
  331.       implicit none
  333.       include 'param.f'
  335.       integer in, nx, nz, r, c
  336.       real Uz(MAXCOLS,MAXROWS),Uzo(MAXCOLS,MAXROWS),
  337.      &     bvp(MAXCOLS),bvs(MAXCOLS),
  338.      &     lvp(MAXROWS),lvs(MAXROWS),rvp(MAXROWS),rvs(MAXROWS)
  339.          
  340.       do r= 1 , nz-2   
  341.          Uz(1,r) = Uz(1,r) - lvp(r) * Uz(2,r) - lvs(r) * Uz(3,r)
  342.          Uz(nx,r) = Uz(nx,r) - rvp(r)*Uz(nx-1,r) - rvs(r)*Uz(nx-2,r)
  343.       enddo
  344.       in = nz-1
  345.       do c=2, nx-1 
  346.          Uz(c,in) = Uz(c,in) - bvp(c)*Uz(c,in-1) - bvs(c)*Uz(c,in-2)
  347.       enddo
  348.       return
  349.       end
  350. C============================================================================