! $Id$
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of special_dummies.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lspecial = .true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 0
!
!***************************************************************
!
module Special
!
  use Cparam
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
!
  implicit none
!
  include '../special.h'
!
  !
  ! Slice precalculation buffers
  !
!
  character(len=24) :: initspecial='nothing'
  logical :: first_time=.true.
!
!  Variables to be used when getting timevarying inlet from file
!
  real, allocatable, dimension(:,:,:,:) :: f_in
  real, allocatable, dimension(:) :: x_in
  real, allocatable, dimension(:) :: y_in
  real, allocatable, dimension(:) :: z_in
  character :: prec_in
  real :: t_in,dx_in,dy_in,dz_in
  integer :: mx_in,my_in,mz_in,nv_in
  integer :: l1_in, nx_in, ny_in, nz_in
  integer :: mvar_in,maux_in,mglobal_in
  integer :: nghost_in
  integer :: m1_in
  integer :: n1_in
  integer :: l2_in
  integer :: m2_in
  integer :: n2_in
  real :: Lx_in
  real :: Ly_in
  real :: Lz_in
  character (len=fnlen) :: turbfile
  character(len=fnlen) :: turb_inlet_dir=''
  logical :: proc_at_inlet
  integer :: ipx_in, ipy_in, ipz_in, iproc_in, nprocx_in, nprocy_in, nprocz_in
  character (len=fnlen) :: directory_in
  real, dimension(2) :: radius=(/0.0182,0.0364/)
  real, dimension(2) :: momentum_thickness=(/0.014,0.0182/)
  real, dimension(2) :: jet_center=(/0.,0./)
  real :: u_t=5.,velocity_ratio=3.3
!
!
!
!
!!  character, len(50) :: initcustom
!
! input parameters
  namelist /jet_init_pars/ &
      initspecial,turb_inlet_dir,u_t,velocity_ratio,radius,momentum_thickness,&
      jet_center
  ! run parameters
  namelist /jet_run_pars/  &
       turb_inlet_dir
!
  contains
!
!***********************************************************************
    subroutine register_special()
!
!  Configure pre-initialised (i.e. before parameter read) variables
!  which should be know to be able to evaluate
!
!  6-oct-03/tony: coded
!
      if (lroot) call svn_id( &
           "$Id$")
!
    endsubroutine register_special
!***********************************************************************
    subroutine initialize_special(f)
!
!  called by run.f90 after reading parameters, but before the time loop
!
!  19-jan-10/nils: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      call keep_compiler_quiet(f)
!
    endsubroutine initialize_special
!***********************************************************************
    subroutine init_special(f)
!
!  initialise special condition; called from start.f90
!  06-oct-2003/tony: coded
!
      use Mpicomm
      use Sub
      use Initcond
!
      real, dimension (mx,my,mz,mfarray) :: f
      integer :: i,j,jjj,kkk
      real, dimension(3) :: velo,tmp
      real :: radius_mean,An,rad
      logical :: non_zero_transveral_velo
!
      intent(inout) :: f
      !
      ! Select case
      !
      select case (initspecial)
      case ('nothing'); if (lroot) print*,'init_special: nothing'
      case ('coaxial_jet')
        if (lroot) print*,'init_special: coaxial_jet'
        velo(1)=u_t
        velo(2)=velo(1)*velocity_ratio
        velo(3)=0.04*velo(2)
        radius_mean=(radius(1)+radius(2))/2.
!
! Set velocity profiles
!
        do jjj=1,ny
          do kkk=1,nz
            rad=sqrt(&
                (y(jjj+m1-1)-jet_center(1))**2+&
                (z(kkk+n1-1)-jet_center(2))**2)
            ! Add mean velocity profile
            if (rad < radius_mean) then
              f(:,jjj+m1-1,kkk+n1-1,1)&
                  =(velo(1)+velo(2))/2&
                  +(velo(2)-velo(1))/2*tanh((rad-radius(1))/&
                  (2*momentum_thickness(1)))
            else
              f(:,jjj+m1-1,kkk+n1-1,1)&
                  =(velo(2)+velo(3))/2&
                  +(velo(3)-velo(2))/2*tanh((rad-radius(2))/&
                  (2*momentum_thickness(2)))
            endif
          enddo
        enddo
        f(:,:,:,iuy:iuz)=0
!
        case ('single_jet')
          if (lroot) print*,'init_special: single_jet'
          velo(1)=u_t
          velo(2)=velo(1)/velocity_ratio
