subroutine wrwaqfil  ( out_prefix , np            , ns           , ne      , nolay  ,        &
     &                       dp         , area_cell     , isidenode    , distj   , is     ,        &
     &                       x          , y             , xcj          , ycj     , nm     ,        &
     &                       dist_dw    , nopenbndnodes , openbndnodes )

!     Deltares Software Centre

!>\file
!>        Writes static Delwaq information out of Selfe data

!     Created            : January 2013  By Leo Postma

!     Modified           :

!     Subroutines called : wrwaqhyd

!     Logical units      :
!
!                          536 = unit Delwaq  dps, depths below ref plane in nodes
!                          537 = unit Delwaq  srf, horizontalk surface areas around nodes

      implicit none

!     Parameters         :

!     kind           function         name                 Descriptipon

      character( *), intent(in   ) :: out_prefix       !<  Output files prefix string
      integer  ( 4), intent(in   ) :: np               !<  Number of nodes (points)
      integer  ( 4), intent(in   ) :: ns               !<  Number of sides (edges)
      integer  ( 4), intent(in   ) :: ne               !<  Number of elements
      integer  ( 4), intent(in   ) :: nolay            !<  Number of layers
      real     ( 4), intent(in   ) :: dp       (np)    !<  Depth values at the points
      real     ( 4), intent(in   ) :: area_cell(np)    !<  Horizontal surface areas
      integer  ( 4), intent(in   ) :: isidenode(ns,2)  !<  Number of sides (edges)
      real     ( 4), intent(in   ) :: distj    (ns)    !<  Distances between the nodes of sides
      integer  ( 4), intent(in   ) :: is       (ns,2)  !<  Two elements sharing a side
      real     ( 4), intent(in   ) :: x        (np)    !<  x of the nodes
      real     ( 4), intent(in   ) :: y        (np)    !<  y of the nodes
      real     ( 4), intent(in   ) :: xcj      (ns)    !<  x of the side center
      real     ( 4), intent(in   ) :: ycj      (ns)    !<  y of the side center
      integer  ( 4), intent(in   ) :: nm       (ne,3)  !<  nodes of an element
      real     ( 4), intent(  out) :: dist_dw  (ns)    !<  distance perpendicular to the side
      integer  ( 4), intent(in   ) :: nopenbndnodes  !<  Number of open boundary nodes
      integer  ( 4), intent(in   ) :: openbndnodes(nopenbndnodes)  !< Their nodenumbers

!     Locals

      character(128)  outstrng                  !  Help string for output
      integer  ( 4 )  i, k                      !  Loop and help variables
      real     ( 4 ), allocatable :: xel(:)     !  x of element center
      real     ( 4 ), allocatable :: yel(:)     !  y of element center
      real     ( 4 )  x1, y1, x2, y2            !  help variables

      outstrng = out_prefix
      i = len_trim(outstrng)
      outstrng(i:i+3) = '.dps'
      open  ( 536, file=outstrng, form='unformatted', convert='big_endian' )
      outstrng = out_prefix
      i = len_trim(outstrng)
      outstrng(i:i+3) = '.srf'
      open  ( 537, file=outstrng, form='unformatted', convert='big_endian' )
      outstrng = out_prefix
      i = len_trim(outstrng)
      outstrng(i:i+3) = '.poi'
      open  ( 538, file=outstrng, form='binary' )
      outstrng = out_prefix
      i = len_trim(outstrng)
      outstrng(i:i+3) = '.len'
      open  ( 539, file=outstrng, form='binary' )
      outstrng = out_prefix
      i = len_trim(outstrng)
      outstrng(i:i+3) = '.vol'
      open  ( 540, file=outstrng, form='binary' )
      outstrng = out_prefix
      i = len_trim(outstrng)
      outstrng(i:i+3) = '.flo'
      open  ( 541, file=outstrng, form='binary' )
      outstrng = out_prefix
      i = len_trim(outstrng)
      outstrng(i:i+3) = '.are'
      open  ( 542, file=outstrng, form='binary' )
      outstrng = out_prefix
      i = len_trim(outstrng)
      outstrng(i:i+5) = '-u.vel'
      open  ( 543, file=outstrng, form='binary' )
      outstrng = out_prefix
      i = len_trim(outstrng)
      outstrng(i:i+5) = '-v.vel'
      open  ( 544, file=outstrng, form='binary' )
      outstrng = out_prefix
      i = len_trim(outstrng)
      outstrng(i:i+3) = '.kmk'
      open  ( 545, file=outstrng, form='binary' )

      write ( 536 ) np, 0, np, np, np, 0
      write ( 536 ) dp
      write ( 537 ) np, 0, np, np, np, 0
      write ( 537 ) area_cell
      do k = 0, nolay-1
         write ( 538 )  ( ( -i -k*nopenbndnodes, openbndnodes(i) + k   *np, 0, 0 ) , i=1,nopenbndnodes )
         write ( 538 )  ( ( isidenode(i,1)+k*np,   isidenode(i,2)+ k   *np, 0, 0 ) , i=1,ns            )
      enddo
      write ( 538 ) ( ( (  i+k*np, i+(k+1)*np, 0, 0 ) , i = 1,np ), k=0,nolay-2 )
      write ( 539 ) 0.0E+00
      do k = 0, nolay-1
         write ( 539 ) ( (   100.0     ,    100.0    ) , i=1,nopenbndnodes )
         write ( 539 ) ( ( 0.5*distj(i),0.5*distj(i) ) , i=1,ns            )
      enddo
      write ( 539 ) ( ( ( 0.0E+00       , 0.0E+00      ) ,i=1,np )   , k=0, nolay-2   )

      allocate ( xel(ne), yel(ne) )
      do i = 1, ne
         xel(i) = ( x(nm(i,1))+x(nm(i,2))+x(nm(i,3)) ) / 3.0
         yel(i) = ( y(nm(i,1))+y(nm(i,2))+y(nm(i,3)) ) / 3.0
      enddo
      dist_dw = 0.0
      do i = 1, ns
         x1 = xcj(i) ; y1 = ycj(i)
         x2 = xcj(i) ; y2 = ycj(i)
         if ( is(i,1) .gt. 0 ) then
            x1 = xel(is(i,1))
            y1 = yel(is(i,1))
         endif
         if ( is(i,2) .gt. 0 ) then
            x2 = xel(is(i,2))
            y2 = yel(is(i,2))
         endif
         dist_dw(i) = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))
      enddo

      return
      end