! $Id$
!
!  This modules solves mean-field contributions to both the
!  induction and the momentum equations.
!
!** 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 :: lmagn_mf = .true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 0
!
! PENCILS PROVIDED mf_EMF(3); mf_EMFdotB; jxb_mf(3); jxbr_mf(3); chiB_mf
! PENCILS PROVIDED mf_qp; mf_Beq21
!
!***************************************************************
module Magnetic_meanfield
!
  use Cparam
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages, only: fatal_error,inevitably_fatal_error,svn_id
  use Magnetic_meanfield_demfdt
!
  implicit none
!
  include 'meanfield.h'
!
!  array for inputting alpha profile
!
  real, dimension (mx,my) :: alpha_input
  real, pointer :: B_ext2
  logical, pointer :: lweyl_gauge
!
  real, dimension (nx) :: kf_x, kf_x1
  real, dimension (my) :: kf_y
  real, dimension (mz) :: kf_z
!
! Parameters
!
  character (len=labellen) :: meanfield_kf_profile='const'
  character (len=labellen) :: meanfield_etat_profile='const'
  character (len=labellen) :: meanfield_Beq_profile='const',qp_model='atan'
  character (len=labellen) :: Omega_profile='nothing', alpha_profile='const'
  character (len=labellen) :: EMF_profile='nothing', delta_profile='const'
!
! Input parameters
!
  real :: Omega_ampl=0.0, dummy=0.0
  real :: alpha_effect=0.0, alpha_quenching=0.0, delta_effect=0.0, alpha_zz=0.
  real :: gamma_effect=0.0, gamma_quenching=0.0
  real :: chit_quenching=0.0, chi_t0=0.0
  real :: meanfield_etat=0.0, meanfield_etat_height=1., meanfield_pumping=1.
  real :: meanfield_delta_height=10., meanfield_delta_width=.05
  real :: meanfield_Beq=1.0,meanfield_Beq_height=0., meanfield_Beq2_height=0.
  real :: meanfield_etat_exp=1.0, uturb=.1
  real :: alpha_eps=0.0, alpha_pom0=0.0
  real :: x_surface=0., x_surface2=0., z_surface=0., qp_width=impossible, qpx_width=impossible
  real :: alpha_equator=impossible, alpha_equator_gap=0.0, alpha_gap_step=0.0
  real :: alpha_rmin=0.
  real :: alpha_cutoff_up=0.0, alpha_cutoff_down=0.0
  real :: meanfield_qs=0.0, meanfield_qp=0.0, meanfield_qe=0.0, meanfield_qa=0.0
  real :: meanfield_Bs=1.0, meanfield_Bp=1.0, meanfield_Be=1.0, meanfield_Ba=1.0
  real :: meanfield_qp1, meanfield_qs1, meanfield_qe1, meanfield_qa1
  real :: meanfield_etaB=0.0
  real :: x1_alp=0., x2_alp=0., y1_alp=0., y2_alp=0.
  logical :: lOmega_effect=.false., lalpha_Omega_approx=.false.
  logical :: lmeanfield_noalpm=.false., lmeanfield_pumping=.false.
  logical :: lmeanfield_jxb=.false., lmeanfield_jxb_with_vA2=.false.
  logical :: lmeanfield_chitB=.false., lignore_gradB2_inchiB=.false.
  logical :: lchit_with_glnTT=.false., lrho_chit=.true., lchit_Bext2_equil=.false.
  logical :: lturb_temp_diff=.false., lqp_profile=.false., lqpx_profile=.false.
!
  namelist /magn_mf_init_pars/ &
      x1_alp, x2_alp, y1_alp, y2_alp
!
! Run parameters
!
  real :: alpha_rmax=0.0, alpha_width=0.0, alpha_width2=0.0, alpha_exp=0.
  real :: meanfield_etat_width=0.0, meanfield_Beq_width=0.0
  real :: meanfield_kf=1.0, meanfield_kf_width=0.0, meanfield_kf_width2=0.0
  real :: Omega_rmax=0.0, Omega_rwidth=0.0
  real :: rhs_term_kx=0.0, rhs_term_ampl=0.0
  real :: rhs_term_amplz=0.0, rhs_term_amplphi=0.0
  real :: mf_qJ2=0.0, qp_aniso_factor=1.0
  real :: kx_alpha=1.
  real, dimension(3) :: alpha_aniso=0.
  real, dimension(3,3) :: alpha_tensor=0., eta_tensor=0.
  real, dimension(ny,3,3) :: alpha_tensor_y=0., eta_tensor_y=0.
  real, dimension(nz,2,2) :: alpha_tensor_z=0., eta_tensor_z=0.
  real, dimension(nx) :: rhs_termz, rhs_termy
  real, dimension(nx) :: etat_x, detat_x, rhs_term
  real, dimension(my) :: etat_y, detat_y
  real, dimension(mz) :: etat_z, detat_z, qp_profile, qp_profder !, qe_profder
  real, dimension(nx) :: qpx_profile, qpx_profder
  logical :: llarge_scale_velocity=.false.
  logical :: lEMF_profile=.false.
  logical :: lalpha_profile_total=.false., lalpha_aniso=.false.
  logical :: ldelta_profile=.false., lalpha_tensor=.false., leta_tensor=.false.
  logical :: lrhs_term=.false., lrhs_term2=.false.
  logical :: lqpcurrent=.false., lNEMPI_correction=.true.
  logical :: lNEMPI_correction_qp_profile=.true.
  logical :: lqp_aniso_factor=.false.
  logical :: lread_alpha_tensor_z=.false., lread_alpha_tensor_z_as_y=.false.
  logical :: lread_eta_tensor_z=.false., lread_eta_tensor_z_as_y=.false.
!
  namelist /magn_mf_run_pars/ &
      alpha_effect, alpha_quenching, alpha_rmax, alpha_exp, alpha_zz, &
      gamma_effect, gamma_quenching, &
      alpha_eps, alpha_pom0, alpha_width, alpha_width2, alpha_aniso, &
      alpha_tensor, eta_tensor, &
      lalpha_profile_total, lmeanfield_noalpm, alpha_profile, &
      chit_quenching, chi_t0, lqp_profile, lqpx_profile, qp_width, qpx_width, &
      x_surface, x_surface2, z_surface, &
      alpha_rmin, kx_alpha, &
      qp_model,&
      ldelta_profile, delta_effect, delta_profile, &
      meanfield_etat, meanfield_etat_height, meanfield_etat_profile, &
      meanfield_etat_width, meanfield_etat_exp, meanfield_Beq_width, &
      meanfield_kf, meanfield_kf_profile, &
      meanfield_kf_width, meanfield_kf_width2, &
      meanfield_Beq, meanfield_Beq_height, meanfield_Beq2_height, &
      meanfield_Beq_profile, uturb, &
      meanfield_delta_height, meanfield_delta_width, &
      lmeanfield_pumping, meanfield_pumping, &
      lmeanfield_jxb, lmeanfield_jxb_with_vA2, &
      lmeanfield_chitB, lchit_with_glnTT, lrho_chit, lchit_Bext2_equil, &
      lturb_temp_diff, lignore_gradB2_inchiB, &
      meanfield_qs, meanfield_qp, meanfield_qe, meanfield_qa, &
      meanfield_Bs, meanfield_Bp, meanfield_Be, meanfield_Ba, &
      lqpcurrent,mf_qJ2, lNEMPI_correction, lNEMPI_correction_qp_profile, &
      lqp_aniso_factor, qp_aniso_factor, &
      alpha_equator, alpha_equator_gap, alpha_gap_step, &
      alpha_cutoff_up, alpha_cutoff_down, &
      lalpha_Omega_approx, lOmega_effect, Omega_profile, Omega_ampl, &
      llarge_scale_velocity, EMF_profile, lEMF_profile, &
      lrhs_term, lrhs_term2, rhs_term_amplz, rhs_term_amplphi, rhs_term_ampl, &
      Omega_rmax, Omega_rwidth, lread_alpha_tensor_z, lread_eta_tensor_z, &
      lread_alpha_tensor_z_as_y, lread_eta_tensor_z_as_y, &
      x1_alp, x2_alp, y1_alp, y2_alp
!
! Diagnostic variables (need to be consistent with reset list below)
!
  integer :: idiag_qsm=0        ! DIAG_DOC: $\left<q_p(\overline{B})\right>$
  integer :: idiag_qpm=0        ! DIAG_DOC: $\left<q_p(\overline{B})\right>$
  integer :: idiag_qem=0        ! DIAG_DOC: $\left<q_e(\overline{B})\right>$,
                                ! DIAG_DOC: in the paper referred to as
                                ! DIAG_DOC: $\left<q_g(\overline{B})\right>$
  integer :: idiag_qam=0        ! DIAG_DOC: $\left<q_a(\overline{B})\right>$
  integer :: idiag_alpm=0       ! DIAG_DOC: $\left<\alpha\right>$
  integer :: idiag_etatm=0      ! DIAG_DOC: $\left<\eta_{\rm t}\right>$
  integer :: idiag_EMFmz1=0     ! DIAG_DOC: $\left<{\cal E}\right>_{xy}|_x$
  integer :: idiag_EMFmz2=0     ! DIAG_DOC: $\left<{\cal E}\right>_{xy}|_y$
  integer :: idiag_EMFmz3=0     ! DIAG_DOC: $\left<{\cal E}\right>_{xy}|_z$
  integer :: idiag_EMFdotBm=0   ! DIAG_DOC: $\left<{\cal E}\cdot\Bv \right>$
  integer :: idiag_EMFdotB_int=0! DIAG_DOC: $\int{\cal E}\cdot\Bv dV$
  integer :: idiag_peffmxz=0    ! YAVG_DOC: $\left<{\cal P}_{\rm eff}\right>_{y}$
  integer :: idiag_alpmxz=0     ! YAVG_DOC: $\left<\alpha\right>_{y}$
!
! xy averaged diagnostics given in xyaver.in
!
  integer :: idiag_qpmz=0       ! XYAVG_DOC: $\left<q_p\right>_{xy}$
!
  contains
!***********************************************************************
    subroutine register_magn_mf()
!
!  No additional variables to be initialized here, but other mean-field
!  modules that are being called from here may need additional variables.
!
!  6-jan-11/axel: adapted from magnetic
!
!  Identify version number.
!
      if (lroot) call svn_id( &
          "$Id$")
!
!  Register secondary mean-field modules.
!
      if (lmagn_mf_demfdt) call register_magn_mf_demfdt()
!
    endsubroutine register_magn_mf
!***********************************************************************
    subroutine initialize_magn_mf(f)
