! $Id$
!
! This special module will contribute a mean electromotive
! force to the vector potential differential equation.
!
! Calculation is done term by term, where each contribution
! is calculated using tensors supplied in a HDF5 file emftensors.h5.
!
! Simo Tuomisto, simo.tuomisto@aalto.fi
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of special_dummies.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lspecial = .true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 0
!
! PENCILS PROVIDED alpha_coefs(3,3); beta_coefs(3,3); gamma_coefs(3)
! PENCILS PROVIDED delta_coefs(3); kappa_coefs(3,3,3); umean_coefs(3)
! PENCILS PROVIDED acoef_coefs(3,3); bcoef_coefs(3,3,3)
! PENCILS PROVIDED alpha_emf(3); beta_emf(3); gamma_emf(3)
! PENCILS PROVIDED delta_emf(3); kappa_emf(3); umean_emf(3)
! PENCILS PROVIDED acoef_emf(3); bcoef_emf(3)
! PENCILS PROVIDED bij_symm(3,3)
! PENCILS PROVIDED emf(3)
!
!***************************************************************
!
module Special
!
  use Cparam
  use Cdata
  use Diagnostics
  use General, only: keep_compiler_quiet,numeric_precision,itoa
  use Messages, only: svn_id, fatal_error, fatal_error_local, warning
  use Mpicomm, only: mpibarrier,MPI_COMM_WORLD,MPI_INFO_NULL,mpireduce_min, mpireduce_max,mpibarrier

  use Sub, only: dot_mn, dot_mn_vm, curl_mn, cross_mn, vec_dot_3tensor,dot2_mn
  use HDF5
  use File_io, only: parallel_unit
!
  implicit none
!
  include '../special.h'

  real, dimension(nx) :: diffus_special=0., advec_special=0.

  ! HDF debug parameters:

  integer :: hdferr
  logical :: hdf_grid_exists

  ! Main HDF object ids and variables

  character (len=fnlen) :: hdf_emftensors_filename, emftensors_file='emftensors.h5'
  integer(HID_T) :: hdf_memtype, &
                    hdf_emftensors_file,  &
                    hdf_emftensors_plist, &
                    hdf_emftensors_group, &
                    hdf_grid_group
!
!  dimensions of the stored coefficients, needs in general to be learned from the input HDF5-file.
!  nx_stored, ny_stored at the moment assumed to coincide with nx, ny.
!
  integer :: nx_stored=nx, ny_stored=ny, nz_stored=1
  integer, parameter :: ntensors = 8

  ! Dataset HDF object ids

  integer, parameter :: id_record_SPECIAL_ILOAD=1
  integer(HID_T),   dimension(ntensors) :: tensor_id_G, tensor_id_D, &
                                           tensor_id_S, tensor_id_memS
  integer(HSIZE_T), dimension(ntensors,10) :: tensor_dims, &
                                              tensor_offsets, tensor_counts, &
                                              tensor_memdims, &
                                              tensor_memoffsets, tensor_memcounts
  integer, parameter :: nscalars = 4
  integer(HID_T),   dimension(nscalars) :: scalar_id_D, scalar_id_S
  integer(HSIZE_T), dimension(nscalars) :: scalar_dims, &
                                           scalar_offsets, scalar_counts, &
                                           scalar_memdims, &
                                           scalar_memoffsets, scalar_memcounts

  ! Actual datasets

  real, dimension(:,:,:,:,:,:)  , allocatable :: alpha_data, beta_data, acoef_data
  real, dimension(:,:,:,:,:)    , allocatable :: gamma_data, delta_data, umean_data
  real, dimension(:,:,:,:,:,:,:), allocatable :: kappa_data, bcoef_data

  real, dimension(:), allocatable :: tensor_times

  logical, dimension(3,3)::lalpha_sym=reshape((/.false.,.true. ,.false., &
                                                .true. ,.false.,.true. , &
                                                .false.,.true. ,.false. /), shape(lalpha_sym)), &
                           lbeta_sym =reshape((/.true. ,.false.,.true. , &
                                                .false.,.true. ,.false., &
                                                .true. ,.false.,.true.  /), shape(lbeta_sym))
  logical, dimension(3) :: lgamma_sym  =[.true. ,.false.,.true. ], &
                           ldelta_sym  =[.false.,.true. ,.false.], &
                           lumean_sym  =[.true. ,.false.,.true. ]
  !logical, dimension(3,3,3)::lkappa_sym=[[.false.,.true. ,.false.],[.true. ,.false.,.true. ],[.false.,.true.,.false.], &
  !                                       [.true. ,.false.,.true. ],[.false.,.true. ,.false.],[.true. ,.false.,.true.], &
  !                                       [.false.,.true. ,.false.],[.true. ,.false.,.true. ],[.false.,.true.,.false.]]
  logical, dimension(3,3,3)::lkappa_sym=reshape((/.false.,.true. ,.false.,.true. ,.false.,.true. ,.false.,.true.,.false., &
                                                  .true. ,.false.,.true. ,.false.,.true. ,.false.,.true. ,.false.,.true., &
                                                  .false.,.true. ,.false.,.true. ,.false.,.true. ,.false.,.true.,.false./), &
                                                shape(lkappa_sym))
  ! Dataset mappings

  integer,parameter :: alpha_id=1, beta_id=2,    &
                       gamma_id=3, delta_id=4,   &
                       kappa_id=5, umean_id=6, &
                       acoef_id=7, bcoef_id=8
  integer,parameter :: time_id=1, z_id=2, y_id=3, x_id=4

  integer, dimension(ntensors),parameter :: tensor_ndims = (/ 6, 6, 5, 5, 7, 5, 6, 7 /)
  integer, dimension(nscalars),parameter :: scalar_ndims = (/ 1, 1, 1, 1 /)

  ! Dataset logical variables

  logical, dimension(3,3)  :: lalpha_arr, lbeta_arr, lacoef_arr
  logical, dimension(3)    :: lgamma_arr, ldelta_arr, lumean_arr
  logical, dimension(3,3,3):: lkappa_arr, lbcoef_arr
  logical, dimension(6)    :: lalpha_c, lbeta_c, lacoef_c
  logical, dimension(3)    :: lgamma_c, ldelta_c, lumean_c
  logical, dimension(3,6)  :: lkappa_c, lbcoef_c
  logical :: lalpha=.false., lbeta=.false., lgamma=.false., ldelta=.false., lkappa=.false.
  logical :: lumean=.false., lacoef=.false., lbcoef=.false., lusecoefs=.false.
  logical :: lread_datasets=.true., lread_time_series=.false., lloop=.false., lbblimit=.false.
  real :: alpha_scale, beta_scale, gamma_scale, delta_scale, kappa_scale, utensor_scale, umean_scale, acoef_scale, bcoef_scale

  character (len=fnlen) :: dataset, dataset_type, dataset_name

  character (len=fnlen) :: alpha_name, beta_name,    &
                           gamma_name, delta_name,   &
                           kappa_name, umean_name, &
                           acoef_name, bcoef_name
  character (len=fnlen), dimension(ntensors) :: tensor_names
  character (len=fnlen), dimension(nscalars) :: scalar_names=(/'t','z','y','x'/)
  real, dimension(ntensors) :: tensor_scales, tensor_maxvals, tensor_minvals

  ! Diagnostic variables
  ! Max diagnostics

  integer :: idiag_alphaxmax=0
  integer :: idiag_alphaymax=0
  integer :: idiag_alphazmax=0
  integer :: idiag_betaxmax=0
  integer :: idiag_betaymax=0
  integer :: idiag_betazmax=0
  integer :: idiag_gammaxmax=0
  integer :: idiag_gammaymax=0
  integer :: idiag_gammazmax=0
  integer :: idiag_deltaxmax=0
  integer :: idiag_deltaymax=0
  integer :: idiag_deltazmax=0
  integer :: idiag_kappaxmax=0
  integer :: idiag_kappaymax=0
  integer :: idiag_kappazmax=0
  integer :: idiag_umeanxmax=0
  integer :: idiag_umeanymax=0
  integer :: idiag_umeanzmax=0
  integer :: idiag_acoefxmax=0
  integer :: idiag_acoefymax=0
  integer :: idiag_acoefzmax=0
  integer :: idiag_bcoefxmax=0
  integer :: idiag_bcoefymax=0
  integer :: idiag_bcoefzmax=0
  integer :: idiag_emfxmax=0
  integer :: idiag_emfymax=0
  integer :: idiag_emfzmax=0
  integer :: idiag_emfcoefxmax=0
  integer :: idiag_emfcoefymax=0
  integer :: idiag_emfcoefzmax=0
  integer :: idiag_emfdiffmax=0
  integer :: idiag_emfxdiffmax=0
  integer :: idiag_emfydiffmax=0
  integer :: idiag_emfzdiffmax=0
! RMS diagnostics
  integer :: idiag_alpharms=0
  integer :: idiag_betarms=0
  integer :: idiag_gammarms=0
  integer :: idiag_deltarms=0
  integer :: idiag_kapparms=0
  integer :: idiag_umeanrms=0
  integer :: idiag_acoefrms=0
  integer :: idiag_bcoefrms=0
  integer :: idiag_emfrms=0
  integer :: idiag_emfcoefrms=0
  integer :: idiag_emfdiffrms=0
! timestep diagnostics
  integer :: idiag_dtemf_ave=0
  integer :: idiag_dtemf_dif=0
!
! Regularization diagnostics.
!
  integer :: idiag_nkappareg=0
!
!  2D diagnostics
!
  integer :: idiag_emfxmxy=0, idiag_emfymxy=0, idiag_emfzmxy=0
  integer :: idiag_emfcoefxmxy=0, idiag_emfcoefymxy=0, idiag_emfcoefzmxy=0
  integer :: idiag_alphaxxmxy=0, idiag_alphayymxy=0, idiag_alphazzmxy=0
  integer :: idiag_alphaxymxy=0, idiag_alphaxzmxy=0, idiag_alphayzmxy=0
  integer :: idiag_betaxxmxy=0, idiag_betayymxy=0, idiag_betazzmxy=0
  integer :: idiag_betaxymxy=0, idiag_betaxzmxy=0, idiag_betayzmxy=0
  integer :: idiag_gammaxmxy=0, idiag_gammaymxy=0, idiag_gammazmxy=0
  integer :: idiag_deltaxmxy=0, idiag_deltaymxy=0, idiag_deltazmxy=0
  integer :: idiag_umeanxmxy=0, idiag_umeanymxy=0, idiag_umeanzmxy=0
  integer :: idiag_kappaxxxmxy=0, idiag_kappayxxmxy=0, idiag_kappazxxmxy=0
  integer :: idiag_kappaxxymxy=0, idiag_kappayxymxy=0, idiag_kappazxymxy=0
  integer :: idiag_kappaxxzmxy=0, idiag_kappayxzmxy=0, idiag_kappazxzmxy=0
  integer :: idiag_kappaxyymxy=0, idiag_kappayyymxy=0, idiag_kappazyymxy=0
  integer :: idiag_kappaxyzmxy=0, idiag_kappayyzmxy=0, idiag_kappazyzmxy=0

!
! Interpolation parameters
!
  character (len=fnlen) :: interpname
  integer :: dataload_len=1, tensor_times_len=-1, iload=0
  real, dimension(nx)     :: tmpline
  real, dimension(nx,3)   :: tmppencil,emftmp
  real, dimension(nx,3,3) :: tmptensor
!
! special symmetries
!
  logical :: lsymmetrize=.false.
  integer :: field_symmetry=0
!
! smoothing near boundaries
!
  integer :: nsmooth_rbound=0, nsmooth_thbound=0
!
! Regularization of beta and kappa.
!
  logical :: lregularize_beta=.false., lregularize_kappa=.false., lregularize_kappa_simple=.false., &
             lreconstruct_tensors=.false., &
             lalt_decomp=.false.,&
             lremove_beta_negativ=.false.
  real :: rel_eta=1e-3, jthreshold=0.3, rel_kappa=1e-3    ! must be > 0
  real, pointer :: eta
  real :: kappa_floor=-1e-5
  real, dimension(nx) :: kappa_mask

! Input dataset name

  namelist /special_init_pars/ &
      emftensors_file, &
      lalpha,   lalpha_c,   alpha_name,   alpha_scale, &
      lbeta,    lbeta_c,    beta_name,    beta_scale, &
      lgamma,   lgamma_c,   gamma_name,   gamma_scale, &
      ldelta,   ldelta_c,   delta_name,   delta_scale, &
      lkappa,   lkappa_c,   kappa_name,   kappa_scale, &
      lumean,   lumean_c,   umean_name,   umean_scale, &
      lacoef,   lacoef_c,   acoef_name,   acoef_scale, &
      lbcoef,   lbcoef_c,   bcoef_name,   bcoef_scale, &
      interpname, dataset, dataset_name, dataset_type, &
      lusecoefs, lloop
  namelist /special_run_pars/ &
      emftensors_file, &
      lalpha,   lalpha_c,   alpha_name,   alpha_scale, &
      lbeta,    lbeta_c,    beta_name,    beta_scale, &
      lgamma,   lgamma_c,   gamma_name,   gamma_scale, &
      ldelta,   ldelta_c,   delta_name,   delta_scale, &
      lkappa,   lkappa_c,   kappa_name,   kappa_scale, &
      lumean,   lumean_c,   umean_name,   umean_scale, &
      lacoef,   lacoef_c,   acoef_name,   acoef_scale, &
      lbcoef,   lbcoef_c,   bcoef_name,   bcoef_scale, &
      interpname, dataset, dataset_name, dataset_type, &
      lusecoefs, lloop, lsymmetrize, field_symmetry, &
      nsmooth_rbound, nsmooth_thbound, lregularize_beta, lreconstruct_tensors, &
      lregularize_kappa, lregularize_kappa_simple, lalt_decomp, lremove_beta_negativ, &
      rel_eta, rel_kappa, jthreshold, kappa_floor, lbblimit

  interface loadDataset
    module procedure loadDataset_rank1
    module procedure loadDataset_rank2
    module procedure loadDataset_rank3
  end interface

  interface symmetrize
    module procedure symmetrize_3d
    module procedure symmetrize_4d
  end interface

  interface smooth_thbound
    module procedure smooth_thbound_4d
  end interface

  interface smooth_rbound
    module procedure smooth_rbound_4d
  end interface

  contains
!***********************************************************************
    subroutine register_special
!
!  Set up indices for variables in special modules.
!
!  6-oct-03/tony: coded
!
      integer :: i
      logical :: hdf_exists
!
      if (lroot) call svn_id( &
           "$Id$")