!
! Set velocity profiles
!
          do jjj=1,ny
            do kkk=1,nz
              rad=sqrt(&
                  (y(jjj+m1-1)-jet_center(1))**2+&
                  (z(kkk+n1-1)-jet_center(1))**2)
              !Add velocity profile
              f(:,jjj+m1-1,kkk+n1-1,1)&
                  =velo(1)*(1-tanh((rad-radius(1))/&
                  momentum_thickness(1)))*0.5+velo(2)
            enddo
          enddo
          f(:,:,:,iuy:iuz)=0
!
        case ('single_laminar_wall_jet')
          if (lroot) print*,'init_special: single_laminar_wall_jet'
!
! Set velocity profiles
!
!!$          do jjj=1,ny
!!$            do kkk=1,nz
!!$              rad=sqrt(&
!!$                  (y(jjj+m1-1)-jet_center(1))**2+&
!!$                  (z(kkk+n1-1)-jet_center(1))**2)
!!$              !Add velocity profile
!!$              if (rad < radius(1)) then
!!$                f(:,jjj+m1-1,kkk+n1-1,1)=u_t*(1-(rad/radius(1))**2)
!!$              else
!!$                f(:,jjj+m1-1,kkk+n1-1,1)=0.
!!$              endif
!!$            enddo
!!$          enddo
!
!
!
          do jjj=1,my
            do kkk=1,mz
              rad=sqrt(&
                  (y(jjj)-jet_center(1))**2+&
                  (z(kkk)-jet_center(1))**2)
              !Add velocity profile
              if (rad < radius(1)) then
                f(:,jjj,kkk,1)=u_t*(1-(rad/radius(1))**2)
              else
                f(:,jjj,kkk,1)=0.
              endif
            enddo
          enddo
!
          f(:,:,:,iuy:iuz)=0
!
      case default
        !
        !  Catch unknown values
        !
        if (lroot) print*,'init_special: No such value for initspecial: ', &
             trim(initspecial)
        call stop_it("")
      endselect
      !
    endsubroutine init_special
!***********************************************************************
    subroutine pencil_criteria_special()
!
!  All pencils that this special module depends on are specified here.
!
!  18-07-06/tony: coded
!
    endsubroutine pencil_criteria_special
!***********************************************************************
    subroutine pencil_interdep_special(lpencil_in)
!
!  Interdependency among pencils provided by this module are specified here.
!
!  18-07-06/tony: coded
!
      logical, dimension(npencils) :: lpencil_in
!
      call keep_compiler_quiet(lpencil_in)
!
    endsubroutine pencil_interdep_special
!***********************************************************************
    subroutine calc_pencils_special(f,p)
!
!  Calculate Hydro pencils.
!  Most basic pencils should come first, as others may depend on them.
!
!   24-nov-04/tony: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
!
      intent(in) :: f
      intent(inout) :: p
!
      call keep_compiler_quiet(f)
      call keep_compiler_quiet(p)
!
    endsubroutine calc_pencils_special
!***********************************************************************
    subroutine dspecial_dt(f,df,p)
!
!  calculate right hand side of ONE OR MORE extra coupled PDEs
!  along the 'current' Pencil, i.e. f(l1:l2,m,n) where
!  m,n are global variables looped over in equ.f90
!
!  Due to the multi-step Runge Kutta timestepping used one MUST always
!  add to the present contents of the df array.  NEVER reset it to zero.
!
!  several precalculated Pencils of information are passed if for
!  efficiency.
!
!   06-oct-03/tony: coded
!
      use Diagnostics
      use Mpicomm
      use Sub
      use Deriv, only: der_pencil
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      real, dimension (3) :: meanx_oo
      real, dimension (3) :: meanx_uu
      real, dimension (nx,3) :: ufluct
      real, dimension (nx) :: ufluct2
      type (pencil_case) :: p
      integer :: i,j
      real, dimension (my) :: tmp,du_mean_dy
      real :: tau_tmp,nu
!
      intent(in) :: f,p
      intent(inout) :: df
!
!
    endsubroutine dspecial_dt
!***********************************************************************
    subroutine read_special_init_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=jet_init_pars, IOSTAT=iostat)
!
    endsubroutine read_special_init_pars
!***********************************************************************
    subroutine write_special_init_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=jet_init_pars)
!
    endsubroutine write_special_init_pars
