CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
liquidchannel.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine liquidchannel (node1, node2, nodem, nelem, lakon, nactdog, identity, ielprop, prop, iflag, v, xflow, f, nodef, idirf, df, rho, g, co, dvi, numf, mi, ipkon, kon)
 

Function/Subroutine Documentation

◆ liquidchannel()

subroutine liquidchannel ( integer  node1,
integer  node2,
integer  nodem,
integer  nelem,
character*8, dimension(*)  lakon,
integer, dimension(0:3,*)  nactdog,
logical  identity,
integer, dimension(*)  ielprop,
real*8, dimension(*)  prop,
integer  iflag,
real*8, dimension(0:mi(2),*)  v,
real*8  xflow,
real*8  f,
integer, dimension(*)  nodef,
integer, dimension(*)  idirf,
real*8, dimension(*)  df,
real*8  rho,
real*8, dimension(3)  g,
real*8, dimension(3,*)  co,
real*8  dvi,
integer  numf,
integer, dimension(*)  mi,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon 
)
22 !
23 ! open channel for incompressible media
24 !
25 ! SG: sluice gate
26 ! SO: sluice opening
27 ! WE: weir
28 ! WO: weir opening
29 ! DS: discontinuous slope
30 ! DO: discontinuous slope opening
31 ! : default channel mit linearly varying trapezoid cross
32 ! section
33 !
34  implicit none
35 !
36  logical identity,bresse,jump
37  character*8 lakon(*)
38 !
39  integer nelem,nactdog(0:3,*),node1,node2,nodem,indexup,i,
40  & ielprop(*),nodef(*),idirf(*),index,iflag,mi(*),nsol,
41  & inv,numf,nodesg,nelemdown,nelemup,node0,kon(*),ipkon(*)
42 !
43  real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),b,d,c,p,
44  & h1,h2,rho,dvi,friction,reynolds,dg,b1,b2,
45  & g(3),dl,xks,z1,z2,co(3,*),xflow2,dyg3dbj,dyg4dbj,
46  & s0,sqrts0,hk,form_fact,h1ns,h2ns,h0,dyg3deta,dyg4deta,
47  & dh3dh1,dh4dh2,dh3dm,dh4dm,eta,da3deta,da4deta,bj,
48  & theta,cth,tth,um1,um2,a1,a2,p1,p2,d1,d2,da1dh1,da2dh2,
49  & dp1dh1,dp2dh2,dd1dh1,dd2dh2,h3,h4,dh3deta,xn1,xn2,xt1,xt2,
50  & dh4deta,yg3,yg4,dyg3dh3,dyg4dh4,a3,a4,da3dh3,da4dh4,
51  & dum1dh1,dum2dh2,c1,c2,dbds,dbjdeta,e0,e1,e2,e3,
52  & dyg3dm,dyg4dm,da3dm,da4dm,dyg3dh1,dyg4dh2,
53  & da3dh1,da4dh2,solreal(3),solimag(3),dist
54 !
55 ! iflag=0: check whether all parameters in the element equation
56 ! are known => equation is not needed
57 ! iflag=1: calculation of the initial flux
58 ! iflag=2: evaluate the element equation and all derivatives
59 ! iflag=3: correct the channel depth in order to move a jump
60 !
61  if (iflag.eq.0) then
62  identity=.true.
63 !
64  if(nactdog(2,node1).ne.0)then
65  identity=.false.
66  elseif(nactdog(2,node2).ne.0)then
67  identity=.false.
68  elseif(nactdog(1,nodem).ne.0)then
69  identity=.false.
70  endif
71 !
72  elseif((iflag.eq.1).or.(iflag.eq.2))then
73 !
74  index=ielprop(nelem)
75 !
76  h1=v(2,node1)
77  h2=v(2,node2)
78 !
79  z1=-g(1)*co(1,node1)-g(2)*co(2,node1)-g(3)*co(3,node1)
80  z2=-g(1)*co(1,node2)-g(2)*co(2,node2)-g(3)*co(3,node2)
81 !
82  dg=dsqrt(g(1)*g(1)+g(2)*g(2)+g(3)*g(3))
83 !
84  if(iflag.eq.1) then
85 !
86 ! in a first call of liquidchannel the flow is determined,
87 ! in a second call the channel depth is calculated
88 !
89  if(lakon(nelem)(6:7).eq.'SG') then
90 !
91 ! sluice gate
92 !
93  b=prop(index+1)
94  s0=prop(index+2)
95  if(s0.lt.-1.d0) then
96  s0=dasin((z1-z2)/dl)
97  endif
98  sqrts0=dsqrt(1.d0-s0*s0)
99  theta=0.d0
100  h2=prop(index+3)
101 !
102  if(dabs(xflow).lt.1.d-30) then
103 !
104 ! determine initial mass flow
105 !
106  if(nactdog(2,node1).ne.0) then
107 !
108 ! upstream level not known
109 !
110  xflow=0.d0
111  else
112  xflow=2.d0*dg*(rho*b*h2)**2*(h1-h2*sqrts0)
113  if(xflow.lt.0.d0) then
114  write(*,*)'*ERROR in liquidchannel: water level'
115  write(*,*) ' upstream of sluice gate is '
116  write(*,*) ' smaller than downstream heigh
117  &t'
118  call exit(201)
119  else
120  xflow=dsqrt(xflow)
121  endif
122  endif
123  else
124 !
125 ! determine the downstream depth
126 ! and the upstream depth if not defined as BC
127 !
128  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
129  v(3,node2)=hk
130  if(h2.gt.hk) then
131 !
132 ! for initial conditions
133 !
134  if(nactdog(2,node1).ne.0) v(2,node1)=3.d0*hk/2.d0
135  v(2,node2)=hk
136  else
137 !
138 ! for initial conditions
139 !
140  if(nactdog(2,node1).ne.0) v(2,node1)=
141  & xflow**2/(2.d0*dg*(rho*b*h2)**2)+h2*sqrts0
142  v(2,node2)=h2
143  endif
144  endif
145  elseif(lakon(nelem)(6:7).eq.'WE') then
146 !
147 ! weir
148 !
149  b=prop(index+1)
150  p=prop(index+2)
151  c=prop(index+3)
152  sqrts0=1.d0
153  theta=0.d0
154 !
155  if(dabs(xflow).lt.1.d-30) then
156 !
157 ! determine initial mass flow
158 !
159  if(nactdog(2,node1).ne.0) then
160 !
161 ! upstream level unknown
162 !
163  xflow=0.d0
164  else
165  if(h1.le.p) then
166  write(*,*) '*ERROR in liquidchannel'
167  write(*,*) ' weir height exceeds'
168  write(*,*) ' upstream level'
169  call exit(201)
170  endif
171  xflow=rho*c*b*(h1-p)**(1.5d0)
172  endif
173  else
174 !
175 ! determine the downstream depth
176 ! and the upstream depth if not defined as BC
177 !
178  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
179  v(3,node2)=hk
180 !
181 ! for initial conditions
182 !
183  if(nactdog(2,node1).ne.0) v(2,node1)=p+3.d0*hk/2.d0
184 !
185 ! next value is used for downstream initial values
186 !
187  v(2,node2)=hk
188  endif
189 !
190  elseif(lakon(nelem)(6:7).eq.'DS') then
191  if(dabs(xflow).lt.1.d-30) then
192 !
193 ! initial mass flow cannot be determined for this
194 ! type of element
195 !
196  xflow=0.d0
197  else
198 !
199 ! determine the downstream depth
200 !
201  b=prop(index+1)
202  s0=prop(index+2)
203  if(s0.lt.-1.d0) then
204  s0=dasin((z1-z2)/dl)
205  endif
206  sqrts0=dsqrt(1.d0-s0*s0)
207  theta=prop(index+4)
208 !
209  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
210  v(3,node2)=hk
211 !
212 ! initial condition for fluid depth
213 ! supercritical value
214 !
215  v(2,node2)=hk/2.d0
216  endif
217 !
218  endif
219  else
220 !
221 ! calculating f and its derivatives
222 !
223  bresse=.false.
224  jump=.false.
225 !
226  xflow2=xflow*xflow
227 !
228 ! element properties
229 !
230  if((lakon(nelem)(6:7).eq.'SG').or.
231  & (lakon(nelem)(6:7).eq.'SO').or.
232  & (lakon(nelem)(6:7).eq.'WO').or.
233  & (lakon(nelem)(6:7).eq.'RE').or.
234  & (lakon(nelem)(6:7).eq.' ').or.
235  & (lakon(nelem)(6:7).eq.'DS').or.
236  & (lakon(nelem)(6:7).eq.'DO')) then
237  b=prop(index+1)
238  s0=prop(index+2)
239  if(s0.lt.-1.d0) then
240  s0=dasin((z1-z2)/dl)
241  endif
242  sqrts0=dsqrt(1.d0-s0*s0)
243  if(lakon(nelem)(6:7).ne.'SG') then
244  dl=prop(index+3)
245  theta=prop(index+4)
246  xks=prop(index+5)
247  if(dl.le.0.d0) then
248  dl=dsqrt((co(1,node2)-co(1,node1))**2+
249  & (co(2,node2)-co(2,node1))**2+
250  & (co(3,node2)-co(3,node1))**2)
251  endif
252  else
253  theta=0.d0
254  endif
255  elseif(lakon(nelem)(6:7).eq.'WE') then
256  b=prop(index+1)
257  p=prop(index+2)
258  c=prop(index+3)
259  sqrts0=1.d0
260  theta=0.d0
261  elseif((lakon(nelem)(6:7).eq.'CO').or.
262  & (lakon(nelem)(6:7).eq.'EL')) then
263  b1=prop(index+1)
264 !
265  s0=prop(index+2)
266  if(s0.lt.-1.d0) then
267  s0=0.d0
268  endif
269  sqrts0=dsqrt(1.d0-s0*s0)
270 !
271  dl=prop(index+3)
272  if(dl.le.0.d0) then
273  dl=dsqrt((co(1,node2)-co(1,node1))**2+
274  & (co(2,node2)-co(2,node1))**2+
275  & (co(3,node2)-co(3,node1))**2)
276  endif
277 !
278  b2=prop(index+4)
279  b=(b1+b2)/2.d0
280  theta=0.d0
281  xks=0.d0
282  elseif((lakon(nelem)(6:7).eq.'ST').or.
283  & (lakon(nelem)(6:7).eq.'DR')) then
284  b=prop(index+1)
285 !
286  s0=prop(index+2)
287  if(s0.lt.-1.d0) then
288  s0=0.d0
289  endif
290  sqrts0=dsqrt(1.d0-s0*s0)
291 !
292  dl=prop(index+3)
293  if(dl.le.0.d0) then
294  dl=dsqrt((co(1,node2)-co(1,node1))**2+
295  & (co(2,node2)-co(2,node1))**2+
296  & (co(3,node2)-co(3,node1))**2)
297  endif
298 !
299  d=prop(index+4)
300  b1=b
301  b2=b
302  theta=0.d0
303  xks=0.d0
304  endif
305 !
306  if(xflow.ge.0.d0) then
307  inv=1
308  else
309  inv=-1
310  endif
311 !
312 ! standard element equation: unknowns are the mass flow
313 ! and the depth upstream and downstream
314 !
315  numf=3
316  nodef(1)=node1
317  nodef(2)=nodem
318  nodef(3)=node2
319  idirf(1)=2
320  idirf(2)=1
321  idirf(3)=2
322 !
323  if(lakon(nelem)(6:7).eq.'SG') then
324 !
325 ! sluice gate
326 ! 1-SG-2-SO-3
327 !
328 ! h2 cannot exceed HKmax
329 !
330  h2=prop(index+3)
331  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
332  v(3,node2)=hk
333  if(h2.gt.hk) h2=hk
334 !
335  nelemdown=nint(prop(index+5))
336  h3=v(2,kon(ipkon(nelemdown)+3))
337  call hns(b,theta,rho,dg,sqrts0,xflow,h2,h2ns)
338  if(h3.lt.h2ns) then
339 !
340 ! Q=f_SG(h1,h2): sluice gate equation between
341 ! 1 and 2
342 !
343 ! next line for output only
344 !
345  v(2,node2)=h2
346 c write(30,*) 'SG: sluice gate equation '
347 c write(30,*)'h1= ',h1,'h2= ',h2,'h3= ',h3,'h2ns= ',h2ns
348  df(1)=2.d0*dg*(rho*b*h2)**2
349  df(2)=-2.d0*xflow
350  f=df(1)*(h1-h2*sqrts0)
351  df(3)=2.d0*f/h2-df(1)*sqrts0
352  f=f-xflow2
353  else
354 !
355 ! fake equation
356 !
357 c write(30,*) 'SG: fake equation '
358 c write(30,*)'h1= ',h1,'h2= ',h2,'h3= ',h3,'h2ns= ',h2ns
359  numf=1
360  nodef(1)=nodem
361  idirf(1)=3
362  f=prop(index+4)-0.5d0
363  df(1)=1.d0
364  endif
365  elseif(lakon(nelem)(6:7).eq.'SO') then
366 !
367 ! sluice opening (element streamdown of sluice gate)
368 ! 0-SG-1-SO-2
369 !
370  nelemup=nint(prop(index+6))
371  node0=kon(ipkon(nelemup)+1)
372  h0=v(2,node0)
373  h1=prop(ielprop(nelemup)+3)
374 !
375 ! h1 cannot exceed HKmax
376 !
377  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
378  v(3,node2)=hk
379  if(h1.gt.hk) h1=hk
380 !
381  call hns(b,theta,rho,dg,sqrts0,xflow,h1,h1ns)
382  if(h2.lt.h1ns) then
383 !
384 ! bresse (frontwater)
385 !
386 c write(30,*) 'SO: Bresse equation '
387 c write(30,*)'h0= ',h0,'h1= ',h1,'h2= ',h2,'h1ns= ',h1ns
388  bresse=.true.
389  else
390 !
391 ! Q=f_SG(h0,h2): sluice gate equation between 0 and 2
392 ! (backwater)
393 !
394 ! reset gate height
395 !
396  h1=prop(ielprop(nelemup)+3)
397 !
398 c write(30,*) 'SO: Sluice gate eqn. between 0 and 2 '
399 c write(30,*)'h0= ',h0,'h1= ',h1,'h2= ',h2,'h1ns= ',h1ns
400  numf=4
401  nodef(4)=node0
402  idirf(4)=2
403 !
404  if(h2.gt.h1) then
405 !
406 ! gate flow (water touches gate)
407 ! section = b * h1
408 !
409 ! next line for output only
410 !
411  v(2,node1)=h1
412  df(4)=2.d0*dg*(rho*b*h1)**2
413  df(3)=-df(4)*sqrts0
414  df(2)=-2.d0*xflow
415  f=df(4)*(h0-h2*sqrts0)
416  df(1)=2.d0*f/h1
417  else
418 !
419 ! incomplete inflexion (water does not touch gate)
420 ! section = b * h2
421 !
422 ! next line for output only
423 !
424  v(2,node1)=h2
425  df(4)=2.d0*dg*(rho*b*h2)**2
426  df(3)=-df(4)*sqrts0
427  df(2)=-2.d0*xflow
428  f=df(4)*(h0-h2*sqrts0)
429  df(3)=df(3)+2.d0*f/h2
430  df(1)=0.d0
431  endif
432  f=f-xflow2
433  endif
434  elseif(lakon(nelem)(6:7).eq.'WE') then
435 !
436 ! weir
437 ! 1-WE-2-WO-3
438 !
439  nelemdown=nint(prop(index+5))
440  h3=v(2,kon(ipkon(nelemdown)+3))
441 !
442 ! default depth for weir is hk
443 !
444  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
445  v(3,node2)=hk
446 !
447  if(h3.lt.p+hk) then
448 !
449 ! output only
450 !
451  v(2,node2)=p+hk
452 !
453 ! Q=f_WE(h1): weir equation
454 !
455 c write(30,*) 'WE: weir equation '
456 c write(30,*)'h1= ',h1,'h2= ',h2,'h3= ',h3,'hk= ',hk
457  f=rho*c*b*(h1-p)**(1.5d0)
458  df(1)=3.d0*f/(2.d0*(h1-p))
459  f=f-xflow
460  df(2)=-1.d0
461  df(3)=0.d0
462  else
463 !
464 ! fake equation
465 !
466 c write(30,*) 'WE: weir equation '
467 c write(30,*)'h1= ',h1,'h2= ',h2,'h3= ',h3,'hk= ',hk
468  numf=1
469  nodef(1)=nodem
470  idirf(1)=3
471  f=prop(index+4)-0.5d0
472  df(1)=1.d0
473  endif
474  elseif(lakon(nelem)(6:7).eq.'WO') then
475 !
476 ! weir opening (element streamdown of weir)
477 ! 0-WE-1-WO-2
478 !
479  nelemup=nint(prop(index+6))
480  node0=kon(ipkon(nelemup)+1)
481  h0=v(2,node0)
482 !
483  p=prop(ielprop(nelemup)+2)
484 !
485 ! default depth for weir is hk
486 !
487  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
488  v(3,node2)=hk
489 !
490  if(h2.lt.p+hk) then
491 !
492 ! bresse between 1 and 2
493 !
494  h1=hk
495 c write(30,*) 'WO: Bresse equation '
496 c write(30,*)'h0= ',h0,'h1= ',h1,'h2= ',h2,'hk= ',hk
497  p=prop(ielprop(nelemup)+2)
498  s0=dasin(p/dsqrt(dl**2+p**2))
499 c write(*,*) 's0=',p,dl,s0
500  sqrts0=dsqrt(1.d0-s0*s0)
501  bresse=.true.
502  else
503 !
504 ! output only
505 !
506  v(2,node1)=h2
507 !
508 ! bresse between 0 and 2
509 !
510 c write(30,*) 'WO: Bresse eqn. between 0 and 2 '
511 c write(30,*)'h0= ',h0,'h1= ',h1,'h2= ',h2,'hk= ',hk
512  nodef(1)=node0
513  h1=h0
514  bresse=.true.
515  endif
516  elseif(lakon(nelem)(6:7).eq.'DS') then
517 !
518 ! discontinuous slope
519 ! 1-DS-2-DO-3
520 !
521  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
522  v(3,node2)=hk
523 !
524  if(h1.gt.hk) then
525  nelemdown=nint(prop(index+8))
526  h3=v(2,kon(ipkon(nelemdown)+3))
527  if(h3.le.hk) then
528 !
529 ! upstream: backwater curve
530 ! downstream: frontwater curve
531 !
532  h2=hk
533  bresse=.true.
534 c write(30,*) 'DS: back/front bresse'
535 c write(30,*)'h1= ',h1,'h2= ',h2,'h3= ',h3
536 !
537 ! for output purposes
538 !
539  v(2,node2)=h2
540  else
541 !
542 ! both curves are backwater curves
543 ! fake equation
544 !
545 c write(30,*) 'DS: back/back fake equation '
546 c write(30,*)'h1= ',h1,'h2= ',h2,'h3= ',h3
547  numf=1
548  nodef(1)=nodem
549  idirf(1)=3
550  f=prop(index+7)-0.5d0
551  df(1)=1.d0
552  endif
553  else
554 !
555 ! both curves are frontwater curves
556 ! fake equation
557 !
558 c write(30,*) 'DS: front/front fake equation '
559 c write(30,*)'h1= ',h1,'h2= ',h2
560  nelemup=nint(prop(index+6))
561  numf=1
562  nodef(1)=kon(ipkon(nelemup)+2)
563  idirf(1)=3
564  f=prop(index+7)-0.5d0
565  df(1)=1.d0
566  endif
567  elseif(lakon(nelem)(6:7).eq.'DO') then
568 !
569 ! discontinuous slope opening
570 ! (element streamdown of discontinuous slope)
571 ! 0-DS-1-DO-2
572 !
573  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
574  v(3,node2)=hk
575 !
576  nelemup=nint(prop(index+6))
577  node0=kon(ipkon(nelemup)+1)
578  h0=v(2,node0)
579 !
580  if(h0.gt.hk) then
581  if(h2.le.hk) then
582 !
583 ! upstream: backwater curve
584 ! downstream: frontwater curve
585 ! bresse between 1 and 2
586 !
587  h1=hk
588 c write(30,*) 'DO: back/front bresse 1-2'
589 c write(30,*)'h0= ',h0,'h1= ',h1,'h2= ',h2
590  bresse=.true.
591  else
592 !
593 ! both curves are backwater curves
594 ! bresse between 0 and 2
595 !
596 c write(30,*) 'DO: back/back bresse 0-2'
597 c write(30,*)'h0= ',h0,'h1= ',h1,'h2= ',h2
598  nodef(1)=node0
599  h1=h0
600  bresse=.true.
601 !
602 ! output purposes
603 !
604  v(2,node1)=(h0+h2)/2.d0
605  endif
606  else
607 !
608 ! both curves are frontwater curves
609 ! bresse between 0 and 2
610 !
611 c write(30,*) 'DO: front/front bresse 0-2'
612 c write(30,*)'h0= ',h0,'h1= ',h1,'h2= ',h2
613  nodef(1)=node0
614  h1=h0
615  bresse=.true.
616 !
617 ! output purposes
618 !
619  v(2,node1)=(h0+h2)/2.d0
620  endif
621  elseif(lakon(nelem)(6:7).eq.'RE') then
622 !
623 ! element upstream of a reservoir
624 ! calculating the critical depth
625 !
626  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
627  v(3,node2)=hk
628  if(h1.ge.hk) then
629 !
630 ! backwater curve
631 !
632  if(h2.lt.hk) h2=hk
633 c write(30,*) 'RE: Bresse downstream equation '
634 c write(30,*) 'h1= ',h1,'h2= ',h2,'hk= ',hk
635  bresse=.true.
636  else
637 !
638 ! frontwater curve
639 !
640  call hns(b,theta,rho,dg,sqrts0,xflow,h1,h1ns)
641  if(h2.le.h1ns) then
642 c write(30,*) 'RE: fake equation '
643 c write(30,*) 'h1= ',h1,'h2= ',h2,'h1ns= ',h1ns
644 !
645 ! fake equation
646 !
647  nelemup=nint(prop(index+6))
648  nodesg=kon(ipkon(nelemup)+2)
649  numf=1
650  nodef(1)=nodesg
651  idirf(1)=3
652 !
653 ! retrieving previous value of eta
654 !
655  index=ielprop(nelemup)
656  if(lakon(nelemup)(6:7).eq.'SG') then
657  f=prop(index+4)-0.5d0
658  elseif(lakon(nelemup)(6:7).eq.'WE') then
659  f=prop(index+4)-0.5d0
660  elseif(lakon(nelemup)(6:7).eq.'DS') then
661  f=prop(index+7)-0.5d0
662  endif
663  df(1)=1.d0
664  else
665 c write(30,*) 'RE: Bresse downstream equation '
666 c write(30,*) 'h1= ',h1,'h2= ',h2,'h1ns= ',h1ns
667  bresse=.true.
668  endif
669  endif
670  elseif(lakon(nelem)(6:7).eq.'CO') then
671 c write(30,*) 'CO: contraction '
672 c write(30,*)'h1= ',h1,'h2= ',h2
673 !
674  call hcrit(xflow,rho,b2,theta,dg,sqrts0,hk)
675  v(3,node2)=hk
676 !
677  if(inv.eq.-1) then
678  if((h1.gt.hk).and.(h2.lt.hk)) then
679  jump=.true.
680  endif
681  else
682  if((h1.lt.hk).and.(h2.gt.hk)) then
683  jump=.true.
684  endif
685  endif
686 !
687 c write(*,*) 'CO ',jump
688 !
689  if(.not.jump) then
690  c1=rho*rho*dg
691  c2=b1*b2*h1*h2
692  df(1)=b1*(2.d0*xflow2+c1*b1*b2*h2**3)
693  df(3)=b2*(2.d0*xflow2+c1*b1*b1*h1**3)
694  f=h1*df(1)-h2*df(3)
695  df(1)=df(1)-3.d0*c1*c2*b1*h1
696  df(3)=3.d0*c1*c2*b1*h2-df(3)
697  df(2)=4.d0*(b1*h1-b2*h2)*xflow
698  endif
699  elseif(lakon(nelem)(6:7).eq.'EL') then
700 c write(30,*) 'EL: enlargement '
701 c write(30,*)'h1= ',h1,'h2= ',h2
702 !
703  call hcrit(xflow,rho,b2,theta,dg,sqrts0,hk)
704  v(3,node2)=hk
705 !
706  if(inv.eq.-1) then
707  if((h1.gt.hk).and.(h2.lt.hk)) then
708  jump=.true.
709  endif
710  else
711  if((h1.lt.hk).and.(h2.gt.hk)) then
712  jump=.true.
713  endif
714  endif
715 !
716 c write(*,*) 'EL ',jump
717 !
718  if(.not.jump) then
719  c1=rho*rho*dg
720  c2=b1*b2*h1*h2
721  df(1)=b1*(2.d0*xflow2+c1*b2*b2*h2**3)
722  df(3)=b2*(2.d0*xflow2+c1*b1*b2*h1**3)
723  f=h1*df(1)-h2*df(3)
724  df(1)=df(1)-3.d0*c1*c2*b2*h1
725  df(3)=3.d0*c1*c2*b2*h2-df(3)
726  df(2)=4.d0*(b1*h1-b2*h2)*xflow
727  endif
728  elseif(lakon(nelem)(6:7).eq.'DR') then
729 c write(30,*) 'DR: drop '
730 c write(30,*)'h1= ',h1,'h2= ',h2
731 !
732  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
733  v(3,node2)=hk
734 !
735  if(inv.eq.-1) then
736  if((h1.gt.hk).and.(h2.lt.hk)) then
737  jump=.true.
738  endif
739  else
740  if((h1.lt.hk).and.(h2.gt.hk)) then
741  jump=.true.
742  endif
743  endif
744 !
745  if(.not.jump) then
746  c1=rho*rho*dg
747  df(1)=2.d0*xflow2+c1*b*b*h2**3
748  df(3)=2.d0*xflow2+c1*b*b*h1*(h1+d)**2
749  f=h1*df(1)-h2*df(3)
750  df(1)=df(1)-c1*b*b*h2*(3.d0*h1+d)*(h1+d)
751  df(3)=3.d0*c1*b*b*h1*h2*h2-df(3)
752  df(2)=4.d0*(h1-h2)*xflow
753  endif
754  elseif(lakon(nelem)(6:7).eq.'ST') then
755 c write(30,*) 'ST: step '
756 c write(30,*)'h1= ',h1,'h2= ',h2
757 !
758  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
759  v(3,node2)=hk
760 !
761  if(inv.eq.-1) then
762  if((h1.gt.hk).and.(h2.lt.hk)) then
763  jump=.true.
764  endif
765  else
766  if((h1.lt.hk).and.(h2.gt.hk)) then
767  jump=.true.
768  endif
769  endif
770 !
771  if(.not.jump) then
772  c1=rho*rho*dg
773  df(1)=2.d0*xflow2+c1*b*b*h2*(h2+d)**2
774  df(3)=2.d0*xflow2+c1*b*b*h1**3
775  f=h1*df(1)-h2*df(3)
776  df(1)=df(1)-3.d0*c1*b*b*h1*h1*h2
777  df(3)=c1*b*b*h1*(3.d0*h2+d)*(h2+d)-df(3)
778  df(2)=4.d0*(h1-h2)*xflow
779  endif
780  elseif(lakon(nelem)(6:7).eq.' ') then
781  bresse=.true.
782 c write(30,*) 'straight: Bresse equation '
783 c write(30,*) 'h1= ',h1,'h2= ',h2
784  endif
785 !
786 ! bresse equation
787 !
788  if((bresse).or.(jump)) then
789 !
790  if(xks.gt.0.d0) then
791 !
792 ! White-Coolebrook
793 !
794 ! hydraulic diameter
795 !
796  d=2.d0*(h1+h2)
797  reynolds=4.d0*xflow/(b*dvi)
798  form_fact=1.d0
799  call friction_coefficient(dl,d,xks,reynolds,form_fact,
800  & friction)
801  endif
802 !
803  if(bresse) then
804  call hcrit(xflow,rho,b,theta,dg,sqrts0,hk)
805  v(3,node2)=hk
806  if(inv.eq.-1) then
807  if((h1.gt.hk).and.(h2.lt.hk)) then
808  jump=.true.
809  endif
810  else
811  if((h1.lt.hk).and.(h2.gt.hk)) then
812  jump=.true.
813  endif
814  endif
815  b1=b
816  b2=b
817  endif
818 !
819 ! geometric data
820 !
821  cth=dcos(theta)
822  tth=dtan(theta)
823 !
824 ! nonprismatic cross section
825 !
826  if(lakon(nelem)(6:7).eq.' ') then
827  dbds=prop(index+7)
828  else
829  dbds=0.d0
830  endif
831 !
832 ! width at water surface
833 !
834  dd1dh1=2.d0*tth
835  dd2dh2=dd1dh1
836  d1=b1+h1*dd1dh1
837  d2=b2+dl*dbds+h2*dd2dh2
838 !
839 ! cross section
840 !
841  a1=h1*(b1+h1*tth)
842  a2=h2*(b2+dl*dbds+h2*tth)
843  da1dh1=d1
844  da2dh2=d2
845 !
846 ! perimeter
847 !
848  p1=b1+2.d0*h1/cth
849  p2=b2+dl*dbds+2.d0*h2/cth
850  dp1dh1=2.d0/cth
851  dp2dh2=dp1dh1
852 !
853 ! factor for friction
854 !
855  if(xks.gt.0.d0) then
856 ! White-Coolebrook
857  um1=friction/8.d0
858  um2=um1
859  dum1dh1=0.d0
860  dum2dh2=0.d0
861  else
862 ! Manning
863  um1=xks*xks*dg*(p1/a1)**(1.d0/3.d0)
864  um2=xks*xks*dg*(p2/a2)**(1.d0/3.d0)
865  dum1dh1=xks*xks*dg*
866  & (p1**(-2.d0/3.d0)*dp1dh1*a1**(1.d0/3.d0)-
867  & a1**(-2.d0/3.d0)*da1dh1*p1**(1.d0/3.d0))/
868  & (3.d0*a1**(2.d0/3d0))
869  dum2dh2=xks*xks*dg*
870  & (p2**(-2.d0/3.d0)*dp2dh2*a2**(1.d0/3.d0)-
871  & a2**(-2.d0/3.d0)*da2dh2*p2**(1.d0/3.d0))/
872  & (3.d0*a2**(2.d0/3d0))
873  endif
874 !
875 ! constants
876 !
877  c1=rho*rho*dg
878  c2=c1*sqrts0
879  c1=c1*s0
880 !
881 ! hydraulic jump
882 !
883  if(jump) then
884 c write(30,*)
885 c & 'liquidchannel: jump in element,hk ',nelem,hk
886  nelemup=prop(index+6)
887  indexup=ielprop(nelemup)
888  if(lakon(nelemup)(6:7).eq.'SG') then
889  eta=prop(indexup+4)
890  prop(indexup+7)=nelem
891  elseif(lakon(nelemup)(6:7).eq.'WE') then
892  eta=prop(indexup+4)
893  prop(indexup+7)=nelem
894  elseif(lakon(nelemup)(6:7).eq.'DS') then
895  eta=prop(indexup+7)
896  prop(indexup+9)=nelem
897  endif
898 !
899 ! determining h3, h4 and derivatives
900 !
901 ! numerator
902 !
903  xt1=c1*a1**3+(h1*dbds-um1*p1)*xflow2
904  xt2=c1*a2**3+(h2*dbds-um2*p2)*xflow2
905 !
906 ! denominator
907 !
908  xn1=c2*a1**3-d1*xflow2
909  xn2=c2*a2**3-d2*xflow2
910 !
911 ! h3 and h4
912 !
913  h3=h1+dl*xt1/xn1*eta
914  h4=h2-dl*xt2/xn2*(1.d0-eta)
915 c write(30,*)
916 c & 'liquidchannel: h3,h4,eta ',h3,h4,eta
917 !
918  if(bresse) then
919 !
920 ! width at jump
921 !
922  bj=b+dbds*eta*dl
923 !
924 ! cross sections and derivatives
925 !
926  a3=h3*(bj+h3*tth)
927  a4=h4*(bj+h4*tth)
928  da3dh3=bj+2.d0*h3*tth
929  da4dh4=bj+2.d0*h4*tth
930 !
931 ! center of gravity and derivatives
932 !
933  yg3=h3*(3.d0*bj+2.d0*h3*tth)/(6.d0*(bj+h3*tth))
934  yg4=h4*(3.d0*bj+2.d0*h4*tth)/(6.d0*(bj+h4*tth))
935  dyg3dh3=((3.d0*bj+4.d0*h3*tth)*(bj+tth)
936  & -tth*h3*(3.d0*bj+2.d0*h3*tth))/
937  & (6.d0*(bj+h3*tth)**2)
938  dyg4dh4=((3.d0*bj+4.d0*h4*tth)*(bj+tth)
939  & -tth*h4*(3.d0*bj+2.d0*h4*tth))/
940  & (6.d0*(bj+h4*tth)**2)
941  dyg3dbj=h3*h3*tth/(6.d0*(bj+h3*tth)**2)
942  dyg4dbj=h4*h4*tth/(6.d0*(bj+h4*tth)**2)
943  endif
944 !
945 ! derivative of h3 w.r.t. h1 and of h4 w.r.t. h2
946 !
947  dh3dh1=1.d0+((3.d0*c1*a1*a1*da1dh1
948  & +(dbds-dum1dh1*p1-um1*dp1dh1)*xflow2)*xn1
949  & -(3.d0*c2*a1*a1*da1dh1-dd1dh1*xflow2)*xt1)/
950  & (xn1*xn1)*eta*dl
951  dh4dh2=1.d0-((3.d0*c1*a2*a2*da2dh2
952  & +(dbds-dum2dh2*p2-um2*dp2dh2)*xflow2)*xn2
953  & -(3.d0*c2*a2*a2*da2dh2-dd2dh2*xflow2)*xt2)/
954  & (xn2*xn2)*(1.d0-eta)*dl
955 !
956  if(bresse) then
957  da3dh1=da3dh3*dh3dh1
958  da4dh2=da4dh4*dh4dh2
959  dyg3dh1=dyg3dh3*dh3dh1
960  dyg4dh2=dyg4dh4*dh4dh2
961  endif
962 !
963 ! derivative of h3 and h4 w.r.t. the mass flow
964 !
965  dh3dm=((dbds*h1-um1*p1)*xn1+d1*xt1)*2.d0*xflow/
966  & (xn1*xn1)*eta*dl
967  dh4dm=-((dbds*h2-um2*p2)*xn2+d2*xt2)*2.d0*xflow/
968  & (xn2*xn2)*(1.d0-eta)*dl
969 !
970  if(bresse) then
971  da3dm=da3dh3*dh3dm
972  da4dm=da4dh4*dh4dm
973  dyg3dm=dyg3dh3*dh3dm
974  dyg4dm=dyg4dh4*dh4dm
975  endif
976 !
977 ! derivative of h3 and h4 w.r.t. eta
978 !
979  dh3deta=dl*xt1/xn1
980  dh4deta=dl*xt2/xn2
981 !
982  if(bresse) then
983  dbjdeta=dbds*dl
984 !
985 ! derivative of A3, A4, yg3 and yg4 w.r.t. eta
986 !
987  da3deta=da3dh3*dh3deta+h3*dbjdeta
988  da4deta=da4dh4*dh4deta+h4*dbjdeta
989  dyg3deta=dyg3dh3*dh3deta+dyg3dbj*dbjdeta
990  dyg4deta=dyg4dh4*dh4deta+dyg4dbj*dbjdeta
991  endif
992 !
993  numf=4
994  nodef(4)=kon(ipkon(nelemup)+2)
995  idirf(4)=3
996 !
997  if(bresse) then
998  f=a4*xflow2+c2*(a3*a3*a4*yg3-a3*a4*a4*yg4)
999  & -a3*xflow2
1000  df(1)=c2*(2.d0*a3*da3dh1*a4*yg3+a3*a3*a4*dyg3dh1
1001  & -da3dh1*a4*a4*yg4)-da3dh1*xflow2
1002  df(2)=2.d0*xflow*(a4-a3)+
1003  & (c2*(2.d0*a3*a4*yg3-a4*a4*yg4)-xflow2)*da3dm+
1004  & (c2*(a3*a3*yg3-2.d0*a3*a4*yg4)+xflow2)*da4dm+
1005  & c2*a3*a3*a4*dyg3dm-c2*a3*a4*a4*dyg4dm
1006  df(3)=c2*(a3*a3*da4dh2*yg3-2.d0*a3*a4*da4dh2*yg4
1007  & -a3*a4*a4*dyg4dh2)+da4dh2*xflow2
1008  df(4)=da4deta*xflow2+
1009  & c2*(2.d0*a3*da3deta*a4*yg3+a3*a3*da4deta*yg3
1010  & +a3*a3*a4*dyg3deta-da3deta*a4*a4*yg4
1011  & -a3*2.d0*a4*da4deta*yg4-a3*a4*a4*dyg4deta)
1012  & -da3deta*xflow2
1013  elseif(lakon(nelem)(6:7).eq.'CO') then
1014  f=b2*h4*(2.d0*xflow2+c2*b1*b1*h3**3)-
1015  & b1*h3*(2.d0*xflow2+c2*b1*b2*h4**3)
1016 ! dfdh3
1017  df(1)=3.d0*b2*h4*c2*b1*b1*h3*h3-
1018  & b1*(2.d0*xflow2+c2*b1*b2*h4**3)
1019 ! dfdh4
1020  df(3)=b2*(2.d0*xflow2+c2*b1*b1*h3**3)-
1021  & 3.d0*b1*h3*c2*b1*b2*h4*h4
1022 ! dfdm
1023  df(2)=4.d0*xflow*(b2*h4-b1*h3)+
1024  & df(1)*dh3dm+df(3)*dh4dm
1025 ! dfdeta
1026  df(4)=df(1)*dh3deta+df(3)*dh4deta
1027 ! dfdh1
1028  df(1)=df(1)*dh3dh1
1029 ! dfdh2
1030  df(3)=df(3)*dh4dh2
1031  elseif(lakon(nelem)(6:7).eq.'EL') then
1032  f=b2*h4*(2.d0*xflow2+c2*b1*b2*h3**3)-
1033  & b1*h3*(2.d0*xflow2+c2*b2*b2*h4**3)
1034 ! dfdh3
1035  df(1)=3.d0*b2*h4*c2*b1*b2*h3*h3-
1036  & b1*(2.d0*xflow2+c2*b2*b2*h4**3)
1037 ! dfdh4
1038  df(3)=b2*(2.d0*xflow2+c2*b1*b2*h3**3)-
1039  & 3.d0*b1*h3*c2*b2*b2*h4*h4
1040 ! dfdm
1041  df(2)=4.d0*xflow*(b2*h4-b1*h3)+
1042  & df(1)*dh3dm+df(3)*dh4dm
1043 ! dfdeta
1044  df(4)=df(1)*dh3deta+df(3)*dh4deta
1045 ! dfdh1
1046  df(1)=df(1)*dh3dh1
1047 ! dfdh2
1048  df(3)=df(3)*dh4dh2
1049  elseif(lakon(nelem)(6:7).eq.'DR') then
1050  f=h4*(2.d0*xflow2+c2*b*b*h3*(h3+d)**2)-
1051  & h3*(2.d0*xflow2+c2*b*b*h4**3)
1052 ! dfdh3
1053  df(1)=h4*c2*b*b*(3.d0*h3+d)*(h3+d)-
1054  & (2.d0*xflow2+c2*b*b*h4**3)
1055 ! dfdh4
1056  df(3)=(2.d0*xflow2+c2*b*b*h3*(h3+d)**2)-
1057  & 3.d0*h3*c2*b*b*h4*h4
1058 ! dfdm
1059  df(2)=4.d0*xflow*(h4-h3)+
1060  & df(1)*dh3dm+df(3)*dh4dm
1061 ! dfdeta
1062  df(4)=df(1)*dh3deta+df(3)*dh4deta
1063 ! dfdh1
1064  df(1)=df(1)*dh3dh1
1065 ! dfdh2
1066  df(3)=df(3)*dh4dh2
1067  elseif(lakon(nelem)(6:7).eq.'ST') then
1068  f=h4*(2.d0*xflow2+c2*b*b*h3**3)-
1069  & h3*(2.d0*xflow2+c2*b*b*h4*(h4+d)**2)
1070 ! dfdh3
1071  df(1)=3.d0*h4*c2*b*b*h3*h3-
1072  & (2.d0*xflow2+c2*b*b*h4*(h4+d)**2)
1073 ! dfdh4
1074  df(3)=(2.d0*xflow2+c2*b*b*h3**3)-
1075  & h3*c2*b*b*(3.d0*h4+d)*(h4+d)
1076 ! dfdm
1077  df(2)=4.d0*xflow*(h4-h3)+
1078  & df(1)*dh3dm+df(3)*dh4dm
1079 ! dfdeta
1080  df(4)=df(1)*dh3deta+df(3)*dh4deta
1081 ! dfdh1
1082  df(1)=df(1)*dh3dh1
1083 ! dfdh2
1084  df(3)=df(3)*dh4dh2
1085  endif
1086  else
1087 !
1088 ! regular Bresse equation
1089 !
1090  f=c2*(a1**3+a2**3)-xflow2*(d1+d2)
1091  df(1)=-f+(h2-h1)*(c2*da1dh1*3.d0*a1*a1-xflow2*dd1dh1)
1092  & -dl*(c1*3.d0*a1*a1*da1dh1
1093  & -(dum1dh1*p1+um1*dp1dh1-dbds)*xflow2)
1094  df(2)=(-(h2-h1)*(d1+d2)
1095  & +dl*(um1*p1+um2*p2-(h1+h2)*dbds))*2.d0*xflow
1096  df(3)=f+(h2-h1)*(c2*da2dh2*3.d0*a2*a2-xflow2*dd2dh2)
1097  & -dl*(c1*3.d0*a2*a2*da2dh2
1098  & -(dum2dh2*p2+um2*dp2dh2-dbds)*xflow2)
1099  f=(h2-h1)*f-dl*(c1*(a1**3+a2**3)
1100  & -(um1*p1+um2*p2-(h1+h2)*dbds)*xflow2)
1101  endif
1102  endif
1103  endif
1104  elseif(iflag.eq.3) then
1105 !
1106 ! only if called from resultnet in case the element contains
1107 ! a hydraulic jump and eta<0 or eta>1. This means that the
1108 ! jump does not take place in the element itself. By adjusting
1109 ! h1 or h2 the jump is forced into a neighboring element
1110 !
1111  index=ielprop(nelem)
1112 c write(30,*) 'iflag=3, nelem',nelem,lakon(nelem)
1113 !
1114  h1=v(2,node1)
1115  h2=v(2,node2)
1116 !
1117  z1=-g(1)*co(1,node1)-g(2)*co(2,node1)-g(3)*co(3,node1)
1118  z2=-g(1)*co(1,node2)-g(2)*co(2,node2)-g(3)*co(3,node2)
1119 !
1120  dg=dsqrt(g(1)*g(1)+g(2)*g(2)+g(3)*g(3))
1121 !
1122  xflow2=xflow*xflow
1123 !
1124 ! determine eta (present location of jump)
1125 !
1126  nelemup=prop(index+6)
1127  indexup=ielprop(nelemup)
1128  if(lakon(nelemup)(6:7).eq.'SG') then
1129  eta=prop(indexup+4)
1130  prop(indexup+4)=0.5d0
1131  prop(indexup+7)=0.d0
1132  elseif(lakon(nelemup)(6:7).eq.'WE') then
1133  eta=prop(indexup+4)
1134  prop(indexup+4)=0.5d0
1135  prop(indexup+7)=0.d0
1136  elseif(lakon(nelemup)(6:7).eq.'DS') then
1137  eta=prop(indexup+7)
1138  prop(indexup+7)=0.5d0
1139  prop(indexup+9)=0.d0
1140  endif
1141 !
1142 ! element properties
1143 !
1144  if((lakon(nelem)(6:7).eq.'SG').or.
1145  & (lakon(nelem)(6:7).eq.'SO').or.
1146  & (lakon(nelem)(6:7).eq.'RE').or.
1147  & (lakon(nelem)(6:7).eq.' ').or.
1148  & (lakon(nelem)(6:7).eq.'DS').or.
1149  & (lakon(nelem)(6:7).eq.'DO')) then
1150  b=prop(index+1)
1151  s0=prop(index+2)
1152  if(s0.lt.-1.d0) then
1153  s0=dasin((z1-z2)/dl)
1154  endif
1155  sqrts0=dsqrt(1.d0-s0*s0)
1156  if(lakon(nelem)(6:7).ne.'SG') then
1157  dl=prop(index+3)
1158  theta=prop(index+4)
1159  xks=prop(index+5)
1160  if(dl.le.0.d0) then
1161  dl=dsqrt((co(1,node2)-co(1,node1))**2+
1162  & (co(2,node2)-co(2,node1))**2+
1163  & (co(3,node2)-co(3,node1))**2)
1164  endif
1165  else
1166  theta=0.d0
1167  endif
1168  elseif(lakon(nelem)(6:7).eq.'WE') then
1169  b=prop(index+1)
1170  p=prop(index+2)
1171  c=prop(index+3)
1172  elseif((lakon(nelem)(6:7).eq.'CO').or.
1173  & (lakon(nelem)(6:7).eq.'EL')) then
1174  b1=prop(index+1)
1175  s0=prop(index+2)
1176  if(s0.lt.-1.d0) then
1177  s0=dasin((z1-z2)/dl)
1178  endif
1179  sqrts0=dsqrt(1.d0-s0*s0)
1180  b2=prop(index+4)
1181  elseif((lakon(nelem)(6:7).eq.'DR').or.
1182  & (lakon(nelem)(6:7).eq.'ST'))then
1183  b=prop(index+1)
1184  s0=prop(index+2)
1185  if(s0.lt.-1.d0) then
1186  s0=dasin((z1-z2)/dl)
1187  endif
1188  sqrts0=dsqrt(1.d0-s0*s0)
1189  d=prop(index+4)
1190  endif
1191 !
1192 ! contraction, enlargement, drop and step:
1193 ! adjust h1 or h2 by solving the appropriate
1194 ! momentum equation
1195 !
1196  if((lakon(nelem)(6:7).eq.'CO').or.
1197  & (lakon(nelem)(6:7).eq.'EL').or.
1198  & (lakon(nelem)(6:7).eq.'DR').or.
1199  & (lakon(nelem)(6:7).eq.'ST'))then
1200  c2=rho*rho*dg*sqrts0
1201 !
1202  if(eta.gt.1.d0) then
1203 !
1204 ! h1 is given, h2 is unknown
1205 !
1206  if(lakon(nelem)(6:7).eq.'CO') then
1207  e3=b1*h1*c2*b1*b2
1208  e0=2.d0*b1*h1*xflow2/e3
1209  e1=-(2.d0*xflow2+c2*b1*b1*h1**3)*b2/e3
1210  e2=0.d0
1211  elseif(lakon(nelem)(6:7).eq.'EL') then
1212  e3=b1*h1*c2*b2*b2
1213  e0=2.d0*b1*h1*xflow2/e3
1214  e1=-(2.d0*xflow2+c2*b1*b2*h1**3)*b2/e3
1215  e2=0.d0
1216  elseif(lakon(nelem)(6:7).eq.'DR') then
1217  e3=h1*c2*b*b
1218  e0=h1*2.d0*xflow2/e3
1219  e1=-(2.d0*xflow2+c2*b*b*h1*(h1+d)**2)/e3
1220  e2=0.d0
1221  elseif(lakon(nelem)(6:7).eq.'ST') then
1222  e3=h1*c2*b*b
1223  e0=h1*2.d0*xflow2/e3
1224  e1=(h1*c2*b*b*d*d-(2.d0*xflow2+c2*b*b*h1**3))/e3
1225  e2=h1*c2*b*b*2.d0*d/e3
1226  endif
1227 !
1228 ! solve the cubic equation
1229 !
1230  call cubic(e0,e1,e2,solreal,solimag,nsol)
1231 !
1232 ! determine the real solution closest to h1
1233 !
1234  dist=1.d30
1235  do i=1,nsol
1236  if(dabs(solreal(i)-h1).lt.dist) then
1237  dist=dabs(solreal(i)-h1)
1238  h2=solreal(i)
1239  endif
1240  enddo
1241  if(nactdog(2,node2).ne.0) v(2,node2)=h2
1242  elseif(eta.lt.0.d0) then
1243 !
1244 ! h2 is given, h1 is unknown
1245 !
1246  if(lakon(nelem)(6:7).eq.'CO') then
1247  e3=c2*b1*b1*b2*h2
1248  e0=2.d0*xflow2*b2*h2/e3
1249  e1=-b1*(2.d0*xflow2+c2*b1*b2*h2**3)/e3
1250  e2=0.d0
1251  elseif(lakon(nelem)(6:7).eq.'EL') then
1252  e3=c2*b1*b2*b2*h2
1253  e0=2.d0*xflow2*b2*h2/e3
1254  e1=-b1*(2.d0*xflow2+c2*b2*b2*h2**3)/e3
1255  e2=0.d0
1256  elseif(lakon(nelem)(6:7).eq.'DR') then
1257  e3=c2*b*b*h2
1258  e0=2.d0*xflow2*h2/e3
1259  e1=(c2*b*b*d*d*h2-(2.d0*xflow2+c2*b*b*h2**3))/e3
1260  e2=c2*b*b*2.d0*d*h2/e3
1261  elseif(lakon(nelem)(6:7).eq.'ST') then
1262  e3=c2*b*b*h2
1263  e0=2.d0*xflow2*h2/e3
1264  e1=-(2.d0*xflow2+c2*b*b*h2*(h2+d)**2)/e3
1265  e2=0.d0
1266  endif
1267 !
1268 ! solve the cubic equation
1269 !
1270  call cubic(e0,e1,e2,solreal,solimag,nsol)
1271 c write(30,*) 'check ',solreal(1)**3+e1*solreal(1)+e0
1272 !
1273 c write(30,*) 'nsol',nsol
1274 c write(30,*) 'solreal',(solreal(i),i=1,3)
1275 c write(30,*) 'solimag',(solimag(i),i=1,3)
1276 !
1277 ! determine the real solution closest to h2
1278 !
1279  dist=1.d30
1280  do i=1,nsol
1281  if(dabs(solreal(i)-h2).lt.dist) then
1282  dist=dabs(solreal(i)-h2)
1283  h1=solreal(i)
1284  endif
1285  enddo
1286  if(nactdog(2,node1).ne.0) v(2,node1)=h1
1287  endif
1288  return
1289  endif
1290 !
1291  if(xks.gt.0.d0) then
1292 !
1293 ! White-Coolebrook
1294 !
1295 ! hydraulic diameter
1296 !
1297  d=2.d0*(h1+h2)
1298  reynolds=4.d0*xflow/(b*dvi)
1299  form_fact=1.d0
1300  call friction_coefficient(dl,d,xks,reynolds,form_fact,
1301  & friction)
1302  endif
1303 !
1304 ! geometric data
1305 !
1306  cth=dcos(theta)
1307  tth=dtan(theta)
1308 !
1309 ! nonprismatic cross section
1310 !
1311  if(lakon(nelem)(6:7).eq.' ') then
1312  dbds=prop(index+7)
1313  else
1314  dbds=0.d0
1315  endif
1316 !
1317 ! width at water surface
1318 !
1319  dd1dh1=2.d0*tth
1320  dd2dh2=dd1dh1
1321  d1=b+h1*dd1dh1
1322  d2=b+dl*dbds+h2*dd2dh2
1323 !
1324 ! cross section
1325 !
1326  a1=h1*(b+h1*tth)
1327  a2=h2*(b+dl*dbds+h2*tth)
1328 !
1329 ! perimeter
1330 !
1331  p1=b+2.d0*h1/cth
1332  p2=b+dl*dbds+2.d0*h2/cth
1333 !
1334 ! factor for friction
1335 !
1336  if(xks.gt.0.d0) then
1337 ! White-Coolebrook
1338  um1=friction/8.d0
1339  um2=um1
1340  else
1341 ! Manning
1342  um1=xks*xks*dg*(p1/a1)**(1.d0/3.d0)
1343  um2=xks*xks*dg*(p2/a2)**(1.d0/3.d0)
1344  endif
1345 !
1346 ! constants
1347 !
1348  c1=rho*rho*dg
1349  c2=c1*sqrts0
1350  c1=c1*s0
1351 !
1352  if(eta.gt.1.d0) then
1353  xt1=c1*a1**3+(h1*dbds-um1*p1)*xflow2
1354  xn1=c2*a2**3-d2*xflow2
1355  if(nactdog(2,node2).ne.0) v(2,node2)=h1+dl*xt1/xn1
1356 c write(30,*) 'move jump: h1 h2,h2new ',h1,h2,v(2,node2)
1357  elseif(eta.lt.0.d0) then
1358  xt2=c1*a2**3+(h2*dbds-um2*p2)*xflow2
1359  xn2=c2*a2**3-d2*xflow2
1360  if(nactdog(2,node1).ne.0)
1361  & v(2,node1)=h2-dl*xt2/xn2
1362 c write(30,*) 'move jump: h1 h1new h2 ',h1,v(2,node1),h2
1363  endif
1364  endif
1365 !
1366  return
z1
static double * z1
Definition: filtermain.c:48
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
c1
static double * c1
Definition: mafillvcompmain.c:30
hns
subroutine hns(b, theta, rho, dg, sqrts0, xflow, h1, h2)
Definition: hns.f:20
friction_coefficient
subroutine friction_coefficient(l, d, ks, reynolds, form_fact, lambda)
Definition: friction_coefficient.f:26
dist
static double * dist
Definition: radflowload.c:42
cubic
subroutine cubic(a0, a1, a2, solreal, solimag, n)
Definition: cubic.f:20
hcrit
subroutine hcrit(xflow, rho, b, theta, dg, sqrts0, hk)
Definition: hcrit.f:20
b1
static double * b1
Definition: mafillkmain.c:30
Hosted by OpenAircraft.com, (Michigan UAV, LLC)