!
      if (lrun) then

        call H5open_F(hdferr)                                              ! Initializes HDF5 library.

        if (numeric_precision() == 'S') then
          if (lroot) write(*,*) 'register_special: loading data as single precision'
          hdf_memtype = H5T_NATIVE_REAL
        else
          if (lroot) write(*,*) 'register_special: loading data as double precision'
          hdf_memtype = H5T_NATIVE_DOUBLE
        end if

        if (lroot) write (*,*) 'register_special: setting parallel HDF5 IO for data file reading'   !MR: Why parallel?
        call H5Pcreate_F(H5P_FILE_ACCESS_F, hdf_emftensors_plist, hdferr)   ! Creates porperty list for HDF5 file.

        if (lmpicomm) &     !MR: doesn't work for nompicomm
          call H5Pset_fapl_mpio_F(hdf_emftensors_plist, MPI_COMM_WORLD, MPI_INFO_NULL, hdferr)

        hdf_emftensors_filename = trim(datadir_snap)//'/'//trim(emftensors_file)

        if (lroot) print *, 'register_special: opening datafile '//trim(hdf_emftensors_filename)// &
                            ' and loading relevant fields into memory...'

        inquire(file=hdf_emftensors_filename, exist=hdf_exists)

        if (.not. hdf_exists) then                                        ! If HDF5 file doesn't exist:
          call H5Pclose_F(hdf_emftensors_plist, hdferr)                   ! Terminates access to property list.
          call H5close_F(hdferr)                                          ! Frees resources used by library.
          call fatal_error('register_special','File '//trim(hdf_emftensors_filename)//' does not exist!')
        end if
!
! Opens HDF5 file for read access only, returns file identifier hdf_emftensors_file.
!
        call H5Fopen_F(hdf_emftensors_filename, H5F_ACC_RDONLY_F, hdf_emftensors_file, hdferr, access_prp = hdf_emftensors_plist)
!
! Checks whether in HDF5 file there is a group /emftensor/.
!
        call H5Lexists_F(hdf_emftensors_file,'/emftensor/', hdf_exists, hdferr)

        if (.not. hdf_exists) then                                         ! If link '/emftensor/' doesn't exist:
          call H5Fclose_F(hdf_emftensors_file, hdferr)                     ! Terminates access to HDF5 file.
          call H5Pclose_F(hdf_emftensors_plist, hdferr)
          call H5close_F(hdferr)
          call fatal_error('register_special','group /emftensor/ does not exist!')
        end if

        call H5Gopen_F(hdf_emftensors_file, 'emftensor', hdf_emftensors_group, hdferr) ! Opens group emftensor in HDF5 file.

      endif
!!      call farray_register_pde('special',ispecial)
!!      call farray_register_auxiliary('specaux',ispecaux)
!!      call farray_register_auxiliary('specaux',ispecaux,communicated=.true.)
!
    endsubroutine register_special
!***********************************************************************
    subroutine initialize_special(f)
!
!  Called after reading parameters, but before the time loop.
!
!  06-oct-03/tony: coded
!
      use Messages, only: information
      use SharedVariables, only: get_shared_variable

      real, dimension (mx,my,mz,mfarray) :: f
      integer :: i,j
      integer(HID_T) :: datatype_id
      logical :: flag
!
      call keep_compiler_quiet(f)

      if (lrun) then

        call get_shared_variable('eta', eta)
        if (trim(dataset_type) == 'time-series' .or. trim(dataset_type) == 'time-crop') then
          lread_time_series=.true.
        elseif (trim(dataset_type) /= 'mean') then
          call fatal_error('initialize_special','Unknown dataset type chosen!')
        endif

        if (lalt_decomp.or.lreconstruct_tensors) then
          lacoef=.true.; lacoef=.true.
        endif

        if (.not.lreloading) then

          call H5Lexists_F(hdf_emftensors_file,'/grid/', hdf_grid_exists, hdferr)
          if (.not. hdf_grid_exists) then
            call warning('initialize_special','error while opening /grid/ - time not available')
            if (lread_time_series) then
              lread_time_series=.false.
              call warning('initialize_special','time-series cannot be read')
            endif
          else

            call H5Gopen_F(hdf_emftensors_file, 'grid', hdf_grid_group, hdferr)
            if (hdferr /= 0) call fatal_error('initialize special','cannot open group "grid"')

            call openDataset_grid(x_id)
            if (scalar_dims(x_id)/=nxgrid) &
              call warning('initialize_special', &
              'x extent '//trim(itoa(int(scalar_dims(x_id))))//' in dataset "grid" different from nxgrid')

            call openDataset_grid(y_id)
            if (scalar_dims(y_id)/=nygrid) &
              call warning('initialize_special', &
              'y extent '//trim(itoa(int(scalar_dims(y_id))))//' in dataset "grid" different from nygrid')

            call openDataset_grid(time_id)

            call H5Dget_type_F(scalar_id_D(time_id), datatype_id, hdferr)
            call H5Tequal_f(datatype_id, hdf_memtype, flag, hdferr)
            if (.not.flag.and.lroot) &
              call information('initialize_special', &
              'Type of stored HDF5 data different from type in memory - converting while reading!')

            if (lread_time_series) then

              allocate(tensor_times(scalar_dims(time_id)))
              call H5Dread_F(scalar_id_D(time_id), hdf_memtype, tensor_times, &
                             [scalar_dims(time_id)], hdferr)
              if (hdferr /= 0) call fatal_error('initialize special','cannot read grid/t')

              t = tensor_times(1)
              if (lroot) then
                print*, 'min time array=',minval(tensor_times)
                print*, 'max time array=',maxval(tensor_times)
              endif

            endif
          endif
        endif

        ! Preset dataset offsets, counts and dimensions
        tensor_offsets = 0
        tensor_counts(:,5:) = 1
        tensor_dims(:,5:) = 3
        tensor_memoffsets = 0
        tensor_memcounts = 1
        tensor_memdims = 3

        if (lalpha) then
          if (.not.allocated(alpha_data)) then
            call openDataset((/'alpha'/),alpha_id)
            allocate(alpha_data(dataload_len,nx_stored,ny_stored,nz_stored,3,3))
          endif
          alpha_data = 0.
        elseif (allocated(alpha_data)) then
          deallocate(alpha_data)
          call closeDataset(alpha_id)
        endif

        if (lbeta) then
          if (.not.allocated(beta_data)) then
            allocate(beta_data(dataload_len,nx_stored, ny_stored,nz_stored,3,3))
            call openDataset((/'beta'/),beta_id)
          endif
          beta_data = 0.
        elseif (allocated(beta_data)) then
          deallocate(beta_data)
          call closeDataset(beta_id)
        endif

        if (lgamma) then
          if (.not.allocated(gamma_data)) then
            allocate(gamma_data(dataload_len,nx_stored, ny_stored,nz_stored,3))
            call openDataset((/'gamma'/),gamma_id)
          endif
          gamma_data = 0.
        elseif (allocated(gamma_data)) then
          deallocate(gamma_data)
          call closeDataset(gamma_id)
        endif

        if (ldelta) then
          if (.not.allocated(delta_data)) then
            allocate(delta_data(dataload_len,nx_stored, ny_stored,nz_stored,3))
            call openDataset((/'delta'/),delta_id)
          endif
          delta_data = 0.
        elseif (allocated(delta_data)) then
          deallocate(delta_data)
          call closeDataset(delta_id)
        endif

        if (lkappa) then
          if (.not.allocated(kappa_data)) then
            allocate(kappa_data(dataload_len,nx_stored, ny_stored,nz_stored,3,3,3))
            call openDataset((/'kappa'/),kappa_id)
          endif
          kappa_data = 0.
        elseif (allocated(kappa_data)) then
          deallocate(kappa_data)
          call closeDataset(kappa_id)
        endif

        if (lumean) then
          if (.not.allocated(umean_data)) then
            allocate(umean_data(dataload_len,nx_stored, ny_stored,nz_stored,3))
            call openDataset((/'umean  ','utensor'/),umean_id)
          endif
          umean_data = 0.
        elseif (allocated(umean_data)) then
          deallocate(umean_data)
          call closeDataset(umean_id)
        endif

        if (lacoef) then
          if (.not.allocated(acoef_data)) then
            allocate(acoef_data(dataload_len,nx_stored, ny_stored,nz_stored,3,3))
            call openDataset((/'acoef'/), acoef_id)
          endif
          acoef_data = 0.
        elseif (allocated(acoef_data)) then
          deallocate(acoef_data)
          call closeDataset(acoef_id)
        endif

        if (lbcoef) then
          if (.not.allocated(bcoef_data)) then
            allocate(bcoef_data(dataload_len,nx_stored,ny_stored,nz_stored,3,3,3))
            call openDataset((/'bcoef'/), bcoef_id)
          endif
          bcoef_data = 0.
        elseif (allocated(bcoef_data)) then
          deallocate(bcoef_data)
          call closeDataset(bcoef_id)
        endif

        if (.not.lusecoefs) then
          idiag_emfcoefxmxy=0; idiag_emfcoefymxy=0; idiag_emfcoefzmxy=0
        endif

        if (lroot) then
          if (lalpha) then
            write (*,*) 'Alpha scale:   ', alpha_scale
            write (*,*) 'Alpha components used: '
            do i=1,3
              write (*,'(A3,3L3,A3)') '|', lalpha_arr(:,i), '|'
            end do
          endif
          if (lbeta) then
            write (*,*) 'Beta scale:   ', beta_scale
            write (*,*) 'Beta components used: '
            do i=1,3
              write (*,'(A3,3L3,A3)') '|', lbeta_arr(:,i), '|'
            end do
          endif
          if (lgamma) then
            write (*,*) 'Gamma scale:   ', gamma_scale
            write (*,*) 'Gamma components used: '
            write (*,'(A3,3L3,A3)') '|', lgamma_arr, '|'
          endif
          if (ldelta) then
            write (*,*) 'Delta scale:   ', delta_scale
            write (*,*) 'Delta components used: '
            write (*,'(A3,3L3,A3)') '|', ldelta_arr, '|'
          endif
          if (lkappa) then
            write (*,*) 'Kappa scale:   ', kappa_scale
            write (*,*) 'Kappa components used: '
            do j=1,3
              write(*,*) '|'
              do i=1,3
                write (*,'(A3,3L3,A3)') '||', lkappa_arr(:,j,i), '||'
              end do
            end do
            write(*,*) '|'
          endif
          if (lumean) then
            write (*,*) 'U-mean scale:  ', umean_scale
            write (*,*) 'U-mean components used: '
            write (*,'(A3,3L3,A3)') '|', lumean_arr, '|'
          endif
          if (lacoef) then
            write (*,*) 'acoef scale:   ', acoef_scale
            write (*,*) 'acoef components used: '
            do i=1,3
              write (*,'(A3,3L3,A3)') '|', lacoef_arr(:,i), '|'
            end do
          endif
          if (lbcoef) then
            write (*,*) 'bcoef scale:   ', bcoef_scale
            write (*,*) 'bcoef components used: '
            do j=1,3
              write(*,*) '|'
              do i=1,3
                write (*,'(A3,3L3,A3)') '||', lbcoef_arr(:,j,i), '||'
              end do
            end do
          endif
        end if

        lread_datasets=.true.

      endif
!
    endsubroutine initialize_special
!***********************************************************************
    subroutine smooth_thbound_4d(arr,nsmooth)
!
!  Simple smothing at theta boundaries: last nsmooth points follow linear trend,
!  defined by last unsmoothed point and zero in boundary point.
!
!  20-jan-2019/MR: coded
!
      real, dimension(:,:,:,:), intent(INOUT) :: arr
      integer, intent(IN) :: nsmooth

      integer :: len_theta,len_r,ir,is
      real, dimension(size(arr,1),size(arr,4)) :: del

      len_theta=size(arr,3); len_r=size(arr,2)

      do ir=1,len_r

        if (lfirst_proc_y) then
          del=arr(:,ir,nsmooth+1,:)/nsmooth
          do is=1,nsmooth
            arr(:,ir,nsmooth+1-is,:)=arr(:,ir,nsmooth+1,:)-is*del
          enddo
          arr(:,ir,1,:)=0.
        endif

        if (llast_proc_y) then
          del=arr(:,ir,len_theta-nsmooth,:)/nsmooth
          do is=1,nsmooth
            arr(:,ir,len_theta-nsmooth+is,:)=arr(:,ir,len_theta-nsmooth,:)-is*del
          enddo
          arr(:,ir,len_theta,:)=0.
        endif

      enddo

    endsubroutine smooth_thbound_4d
!***********************************************************************
    subroutine smooth_rbound_4d(arr,nsmooth)
!
!  Simple smothing at r boundaries: last nsmooth points follow linear trend,
!  defined by last unsmoothed point and zero in boundary point.
!
!  20-jan-2019/MR: coded
!
      real, dimension(:,:,:,:), intent(INOUT) :: arr
      integer, intent(IN) :: nsmooth

      integer :: len_theta,len_r,ith,is
      real, dimension(size(arr,1),size(arr,4)) :: del

      len_theta=size(arr,3); len_r=size(arr,2)

      do ith=1,len_theta

        if (lfirst_proc_x) then
          del=arr(:,nsmooth+1,ith,:)/nsmooth
          do is=1,nsmooth
            arr(:,nsmooth+1-is,ith,:)=arr(:,nsmooth+1,ith,:)-is*del
          enddo
!          arr(:,ir,1,:)=0.
        endif

        if (llast_proc_x) then
          del=arr(:,len_r-nsmooth,ith,:)/nsmooth
          do is=1,nsmooth
            arr(:,len_theta-nsmooth+is,ith,:)=arr(:,len_theta-nsmooth,ith,:)-is*del
          enddo
!          arr(:,ir,len_theta,:)=0.
        endif

      enddo

    endsubroutine smooth_rbound_4d
!***********************************************************************
    subroutine symmetrize_4d(arr,lsym)

      use Mpicomm, only: mpisendrecv_real,mpibarrier,MPI_ANY_TAG
      use General, only: find_proc

      real, dimension(:,:,:,:), intent(INOUT) :: arr
      logical                 , intent(IN)    :: lsym

      integer :: len_theta,len_theta_h,symthproc
      integer, dimension(4) :: sz
      logical :: lmiddle
      real, dimension(:,:,:,:), allocatable :: buffer

      len_theta=size(arr,3); len_theta_h=floor(len_theta/2.)
      lmiddle=mod(nprocy,2)/=0.and.ipy==floor(nprocy/2.)
      sz=(/size(arr,1),size(arr,2),len_theta,size(arr,4)/)

      if (lmiddle) then

        if (lsym) then
          arr(:,:,:len_theta_h,:) = 0.5*(arr(:,:,:len_theta_h,:)+arr(:,:,len_theta:len_theta_h:-1,:))
          arr(:,:,len_theta:len_theta_h:-1,:) = arr(:,:,:len_theta_h,:)
        else
          arr(:,:,:len_theta_h,:) = 0.5*(arr(:,:,:len_theta_h,:)-arr(:,:,len_theta:len_theta_h:-1,:))
          arr(:,:,len_theta:len_theta_h:-1,:) = -arr(:,:,:len_theta_h,:)
        endif

      else

        allocate(buffer(sz(1),sz(2),sz(3),sz(4)))
        symthproc=find_proc(ipx,nprocy-1-ipy,ipz)
        call mpisendrecv_real(arr,sz,symthproc,iproc,buffer,symthproc,symthproc)

        if (lsym) then
          arr = 0.5*(arr+buffer(:,:,len_theta:1:-1,:))
        else
          arr = 0.5*(arr-buffer(:,:,len_theta:1:-1,:))
        endif

      endif
      call mpibarrier

    endsubroutine symmetrize_4d
!***********************************************************************
    subroutine symmetrize_3d(arr,lsym)

      use Mpicomm, only: mpisendrecv_real,mpibarrier,MPI_ANY_TAG
      use General, only: find_proc

      real, dimension(:,:,:), intent(INOUT) :: arr
      logical               , intent(IN)    :: lsym

      integer :: len_theta,len_theta_h,symthproc
      integer, dimension(3) :: sz
      logical :: lmiddle
      real, dimension(:,:,:), allocatable :: buffer

      len_theta=size(arr,2); len_theta_h=floor(len_theta/2.)
      lmiddle=mod(nprocy,2)/=0.and.ipy==floor(nprocy/2.)
      sz=(/size(arr,1),len_theta,size(arr,3)/)

      if (lmiddle) then

        if (lsym) then
          arr(:,:len_theta_h,:) = 0.5*(arr(:,:len_theta_h,:)+arr(:,len_theta:len_theta_h:-1,:))
          arr(:,len_theta:len_theta_h:-1,:) = arr(:,:len_theta_h,:)
        else
          arr(:,:len_theta_h,:) = 0.5*(arr(:,:len_theta_h,:)-arr(:,len_theta:len_theta_h:-1,:))
          arr(:,len_theta:len_theta_h:-1,:) = -arr(:,:len_theta_h,:)
        endif

      else

        allocate(buffer(sz(1),sz(2),sz(3)))
        symthproc=find_proc(ipx,nprocy-1-ipy,ipz)
        call mpisendrecv_real(arr,sz,symthproc,iproc,buffer,symthproc,MPI_ANY_TAG)  ! symthproc

        if (lsym) then
          arr = 0.5*(arr+buffer(:,len_theta:1:-1,:))
        else
          arr = 0.5*(arr-buffer(:,len_theta:1:-1,:))
        endif

      endif
      call mpibarrier

    endsubroutine symmetrize_3d
!***********************************************************************
    subroutine finalize_special(f)
!
!    Called right before exiting.
!
!    14-aug-2011/Bourdin.KIS: coded
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
      call keep_compiler_quiet(f)

      if (lrun) then

        if (lroot.and.(headtt.or.ldebug)) print *,'finalize_special: Deallocating coefficients'

        ! Deallocate data
        if (allocated(alpha_data)) deallocate(alpha_data)
        if (allocated(beta_data)) deallocate(beta_data)
        if (allocated(gamma_data)) deallocate(gamma_data)
        if (allocated(delta_data)) deallocate(delta_data)
        if (allocated(kappa_data)) deallocate(kappa_data)
        if (allocated(umean_data)) deallocate(umean_data)

        if (lroot.and.(headtt.or.ldebug)) print *,'finalize_special: Closing datasets'

        if (lalpha) call closeDataset(alpha_id)
        if (lbeta)  call closeDataset(beta_id)
        if (lgamma) call closeDataset(gamma_id)
        if (ldelta) call closeDataset(delta_id)
        if (lkappa) call closeDataset(kappa_id)
        if (lumean) call closeDataset(umean_id)
        if (lacoef) call closeDataset(acoef_id)
        if (lbcoef) call closeDataset(bcoef_id)


        if (hdf_grid_exists) then
          if (lroot.and.(headtt.or.ldebug)) print *,'finalize_special: Closing grid'
          call closeDataset_grid(time_id)
          call closeDataset_grid(x_id)
          call closeDataset_grid(y_id)
          call H5Gclose_F(hdf_grid_group, hdferr)
          if (lread_time_series) then
            if (allocated(tensor_times)) deallocate(tensor_times)
          endif
        endif

        if (lroot.and.(headtt.or.ldebug)) print *,'finalize_special: Closing emftensors group'
        call H5Gclose_F(hdf_emftensors_group, hdferr)
        if (lroot.and.(headtt.or.ldebug)) print *,'finalize_special: Closing emftensors file'
        call H5Fclose_F(hdf_emftensors_file, hdferr)
        if (lroot.and.(headtt.or.ldebug)) print *,'finalize_special: Closing emftensors parameters'
        call H5Pclose_F(hdf_emftensors_plist, hdferr)
        if (lroot.and.(headtt.or.ldebug)) print *,'finalize_special: Closing HDF5'
        call H5close_F(hdferr)
        call mpibarrier

      end if
  !
    endsubroutine finalize_special
!***********************************************************************
    subroutine special_before_boundary(f)
!
!  Possibility to modify the f array before the boundaries are
!  communicated.
!
!  Some precalculated pencils of data are passed in for efficiency
!  others may be calculated directly from the f array
!
!  06-jul-06/tony: coded
!  20-may-19/MR: added beta regularization, reconstruction of alpha and gamma from acoef and bcoef,
!                alternative decomposition of tensors
!
      use Mpicomm, only: mpireduce_sum_int, mpiallreduce_sum_int, mpireduce_min
      use PolynomialRoots, only: cubicroots
      use Sub, only: dyadic2_other

      real, dimension (mx,my,mz,mfarray), intent(INOUT) :: f
!
      real :: delt,minbeta,minbeta_,thmin,thmax,rmin,rmax,trmin,det,oldtrace,newtrace
      integer :: i,j,k,numzeros,numcomplex,numzeros_,mm,ll,iv,iv0,ik,icheck
      integer, dimension(:,:,:,:,:), allocatable :: beta_mask
      real, dimension(4) :: polcoeffs
      complex, dimension(3) :: eigenvals
      real, dimension(3,3) :: ev, beta_sav

      call keep_compiler_quiet(f)

      if (lfirst) then
        if (lread_time_series) then
          if (t >= tensor_times(iload+1)) then
            lread_datasets=.true.
            iload = iload + 1
          endif
        else
          iload=1
        endif

        if (lread_datasets) then
          if (iload > tensor_times_len) then
            if (lloop) then
              iload=1
              delt=2*tensor_times(tensor_times_len)-tensor_times(1)-tensor_times(tensor_times_len-1)
              tensor_times = tensor_times+delt
            else
              call fatal_error('special_before_boundary', 'no more data to load')
            endif
          endif
!if (lread_time_series) &
!print *, 'loading: t,iload=',tensor_times(iload),iload
          ! Load datasets
          if (lalpha) call loadDataset(alpha_data, lalpha_arr, alpha_id, iload-1,'Alpha')
          if (lbeta)  call loadDataset(beta_data,  lbeta_arr,  beta_id,  iload-1,'Beta')
          if (lgamma) call loadDataset(gamma_data, lgamma_arr, gamma_id, iload-1,'Gamma')
          if (ldelta) call loadDataset(delta_data, ldelta_arr, delta_id, iload-1,'Delta')
          if (lkappa) call loadDataset(kappa_data, lkappa_arr, kappa_id, iload-1,'Kappa')
          if (lumean) call loadDataset(umean_data, lumean_arr, umean_id, iload-1,'Umean')
          if (lacoef) call loadDataset(acoef_data, lacoef_arr, acoef_id, iload-1,'Acoef')
          if (lbcoef) call loadDataset(bcoef_data, lbcoef_arr, bcoef_id, iload-1,'Bcoef')
          lread_datasets=.false.

          if (lreconstruct_tensors.and.lacoef.and.lbcoef) then                 ! beta, delta, kappa supposed to be correct
            if (lroot) print*, 'Reconstruct alpha and gamma from raw tensors.'
            do mm=1,ny
              alpha_data(1,:,mm,1,1,1)=     acoef_data(1,:,mm,1,1,1)-bcoef_data(1,:,mm,1,1,2,2)/x(l1:l2)
              alpha_data(1,:,mm,1,1,2)=0.5*(acoef_data(1,:,mm,1,1,2)+bcoef_data(1,:,mm,1,1,1,2)/x(l1:l2)  &
                                           +acoef_data(1,:,mm,1,2,1)-bcoef_data(1,:,mm,1,2,2,2)/x(l1:l2))
              alpha_data(1,:,mm,1,2,2)=     acoef_data(1,:,mm,1,2,2)+bcoef_data(1,:,mm,1,2,1,2)/x(l1:l2)
              alpha_data(1,:,mm,1,1,3)=0.5*(acoef_data(1,:,mm,1,1,3)+acoef_data(1,:,mm,1,3,1) &
                                           -bcoef_data(1,:,mm,1,3,2,2)/x(l1:l2))
              alpha_data(1,:,mm,1,2,3)=0.5*(acoef_data(1,:,mm,1,2,3)+acoef_data(1,:,mm,1,3,2) &
                                           +bcoef_data(1,:,mm,1,3,1,2)/x(l1:l2))

              gamma_data(1,:,mm,1,1)=0.5*(acoef_data(1,:,mm,1,3,2)-acoef_data(1,:,mm,1,2,3) &
                                         +bcoef_data(1,:,mm,1,3,1,2)/x(l1:l2))
              gamma_data(1,:,mm,1,2)=0.5*(acoef_data(1,:,mm,1,1,3)-acoef_data(1,:,mm,1,3,1) &
                                         +bcoef_data(1,:,mm,1,3,2,2)/x(l1:l2))
              gamma_data(1,:,mm,1,3)=0.5*(acoef_data(1,:,mm,1,2,1)-acoef_data(1,:,mm,1,1,2) &
                                         -bcoef_data(1,:,mm,1,1,1,2)/x(l1:l2)-bcoef_data(1,:,mm,1,2,2,2)/x(l1:l2))
            enddo

            alpha_data(:,:,:,:,2,1)=alpha_data(:,:,:,:,1,2)
            alpha_data(:,:,:,:,3,2)=alpha_data(:,:,:,:,2,3)
            alpha_data(:,:,:,:,3,1)=alpha_data(:,:,:,:,1,3)
            alpha_data=alpha_data*alpha_scale
            gamma_data=gamma_data*gamma_scale

          endif

          if (lalt_decomp.and.lacoef.and.lbcoef) then                    ! kappa supposed to be correct.
            if (lroot) print*, 'Construct alternative decomposition from raw tensors.'
            do mm=1,ny

              alpha_data(1,:,mm,1,1,1)=      acoef_data(1,:,mm,1,1,1)-bcoef_data(1,:,mm,1,1,2,2)/x(l1:l2)
              alpha_data(1,:,mm,1,1,2)=0.5*( acoef_data(1,:,mm,1,1,2)+bcoef_data(1,:,mm,1,1,1,2)/x(l1:l2)  &
                                            +acoef_data(1,:,mm,1,2,1)-bcoef_data(1,:,mm,1,2,2,2)/x(l1:l2))
              alpha_data(1,:,mm,1,2,2)=      acoef_data(1,:,mm,1,2,2)+bcoef_data(1,:,mm,1,2,1,2)/x(l1:l2)

              alpha_data(1,:,mm,1,1,3)=0.5*( acoef_data(1,:,mm,1,1,3)+acoef_data(1,:,mm,1,3,1) &
                                            -(                  bcoef_data(1,:,mm,1,3,2,2) &
                                              +                 bcoef_data(1,:,mm,1,1,3,1) &
                                              +cotth(mm+nghost)*bcoef_data(1,:,mm,1,1,3,2))/x(l1:l2))

              alpha_data(1,:,mm,1,2,3)=0.5*( acoef_data(1,:,mm,1,2,3)+acoef_data(1,:,mm,1,3,2) &
                                            -(                  bcoef_data(1,:,mm,1,2,3,1) &
                                              -                 bcoef_data(1,:,mm,1,3,1,2) &
                                              +cotth(mm+nghost)*bcoef_data(1,:,mm,1,2,3,2))/x(l1:l2))

              alpha_data(1,:,mm,1,3,3)= acoef_data(1,:,mm,1,3,3) - (                  bcoef_data(1,:,mm,1,3,3,1) &
                                                                    +cotth(mm+nghost)*bcoef_data(1,:,mm,1,3,3,2))/x(l1:l2)

              gamma_data(1,:,mm,1,1)=0.5*( acoef_data(1,:,mm,1,3,2)-acoef_data(1,:,mm,1,2,3) &
                                          +(                  bcoef_data(1,:,mm,1,2,3,1) &
                                            +                 bcoef_data(1,:,mm,1,3,1,2) &
                                            +cotth(mm+nghost)*bcoef_data(1,:,mm,1,2,3,2))/x(l1:l2))

              gamma_data(1,:,mm,1,2)=0.5*( acoef_data(1,:,mm,1,1,3)-acoef_data(1,:,mm,1,3,1) &
                                          -(                  bcoef_data(1,:,mm,1,1,3,1) &
                                            -                 bcoef_data(1,:,mm,1,3,2,2) &
                                            +cotth(mm+nghost)*bcoef_data(1,:,mm,1,1,3,2))/x(l1:l2))

              gamma_data(1,:,mm,1,3)=0.5*( acoef_data(1,:,mm,1,2,1)-acoef_data(1,:,mm,1,1,2) &
                                         -(bcoef_data(1,:,mm,1,1,1,2)+bcoef_data(1,:,mm,1,2,2,2))/x(l1:l2))

              delta_data(1,:,mm,1,1)=0.25*(bcoef_data(1,:,mm,1,2,2,1)-bcoef_data(1,:,mm,1,2,1,2)+2.*bcoef_data(1,:,mm,1,3,3,1))

              delta_data(1,:,mm,1,2)=0.25*(bcoef_data(1,:,mm,1,1,1,2)-bcoef_data(1,:,mm,1,1,2,1)+2.*bcoef_data(1,:,mm,1,3,3,2))

              delta_data(1,:,mm,1,3)=-0.5*(bcoef_data(1,:,mm,1,1,3,1)+bcoef_data(1,:,mm,1,2,3,2))

              beta_data(1,:,mm,1,1,1)=-bcoef_data(1,:,mm,1,1,3,2)
              beta_data(1,:,mm,1,2,2)= bcoef_data(1,:,mm,1,2,3,1)
              beta_data(1,:,mm,1,3,3)=0.5*(-bcoef_data(1,:,mm,1,3,2,1)+bcoef_data(1,:,mm,1,3,1,2))

              beta_data(1,:,mm,1,1,2)=0.5*(-bcoef_data(1,:,mm,1,2,3,2)+bcoef_data(1,:,mm,1,1,3,1))
              beta_data(1,:,mm,1,1,3)=0.25*( -2.*bcoef_data(1,:,mm,1,3,3,2)+bcoef_data(1,:,mm,1,1,1,2)-bcoef_data(1,:,mm,1,1,2,1))
              beta_data(1,:,mm,1,2,3)=0.25*(2.*bcoef_data(1,:,mm,1,3,3,1)+bcoef_data(1,:,mm,1,2,1,2)-bcoef_data(1,:,mm,1,2,2,1))

              do ik=1,3
                kappa_data(1,:,mm,1,ik,:,3)=0.; kappa_data(1,:,mm,1,ik,3,:)=0.
              enddo
            enddo

            alpha_data(:,:,:,:,2,1)=alpha_data(:,:,:,:,1,2)
            alpha_data(:,:,:,:,3,1)=alpha_data(:,:,:,:,1,3)
            alpha_data(:,:,:,:,3,2)=alpha_data(:,:,:,:,2,3)

            beta_data(:,:,:,:,2,1)=beta_data(:,:,:,:,1,2)
            beta_data(:,:,:,:,3,1)=beta_data(:,:,:,:,1,3)
            beta_data(:,:,:,:,3,2)=beta_data(:,:,:,:,2,3)

            alpha_data=alpha_data*alpha_scale
            gamma_data=gamma_data*gamma_scale
            delta_data=delta_data*delta_scale
            beta_data =beta_data*beta_scale

          endif

          if (nsmooth_rbound/=0) then
            do i=1,3
              if (lgamma) call smooth_rbound(gamma_data(:,:,:,:,i),nsmooth_rbound)
              if (ldelta) call smooth_rbound(delta_data(:,:,:,:,i),nsmooth_rbound)
              if (lumean) call smooth_rbound(umean_data(:,:,:,:,i),nsmooth_rbound)
              do j=1,3
                if (lacoef) call smooth_rbound(acoef_data(:,:,:,:,i,j),nsmooth_rbound)
                if (lalpha) call smooth_rbound(alpha_data(:,:,:,:,i,j),nsmooth_rbound)
                if (lbeta) call smooth_rbound(beta_data(:,:,:,:,i,j),nsmooth_rbound)
                do k=1,3
                  if (lbcoef) call smooth_rbound(bcoef_data(:,:,:,:,i,j,k),nsmooth_rbound)
                  if (lkappa) call smooth_rbound(kappa_data(:,:,:,:,i,j,k),nsmooth_rbound)
                enddo
              enddo
            enddo
          endif

          if (nsmooth_thbound/=0) then
            do i=1,3
              if (lgamma) call smooth_thbound(gamma_data(:,:,:,:,i),nsmooth_thbound)
              if (ldelta) call smooth_thbound(delta_data(:,:,:,:,i),nsmooth_thbound)
              if (lumean) call smooth_thbound(umean_data(:,:,:,:,i),nsmooth_thbound)
              do j=1,3
                if (lacoef) call smooth_thbound(acoef_data(:,:,:,:,i,j),nsmooth_thbound)
                if (lalpha) call smooth_thbound(alpha_data(:,:,:,:,i,j),nsmooth_thbound)
                if (lbeta) call smooth_thbound(beta_data(:,:,:,:,i,j),nsmooth_thbound)
                do k=1,3
                  if (lbcoef) call smooth_thbound(bcoef_data(:,:,:,:,i,j,k),nsmooth_thbound)
                  if (lkappa) call smooth_thbound(kappa_data(:,:,:,:,i,j,k),nsmooth_thbound)
                enddo
              enddo
            enddo
          endif

!if (lroot.and.lbeta) write(200,*) beta_data(1,:,:,1,:,:)

          if (lbeta.and.lremove_beta_negativ) then
            do i=1,3
              where(beta_data(:,:,:,:,i,i)<eta*rel_eta) beta_data(:,:,:,:,i,i)=eta*rel_eta
            enddo
          endif

!if (lroot.and.lbeta) write(300,*) beta_data(1,:,:,1,:,:)

          if (lbeta.and.lregularize_beta) then

            allocate(beta_mask(dataload_len,nx_stored,ny_stored,nz_stored,3))
            beta_mask=0

            do i=1,3
!
!  Look for vanishing diagonal elements of beta.
!
              where(beta_data(:,:,:,:,i,i)+eta<=0.) beta_mask(:,:,:,:,i)=1
              minbeta=minval(beta_data(:,:,:,:,i,i),beta_mask(:,:,:,:,i)==1)
              !where(beta_mask(:,:,:,:,i)==1) beta_data(:,:,:,:,i,i)=(-1.+rel_eta)*eta

              numzeros=sum(beta_mask(:,:,:,:,i))
              call mpiallreduce_sum_int(numzeros,numzeros_)
              if (numzeros_>0) then
                call mpireduce_min(minbeta,minbeta_)
                if (lroot) then
                  print'(a,i1,a,i1,a,i15,a$)', 'beta(',i,',',i,')+eta<=0 at ', numzeros_, ' positions'
                  print'(a,i1,a,i1,a,e12.5)', ', min(beta(',i,',',i,')) = ', minbeta_
                endif
              endif
            enddo
!
!  Look for locations of negative definite beta.
!
            icheck=0
            do

              numzeros=0; numcomplex=0; rmin=x(l2); rmax=x(l1); thmin=y(m2); thmax=y(m1); trmin=impossible
              do mm=1,ny_stored; do ll=1,nx_stored

                oldtrace=beta_data(1,ll,mm,1,1,1)+beta_data(1,ll,mm,1,2,2)+beta_data(1,ll,mm,1,3,3)
                beta_sav=beta_data(1,ll,mm,1,:,:)

                polcoeffs=(/2.*beta_data(1,ll,mm,1,1,2)*beta_data(1,ll,mm,1,1,3)*beta_data(1,ll,mm,1,2,3) &
                             - beta_data(1,ll,mm,1,1,2)**2*beta_data(1,ll,mm,1,3,3) &
                             - beta_data(1,ll,mm,1,1,3)**2*beta_data(1,ll,mm,1,2,2) &
                             - beta_data(1,ll,mm,1,2,3)**2*beta_data(1,ll,mm,1,1,1) &
                             + beta_data(1,ll,mm,1,1,1)*beta_data(1,ll,mm,1,2,2)*beta_data(1,ll,mm,1,3,3), &
!
                               beta_data(1,ll,mm,1,1,2)**2+beta_data(1,ll,mm,1,1,3)**2+beta_data(1,ll,mm,1,2,3)**2 &
                             - beta_data(1,ll,mm,1,1,1)*beta_data(1,ll,mm,1,2,2) &
                             - beta_data(1,ll,mm,1,1,1)*beta_data(1,ll,mm,1,3,3) &
                             - beta_data(1,ll,mm,1,2,2)*beta_data(1,ll,mm,1,3,3), &
!
                               beta_data(1,ll,mm,1,1,1)+beta_data(1,ll,mm,1,2,2)+beta_data(1,ll,mm,1,3,3), &
                             -1. /)
                call cubicroots(polcoeffs, eigenvals)
                if (any(abs(imag(eigenvals))>0.e-9*abs(real(eigenvals)))) numcomplex=numcomplex+1

                newtrace=sum(real(eigenvals))
                if (abs(oldtrace-newtrace)>1.e-9) print*, 'll,mm,oldtrace-newtrace (1)=', &
                  ll,mm,abs(oldtrace-newtrace),oldtrace,newtrace

                if (any(real(eigenvals)<-eta)) then
!
!  If there are negative eigenvalues of beta,
!
!write(100,*) iproc, ll,mm,real(eigenvals)  !, &
!sum(polcoeffs*(/1.d0,real(eigenvals(1)),real(eigenvals(1))**2,real(eigenvals(1))**3/)), &
!sum(polcoeffs*(/1.d0,real(eigenvals(2)),real(eigenvals(2))**2,real(eigenvals(2))**3/)), &
!sum(polcoeffs*(/1.d0,real(eigenvals(3)),real(eigenvals(3))**2,real(eigenvals(3))**3/))
                  numzeros=numzeros+1
                  trmin=min(trmin,sum(real(eigenvals)))
                  rmin=min(rmin,x(ll+nghost)); rmax=max(rmax,x(ll+nghost))
                  thmin=min(thmin,y(mm+nghost)); thmax=max(thmax,y(mm+nghost))
!
!  calculate eigenvectors (v1,v2,-1) for all three eigenvalues.
!
                  ev(:,3)=-1.
                  do iv=1,3
!
! Eigenvalues < -eta are set to (-1.+rel_eta)*eta with rel_eta>0.
!
! output here: eigenvalues JOERN
!
                    if (real(eigenvals(iv))<-eta) eigenvals(iv)=cmplx((-1.+rel_eta)*eta,0.)

                    det = (beta_data(1,ll,mm,1,1,1)-real(eigenvals(iv)))*(beta_data(1,ll,mm,1,2,2)-real(eigenvals(iv))) &
                          -beta_data(1,ll,mm,1,1,2)**2
                    ev(iv,1)= (beta_data(1,ll,mm,1,1,3)*(beta_data(1,ll,mm,1,2,2)-real(eigenvals(iv))) &
                              -beta_data(1,ll,mm,1,2,3)*beta_data(1,ll,mm,1,1,2))/det
                    ev(iv,2)=((beta_data(1,ll,mm,1,1,1)-real(eigenvals(iv)))*beta_data(1,ll,mm,1,2,3) &
                              -beta_data(1,ll,mm,1,1,2)*beta_data(1,ll,mm,1,1,3))/det
                    ev(iv,:)=ev(iv,:)/sqrt(sum(ev(iv,:)**2))    ! normalization
                  enddo

                  beta_data(1,ll,mm,1,:,:)=0.
!
!  Transform back with non-negative eigenvalues.
!
                  oldtrace=0.
                  do iv=1,3
                    beta_data(1,ll,mm,1,:,:)=beta_data(1,ll,mm,1,:,:)+real(eigenvals(iv))*dyadic2_other(ev(iv,:))
                    oldtrace=oldtrace+real(eigenvals(iv))
                  enddo
                  newtrace=beta_data(1,ll,mm,1,1,1)+beta_data(1,ll,mm,1,2,2)+beta_data(1,ll,mm,1,3,3)
                  if (abs(oldtrace-newtrace)>1.e-9) print*, 'll,mm,oldtrace-newtrace (2) =', ll,mm,abs(oldtrace-newtrace)

                else
                  if (icheck==0.and.any(beta_mask(1,ll,mm,1,:)==1)) &
                    call warning('meanfield_e_tensor','beta regularization discrepancy at iproc,l,m='//trim(itoa(iproc))//','// &
                                                      trim(itoa(ll))//','//trim(itoa(mm)),iproc_world)
                endif
111 do i=1,3; do j=1,3
    if (beta_sav(i,j)/=0.) then
!     if (abs(beta_sav(i,j)-beta_data(1,ll,mm,1,i,j))/beta_sav(i,j) > 1e-1)  &
!       print*, 'JOERN', ll, mm, abs((beta_sav(i,j)-beta_data(1,ll,mm,1,i,j))/beta_sav(i,j)),beta_sav(i,j)
    endif
enddo; enddo
              enddo; enddo

              call mpireduce_sum_int(numcomplex,numzeros_)
              if (lroot.and.numzeros_>0) print'(a,i15,a)', 'beta has complex eigenvalues at', numzeros_,' positions.'

              call mpiallreduce_sum_int(numzeros,numzeros_)
              if (numzeros_>0) then
                if (lroot) print'(a,i15,a)', 'beta+eta is negative definite at', numzeros_,' positions.'
                if (numzeros>0) print'(4(a,f6.3),a,i4)', &
                                'in r-theta region (',rmin,',',rmax,')x(',thmin,',',thmax,') of proc ',iproc,'.'
                call mpireduce_min(trmin,minbeta_)
                if (lroot) print'(a,e12.5)', 'minimal trace = ', minbeta_
              else
                exit
              endif

              icheck=icheck+1
            enddo

            deallocate(beta_mask)

          endif

          if (lkappa) then
            if (lregularize_kappa_simple) then
!
              lregularize_kappa=.false.
!
! Setting kappa_floor for kappa_phi,r,theta and kappa_phi,theta,r by hand.
!
              where(kappa_data(1,:,:,1,3,2,1) < kappa_floor) &
                kappa_data(1,:,:,1,3,2,1)= kappa_floor
              where(kappa_data(1,:,:,1,3,1,2) < kappa_floor) &
                kappa_data(1,:,:,1,3,1,2)= kappa_floor
            endif
          endif

!if (lroot.and.lbeta) write(100,*) beta_data(1,:,:,1,:,:)

          if (lsymmetrize) then
            do i=1,3
              if (lgamma) call symmetrize(gamma_data(:,:,:,:,i),lgamma_sym(i))
              if (ldelta) call symmetrize(delta_data(:,:,:,:,i),ldelta_sym(i))
              if (lumean) call symmetrize(umean_data(:,:,:,:,i),lumean_sym(i))
              do j=1,3
                if (lalpha) call symmetrize(alpha_data(:,:,:,:,i,j),lalpha_sym(i,j))
                if (lbeta) call symmetrize(beta_data(:,:,:,:,i,j),lbeta_sym(i,j))
                do k=1,3
                  if (lkappa) call symmetrize(kappa_data(:,:,:,:,i,j,k),lkappa_sym(i,j,k))
                enddo
              enddo
            enddo
          endif

        end if
      end if

      if (field_symmetry==1) then
        call symmetrize_3d(f(:,:,:,iax),.false.)
        call symmetrize_3d(f(:,:,:,iay),.true.)
        call symmetrize_3d(f(:,:,:,iaz),.false.)
      elseif (field_symmetry==-1) then
        call symmetrize_3d(f(:,:,:,iax),.true.)
        call symmetrize_3d(f(:,:,:,iay),.false.)
        call symmetrize_3d(f(:,:,:,iaz),.true.)
      endif

    endsubroutine special_before_boundary
!***********************************************************************
    subroutine pencil_criteria_special
!
!  All pencils that this special module depends on are specified here.
!
      lpenc_requested(i_bb)=.true.
      if (lactive_dimension(3)) then    ! 3D-case
        lpenc_requested(i_bij)=.true.
      else                              ! 2D-case
        lpenc_requested(i_bijtilde)=.true.
      endif
      lpenc_requested(i_jj)=.true.
      !if (lbcoef.and.lusecoefs) lpenc_requested(i_bijtilde)=.true.
!
    endsubroutine pencil_criteria_special
!***********************************************************************
    subroutine calc_pencils_special(f,p)
!
!  Calculate Special pencils.
!  Most basic pencils should come first, as others may depend on them.
!
!  24-nov-04/tony: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
!
      intent(in) :: f
      intent(inout) :: p
!
      integer :: i,j,k,iii,nind
      real, dimension(nx) :: jrt,jtr
      character(LEN=80) :: mess
      character(LEN=20) :: cbuf
      real :: limit=1e30
!
      call keep_compiler_quiet(f)
!
! Check if p%bb is not already to large
!
      if (lbblimit) then
        if (maxval(abs(p%bb)) > limit) call fatal_error_local('calc_pencils_special', &
          'magnetic field too big')
      endif

      nind=1     !n-nghost if coefficients would be z-dependent
!
! Calculate emf pencil
!
      p%emf = 0
!
      if (lacoef) then
        ! Calculate acoef B
        do j=1,3; do i=1,3
          if (lacoef_arr(i,j)) then
            p%acoef_coefs(:,i,j)=emf_interpolate(acoef_data(1:dataload_len,:,m-nghost,nind,i,j))
          else
            p%acoef_coefs(:,i,j)=0
          end if
        end do; end do
        call dot_mn_vm(p%bb,p%acoef_coefs,p%acoef_emf)
      end if
!
      if (lbcoef) then
        ! Calculate bcoef (grad B)
        do k=1,3; do j=1,3; do i=1,3
          if (lbcoef_arr(i,j,k)) then
            p%bcoef_coefs(:,i,j,k)=emf_interpolate(bcoef_data(1:dataload_len,:,m-nghost,nind,i,j,k))
          else
            p%bcoef_coefs(:,i,j,k)=0
          end if
        end do; end do; end do
        p%bcoef_emf = 0
!
!  Use partial (non-covariant) derivatives of B in the form \partial B_{r,theta,phi}/\partial r,
!  \partial B_{r,theta,phi}/(r \partial theta).
!
        do k=1,2; do j=1,3; do i=1,3
          if (lbcoef_arr(i,j,k)) &
            p%bcoef_emf(:,i)=p%bcoef_emf(:,i)+p%bcoef_coefs(:,i,j,k)*p%bijtilde(:,j,k)    !!! tb done: 3D case
        end do; end do; end do
      end if
!
      if (lalpha) then
        ! Calculate alpha B
        do j=1,3; do i=1,3
          if (lalpha_arr(i,j)) then
            p%alpha_coefs(:,i,j)=emf_interpolate(alpha_data(1:dataload_len,:,m-nghost,nind,i,j))
          else
            p%alpha_coefs(:,i,j)=0
          end if
        end do; end do

        call dot_mn_vm(p%bb,p%alpha_coefs,p%alpha_emf)
        p%emf = p%emf + p%alpha_emf
      end if
!
      if (lbeta) then
        ! Calculate beta (curl B)
        do j=1,3; do i=1,3
          if (lbeta_arr(i,j)) then
            p%beta_coefs(:,i,j)=emf_interpolate(beta_data(1:dataload_len,:,m-nghost,nind,i,j))
          else
            p%beta_coefs(:,i,j)=0
          end if
        end do; end do
        call dot_mn_vm(p%jj,p%beta_coefs,p%beta_emf)
        p%emf = p%emf - p%beta_emf
      end if
!
      if (lgamma) then
        ! Calculate gamma x B
        do i=1,3
          if (lgamma_arr(i)) then
            p%gamma_coefs(:,i)=emf_interpolate(gamma_data(1:dataload_len,:,m-nghost,nind,i))
          else
            p%gamma_coefs(:,i)=0
          end if
        end do

        call cross_mn(p%gamma_coefs,p%bb,p%gamma_emf)
        p%emf = p%emf + p%gamma_emf
      end if
!
      if (ldelta) then
        ! Calculate delta x (curl B)
        do i=1,3
          if (ldelta_arr(i)) then
            p%delta_coefs(:,i)=emf_interpolate(delta_data(1:dataload_len,:,m-nghost,nind,i))
          else
            p%delta_coefs(:,i)=0
          end if
        end do
        call cross_mn(p%delta_coefs,p%jj,p%delta_emf)
        p%emf = p%emf - p%delta_emf
      end if

      if (lkappa) then
        ! Calculate (grad B)_symm
        if (lactive_dimension(3)) then
          do j=1,3; do i=1,3
            p%bij_symm(:,i,j)=0.5*(p%bij(:,i,j) + p%bij(:,j,i))
          end do; end do
        else
          do j=1,3; do i=1,3
            p%bij_symm(:,i,j)=0.5*(p%bijtilde(:,i,j)+p%bij_cov_corr(:,i,j) + &
                                   p%bijtilde(:,j,i)+p%bij_cov_corr(:,j,i))
          end do; end do
        endif

        do k=1,3; do j=1,3; do i=1,3
          if (lkappa_arr(i,j,k)) then
            p%kappa_coefs(:,i,j,k)=emf_interpolate(kappa_data(1:dataload_len,:,m-nghost,nind,i,j,k))
          else
            p%kappa_coefs(:,i,j,k)=0
          endif
        end do; end do; end do

        ! Calculate kappa (grad B)_symm
        p%kappa_emf = 0; kappa_mask=0.
        do k=1,3; do j=1,3; do i=1,3
          if (lkappa_arr(i,j,k)) then
            if (lregularize_kappa) then
              if (i==3.and.j+k==3) then
!
! For cases where |B_x;y| << |B_y;x| or |B_y;x| << |B_x;y|: ensure that beta+/-kappa are positive definit.
!
                if (lactive_dimension(3)) then
                  jrt=p%bij(:,1,2)
                  jtr=p%bij(:,2,1)
                else             
                  jrt=p%bijtilde(:,1,2)+p%bij_cov_corr(:,1,2)
                  jtr=p%bijtilde(:,2,1)+p%bij_cov_corr(:,2,1)
                endif

                where (abs(jrt)<jthreshold*abs(jtr) .and. &
                       !!p%kappa_coefs(:,3,j,k)<-eta-p%beta_coefs(:,3,3) )
                       (1.-jthreshold**2)*p%kappa_coefs(:,3,j,k)<-(1.-jthreshold)**2*(eta+p%beta_coefs(:,3,3)) )
                       !!!(jtr**2-jrt**2)*p%kappa_coefs(:,3,j,k)+(jrt-jtr)**2*(eta+p%beta_coefs(:,3,3))<0. )
                  p%kappa_coefs(:,3,j,k)=(-1.+rel_kappa)*(eta+p%beta_coefs(:,3,3))
                  kappa_mask=1.
                elsewhere (abs(jtr)<jthreshold*abs(jrt) .and. &
                           !!p%kappa_coefs(:,3,j,k)>eta+p%beta_coefs(:,3,3))
                           (1.-jthreshold**2)*p%kappa_coefs(:,3,j,k)>(1.-jthreshold)**2*(eta+p%beta_coefs(:,3,3)))
                           !!!(jrt**2-jtr**2)*p%kappa_coefs(:,3,j,k)-(jtr-jrt)**2*(eta+p%beta_coefs(:,3,3))>0.)
                  p%kappa_coefs(:,3,j,k)=(1.-rel_kappa)*(eta+p%beta_coefs(:,3,3))
                  kappa_mask=-1.
                endwhere
                !where ( (eta+p%beta_coefs(:,3,3))*(jtr-jrt)**2 + (jtr**2 - jrt**2)*p%kappa_coefs(:,3,j,k) < 0.)
                !  p%kappa_coefs(:,3,j,k)=0.
                !  kappa_mask=1.
                !endwhere