!
!  Perform any post-parameter-read initialization
!
!  20-may-03/axel: reinitialize_aa added
!
      use Sub, only: erfunc
      use Mpicomm, only: stop_it
      use SharedVariables, only: put_shared_variable, get_shared_variable
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: kf_x_tmp, kf_x1_tmp, prof_tmp
!      character (len=linelen) :: dummy
      integer :: ierr, i, j
!
!  check for alpha profile
!
      if (alpha_profile=='read') then
        print*,'read alpha profile'
        open(1,file='alpha_input.dat',form='unformatted')
        read(1) alpha_input
        close(1)
      endif
!
!  check for (simple-minded) anisotropy
!
      if (any(alpha_aniso/=0.)) then
        lalpha_aniso=.true.
      else
        lalpha_aniso=.false.
      endif
      if (lroot) print*,'lalpha_aniso=',lalpha_aniso
!
!  check for (slightly less-simple minded) anisotropy
!  Put alpha_zz=1 to get alpha_ij = delta_ij - OiOj.
!
      if (any(alpha_tensor/=0.) .or. alpha_zz/=0.) then
        lalpha_tensor=.true.
        if (alpha_zz /= 0.) then
          alpha_tensor(1,1)=1.
          alpha_tensor(2,2)=1.
          alpha_tensor(3,3)=1.
          alpha_tensor_y(:,1,1)=-alpha_zz*cos(y(m1:m2))**2
          alpha_tensor_y(:,2,2)=-alpha_zz*sin(y(m1:m2))**2
          alpha_tensor_y(:,1,2)=+alpha_zz*sin(y(m1:m2))*cos(y(m1:m2))
          alpha_tensor_y(:,2,1)=alpha_tensor_y(:,1,2)
        elseif (lread_alpha_tensor_z.or.lread_alpha_tensor_z_as_y) then
!
!  Read in alpha tensor (z-dependent, and only 2x2)
!
          open(1,file='alpha_tensor_z.dat',form='unformatted')
          read(1) alpha_tensor_z
          close(1)
!
!  Each component can be scaled separately.
!  Conversely, they must all be turned on separately.
!
          do i=1,2; do j=1,2
            alpha_tensor_z(:,i,j)=alpha_tensor(i,j)*alpha_tensor_z(:,i,j)
          enddo; enddo
!
!  if lread_alpha_tensor_z_as_y: swap back (x,y,z) <-- (y,z,x)
!
          if (lread_alpha_tensor_z_as_y) then
            if (ny==nz) then
              alpha_tensor_y(n1:n2,3,3)=alpha_tensor_z(:,1,1)
              alpha_tensor_y(n1:n2,3,1)=alpha_tensor_z(:,1,2)
              alpha_tensor_y(n1:n2,1,3)=alpha_tensor_z(:,2,1)
              alpha_tensor_y(n1:n2,1,1)=alpha_tensor_z(:,2,2)
            else
              call fatal_error('initialize_magn_mf', &
                'lread_alpha_tensor_z_as_y works only when ny=nz')
            endif
          endif
        endif
      else
        lalpha_tensor=.false.
      endif
      if (lroot) print*,'lalpha_tensor=',lalpha_tensor
!
!  check for (slightly less-simple minded) anisotropy
!  Put alpha_zz=1 to get alpha_ij = delta_ij - OiOj.
!
      if (any(eta_tensor/=0.)) then
        leta_tensor=.true.
        if (lread_eta_tensor_z.or.lread_eta_tensor_z_as_y) then
!
!  Read in eta tensor (z-dependent, and only 2x2)
!
          open(1,file='eta_tensor_z.dat',form='unformatted')
          read(1) eta_tensor_z
          close(1)
!
!  Each component can be scaled separately.
!  Conversely, they must all be turned on separately.
!
          do i=1,2; do j=1,2
            eta_tensor_z(:,i,j)=eta_tensor(i,j)*eta_tensor_z(:,i,j)
          enddo; enddo
!
!  if lread_alpha_tensor_z_as_y: swap back (x,y,z) <-- (y,z,x)
!
          if (lread_eta_tensor_z_as_y) then
            if (ny==nz) then
              eta_tensor_y(n1:n2,3,3)=eta_tensor_z(:,1,1)
              eta_tensor_y(n1:n2,3,1)=eta_tensor_z(:,1,2)
              eta_tensor_y(n1:n2,1,3)=eta_tensor_z(:,2,1)
              eta_tensor_y(n1:n2,1,1)=eta_tensor_z(:,2,2)
            else
              call fatal_error('initialize_magn_mf', &
                'lread_eta_tensor_z_as_y works only when ny=nz')
            endif
          endif
        endif
      else
        leta_tensor=.false.
      endif
      if (lroot) print*,'leta_tensor=',leta_tensor
!
!  write profile (uncomment for debugging)
!
!     if (lroot) then
!       do n=n1,n2
!         print*,z(n),eta_z(n)
!       enddo
!     endif
!
!  Possibility of adding a coskx term to the rhs of the dAz/dt equation.
!
      if (lrhs_term) then
        !rhs_term=(rhs_term_ampl/rhs_term_kx)*cos(rhs_term_kx*x(l1:l2))
        !rhs_term=rhs_term_ampl*(.25*x(l1:l2)**2-.0625*x(l1:l2)**4)
        rhs_term=rhs_term_ampl*(1.-x(l1:l2)**2)
        !rhs_term=rhs_term_ampl
      endif
!
      if (lrhs_term2) then
        !rhs_termz=rhs_term_amplz*(-x(l1:l2)**2)/4.      !kk
        !rhs_termy=rhs_term_amplphi*x(l1:l2)/2.          !kk
        rhs_termy=.50*rhs_term_amplphi*(1.-x(l1:l2))*x(l1:l2)
        rhs_termz=.25*rhs_term_amplz*(.75-(1.-.25*x(l1:l2)**2)*x(l1:l2)**2)
      endif
!
!  if meanfield theory is invoked, we want to send meanfield_etat to
!  other subroutines
!
      !if (lmagn_mf .and. lrun) then
      if (lrun) then
        call put_shared_variable('meanfield_etat',meanfield_etat,ierr)
        if (ierr/=0) call fatal_error('initialize_magn_mf', &
            'there was a problem when putting meanfield_etat')
!
!  Quenching parameters: lmeanfield_chitB and chi_t0.
!
        call put_shared_variable('lmeanfield_chitB',lmeanfield_chitB,ierr)
        if (ierr/=0) call fatal_error('initialize_magn_mf', &
            'there was a problem when putting lmeanfield_chitB')
        if (lmeanfield_chitB) then
        call put_shared_variable('chi_t0',chi_t0,ierr)
        if (ierr/=0) call fatal_error('initialize_magn_mf', &
            'there was a problem when putting chi_t0')
        endif
!
!  chit_quenching for flux boundary condition
!
        !dummy=meanfield_Beq_profile
        !call put_shared_variable('meanfield_Beq_profile',dummy,ierr)
        !if (ierr/=0) call stop_it("meanfield_Beq_profile: "//&
        !     "there was a problem when putting meanfield_Beq_profile")
        call put_shared_variable('meanfield_Beq',meanfield_Beq,ierr)
        if (ierr/=0) call stop_it("meanfield_Beq: "//&
             "there was a problem when putting meanfield_Beq")
        call put_shared_variable('chit_quenching',chit_quenching,ierr)
        if (ierr/=0) call stop_it("chit_quenching: "//&
             "there was a problem when putting chit_quenching")
        call put_shared_variable('uturb',uturb,ierr)
        if (ierr/=0) call stop_it("uturb: "//&
             "there was a problem when putting uturb")
!
      endif
!
!  Compute etat profile and share with other routines.
!  Here we also set the etat_i and getat_i profiles.
!
      if (meanfield_kf/=0.0) then
        select case (meanfield_kf_profile)
        case ('const')
          kf_x1=1./meanfield_kf
          kf_x=meanfield_kf
          kf_y=1.
          kf_z=1.
        case ('sqrt_surface_x1')
          prof_tmp=.5*(1.+erfunc((x(l1:l2)-x_surface2)/meanfield_kf_width2))
          kf_x1_tmp=max(1.-x(l1:l2)/x_surface,tini)/meanfield_kf
          kf_x1=prof_tmp*kf_x1_tmp
          kf_x_tmp=1./kf_x1_tmp
          kf_x=prof_tmp*kf_x_tmp*exp(-(meanfield_kf_width*kf_x_tmp)**2)
          kf_y=1.
          kf_z=1.
        case default;
          call inevitably_fatal_error('initialize_magnetic', &
          'no such meanfield_kf_profile profile')
        endselect
      endif
!
!  Compute etat profile and share with other routines.
!  Here we also set the etat_i and getat_i profiles.
!
      if (meanfield_etat/=0.0) then
        select case (meanfield_etat_profile)
        case ('const')
          etat_x=1.
          etat_y=1.
          etat_z=meanfield_etat
          detat_x=0.
          detat_y=0.
          detat_z=0.
        case ('exp(z/H)')
          etat_x=1.
          etat_y=meanfield_etat
          etat_z=exp(z/meanfield_etat_height)
          detat_x=0.
          detat_y=0.
          detat_z=etat_z/meanfield_etat_height
        case ('surface_x1'); etat_x=0.5 &
          *(1.+erfunc((x(l1:l2)-x_surface2)/meanfield_etat_width))
          etat_y=meanfield_etat
          etat_z=1.
          detat_x=1./(meanfield_etat_width*sqrtpi) &
            *exp(-((x(l1:l2)-x_surface2)/meanfield_etat_width)**2)
          detat_y=0.
          detat_z=0.
        case ('surface_x2'); etat_x=0.25 &
          *(1.-erfunc((x(l1:l2)-x_surface)/meanfield_etat_width)) &
          *(1.+erfunc((x(l1:l2)-x_surface2)/meanfield_etat_width))
          etat_y=meanfield_etat
          etat_z=1.
          detat_x=.5/(meanfield_etat_width*sqrtpi)*( &
             exp(-((x(l1:l2)-x_surface)/meanfield_etat_width)**2) &
            *(1.+erfunc((x(l1:l2)-x_surface2)/meanfield_etat_width)) &
            +(1.+erfunc((x(l1:l2)-x_surface)/meanfield_etat_width)) &
            *exp(-((x(l1:l2)-x_surface2)/meanfield_etat_width)**2))
          detat_y=0.
          detat_z=0.
        case ('sqrt_surface_x1')
          etat_x=sqrt(kf_x1)
          etat_y=meanfield_etat
          etat_z=1.
          if (ip<=6) print*,'etat(sqrt_surface_x1)=',etat_x
          detat_x=1./(meanfield_etat_width*sqrtpi) &
            *exp(-((x(l1:l2)-x_surface2)/meanfield_etat_width)**2)
          detat_y=0.
          detat_z=0.
        case ('sin2y')
          etat_x=meanfield_etat
          etat_y=sin(y)**2
          etat_z=1.
          detat_x=0.
          detat_y=2.*sin(y)*cos(y)
          detat_z=0.
        case ('sin4y')
          etat_x=meanfield_etat
          etat_y=sin(y)**4
          etat_z=1.
          detat_x=0.
          detat_y=4.*sin(y)**3*cos(y)
          detat_z=0.
