! $Id$
!
!  This modules solves two reactive scalar advection equations
!  This is used for modeling the spatial evolution of left and
!  right handed aminoacids.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lchiral = .true.
!
! MVAR CONTRIBUTION 2
! MAUX CONTRIBUTION 0
!
!***************************************************************
module Chiral
!
  use Cparam
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
!
  implicit none
!
  include 'chiral.h'
!
  integer :: iXX_chiral=0,iYY_chiral=0
  character (len=labellen) :: initXX_chiral='zero',initYY_chiral='zero'
  logical :: llorentzforceEP=.false.
  logical :: linitialize_aa_from_EP=.false.
  logical :: linitialize_aa_from_EP_alpgrad=.false.
  logical :: linitialize_aa_from_EP_betgrad=.false.
  real :: tinitialize_aa_from_EP=0.
  real :: amplXX_chiral=.1, widthXX_chiral=.5
  real :: amplYY_chiral=.1, widthYY_chiral=.5
  real :: kx_XX_chiral=1.,ky_XX_chiral=1.,kz_XX_chiral=1.,radiusXX_chiral=0.
  real :: kx_YY_chiral=1.,ky_YY_chiral=1.,kz_YY_chiral=1.,radiusYY_chiral=0.
  real :: xposXX_chiral=0.,yposXX_chiral=0.,zposXX_chiral=0.
  real :: xposYY_chiral=0.,yposYY_chiral=0.,zposYY_chiral=0.
!
  namelist /chiral_init_pars/ &
       initXX_chiral,amplXX_chiral,kx_XX_chiral,ky_XX_chiral,kz_XX_chiral, &
       initYY_chiral,amplYY_chiral,kx_YY_chiral,ky_YY_chiral,kz_YY_chiral, &
       radiusXX_chiral,widthXX_chiral, &
       radiusYY_chiral,widthYY_chiral, &
       xposXX_chiral,yposXX_chiral,zposXX_chiral, &
       xposYY_chiral,yposYY_chiral,zposYY_chiral
!
  real :: chiral_diffXX=impossible, chiral_diff=0., chiral_crossinhibition=1.,chiral_fidelity=1.
  real :: chiral_fishernu=0., chiral_fisherK=1. !fishers equation growth rate, carrying capac.
  real :: chiral_fisherR=0. !(reinfection)
  real, dimension(3) :: gradX0=(/0.0,0.0,0.0/), gradY0=(/0.0,0.0,0.0/)
  character (len=labellen) :: chiral_reaction='BAHN_model'
  logical :: limposed_gradient=.false.
  logical :: lupw_chiral=.false.
!
  namelist /chiral_run_pars/ &
       chiral_diffXX, chiral_diff, chiral_crossinhibition, chiral_fidelity, &
       chiral_reaction, limposed_gradient, gradX0, gradY0, &
       chiral_fishernu, chiral_fisherK, chiral_fisherR, &
       lupw_chiral, llorentzforceEP, &
       linitialize_aa_from_EP,tinitialize_aa_from_EP, &
       linitialize_aa_from_EP_alpgrad,linitialize_aa_from_EP_betgrad
!
  integer :: idiag_XX_chiralmax=0, idiag_XX_chiralm=0
  integer :: idiag_YY_chiralmax=0, idiag_YY_chiralm=0
  integer :: idiag_QQm_chiral=0, idiag_QQ21m_chiral=0, idiag_QQ21QQm_chiral=0
  integer :: idiag_brmsEP=0, idiag_bmaxEP=0
  integer :: idiag_jrmsEP=0, idiag_jmaxEP=0
  integer :: idiag_jbmEP=0
!
  contains
!***********************************************************************
    subroutine register_chiral()
!
!  Initialise variables which should know that we solve for passive
!  scalar: iXX_chiral and iYY_chiral; increase nvar accordingly
!
!  28-may-04/axel: adapted from pscalar
!
      use FArrayManager
!
      call farray_register_pde('XX_chiral',iXX_chiral)
      call farray_register_pde('YY_chiral',iYY_chiral)
