!    Copyright (C) 2006 Imperial College London and others.
!    
!    Please see the AUTHORS file in the main source directory for a full list
!    of copyright holders.
!
!    Prof. C Pain
!    Applied Modelling and Computation Group
!    Department of Earth Science and Engineering
!    Imperial College London
!
!    amcgsoftware@imperial.ac.uk
!    
!    This library is free software; you can redistribute it and/or
!    modify it under the terms of the GNU Lesser General Public
!    License as published by the Free Software Foundation,
!    version 2.1 of the License.
!
!    This library 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
!    Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public
!    License along with this library; if not, write to the Free Software
!    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
!    USA

#include "fdebug.h"

module synthetic_bc

use fldebug
use spud
use global_parameters, only: dt, option_path_len, FIELD_NAME_LEN
use elements
use parallel_tools
use boundary_conditions
use transform_elements
use fetools
use fields
use state_module
use sparse_tools, only : csr_sparsity

implicit none

private
public synthetic_eddy_method, add_sem_bc, initialise_sem_memory, add_dfm_bc

  
type eddy
   !!< Store eddy info, for synthetic eddy method
   real, dimension(:),pointer:: xeddy
   real, dimension(:),pointer:: yeddy
   real, dimension(:),pointer:: zeddy

   integer, dimension(:),pointer:: eu
   integer, dimension(:),pointer:: ev
   integer, dimension(:),pointer:: ew
end type eddy


integer, save:: sem_bc_count=0
type(eddy),dimension(:),pointer,save:: eddies

contains

  subroutine synthetic_eddy_method(surface_field, surface_field1,  & 
       surface_field2, surface_field3, bc_position, bc_path_i, ns)

      ! declarations
      type(vector_field),intent(inout)     :: surface_field
      type(vector_field),intent(in)        :: surface_field1,surface_field2
      type(vector_field),intent(in)        :: surface_field3
      type(vector_field),intent(in)        :: bc_position
      type(element_type),pointer           :: x_shape,vel_shape
      real,allocatable,dimension(:)        :: detwei,vel_int
      real,allocatable,dimension(:,:),save :: uav
      real,allocatable,dimension(:),save   :: area
      real                                 :: xmin,xmax,ymin,ymax,zmin,zmax,tarea
      real                                 :: lxmin,lxmax,lymin,lymax,lzmin,lzmax
      real                                 :: xminb,xmaxb,yminb,ymaxb,zminb,zmaxb
      real                                 :: dxb,dyb,dzb
      real                                 :: rdmax,rts_x,rts_y,rts_z
      real                                 :: dx,dy,dz,vol
      real                                 :: dxp,dyp,dzp
      integer                              :: nsem,ns,nots,bcnod,i,j,ele,ratio
      real,allocatable,dimension(:,:)      :: cords
      real,dimension(2)                    :: tmpcord
      real,dimension(:),pointer            :: x,y,z
      real,dimension(:),pointer            :: lx,ly,lz
      real                                 :: uf,vf,wf,ff
      real,allocatable,dimension(:)        :: ufl,vfl,wfl
      real                                 :: a11,a21,a22,a31,a32,a33
      real                                 :: resuu,resvv,resww
      real,dimension(3)                    :: vl
      real                                 :: minum,minvm,minwm,maxum,maxvm,maxwm
      real                                 :: reuumn,reuumx,revvmn,revvmx,rewwmn,rewwmx
      real                                 :: uflmn,uflmx,vflmn,vflmx,wflmn,wflmx
      logical,save                         :: initz=.false.
      logical,allocatable,dimension(:),save:: initeddy
      character(len=OPTION_PATH_LEN)       :: bc_path_i

#ifdef HAVE_MPI
      INCLUDE 'mpif.h'
      INTEGER IERR
