SPEED
GET_NODE_DEPTH_FROM_ALLUVIAL.f90
Go to the documentation of this file.
1! Copyright (C) 2012 The SPEED FOUNDATION
2! Author: Ilario Mazzieri
3!
4! This file is part of SPEED.
5!
6! SPEED is free software; you can redistribute it and/or modify it
7! under the terms of the GNU Affero General Public License as
8! published by the Free Software Foundation, either version 3 of the
9! License, or (at your option) any later version.
10!
11! SPEED is distributed in the hope that it will be useful, but
12! WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14! Affero General Public License for more details.
15!
16! You should have received a copy of the GNU Affero General Public License
17! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
18
46
47
48 subroutine get_node_depth_from_alluvial(loc_n_num, n_elev, nn_elem, &
49 x_elev, y_elev, z_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, zz_elevation, &
53 tagmat, max_as, tol)
54
55 implicit none
56
57 integer*4 :: n_elev,nn_elem,cs_nnz_loc,nm,ne,nn_s
58 integer*4 :: im,ie,i,j,k,nn,ip,isn,ic
59 integer*4 :: h
60 integer*4 :: tagmat
61
62 integer*4, dimension(nn_s) :: loc_n_num
63 integer*4, dimension(nn_elem) :: node1_elem,node2_elem,node3_elem
64 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
65 integer*4, dimension(nm) :: tm
66 integer*4, dimension(nm) :: sd
67
68 real*8 :: dx,dy,dz,tol
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_as
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 :: z_elev_min
81
82 real*8, dimension(:), allocatable :: ct,ww
83 real*8, dimension(n_elev) :: x_elev,y_elev,z_elev
84 real*8, dimension(nn_s) :: xx_s,yy_s,zz_s
85 real*8, dimension(nn_s) :: zz_elevation
86
87 real*8, dimension(:,:), allocatable :: dd
88
89
90 d2min = (5 * max_as)**2
91
92 z_elev_min = z_elev(1)
93 do i = 1,n_elev
94 if (z_elev(i).lt.z_elev_min) then
95 z_elev_min = z_elev(i)
96 endif
97 enddo
98
99 ne = cs_loc(0) - 1
100
101
102! nn = 2
103! allocate(ct(nn),ww(nn),dd(nn,nn))
104! call MAKE_LGL_NW(nn,ct,ww,dd)
105
106! do im = 1,nm
107! if ((sd(im) +1).ne.nn) then
108! deallocate(ct,ww,dd)
109! nn = sd(im) +1
110! allocate(ct(nn),ww(nn),dd(nn,nn))
111! call MAKE_LGL_NW(nn,ct,ww,dd)
112! endif
113
114 do ie = 1,ne
115 im = cs_loc(cs_loc(ie -1) +0);
116 nn = sd(im) +1
117
118 if (im .eq. tagmat) then
119
120
121 allocate(ct(nn),ww(nn),dd(nn,nn))
122 call make_lgl_nw(nn,ct,ww,dd)
123
124 do k = 1,nn
125 do j = 1,nn
126 do i = 1,nn
127
128 ip = nn*nn*(k -1) +nn*(j -1) +i
129 isn = cs_loc(cs_loc(ie -1) + ip)
130 ic = isn
131
132
133
134 if (zz_elevation(ic) .eq. -1.0e+30) then
135
136 do h = 1,nn_elem
137 if (zz_s(ic).lt.z_elev_min) exit
138
139 x1 = x_elev(node1_elem(h))
140 y1 = y_elev(node1_elem(h))
141 z1 = z_elev(node1_elem(h))
142
143 if (((x1 - xx_s(ic))**2 + (y1 - yy_s(ic))**2).le.d2min) then
144
145 x2 = x_elev(node2_elem(h))
146 y2 = y_elev(node2_elem(h))
147 z2 = z_elev(node2_elem(h))
148
149 x3 = x_elev(node3_elem(h))
150 y3 = y_elev(node3_elem(h))
151 z3 = z_elev(node3_elem(h))
152
153 if (zz_s(ic).ge.min(z1,z2,z3)) then
154
155 !+ Point in triangle test
156
157 ! P = (xx_s(ic) yy_s(ic))
158 ! A = (X1 Y1)
159 ! B = (X2 Y2)
160 ! C = (X3 Y3)
161 ! Compute vectors
162 ! v0 = C - A
163 v0x=(x3 - x1)
164 v0y=(y3 - y1)
165
166 ! v1 = B - A
167 v1x=(x2 - x1)
168 v1y=(y2 - y1)
169
170 ! v2 = P - A
171 v2x=(xx_s(ic) - x1)
172 v2y=(yy_s(ic) - y1)
173
174 ! Compute dot products
175 ! [u].[v] = ux * vx + uy * vy
176 ! dot([u],[v])
177
178 !dot00 = dot(v0, v0)
179 dot00 = v0x * v0x + v0y * v0y
180
181 !dot01 = dot(v0, v1)
182 dot01 = v0x * v1x + v0y * v1y
183
184 !dot02 = dot(v0, v2)
185 dot02 = v0x * v2x + v0y * v2y
186
187 !dot11 = dot(v1, v1)
188 dot11 = v1x * v1x + v1y * v1y
189
190 !dot12 = dot(v1, v2)
191 dot12 = v1x * v2x + v1y * v2y
192
193 ! Compute barycentric coordinates
194 invdenom = 1 / (dot00 * dot11 - dot01 * dot01)
195 u = (dot11 * dot02 - dot01 * dot12) * invdenom
196 v = (dot00 * dot12 - dot01 * dot02) * invdenom
197 !Point in triangle test
198 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
199
200 if ( (u.ge.0.0d0).and.(v.ge.0.0d0).and.((u + v).le.1.0d0) ) then
201
202 ! Build up the plane passing through the points P1, P2 and P3
203
204 ux=(x1-x2)/sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
205 uy=(y1-y2)/sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
206 uz=(z1-z2)/sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)
207 vx=(x3-x2)/sqrt((x3-x2)**2+(y3-y2)**2+(z3-z2)**2)
208 vy=(y3-y2)/sqrt((x3-x2)**2+(y3-y2)**2+(z3-z2)**2)
209 vz=(z3-z2)/sqrt((x3-x2)**2+(y3-y2)**2+(z3-z2)**2)
210
211 a = uy * vz - uz * vy
212 b = uz * vx - ux * vz
213 c = ux * vy - uy * vx
214
215 zz_interp = -a/c * (xx_s(ic)-x1) -b/c * (yy_s(ic)-y1) + z1
216 zz_elevation(ic) = -1 * ( zz_interp - zz_s(ic) )
217 if (abs(zz_elevation(ic)).lt.tol) then
218 zz_elevation(ic) = 0.0d0
219 endif
220
221 endif !if ( (u.ge.0.0d0).and.(v.ge.0.0d0).and.((u + v).le.1.0d0) ) then
222
223 if ( (u.ge.0.0d0).and.(v.ge.0.0d0).and.((u + v).le.1.0d0) ) exit
224
225 endif !if (zz_s(ic).ge.min(Z1,Z2,Z3)) then
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
231 endif !if (zz_elevation(ic).le.0.0d0) then
232
233
234
235 enddo !do k = 1,nn
236 enddo !do j = 1,nn
237 enddo !do i = 1,nn
238
239 deallocate(ct,ww,dd)
240
241 endif
242 enddo
243! enddo
244
245 return
246
247 end subroutine get_node_depth_from_alluvial
248
get_node_depth_from_alluvial
subroutine get_node_depth_from_alluvial(loc_n_num, n_elev, nn_elem
Computes elevation from alluvial basin (ALL.out).
Definition GET_NODE_DEPTH_FROM_ALLUVIAL.f90:49
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