!----- 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: wrwaq.F90 43424 2015-12-04 17:30:45Z kernkam $ ! $HeadURL: https://repos.deltares.nl/repos/ds/trunk/additional/unstruc/src/wrwaq.F90 $ module wrwaq !include preprocessing flags from autotools #ifdef HAVE_CONFIG_H #include "config.h" #endif use unstruc_files implicit none !!--copyright------------------------------------------------------------------- ! Copyright (c) 2007, 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. !!--description----------------------------------------------------------------- ! This module contains the routines for writing WAQ input files !!--declarations---------------------------------------------------------------- contains function openwaqbinfile(filename) result (lun) character(len=*), intent(in) :: filename !< Output filename. integer :: lun integer, external :: numuni ! ! WARNING: WAQ input files must be written using form=binary ! instead of unformatted. ! Although it is not standard Fortran ! lun = numuni() #ifdef HAVE_FC_FORM_BINARY open ( lun , file=filename , form = 'binary', status = 'replace') #else ! standardized way if binary is not available open ( lun , file=filename , form = 'unformatted' , access='stream', status = 'replace') #endif call reg_file_open(lun, filename) call mess(LEVEL_INFO, 'Opened file : ', filename) end function openwaqbinfile function openasciifile(filename) result (lun) character(len=*), intent(in) :: filename integer :: lun integer, external :: numuni ! ! NOTE: Opens a simple ASCII-file. Function is intended only ! to isolate newlun_nogdp dependency. ! lun = numuni() #ifdef HAVE_FC_FORM_BINARY open ( lun , file=filename , SHARED ) #else open ( lun , file=filename , access='stream') #endif call reg_file_open(lun, filename) call mess(LEVEL_INFO, 'Opened file : ', filename) end function openasciifile ! !------------------------------------------------------------------------------ !> Write ASCII or binary pointer file for WAQ. subroutine wrwaqpoi(ifrmto, noq, filename, ascii) implicit none integer , intent(in) :: noq !< Nr. of linkages (pointers) between computational cells. integer, dimension(:,:), intent(in) :: ifrmto !< Pointer table with all linkages. !! ifrmto(1,:) = 'from' cell number !! ifrmto(2,:) = 'to' cell number !! ifrmto(3,:) = 'from-1' cell number !! ifrmto(4,:) = 'to+1' cell number logical , intent(in) :: ascii !< Produce ascii file or not (then binary). character(*) , intent(in) :: filename !< Name for output pointer file. ! ! Local variables ! integer :: i integer :: lunout integer :: q ! !! executable statements ------------------------------------------------------- ! if (ascii) then ! ! ascii output ! lunout = openasciifile(filename) do q = 1,noq write(lunout,'(4i10)') ( ifrmto(i,q), i=1,4 ) enddo close(lunout) else ! ! binary output ! lunout = openwaqbinfile(filename) write(lunout) ( ( ifrmto(i,q), i=1,4 ), q=1,noq ) close(lunout) endif end subroutine wrwaqpoi ! !------------------------------------------------------------------------------ !> Write (binary) from/to length file for DelWAQ. subroutine wrwaqlen(noq, lenex, filename, ascii) use precision ! implicit none ! ! Global variables ! integer , intent(in) :: noq !< Nr. of linkages (pointers) between computational cells. double precision , intent(in) :: lenex(2, noq) !< Dispersion half-lengths of computational cells, segment !! centre to exchange point. (2 values: from/to direction) logical , intent(in) :: ascii !< Produce ascii file or not (then binary). character(*) , intent(in) :: filename !< Output filename. ! ! Local variables ! integer :: i integer :: lunout integer :: q ! !! executable statements ------------------------------------------------------- ! if (ascii) then ! ! ascii output ! lunout = openasciifile(filename) do q = 1,noq write(lunout,'(i10,2f18.8)') q, ( lenex(i,q), i=1,2 ) enddo close(lunout) else ! ! binary output ! lunout = openwaqbinfile(filename) write(lunout) noq write(lunout) (( real(lenex(i,q),sp), i=1,2 ), q=1,noq ) close(lunout) endif end subroutine wrwaqlen ! !------------------------------------------------------------------------------ !> Write (binary) horizontal surface file for DelWAQ. subroutine wrwaqsrf(srf, nosegl, nolay, filename, ascii) use precision implicit none ! ! Global variables ! integer , intent(in) :: nosegl !< Nr. of segment per layer. integer , intent(in) :: nolay !< Nr. of layers. double precision, dimension(:) , intent(in) :: srf !< Horizontal surfaces of computational cells. (size nosegl) logical , intent(in) :: ascii !< Produce ascii file or not (then binary). character(*) , intent(in) :: filename !< Output filename. character(256) :: fileold !< Old output filename. ! ! Local variables ! integer :: is, il integer :: lunout ! !! executable statements ------------------------------------------------------- ! if (ascii) then ! ! ascii output ! lunout = openasciifile(filename) ! nmax, mmax, nosegl, nosegl, nosegl, 0.0 ! noseg in all places is a dummy: no aggr yet, and no three m-n-k dimensions write (lunout, '(i10)') 0 do il = 1,nolay do is = 1,nosegl write(lunout,'(i10,f18.8)') is + (il - 1) * nosegl , srf(is) enddo enddo close(lunout) else ! ! binary output ! ! new style surf file for NGHS WAQ UI lunout = openwaqbinfile(filename) write(lunout) 0 write(lunout) ((real(srf(is),sp), is=1,nosegl ), il=1,nolay) close(lunout) ! temporaraly keep writing then old style surf file for WAQ_GUI fileold = trim(filename)//'old' lunout = openwaqbinfile(fileold) write(lunout) nosegl, 1, nosegl, nosegl, nosegl, 0.0 write(lunout) (real(srf(is),sp), is=1,nosegl ) close(lunout) endif end subroutine wrwaqsrf ! !------------------------------------------------------------------------------ !> Write ASCII attributes file for WAQ. subroutine wrwaqatr(nosegl, nolay, kmk1, kmk2, filename) implicit none integer , intent(in) :: nosegl !< Nr. of segments per layer integer , intent(in) :: nolay !< Nr. of layers integer, dimension(:), intent(in) :: kmk1 !< First WAQ segment features at start of calculation integer, dimension(:), intent(in) :: kmk2 !< Second WAQ segment features at start of calculation character(*) , intent(in) :: filename !< Name for output pointer file. ! ! Local variables ! integer :: il, is integer :: lunatr character( 2 ) kenout(nosegl) !! this is now allocated on the stack !!! ! !! executable statements ------------------------------------------------------- ! call newfil(lunatr, trim(filename)) write ( lunatr , '(a)' ) ' ; DELWAQ_COMPLETE_ATTRIBUTES' write ( lunatr , '(a)' ) ' 2 ; two blocks with input ' write ( lunatr , '(a)' ) ' 1 ; number of attributes, they are :' write ( lunatr , '(a)' ) ' 1 ; ''1'' is active ''0'' is not' write ( lunatr , '(a)' ) ' 1 ; data follows in this file ' write ( lunatr , '(a)' ) ' 1 ; all data is given without defaults' do il = 1,nolay write ( lunatr , * ) ' ; layer: ',il do is = 1, nosegl kenout(is) = ' ' write( kenout(is), '(I2)' ) kmk1( is + (il - 1) * nosegl ) enddo write ( lunatr, '(500a2)' ) kenout enddo write ( lunatr , '(a)' ) ' 1 ; number of attributes, they are :' write ( lunatr , '(a)' ) ' 2 ; ''1'' has surface ''3'' has bottom' write ( lunatr , '(a)' ) ' ; ''0'' has both ''2'' has none ' write ( lunatr , '(a)' ) ' 1 ; data follows in this file ' write ( lunatr , '(a)' ) ' 1 ; all data is given without defaults' do il = 1,nolay write ( lunatr , * ) ' ; layer: ',il do is = 1, nosegl kenout(is) = ' ' write( kenout(is), '(I2)' ) kmk2( is + (il - 1) * nosegl ) enddo write ( lunatr, '(500a2)' ) kenout enddo write ( lunatr , '(a)' ) ' 0 ; no time dependent attributes' close(lunatr) end subroutine wrwaqatr ! !------------------------------------------------------------------------------ !> Write (binary) LGA-numbering file for DelWAQ. !! Mapping from flow nodes to WAQ segments. subroutine wrwaqlga(iapnt, nmax, mmax, nosegl, nolay, noq1, noq2, noq3, & filename, ascii) implicit none ! ! Global variables ! integer, intent(in) :: iapnt(nmax*mmax) !< Mapping from flow nodes to WAQ segments. integer, intent(in) :: nmax, mmax !< Nr. of rows and columns in flow grid. integer, intent(in) :: nosegl !< Nr. of WAQ segments per layer. integer, intent(in) :: nolay !< Nr. of WAQ layers. integer, intent(in) :: noq1, noq2, noq3 !< Nr. of exchange points in column, row and layer direction. logical, intent(in) :: ascii !< Produce ascii file or not (then binary). character(*), intent(in) :: filename !< Output filename. ! ! Global variables ! integer :: lunout, i ! !! executable statements ------------------------------------------------------- ! if (ascii) then ! ! ascii output ! lunout = openasciifile(filename) write (lunout, '(a,i4,a,i4,a,i3)') 'Flow: nmax=', nmax, ', mmax=', mmax, ', nosegl=', nosegl write (lunout, '(a,i3,a,i4,a,i4,a,i3)') 'WAQ: nolay=', nolay, ', noq1=', noq1, ', noq2=', noq2, ', noq3=', noq3 do i = 1,nmax write(lunout, '(i10,i10)') i, iapnt(i) enddo close(lunout) else ! ! binary output ! lunout = openwaqbinfile(filename) write ( lunout ) nmax, mmax, nosegl, nolay, noq1, noq2, noq3 write ( lunout ) iapnt(1:nmax) close ( lunout ) end if end subroutine wrwaqlga ! !------------------------------------------------------------------------------ !> Write (binary) exchange file(s) for DelWAQ: area and fluxes. subroutine wrwaqbin(itim, quant, nquant, filename, ascii, lunout) use precision implicit none ! ! Global variables ! integer , intent(in) :: itim !< Time for new data block integer , intent(in) :: nquant !< Size of quant(ity) array. double precision, dimension(:) , intent(in) :: quant !< Quantity array to be written. logical , intent(in) :: ascii !< Produce ascii file or not (then binary). character(*) , intent(in) :: filename !< Output filename (only used if lunout not connected yet). integer , intent(inout) :: lunout !< File pointer for output file. Used if already connected, !! or set to new value for filename. ! ! Local variables ! integer :: q ! !! executable statements ------------------------------------------------------- ! if (ascii) then ! ! ascii output ! if (lunout<0) then lunout = openasciifile(filename) endif write(lunout,'(a,i10)') 'Time = ', itim do q = 1,nquant write(lunout,'(i10,f18.8)') q, quant(q) enddo !close(lunout) else ! ! binary output ! if (lunout<0) then lunout = openwaqbinfile(filename) endif write(lunout) itim write(lunout) ( real(quant(q),sp), q=1,nquant ) !close(lunout) endif end subroutine wrwaqbin ! !------------------------------------------------------------------------------ !!! routines below not used yet in Unstruc, copied from morphology branch of Bert Jagers. !> Write monitoring segments file for DelWAQ (each segment is a monitoring area). subroutine wrmonseg(noseg, filename) implicit none ! ! Global variables ! integer , intent(in) :: noseg character(*) , intent(in) :: filename !< Output filename. ! ! Local variables ! integer :: lunout integer :: s ! !! executable statements ------------------------------------------------------- ! lunout = openasciifile(filename) write(lunout,'(i5)') noseg do s = 1,noseg write(lunout,'(a,i4,a,i5)') '''Segment ',s,''' 1 ',s enddo close(lunout) end subroutine wrmonseg ! !------------------------------------------------------------------------------ !> Write NROFSEGM.DAT file for DelWAQ. subroutine wr_nrofseg(noseg, filename) implicit none ! ! Global variables ! integer , intent(in) :: noseg character(*) , intent(in) :: filename !< Output filename. ! ! Local variables ! integer :: lunout ! !! executable statements ------------------------------------------------------- ! lunout = openasciifile(filename) write(lunout,'(i12,a)') noseg,' ; number of segments' close(lunout) end subroutine wr_nrofseg ! !------------------------------------------------------------------------------ !> Write NROFEXCH.DAT file for DelWAQ. subroutine wr_nrofexch(noq1, noq2, noq3, filename) implicit none ! ! Global variables ! integer , intent(in) :: noq1 integer , intent(in) :: noq2 integer , intent(in) :: noq3 character(*) , intent(in) :: filename !< Output filename. ! ! Local variables ! integer :: lunout ! !! executable statements ------------------------------------------------------- ! lunout = openasciifile(filename) write(lunout,'(3i12,a)') noq1, noq2, noq3, ' ; number of exchanges in three directions' close(lunout) end subroutine wr_nrofexch ! !------------------------------------------------------------------------------ end module wrwaq !> Module for coupling with WAQ. !! Currently only writing of WAQ-files. module waq use unstruc_messages implicit none logical, parameter :: waq_format_ascii = .false. !< For debugging: .true. produces ascii files. .false. the binary files used for DELWAQ input. type gd_waqpar ! ! ! ! doubles ! ! ! ! ! ! reals ! ! ! real(fp) :: mtimstep ! maximum explicit time step for whole simulation ! ! ! ! integers ! ! integer :: aggre ! 0: no aggregation (=active cells only in FM), 1: aggregation according to content of flaggr integer :: aggrel ! 0: no layer aggregation, 1: layer aggregation ! integer :: itwqff ! start time for writing binary waq files ! integer :: itwqfi ! time step for writing binary waq files ! integer :: itwqfl ! end time for writing binary waq files ! integer :: itim ! last writen time integer :: lunvol ! file unit number to an output file integer :: lunare ! file unit number to an output file integer :: lunflo ! file unit number to an output file integer :: lunsal ! file unit number to an output file ! integer :: luntem ! file unit number to an output file integer :: lunvdf ! file unit number to an output file integer :: luntau ! file unit number to an output file ! integer :: lunsrctmp ! file unit number to an output file ! integer :: lunwlk ! file unit number to an output file ! integer :: lunsrc ! file unit number to an output file ! integer :: lunkmk ! file unit number to an output file integer :: noseg ! number of WAQ segments integer :: nosegl ! number of WAQ segments per layer integer :: noq ! total number of WAQ exchanges integer :: noql ! number of horizontal WAQ exchanges per layer integer :: noq12 ! number of horizontal WAQ exchanges (excluding sink/sources) integer :: noq12s ! number of horizontal WAQ exchanges (including sink/sources) integer :: numsrc ! number of sinks/sources integer :: numsrcbnd ! number of sinks/sources that are boundary conditions integer :: numsrcwaq ! number of adition sources/sinks exchanges in waq (based on posible combinations) integer :: kmxnx ! maximum number of active layers integer :: kmxnxa ! maximum number of aggregated layers integer :: ndkxi ! nr of internal flowcells (3D) ! integer :: nobrk ! number of breakpoints in loads file ! integer :: nowalk ! number of walking discharges ! ! ! integer :: cfoutset ! 'DataCFtoWQ' ! integer :: wqinset ! 'DataWQtoCF' ! integer :: wqiinset ! 'DataWQtoWQI' ! integer :: wqioutset ! 'DataWQItoWQ' ! ! pointers ! ! integer , dimension(:) , pointer :: iwlk ! walkings integer , allocatable :: ilaggr(:) ! layer aggregation pointer integer , allocatable :: ilaggrhyd(:) ! layer aggregation table for hyd-file integer , allocatable :: isaggr(:) ! segment aggregation pointer integer , allocatable :: iapnt (:) ! flow-to-waq pointer integer , allocatable :: iqaggr(:) ! exchange aggregation pointer integer , allocatable :: iqwaggr(:) ! exchange aggregation pointer for the vertical integer , allocatable :: ifrmto(:,:) ! from-to pointer table integer , allocatable :: ifrmtosrc(:,:) ! from-to pointer table for sources integer , allocatable :: ifsmin(:) ! first active layer (counting from the top) (z-model) integer , allocatable :: ifsmax(:) ! maximum active layer (counting from the top) (z-model) integer , allocatable :: nosega(:) ! no of segments aggregated into WAQ segments integer , allocatable :: kmk1(:) ! First WAQ segment features at start of calculation (1 is active 0 is not) integer , allocatable :: kmk2(:) ! Second WAQ segment features at start of calculation (1 surface, 3 bottom, 0 both, 2 neither) ! integer , dimension(:) , pointer :: ksrwaq ! stored value of nr of layers source locations ! integer , dimension(:) , pointer :: lunsed ! file unit numbers to sediment output files ! real(fp), dimension(:,:), pointer :: quwaq ! Cumulative qxk ! real(fp), dimension(:,:), pointer :: qvwaq ! Cumulative qyk ! real(fp), dimension(:,:), pointer :: qwwaq ! Cumulative qzk ! real(fp), dimension(:) , pointer :: discumwaq ! Cumulated sources m3/s*nstep double precision, allocatable :: horsurf(:) ! horizontal surfaces of segments double precision, allocatable :: vol(:) ! WAQ (aggregated) volumes double precision, allocatable :: sal(:) ! WAQ (aggregated) salinity double precision, allocatable :: tau(:) ! WAQ (aggregated) taus double precision, allocatable :: vdf(:) ! WAQ (aggregated) vertical diffusion ! real(sp), allocatable :: vol2 ! WAQ volume at end of step ! real(sp), allocatable :: sag ! WAQ segment aggregator ! real(sp), allocatable :: sag2 ! WAQ segment aggregator2 double precision, allocatable :: qag(:) ! WAQ (aggregated) flux double precision, allocatable :: area(:) ! WAQ (aggregated) exchange areas ! real(sp), dimension(:) , pointer :: loads ! Value of the loads at last step ! ! ! ! logicals ! ! ! logical :: first_cf ! first time of writing to binary waq files online ! logical :: firsttime ! first time of writing to binary waq files ! logical :: firstwaq ! skip the first time of getting bedlevel from WAQ ! logical :: waqfil ! Write binary waq files ! logical :: waqol = .false. ! flag for online coupling via WAQ files ! ! ! ! characters ! ! character(256) :: flaggr ! Name of input aggregation file end type gd_waqpar type(gd_waqpar) :: waqpar !-- Start subroutines ------------- contains subroutine reset_waq() ! !! executable statements ------------------------------------------------------- ! implicit none call close_and_reset(waqpar%lunvol) call close_and_reset(waqpar%lunsal) call close_and_reset(waqpar%lunare) call close_and_reset(waqpar%lunflo) call close_and_reset(waqpar%luntau) call close_and_reset(waqpar%lunvdf) contains subroutine close_and_reset(lun) implicit none integer, intent(inout) :: lun if (lun > 0) then call doclose(lun) end if lun = -1 end subroutine close_and_reset end subroutine reset_waq ! !------------------------------------------------------------------------------ !> Write the coupled hydrodynamics information file (.hyd) subroutine waq_wri_hyd() use unstruc_version_module, only: unstruc_version_full use unstruc_files use m_flowparameters use m_flowtimes use m_flow use m_flowexternalforcings use m_flowgeom use unstruc_model implicit none ! ! Local variables ! integer :: lunhyd, lunattr, ierr character(len=255) :: filename, stmp character tex*80, datetime*20 integer :: i, ibnd, isrc, kk1, kk2 integer :: itdate, julday, idatum, itijd, iyea,imon,iday, ihou,imin,isec double precision :: anl double precision :: x1, y1, x2, y2 double precision, parameter :: rmissval = -999.0d0 ! !! executable statements ------------------------------------------------------- ! filename = defaultFilename('hyd') call newfil(lunhyd, trim(filename)) call dateandtimenow(iyea,imon,iday,ihou,imin,isec) write(datetime ,'(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') iyea, '-', imon, '-', iday, ', ', ihou, ':', imin, ':', isec tex = 'file-creation-date ' //datetime write(lunhyd,'(a/a)') 'file-created-by '//trim(unstruc_version_full), tex write (lunhyd, '(a,a)') 'task ','full-coupling' if (layertype==LAYTP_SIGMA) then ! all sigma layers write (lunhyd, '(a,a)') 'geometry ','unstructured' elseif (layertype==LAYTP_Z) then ! all z layers write (lunhyd, '(a,a)') 'geometry ','unstructured z-layers' elseif (layertype==LAYTP_LEFTSIGMA) then ! mixed sigma/z layers write (lunhyd, '(a,a)') 'geometry ','unstructured left-sigma-layers' elseif (layertype==LAYTP_LEFTZ) then ! mixed sigma/z layers write (lunhyd, '(a,a)') 'geometry ','unstructured left-z-layers' else ! other? write (lunhyd, '(a,a)') 'geometry ','unstructured other' endif !if ( ! TODO: include check on aggregation if ( waqpar%aggre <= 0 ) then write (lunhyd, '(a,a)') 'horizontal-aggregation ','no' else write (lunhyd, '(a,a)') 'horizontal-aggregation ','yes' end if !! Time block read (refdat, '(i8)') itdate call juldat ( itdate, julday ) write ( lunhyd , '(A,I8,A )' ) 'reference-time ''',itdate,'000000''' call timdat( julday, tstart_user, idatum, itijd ) write ( lunhyd , '(A,I8,I6.6,A)') 'hydrodynamic-start-time ''',idatum,itijd,'''' call timdat( julday, tstop_user, idatum, itijd ) write ( lunhyd , '(A,I8,I6.6,A)') 'hydrodynamic-stop-time ''',idatum,itijd,'''' call timdat( julday, ti_waq, idatum, itijd ) ! TODO: niet ti_waq maar een soort comp.dt write ( lunhyd , '(A,I6.6,A )' ) 'hydrodynamic-timestep ''00000000',itijd,'''' write ( lunhyd , '(A,I8,A )' ) 'conversion-ref-time ''',itdate,'000000''' call timdat( julday, tstart_user, idatum, itijd ) write ( lunhyd , '(A,I8,I6.6,A)') 'conversion-start-time ''',idatum,itijd,'''' call timdat( julday, tstop_user, idatum, itijd ) write ( lunhyd , '(A,I8,I6.6,A)') 'conversion-stop-time ''',idatum,itijd,'''' call timdat( julday, ti_waq, idatum, itijd ) idatum = idatum - itdate write ( lunhyd , '(A,I8.8,I6.6,A )' ) 'conversion-timestep ''',idatum,itijd,'''' !! Grid dimensions block write (lunhyd, '(a,i10)') 'grid-cells-first-direction ', ndxi write (lunhyd, '(a,i10)') 'grid-cells-second-direction ', 1 write (lunhyd, '(a,i10)') 'number-hydrodynamic-layers ', waqpar%kmxnx write (lunhyd, '(a,i10)') 'number-horizontal-exchanges ', waqpar%noq12s ! Including open boundaries if (waqpar%kmxnxa > 1) then write (lunhyd, '(a,i10)') 'number-vertical-exchanges ', waqpar%noq - waqpar%noq12s else write (lunhyd, '(a,i10)') 'number-vertical-exchanges ', 0 end if write (lunhyd, '(a,i10)') 'number-water-quality-segments-per-layer ', waqpar%nosegl write (lunhyd, '(a,i10)') 'number-water-quality-layers ', waqpar%kmxnxa write ( lunhyd , '(A,A )' ) 'hydrodynamic-file ', ''''//trim(md_ident)//'''' if ( waqpar%aggre <= 0 ) then write ( lunhyd , '(A,A )' ) 'aggregation-file ', 'none' else write ( lunhyd , '(A,A )' ) 'aggregation-file ', ''''//trim(waqpar%flaggr)//'''' endif ! Grid and boundary conditions write (lunhyd, '(a,a)') 'boundaries-file ', ''''//trim(defaultFilename('bnd', prefixWithDirectory=.false.))//'''' write (lunhyd, '(a,a)') 'waqgeom-file ', ''''//trim(defaultFilename('waqgeom', prefixWithDirectory=.false.))//'''' !! MJ: strange, the real grid-indices-file is the lga file!!, not the bnd-file.... 'com-tut_fti_waq.lga' !! MJ: keep this until the old WAQ_GUI is updated write (lunhyd, '(a,a)') 'grid-indices-file ', ''''//trim(defaultFilename('bnd', prefixWithDirectory=.false.))//'''' write (lunhyd, '(a,a)') 'grid-coordinates-file ', ''''//trim(defaultFilename('waqgeom', prefixWithDirectory=.false.))//'''' !! Coupling files block write (lunhyd, '(a,a)') 'volumes-file ', ''''//trim(defaultFilename('vol', prefixWithDirectory=.false.))//'''' write (lunhyd, '(a,a)') 'areas-file ', ''''//trim(defaultFilename('are', prefixWithDirectory=.false.))//'''' write (lunhyd, '(a,a)') 'flows-file ', ''''//trim(defaultFilename('flo', prefixWithDirectory=.false.))//'''' write (lunhyd, '(a,a)') 'pointers-file ', ''''//trim(defaultFilename('poi', prefixWithDirectory=.false.))//'''' write (lunhyd, '(a,a)') 'lengths-file ', ''''//trim(defaultFilename('len', prefixWithDirectory=.false.))//'''' stmp = ' ' if (jasal > 0) then stmp = ''''//trim(defaultFilename('sal', prefixWithDirectory=.false.))//'''' else stmp = 'none' end if write (lunhyd, '(a,a)') 'salinity-file ', trim(stmp) write (lunhyd, '(a,a)') 'temperature-file ', 'none' if (waqpar%kmxnx == 1) then write (lunhyd, '(a,a)') 'vert-diffusion-file ', 'none' else write (lunhyd, '(a,a)') 'vert-diffusion-file ', ''''//trim(defaultFilename('vdf', prefixWithDirectory=.false.))//'''' endif ! new style surf file for NGHS WAQ UI write (lunhyd, '(a,a)') 'horizontal-surfaces-file ', ''''//trim(defaultFilename('srf', prefixWithDirectory=.false.))//'''' ! temporaraly keep writing then old style surf file for WAQ_GUI write (lunhyd, '(a,a)') 'surfaces-file ', ''''//trim(defaultFilename('srf', prefixWithDirectory=.false.))//'old''' !discharges-file 'com-f34.src' !chezy-coefficients-file 'com-f34.chz' write (lunhyd, '(a,a)') 'shear-stresses-file ', ''''//trim(defaultFilename('tau', prefixWithDirectory=.false.))//'''' !walking-discharges-file 'com-f34.wlk' !minimum-vert-diffusion ! upper-layer 0.0000E+00 ! lower-layer 0.0000E+00 ! interface-depth 0.0000E+00 !end-minimum-vert-diffusion !constant-dispersion ! first-direction 0.0000E+00 ! second-direction 0.0000E+00 ! third-direction 0.0000E+00 !end-constant-dispersion write (lunhyd, '(a,a)') 'attributes-file ', ''''//trim(defaultFilename('atr', prefixWithDirectory=.false.))//'''' if (layertype==LAYTP_Z) then ! all z layers ! do j = 1, mxlaydefs ! do k = 0, laymx(j) ! write ( 1972, '(I,I,F15.6)') j, k, zslay(k, j) ! enddo ! enddo if (mxlaydefs == 1) then write ( lunhyd , '(A, F15.6 )' ) 'z-layers-ztop ', zslay(laymx(1),1) write ( lunhyd , '(A, F15.6 )' ) 'z-layers-zbot ', zslay(0,1) else ! we do not yet know what to do with multiple defintions of z-layers write ( lunhyd , '(A )' ) 'z-layers-ztop multi-layerdefinition ' write ( lunhyd , '(A )' ) 'z-layers-zbot multi-layerdefinition ' endif elseif (layertype/=LAYTP_SIGMA) then ! all other options ! we do not yet know what to do with z/sigma combinations write ( lunhyd , '(A )' ) 'z-layers-ztop mixed-layerdefinition ' write ( lunhyd , '(A )' ) 'z-layers-zbot mixed-layerdefinition ' endif write ( lunhyd , '(A )' ) 'hydrodynamic-layers' do i = 1, waqpar%kmxnx write ( lunhyd , '(F15.6 )' ) 1.0e0 / real(waqpar%kmxnx) enddo write ( lunhyd , '(A )' ) 'end-hydrodynamic-layers' write ( lunhyd , '(A )' ) 'water-quality-layers ' anl = 1.000 do i = 1,waqpar%kmxnxa write ( lunhyd , '(F15.6 )' ) real(waqpar%ilaggrhyd(i)) enddo write ( lunhyd , '(A )' ) 'end-water-quality-layers' !discharges ! 2 14 1 '(14,2)' !end-discharges if(numsrc > 0) then ibnd = 0 if (nopenbndsect>0) ibnd = nopenbndlin(nopenbndsect) write ( lunhyd , '(A )' ) 'sink-sources' do isrc = 1,numsrc kk1 = ksrc(1,isrc) kk2 = ksrc(4,isrc) if (kk1 > 0) then x1 = xz(kk1) y1 = yz(kk1) else ibnd = ibnd + 1 kk1 = - ibnd x1 = rmissval y1 = rmissval endif if (kk2 > 0) then x2 = xz(kk2) y2 = yz(kk2) else ibnd = ibnd + 1 kk2 = - ibnd x2 = rmissval y2 = rmissval endif write ( lunhyd , '(3I10,4F18.6,2X,A)' ) isrc, kk1, kk2, x1, y1, x2, y2, trim(srcname(isrc)) enddo write ( lunhyd , '(A )' ) 'end-sink-sources' endif call doclose(lunhyd) end subroutine waq_wri_hyd ! !------------------------------------------------------------------------------ !> Write flow geometry information for use in WAQ-GUI. subroutine waq_wri_geom() use unstruc_netcdf use unstruc_files use unstruc_model, only: md_ident use m_flowtimes, only: rundat2 implicit none ! ! Local variables ! character(len=255) :: filename ! !! executable statements ------------------------------------------------------- ! filename = trim('DFM_DELWAQ_'//trim(md_ident)//trim(rundat2)//defaultFilename('waqgeom', prefixWithDirectory=.false.)) call unc_write_net_flowgeom(filename) ! TODO: AvD: this overall mixing up of DFM_DELWAQ_ and rundat2 is a mess, centralize this! ! add segment aggregation table to this file! label with delwaq_role = "segment_aggregation_table" and a "segdim" ! ierr = unc_open(filename, nf90_write, inetfile); call check_error(ierr, 'file '''//trim(filename)//'''') ! if (nerr_ > 0) return ! write waqpar%iapnt(1:ndxi) to a new table end subroutine waq_wri_geom ! !------------------------------------------------------------------------------ subroutine wrwaqgeom (filename, version_full, sferic, epsg, nr_nodes, xk, yk, zk, max_vertex, nr_elems, netelem, nr_edges, netlink, & nr_flowlinks, flowlink, xu, yu, ierr) ! !=============================================================================== ! Write the waqgeom netcdf file !=============================================================================== ! use netcdf use io_ugrid implicit none character(len=*) :: filename ! sobek_waqgeom.nc character(len=*) :: version_full ! SOBEK Version logical :: sferic ! No Sferic here, .false. integer :: epsg ! Projection type, 0 integer :: nr_nodes ! Total number of polypoints real(8) :: xk(nr_nodes), yk(nr_nodes), zk(nr_nodes) ! Poly Points integer :: max_vertex ! Maximaal aantal polypoints in polygon (here all 6) integer :: nr_elems ! Number of waq segments (isegtotal) integer :: netelem(max_vertex, nr_elems) ! idxPolygons(6, isegtotal) integer :: nr_edges ! Number of all line pieces to be drawn. integer :: netlink(2, nr_edges) ! From ploypoint to polypoint integer :: nr_flowlinks integer :: flowlink(2, nr_flowlinks) real(8) :: xu(nr_flowlinks), yu(nr_flowlinks) integer :: lundia = 1999 ! ! Local variables ! integer :: ierr integer :: igeomfile integer :: id_netelemdim ! netcdf id for mesh element dimension integer :: id_netnodedim ! netcdf id for node dimension integer :: id_netelemmaxnodedim ! netcdf id maximum number of vertices for an element (== 4 in Delft3D-FLOW) integer :: id_netnodex ! netcdf id for x-coordinate of node integer :: id_netnodey ! netcdf id for y-coordinate of node integer :: id_netnodez ! netcdf id for z-coordinate of node integer :: id_netlink, id_netelem, id_netlinkdim, id_netlinkptsdim integer :: id_flowlink, id_flowlinkdim, id_flowlinkptsdim, id_flowlinktype, id_flowlinkxu, id_flowlinkyu integer :: id_cfdim, id_cfmesh type(t_crs) :: crs character(len=20) :: datetime integer :: iyea, imon, iday, ihou, imin, isec ! !! executable statements ------------------------------------------------------- call iosdate(iyea,imon,iday) call iostime(ihou,imin,isec) write(datetime ,'(i4.4,a,i2.2,a,i2.2,a,i2.2,a,i2.2,a,i2.2)') iyea, '-', imon, '-', iday, ', ', ihou, ':', imin, ':', isec ! ierr = 0 crs%is_spherical = sferic crs%varvalue = epsg ! ! create or open the file ! ierr = nf90_create(filename, 0, igeomfile) !; call nc_check_err(lundia, ierr, "creating file", filename) if (ierr/=0) goto 9999 ! ! global attributes ! ierr = ug_addglobalatts(igeomfile) !, trim(version_full)) ierr = nf90_def_dim(igeomfile, 'dim' , 1, id_cfdim) ierr = nf90_def_var(igeomfile, 'mesh', nf90_int, id_cfdim, id_cfmesh) ierr = nf90_put_att(igeomfile, id_cfmesh, 'long_name' , 'Delft3D FM aggregated mesh') ierr = nf90_put_att(igeomfile, id_cfmesh, 'cf_role' , 'mesh_topology') ierr = nf90_put_att(igeomfile, id_cfmesh, 'topology_dimension' , '2 d') ierr = nf90_put_att(igeomfile, id_cfmesh, 'node_coordinates' , 'NetNode_x NetNode_y') ierr = nf90_put_att(igeomfile, id_cfmesh, 'face_node_connectivity', 'NetElemNode') ierr = nf90_put_att(igeomfile, id_cfmesh, 'edge_node_connectivity', 'NetLink') ierr = nf90_put_att(igeomfile, id_cfmesh, 'edge_face_connectivity', 'FlowLink') ! ! Dimensions ! ierr = nf90_def_dim(igeomfile, 'nNetNode' , nr_nodes, id_netnodedim) ierr = nf90_def_dim(igeomfile, 'nNetLink' , nr_edges, id_netlinkdim) ierr = nf90_def_dim(igeomfile, 'nNetLinkPts' , 2, id_netlinkptsdim) ! each edges has only a begin and end point ierr = nf90_def_dim(igeomfile, 'nNetElem' , nr_elems, id_netelemdim) ! number of elements ierr = nf90_def_dim(igeomfile, 'nNetElemMaxNode', max_vertex, id_netelemmaxnodedim) ! each element has exactly four vertices ! ierr = nf90_def_dim(igeomfile, 'nFlowLink' , nr_flowlinks, id_flowlinkdim) ierr = nf90_def_dim(igeomfile, 'nFlowLinkPts' , 2, id_flowlinkptsdim) ! each flow link has only a begin and end point ! ! Coordinates ! ierr = ug_add_coordmapping(igeomfile,crs) ! ierr = nf90_def_var(igeomfile, 'NetNode_x', nf90_double, id_netnodedim, id_netnodex) ierr = nf90_def_var(igeomfile, 'NetNode_y', nf90_double, id_netnodedim, id_netnodey) ierr = ug_addcoordatts(igeomfile, id_netnodex, id_netnodey, crs) ! ierr = nf90_def_var(igeomfile, 'NetNode_z', nf90_double, id_netnodedim, id_netnodez) ierr = nf90_put_att(igeomfile, id_netnodez, 'units', 'm') ierr = nf90_put_att(igeomfile, id_netnodez, 'positive', 'up') ierr = nf90_put_att(igeomfile, id_netnodez, 'standard_name', 'sea_floor_depth') ierr = nf90_put_att(igeomfile, id_netnodez, 'long_name', 'Bottom level at net nodes (flow element''s corners)') ! ierr = nf90_put_att(igeomfile, id_netnodez, 'coordinates', 'NetNode_x NetNode_y') ! ! Netlinks ! ierr = nf90_def_var(igeomfile, 'NetLink', nf90_int, (/ id_netlinkptsdim, id_netlinkdim /) , id_netlink) ierr = nf90_put_att(igeomfile, id_netlink, 'long_name' , 'link between two netnodes') ierr = nf90_put_att(igeomfile, id_netlink, 'start_index', 1) ! NetElements ierr = nf90_def_var(igeomfile, 'NetElemNode', nf90_int, (/ id_netelemmaxnodedim, id_netelemdim /) , id_netelem) ierr = nf90_put_att(igeomfile, id_netelem, 'long_name' , 'Net element defined by nodes') ierr = nf90_put_att(igeomfile, id_netelem, 'start_index', 1) ierr = nf90_put_att(igeomfile, id_netelem, '_FillValue' , 0) ! FLowlinks ierr = nf90_def_var(igeomfile, 'FlowLink', nf90_int, (/ id_flowlinkptsdim, id_flowlinkdim /) , id_flowlink) ierr = nf90_put_att(igeomfile, id_flowlink, 'long_name' , 'link/interface between two flow elements') ierr = nf90_put_att(igeomfile, id_flowlink, 'start_index', 1) ! ierr = nf90_def_var(igeomfile, 'FlowLinkType', nf90_int, (/ id_flowlinkdim /) , id_flowlinktype) ierr = nf90_put_att(igeomfile, id_flowlinktype, 'long_name' , 'type of flowlink') ierr = nf90_put_att(igeomfile, id_flowlinktype, 'valid_range' , (/ 1, 2 /)) ierr = nf90_put_att(igeomfile, id_flowlinktype, 'flag_values' , (/ 1, 2 /)) ierr = nf90_put_att(igeomfile, id_flowlinktype, 'flag_meanings', 'link_between_1D_flow_elements link_between_2D_flow_elements') ! ierr = nf90_def_var(igeomfile, 'FlowLink_xu', nf90_double, (/ id_flowlinkdim /) , id_flowlinkxu) ierr = nf90_def_var(igeomfile, 'FlowLink_yu', nf90_double, (/ id_flowlinkdim /) , id_flowlinkyu) ierr = ug_addcoordatts(igeomfile, id_flowlinkxu, id_flowlinkyu, crs) ierr = nf90_put_att(igeomfile, id_flowlinkxu, 'long_name' , 'x-Coordinate of velocity point on flow link.') ierr = nf90_put_att(igeomfile, id_flowlinkyu, 'long_name' , 'y-Coordinate of velocity point on flow link.') ! ierr = nf90_enddef(igeomfile) ! !=================================================================================================================== ! ! write data ! ierr = nf90_put_var(igeomfile, id_netnodex, xk(1:nr_nodes)) ierr = nf90_put_var(igeomfile, id_netnodey, yk(1:nr_nodes)) ierr = nf90_put_var(igeomfile, id_netnodez, zk(1:nr_nodes)) ! ierr = nf90_put_var(igeomfile, id_netlink, netlink, count=(/ 2, nr_edges /)) ierr = nf90_put_var(igeomfile, id_netelem, netelem, count=(/ max_vertex, nr_elems /)) ! ierr = nf90_put_var(igeomfile, id_flowlink, flowlink(:,1:nr_flowlinks)) ierr = nf90_put_var(igeomfile, id_flowlinkxu, xu(1:nr_flowlinks)) ierr = nf90_put_var(igeomfile, id_flowlinkyu, yu(1:nr_flowlinks)) ! ierr = nf90_sync(igeomfile) !; call nc_check_err(lundia, ierr, "sync file", filename) if (ierr/=0) goto 9999 ierr = nf90_close(igeomfile) !; call nc_check_err(lundia, ierr, "closing file", filename) 9999 continue ! end subroutine wrwaqgeom !> Write boundary information for use in WAQ-GUI. subroutine waq_wri_bnd() use m_flowgeom use network_data use m_flowexternalforcings use unstruc_files implicit none ! ! Local variables ! integer :: LL, L, Lf, n, i, istart, nci, k, n1, n2 integer :: ibnd, isrc, kk integer :: lunbnd, ierr character(len=255) :: filename double precision :: x1, y1, x2, y2, xn, yn integer, parameter :: waqmaxnamelen = 20 integer :: namelen ! !! executable statements ------------------------------------------------------- ! filename = defaultFilename('bnd') call newfil(lunbnd, trim(filename)) write(lunbnd, '(i8)') nopenbndsect + waqpar%numsrcbnd ! Nr of open boundary sections and sink sources. istart = 0 do i=1,nopenbndsect namelen = len_trim(openbndname(i)) write(lunbnd, '(a)') trim(openbndname(i)(1:min(namelen, waqmaxnamelen))) ! Section name write(lunbnd, '(i8)') nopenbndlin(i)-istart ! Nr of lins in section do LL=istart+1,nopenbndlin(i) L = openbndlin(LL) Lf = lne2ln(L) if (Lf <= 0 .or. Lf > lnx) then n = 0 x1 = 0d0 y1 = 0d0 x2 = 0d0 y2 = 0d0 else n = ln(1,Lf) if (kn(3,L) == 1) then ! 1D link ! n1 = abs(lne(1,L)) ! external 1D flow node ! n2 = abs(lne(2,L)) ! internal 1D flow node n1 = abs(ln(1,Lf)) ! external 1D flow node n2 = abs(ln(2,Lf)) ! internal 1D flow node ! For 1D links: produce fictious 'cross/netlink' call normalout(xz(n1), yz(n1), xz(n2), yz(n2), xn, yn) xn = wu(Lf)*xn yn = wu(Lf)*yn x1 = .5d0*(xz(n1)+xz(n2)) - .5d0*xn y1 = .5d0*(yz(n1)+yz(n2)) - .5d0*yn x2 = .5d0*(xz(n1)+xz(n2)) + .5d0*xn y2 = .5d0*(yz(n1)+yz(n2)) + .5d0*yn else x1 = xk(kn(1,L)) y1 = yk(kn(1,L)) x2 = xk(kn(2,L)) y2 = yk(kn(2,L)) end if end if ! Bnd flow nodes as delwaq bnd segments (<0) ! and 2 points of net link write(lunbnd, '(i8,4f18.8)') -(n-ndxi), x1, y1, x2, y2 end do istart = nopenbndlin(i) end do ibnd = ndx - ndxi do isrc = 1, numsrc if (ksrc(1,isrc) == 0 .or. ksrc(4,isrc) == 0) then ! This is a boundary condition ibnd = ibnd + 1 if (ksrc(1,isrc) /= 0) then kk = ksrc(1,isrc) else kk = ksrc(4,isrc) endif write(lunbnd, '(a)') trim(srcname(isrc)) ! Section name write(lunbnd, '(i8)') 1 ! Nr of lins in section write(lunbnd, '(i8,4f18.8)') -(ibnd), xz(kk), yz(kk), xz(kk), yz(kk) endif enddo call doclose(lunbnd) end subroutine waq_wri_bnd ! !------------------------------------------------------------------------------ !> Writes all necessary time-independent model schematisation files for DelWAQ. !! (.len, .poi, .srf) subroutine waq_wri_model_files() use m_flowgeom use m_flow, only : Lnkx use unstruc_files, only : defaultFilename implicit none integer :: ierr ! ! Local variables ! ! !! executable statements ------------------------------------------------------- ! write(msgbuf, '(a)') 'In waq_wri_model_files().' call msg_flush() ! Automatically prepare aggration and mapping tables. call waq_prepare_aggr() ! Hyd file (hydrodynamic coupling) call waq_wri_hyd() ! (Unstructured) grid file call waq_wri_geom() ! Boundaries file call waq_wri_bnd() ! Pointer file (linkages) call waq_wri_poi(defaultFilename('poi')) ! Length file (dispersion lengths) call waq_wri_len(lnx, dx, acl, defaultFilename('len')) ! Surface file (horizontal surfaces) call waq_wri_srf(ndx2D, ndxi, ndx, ba, defaultFilename('srf')) ! Attributes file call waq_wri_atr(defaultFilename('atr')) end subroutine waq_wri_model_files ! !------------------------------------------------------------------------------ !> Writes all necessary time-dependent couple files for DelWAQ. !! (.are, .flo, .vol) !! !! Note: flow-related files (.are and .flo) are not written for first !! time. Thereafter, accumulated flux is associated with previous !! timestep. At the last time a dummy record is written in .are and .flo. subroutine waq_wri_couple_files(time) use m_flowtimes use m_flowgeom use m_flow use m_flowexternalforcings use unstruc_files, only: defaultFilename implicit none ! ! Global variables ! double precision, intent(in) :: time !< Current simulation time ! ! Local variables ! integer :: itim !< Time (seconds) since simulation start. integer, save :: itim_prev = -1 ! integer :: i, n, L, LL, nn, ifrctyp ! double precision :: wrc, cf, cfn, cz, frcn, ar, wa ! !! executable statements ------------------------------------------------------- ! write (msgbuf,*) 'In waq_wri_couple_files: #', it_waq, ', time: ', time it_waq = it_waq+1 call msg_flush() itim = nint(time) ! Volume file (volumes of computational cells) call waq_wri_vol(itim, defaultFilename('vol'), waqpar%lunvol) ! TODO: AvD: add a 'mode' 0/1/2 similar to Delft3D, so that we can write some quantities ! at *start* of *next* timestep instead of currently at the *end* of *current* timestep. ! Salinty file (salinity of computational cells) if (jasal > 0) then call waq_wri_sal(itim, defaultFilename('sal'), waqpar%lunsal) end if ! Taus file (contains taus at the bottom of computational cells) call gettaus(1) call waq_wri_tau(itim, defaultFilename('tau'), waqpar%luntau) !> Vertical diffusion file (contains vertical diffusion of computational cells) if (waqpar%kmxnxa > 1) then call waq_wri_vdf(itim, defaultFilename('vdf'), waqpar%lunvdf) end if ! For first step, do not write any flux-related quantities (note the return). if (it_waq == 1) then itim_prev = itim return end if ! Area file (flow areas) call waq_wri_are(itim_prev, defaultFilename('are'), waqpar%lunare) ! AvD: NOTE: A bit strange, au is *not* accumulated, defined at time1, but still printed ! in file with itim_prev. ! Flow file (discharges) call waq_wri_flo(itim_prev, int(ti_waq), defaultFilename('flo'), waqpar%lunflo) ! Write a dummy last record in area and flow file to make them complete. if (time == tstop_user) then au = 0d0 q1waq = 0d0 if (kmx > 0) then qwwaq = 0d0 end if if (numsrc > 0) then qsrcwaq = 0d0 ! Reset accumulated discharges end if ! Dummy area record call waq_wri_are(itim, defaultFilename('are'), waqpar%lunare) ! Dummy flow record call waq_wri_flo(itim, int(ti_waq), defaultFilename('flo'), waqpar%lunflo) end if q1waq = 0d0 ! Reset accumulated discharges if (kmx > 0) then qwwaq = 0d0 ! Reset accumulated discharges end if if (numsrc > 0) then qsrcwaq = 0d0 ! Reset accumulated discharges end if itim_prev = itim end subroutine waq_wri_couple_files ! !------------------------------------------------------------------------------ !> Prepare aggregation and mapping tables between flow grid and WAQ cells. !! Currently, default: one-to-one. subroutine waq_prepare_aggr() use m_flowgeom use m_flow use m_flowexternalforcings use m_alloc implicit none integer :: i, ilay, kb, ktx, vaglay integer :: lunvag, istat, ierr logical :: test_aggr, test_layeraggr call realloc(waqpar%iapnt, ndx, keepExisting=.false., fill=0) call realloc(waqpar%isaggr, ndkx, keepExisting=.false., fill=0) call realloc(waqpar%iqaggr, lnkx, keepExisting=.false., fill=0) call realloc(waqpar%iqwaggr, ndkx, keepExisting=.false., fill=0) ! Is there a DIDO aggregation? Otherwise set default aggregation inquire (file = "waqtest.dwq", exist = test_aggr)! no interface supported yet if (test_aggr) then ! Change to true to test DIDO aggregation waqpar%aggre = 1 waqpar%flaggr = "waqtest.dwq" call waq_read_dwq(ndxi, ndx, waqpar%iapnt, waqpar%flaggr) waqpar%nosegl = maxval(waqpar%iapnt) else ! no aggregation, create default one to one aggregation waqpar%aggre = 0 waqpar%flaggr = " " do i = 1, ndxi waqpar%iapnt(i) = i end do waqpar%nosegl = ndxi ! add pointer for boundary nodes do i = ndxi+1,ndx waqpar%iapnt(i) = -(i-ndxi) enddo end if ! Maximum nr of layers per node in flow ! MJ: this might not yet be the absolute total number of layers in a z-model, but we might not need it... if(kmx == 0) then waqpar%kmxnx = 1 else waqpar%kmxnx = maxval(kmxn) end if call realloc(waqpar%ilaggrhyd, waqpar%kmxnx, keepExisting=.false., fill=1) call realloc(waqpar%ilaggr, waqpar%kmxnx, keepExisting=.false., fill=0) waqpar%aggrel = 0 waqpar%kmxnxa = waqpar%kmxnx ! Is there a layer aggregation file? inquire (file = "waqtest.vag", exist = test_layeraggr)! no interface supported yet if (test_layeraggr) then call mess(LEVEL_INFO, 'Try to read vertical aggregation file for DELWAQ.') call oldfil(lunvag, "waqtest.vag") read (lunvag, *, iostat=istat) vaglay if (vaglay == -1) then call mess(LEVEL_INFO, 'Found -1 as number of layers in waqtest.vag. Will aggregate DELWAQ output to 2D.') waqpar%kmxnxa = 1 waqpar%ilaggrhyd(1) = waqpar%kmxnx elseif (vaglay /= waqpar%kmxnx) then call mess(LEVEL_WARN, 'Mismatch in number of layers in DELWAQ vag layer aggregation file in '''//trim("waqtest.vag")//'''.') call mess(LEVEL_WARN, 'no layer aggregation applied.') else read (lunvag, *, iostat=istat) waqpar%ilaggr(1:vaglay) if ( istat /= 0) then call mess(LEVEL_WARN, 'unable to read the layer aggregation pointer in '''//trim("waqtest.vag")//'''.') call mess(LEVEL_WARN, 'no layer aggregation applied.') else ! Check validity of the layer aggregation and create waqpar%ilaggrhyd ierr = 0 if (waqpar%ilaggr(1) == 1 ) then waqpar%kmxnxa = 1 do i = 2,vaglay if(waqpar%ilaggr(i-1) == waqpar%ilaggr(i)) then waqpar%ilaggrhyd(waqpar%kmxnxa) = waqpar%ilaggrhyd(waqpar%kmxnxa) + 1 elseif (waqpar%ilaggr(i) == waqpar%ilaggr(i-1) + 1) then waqpar%kmxnxa = waqpar%kmxnxa + 1 else ierr = 1 exit end if enddo else ierr = 1 endif if (ierr == 1) then ! No correct layer aggregation, default: no layer aggragation call mess(LEVEL_WARN, 'layer aggregation is incorrect. no layer aggregation applied.') waqpar%ilaggrhyd = 1 waqpar%kmxnxa = waqpar%kmxnx else waqpar%aggrel = 1 end if endif endif else waqpar%aggrel = 0 waqpar%kmxnxa = waqpar%kmxnx end if if (waqpar%aggrel == 1 .and. waqpar%kmxnxa == 1) then ! If layer aggregation results in one layer left, then waq output will be 2D from 2D FM data directly waqpar%aggrel = 0 end if if (waqpar%aggrel == 0) then do i=1, waqpar%kmxnxa waqpar%ilaggr(i) = i enddo endif waqpar%noseg = waqpar%nosegl * waqpar%kmxnxa call realloc(waqpar%nosega, waqpar%noseg, keepExisting=.false., fill=0) call realloc(waqpar%vol, waqpar%noseg, keepExisting=.false., fill=0d0) call realloc(waqpar%sal, waqpar%noseg, keepExisting=.false., fill=0d0) call realloc(waqpar%tau, waqpar%noseg, keepExisting=.false., fill=0d0) call realloc(waqpar%vdf, waqpar%noseg, keepExisting=.false., fill=0d0) call realloc(waqpar%kmk1, waqpar%noseg, keepExisting=.false., fill=0) call realloc(waqpar%kmk2, waqpar%noseg, keepExisting=.false., fill=0) call getkbotktopmax(ndxi,kb,ktx) waqpar%ndkxi = ktx ! Maximum internal 3D node call waq_make_aggr_seg() ! Prepare arrays for sinks and sources. call waq_prepare_src() ! allocate maximum possible number of exchanges before aggregation waqpar%noq12 = lnx * waqpar%kmxnxa if (waqpar%kmxnxa > 1) then waqpar%noq = waqpar%noq12 + waqpar%numsrcwaq + ndxi * waqpar%kmxnxa else waqpar%noq = waqpar%noq12 + numsrc end if call realloc(waqpar%ifrmto, (/ 4, waqpar%noq /), keepExisting=.false., fill = 0) call waq_make_aggr_lnk() call realloc(waqpar%ifrmto, (/ 4, waqpar%noq /), keepExisting=.true., fill = 0) call realloc(waqpar%qag, waqpar%noq, keepExisting=.false., fill = 0d0) call realloc(waqpar%area, waqpar%noq, keepExisting=.false., fill = 0d0) waqpar%noql = waqpar%noq12 / waqpar%kmxnxa end subroutine waq_prepare_aggr ! !------------------------------------------------------------------------------ !> Creates aggregation and mapping tables to WAQ cells. subroutine waq_make_aggr_seg() use m_flowgeom use m_flow use wrwaq implicit none integer, parameter :: kmktopbot = 0 !< Segment is both the highest and lowest layer at this node integer, parameter :: kmktop = 1 !< Segment is at the highest layer at this node integer, parameter :: kmkmiddle = 2 !< Segment is an intermediate layer at this node integer, parameter :: kmkbot = 3 !< Segment is at the highest layer at this node integer :: k, kk, kb, ktx, ibnd, iseg ! clear segment aggregation pointer waqpar%isaggr = 0 waqpar%kmk1 = 0 waqpar%kmk2 = 0 ! 2D do iseg = 1, ndx waqpar%isaggr(iseg) = waqpar%iapnt(iseg) if(waqpar%isaggr(iseg) > 0) then waqpar%kmk1(iseg) = 1 !< segment is active (note MJ: not tested yet for dry cells...) endif enddo ! 3D if (waqpar%kmxnxa > 1) then do k = 1, ndxi call getkbotktopmax(k,kb,ktx) do kk = ktx, kb, -1 iseg = waqpar%iapnt(k) + (waqpar%ilaggr(ktx - kk + 1) - 1) * waqpar%nosegl waqpar%isaggr(kk) = iseg waqpar%nosega(iseg) = waqpar%nosega(iseg) + 1 waqpar%kmk2(iseg) = kmkmiddle if (ktx == kb) then waqpar%kmk2(iseg) = kmktopbot else if (kk == ktx) then waqpar%kmk2(iseg) = kmktop else if (kk == kb) then waqpar%kmk2(iseg) = kmkbot endif waqpar%kmk1(iseg) = 1 !< segment is active enddo end do ! also aggregate boundary nodes in 3D! do k = ndxi + 1, ndx call getkbotktopmax(k,kb,ktx) do kk = ktx, kb, -1 waqpar%isaggr(kk) = -((k - ndxi) + (waqpar%ilaggr(ktx - kk + 1) - 1) * (ndx - ndxi)) enddo end do end if end subroutine waq_make_aggr_seg ! !------------------------------------------------------------------------------ !> Creates aggregation and mapping tables to WAQ links. subroutine waq_make_aggr_lnk() use m_flowgeom use m_flow use wrwaq implicit none integer :: dseg, dbnd, isrc integer :: L, LL, Lb, Ltx, kxL, ip, ip1, ip2, ip3, ip4, iq integer :: k, kk, kb, ktx ! first 2D waqpar%noq12 = 0 do L=1,lnx ip1 = waqpar%iapnt(ln(1,L)) ip2 = waqpar%iapnt(ln(2,L)) ip3 = 0 if(klnup(1,L) /= 0) ip3 = waqpar%iapnt(abs(klnup(1,L))) ! abs to correct single value weighting flag. ip4 = 0 if(klnup(4,L) /= 0) ip4 = waqpar%iapnt(abs(klnup(4,L))) ! abs to correct single value weighting flag. if ( (ip1 > 0 .or. ip2 > 0) .and. ip1 /= ip2 ) then ip = 0 if (waqpar%aggre == 1 .or. waqpar%aggrel == 1) then do iq = 1 , waqpar%noq12 if (waqpar%ifrmto(1,iq) == ip1 .and. waqpar%ifrmto(2,iq) == ip2) then ip = iq waqpar%ifrmto(3,iq) = 0 waqpar%ifrmto(4,iq) = 0 exit endif if (waqpar%ifrmto(1,iq) == ip2 .and. waqpar%ifrmto(2,iq) == ip1) then ip = -iq waqpar%ifrmto(3,iq) = 0 waqpar%ifrmto(4,iq) = 0 exit endif enddo endif if ( ip == 0 ) then waqpar%noq12 = waqpar%noq12 + 1 waqpar%ifrmto(1,waqpar%noq12) = ip1 waqpar%ifrmto(2,waqpar%noq12) = ip2 waqpar%ifrmto(3,waqpar%noq12) = ip3 waqpar%ifrmto(4,waqpar%noq12) = ip4 waqpar%iqaggr(L) = waqpar%noq12 else waqpar%iqaggr(L) = ip endif end if end do if (waqpar%kmxnxa == 1) then ! In 2D total number of equals number of vertical exchanges ! Add links from sink source here!! if (waqpar%numsrcwaq > 0) then do isrc = 1, waqpar%numsrcwaq waqpar%ifrmto(1,waqpar%noq12 + isrc) = waqpar%ifrmtosrc(1, isrc) waqpar%ifrmto(2,waqpar%noq12 + isrc) = waqpar%ifrmtosrc(2, isrc) waqpar%ifrmto(3,waqpar%noq12 + isrc) = 0 waqpar%ifrmto(4,waqpar%noq12 + isrc) = 0 enddo endif waqpar%noq12s = waqpar%noq12 + waqpar%numsrcwaq waqpar%noq = waqpar%noq12s else ! In 3D copy aggregated 2D exchanges to all layers do L=1,lnx call getLbotLtopmax(L,Lb,Ltx) do LL = Ltx, Lb, -1 ip = waqpar%iqaggr(L) if (ip == 0) then waqpar%iqaggr(LL) = 0 else waqpar%iqaggr(LL) = ip + sign((waqpar%ilaggr(Ltx - LL + 1) - 1) * waqpar%noq12,ip) iq = abs(waqpar%iqaggr(LL)) dseg = (waqpar%ilaggr(Ltx - LL + 1) - 1) * waqpar%nosegl dbnd = (waqpar%ilaggr(Ltx - LL + 1) - 1) * (ndx - ndxi + waqpar%numsrcbnd) ! current number of external links in FM, account for sinks sources here too! if (waqpar%ifrmto(1,iq) == 0) then if (waqpar%ifrmto(1,ip) > 0) waqpar%ifrmto(1,iq) = waqpar%ifrmto(1,ip) + dseg if (waqpar%ifrmto(1,ip) < 0) waqpar%ifrmto(1,iq) = waqpar%ifrmto(1,ip) - dbnd if (waqpar%ifrmto(2,ip) > 0) waqpar%ifrmto(2,iq) = waqpar%ifrmto(2,ip) + dseg if (waqpar%ifrmto(2,ip) < 0) waqpar%ifrmto(2,iq) = waqpar%ifrmto(2,ip) - dbnd if (waqpar%ifrmto(3,ip) > 0) waqpar%ifrmto(3,iq) = waqpar%ifrmto(3,ip) + dseg if (waqpar%ifrmto(4,ip) > 0) waqpar%ifrmto(4,iq) = waqpar%ifrmto(4,ip) + dseg end if end if end do end do waqpar%noq12 = waqpar%noq12 * waqpar%kmxnxa waqpar%noq = waqpar%noq12 ! Add links from sink source here!! ! No checking on doubles due to aggregation yet!! if (waqpar%numsrcwaq > 0) then do isrc = 1, waqpar%numsrcwaq waqpar%ifrmto(1,waqpar%noq12 + isrc) = waqpar%ifrmtosrc(1, isrc) waqpar%ifrmto(2,waqpar%noq12 + isrc) = waqpar%ifrmtosrc(2, isrc) waqpar%ifrmto(3,waqpar%noq12 + isrc) = 0 waqpar%ifrmto(4,waqpar%noq12 + isrc) = 0 enddo endif waqpar%noq12s = waqpar%noq12 + waqpar%numsrcwaq waqpar%noq = waqpar%noq12s ! And now the vertical exchanges and pointer (note: upward flow direction in FM is reversed for WAQ!) do k = 1, waqpar%nosegl do kk = 1, waqpar%kmxnxa - 1 iq = waqpar%noq + k + (kk - 1) * waqpar%nosegl waqpar%ifrmto(1,iq) = k + (kk - 1) * waqpar%nosegl waqpar%ifrmto(2,iq) = k + kk * waqpar%nosegl waqpar%ifrmto(3,iq) = max(k + (kk - 2) * waqpar%nosegl, 0) waqpar%ifrmto(4,iq) = 0 if(kk < waqpar%kmxnxa - 1) waqpar%ifrmto(4,iq) = k + (kk + 1) * waqpar%nosegl end do end do waqpar%noq = waqpar%noq + waqpar%nosegl * (waqpar%kmxnxa - 1) do k = 1, ndxi call getkbotktopmax(k,kb,ktx) do kk = ktx - 1, kb, -1 if(waqpar%ilaggr(ktx - kk + 1) /= waqpar%ilaggr(ktx - kk)) then waqpar%iqwaggr(kk) = waqpar%noq12s + waqpar%iapnt(k) + (waqpar%ilaggr(ktx - kk) - 1) * waqpar%nosegl end if enddo end do end if end subroutine waq_make_aggr_lnk ! !------------------------------------------------------------------------------ !> Prepare additional exchanges for waq to store sink/source discharges !! Currently, we allocate all posible combinations in an inlet-outlet situation !! since we doe not know which one are used subroutine waq_prepare_src() use m_flowgeom use m_flow use m_flowexternalforcings use m_alloc implicit none integer :: ibnd, nbnd, isrc, K, K1, K2, kk integer :: kmax1, kk1,kb1,ktx1 integer :: kmax2, kk2,kb2,ktx2 waqpar%numsrcbnd = 0 waqpar%numsrcwaq = 0 if (numsrc==0) return ! skip is no resources call realloc (ksrcwaq , numsrc, keepexisting=.true., fill=0 ) ! First determine the number of external sink/sources and the allocations needed do isrc = 1, numsrc ksrcwaq (isrc) = waqpar%numsrcwaq if (ksrc(1,isrc) == 0 .or. ksrc(4,isrc) == 0) then ! This is a boundary condition waqpar%numsrcbnd = waqpar%numsrcbnd + 1 waqpar%numsrcwaq = waqpar%numsrcwaq + waqpar%kmxnxa else ! This is a sink/source combination waqpar%numsrcwaq = waqpar%numsrcwaq + waqpar%kmxnxa * waqpar%kmxnxa endif enddo call realloc (waqpar%ifrmtosrc, (/ 2,waqpar%numsrcwaq /), keepexisting=.true., fill=0 ) call realloc (qsrcwaq, waqpar%numsrcwaq , keepexisting=.true., fill=0.0D0 ) nbnd = ndx - ndxi + waqpar%numsrcbnd ! total number of boudaries ibnd = ndx - ndxi ! starting number for sink source boundaries ! Create additional pointer for sink/sources do isrc = 1, numsrc if (ksrc(1,isrc) == 0 .or. ksrc(4,isrc) == 0) then ! This is a boundary condition ibnd = ibnd + 1 if (ksrc(1,isrc) /= 0) then kk1 = ksrc(1,isrc) do K=1,waqpar%kmxnxa waqpar%ifrmtosrc(1, ksrcwaq (isrc) + K) = waqpar%iapnt(kk1) + (K - 1) * waqpar%nosegl waqpar%ifrmtosrc(2, ksrcwaq (isrc) + K) = -ibnd - nbnd * (K - 1) enddo else kk2 = ksrc(4,isrc) do K=1,waqpar%kmxnxa waqpar%ifrmtosrc(1, ksrcwaq (isrc) + K) = - ibnd - nbnd * (K - 1) waqpar%ifrmtosrc(2, ksrcwaq (isrc) + K) = waqpar%iapnt(kk2) + (K - 1) * waqpar%nosegl enddo endif else ! This is a sink/source combination kk1 = ksrc(1,isrc) kk2 = ksrc(4,isrc) if (waqpar%kmxnxa > 1) then do K1=1,waqpar%kmxnxa do K2=1,waqpar%kmxnxa kk = ksrcwaq (isrc) + K1 + (K2 - 1) * waqpar%kmxnxa waqpar%ifrmtosrc(1, kk) = waqpar%iapnt(kk1) + (K1 - 1) * waqpar%nosegl waqpar%ifrmtosrc(2, kk) = waqpar%iapnt(kk2) + (K2 - 1) * waqpar%nosegl enddo enddo else waqpar%ifrmtosrc(1, ksrcwaq (isrc) + 1) = waqpar%iapnt(kk1) waqpar%ifrmtosrc(2, ksrcwaq (isrc) + 1) = waqpar%iapnt(kk2) endif endif enddo end subroutine waq_prepare_src ! !------------------------------------------------------------------------------ !> Write WAQ pointer file. subroutine waq_wri_poi(filename) use m_alloc use wrwaq implicit none ! ! Global variables ! character(len=*), intent(in) :: filename !< Output filename. ! !! executable statements ------------------------------------------------------- ! ! Call the waq-poi file writer call wrwaqpoi(waqpar%ifrmto, waqpar%noq, filename, waq_format_ascii) end subroutine waq_wri_poi ! !------------------------------------------------------------------------------ !> Write WAQ len file. !! (contains dispersion lengths of computational cells) subroutine waq_wri_len(lnx, dx, acl, filename) use m_alloc use wrwaq implicit none ! ! Global variables ! integer, intent(in) :: lnx !< nr of flow links (internal + boundary) double precision, intent(in) :: dx(lnx) !< link length (m) double precision, intent(in) :: acl(lnx) !< left dx fraction, 0<=alfacl<=1 character(len=*), intent(in) :: filename !< Output filename. ! ! Local variables ! integer :: L, ip, kk integer, allocatable :: noqa(:) double precision, allocatable :: lenex(:,:) !< Length table: 'half' dx length from cell center to interface. !! lenex(1,:) = dx for left/1st cell to interface !! lenex(2,:) = dx for right/2nd cell to interface ! !! executable statements ------------------------------------------------------- ! call realloc(lenex, (/ 2, waqpar%noq /), keepExisting=.false., fill = 0d0) call realloc(noqa, waqpar%noql, keepExisting=.false., fill = 0) do L=1,lnx ip = waqpar%iqaggr(L) if (ip > 0) then ! MJ: TODO for now a simple average of the dispersion lengths, may be better to weight by wu (link initial width) lenex(1,ip) = lenex(1,ip) + dx(L)*acl(L) lenex(2,ip) = lenex(2,ip) + dx(L)*(1d0-acl(L)) noqa(ip) = noqa(ip) + 1 else if (ip < 0) then lenex(1,-ip) = lenex(1,-ip) + dx(L)*(1d0-acl(L)) lenex(2,-ip) = lenex(2,-ip) + dx(L)*acl(L) noqa(-ip) = noqa(-ip) + 1 end if end do do ip = 1, waqpar%noql if (waqpar%aggre == 1) then lenex(1,ip) = lenex(1,ip) / dble(noqa(ip)) lenex(2,ip) = lenex(2,ip) / dble(noqa(ip)) end if ! Copy lenghts to other layers do kk = 1, waqpar%kmxnxa - 1 lenex(1,ip + kk * waqpar%noql) = lenex(1,ip) lenex(2,ip + kk * waqpar%noql) = lenex(2,ip) end do end do ! dummy lengthes for sinks/sources do ip = waqpar%noq12 + 1, waqpar%noq12s lenex(1,ip) = 1d5 lenex(2,ip) = 1d5 enddo ! dummy lengthes in third direction for all layers (will be calculated by WAQ from volume and surface) do ip = waqpar%noq12s + 1, waqpar%noq lenex(1,ip) = 1d0 lenex(2,ip) = 1d0 end do ! Call the waq-len file writer call wrwaqlen(waqpar%noq, lenex, filename, waq_format_ascii) deallocate(lenex) end subroutine waq_wri_len ! !------------------------------------------------------------------------------ !> Write WAQ srf file. !! (contains horizontal surface areas of computational cells) subroutine waq_wri_srf(ndx2D, ndxi, ndx, ba, filename) use m_alloc use wrwaq implicit none ! ! Global variables ! integer, intent(in) :: ndx2D !< nr of 2D flow cells integer, intent(in) :: ndxi !< nr of internal flowcells (internal = 2D + 1D) integer, intent(in) :: ndx !< nr of flow nodes (internal + boundary) double precision, intent(in) :: ba(ndx) !< bottom area (m2), if < 0 use table in node type character(len=*), intent(in) :: filename !< Output filename. ! ! Local variables ! integer :: k, i ! !! executable statements ------------------------------------------------------- ! call realloc(waqpar%horsurf, waqpar%noseg, keepExisting=.false., fill=0d0) ! !! executable statements ------------------------------------------------------- ! ! AvD: TODO: What if ba(..) < 0. do k = 1, ndxi waqpar%horsurf(waqpar%iapnt(k)) = waqpar%horsurf(waqpar%iapnt(k)) + max(ba(k), 0d0) end do ! Copy to all layers if (waqpar%kmxnxa > 1) then do k = 1, waqpar%kmxnxa - 1 do i = 1, waqpar%nosegl waqpar%horsurf(k * waqpar%nosegl + i) = waqpar%horsurf(i) end do end do end if ! Call the waq-srf file writer call wrwaqsrf(waqpar%horsurf, waqpar%nosegl, waqpar%kmxnxa, filename, waq_format_ascii) end subroutine waq_wri_srf ! !------------------------------------------------------------------------------ !> Write WAQ atr file. !! (contains atributes of computational cells) subroutine waq_wri_atr(filename) use m_alloc use wrwaq implicit none ! ! Global variables ! character(len=*), intent(in) :: filename !< Output filename. ! !! executable statements ------------------------------------------------------- ! ! Call the waq-atr file writer call wrwaqatr(waqpar%nosegl, waqpar%kmxnxa, waqpar%kmk1, waqpar%kmk2, filename) end subroutine waq_wri_atr ! !------------------------------------------------------------------------------ !> Write WAQ vol file. !! (contains volumes of computational cells) subroutine waq_wri_vol(itim, filenamevol, lunvol) use m_flowgeom use m_flow use wrwaq implicit none ! ! Global variables ! integer, intent(in) :: itim !< time (seconds) since simulation start character(len=*), intent(in) :: filenamevol !< Output filename for volumes (only used when lunvol < 0). integer, intent(inout) :: lunvol !< File pointer for output vol-file (opened upon first call). ! ! Local variables ! integer :: i, k, kb, ktx, kk ! !! executable statements ------------------------------------------------------- ! waqpar%vol = 0d0 if (waqpar%aggre == 0 .and. waqpar%kmxnxa == 1) then do i = 1, ndxi waqpar%vol(i) = vol1(i) end do else if (waqpar%aggre == 0 .and. waqpar%aggrel == 0) then do k = 1, ndxi call getkbotktopmax(k,kb,ktx) do kk = kb, ktx waqpar%vol(waqpar%isaggr(kk)) = vol1(kk) end do end do else do k = 1, ndxi call getkbotktopmax(k,kb,ktx) do kk = kb, ktx waqpar%vol(waqpar%isaggr(kk)) = waqpar%vol(waqpar%isaggr(kk)) + vol1(kk) end do end do end if ! Call the waq-vol file writer call wrwaqbin(itim, waqpar%vol, waqpar%noseg, filenamevol, waq_format_ascii, lunvol) end subroutine waq_wri_vol ! !------------------------------------------------------------------------------ !> Write WAQ sal file. !! (contains salinity of computational cells) subroutine waq_wri_sal(itim, filenamesal, lunsal) use m_flowgeom use m_flow use wrwaq implicit none ! ! Global variables ! integer, intent(in) :: itim !< time (seconds) since simulation start character(len=*), intent(in) :: filenamesal !< Output filename for salinity (only used when lunsal < 0 and do_sal==1). integer, intent(inout) :: lunsal !< File pointer for output sal-file (opened upon first call, if do_sal==1). ! ! Local variables ! integer :: i, k, kb, ktx, kk ! !! executable statements ------------------------------------------------------- ! waqpar%sal = 0d0 if (waqpar%aggre == 0 .and. waqpar%kmxnxa == 1) then do i=1,ndxi if ( vol1(i) > 1d-25 ) then waqpar%sal(i) = sa1(i) end if end do else if (waqpar%aggre == 0 .and. waqpar%aggrel == 0) then do k = 1, ndxi call getkbotktopmax(k,kb,ktx) do kk = kb, ktx if ( vol1(kk) > 1d-25 ) then waqpar%sal(waqpar%isaggr(kk)) = sa1(kk) end if end do end do else ! Salinity is aggregated volume weighted do k = 1, ndxi call getkbotktopmax(k,kb,ktx) do kk = kb, ktx if ( vol1(kk) > 1d-25 ) then waqpar%sal(waqpar%isaggr(kk)) = waqpar%sal(waqpar%isaggr(kk)) + sa1(kk) * vol1(kk) end if end do end do do i = 1, waqpar%noseg if ( waqpar%vol(i) > 1d-25 ) then waqpar%sal(i) = waqpar%sal(i) / waqpar%vol(i) else waqpar%sal(i) = 0d0 end if end do end if ! Call the waq-vol file writer for salinity call wrwaqbin(itim, waqpar%sal, waqpar%noseg, filenamesal, waq_format_ascii, lunsal) end subroutine waq_wri_sal ! !------------------------------------------------------------------------------ !> Write WAQ tau file. !! (contains taus at the bottom of computational cells) subroutine waq_wri_tau(itim, filenametau, luntau) use m_flowgeom use m_flow use wrwaq implicit none ! ! Global variables ! integer, intent(in) :: itim !< time (seconds) since simulation start character(len=*), intent(in) :: filenametau !< Output filename for tau (only used when luntau < 0). integer, intent(inout) :: luntau !< File pointer for output tau-file (opened upon first call). ! ! Local variables ! integer :: i, iseg, k, kb, ktx, kk ! !! executable statements ------------------------------------------------------- ! waqpar%tau = 0d0 if (waqpar%aggre == 0 .and. waqpar%kmxnxa == 1) then do i = 1, ndxi waqpar%tau(i) = taus(i) end do else if (waqpar%aggre == 0 .and. waqpar%aggrel == 0) then do k = 1, ndxi call getkbotktopmax(k,kb,ktx) do kk = kb, ktx waqpar%tau(waqpar%isaggr(kk)) = taus(k) end do end do else ! Taus are aggregated horisontal surface weighted do k = 1, ndxi waqpar%tau(waqpar%isaggr(k)) = waqpar%tau(waqpar%isaggr(k)) + taus(k) * max(ba(k), 0d0) end do do i = 1, waqpar%nosegl if (waqpar%horsurf(i) > 1d-25) then waqpar%tau(i) = waqpar%tau(i) / waqpar%horsurf(i) end if end do do i = 1, waqpar%nosegl do k = 1, waqpar%kmxnxa - 1 waqpar%tau(i + k * waqpar%nosegl) = waqpar%tau(i) end do end do end if ! Call the waq-vol file writer for tau call wrwaqbin(itim, waqpar%tau, waqpar%noseg, filenametau, waq_format_ascii, luntau) end subroutine waq_wri_tau ! !------------------------------------------------------------------------------ !> Write WAQ vdf file. !! (contains vertical diffusion of computational cells) subroutine waq_wri_vdf(itim, filenamevdf, lunvdf) use m_flowgeom use m_flow use wrwaq implicit none ! ! Global variables ! integer, intent(in) :: itim !< time (seconds) since simulation start character(len=*), intent(in) :: filenamevdf !< Output filename for vertical diffusion (only used when there's more than one layer). integer, intent(inout) :: lunvdf !< File pointer for output tau-file (opened upon first call). ! ! Local variables ! integer :: i, k, kb, ktx, kk ! !! executable statements ------------------------------------------------------- ! waqpar%vdf = 0d0 if (waqpar%aggre == 0 .and. waqpar%aggrel == 0) then do k = 1, ndxi call getkbotktopmax(k,kb,ktx) do kk = kb, ktx if ( vol1(kk) > 1d-25 ) then waqpar%vdf(waqpar%isaggr(kk)) = dble(vicwws(kk)) end if end do end do else do k = 1, ndxi call getkbotktopmax(k,kb,ktx) do kk = kb, ktx if ( vol1(kk) > 1d-25 ) then waqpar%vdf(waqpar%isaggr(kk)) = waqpar%vdf(waqpar%isaggr(kk)) + dble(vicwws(kk)) * vol1(kk) end if end do end do do i = 1, waqpar%noseg if ( waqpar%vol(i) > 1d-25 ) then waqpar%vdf(i) = waqpar%vdf(i) / waqpar%vol(i) else waqpar%vdf(i) = 0d0 end if end do end if ! Call the waq-vol file writer for vertical diffusion call wrwaqbin(itim, waqpar%vdf, waqpar%noseg, filenamevdf, waq_format_ascii, lunvdf) end subroutine waq_wri_vdf ! !------------------------------------------------------------------------------ !> Write WAQ are file. !! (contains areas at flow interfaces) subroutine waq_wri_are(itim, filename, lun) use m_flowgeom use m_flow use wrwaq implicit none ! ! Global variables ! integer, intent(in) :: itim !< time (seconds) since simulation start character(len=*), intent(in) :: filename !< Output filename (only used when lun < 0). integer, intent(inout) :: lun !< File pointer for output are-file. ! ! Local variables ! integer :: i, L, Lb, Ltx, kxL, LL, ip ! !! executable statements ------------------------------------------------------- ! waqpar%area = 0d0 if (waqpar%aggre == 0 .and. waqpar%kmxnxa == 1) then do i = 1, lnx waqpar%area(i) = au(i) end do else if (waqpar%aggre == 0 .and. waqpar%aggrel == 0) then do L = 1, lnx call getLbotLtopmax(L,Lb,Ltx) do LL = Lb, Ltx waqpar%area(waqpar%iqaggr(LL)) = au(LL) end do end do else do L = 1, lnx call getLbotLtopmax(L,Lb,Ltx) if (abs(waqpar%iqaggr(Lb)) > 0) then do LL = Lb, Ltx ip = abs(waqpar%iqaggr(LL)) waqpar%area(ip) = waqpar%area(ip) + au(LL) end do end if end do end if ! dummy areas for sink/sources do i = waqpar%noq12 + 1, waqpar%noq12s waqpar%area(i) = 0.1D0 enddo ! Add area of the vertical exchanges if (waqpar%kmxnxa > 1) then do i = 1, waqpar%noseg - waqpar%nosegl waqpar%area(waqpar%noq12s + i) = waqpar%horsurf(i) end do end if ! Call the waq-flo file writer call wrwaqbin(itim, waqpar%area, waqpar%noq, filename, waq_format_ascii, lun) end subroutine waq_wri_are ! !------------------------------------------------------------------------------ !> Write WAQ flo file. !! (contains discharges at flow interfaces) subroutine waq_wri_flo(itim, ti_waq, filename, lun) use m_flowgeom use m_flow use wrwaq implicit none ! ! Global variables ! integer, intent(in) :: itim !< time (seconds) since simulation start integer, intent(in) :: ti_waq !< time step (seconds) between two waq outputs character(len=*), intent(in) :: filename !< Output filename (only used when lun < 0). integer, intent(inout) :: lun !< File pointer for output flo-file. ! ! Local variables ! integer :: isrc integer :: i, L, LL, Lb, kxL, Ltx, ip integer :: k, kk, kb, ktx integer, save :: itim_prev = -1 ! !! executable statements ------------------------------------------------------- ! waqpar%qag = 0d0 ! Average the accumulated discharges. if (waqpar%aggre == 0 .and. waqpar%kmxnxa == 1) then do i = 1, lnx waqpar%qag(i) = q1waq(i) / dble(ti_waq) end do else if (waqpar%aggre == 0 .and. waqpar%aggrel == 0) then do L = 1, lnx call getLbotLtopmax(L,Lb,Ltx) do LL = Lb, Ltx waqpar%qag(waqpar%iqaggr(LL)) = q1waq(LL) / dble(ti_waq) end do end do else do L = 1, lnx call getLbotLtopmax(L,Lb,Ltx) if (abs(waqpar%iqaggr(Lb)) > 0) then do LL = Lb, Ltx ip = abs(waqpar%iqaggr(LL)) if(waqpar%iqaggr(LL) > 0) then waqpar%qag(ip) = waqpar%qag(ip) + q1waq(LL) / dble(ti_waq) else waqpar%qag(ip) = waqpar%qag(ip) - q1waq(LL) / dble(ti_waq) end if end do end if end do end if ! Add sink/source dicharges !! TODO: write out discharges to a separe (ascii) file for additional wasteloads? if(waqpar%numsrcwaq > 0) then do isrc = 1, waqpar%numsrcwaq waqpar%qag(waqpar%noq12 + isrc) = qsrcwaq(isrc) / dble(ti_waq) enddo endif ! Add discharge of the vertical exchanges (note: upward flow direction in FM is reversed for WAQ!) if (waqpar%kmxnxa > 1) then if (waqpar%aggre == 0 .and. waqpar%aggrel == 0) then do k = 1, ndxi call getkbotktopmax(k,kb,ktx) do kk = kb, ktx - 1 if (waqpar%iqwaggr(kk)>0) then waqpar%qag(waqpar%iqwaggr(kk)) = - qwwaq(kk) / dble(ti_waq) end if end do end do else do k = 1, ndxi call getkbotktopmax(k,kb,ktx) do kk = kb, ktx if (waqpar%iqwaggr(kk)>0) then waqpar%qag(waqpar%iqwaggr(kk)) = waqpar%qag(waqpar%iqwaggr(kk)) - qwwaq(kk) / dble(ti_waq) end if end do end do end if end if ! Call the waq-flo file writer call wrwaqbin(itim, waqpar%qag, waqpar%noq, filename, waq_format_ascii, lun) end subroutine waq_wri_flo ! !------------------------------------------------------------------------------ !> Read an aggregation file (.dwq) into the global aggregation table. subroutine waq_read_dwq(ndxi, ndx, iapnt, filename) use unstruc_files implicit none ! ! Global variables ! integer, intent(in) :: ndxi !< Nr. of (internal) flow cells expected integer, intent(in) :: ndx !< Nr. of total flow cells (incl boundaries) integer, intent(out) :: iapnt(ndx) !< Mapping table flow cell -> waq cell (size ndx) character(len=*), intent(in) :: filename !< Input filename ! ! Local variables ! logical :: aggregate integer :: lundwq, istat, i, l1 integer :: headervals(5) character(len=4) :: fmt ! !! executable statements ------------------------------------------------------- ! aggregate = .true. l1 = len_trim(filename) if (l1 == 0) then aggregate = .false. else call oldfil(lundwq, trim(filename)) ! when the file does not exist the program stops(!) if ( lundwq == 0 ) then call mess(LEVEL_WARN, 'unable to open aggregation file '''//trim(filename)//'''.') aggregate = .false. else read (lundwq, *, iostat=istat) headervals if ( istat /= 0) then call mess(LEVEL_WARN, 'unable to read dimensions in aggregation file '''//trim(filename)//'''.') aggregate = .false. else if (headervals(3) /= ndxi) then call mess(LEVEL_WARN, 'mmax*nmax in dwqfile '''//trim(filename)//''' does not match ndxi: ', headervals(3), ndxi) aggregate = .false. else ! read the aggregation array read ( lundwq , *, iostat=istat) iapnt(1:ndxi) if ( istat /= 0) then call mess(LEVEL_WARN, 'unable to read aggregation pointer in '''//trim(filename)//'''.') aggregate = .false. end if end if end if close ( lundwq ) end if end if if (.not. aggregate) then ! if no aggregation could be read correctly, create default one to one aggregation call mess(LEVEL_WARN, 'no aggregation will be used '''//trim(filename)//'''.') do i = 1, ndxi iapnt(i) = i end do end if ! add pointer for boundary nodes do i = ndxi+1,ndx iapnt(i) = -(i-ndxi) enddo end subroutine waq_read_dwq ! !------------------------------------------------------------------------------ end module waq subroutine gettaus(typout) use m_flowgeom use m_flow use m_alloc implicit none integer, intent (in) :: typout !< type of setting, 1: set czs and taus, 2: just set czs: double precision :: taucurc !< local variable for taucurrent double precision :: czc !< local variable for chezy integer :: ierr !< Error code integer :: n !< Counter if (.not. allocated(czs) ) then call realloc(czs, ndxi, keepExisting = .false., fill = 0d0, stat = ierr) else if (size(czs) < ndxi) then call realloc(czs, ndxi, keepExisting = .false., fill = 0d0, stat = ierr) endif if (typout == 1) then if (.not. allocated(taus) ) then call realloc(taus, ndxi, keepExisting = .false., fill = 0d0, stat = ierr) else if (size(taus) < ndxi) then call realloc(taus, ndxi, keepExisting = .false., fill = 0d0, stat = ierr) endif endif do n = 1,ndxi call gettau(n,taucurc,czc) czs(n) = czc if (typout == 1) then taus(n) = taucurc endif enddo end subroutine gettaus subroutine gettau(n,taucurc,czc) use m_flowgeom use m_flow implicit none integer, intent(in) :: n !< Flow node number double precision, intent(out) :: taucurc !< Bed shear stress from current (taucurrent) double precision, intent(out) :: czc !< Chezy at flow node (taucurrent) ! ! Local variables ! integer :: L, LL, Lb, Lt, nn !< Local link counters integer :: kb, kt !< node counters double precision :: cf, cfn, cz, frcn, ar, wa, ust !< Local intermediate variables taucurc = 0d0 czc = 0d0 cfn = 0d0 wa = 0d0 ust = 0d0 do nn = 1,nd(n)%lnx LL = abs( nd(n)%ln(nn) ) frcn = frcu(LL) if (frcn > 0 .and. hu(LL) > 0) then call getcz(hu(LL), frcn, ifrcutp(LL), cz) cf = ag/(cz*cz) ar = au(LL)*dx(LL) wa = wa + ar ! area weigthed cfn = cfn + cf*ar if (kmx > 0) then ust = ust + ustb(LL)*ar endif endif enddo if (wa > 0) then cfn = cfn / wa endif if (cfn > 0) then czc = sqrt(ag/cfn) endif if (kmx == 0) then taucurc = rhomean*cfn*(ucx(n)*ucx(n) + ucy(n)*ucy(n)) else if (wa > 0) then ust = ust / wa endif taucurc = rhomean*ust*ust endif end subroutine gettau