! $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, iZZ_chiral=0
  integer :: iXX2_chiral=0, iYY2_chiral=0
  character (len=labellen) :: initXX_chiral='zero', initYY_chiral='zero', initZZ_chiral='zero'
  logical :: llorentzforceEP=.false., lZZ_chiral=.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 :: amplZZ_chiral=.1, widthZZ_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.
  real :: initpowerYY_chiral=2., kpeakYY_chiral=0., cutoffYY_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, &
       initZZ_chiral, amplZZ_chiral, lZZ_chiral, &
       radiusXX_chiral, widthXX_chiral, &
       radiusYY_chiral, widthYY_chiral, &
       xposXX_chiral, yposXX_chiral, zposXX_chiral, &
       xposYY_chiral, yposYY_chiral, zposYY_chiral, &
       initpowerYY_chiral, kpeakYY_chiral, cutoffYY_chiral
!
  real :: chiral_diffXX=impossible, chiral_diff=0., chiral_crossinhibition=1.,chiral_fidelity=1.
  real :: chiral_diffZZ=impossible
  real :: chiral_fishernu=0., chiral_fisherK=1. !fishers equation growth rate, carrying capac.
  real :: chiral_fishermu=0., chiral_fisherH=1. !fishers equation growth rate, carrying capac.
  real :: chiral_fisherR=0., chiral_fisherR2=0. !(reinfections, second model corresponds to SIRS)
  real :: chiral_fisherR2_tend=0., chiral_fisherR2_tstart=0. !(second model, tend and start time)
  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_diffZZ, chiral_crossinhibition, chiral_fidelity, &
       chiral_reaction, limposed_gradient, gradX0, gradY0, &
       chiral_fishernu, chiral_fisherK, chiral_fisherR, chiral_fisherR2, &
       chiral_fishermu, chiral_fisherH, &
       lupw_chiral, llorentzforceEP, &
       linitialize_aa_from_EP,tinitialize_aa_from_EP, &
       linitialize_aa_from_EP_alpgrad,linitialize_aa_from_EP_betgrad, &
       chiral_fisherR2_tend, chiral_fisherR2_tstart
!
  integer :: idiag_XX_chiralmax=0, idiag_XX_chiralm=0
  integer :: idiag_YY_chiralmax=0, idiag_YY_chiralm=0
  integer :: idiag_ZZ_chiralmax=0, idiag_ZZ_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
  integer :: idiag_R2tdep=0
!
! Auxiliary variables.
!
  real, dimension (nx,3) :: gXX_chiral, gYY_chiral, gZZ_chiral
  real, dimension (nx) :: del2XX_chiral, del2YY_chiral,del2ZZ_chiral
  real :: chiral_fisherR2_tdep
!
  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)
      if (lZZ_chiral) call farray_register_pde('ZZ_chiral',iZZ_chiral)
!
!  extra fields for Higgs
!
      if (initYY_chiral=='power_randomphase_higgs') then
        call farray_register_pde('XX2_chiral',iXX2_chiral)
        call farray_register_pde('YY2_chiral',iYY2_chiral)
      endif
!
!  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 (chiral_diffZZ==impossible) chiral_diffZZ=chiral_diff
      if (headt) print*,'chiral_diffXX, chiral_diff, chiral_diffZZ =', &
                         chiral_diffXX, chiral_diff, chiral_diffZZ
!
!  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 Sub
      use General
      use Initcond
      use InitialCondition, only: initial_condition_chiral
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz) :: r
      real, dimension (mx,my,mz) :: norm1