!
!  general expression, includes previous cases with meanfield_etat_exp=2 and 4
!
        case ('siny**n')
          etat_x=meanfield_etat
          etat_y=sin(y)**meanfield_etat_exp
          etat_z=1.
          detat_x=0.
          detat_y=meanfield_etat_exp*sin(y)**(meanfield_etat_exp-1.)*cos(y)
          detat_z=0.
        case ('Jouve-2008-benchmark')
          etat_x = meanfield_etat*(0.01 + 0.5*(1.-0.01)*(1.0+erfunc((x(l1:l2)-0.7)/0.02)))
          etat_y = 1.
          etat_z = 1.
          detat_x= meanfield_etat*0.5*(1.-0.01)*exp(-((x(l1:l2)-0.7)/0.02)**2)
          detat_y= 0.
          detat_z= 0.
        case ('sinx**2*exp(-z/H)')
          etat_x=sin(x(l1:l2))**2
          etat_y=meanfield_etat
          etat_z=exp(-z/meanfield_etat_height)
          detat_x=2.*sin(x(l1:l2))*cos(x(l1:l2))
          detat_y=0.
          detat_z=-etat_z/meanfield_etat_height
        case ('surface_z'); etat_z=0.5 &
          *(1.-erfunc((z-z_surface)/meanfield_etat_width))
          etat_y=meanfield_etat
          etat_x=1.
          detat_z=-1./(meanfield_etat_width*sqrtpi) &
            *exp(-((z-z_surface)/meanfield_etat_width)**2)
          detat_y=0.
          detat_x=0.
        case default;
          call inevitably_fatal_error('initialize_magnetic', &
          'no such meanfield_etat_profile profile')
        endselect
!
!  share etat profile with viscosity module
!
        if (lviscosity) then
          call put_shared_variable('etat_x',etat_x,ierr)
          call put_shared_variable('etat_y',etat_y,ierr)
          call put_shared_variable('etat_z',etat_z,ierr)
          call put_shared_variable('detat_x',detat_x,ierr)
          call put_shared_variable('detat_y',detat_y,ierr)
          call put_shared_variable('detat_z',detat_z,ierr)
          print*,'ipz,z(n),etat_z(n),detat_z(n)'
          do n=n1,n2
            print*,ipz,z(n),etat_z(n),detat_z(n)
          enddo
          print*
        endif
      endif
!
!  define inverse of meanfield_qp
!
      if (meanfield_qp==0.) then
        meanfield_qp1=0.
      else
        meanfield_qp1=1./meanfield_qp
      endif
!
!  define inverse of meanfield_qs
!
      if (meanfield_qs==0.) then
        meanfield_qs1=0.
      else
        meanfield_qs1=1./meanfield_qs
      endif
!
!  define inverse of meanfield_qe
!
      if (meanfield_qe==0.) then
        meanfield_qe1=0.
      else
        meanfield_qe1=1./meanfield_qe
      endif
!
!  define inverse of meanfield_qa
!
      if (meanfield_qa==0.) then
        meanfield_qa1=0.
      else
        meanfield_qa1=1./meanfield_qa
      endif
!
!  define meanfield_qp_profile
!
      if (lqp_profile) then
        qp_profile=0.5*(1.-erfunc((z-z_surface)/qp_width))
        !qe_profile=0.5*(1.-erfunc((z-z_surface)/qe_width))
        qp_profder=-exp(-((z-z_surface)/qp_width)**2)/(qp_width*sqrtpi)
        !qe_profder=-exp(-((z-z_surface)/qe_width)**2)/(qe_width*sqrtpi)
      else
        qp_profile=1.
        qp_profder=0.
      endif
!
!  define meanfield_qpx_profile (x direction)
!
      if (lqpx_profile) then
        qpx_profile=exp(-(x(l1:l2)/qpx_width)**2)
        qpx_profder=-2.*x(l1:l2)/qpx_width**2*exp(-(x(l1:l2)/qpx_width)**2)
      else
        qpx_profile=1.
        qpx_profder=0.
      endif
!
!  Initialize module variables which are parameter dependent
!  wave speed of gauge potential
!
      if (lrun) then
        call get_shared_variable('lweyl_gauge',lweyl_gauge,ierr)
        if (ierr/=0) &
            call fatal_error("initialize_special: ", "cannot get lweyl_gauge")
        if (lroot) print*,'initialize_magn_mf: lweyl_gauge=',lweyl_gauge
!       if (.not.lweyl_gauge) then
!         call get_shared_variable('eta',eta,ierr)
!         if (ierr/=0) &
!             call fatal_error("initialize_special: ", "cannot get shared eta")
!       endif
      endif
!
!  Initialize secondary mean-field modules:
!
      if (lmagn_mf_demfdt .or. lalpm .or. lalpm_alternate ) then
        call put_shared_variable('kf_x',kf_x,ierr)
        call put_shared_variable('kf_y',kf_y,ierr)
        call put_shared_variable('kf_z',kf_z,ierr)
        call put_shared_variable('kf_x1',kf_x1,ierr)
        call put_shared_variable('etat_x',etat_x,ierr)
        call put_shared_variable('etat_y',etat_y,ierr)
        call put_shared_variable('etat_z',etat_z,ierr)
        if (ierr/=0) call fatal_error('initialize_magnetic',&
            'there was a problem when sharing etat_xyz')
        call initialize_magn_mf_demfdt(f)
      endif
!
!  Get B_ext2 from magnetic module.
!
      call get_shared_variable('B_ext2',B_ext2,ierr)
      if (ierr/=0) &
          call fatal_error("initialize_magn_mf:", "cannot get B_ext2")
!
    endsubroutine initialize_magn_mf
!***********************************************************************
    subroutine init_aa_mf(f)
!
!  Initialise mean-field related magnetic field; called from magnetic.f90
!  At the moment, no own initial conditions are allowed, but we need this
!  to call secondary modules
!
!   6-jan-2011/axel: adapted from magnetic
!
      real, dimension (mx,my,mz,mfarray) :: f
!
!  Initialize secondary mean-field modules.
!
      if (lmagn_mf_demfdt) call init_aa_mf_demfdt(f)
!
    endsubroutine init_aa_mf
!***********************************************************************
    subroutine pencil_criteria_magn_mf()
!
!   All pencils that the magnetic mean-field module depends on
!   are specified here.
!
!  28-jul-10/axel: adapted from magnetic
!
      lpenc_requested(i_bb)=.true.
!
      if (meanfield_etat/=0.0.or.ietat/=0) &
          lpenc_requested(i_del2a)=.true.
!
!  In mean-field theory, with variable etat, need divA for resistive gauge.
!
      if (meanfield_etat_profile/='const') then
        lpenc_requested(i_diva)=.true.
      endif
!
!  In spherical coordinates, we also need grad(divA)
!
      if (lspherical_coords) then
        lpenc_requested(i_graddiva)=.true.
        lpenc_requested(i_x_mn)=lOmega_effect
      endif
!
!  For mean-field modelling in momentum equation:
!
      if (lmeanfield_jxb) then
        lpenc_requested(i_b2)=.true.
        lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_jxbr_mf)=.true.
        if (lmeanfield_jxb_with_vA2) lpenc_requested(i_va2)=.true.
      endif
!
!  For mean-field modelling in entropy equation:
!
      if (lmeanfield_chitB) then
        lpenc_requested(i_b2)=.true.
        lpenc_requested(i_glnrho)=.true.
        lpenc_requested(i_chiB_mf)=.true.
        if (lrho_chit) lpenc_requested(i_rho1)=.true.
        if (lturb_temp_diff) then
          lpenc_requested(i_cp)=.true.
          lpenc_requested(i_glnTT)=.true.
          lpenc_requested(i_del2lnTT)=.true.
        else
          lpenc_requested(i_gss)=.true.
          lpenc_requested(i_del2ss)=.true.
        endif
      endif
!
!  Turbulent diffusivity
!
      if (meanfield_etat/=0.0 .or. ietat/=0 .or. &
          alpha_effect/=0.0 .or. delta_effect/=0.0 .or. &
          gamma_effect/=0.0 .or. lread_alpha_tensor_z .or. &
          lread_eta_tensor_z .or. lread_alpha_tensor_z_as_y .or. &
          lread_eta_tensor_z_as_y) &
          lpenc_requested(i_mf_EMF)=.true.
      if (delta_effect/=0.0) lpenc_requested(i_oxJ)=.true.
!
!  J is needed when the eta tensor is invoked.
!  (But there should be other such requests that have apparently
!  not yet caused pencil checks to fail yet...)
!
      if (leta_tensor) then
        lpenc_requested(i_jj)=.true.
      endif
!
      if (idiag_EMFdotBm/=0.or.idiag_EMFdotB_int/=0) lpenc_diagnos(i_mf_EMFdotB)=.true.
!
! If large_scale_velocity is included we need this pencil
!
      if (llarge_scale_velocity) lpenc_requested(i_uxb)=.true.
!
!  Pencil criteria for secondary modules
!
      if (lmagn_mf_demfdt) call pencil_criteria_magn_mf_demfdt()
!
!  Nempi model C requires currents
!
      if (lqpcurrent) then
        lpenc_requested(i_jj)=.true.
        lpenc_requested(i_j2)=.true.
        lpenc_requested(i_jij)=.true.
      endif
!
      if (lNEMPI_correction) then
        lpenc_requested(i_glnrho)=.true.
        lpenc_requested(i_rho1)=.true.
        lpenc_requested(i_b2)=.true.
      endif
!
    endsubroutine pencil_criteria_magn_mf
!***********************************************************************
    subroutine pencil_interdep_magn_mf(lpencil_in)
!
!  Interdependency among pencils from the Magnetic module is specified here.
!
!  28-jul-10/axel: adapted from magnetic
!
      logical, dimension(npencils) :: lpencil_in
!
!  exa
!
      if (lpencil_in(i_exa)) then
        lpencil_in(i_aa)=.true.
        lpencil_in(i_mf_EMF)=.true.
      endif