!if (ldiagnos.and.any(kappa_mask/=0.)) then
if (.false.) then
  mess=''
  !mess='it,iprocs, m='//trim(itoa(it))//' '//trim(itoa(ipx))//' '//trim(itoa(ipy))//' '//trim(itoa(m))
  write(mess,'(e12.4)') y(m)

  do iii=1,nx
    if (kappa_mask(iii)/=0.) then
      mess=trim(mess)//'|'//trim(itoa(iii))
      write(cbuf,'(e12.4)') x(iii)
      write(20+iproc,*) trim(mess)//' '//trim(cbuf)
    endif
  enddo
endif
              endif
            endif
            p%kappa_emf(:,i)=p%kappa_emf(:,i)+p%kappa_coefs(:,i,j,k)*p%bij_symm(:,j,k)
          endif
        end do; end do; end do
        p%emf = p%emf - p%kappa_emf
      end if
!
      if (lumean) then
        ! Calculate Umean x B
        do i=1,3
          if (lumean_arr(i)) then
            p%umean_coefs(:,i)=emf_interpolate(umean_data(1:dataload_len,:,m-nghost,nind,i))
          else
            p%umean_coefs(:,i)=0
          end if
        end do
        call cross_mn(p%umean_coefs,p%bb,p%umean_emf)
        p%emf = p%emf + p%umean_emf
      end if