!***********************************************************************
    subroutine read_special_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=jet_run_pars, IOSTAT=iostat)
!
    endsubroutine read_special_run_pars
!***********************************************************************
    subroutine write_special_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=jet_run_pars)
!
    endsubroutine write_special_run_pars
!***********************************************************************
    subroutine rprint_special(lreset,lwrite)
!
!  reads and registers print parameters relevant to special
!
!   06-oct-03/tony: coded
!
      use Diagnostics
!!$      use FArrayManager, only: farray_index_append
      use Sub
!
      integer :: iname
      logical :: lreset,lwr
      logical, optional :: lwrite
!
!!$      lwr = .false.
!!$      if (present(lwrite)) lwr=lwrite
!!$!
!!$!  reset everything in case of reset
!!$!  (this needs to be consistent with what is defined above!)
!!$!
!!$      if (lreset) then
!!$        idiag_turbint=0
!!$        idiag_tau_w=0
!!$        idiag_uxm_central=0
!!$      endif
!!$!
!!$      do iname=1,nname
!!$        call parse_name(iname,cname(iname),cform(iname),'turbint',idiag_turbint)
!!$        call parse_name(iname,cname(iname),cform(iname),'tau_w',idiag_tau_w)
!!$        call parse_name(iname,cname(iname),cform(iname),'uxm_central',idiag_uxm_central)
!!$      enddo
!!$!
!!$!  write column where which magnetic variable is stored
!!$      if (lwr) then
!!$        call farray_index_append('i_turbint',idiag_turbint)
!!$        call farray_index_append('i_tau_w',idiag_tau_w)
!!$        call farray_index_append('i_uxm_central',idiag_uxm_central)
!!$      endif
!
    endsubroutine rprint_special
!***********************************************************************
    subroutine get_slices_special(f,slices)
!
!  Write slices for animation of special variables.
!
!  26-jun-06/tony: dummy
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (slice_data) :: slices
      !
      !  Loop over slices
      !
!!$      select case (trim(slices%name))
!!$        !
!!$        !  Vorticity (derived variable)
!!$        !
!!$      case ('oo_meanx')
!!$        if (slices%index == 3) then
!!$          slices%ready = .false.
!!$        else
!!$          slices%index = slices%index+1
!!$          slices%xy=>oo_xy_meanx(:,:,slices%index)
!!$          if (slices%index < 3) slices%ready = .true.
!!$        endif
!!$      case ('uu_meanx')
!!$        if (slices%index >= 3) then
!!$          slices%ready = .false.
!!$        else
!!$          slices%index = slices%index+1
!!$          slices%xy=uu_xy_meanx(:,:,slices%index)
!!$          if (slices%index < 3) slices%ready = .true.
!!$        endif
!!$      endselect
!
    endsubroutine get_slices_special
!***********************************************************************
    subroutine special_after_boundary(f)
!
!  Mean flow velocitites
!
!  14-mar-08/nils: coded
!
      use Sub
      use Mpicomm, only: mpireduce_sum, mpibcast_real
!
      real, dimension (mx,my,mz,mfarray), intent(in) :: f
      real, dimension(nygrid,3) :: mean_u_tmp
      real :: faq
      integer :: j,k
!!$!
!!$!  calculate mean of velocity in xz planes
!!$!
!!$      if (lvideo.and.lfirst .or. ldiagnos) then
!!$        mean_u_tmp=0
!!$        faq=nxgrid*nzgrid
!!$        do j=m1,m2
!!$          do k=1,3
!!$            mean_u_tmp(j+ny*ipy-nghost,k)=sum(f(l1:l2,j,n1:n2,k+iux-1))/faq
!!$          enddo
!!$        enddo
!!$        do k=1,3
!!$          call mpireduce_sum(mean_u_tmp(:,k),mean_u(:,k),nygrid)
!!$          call mpibcast_real(mean_u(:,k),nygrid)
!!$        enddo
!!$      endif
!
    endsubroutine special_after_boundary
!***********************************************************************
    subroutine special_boundconds(f,bc)
!
!   Some precalculated pencils of data are passed in for efficiency
!   others may be calculated directly from the f array
!
!   2008-06-19/nils: coded
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
      type (boundary_condition) :: bc
      character (len=3) :: topbot
!
!
!
      topbot='top'
      if (bc%location==-1) topbot='bot'