!
!  mf_EMFdotB
!
      if (lpencil_in(i_mf_EMFdotB)) then
        lpencil_in(i_mf_EMF)=.true.
        lpencil_in(i_bb)=.true.
      endif
!
!  mf_EMFdotB
!
      if (lpencil_in(i_mf_EMF)) then
        if (lspherical_coords) then
          lpencil_in(i_jj)=.true.
          lpencil_in(i_graddivA)=.true.
        endif
        lpencil_in(i_b2)=.true.
!
!  compute del2a regardless of gauge
!
        if (meanfield_etat/=0.0 .or. ietat/=0) then
!         if (lweyl_gauge) then
!           lpencil_in(i_jj)=.true.
!         else
            lpencil_in(i_del2a)=.true.
!         endif
        endif
      endif
!
!  oxJ effect
!
        if (lpencil_in(i_oxJ)) lpencil_in(i_jj)=.true.
!
!  ??
!
!     if (lpencil_in(i_del2A)) then
!       if (lspherical_coords) then
!       endif
!     endif
!
!  Mean-field Lorentz force: jxb_mf
!
      if (lpencil_in(i_jxbr_mf)) lpencil_in(i_jxb_mf)=.true.
      if (lpencil_in(i_jxb_mf)) lpencil_in(i_jxb)=.true.
!
!  Pencil criteria for secondary modules
!
      if (lmagn_mf_demfdt) call pencil_interdep_magn_mf_demfdt(lpencil_in)
!
    endsubroutine pencil_interdep_magn_mf
!***********************************************************************
    subroutine calc_pencils_magn_mf(f,p)
!
!  Calculate Magnetic mean-field pencils.
!  Most basic pencils should come first, as others may depend on them.
!
!  28-jul-10/axel: adapted from magnetic
!
      use Sub
      use General, only: bessj
      use Diagnostics, only: sum_mn_name, xysum_mn_name_z, ysum_mn_name_xz
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
!
      real, dimension (nx) :: alpha_total, rr, pp, pp0, spiral
      real, dimension (nx) :: meanfield_etat_tmp
      real, dimension (nx) :: alpha_tmp, alpha_quenching_tmp, delta_tmp
      real, dimension (nx) :: kf_tmp, EMF_prof, alpm, prefact, Beq21
      real, dimension (nx) :: meanfield_qs_func, meanfield_qp_func, meanfield_qa_func
      real, dimension (nx) :: meanfield_qe_func, meanfield_qs_der, meanfield_qa_der
      real, dimension (nx) :: meanfield_qp_der, meanfield_qe_der, BiBk_Bki
      real, dimension (nx) :: meanfield_Bs21, meanfield_Bp21, meanfield_Be21, meanfield_Ba21
      real, dimension (nx) :: meanfield_etaB2, g2, chit_prof !, quench_chiB
      real, dimension (nx) :: oneQbeta02, oneQbeta2, XXj_BBj, BjBzj
      real, dimension (nx,3) :: Bk_Bki, exa_meanfield, glnchit_prof, glnchit, XXj
      real, dimension (nx,3) :: meanfield_getat_tmp, getat_cross_B_tmp, B2glnrho, glnchit2
      real :: kx,fact
      integer :: i,j,l
!
      intent(inout) :: f,p
!
      save :: chit_prof, glnchit_prof
!
!  mean-field Lorentz force
!
      if (lpencil(i_jxbr_mf)) then
!
!  The following 9 lines have not been used for any publications so far.
!
!        if (lmeanfield_jxb_with_vA2) then
!          call inevitably_fatal_error('calc_pencils_magnetic', &
!            'lmeanfield_jxb_with_vA2 should be checked')
!--
!         meanfield_urms21=1./(3.*meanfield_kf*meanfield_etat)**2
!         meanfield_qs_func=meanfield_qs*(1.-2*pi_1*atan(p%vA2*meanfield_urms21))
!         meanfield_qp_func=meanfield_qp*(1.-2*pi_1*atan(p%vA2*meanfield_urms21))
!         meanfield_qe_func=meanfield_qe*(1.-2*pi_1*atan(p%b2*meanfield_Be21))
!         meanfield_qs_der=2*pi_1*meanfield_qs/(1.+(p%vA2*meanfield_urms21)**2)
!         meanfield_qp_der=2*pi_1*meanfield_qp/(1.+(p%vA2*meanfield_urms21)**2)
!         meanfield_qe_der=2*pi_1*meanfield_qe*meanfield_Be21/(1.+(p%b2*meanfield_Be21)**2)
!         call multsv_mn(meanfield_qs_func,p%jxb,tmp_jxb); p%jxb_mf=tmp_jxb
!--
!
!  The following (not lmeanfield_jxb_with_vA2) has been used for the
!  various publications so far.
!
!        else
          select case (meanfield_Beq_profile)
          case ('exp(z/H)');
!           H = -2* density scale height, old version
            Beq21=1./meanfield_Beq**2*exp(-2*z(n)/meanfield_Beq_height)
!
!  Beq ~ exp(-z/2H), where H=Hrho is the density scale height.
!  Give here meanfield_Beq2_height, which is the scale height for Beq^2.
!
          case ('exp(z/2H)');
!           H = density scale height, more straightforward
            Beq21=1./meanfield_Beq**2*exp(z(n)/meanfield_Beq2_height)
!
!  Beq profile for constant uturb
!
          case ('uturbconst');
!
!  allows to take into account a shifted equilibrium solution
!  caused by the additional pressure contribution
!
            Beq21=mu01*p%rho1/(uturb**2)
!
!  Beq21 becomes large in the corona
!
          case ('uturb_surface');
            Beq21=mu01*p%rho1/(uturb**2)* &
              0.5*(1.-erfunc((z(n)-z_surface)/meanfield_Beq_width))
!
!  default
!
          case default;
            Beq21=1./meanfield_Beq**2
          endselect
          p%mf_Beq21=Beq21
!
!  Here, p%b2*meanfield_Bp21 = beta^2/beta_p^2, etc.
!
          meanfield_Bs21=Beq21/meanfield_Bs**2
          meanfield_Bp21=Beq21/meanfield_Bp**2
          meanfield_Be21=Beq21/meanfield_Be**2
          meanfield_Ba21=Beq21/meanfield_Ba**2
!
!  Compute qp(B^2), qs(B^2), qe(B^2), and their derivatives for 2 different parameterizations
!  qp=qp0/(1+B^2*meanfield_Bp21)
!  qp'=-qp0/(1+B^2*meanfield_Bp21)^2*meanfield_Bp21=-qp^2*meanfield_Bp21/qp0
!
          if(qp_model=='rational') then
            meanfield_qp_func=meanfield_qp/(1.+p%b2*meanfield_Bp21)*qp_profile(n)*qpx_profile
            meanfield_qs_func=meanfield_qs/(1.+p%b2*meanfield_Bs21)
            meanfield_qe_func=meanfield_qe/(1.+p%b2*meanfield_Be21)
            meanfield_qa_func=meanfield_qa/(1.+p%b2*meanfield_Ba21)
            meanfield_qp_der=-meanfield_qp_func**2*meanfield_Bp21*meanfield_qp1
            meanfield_qs_der=-meanfield_qs_func**2*meanfield_Bs21*meanfield_qs1
            meanfield_qe_der=-meanfield_qe_func**2*meanfield_Be21*meanfield_qe1
            meanfield_qa_der=-meanfield_qa_func**2*meanfield_Ba21*meanfield_qa1
          else
            meanfield_qp_func=meanfield_qp*(1.-2*pi_1*atan(p%b2*meanfield_Bp21))*qp_profile(n)
            meanfield_qs_func=meanfield_qs*(1.-2*pi_1*atan(p%b2*meanfield_Bs21))
            meanfield_qe_func=meanfield_qe*(1.-2*pi_1*atan(p%b2*meanfield_Be21))
            meanfield_qp_der=-2*pi_1*meanfield_qp*meanfield_Bp21/(1.+(p%b2*meanfield_Bp21)**2)
            meanfield_qs_der=-2*pi_1*meanfield_qs*meanfield_Bs21/(1.+(p%b2*meanfield_Bs21)**2)
            meanfield_qe_der=-2*pi_1*meanfield_qe*meanfield_Be21/(1.+(p%b2*meanfield_Be21)**2)
          endif
          p%mf_qp=meanfield_qp_func
!
!  Add (1/2)*grad[qp*B^2]. This initializes p%jxb_mf.
!  Note: p%jij is not J_i,j; omit extra term for the time being.
!  Also: "p%b2*meanfield_qp_der" is the same as dqp/dlnbeta^2.
!
          call multmv_transp(p%bij,p%bb,Bk_Bki) !=1/2 grad B^2
          if (lqpcurrent) then
            call multsv_mn(mu01*(1+p%j2*mf_qJ2)* &
              (meanfield_qp_func+p%b2*meanfield_qp_der), &
              Bk_Bki,p%jxb_mf)
          else
            call multsv_mn(mu01*(meanfield_qp_func+p%b2*meanfield_qp_der),Bk_Bki,p%jxb_mf)
            if (lNEMPI_correction) then
              call multsv_mn_add(-.5*mu01*p%b2**2*meanfield_qp_der,p%glnrho,p%jxb_mf)
            endif
!
!  Allow for qp profile
!
            if (lNEMPI_correction_qp_profile) then
              p%jxb_mf(:,3)=p%jxb_mf(:,3)+.5*mu01*qp_profder(n)*meanfield_qp_func*p%b2
            endif
          endif
!
!  Allow for simple-minded anisotropy
!
          if (lqp_aniso_factor) p%jxb_mf(:,3)=p%jxb_mf(:,3)*qp_aniso_factor
!
!  Add -B.grad[qs*B_i]. This term does not promote instability.
!
          call multsv_mn_add(-meanfield_qs_func,p%jxb+mu01*Bk_Bki,p%jxb_mf)
          call dot(Bk_Bki,p%bb,BiBk_Bki)
          call multsv_mn_add(-2.*mu01*meanfield_qs_der*BiBk_Bki,p%bb,p%jxb_mf)
!
!  Add e_z*grad(qe*B^2). This has not yet been found to promote instability.
!  Added below (in commented out form, how I think it should be if used)
!
          p%jxb_mf(:,3)=p%jxb_mf(:,3)+2.*mu01*(meanfield_qe_der*p%b2+meanfield_qe_func)*Bk_Bki(:,3)
          !p%jxb_mf(:,3)=p%jxb_mf(:,3)+mu01*(meanfield_qe_der*p%b2+meanfield_qe_func)*Bk_Bki(:,3) &
          !                        -.5*mu01*p%b2**2*meanfield_qe_der,p%glnrho(:,3) &
          !                        +.5*mu01*qe_profder(n)*meanfield_qe_func*p%b2