!
if (.false.) then
if (maxval(abs(p%bcoef_coefs(:,2,2,1)-p%bcoef_coefs(:,2,1,2)-2.*(p%beta_coefs(:,2,3)-p%delta_coefs(:,1))))>1.e-9) &
  print*, '2,3,m,n=', m,n
if (maxval(abs(p%bijtilde(:,2,1)-p%bijtilde(:,1,2)+p%bb(:,2)/x(l1:l2)-p%jj(:,3)))>1.e-10) print*, 'jphi,m=', m, &
maxval(abs(p%jj(:,3))), maxval(abs(p%bijtilde(:,2,1)-p%bijtilde(:,1,2)+p%bb(:,2)/x(l1:l2)))
endif
!
    endsubroutine calc_pencils_special
!***********************************************************************
    subroutine dspecial_dt(f,df,p)
!
!  calculate right hand side of ONE OR MORE extra coupled PDEs
!  along the 'current' Pencil, i.e. f(l1:l2,m,n) where
!  m,n are global variables looped over in equ.f90
!
!  Due to the multi-step Runge Kutta timestepping used one MUST always
!  add to the present contents of the df array.  NEVER reset it to zero.
!
!  Several precalculated Pencils of information are passed for
!  efficiency.
!
!  06-oct-03/tony: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
!
      intent(in) :: f,p
      intent(inout) :: df

      integer :: i, j