!
      select case (bc%bcname)
      case ('tur')
        select case (bc%location)
        case (iBC_X_TOP)
          call bc_turb(f,bc%value1,'top',1,bc%ivar)
        case (iBC_X_BOT)
          call bc_turb(f,bc%value1,'bot',1,bc%ivar)
          bc%done=.true.
        end select
      case ('wi')
        call bc_wi_x(f,+1,topbot,bc%ivar,val=bc%value1)
        bc%done=.true.
      case ('wip')
        call bc_wip_x(f,+1,topbot,bc%ivar,val=bc%value1)
        bc%done=.true.
      end select
!
    endsubroutine special_boundconds
!***********************************************************************
    subroutine bc_turb(f,u_t,topbot,j,ivar)
!
! Use a prerun simulation as inlet condition
!
! 2010-10-15/nils: coded
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
      real :: shift, grid_shift, weight, round
      integer :: iround,lower,upper,ii,j,imin,imax,ivar
      character (len=3) :: topbot
      real :: T_t,u_t
!
! Read data from file only initially
! At later times this is stored in processor memory.
!
      if (first_time) then
        call read_turbulent_data(topbot,j)
      end if
!
! Set all ghost points at the same time
!
      if (ivar==iux) then
!
! Set which ghost points to update
!
        if (topbot=='bot') then
          imin=1
          imax=3
        else
          imin=l1+1
          imax=mx
        endif
!
! Set data at ghost points
!
        if (Lx_in == 0) call fatal_error('bc_nscbc_prf_x',&
            'Lx_in=0. Check that the precisions are the same.')
        round=t*u_t/Lx_in
        iround=int(round)
        shift=round-iround
        grid_shift=shift*nx_in
        lower=l1_in+int(grid_shift)
        upper=lower+1
        weight=grid_shift-int(grid_shift)
        f(imin:imax,m1:m2,n1:n2,iux:iuz)&
            =f_in(lower-3:lower-1,m1_in:m2_in,n1_in:n2_in,iux:iuy)*(1-weight)&
            +f_in(upper-3:upper-1,m1_in:m2_in,n1_in:n2_in,iux:iuy)*weight
        f(imin:imax,m1:m2,n1:n2,ilnrho)&
            =f_in(lower-3:lower-1,m1_in:m2_in,n1_in:n2_in,ilnrho)*(1-weight)&
            +f_in(upper-3:upper-1,m1_in:m2_in,n1_in:n2_in,ilnrho)*weight
!
! Add mean flow velocity on top of the turubulence
!
        f(imin:imax,m1:m2,n1:n2,iux)=f(imin:imax,m1:m2,n1:n2,iux)+u_t
      endif
!
    end subroutine bc_turb
!***********************************************************************
    subroutine special_before_boundary(f)
!
!   Possibility to modify the f array before the boundaries are
!   communicated.
!
!   Some precalculated pencils of data are passed in for efficiency
!   others may be calculated directly from the f array
!
!   06-jul-06/tony: coded
!
      real, dimension (mx,my,mz,mfarray), intent(in) :: f
!
      call keep_compiler_quiet(f)
!
    endsubroutine special_before_boundary
!***********************************************************************
    subroutine read_turbulent_data(topbot,j)
!
!   This subroutine will read data from pre-run isotropic box
!
!   19-jan-10/nils: coded
!
      use General, only: safe_character_assign, itoa
      use Sub, only : rdim
!
      character (len=3), intent(in) :: topbot
      integer, intent(in) :: j
      integer :: i,stat
      logical :: exist
      character (len=fnlen) :: file
!
! Read the size of the data to be found on the file.
!
      file=trim(turb_inlet_dir)//'/data/proc0/dim.dat'
      inquire(FILE=trim(file),EXIST=exist)
      if (exist) then
        call rdim(file,&
            mx_in,my_in,mz_in,mvar_in,maux_in,mglobal_in,prec_in,&
            nghost_in,ipx_in, ipy_in, ipz_in)
        nv_in=mvar_in+maux_in+mglobal_in
      else
        print*,'file=',file
        call fatal_error('bc_turb','Could not find file!')
      endif