!
!  Add div[qa*Bz*(Bi*gj+Bj*gi)]
!  Begin by initializing auxiliary vector X_j (see notes)
!
          XXj=2.*Bk_Bki
          call multsv_mn_add(-mu01*p%rho1*p%b2,p%grho,XXj)
!
!  Do dqa/dB2 * Bz * (BiX3+z3B.X) term
!
          call dot(XXj,p%bb,XXj_BBj)
          p%jxb_mf(:,3)=p%jxb_mf(:,3)+meanfield_qa_der*p%bb(:,3)*XXj_BBj
          call multsv_mn_add(meanfield_qa_der*p%bb(:,3)*XXj(:,3),p%bb,p%jxb_mf)
!
!  Do qa*[Bz,z*(Bi+Bj*Bz,j*zi)] term
!
          call dot(p%bb,p%bij(:,3,:),BjBzj)
          p%jxb_mf(:,3)=p%jxb_mf(:,3)+mu01*meanfield_qa_func*p%bij(:,3,3)*BjBzj
          call multsv_mn_add(mu01*meanfield_qa_func*p%bij(:,3,3),p%bb,p%jxb_mf)
!
!  Do qa*Bz*B_i,z
!
          call multsv_mn_add(mu01*meanfield_qa_func*p%bb(:,3),p%bij(:,:,3),p%jxb_mf)
!
!        endif
        call multsv_mn(p%rho1,p%jxb_mf,p%jxbr_mf)
      endif
!
!  compute alpha effect for EMF
!
!  mf_EMF for alpha effect dynamos
!
      if (lpencil(i_mf_EMF)) then
!
!  compute alpha profile (alpha_tmp)
!
        if (nxgrid/=1) kx=2*pi/Lx
        select case (alpha_profile)
        case ('const'); alpha_tmp=1.
        case ('box')
          if (y(m) >= y1_alp .and. y(m) <= y2_alp) then
            where(x(l1:l2)>=x1_alp .and. x(l1:l2)<=x2_alp)
              alpha_tmp=1.
            elsewhere
              alpha_tmp=0.
            endwhere
          else
            alpha_tmp=0.
          endif
        case ('coskx'); alpha_tmp=sqrt2*cos(kx_alpha*x(l1:l2))
        case ('siny'); alpha_tmp=sin(y(m))
        case ('sinz'); alpha_tmp=sin(z(n))
        case ('cos(z/2)'); alpha_tmp=cos(.5*z(n))
        case ('cos(z/2)_with_halo'); alpha_tmp=max(cos(.5*z(n)),0.)
        case ('sphere')
          if (lspherical_coords) then
            rr=x(l1:l2)
          elseif (lcylindrical_coords) then
            rr=sqrt(x(l1:l2)**2+z(n)**2)
          else
            rr=sqrt(x(l1:l2)**2+y(m)**2+z(n)**2)
          endif
          alpha_tmp=.5*(1.-erfunc((rr-alpha_rmax)/alpha_width))
        case ('z')
          if (alpha_rmax/=0.) then
            rr=sqrt(x(l1:l2)**2+y(m)**2+z(n)**2)
            pp=.5*(1.-erfunc((rr-alpha_rmax)/alpha_width))
            alpha_tmp=z(n)*pp
          else
            alpha_tmp=z(n)
          endif
        case ('z/r')
          rr=sqrt(x(l1:l2)**2+y(m)**2+z(n)**2)
          if (alpha_rmax/=0.) then
            pp=.5*(1.-erfunc((rr-alpha_rmax)/alpha_width))
            alpha_tmp=z(n)*pp/rr
          else
            alpha_tmp=z(n)/rr
          endif
        case ('1-tanhz'); alpha_tmp=.5*(1.-tanh(z(n)/alpha_width))
        case ('z/H'); alpha_tmp=z(n)/xyz1(3)
        case ('z/H_0'); alpha_tmp=z(n)/xyz1(3); if (z(n)==xyz1(3)) alpha_tmp=0.
        case ('y/H'); alpha_tmp=y(m)/xyz1(3)
        case ('J0x'); do l=l1,l2; alpha_tmp(l-l1+1)=bessj(0,k1bessel0*x(l)); enddo
        case ('J0x_2nd'); do l=l1,l2; alpha_tmp(l-l1+1)=bessj(0,k2bessel0*x(l)); enddo
        case ('cosy'); alpha_tmp=cos(y(m))
        case ('cosy*sin2y'); alpha_tmp=cos(y(m))*sin(y(m))**2
        case ('cosy*sin4y'); alpha_tmp=cos(y(m))*sin(y(m))**4
        case ('cosy*sin6y'); alpha_tmp=cos(y(m))*sin(y(m))**6
        case ('cosy*sin8y'); alpha_tmp=cos(y(m))*sin(y(m))**8
        case ('cosy*siny**n'); alpha_tmp=cos(y(m))*sin(y(m))**alpha_exp
        case ('surface_x*cosy'); alpha_tmp=0.5*(1.-erfunc((x(l1:l2)-x_surface)/alpha_width))*cos(y(m))
        case ('surface_x2*cosy'); alpha_tmp=0.25 &
          *(1.-erfunc((x(l1:l2)-x_surface)/alpha_width)) &
          *(1.+erfunc((x(l1:l2)-x_surface2)/alpha_width2))*cos(y(m))
        case ('sqrt_surface_x2*cosy')
          !alpha_tmp=sqrt(kf_x)*cos(y(m))
          alpha_tmp=kf_x1*cos(y(m))
          if (ip<=6.and.m==m1) print*,'alpha(sqrt_surface_x2)=',alpha_tmp
        case ('y*(1+eps*sinx)'); alpha_tmp=y(m)*(1.+alpha_eps*sin(kx*x(l1:l2)))
        case ('step-nhemi'); alpha_tmp=-tanh((y(m)-pi/2)/alpha_gap_step)
        case ('stepy'); alpha_tmp=-tanh((y(m)-yequator)/alpha_gap_step)
        case ('stepz'); alpha_tmp=-tanh((z(n)-zequator)/alpha_gap_step)
        case ('ystep-xcutoff')
           alpha_tmp=-tanh((y(m)-pi/2)/alpha_gap_step)&
             *(1+stepdown(x(l1:l2),alpha_rmax,alpha_width))
        case ('cosy*sin**n-xcutoff')
           alpha_tmp=(cos(y(m))*sin(y(m))**alpha_exp)&
             *(1+stepdown(x(l1:l2),alpha_rmax,alpha_width))
        case ('step-drop'); alpha_tmp=(1. &
                -step(y(m),pi/2.-alpha_equator_gap,alpha_gap_step) &
                -step(y(m),pi/2.+alpha_equator_gap,alpha_gap_step) &
                -step(alpha_cutoff_up,y(m),alpha_gap_step) &
                +step(y(m),alpha_cutoff_down,alpha_gap_step))
        case ('surface_z'); alpha_tmp=0.5*(1.-erfunc((z(n)-z_surface)/alpha_width))
        case ('above_x0'); alpha_tmp=0.5*(1.+erfunc((x(l1:l2)-alpha_rmin)/alpha_width))
        case ('z/H+surface_z'); alpha_tmp=(z(n)/z_surface)*0.5*(1.-erfunc((z(n)-z_surface)/alpha_width))
          if (headtt) print*,'alpha_profile=z/H+surface_z: z_surface,alpha_width=',z_surface,alpha_width
!
!  Galactic alpha effect profile
!
        case ('z/H*exp(-z2/2h2)')
          if (lcylindrical_coords) then
            if (headtt) print*,'z/H*exp(-z2/2h2) profile (cylindrical coords)',alpha_rmax
            if (alpha_rmax/=0.) then
              rr=x(l1:l2)
              pp=.5*(1.-erfunc((rr-alpha_rmax)/alpha_width2))
              alpha_tmp=z(n)/alpha_width*exp(-.5*(z(n)/alpha_width)**2)*pp
              !alpha_tmp=exp(-.5*(z(n)/alpha_width)**2)*pp
              if (alpha_pom0/=0.) alpha_tmp=alpha_tmp/(1.+(rr/alpha_pom0)**4)**.25
            else
              alpha_tmp=z(n)/alpha_width*exp(-.5*(z(n)/alpha_width)**2)
            endif
          else
            if (headtt) print*,'z/H*exp(-z2/2h2) profile (Cartesian)',alpha_rmax
            if (alpha_rmax/=0.) then
              rr=sqrt(x(l1:l2)**2+y(m)**2)
              pp=.5*(1.-erfunc((rr-alpha_rmax)/alpha_width2))
              alpha_tmp=z(n)/alpha_width*exp(-.5*(z(n)/alpha_width)**2)*pp
              if (alpha_pom0/=0.) alpha_tmp=alpha_tmp/(1.+(rr/alpha_pom0)**4)**.25
            else
              alpha_tmp=z(n)/alpha_width*exp(-.5*(z(n)/alpha_width)**2)
            endif
          endif
!
!  Galactic alpha effect profile with spiral
!
        case ('z/H*exp(-z2/2h2)*spiral')
          rr=sqrt(x(l1:l2)**2+y(m)**2)
          pp=atan2(y(m),x(l1:l2))
          pp0=2.*alog(rr)
          spiral=sin((pp-pp0))**2
          alpha_tmp=z(n)/alpha_width*exp(-.5*(z(n)/alpha_width)**2)*spiral
        case ('z/H*erfunc(H-z)'); alpha_tmp=z(n)/xyz1(3)*erfunc((xyz1(3)-abs(z(n)))/alpha_width)
        case ('read'); alpha_tmp=alpha_input(l1:l2,m)
        case ('Jouve-2008-benchmark')
          alpha_tmp=(3.*sqrt(3.)/4.)*(sin(y(m)))**2*cos(y(m))*(1.0+erfunc((x(l1:l2)-0.7)/0.02))
        case ('nothing');
          call inevitably_fatal_error('calc_pencils_magnetic', &
            'alpha_profile="nothing" has been renamed to "const", please update your run.in')
        case default;
          call inevitably_fatal_error('calc_pencils_magnetic', &
            'alpha_profile no such alpha profile')
        endselect
!
!  compute delta effect profile (delta_tmp)
!
        select case (delta_profile)
        case ('const'); delta_tmp=1.
        case ('cos(z/2)_with_halo'); delta_tmp=max(cos(.5*z(n)),0.)
        case ('sincos(z/2)_with_halo'); delta_tmp=max(cos(.5*z(n)),0.)*sin(.5*z(n))
        case ('cos(theta)'); delta_tmp=cos(y(m))
        case ('sin(theta)'); delta_tmp=sin(y(m))
        case ('cosy*sin4y'); delta_tmp=cos(y(m))*sin(y(m))**4
        case ('tanhz'); delta_tmp=.5*(1.-tanh((abs(z(n))-meanfield_delta_height)/meanfield_delta_width))
        case default;
          call inevitably_fatal_error('calc_pencils_magnetic', &
            'delta_profile no such delta profile')
        endselect