#endif

      ! time varying turbulent inlet conditions
      ! input:: mean flow, Reynolds stress tensor, turbulence length-scale
      ! synthetic eddy method (see Jarrin et al. 2006)
      ! can be used to generate turbulent inlet bcs for any given surface
      ! aligned with the cartesian coordinate system

      ! known issues::
      ! 1.there is a bug when calculating the inflow plane surface area in parallel-tarea
      ! 2.there is a bug when calculating the inflow bulk velocity in parallel-uav
      ! both because of the halo

      ! field names::
      ! surface_field  - inlet boundary condition (calculated here)
      ! surface_field1 - mean velocity
      ! surface_field2 - Re_ij (uu,vv,ww)
      ! surface_field3 - turbulence length scale
      ! bc_position    - inlet plane mesh

      ewrite(3,*) 'setting a turbulent inlet boundary condition',ns
      nsem = sem_bc_count
      ! get additional turbulence related options
      call get_option(trim(bc_path_i),nots)

      ! calculate min&max turbulence lengthscale
      lx => surface_field3%val(1,:)  
      ly => surface_field3%val(2,:)
      lz => surface_field3%val(3,:)

      lxmin = minval(lx); lxmax = maxval(lx)
      lymin = minval(ly); lymax = maxval(ly)
      lzmin = minval(lz); lzmax = maxval(lz)

      call allmin(lxmin); call allmax(lxmax)
      call allmin(lymin); call allmax(lymax)
      call allmin(lzmin); call allmax(lzmax)

      ewrite(3,*) 'number of eddies',nots
      ewrite(3,*) 'turbulence lengthscale min&max',lxmin,lxmax,':',lymin,lymax,':',lzmin,lzmax

      ! number of nodes on boundary surface
      bcnod=surface_field%mesh%nodes

      ! coordinates of nodes
      x => bc_position%val(1,:)  
      y => bc_position%val(2,:)
      z => bc_position%val(3,:)

      ! work out min & max of boundary surface and orientation
      xmin = minval(x); xmax = maxval(x)
      ymin = minval(y); ymax = maxval(y)
      zmin = minval(z); zmax = maxval(z)

      call allmin(xmin)
      call allmin(ymin)
      call allmin(zmin)

      call allmax(xmax)
      call allmax(ymax)
      call allmax(zmax)

      dxp = xmax-xmin
      dyp = ymax-ymin
      dzp = zmax-zmin

      ! generate bounding box
      if (dxp<dyp .and. dxp<dzp) then
         ewrite(3,*) 'using a yz plane'

         xminb=-lxmax+xmax
         xmaxb= lxmax+xmax
         yminb= ymin-lymax
         ymaxb= ymax+lymax
         zminb= zmin-lzmax
         zmaxb= zmax+lzmax

         rdmax=lxmax

      elseif (dyp<dxp .and. dyp<dzp) then
         ewrite(3,*) 'using a xz plane'

         xminb= xmin-lxmax
         xmaxb= xmax+lxmax
         yminb=-lymax+ymax
         ymaxb= lymax+ymax
         zminb= zmin-lzmax
         zmaxb= zmax+lzmax

         rdmax=lymax

      elseif (dzp<dxp .and. dzp<dyp) then
         ewrite(3,*) 'using a xy plane'

         xminb= xmin-lxmax
         xmaxb= xmax+lxmax
         yminb= ymin-lymax
         ymaxb= ymax+lymax
         zminb=-lzmax+zmax
         zmaxb= lzmax+zmax

         rdmax=lzmax
      end if

      dxb= xmaxb-xminb
      dyb= ymaxb-yminb
      dzb= zmaxb-zminb

      ! generate initial coordinades, signs of turbulent spots
      ! calculate convection velocity
      if (.not.initz) then
         allocate(initeddy(nsem))
         initeddy = .true.

         allocate(uav(nsem,3))
         allocate(area(nsem))
         uav=0.; area=0.

         initz=.true.
      endif

      if(initeddy(ns))then
         if(GetRank()==0) then
            allocate(cords(nots,3))
            call random_seed()
            call random_number(cords)
            
            do i=1,nots
               eddies(ns)%xeddy(i)=xminb+(xmaxb-xminb)*cords(i,1)
               eddies(ns)%yeddy(i)=yminb+(ymaxb-yminb)*cords(i,2)
               eddies(ns)%zeddy(i)=zminb+(zmaxb-zminb)*cords(i,3)

               eddies(ns)%eu(i)=esign()
               eddies(ns)%ev(i)=esign()
               eddies(ns)%ew(i)=esign()
            enddo
            deallocate(cords)
         end if

