! $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'
!
!   n_pivot = number of pivot points
!   udrive = strength of the forcing velocity
!   xp = pivot point
!   yp = pivot point
!   rad = distance of the footpoint to the pivot point
!   lam_twist = Gaussian RMS width
!   lam_u = coefficient for the exponential saturation for the velocity
!   r_profile = velocity profile in the radial direction
!   z_profile = velocity profile along the z-direction
!   lam_z = coefficient for the z-profiles
!   blinking = flag for the 'vortex blinking'
!   delay_blink = waiting time for the blinking process to commence
!   t_blink_up = time period for the driver to be positive
!   t_blink_0u = time period for the driver to be zero after positive driving
!   t_blink_down = time period for the driver to be negative
!   t_blink_0d= time period for the driver to be zero after negative driving
!
  integer :: n_pivot = 1
  real, dimension(6) :: udrive = (/0.1,0.,0.,0.,0.,0./)
  real, dimension(6) :: xp = (/0.,0.,0.,0.,0.,0./)
  real, dimension(6) :: yp = (/0.,0.,0.,0.,0.,0./)
  real, dimension(6) :: rad = (/2.,0.,0.,0.,0.,0./)
  real, dimension(6) :: lam_twist = (/0.17,0.,0.,0.,0.,0./)
  real :: lam_u = 1
  character (len=labellen) :: r_profile = 'gaussian'
  character (len=labellen) :: z_profile = 'exp'
  real :: lam_z = 1
  logical, save :: lblink = .False.
  real :: t_e3 = 8.
  real, dimension(6) :: delay_blink = (/0.,0.,0.,0.,0.,0./)
  real, dimension(6) :: t_blink_up = (/1.,0.,0.,0.,0.,0./)
  real, dimension(6) :: t_blink_0u = (/0.,0.,0.,0.,0.,0./)
  real, dimension(6) :: t_blink_down = (/0.,0.,0.,0.,0.,0./)
  real, dimension(6) :: t_blink_0d = (/0.,0.,0.,0.,0.,0./)
!    
! input parameters
  namelist /special_init_pars/ n_pivot
!
! run parameters
  namelist /special_run_pars/ &
    n_pivot, udrive, xp, yp, rad, lam_twist, lam_u, r_profile, z_profile, lam_z, &
    t_e3, lblink, delay_blink, t_blink_up, t_blink_0u, t_blink_down, t_blink_0d
!
  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 with lreloading indicating a RELOAD
!
!  07-may-2015/iomsn (Simon Candelaresi): coded
!
      real, dimension(mx,my,mz,mfarray) :: f
!
      call keep_compiler_quiet(f)
!
    endsubroutine initialize_special
!***********************************************************************
    subroutine init_special(f)
!
!  Initialize special condition; called by start.f90.
!
!  07-may-2015/iomsn (Simon Candelaresi): coded
!      
      real, dimension(mx,my,mz,mfarray), intent(out) :: f

      call keep_compiler_quiet(f)
!
    endsubroutine init_special
!***********************************************************************
    subroutine finalize_special(f)
!
!  Called right before exiting.
!
!  07-may-2015/iomsn (Simon Candelaresi): coded
!
      real, dimension(mx,my,mz,mfarray) :: f
!
      call keep_compiler_quiet(f)
!
    endsubroutine finalize_special
!***********************************************************************
    subroutine read_special_init_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=special_init_pars, IOSTAT=iostat)
!
    endsubroutine read_special_init_pars
!***********************************************************************
    subroutine write_special_init_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=special_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=special_run_pars, IOSTAT=iostat)
!
    endsubroutine read_special_run_pars
!***********************************************************************
    subroutine write_special_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=special_run_pars)
!
    endsubroutine write_special_run_pars
!***********************************************************************
    subroutine rprint_special(lreset,lwrite)
!
!  reads and registers print parameters relevant to special
!
!  07-may-2015/iomsn (Simon Candelaresi): coded
!
      logical :: lreset
      logical, optional :: lwrite
!
      call keep_compiler_quiet(lreset)
      call keep_compiler_quiet(lwrite)
!
    endsubroutine rprint_special
!***********************************************************************
    subroutine special_calc_hydro(f,df,p)
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
      real, dimension(mx,my,mz,mvar), intent(inout) :: df
      type (pencil_case), intent(in) :: p
!
      ! Apply driving of velocity field.
      if (.not. lpencil_check_at_work) &
          call vel_driver (f, df)
!
    endsubroutine special_calc_hydro
!***********************************************************************
  subroutine vel_driver(f, df)
!
! Drive bottom boundary horizontal velocities with given velocity field.
!
!  07-may-2015/iomsn (Simon Candelaresi): coded
!
      real, dimension(mx,my,mz,mfarray), intent(inout) :: f
      real, dimension(mx,my,mz,mvar), intent(inout) :: df
      integer :: ip
      real, dimension(mx) :: dist       ! distance of point to pivot point
      real, dimension(mx) :: dist0      ! distance for the normalization
      real, dimension(mx) :: ux, uy     ! velocity in x and y direction      
      real :: z_factor                  ! multiplication factor for the z-dependence
      integer :: blink                  ! either -1, 0 or 1, depending on the blinking stage
      real :: t_blink_tot               ! total time for a blink switching on, off and negative
      real :: xc, yc, kc                ! variable sofr the e3 braid