!
!  check first for initXX_chiral
!
      select case (initXX_chiral)
        case ('zero'); f(:,:,:,iXX_chiral)=0.
        case ('const'); f(:,:,:,iXX_chiral)=amplXX_chiral
        case ('random_minus1_to1'); call random_number_wrapper(r); f(:,:,:,iXX_chiral)=2.*r-1.
        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 fatal_error('init_chiral','no such initXX_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 ('random_0to_pi'); call random_number_wrapper(r); f(:,:,:,iYY_chiral)=2.*pi*r
        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,'chiral_list')
        case ('power_randomphase')
          call power_randomphase(amplYY_chiral,initpowerYY_chiral,0.,kpeakYY_chiral,cutoffYY_chiral,&
            f,iXX_chiral,iXX_chiral,lscale_tobox=.false.)
          call power_randomphase(amplYY_chiral,initpowerYY_chiral,0.,kpeakYY_chiral,cutoffYY_chiral,&
            f,iYY_chiral,iYY_chiral,lscale_tobox=.false.)
        case ('power_randomphase_higgs')
          call power_randomphase(amplYY_chiral,initpowerYY_chiral,0.,kpeakYY_chiral,cutoffYY_chiral,&
            f,iXX_chiral,iXX_chiral,lscale_tobox=.false.)
          call power_randomphase(amplYY_chiral,initpowerYY_chiral,0.,kpeakYY_chiral,cutoffYY_chiral,&
            f,iYY_chiral,iYY_chiral,lscale_tobox=.false.)
          call power_randomphase(amplYY_chiral,initpowerYY_chiral,0.,kpeakYY_chiral,cutoffYY_chiral,&
            f,iXX2_chiral,iXX2_chiral,lscale_tobox=.false.)
          call power_randomphase(amplYY_chiral,initpowerYY_chiral,0.,kpeakYY_chiral,cutoffYY_chiral,&
            f,iYY2_chiral,iYY2_chiral,lscale_tobox=.false.)
          r=1./sqrt(f(:,:,:,iXX_chiral)**2 &
                   +f(:,:,:,iYY_chiral)**2 &
                   +f(:,:,:,iXX2_chiral)**2 &
                   +f(:,:,:,iYY2_chiral)**2)
          f(:,:,:,iXX_chiral) =r*f(:,:,:,iXX_chiral)
          f(:,:,:,iYY_chiral) =r*f(:,:,:,iYY_chiral)
          f(:,:,:,iXX2_chiral)=r*f(:,:,:,iXX2_chiral)
          f(:,:,:,iYY2_chiral)=r*f(:,:,:,iYY2_chiral)
        case default; call fatal_error('init_chiral','no such initYY_chiral: '//trim(initYY_chiral))
      endselect
!
!  check next for initZZ_chiral
!
      if (lZZ_chiral) then
        select case (initZZ_chiral)
          case ('zero'); f(:,:,:,iZZ_chiral)=0.
          case ('const'); f(:,:,:,iZZ_chiral)=amplZZ_chiral
          case ('chiral_list'); call chiral_list(amplZZ_chiral,f,iZZ_chiral,'chiral_listZZ')
          case default; call fatal_error('init_chiral','no such initZZ_chiral: '//trim(initZZ_chiral))
        endselect
      endif
!
!  Interface for user's own initial condition
!
      if (linitial_condition) call initial_condition_chiral(f)
!
    endsubroutine init_chiral
!***********************************************************************
    subroutine chiral_list(fact,f,i,filename)
!
!  Read intial conditions from file
!
!  13-apr-20/axel: coded
!
      integer :: nrow, irow, ix, iy, i, l, m, n=1
      character (len=*) :: filename
      real, dimension (mx,my,mz,mfarray) :: f
      real :: ampl, fact
!
      open(1,file=trim(filename)//'.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 Sub
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
!
      intent(in)  :: p
      intent(inout) :: f
      intent(inout) :: df
!
      real, dimension (nx,3) :: bbEP, jjEP, jxbEP
      real, dimension (nx) :: XX_chiral, YY_chiral, ZZ_chiral
      real, dimension (nx) :: ugXX_chiral, dXX_chiral
      real, dimension (nx) :: ugYY_chiral, dYY_chiral
      real, dimension (nx) :: ugZZ_chiral, dZZ_chiral
      real, dimension (nx) :: RRXX_chiral,XX2_chiral
      real, dimension (nx) :: RRYY_chiral,YY2_chiral
      real, dimension (nx) :: RR21_chiral
      real, dimension (nx) :: diffus_chiral
      real :: pp,qq
      integer :: j
!
!  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.
!
      if (lhydro) then
        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
      endif
!
!  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
!
!  Advection & diffusion for the containment in the SIR model.
!
      if (lZZ_chiral) then
        if (lhydro) then
          call u_dot_grad(f,iZZ_chiral,gZZ_chiral,p%uu,ugZZ_chiral,UPWIND=lupw_chiral)
          df(l1:l2,m,n,iZZ_chiral)=df(l1:l2,m,n,iZZ_chiral)-ugZZ_chiral
        endif
        call del2(f,iZZ_chiral,del2ZZ_chiral)
        df(l1:l2,m,n,iZZ_chiral)=df(l1:l2,m,n,iZZ_chiral)+chiral_diffZZ*del2ZZ_chiral
      endif
!
!  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 &
        +chiral_fisherR2_tdep*(1.-(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
!
!  Containment in the SIR model.
!
      if (lZZ_chiral) then
        ZZ_chiral=f(l1:l2,m,n,iZZ_chiral)
        dXX_chiral=-chiral_fishermu*XX_chiral*ZZ_chiral
        dZZ_chiral=+chiral_fishermu*XX_chiral*ZZ_chiral-chiral_fisherH*ZZ_chiral
        df(l1:l2,m,n,iXX_chiral)=df(l1:l2,m,n,iXX_chiral)+dXX_chiral
        df(l1:l2,m,n,iZZ_chiral)=df(l1:l2,m,n,iZZ_chiral)+dZZ_chiral
      endif
!
!  Alternative when nothing is set.
!
      case ('nothing')
        if (lroot.and.ip<=5) print*,"chiral_reaction='nothing'"
!
      case default
        call fatal_error('chiral_reaction','no such chiral_reaction: '//trim(chiral_reaction))
      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

      call calc_diagnostics_chiral(f,p)

    endsubroutine dXY_chiral_dt
!***********************************************************************
    subroutine calc_diagnostics_chiral(f,p)
!
!  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>
!
      use Diagnostics
      use Sub

      real, dimension (mx,my,mz,mfarray) :: f
      type (pencil_case) :: p

      real, dimension (nx) :: XX_chiral,YY_chiral,ZZ_chiral,QQ_chiral,QQ21QQ_chiral,bbEP2,jjEP2,jbEP
      real, dimension (nx,3) :: bbEP, jjEP
      real, dimension (nx,3,3) :: gXXij_chiral,gYYij_chiral
      real :: lamchiral
      integer :: j

      if (ldiagnos) then

        XX_chiral=f(l1:l2,m,n,iXX_chiral)
        YY_chiral=f(l1:l2,m,n,iYY_chiral)

        call max_mn_name(XX_chiral,idiag_XX_chiralmax)
        call max_mn_name(YY_chiral,idiag_YY_chiralmax)
        call sum_mn_name(XX_chiral,idiag_XX_chiralm)
        call sum_mn_name(YY_chiral,idiag_YY_chiralm)
        if (lZZ_chiral) then
          ZZ_chiral=f(l1:l2,m,n,iZZ_chiral)
          call sum_mn_name(ZZ_chiral,idiag_ZZ_chiralm)
          call max_mn_name(ZZ_chiral,idiag_ZZ_chiralmax)
        endif
        call save_name(chiral_fisherR2_tdep,idiag_R2tdep)
!
        if (idiag_QQm_chiral/=0 .or. idiag_QQ21m_chiral/=0 .or. idiag_QQ21QQm_chiral/=0) then
          QQ_chiral=XX_chiral-YY_chiral
          call sum_mn_name(QQ_chiral,idiag_QQm_chiral)
          if (idiag_QQ21m_chiral/=0) call sum_mn_name(1.-QQ_chiral**2,idiag_QQ21m_chiral)
          if (idiag_QQ21QQm_chiral/=0) then
            lamchiral=2.*chiral_fidelity-1.
            QQ21QQ_chiral=(lamchiral-QQ_chiral**2)/(1.+QQ_chiral**2)*QQ_chiral
            call sum_mn_name(QQ21QQ_chiral,idiag_QQ21QQm_chiral)
          endif
        endif
!
!  Calculate BB = gXX_chiral x gYY_chiral
!
        if (idiag_brmsEP/=0.or.idiag_bmaxEP/=0.or.idiag_jbmEP/=0) then
          call cross(gXX_chiral,gYY_chiral,bbEP)
          if (idiag_brmsEP/=0.or.idiag_bmaxEP/=0) then
            call dot2(bbEP,bbEP2)
            call sum_mn_name(bbEP2,idiag_brmsEP,lsqrt=.true.)
            call max_mn_name(bbEP2,idiag_bmaxEP,lsqrt=.true.)
          endif
        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)
          call sum_mn_name(jjEP2,idiag_jrmsEP,lsqrt=.true.)
          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
            call dot(jjEP,bbEP,jbEP)
            call sum_mn_name(jbEP,idiag_jbmEP)
          endif
        endif
      endif
!
    endsubroutine calc_diagnostics_chiral
!***********************************************************************
    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,gXX2_chiral,gYY2_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 (initYY_chiral=='power_randomphase_higgs') then
                call grad(f,iXX2_chiral,gXX2_chiral)
                call grad(f,iYY2_chiral,gYY2_chiral)
              endif
              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))
                  if (initYY_chiral=='power_randomphase_higgs') then
                    f(l1:l2,m,n,j+iaa-1)=f(l1:l2,m,n,j+iaa-1) &
                      +.5*( f(l1:l2,m,n,iXX2_chiral)*gYY2_chiral(:,j) &
                           -f(l1:l2,m,n,iYY2_chiral)*gXX2_chiral(:,j))
                    if (lroot) print*,'initialize_aa_from_Higgs at t=',t
                  endif
                endif
              enddo
            enddo
            enddo
          endif
        endif
      endif
!
!  time-dependent reinfection rate
!
      if (chiral_fisherR2_tend==0.) then
        chiral_fisherR2_tdep=chiral_fisherR2
      else
        chiral_fisherR2_tdep = chiral_fisherR2*(max(0.,1.-(max(0.,real(t)-chiral_fisherR2_tstart)/ &
                               (chiral_fisherR2_tend-chiral_fisherR2_tstart))**2)**2)
      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_ZZ_chiralmax=0; idiag_ZZ_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; idiag_R2tdep=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),&
            'ZZm',idiag_ZZ_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),&
            'ZZmax',idiag_ZZ_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)
        call parse_name(iname,cname(iname),cform(iname),&
            'R2tdep',idiag_R2tdep)
      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=='ZZ_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: ZZ
!
        case ('ZZ_chiral'); call assign_slices_scal(slices,f,iZZ_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