!
!  Compute diffusion term.
!  This sets the meanfield_etat_tmp term at each time step.
!
        if (meanfield_etat/=0.0) then
          meanfield_etat_tmp=etat_x*etat_y(m)*etat_z(n)
          meanfield_getat_tmp(:,1)=detat_x*etat_y(m)*etat_z(n)
          meanfield_getat_tmp(:,2)=etat_x*detat_y(m)*etat_z(n)
          meanfield_getat_tmp(:,3)=etat_x*etat_y(m)*detat_z(n)
!
!  1/r correction for spherical symmetry (assuming axisymmetric profiles,
!  i.e.  meanfield_getat_tmp(:,3)=0.)
!
          if (lspherical_coords) then
            meanfield_getat_tmp(:,2)=r1_mn*meanfield_getat_tmp(:,2)
          endif
!
!  Magnetic etat quenching (contribution to pumping currently ignored)
!
          if (meanfield_etaB/=0.0) then
            meanfield_etaB2=meanfield_etaB**2
            meanfield_etat_tmp=meanfield_etat_tmp/sqrt(1.+p%b2/meanfield_etaB2)
!           call quench_simple(p%b2,meanfield_etaB2,meanfield_etat_tmp)
          endif
        endif
!
!  Here we initialize alpha_total.
!  Allow for the possibility of dynamical alpha.
!  By default, lmeanfield_noalpm=F, so normally treat
!  alpha_tmp as a general profile in front of the total alpha.
!  Invoke lmeanfield_noalpm to treat magnetic alpha separately
!  (which makes physically more sense!)
!
        if (lalpm .and. .not. lmeanfield_noalpm) then
          if (lalpha_profile_total) then
            alpha_total=(alpha_effect+f(l1:l2,m,n,ialpm))*alpha_tmp
            if (headtt) print*,'use alp=(alpK+alpM)*profile'
          else
            alpha_total=alpha_effect*alpha_tmp+f(l1:l2,m,n,ialpm)
            if (headtt) print*,'use alp=alpK*profile+alpM'
          endif
!
!  Possibility of *alternate* dynamical alpha.
!  Here we initialize alpha_total.
!
        elseif (lalpm_alternate .and. .not. lmeanfield_noalpm) then
          kf_tmp=kf_x*kf_y(m)*kf_z(n)
          prefact=meanfield_etat_tmp*(kf_tmp/meanfield_Beq)**2
          alpm=prefact*(f(l1:l2,m,n,ialpm)-p%ab)
          if (lalpha_profile_total) then
            alpha_total=(alpha_effect+alpm)*alpha_tmp
          else
            alpha_total=alpha_effect*alpha_tmp+alpm
          endif
        else
          alpha_total=alpha_effect*alpha_tmp
        endif
!
!  Possibility of conventional alpha quenching (rescales alpha_total).
!  Initialize EMF with alpha_total*bb.
!  Here we initialize p%mf_EMF.
!  NOTE: the following with alpha_quenching*sqrt(kf_x) was a hardwired fix
!  and and is now dealt with using meanfield_Beq_profile='sqrt(kf_x)'.
!  NOTE2: the calculation of Beq21 is here repeated, which should be avoided.
!
        if (alpha_quenching/=0.0) then
          select case (meanfield_Beq_profile)
          case ('uturbconst');
            Beq21=mu01*p%rho1/(uturb**2)
          case ('sqrt(kf_x)');
            Beq21=sqrt(kf_x)/meanfield_Beq**2
          case default;
            Beq21=1./meanfield_Beq**2
          endselect
          alpha_quenching_tmp=alpha_quenching
          alpha_total=alpha_total/(1.+alpha_quenching_tmp*p%b2*Beq21)
        endif
!
!  Apply alpha effect; allow for anisotropy
!
        if (lalpha_aniso) then
          do j=1,3
            p%mf_EMF(:,j)=alpha_total*alpha_aniso(j)*p%bb(:,j)
          enddo
        elseif (lalpha_tensor) then
          do i=1,3
            p%mf_EMF(:,i)=0.
            do j=1,3
              if (alpha_zz /= 0.) then
                p%mf_EMF(:,i)=p%mf_EMF(:,i)+alpha_total &
                  *(alpha_tensor(i,j)+alpha_tensor_y(m-m1+1,i,j))*p%bb(:,j)
              elseif (lread_alpha_tensor_z_as_y) then
                p%mf_EMF(:,i)=p%mf_EMF(:,i)+ &
                  alpha_tensor_y(m-m1+1,i,j)*p%bb(:,j)
              elseif (lread_alpha_tensor_z) then
                if (i<=2.and.j<=2) then
                  !p%mf_EMF(:,i)=p%mf_EMF(:,i)+alpha_total &
!AB: working with alpha_total is questionable in this case
                  p%mf_EMF(:,i)=p%mf_EMF(:,i)+1. &
                      *alpha_tensor_z(n-n1+1,i,j)*p%bb(:,j)
                endif
              else
                p%mf_EMF(:,i)=p%mf_EMF(:,i)+alpha_total*alpha_tensor(i,j)*p%bb(:,j)
              endif
            enddo
          enddo
        else
          call multsv_mn(alpha_total,p%bb,p%mf_EMF)
        endif
!
!  Optionally, run the alpha-Omega approximation, i.e. apply
!  alpha effect only to the toroidal component. Since p%mf_EMF
!  was initialized only in the previous line, we can just set
!  the r and theta components to zero (in spherical coordinates).
!
        if (lalpha_Omega_approx) then
          p%mf_EMF(:,1:2)=0.
          call fatal_error("calc_pencils_magn_mf: ", &
              "lalpha_Omega_approx not implemented for this case")
        endif
!
!  Apply eta tensor, but subtract part from etat for stability reasons.
!  In other words, any constant meanfield_etat should have formally no
!  effect, but is always needed because the pure Weyl gauge is unstable.
!
        if (leta_tensor) then
          if (lread_eta_tensor_z_as_y) then
            do i=1,3
              do j=1,3
                p%mf_EMF(:,i)=p%mf_EMF(:,i)-eta_tensor_y(m-m1+1,i,j)*p%jj(:,j)
              enddo
              p%mf_EMF(:,i)=p%mf_EMF(:,i)+meanfield_etat*p%jj(:,i)
            enddo
          else
            do i=1,2
              do j=1,2
                p%mf_EMF(:,i)=p%mf_EMF(:,i)-eta_tensor_z(n-n1+1,i,j)*p%jj(:,j)
              enddo
              p%mf_EMF(:,i)=p%mf_EMF(:,i)+meanfield_etat*p%jj(:,i)
            enddo
          endif
        endif
!
!  Add possible delta x J effect and turbulent diffusion to EMF.
!
        if (ldelta_profile) then
          p%mf_EMF=p%mf_EMF+delta_effect*p%oxJ
        else
          if (lspherical_coords) then
! assuming 1D in theta, delta is only delta_theta, J is only J_theta
            p%mf_EMF(:,1)=p%mf_EMF(:,1)+delta_effect*delta_tmp*p%jj(:,3)
            p%mf_EMF(:,3)=p%mf_EMF(:,3)-delta_effect*delta_tmp*p%jj(:,1)
          else
            p%mf_EMF(:,1)=p%mf_EMF(:,1)-delta_effect*delta_tmp*p%jj(:,2)
            p%mf_EMF(:,2)=p%mf_EMF(:,2)+delta_effect*delta_tmp*p%jj(:,1)
          endif
        endif
!
!  apply pumping effect: EMF=...-.5*grad(etat) x B
!
        if (lmeanfield_pumping) then
          fact=.5*meanfield_pumping
          call cross_mn(meanfield_getat_tmp, p%bb, getat_cross_B_tmp)
          p%mf_EMF=p%mf_EMF-fact*getat_cross_B_tmp
        endif
!
!  apply pumping effect in spherical coordinates
!
        if (gamma_effect/=0.) then
          fact=gamma_effect
          p%mf_EMF(:,2)=p%mf_EMF(:,2)-fact*p%bb(:,3)
          p%mf_EMF(:,3)=p%mf_EMF(:,3)+fact*p%bb(:,2)
        endif
!
!  Apply diffusion term: simple in Weyl gauge, which is not the default!
!  In diffusive gauge, add (divA) grad(etat) term.
!
        if (lweyl_gauge) then
          call multsv_mn_add(-meanfield_etat_tmp,p%jj,p%mf_EMF)
        else
          call multsv_mn_add(meanfield_etat_tmp,p%del2a,p%mf_EMF)
          call multsv_mn_add(p%diva,meanfield_getat_tmp,p%mf_EMF)
        endif
!
!  Allow for possibility of variable etat from the f-array.
!
        if (ietat/=0) then
          call multsv_mn_add(-f(l1:l2,m,n,ietat),p%jj,p%mf_EMF)
        endif
!
!  Possibility of adding contribution from large-scale velocity.
!
        if (llarge_scale_velocity) p%mf_EMF=p%mf_EMF+p%uxb
!
!  Possibility of turning EMF to zero in a certain region.
!
        if (lEMF_profile) then
          select case (EMF_profile)
          case ('xcutoff');
            EMF_prof= 1+stepdown(x(l1:l2),alpha_rmax,alpha_width)
          case ('surface_z');
            EMF_prof=0.5*(1.-erfunc((z(n)-z_surface)/alpha_width))
          case ('nothing');
          call inevitably_fatal_error('calc_pencils_magnetic', &
            'lEMF_profile=T, but no profile selected !')
          endselect
            p%mf_EMF(:,1)=p%mf_EMF(:,1)*EMF_prof(:)
            p%mf_EMF(:,2)=p%mf_EMF(:,2)*EMF_prof(:)
            p%mf_EMF(:,3)=p%mf_EMF(:,3)*EMF_prof(:)
        endif
      endif
!
!  Evaluate exa, but do this at the end after mf_EMF has been
!  fully assembled.
!
      if (lpencil(i_exa)) then
        if (lmagn_mf) then
          call cross_mn(-p%mf_EMF,p%aa,exa_meanfield)
          p%exa=p%exa+exa_meanfield
        endif
      endif
!
!  EMFdotB
!
      if (lpencil(i_mf_EMFdotB)) call dot_mn(p%mf_EMF,p%bb,p%mf_EMFdotB)