!
! Allocate array for data to be used at the inlet.
! For now every processor reads all the data - this is clearly an overkill,
! but we leav it like this during the development of this feature.
!
      allocate( f_in(mx_in,my_in,mz_in,nv_in),STAT=stat)
      if (stat>0) call fatal_error('bc_turb',&
          "Couldn't allocate memory for f_in ")
      allocate( x_in(mx_in),STAT=stat)
      if (stat>0) call fatal_error('bc_turb',&
          "Couldn't allocate memory for x_in ")
      allocate( y_in(my_in),STAT=stat)
      if (stat>0) call fatal_error('bc_turb',&
          "Couldn't allocate memory for y_in ")
      allocate( z_in(mz_in),STAT=stat)
      if (stat>0) call fatal_error('bc_turb',&
          "Couldn't allocate memory for z_in ")
!
! Check which processor we want to read from.
! In the current implementation it is required that:
!   1) The number of mesh points and processors at the interface between
!      the two computational domains are equal. The two comp. domains I
!      am refering to here is the domain of the current simulation and the
!      domain of the pre-run isotropic turbulence simulation defining the
!      turbulence at the inlet.
!   2) The pre-run simulaion can not have multiple processors in the flow
!      direction of the current simulation.
!
      if (lprocz_slowest) then
        ipx_in=ipx
        ipy_in=ipy
        ipz_in=ipz
        nprocx_in=nprocx
        nprocy_in=nprocy
        nprocz_in=nprocz
        if (j==1) then
          if ((topbot=='bot'.and.lfirst_proc_x).or.&
              (topbot=='top'.and.llast_proc_x)) then
            proc_at_inlet=.true.
            ipx_in=0
            nprocx_in=1
          endif
        elseif (j==2) then
          if ((topbot=='bot'.and.lfirst_proc_y).or.&
              (topbot=='top'.and.llast_proc_y)) then
            proc_at_inlet=.true.
            ipy_in=0
            nprocy_in=1
          endif
        elseif (j==3) then
          if ((topbot=='bot'.and.lfirst_proc_z).or.&
              (topbot=='top'.and.llast_proc_z)) then
            proc_at_inlet=.true.
            ipz_in=0
            nprocz_in=1
          endif
        else
          call fatal_error("bc_turb_x",'No such direction!')
        endif
        iproc_in=ipz_in*nprocy_in*nprocx_in+ipy_in*nprocx_in+ipx_in
      else
        call fatal_error("bc_turb_x",&
            'lprocz_slowest=F not implemeted for inlet from file!')
      endif
