GET_NODE_DEPTH_AND_VS.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 get_node_depth_and_vs(loc_n_num, nn_elev, nn_elem, &
                                            xx_elev, yy_elev, zz_elev, vs_elev, sediments, &
                                            node1_elem, node2_elem, node3_elem, &
                                            cs_nnz_loc, cs_loc, nm, tm, sd, &
                                            nn_s, xx_s, yy_s, zz_s, &
                                            zz_elevation, zz_alluvial, vs, thickness, &
                                            tagmat, max_es,tol)      
 
      implicit none
      
      integer*4 :: im,ie,i,j,k,nn,ip,isn,ic
      integer*4 :: nn_elev, nn_elem, cs_nnz_loc, nm, ne, nn_s
      integer*4 :: h, h_sel
      integer*4 :: tagmat
 
      integer*4, dimension(nn_s) :: loc_n_num                                                         
      integer*4, dimension(nn_elem) :: node1_elem,node2_elem,node3_elem
      integer*4, dimension(nm) :: tm
      integer*4, dimension(nm) :: sd
      integer*4, dimension(0:cs_nnz_loc) :: cs_loc
 
      real*8 :: x1,y1,z1                                
      real*8 :: x2,y2,z2                                
      real*8 :: x3,y3,z3                                
      real*8 :: ux,uy,uz,vx,vy,vz                
      real*8 :: a,b,c                                        
      real*8 :: max_es                                        
      real*8 :: zz_interp                                
      real*8 :: v0x,v0y,v1x,v1y,v2x,v2y                        
      real*8 :: dot00,dot01,dot02,dot11,dot12        
      real*8 :: invdenom,u,v                                        
      real*8 :: d2min 
      real*8 :: zz_elev_min 
      real*8 :: dx,dy,dz,tol, dist, dist_min
 
      real*8, dimension(:), allocatable :: ct,ww
      real*8, dimension(nn_elev) :: xx_elev,yy_elev,zz_elev
      real*8, dimension(nn_elem) :: vs_elev,sediments 
      real*8, dimension(nn_s) :: xx_s,yy_s,zz_s
      real*8, dimension(nn_s) :: zz_elevation
      real*8, dimension(nn_s) :: zz_alluvial
      real*8, dimension(nn_s) :: vs, thickness
      
      real*8, dimension(:,:), allocatable :: dd
 
      d2min = (5 * max_es)**2 
          
      zz_elev_min = zz_elev(1)
      do i = 1,nn_elev
         if (zz_elev(i).lt.zz_elev_min) then
                 zz_elev_min = zz_elev(i)
          endif
      enddo
      !write(*,*) d2min, max_es
 
      
      ne = cs_loc(0) - 1
      
      do ie = 1,ne
         
         im = cs_loc(cs_loc(ie -1) +0)
         
         if (im .eq. tagmat) then
         
             nn = sd(im) +1
             allocate(ct(nn),ww(nn),dd(nn,nn))
             call make_lgl_nw(nn,ct,ww,dd)  
 
                            
               do k = 1,nn
                  do j = 1,nn
                     do i = 1,nn
                        
                        ip = nn*nn*(k -1) +nn*(j -1) +i
                        isn = cs_loc(cs_loc(ie -1) + ip)
                        ic = isn
                        
                        dist_min = 1.d+30;
                        h_sel = 0;
                        do h = 1,nn_elem
                            
                            
                            x1 = xx_elev(node1_elem(h)) 
                            y1 = yy_elev(node1_elem(h)) 
                            z1 = zz_elev(node1_elem(h)) 
                            
                            dist = (x1 - xx_s(ic))**2 + (y1 - yy_s(ic))**2
                            if (dist .le. dist_min) then
                               h_sel = h;
                               dist_min = dist;
                            endif          
                        enddo
                                     
                            !write(*,*) h_sel 
                            !read(*,*)                            
                            x1 = xx_elev(node1_elem(h_sel)) 
                            y1 = yy_elev(node1_elem(h_sel)) 
                            z1 = zz_elev(node1_elem(h_sel)) 
 
                            x2 = xx_elev(node2_elem(h_sel)) 
                            y2 = yy_elev(node2_elem(h_sel)) 
                            z2 = zz_elev(node2_elem(h_sel)) 
 
                            x3 = xx_elev(node3_elem(h_sel)) 
                            y3 = yy_elev(node3_elem(h_sel)) 
                            z3 = zz_elev(node3_elem(h_sel)) 
 
                           
                            zz_interp = (z1 + z2 + z3)/3.d0
                            zz_elevation(ic) = ( zz_interp - zz_s(ic) )
                            vs(ic) = vs_elev(h_sel)
                            thickness(ic) = sediments(h_sel)
                                                        
                            if(zz_elevation(ic) .le. tol) then
                               zz_elevation(ic) = 0.0d0
                               vs(ic) = vs_elev(h_sel)
                               thickness(ic) = sediments(h_sel)
                            endif
                                                        
                                                                                
 
 
              enddo !do k = 1,nn
            enddo !do j = 1,nn
          enddo !do i = 1,nn
 
          deallocate(ct,ww,dd)
                    
        endif !if im.eq.tagmat
      enddo
 
      
      return
      
      subroutine get_node_depth_and_vs(loc_n_num, nn_elev, nn_elem, & …
      end subroutine get_node_depth_and_vs