READ_SISM.f90

Go to the documentation of this file.

!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.
 
 
      subroutine read_sism(Xipo,Yipo,Zipo, &                
                                 X1,Y1,Z1, &                        
                                 X2,Y2,Z2, &                        
                                 X3,Y3,Z3, &                        
                                 nnod,xs,ys,zs, &                
                                 num_ns,sour_ns,i_sism,&        
                                 dist_sour_ns,nl_sism,&                
                                 max_num_ns,loc_n_num, nn_loc, &
                                 pos_sour_nx,pos_sour_ny,pos_sour_nz)                
 
 
      implicit none
      
      integer*4 :: i_sism,isn,i,nn_loc
      integer*4 :: nnod,node_sism,num_ns,nl_sism,max_num_ns
      
      integer*4, dimension(nn_loc) :: loc_n_num(nn_loc)
 
      integer*4, dimension(max_num_ns,nl_sism) :: sour_ns
 
      real*8 :: xipo,yipo,zipo,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,ux,uy,uz,vx,vy,vz,tol
      real*8 :: p1x,p1y,p1z,p2x,p2y,p2z 
      real*8 :: epsilon, twopi, anglesum, costheta,m1,m2
 
      real*8, dimension(nn_loc) :: xs,ys,zs
      real*8, dimension(4) :: x,y,z
 
      real*8, dimension(max_num_ns,nl_sism) :: dist_sour_ns
      real*8, dimension(max_num_ns,nl_sism) :: pos_sour_nx,pos_sour_ny,pos_sour_nz
 
      node_sism = 0
 
      tol = 6.28
      epsilon = 1.d-8
      twopi = 6.283185307179586476925287
                    
      x(1) = x1
      y(1) = y1
      z(1) = z1
 
      x(2) = x2
      y(2) = y2
      z(2) = z2 
      
      x(3) = x3
      y(3) = y3
      z(3) = z3
 
      x(4) = x1
      y(4) = y1
      z(4) = z1
 
 
      do isn = 1,nnod
 
          anglesum = 0.0d0
          do i = 1,3
                   p1x = x(i) - xs(isn)
                  p1y = y(i) - ys(isn)
                  p1z = z(i) - zs(isn)
                  p2x = x(i+1) - xs(isn)
                  p2y = y(i+1) - ys(isn)
                  p2z = z(i+1) - zs(isn)
 
                  m1 = dsqrt(p1x*p1x + p1y*p1y + p1z*p1z)
                  m2 = dsqrt(p2x*p2x + p2y*p2y + p2z*p2z)
      
                  if ((m1*m2).le.epsilon) then
                        anglesum = twopi
                  else
                        costheta = (p1x*p2x + p1y*p2y + p1z*p2z)/(m1*m2)
                        anglesum = anglesum + dacos(costheta)
                  endif
          enddo
 
          if (anglesum.ge.tol) then
 
              node_sism = node_sism + 1
              sour_ns(node_sism,i_sism) = loc_n_num(isn) 
              dist_sour_ns(node_sism,i_sism) = dsqrt((xipo - xs(isn))**2 +(yipo - ys(isn))**2 +(zipo - zs(isn))**2)
              pos_sour_nx(node_sism,i_sism) = xs(isn)
              pos_sour_ny(node_sism,i_sism) = ys(isn)              
              pos_sour_nz(node_sism,i_sism) = zs(isn)
 
          endif
 
      enddo
 
     
      return
      end subroutine read_sism