!
!  Read data only if required, i.e. if we are at a processor handling inlets
!
      if (proc_at_inlet) then
        print*,'datadir=',datadir
        call safe_character_assign(directory_in,&
            trim(turb_inlet_dir)//'/data/proc'//itoa(iproc_in))
        print*,'directory_in=',directory_in
        call safe_character_assign(turbfile,&
            trim(directory_in)//'/var.dat')
        print*,'turbfile=',turbfile
        open(1,FILE=turbfile,FORM='unformatted')
        if (ip<=8) print*,'input: open, mx_in,my_in,mz_in,nv_in=',&
            mx_in,my_in,mz_in,nv_in
        read(1) f_in
        read(1) t_in,x_in,y_in,z_in,dx_in,dy_in,dz_in
        nx_in=mx_in-2*nghost
        ny_in=my_in-2*nghost
        nz_in=mz_in-2*nghost
        l1_in=nghost+1
        m1_in=nghost+1
        n1_in=nghost+1
        l2_in=mx_in-nghost
        m2_in=my_in-nghost
        n2_in=mz_in-nghost
        Lx_in=x_in(l2_in+1)-x_in(l1_in)
        Ly_in=y_in(m2_in+1)-y_in(m1_in)
        Lz_in=z_in(n2_in+1)-z_in(n1_in)
        first_time=.false.
        close(1)
      endif
!
    end subroutine read_turbulent_data
!***********************************************************************
    subroutine bc_wip_x(f,sgn,topbot,j,rel,val)
!
!
!  23-may-10/nils+marianne: coded
!
      character (len=3) :: topbot
      real, dimension (mx,my,mz,mfarray) :: f
      real :: val
      integer :: sgn,i,j
      logical, optional :: rel
      logical :: relative
!
      real :: radius, rad,z0,y0,yp0,zp0,dyp,minrad,distance,yp1
      integer :: jj,kk,nrad,irad,nb_lines,nb_colums,iz,iy
      real, allocatable, dimension(:,:,:) :: center
!
! Allocate array
!
      dyp=0.01
      radius=0.025
      distance=(dyp+radius)*2.
      yp0=xyz0(2)+radius*1.1
      yp1=xyz1(2)-radius*1.1
      nb_lines =floor((yp1-yp0)/(2*radius+dyp))
      nb_colums=max(floor(Lxyz(3)/(2*radius+dyp)),1)
      allocate(center(nb_lines,nb_colums,2))
!
!
      if (Lxyz(3) > 2*radius) then
        zp0=xyz0(3)+radius*1.3
      else
        zp0=xyz0(3)
      endif
      nrad=Lxyz(2)/(2*radius+dyp)
      z0=0.
!
! Find the center position of all the holes and put in array
!
      do iy=1,Nb_lines
        do iz=1,Nb_colums
          if (mod(iz,2)==0) then
            center(iy,iz,1)=yp0-distance/sqrt(2.)+(iy-1)*distance
          else
            center(iy,iz,1)=yp0+distance*(iy-1)
          endif
          center(iy,iz,2)=zp0-distance/sqrt(2.)*(iz-1)
        enddo
      enddo
!
!  Loop over all grid points
!
        do jj=1,my
          do kk=1,mz
            minrad=impossible
            do iy=1,Nb_lines
              do iz=1,Nb_colums
                rad=sqrt((y(jj)-center(iy,iz,1))**2+(z(kk)-center(iy,iz,2))**2)
                if (rad < minrad) minrad=rad
              enddo
            enddo
!
!  Zero derivative for density
!
            if (j == ilnrho) then
              do i=1,nghost
                f(l1-i,jj,kk,j)= f(l1+i,jj,kk,j)
              enddo
!
!  Constant temperature
!
            elseif (j == ilnTT) then
              f(l1,jj,kk,j)=val
              do i=1,nghost
                f(l1-i,jj,kk,j)=2*val-f(l1+i,jj,kk,j)
              enddo
            else
!
! Check if we are inside the radius if the inlet
!
              if (minrad > radius) then
            ! Solid wall
                do i=1,nghost
                  if (j <= iuz) then
                    f(l1,jj,kk,j)=0.
                    f(l1-i,jj,kk,j)=-f(l1+i,jj,kk,j)
                  else
                    f(l1-i,jj,kk,j)= f(l1+i,jj,kk,j)
                  endif
                enddo
              else
              ! Inlet
                if (j == iux) then
                  f(l1,jj,kk,j)=val*(1-(minrad/radius)**2)
                  do i=1,nghost
                    f(l1-i,jj,kk,j)=2*val*(1-(minrad/radius)**2)-f(l1+i,jj,kk,j)
                  enddo
                elseif (j == iuy .or. j == iuz) then
                  f(l1,jj,kk,j)=0.
                  do i=1,nghost
                    f(l1-i,jj,kk,j)=-f(l1+i,jj,kk,j)
                  enddo
                else
                  do i=1,nghost
                    f(l1-i,jj,kk,j)=2*val-f(l1+i,jj,kk,j)
                  enddo
                endif
              endif
            endif
         enddo
        enddo
!
  endsubroutine bc_wip_x
!***********************************************************************
    subroutine bc_wi_x(f,sgn,topbot,j,rel,val)
!
!
!  23-may-10/nils+marianne: coded
!
      character (len=3) :: topbot
      real, dimension (mx,my,mz,mfarray) :: f
      real  :: val
      integer :: sgn,i,j
      logical, optional :: rel
      logical :: relative
      real :: radius, rad,z0,y0
      integer :: jj,kk
!
!
!
      radius=0.015
      z0=0.
      y0=0.
!
!  Loop over all grid points
!
        do jj=1,my
          do kk=1,mz
            rad=sqrt((y(jj)-y0)**2+(z(kk)-z0)**2)
!
!  Zero derivative for density
!
            if (j == ilnrho) then
              do i=1,nghost
                f(l1-i,jj,kk,j)= f(l1+i,jj,kk,j)
              enddo
!
!  Constant temperature
!
            elseif (j == ilnTT) then
              f(l1,jj,kk,j)=val
              do i=1,nghost
                f(l1-i,jj,kk,j)=2*val-f(l1+i,jj,kk,j)
              enddo
            else
!
! Check if we are inside the radius if the inlet
!
              if (rad > radius) then
            ! Solid wall
                do i=1,nghost
                  if (j <= iuz) then
                    f(l1,jj,kk,j)=0.
                    f(l1-i,jj,kk,j)=-f(l1+i,jj,kk,j)
                  else
                    f(l1-i,jj,kk,j)= f(l1+i,jj,kk,j)
                  endif
                enddo
              else
              ! Inlet
                if (j == iux) then
                  f(l1,jj,kk,j)=val*(1-(rad/radius)**2)
                  do i=1,nghost
                    f(l1-i,jj,kk,j)=2*val*(1-(rad/radius)**2)-f(l1+i,jj,kk,j)
                  enddo
                elseif (j == iuy .or. j == iuz) then
                  f(l1,jj,kk,j)=0.
                  do i=1,nghost
                    f(l1-i,jj,kk,j)=-f(l1+i,jj,kk,j)
                  enddo
                else
                  do i=1,nghost
                    f(l1-i,jj,kk,j)=2*val-f(l1+i,jj,kk,j)
                  enddo
                endif
              endif
            endif
         enddo
        enddo
!
  endsubroutine bc_wi_x
!!$!***********************************************************************
!!$    subroutine bc_wo_x(f,sgn,topbot,j,rel,val)
!!$!
!!$!  23-may-10/nils: coded
!!$!
!!$      use EquationOfState
!!$      use chemistry
!!$!
!!$      character (len=3) :: topbot
!!$      real, dimension (mx,my,mz,mfarray) :: f
!!$      real :: val
!!$      integer :: sgn,i,j
!!$      logical, optional :: rel
!!$      logical :: relative
!!$
!!$      real :: radius, rad,z0,y0,T0,P0,r,mu
!!$      integer :: jj,kk
!!$!
!!$!
!!$!
!!$      radius=0.0
!!$      z0=0.
!!$      y0=0.
!!$!
!!$! First of all we set one-sided derivatives for all variables over the full
!!$! boundary. This will then mostly be overwritten later.
!!$
!!$!NILS: COmmented out the below statements for now in order to make the
!!$!NILS: code compile after I moved this from the boundary module.
!!$!NILS: will have to duplicate whatever bc i need here since this
!!$!NILS: module is used by the boundary module.....
!!$!
!!$!      call bc_onesided_x(f,topbot,j)
!!$!      call bc_extrap_2_1(f,topbot,j)
!!$!
!!$!  Loop over all grid points
!!$!
!!$        do jj=1,my
!!$          do kk=1,mz
!!$            rad=sqrt((y(jj)-y0)**2+(z(kk)-z0)**2)
!!$!
!!$!  Zero derivative for density
!!$!
!!$            if (j == ilnrho) then
!!$              do i=1,nghost
!!$                f(l2+i,jj,kk,j)= f(l2-i,jj,kk,j)
!!$              enddo
!!$            else
!!$!
!!$! Check if we are inside the radius if the inlet
!!$!
!!$              if (rad > radius) then
!!$                ! Solid wall
!!$                do i=1,nghost
!!$                  if (j <= iuz) then
!!$                    f(l2,jj,kk,j)=0.
!!$                    f(l2+i,jj,kk,j)=-f(l2-i,jj,kk,j)
!!$                  elseif (j == ilnTT) then
!!$                    f(l2,jj,kk,j)=val
!!$                    f(l2+i,jj,kk,j)=2*val-f(l2-i,jj,kk,j)
!!$                  else
!!$                    f(l2+i,jj,kk,j)= f(l2-i,jj,kk,j)
!!$                  endif
!!$                enddo
!!$              else
!!$                ! Outlet
!!$                if (j == ilnTT) then
!!$                  call getmu(f,mu,l2,jj,kk)
!!$                  r=Rgas*mu
!!$                  P0=1.013e6
!!$                  T0=min(P0/(exp(f(l2,jj,kk,ilnrho))*r),2500.)
!!$                 f(l2,jj,kk,j)=log(T0)
!!$                  do i=1,nghost
!!$                    f(l2+i,jj,kk,j)=2*log(T0)-f(l2-i,jj,kk,j)
!!$                  enddo
!!$                else
!!$                  ! Do nothing because one sided conditions has already
!!$                  ! been set in the top of the rutine
!!$                endif
!!$              endif
!!$            endif
!!$         enddo
!!$        enddo
!!$!
!!$  endsubroutine bc_wo_x
!***********************************************************************
!
!********************************************************************
!************        DO NOT DELETE THE FOLLOWING       **************
!********************************************************************
!**  This is an automatically generated include file that creates  **
!**  copies dummy routines from nospecial.f90 for any Special      **
!**  routines not implemented in this file                         **
!**                                                                **
    include '../special_dummies.inc'
!********************************************************************
!
endmodule Special