!
!  Identify module and boundary conditions.
!
      if (lroot.and.(headtt.or.ldebug)) print*,'dspecial_dt: SOLVE dspecial_dt'
!!      if (headtt) call identify_bcs('special',ispecial)
!
!
!! SAMPLE DIAGNOSTIC IMPLEMENTATION
!!
!!      if (ldiagnos) then
!!        if (idiag_SPECIAL_DIAGNOSTIC/=0) then
!!          call sum_mn_name(MATHEMATICAL EXPRESSION,idiag_SPECIAL_DIAGNOSTIC)
!!! see also integrate_mn_name
!!        endif
!!      endif
!      emftmp=0
!
      if (ldiagnos) then

        tmppencil = emftmp - p%emf
        !
        if (idiag_alphaxmax/=0) call max_mn_name(p%alpha_emf(:,1),idiag_alphaxmax)
        if (idiag_alphaymax/=0) call max_mn_name(p%alpha_emf(:,2),idiag_alphaymax)
        if (idiag_alphazmax/=0) call max_mn_name(p%alpha_emf(:,3),idiag_alphazmax)
        if (idiag_betaxmax/=0) call max_mn_name(p%beta_emf(:,1),idiag_betaxmax)
        if (idiag_betaymax/=0) call max_mn_name(p%beta_emf(:,2),idiag_betaymax)
        if (idiag_betazmax/=0) call max_mn_name(p%beta_emf(:,3),idiag_betazmax)
        if (idiag_gammaxmax/=0) call max_mn_name(p%gamma_emf(:,1),idiag_gammaxmax)
        if (idiag_gammaymax/=0) call max_mn_name(p%gamma_emf(:,2),idiag_gammaymax)
        if (idiag_gammazmax/=0) call max_mn_name(p%gamma_emf(:,3),idiag_gammazmax)
        if (idiag_deltaxmax/=0) call max_mn_name(p%delta_emf(:,1),idiag_deltaxmax)
        if (idiag_deltaymax/=0) call max_mn_name(p%delta_emf(:,2),idiag_deltaymax)
        if (idiag_deltazmax/=0) call max_mn_name(p%delta_emf(:,3),idiag_deltazmax)
        if (idiag_kappaxmax/=0) call max_mn_name(p%kappa_emf(:,1),idiag_kappaxmax)
        if (idiag_kappaymax/=0) call max_mn_name(p%kappa_emf(:,2),idiag_kappaymax)
        if (idiag_kappazmax/=0) call max_mn_name(p%kappa_emf(:,3),idiag_kappazmax)
        if (idiag_umeanxmax/=0) call max_mn_name(p%umean_emf(:,1),idiag_umeanxmax)
        if (idiag_umeanymax/=0) call max_mn_name(p%umean_emf(:,2),idiag_umeanymax)
        if (idiag_umeanzmax/=0) call max_mn_name(p%umean_emf(:,3),idiag_umeanzmax)
        if (idiag_acoefxmax/=0) call max_mn_name(p%acoef_emf(:,1),idiag_acoefxmax)
        if (idiag_acoefymax/=0) call max_mn_name(p%acoef_emf(:,2),idiag_acoefymax)
        if (idiag_acoefzmax/=0) call max_mn_name(p%acoef_emf(:,3),idiag_acoefzmax)
        if (idiag_bcoefxmax/=0) call max_mn_name(p%bcoef_emf(:,1),idiag_bcoefxmax)
        if (idiag_bcoefymax/=0) call max_mn_name(p%bcoef_emf(:,2),idiag_bcoefymax)
        if (idiag_bcoefzmax/=0) call max_mn_name(p%bcoef_emf(:,3),idiag_bcoefzmax)
        if (idiag_emfxmax/=0) call max_mn_name(p%emf(:,1),idiag_emfxmax)
        if (idiag_emfymax/=0) call max_mn_name(p%emf(:,2),idiag_emfymax)
        if (idiag_emfzmax/=0) call max_mn_name(p%emf(:,3),idiag_emfzmax)
        if (idiag_emfcoefxmax/=0) call max_mn_name(emftmp(:,1),idiag_emfcoefxmax)
        if (idiag_emfcoefymax/=0) call max_mn_name(emftmp(:,2),idiag_emfcoefymax)
        if (idiag_emfcoefzmax/=0) call max_mn_name(emftmp(:,3),idiag_emfcoefzmax)
        if (idiag_emfdiffmax/=0) then
          call dot2_mn(tmppencil,tmpline)
          call max_mn_name(tmpline,idiag_emfdiffmax,lsqrt=.true.)
        end if
        if (idiag_emfxdiffmax/=0) &
          call max_mn_name(abs(tmppencil(:,1)),idiag_emfxdiffmax)
        if (idiag_emfydiffmax/=0) &
          call max_mn_name(abs(tmppencil(:,2)),idiag_emfydiffmax)
        if (idiag_emfzdiffmax/=0) &
          call max_mn_name(abs(tmppencil(:,3)),idiag_emfzdiffmax)
        if (idiag_alpharms/=0) then
          call dot2_mn(p%alpha_emf,tmpline)
          call sum_mn_name(tmpline,idiag_alpharms,lsqrt=.true.)
        end if
        if (idiag_betarms/=0) then
          call dot2_mn(p%beta_emf,tmpline)
          call sum_mn_name(tmpline,idiag_betarms,lsqrt=.true.)
        end if
        if (idiag_gammarms/=0) then
          call dot2_mn(p%gamma_emf,tmpline)
          call sum_mn_name(tmpline,idiag_gammarms,lsqrt=.true.)
        end if
        if (idiag_deltarms/=0) then
          call dot2_mn(p%delta_emf,tmpline)
          call sum_mn_name(tmpline,idiag_deltarms,lsqrt=.true.)
        end if
        if (idiag_kapparms/=0) then
          call dot2_mn(p%kappa_emf,tmpline)
          call sum_mn_name(tmpline,idiag_kapparms,lsqrt=.true.)
        end if
        if (idiag_umeanrms/=0) then
          call dot2_mn(p%umean_emf,tmpline)
          call sum_mn_name(tmpline,idiag_umeanrms,lsqrt=.true.)
        end if
        if (idiag_acoefrms/=0) then
          call dot2_mn(p%acoef_emf,tmpline)
          call sum_mn_name(tmpline,idiag_acoefrms,lsqrt=.true.)
        end if
        if (idiag_bcoefrms/=0) then
          call dot2_mn(p%bcoef_emf,tmpline)
          call sum_mn_name(tmpline,idiag_bcoefrms,lsqrt=.true.)
        end if
        if (idiag_emfrms/=0) then
          call dot2_mn(p%emf,tmpline)
          call sum_mn_name(tmpline,idiag_emfrms,lsqrt=.true.)
        end if
        if (idiag_emfcoefrms/=0) then
          call dot2_mn(emftmp,tmpline)
          call sum_mn_name(tmpline,idiag_emfcoefrms,lsqrt=.true.)
        end if
        if (idiag_emfdiffrms/=0) then
          call dot2_mn(tmppencil,tmpline)
          call sum_mn_name(tmpline,idiag_emfdiffrms,lsqrt=.true.)
        end if

        if (idiag_nkappareg/=0) call sum_mn_name(abs(kappa_mask),idiag_nkappareg,lplain=.true.)
      end if
!
      if (l2davgfirst) then
        call zsum_mn_name_xy(p%emf(:,1),idiag_emfxmxy)
        call zsum_mn_name_xy(p%emf(:,2),idiag_emfymxy)
        call zsum_mn_name_xy(p%emf(:,3),idiag_emfzmxy)
        call zsum_mn_name_xy(emftmp(:,1),idiag_emfcoefxmxy)
        call zsum_mn_name_xy(emftmp(:,2),idiag_emfcoefymxy)
        call zsum_mn_name_xy(emftmp(:,3),idiag_emfcoefzmxy)
        call zsum_mn_name_xy(p%alpha_coefs(:,1,1),idiag_alphaxxmxy)
        call zsum_mn_name_xy(p%alpha_coefs(:,1,2),idiag_alphaxymxy)
        call zsum_mn_name_xy(p%alpha_coefs(:,1,3),idiag_alphaxzmxy)
        call zsum_mn_name_xy(p%alpha_coefs(:,2,2),idiag_alphayymxy)
        call zsum_mn_name_xy(p%alpha_coefs(:,2,3),idiag_alphayzmxy)
        call zsum_mn_name_xy(p%alpha_coefs(:,3,3),idiag_alphazzmxy)
        call zsum_mn_name_xy(p%beta_coefs(:,1,1),idiag_betaxxmxy)
        call zsum_mn_name_xy(p%beta_coefs(:,1,2),idiag_betaxymxy)
        call zsum_mn_name_xy(p%beta_coefs(:,1,3),idiag_betaxzmxy)
        call zsum_mn_name_xy(p%beta_coefs(:,2,2),idiag_betayymxy)
        call zsum_mn_name_xy(p%beta_coefs(:,2,3),idiag_betayzmxy)
        call zsum_mn_name_xy(p%beta_coefs(:,3,3),idiag_betazzmxy)
        call zsum_mn_name_xy(p%gamma_coefs(:,1),idiag_gammaxmxy)
        call zsum_mn_name_xy(p%gamma_coefs(:,2),idiag_gammaymxy)
        call zsum_mn_name_xy(p%gamma_coefs(:,3),idiag_gammazmxy)
        call zsum_mn_name_xy(p%delta_coefs(:,1),idiag_deltaxmxy)
        call zsum_mn_name_xy(p%delta_coefs(:,2),idiag_deltaymxy)
        call zsum_mn_name_xy(p%delta_coefs(:,3),idiag_deltazmxy)
        call zsum_mn_name_xy(p%umean_coefs(:,1),idiag_umeanxmxy)
        call zsum_mn_name_xy(p%umean_coefs(:,2),idiag_umeanymxy)
        call zsum_mn_name_xy(p%umean_coefs(:,3),idiag_umeanzmxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,1,1,1),idiag_kappaxxxmxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,2,1,1),idiag_kappayxxmxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,3,1,1),idiag_kappazxxmxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,1,1,2),idiag_kappaxxymxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,2,1,2),idiag_kappayxymxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,3,1,2),idiag_kappazxymxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,1,1,3),idiag_kappaxxzmxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,2,1,3),idiag_kappayxzmxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,3,1,3),idiag_kappazxzmxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,1,2,2),idiag_kappaxyymxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,2,2,2),idiag_kappayyymxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,3,2,2),idiag_kappazyymxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,1,2,3),idiag_kappaxyzmxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,2,2,3),idiag_kappayyzmxy)
        call zsum_mn_name_xy(p%kappa_coefs(:,3,2,3),idiag_kappazyzmxy)
      endif

      call keep_compiler_quiet(f,df)