!
!  Identify version number.
!
      if (lroot) call svn_id( &
          "$Id$")
!
    endsubroutine register_chiral
!***********************************************************************
    subroutine initialize_chiral(f)
!
!  Perform any necessary post-parameter read initialization
!  Dummy routine
!
!  28-may-04/axel: adapted from pscalar
!
      real, dimension (mx,my,mz,mfarray) :: f
!
!  default: chiral_diffXX=chiral_diff
!
      if (chiral_diffXX==impossible) chiral_diffXX=chiral_diff
      if (headt) print*,'chiral_diffXX,chiral_diff=',chiral_diffXX,chiral_diff
!
!  set f to zero and then call the same initial condition
!  that was used in start.csh
!
      call keep_compiler_quiet(f)
!
    endsubroutine initialize_chiral
!***********************************************************************
    subroutine init_chiral(f)
!
!  initialise passive scalar field; called from start.f90
!
!  28-may-04/axel: adapted from pscalar
!
      use Mpicomm
      use Sub
      use Initcond
      use InitialCondition, only: initial_condition_chiral
!
      real, dimension (mx,my,mz,mfarray) :: f
!
!  check first for initXX_chiral
!
      select case (initXX_chiral)
        case ('zero'); f(:,:,:,iXX_chiral)=0.
        case ('const'); f(:,:,:,iXX_chiral)=amplXX_chiral
        case ('blob'); call blob(amplXX_chiral,f,iXX_chiral,radiusXX_chiral,xposXX_chiral,yposXX_chiral,zposXX_chiral)
        case ('hat-x'); call hat(amplXX_chiral,f,iXX_chiral,widthXX_chiral,kx=kx_XX_chiral)
        case ('hat-y'); call hat(amplXX_chiral,f,iXX_chiral,widthXX_chiral,ky=ky_XX_chiral)
        case ('hat-z'); call hat(amplXX_chiral,f,iXX_chiral,widthXX_chiral,kz=kz_XX_chiral)
        case ('gaussian-x'); call gaussian(amplXX_chiral,f,iXX_chiral,kx=kx_XX_chiral)
        case ('gaussian-y'); call gaussian(amplXX_chiral,f,iXX_chiral,ky=ky_XX_chiral)
        case ('gaussian-z'); call gaussian(amplXX_chiral,f,iXX_chiral,kz=kz_XX_chiral)
        case ('gaussian-noise'); call gaunoise(amplXX_chiral,f,iXX_chiral)
        case ('positive-noise'); call posnoise(amplXX_chiral,f,iXX_chiral)
        case ('wave-x'); call wave(amplXX_chiral,f,iXX_chiral,kx=kx_XX_chiral)
        case ('wave-y'); call wave(amplXX_chiral,f,iXX_chiral,ky=ky_XX_chiral)
        case ('wave-z'); call wave(amplXX_chiral,f,iXX_chiral,kz=kz_XX_chiral)
        case ('cos(x-cosz)'); call cosxz_cosz(amplXX_chiral,f,iXX_chiral,kx_XX_chiral,kz_XX_chiral)
        case ('cosy_sinz'); call cosy_sinz(amplXX_chiral,f,iXX_chiral,ky_XX_chiral,kz_XX_chiral)
        case ('cosx_cosz'); call cosx_cosz(amplXX_chiral,f,iXX_chiral,kx_XX_chiral,kz_XX_chiral)
        case ('cosx_cosy_cosz'); call cosx_cosy_cosz(amplXX_chiral,f,iXX_chiral,kx_XX_chiral,ky_XX_chiral,kz_XX_chiral)
        case ('cosx_siny_cosz'); call cosx_siny_cosz(amplXX_chiral,f,iXX_chiral,kx_XX_chiral,ky_XX_chiral,kz_XX_chiral)
        case default; call stop_it('init_chiral: bad init_chiral='//trim(initXX_chiral))
      endselect
!
!  check next for initYY_chiral
!
      select case (initYY_chiral)
        case ('zero'); f(:,:,:,iYY_chiral)=0.
        case ('const'); f(:,:,:,iYY_chiral)=amplYY_chiral
        case ('blob'); call blob(amplYY_chiral,f,iYY_chiral,radiusYY_chiral,xposYY_chiral,yposYY_chiral,zposYY_chiral)
        case ('hat-x'); call hat(amplYY_chiral,f,iYY_chiral,widthYY_chiral,kx=kx_YY_chiral)
        case ('hat-y'); call hat(amplYY_chiral,f,iYY_chiral,widthYY_chiral,ky=ky_YY_chiral)
        case ('hat-z'); call hat(amplYY_chiral,f,iYY_chiral,widthYY_chiral,kz=kz_YY_chiral)
        case ('gaussian-x'); call gaussian(amplYY_chiral,f,iYY_chiral,kx=kx_YY_chiral)
        case ('gaussian-y'); call gaussian(amplYY_chiral,f,iYY_chiral,ky=ky_YY_chiral)
        case ('gaussian-z'); call gaussian(amplYY_chiral,f,iYY_chiral,kz=kz_YY_chiral)
        case ('gaussian-noise'); call gaunoise(amplYY_chiral,f,iYY_chiral)
        case ('positive-noise'); call posnoise(amplYY_chiral,f,iYY_chiral)
        case ('wave-x'); call wave(amplYY_chiral,f,iYY_chiral,kx=kx_YY_chiral)
        case ('wave-y'); call wave(amplYY_chiral,f,iYY_chiral,ky=ky_YY_chiral)
        case ('wave-z'); call wave(amplYY_chiral,f,iYY_chiral,kz=kz_YY_chiral)
        case ('cos(y-sinz)'); call cosyz_sinz(amplYY_chiral,f,iYY_chiral,ky_YY_chiral,kz_YY_chiral)
        case ('cosy_sinz'); call cosy_sinz(amplYY_chiral,f,iYY_chiral,ky_YY_chiral,kz_YY_chiral)
        case ('cosy_cosz'); call cosy_cosz(amplYY_chiral,f,iYY_chiral,ky_YY_chiral,kz_YY_chiral)
        case ('cosx_cosy_cosz'); call cosx_cosy_cosz(amplYY_chiral,f,iYY_chiral,kx_YY_chiral,ky_YY_chiral,kz_YY_chiral)
        case ('cosx_siny_cosz'); call cosx_siny_cosz(amplYY_chiral,f,iYY_chiral,kx_YY_chiral,ky_YY_chiral,kz_YY_chiral)
        case ('chiral_list'); call chiral_list(amplYY_chiral,f,iYY_chiral)
        case default; call stop_it('init_chiral: bad init_chiral='//trim(initYY_chiral))
      endselect
!
!  Interface for user's own initial condition
!
      if (linitial_condition) call initial_condition_chiral(f)
!
    endsubroutine init_chiral
!***********************************************************************
    subroutine chiral_list(fact,f,i)
!
!  Read intial conditions from file
!
!  13-apr-20/axel: coded
!
      integer :: nrow, irow, ix, iy, i, l, m, n=1
      real, dimension (mx,my,mz,mfarray) :: f
      real :: ampl, fact
!
      open(1,file='chiral_list.dat')
      read(1,*) nrow
      do irow=1,nrow
        read(1,*) ix,iy,ampl
        if (ix/nx==ipx.and.iy/ny==ipy) then
          l=l1+mod(ix,nx)
          m=m1+mod(iy,ny)
          n=n1
          print*,'AXEL: ',ix,iy,ipx,ipy,l,m
          f(l,m,n,i)=f(l,m,n,i)+fact*ampl
        endif
      enddo
      close(1)
!
    endsubroutine chiral_list
!***********************************************************************
    subroutine pencil_criteria_chiral()
!
!  All pencils that the Chiral module depends on are specified here.
!
!  21-nov-04/anders: coded
!
      lpenc_requested(i_uu)=.true.
      if (llorentzforceEP) lpenc_requested(i_rho1)=.true.
!
    endsubroutine pencil_criteria_chiral
!***********************************************************************
    subroutine pencil_interdep_chiral(lpencil_in)
!
!  Interdependency among pencils provided by the Chiral module
!  is specified here.
!
!  21-11-04/anders: coded
!
      logical, dimension(npencils) :: lpencil_in
!
      call keep_compiler_quiet(lpencil_in)
!
    endsubroutine pencil_interdep_chiral
!***********************************************************************
    subroutine calc_pencils_chiral(f,p)
!
!  Calculate Chiral pencils.
!  Most basic pencils should come first, as others may depend on them.
!
!  21-11-04/anders: coded
!
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
!
      intent(in) :: f,p
!
      call keep_compiler_quiet(f)
      call keep_compiler_quiet(p)
!
    endsubroutine calc_pencils_chiral
!***********************************************************************
    subroutine dXY_chiral_dt(f,df,p)
!
!  passive scalar evolution
!  calculate chirality equations in reduced form; see q-bio/0401036
!
!  28-may-04/axel: adapted from pscalar
!   1-jul-09/axel: included gradX, gradY, and allowed for different reactions
!
      use Diagnostics
      use Sub
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
!
      real, dimension (nx,3,3) :: gXXij_chiral,gYYij_chiral
      real, dimension (nx,3) :: gXX_chiral,gYY_chiral,bbEP,jjEP,jxbEP
      real, dimension (nx) :: bbEP2,jjEP2,jbEP
      real, dimension (nx) :: XX_chiral,ugXX_chiral,del2XX_chiral,dXX_chiral
      real, dimension (nx) :: YY_chiral,ugYY_chiral,del2YY_chiral,dYY_chiral
      real, dimension (nx) :: RRXX_chiral,XX2_chiral
      real, dimension (nx) :: RRYY_chiral,YY2_chiral
      real, dimension (nx) :: RR21_chiral
      real, dimension (nx) :: QQ_chiral,QQ21_chiral,QQ21QQ_chiral
      real, dimension (nx) :: diffus_chiral
      real :: pp,qq,lamchiral
      integer :: j
!
      intent(in)  :: p
      intent(inout) :: f
      intent(inout) :: df
!
!  identify module and boundary conditions
!
      if (headtt.or.ldebug) print*,'SOLVE dXY_dt'
      if (headtt) call identify_bcs('XX_chiral',iXX_chiral)
      if (headtt) call identify_bcs('YY_chiral',iYY_chiral)
!
!  gradient of passive scalar
!
      call grad(f,iXX_chiral,gXX_chiral)
      call grad(f,iYY_chiral,gYY_chiral)
!
!  Add diffusion of imposed spatially constant gradient of X or Y.
!  This makes sense mainly for periodic boundary conditions.
!
      if (limposed_gradient) then
        do j=1,3
          gXX_chiral(:,j)=gXX_chiral(:,j)+gradX0(j)
          gYY_chiral(:,j)=gYY_chiral(:,j)+gradY0(j)
        enddo
      endif
!
!  advection term
!
      !call dot_mn(p%uu,gXX_chiral,ugXX_chiral)
      !call dot_mn(p%uu,gYY_chiral,ugYY_chiral)
      call u_dot_grad(f,iXX_chiral,gXX_chiral,p%uu,ugXX_chiral,UPWIND=lupw_chiral)
      call u_dot_grad(f,iYY_chiral,gYY_chiral,p%uu,ugYY_chiral,UPWIND=lupw_chiral)
      df(l1:l2,m,n,iXX_chiral)=df(l1:l2,m,n,iXX_chiral)-ugXX_chiral
      df(l1:l2,m,n,iYY_chiral)=df(l1:l2,m,n,iYY_chiral)-ugYY_chiral
!
!  diffusion term
!
      call del2(f,iXX_chiral,del2XX_chiral)
      call del2(f,iYY_chiral,del2YY_chiral)
      df(l1:l2,m,n,iXX_chiral)=df(l1:l2,m,n,iXX_chiral)+chiral_diffXX*del2XX_chiral
      df(l1:l2,m,n,iYY_chiral)=df(l1:l2,m,n,iYY_chiral)+chiral_diff  *del2YY_chiral
!
!  For Euler Potentials, possibility to add Lorentz force
!
      if (llorentzforceEP) then
        if (lhydro) then
          call cross(gXX_chiral,gYY_chiral,bbEP)
          do j=1,3
            jjEP(:,j)=gXX_chiral(:,j)*del2YY_chiral &
                     -gYY_chiral(:,j)*del2XX_chiral
          enddo
          call cross(jjEP,bbEP,jxbEP)
          if (ldensity) call multsv_mn(p%rho1,jxbEP,jxbEP)
          df(l1:l2,m,n,iux:iuz)=df(l1:l2,m,n,iux:iuz)+jxbEP
        endif
      endif
!
!  selection of different reaction terms
!
      select case (chiral_reaction)
!
!  BAHN model
!
      case ('BAHN_model')
!
!  reaction terms
!  X^2/Rtilde^2 - X*R
!  Y^2/Rtilde^2 - Y*R
!
!  for finite crossinhibition (=kI/kS) and finite fidelity (=f) we have
!  R --> RX=X+Y*kI/kS, R --> RY=Y+X*kI/kS, and
!  X2tilde=X^2/2RX, Y2tilde=Y^2/2RY.
!
      XX_chiral=f(l1:l2,m,n,iXX_chiral)
      YY_chiral=f(l1:l2,m,n,iYY_chiral)
      RRXX_chiral=XX_chiral+YY_chiral*chiral_crossinhibition
      RRYY_chiral=YY_chiral+XX_chiral*chiral_crossinhibition
!
!  abbreviations for quadratic quantities
!
      XX2_chiral=.5*XX_chiral**2/max(RRXX_chiral, tini)
      YY2_chiral=.5*YY_chiral**2/max(RRYY_chiral, tini)
      RR21_chiral=1./max(XX2_chiral+YY2_chiral, tini)
!
!  fidelity factor
!
      pp=.5*(1.+chiral_fidelity)
      qq=.5*(1.-chiral_fidelity)
!
!  final reaction equation
!
      dXX_chiral=(pp*XX2_chiral+qq*YY2_chiral)*RR21_chiral-XX_chiral*RRXX_chiral
      dYY_chiral=(pp*YY2_chiral+qq*XX2_chiral)*RR21_chiral-YY_chiral*RRYY_chiral
      df(l1:l2,m,n,iXX_chiral)=df(l1:l2,m,n,iXX_chiral)+dXX_chiral
      df(l1:l2,m,n,iYY_chiral)=df(l1:l2,m,n,iYY_chiral)+dYY_chiral
!
      case('fisher')
      if (headtt) print*,"chiral_reaction='fishers equation'"
      if (headtt) print*,"growth rate=", chiral_fishernu,"carrying capacity=",chiral_fisherK
      XX_chiral=f(l1:l2,m,n,iXX_chiral)
      YY_chiral=f(l1:l2,m,n,iYY_chiral)
      dXX_chiral=(1.-XX_chiral/chiral_fisherK)*XX_chiral*chiral_fishernu
      dYY_chiral=(1.-YY_chiral/chiral_fisherK)*YY_chiral*chiral_fishernu
      df(l1:l2,m,n,iXX_chiral)=df(l1:l2,m,n,iXX_chiral)+dXX_chiral
      df(l1:l2,m,n,iYY_chiral)=df(l1:l2,m,n,iYY_chiral)+dYY_chiral
!
      case('SIR')
      if (headtt) print*,"chiral_reaction='SIR equation'"
      if (headtt) print*,"growth rate=", chiral_fishernu,"carrying capacity=",chiral_fisherK
      XX_chiral=f(l1:l2,m,n,iXX_chiral)
      YY_chiral=f(l1:l2,m,n,iYY_chiral)
      dXX_chiral=-chiral_fishernu*XX_chiral*YY_chiral
      dYY_chiral=+chiral_fishernu*XX_chiral*YY_chiral-chiral_fisherK*YY_chiral &
        +chiral_fisherR*(1.-(XX_chiral+YY_chiral))
      df(l1:l2,m,n,iXX_chiral)=df(l1:l2,m,n,iXX_chiral)+dXX_chiral
      df(l1:l2,m,n,iYY_chiral)=df(l1:l2,m,n,iYY_chiral)+dYY_chiral
!
      case ('nothing')
        if (lroot.and.ip<=5) print*,"chiral_reaction='nothing'"
!
      case default
        write(unit=errormsg,fmt=*) &
             'dXY_chiral_dt: No such value for chiral_reaction: ', &
             trim(chiral_reaction)
        call fatal_error('chiral_reaction',errormsg)
      endselect
!
!  For the timestep calculation, need maximum diffusion
!
      if (lfirst.and.ldt) then
        diffus_chiral=max(chiral_diffXX,chiral_diff)*dxyz_2
        if (headtt.or.ldebug) print*,'dXY_chiral_dt: max(diffus_chiral) =', &
                                      maxval(diffus_chiral)
        maxdiffus=max(maxdiffus,diffus_chiral)
      endif
!
!  diagnostics
!
!  output for double and triple correlators (assume z-gradient of cc)
!  <u_k u_j d_j c> = <u_k c uu.gradXX_chiral>
!
      if (ldiagnos) then
        if (idiag_XX_chiralmax/=0) &
            call max_mn_name(XX_chiral,idiag_XX_chiralmax)
        if (idiag_YY_chiralmax/=0) &
            call max_mn_name(YY_chiral,idiag_YY_chiralmax)
        if (idiag_XX_chiralm/=0) call sum_mn_name(XX_chiral,idiag_XX_chiralm)
        if (idiag_YY_chiralm/=0) call sum_mn_name(YY_chiral,idiag_YY_chiralm)
!
!  extra diagnostics
!
        lamchiral=2.*chiral_fidelity-1.
        QQ_chiral=XX_chiral-YY_chiral
        QQ21_chiral=1.-QQ_chiral**2
        QQ21QQ_chiral=(lamchiral-QQ_chiral**2)/(1.+QQ_chiral**2)*QQ_chiral
        if (idiag_QQm_chiral/=0) call sum_mn_name(QQ_chiral,idiag_QQm_chiral)
        if (idiag_QQ21m_chiral/=0) &
            call sum_mn_name(QQ21_chiral,idiag_QQ21m_chiral)
        if (idiag_QQ21QQm_chiral/=0) &
            call sum_mn_name(QQ21QQ_chiral,idiag_QQ21QQm_chiral)
!
!  Calculate BB = gXX_chiral x gYY_chiral
!
        if (idiag_brmsEP/=0.or.idiag_bmaxEP/=0) then
          call cross(gXX_chiral,gYY_chiral,bbEP)
          call dot2(bbEP,bbEP2)
          if (idiag_brmsEP/=0) &
              call sum_mn_name(bbEP2,idiag_brmsEP,lsqrt=.true.)
          if (idiag_bmaxEP/=0) &
              call max_mn_name(bbEP2,idiag_bmaxEP,lsqrt=.true.)
        endif
!
!  Calculate Ji = Xi*del2Y - Yi*del2X + Xij Yj - Yij Xj
!  and then the max and rms values of that
!
        if (idiag_jrmsEP/=0.or.idiag_jmaxEP/=0.or.idiag_jbmEP/=0) then
          do j=1,3
            jjEP(:,j)=gXX_chiral(:,j)*del2YY_chiral &
                     -gYY_chiral(:,j)*del2XX_chiral
          enddo
          call g2ij(f,iXX_chiral,gXXij_chiral)
          call g2ij(f,iYY_chiral,gYYij_chiral)
          call multmv_mn(+gXXij_chiral,gYY_chiral,jjEP,ladd=.true.)
          call multmv_mn(-gYYij_chiral,gXX_chiral,jjEP,ladd=.true.)
          call dot2(jjEP,jjEP2)
          if (idiag_jrmsEP/=0) &
              call sum_mn_name(jjEP2,idiag_jrmsEP,lsqrt=.true.)
          if (idiag_jmaxEP/=0) &
              call max_mn_name(jjEP2,idiag_jmaxEP,lsqrt=.true.)
!
!  For J.B, we need B, but only if not already calculated.
!
          if (idiag_jbmEP/=0) then
            if (idiag_brmsEP/=0.or.idiag_bmaxEP/=0) then
              call cross(gXX_chiral,gYY_chiral,bbEP)
            endif
            call dot(jjEP,bbEP,jbEP)
            call sum_mn_name(jbEP,idiag_jbmEP)
          endif
        endif
      endif
!
    endsubroutine dXY_chiral_dt
!***********************************************************************
    subroutine chiral_before_boundary(f)
!
!  initialize aa from Euler potentials (EP), but only once when t=0
!
!   4-jul-09/axel: coded
!
      use Cdata
      use Sub
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: gXX_chiral,gYY_chiral
      integer :: j
!
      intent(inout) :: f
!
!  Possibility to initialize magnetic vector potential in terms
!  of Euler potentials X and Y.
!
      if (linitialize_aa_from_EP) then
        if (t>=tinitialize_aa_from_EP) then
          linitialize_aa_from_EP=.false.
          if (lmagnetic) then
            if (lroot) print*,'initialize_aa_from_EP at t=',t
            do n=n1,n2
            do m=m1,m2
              call grad(f,iXX_chiral,gXX_chiral)
              call grad(f,iYY_chiral,gYY_chiral)
              if (limposed_gradient) then
                do j=1,3
                  gXX_chiral(:,j)=gXX_chiral(:,j)+gradX0(j)
                  gYY_chiral(:,j)=gYY_chiral(:,j)+gradY0(j)
                enddo
              endif
              do j=1,3
                if (linitialize_aa_from_EP_alpgrad) then
                  f(l1:l2,m,n,j+iaa-1)=f(l1:l2,m,n,iXX_chiral)*gYY_chiral(:,j)
                elseif (linitialize_aa_from_EP_betgrad) then
                  f(l1:l2,m,n,j+iaa-1)=-f(l1:l2,m,n,iYY_chiral)*gXX_chiral(:,j)
                else
                  f(l1:l2,m,n,j+iaa-1)=.5*( &
                    f(l1:l2,m,n,iXX_chiral)*gYY_chiral(:,j) &
                   -f(l1:l2,m,n,iYY_chiral)*gXX_chiral(:,j))
                endif
              enddo
            enddo
            enddo
          endif
        endif
      endif
!
    endsubroutine chiral_before_boundary
!***********************************************************************
    subroutine read_chiral_init_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=chiral_init_pars, IOSTAT=iostat)
!
    endsubroutine read_chiral_init_pars
!***********************************************************************
    subroutine write_chiral_init_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=chiral_init_pars)