!
!  Compute chiB_term in specific entropy equation.
!  Ds/Dt = ... (1/rho*T)*div[rho*T*chit(B)*grads]
!        = ... chit(B)*{[gradln(rho*T)+gradln(chit(B))].gs+del2s}
!        = ... chit(B)*[gradln(rho*T).grads+del2s] + chit(B)*gradln(chit(B)).gs
!
      if (lpencil(i_chiB_mf)) then
        select case (meanfield_Beq_profile)
        case ('uturbconst');
          Beq21=mu01*p%rho1/(uturb**2)
        case default;
          Beq21=1./meanfield_Beq**2
        endselect
!
!  Choice between 2 versions:
!
        oneQbeta2=1.+chit_quenching*p%b2*Beq21
!
!  Currently no additional profile, so zero gradient.
!
        chit_prof=1.
        glnchit_prof=0.
!
!  Add contribution from gradient of chiB*B^2/Beq^2 term.
!  The B2glnrho term is critical for the magneto-thermodiffusive instability.
!
        call multsv_mn(p%b2,p%glnrho,B2glnrho)
        if (lignore_gradB2_inchiB) then
          call multsv_mn(-chit_quenching*Beq21/oneQbeta2,-B2glnrho,glnchit)
        else
          call multmv_transp(p%bij,p%bb,Bk_Bki) !=1/2 grad B^2
          call multsv_mn(-chit_quenching*Beq21/oneQbeta2,2.*Bk_Bki-B2glnrho,glnchit)
        endif
!
!  In the presence of a uniform imposed field, chit or calKt are no
!  longer constant. To prevent this, we compensate by a 1+Q*B0^2 factor.
!  This leads to an additional contribution in the glnchit expression.
!
        if (lchit_Bext2_equil) then
          oneQbeta02=1.+chit_quenching*B_ext2*Beq21
          call multsv_mn(-chit_quenching*B_ext2*Beq21/oneQbeta02,p%glnrho,glnchit2)
          glnchit=glnchit+glnchit2
        endif
!
!       if (lchit_with_glnTT) then
!         call dot(p%glnrho+p%glnTT+glnchit_prof+glnchit,p%gss,g2)
!         p%chiB_mf=chi_t0*quench_chiB*(g2+p%del2ss)
!       else
!
!  The lrho_chit option corresponds to solving the equation
!  Ds/Dt = (calKt/rho)*(glncalKt.grads + del2s). Otherwise we have
!  Ds/Dt = chit*[(glnrho+glncalKt).grads + del2s].
!  This corresponds to turbulent diffusion with entropy gradient.
!
          if (lrho_chit) then
            call dot(glnchit_prof+glnchit,p%gss,g2)
            if (lchit_Bext2_equil) then
              p%chiB_mf=p%rho1*chi_t0*oneQbeta02/oneQbeta2*(g2+p%del2ss)
            else
              p%chiB_mf=p%rho1*chi_t0/oneQbeta2*(g2+p%del2ss)
            endif
!
!  turbulent diffusion with temperature gradient
!
          elseif (lturb_temp_diff) then
            call dot(p%glnrho+p%glnTT+glnchit_prof+glnchit,p%glnTT,g2)
            if (lchit_Bext2_equil) then
              p%chiB_mf=p%cp*chi_t0*oneQbeta02/oneQbeta2*(g2+p%del2lnTT)
            else
              p%chiB_mf=p%cp*chi_t0/oneQbeta2*(g2+p%del2lnTT)
            endif
!
!  Turbulent diffusion with entropy gradient and coefficient
!  proportional to chi_t*rho, without temperature factor.
!
          else
            call dot(p%glnrho+glnchit_prof+glnchit,p%gss,g2)
            if (lchit_Bext2_equil) then
              p%chiB_mf=chi_t0*oneQbeta02/oneQbeta2*(g2+p%del2ss)
            else
              p%chiB_mf=chi_t0/oneQbeta2*(g2+p%del2ss)
            endif
          endif
!       endif
!?? end of wrong indentation??
!?? end of chit and chiB apparently (!AB)
      endif
!
!  Calculate diagnostics.
!
      if (ldiagnos) then
        if (idiag_qsm/=0) call sum_mn_name(meanfield_qs_func,idiag_qsm)
        if (idiag_qpm/=0) call sum_mn_name(meanfield_qp_func,idiag_qpm)
        if (idiag_qem/=0) call sum_mn_name(meanfield_qe_func,idiag_qem)
      endif
!
!  1d-averages. Happens at every it1d timesteps, NOT at every it1.
!
      if (l1davgfirst .or. (ldiagnos .and. ldiagnos_need_zaverages)) then
        call xysum_mn_name_z(meanfield_qp_func,idiag_qpmz)
      endif
!
!  2-D averages.
!  y-averaged alpha effect
!
      if (l2davgfirst) then
        if (idiag_alpmxz/=0)  then
          call ysum_mn_name_xz(alpha_total,idiag_alpmxz)
        endif
      endif
!
    endsubroutine calc_pencils_magn_mf
!***********************************************************************
    subroutine daa_dt_meanfield(f,df,p)
!
!  add mean-field evolution to magnetic field.
!
!  27-jul-10/axel: coded
!
      use Diagnostics
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
!
      intent(in) :: p
      intent(inout) :: df,f
!
      real, dimension(nx) :: diffus_eta
!
!  Identify module and boundary conditions.
!
      if (headtt.or.ldebug) print*,'daa_dt_meanfield: SOLVE'
!
!  Add jxb/rho to momentum equation.
!
      if (lhydro) then
        if (lmeanfield_jxb) df(l1:l2,m,n,iux:iuz)=df(l1:l2,m,n,iux:iuz)+p%jxbr_mf
      endif
!
!  Add div(chit(B)*rho*T*grads) to entropy equation
!AB: the following is not correct, so disable it for now!
!
!      if (lentropy) then
!        if (lmeanfield_chitB) df(l1:l2,m,n,iss)=df(l1:l2,m,n,iss)+p%chiB_mf
!      endif
!
!  Multiply resistivity by Nyquist scale, for resistive time-step.
!  We include possible contribution from meanfield_etat, which is however
!  only invoked in mean field models.
!  Allow for variable etat (mean field theory).
!
      if (lfirst.and.ldt) then
        diffus_eta=meanfield_etat*dxyz_2
        if (headtt.or.ldebug) &
          print*, 'daa_dt_meanfield: max(diffus_eta)  =', maxval(diffus_eta)
        maxdiffus=max(maxdiffus,diffus_eta)
      endif
!
!  Alpha effect.
!  Additional terms if Mean Field Theory is included.
!
      if (lmagn_mf .and. &
        (meanfield_etat/=0.0 .or. ietat/=0 .or. &
        alpha_effect/=0.0 .or. delta_effect/=0.0) .or. &
        gamma_effect/=0.0 .or. lread_alpha_tensor_z .or. &
        lread_eta_tensor_z .or. lread_alpha_tensor_z_as_y .or. &
        lread_eta_tensor_z_as_y) then
!
!  Apply p%mf_EMF only if .not.lmagn_mf_demfdt; otherwise postpone.
!
        if (.not.lmagn_mf_demfdt) then
          df(l1:l2,m,n,iax:iaz)=df(l1:l2,m,n,iax:iaz)+p%mf_EMF
        endif
!
!  Apply Omega effect. This is currently applied directly to A
!  and is therefore not part of the pencil p%mf_EMF.
!
        if (lOmega_effect) call Omega_effect(f,df,p)
      endif
!
!  Impose a By(x)=B0*sinkx field by adding a dAz/dt= ... + (B/tau*k)*coskx term.
!
      if (lrhs_term) df(l1:l2,m,n,iaz)=df(l1:l2,m,n,iaz)+rhs_term
!
!  2-component external emf
!
      if (lrhs_term2) then
        df(l1:l2,m,n,iay)=df(l1:l2,m,n,iay)+rhs_termy
        df(l1:l2,m,n,iaz)=df(l1:l2,m,n,iaz)+rhs_termz
      endif
!
!  Calculate diagnostic quantities.
!  Diagnostic output for mean field dynamos.
!
      if (ldiagnos) then
        if (idiag_EMFdotBm/=0) call sum_mn_name(p%mf_EMFdotB,idiag_EMFdotBm)
        if (idiag_EMFdotB_int/=0) call integrate_mn_name(p%mf_EMFdotB,idiag_EMFdotB_int)
      endif ! endif (ldiagnos)
!
!  1d-averages. Happens at every it1d timesteps, NOT at every it1.
!
      if (l1davgfirst .or. (ldiagnos .and. ldiagnos_need_zaverages)) then
!
!  Calculate magnetic helicity flux (ExA contribution).
!
        if (idiag_EMFmz1/=0) call xysum_mn_name_z(p%mf_EMF(:,1),idiag_EMFmz1)
        if (idiag_EMFmz2/=0) call xysum_mn_name_z(p%mf_EMF(:,2),idiag_EMFmz2)
        if (idiag_EMFmz3/=0) call xysum_mn_name_z(p%mf_EMF(:,3),idiag_EMFmz3)
      endif
!
!  2-D averages.
!  y-averaged effective magnetic pressure,
!  only makes sense for the 'uturbconst' profile
!
      if (l2davgfirst) then
        if (idiag_peffmxz/=0)  then
          call ysum_mn_name_xz(.5*(1.-p%mf_qp)*p%b2*p%mf_Beq21,idiag_peffmxz)
        endif
      endif
!  Note that this does not necessarily happen with ldiagnos=.true.
!
!     if (l2davgfirst) then
!       if (idiag_Ezmxz/=0) call ysum_mn_name_xz(p%uxb(:,3),idiag_Ezmxz)
!     else
!
!  We may need to calculate bxmxy without calculating bmx. The following
!  if condition was messing up calculation of bmxy_rms
!
!       if (ldiagnos) then
!         if (idiag_bxmxy/=0) call zsum_mn_name_xy(p%bb(:,1),idiag_bxmxy)
!       endif
!     endif
!
!  Time-advance of secondary mean-field modules.
!
      if (lmagn_mf_demfdt) call demf_dt_meanfield(f,df,p)
!
    endsubroutine daa_dt_meanfield
!***********************************************************************
    subroutine meanfield_chitB(rho,b2,quench)
!
!  Calculate magnetic properties needed for z boundary conditions.
!  This routine contails calls to more specialized routines.
!
!   8-jun-13/axel: coded
!
      real, dimension (nx) :: rho,b2,Beq21,quench
!
!  compute Beq21 = 1/Beq^2
!
      select case (meanfield_Beq_profile)
      case ('uturbconst');
        Beq21=mu01/(rho*uturb**2)
      case default;
        Beq21=1./meanfield_Beq**2
      endselect