!
    endsubroutine dspecial_dt
!***********************************************************************
    subroutine read_special_init_pars(iostat)
!
      integer, intent(out) :: iostat
!
      iostat = 0
      call setParameterDefaults
      read(parallel_unit, NML=special_init_pars, IOSTAT=iostat)
      if (lroot) write (*,*) 'read_special_init_pars parameters read...'
      call parseParameters
      if (lroot) write (*,*) 'read_special_init_pars parameters parsed...'
!
    endsubroutine read_special_init_pars
!***********************************************************************
    subroutine write_special_init_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=special_init_pars)

   endsubroutine write_special_init_pars
!***********************************************************************
    subroutine read_special_run_pars(iostat)
!
      integer, intent(out) :: iostat
!
      iostat = 0

      call setParameterDefaults
      if (lroot) write (*,*) 'read_special_run_pars parameters read...'
      read(parallel_unit, NML=special_run_pars, IOSTAT=iostat)
      call parseParameters
      if (lroot) write (*,*) 'read_special_run_pars parameters parsed...'
!
    endsubroutine read_special_run_pars
!***********************************************************************
    subroutine write_special_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=special_run_pars)
!
    endsubroutine write_special_run_pars
!***********************************************************************
    subroutine rprint_special(lreset,lwrite)
!
!  Reads and registers print parameters relevant to special.
!
!  06-oct-03/tony: coded
!
!!      use FArrayManager, only: farray_index_append
!
      integer :: iname
      logical :: lreset,lwr
      logical, optional :: lwrite
!
      lwr = .false.
      if (present(lwrite)) lwr=lwrite
!
!  reset everything in case of reset
!  (this needs to be consistent with what is defined above!)
!
      if (lreset) then
        ! Max diagnostics
        idiag_alphaxmax=0
        idiag_alphaymax=0
        idiag_alphazmax=0
        idiag_betaxmax=0
        idiag_betaymax=0
        idiag_betazmax=0
        idiag_gammaxmax=0
        idiag_gammaymax=0
        idiag_gammazmax=0
        idiag_deltaxmax=0
        idiag_deltaymax=0
        idiag_deltazmax=0
        idiag_kappaxmax=0
        idiag_kappaymax=0
        idiag_kappazmax=0
        idiag_umeanxmax=0
        idiag_umeanymax=0
        idiag_umeanzmax=0
        idiag_acoefxmax=0
        idiag_acoefymax=0
        idiag_acoefzmax=0
        idiag_bcoefxmax=0
        idiag_bcoefymax=0
        idiag_bcoefzmax=0
        idiag_emfxmax=0
        idiag_emfymax=0
        idiag_emfzmax=0
        idiag_emfcoefxmax=0
        idiag_emfcoefymax=0
        idiag_emfcoefzmax=0
        idiag_emfdiffmax=0
        idiag_emfxdiffmax=0
        idiag_emfydiffmax=0
        idiag_emfzdiffmax=0
!     RMS diagnostics
        idiag_alpharms=0
        idiag_betarms=0
        idiag_gammarms=0
        idiag_deltarms=0
        idiag_kapparms=0
        idiag_umeanrms=0
        idiag_acoefrms=0
        idiag_bcoefrms=0
        idiag_emfrms=0
        idiag_emfcoefrms=0
        idiag_emfdiffrms=0
!    timestep diagnostics
        idiag_dtemf_ave=0
        idiag_dtemf_dif=0
        idiag_nkappareg=0
!    2D diagnostics
        idiag_emfxmxy=0; idiag_emfymxy=0; idiag_emfzmxy=0
        idiag_emfcoefxmxy=0; idiag_emfcoefymxy=0; idiag_emfcoefzmxy=0
        idiag_alphaxxmxy=0; idiag_alphayymxy=0; idiag_alphazzmxy=0;
        idiag_alphaxymxy=0; idiag_alphaxzmxy=0; idiag_alphayzmxy=0;
        idiag_betaxxmxy=0; idiag_betayymxy=0; idiag_betazzmxy=0;
        idiag_betaxymxy=0; idiag_betaxzmxy=0; idiag_betayzmxy=0;
        idiag_gammaxmxy=0; idiag_gammaymxy=0; idiag_gammazmxy=0;
        idiag_deltaxmxy=0; idiag_deltaymxy=0; idiag_deltazmxy=0;
        idiag_umeanxmxy=0; idiag_umeanymxy=0; idiag_umeanzmxy=0;
        idiag_kappaxxxmxy=0; idiag_kappayxxmxy=0; idiag_kappazxxmxy=0;
        idiag_kappaxxymxy=0; idiag_kappayxymxy=0; idiag_kappazxymxy=0;
        idiag_kappaxxzmxy=0; idiag_kappayxzmxy=0; idiag_kappazxzmxy=0;
        idiag_kappaxyymxy=0; idiag_kappayyymxy=0; idiag_kappazyymxy=0;
        idiag_kappaxyzmxy=0; idiag_kappayyzmxy=0; idiag_kappazyzmxy=0;
      endif
!
      do iname=1,nname
        ! Maximum values of emf terms
        call parse_name(iname,cname(iname),cform(iname),'alphaxmax',idiag_alphaxmax)
        call parse_name(iname,cname(iname),cform(iname),'alphaymax',idiag_alphaymax)
        call parse_name(iname,cname(iname),cform(iname),'alphazmax',idiag_alphazmax)
        call parse_name(iname,cname(iname),cform(iname),'betaxmax',idiag_betaxmax)
        call parse_name(iname,cname(iname),cform(iname),'betaymax',idiag_betaymax)
        call parse_name(iname,cname(iname),cform(iname),'betazmax',idiag_betazmax)
        call parse_name(iname,cname(iname),cform(iname),'gammaxmax',idiag_gammaxmax)
        call parse_name(iname,cname(iname),cform(iname),'gammaymax',idiag_gammaymax)
        call parse_name(iname,cname(iname),cform(iname),'gammazmax',idiag_gammazmax)
        call parse_name(iname,cname(iname),cform(iname),'deltaxmax',idiag_deltaxmax)
        call parse_name(iname,cname(iname),cform(iname),'deltaymax',idiag_deltaymax)
        call parse_name(iname,cname(iname),cform(iname),'deltazmax',idiag_deltazmax)
        call parse_name(iname,cname(iname),cform(iname),'kappaxmax',idiag_kappaxmax)
        call parse_name(iname,cname(iname),cform(iname),'kappaymax',idiag_kappaymax)
        call parse_name(iname,cname(iname),cform(iname),'kappazmax',idiag_kappazmax)
        call parse_name(iname,cname(iname),cform(iname),'umeanxmax',idiag_umeanxmax)
        call parse_name(iname,cname(iname),cform(iname),'umeanymax',idiag_umeanymax)
        call parse_name(iname,cname(iname),cform(iname),'umeanzmax',idiag_umeanzmax)
        call parse_name(iname,cname(iname),cform(iname),'acoefxmax',idiag_acoefxmax)
        call parse_name(iname,cname(iname),cform(iname),'acoefymax',idiag_acoefymax)
        call parse_name(iname,cname(iname),cform(iname),'acoefzmax',idiag_acoefzmax)
        call parse_name(iname,cname(iname),cform(iname),'bcoefxmax',idiag_bcoefxmax)
        call parse_name(iname,cname(iname),cform(iname),'bcoefymax',idiag_bcoefymax)
        call parse_name(iname,cname(iname),cform(iname),'bcoefzmax',idiag_bcoefzmax)
        call parse_name(iname,cname(iname),cform(iname),'emfxmax',idiag_emfxmax)
        call parse_name(iname,cname(iname),cform(iname),'emfymax',idiag_emfymax)
        call parse_name(iname,cname(iname),cform(iname),'emfzmax',idiag_emfzmax)
        call parse_name(iname,cname(iname),cform(iname),'emfcoefxmax',idiag_emfcoefxmax)
        call parse_name(iname,cname(iname),cform(iname),'emfcoefymax',idiag_emfcoefymax)
        call parse_name(iname,cname(iname),cform(iname),'emfcoefzmax',idiag_emfcoefzmax)
        call parse_name(iname,cname(iname),cform(iname),'emfdiffmax',idiag_emfdiffmax)
        call parse_name(iname,cname(iname),cform(iname),'emfxdiffmax',idiag_emfxdiffmax)
        call parse_name(iname,cname(iname),cform(iname),'emfydiffmax',idiag_emfydiffmax)
        call parse_name(iname,cname(iname),cform(iname),'emfzdiffmax',idiag_emfzdiffmax)
        ! RMS values of emf terms
        call parse_name(iname,cname(iname),cform(iname),'alpharms',idiag_alpharms)
        call parse_name(iname,cname(iname),cform(iname),'betarms',idiag_betarms)
        call parse_name(iname,cname(iname),cform(iname),'gammarms',idiag_gammarms)
        call parse_name(iname,cname(iname),cform(iname),'deltarms',idiag_deltarms)
        call parse_name(iname,cname(iname),cform(iname),'kapparms',idiag_kapparms)
        call parse_name(iname,cname(iname),cform(iname),'umeanrms',idiag_umeanrms)
        call parse_name(iname,cname(iname),cform(iname),'acoefrms',idiag_acoefrms)
        call parse_name(iname,cname(iname),cform(iname),'bcoefrms',idiag_bcoefrms)
        call parse_name(iname,cname(iname),cform(iname),'emfrms',idiag_emfrms)
        call parse_name(iname,cname(iname),cform(iname),'emfcoefrms',idiag_emfcoefrms)
        call parse_name(iname,cname(iname),cform(iname),'emfdiffrms',idiag_emfdiffrms)
!    timestep diagnostics
        call parse_name(iname,cname(iname),cform(iname),'dtemf_ave',idiag_dtemf_ave)
        call parse_name(iname,cname(iname),cform(iname),'dtemf_dif',idiag_dtemf_dif)
        call parse_name(iname,cname(iname),cform(iname),'nkappareg',idiag_nkappareg)
      enddo

      do iname=1,nnamexy
        call parse_name(iname,cnamexy(iname),cformxy(iname),'EMFxmxy',idiag_emfxmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'EMFymxy',idiag_emfymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'EMFzmxy',idiag_emfzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'EMFcoefxmxy',idiag_emfcoefxmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'EMFcoefymxy',idiag_emfcoefymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'EMFcoefzmxy',idiag_emfcoefzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'alphaxxmxy',idiag_alphaxxmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'alphayymxy',idiag_alphayymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'alphazzmxy',idiag_alphazzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'alphaxymxy',idiag_alphaxymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'alphaxzmxy',idiag_alphaxzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'alphayzmxy',idiag_alphayzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'betaxxmxy',idiag_betaxxmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'betayymxy',idiag_betayymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'betazzmxy',idiag_betazzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'betaxymxy',idiag_betaxymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'betaxzmxy',idiag_betaxzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'betayzmxy',idiag_betayzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'umeanxmxy',idiag_umeanxmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'umeanymxy',idiag_umeanymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'umeanzmxy',idiag_umeanzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'deltaxmxy',idiag_deltaxmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'deltaymxy',idiag_deltaymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'deltazmxy',idiag_deltazmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'gammaxmxy',idiag_gammaxmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'gammaymxy',idiag_gammaymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'gammazmxy',idiag_gammazmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappaxxxmxy',idiag_kappaxxxmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappayxxmxy',idiag_kappayxxmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappazxxmxy',idiag_kappazxxmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappaxxymxy',idiag_kappaxxymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappayxymxy',idiag_kappayxymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappazxymxy',idiag_kappazxymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappaxxzmxy',idiag_kappaxxzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappayxzmxy',idiag_kappayxzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappazxzmxy',idiag_kappazxzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappaxyymxy',idiag_kappaxyymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappayyymxy',idiag_kappayyymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappazyymxy',idiag_kappazyymxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappaxyzmxy',idiag_kappaxyzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappayyzmxy',idiag_kappayyzmxy)
        call parse_name(iname,cnamexy(iname),cformxy(iname),'kappazyzmxy',idiag_kappazyzmxy)
      enddo

    endsubroutine rprint_special
!***********************************************************************
    subroutine special_calc_magnetic(f,df,p)
!
!  Calculate an additional 'special' term on the right hand side of the
!  induction equation.
!
!  Some precalculated pencils of data are passed in for efficiency
!  others may be calculated directly from the f array.
!
!  06-oct-03/tony: coded
!
      real, dimension (mx,my,mz,mfarray), intent(in) :: f
      real, dimension (mx,my,mz,mvar), intent(inout) :: df
      real :: diffus_tmp
      type (pencil_case), intent(in) :: p
      integer :: i,j,k
!
      call keep_compiler_quiet(f)
!
! Overwrite with a and b coefs if needed
!
      if (lusecoefs) then
        emftmp=0
        if (lacoef) emftmp = emftmp + p%acoef_emf
        if (lbcoef) emftmp = emftmp + p%bcoef_emf
        if (lumean) emftmp = emftmp + p%umean_emf
      else
        emftmp = p%emf
      end if

      df(l1:l2,m,n,iax:iaz)=df(l1:l2,m,n,iax:iaz)+emftmp
!
      if (lfirst.and.ldt) then
!
! Calculate advec_special
!
        advec_special=0.
        if (lalpha) then
          call dot_mn_vm(dline_1, abs(p%alpha_coefs), tmppencil)
          advec_special=advec_special+sum(tmppencil,2)
        end if
        if (lgamma) then
          call dot_mn(dline_1, abs(p%gamma_coefs), tmpline)
          advec_special=advec_special+tmpline
        end if
        if (lumean) then
          call dot_mn(dline_1, abs(p%umean_coefs), tmpline)
          advec_special=advec_special+tmpline
        end if
!
        maxadvec=maxadvec+advec_special
!
        if (ldiagnos.and.idiag_dtemf_ave/=0) then
          call max_mn_name(advec_special/cdt,idiag_dtemf_ave,l_dt=.true.)
        endif
!
! Calculate diffus_special
!
        diffus_special=0.
        if (lbeta) then
          call dot_mn_vm(dline_1,abs(p%beta_coefs), tmppencil)
          call dot_mn(dline_1, tmppencil, tmpline)
          diffus_special=diffus_special+tmpline
        end if
!
        if (ldelta) then
          call cross_mn(dline_1,abs(p%delta_coefs), tmppencil)
          call dot_mn(dline_1,tmppencil,tmpline)
          diffus_special=diffus_special+tmpline
        end if
!
        if (lkappa) then
          call vec_dot_3tensor(dline_1, abs(p%kappa_coefs), tmptensor)
          call dot_mn_vm(dline_1,tmptensor, tmppencil)
          diffus_special=diffus_special+sum(tmppencil,2)
        end if

        maxdiffus=max(maxdiffus,diffus_special)
!
        if (ldiagnos.and.idiag_dtemf_dif/=0) then
          call max_mn_name(diffus_special/cdtv,idiag_dtemf_dif,l_dt=.true.)
        endif
!
      end if