!
      if (r_profile == 'e3') n_pivot = 1
        
      do ip = 1, n_pivot
        dist = sqrt((x(l1:l2)-xp(ip))**2 + (y(m)-yp(ip))**2)
        
        select case (r_profile)
        case ('e3')
            yc = 0
            if (mod(int(t/t_e3),2) == 0) then
                xc = -1
                kc = -1
            else
                xc = 1
                kc = 1
            endif
            ux = -udrive(ip)*sqrt(2.)*kc*exp((-(x(l1:l2)-xc)**2-(y(m)-yc)**2)/2-(mod(t,t_e3)**2)/(t_e3/4)**2)*(-y(m)+yc)
            uy = -udrive(ip)*sqrt(2.)*kc*exp((-(x(l1:l2)-xc)**2-(y(m)-yc)**2)/2-(mod(t,t_e3)**2)/(t_e3/4)**2)*(x(l1:l2)-xc)
        case ('gaussian')
            ux = exp(-(dist-rad(ip))**2/(2*lam_twist(ip)**2))*udrive(ip)*(-y(m)+yp(ip))
            uy = exp(-(dist-rad(ip))**2/(2*lam_twist(ip)**2))*udrive(ip)*(x(l1:l2)-xp(ip))
            ! normalize
            dist0 = (rad(ip) + sqrt(rad(ip)**2 + 4*lam_twist(ip)**2))/2
            ux = ux/(exp(-(dist0-rad(ip))**2/(2*lam_twist(ip)**2))*dist0)
            uy = uy/(exp(-(dist0-rad(ip))**2/(2*lam_twist(ip)**2))*dist0)
        case ('linear_exp')
            ux = exp(-abs(dist-rad(ip))/lam_twist(ip))*udrive(ip)*(-y(m)+yp(ip))
            uy = exp(-abs(dist-rad(ip))/lam_twist(ip))*udrive(ip)*(x(l1:l2)-xp(ip))
            ! normalize
            ux = ux*exp(1/sqrt(lam_twist(ip)))/sqrt(lam_twist(ip)) 
            uy = uy*exp(1/sqrt(lam_twist(ip)))/sqrt(lam_twist(ip))
         case default
            ux = exp(-(dist-rad(ip))**2/(2*lam_twist(ip)**2))*udrive(ip)*(-y(m)+yp(ip))/dist
            uy = exp(-(dist-rad(ip))**2/(2*lam_twist(ip)**2))*udrive(ip)*(x(l1:l2)-xp(ip))/dist
        end select
        
        ! add z-dependence
        select case (z_profile)
        case ('sharp')
            z_factor = 0
            if (z(n) == xyz0(nghost)) z_factor = 1
        case ('exp')
            z_factor = exp(-(z(n)-xyz0(nghost))*lam_z)
            if (z(n) <= xyz0(nghost)) z_factor = 1
        case ('erf')
            z_factor = 1-erf((z(n)-xyz0(nghost))*lam_z)
            if (z(n) <= xyz0(nghost)) z_factor = 1
        case default
            z_factor = exp(-(z(n)-xyz0(nghost))*lam_z)
            if (z(n) <= xyz0(nghost)) z_factor = 1
        end select
        ux = z_factor*ux
        uy = z_factor*uy

        ! add 'vortex blinking' by switching the sign in time
        blink = 1
        t_blink_tot = t_blink_up(ip) + t_blink_0u(ip) + t_blink_down(ip) + t_blink_0d(ip)
        if (lblink) blink = 0
        if ((lblink) .and. (t >= delay_blink(ip))) then
            if (((t-delay_blink(ip))/t_blink_tot-floor((t-delay_blink(ip))/t_blink_tot))*t_blink_tot < t_blink_up(ip)) then
                blink = 1
            elseif ((((t-delay_blink(ip))/t_blink_tot-floor((t-delay_blink(ip))/t_blink_tot))*t_blink_tot >= t_blink_up(ip)) .and. &
                ((t-delay_blink(ip))/t_blink_tot-floor((t-delay_blink(ip))/t_blink_tot))*t_blink_tot < &
                (t_blink_up(ip)+t_blink_0u(ip))) then
                blink = 0
            elseif ((((t-delay_blink(ip))/t_blink_tot-floor((t-delay_blink(ip))/t_blink_tot))*t_blink_tot >= &
                (t_blink_up(ip)+t_blink_0u(ip))) .and. &
                ((t-delay_blink(ip))/t_blink_tot-floor((t-delay_blink(ip))/t_blink_tot))*t_blink_tot < &
                (t_blink_up(ip)+t_blink_0u(ip)+t_blink_down(ip))) then
                blink = -1
             else
                blink = 0
            endif
        endif
        
        df(l1:l2,m,n,iux) = df(l1:l2,m,n,iux) - (f(l1:l2,m,n,iux) - blink*ux)/lam_u*z_factor*dt
        df(l1:l2,m,n,iuy) = df(l1:l2,m,n,iuy) - (f(l1:l2,m,n,iuy) - blink*uy)/lam_u*z_factor*dt
      enddo
!
    endsubroutine vel_driver
!***********************************************************************
!************        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