#ifdef HAVE_MPI
         if(IsParallel())then
            call mpi_bcast(eddies(ns)%xeddy(1:nots),nots,GetPREAL(),0,MPI_COMM_FEMTOOLS,ierr)
            call mpi_bcast(eddies(ns)%yeddy(1:nots),nots,GetPREAL(),0,MPI_COMM_FEMTOOLS,ierr)
            call mpi_bcast(eddies(ns)%zeddy(1:nots),nots,GetPREAL(),0,MPI_COMM_FEMTOOLS,ierr)

            call mpi_bcast(eddies(ns)%eu(1:nots),nots,mpi_integer,0,MPI_COMM_FEMTOOLS,ierr)
            call mpi_bcast(eddies(ns)%ev(1:nots),nots,mpi_integer,0,MPI_COMM_FEMTOOLS,ierr)
            call mpi_bcast(eddies(ns)%ew(1:nots),nots,mpi_integer,0,MPI_COMM_FEMTOOLS,ierr)
         endif
#endif

         x_shape   => ele_shape(bc_position,1)
         vel_shape => ele_shape(surface_field1,1)

         allocate(detwei(1:ele_ngi(bc_position,1)))
         allocate(vel_int(1:ele_loc(surface_field1,1)))

         tarea=0.
         do ele=1,element_count(surface_field1)
            call transform_to_physical(bc_position,ele,detwei=detwei)
            vel_int=shape_rhs(vel_shape, detwei)
            tarea=tarea+sum(vel_int)
            do i=1,surface_field1%dim
               uav(ns,i)=uav(ns,i)+sum(ele_val(surface_field1,i,ele)*vel_int)
            enddo
         enddo

         deallocate(detwei,vel_int)

         call allsum(tarea)
         call allsum(uav(ns,:))

         do i=1,surface_field1%dim
            uav(ns,i)=uav(ns,i)/tarea
         enddo

         area(ns)=tarea
         initeddy(ns)=.false.
      endif
      
      ewrite(3,*) 'surface info'
      ewrite(3,*) 'x min&max',xmin,xmax,dxp
      ewrite(3,*) 'y min&max',ymin,ymax,dyp
      ewrite(3,*) 'z min&max',zmin,zmax,dzp
      ewrite(3,*) 'surface area',area(ns)
      ewrite(3,*) 'number of nodes',bcnod
      ewrite(3,*) 'convection velocity:',uav(ns,1),uav(ns,2),uav(ns,3)

      ! bounding box volume
      vol=area(ns)*2*rdmax

      ewrite(3,*) 'bounding box info'
      ewrite(3,*) 'x min&max',xminb,xmaxb
      ewrite(3,*) 'y min&max',yminb,ymaxb
      ewrite(3,*) 'z min&max',zminb,zmaxb
      ewrite(3,*) 'box volume',vol

      if(GetRank()==0) then
         ! convect turbulent spots
         do i=1,nots
            eddies(ns)%xeddy(i)=eddies(ns)%xeddy(i)+uav(ns,1)*dt
            eddies(ns)%yeddy(i)=eddies(ns)%yeddy(i)+uav(ns,2)*dt
            eddies(ns)%zeddy(i)=eddies(ns)%zeddy(i)+uav(ns,3)*dt
         enddo

         ! regenerate
         do i=1,nots

            if(eddies(ns)%xeddy(i)>xmaxb)then
               call random_number(tmpcord)

               ratio= int((abs(eddies(ns)%xeddy(i)-xminb))/dxb)

               eddies(ns)%xeddy(i)=eddies(ns)%xeddy(i)-ratio*dxb
               eddies(ns)%yeddy(i)=yminb+(ymaxb-yminb)*tmpcord(1)
               eddies(ns)%zeddy(i)=zminb+(zmaxb-zminb)*tmpcord(2)

               eddies(ns)%eu(i)=esign()
               eddies(ns)%ev(i)=esign()
               eddies(ns)%ew(i)=esign()

            elseif(eddies(ns)%xeddy(i)<xminb)then
               call random_number(tmpcord)

               ratio= int((abs(eddies(ns)%xeddy(i)-xmaxb))/dxb)

               eddies(ns)%xeddy(i)=eddies(ns)%xeddy(i)+ratio*dxb
               eddies(ns)%yeddy(i)=yminb+(ymaxb-yminb)*tmpcord(1)
               eddies(ns)%zeddy(i)=zminb+(zmaxb-zminb)*tmpcord(2)

               eddies(ns)%eu(i)=esign()
               eddies(ns)%ev(i)=esign()
               eddies(ns)%ew(i)=esign()

            elseif(eddies(ns)%yeddy(i)>ymaxb)then
               call random_number(tmpcord)

               ratio= int((abs(eddies(ns)%yeddy(i)-yminb))/dyb)

               eddies(ns)%xeddy(i)=xminb+(xmaxb-xminb)*tmpcord(1)
               eddies(ns)%yeddy(i)=eddies(ns)%yeddy(i)-ratio*dyb
               eddies(ns)%zeddy(i)=zminb+(xmaxb-xminb)*tmpcord(2)

               eddies(ns)%eu(i)=esign()
               eddies(ns)%ev(i)=esign()
               eddies(ns)%ew(i)=esign()

            elseif(eddies(ns)%yeddy(i)<yminb)then
               call random_number(tmpcord)

               ratio= int((abs(eddies(ns)%yeddy(i)-ymaxb))/dyb)

               eddies(ns)%xeddy(i)=xmaxb+(xmaxb-xminb)*tmpcord(1)
               eddies(ns)%yeddy(i)=eddies(ns)%yeddy(i)+ratio*dyb
               eddies(ns)%zeddy(i)=zminb+(zmaxb-zminb)*tmpcord(2)

               eddies(ns)%eu(i)=esign()
               eddies(ns)%ev(i)=esign()
               eddies(ns)%ew(i)=esign()

            elseif(eddies(ns)%zeddy(i)>zmaxb)then
               call random_number(tmpcord)

               ratio= int((abs(eddies(ns)%zeddy(i)-zminb))/dzb)

               eddies(ns)%xeddy(i)=xminb+(xmaxb-xminb)*tmpcord(1)
               eddies(ns)%yeddy(i)=yminb+(ymaxb-yminb)*tmpcord(2)
               eddies(ns)%zeddy(i)=eddies(ns)%zeddy(i)-ratio*dzb

               eddies(ns)%eu(i)=esign()
               eddies(ns)%ev(i)=esign()
               eddies(ns)%ew(i)=esign()

            elseif(eddies(ns)%zeddy(i)<zminb)then
               call random_number(tmpcord)

               ratio= int((abs(eddies(ns)%zeddy(i)-zmaxb))/dzb)

               eddies(ns)%xeddy(i)=xmaxb+(xmaxb-xminb)*tmpcord(1)
               eddies(ns)%yeddy(i)=yminb+(ymaxb-yminb)*tmpcord(2)
               eddies(ns)%zeddy(i)=eddies(ns)%zeddy(i)+ratio*dzb
               
               eddies(ns)%eu(i)=esign()
               eddies(ns)%ev(i)=esign()
               eddies(ns)%ew(i)=esign()
            endif
         enddo
      end if

