SPEED
Functions/Subroutines
GET_NODE_DEPTH_FROM_CMPLX.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine get_node_depth_from_cmplx (loc_n_num, nn_elev, nn_elem,
 Computes elevation from topography (XYZ.out).
 

Function/Subroutine Documentation

◆ get_node_depth_from_cmplx()

subroutine get_node_depth_from_cmplx ( integer*4, dimension(nn_s)  loc_n_num,
integer*4  nn_elev,
integer*4  nn_elem 
)

Computes elevation from topography (XYZ.out).

Note
See MAKE_NOT_HONORING.f90
Author
Ilario Mazzieri
Date
September, 2013
Version
1.0
Parameters
[in]nn_snumber of local nodes
[in]loc_n_numlocal numeration vector
[in]nn_elevnumber nodes in the triangular grid
[in]x_elevelevation values of local nodes
[in]y_elevelevation values of local nodes
[in]z_elevelevation values of local nodes
[in]nn_elemnumber of triangular elements
[in]node1_elemindex triangle vertex
[in]node2_elemindex triangle vertex
[in]node3_elemindex triangle vertex
[in]cs_nnx_loclength cs_loc
[in]cs_loclocal connectivity vector
[in]xx_svertex x- coordinate of local nodes
[in]yy_svertex y- coordinate of local nodes
[in]zz_svertex z- coordinate of local nodes
[in]nmnumber of materials
[in]tmlabels for material vector
[in]sdpolynomial degree vector
[in]tagmatspecific material tag given in CASE option
[in]max_esmax topography spacing
[in]toltolerance given in CASE option
[in]zz_alluvialelevation of the nodes from alluvial soil
[out]zz_elevationelevation of the nodes from complex topography

Definition at line 48 of file GET_NODE_DEPTH_FROM_CMPLX.f90.

49 xx_elev, yy_elev, zz_elev, &
50 node1_elem, node2_elem, node3_elem, &
51 cs_nnz_loc, cs_loc, nm, tm, sd, &
52 nn_s, xx_s, yy_s, zz_s, &
53 zz_elevation, zz_alluvial, &
54 tagmat, max_es,tol)
55
56 implicit none
57
58 integer*4 :: im,ie,i,j,k,nn,ip,isn,ic
59 integer*4 :: nn_elev, nn_elem, cs_nnz_loc, nm, ne, nn_s
60 integer*4 :: h
61 integer*4 :: tagmat
62
63 integer*4, dimension(nn_s) :: loc_n_num
64 integer*4, dimension(nn_elem) :: node1_elem,node2_elem,node3_elem
65 integer*4, dimension(nm) :: tm
66 integer*4, dimension(nm) :: sd
67 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
68
69 real*8 :: x1,y1,z1
70 real*8 :: x2,y2,z2
71 real*8 :: x3,y3,z3
72 real*8 :: ux,uy,uz,vx,vy,vz
73 real*8 :: a,b,c
74 real*8 :: max_es
75 real*8 :: zz_interp
76 real*8 :: v0x,v0y,v1x,v1y,v2x,v2y
77 real*8 :: dot00,dot01,dot02,dot11,dot12
78 real*8 :: invdenom,u,v
79 real*8 :: d2min
80 real*8 :: zz_elev_min
81 real*8 :: dx,dy,dz,tol
82
83 real*8, dimension(:), allocatable :: ct,ww
84 real*8, dimension(nn_elev) :: xx_elev,yy_elev,zz_elev
85 real*8, dimension(nn_s) :: xx_s,yy_s,zz_s
86 real*8, dimension(nn_s) :: zz_elevation
87 real*8, dimension(nn_s) :: zz_alluvial
88
89 real*8, dimension(:,:), allocatable :: dd
90
91 d2min = (5 * max_es)**2
92
93 zz_elev_min = zz_elev(1)
94 do i = 1,nn_elev
95 if (zz_elev(i).lt.zz_elev_min) then
96 zz_elev_min = zz_elev(i)
97 endif
98 enddo
99
100
101 nn = 2
102 allocate(ct(nn),ww(nn),dd(nn,nn))
103 call make_lgl_nw(nn,ct,ww,dd)
104
105 ne = cs_loc(0) - 1
106
107 do im = 1,nm
108 if ((sd(im) +1).ne.nn) then
109 deallocate(ct,ww,dd)
110 nn = sd(im) +1
111 allocate(ct(nn),ww(nn),dd(nn,nn))
112 call make_lgl_nw(nn,ct,ww,dd)
113 endif
114
115 do ie = 1,ne
116 if (cs_loc(cs_loc(ie -1) +0).eq.tagmat) then
117
118 do k = 1,nn
119 do j = 1,nn
120 do i = 1,nn
121
122 ip = nn*nn*(k -1) +nn*(j -1) +i
123 isn = cs_loc(cs_loc(ie -1) + ip)
124 ic = isn
125
126
127
128! if ((zz_elevation(ic).eq.-1.0e+30).and.(zz_alluvial(ic).ge.0.0d0)) then
129 if (zz_elevation(ic).eq.-1.0e+30) then
130 do h = 1,nn_elem
131
132 x1 = xx_elev(node1_elem(h))
133 y1 = yy_elev(node1_elem(h))
134 z1 = zz_elev(node1_elem(h))
135
136 if (((x1 - xx_s(ic))**2 + (y1 - yy_s(ic))**2).le.d2min) then
137
138 x2 = xx_elev(node2_elem(h))
139 y2 = yy_elev(node2_elem(h))
140 z2 = zz_elev(node2_elem(h))
141
142 x3 = xx_elev(node3_elem(h))
143 y3 = yy_elev(node3_elem(h))
144 z3 = zz_elev(node3_elem(h))
145
146 !Point in triangle test
147
148 ! P = (xx_s(ic) yy_s(ic))
149 ! A = (X1 Y1)
150 ! B = (X2 Y2)
151 ! C = (X3 Y3)
152
153 ! Compute vectors
154 ! v0 = C - A
155 v0x=(x3 - x1)
156 v0y=(y3 - y1)
157
158 !v0x=(X3 - X1)/sqrt((X3 - X1)**2+(Y3 - Y1)**2)
159 !v0y=(Y3 - Y1)/sqrt((X3 - X1)**2+(Y3 - Y1)**2)
160
161 ! v1 = B - A
162 v1x=(x2 - x1)
163 v1y=(y2 - y1)
164
165 !v1x=(X2 - X1)/sqrt((X2 - X1)**2+(Y2 - Y1)**2)
166 !v1y=(Y2 - Y1)/sqrt((X2 - X1)**2+(Y2 - Y1)**2)
167
168 ! v2 = P - A
169 v2x=(xx_s(ic) - x1)
170 v2y=(yy_s(ic) - y1)
171
172 !v2x=(xx_s(ic) - X1)/sqrt((xx_s(ic) - X1)**2+(yy_s(ic) - Y1)**2)
173 !v2y=(yy_s(ic) - Y1)/sqrt((xx_s(ic) - X1)**2+(yy_s(ic) - Y1)**2)
174
175 ! Compute dot products
176 ! [u].[v] = ux * vx + uy * vy
177 ! dot([u],[v])
178
179 !dot00 = dot(v0, v0)
180 dot00 = v0x * v0x + v0y * v0y
181
182 !dot01 = dot(v0, v1)
183 dot01 = v0x * v1x + v0y * v1y
184
185 !dot02 = dot(v0, v2)
186 dot02 = v0x * v2x + v0y * v2y
187
188 !dot11 = dot(v1, v1)
189 dot11 = v1x * v1x + v1y * v1y
190
191 !dot12 = dot(v1, v2)
192 dot12 = v1x * v2x + v1y * v2y
193
194 ! Compute barycentric coordinates
195 invdenom = 1 / (dot00 * dot11 - dot01 * dot01)
196 u = (dot11 * dot02 - dot01 * dot12) * invdenom
197 v = (dot00 * dot12 - dot01 * dot02) * invdenom
198
199 !Point in triangle test
200 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
201
202 if ( (u.ge.0.0d0).and.(v.ge.0.0d0).and.((u + v).le.1.0d0) ) then
203
204 ! Build up the plane passing through the points P1, P2 and P3
205
206 ux=(x1-x2)/sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
207 uy=(y1-y2)/sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
208 uz=(z1-z2)/sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
209 vx=(x3-x2)/sqrt((x3-x2)**2+(y3-y2)**2+(z3-z2)**2)
210 vy=(y3-y2)/sqrt((x3-x2)**2+(y3-y2)**2+(z3-z2)**2)
211 vz=(z3-z2)/sqrt((x3-x2)**2+(y3-y2)**2+(z3-z2)**2)
212
213 a = uy * vz - uz * vy
214 b = uz * vx - ux * vz
215 c = ux * vy - uy * vx
216
217 zz_interp = -a/c * (xx_s(ic)-x1) -b/c * (yy_s(ic)-y1) + z1
218 zz_elevation(ic) = ( zz_interp - zz_s(ic) )
219
220 if (abs(zz_elevation(ic)).lt.tol) then
221 zz_elevation(ic) = 0.0d0
222 endif
223 endif !if ( (u.ge.0.0d0).and.(v.ge.0.0d0).and.((u + v).le.1.0d0) ) then
224
225 if ( (u.ge.0.0d0).and.(v.ge.0.0d0).and.((u + v).le.1.0d0) ) exit
226
227 endif !if (((X1 - xx_s(ic))**2 + (Y1 - yy_s(ic))**2).le.d2min) then
228
229 enddo !do h = 1,nn_elem
230 endif !if ((zz_elevation(ic).eq.-1.0e+30).and.(zz_alluvial(ic).ge.0.0d0)) then
231 enddo !do k = 1,nn
232 enddo !do j = 1,nn
233 enddo !do i = 1,nn
234
235
236 endif
237 enddo
238 enddo
239
240 return
241
make_lgl_nw
subroutine make_lgl_nw(nb_pnt, xq, wq, dd)
Makes Gauss-Legendre-Lobatto nodes, weigths and spectral derivatives.
Definition MAKE_LGL_NW.f90:29

References make_lgl_nw().

Referenced by make_nothonoring().

Here is the call graph for this function:
Here is the caller graph for this function: