!----- AGPL -------------------------------------------------------------------- ! ! Copyright (C) Stichting Deltares, 2015. ! ! This file is part of Delft3D (D-Flow Flexible Mesh component). ! ! Delft3D 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 version 3. ! ! Delft3D 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 Delft3D. If not, see . ! ! contact: delft3d.support@deltares.nl ! Stichting Deltares ! P.O. Box 177 ! 2600 MH Delft, The Netherlands ! ! All indications and logos of, and references to, "Delft3D", ! "D-Flow Flexible Mesh" and "Deltares" are registered trademarks of Stichting ! Deltares, and remain the property of Stichting Deltares. All rights reserved. ! !------------------------------------------------------------------------------- ! $Id: manholes.f90 42642 2015-10-21 11:34:20Z dam_ar $ ! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/manholes.f90 $ !> Manholes connect 1D networks with 2D/3D grids or other 1D networks. !! !! A manhole has an x,y-location, which determines the unique 2D flow cell !! that it lies in. A manhole can also be 'pressurized' (closed lid), then !! no 2D cell is associated, only 1D netlinks can be connected then. module m_manholes use properties use unstruc_messages implicit none integer, parameter :: MANHOLE_CLOSED = 1 integer, parameter :: MANHOLE_RESERVOIR = 2 integer, parameter :: MANHOLE_OPEN_MOMENTUM = 3 integer, parameter :: MANHOLE_OPEN_NOMOMENTUM = 4 type manhole_t double precision :: x !< exact x-coordinate double precision :: y !< exact y-coordinate character(len=40):: name !< name of manhole integer :: itype !< type of manhole (see module parameters) ! More attributes... end type manhole_t type(manhole_t), allocatable :: manholes(:) !< Global array of manholes integer :: nummh !< Current number of manholes in array. contains ! !------------------------------------------------------------------------------ subroutine init_manholes() nummh = 0 end subroutine init_manholes ! !------------------------------------------------------------------------------ !> Adds a new manhole at the end of the global manhole list. subroutine add_manhole(x, y, name, itype) double precision, intent(in) :: x, y !< Location coordinates character(len=*), intent(in) :: name !< Name of manhole integer, intent(in) :: itype !< Type of manhole (see module parameters) integer :: maxmh maxmh = size(manholes) if (nummh >= maxmh) then maxmh = ceiling(1.2*max(1,nummh)) call realloc_manholes(manholes, maxmh, .true.) end if nummh = nummh + 1 manholes(nummh)%x = x manholes(nummh)%y = y manholes(nummh)%name = name manholes(nummh)%itype = itype end subroutine add_manhole ! !------------------------------------------------------------------------------ !> Resized an allocatable array of manholes. subroutine realloc_manholes(mhs, n, keepExisting) type(manhole_t), allocatable, intent(inout) :: mhs(:) !< Array of manholes. integer, intent(in) :: n !< new size logical, intent(in) :: keepExisting !< Preserve existing manholes in array or not (set to .false. to save on copy overhead). type(manhole_t), allocatable :: mhb(:) integer :: nold, ncopy logical :: equalSize ! 1. Backup if necessary. if (allocated(mhs)) then nold = ubound(mhs,1) equalSize = nold == n if (equalSize .and. keepExisting) return ! output=input ! if (keepExisting) then ncopy = min(n, nold) allocate (mhb(1:ncopy)) mhb(1:ncopy) = mhs(1:ncopy) ! TODO: If manhole_t will contain allocatables/arrays itself, make allocate and copy more sophisticated. endif if (.not.equalSize) deallocate(mhs) endif ! Reallocate if necessary. if (.not. allocated(mhs)) then allocate(mhs(1:n)) end if ! Put back backup if necessary. if (allocated(mhb)) then mhs(1:ncopy) = mhb(1:ncopy) deallocate(mhb) endif end subroutine realloc_manholes ! !------------------------------------------------------------------------------ !> Deletes all manholes and deallocates global manhole array. subroutine delete_manholes() if (allocated(manholes)) then deallocate(manholes) end if nummh = 0 end subroutine delete_manholes ! !------------------------------------------------------------------------------ !> Reads manholes from file. subroutine load_manholes(filename, jadoorladen) character(len=*), intent(in) :: filename integer, intent(in) :: jadoorladen !< Append to existing manholespoints or not type(tree_data), pointer :: mhs_ptr, mh_ptr integer :: istat, i, itype logical :: success double precision :: x, y call tree_create(trim(filename), mhs_ptr) call prop_file('ini',trim(filename),mhs_ptr,istat) if (istat /= 0) then call mess(LEVEL_ERROR, 'Manhole file '''//trim(filename)//''' not read. Code: ', istat) return else call mess(LEVEL_DEBUG, 'Opened manhole file : ', trim(filename) ) endif do i = 1,size(mhs_ptr%child_nodes) mh_ptr => mhs_ptr%child_nodes(i)%node_ptr if (trim(tree_get_name(mh_ptr)) == 'manhole') then call prop_get_double ( mh_ptr, '*', 'x', x, success) if (.not. success) then cycle end if call prop_get_double ( mh_ptr, '*', 'y', y, success) if (.not. success) then cycle end if call prop_get_integer ( mh_ptr, '*', 'Type', itype, success) if (.not. success) then cycle end if end if call add_manhole(x, y, "Manhole", itype) end do end subroutine load_manholes ! !------------------------------------------------------------------------------ end module m_manholes subroutine flgsfurufm(formno, m, teken, husb, hdsb, velhght, zs, ds, dg, dc, wstr,& & cwfa, cwd, mugfa, cgfa, cgda, strdamf, lambda) !!--copyright------------------------------------------------------------------- ! Copyright (c) 2003, Deltares. All rights reserved. !!--disclaimer------------------------------------------------------------------ ! This code is part of the Delft3D software system. Deltares has ! developed c.q. manufactured this code to its best ability and according to the ! state of the art. Nevertheless, there is no express or implied warranty as to ! this software whether tangible or intangible. In particular, there is no ! express or implied warranty as to the fitness for a particular purpose of this ! software, whether tangible or intangible. The intellectual property rights ! related to this software code remain with Deltares at all times. ! For details on the licensing agreement, we refer to the Delft3D software ! license and any modifications to this license, if applicable. These documents ! are available upon request. !!--version information--------------------------------------------------------- ! $Author$ ! $Date$ ! $Revision$ !!--description----------------------------------------------------------------- ! NONE !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- ! use m_GlobalParameters ! use cpluv use m_flowgeom, only : dx use m_flow use m_flowtimes implicit none ! ! Local parameters ! double precision, parameter :: relax = 0.0D0, alfa = 0.9D0 ! ! Global variables ! integer, intent(in) :: formno integer, intent(in) :: m double precision, intent(in) :: cgda double precision, intent(in) :: cgfa double precision, intent(in) :: cwd double precision, intent(in) :: cwfa double precision, intent(in) :: dc double precision, intent(in) :: dg double precision, intent(in) :: ds double precision, intent(in) :: hdsb double precision, intent(in) :: husb double precision, intent(in) :: mugfa double precision, intent(in) :: strdamf, lambda double precision, intent(in) :: teken double precision, intent(in) :: velhght double precision, intent(in) :: wstr double precision, intent(in) :: zs ! ! ! Local variables ! integer :: itgenstr logical :: again double precision :: cu double precision :: dh double precision :: dsqrt double precision :: dxdt double precision :: hs1 double precision :: mu double precision :: rhsc double precision :: ustru double precision :: su double precision :: sd logical, external :: iterfurufm ! ! !! executable statements ------------------------------------------------------- ! ! !======================================================================= ! Deltares ! One-Two Dimensional Modelling System ! S O B E K ! ! Subsystem: Flow Module ! ! Programmer: J.Kuipers ! ! Module: FLGSFURU (FLow General Structure ! calculate FU and RU) ! ! Module description: The linearization coefficients FU and RU are ! calculated for the general structure ! ! The stage of the flow was already determined. ! ! ! Parameters: ! NR NAME IO DESCRIPTION ! 16 cgda I Contraction coefficient for drowned gate flow ! (adapted) ! 15 cgfa I Contraction coefficient for gate flow ! (adapted) ! 13 cwd I Contraction coefficient for drowned weir flow. ! 12 cwfa I Contraction coefficient for free weir flow. ! (adapted) ! 10 dc I Critical water level (free gate flow) ! 9 dg I Gate opening height. ! 8 ds I Water level immediately downstream the gate. ! 1 formno I Flow condition of general structure: ! 0 : closed or dry ! 1 : free weir flow ! 2 : drowned weir flow ! 3 : free gate flow ! 4 : drowned gate flow ! 5 hdsb I Downstream water level. ! 4 husb I Upstream water level. ! 2 m I Grid index of structure ! 14 mugfa I Vertical contraction coefficient for free gate ! flow (adapted) ! 3 teken I Flow direction (+1/-1). ! 6 velhght I Velocity height ! 11 wstr I Width at centre of structure. ! 7 zs I Bed level at centre of structure. !======================================================================= ! Include Pluvius data space ! Declaration of parameters: ! Declaration of local variables: ! if (formno==0) then ! closed or dry ! hu(m) = 0d0 au(m) = 0d0 ; fu(m) = 0d0 ; ru(m) = 0d0 else again = .true. endif ! ! Calculate upstream energy level w.r.t sill ! hs1 = husb + velhght - zs ! itgenstr = 0 dxdt = strdamf*dx(m)/dts do while (again) itgenstr = itgenstr + 1 if (formno==1) then ! free weir flow cu = cwfa**2*ag /1.5D0 !TEM WRITE (11,*) cu,cwfa au(m) = wstr*hs1*2.0D0/3.0D0 ustru = cwfa*sqrt(ag*2.0D0/3.0D0*hs1) rhsc = cu*(hdsb + velhght - zs)*teken elseif (formno==2) then ! drowned weir flow cu = cwd**2*2.0D0*ag au(m) = wstr*ds dh = max(hs1 - ds, 0.D0) ustru = cwd*dsqrt(ag*2.0D0*dh) rhsc = cu*(hdsb + velhght - (ds + zs))*teken elseif (formno==3) then ! free gate flow mu = mugfa*cgfa cu = mu**2*2.0D0*ag au(m) = wstr*dg dh = max(hs1 - dc, 0.D0) ustru = mu*dsqrt(ag*2.0D0*dh) rhsc = cu*(hdsb + velhght - (dc + zs))*teken elseif (formno==4) then ! drowned gate flow mu = mugfa*cgda cu = mu**2*2.0D0*ag au(m) = wstr*dg dh = max(hs1 - ds, 0.D0) ustru = mu*dsqrt(ag*2.0D0*dh) rhsc = cu*(hdsb + velhght - (ds + zs))*teken else endif ! dads(m) = wstr if (teken>0) then su = husb sd = hdsb else sd = husb su = hdsb endif again = iterfurufm(m, su, sd, ustru, cu, rhsc, dxdt, lambda) .and. itgenstr<100 enddo q1(m) = au(m)*u1(m) ! this may be done without end subroutine flgsfurufm logical function iterfurufm(m, su, sd, ustru, cu, rhsc, dxdt, lambda) !======================================================================= ! Deltares ! One-Two Dimensional Modelling System ! S O B E K ! ! Subsystem: Flow Module ! ! Programmer: J. Noort, Guus Stelling ! ! Module: iterfurufm (ITERFURU) ! ! Module description: coefficients for momentum equation in wet weir point ! ! ! update information ! person date ! ! ! ! Include Pluvius data space ! ! use m_GlobalParameters ! use cpluv use m_strucs use m_flow use m_flowgeom, only : dx implicit none ! ! Global variables ! ! integer, intent(in) :: m double precision, intent(in) :: ustru, lambda double precision, intent(in) :: cu double precision, intent(in) :: rhsc double precision :: su ! not s(up) but s(k1) double precision :: sd ! not s(do) but s(k2), see switch in calling routine double precision :: dxdt ! ! Local variables ! ! double precision, parameter :: relax = 0d0 double precision :: bu double precision :: du double precision :: u1mi, dxfrL ! !! executable statements -------------------------------- ----------------------- ! dxfrL = 0d0 if (lambda == 0) then ! if structure defined friction == 0, use standard friction dxfrL = dx(m)*cfuhi(m) endif bu = dxdt + (1.0 + relax + dxfrL)*ustru du = (strucalfa*q1(m)/max(au(m),1d-4) + (1-strucalfa)*u0(m))*dxdt + relax*ustru*u1(m) + rhsc fu(m) = cu/bu ru(m) = du/bu u1mi = u1(m) u1(m) = ru(m) + fu(m)*(su - sd) if (relax == 0.0) then iterfurufm = .false. else if (abs(u1mi - u1(m))>1.0D-6) then iterfurufm = .true. else iterfurufm = .false. endif end function iterfurufm subroutine flgsfm( n, ng, L, firstiter, jarea) use m_flowgeom !!--copyright------------------------------------------------------------------- ! Copyright (c) 2003, Deltares. All rights reserved. !!--disclaimer------------------------------------------------------------------ ! This code is part of the Delft3D software system. Deltares has ! developed c.q. manufactured this code to its best ability and according to the ! state of the art. Nevertheless, there is no express or implied warranty as to ! this software whether tangible or intangible. In particular, there is no ! express or implied warranty as to the fitness for a particular purpose of this ! software, whether tangible or intangible. The intellectual property rights ! related to this software code remain with Deltares at all times. ! For details on the licensing agreement, we refer to the Delft3D software ! license and any modifications to this license, if applicable. These documents ! are available upon request. !!--version information--------------------------------------------------------- ! $Author$ ! $Date$ ! $Revision$ !!--description----------------------------------------------------------------- ! NONE !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- ! use cpluv ! use m_strucs ! use ident use m_strucs use m_flow implicit none ! ! Global variables ! integer, intent(in) :: n !< general structure point n integer, intent(in) :: ng !< is a member of general structure sigal ng integer, intent(in) :: L !< Flow link number, signed! If L < 0 then flow link is in opposite direction than structure left-right orientation. logical, intent(in) :: firstiter logical :: jarea ! ! ! Local variables ! integer :: il, ir, k1, k2, kL, kR, m, Lf integer :: L0 integer :: errornumber integer, save :: m19=0, m48=0 integer :: ker logical :: velheight real :: abran character(8) :: string1 double precision :: cgd, auL double precision :: cgf double precision :: crest double precision :: cwd double precision :: cwf double precision :: dg double precision :: ds double precision :: ds1 double precision :: ds2 double precision :: hdsb double precision :: husb double precision :: lambda double precision :: mugf double precision :: relax double precision :: rhoast double precision :: rholeft double precision :: rhoright double precision :: strdamf double precision :: teken, tekenstr double precision :: ud double precision :: uu double precision :: w2 double precision :: wsd double precision :: wstr double precision :: zb2 double precision :: zs, gateloweredgelevel, gatedoorheight double precision :: Fusa, Rusa, Ausa, DsL ! !! executable statements ------------------------------------------------------- ! ! !======================================================================= ! Deltares ! One-Two Dimensional Modelling System ! S O B E K ! ! Subsystem: Flow Module ! ! Programmer: J.Kuipers ! ! Module: FLGS (FLow General Structure) ! ! Module description: In subroutine FLGS the QH-relationship for a ! general structure will be transformed to a ! linearized equation ! ! ! Parameters: ! NR NAME IO DESCRIPTION ! 2 il I Grid point on left side of structure (lower ! index). ! 3 ir I Grid point on right side of structure (upper ! index). ! 8 jarea I If True then claculate only area ! 1 m I Grid index of structure ! 4 istru I Number of structure. ! 7 firstiter I True in case of first iteration step. ! Subprogram calls: ! NAME DESCRIPTION ! flgtar FLow get General sTructure ARguments ! flupdg FLow UP/Downstream near General structure ! errmsg generate ERRer MeSsaGe to log file ! flqhgs FLow QH relation for General Structure !======================================================================= ! Include Pluvius data space ! Include identifiers of objects ! ! Declaration of parameters: ! ! ! Declaration of local variables: ! ! ! !TEM WRITE (11,*) 'Call structure',istru,'(',m,il,ir,istru,')' Lf = abs(L) ! NOTE: Since a single general structure may be crossed by multiple flow links, ! pay attention to proper directions: structure parameters are typically determined ! by the structure's left-right direction, whereas upwinding and furu-computations ! are typically in the flow link's 1-2 direction. k1 = ln(1,Lf) ; k2 = ln(2,Lf) ! 1 -> 2 flow link direction kL = kcgen(1,n) ; kR = kcgen(2,n) ! L -> R structure direction m = L il = k1 ir = k2 L0 = n - L1cgensg(ng)+1 zs = min ( bob(1,Lf), bob(2,Lf) ) ! == zcgen(3*ng - 2) crest/silllevel gateloweredgelevel = generalstruc(ng)%gateheightonlink(L0) ! == zcgen(3*ng - 1) under gate door and infinity in open part. ! TODO: RTC: AvD/Herman: hier ook wu'tjes en zb1-tjes etc gaan zetten, voordat we de flupd/flgtar-subroutines gaan callen? ! Velheight is always true for river structures ! velheight = istrtyp(7, istru)==1 velheight = .true. relax = 1.0D0 au(Lf) = 0d0 ; fu(Lf) = 0d0 ; ru(Lf) = 0d0 ! ! ng instead of istru dg = gateloweredgelevel - zs call flupdofm(m, il, ir, ng, velheight, rholeft, rhoright, crest, husb, hdsb, & uu, ud, teken, relax) gatedoorheight = 0d0 tekenstr = teken*sign(1, L) ! if flow link abs(L) is in opposite orientation to the structure's orientation, then negate the just computed upwind (flow) teken. if (husb > zs) then call flgtarfm(ng, L0, wu(Lf), bl(kL), bl(kR), tekenstr, zs, wstr, w2, wsd, zb2, dg, ds1, ds2, cgf, cgd, & cwf, cwd, mugf, lambda, strdamf, gatedoorheight) DsL = s1(k2) - s1(k1) u1(Lf) = Rusav(1,n) - Fusav(1,n)*DsL ; u0(Lf) = u1(Lf) ; q1(Lf) = Ausav(1,n)*u1(Lf) call flqhgsfm(Lf, teken, husb, hdsb, uu, zs, wstr, w2, wsd, zb2, ds1, ds2, dg, & cgf, cgd, cwf, cwd, mugf, lambda, strdamf, jarea, ds) fusav(1,n) = fu(Lf) ; Rusav(1,n) = ru(Lf) ; Ausav(1,n) = au(Lf) endif if (gatedoorheight > 0d0) then ! now add water overflowing top of gate zs = gateloweredgelevel + gatedoorheight if (husb > zs) then ! husb = upwind waterlevel instead of height dg = 1d9 ! sky is the limit, this gate fully open u1(Lf) = rusav(2,n) - Fusav(2,n)*dsL ; u0(Lf) = u1(Lf) ; q1(Lf) = Ausav(2,n)*u1(Lf) call flgtarfm(ng, L0, wu(Lf), bl(kL), bl(kR), tekenstr, zs, wstr, w2, wsd, zb2, dg, ds1, ds2, cgf, cgd, & cwf, cwd, mugf, lambda, strdamf, gatedoorheight) call flqhgsfm(Lf, teken, husb, hdsb, uu, zs, wstr, w2, wsd, zb2, ds1, ds2, dg, & cgf, cgd, cwf, cwd, mugf, lambda, strdamf, jarea, ds) fusav(2,n) = fu(Lf) ; rusav(2,n) = ru(Lf) ; ausav(2,n) = au(Lf) au(Lf) = ausav(1,n) + ausav(2,n) fu(Lf) = (fusav(1,n)*ausav(1,n) + fusav(2,n)*ausav(2,n) ) / au(Lf) ru(Lf) = (rusav(1,n)*ausav(1,n) + rusav(2,n)*ausav(2,n) ) / au(Lf) else fusav(2,n) = 0d0 ; rusav(2,n) = 0d0 ; ausav(2,n) = 0d0 endif endif if (au(Lf) == 0d0) then hu(Lf) = 0d0 endif ! TEMP = laatste statement ! strhis(15, istru) = ds + crest ! waterlevel on crest end subroutine flgsfm subroutine flgsareafm(formno, m, husb, velhght, zs, ds, dg, wstr) !!--copyright------------------------------------------------------------------- ! Copyright (c) 2003, Deltares. All rights reserved. !!--disclaimer------------------------------------------------------------------ ! This code is part of the Delft3D software system. Deltares has ! developed c.q. manufactured this code to its best ability and according to the ! state of the art. Nevertheless, there is no express or implied warranty as to ! this software whether tangible or intangible. In particular, there is no ! express or implied warranty as to the fitness for a particular purpose of this ! software, whether tangible or intangible. The intellectual property rights ! related to this software code remain with Deltares at all times. ! For details on the licensing agreement, we refer to the Delft3D software ! license and any modifications to this license, if applicable. These documents ! are available upon request. !!--version information--------------------------------------------------------- ! $Author$ ! $Date$ ! $Revision$ !!--description----------------------------------------------------------------- ! NONE !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- ! use cpluv use m_flow, only: au implicit none ! ! Global variables ! integer, intent(in) :: formno integer, intent(in) :: m double precision, intent(in) :: dg double precision, intent(in) :: ds double precision, intent(in) :: husb double precision, intent(in) :: velhght double precision, intent(in) :: wstr double precision, intent(in) :: zs ! ! ! Local variables ! double precision :: hs1 ! ! !! executable statements ------------------------------------------------------- ! ! !======================================================================= ! Deltares ! One-Two Dimensional Modelling System ! S O B E K ! ! Subsystem: Flow Module ! ! Programmer: J.Kuipers ! ! Module: FLGSAREA (FLow General Structure ! calculate AREA thru structure) ! ! Module description: The area through the general structure will ! be deermined. ! ! The stage of the flow was already determined. ! ! ! Parameters: ! NR NAME IO DESCRIPTION ! 7 dg I Gate opening height. ! 6 ds I Water level immediately downstream the gate. ! 1 formno I Flow condition of general structure: ! 0 : closed or dry ! 1 : free weir flow ! 2 : drowned weir flow ! 3 : free gate flow ! 4 : drowned gate flow ! 3 husb I Upstream water level. ! 2 m I Grid index of structure ! 4 velhght I Velocity height ! 8 wstr I Width at centre of structure. ! 5 zs I Bed level at centre of structure. !======================================================================= ! Include Pluvius data space ! Declaration of parameters: ! Declaration of local variables: ! if (formno==0) then ! closed or dry ! au(m) = 0.0 ! kfu(m) = 0 else ! ! Calculate upstream energy level w.r.t sill ! hs1 = husb + velhght - zs ! kfu(m) = 1 ! if (formno==1) then ! free weir flow au(m) = wstr*hs1*2.0D0/3.0D0 elseif (formno==2) then ! drowned weir flow au(m) = wstr*ds elseif (formno==3) then ! free gate flow au(m) = wstr*dg elseif (formno==4) then ! drowned gate flow au(m) = wstr*dg else endif endif end subroutine flgsareafm subroutine flgsd2fm(wsd, wstr, zs, w2, zb2, dg, ds1, ds2, elu, hd, rhoast, & & cgd, imag, ds, lambda) !!--copyright------------------------------------------------------------------- ! Copyright (c) 2003, Deltares. All rights reserved. !!--disclaimer------------------------------------------------------------------ ! This code is part of the Delft3D software system. Deltares has ! developed c.q. manufactured this code to its best ability and according to the ! state of the art. Nevertheless, there is no express or implied warranty as to ! this software whether tangible or intangible. In particular, there is no ! express or implied warranty as to the fitness for a particular purpose of this ! software, whether tangible or intangible. The intellectual property rights ! related to this software code remain with Deltares at all times. ! For details on the licensing agreement, we refer to the Delft3D software ! license and any modifications to this license, if applicable. These documents ! are available upon request. !!--version information--------------------------------------------------------- ! $Author$ ! $Date$ ! $Revision$ !!--description----------------------------------------------------------------- ! NONE !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- implicit none ! ! Local parameters ! double precision, parameter :: c23 = 2.0D0/3.0D0, c13 = 1.0D0/3.0D0 ! ! Global variables ! logical, intent(out) :: imag double precision, intent(in) :: cgd double precision, intent(in) :: dg double precision, intent(out) :: ds double precision, intent(in) :: ds1 double precision, intent(in) :: ds2 double precision, intent(in) :: elu double precision, intent(in) :: hd double precision, intent(in) :: lambda double precision, intent(in) :: rhoast double precision, intent(in) :: w2 double precision, intent(in) :: wsd double precision, intent(in) :: wstr double precision, intent(in) :: zb2 double precision, intent(in) :: zs ! ! ! Local variables ! double precision :: ag double precision :: bg double precision :: cg double precision :: d2 double precision :: det double precision :: hsl double precision :: terma double precision :: termb ! ! !! executable statements ------------------------------------------------------- ! ! !======================================================================= ! Deltares ! One-Two Dimensional Modelling System ! S O B E K ! ! Subsystem: Flow Module ! ! Programmer: J.Brouwer/J.Kuipers ! ! Module: FLGSD2 (FLow Gen. Struct. Depth sill 2nd ord. eq.) ! ! Module description: Compute water depth ds at the sill by a second ! order algebraic equation. ! ! In case of drowned gate flow the water level at ! the sill is required. The water depth is calcu- ! lated in this routine. ! ! ! Parameters: ! NR NAME IO DESCRIPTION ! 12 cgd I Correction coefficient for drowned gate flow. ! 6 dg I Gate opening height. ! 7 ds1 I Delta s1 general structure. ! 8 ds2 I Delta s2 general structure. ! 14 ds IO Water level immediately downstream the gate. ! 9 elu I Upstream energy level. ! 10 hd I Downstream water level. ! 13 imag O Logical indicator, = TRUE when determinant of ! second order algebraic equation less than ! zero. ! 15 lambda I Extra resistance in general structure. ! 11 rhoast I Downstream water density divided by upstream ! water density. ! 4 w2 I Width at right side of structure. ! 1 wsd I Width structure right or left side. ! 2 wstr I Width at centre of structure. ! 5 zb2 I Bed level at right side of structure. ! 3 zs I Bed level at centre of structure. !======================================================================= ! ! Declaration of parameters: ! ! ! Declaration of local variables: ! ! !JK LOGICAL uitput !JK COMMON /UITPUT/uitput ! ! Calculate Ag, Bg and Cg according to appendix C of ! the design document River Rural integratietraject deel 3. !JK WRITE (11,*) 'IN FLGSD2 ----' ! ag = (1.0D0 - rhoast)*(w2/12.0D0 + wsd/4.0D0) + 0.5D0*(rhoast + 1.0D0) & & *(c13*w2 + c23*wsd) d2 = hd - zb2 ! terma = (4.0D0*rhoast*cgd*cgd*dg*dg*wstr*wstr)/(w2*d2)*(1.0D0 + lambda/d2) termb = 4.0D0*cgd*dg*wstr ! bg = (1.0D0 - rhoast)*((d2 + ds1)*(w2 + wsd)/6.D0 + ds1*wsd*c13) & & + 0.5D0*(rhoast + 1.0D0) & & *((ds1 + ds2 - d2)*(c13*w2 + c23*wsd) + (c23*d2 + c13*ds1) & & *w2 + (c13*d2 + c23*ds1)*wsd) + terma - termb ! hsl = elu - zs ! cg = (1.0D0 - rhoast)*((d2 + ds1)**2*(w2 + wsd)/12.D0 + ds1**2*wsd/6.0D0) & & + 0.5D0*(rhoast + 1.0D0)*(ds1 + ds2 - d2) & & *((c23*d2 + c13*ds1)*w2 + (c13*d2 + c23*ds1)*wsd) - terma*hsl + & & termb*hsl ! det = bg*bg - 4.0D0*ag*cg if (det<0.0D0) then imag = .true. !JK WRITE (11,*) 'Det=',det else imag = .false. ds = ( - bg + sqrt(det))/(2.0D0*ag) endif end subroutine flgsd2fm subroutine flgsd3fm(wsd, wstr, zs, w2, zb2, ds1, ds2, elu, hd, rhoast, cwd, & & ds, lambda) !!--copyright------------------------------------------------------------------- ! Copyright (c) 2003, Deltares. All rights reserved. !!--disclaimer------------------------------------------------------------------ ! This code is part of the Delft3D software system. Deltares has ! developed c.q. manufactured this code to its best ability and according to the ! state of the art. Nevertheless, there is no express or implied warranty as to ! this software whether tangible or intangible. In particular, there is no ! express or implied warranty as to the fitness for a particular purpose of this ! software, whether tangible or intangible. The intellectual property rights ! related to this software code remain with Deltares at all times. ! For details on the licensing agreement, we refer to the Delft3D software ! license and any modifications to this license, if applicable. These documents ! are available upon request. !!--version information--------------------------------------------------------- ! $Author$ ! $Date$ ! $Revision$ !!--description----------------------------------------------------------------- ! NONE !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- implicit none ! ! Local parameters ! double precision, parameter :: c23 = 2.0D0/3.0D0, c13 = 1.0D0/3.0D0 ! ! Global variables ! double precision, intent(in) :: cwd double precision, intent(out) :: ds double precision, intent(in) :: ds1 double precision, intent(in) :: ds2 double precision, intent(in) :: elu double precision, intent(in) :: hd double precision, intent(in) :: lambda double precision, intent(in) :: rhoast double precision, intent(in) :: w2 double precision, intent(in) :: wsd double precision, intent(in) :: wstr double precision, intent(in) :: zb2 double precision, intent(in) :: zs ! ! ! Local variables ! double precision :: aw double precision :: bw double precision :: cw double precision :: d2 double precision :: fac double precision :: h2a double precision :: h2b double precision :: h2c double precision :: hsl double precision :: hulp double precision :: hulp1 double precision :: p double precision :: phi double precision :: q double precision :: r60 double precision :: term double precision :: u double precision :: v ! ! !! executable statements ------------------------------------------------------- ! ! !======================================================================= ! Deltares ! One-Two Dimensional Modelling System ! S O B E K ! ! Subsystem: Flow Module ! ! Programmer: J.Brouwer/J.Kuipers ! ! Module: FLGSD3 (FLow Gen. Struct. Depth sill 3rd ord. eq.) ! ! Module description: Compute water depth ds at the sill by solving a ! third order algebraic equation. ! ! In case of drowned weir flow the water level at ! the sill is required. The water depth is calcu- ! lated in this routine. ! ! ! Parameters: ! NR NAME IO DESCRIPTION ! 11 cwd I Correction coefficient for drowned weir flow. ! 6 ds1 I Delta s1 general structure. ! 7 ds2 I Delta s2 general structure. ! 12 ds IO Water level immediately downstream the gate. ! 8 elu I Upstream energy level. ! 9 hd I Downstream water level. ! 13 lambda I Extra resistance in general structure. ! 10 rhoast I Downstream water density divided by upstream ! water density. ! 4 w2 I Width at right side of structure. ! 1 wsd I Width structure right or left side. ! 2 wstr I Width at centre of structure. ! 5 zb2 I Bed level at right side of structure. ! 3 zs I Bed level at centre of structure. !======================================================================= ! ! Declaration of parameters: ! ! ! Declaration of local variables: ! ! !JK LOGICAL uitput !JK COMMON /UITPUT/uitput ! ! Calculate Dw (=term), Aw, Bw and Cw according to appendix C of ! the design document River Rural integratietraject deel 3. ! !JK WRITE (11,*) 'FLGSD3' !JK WRITE (11,*) 'wsd,wstr,zs ,w2 ,zb2,ds1 ,ds2 ,elu ,hd' , !JK + wsd,wstr,zs ,w2 ,zb2,ds1 ,ds2 ,elu ,hd d2 = hd - zb2 hsl = elu - zs !JK WRITE (11,*) 'hsl',hsl term = ((4.0D0*cwd*cwd*rhoast*wstr*wstr)/(w2*d2))*(1.0D0 + lambda/d2) ! aw = ( - term*hsl - 4.0D0*cwd*wstr + (1.0D0 - rhoast) & & *(w2/12.0D0 + wsd/4.0D0) + 0.5D0*(rhoast + 1.0D0)*(c13*w2 + c23*wsd)) & & /term ! bw = (4.0D0*cwd*wstr*hsl + (1.0D0 - rhoast) & & *((d2 + ds1)*(w2 + wsd)/6.D0 + ds1*wsd*c13) + 0.5D0*(rhoast + 1.0D0) & & *((ds1 + ds2 - d2)*(c13*w2 + c23*wsd) + (c23*d2 + c13*ds1) & & *w2 + (c13*d2 + c23*ds1)*wsd))/term ! cw = ((1.0D0 - rhoast)*((d2 + ds1)**2*(w2 + wsd)/12.D0 + ds1**2*wsd/6.0D0) & & + 0.5D0*(rhoast + 1.0D0)*(ds1 + ds2 - d2) & & *((c23*d2 + c13*ds1)*w2 + (c13*d2 + c23*ds1)*wsd))/term ! ! Solve the equation ds**3 + aw*ds**2 + bw*ds +cw to get the water ! level at the sill ! p = bw/3.0D0 - aw*aw/9.0D0 q = aw*aw*aw/27.0D0 - aw*bw/6.0D0 + cw/2.0D0 hulp = q*q + p*p*p ! if (hulp<0.0D0) then p = abs(p) phi = acos(abs(q)/p/sqrt(p))/3.0D0 r60 = acos(0.5D0) fac = sign(2.D0, q)*sqrt(p) h2a = -fac*cos(phi) h2b = fac*cos(r60 - phi) h2c = fac*cos(r60 + phi) ds = max(h2a, h2b, h2c) - aw/3.0D0 else hulp = sqrt(hulp) hulp1 = -q + hulp if (abs(hulp1)<1E-6) then u = 0 ; v = 0 else ! hk: ook fix for Erwin, ARS 15132 u = abs(hulp1)**c13*sign(1.0D0, hulp1) hulp1 = -q - hulp v = abs(hulp1)**c13*sign(1.0D0, hulp1) endif ds = u + v - aw/3.0D0 endif end subroutine flgsd3fm subroutine flgtarfm(ng, L0, wuL, bl1, bl2, teken, zs, wstr, w2, wsd, zb2, dg, ds1, ds2, cgf, & ! fromgeneral cgd, cwf, cwd, mugf, lambda, strdamf, gatedoorheight) !!--copyright------------------------------------------------------------------- ! Copyright (c) 2003, Deltares. All rights reserved. !!--disclaimer------------------------------------------------------------------ ! This code is part of the Delft3D software system. Deltares has ! developed c.q. manufactured this code to its best ability and according to the ! state of the art. Nevertheless, there is no express or implied warranty as to ! this software whether tangible or intangible. In particular, there is no ! express or implied warranty as to the fitness for a particular purpose of this ! software, whether tangible or intangible. The intellectual property rights ! related to this software code remain with Deltares at all times. ! For details on the licensing agreement, we refer to the Delft3D software ! license and any modifications to this license, if applicable. These documents ! are available upon request. !!--version information--------------------------------------------------------- ! $Author$ ! $Date$ ! $Revision$ !!--description----------------------------------------------------------------- ! NONE !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- ! use cpluv use m_strucs use m_missing implicit none ! ! Global variables ! integer :: ng integer, intent(in) :: L0 !< counter for the current flow link under genstru #ng (1:ncgen for each separate genstru) double precision, intent(in) :: wuL !< wu of this flow link. double precision, intent(in) :: bl1 !< bl of nod1 double precision, intent(in) :: bl2 !< bl of nod2 double precision, intent(out) :: cgd double precision, intent(out) :: cgf double precision, intent(out) :: cwd double precision, intent(out) :: cwf double precision :: dg double precision, intent(out) :: ds1 double precision, intent(out) :: ds2 double precision :: lambda double precision, intent(out) :: mugf double precision :: strdamf double precision, intent(in) :: teken !< Flow direction, w.r.t. the structure's orientation. So: based on both upwind *and* flow-link<-->str-pli crossing. double precision :: w2 double precision, intent(out) :: wsd double precision, intent(out) :: wstr double precision, intent(out) :: gatedoorheight double precision :: zb2 double precision :: zs ! ! ! Local variables ! double precision :: help double precision :: w1 double precision :: wsdl double precision :: wsdr double precision :: zb1 double precision :: zbsl double precision :: zbsr ! ! !! executable statements ------------------------------------------------------- ! ! !======================================================================= ! Deltares ! One-Two Dimensional Modelling System ! S O B E K ! ! Subsystem: Flow Module ! ! Programmer: J.Brouwer ! ! Module: FLGTAR (FLow get General sTructure ARguments) ! ! Module description: Parameters for the general structure are extracted ! from the structures module. ! ! ! Parameters: ! NR NAME IO DESCRIPTION ! 13 cgd O Correction coefficient for drowned gate flow. ! 12 cgf O Correction coefficient for free gate flow. ! 15 cwd O Correction coefficient for drowned weir flow. ! 14 cwf O Correction coefficient for free weir flow. ! 9 dg O Gate opening height. ! 10 ds1 O Delta s1 general structure. ! 11 ds2 O Delta s2 general structure. ! 1 istru I Number of structure. ! 17 lambda O Extra resistance ! 16 mugf O Contraction coefficient for free gate flow. ! 3 tekenstr I Flow direction, w.r.t. structure orientation (+/-). ! 6 w2 O Width at right side of structure. ! 7 wsd O Width structure right or left side. ! 5 wstr O Width at centre of structure. ! 8 zb2 O Bed level at right side of structure. ! 4 zs O Bed level at centre of structure. !======================================================================= ! Include Pluvius data space ! ! Declaration of parameters: ! ! ! Declaration of local variables: ! ! ! Fetch parameters from structure info array ! if (generalstruc(ng)%numlinks <= 1) then ! Structure crosses just one link, use user-specified widths w1 = min(wuL, generalstruc(ng)%widthleftW1) wsdl = min(wuL, generalstruc(ng)%widthleftWsdl) ! wstr = generalstruc(ng)%widthcenter wstr = min(wuL, generalstruc(ng)%widthcenteronlink(L0)) ! Possible realtime-controlled, even when crossing one link, so use widthcenteronlink here already. wsdr = min(wuL, generalstruc(ng)%widthrightWsdr) w2 = min(wuL, generalstruc(ng)%widthrightW2) else ! Structure crosses more than one link: nonsensible to use single width left/right etc. same for all links. Use center linkwidth instead (i.e., typically wu(Lf)) ! TODO: UNST-695: Support sideways movement/closing of a gate. w1 = min(wuL, generalstruc(ng)%widthcenteronlink(L0)) !widthleftW1 wsdl = min(wuL, generalstruc(ng)%widthcenteronlink(L0)) !widthleftWsdl ! wstr = generalstruc(ng)%widthcenter wstr = min(wuL, generalstruc(ng)%widthcenteronlink(L0)) wsdr = min(wuL, generalstruc(ng)%widthcenteronlink(L0)) !widthrightWsdr w2 = min(wuL, generalstruc(ng)%widthcenteronlink(L0)) !widthrightW2 end if ! zs = generalstruc(ng)%levelcenter ! comes from ec zb1 = max(bl1, generalstruc(ng)%levelleftZb1) zbsl = max(bl1, generalstruc(ng)%levelleftZbsl) zbsr = max(bl2, generalstruc(ng)%levelrightZbsr) zb2 = max(bl2, generalstruc(ng)%levelrightZb2) ! dg = generalstruc(ng)%gateheight - zs ! also comes from ec lambda = generalstruc(ng)%extraresistance strdamf = generalstruc(ng)%dynstructext gatedoorheight = generalstruc(ng)%gatedoorheight !if (strdamf< - 0.5D0) strdamf = dynstructext !if (lambda < - 0.5D0) lambda = extra_resist_genstruc ! ! Determine cgf, cgd, cwf, cwd, mugf ! (flow direction dependent) ! if (teken>0.0D0) then cgf = generalstruc(ng)%pos_freegateflowcoeff cgd = generalstruc(ng)%pos_drowngateflowcoeff cwf = generalstruc(ng)%pos_freeweirflowcoeff cwd = generalstruc(ng)%pos_drownweirflowcoeff mugf = generalstruc(ng)%pos_contrcoeffreegate else cgf = generalstruc(ng)%neg_freegateflowcoeff cgd = generalstruc(ng)%neg_drowngateflowcoeff cwf = generalstruc(ng)%neg_freeweirflowcoeff cwd = generalstruc(ng)%neg_drownweirflowcoeff mugf = generalstruc(ng)%neg_contrcoeffreegate endif ! ! Determine flow direction dependent parameters ! if (teken>0.0D0) then wsd = wsdr ds1 = zs - zbsr ds2 = zbsr - zb2 else wsd = wsdl ds1 = zs - zbsl ds2 = zbsl - zb1 help = w1 w1 = w2 w2 = help help = zb1 zb1 = zb2 zb2 = help endif end subroutine flgtarfm subroutine togeneral(ng, hulp, ngen, widths) use m_strucs use m_alloc implicit none integer, intent(in) :: ng !< Index of this general structure in the generalstruc(:) array double precision, intent(in) :: hulp(25) !< genstru params read from file integer, intent(in) :: ngen !< Number of flow links crossed by this single general structure double precision, intent(in) :: widths(ngen) !< wu(L) values for all links crossed by this single general structure generalstruc(ng)%widthleftW1 = hulp( 1) !< this and following: see Sobek manual generalstruc(ng)%levelleftZb1 = hulp( 2) generalstruc(ng)%widthleftWsdl = hulp( 3) generalstruc(ng)%levelleftZbsl = hulp( 4) generalstruc(ng)%widthcenter = hulp( 5) generalstruc(ng)%levelcenter = hulp( 6) generalstruc(ng)%widthrightWsdr = hulp( 7) generalstruc(ng)%levelrightZbsr = hulp( 8) generalstruc(ng)%widthrightW2 = hulp( 9) generalstruc(ng)%levelrightZb2 = hulp(10) generalstruc(ng)%gateheight = hulp(11) generalstruc(ng)%gateheightintervalcntrl = hulp(12) generalstruc(ng)%pos_freegateflowcoeff = hulp(13) generalstruc(ng)%pos_drowngateflowcoeff = hulp(14) generalstruc(ng)%pos_freeweirflowcoeff = hulp(15) generalstruc(ng)%pos_drownweirflowcoeff = hulp(16) generalstruc(ng)%pos_contrcoeffreegate = hulp(17) generalstruc(ng)%neg_freegateflowcoeff = hulp(18) generalstruc(ng)%neg_drowngateflowcoeff = hulp(19) generalstruc(ng)%neg_freeweirflowcoeff = hulp(20) generalstruc(ng)%neg_drownweirflowcoeff = hulp(21) generalstruc(ng)%neg_contrcoeffreegate = hulp(22) generalstruc(ng)%extraresistance = hulp(23) ! lambda = L*g/ (C*C) generalstruc(ng)%dynstructext = hulp(24) if (hulp(25) > 0d0) then generalstruc(ng)%gatedoorheight = hulp(25) endif generalstruc(ng)%stabilitycounter = 0d0 ! hulp(25) call realloc(generalstruc(ng)%widthcenteronlink, ngen) generalstruc(ng)%widthcenteronlink(1:ngen) = widths(1:ngen) call realloc(generalstruc(ng)%gateheightonlink, ngen) generalstruc(ng)%gateheightonlink(1:ngen) = generalstruc(ng)%gateheight generalstruc(ng)%numlinks = ngen end subroutine togeneral !> Determines flow link' upwind/downwind parameters based on current velocities and water levels. !! NOTE that this is purely for this flow link, independent of left-right orientation of the structure itself. !! (Motivation: a single structure in 2D may be crossed by multiple flow links, with varying 1->2 orientation.) subroutine flupdofm(m, il, ir, istru, velheight, rholeft, rhoright, crest, & husb, hdsb, uu, ud, teken, relax) use m_strucs use m_flowgeom use m_flow implicit none ! ! Global variables ! integer, intent(in) :: m !< Flow link number, signed! If m < 0 then flow link is in opposite direction than structure left-right orientation. integer, intent(in) :: il,ir,istru logical, intent(in) :: velheight double precision, intent(in) :: crest, relax double precision :: hdsb double precision :: husb double precision :: rholeft double precision :: rhoright double precision :: teken double precision :: ud double precision :: uu double precision :: tem, ucxku, ucyku integer :: L, k, LL,iflip ! Parameters: ! NR NAME IO DESCRIPTION ! 7 crest I crest. ! 10 hd O Downstream water level at t=(n+1). ! 9 husb O Upstream water level at t=(n+1). ! 2 il I Grid point on left side of structure (lower ! index). ! 3 ir I Grid point on right side of structure (upper ! index). ! 4 istru I structure number ! 1 m I Index of velocity point, signed! If m < 0 then flow link is in opposite direction than structure left-right orientation. ! 5 rholeft O Density of diluted water on left grid point. ! 6 rhoright O Density of diluted water on left grid point. ! 13 teken O Flow direction (+1/-1). ! 12 ud O Downstream velocity. ! 11 uu O Upstream velocity. ! 8 velheight I True if velicity height will be taken into account L = abs(m) iflip = max(0, sign(1,m)) ! iflip: 0 if flow link has same orientation as structure's left-right, 1 if opposite (because then the 'left' str point == ln(2,Lf)) if (relax .ne. 1d0) then husb = s1(il)*relax + (1.D0 - relax)*strhis2(9+iflip,istru) ! TODO: HK: strhis2 is not yet filled anywhere (no relaxation possible) hdsb = s1(ir)*relax + (1.D0 - relax)*strhis2(10-iflip,istru) else husb = s1(il) hdsb = s1(ir) endif uu = 0.D0 ud = 0.D0 if (velheight) then uu = 0d0 ; ud = 0d0 do k = 1,nd(il)%lnx LL = iabs( nd(il)%ln(k) ) if (iadv(LL) .ne. 22) then ! any non-structure point uu = max(uu, abs(u1(LL)) ) endif enddo do k = 1,nd(ir)%lnx LL = iabs( nd(ir)%ln(k) ) if (iadv(LL) .ne. 22) then ! any non-structure point ud = max(ud, abs(u1(LL)) ) endif enddo !uu = call getucxucynoweirs(il, ucxku, ucyku, 6) !uu = csu(L)*ucxku + snu(L)*ucyku !call getucxucynoweirs(ir, ucxku, ucyku, 6) !ud = csu(L)*ucxku + snu(L)*ucyku endif if (u1(L) > 0d0) then teken = 1d0 else if (u1(L) < 0d0) then teken = -1d0 else if (s1(iL) > s1(ir) ) then teken = 1d0 elseif (s1(iL) < s1(ir) ) then teken = -1d0 else ! s1(iL) == s1(ir) teken = -dble(sign(1,m)) ! account for orientation of flow link w.r.t. structure endif if (teken < 0) then tem = hdsb ; hdsb = husb ; husb = tem tem = ud ; ud = uu ; uu = tem endif end subroutine flupdofm subroutine flqhgsfm(m, teken, husb, hdsb, uu, zs, wstr, w2, wsd, zb2, ds1, ds2, & dg, cgf, cgd, cwf, cwd, mugf, lambda, strdamf, jarea, ds) use m_flow, only : au, fu, ru use m_physcoef, only : ag implicit none ! ! Global variables ! integer :: m logical, intent(in) :: jarea double precision :: cgd double precision :: cgf double precision :: cwd double precision, intent(in) :: cwf double precision :: dg double precision :: ds double precision :: ds1 double precision :: ds2 double precision :: hdsb double precision :: husb double precision :: lambda double precision :: mugf double precision :: rhoast=1d0 double precision :: strdamf double precision :: teken double precision, intent(in) :: uu double precision :: w2 double precision :: wsd double precision :: wstr double precision :: zb2 double precision :: zs ! ! ! Local variables ! integer :: formno logical :: dpsequfm logical :: imag double precision :: cgd2 double precision :: cgda double precision :: cgfa double precision :: cwfa double precision :: dc double precision :: dlim double precision :: elu double precision :: hd1 double precision :: hs1 double precision :: mugfa double precision :: velhght, tr ! ! !! executable statements ------------------------------------------------------- ! ! ! !======================================================================= ! Deltares ! One-Two Dimensional Modelling System ! S O B E K ! ! Subsystem: Flow Module ! ! Programmer: J.Kuipers ! ! Module: FLQHGS (FLow QH relation for General Structure) ! ! Module description: The QH-relationship for a general structure ! will be transformed to a linearized equation ! ! In this subroutine for given upstream and down- ! stream water levels and upstream velocity ! the flow condition will be determined. ! ! ! Parameters: ! NR NAME IO DESCRIPTION ! 16 cgd I Contraction coefficient for drowned gate flow ! 15 cgf I Contraction coefficient for gate flow ! 18 cwd I Contraction coefficient for drowned weir flow. ! 17 cwf I Contraction coefficient for free weir flow. ! 13 dg I Gate opening height. ! 22 ds O Water level immediately downstream the gate. ! 11 ds1 I Delta s1 general structure. ! 12 ds2 I Delta s2 general structure. ! 4 hdsb I Downstream water level. ! 3 husb I Upstream water level. ! 21 jarea I If True then claculate only area ! 1 m I Grid index of structure ! 20 lambda I Extra resistance ! 19 mugf I Vertical contraction coefficient for free ! gate flow. ! 14 rhoast I Ratio of density right and left of structure ! 2 teken I Flow direction (+1/-1). ! 5 uu I Upstream velocity. ! 8 w2 I Width at right side of structure. ! 9 wsd I Width structure right or left side. ! 7 wstr I Width at centre of structure. ! 10 zb2 I Bed level at right side of structure. ! 6 zs I Bed level at centre of structure. ! ! Subprogram calls: ! NAME DESCRIPTION ! flccgs FLow contraction coefficients for general structure ! flgsd2 FLow general structure depth sill 2nd order equation ! flgsd3 FLow general structure depth sill 3rd order equation ! flgsfuru FLow general structure calculate FU and RU !======================================================================= ! Include Pluvius data space ! Declaration of parameters: ! Declaration of local variables: ! Function declaration: ! ! Compute upstream velocity height and energy level ! velhght = uu*uu/(2.0D0*ag) elu = husb + velhght hs1 = elu - zs ! tr = 1d-4 if (hs1 < tr .or. wstr < tr .or. dg < tr) THEN ! .or. min(cgf, cgd, cwf, cwd)<=0. .or. & ! & dg<.0001) then !hk: or gate closed formno = 0 ; return else ! ! Compute critical water depth at the ! sill, dc and water depth at the sill,ds ! dlim = hs1*(wstr/w2*2./3.*sqrt(2./3.))**(2.0/3.0) hd1 = max(hdsb, zb2 + dlim*0.9D0) ! dc = 2.0D0/3.0D0*hs1 ! ! Calculate ds by solving third order algebraic equation ! call flgsd3fm(wsd, wstr, zs, w2, zb2, ds1, ds2, elu, hd1, rhoast, cwd, ds, & & lambda) ! if (ds>=dc) then ! waterheight on crest larger than critical height on crest if (dg>=ds) then ! ! - drowned weir - ! formno = 2 else ! ! - gate flow - ! formno = 3 ! ! adapt coefficients on basis of Ds & Cwd ! call flccgsfm(dg, ds, cgd, cgf, cwd, mugf, cgda, cgfa, mugfa) endif else ! ! Adapt Cwf coefficient ! if (cwf0.0D0) then cwfa = min(dc/ds*cwd, cwf) else cwfa = cwf endif ! if (dg>=dc) then ! ! - free weir - ! formno = 1 ds = dc else ! ! - gate flow - ! formno = 3 ! ! adapt coefficients on basis of Dc & Cwf ! call flccgsfm(dg, dc, cgd, cgf, cwfa, mugf, cgda, cgfa, mugfa) endif endif ! ! In case of gate flow determine type of gate flow ! (drowned or free) ! if (formno==3) then dc = mugfa*dg ! ! Cgd for second order equation = Cgd' * Mu' ! cgd2 = cgda*mugfa ! call flgsd2fm(wsd, wstr, zs, w2, zb2, dg, ds1, ds2, elu, hd1, rhoast, & & cgd2, imag, ds, lambda) ! if (imag) then ! ! - free gate - ! formno = 3 ds = dc elseif (ds<=dc) then ! ! - free gate - ! formno = 3 ! ! Adapt coefficients ! if (cgda>cgfa) then if (.not.dpsequfm(dc, 0.0D0, 1.0D-20)) then cgfa = max(ds/dc*cgda, cgfa) endif elseif (ds>0.0D0) then cgfa = min(dc/ds*cgda, cgfa) else endif ds = dc !TEM WRITE (11,*) 'cgfa,mugfa',cgfa,mugfa else ! ! - drowned gate - ! formno = 4 endif endif ! ! endif ! !TEM WRITE (11,*) 'formno,ds,dc,dg',formno,ds,dc,dg ! ! The flowe condition is known so calculate ! the linearization coefficients FU and RU ! if (jarea) then call flgsareafm(formno, m, husb, velhght, zs, ds, dg, wstr) else call flgsfurufm(formno, m, teken, husb, hdsb, velhght, zs, ds, dg, dc, wstr, & cwfa, cwd, mugfa, cgfa, cgda, strdamf, lambda) endif end subroutine flqhgsfm subroutine flccgsfm(dg, dsc, cgd, cgf, cw, mugf, cgda, cgfa, mugfa) !!--copyright------------------------------------------------------------------- ! Copyright (c) 2003, Deltares. All rights reserved. !!--disclaimer------------------------------------------------------------------ ! This code is part of the Delft3D software system. Deltares has ! developed c.q. manufactured this code to its best ability and according to the ! state of the art. Nevertheless, there is no express or implied warranty as to ! this software whether tangible or intangible. In particular, there is no ! express or implied warranty as to the fitness for a particular purpose of this ! software, whether tangible or intangible. The intellectual property rights ! related to this software code remain with Deltares at all times. ! For details on the licensing agreement, we refer to the Delft3D software ! license and any modifications to this license, if applicable. These documents ! are available upon request. !!--version information--------------------------------------------------------- ! $Author$ ! $Date$ ! $Revision$ !!--description----------------------------------------------------------------- ! NONE !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- implicit none ! ! Global variables ! double precision, intent(in) :: cgd double precision, intent(out) :: cgda double precision, intent(in) :: cgf double precision, intent(out) :: cgfa double precision, intent(in) :: cw double precision :: dg double precision :: dsc double precision, intent(in) :: mugf double precision, intent(out) :: mugfa ! ! ! Local variables ! logical :: dpsequfm ! ! !! executable statements ------------------------------------------------------- ! ! !======================================================================= ! Deltares ! One-Two Dimensional Modelling System ! S O B E K ! ! Subsystem: Flow Module ! ! Programmer: J.Brouwer/J.Kuipers ! ! Module: FLCCGS (FLow Corr. Coefficients for General Structure) ! ! Module description: Correct coefficients for gate flow ! ! In the formulas for the gate and weir several ! coefficients are applied. To avoid discontinuities ! in the transition from weir to gate flow, the ! correction coefficient cgd should be corrected. ! ! ! Parameters: ! NR NAME IO DESCRIPTION ! 3 cgd I Correction coefficient for drowned gate flow. ! 7 cgda O Adapted correction coefficient for drowned ! gate flow. ! 4 cgf I Correction coefficient for free gate flow. ! 8 cgfa O Adapted correction coefficient for free gate ! flow. ! 5 cw I Correction coefficient for weir flow. ! 1 dg I Gate opening height. ! 2 dsc I Depth at sill or critical depth. ! 6 mugf I Contraction coefficient for free gate flow. ! 9 mugfa O Adapted contraction coefficient for free gate ! flow. !======================================================================= ! ! Declaration of parameters: ! ! ! Logical function ! ! ! dsc contains ds or dc ! if (.not.dpsequfm(dsc, 0.0D0, 1.D-20)) then ! if (dg/dsc>mugf) then mugfa = dg/dsc else mugfa = mugf endif ! if (cgd>cw) then if (dpsequfm(dg, 0.0D0, 1.0D-20)) then cgda = cgd else cgda = min(dsc/dg*cw, cgd) endif else cgda = max(dg/dsc*cw, cgd) endif ! if (cgf>cw) then if (dpsequfm(dg, 0.0D0, 1.0D-20)) then cgfa = cgf else cgfa = min(dsc/dg*cw, cgf) endif else cgfa = max(dg/dsc*cw, cgf) endif ! else mugfa = mugf cgda = cgd cgfa = cgf endif end subroutine flccgsfm function dpsequfm(dvar1, dvar2, eps) ! equal within eps? implicit none logical :: dpsequfm double precision, intent(in) :: dvar1, dvar2, eps dpsequfm = abs(dvar1 - dvar2) 0) then if (allocated (generalstruc) ) deallocate (generalstruc) allocate (generalstruc(mxgeneral), stat=ierr) endif if (mxuniversal > 0) then if (allocated (universalstruc) ) deallocate (universalstruc) allocate (universalstruc(mxuniversal), stat=ierr) endif end subroutine readandallocstructures subroutine furusobekstructures() use m_flow use m_flowgeom use m_strucs implicit none integer :: istru, ng, n, numgen, L, Ls double precision :: zup, bup, a, fac logical :: firstiter=.true. , jarea= .false. firstiter = .true. jarea = .false. do ng = 1, ncgensg ! loop over generalstruc signals, sethu do n = L1cgensg(ng), L2cgensg(ng) L = kcgen(3,n) if (kcgen(1,n) == ln(2,L)) then Ls = -L ! Flow link has opposite orientation to structure's orientation. else Ls = L end if if (hu(L) > 0d0) then ! hu is above lowest sill call flgsfm( n, ng, Ls, firstiter , jarea ) endif enddo enddo do ng = 1, ngatesg zup = zgate(ng) do n = L1gatesg(ng), L2gatesg(ng) L = kgate(3,n) if (hu(L) > 0d0) then bup = min(bob(1,L), bob(2,L) ) a = zup - bup fac = min (1d0, max(1d-2, a*huvli(L) ) ) fu(L) = fac*fu(L) endif enddo enddo end subroutine furusobekstructures !> Computes and sets the widths and gate lower edge levels on each of the flow links !! crossed by a general structure (gate/weir/true genstru). !! This is now an extended version of SOBEK's setLineStructure, because it also enables !! a sideways closing gate with two doors from the left and right side, where the partially !! closed portions have gate flow, and the center open portion still only has normal weir !! flow across the sill. subroutine update_zcgen_widths_and_heights() use m_strucs use m_flowexternalforcings use m_flowgeom use m_structures implicit none double precision :: crestwidth, totalWidth, closedWidth, closedGateWidthL, closedGateWidthR, help integer :: ng, L, L0, Lf do ng=1,ncgensg ! Loop over general structures ! 1: First determine total width of all genstru links (TODO: AvD: we should not recompute this every user time step) totalWidth = 0d0 if (generalstruc(ng)%numlinks == 0) then cycle ! Only upon invalid input (see warnings in log about missing structure params) end if do L=L1cgensg(ng),L2cgensg(ng) L0 = L-L1cgensg(ng)+1 Lf = kcgen(3,L) generalstruc(ng)%widthcenteronlink(L0) = wu(Lf) totalWidth = totalWidth + wu(Lf) end do ! 2a: the desired crest width for this overall structure (hereafter, the open links for this genstru should add up to this width) ! Also: only for gates, the desired door opening width for this overall structure ! (should be smaller than crestwidth, and for this portion the open gate door is emulated by dummy very high lower edge level) if (cgen_type(ng) == ICGENTP_WEIR) then crestwidth = zcgen((ng-1)*3+3) closedGateWidthL = 0d0 closedGateWidthR = 0d0 else if (cgen_type(ng) == ICGENTP_GENSTRU) then crestwidth = totalWidth ! No crest/sill-width setting for true general structure yet (not old ext, nor new ext) closedGateWidthL = max(0d0, .5d0*(totalWidth - zcgen((ng-1)*3+3))) ! Default symmetric opening closedGateWidthR = max(0d0, .5d0*(totalWidth - zcgen((ng-1)*3+3))) else if (cgen_type(ng) == ICGENTP_GATE) then ! For a gate: zcgen(3,ng) is limited to the door opening width, but we want to open all links ! *underneath* the two doors as well, (if lower_edge_level is still high enough above sill_level) crestwidth = gates(cgen2str(ng))%sill_width if (gates(cgen2str(ng))%opening_direction == IOPENDIR_FROMLEFT) then closedGateWidthL = max(0d0, totalWidth - zcgen((ng-1)*3+3)) closedGateWidthR = 0d0 else if (gates(cgen2str(ng))%opening_direction == IOPENDIR_FROMRIGHT) then closedGateWidthL = 0d0 closedGateWidthR = max(0d0, totalWidth - zcgen((ng-1)*3+3)) else ! IOPENDIR_SYMMETRIC closedGateWidthL = max(0d0, .5d0*(totalWidth - zcgen((ng-1)*3+3))) closedGateWidthR = max(0d0, .5d0*(totalWidth - zcgen((ng-1)*3+3))) end if generalstruc(ng)%gateheightonlink(1:generalstruc(ng)%numlinks) = huge(1d0) ! As a start, gate door is open everywhere. Below, we will close part of the gate doors. end if ! 2a: Determine the width that needs to be closed on 'left' side ! close the line structure from the outside to the inside: first step increasing increments closedWidth = max(0d0, totalWidth - crestwidth)/2d0 do L=L1cgensg(ng),L2cgensg(ng) L0 = L-L1cgensg(ng)+1 Lf = kcgen(3,L) if (closedWidth > 0d0) then help = min (wu(Lf), closedWidth) generalstruc(ng)%widthcenteronlink(L0) = wu(Lf) - help ! 0d0 if closed closedWidth = closedWidth - help end if if (cgen_type(ng) == ICGENTP_GATE .and. closedGateWidthL > 0d0 .or. cgen_type(ng) == ICGENTP_GENSTRU ) then !if (closedGateWidthL > .5d0*wu(Lf)) then generalstruc(ng)%gateheightonlink(L0) = zcgen((ng-1)*3+2) help = min (wu(Lf), closedGateWidthL) closedGateWidthL = closedGateWidthL - help !end if else end if if (closedWidth <= 0d0 .and. closedGateWidthL <= 0d0) then ! finished exit endif enddo ! 2b: Determine the width that needs to be closed on 'right' side ! second step decreasing increments closedWidth = max(0d0, totalWidth - crestwidth)/2d0 do L=L2cgensg(ng),L1cgensg(ng),-1 L0 = L-L1cgensg(ng)+1 Lf = kcgen(3,L) if (closedWidth > 0d0) then help = min (wu(Lf), closedWidth) ! generalstruc(ng)%widthcenteronlink(L0) = wu(Lf) - help ! 0d0 if closed generalstruc(ng)%widthcenteronlink(L0) = generalstruc(ng)%widthcenteronlink(L0) - help ! 0d0 if closed closedWidth = closedWidth - help end if if (cgen_type(ng) == ICGENTP_GATE .and. closedGateWidthR > 0d0 .or. cgen_type(ng) == ICGENTP_GENSTRU) then !if (closedGateWidthL > .5d0*wu(Lf)) then generalstruc(ng)%gateheightonlink(L0) = zcgen((ng-1)*3+2) help = min (wu(Lf), closedGateWidthR) closedGateWidthR = closedGateWidthR - help !end if end if if (closedWidth <= 0d0 .and. closedGateWidthR <= 0d0) then ! finished exit endif enddo end do ! 1,ngensg end subroutine update_zcgen_widths_and_heights subroutine enloss(ag ,d1 ,eweir ,hkruin ,hov , & & qunit ,qvolk ,toest ,vov , & & ewbov ,ewben ,wsben ,dte ,dtefri ) !----- GPL --------------------------------------------------------------------- ! ! Copyright (C) Stichting Deltares, 2011-2014. ! ! This program is free software: you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by ! the Free Software Foundation version 3. ! ! This program 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 General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with this program. If not, see . ! ! contact: delft3d.support@deltares.nl ! Stichting Deltares ! P.O. Box 177 ! 2600 MH Delft, The Netherlands ! ! All indications and logos of, and references to, "Delft3D" and "Deltares" ! are registered trademarks of Stichting Deltares, and remain the property of ! Stichting Deltares. All rights reserved. ! !------------------------------------------------------------------------------- ! Original URL: https://svn.oss.deltares.nl/repos/delft3d/trunk/src/engines_gpl/flow2d3d/packages/kernel/src/compute/enloss.f90 !!--description----------------------------------------------------------------- ! ! Function: Determines additional energy loss due to weir. ! Energy loss dependent on velocity perpendicular ! to weir (Schonfeld, 1955). ! Subroutine based on subroutine ENLOSS in WAQUA. ! Energyloss based on Carnot relation for ! decellerated flow VOV < 0.25 m/s ! or on "Tables" VOV > 0.5 m/s or ! average for 0.25 < VOV < 0.50 m/s. ! In 2003 the improved formulations described in ! "Beoordeling nieuwe overlaatroutines WAQUA" ! WL-report Z3063, have been implemented and tested. ! !!--pseudo code and references-------------------------------------------------- ! ! "Weergave van extra energieverlies in RIVCUR." ! (J.H.A. Wybenga, report Deltares Q779 voor RWS/RIZA, 1989). ! "Energieverliezen door overlaten: een gewijzigde berekeningsprocedure voor ! WAQUA-rivieren versie." ! (H. Vermaas, verslag onderzoek Deltares Q92, 1987) ! "Beoordeling nieuwe overlaatroutines WAQUA" ! (J. van Kester, verslag Z3063, juli 2001) ! !!--declarations---------------------------------------------------------------- use precision implicit none ! ! Global variables ! real(fp) , intent(in) :: ag ! Description and declaration in esm_alloc_real.f90 real(fp) :: d1 !! Distance between crest and downstream depth real(fp) , intent(out):: dte !! Subgrid energy loss due to weir real(fp) , intent(in) :: dtefri !! Energy loss due to friction real(fp) , intent(in) :: ewben !! Excess waterdepth downstream real(fp) , intent(in) :: ewbov !! Excess waterdepth upstream real(fp) :: eweir !! Energy level at weir real(fp) , intent(in) :: hkruin !! Crest height (downward positive). real(fp) :: hov !! Total water depth at crest weir real(fp) :: qunit !! Discharge at weir crest real(fp) :: qvolk !! Maximum discharge (super critical flow) real(fp) , intent(in) :: vov !! Velocity at crest of weir real(fp) , intent(in) :: wsben !! Downstream water level character(4) :: toest !! State weir: !! volk = perfect weir !! onvo = imperfect weir ! ! Local variables ! real(fp) :: dtecar real(fp) :: dteonv real(fp) :: dtetab real(fp) :: dtevol real(fp) :: qqv double precision :: tabellenboek real(fp) :: theta real(fp) :: vben real(fp) :: wsov ! !! executable statements ------------------------------------------------------- ! ! ! Determine energy loss ! qqv = qunit/qvolk ! ! Imperfect weir (sub critical flow) ! if ((wsben + hkruin + d1)==0.0) then ! ! Dry bed downstream, could perhaps also check on qunit==0.0 ! dteonv = 0.0 elseif (abs(vov)<=0.25) then ! ! Energy loss according Carnot law ! WSBEN+HKRUIN+D1 := H0 (S0+DPS) ! dteonv = (vov - qunit/(wsben + hkruin + d1))**2/(2.*ag) elseif (abs(vov)>0.25 .and. abs(vov)<0.5) then ! ! Weigthing between Carnot and Tables ! WSBEN+HKRUIN+D1 := H0 (S0+DPS) ! dtecar = (vov - qunit/(wsben + hkruin + d1))**2/(2.*ag) dtetab = tabellenboek(d1, eweir, qunit, qvolk) theta = (abs(vov) - 0.25)/0.25 dteonv = (1 - theta)*dtecar + theta*dtetab elseif (abs(vov)>=0.5) then ! ! Energy loss according to Tables from Vermaas ! dteonv = tabellenboek(d1, eweir, qunit, qvolk) else endif ! ! Determine energy loss for free weir flow ! dtevol = ewbov - ewben - dtefri ! ! Determine if it is free weir flow or submerged weir flow ! if (dtevol*qqv**2>=dteonv) then ! ! It is a free weir flow ! toest = 'volk' else ! ! It is a submerged weir flow ! toest = 'onvo' endif ! ! Energy loss ! if (toest=='volk') then ! ! It is a free weir flow ! dte = dtevol elseif (toest=='onvo') then ! ! It is a submerged weir flow ! dte = dteonv else endif end subroutine enloss double precision function tabellenboek(d1 ,eweir ,qunit ,qvolk ) !----- GPL --------------------------------------------------------------------- ! ! Copyright (C) Stichting Deltares, 2011-2014. ! ! This program is free software: you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by ! the Free Software Foundation version 3. ! ! This program 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 General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with this program. If not, see . ! ! contact: delft3d.support@deltares.nl ! Stichting Deltares ! P.O. Box 177 ! 2600 MH Delft, The Netherlands ! ! All indications and logos of, and references to, "Delft3D" and "Deltares" ! are registered trademarks of Stichting Deltares, and remain the property of ! Stichting Deltares. All rights reserved. ! !------------------------------------------------------------------------------- ! $Id: manholes.f90 42642 2015-10-21 11:34:20Z dam_ar $ ! Original URL: https://svn.oss.deltares.nl/repos/delft3d/trunk/src/engines_gpl/flow2d3d/packages/kernel/src/compute/tabel.f90 !!--description----------------------------------------------------------------- ! ! Function: Determines additional energy loss due to weir. ! Energy loss dependent on velocity perpendicular ! to weir (Schonfeld, 1955). ! Subroutine based on subroutine TABEL in WAQUA. ! Energyloss based on Carnot relation for ! decellerated flow, or on "Tabellenboek". ! This subroutine is also called for supercritical ! flow (qqv = 1.0), see the formulations in ! "Beoordeling nieuwe overlaatroutines WAQUA" ! WL-report Z3063. ! Method used: Reference : Weergave van extra energieverlies ! in RIVCUR. (J.H.A. Wybenga, report ! Deltares Q779 voor RWS/RIZA, 1989). ! Energieverliezen door overlaten: een gewijzigde ! berekeningsprocedure voor WAQUA-rivieren versie. ! (H. Vermaas, verslag onderzoek Deltares ! Q92, 1987) ! !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- use precision implicit none ! ! Global variables ! real(fp) , intent(in) :: d1 !! Distance between crest and downstream depth real(fp) , intent(in) :: eweir !! Energy level at weir real(fp) , intent(in) :: qunit !! Discharge at weir crest real(fp) , intent(in) :: qvolk !! Maximum discharge (super critical flow) ! ! ! Local variables ! real(fp) :: f1 real(fp) :: f1low real(fp) :: f1up real(fp) :: f2 real(fp) :: f2low real(fp) :: f2up real(fp) :: qqv real(fp) :: qqvlow real(fp) :: qqvup real(fp) :: theta ! ! !! executable statements ------------------------------------------------------- ! ! qqv = min(qunit/qvolk, 1.0_fp) ! !-----Calculate energy loss according to tables (based on experinments ! in the delta flums) ! if (qqv<0.3) then qqvlow = 0. qqvup = 0.3 if (eweir<0.5) then f1low = 0.0 f1up = 0.0092*eweir f2low = 0.0 f2up = 0.328 else f1low = 0.0 f1up = 0.0057 + 0.002*log10(eweir - 0.21) f2low = 0.0 f2up = 0.314 + 0.0237*eweir**0.808 endif ! elseif (qqv<0.6) then qqvlow = 0.3 qqvup = 0.6 if (eweir<0.5) then f1low = 0.0092*eweir f1up = 0.0211*eweir f2low = 0.328 f2up = 0.316 else f1low = 0.0057 + 0.002*log10(eweir - 0.21) f1up = -0.0017 + 0.033*log10(eweir + 1.85) f2low = 0.314 + 0.0237*eweir**0.808 f2up = 0.283 + 0.0531*eweir**0.702 endif ! elseif (qqv<0.8) then qqvlow = 0.6 qqvup = 0.8 if (eweir<0.5) then f1low = 0.0211*eweir f1up = 0.0463*eweir f2low = 0.316 f2up = 0.295 else f1low = -0.0017 + 0.033*log10(eweir + 1.85) f1up = -0.0022 + 0.064*log10(eweir + 1.99) f2low = 0.283 + 0.0531*eweir**0.702 f2up = 0.280 + 0.0272*eweir**0.912 endif ! elseif (qqv<0.9) then qqvlow = 0.8 qqvup = 0.9 if (eweir<0.5) then f1low = 0.0463*eweir f1up = 0.0657*eweir f2low = 0.295 f2up = 0.238 else f1low = -0.0022 + 0.064*log10(eweir + 1.99) f1up = 0.039 + 0.076*log10(eweir + 0.33) f2low = 0.280 + 0.0272*eweir**0.912 f2up = 0.224 + 0.0256*eweir**0.869 endif ! elseif (qqv<0.95) then qqvlow = 0.9 qqvup = 0.95 if (eweir<0.5) then f1low = 0.0657*eweir f1up = 0.0772*eweir f2low = 0.238 f2up = 0.206 else f1low = 0.039 + 0.076*log10(eweir + 0.33) f1up = 0.065 + 0.116*log10(eweir + 0.092) f2low = 0.224 + 0.0256*eweir**0.869 f2up = 0.124 + 0.11*eweir**0.422 endif ! elseif (qqv<0.99) then qqvlow = 0.95 qqvup = 0.99 if (eweir<0.5) then f1low = 0.0772*eweir f2low = 0.206 else f1low = 0.065 + 0.116*log10(eweir + 0.092) f2low = 0.124 + 0.11*eweir**0.422 endif if (eweir<1.0) then f1up = 0.115*eweir**0.9 f2up = 0.242 else f1up = 0.133 + 0.213*log10(eweir - 0.173) f2up = 0.133 + 0.109*eweir**0.619 endif ! elseif (qqv>=0.99) then qqvlow = 0.99 qqvup = 1.0 if (eweir<1.0) then f1low = 0.115*eweir**0.9 f2low = 0.242 else f1low = 0.133 + 0.213*log10(eweir - 0.173) f2low = 0.133 + 0.109*eweir**0.619 endif if (eweir<2.0) then f1up = 0.156*eweir**0.75 f2up = 0.343 else f1up = 0.244 + 0.172*log10(eweir - 0.718) f2up = 0.075 + 0.176*eweir**0.609 endif ! else endif ! !-----Calculate terms for energy loss ! theta = (qqv - qqvlow)/(qqvup - qqvlow) f1 = (1 - theta)*f1low + theta*f1up f2 = (1 - theta)*f2low + theta*f2up tabellenboek = f1*d1**f2 end function tabellenboek ! subroutine usrbrl2d(icx ,icy ,nmmax ,kmax ,kfu , & ! & kspu ,guu ,gvu ,qxk ,bbk , & ! & ddk ,ubrlsu ,dps ,hkru ,s0 , & ! & hu ,umean ,thick ,dteu ,taubpu , & ! & gdp ) !----- GPL --------------------------------------------------------------------- ! ! ! ! real(fp) , pointer :: eps ! real(fp) , pointer :: ag ! real(fp) , pointer :: thetaw ! ! Global variables ! ! integer , intent(in) :: icx !! Increment in the X-dir., if ICX= NMAX !! then computation proceeds in the X- !! dir. If icx=1 then computation pro- !! ceeds in the Y-dir. ! integer , intent(in) :: icy !! Increment in the Y-dir. (see ICX) ! integer , intent(in) :: kmax ! Description and declaration in esm_alloc_int.f90 ! integer , intent(in) :: nmmax ! Description and declaration in dimens.igs ! integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfu ! Description and declaration in esm_alloc_int.f90 ! integer , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax), intent(in) :: kspu ! Description and declaration in esm_alloc_int.f90 ! real(prec), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: dps ! Description and declaration in esm_alloc_real.f90 ! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: dteu ! Description and declaration in esm_alloc_real.f90 ! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: guu ! Description and declaration in esm_alloc_real.f90 ! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: gvu ! Description and declaration in esm_alloc_real.f90 ! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: hkru ! Description and declaration in esm_alloc_real.f90 ! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: hu ! Description and declaration in esm_alloc_real.f90 ! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: s0 ! Description and declaration in esm_alloc_real.f90 ! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: umean ! Description and declaration in esm_alloc_real.f90 ! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: bbk !! Internal work array, coefficient la- !! yer velocity in (N,M,K) implicit part ! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: ddk !! Internal work array, diagonal space !! at (N,M,K) ! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: qxk ! Description and declaration in esm_alloc_real.f90 ! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: taubpu ! Description and declaration in esm_alloc_real.f90 ! real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: ubrlsu ! Description and declaration in esm_alloc_real.f90 ! real(fp) , dimension(kmax) , intent(in) :: thick ! Description and declaration in esm_alloc_real.f90 ! ! Local variables ! ! integer :: itel ! Loop counter ! integer :: k ! Loop counter over KMAX ! integer :: ksp ! Value for structure point KFU(NM)*ABS (KSPU(NM,0))*KSPU(NM,K) ! integer :: ndm ! NM-ICY ! integer :: ndmu ! NMU-ICY ! integer :: nm ! Loop counter over NMMAX ! integer :: nmu ! NM+ICX ! real(fp) :: avolk ! help constant in computation of qvolk ! real(fp) :: d1 ! Distance between crest and downstream depth ! real(fp) :: dd ! Downstream depth ! real(fp) :: dte0 ! Local backup of DTEU value ! real(fp) :: dtefri ! Energy loss due to friction ! real(fp) :: ewben ! Excess waterdepth downstream ! real(fp) :: ewbov ! Excess waterdepth upstream ! real(fp) :: eweir ! Energy level at weir ! real(fp) :: hkruin ! Crest height (downward positive). ! real(fp) :: hov ! Total water depth at crest weir ! real(fp) :: hvolk ! Water depth associated with critical flow ! real(fp) :: qov ! real(fp) :: qunit ! Discharge at weir crest ! real(fp) :: qvolk ! Maximal discharge (super critical flow ) ! real(fp) :: vbov ! real(fp) :: vov ! real(fp) :: wsben ! Water level d/s ! real(fp) :: wsbov ! Water level u/s ! character(4) :: toest ! Status weir: volk = perfect weir onvo = imperfect weir ! !! executable statements ------------------------------------------------------- ! ! eps => gdp%gdconst%eps ! ag => gdp%gdphysco%ag ! thetaw => gdp%gdrivpro%thetaw ! ! avolk = 2.0*sqrt(2.0*ag/3.0)/3.0 ! ! ndm = -icy ! nmu = icx ! ndmu = -icy + icx ! do nm = 1, nmmax ! ndm = ndm + 1 ! nmu = nmu + 1 ! ndmu = ndmu + 1 ! ! 2D weir, quadratic friction ! ! ksp = kfu(nm)*abs(kspu(nm, 0)) ! hkruin = hkru(nm) ! if ((ksp==9) .and. (real(dps(nm) ,fp)>hkruin) & ! & .and. (real(dps(nmu),fp)>hkruin) ) then ! ! Initial values ! ! qunit = umean(nm)*hu(nm) ! vbov = umean(nm) ! ! Define flow direction and downstream resp. upstream points ! ! if (vbov>=0.001) then ! dd = real(dps(nmu),fp) ! wsbov = s0(nm) ! wsben = s0(nmu) ! elseif (vbov<= - 0.001) then ! dd = real(dps(nm),fp) ! wsbov = s0(nmu) ! wsben = s0(nm) ! elseif (s0(nmu)>s0(nm)) then ! dd = real(dps(nm),fp) ! wsbov = s0(nmu) ! wsben = s0(nm) ! else ! dd = real(dps(nmu),fp) ! wsbov = s0(nm) ! wsben = s0(nmu) ! endif ! ! Initial values ! ! qunit = abs(qunit) ! eweir = (max(1e-6_fp, wsbov + hkruin)) + vbov**2/(2.*ag) ! ewbov = wsbov + hkruin ! ewben = wsben + hkruin ! ! ! qvolk = avolk*eweir**1.5 ! ! ! ! calculate the maximum velocity above the weir crest ! ! reduced wet surface ! ! ! hov = max(1e-6_fp, wsbov + hkruin) ! vov = qunit/(max(1e-6_fp, wsbov + hkruin)) ! itel = 0 ! if (vov<0.5) then ! ! ! ! this is only interesting for velocity less than 0.5 m/s ! ! vov will only be used in cases vov is less than 0.5 m/s ! ! hvolk = (2./3.)*eweir ! ! ! 100 continue ! itel = itel + 1 ! vov = qunit/hov ! hov = max(hvolk, eweir - vov**2/(2.*ag)) ! qov = vov*hov ! ! ! if ((abs(qunit - qov))>0.001*max(0.0001_fp, qunit) .and. (itel<101)) & ! & goto 100 ! endif ! ! ! dte0 = dteu(nm) ! d1 = dd - hkruin ! dtefri = taubpu(nm)*abs(umean(nm))*gvu(nm)/(ag*hu(nm)) ! ! ! ! Determine subgrid energy loss (DTEU) ! ! call enloss(ag ,d1 ,eweir ,hkruin ,hov , & ! & qunit ,qvolk ,toest ,vov , & ! & ewbov ,ewben ,wsben ,dteu(nm) ,dtefri ) ! ! ! ! Use overrelaxation (following experience Stelling '98) ! ! Notice the unusual thetaw definition (thetaw=0: dteu, thetaw=1: dte0) ! ! ! dteu(nm) = (1.0_fp - thetaw)*dteu(nm) + thetaw*dte0 ! if (toest=='volk') vbov = qvolk/max(hu(nm), 1E-6_fp) ! if (vbov>eps) then ! do k = 1, kmax ! bbk(nm, k) = bbk(nm, k) + (ag*dteu(nm)/vbov)*ubrlsu(nm, k) & ! & /(gvu(nm)*thick(k)) ! enddo ! endif ! endif ! enddo ! end subroutine usrbrl2d