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