subroutine culver(icx ,icy ,kmax ,nsrc ,kfs , & & kfsmax ,kfsmin ,mnksrc ,disch ,dps , & & s0 ,zk ,thick ,voldis ,timsec , & & sumrho ,gdp ) !----- 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$ ! $HeadURL$ !!--description----------------------------------------------------------------- ! Computes the discharge relation through a culvert. ! !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- use precision ! use globaldata ! implicit none ! type(globdat),target :: gdp ! ! The following list of pointer parameters is used to point inside the gdp structure ! real(fp) , dimension(:) , pointer :: arcul real(fp) , dimension(:) , pointer :: calfa real(fp) , dimension(:) , pointer :: clcul real(fp) , dimension(:) , pointer :: cleng real(fp) , dimension(:,:) , pointer :: closs1 real(fp) , dimension(:,:) , pointer :: closs2 real(fp) , dimension(:,:) , pointer :: closs3 real(fp) , dimension(:) , pointer :: cmann real(fp) , dimension(:) , pointer :: htcul integer , dimension(:) , pointer :: numrel1 integer , dimension(:) , pointer :: numrel2 integer , dimension(:) , pointer :: numrel3 real(fp) , dimension(:) , pointer :: poscul real(fp) , dimension(:,:) , pointer :: wetar1 real(fp) , dimension(:,:) , pointer :: wetar2 real(fp) , dimension(:,:) , pointer :: wetar3 real(fp) , dimension(:) , pointer :: wtcul integer , pointer :: lundia real(fp) , pointer :: hdt real(fp) , pointer :: ag real(fp) , pointer :: rhow logical , pointer :: culvert logical , pointer :: zmodel ! ! 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) :: nsrc ! Description and declaration in dimens.igs integer, dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfs ! Description and declaration in esm_alloc_int.f90 integer, dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfsmax ! Description and declaration in esm_alloc_int.f90 integer, dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfsmin ! Description and declaration in esm_alloc_int.f90 integer , intent(in) :: kmax integer, dimension(7, nsrc) :: mnksrc ! Description and declaration in r-i-ch.igs real(fp) , intent(in) :: timsec ! Time in seconds since reference date real(fp), dimension(nsrc) , intent(out) :: disch ! Description and declaration in esm_alloc_real.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) , intent(in) :: s0 ! Description and declaration in esm_alloc_real.f90 real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax), intent(in) :: sumrho ! Description and declaration in esm_alloc_real.f90 real(fp), dimension(kmax) , intent(in) :: thick ! Description and declaration in esm_alloc_real.f90 real(fp), dimension(nsrc) :: voldis ! Description and declaration in esm_alloc_real.f90 real(fp), dimension(0:kmax) , intent(in) :: zk ! See description and declaration of sig in esm_alloc_real.f90 ! ! Local variables ! integer :: ddb integer :: i integer :: icxy integer :: iflow integer, external :: intlay integer :: isrc integer :: iswitch integer :: k integer :: nmin integer :: nmout integer :: kin integer :: kout integer :: swapval real(fp) :: area real(fp) :: cc1 real(fp) :: cd1 real(fp) :: cd2 real(fp) :: cd3 real(fp) :: coefl real(fp) :: delpres ! Water pressure gradient real(fp) :: delzet ! Water level gradient real(fp) :: h0 real(fp) :: height real(fp) :: hin real(fp) :: hout real(fp) :: pos real(fp) :: rdis real(fp) :: rhorefin real(fp) :: rhorefout real(fp) :: zz1 real(fp) :: zz2 real(fp) :: zdown real(fp) :: zup ! Interface to dll is in High precision! ! real(hp) :: disch_dll real(hp) :: pos1_dll real(hp) :: pos2_dll real(hp) :: rmissval integer(pntrsize) :: error_ptr integer(pntrsize), external :: perf_function_culvert character(256) :: errmsg character(256) :: message ! Contains message from ! integer , pointer :: max_integers integer , pointer :: max_reals integer , pointer :: max_strings character(256) , dimension(:), pointer :: dll_function character(256) , dimension(:), pointer :: dll_usrfil integer(pntrsize), dimension(:), pointer :: dll_handle integer , dimension(:), pointer :: dll_integers real(hp) , dimension(:), pointer :: dll_reals character(256) , dimension(:), pointer :: dll_strings ! !! executable statements ------------------------------------------------------- ! arcul => gdp%gdculver%arcul calfa => gdp%gdculver%calfa clcul => gdp%gdculver%clcul cleng => gdp%gdculver%cleng closs1 => gdp%gdculver%closs1 closs2 => gdp%gdculver%closs2 closs3 => gdp%gdculver%closs3 cmann => gdp%gdculver%cmann htcul => gdp%gdculver%htcul numrel1 => gdp%gdculver%numrel1 numrel2 => gdp%gdculver%numrel2 numrel3 => gdp%gdculver%numrel3 poscul => gdp%gdculver%poscul wetar1 => gdp%gdculver%wetar1 wetar2 => gdp%gdculver%wetar2 wetar3 => gdp%gdculver%wetar3 wtcul => gdp%gdculver%wtcul lundia => gdp%gdinout%lundia hdt => gdp%gdnumeco%hdt ag => gdp%gdphysco%ag rhow => gdp%gdphysco%rhow culvert => gdp%gdprocs%culvert zmodel => gdp%gdprocs%zmodel ! max_integers => gdp%gdculver%max_integers max_reals => gdp%gdculver%max_reals max_strings => gdp%gdculver%max_strings dll_function => gdp%gdculver%dll_function dll_usrfil => gdp%gdculver%dll_usrfil dll_handle => gdp%gdculver%dll_handle dll_integers => gdp%gdculver%dll_integers dll_reals => gdp%gdculver%dll_reals dll_strings => gdp%gdculver%dll_strings ! ! Statics ! ddb = gdp%d%ddbound icxy = max(icx, icy) rmissval = -1.0e20_hp ! ! Discharge relation through culvert ! MNKSRC(7,.) = 3: one-way culvert "C" ! MNKSRC(7,.) = 4: one-way culvert especially for "WL Borgerhout" "E" ! MNKSRC(7,.) = 5: two-way culvert especially for "WL Borgerhout" "D" ! MNKSRC(7,.) = 6: power station; not handled here "Q" ! MNKSRC(7,.) = 7: user defined culvert in a dll "U" ! MNKSRC(7,.) = 8: two-way culvert "F" ! do isrc = 1, nsrc ! ! Calculate the loss coefficient ! Only for type 4 and 5 ! if ( (mnksrc(7,isrc) == 4) .or. (mnksrc(7,isrc) == 5) ) then call n_and_m_to_nm(mnksrc(2,isrc), mnksrc(1,isrc), nmin, gdp) call n_and_m_to_nm(mnksrc(5,isrc), mnksrc(4,isrc), nmout, gdp) ! ! for intake point: ! mnksrc(3,isrc) = intlay(lundia ,zmodel ,kfsmin(nmin) ,kfsmax(nmin) , & & poscul(isrc),zk ,dps(nmin) ,s0(nmin) , & & kmax ,thick ,isrc ,'intake ' ) ! ! for outfall point: ! mnksrc(6,isrc) = intlay(lundia ,zmodel ,kfsmin(nmout),kfsmax(nmout), & & poscul(isrc),zk ,dps(nmout) ,s0(nmout) , & & kmax ,thick ,isrc ,'outfall' ) ! ! compute loss coefficient: ! area = 0.0 coefl = 0.0 if (kfs(nmin) == 1 .or. kfs(nmout) == 1) then hin = max (0.0_fp,s0(nmin )-poscul(isrc)) hout = max (0.0_fp,s0(nmout)-poscul(isrc)) height = 0.5 * ( hin + hout) height = min ( height , htcul(isrc) ) area = height * wtcul(isrc) ! ! compute cd1: ! if ( area < wetar1(isrc,1) ) cd1 = closs1(isrc,1) if ( area >= wetar1(isrc,numrel1(isrc)) ) & & cd1 = closs1(isrc,numrel1(isrc)) do i =2,numrel1(isrc) if ( area > wetar1(isrc,i-1) .and. & & area <= wetar1(isrc,i) ) then zz1 = area - wetar1(isrc,i-1) zz2 = wetar1(isrc,i) - wetar1(isrc,i-1) cc1 = closs1(isrc,i) - closs1(isrc,i-1) cd1 = closs1(isrc,i-1) + zz1 / zz2 * cc1 endif enddo ! ! compute cd2: ! if ( area <= wetar2(isrc,1) ) cd2 = closs2(isrc,1) if ( area >= wetar2(isrc,numrel2(isrc)) ) & & cd2 = closs2(isrc,numrel2(isrc)) do i =2,numrel2(isrc) if ( area > wetar2(isrc,i-1) .and. & & area <= wetar2(isrc,i) ) then zz1 = area - wetar2(isrc,i-1) zz2 = wetar2(isrc,i) - wetar2(isrc,i-1) cc1 = closs2(isrc,i) - closs2(isrc,i-1) cd2 = closs2(isrc,i-1) + zz1 / zz2 * cc1 endif enddo ! ! compute cd3: ! if ( area <= wetar3(isrc,1) ) cd3 = closs3(isrc,1) if ( area >= wetar3(isrc,numrel3(isrc)) ) & & cd3 = closs3(isrc,numrel3(isrc)) do i =2,numrel3(isrc) if ( area > wetar3(isrc,i-1) .and. & & area <= wetar3(isrc,i) ) then zz1 = area - wetar3(isrc,i-1) zz2 = wetar3(isrc,i) - wetar3(isrc,i-1) cc1 = closs3(isrc,i) - closs3(isrc,i-1) cd3 = closs3(isrc,i-1) + zz1 / zz2 * cc1 endif enddo endif endif ! ! Compute discharge for each types (3, 4, 5, 7, 8) ! if (mnksrc(7,isrc) == 3) then ! ! Compute discharge for type 3 ! call n_and_m_to_nm(mnksrc(2,isrc), mnksrc(1,isrc), nmin, gdp) call n_and_m_to_nm(mnksrc(5,isrc), mnksrc(4,isrc), nmout, gdp) if (kfs(nmin)==1) then delzet = s0(nmin) - s0(nmout) if (delzet>0.0) then disch(isrc) = clcul(isrc)*arcul(isrc)*sqrt(2.*ag*delzet) else disch(isrc) = 0.0 endif endif endif if (mnksrc(7,isrc) == 4) then ! ! Compute discharge for type 4 ! if (kfs(nmin) == 1) then if (s0(nmin) >= s0(nmout)) then iswitch = 0 call cptdis(lundia ,ag ,area ,calfa(isrc), & & cd1 ,cd2 ,cd3 ,cleng(isrc), & & cmann(isrc),coefl ,disch(isrc),htcul(isrc), & & iflow ,iswitch ,poscul(isrc),rdis , & & s0(nmin) ,s0(nmout) ,wtcul(isrc) ,gdp ) disch(isrc) = rdis else disch(isrc) = 0.0 endif else write(lundia,'(a,i2,a)') & & 'intake point of one-way culvert nr.',isrc, & & ' is dry' endif endif if (mnksrc(7,isrc) == 5) then ! ! Compute discharge for type 5 ! iswitch = 0 if (s0(nmin) < s0(nmout)) then ! ! intake and outfall exchange; recompute area and energy loss ! iswitch = 1 swapval = nmin nmin = nmout nmout = swapval endif if (kfs(nmin) == 1) then call cptdis(lundia ,ag ,area ,calfa(isrc), & & cd1 ,cd2 ,cd3 ,cleng(isrc), & & cmann(isrc),coefl ,disch(isrc),htcul(isrc), & & iflow ,iswitch ,poscul(isrc),rdis , & & s0(nmin) ,s0(nmout) ,wtcul(isrc) ,gdp ) disch(isrc) = rdis else write(lundia,'(a,i2,a)') & & 'intake point of two-way culvert nr.',isrc, & & ' is dry' endif endif if (mnksrc(7,isrc) == 7) then ! ! compute loss coefficient and discharge for type 7 ! call n_and_m_to_nm(mnksrc(2,isrc), mnksrc(1,isrc), nmin, gdp) call n_and_m_to_nm(mnksrc(5,isrc), mnksrc(4,isrc), nmout, gdp) if (kfs(nmin) == 1 .or. kfs(nmout) == 1) then ! ! User defined formula in DLL ! Input parameters are passed via dll_reals/integers/strings-arrays ! if (max_reals < 7) then write(errmsg,'(a,a,a)') 'Insufficient space to pass real values to culvert library.' call prterr (lundia,'U021', trim(errmsg)) call d3stop(1, gdp) endif dll_reals( 1) = real(timsec ,hp) dll_reals( 2) = real(s0(nmin) ,hp) dll_reals( 3) = real(s0(nmout),hp) dll_reals( 4) = real(-dps(nmin) ,hp) dll_reals( 5) = real(-dps(nmout),hp) dll_reals( 6) = real(disch(isrc),hp) dll_reals( 7) = real(ag,hp) ! if (max_integers < 8) then write(errmsg,'(a,a,a)') 'Insufficient space to pass integer values to culvert library.' call prterr (lundia,'U021', trim(errmsg)) call d3stop(1, gdp) endif dll_integers( 1) = nmin dll_integers( 2) = mnksrc(2,isrc) dll_integers( 3) = mnksrc(1,isrc) dll_integers( 4) = kfs(nmin) dll_integers( 5) = nmout dll_integers( 6) = mnksrc(5,isrc) dll_integers( 7) = mnksrc(4,isrc) dll_integers( 8) = kfs(nmout) ! if (max_strings < 2) then write(errmsg,'(a,a,a)') 'Insufficient space to pass strings to culvert library.' call prterr (lundia,'U021', trim(errmsg)) call d3stop(1, gdp) endif dll_strings( 1) = gdp%runid dll_strings( 2) = dll_usrfil(isrc) ! ! Initialisation of output variables of user defined transport formulae ! disch_dll = 0.0_hp pos1_dll = rmissval pos2_dll = rmissval message = ' ' call psemlun error_ptr = 0 error_ptr = perf_function_culvert(dll_handle(isrc), dll_function(isrc), & dll_integers , max_integers , & dll_reals , max_reals , & dll_strings , max_strings , & disch_dll , pos1_dll , & pos2_dll , message) call vsemlun if (error_ptr /= 0) then write(errmsg,'(a,a,a)') 'Cannot find function "',trim(dll_function(isrc)),'" in dynamic library.' call prterr (lundia,'U021', trim(errmsg)) call d3stop(1, gdp) endif if (message /= ' ') then write (lundia,'(a,a,a)') '*** ERROR Message from user defined culvert formula ',trim(dll_function(isrc)),' :' write (lundia,'(a,a )') ' ', trim(message) write (lundia,'(a )') ' ' call d3stop(1, gdp) endif ! ! Output parameters ! disch(isrc) = real(disch_dll,fp) ! ! Layer for intake point in specific layer or default over full depth (layer = 0) ! if (comparereal(pos1_dll,rmissval) /= 0) then pos = real(pos1_dll,fp) mnksrc(3,isrc) = intlay(lundia ,zmodel ,kfsmin(nmin) ,kfsmax(nmin) , & & pos ,zk ,dps(nmin) ,s0(nmin) , & & kmax ,thick ,isrc ,'intake ' ) else mnksrc(3,isrc) = 0 endif ! ! Layer for outfall point in specific layer or default over full depth (layer = 0) ! if (comparereal(pos2_dll,rmissval) /= 0) then pos = real(pos2_dll,fp) mnksrc(6,isrc) = intlay(lundia ,zmodel ,kfsmin(nmout),kfsmax(nmout), & & pos ,zk ,dps(nmout) ,s0(nmout) , & & kmax ,thick ,isrc ,'outfall' ) else mnksrc(6,isrc) = 0 endif endif endif if (mnksrc(7, isrc) == 8) then ! ! Compute discharge for type 8 ! call n_and_m_to_nm(mnksrc(2,isrc), mnksrc(1,isrc), nmin, gdp) call n_and_m_to_nm(mnksrc(5,isrc), mnksrc(4,isrc), nmout, gdp) hin = max (0.0_fp,s0(nmin )-poscul(isrc)) hout = max (0.0_fp,s0(nmout)-poscul(isrc)) height = 0.5_fp * ( hin + hout) height = min ( height , htcul(isrc) ) area = height * wtcul(isrc) ! ! vertical position for intake point: ! kin = intlay(lundia ,zmodel ,kfsmin(nmin) ,kfsmax(nmin) , & & poscul(isrc),zk ,dps(nmin) ,s0(nmin) , & & kmax ,thick ,isrc ,'intake ' ) mnksrc(3,isrc) = kin ! ! vertical position for outfall point: ! kout = intlay(lundia ,zmodel ,kfsmin(nmout),kfsmax(nmout), & & poscul(isrc),zk ,dps(nmout) ,s0(nmout) , & & kmax ,thick ,isrc ,'outfall' ) mnksrc(6,isrc) = kout rhorefin = 0.5_fp * thick(1) * rhow rhorefout = rhorefin do k = 2, kin rhorefin = rhorefin + 0.5_fp*(thick(k-1)+thick(k)) * rhow enddo do k = 2, kout rhorefout = rhorefout + 0.5_fp*(thick(k-1)+thick(k)) * rhow enddo delpres = sumrho(nmin,kin)*hin/rhorefin - sumrho(nmout,kout)*hout/rhorefout disch(isrc) = clcul(isrc)*area*sqrt(2.0_fp*ag*abs(delpres)) if (delpres < 0.0_fp) then ! ! intake and outfall exchange; recompute area and energy loss ! swapval = nmout nmout = nmin nmin = swapval swapval = kout kout = kin kin = swapval hin = max (0.0_fp,s0(nmin )-poscul(isrc)) hout = max (0.0_fp,s0(nmout)-poscul(isrc)) height = min ( hin , htcul(isrc) ) area = height * wtcul(isrc) mnksrc(3,isrc) = kin mnksrc(6,isrc) = kout disch(isrc) = - disch(isrc) endif if (kfs(nmin) /= 1) then disch(isrc) = 0.0_fp endif endif ! ! For all types: ! voldis(isrc) = voldis(isrc) + disch(isrc)*hdt enddo end subroutine culver integer function intlay (lundia ,zmodel ,kfsmin ,kfsmax , & & pos ,zk ,dps ,s0 , & & kmax ,thick ,isrc ,intake ) !----- 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$ ! $HeadURL$ !!--description----------------------------------------------------------------- ! Computes the discharge relation through a culvert. ! !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- use precision ! implicit none ! ! Function arguments ! logical , intent(in) :: zmodel ! True = zmodel, False = otherwise (sigma) integer , intent(in) :: isrc ! Number of culvert integer , intent(in) :: kfsmax ! Description and declaration in esm_alloc_int.f90 integer , intent(in) :: kfsmin ! Description and declaration in esm_alloc_int.f90 integer , intent(in) :: kmax ! Number of layers integer , intent(in) :: lundia ! Diagnostic file real(prec) , intent(in) :: dps ! Bathymetry [m below ref. level] real(fp) , intent(in) :: pos ! Level of culvert [m above ref. level] real(fp) , intent(in) :: s0 ! Water level [m above ref. level] real(fp), dimension(kmax) , intent(in) :: thick ! Description and declaration in esm_alloc_real.f90 real(fp), dimension(0:kmax) , intent(in) :: zk ! See description and declaration of sig in esm_alloc_real.f90 character(7) , intent(in) :: intake ! 'intake' or 'outfall' ! ! Local variables ! integer :: k real(fp) :: h0 real(fp) :: zdown real(fp) :: zup ! !! executable statements ------------------------------------------------------- ! intlay=0 if (zmodel) then do k=kfsmin, kfsmax if (pos >= zk(k-1) .and. pos < zk(k)) intlay = k enddo else h0 = real(dps,fp)+s0 zup = s0 do k=1,kmax zdown = zup - h0 *thick(k) if (pos >= zdown .and. pos < zup) intlay = k zup = zdown enddo endif if (pos < -real(dps,fp) ) then write(lundia,'(a,a,a,i2,a,f5.2,a,f5.2)') & & 'vertical position of culvert ',trim(intake),' nr. ',isrc, & & ' (',pos,') is below the bottom ',-real(dps,fp) endif end function intlay