!
!  compute chit_quenching
!
      quench=1./(1.+chit_quenching*b2*Beq21)
!
    endsubroutine meanfield_chitB
!***********************************************************************
    subroutine Omega_effect(f,df,p)
!
!  Omega effect coded (normally used in context of mean-field theory)
!  Can do uniform shear (0,Sx,0), and the cosx*cosz profile (solar CZ).
!  In most cases the Omega effect can be modeled using hydro_kinematic,
!  but this is not possible when the flow varies in a direction that
!  is not a coordinate direction, e.g. when U=(0,Sx,0) and A=A(z,t).
!  In such cases the Omega effect must be rewritten in terms of
!  velocity gradients operating on A, i.e. (gradU)^T.A, instead of UxB.
!
!  30-apr-05/axel: coded
!
      use Sub, only: stepdown
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
      real kx
!
      intent(in) :: f,p
      intent(inout) :: df
!
!  use gauge transformation, uxB = -Ay*grad(Uy) + gradient-term
!
      select case (Omega_profile)
      case ('nothing'); print*,'Omega_profile=nothing'
      case ('(0,Sx,0)')
        if (headtt) print*,'Omega_effect: uniform shear in x, S=',Omega_ampl
        df(l1:l2,m,n,iax)=df(l1:l2,m,n,iax)-Omega_ampl*f(l1:l2,m,n,iay)
      case ('(0,0,Sx)')
        if (headtt) print*,'Omega_effect: uniform shear in x, S=',Omega_ampl
        df(l1:l2,m,n,iax)=df(l1:l2,m,n,iax)-Omega_ampl*f(l1:l2,m,n,iaz)
      case ('(0,0,pomSx)')
        if (headtt) print*,'Omega_effect: uniform shear in radius, S=',Omega_ampl
        df(l1:l2,m,n,iax)=df(l1:l2,m,n,iax)-Omega_ampl*f(l1:l2,m,n,iaz)*p%x_mn
      case ('(Sz,0,0)')
        if (headtt) print*,'Omega_effect: uniform shear in z, S=',Omega_ampl
        df(l1:l2,m,n,iaz)=df(l1:l2,m,n,iaz)-Omega_ampl*f(l1:l2,m,n,iax)
        if (lhydro) df(l1:l2,m,n,iux)=df(l1:l2,m,n,iux)-Omega_ampl*f(l1:l2,m,n,iuz)
      case ('(0,cosx*cosz,0)')
        if (headtt) print*,'Omega_effect: solar shear, S=',Omega_ampl
        df(l1:l2,m,n,iax)=df(l1:l2,m,n,iax)+Omega_ampl*f(l1:l2,m,n,iay) &
            *sin(x(l1:l2))*cos(z(n))
        df(l1:l2,m,n,iaz)=df(l1:l2,m,n,iaz)+Omega_ampl*f(l1:l2,m,n,iay) &
            *cos(x(l1:l2))*sin(z(n))
      case ('(0,0,cosx)')
        kx=2*pi/Lx
        if (headtt) print*,'Omega_effect: (0,0,cosx), S,kx=',Omega_ampl,kx
        df(l1:l2,m,n,iax)=df(l1:l2,m,n,iax)+Omega_ampl*f(l1:l2,m,n,iaz) &
            *kx*sin(kx*x(l1:l2))
      case ('(0,0,siny)')
        if (headtt) print*,'Omega_effect: (0,0,siny), Omega_ampl=',Omega_ampl
        df(l1:l2,m,n,iax)=df(l1:l2,m,n,iax)+Omega_ampl*f(l1:l2,m,n,iaz) &
            *sin(y(m))
      case('rcutoff_sin_theta')
        if (headtt) print*,'Omega_effect: r cutoff sin(theta), Omega_ampl=',Omega_ampl
!        df(l1:l2,m,n,iax)=Omega_ampl*df(l1:l2,m,n,iax)
        df(l1:l2,m,n,iax)=df(l1:l2,m,n,iax)-Omega_ampl*p%bb(:,2) &
            *sin(y(m))*x(l1:l2)*(1+stepdown(x(l1:l2),Omega_rmax,Omega_rwidth))
        df(l1:l2,m,n,iay)=df(l1:l2,m,n,iay)+Omega_ampl*p%bb(:,1) &
            *sin(y(m))*x(l1:l2)*(1+stepdown(x(l1:l2),Omega_rmax,Omega_rwidth))
      case default; print*,'Omega_profile=unknown'
      endselect
!
    endsubroutine Omega_effect
!***********************************************************************
    subroutine read_magn_mf_init_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=magn_mf_init_pars, IOSTAT=iostat)
!
!  read namelist for secondary modules in mean-field theory (if invoked)
!
      if (lmagn_mf_demfdt) call read_magn_mf_demfdt_init_pars(iostat)
!
    endsubroutine read_magn_mf_init_pars
!***********************************************************************
    subroutine write_magn_mf_init_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=magn_mf_init_pars)
!
!  write namelist for secondary modules in mean-field theory (if invoked)
!
      if (lmagn_mf_demfdt) call write_magn_mf_demfdt_init_pars(unit)
!
    endsubroutine write_magn_mf_init_pars
!***********************************************************************
    subroutine read_magn_mf_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=magn_mf_run_pars, IOSTAT=iostat)
!
!  read namelist for secondary modules in mean-field theory (if invoked)
!
      if (lmagn_mf_demfdt) call read_magn_mf_demfdt_run_pars(iostat)
!
    endsubroutine read_magn_mf_run_pars
!***********************************************************************
    subroutine write_magn_mf_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=magn_mf_run_pars)
!
!  write namelist for secondary modules in mean-field theory (if invoked)
!
      if (lmagn_mf_demfdt) call write_magn_mf_demfdt_run_pars(unit)
!
    endsubroutine write_magn_mf_run_pars
!***********************************************************************
    subroutine rprint_magn_mf(lreset,lwrite)
!
!  Reads and registers print parameters relevant for magnetic fields.
!
!   3-may-02/axel: coded
!  27-may-02/axel: added possibility to reset list
!
      use Diagnostics, only: parse_name
!
      integer :: iname,inamey,inamez,ixy,ixz,irz,inamer
!      integer :: inamex
      logical :: lreset,lwr
      logical, optional :: lwrite
!
      lwr = .false.
      if (present(lwrite)) lwr=lwrite
!
!  Reset everything in case of RELOAD.
!  (this needs to be consistent with what is defined above!)
!
      if (lreset) then
        idiag_qsm=0; idiag_qpm=0; idiag_qem=0;
        idiag_EMFmz1=0; idiag_EMFmz2=0; idiag_EMFmz3=0; idiag_peffmxz=0
        idiag_qpmz=0; idiag_alpmxz=0
      endif
!
!  Check for those quantities that we want to evaluate online.
!
      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),'qsm',idiag_qsm)
        call parse_name(iname,cname(iname),cform(iname),'qpm',idiag_qpm)
        call parse_name(iname,cname(iname),cform(iname),'qem',idiag_qem)
      enddo
!
!  Check for those quantities for which we want xy-averages.
!
!     do inamex=1,nnamex
!    enddo
!
!  Check for those quantities for which we want xz-averages.
!
      do inamey=1,nnamey
      enddo
!
!  Check for those quantities for which we want yz-averages.
!
      do inamez=1,nnamez
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'EMFmz1',idiag_EMFmz1)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'EMFmz2',idiag_EMFmz2)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'EMFmz3',idiag_EMFmz3)
        call parse_name(inamez,cnamez(inamez),cformz(inamez),'qpmz',idiag_qpmz)
      enddo
!
!  Check for those quantities for which we want y-averages.
!
      do ixz=1,nnamexz
        call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'peffmxz',idiag_peffmxz)
        call parse_name(ixz,cnamexz(ixz),cformxz(ixz),'alpmxz',idiag_alpmxz)
      enddo
!
!  Check for those quantities for which we want z-averages.
!
      do ixy=1,nnamexy
      enddo
!
!  Check for those quantities for which we want phi-averages.
!
      do irz=1,nnamerz
      enddo
!
!  Check for those quantities for which we want phiz-averages.
!
      do inamer=1,nnamer
      enddo
!
!  Write column, idiag_XYZ, where our variable XYZ is stored.
!
      if (lwr) then
      endif
!
    endsubroutine rprint_magn_mf
!***********************************************************************
    subroutine pc_aasb_const_alpha(f,topbot,j)
!
!  Perfect-conductor BC for mean field in 1D (only z dependent) with constant alpha and etaT,
!  IF A IS CONSIDERED AS B. 
!  Only implemented at lower z boundary.
!
!  2-jan-17/MR: coded
!
    real, dimension(:,:,:,:) :: f
    integer :: j
    character (LEN=*) :: topbot

    real :: g0, D
    integer :: k

      if (topbot=='bot') then

        if (j/=iax) return

        g0=-147./(60.*dz); D=(alpha_effect**2+g0**2)*60.*dz
        k=n1
        f(:,:,n1,iax) = ( 360.*(-g0*f(:,:,k+1,iax)-alpha_effect*f(:,:,k+1,iay)) &
                         -450.*(-g0*f(:,:,k+2,iax)-alpha_effect*f(:,:,k+2,iay)) &
                         +400.*(-g0*f(:,:,k+3,iax)-alpha_effect*f(:,:,k+3,iay)) &
                         -225.*(-g0*f(:,:,k+4,iax)-alpha_effect*f(:,:,k+4,iay)) &
                         + 72.*(-g0*f(:,:,k+5,iax)-alpha_effect*f(:,:,k+5,iay)) &
                         - 10.*(-g0*f(:,:,k+6,iax)-alpha_effect*f(:,:,k+6,iay))  )/D

        f(:,:,n1,iay) = ( 360.*(-g0*f(:,:,k+1,iay)+alpha_effect*f(:,:,k+1,iax)) &
                         -450.*(-g0*f(:,:,k+2,iay)+alpha_effect*f(:,:,k+2,iax)) &
                         +400.*(-g0*f(:,:,k+3,iay)+alpha_effect*f(:,:,k+3,iax)) &
                         -225.*(-g0*f(:,:,k+4,iay)+alpha_effect*f(:,:,k+4,iax)) &
                         + 72.*(-g0*f(:,:,k+5,iay)+alpha_effect*f(:,:,k+5,iax)) &
                         - 10.*(-g0*f(:,:,k+6,iay)+alpha_effect*f(:,:,k+6,iax))  )/D
      else 
        stop
      endif
      
    endsubroutine pc_aasb_const_alpha
!***********************************************************************
endmodule Magnetic_meanfield