#ifdef HAVE_MPI
      if(IsParallel())then
         call mpi_bcast(eddies(ns)%xeddy(1:nots),nots,GetPREAL(),0,MPI_COMM_FEMTOOLS,ierr)
         call mpi_bcast(eddies(ns)%yeddy(1:nots),nots,GetPREAL(),0,MPI_COMM_FEMTOOLS,ierr)
         call mpi_bcast(eddies(ns)%zeddy(1:nots),nots,GetPREAL(),0,MPI_COMM_FEMTOOLS,ierr)

         call mpi_bcast(eddies(ns)%eu(1:nots),nots,mpi_integer,0,MPI_COMM_FEMTOOLS,ierr)
         call mpi_bcast(eddies(ns)%ev(1:nots),nots,mpi_integer,0,MPI_COMM_FEMTOOLS,ierr)
         call mpi_bcast(eddies(ns)%ew(1:nots),nots,mpi_integer,0,MPI_COMM_FEMTOOLS,ierr)        
      endif
#endif 

      !ewrite(3,*) 'eddy coordinates & signs'
      !do i=1,nots
      !   ewrite(3,*) i,                                                    &
      !        eddies(ns)%xeddy(i),eddies(ns)%yeddy(i),eddies(ns)%zeddy(i), & 
      !        eddies(ns)%eu(i),eddies(ns)%ev(i),eddies(ns)%ew(i)
      !enddo
      
      ! calculate fluctuating component
      allocate(ufl(bcnod))
      allocate(vfl(bcnod))
      allocate(wfl(bcnod))

      ufl=0.; vfl=0.; wfl=0.

      do i=1,bcnod

         uf=0.; vf=0.; wf=0.
         rts_x=node_val(surface_field3,1,i)
         rts_y=node_val(surface_field3,2,i)
         rts_z=node_val(surface_field3,3,i)

         do j=1,nots

            dx=abs(x(i)-eddies(ns)%xeddy(j))
            dy=abs(y(i)-eddies(ns)%yeddy(j))
            dz=abs(z(i)-eddies(ns)%zeddy(j))

            if (dx<rts_x .and. dy<rts_y .and. dz<rts_z) then

               ff=(1.-dx/rts_x)*(1.-dy/rts_y)*(1.-dz/rts_z)
               ff=ff/((sqrt(3./2.*rts_x))*(sqrt(3./2.*rts_y))*(sqrt(3./2.*rts_z)))

               uf=uf+eddies(ns)%eu(j)*ff
               vf=vf+eddies(ns)%ev(j)*ff
               wf=wf+eddies(ns)%ew(j)*ff

            end if
         end do

         uf=uf*sqrt(vol/nots)
         vf=vf*sqrt(vol/nots)
         wf=wf*sqrt(vol/nots)

         ! Cholesky decomposition of the Re_ij
         resuu=node_val(surface_field2,1,i)
         resvv=node_val(surface_field2,2,i)
         resww=node_val(surface_field2,3,i)

         if (resuu<0.)then
            ewrite(3,*) 'WARNING: there is a negative value in Re_uu'
            resuu=1.e-5
         end if
         if (resvv<0.)then
            ewrite(3,*) 'WARNING: there is a negative value in Re_vv'
            resvv=1.e-5
         end if
         if (resww<0.)then
            ewrite(3,*) 'WARNING: there is a negative value in Re_ww'
            resww=1.e-5
         end if
         
         a11 = sqrt(resuu)
         a21 = 0.
         a22 = sqrt(resvv)
         a31 = 0.
         a32 = 0.
         a33 = sqrt(resww)

         ufl(i) = a11* uf
         vfl(i) = a21*uf + a22*vf
         wfl(i) = a31*uf + a32*vf + a33*wf
      end do

      ! get min&max of mean velocities
      minum=minval(surface_field1%val(1,:)); maxum=maxval(surface_field1%val(1,:))
      minvm=minval(surface_field1%val(2,:)); maxvm=maxval(surface_field1%val(2,:))
      minwm=minval(surface_field1%val(3,:)); maxwm=maxval(surface_field1%val(3,:))

      call allmin(minum); call allmax(maxum)
      call allmin(minvm); call allmax(maxvm)
      call allmin(minwm); call allmax(maxwm)

      ! calculate final velocity signal
      do i = 1, bcnod
         vl(1) = node_val(surface_field1,1,i)+ufl(i)
         vl(2) = node_val(surface_field1,2,i)+vfl(i)
         vl(3) = node_val(surface_field1,3,i)+wfl(i)
         call set(surface_field,i,vl)
      end do

      ! calculate min&max of Re_ij
      reuumn = minval(surface_field2%val(1,:)); reuumx=maxval(surface_field2%val(1,:))
      revvmn = minval(surface_field2%val(2,:)); revvmx=maxval(surface_field2%val(2,:))
      rewwmn = minval(surface_field2%val(3,:)); rewwmx=maxval(surface_field2%val(3,:))

      call allmin(reuumn); call allmax(reuumx)
      call allmin(revvmn); call allmax(revvmx)
      call allmin(rewwmn); call allmax(rewwmx)

      ! calculate min&max of fluctuating component
      uflmn=minval(ufl); uflmx=maxval(ufl)
      vflmn=minval(vfl); vflmx=maxval(vfl)
      wflmn=minval(wfl); wflmx=maxval(wfl)

      call allmin(uflmn); call allmax(uflmx)
      call allmin(vflmn); call allmax(vflmx)
      call allmin(wflmn); call allmax(wflmx)

      ewrite(3,*) 'inlet condition diagnostics'
      ewrite(3,*) 'Re stresses min&max'
      ewrite(3,*) 'Re_uu',reuumn,':',reuumx
      ewrite(3,*) 'Re_vv',revvmn,':',revvmx
      ewrite(3,*) 'Re_ww',rewwmn,':',rewwmx
      ewrite(3,*) 'mean:fluctuating min&max'
      ewrite(3,*) 'u',minum,maxum,':',uflmn,uflmx
      ewrite(3,*) 'v',minvm,maxvm,':',vflmn,vflmx
      ewrite(3,*) 'w',minwm,maxwm,':',wflmn,wflmx

      deallocate(ufl,vfl,wfl)

      ewrite(3,*) 'leaving turbulent inlet boundary conditions routine'
    end subroutine synthetic_eddy_method

    !----------------------------------------------------------

    function esign()

      integer :: esign
      real    :: i

      call random_number(i)

      if (i<0.5) then
        esign = -1        
      else
        esign = +1
      endif

      return 
     end function esign
     
     !----------------------------------------------------------
     
     subroutine initialise_sem_memory(ns,nots)

       logical,save::initialise_memory=.false.
       logical,allocatable,dimension(:),save::initeddymem
       integer:: ns, nots, nsem

       nsem=sem_bc_count

       if(.not.initialise_memory)then
          allocate(initeddymem(nsem))
          initeddymem=.true.
          allocate(eddies(nsem))
          initialise_memory=.true.
       endif
       
       if(initeddymem(ns))then
          allocate(eddies(ns)%xeddy(nots));allocate (eddies(ns)%yeddy(nots));allocate(eddies(ns)%zeddy(nots))
          allocate(eddies(ns)%eu(nots))   ;allocate(eddies(ns)%ev(nots))    ;allocate(eddies(ns)%ew(nots))
          initeddymem(ns)=.false.
       endif
       
     end subroutine initialise_sem_memory
     
     !----------------------------------------------------------
     
     subroutine add_sem_bc(have_sem_bc)

       logical:: have_sem_bc

       if (have_sem_bc) then
          sem_bc_count=sem_bc_count+1
       end if
       
     end subroutine add_sem_bc
     
  subroutine add_dfm_bc(field, surface_field, surface_element_list, bc_position, surface_mesh, bc_name, bc_path, timestep)
  
    implicit none
    type(vector_field), intent(in) :: field
    type(vector_field), intent(inout) :: surface_field
    type(vector_field), intent(in) :: bc_position
    type(mesh_type), intent(inout) :: surface_mesh
    real, intent(in) :: timestep
    character(len=FIELD_NAME_LEN), intent(in) :: bc_path, bc_name
    integer, dimension(:), pointer, intent(in):: surface_element_list
    
    integer :: nlevel, rank, nneigh, node, neigh, ele, i, j, k, l, m, n, e
    integer, dimension(100) :: neighbours_list
    real :: dist, ui_old, bi, ai, ff, sum_bi, intt
    real, dimension(bc_position%dim) :: intl, reys, ui
    real, dimension(bc_position%dim) :: X
    integer, dimension(:), pointer :: nodes
    
    type(scalar_field) :: neighbour_flags, node_flags, filter_coeff
    type(vector_field), pointer :: length_scale, time_scale, reynolds_stress, fluctuations, old_fluctuations
    type(vector_field) :: rand
    type(csr_sparsity), pointer :: nelist
    real, parameter :: pi=3.14159265359
    
    ewrite(1,*) 'In add_dfm_bc'

    call get_option(trim(bc_path)//"/digital_filter_turbulence/number_of_layers", nlevel)
    
    call allocate (filter_coeff, surface_mesh, name="FilterCoefficients")
    call allocate (rand, surface_field%dim, surface_mesh, name="Random")
    call allocate (neighbour_flags, surface_mesh, name="NeighbourFlags")
    call allocate (node_flags, surface_mesh, name="NodeFlags")
    
    fluctuations => extract_surface_field(field,bc_name,"FluctuationField")
    old_fluctuations => extract_surface_field(field,bc_name,"OldFluctuationField")
    length_scale => extract_surface_field(field,bc_name,"TurbulenceLengthScale")
    time_scale => extract_surface_field(field,bc_name,"TurbulenceTimeScale")
    reynolds_stress => extract_surface_field(field,bc_name,"ReynoldsStress")
        
   ! Generate pseudo-random numbers
    call random_seed()
    do node = 1, node_count(surface_field)
       call random_number(X)
       X = 2.*X-1.
       call set(rand,node,X)
    enddo
    
    ! Loop over elements
    call zero(node_flags)
    nelist => extract_nelist(surface_mesh)
    do ele = 1, ele_count(surface_field)
       nodes => ele_nodes(surface_mesh,ele)
       
       do n = 1, size(nodes)
         node = nodes(n)
	 if (node_val(node_flags,node) == 0.) then
       
         ! First, calculate filter coefficients
         intl = node_val(length_scale,node)
         intt = node_val(time_scale,1,node)
         reys = node_val(reynolds_stress,node)	 
	 
         ! Find elements within filter width
	 call zero(neighbour_flags)
         call zero(filter_coeff)
	 neighbours_list=0
	 nneigh=0
	 rank=0
	 
	 call find_neighbours (surface_mesh, node, neighbour_flags, neighbours_list, nneigh, rank)
	 
	 ! Calculate filter coefficients
         sum_bi= 0.
	 do i = 1, nneigh
	   neigh = neighbours_list(i)
	   
	   bi = 1.
	   do k = 1, surface_field%dim
	     dist=node_val(bc_position,k,neigh)-node_val(bc_position,k,node)
	     bi=bi*exp(-pi/2.*dist**2./intl(k)**2.)
	   end do
	   sum_bi=sum_bi+bi*bi
	   
	   call set(filter_coeff,neigh,bi)
	 enddo
	 call scale(filter_coeff, 1./sqrt(sum_bi))
	 
         ! Calculate rescaled fluctuations
	 ui = 0.
         do i = 1, nneigh
           neigh = neighbours_list(i) 
           do j = 1, surface_field%dim
	     ui(j) = ui(j) + node_val(filter_coeff,neigh)*node_val(rand,j,neigh)
           enddo
         enddo
         ui = ui*sqrt(reys)
       
         ! Time correlation
         ai = exp(-pi/2.*timestep**2./intt**2.)
         ff = 1./sqrt(1. + ai*ai)
         do j = 1, surface_field%dim
	   ui_old = node_val(old_fluctuations,j,node)
	   ui(j) = ff*ui(j) + ai*ff*ui_old
         enddo
	 
	 call set(fluctuations,node,ui)
	 call set(node_flags,node,1.)
	 
	 endif
       end do
       
    enddo
      
    call addto(surface_field,fluctuations)
    call set(old_fluctuations,fluctuations)
    
    call deallocate (filter_coeff)
    call deallocate (rand)
    call deallocate (neighbour_flags)
    call deallocate (node_flags)

  contains
    
    recursive subroutine find_neighbours (surface_mesh, node, neighbour_flags, neighbours_list, nneigh, rank)
      
      integer :: node, nneigh, rank
      type(mesh_type) :: surface_mesh
      type(scalar_field) :: neighbour_flags  
      integer, dimension(:) :: neighbours_list
         
      integer :: l, m
      integer, dimension(:), pointer :: ele_neighbours, neighbours
           
      ! Loop over neighbours
      rank = rank+1
      ele_neighbours => node_neigh(surface_mesh,node)
      do l = 1, size(ele_neighbours)

         neighbours => ele_nodes(surface_mesh,ele_neighbours(l))
         do m = 1, size(neighbours)
    	   if (node_val(neighbour_flags,neighbours(m)) == 0.) then
    	     nneigh=nneigh+1
    	     neighbours_list(nneigh)=neighbours(m)
    	     call set(neighbour_flags, neighbours(m), 1.)	     
    	   endif
	   
	   if (rank < nlevel) &
	      call find_neighbours (surface_mesh, neighbours(m), neighbour_flags, neighbours_list, nneigh, rank)
         enddo
       
      enddo
    
    end subroutine find_neighbours

  end subroutine add_dfm_bc
     
   end module synthetic_bc
   