!
    endsubroutine write_chiral_init_pars
!***********************************************************************
    subroutine read_chiral_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=chiral_run_pars, IOSTAT=iostat)
!
    endsubroutine read_chiral_run_pars
!***********************************************************************
    subroutine write_chiral_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=chiral_run_pars)
!
    endsubroutine write_chiral_run_pars
!***********************************************************************
    subroutine rprint_chiral(lreset,lwrite)
!
!  reads and registers print parameters relevant for chirality
!
!  28-may-04/axel: adapted from pscalar
!
      use Diagnostics
!
      integer :: iname
      logical :: lreset
      logical, optional :: lwrite
!
!  reset everything in case of reset
!  (this needs to be consistent with what is defined above!)
!
      if (lreset) then
        idiag_XX_chiralmax=0; idiag_XX_chiralm=0
        idiag_YY_chiralmax=0; idiag_YY_chiralm=0
        idiag_QQm_chiral=0; idiag_QQ21m_chiral=0; idiag_QQ21QQm_chiral=0
        idiag_brmsEP=0; idiag_bmaxEP=0
        idiag_jrmsEP=0; idiag_jmaxEP=0
        idiag_jbmEP=0
      endif
!
!  check for those quantities that we want to evaluate online
!
      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),&
            'XXm',idiag_XX_chiralm)
        call parse_name(iname,cname(iname),cform(iname),&
            'YYm',idiag_YY_chiralm)
        call parse_name(iname,cname(iname),cform(iname),&
            'XXmax',idiag_XX_chiralmax)
        call parse_name(iname,cname(iname),cform(iname),&
            'YYmax',idiag_YY_chiralmax)
        call parse_name(iname,cname(iname),cform(iname),&
            'QQm',idiag_QQm_chiral)
        call parse_name(iname,cname(iname),cform(iname),&
            'QQ21m',idiag_QQ21m_chiral)
        call parse_name(iname,cname(iname),cform(iname),&
            'QQ21QQm',idiag_QQ21QQm_chiral)
        call parse_name(iname,cname(iname),cform(iname),&
            'brmsEP',idiag_brmsEP)
        call parse_name(iname,cname(iname),cform(iname),&
            'bmaxEP',idiag_bmaxEP)
        call parse_name(iname,cname(iname),cform(iname),&
            'jrmsEP',idiag_jrmsEP)
        call parse_name(iname,cname(iname),cform(iname),&
            'jmaxEP',idiag_jmaxEP)
        call parse_name(iname,cname(iname),cform(iname),&
            'jbmEP',idiag_jbmEP)
      enddo