!
    endsubroutine special_calc_magnetic
!***********************************************************************
    subroutine openDataset(datagroup_,tensor_id)

      ! Open a dataset e.g. /emftensor/alpha/data and auxillary dataspaces

      character(len=*), dimension(:), intent(in) :: datagroup_     ! name of data group
      integer, intent(in)             :: tensor_id
!
      integer(HSIZE_T), dimension(10) :: dimsizes, maxdimsizes
      integer :: ndims, i
      integer(HSIZE_T) :: num
      logical :: hdf_exists, lok
      character(len=fnlen)            :: dataset
      character(len=len(datagroup_))  :: datagroup

      dataset = tensor_names(tensor_id)               ! name of dataset.
      ! Check that datagroup e.g. /emftensor/alpha exists

      lok=.false.
      do i=1,size(datagroup_)
        datagroup=datagroup_(i)
        call H5Lexists_F(hdf_emftensors_group, datagroup, hdf_exists, hdferr)
        if (hdf_exists) then
          lok=.true.
          exit
        endif
      enddo
      if (.not. lok) &
          call fatal_error('openDataset','/emftensor/'//trim(datagroup)// ' does not exist')

      ! Open datagroup, returns identifier in tensor_id_G.
      call H5Gopen_F(hdf_emftensors_group, datagroup, tensor_id_G(tensor_id),hdferr)
      if (hdferr /= 0) then
        call fatal_error('openDataset','Error opening /emftensor/'//trim(datagroup))
      end if
      ! Check that dataset e.g. /emftensor/alpha/mean exists
      call H5Lexists_F(tensor_id_G(tensor_id), dataset, hdf_exists, hdferr)
      if (.not. hdf_exists) then
        call H5Gclose_F(tensor_id_G(tensor_id), hdferr)
        call fatal_error('openDataset','/emftensor/'//trim(datagroup)// &
                          '/'//trim(dataset)//' does not exist')
      end if
      ! Open dataset <dataset> in group, returns identifier in tensor_id_D.
      call H5Dopen_F(tensor_id_G(tensor_id), dataset, tensor_id_D(tensor_id),hdferr)
      if (hdferr /= 0) then
        call H5Gclose_F(tensor_id_G(tensor_id), hdferr)
        call fatal_error('openDataset','Error opening /emftensor/'// &
                          trim(datagroup)//'/'//trim(dataset))
      end if
      ! Get dataspace of dataset (identifier of copy of dataspace in tensor_id_S).
      call H5Dget_space_F(tensor_id_D(tensor_id), tensor_id_S(tensor_id),hdferr)
      if (hdferr /= 0) then
        call H5Dclose_F(tensor_id_D(tensor_id), hdferr)
        call H5Gclose_F(tensor_id_G(tensor_id), hdferr)
        call fatal_error('openDataset','Error opening dataspace for/emftensor/'// &
                          trim(datagroup)//'/'//trim(dataset))
      end if
      ! Get dataspace dimensions in tensor_dims.
      ndims = tensor_ndims(tensor_id)
      call H5Sget_simple_extent_dims_F(tensor_id_S(tensor_id), &
                                       dimsizes(1:ndims), &
                                       maxdimsizes(1:ndims), &
                                       hdferr)                 !MR: hdferr/=0!
!print*, 'from H5Sget_simple_extent_dims_F, line 1330: hdferr=', hdferr
      call H5Sget_simple_extent_npoints_F(tensor_id_S(tensor_id),num,hdferr) ! This is to mask the error of the preceding call.

      tensor_dims(tensor_id,1:ndims)=dimsizes(1:ndims)
      if (tensor_dims(tensor_id,2)/=nxgrid) &
        call fatal_error('initialize_special', &
        'x extent in dataset "'//trim(datagroup)//'" '//trim(itoa(int(tensor_dims(tensor_id,2))))//' different from nxgrid')
      if (tensor_dims(tensor_id,3)/=nygrid) &
        call fatal_error('initialize_special', &
        'y extent in dataset "'//trim(datagroup)//'" '//trim(itoa(int(tensor_dims(tensor_id,3))))//' different from nygrid')

      nx_stored=tensor_dims(tensor_id,2)/nprocx
      ny_stored=tensor_dims(tensor_id,3)/nprocy

      ! Set dataset counts, memory offsets, counts and dimensions

      tensor_counts(tensor_id,1:4) = [ dataload_len, nx_stored, ny_stored, nz_stored ]
      tensor_memcounts(tensor_id,1:4) = [ dataload_len, nx_stored, ny_stored, nz_stored ]
      tensor_memdims(tensor_id,1:4) = [ dataload_len, nx_stored, ny_stored, nz_stored ]
      if (nz_stored == 1) then
        tensor_offsets(tensor_id,1:4) = [ 0, ipx*nx_stored, ipy*ny_stored, 0 ]
      else
        tensor_offsets(tensor_id,1:4) = [ 0, ipx*nx_stored, ipy*ny_stored, ipz*nz_stored ]
      endif

      if (tensor_times_len==-1) then
        tensor_times_len=dimsizes(1)
      elseif (tensor_times_len/=dimsizes(1)) then
        call fatal_error('openDataset','dataset emftensor/'//trim(datagroup)//'/'//trim(dataset)//' has deviating time extent')
      endif

      if (hdferr /= 0) then
        call H5Sclose_F(tensor_id_S(tensor_id), hdferr)
        call H5Dclose_F(tensor_id_D(tensor_id), hdferr)
        call H5Gclose_F(tensor_id_G(tensor_id), hdferr)
        call fatal_error('openDataset','Cannot get dimensions extent '// &
                          'for /emftensor/'//trim(datagroup)//'/'//trim(dataset))
      end if

      ! Create a memory space mapping for input data (identifier in tensor_id_memS).
      call H5Screate_simple_F(ndims, tensor_memdims(tensor_id,1:ndims), &
                              tensor_id_memS(tensor_id), hdferr)
      if (hdferr /= 0) then
        call H5Sclose_F(tensor_id_S(tensor_id), hdferr)
        call H5Dclose_F(tensor_id_D(tensor_id), hdferr)
        call H5Gclose_F(tensor_id_G(tensor_id), hdferr)
        call fatal_error('openDataset','Error creating memory mapping '// &
                          'for /emftensor/'//trim(datagroup)//'/'//trim(dataset))
      end if

      if (lroot) write(*,*) 'Using dataset /emftensor/'//trim(datagroup)//'/'//trim(dataset)//'.'

    end subroutine openDataset
!***********************************************************************
    subroutine openDataset_grid(scalar_id)

      ! Open a dataset in group grid, e.g. /grid/t

      integer, intent(in)             :: scalar_id
!
      integer(HSIZE_T), dimension(10) :: maxdimsizes
      integer(HSIZE_T) :: num
      integer :: ndims
      logical :: hdf_exists
      character(len=fnlen)            :: dataset

      dataset = scalar_names(scalar_id)               ! name of dataset

      ! Check that dataset exists
      call H5Lexists_F(hdf_grid_group, dataset, hdf_exists, hdferr)
      if (.not. hdf_exists) then
        call H5Gclose_F(hdf_grid_group, hdferr)
        call fatal_error('openDataset_grid','/grid/'//trim(dataset)//' does not exist')
      end if
      ! Open dataset dataset in group, returns identifier in scalar_id_D.
      call H5Dopen_F(hdf_grid_group, dataset, scalar_id_D(scalar_id),hdferr)
      if (hdferr /= 0) then
        call H5Gclose_F(hdf_grid_group, hdferr)
        call fatal_error('openDataset_grid','Error opening /grid/'//trim(dataset))
      end if
      ! Get dataspace of dataset (identifier of copy of dataspace in scalar_id_S).
      call H5Dget_space_F(scalar_id_D(scalar_id), scalar_id_S(scalar_id),hdferr)
      if (hdferr /= 0) then
        call H5Dclose_F(scalar_id_D(scalar_id), hdferr)
        call H5Gclose_F(hdf_grid_group, hdferr)
        call fatal_error('openDataset_grid','Error opening dataspace for/grid/'//trim(dataset))
      end if

      ! Get dataspace dimensions in scalar_dims.
      call H5Sget_simple_extent_dims_F(scalar_id_S(scalar_id), &
                                       scalar_dims(scalar_id), &
                                       maxdimsizes(1:ndims), &
                                       hdferr)                 !MR: hdferr/=0!
!This is to mask the error of the preceding call.
      call H5Sget_simple_extent_npoints_F(scalar_id_S(scalar_id),num,hdferr)

      if (hdferr /= 0) then
        call H5Sclose_F(scalar_id_S(scalar_id), hdferr)
        call H5Dclose_F(scalar_id_D(scalar_id), hdferr)
        call H5Gclose_F(hdf_grid_group, hdferr)
        call fatal_error('openDataset_grid','Error in getting dimensions for grid/'//trim(dataset))
      end if

    end subroutine openDataset_grid
!***********************************************************************
    subroutine closeDataset_grid(id)

      ! Close opened dataspaces, dataset and group

      integer :: id

      call H5Sclose_F(scalar_id_S(id), hdferr)
      call H5Dclose_F(scalar_id_D(id), hdferr)

    end subroutine closeDataset_grid
!***********************************************************************
    subroutine closeDataset(tensor_id)

      ! Close opened dataspaces, dataset and group

      integer :: tensor_id

      call H5Sclose_F(tensor_id_memS(tensor_id), hdferr)
      call H5Sclose_F(tensor_id_S(tensor_id), hdferr)
      call H5Dclose_F(tensor_id_D(tensor_id), hdferr)
      call H5Gclose_F(tensor_id_G(tensor_id), hdferr)

    end subroutine closeDataset
!***********************************************************************
    subroutine loadDataset_rank1(dataarray, datamask, tensor_id, loadstart,name)

      ! Load a chunk of data for a vector, beginning at loadstart

      use General, only: itoa

      real, dimension(:,:,:,:,:), intent(inout) :: dataarray
      logical, dimension(3), intent(in) :: datamask
      integer, intent(in) :: tensor_id
      integer, intent(in) :: loadstart
      character(*), optional :: name

      integer :: ndims,i
      integer(HSIZE_T) :: mask_i
      real :: globmin, globmax, sum, rms ! output for diagnostics

      ndims = tensor_ndims(tensor_id)
      tensor_offsets(tensor_id,1) = loadstart
      call H5Sselect_none_F(tensor_id_S(tensor_id), hdferr)     ! resets selection region.
      call H5Sselect_none_F(tensor_id_memS(tensor_id), hdferr)  ! perhaps dispensable when
                                                                ! H5S_SELECT_SET_F is used below.
      do mask_i=1,3
        ! Load only wanted datasets
        if (datamask(mask_i)) then
          ! Set the new offset for data reading
          tensor_offsets(tensor_id,ndims)    = mask_i-1
          tensor_memoffsets(tensor_id,ndims) = mask_i-1
          ! Select hyperslab for data.
!          print '(a,a,5(1x,i3))', 'tensor offset',name,tensor_offsets(tensor_id,1:ndims)
!          print '(a,a,5(1x,i3))', 'tensor memoffset',name,tensor_memoffsets(tensor_id,1:ndims)
!          print '(a,a,5(1x,i3))', 'tensor counts',name,tensor_counts(tensor_id,:ndims)
!          print '(a,a,5(1x,i3))', 'tensor memcounts',name,tensor_memcounts(tensor_id,:ndims)

          call H5Sselect_hyperslab_F(tensor_id_S(tensor_id), H5S_SELECT_OR_F, &
                                     tensor_offsets(tensor_id,1:ndims),       &
                                     tensor_counts(tensor_id,1:ndims),        &
                                     hdferr)
           if (hdferr /= 0) then
             call fatal_error('loadDataset_rank1','Error creating File mapping '// &
                           'for /grid/'//name)
           end if
          ! Select hyperslab for memory.
          call H5Sselect_hyperslab_F(tensor_id_memS(tensor_id), H5S_SELECT_OR_F, &
                                     tensor_memoffsets(tensor_id,1:ndims),       &
                                     tensor_memcounts(tensor_id,1:ndims),        &
                                     hdferr)
          if (hdferr /= 0) then
            call fatal_error('loadDataset_rank1','Error creating memory mapping '// &
                           'for /grid/'//name)
          end if
        end if
      end do
      ! Read data into memory.
      tensor_dims(tensor_id,ndims)=3
      call H5Dread_F(tensor_id_D(tensor_id), hdf_memtype, dataarray, &
                     tensor_dims(tensor_id,1:ndims), hdferr, &
                     tensor_id_memS(tensor_id), tensor_id_S(tensor_id))
      if (hdferr /= 0) &
        call fatal_error('loadDataset_rank1','Error reading dataset '// &
                         'for /grid/'//name)
      sum = 0.; rms = 0.
      do i=1,3
        tensor_maxvals(tensor_id) = maxval(dataarray(:,:,:,:,i))
        tensor_minvals(tensor_id) = minval(dataarray(:,:,:,:,i))
        if (present(name)) then
          call mpireduce_min(tensor_minvals(tensor_id),globmin)
          call mpireduce_max(tensor_maxvals(tensor_id),globmax)
          if (lroot) write (*,*) trim(name)//'(', trim(itoa(i)), ')  min/max:', globmin, globmax
        endif
      sum = sum + tensor_maxvals(tensor_id)*tensor_maxvals(tensor_id)
      enddo
      rms = sqrt(sum)
      if (present(name) .and. lroot) write (*,*) trim(name)//' (rms):', rms

      dataarray = tensor_scales(tensor_id) * dataarray

    end subroutine loadDataset_rank1
!***********************************************************************
    subroutine loadDataset_rank2(dataarray, datamask, tensor_id, loadstart,name)

      ! Load a chunk of data for a 2-rank tensor, beginning at loadstart

      use General, only: itoa

      real, dimension(:,:,:,:,:,:), intent(inout) :: dataarray
      logical, dimension(3,3), intent(in) :: datamask
      integer, intent(in) :: tensor_id
      integer, intent(in) :: loadstart
      character(*), optional :: name

      integer :: ndims, i, j
      integer(HSIZE_T) :: mask_i, mask_j
      real :: globmin, globmax, sum, rms ! output for diagnostics

      ndims = tensor_ndims(tensor_id)
      tensor_offsets(tensor_id,1) = loadstart
      call H5Sselect_none_F(tensor_id_S(tensor_id), hdferr)
      call H5Sselect_none_F(tensor_id_memS(tensor_id), hdferr)

      do mask_j=1,3; do mask_i=1,3

        ! Load only wanted datasets
        if (datamask(mask_i,mask_j)) then

          ! Set the new offset for data reading
          tensor_offsets(tensor_id,ndims-1)    = mask_i-1
          tensor_offsets(tensor_id,ndims)      = mask_j-1
          tensor_memoffsets(tensor_id,ndims-1) = mask_i-1
          tensor_memoffsets(tensor_id,ndims)   = mask_j-1

          ! Hyperslab for data
          !if (mask_i==1.and.mask_j==1) print '(a,i2,a,6(1x,i3))', 'iproc, tensor offset',iproc,name,tensor_offsets(tensor_id,1:ndims-2)
!          print '(a,i2,a,6(1x,i3))', 'iproc, tensor memoffset',iproc,name,tensor_memoffsets(tensor_id,1:ndims)
!          print '(a,a,6(1x,i3))', 'tensor counts',name,tensor_counts(tensor_id,:ndims)
!          print '(a,a,6(1x,i3))', 'tensor memcounts',name,tensor_memcounts(tensor_id,:ndims)

!print*, 'before H5Sselect_hyperslab_F for file'
          call H5Sselect_hyperslab_F(tensor_id_S(tensor_id), H5S_SELECT_OR_F, &
                                     tensor_offsets(tensor_id,1:ndims),       &
                                     tensor_counts(tensor_id,1:ndims),        &
                                     hdferr)
          if (hdferr /= 0) &
            call fatal_error('loadDataset_rank2','Error creating File mapping '// &
                             'for /emftensor/'//name)
          ! Hyperslab for memory
!print*, 'before H5Sselect_hyperslab_F for memory'
          call H5Sselect_hyperslab_F(tensor_id_memS(tensor_id), H5S_SELECT_OR_F, &
                                     tensor_memoffsets(tensor_id,1:ndims),        &
                                     tensor_memcounts(tensor_id,1:ndims),         &
                                     hdferr)
           if (hdferr /= 0) then
             call fatal_error('loadDataset_rank2','Error creating memory mapping '// &
                           'for /emftensor/'//name)
           end if
        end if
      end do; end do


      ! Read data into memory
      tensor_dims(tensor_id,ndims-1:ndims)=3
!print*, 'before H5Dread_F'
      call H5Dread_F(tensor_id_D(tensor_id), hdf_memtype, dataarray, &
                     tensor_dims(tensor_id,1:ndims), hdferr, &
                     tensor_id_memS(tensor_id), tensor_id_S(tensor_id))
      if (hdferr /= 0) &
        call fatal_error('loadDataset_rank2','Error reading dataset '// &
                         'for /emftensor/'//name)
      sum = 0.; rms = 0.
      do i=1,3 ; do j=1,3
        tensor_maxvals(tensor_id) = maxval(dataarray(:,:,:,:,i,j))
        tensor_minvals(tensor_id) = minval(dataarray(:,:,:,:,i,j))
        if (present(name)) then
          call mpireduce_min(tensor_minvals(tensor_id),globmin)
          call mpireduce_max(tensor_maxvals(tensor_id),globmax)
          if (i == j .and. lroot) write (*,*) trim(name)// &
            '(', trim(itoa(i)), ',', trim(itoa(j)),') min/max: ', globmin, globmax
        endif
        sum = sum + tensor_maxvals(tensor_id)*tensor_maxvals(tensor_id)
      enddo ; enddo
      rms=sqrt(sum)
      if (present(name) .and. lroot) write (*,*) trim(name)//' (rms):', rms

      dataarray = tensor_scales(tensor_id) * dataarray

    end subroutine loadDataset_rank2
!***********************************************************************
    subroutine loadDataset_rank3(dataarray, datamask, tensor_id, loadstart,name)

      ! Load a chunk of data for a 3-rank tensor, beginning at loadstart

      real, dimension(:,:,:,:,:,:,:), intent(inout) :: dataarray
      logical, dimension(3,3,3), intent(in) :: datamask
      integer, intent(in) :: tensor_id
      integer, intent(in) :: loadstart
      character(*), optional :: name

      integer :: ndims
      integer(HSIZE_T) :: mask_i, mask_j, mask_k
      real :: globmin, globmax ! output for diagnostics

      ndims = tensor_ndims(tensor_id)
      tensor_offsets(tensor_id,1) = loadstart
      call H5Sselect_none_F(tensor_id_S(tensor_id), hdferr)
      call H5Sselect_none_F(tensor_id_memS(tensor_id), hdferr)

      do mask_k=1,3; do mask_j=1,3; do mask_i=1,3
        ! Load only wanted datasets
        if (datamask(mask_i,mask_j,mask_k)) then
          ! Set the new offset for data reading
          tensor_offsets(tensor_id,ndims-2)    = mask_i-1
          tensor_offsets(tensor_id,ndims-1)    = mask_j-1
          tensor_offsets(tensor_id,ndims)      = mask_k-1
          tensor_memoffsets(tensor_id,ndims-2) = mask_i-1
          tensor_memoffsets(tensor_id,ndims-1) = mask_j-1
          tensor_memoffsets(tensor_id,ndims)   = mask_k-1
          ! Hyperslab for data
!          print '(a,a,7(1x,i3))', 'tensor offset',name,tensor_offsets(tensor_id,1:ndims)
!          print '(a,a,7(1x,i3))', 'tensor memoffset',name,tensor_memoffsets(tensor_id,1:ndims)
!          print '(a,a,7(1x,i3))', 'tensor counts',name,tensor_counts(tensor_id,:ndims)
!          print '(a,a,7(1x,i3))', 'tensor memcounts',name,tensor_memcounts(tensor_id,:ndims)

          call H5Sselect_hyperslab_F(tensor_id_S(tensor_id), H5S_SELECT_OR_F, &
                                     tensor_offsets(tensor_id,1:ndims),       &
                                     tensor_counts(tensor_id,1:ndims),        &
                                     hdferr)
          if (hdferr /= 0) &
            call fatal_error('loadDataset_rank3','Error creating File mapping '// &
                             'for /emftensor/'//name)
          ! Hyperslab for memory
          call H5Sselect_hyperslab_F(tensor_id_memS(tensor_id), H5S_SELECT_OR_F, &
                                     tensor_memoffsets(tensor_id,1:ndims),        &
                                     tensor_memcounts(tensor_id,1:ndims),         &
                                     hdferr)
           if (hdferr /= 0) &
             call fatal_error('loadDataset_rank3','Error creating memory mapping '// &
                              'for /emftensor/'//name)
        end if
      end do; end do; end do
      ! Read data into memory
      tensor_dims(tensor_id,ndims-1:ndims)=3
      call H5Dread_F(tensor_id_D(tensor_id), hdf_memtype, dataarray, &
                     tensor_dims(tensor_id,1:ndims), hdferr, &
                     tensor_id_memS(tensor_id), tensor_id_S(tensor_id))
      if (hdferr /= 0) &
        call fatal_error('loadDataset_rank3','Error reading dataset '// &
                         'for /emftensor/'//name)
      tensor_maxvals(tensor_id) = maxval(dataarray)
      tensor_minvals(tensor_id) = minval(dataarray)
      if (present(name)) then
        call mpireduce_min(tensor_minvals(tensor_id),globmin)
        call mpireduce_max(tensor_maxvals(tensor_id),globmax)
        if (lroot) write (*,*) trim(name)//' min/max: ', globmin, globmax
      endif
      dataarray = tensor_scales(tensor_id) * dataarray

    end subroutine loadDataset_rank3
!***********************************************************************
    function emf_interpolate(dataarray) result(interp_data)

      real, intent(in), dimension(dataload_len,nx) :: dataarray
      real, dimension(nx) :: interp_data

     ! interp_data=dataarray(:,1)
      interp_data=dataarray(1,:)

    end function emf_interpolate
!***********************************************************************
    subroutine setParameterDefaults

      ! alpha
      lalpha=.false.
      lalpha_c=.false.
      lalpha_arr=.false.
      alpha_scale=1.0
      alpha_name='mean'
      ! beta
      lbeta=.false.
      lbeta_c=.false.
      lbeta_arr=.false.
      beta_scale=1.0
      beta_name='mean'
      ! gamma
      lgamma=.false.
      lgamma_c=.false.
      lgamma_arr=.false.
      gamma_scale=1.0
      gamma_name='mean'
      ! delta
      ldelta=.false.
      ldelta_c=.false.
      ldelta_arr=.false.
      delta_scale=1.0
      delta_name='mean'
      ! kappa
      lkappa=.false.
      lkappa_c=.false.
      lkappa_arr=.false.
      kappa_scale=1.0
      kappa_name='mean'
      ! umean
      lumean=.false.
      lumean_c=.false.
      lumean_arr=.false.
      umean_scale=1.0
      umean_name='mean'
      ! acoef
      lacoef=.false.
      lacoef_c=.false.
      lacoef_arr=.false.
      acoef_scale=1.0
      acoef_name='mean'
      ! bcoef
      lbcoef=.false.
      lbcoef_c=.false.
      lbcoef_arr=.false.
      bcoef_scale=1.0
      bcoef_name='mean'
      ! other
      interpname  = 'none'
      dataset = ''
      dataset_type = 'mean'
      dataset_name = ''
      tensor_maxvals=0.0
      tensor_minvals=0.0
      lusecoefs = .false.
      lloop = .false.

    end subroutine setParameterDefaults
!***********************************************************************
    subroutine parseParameters

    integer :: i
!
! Load boolean array for alpha
!
      if (any(lalpha_c)) then
        lalpha = .true.
        lalpha_arr(1,1) = lalpha_c(1)
        lalpha_arr(2,1) = lalpha_c(2)
        lalpha_arr(1,2) = lalpha_c(2)
        lalpha_arr(3,1) = lalpha_c(3)
        lalpha_arr(1,3) = lalpha_c(3)
        lalpha_arr(2,2) = lalpha_c(4)
        lalpha_arr(2,3) = lalpha_c(5)
        lalpha_arr(3,2) = lalpha_c(5)
        lalpha_arr(3,3) = lalpha_c(6)
      else if (lalpha) then
        lalpha_arr = .true.
      end if
!
! Load boolean array for beta
!
      if (any(lbeta_c)) then
        lbeta = .true.
        lbeta_arr(1,1) = lbeta_c(1)
        lbeta_arr(2,1) = lbeta_c(2)
        lbeta_arr(1,2) = lbeta_c(2)
        lbeta_arr(3,1) = lbeta_c(3)
        lbeta_arr(1,3) = lbeta_c(3)
        lbeta_arr(2,2) = lbeta_c(4)
        lbeta_arr(2,3) = lbeta_c(5)
        lbeta_arr(3,2) = lbeta_c(5)
        lbeta_arr(3,3) = lbeta_c(6)
      else if (lbeta) then
        lbeta_arr = .true.
      end if
!
! Load boolean array for gamma
!
      if (any(lgamma_c)) then
        lgamma = .true.
        lgamma_arr  = lgamma_c
      else if (lgamma) then
        lgamma_arr  = .true.
      end if
!
! Load boolean array for delta
!
      if (any(ldelta_c)) then
        ldelta = .true.
        ldelta_arr  = ldelta_c
      else if (ldelta) then
        ldelta_arr  = .true.
      end if
!
! Load boolean array for kappa
!
      if (any(lkappa_c)) then
        lkappa = .true.
        do i=1,3
          lkappa_arr(i,1,1) = lkappa_c(i,1)
          lkappa_arr(i,2,1) = lkappa_c(i,2)
          lkappa_arr(i,1,2) = lkappa_c(i,2)
          lkappa_arr(i,3,1) = lkappa_c(i,3)
          lkappa_arr(i,1,3) = lkappa_c(i,3)
          lkappa_arr(i,2,2) = lkappa_c(i,4)
          lkappa_arr(i,2,3) = lkappa_c(i,5)
          lkappa_arr(i,3,2) = lkappa_c(i,5)
          lkappa_arr(i,3,3) = lkappa_c(i,6)
        enddo
      elseif (lkappa) then
        lkappa_arr = .true.
      end if
!
! Load boolean array for acoef
!
      if (any(lacoef_c)) then
        lacoef=.true.
        lacoef_arr(1,1) = lacoef_c(1)
        lacoef_arr(2,1) = lacoef_c(2)
        lacoef_arr(1,2) = lacoef_c(2)
        lacoef_arr(3,1) = lacoef_c(3)
        lacoef_arr(1,3) = lacoef_c(3)
        lacoef_arr(2,2) = lacoef_c(4)
        lacoef_arr(2,3) = lacoef_c(5)
        lacoef_arr(3,2) = lacoef_c(5)
        lacoef_arr(3,3) = lacoef_c(6)
        if (any([lalpha,lgamma])) then
          if (lroot) call warning('initialize_special', &
            'any lacoef_c=T overrides settings of lalpha and lgamma')
          lalpha=.false.; lgamma=.false.
          lalpha_arr = .false.; lgamma_arr = .false.
        endif
      elseif (lacoef) then
        lacoef_arr = .true.
      end if
!
! Load boolean array for bcoef
!
      if (any(lbcoef_c)) then
        lbcoef = .true.
        do i=1,3
          lbcoef_arr(i,1,1) = lbcoef_c(i,1)
          lbcoef_arr(i,2,1) = lbcoef_c(i,2)
          lbcoef_arr(i,1,2) = lbcoef_c(i,2)
          lbcoef_arr(i,3,1) = lbcoef_c(i,3)
          lbcoef_arr(i,1,3) = lbcoef_c(i,3)
          lbcoef_arr(i,2,2) = lbcoef_c(i,4)
          lbcoef_arr(i,2,3) = lbcoef_c(i,5)
          lbcoef_arr(i,3,2) = lbcoef_c(i,5)
          lbcoef_arr(i,3,3) = lbcoef_c(i,6)
        enddo
        if (any([lbeta,ldelta,lkappa])) then
          if (lroot) call warning('initialize_special', &
            'any lbcoef_c=T overrides settings of lbeta,ldelta,lkappa')
          lbeta=.false.; lbeta_arr = .false.
          ldelta=.false.; ldelta_arr = .false.
          lkappa=.false.; lkappa_arr = .false.
        endif
      elseif (lbcoef) then
        lbcoef_arr = .true.
      end if
!
! Load boolean array for umean
!
      if (any(lumean_c)) then
        lumean = .true.
        lumean_arr  = lumean_c
      else if (lumean) then
        lumean_arr  = .true.
      end if
!
! Store scales
!
      tensor_scales(alpha_id) = alpha_scale
      tensor_scales(beta_id)  = beta_scale
      tensor_scales(gamma_id) = gamma_scale
      tensor_scales(delta_id) = delta_scale
      tensor_scales(kappa_id) = kappa_scale
      tensor_scales(umean_id) = umean_scale
      tensor_scales(acoef_id) = acoef_scale
      tensor_scales(bcoef_id) = bcoef_scale
!
! Store names
!
      if (len_trim(dataset_name) > 0) then
          if (lroot) call warning('initialize_special', &
            '"dataset_name" given, overwriting individual dataset names.')
        alpha_name  = dataset_name
        beta_name   = dataset_name
        gamma_name  = dataset_name
        delta_name  = dataset_name
        kappa_name  = dataset_name
        umean_name  = dataset_name
        acoef_name  = dataset_name
        bcoef_name  = dataset_name
      end if
      tensor_names(alpha_id)  = alpha_name
      tensor_names(beta_id)   = beta_name
      tensor_names(gamma_id)  = gamma_name
      tensor_names(delta_id)  = delta_name
      tensor_names(kappa_id)  = kappa_name
      tensor_names(umean_id)  = umean_name
      tensor_names(acoef_id)  = acoef_name
      tensor_names(bcoef_id)  = bcoef_name

      ! For backwards compatibility, if dataset has been set,
      ! overwrite dataset_type and give a warning.
      if (len_trim(dataset) > 0) then
          if (lroot) call warning('initialize_special', &
            'parameter "dataset" overwrites "dataset_type"')
          dataset_type = dataset
      endif

    end subroutine parseParameters
!*****************************************************************************
    logical function output_persistent_special()
!
!  Writes out the time of the next SNI
!
!  13-Dec-2011/Bourdin.KIS: reworked
!  14-jul-2015/fred: removed obsolete Remnant persistant variable from current
!  write and added new cluster variables. All now consistent with any io
!
      use IO, only: write_persist
!
!      if (lcollective_IO) call fatal_error ('output_persistent_interstellar', &
!          "The interstellar persistent variables can't be written
!          collectively!")
!
      output_persistent_special = .true.
!
      if (write_persist ('SPECIAL_ILOAD', id_record_SPECIAL_ILOAD, iload)) return
!
      output_persistent_special = .false.
!
    endfunction output_persistent_special
!*****************************************************************************
!********************************************************************
!********************************************************************
!************        DO NOT DELETE THE FOLLOWING        *************
!********************************************************************
!**  This is an automatically generated include file that creates  **
!**  copies dummy routines from nospecial.f90 for any Special      **
!**  routines not implemented in this file                         **
!**                                                                **
    include '../special_dummies.inc'
!***********************************************************************
endmodule Special