!
!  check for those quantities for which we want video slices
!
      if (lwrite_slices) then
        where(cnamev=='XX_chiral'.or.cnamev=='YY_chiral'.or.cnamev=='DQ_chiral') &
          cformv='DEFINED'
      endif
!
    endsubroutine rprint_chiral
!***********************************************************************
    subroutine get_slices_chiral(f,slices)
!
!  Write slices for animation of Chiral variables.
!
!  26-jul-06/tony: coded
!
      use Slices_methods, only:assign_slices_scal

      real, dimension (mx,my,mz,mfarray) :: f
      type (slice_data) :: slices
!
!  Loop over slices
!
      select case (trim(slices%name))
!
!  Chirality fields: XX
!
        case ('XX_chiral'); call assign_slices_scal(slices,f,iXX_chiral)
!
!  Chirality fields: YY
!
        case ('YY_chiral'); call assign_slices_scal(slices,f,iYY_chiral)
!
!  Chirality fields: DQ
!
        case ('DQ_chiral')
          if (lwrite_slice_yz ) &
            call QQ_chiral(f(ix_loc,m1:m2,n1:n2,iXX_chiral)-f(ix_loc,m1:m2,n1:n2,iYY_chiral),slices%yz)
          if (lwrite_slice_xz ) &
            call QQ_chiral(f(l1:l2,iy_loc,n1:n2,iXX_chiral)-f(l1:l2,iy_loc,n1:n2,iYY_chiral),slices%xz)
          if (lwrite_slice_xy ) &
            call QQ_chiral(f(l1:l2,m1:m2,iz_loc,iXX_chiral)-f(l1:l2,m1:m2,iz_loc,iYY_chiral),slices%xy)
          if (lwrite_slice_xy2) &
            call QQ_chiral(f(l1:l2,m1:m2,iz2_loc,iXX_chiral)-f(l1:l2,m1:m2,iz2_loc,iYY_chiral),slices%xy2)
          if (lwrite_slice_xy3) &
            call QQ_chiral(f(l1:l2,m1:m2,iz3_loc,iXX_chiral)-f(l1:l2,m1:m2,iz3_loc,iYY_chiral),slices%xy3)
          if (lwrite_slice_xy4) &
            call QQ_chiral(f(l1:l2,m1:m2,iz4_loc,iXX_chiral)-f(l1:l2,m1:m2,iz4_loc,iYY_chiral),slices%xy4)
          if (lwrite_slice_xz2) &
            call QQ_chiral(f(l1:l2,iy2_loc,n1:n2,iXX_chiral)-f(l1:l2,iy2_loc,n1:n2,iYY_chiral),slices%xz2)
          slices%ready=.true.
!
      endselect
!
contains
      subroutine QQ_chiral(source,dest)
!
        real, dimension(:,:) :: source,dest
!
        dest=source**2 
        dest=source*(1.-dest)/(1.+dest)
!
      endsubroutine QQ_chiral

    endsubroutine get_slices_chiral
!***********************************************************************
endmodule Chiral