! $Id$
!
!  This module contains routines both for delta-correlated
!  and continuous forcing. The fcont pencil is only provided
!  for continuous forcing.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lforcing = .true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 0
!
! PENCILS PROVIDED fcont(3,n_forcing_cont_max)
!
!***************************************************************
module Forcing
!
  use Cdata
  use General, only: keep_compiler_quiet
  use Messages
  use Geometrical_types
!
  implicit none
!
  include 'forcing.h'
  include 'record_types.h'
!
  interface input_persistent_forcing
     module procedure input_persist_forcing_id
     module procedure input_persist_forcing
  endinterface
!
  real :: force=0.,force2=0., force_double=0., force1_scl=1., force2_scl=1.
  real :: relhel=1., height_ff=0., r_ff=0., r_ff_hel=0., rcyl_ff=0., rel_zcomp=1.
  real :: Bconst=1., Bslope=0.
  real :: fountain=1.,width_ff=.5,nexp_ff=1.,n_hel_sin_pow=0.
  real :: crosshel=0.
  real :: radius_ff=0., k1_ff=1., kx_ff=1., ky_ff=1., kz_ff=1., z_center=0.
  real :: slope_ff=0., work_ff=0., omega_ff=0., omega_double_ff=0., n_equator_ff=1.
  real :: ampl_double_ff=0.
  real :: tforce_stop=impossible,tforce_stop2=impossible
  real :: tforce_start=0.,tforce_start2=0.
  real :: wff_ampl=0.,  xff_ampl=0.,  yff_ampl=0.,  zff_ampl=0.
  real :: wff2_ampl=0., xff2_ampl=0., yff2_ampl=0., zff2_ampl=0.
  real :: zff_hel=0.,max_force=impossible
  real :: dtforce=0., dtforce_ampl=.5, dtforce_duration=-1.0, force_strength=0.
  real :: b0_mode=0.
  real(KIND=rkind8) :: tforce_ramp_down=1.1, tauforce_ramp_down=1.
  real, dimension(3) :: force_direction=(/0.,0.,0./)
  real, dimension(3,2) :: location_fixed=0.
  real, dimension(nx) :: profx_ampl=1.,profx_hel=1., profx_ampl1=0.
  real, dimension(my) :: profy_ampl=1.,profy_hel=1.
  real, dimension(mz) :: profz_ampl=1.,profz_hel=1.,qdouble_profile=1.
  integer :: kfountain=5,iff,ifx,ify,ifz,ifff,iffx,iffy,iffz,i2fff,i2ffx,i2ffy,i2ffz,iff_aux
  integer :: kzlarge=1
  integer :: iforcing_zsym=0, nlocation=1
  logical :: lwork_ff=.false.,lmomentum_ff=.false.
  logical :: lhydro_forcing=.true., lneutral_forcing=.false., lmagnetic_forcing=.false.
  logical :: lcrosshel_forcing=.false.,ltestfield_forcing=.false.,ltestflow_forcing=.false.
  logical :: lxxcorr_forcing=.false., lxycorr_forcing=.false.
  logical :: lhelical_test=.false., lrandom_location=.true., lrandom_time=.false.
  logical :: lwrite_psi=.false., lforce_ramp_down=.false.
  logical :: lscale_kvector_tobox=.false.,lwrite_gausspot_to_file=.false.
  logical :: lwrite_gausspot_to_file_always=.false.
  logical :: lscale_kvector_fac=.false.
  logical :: lforce_peri=.false., lforce_cuty=.false.
  logical :: lforcing2_same=.false., lforcing2_curl=.false.
  logical :: lff_as_aux = .false.
  logical :: lforcing_osc = .false., lforcing_osc2 = .false., lforcing_osc_double = .false.
  real :: scale_kvectorx=1.,scale_kvectory=1.,scale_kvectorz=1.
  logical :: old_forcing_evector=.false., lforcing_coefs_hel_double=.false.
  character (len=labellen) :: iforce='zero', iforce2='zero'
  character (len=labellen) :: iforce_profile='nothing'
  character (len=labellen) :: iforce_tprofile='nothing'
  real :: equator=0.
  real :: kx_2df=0.,ky_2df=0.,xminf=0.,xmaxf=0.,yminf=0.,ymaxf=0.
! For helical forcing in spherical polar coordinate system
  real, allocatable, dimension(:,:,:) :: psif
  real, allocatable, dimension(:,:) :: cklist
  logical :: lfastCK=.false.,lsamesign=.true.
  logical :: lisotropize_SR=.false., SR_alpha_shift=.false., SR_ell_shift=.false., lfastSR=.false.
! allocated only if we have lfastCK=T
  real,allocatable,dimension(:,:,:) :: Zpsi_list
  real,allocatable,dimension(:,:,:) :: RYlm_list,IYlm_list
! allocated only if we have lfastSR=T
  real,allocatable,dimension(:,:,:,:) :: RYlm_rot,IYlm_rot
  integer :: nlist_ck=25, ell_shift=0, helsign=0., ncol=7, n_axisrot_angles=10
  real :: fpre = 1.0,ck_equator_gap=0.,ck_gap_step=0.,alpha_shift=5.
  integer :: icklist,jtest_aa0=5,jtest_uu0=1
! For random forcing
  logical :: lavoid_xymean=.false., lavoid_ymean=.false., lavoid_zmean=.false., &
             ldiscrete_phases=.false.
! For random forcing in 2d
  integer,allocatable, dimension (:,:) :: random2d_kmodes
  integer :: random2d_nmodes
  integer :: random2d_kmin=0, random2d_kmax=0
  integer :: channel_force=1
! continuous 2d forcing
  integer :: k2d
  logical :: l2dxz,l2dyz
! anisotropic forcing
  logical :: laniso_forcing_old=.true.
! Persistent stuff
  real :: tsforce=-10.
  real, dimension (3) :: location, location2
!
!  continuous forcing variables
!
  integer :: n_forcing_cont=n_forcing_cont_max
  logical :: lembed=.false.,lshearing_adjust_old=.false.
  logical, dimension(n_forcing_cont_max) :: lgentle=.false.
  character (len=labellen), dimension(n_forcing_cont_max) :: iforcing_cont='nothing'
  real, dimension(n_forcing_cont_max) :: ampl_ff=1., ampl1_ff=0., width_fcont=1., x1_fcont=0., x2_fcont=0.
  real, dimension(n_forcing_cont_max) :: kf_fcont=impossible, kf_fcont_x=impossible, kf_fcont_y=impossible, kf_fcont_z=impossible
  real, dimension(n_forcing_cont_max) :: omega_fcont=0., omegay_fcont=0., omegaz_fcont=0.
  real, dimension(n_forcing_cont_max) :: eps_fcont=0., tgentle=0., z_center_fcont=0.
  real, dimension(n_forcing_cont_max) :: ampl_bb=5.0e-2,width_bb=0.1,z_bb=0.1,eta_bb=1.0e-4
  real, dimension(n_forcing_cont_max) :: fcont_ampl=1., ABC_A=1., ABC_B=1., ABC_C=1., theta_TG=0.
  real :: ampl_diffrot=1.0,omega_exponent=1.0
  real :: omega_tidal=1.0, R0_tidal=1.0, phi_tidal=1.0, Omega_vortex=0.
  real :: cs0eff=impossible
  type(torus_rect), save :: torus    ! save needed due to internal initialization
  ! GP_TC13 forcing
  real :: tcor_GP=1.,kmin_GP=1.,kmax_GP=2.,beta_GP=1.3333
  integer :: nk_GP=2
!
!  auxiliary functions for continuous forcing function
!
  real, dimension (my,n_forcing_cont_max) :: phi1_ff
  real, dimension (mx,n_forcing_cont_max) :: phi2_ff
  real, dimension (mx,n_forcing_cont_max) :: sinx,cosx,sinxt,cosxt,embedx,expmk2x2
  real, dimension (my,n_forcing_cont_max) :: siny,cosy,sinyt,cosyt,embedy,expmk2y2
  real, dimension (mz,n_forcing_cont_max) :: sinz,cosz,sinzt,coszt,embedz
  real, dimension (100,n_forcing_cont_max) :: xi_GP,eta_GP
  real, allocatable, dimension (:,:,:,:) :: fcont_from_file
!
  namelist /forcing_run_pars/ &
       tforce_start,tforce_start2,&
       iforce,force,force_double,relhel,crosshel,height_ff,r_ff,r_ff_hel, &
       rcyl_ff,width_ff,nexp_ff,lff_as_aux,Bconst,Bslope, rel_zcomp, &
       iforce2, force2, force1_scl, force2_scl, iforcing_zsym, &
       kfountain, fountain, tforce_stop, tforce_stop2, &
       radius_ff,k1_ff,kx_ff,ky_ff,kz_ff,slope_ff,work_ff,lmomentum_ff, &
       omega_ff, omega_double_ff, n_equator_ff, ampl_double_ff, &
       location_fixed, lrandom_location, nlocation, &
       lwrite_gausspot_to_file,lwrite_gausspot_to_file_always, &
       wff_ampl, xff_ampl, yff_ampl, zff_ampl, zff_hel, &
       wff2_ampl, xff2_ampl,yff2_ampl, zff2_ampl, &
       lhydro_forcing, lneutral_forcing, lmagnetic_forcing, &
       ltestfield_forcing, ltestflow_forcing, &
       lcrosshel_forcing, lxxcorr_forcing, lxycorr_forcing, &
       jtest_aa0,jtest_uu0, &
       max_force,dtforce,dtforce_duration,old_forcing_evector, &
       lforcing_coefs_hel_double, dtforce_ampl, &
       iforce_profile, iforce_tprofile, lscale_kvector_tobox, &
       force_direction, force_strength, lhelical_test, &
       lfastCK,lisotropize_SR,fpre,helsign,nlist_ck,lwrite_psi,&
       ck_equator_gap,ck_gap_step,&
       ABC_A, ABC_B, ABC_C, &
       lforcing_cont,iforcing_cont, z_center_fcont, z_center, &
       lembed, k1_ff, ampl_ff, ampl1_ff, width_fcont, x1_fcont, x2_fcont, &
       kf_fcont, omega_fcont, omegay_fcont, omegaz_fcont, eps_fcont, &
       lsamesign, lshearing_adjust_old, equator, &
       lscale_kvector_fac,scale_kvectorx,scale_kvectory,scale_kvectorz, &
       lforce_peri,lforce_cuty,lforcing2_same,lforcing2_curl, &
       tgentle,random2d_kmin,random2d_kmax,l2dxz,l2dyz,k2d, &
       z_bb,width_bb,eta_bb,fcont_ampl, &
       ampl_diffrot,omega_exponent,kx_2df,ky_2df,xminf,xmaxf,yminf,ymaxf, &
       lavoid_xymean,lavoid_ymean,lavoid_zmean, ldiscrete_phases, &
       omega_tidal, R0_tidal, phi_tidal, omega_vortex, &
       lforce_ramp_down, tforce_ramp_down, tauforce_ramp_down, &
       n_hel_sin_pow, kzlarge, cs0eff, channel_force, torus, Omega_vortex, &
       lrandom_time,laniso_forcing_old, &
       lforcing_osc, lforcing_osc2, lforcing_osc_double, &
       tcor_GP, kmin_GP,kmax_GP,nk_GP,beta_GP, n_axisrot_angles, &
       SR_alpha_shift, alpha_shift, ell_shift, SR_ell_shift, ncol, lfastSR
!
! other variables (needs to be consistent with reset list below)
!
  integer :: idiag_rufm=0       ! DIAG_DOC: $\left<\rho\fv\cdot\uv\right>$
  integer :: idiag_rufint=0
  integer :: idiag_ufm=0        ! DIAG_DOC: $\left<\fv\cdot\uv\right>$
  integer :: idiag_ofm=0        ! DIAG_DOC: $\left<\ov\cdot\fv\right>$
  integer :: idiag_qfm=0
  integer :: idiag_ffm=0
  integer :: idiag_ruxfxm=0
  integer :: idiag_ruyfym=0
  integer :: idiag_ruzfzm=0
  integer :: idiag_ruxfym=0
  integer :: idiag_ruyfxm=0
  integer :: idiag_fxbxm=0
  integer :: idiag_fxbym=0
  integer :: idiag_fxbzm=0
!
! Auxiliaries
!
  real :: cs0
  real, dimension(:,:), pointer :: reference_state
  real, dimension(3) :: k1xyz=0.
  real, dimension(mz) :: profz_k
  logical :: lmhd_forcing
!
! Data from k.dat.
!
  integer :: nk, nk2
  real :: kav, kav2
  real, dimension(:), allocatable :: kkx,kky,kkz
  real, dimension(:), allocatable :: kkx2,kky2,kkz2
!
  real, allocatable, dimension (:,:) :: KS_k,KS_A,KS_B !or through whole field for each wavenumber?
  real, allocatable, dimension (:) :: KS_omega !or through whole field for each wavenumber?
  integer :: KS_modes = 25
!
  contains
!
!***********************************************************************
    subroutine register_forcing
!
!  add forcing in timestep
!  11-may-2002/wolf: coded

      use FarrayManager, only: farray_register_auxiliary
!
!  identify version number
!
      if (lroot) call svn_id( &
           "$Id$")

      if (iforce=='spherical_radial') then
        if (lisotropize_SR.and.lfastSR) call farray_register_auxiliary('force',iff_aux)
      endif
!
    endsubroutine register_forcing
!***********************************************************************
    subroutine initialize_forcing
!
!  read seed field parameters
!  nothing done from start.f90
!
!  25-sep-2014/MR: determine n_forcing_cont according to the actual selection
!  23-jan-2015/MR: reference state now fetched here and stored in module variable
!  15-feb-2015/MR: returning before entering continuous forcing section when
!                  no such forcing is requested
!  18-dec-2015/MR: minimal wavevectors k1xyz moved here from grid
!  14-Jun-2016/MR+NS: added forcing sinx*exp(-z^2)
!  11-May-2017/NS: added forcing Aycont_z
!  08-Aug-2019/MR: moved reading of k.dat, k_double.dat from individual subroutines
!                  -> nk, kav, kk[xyz] and nk2, kav2, kk2[xyz], respectively, now module 
!                  variables.
!  16-Aug-2023/MR+ABa: Added a non-helical spherical forcing acting only in the 
!                      radial direction that consists of Bessel functions in the 
!                      radial and randomized spherical Harmonics around the x-axis
!                      in the horizontal direction.  
!
      use General, only: bessj,itoa
      use SharedVariables, only: get_shared_variable
      use Sub, only: step,erfunc,stepdown,register_report_aux
      use EquationOfState, only: cs20
!
      real :: zstar,rmin,rmax,a_ell,anum,adenom,jlm_ff,ylm_ff,alphar,Balpha,RYlm,IYlm,intv_rotang
      real :: ang_intv,sthphase,cthphase,costhprime,phprime

      integer :: l,m,n,i,ilread,ilm,ckno,ilist,emm,aindex,Legendrel,iangle
      logical :: lk_dot_dat_exists

      cs0=sqrt(cs20)
!
      if (lstart) then
        if (ip<4) print*,'initialize_forcing: not needed in start'
      else
!
!  Check whether we want ordinary hydro forcing or magnetic forcing
!  Currently there is also the option of forcing just the neutrals,
!  but in future it could be combined with other forcings instead.
!
        if (lcrosshel_forcing) then
          if (.not.(lmagnetic.and.lhydro)) &
            call fatal_error('initialize_forcing','cross helicity forcing requires an MHD run')
          lmagnetic_forcing=.true.
          lhydro_forcing=.true.
          lmhd_forcing=.false.
        else
          lmhd_forcing=lmagnetic_forcing.and.lhydro_forcing .or. &
                       ltestfield_forcing.and.ltestflow_forcing
        endif 
        if (lmagnetic_forcing) then
          ifff=iaa; iffx=iax; iffy=iay; iffz=iaz
        endif 
        if (lneutral_forcing) then
          if (iuun==0) call fatal_error("lneutral_forcing","uun==0")
          ifff=iuun; iffx=iunx; iffy=iuny; iffz=iunz
        endif 
        if (lhydro_forcing) then
          if (lmagnetic_forcing) then
            i2fff=iuu; i2ffx=iux; i2ffy=iuy; i2ffz=iuz
          else
            if (iphiuu==0) then
              ifff=iuu; iffx=iux; iffy=iuy; iffz=iuz
            else
              ifff=iphiuu; iffx=iphiuu; iffy=0; iffz=0
            endif
          endif
        endif
        if (.not.(lmagnetic_forcing.or.lhydro_forcing.or.lneutral_forcing.or. &
                  ltestfield_forcing.or.ltestflow_forcing)) &
          call fatal_error("initialize_forcing","No forcing function set")
        
        if (ldebug) print*,'initialize_forcing: ifff=',ifff
!
!  check whether we want constant forcing at each timestep,
!  in which case lwork_ff is set to true.
!
        if (work_ff/=0.) then
          force=1.
          lwork_ff=.true.
          if (lroot) print*,'initialize_forcing: reset force=1., because work_ff is set'
        endif
      endif
!
!  initialize location to location_fixed
!
      location=location_fixed(:,1)
      location2=location_fixed(:,2)
!
!  vertical profiles for amplitude and helicity of the forcing
!  default is constant profiles for rms velocity and helicity.
!
      if (iforce_profile=='nothing') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.; profz_hel=1.
!
!  sine profile for amplitude of forcing 
!
      elseif (iforce_profile=='sinz') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_hel=1.
        do n=1,mz
          profz_ampl(n)=sin(kzlarge*z(n))
        enddo
!
      elseif (iforce_profile=='equator') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.
        do n=1,mz
          profz_hel(n)=sin(2*pi*z(n)*n_equator_ff/Lz)
        enddo
!
      elseif (iforce_profile=='equator_step') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.
        do n=1,mz
          profz_hel(n)= -1.+2.*step(z(n),equator-ck_equator_gap,ck_gap_step)
        enddo
!
!  step function change in intensity of helicity at zff_hel
!
      elseif (iforce_profile=='step_ampl=z') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
                       profz_hel=1.
        do n=1,mz
          profz_ampl(n)=step(z(n),zff_hel,width_ff)
        enddo
!
!  sign change of helicity proportional to cosy
!
      elseif (iforce_profile=='equator_hel=cosy') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.
        do m=1,my
          profy_hel(m)=cos(y(m))
        enddo
        profz_ampl=1.; profz_hel=1.
!
! step function profile
!
      elseif (iforce_profile=='equator_hel=step') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.
        do m=1,my
          profy_hel(m)= -1.+2.*step(y(m),yequator-ck_equator_gap,ck_gap_step)
        enddo
!
!  sign change of helicity about z=0
!
      elseif (iforce_profile=='equator_hel=z') then
        call fatal_error("initialize_forcing","use 'equator_hel=z/L' instead of 'equator_hel=z'")
!
!  Linear profile of helicity, normalized by ztop (or -zbot, if it is bigger)
!
      elseif (iforce_profile=='equator_hel=z/L') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.
        do n=1,mz
          profz_hel(n)=z(n)/max(-xyz0(3),xyz1(3))
        enddo
!
!  cosine profile of helicity about z=0
!
      elseif (iforce_profile=='cos(z)') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.
        do n=1,mz
          profz_hel(n)=cos(z(n))
        enddo
!
!  cosine profile of helicity about z=0
!
      elseif (iforce_profile=='cos(z/2)') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.
        do n=1,mz
          profz_hel(n)=cos(.5*z(n))
        enddo
!
!  Cosine profile of helicity about z=0, with max function.
!  Should be used only if |z| < 1.5*pi.
!
      elseif (iforce_profile=='cos(z/2)_with_halo') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.
        do n=1,mz
          profz_hel(n)=max(cos(.5*z(n)),0.)
        enddo
!
!  helicity profile proportional to z^2, but vanishing on the boundaries
!
      elseif (iforce_profile=='squared') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.
        do n=1,mz
          profz_hel(n)=(z(n)-xyz0(3))*(xyz1(3)-z(n))
        enddo
!
!  helicity profile proportional to z^3, but vanishing on the boundaries
!
      elseif (iforce_profile=='cubic') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.
        do n=1,mz
          profz_hel(n)=z(n)*(z(n)-xyz0(3))*(xyz1(3)-z(n))
        enddo
!
!  power law
!
      elseif (iforce_profile=='(z-z0)^n') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=(abs(z-xyz0(3))/height_ff)**nexp_ff
        profz_hel=1.
        if (lroot) print*,'profz_ampl=',profz_ampl
!
!  power law with offset
!
      elseif (iforce_profile=='1+(z-z0)^n') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        if (height_ff>0) then
          zstar=xyz0(3)
        elseif (height_ff<0) then
          zstar=xyz1(3)
        else
          zstar=impossible
          call fatal_error('initialize_forcing',"needs height_ff/=0 for iforce_profile='1+(z-z0)^n'")
        endif
        profz_ampl=(1.+(z-zstar)/height_ff)**nexp_ff
        profz_hel=1.
        if (lroot) print*,'profz_ampl=',profz_ampl
!
!  power law with offset
!
      elseif (iforce_profile=='1+(z-z0)^n_prof') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        if (height_ff>0) then
          zstar=xyz0(3)
        elseif (height_ff<0) then
          zstar=xyz1(3)
        else
          zstar=impossible
          call fatal_error('initialize_forcing',"needs height_ff/=0 for iforce_profile='1+(z-z0)^n_prof'")
        endif
        profz_ampl=((1.+(z-zstar)/height_ff)**nexp_ff)*.5*(1.+tanh(20.*cos(.55*z)))
        profz_hel=1.
        if (lroot) print*,'profz_ampl=',profz_ampl
!
!  exponential law
!
      elseif (iforce_profile=='exp(z/H)') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=exp(z/width_ff)
        profz_hel=1.
!
      elseif (iforce_profile=='exp(x/H)') then
        profz_ampl=1.; profz_hel=1.
        profy_ampl=1.; profy_hel=1.
        profx_ampl=exp(x(l1:l2)/width_ff)
        profx_hel=1.
!
      elseif (iforce_profile=='exp(-(z/H)^2)') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=exp(-(z/width_ff)**2)
        profz_hel=1.
!
      elseif (iforce_profile=='xybox') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=.5*(1.-erfunc(z/width_ff))
        profz_hel=1.
        profx_ampl= step(x(l1:l2),xminf,width_ff)-step(x(l1:l2),xmaxf,width_ff)
        do m=1,my
          profy_ampl(m)= step(y(m),yminf,width_ff)-step(y(m),ymaxf,width_ff)
        enddo
!
!  turn off forcing intensity above z=0
!
      elseif (iforce_profile=='surface_z') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=.5*(1.-erfunc((z-r_ff)/width_ff))
        profz_hel=1.
!
!  turn off helicity of forcing above x=r_ff
!
      elseif (iforce_profile=='surface_helx') then
        profx_ampl=1.; profx_hel=.5*(1.-erfunc((x(l1:l2)-r_ff)/width_ff))
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.; profz_hel=1.
!
!  turn off helicity of forcing above x=r_ff
!
      elseif (iforce_profile=='surface_helx_cosy') then
        profx_ampl=1.; profx_hel=.5*(1.-erfunc((x(l1:l2)-r_ff)/width_ff))
       profy_ampl=1.; profy_hel=cos(y)
       profz_ampl=1.; profz_hel=1.
!
!  turn off helicity of forcing above x=r_ff
!  used in Jabbari et al. (2015)
!
      elseif (iforce_profile=='surface_helx_cosy*siny**n_hel_sin_pow') then
        profx_ampl=1.; profx_hel=.5*(1.-erfunc((x(l1:l2)-r_ff)/width_ff))
        profy_ampl=1.; profy_hel=cos(y)*sin(y)**n_hel_sin_pow
        profz_ampl=1.; profz_hel=1.
!
!  turn off helicity of forcing above x=r_ff
!  but with step function in radius, as in Warnecke et al. (2011)
!
      elseif (iforce_profile=='surface_stepx_cosy*siny**n_hel_sin_pow') then
        profx_ampl=.5*(1.-erfunc((x(l1:l2)-r_ff)/width_ff))
        profx_hel=1.
        profy_ampl=1.; profy_hel=cos(y)*sin(y)**n_hel_sin_pow
        profz_ampl=1.; profz_hel=1.
!
!  turn off helicity of forcing above z=r_ff
!
      elseif (iforce_profile=='surface_helz') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.
        profz_hel=.5*(1.-erfunc((z-r_ff)/width_ff))
!
!  turn off helicity of forcing above z=0
!
      elseif (iforce_profile=='surface_amplhelz') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=.5*(1.-erfunc((z-r_ff)/width_ff))
        profz_hel=.5*(1.-erfunc((z-r_ff_hel)/width_ff))
!
!  turn off forcing intensity above z=z0, and
!  stepx profile of helicity
!
      elseif (iforce_profile=='surface_z_stepx') then
        profx_ampl=1.
        profy_ampl=1.; profy_hel=1.
        do l=l1,l2
          profx_hel(l-nghost)= -1.+2.*step(x(l),0.,width_ff)
        enddo
        profz_ampl=.5*(1.-erfunc((z-r_ff)/width_ff))
        profz_hel=1.
!
!  turn off forcing intensity above z=z0, and
!  stepy profile of helicity
!
      elseif (iforce_profile=='surface_z_stepy') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.
        do m=1,my
          profy_hel(m)= -1.+2.*step(y(m),0.,width_ff)
        enddo
        profz_ampl=.5*(1.-erfunc((z-r_ff)/width_ff))
        profz_hel=1.
!
! turn on forcing in a certain layer (generalization of 'forced_convection')
!
      elseif (iforce_profile=='xyzbox') then
        profx_ampl=.5*(1.+erfunc((x(l1:l2)-xff_ampl)/wff_ampl)) &
                  -.5*(1.+erfunc((x(l1:l2)-xff2_ampl)/wff2_ampl))
        profy_ampl=.5*(1.+erfunc((y-yff_ampl)/wff_ampl)) &
                  -.5*(1.+erfunc((y-yff2_ampl)/wff2_ampl))
        profz_ampl=.5*(1.+erfunc((z-zff_ampl)/wff_ampl)) &
                  -.5*(1.+erfunc((z-zff2_ampl)/wff2_ampl))
        profx_hel=1.
        profy_hel=1.
        profz_hel=1.
!
! turn on forcing in a certain layer (generalization of 'forced_convection')
!
      elseif (iforce_profile=='zlayer') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=.5*(1.+erfunc((z-zff_ampl)/wff_ampl)) &
                  -.5*(1.+erfunc((z-zff2_ampl)/wff2_ampl))
!        profz_hel=1.
        do n=1,mz
          profz_hel(n)=z(n)/max(-xyz0(3),xyz1(3))
        enddo
!
! turn on forcing in the bulk of the convection zone
!
      elseif (iforce_profile=='forced_convection') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=.5*(1.+ erfunc((z)/width_ff)) - 0.5*(1.+ erfunc((z-1.)/(2.*width_ff)))
        profz_hel=1.
!
!  cosx modulation
!
      elseif (iforce_profile=='cosx2') then
        print*,'--------------------'
        print*,'using iforce profile ', iforce_profile
        print*,'--------------------'
        profx_ampl=cos(kx_ff*x(l1:l2))**2
        profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.; profz_hel=1.
!
!  cosx modulation
!
      elseif (iforce_profile=='cosx^nexp') then
        profx_ampl=fountain*cos(kx_ff*x(l1:l2))**nexp_ff
        profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.; profz_hel=1.
!
!  turn off forcing intensity above x=x0, and
!  cosy profile of helicity
!
      elseif (iforce_profile=='surface_x_cosy') then
        profx_ampl=.5*(1.-erfunc((x(l1:l2)-r_ff)/width_ff))
        profx_hel=1.
        profy_ampl=1.
        do m=1,my
          profy_hel(m)=cos(y(m))
        enddo
        profz_ampl=1.; profz_hel=1.
!
!  turn off forcing intensity above x=x0, and
!  stepy profile of helicity, used in Warnecke et al. (2011)
!
      elseif (iforce_profile=='surface_x_stepy') then
        profx_ampl=.5*(1.-erfunc((x(l1:l2)-r_ff)/width_ff))
        profx_hel=1.
        profy_ampl=1.
        do m=1,my
          profy_hel(m)= -1.+2.*step(y(m),pi/2.,width_ff)
        enddo
        profz_ampl=1.; profz_hel=1.
!
!  turn off forcing intensity above x=x0
!
      elseif (iforce_profile=='surface_x') then
        profx_ampl=.5*(1.-erfunc((x(l1:l2)-r_ff)/width_ff))
        profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.; profz_hel=1.
!
!  turn off forcing intensity above x=r_ff and y=r_ff
!
      elseif (iforce_profile=='surface_xy') then
        profx_ampl=.5*(1.-erfunc((x(l1:l2)-r_ff)/width_ff)); profx_hel=1.
        profy_ampl=.5*(1.-erfunc((y-r_ff)/width_ff)); profy_hel=1.
        profz_ampl=1.; profz_hel=1.
!
!  turn off forcing intensity above x=r_ff and y=r_ff
!
      elseif (iforce_profile=='surface_xz') then
        profx_ampl=.5*(1.-erfunc((x(l1:l2)-r_ff)/width_ff)); profx_hel=1.
        profz_ampl=.5*(1.-erfunc((z-r_ff)/width_ff)); profz_hel=1.
        profy_ampl=1.; profy_hel=1.
!
!  turn on forcing intensity above x=r_ff
!
      elseif (iforce_profile=='above_x0') then
        profx_ampl=.5*(1.+erfunc((x(l1:l2)-r_ff)/width_ff))
        profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_ampl=1.; profz_hel=1.
!
!  turn on forcing intensity above x=x0 with cosy profile
!
      elseif (iforce_profile=='above_x0_cosy') then
        profx_ampl=.5*(1.+erfunc((x(l1:l2)-r_ff)/width_ff))
        profy_ampl=1.
        profx_hel=1.
        do m=1,my
          profy_hel(m)=cos(y(m));
        enddo
        profz_ampl=1.; profz_hel=1.
!
!  just a change in intensity in the z direction
!
      elseif (iforce_profile=='intensity') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        profz_hel=1.
        do n=1,mz
          profz_ampl(n)=.5+.5*cos(z(n))
        enddo
!
!  Galactic profile both for intensity and helicity
!
      elseif (iforce_profile=='galactic-old') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        do n=1,mz
          if (abs(z(n))<zff_ampl) profz_ampl(n)=.5*(1.-cos(z(n)))
          if (abs(z(n))<zff_hel ) profz_hel (n)=.5*(1.+cos(z(n)/2.))
        enddo
!
!  Galactic profile both for intensity and helicity
!  Normally one would put zff_ampl=zff_hel=pi
!
      elseif (iforce_profile=='galactic') then
        profx_ampl=1.; profx_hel=1.
        profy_ampl=1.; profy_hel=1.
        do n=1,mz
          profz_ampl(n)=.0
          profz_hel (n)=.0
          if (abs(z(n))<zff_ampl) profz_ampl(n)=.5*(1.+cos(z(n)))
          if (abs(z(n))<zff_hel ) profz_hel (n)=sin(z(n))
        enddo
!
! Galactic helicity profile for helicity
      elseif (iforce_profile=='gal-wind') then
         profx_ampl=1.; profx_hel=1.
         profy_ampl=1.; profy_hel=1.
         do n=1,mz
            profz_hel(n)=sin(pi*z(n)/Lz)
         enddo
!
      elseif (iforce_profile=='gal-wind2') then
         profx_ampl=1.; profx_hel=1.
         profy_ampl=1.; profy_hel=1.
         do n=1,mz
            profz_hel(n)=sin(2.*pi*z(n)/Lz)
         enddo
!
!  corona
!
      elseif (iforce_profile=='diffrot_corona') then
        profz_ampl=1.; profz_hel=1.
        profy_ampl=1.; profy_hel=1.
        profx_ampl=.5*(1.-tanh((x(l1:l2)-xff_ampl)/wff_ampl))
        profx_hel=1.
!
      else
        call fatal_error('initialize_forcing','no such iforce_profile: '//trim(iforce_profile))
      endif
!
!  Turn on forcing intensity for force_double above z=0.
!  Adding a z-profile is the default. To turn it off, put width_ff=0.
!  If it is turned off, we put qdouble_profile=.5, to give both 50%,
!  if force_double=force.
!
      if (lforcing_coefs_hel_double) then
        if (width_ff==0.) then
          qdouble_profile=.5
        else
          qdouble_profile=.5*(1.+erfunc((z-r_ff)/width_ff))
        endif
!
!  Debug output
!
        if (ip<17.and.lroot) then
          print*,'r_ff,width_ff=',r_ff,width_ff
          print*,'qdouble_profile=',qdouble_profile
        endif
      endif
!
!  at the first step, the sin and cos functions are calculated for all
!  x,y,z points and are then saved and used for all subsequent steps
!  and pencils
!
      if (ip<=6) print*,'forcing_cont:','lforcing_cont=',lforcing_cont,iforcing_cont

      if (lstart) return

      if (iforce=='2'           .or.iforce=='helical'     .or.iforce=='helical_kprof'.or. &
          iforce=='helical_both'.or.iforce=='irrotational'.or.iforce=='hel_smooth'   .or. &
          iforce=='noshear') then

        inquire(FILE="k.dat", EXIST=lk_dot_dat_exists)

        if (lk_dot_dat_exists) then
          if (lroot.and.ip<14) print*,'initialize_forcing: opening k.dat'
          open(9,file='k.dat',status='old')
          read(9,*) nk,kav
          if (lroot.and.ip<14) print*,'initialize_forcing: average k=',kav
          if (allocated(kkx)) deallocate(kkx,kky,kkz)
          allocate(kkx(nk),kky(nk),kkz(nk))
          read(9,*) kkx
          read(9,*) kky
          read(9,*) kkz
          close(9)
        else
          call inevitably_fatal_error('initialize_forcing', 'you must give an input k.dat file')
        endif

        if (lforcing_coefs_hel_double) then
!
!  Read k_double.dat once at first call.
!
          inquire(FILE="k_double.dat", EXIST=lk_dot_dat_exists)
          if (lk_dot_dat_exists) then
            if (lroot.and.ip<14) print*,'initialize_forcing: opening k_double.dat'
            open(9,file='k_double.dat',status='old')
            read(9,*) nk2,kav2
            if (lroot.and.ip<14) print*,'initialize_forcing: average k(2)=',kav2
            if (allocated(kkx2)) deallocate(kkx2,kky2,kkz2)
            allocate(kkx2(nk2),kky2(nk2),kkz2(nk2))
            read(9,*) kkx2
            read(9,*) kky2
            read(9,*) kkz2
            close(9)
          else
            call inevitably_fatal_error ('initialize_forcing','you must provide a k_double.dat file')
          endif
        endif
!
      elseif (iforce=='spherical_radial') then  
!
        if (.not. lspherical_coords) call fatal_error('initialize_forcing', &
                                     'spherical_radial forcing works only in spherical coordinates!')
!
! Read the list of values for emm, ell and alpha from file "alpha_in.dat". 
! The number of given alpha for each lm pair should be ncol-2. 
!
        open(unit=76,file="alpha_in.dat",status="old")
        read(76,*) ckno,rmin,rmax
        if (ckno/=nlist_ck) then
          call warning("initialize_forcing", &
          "Number of entries in alpha_in.dat "//trim(itoa(ckno))//" unequal nlist_ck="//trim(itoa(nlist_ck)))
          nlist_ck=ckno
        endif
        if (.not.allocated(cklist)) allocate(cklist(nlist_ck,ncol))
        do ilread=1,nlist_ck
          read(76,*) (cklist(ilread,ilm),ilm=1,ncol)
        enddo
        close(76)
!
        if (.not.allocated(psif)) allocate(psif(mx,my,mz))
!
        if (lfastSR.and.lisotropize_SR) then
          if (.not.allocated(RYlm_rot)) then
            allocate(RYlm_rot(my,mz,nlist_ck,n_axisrot_angles))
            allocate(IYlm_rot(my,mz,nlist_ck,n_axisrot_angles))
          endif
          ang_intv=pi/n_axisrot_angles
          do ilist=1,nlist_ck
            emm = cklist(ilist,1)
            Legendrel = cklist(ilist,2)
            do iangle=1,n_axisrot_angles
              do n=1,mz
                do m=1,my
                 
                  sthphase=sin(ang_intv*iangle); cthphase=cos(ang_intv*iangle)
                  costhprime=sthphase*cosph(n)*sinth(m)+cthphase*costh(m)
                  if (abs(costhprime)==1.) then
                    phprime=0.
                  else
                    phprime=acos((cthphase*cosph(n)*sinth(m)-sthphase*costh(m))/sqrt(1.-costhprime**2))
                  endif
                
                  call sp_harm_real_costh(RYlm,Legendrel,emm,costhprime,phprime)
                  call sp_harm_imag_costh(IYlm,Legendrel,emm,costhprime,phprime)
                  RYlm_rot(m,n,ilist,iangle)=RYlm
                  IYlm_rot(m,n,ilist,iangle)=IYlm
                enddo
              enddo
            enddo
          enddo
        endif 
 
      elseif (iforce=='chandra_kendall'.or.iforce=='cktest') then 
!
        if (.not. lspherical_coords) call fatal_error('initialize_forcing', &
                        'Chandrasekhar-Kendall forcing works only in spherical coordinates!')
!
        if (abs(helsign)/=1) call fatal_error("initialize_forcing", "CK forcing: helsign must be ±1")
!
        if (iforce=='cktest') then
          lhelical_test=.true.
          icklist=0
        endif

        if (lroot) print*,'Chandrasekhar-Kendall forcing in spherical polar coordinates'
        if (.not.allocated(psif)) allocate(psif(mx,my,mz),cklist(nlist_ck,5))
!
! Read the list of values for emm, ell and alpha from file "alpha_in.dat". This code is designed for 25 such.
!
        open(unit=76,file="alpha_in.dat",status="old")
        read(76,*) ckno,rmin,rmax
        if (ckno/=nlist_ck) call fatal_error("initialize_forcing", &
          "Number of entries in alpha_in.dat "//trim(itoa(ckno))//" unequal nlist_ck="//trim(itoa(nlist_ck)))

        do ilread=1,nlist_ck
          read(76,*) (cklist(ilread,ilm),ilm=1,5)
        enddo
        close(76)

        if (ck_equator_gap/=0) &
          profy_ampl=1.-step(y,pi/2.-ck_equator_gap,ck_gap_step)+step(y,pi/2.+ck_equator_gap,ck_gap_step)

        if (lfastCK) then

          if (.not.allocated(Zpsi_list)) then
            allocate(Zpsi_list(mx,nlist_ck,3))
            allocate(RYlm_list(my,mz,nlist_ck),IYlm_list(my,mz,nlist_ck))
          endif

          do ilist=1,nlist_ck
            emm = cklist(ilist,1)
            Legendrel = cklist(ilist,2)
            do n=1,mz
              do m=1,my
                call sp_harm_real(RYlm,Legendrel,emm,y(m),z(n))
                call sp_harm_imag(IYlm,Legendrel,emm,y(m),z(n))
                RYlm_list(m,n,ilist)=RYlm
                IYlm_list(m,n,ilist)=IYlm
              enddo
            enddo

            do aindex=1,3
              Balpha = cklist(ilist,2+aindex)
              call sp_bessely_l(anum,Legendrel,Balpha*xyz0(1))
              call sp_besselj_l(adenom,Legendrel,Balpha*xyz0(1))
              a_ell = -anum/adenom
              do l=1,mx
                alphar=Balpha*x(l)
                call sp_besselj_l(jlm_ff,Legendrel,alphar)
                call sp_bessely_l(ylm_ff,Legendrel,alphar)
                Zpsi_list(l,ilist,aindex) = (a_ell*jlm_ff+ylm_ff)
              enddo
            enddo
          enddo

        endif 
      elseif (iforce=='twist') then
        if (r_ff==0.) call fatal_error('initialize_forcing',"for iforce='twist', r_ff=0 impossible")
      endif

      if (lff_as_aux) call register_report_aux('ff', iff, ifx, ify, ifz)
!
!  Get reference_state. Requires that density is initialized before forcing.
!
      if (lreference_state) &
        call get_shared_variable('reference_state',reference_state,caller='initialize_forcing')
!
!  Make sure that testfields/testflows are registered.  
!    
      ltestfield_forcing = ltestfield_forcing.and.iaatest>0
      ltestflow_forcing = ltestflow_forcing.and.iuutest>0
!
!  At the moment, cs0 is used for normalization.
!  It is saved in param2.nml.
!
      if (cs0eff==impossible) then
        if (cs0==impossible) then
          cs0eff=1.
          if (lroot) print*,'initialize_forcing: for normalization, use cs0eff=',cs0eff
        else
          cs0eff=cs0
        endif
      endif

      if (.not.lforcing_cont) return
      
      do i=1,n_forcing_cont_max
        if ( iforcing_cont(i)=='nothing' ) then
          n_forcing_cont=i-1
          exit
        else
          if (kf_fcont(i)  ==impossible) kf_fcont(i)  =k1_ff
          if (kf_fcont_x(i)==impossible) kf_fcont_x(i)=kx_ff
          if (kf_fcont_y(i)==impossible) kf_fcont_y(i)=ky_ff
          if (kf_fcont_z(i)==impossible) kf_fcont_z(i)=kz_ff
          if (z_center_fcont(i)==0.) z_center_fcont(i)=z_center
        endif
!
!  all 3 sin and cos functions needed
!
        if (iforcing_cont(i)=='ABC' .or. iforcing_cont(i)=='Straining') then
          if (lroot) print*,'forcing_cont: ABC--calc sinx, cosx, etc'
          sinx(:,i)=sin(kf_fcont(i)*x); cosx(:,i)=cos(kf_fcont(i)*x)
          siny(:,i)=sin(kf_fcont(i)*y); cosy(:,i)=cos(kf_fcont(i)*y)
          sinz(:,i)=sin(kf_fcont(i)*z); cosz(:,i)=cos(kf_fcont(i)*z)
        elseif (iforcing_cont(i)=='RobertsFlow' .or. &
                iforcing_cont(i)=='RobertsFlowII' .or. &
                iforcing_cont(i)=='RobertsFlowMask' .or. &
                iforcing_cont(i)=='RobertsFlow_exact' ) then
          if (lroot) print*,'forcing_cont: Roberts Flow'
          sinx(:,i)=sin(kf_fcont(i)*x); cosx(:,i)=cos(kf_fcont(i)*x)
          siny(:,i)=sin(kf_fcont(i)*y); cosy(:,i)=cos(kf_fcont(i)*y)
          if (iforcing_cont(i)=='RobertsFlowMask') then
            profx_ampl=exp(-.5*x(l1:l2)**2/radius_ff**2)
            profy_ampl=exp(-.5*y**2/radius_ff**2)
          endif
        elseif (iforcing_cont(i)=='RobertsFlow-zdep' .or. &
                iforcing_cont(i)=='Roberts-for-SSD' .or. &
                iforcing_cont(i)=='elevator-flow') then
          if (lroot) print*,'forcing_cont: z-dependent Roberts Flow'
          sinx(:,i)=sin(kf_fcont(i)*x); cosx(:,i)=cos(kf_fcont(i)*x)
          siny(:,i)=sin(kf_fcont(i)*y); cosy(:,i)=cos(kf_fcont(i)*y)
        elseif (iforcing_cont(i)=='nocos') then
          if (lroot) print*,'forcing_cont: nocos flow'
          sinx(:,i)=sin(kf_fcont(i)*x)
          siny(:,i)=sin(kf_fcont(i)*y)
          sinz(:,i)=sin(kf_fcont(i)*z)
        elseif (iforcing_cont(i)=='Schur_helical') then
          if (lroot) print*,'forcing_cont: Schur helical flow'
          sinx(:,i)=sin(kf_fcont(i)*x); cosx(:,i)=cos(kf_fcont(i)*x)
          siny(:,i)=sin(kf_fcont(i)*y); cosy(:,i)=cos(kf_fcont(i)*y)
          sinz(:,i)=sin(kf_fcont(i)*z); cosz(:,i)=cos(kf_fcont(i)*z)
        elseif (iforcing_cont(i)=='Schur_nonhelical') then
          if (lroot) print*,'forcing_cont: Schur nonhelical flow'
          sinx(:,i)=sin(kf_fcont(i)*x); cosx(:,i)=cos(kf_fcont(i)*x)
          siny(:,i)=sin(kf_fcont(i)*y); cosy(:,i)=cos(kf_fcont(i)*y)
          sinz(:,i)=sin(kf_fcont(i)*z); cosz(:,i)=cos(kf_fcont(i)*z)
        elseif (iforcing_cont(i)=='KolmogorovFlow-x') then
          if (lroot) print*,'forcing_cont: Kolmogorov flow'
          cosx(:,i)=cos(kf_fcont(i)*x)
        elseif (iforcing_cont(i)=='KolmogorovFlow-z') then
          if (lroot) print*,'forcing_cont: Kolmogorov flow z'
          cosz(:,i)=cos(kf_fcont(i)*z)
        elseif (iforcing_cont(i)=='TG') then
          if (lroot) print*,'forcing_cont: TG'
          sinx(:,i)=sin(kf_fcont(i)*x); cosx(:,i)=cos(kf_fcont(i)*x)
          siny(:,i)=sin(kf_fcont(i)*y); cosy(:,i)=cos(kf_fcont(i)*y)
          cosz(:,i)=cos(kf_fcont(i)*z)
          if (lembed) then
            embedx(:,i)=.5+.5*tanh(x/width_fcont(i))*tanh((pi-x)/width_fcont(i))
            embedy(:,i)=.5+.5*tanh(y/width_fcont(i))*tanh((pi-y)/width_fcont(i))
            embedz(:,i)=.5+.5*tanh(z/width_fcont(i))*tanh((pi-z)/width_fcont(i))
            sinx(:,i)=embedx(:,i)*sinx(:,i); cosx(:,i)=embedx(:,i)*cosx(:,i)
            siny(:,i)=embedy(:,i)*siny(:,i); cosy(:,i)=embedy(:,i)*cosy(:,i)
            cosz(:,i)=embedz(:,i)*cosz(:,i)
          endif
        elseif (iforcing_cont(i)=='TG-random-nonhel' .or. &
                iforcing_cont(i)=='TG-random-hel') then
          if (lroot) print*,'forcing_cont: TG-random--calc sinx, cosx, etc'
          sinx(:,i)=sin(kf_fcont(i)*x); cosx(:,i)=cos(kf_fcont(i)*x)
          siny(:,i)=sin(kf_fcont(i)*y); cosy(:,i)=cos(kf_fcont(i)*y)
          sinz(:,i)=sin(kf_fcont(i)*z); cosz(:,i)=cos(kf_fcont(i)*z)
        elseif (iforcing_cont(i)=='cosx*cosy*cosz') then
          if (lroot) print*,'forcing_cont: cosx(:,i)*cosy(:,i)*cosz(:,i)'
          sinx(:,i)=sin(kf_fcont(i)*x); cosx(:,i)=cos(kf_fcont(i)*x)
          siny(:,i)=sin(kf_fcont(i)*y); cosy(:,i)=cos(kf_fcont(i)*y)
          sinz(:,i)=sin(kf_fcont(i)*z); cosz(:,i)=cos(kf_fcont(i)*z)
        elseif (iforcing_cont(i)=='sinx') then
          sinx(:,i)=sin(kf_fcont(i)*x)
          if (tgentle(i) > 0.) then
            lgentle(i)=.true.
            if (lroot) print *, 'initialize_forcing: gentle forcing till t = ', tgentle
          endif
        elseif (iforcing_cont(i)=='(0,sinx,0)'.or.iforcing_cont(i)=='(0,sinxsint,0)') then
          sinx(:,i)=sin(kf_fcont(i)*x)
        elseif (iforcing_cont(i)=='(0,0,cosx)') then
          cosx(:,i)=cos(kf_fcont(i)*x)
        elseif (iforcing_cont(i)=='(0,0,cosxcosy)') then
          cosx(:,i)=cos(kf_fcont(i)*x)
          cosy(:,i)=cos(kf_fcont(i)*y)
        elseif (iforcing_cont(i)=='B=(0,0,cosxcosy)') then
          sinx(:,i)=sin(kf_fcont(i)*x)
          siny(:,i)=sin(kf_fcont(i)*y)
          cosx(:,i)=cos(kf_fcont(i)*x)
          cosy(:,i)=cos(kf_fcont(i)*y)
        elseif (iforcing_cont(i)=='(sinz,cosz,0)') then
          sinz(:,i)=sin(kf_fcont(i)*z)
          cosz(:,i)=cos(kf_fcont(i)*z)
        elseif (iforcing_cont(i)=='(0,cosx*cosz,0)' .or. &
                iforcing_cont(i)=='(0,cosx*cosz,0)_Lor') then
          cosx(:,i)=cos(kf_fcont_x(i)*x)
          cosz(:,i)=cos(kf_fcont_z(i)*z)
          sinx(:,i)=sin(kf_fcont_x(i)*x)
          sinz(:,i)=sin(kf_fcont_z(i)*z)
          profy_ampl=exp(-(kf_fcont_y(i)*y)**2)
        elseif (iforcing_cont(i)=='(0,sinx*exp(-z^2),0)') then
          sinx(:,i)=sin(kf_fcont_x(i)*x)
          cosx(:,i)=cos(kf_fcont_x(i)*x)
        elseif (iforcing_cont(i)=='J0_k1x') then
          do l=l1,l2
            profx_ampl (l-l1+1)=ampl_ff(i) *bessj(0,k1bessel0*x(l))
            profx_ampl1(l-l1+1)=ampl1_ff(i)*bessj(1,k1bessel0*x(l))
          enddo
        elseif (iforcing_cont(i)=='gaussian-z') then
          profx_ampl=exp(-.5*x(l1:l2)**2/radius_ff**2)*ampl_ff(i)
          profy_ampl=exp(-.5*y**2/radius_ff**2)
          profz_ampl=exp(-.5*z**2/radius_ff**2)
        elseif (iforcing_cont(i)=='fluxring_cylindrical') then
          if (lroot) print*,'forcing_cont: fluxring cylindrical'
        elseif (iforcing_cont(i)=='vortex') then
          call torus_init(torus)
        elseif (iforcing_cont(i)=='KS') then
        !call periodic_KS_setup(-5./3.) !Kolmogorov spec. periodic KS
        !call random_isotropic_KS_setup(-5./3.,1.,(nxgrid)/2.) !old form
        call random_isotropic_KS_setup_test !Test KS model code with 3 specific
        !modes.
        elseif (iforcing_cont(i)=='exp(-x2-y2)') then
          expmk2x2(:,i)=exp(-kf_fcont(i)**2*x**2)
          expmk2y2(:,i)=exp(-kf_fcont(i)**2*y**2)
        elseif (iforcing_cont(i)=='xz') then
          if (lroot) print*,'forcing_cont: xz'
          profx_ampl=ampl_ff(i)*(x(l1:l2)-xyz0(1))*(xyz1(1)-x(l1:l2))
          profz_ampl=(z-xyz0(3))*(xyz1(3)-z)
        elseif (iforcing_cont(i)=='1-(4/3)(1-r^2/4)*r^2') then
          if (lroot) print*,'forcing_cont: 1-(4/3)(1-r^2/4)*r^2'
          profx_ampl=ampl_ff(i)*(1.-4./3.*(1.-.25*x(l1:l2)**2)*x(l1:l2)**2)
        elseif (iforcing_cont(i)=='tidal') then
          if (lroot) print*, 'forcing_cont: tidal'
          sinx(:,i)=sin(2.*pi*x/Lxyz(1))
          siny(:,i)=sin(2.*pi*y/Lxyz(2))
       elseif (iforcing_cont(i)=='from_file') then
          if (allocated(fcont_from_file)) deallocate(fcont_from_file)
          allocate(fcont_from_file(3,nxgrid,nygrid,nzgrid))
          
          ! To create forcing_cont.dat, see function pc.util.write_forcing_cont in the Python module.
          if (lroot.and.ip<14) print*,'initialize_forcing: opening forcing_cont.dat'
          open(1,file='forcing_cont.dat',status='old')
          read(1,*) fcont_from_file
          close(1)
        endif
      enddo
      if (n_forcing_cont==0) call warning('forcing','no valid continuous iforcing_cont specified')
!
!  minimal wavenumbers
!
      where( Lxyz/=0.) k1xyz=2.*pi/Lxyz
!
      if (r_ff /=0. .or. rcyl_ff/=0.) profz_k = tanh(z/width_ff)
!
    endsubroutine initialize_forcing
!***********************************************************************
    subroutine addforce(f)
!
!  Add forcing at the end of each time step.
!  Since forcing is constant during one time step,
!  this can be added as an Euler 1st order step
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      logical, save :: lfirstforce=.true., lfirstforce2=.true.
      logical, save :: llastforce=.true., llastforce2=.true.
!
!  Turn off forcing if t<tforce_start or t>tforce_stop.
!  This can be useful for producing good initial conditions
!  for turbulent decay experiments.
!
      call timing('addforce','entered')
!
      if ( (t>tforce_stop) .or. (t<tforce_start) ) then
        if ( (t>tforce_stop) .and. llastforce) then
          if (iforce/='zero' .and. lroot) print*, 'addforce: t>tforce_stop; no forcing'
          llastforce=.false.
        endif
      else
        if ( iforce/='zero' .and. lfirstforce .and. lroot ) lfirstforce=.false.
!--         print*, 'addforce: addforce started'
!
!  calculate and add forcing function
!
        select case (iforce)
        case ('2drandom_xy');     call forcing_2drandom_xy(f)
        case ('2drxy_simple');    call forcing_2drandom_xy_simple(f)
        case ('ABC');             call forcing_ABC(f)
        case ('MHD_mode');        call forcing_mhd_mode(f)
        case ('blobs');           call forcing_blobs(f)
        case ('blobHS_random');   call forcing_blobHS_random(f)
        case ('chandra_kendall','cktest'); call forcing_chandra_kendall(f)
        case ('spherical_radial');call forcing_spherical_radial(f)    
        case ('diffrot');         call forcing_diffrot(f,force)
        case ('fountain', '3');   call forcing_fountain(f)
        case ('hillrain');        call forcing_hillrain(f,force)
        case ('gaussianpot');     call forcing_gaussianpot(f,force)
        case ('white_noise');     call forcing_white_noise(f)
        case ('GP');              call forcing_GP(f)
        case ('Galloway-Proctor-92'); call forcing_GP92(f)
        case ('irrotational');    call forcing_irro(f,force)
        case ('helical', '2');    call forcing_hel(f)
        case ('helical_both');    call forcing_hel_both(f)
        case ('helical_kprof');   call forcing_hel_kprof(f)
        case ('hel_smooth');      call forcing_hel_smooth(f)
        case ('horiz-shear');     call forcing_hshear(f)
        case ('nocos');           call forcing_nocos(f)
        case ('TG');              call forcing_TG(f)
        case ('twist');           call forcing_twist(f)
        case ('tidal');           call forcing_tidal(f,force)
        case ('zero'); if (headt.and.ip<10) print*,'addforce: No forcing'
        case default;  call warning('addforce','no such iforce: '//trim(iforce))
        endselect
      endif
!
!  add *additional* forcing function
!
      if ( (t>tforce_stop2) .or. (t<tforce_start2) ) then
        if ( (t>tforce_stop2) .and. llastforce2 .and. lroot) &
            print*,'addforce: t>tforce_stop2; no forcing'
        if (t>tforce_stop2) llastforce2=.false.
      else
        if ( (iforce2/='zero') .and. lfirstforce2 .and. lroot) print*,'addforce: addforce2 started'
        lfirstforce2=.false.
!
        select case (iforce2)
        case ('diffrot');      call forcing_diffrot(f,force2)
        case ('fountain');     call forcing_fountain(f)
        case ('gaussianpot');  call forcing_gaussianpot(f,force2)
        case ('helical');      call forcing_hel(f)
        case ('helical_both'); call forcing_hel_both(f)
        case ('horiz-shear');  call forcing_hshear(f)
        case ('irrotational'); call forcing_irro(f,force2)
        case ('tidal');        call forcing_tidal(f,force2)
        case ('zero');
          if (headtt .and. lroot) print*,'addforce: No additional forcing'
        case default;
          if (lroot) print*,'addforce: No such forcing iforce2=',trim(iforce2)
        endselect
!
        if (headtt.or.ldebug) print*,'addforce: done addforce'
      endif
      call timing('addforce','finished')
!
    endsubroutine addforce
!***********************************************************************
    subroutine forcing_2drandom_xy(f)
!
!  Random force in two dimensions (x and y) limited to bands of wave-vector
!  in space.
!
!  14-feb-2011/ dhruba : coded
!
      use General, only: random_number_wrapper
      use EquationOfState, only: cs20
!
      real, dimension (mx,my,mz,mfarray) :: f
      logical, save :: lfirst_call=.true.
      integer :: ikmodes,iran1,iran2,kx1,ky1,kx2,ky2
      real, dimension(nx,3) :: forcing_rhs
      real, dimension(nx) :: xkx1,xkx2
      real,dimension(4) :: fran
      real :: phase1,phase2,pi_over_Lx,force_norm
!
      if (lfirst_call) then
! If this is the first time this routine is being called select a set of
! k-vectors in two dimensions according to inputparameters:
        if (lroot) print*,'forcing_2drandom_xy: selecting k vectors'
        call get_2dmodes (.true.)
        allocate(random2d_kmodes (2,random2d_nmodes))
        call get_2dmodes (.false.)
        if (lroot) then
! The root processors also write out the forced modes.
          open(unit=10,file=trim(datadir)//'/2drandomk.out',status='unknown')
          do ikmodes = 1, random2d_nmodes
            write(10,*) random2d_kmodes(1,ikmodes),random2d_kmodes(2,ikmodes)
          enddo
        endif
        lfirst_call=.false.
      endif
!
! force = xhat [ cos (k_1 x + \phi_1 ) + cos (k_2 y + \phi_1) ] +
!         yhat [ cos (k_3 x + \phi_2 ) + cos (k_4 y + \phi_2) ]
! where k_1 -- k_2 and k_3 -- k_4 are two randomly chose pairs in the list
! random2d_kmodes and  \phi_1 and \phi_2 are two random phases.
!
      call random_number_wrapper(fran,CHANNEL=channel_force)
      phase1=pi*(2*fran(1)-1.)
      phase2=pi*(2*fran(2)-1.)
      iran1=random2d_nmodes*(.9999*fran(3))+1
      iran2=random2d_nmodes*(.9999*fran(4))+1
!
!  normally we want to use the wavevectors as they are,
!  but in some cases, e.g. when the box is bigger than 2pi,
!  we want to rescale k so that k=1 now corresponds to a smaller value.
!
      if (lscale_kvector_fac) then
        kx1=random2d_kmodes(1,iran1)*scale_kvectorx
        ky1=random2d_kmodes(2,iran1)*scale_kvectory
        kx2=random2d_kmodes(1,iran2)*scale_kvectorx
        ky2=random2d_kmodes(2,iran2)*scale_kvectory
        pi_over_Lx=0.5
      elseif (lscale_kvector_tobox) then
        kx1=random2d_kmodes(1,iran1)*(2.*pi/Lxyz(1))
        ky1=random2d_kmodes(2,iran1)*(2.*pi/Lxyz(2))
        kx2=random2d_kmodes(1,iran2)*(2.*pi/Lxyz(1))
        ky2=random2d_kmodes(2,iran2)*(2.*pi/Lxyz(2))
        pi_over_Lx=pi/Lxyz(1)
      else
        kx1=random2d_kmodes(1,iran1)
        ky1=random2d_kmodes(2,iran1)
        kx2=random2d_kmodes(1,iran2)
        ky2=random2d_kmodes(2,iran2)
        pi_over_Lx=0.5
      endif
!
! Now add the forcing
!
      force_norm = force*cs20*sqrt(dt)
      do n=n1,n2
        do m=m1,m2
          xkx1 = x(l1:l2)*kx1+phase1
          xkx2 = x(l1:l2)*kx2+phase2
          forcing_rhs(:,1) = force_norm*( cos(xkx1) + cos(y(m)*ky1+phase1) )
          forcing_rhs(:,2) = force_norm*( cos(xkx2) + cos(y(m)*ky2+phase2) )
          forcing_rhs(:,3) = 0.
          if (lhelical_test) then
            f(l1:l2,m,n,iuu:iuu+2)=forcing_rhs(:,1:3)
          else
            f(l1:l2,m,n,iuu:iuu+2)=f(l1:l2,m,n,iuu:iuu+2)+forcing_rhs(:,1:3)
          endif
        enddo
      enddo
!
      if (ip<=9) print*,'forcing_2drandom_xy: forcing OK'
!
    endsubroutine forcing_2drandom_xy
!***********************************************************************
    subroutine get_2dmodes (lonly_total)
!
      integer :: ik1,ik2,modk
      integer :: imode=1
      logical :: lonly_total
!
!  15-feb-2011/dhruba: coded
!
      imode=1
      do ik1=0,random2d_kmax
        do ik2 = 0, random2d_kmax
          modk = nint(sqrt ( float(ik1*ik1) + float (ik2*ik2) ))
          if ( (modk<=random2d_kmax).and.(modk>=random2d_kmin)) then
            if (.not.lonly_total) then
              random2d_kmodes(1,imode) = ik1
              random2d_kmodes(2,imode) = ik2
            endif
            imode=imode+1
          endif
        enddo
      enddo
      if (lonly_total) random2d_nmodes = imode
!
    endsubroutine get_2dmodes
!***********************************************************************
    subroutine forcing_2drandom_xy_simple(f)
!
!  Random force in two dimensions (x and y) limited to bands of wave-vector
!  in space.
!
!  14-feb-2011/ dhruba : coded
!
      use General, only: random_number_wrapper, itoa
      use Sub, only: step
      use EquationOfState, only: cs20
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension(nx,3) :: forcing_rhs
      real, dimension(nx) :: xkx
      real,dimension(4) :: fran
      real :: phase1,phase2,force_norm
!
! force = xhat [ cos (k_2 y + \phi_1) ] +
!         yhat [ cos (k_3 x + \phi_2 ) ]
! where k_1 -- k_2 and k_3 -- k_4 are two randomly chose pairs in the list
! random2d_kmodes and  \phi_1 and \phi_2 are two random phases.
!
      call random_number_wrapper(fran,CHANNEL=channel_force)
      phase1=pi*(2*fran(1)-1.)
      phase2=pi*(2*fran(2)-1.)
!
!  normally we want to use the wavevectors as they are,
!  but in some cases, e.g. when the box is bigger than 2pi,
!  we want to rescale k so that k=1 now corresponds to a smaller value.
!
! Now add the forcing
!
      force_norm = force*cs20*sqrt(dt)
      do n=n1,n2
        do m=m1,m2
          xkx = x(l1:l2)*kx_2df+phase1
          forcing_rhs(:,1) = force_norm*(cos(y(m)*ky_2df+phase2) )*profx_ampl*profy_ampl(m)
          forcing_rhs(:,2) = force_norm*(sin(xkx) )*profx_ampl*profy_ampl(m)
          forcing_rhs(:,3) = 0.
          if (lhelical_test) then
            f(l1:l2,m,n,iuu:iuu+2)=forcing_rhs(:,1:3)
          else
            f(l1:l2,m,n,iuu:iuu+2)=f(l1:l2,m,n,iuu:iuu+2)+forcing_rhs(:,1:3)
          endif
        enddo
      enddo
!
      if (ip<=9) print*,'forcing_2drandom_xy_simple: forcing OK'
!
    endsubroutine forcing_2drandom_xy_simple
!***********************************************************************
    subroutine forcing_irro(f,force_ampl)
!
!  add acoustic forcing function, using a set of precomputed wavevectors
!  This forcing drives pressure waves
!
!  10-sep-01/axel: coded
!   6-feb-13/MR: discard wavevectors [0,0,kz] if lavoid_yxmean
!  06-dec-13/nishant: made kkx etc allocatable
!  20-aug-14/MR: discard wavevectors [kx,0,kz] if lavoid_ymean, [kx,ky,0] if lavoid_zmean
!  21-jan-15/MR: changes for use of reference state.
!
      use Diagnostics, only: sum_mn_name
      use General, only: random_number_wrapper
      use Sub, only: del2v_etc,dot,dot_mn
      use Mpicomm, only: mpireduce_sum
!
      real, dimension (mx,my,mz,mfarray) :: f
      real :: kx0,kx,ky,kz,force_ampl,pi_over_Lx
      real :: phase,ffnorm
      real, dimension (2) :: fran
      real, dimension (nx) :: rho1,qf
      real, dimension (nx,3) :: forcing_rhs,curlo
      complex, dimension (mx) :: fx
      complex, dimension (my) :: fy
      complex, dimension (mz) :: fz
      complex, dimension (3) :: ikk
      integer :: ik,j,jf
!
      do
        call random_number_wrapper(fran,CHANNEL=channel_force)
        phase=pi*(2*fran(1)-1.)
        ik=nk*(.9999*fran(2))+1
!
!  if lavoid_xymean=T and wavevector is close enough to [0,0,kz] discard it
!  and look for a new one
!
        if ( lavoid_xymean ) then
          if ( abs(kkx(ik))>.9*k1xyz(1) .or. abs(kky(ik))>.9*k1xyz(2) ) exit
        elseif ( lavoid_ymean ) then
          if ( abs(kky(ik))>.9*k1xyz(2) ) exit
        elseif ( lavoid_zmean ) then
          if ( abs(kkz(ik))>.9*k1xyz(3) ) exit
        else
          exit
        endif
      enddo
!
      if (ip<=6) print*,'forcing_irro: ik,phase,kk=',ik,phase,kkx(ik),kky(ik),kkz(ik),dt
!
!  normally we want to use the wavevectors as they are,
!  but in some cases, e.g. when the box is bigger than 2pi,
!  we want to rescale k so that k=1 now corresponds to a smaller value.
!
      if (lscale_kvector_fac) then
        kx0=kkx(ik)*scale_kvectorx
        ky=kky(ik)*scale_kvectory
        kz=kkz(ik)*scale_kvectorz
        pi_over_Lx=0.5
      elseif (lscale_kvector_tobox) then
        kx0=kkx(ik)*(2.*pi/Lxyz(1))
        ky=kky(ik)*(2.*pi/Lxyz(2))
        kz=kkz(ik)*(2.*pi/Lxyz(3))
        pi_over_Lx=pi/Lxyz(1)
      else
        kx0=kkx(ik)
        ky=kky(ik)
        kz=kkz(ik)
        pi_over_Lx=0.5
      endif
!
!  in the shearing sheet approximation, kx = kx0 - St*k_y.
!  Here, St=-deltay/Lx. However, to stay near kx0, we ignore
!  integer shifts.
!
      if (Sshear==0.) then
        kx=kx0
      else
        if (lshearing_adjust_old) then
          kx=kx0+ky*deltay/Lx
        else
          kx=kx0+mod(ky*deltay/Lx-pi_over_Lx,2.*pi_over_Lx)+pi_over_Lx
        endif
      endif
!
!  Need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference.
!  Divide also by kav, so it would be ffnorm=force*sqrt(kav/dt)*dt/kav,
!  but this can be simplified to give:
!
      ffnorm=force_ampl*sqrt(dt/kav)
!
!  pre-calculate for the contributions e^(ikx*x), e^(iky*y), e^(ikz*z),
!  as well as the phase factor.
!
      
      fx=exp(cmplx(0.,kx*x+phase))*ffnorm
      fy=exp(cmplx(0.,ky*y))
      fz=exp(cmplx(0.,kz*z))
!
!  write i*k as vector
!
      ikk(1)=cmplx(0.,kx)
      ikk(2)=cmplx(0.,ky)
      ikk(3)=cmplx(0.,kz)
!
!  Loop over all directions, but skip over directions with no extent.
! 
      lfirstpoint=.true.
      do n=n1,n2
      do m=m1,m2
!
!  Can force either velocity (default) or momentum (perhaps more physical).
!
        if (lmomentum_ff) then
          if (ldensity_nolog) then
            if (lreference_state) then
              rho1=1./(f(l1:l2,m,n,irho)+reference_state(:,iref_rho))
            else
              rho1=1./f(l1:l2,m,n,irho)
            endif
          else
            rho1=exp(-f(l1:l2,m,n,ilnrho))
          endif
        else
          rho1=1.
        endif
        do j=1,3
          if (lactive_dimension(j)) then
            jf=j+ifff-1
            forcing_rhs(:,j) = profx_ampl*profy_ampl(m)*profz_ampl(n) &
                              *rho1*real(ikk(j)*fx(l1:l2)*fy(m)*fz(n))
            f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)
          endif
        enddo
        if (lout) then
          if (idiag_qfm/=0) then
            call del2v_etc(f,iuu,curlcurl=curlo)
            call dot_mn(curlo,forcing_rhs,qf)
            call sum_mn_name(qf,idiag_qfm)
          endif
        endif
        lfirstpoint=.false.
      enddo
      enddo
!
!  For printouts, irufm needs to be communicated to other processors.
!
      if (lout) then
        if (idiag_qfm/=0) call mpireduce_sum(fname(idiag_qfm)/nwgrid,fname(idiag_qfm))
      endif
!
    endsubroutine forcing_irro
!***********************************************************************
    subroutine forcing_coefs_hel(coef1,coef2,coef3,fda,fx,fy,fz)
!
!  Calculates position-independent and 1D coefficients for helical forcing.
!
!  4-oct-17/MR: outsourced from forcing_hel.
!               Spotted bug: for old_forcing_evector=T, kk and ee remain undefined
!                            - needs to be fixed
      use Sub
!
      real,    dimension (3), intent(out) :: coef1,coef2,coef3,fda
      complex, dimension (mx),intent(out) :: fx
      complex, dimension (my),intent(out) :: fy
      complex, dimension (mz),intent(out) :: fz

      real :: phase, fact 
      real, dimension(3) :: kk
!
!  The routine fconst_coefs_hel generates various coefficients, including the phase.
!
      call fconst_coefs_hel(force,kkx,kky,kkz,nk,kav,coef1,coef2,coef3,kk,phase,fact,fda)
      call fxyz_coefs_hel(coef1,coef2,coef3,kk,phase,fact,fx,fy,fz)

    endsubroutine forcing_coefs_hel
!***********************************************************************
    subroutine forcing_pars_hel(coef1,coef2,coef3,fda,kk,phase,fact)
!
!  Calculates position-independent and 1D coefficients for helical forcing.
!  But this routine doesn't seem to be called from anywhere anymore.
!
!  4-oct-17/MR: outsourced from forcing_hel.
!               Spotted bug: for old_forcing_evector=T, kk and ee remain undefined
!                            - needs to be fixed
      use Sub
!
      real, dimension (3), intent(out) :: coef1,coef2,coef3,fda,kk
      real, intent(out) :: phase, fact 
!
      call fconst_coefs_hel(force,kkx,kky,kkz,nk,kav,coef1,coef2,coef3,kk,phase,fact,fda)

    endsubroutine forcing_pars_hel
!***********************************************************************
    subroutine forcing_coefs_hel2(coef1,coef2,coef3,fda,fx,fy,fz)
!
!  Modified copy of forcing_coefs_hel
!  Calculates position-independent and 1D coefficients for helical forcing.
!
!  10-nov-18/axel: adapted from forcing_coefs_hel to read also k_double.dat.
!
      use Sub
!
      real,    dimension (3), intent(out) :: coef1,coef2,coef3,fda
      complex, dimension (mx),intent(out) :: fx
      complex, dimension (my),intent(out) :: fy
      complex, dimension (mz),intent(out) :: fz
!
      real, dimension (3) :: kk
      real :: phase, fact 
!
!  This one also calculates coefficients, which are new ones, independent
!  of those in subroutine forcing_coefs_hel
!
      call fconst_coefs_hel(force_double,kkx2,kky2,kkz2,nk2,kav2,coef1,coef2,coef3,kk,phase,fact,fda)
      call fxyz_coefs_hel(coef1,coef2,coef3,kk,phase,fact,fx,fy,fz)
!
    endsubroutine forcing_coefs_hel2
!***********************************************************************
    subroutine fconst_coefs_hel(force_fact,kkx,kky,kkz,nk,kav,coef1,coef2,coef3,kk,phase,fact,fda)
!
!  This routine is can be called with any values of kkx,kky,kkz
!  to produce coef1,coef2,coef3,kk,phase,fact and fda.
!
!  08-aug-19/MR: modified to provide kk,phase,fact instead of fx,fy,fz
! 
      use General, only: random_number_wrapper,random_seed_wrapper
      use Sub
!      use Cdata, only: seed, nseed
!
      real,                   intent(in ) :: force_fact, kav
      integer,                intent(in ) :: nk
      real,    dimension (nk),intent(in ) :: kkx,kky,kkz
      real,    dimension (3), intent(out) :: coef1,coef2,coef3,kk,fda
      real,                   intent(out) :: phase,fact
!
      real :: ffnorm
      real, dimension (2) :: fran
      integer :: ik
      real :: kx0,kx,ky,kz,k2,k,pi_over_Lx
      real :: ex,ey,ez,kde,kex,key,kez,kkex,kkey,kkez
      real, dimension(3) :: e1,e2,ee
      real :: norm,phi
      real :: fd,fd2
      integer, save :: ideter=-1
!
!  generate random coefficients -1 < fran < 1
!  ff=force*Re(exp(i(kx+phase)))
!  |k_i| < akmax
!
      do
        if (ideter>=0) then      ! deterministic selection of wavevector and phase
          ik=mod(ideter-1,nk)+1
          phase=pi*(real(ik)/real(nk))
          ideter=ideter+1
!print*, 'itsub,ideter=', itsub,ideter, ik, kkx(ik),kky(ik),kkz(ik),phase
        else                     !random choice (default)
          call random_number_wrapper(fran,CHANNEL=channel_force)
!          call random_seed_wrapper(GET=seed)
!if (lroot) write(20,*) 'forcing: seed=', seed(1:nseed)
          phase=pi*(2*fran(1)-1.)
!AB: to add time-dependence XX
          ik=nk*(.9999*fran(2))+1
        endif
!
!  if lavoid_xymean=T and wavevector is close enough to [0,0,kz] discard it
!  and look for a new one
!
        if ( lavoid_xymean ) then
          if ( abs(kkx(ik))>.9*k1xyz(1) .or. abs(kky(ik))>.9*k1xyz(2) ) exit
        elseif ( lavoid_ymean ) then
          if ( abs(kky(ik))>.9*k1xyz(2) ) exit
        elseif ( lavoid_zmean ) then
          if ( abs(kkz(ik))>.9*k1xyz(3) ) exit
        else
          exit
        endif
      enddo
!
      if (ip<=6) then
        print*,'fconst_coefs_hel: ik,phase=',ik,phase
        print*,'fconst_coefs_hel: kx,ky,kz=',kkx(ik),kky(ik),kkz(ik)
      endif
!
!  normally we want to use the wavevectors as they are,
!  but in some cases, e.g. when the box is bigger than 2pi,
!  we want to rescale k so that k=1 now corresponds to a smaller value.
!
      if (lscale_kvector_fac) then
        kx0=kkx(ik)*scale_kvectorx
        ky=kky(ik)*scale_kvectory
        kz=kkz(ik)*scale_kvectorz
        pi_over_Lx=0.5
      elseif (lscale_kvector_tobox) then
        kx0=kkx(ik)*(2.*pi/Lxyz(1))
        ky=kky(ik)*(2.*pi/Lxyz(2))
        kz=kkz(ik)*(2.*pi/Lxyz(3))
        pi_over_Lx=pi/Lxyz(1)
      else
        kx0=kkx(ik)
        ky=kky(ik)
        kz=kkz(ik)
        pi_over_Lx=0.5
      endif
!
!  in the shearing sheet approximation, kx = kx0 - St*k_y.
!  Here, St=-deltay/Lx. However, to stay near kx0, we ignore
!  integer shifts.
!
      if (Sshear==0.) then
        kx=kx0
      else
        if (lshearing_adjust_old) then
          kx=kx0+ky*deltay/Lx
        else
          kx=kx0+mod(ky*deltay/Lx-pi_over_Lx,2.*pi_over_Lx)+pi_over_Lx
        endif
      endif
!
!  compute k^2 and output wavenumbers
!
      k2=kx**2+ky**2+kz**2
      k=sqrt(k2)
      if (ip<4) then
        open(89,file='forcing_hel_output.dat',position='append')
        write(89,'(6f10.5)') k,kx0,kx,ky,kz,deltay
        close(89)
      endif
!
!  Find e-vector:
!  Start with old method (not isotropic) for now.
!  Pick e1 if kk not parallel to ee1. ee2 else.
!
      if ((ky==0).and.(kz==0)) then
        ex=0; ey=1; ez=0
      else
        ex=1; ey=0; ez=0
      endif
!
      if (old_forcing_evector) then
!!! kk, ee not defined!
      else
!
!  Isotropize ee in the plane perp. to kk by
!  (1) constructing two basis vectors for the plane perpendicular
!      to kk, and
!  (2) choosing a random direction in that plane (angle phi)
!  Need to do this in order for the forcing to be isotropic.
!
        kk = (/kx, ky, kz/)
        ee = (/ex, ey, ez/)
        call cross(kk,ee,e1)
        call dot2(e1,norm); e1=e1/sqrt(norm) ! e1: unit vector perp. to kk
        call cross(kk,e1,e2)
        call dot2(e2,norm); e2=e2/sqrt(norm) ! e2: unit vector perp. to kk, e1
        call random_number_wrapper(phi,CHANNEL=channel_force); phi = phi*2*pi
        ee = cos(phi)*e1 + sin(phi)*e2
        ex=ee(1); ey=ee(2); ez=ee(3)
      endif
!
!  k.e
!
      call dot(kk,ee,kde)
!
!  k x e
!
      kex=ky*ez-kz*ey
      key=kz*ex-kx*ez
      kez=kx*ey-ky*ex
!
!  k x (k x e)
!
      kkex=ky*kez-kz*key
      kkey=kz*kex-kx*kez
      kkez=kx*key-ky*kex
!
!  ik x (k x e) + i*phase
!
!  Normalize ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!
!  This does already include the new sqrt(2) factor (missing in B01).
!  So, in order to reproduce the 0.1 factor mentioned in B01
!  we have to set force=0.07.
!
!  Furthermore, for |relhel| < 1, sqrt(2) should be replaced by
!  sqrt(1.+relhel**2). This is done now (9-nov-02).
!  This means that the previous value of force=0.07 (for relhel=0)
!  should now be replaced by 0.05.
!
!  Note: kav is not to be scaled with k1_ff (forcing should remain
!  unaffected when changing k1_ff).
!
      ffnorm = sqrt(1.+relhel**2)*k*sqrt(k2-kde**2)/sqrt(kav*cs0eff**3)*(k/kav)**slope_ff
      if (ip<=9) then
        print*,'fconst_coefs_hel: k,kde,kav=',k,kde,kav
        print*,'fconst_coefs_hel: cs0eff,ffnorm=',cs0eff,ffnorm
        print*,'fconst_coefs_hel: k*sqrt(k2-kde**2)=',k*sqrt(k2-kde**2)
      endif
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
      fact=force_fact/ffnorm*sqrt(dt)
!
!  prefactor; treat real and imaginary parts separately (coef1 and coef2),
!  so they can be multiplied by different profiles below.
!
      coef1=k*(/kex,key,kez/)
      coef2=relhel*(/kkex,kkey,kkez/)
!
!  possibly multiply forcing by sgn(z) and radial profile
!
      if (rcyl_ff==0.) coef3=crosshel*k*(/kkex,kkey,kkez/)

      !if (ip<=5) print*,'fconst_coefs_hel: coef=',coef1,coef2
!
! An attempt to implement anisotropic forcing using direction
! dependent forcing amplitude. Activated only if force_strength,
! describing the anisotropic part of the forcing, is
! nonzero. force_direction, which is a vector, defines the preferred
! direction of forcing. The expression for the forcing amplitude used
! at the moment is:
!
!  f(i)=f0*[1+epsilon(delta_ij*(k(i)*fd(j))/(|k||fd|))^2*fd(i)/|fd|]
!
! here f0 and fd are shorthand for force and forcing_direction,
! respectively, and epsilon=force_strength/force.
!
! 09-feb-22/hongzhe: possibility of a new implementation. If
!                    laniso_forcing_old=F then fda is the same in all
!                    three directions, and fda=sqrt(1+epsilon*mu^2)
!                    where mu = cos of the angle between kk and fd
!
      if (force_strength/=0.) then
        call dot(force_direction,force_direction,fd2)
        fd=sqrt(fd2)
        if (laniso_forcing_old) then
          fda = 1. + (force_strength/force)*(kk*force_direction/(k*fd))**2 * force_direction/fd
        else
          fda = sqrt(1. + (force_strength/force)*(dot_product(kk,force_direction)/(k*fd))**2)
        endif
      else
        fda = 1.
      endif
!
    endsubroutine fconst_coefs_hel
!**************************************************************************
    subroutine fxyz_coefs_hel(coef1,coef2,coef3,kk,phase,fact,fx,fy,fz)
!
!  08-aug-19/MR: carved out to produce fx,fy,fz from the other parameters.
! 
      use General, only: itoa

      real,    dimension (3), intent(in ) :: coef1,coef2,coef3,kk
      real,                   intent(in ) :: phase, fact
      complex, dimension (mx),intent(out) :: fx
      complex, dimension (my),intent(out) :: fy
      complex, dimension (mz),intent(out) :: fz

      real :: kx, kz, phase_x

      if (ldiscrete_phases) then
        kx=kk(1)*k1_ff
        phase_x=nint(phase/(kx*dx))*dx
        fx=cmplx(cos(kx*(x+phase_x)),sin(kx*(x+phase_x)))*fact
      else
        fx=exp(cmplx(0.,kk(1)*k1_ff*x+phase))*fact
      endif

      fy=exp(cmplx(0.,kk(2)*k1_ff*y))
!
!  symmetry of forcing function about z direction
!
      kz = kk(3)
      select case (iforcing_zsym)
        case(0); fz=exp(cmplx(0.,kz*k1_ff*z))
        case(1); fz=cos(kz*k1_ff*z)
        case(-1); fz=sin(kz*k1_ff*z)
        case default
          call fatal_error('fxyz_coefs_hel','no such iforcing_zsym: '//trim(itoa(iforcing_zsym)))
      endselect
!
!  possibly multiply forcing by z-profile
!  (This stuff is now supposed to be done in initialize; keep for now)
!
!-    if (height_ff/=0.) then
!-      if (lroot) print*,'forcing_hel: include z-profile'
!-      tmpz=(z/height_ff)**2
!-      fz=fz*exp(-tmpz**5/max(1.-tmpz,1e-5))
!-    endif
!
! need to discuss with axel
!
      if (rcyl_ff/=0.) then
!       if (lroot) print*,'fxyz_coefs_hel: applying sgn(z)*xi(r) profile'
        !
        ! only z-dependent part can be done here; radial stuff needs to go
!       ! into the loop
        !
        fz = fz*profz_k
      endif

      if (ip<=5) then
        print*,'fxyz_coefs_hel: fx=',fx
        print*,'fxyz_coefs_hel: fy=',fy
        print*,'fxyz_coefs_hel: fz=',fz
      endif
!
    endsubroutine fxyz_coefs_hel
!***********************************************************************
    subroutine forcing_hel(f)
!
!  Add helical forcing function, using a set of precomputed wavevectors.
!  The relative helicity of the forcing function is determined by the factor
!  sigma, called here also relhel. If it is +1 or -1, the forcing is a fully
!  helical Beltrami wave of positive or negative helicity. For |relhel| < 1
!  the helicity less than maximum. For relhel=0 the forcing is nonhelical.
!  The forcing function is now normalized to unity (also for |relhel| < 1).
!
!  10-apr-00/axel: coded
!   3-sep-02/axel: introduced k1_ff, to rescale forcing function if k1/=1.
!  25-sep-02/axel: preset force_ampl to unity (in case slope is not controlled)
!   9-nov-02/axel: corrected normalization factor for the case |relhel| < 1.
!  23-feb-10/axel: added helicity profile with finite second derivative.
!  13-jun-13/axel: option of symmetry of forcing function about z direction
!  06-dec-13/nishant: made kkx etc allocatable
!  20-aug-14/MR: discard wavevectors analogously to forcing_irrot
!  21-jan-15/MR: changes for use of reference state.
!  26-nov-15/AB: The cs0eff factor used for normalization can now be set.
!   4-oct-17/MR: outsourced forcing_coefs_hel for calculation of coefficients
!                which can be obtained outside an mn-loop. 
!                Moved density calculation into mn-loop where it belongs.
!                Spotted possible bug: for lforcing2_curl=T, calculation of forcing_rhs 
!                seems not to be meaningful.
!  12-jan-18/axel: added periodic forcing for omega_ff /= 0.
!   3-aug-22/axel: added omega_double_ff for second forcing function
!
      use Mpicomm, only: mpireduce_sum
      use Sub
      use EquationOfState, only: rho0
      use DensityMethods, only: getrho1
      use Diagnostics, only: sum_mn_name
!
      real, dimension (mx,my,mz,mfarray), intent(INOUT) :: f

      real, dimension (nx) :: rho1,ruf,rho,force_ampl
      real, dimension (nx,3) :: variable_rhs,forcing_rhs,forcing_rhs2
      real, dimension (nx,3) :: forcing_rhs_old,forcing_rhs2_old
      real, dimension (nx,3) :: force_all
      real, dimension (3), save :: fda, fda2, fda_old, fda2_old
      complex, dimension (mx), save :: fx=0., fx2=0., fx_old=0., fx2_old=0.
      complex, dimension (my), save :: fy=0., fy2=0., fy_old=0., fy2_old=0.
      complex, dimension (mz), save :: fz=0., fz2=0., fz_old=0., fz2_old=0.
      real, dimension (3), save :: coef1,coef2,coef3,coef1b,coef2b,coef3b
      integer :: j,jf,j2f
      complex, dimension (nx) :: fxyz, fxyz2, fxyz_old, fxyz2_old
      real :: profyz, pforce, qforce, aforce, tmp
      real, dimension(3) :: profyz_hel_coef2, profyz_hel_coef2b
!
!  Compute forcing coefficients.
!
      if (t>=tsforce) then
        fx_old=fx
        fy_old=fy
        fz_old=fz
        fda_old=fda
        call forcing_coefs_hel(coef1,coef2,coef3,fda,fx,fy,fz)
!
!  Possibility of reading in data for second forcing function and
!  computing the relevant coefficients (fx2,fy2,fz2,fda2) here.
!  Unlike coef[1-3], where we add the letter b, we add a 2 for the
!  other coefficients for better readibility.
!
        if (lforcing_coefs_hel_double) then
          fx2_old=fx2
          fy2_old=fy2
          fz2_old=fz2
          fda2_old=fda2
          call forcing_coefs_hel2(coef1b,coef2b,coef3b,fda2,fx2,fy2,fz2)
        elseif (lmhd_forcing) then
          fx2_old=fx2
          fy2_old=fy2
          fz2_old=fz2
          fda2_old=fda2
          call forcing_coefs_hel(coef1b,coef2b,coef3b,fda2,fx2,fy2,fz2)
        endif
!
!  By default, dtforce=0, so new forcing is applied at every time step.
!  Alternatively, it can be set to any other time, so the forcing is
!  updated every dtforce.
!
        tsforce=t+dtforce
      endif
!
!  weight factor, pforce is the factor for the new term,
!  and qforce is that for the outgoing term.
!  If no outphasing is involved, pforce=1. and qforce=0.
!
      if (dtforce==0.) then
        pforce=1.
        qforce=0.
      else
        aforce=sin(   pi*(tsforce-t)/dtforce)**2*dtforce_ampl+(1.-dtforce_ampl)
        pforce=cos(.5*pi*(tsforce-t)/dtforce)**2*aforce
        qforce=sin(.5*pi*(tsforce-t)/dtforce)**2*aforce
      endif
!
!  loop the two cases separately, so we don't check for r_ff during
!  each loop cycle which could inhibit (pseudo-)vectorisation
!  calculate energy input from forcing; must use lout (not ldiagnos)
!
!  Here standard case. The case rcyl_ff==0 is to be removed sometime.
!
      if (rcyl_ff == 0) then
        lfirstpoint = .true.
        do n=n1,n2
          do m=m1,m2
!
!  Compute useful shorthands for primary forcing function.
!
            profyz=profy_ampl(m)*profz_ampl(n)
            profyz_hel_coef2=profy_hel(m)*profz_hel(n)*coef2
!
!  Compute the combined complex forcing function fxyz.
!  If lforcing_coefs_hel_double is set, we also add a
!  contribution from fxyz2=fx2(l1:l2)*fy2(m)*fz2(n),
!  which is weighted with factor qdouble_profile(n).
!  qdouble_profile turns on fxyz2 in the upper parts.
!
            fxyz=fx(l1:l2)*fy(m)*fz(n)
            fxyz_old=fx_old(l1:l2)*fy_old(m)*fz_old(n)
            force_ampl=profx_ampl*profyz
!
!  Do the same for secondary forcing function.
!
            if (lforcing_coefs_hel_double.or.lmhd_forcing) then
              profyz_hel_coef2b=profy_hel(m)*profz_hel(n)*coef2b
              fxyz2=fx2(l1:l2)*fy2(m)*fz2(n)
              fxyz2_old=fx2_old(l1:l2)*fy2_old(m)*fz2_old(n)
            endif
!
!  Possibility of compute work done by forcing.
!
            if (lwork_ff) &
              force_ampl=force_ampl*calc_force_ampl(f,fx,fy,fz,profyz*cmplx(coef1,profyz_hel_coef2))
!
!  In the past we always forced the du/dt, but in some cases
!  it may be better to force rho*du/dt (if lmomentum_ff=.true.)
!  For compatibility with earlier results, lmomentum_ff=.false. by default.
!
            if (ldensity) then
              if (lmomentum_ff.or.lout) then
                call getrho1(f(:,m,n,ilnrho),rho1)
                if (lmomentum_ff) force_ampl=force_ampl*rho1
              endif
            endif

            do j=1,3
              if (lactive_dimension(j)) then
!
!  Primary forcing function: assemble here forcing_rhs(:,j).
!  Add here possibility of periodic forcing proportional to cos(om*t).
!  By default, omega_ff=0.
!
                forcing_rhs(:,j) = force_ampl*fda(j)*cos(omega_ff*t) &
                                  *real(cmplx(coef1(j),profx_hel*profyz_hel_coef2(j))*fxyz)

                if (qforce/=0.) &
                  forcing_rhs_old(:,j) = force_ampl*fda_old(j)*cos(omega_ff*t) &
                                        *real(cmplx(coef1(j),profx_hel*profyz_hel_coef2(j))*fxyz_old)

                if (lmhd_forcing) then
                  forcing_rhs2(:,j) = force_ampl*fda2(j)*cos(omega_ff*t) & 
                                     *real(cmplx(coef1b(j),profx_hel*profyz_hel_coef2b(j))*fxyz2)
                  if (qforce/=0.) &
                    forcing_rhs2_old(:,j) = force_ampl*fda2_old(j)*cos(omega_ff*t) &
                                           *real(cmplx(coef1b(j),profx_hel*profyz_hel_coef2b(j))*fxyz2_old)
                endif
!
!  Possibility of adding second forcing function.
! 
                if (lforcing_coefs_hel_double) then
                  forcing_rhs(:,j) = (1.-qdouble_profile(n))*forcing_rhs(:,j) &
                                    +qdouble_profile(n)*force_ampl*fda2(j)*cos(omega_double_ff*t) &
                                    *real(cmplx(coef1b(j),profx_hel*profyz_hel_coef2b(j))*fxyz2)
!
                  if (qforce/=0.) &
                    forcing_rhs_old(:,j) = (1.-qdouble_profile(n))*forcing_rhs_old(:,j) &
                                          +qdouble_profile(n)*force_ampl*fda2_old(j)*cos(omega_double_ff*t) &
                                          *real(cmplx(coef1b(j),profx_hel*profyz_hel_coef2b(j))*fxyz2_old)
                endif
!
!  assemble new combination
!
                if (qforce/=0.) forcing_rhs(:,j) = pforce*forcing_rhs(:,j) + qforce*forcing_rhs_old(:,j)
!
! put force into auxiliary variable, if requested
!
                if (lff_as_aux) f(l1:l2,m,n,iff+j-1) = f(l1:l2,m,n,iff+j-1)+forcing_rhs(:,j)
!
!  Compute additional forcing function (used for velocity if crosshel=1).
!  It can optionally be the same. Alternatively, one has to set crosshel=1.
!
                if (lforcing2_same) then
                  forcing_rhs2(:,j)=forcing_rhs(:,j)
                elseif (lforcing2_curl) then
!
!  Something seems to be missing here as real(cmplx(0.,coef3(j))) is zero.
!
                  forcing_rhs2(:,j) = force_ampl*real(cmplx(0.,coef3(j)))*fda(j) 
                else
                  forcing_rhs2(:,j) = force_ampl*real(cmplx(0.,coef3(j))*fxyz)*fda(j)
                endif
!
!  Choice of different possibilities.
!  Here is the decisive step where forcing is applied.
!  Added possibility of linearly ramping down the forcing in time.
!
                if (ifff/=0) then

                  jf=j+ifff-1
                  j2f=j+i2fff-1
!
                  if (lhelical_test) then
                    f(l1:l2,m,n,jf)=forcing_rhs(:,j)
                  else
                    if (lforce_ramp_down) then
                      tmp=max(0d0,1d0+min(0d0,(tforce_ramp_down-t)/tauforce_ramp_down))
                      f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)*force1_scl*tmp
                    else
                      f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)*force1_scl
                    endif
!
!  Allow here for forcing both in u and in b=curla. In that case one sets
!  lhydro_forcing=T and lmagnetic_forcing=T.
!
                    if (lmhd_forcing) then 
                      f(l1:l2,m,n,j2f)=f(l1:l2,m,n,j2f)+forcing_rhs2(:,j)*force2_scl
                    elseif (lcrosshel_forcing) then
                      f(l1:l2,m,n,j2f)=f(l1:l2,m,n,j2f)+forcing_rhs(:,j)*force2_scl
                    endif
!
!  Forcing with enhanced xx correlation.
!
                    if (j==1.and.lxxcorr_forcing) &
                      f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,1)*force2_scl
                  endif
                endif
!
!  If one of the testfield methods is used, we need to add a forcing term
!  in one of the auxiliary equations. Their location is denoted by jtest_aa0
!  and jtest_uu0 for the testfield and testflow equations, respectively.
!  In the testflow module, jtest_uu0=1 is used, while in the testfield_nonlinear_z
!  module jtest_uu0=5 is used, so for now we give them by hand.
!
                if (ltestfield_forcing) then
                  jf=j+iaatest-1+3*(jtest_aa0-1)
                  f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)*force1_scl
                endif
                if (ltestflow_forcing) then
                  jf=j+iuutest-1+3*(jtest_uu0-1)
                  if (ltestfield_forcing) then
                    f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs2(:,j)*force2_scl
                  else
                    f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)*force2_scl
                  endif
                endif
              endif
            enddo
!
!  Forcing with enhanced xy correlation.
!  Can only come outside the previous j loop.
!  This is apparently not yet applied to testfield or testflow forcing.
!
            if (ifff/=0.and..not.lhelical_test.and.lxycorr_forcing) then
              do j=1,3
                if (lactive_dimension(j)) then
                  jf=j+ifff-1
                  if (j==1) f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf) + forcing_rhs(:,2)*force2_scl
                  if (j==2) f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf) + forcing_rhs(:,1)*force2_scl
                endif
              enddo
            endif
!
!  Sum up.
!
            if (lout) then
              if (idiag_rufm  /=0 .or. idiag_ruxfxm/=0 .or. idiag_ruxfym/=0 .or. &
                  idiag_ruyfxm/=0 .or. idiag_ruyfym/=0 .or. idiag_ruzfzm/=0     ) then
!
!  Compute density.
!
                if (ldensity) then
                  rho=1./rho1
                else
                  rho=rho0
                endif
!
!  Evaluate the various choices.
!
                if (idiag_rufm/=0) then
!
!  Compute rhs.
!
                  variable_rhs=f(l1:l2,m,n,iffx:iffz)          !MR: here already updated!
                  call multsv_mn(rho/dt,forcing_rhs,force_all)
                  call dot_mn(variable_rhs,force_all,ruf)
                  call sum_mn_name(ruf,idiag_rufm)
                endif

                if (idiag_ruxfxm/=0) call sum_mn_name(rho*f(l1:l2,m,n,iux)*forcing_rhs(:,1),idiag_ruxfxm)
                if (idiag_ruxfym/=0) call sum_mn_name(rho*f(l1:l2,m,n,iux)*forcing_rhs(:,2),idiag_ruxfym)
                if (idiag_ruyfxm/=0) call sum_mn_name(rho*f(l1:l2,m,n,iuy)*forcing_rhs(:,1),idiag_ruyfxm)
                if (idiag_ruyfym/=0) call sum_mn_name(rho*f(l1:l2,m,n,iuy)*forcing_rhs(:,2),idiag_ruyfym)
                if (idiag_ruzfzm/=0) call sum_mn_name(rho*f(l1:l2,m,n,iuz)*forcing_rhs(:,3),idiag_ruzfzm)
              endif
            endif
!
!  End of mn loop.
!
            lfirstpoint = .false.
          enddo
        enddo
      else
!
!  Radial profile, but this is old fashioned and probably no longer used.
!
        call fatal_error('forcing_hel','check that radial profile with rcyl_ff/=0. works')
        do j=1,3
          if (lactive_dimension(j)) then
            jf=j+ifff-1
            do n=n1,n2
              do m=m1,m2
                if (ldensity) then
                  call getrho1(f(:,m,n,ilnrho),rho1)
                else
                  rho1=1./rho0
                endif
                forcing_rhs(:,j) = rho1*profx_ampl*profy_ampl(m)*profz_ampl(n) &
                                  *real(cmplx(coef1(j),profy_hel(m)*profz_k(n)*coef2(j)) &
                                        *fx(l1:l2)*fy(m)*fz(n))

                if (lhelical_test) then
                  f(l1:l2,m,n,jf)=forcing_rhs(:,j)
                else
                  f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)
                endif

              enddo
            enddo
          endif
        enddo
      endif
!
!  For printouts, rufm needs to be communicated.
!
      if (lout) then
        if (idiag_rufm/=0)   call mpireduce_sum(fname(idiag_rufm)/nwgrid,fname(idiag_rufm))
        if (idiag_ruxfxm/=0) call mpireduce_sum(fname(idiag_ruxfxm)/nwgrid,fname(idiag_ruxfxm))
        if (idiag_ruxfym/=0) call mpireduce_sum(fname(idiag_ruxfym)/nwgrid,fname(idiag_ruxfym))
        if (idiag_ruyfxm/=0) call mpireduce_sum(fname(idiag_ruyfxm)/nwgrid,fname(idiag_ruyfxm))
        if (idiag_ruyfym/=0) call mpireduce_sum(fname(idiag_ruyfym)/nwgrid,fname(idiag_ruyfym))
        if (idiag_ruzfzm/=0) call mpireduce_sum(fname(idiag_ruzfzm)/nwgrid,fname(idiag_ruzfzm))
      endif
!
      if (ip<=9) print*,'forcing_hel: forcing OK'
!
    endsubroutine forcing_hel
!***********************************************************************
    subroutine forcing_hel_kprof(f)
!
!  Add helical forcing function, using a set of precomputed wavevectors.
!  The relative helicity of the forcing function is determined by the factor
!  sigma, called here also relhel. If it is +1 or -1, the forcing is a fully
!  helical Beltrami wave of positive or negative helicity. For |relhel| < 1
!  the helicity less than maximum. For relhel=0 the forcing is nonhelical.
!  The forcing function is now normalized to unity (also for |relhel| < 1).
!
!  30-jan-11/axel: adapted from forcing_hel and added z-dependent scaling of k
!  06-dec-13/nishant: made kkx etc allocatable
!
      use Diagnostics, only: sum_mn_name
      use EquationOfState, only: rho0
      use General, only: random_number_wrapper
      use Sub, only: del2v_etc,curl,cross,dot,dot2
      use DensityMethods, only: getrho1, getrho
      use Mpicomm, only: mpireduce_sum
!
      real :: phase,ffnorm
      real, dimension (2) :: fran
      real, dimension (nx) :: rho1,ff,uf,of,qf,rho
      real, dimension (mz) :: kfscl
      real, dimension (nx,3) :: variable_rhs,forcing_rhs,forcing_rhs2
      real, dimension (nx,3) :: fda,uu,oo,bb,fxb,curlo
      real, dimension (mx,my,mz,mfarray) :: f
      complex, dimension (mx) :: fx
      complex, dimension (my) :: fy
      complex, dimension (mz) :: fz
      real, dimension (3) :: coef1,coef2,coef3
      integer :: ik,j,jf,j2f
      real :: kx0,kx,ky,kz,k2,k,force_ampl,pi_over_Lx
      real :: ex,ey,ez,kde,sig,fact,kex,key,kez,kkex,kkey,kkez
      real, dimension(3) :: e1,e2,ee,kk
      real :: norm,phi
      real :: fd,fd2
!
!  generate random coefficients -1 < fran < 1
!  ff=force*Re(exp(i(kx+phase)))
!  |k_i| < akmax
!
      call random_number_wrapper(fran,CHANNEL=channel_force)
      phase=pi*(2*fran(1)-1.)
      ik=nk*(.9999*fran(2))+1
      if (ip<=6) then
        print*,'forcing_hel_kprof: ik,phase=',ik,phase
        print*,'forcing_hel_kprof: kx,ky,kz=',kkx(ik),kky(ik),kkz(ik)
        print*,'forcing_hel_kprof: dt=',dt
      endif
!
!  compute z-dependent scaling factor, regard kav as kmax and write
!    kinv=1./kav+(1.-1./kav)*(ztop-z)/(ztop-zbot)
!    kfscl=exp((z(n)-xyz0(3))/Lxyz(3))/kav
!  i.e. divide by kav, so that k=1 when z=zbot.
!
      kfscl=1./(1.+(kav-1.)*(xyz0(3)+Lxyz(3)-z)/Lxyz(3))
      !kfscl=exp((z-xyz0(3))/Lxyz(3))/kav
!
!  Loop over all k, but must have the call to random_number_wrapper
!  outside this loop.
!
      call random_number_wrapper(phi,CHANNEL=channel_force); phi = phi*2*pi
      lfirstpoint = .true.
      do n=n1,n2
!
!  normally we want to use the wavevectors as they are,
!  but in some cases, e.g. when the box is bigger than 2pi,
!  we want to rescale k so that k=1 now corresponds to a smaller value.
!
        if (lscale_kvector_fac) then
          kx0=kkx(ik)*scale_kvectorx
          ky=kky(ik)*scale_kvectory
          kz=kkz(ik)*scale_kvectorz
          pi_over_Lx=0.5
        elseif (lscale_kvector_tobox) then
          kx0=kkx(ik)*(2.*pi/Lxyz(1))
          ky=kky(ik)*(2.*pi/Lxyz(2))
          kz=kkz(ik)*(2.*pi/Lxyz(3))
          pi_over_Lx=pi/Lxyz(1)
        else
          kx0=kkx(ik)
          ky=kky(ik)
          kz=kkz(ik)
          pi_over_Lx=0.5
        endif
!
!  scale all components of k
!
        kx0=kx0*kfscl(n)
        ky=ky*kfscl(n)
        kz=kz*kfscl(n)
!
!  in the shearing sheet approximation, kx = kx0 - St*k_y.
!  Here, St=-deltay/Lx. However, to stay near kx0, we ignore
!  integer shifts.
!
        if (Sshear==0.) then
          kx=kx0
        else
          if (lshearing_adjust_old) then
            kx=kx0+ky*deltay/Lx
          else
            kx=kx0+mod(ky*deltay/Lx-pi_over_Lx,2.*pi_over_Lx)+pi_over_Lx
          endif
        endif
!
!  compute k^2 and output wavenumbers
!
        k2=kx**2+ky**2+kz**2
        k=sqrt(k2)
        if (ip<4) then
          open(89,file='forcing_hel_kprof_output.dat',position='append')
          write(89,'(6f10.5)') k,kx0,kx,ky,kz,deltay
          close(89)
        endif
!
!  Find e-vector:
!  Start with old method (not isotropic) for now.
!  Pick e1 if kk not parallel to ee1. ee2 else.
!
        if ((ky==0).and.(kz==0)) then
          ex=0; ey=1; ez=0
        else
          ex=1; ey=0; ez=0
        endif
        if (.not. old_forcing_evector) then
   !
   !  Isotropize ee in the plane perp. to kk by
   !  (1) constructing two basis vectors for the plane perpendicular
   !      to kk, and
   !  (2) choosing a random direction in that plane (angle phi)
   !  Need to do this in order for the forcing to be isotropic.
   !
          kk = (/kx, ky, kz/)
          ee = (/ex, ey, ez/)
          call cross(kk,ee,e1)
          call dot2(e1,norm); e1=e1/sqrt(norm) ! e1: unit vector perp. to kk
          call cross(kk,e1,e2)
          call dot2(e2,norm); e2=e2/sqrt(norm) ! e2: unit vector perp. to kk, e1
          ee = cos(phi)*e1 + sin(phi)*e2
          ex=ee(1); ey=ee(2); ez=ee(3)
        endif
!
!  k.e
!
        call dot(kk,ee,kde)
!
!  k x e
!
        kex=ky*ez-kz*ey
        key=kz*ex-kx*ez
        kez=kx*ey-ky*ex
!
!  k x (k x e)
!
        kkex=ky*kez-kz*key
        kkey=kz*kex-kx*kez
        kkez=kx*key-ky*kex
!
!  ik x (k x e) + i*phase
!
!  Normalize ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!
!  This does already include the new sqrt(2) factor (missing in B01).
!  So, in order to reproduce the 0.1 factor mentioned in B01
!  we have to set force=0.07.
!
!  Furthermore, for |relhel| < 1, sqrt(2) should be replaced by
!  sqrt(1.+relhel**2). This is done now (9-nov-02).
!  This means that the previous value of force=0.07 (for relhel=0)
!  should now be replaced by 0.05.
!
!  Note: kav is not to be scaled with k1_ff (forcing should remain
!  unaffected when changing k1_ff).
!
        ffnorm=sqrt(1.+relhel**2)*k*sqrt(k2-kde**2)/sqrt(kav*cs0**3)*(k/kav)**slope_ff
        if (ip<=9) then
          print*,'forcing_hel_kprof: k,kde,ffnorm,kav=',k,kde,ffnorm,kav
          print*,'forcing_hel_kprof: k*sqrt(k2-kde**2)=',k*sqrt(k2-kde**2)
        endif
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
        fact=force/ffnorm*sqrt(dt)
!
!  The wavevector is for the case where Lx=Ly=Lz=2pi. If that is not the
!  case one needs to scale by 2pi/Lx, etc.
!
        fx=exp(cmplx(0.,kx*k1_ff*x+phase))*fact
        fy=exp(cmplx(0.,ky*k1_ff*y))
        fz=exp(cmplx(0.,kz*k1_ff*z))
!
!  possibly multiply forcing by sgn(z) and radial profile
!
        if (rcyl_ff/=0.) then
          if (lroot) print*,'forcing_hel_kprof: applying sgn(z)*xi(r) profile'
          !
          ! only z-dependent part can be done here; radial stuff needs to go
          ! into the loop
          !
          fz = fz*profz_k
        endif
!
        if (ip<=5) then
          print*,'forcing_hel_kprof: fx=',fx
          print*,'forcing_hel_kprof: fy=',fy
          print*,'forcing_hel_kprof: fz=',fz
        endif
!
!  prefactor; treat real and imaginary parts separately (coef1 and coef2),
!  so they can be multiplied by different profiles below.
!
        coef1(1)=k*kex; coef2(1)=relhel*kkex; coef3(1)=crosshel*k*kkex
        coef1(2)=k*key; coef2(2)=relhel*kkey; coef3(2)=crosshel*k*kkey
        coef1(3)=k*kez; coef2(3)=relhel*kkez; coef3(3)=crosshel*k*kkez
        if (ip<=5) print*,'forcing_hel_kprof: coef=',coef1,coef2
!
! An attempt to implement anisotropic forcing using direction
! dependent forcing amplitude. Activated only if force_strength,
! describing the anisotropic part of the forcing, is
! nonzero. force_direction, which is a vector, defines the preferred
! direction of forcing. The expression for the forcing amplitude used
! at the moment is:
!
!  f(i)=f0*[1+epsilon(delta_ij*(k(i)*fd(j))/(|k||fd|))^2*fd(i)/|fd|]
!
! here f0 and fd are shorthand for force and forcing_direction,
! respectively, and epsilon=force_strength/force.
!
        if (force_strength/=0.) then
          call dot(force_direction,force_direction,fd2)
          fd=sqrt(fd2)
          do j=1,3
            fda(:,j) = 1. + (force_strength/force)*(kk(j)*force_direction(j)/(k*fd))**2 &
                      *force_direction(j)/fd
          enddo
        else
          fda = 1.
        endif
!
!  loop the two cases separately, so we don't check for r_ff during
!  each loop cycle which could inhibit (pseudo-)vectorisation
!  calculate energy input from forcing; must use lout (not ldiagnos)
!
        force_ampl=1.0
        if (rcyl_ff == 0) then       ! no radial profile
          do m=m1,m2
!
!  In the past we always forced the du/dt, but in some cases
!  it may be better to force rho*du/dt (if lmomentum_ff=.true.)
!  For compatibility with earlier results, lmomentum_ff=.false. by default.
!
            if (ldensity.and.lmomentum_ff) then
              call getrho1(f(:,m,n,ilnrho),rho1)
            else
              rho1=1./rho0
            endif
            if (lwork_ff) force_ampl=calc_force_ampl(f,fx,fy,fz,profy_ampl(m)*profz_ampl(n) &
                                                     *cmplx(coef1,profy_hel(m)*profz_hel(n)*coef2))
            variable_rhs=f(l1:l2,m,n,iffx:iffz)
            do j=1,3
              if (lactive_dimension(j)) then
!
                forcing_rhs(:,j)=rho1*profx_ampl*profy_ampl(m)*profz_ampl(n)*force_ampl &
                  *real(cmplx(coef1(j),profx_hel*profy_hel(m)*profz_hel(n)*coef2(j)) &
                        *fx(l1:l2)*fy(m)*fz(n))*fda(:,j)
!
                forcing_rhs2(:,j)=rho1*profx_ampl*profy_ampl(m)*profz_ampl(n)*force_ampl &
                                  *real(cmplx(0.,coef3(j))*fx(l1:l2)*fy(m)*fz(n))*fda(:,j)
!
                if (ifff/=0) then
!
                  jf=j+ifff-1
                  j2f=j+i2fff-1
!
                  if (lhelical_test) then
                    f(l1:l2,m,n,jf)=forcing_rhs(:,j)
                  else
                    f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)*force1_scl
!
!  allow here for forcing both in u and in b=curla. In that case one sets
!  lhydro_forcing=T and lmagnetic_forcing=T.
!
                    if (lmhd_forcing) then
                      f(l1:l2,m,n,j2f)=f(l1:l2,m,n,j2f)+forcing_rhs2(:,j)*force2_scl
                    elseif (lcrosshel_forcing) then
                      f(l1:l2,m,n,j2f)=f(l1:l2,m,n,j2f)+forcing_rhs(:,j)*force2_scl
                    endif
                  endif
!
                endif
!
!  If one of the testfield methods is used, we need to add a forcing term
!  in one of the auxiliary equations. Their location is denoted by jtest_aa0
!  and jtest_uu0 for the testfield and testflow equations, respectively.
!  In the testflow module, jtest_uu0=1 is used, while in the testfield_nonlinear
!  module jtest_uu0=5 is used, so for now we give them by hand.
!
                if (ltestfield_forcing) then
                  jf=j+iaatest-1+3*(jtest_aa0-1)
                  f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)*force1_scl
                endif
                if (ltestflow_forcing) then
                  jf=j+iuutest-1+3*(jtest_uu0-1)
                  if (ltestfield_forcing) then
                    f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs2(:,j)*force2_scl
                  else
                    f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)*force2_scl
                  endif
                endif
              endif
            enddo
!
! Diagnostics: 
!
            if (lout) then
              if (ldensity.and.idiag_rufm/=0) then

                if (lmomentum_ff) then
                  rho=1./rho1
                else
                  call getrho(f(:,m,n,ilnrho),rho)
                endif

                uu=f(l1:l2,m,n,iux:iuz)
                call dot(uu,forcing_rhs,uf)
                call sum_mn_name(rho*uf,idiag_rufm)

              endif
              if (idiag_ufm/=0) then
                uu=f(l1:l2,m,n,iux:iuz)
                call dot(uu,forcing_rhs,uf)
                call sum_mn_name(uf,idiag_ufm)
              endif
              if (idiag_ofm/=0) then
                call curl(f,iuu,oo)
                call dot(oo,forcing_rhs,of)
                call sum_mn_name(of,idiag_ofm)
              endif
              if (idiag_qfm/=0) then
                call del2v_etc(f,iuu,curlcurl=curlo)
                call dot(curlo,forcing_rhs,qf)
                call sum_mn_name(qf,idiag_qfm)
              endif
              if (idiag_ffm/=0) then
                call dot2(forcing_rhs,ff)
                call sum_mn_name(ff,idiag_ffm)
              endif
              if (lmagnetic) then
!               if (idiag_fxbxm/=0.or.idiag_fxbym/=0.or.idiag_fxbzm/=0) then
!                 call curl(f,iaa,bb)
!                 call cross(forcing_rhs,bb,fxb)
!                 call sum_mn_name(fxb(:,1),idiag_fxbxm)
!                 call sum_mn_name(fxb(:,2),idiag_fxbym)
!                 call sum_mn_name(fxb(:,3),idiag_fxbzm)
!               endif
              endif
            endif
            lfirstpoint = .false.
          enddo

        else
!
!  Radial profile, but this is old fashioned and probably no longer used.
!
          call fatal_error('forcing_hel_kprof','check that radial profile with rcyl_ff works')
          do j=1,3
            if (lactive_dimension(j)) then
              jf=j+ifff-1
              sig = relhel*profz_k(n)
              coef1(1)=cmplx(k*kex,sig*kkex)
              coef1(2)=cmplx(k*key,sig*kkey)
              coef1(3)=cmplx(k*kez,sig*kkez)
              do m=m1,m2
!
!  In the past we always forced the du/dt, but in some cases
!  it may be better to force rho*du/dt (if lmomentum_ff=.true.)
!  For compatibility with earlier results, lmomentum_ff=.false. by default.
!
                if (ldensity.and.lmomentum_ff) then
                  call getrho1(f(:,m,n,ilnrho),rho1)
                else
                  rho1=1./rho0
                endif

                forcing_rhs(:,j)= rho1*profx_ampl*profy_ampl(m)*profz_ampl(n) &
                                 *real(cmplx(coef1(j),profy_hel(m)*coef2(j))*fx(l1:l2)*fy(m)*fz(n))
                if (lhelical_test) then
                  f(l1:l2,m,n,jf)=forcing_rhs(:,j)
                else
                  f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)
                endif
                lfirstpoint = .false.
              enddo
            endif
          enddo
        endif
!
      enddo
!
      if (ip<=9) print*,'forcing_hel_kprof: forcing OK'
!
!  For printouts, rufm needs to be communicated.
!
      if (lout) then
        if (idiag_rufm/=0) call mpireduce_sum(fname(idiag_rufm)/nwgrid,fname(idiag_rufm))
        if (idiag_ufm/=0)  call mpireduce_sum(fname(idiag_ufm)/nwgrid,fname(idiag_ufm))
        if (idiag_ofm/=0)  call mpireduce_sum(fname(idiag_ofm)/nwgrid,fname(idiag_ofm))
        if (idiag_qfm/=0)  call mpireduce_sum(fname(idiag_qfm)/nwgrid,fname(idiag_qfm))
        if (idiag_ffm/=0)  call mpireduce_sum(fname(idiag_ffm)/nwgrid,fname(idiag_ffm))
!       if (lmagnetic) then
!         if (idiag_fxbxm/=0) call mpireduce_sum(fname(idiag_fxbxm)/nwgrid,fname(idiag_fxbxm))
!         if (idiag_fxbym/=0) call mpireduce_sum(fname(idiag_fxbym)/nwgrid,fname(idiag_fxbym))
!         if (idiag_fxbzm/=0) call mpireduce_sum(fname(idiag_fxbzm)/nwgrid,fname(idiag_fxbzm))
!       endif
      endif
!
    endsubroutine forcing_hel_kprof
!***********************************************************************
    subroutine forcing_hel_both(f)
!
!  Add helical forcing function, using a set of precomputed wavevectors.
!  The relative helicity of the forcing function is determined by the factor
!  sigma, called here also relhel. If it is +1 or -1, the forcing is a fully
!  helical Beltrami wave of positive or negative helicity. For |relhel| < 1
!  the helicity less than maximum. For relhel=0 the forcing is nonhelical.
!  The forcing function is now normalized to unity (also for |relhel| < 1).
!  This adds positive helical forcing to the "northern hemisphere" (y above the
!  midplane and negative helical forcing to the "southern hemisphere". The
!  two forcing are merged at the "equator" (midplane) where both are smoothly
!  set to zero.
!
!  22-sep-08/dhruba: adapted from forcing_hel
!   6-oct-09/MR: according to Axel, this routine is now superseded by forcing_hel and should be deleted
!  06-dec-13/nishant: made kkx etc allocatable
!
      use General, only: random_number_wrapper
      use Sub
!
      real :: phase,ffnorm
      real, dimension (2) :: fran
      real, dimension (nx) :: rho1
      real, dimension (nx,3) :: variable_rhs,forcing_rhs
      real, dimension (mx,my,mz,mfarray) :: f
      complex, dimension (mx) :: fx
      complex, dimension (my) :: fy
      complex, dimension (mz) :: fz
      real, dimension (3) :: coef1,coef2
      integer :: ik,j,jf
      real :: kx0,kx,ky,kz,k2,k
      real :: ex,ey,ez,kde,fact,kex,key,kez,kkex,kkey,kkez
      real, dimension(3) :: e1,e2,ee,kk
      real :: norm,phi,pi_over_Lx
!
!  additional stuff for test fields
!
      integer :: jtest
!
!  generate random coefficients -1 < fran < 1
!  ff=force*Re(exp(i(kx+phase)))
!  |k_i| < akmax
!
      call random_number_wrapper(fran,CHANNEL=channel_force)
      phase=pi*(2*fran(1)-1.)
      ik=nk*(.9999*fran(2))+1
      if (ip<=6) then 
        print*,'forcing_hel_both: ik,phase=',ik,phase
        print*,'forcing_hel_both: kx,ky,kz=',kkx(ik),kky(ik),kkz(ik)
        print*,'forcing_hel_both: dt=',dt
      endif
!
!  normally we want to use the wavevectors as they are,
!  but in some cases, e.g. when the box is bigger than 2pi,
!  we want to rescale k so that k=1 now corresponds to a smaller value.
!
      if (lscale_kvector_fac) then
        kx0=kkx(ik)*scale_kvectorx
        ky=kky(ik)*scale_kvectory
        kz=kkz(ik)*scale_kvectorz
        pi_over_Lx=0.5
      elseif (lscale_kvector_tobox) then
        kx0=kkx(ik)*(2.*pi/Lxyz(1))
        ky=kky(ik)*(2.*pi/Lxyz(2))
        kz=kkz(ik)*(2.*pi/Lxyz(3))
        pi_over_Lx=pi/Lxyz(1)
      else
        kx0=kkx(ik)
        ky=kky(ik)
        kz=kkz(ik)
        pi_over_Lx=0.5
      endif
!
!  in the shearing sheet approximation, kx = kx0 - St*k_y.
!  Here, St=-deltay/Lx
!
      if (Sshear==0.) then
        kx=kx0
      else
        kx=kx0+ky*deltay/Lx
      endif
!
      if (headt.or.ip<5) print*, 'forcing_hel_both: kx0,kx,ky,kz=',kx0,kx,ky,kz
      k2=kx**2+ky**2+kz**2
      k=sqrt(k2)
!
! Find e-vector
!
      !
      ! Start with old method (not isotropic) for now.
      ! Pick e1 if kk not parallel to ee1. ee2 else.
      !
      if ((ky==0).and.(kz==0)) then
        ex=0; ey=1; ez=0
      else
        ex=1; ey=0; ez=0
      endif

      if (.not. old_forcing_evector) then
        !
        !  Isotropize ee in the plane perp. to kk by
        !  (1) constructing two basis vectors for the plane perpendicular
        !      to kk, and
        !  (2) choosing a random direction in that plane (angle phi)
        !  Need to do this in order for the forcing to be isotropic.
        !
        kk = (/kx, ky, kz/)
        ee = (/ex, ey, ez/)
        call cross(kk,ee,e1)
        call dot2(e1,norm); e1=e1/sqrt(norm) ! e1: unit vector perp. to kk
        call cross(kk,e1,e2)
        call dot2(e2,norm); e2=e2/sqrt(norm) ! e2: unit vector perp. to kk, e1
        call random_number_wrapper(phi,CHANNEL=channel_force); phi = phi*2*pi
        ee = cos(phi)*e1 + sin(phi)*e2
        ex=ee(1); ey=ee(2); ez=ee(3)
      endif
!
!  k.e
!
      call dot(kk,ee,kde)
!
!  k x e
!
      kex=ky*ez-kz*ey
      key=kz*ex-kx*ez
      kez=kx*ey-ky*ex
!
!  k x (k x e)
!
      kkex=ky*kez-kz*key
      kkey=kz*kex-kx*kez
      kkez=kx*key-ky*kex
!
!  ik x (k x e) + i*phase
!
!  Normalize ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!  For further details on normalization see forcing_hel subroutine.
!
      ffnorm = sqrt(1.+relhel**2)*k*sqrt(k2-kde**2)/sqrt(kav*cs0**3)*(k/kav)**slope_ff
      if (ip<=9) then
        print*,'forcing_hel_both: k,kde,ffnorm,kav=',k,kde,ffnorm,kav
        print*,'forcing_hel_both: k*sqrt(k2-kde**2)=',k*sqrt(k2-kde**2)
      endif
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
      fact=force/ffnorm*sqrt(dt)
!
!  The wavevector is for the case where Lx=Ly=Lz=2pi. If that is not the
!  case one needs to scale by 2pi/Lx, etc.
!
      fx=exp(cmplx(0.,kx*k1_ff*x+phase))*fact
! only sines, to make sure that the force goes to zero at the equator
      fy=cmplx(0.,sin(ky*k1_ff*y))
      fz=exp(cmplx(0.,kz*k1_ff*z))
!
      if (ip<=5) then
        print*,'forcing_hel_both: fx=',fx
        print*,'forcing_hel_both: fy=',fy
        print*,'forcing_hel_both: fz=',fz
      endif
!
!  prefactor; treat real and imaginary parts separately (coef1 and coef2),
!  so they can be multiplied by different profiles below.
!
      coef1(1)=k*kex; coef2(1)=relhel*kkex
      coef1(2)=k*key; coef2(2)=relhel*kkey
      coef1(3)=k*kez; coef2(3)=relhel*kkez
      if (ip<=5) print*,'forcing_hel_both: coef=',coef1,coef2
!
!  loop the two cases separately, so we don't check for r_ff during
!  each loop cycle which could inhibit (pseudo-)vectorisation
!  calculate energy input from forcing; must use lout (not ldiagnos)
!
     rho1=1.
      do n=n1,n2
        do m=m1,m2
          variable_rhs=f(l1:l2,m,n,iffx:iffz)
          do j=1,3
            if (lactive_dimension(j)) then
              jf=j+ifff-1
              if (y(m)>equator)then
                forcing_rhs(:,j) = rho1*real(cmplx(coef1(j),coef2(j))*fx(l1:l2)*fy(m)*fz(n)) &
                                  *(1.-step(y(m),equator-ck_equator_gap,ck_gap_step)+ &
                                    step(y(m),equator+ck_equator_gap,ck_gap_step))
              else
                forcing_rhs(:,j) = rho1*real(cmplx(coef1(j),-coef2(j))*fx(l1:l2)*fy(m)*fz(n)) &
                                  *(1.-step(y(m),equator-ck_equator_gap,ck_gap_step)+ &
                                    step(y(m),equator+ck_equator_gap,ck_gap_step))
              endif
              if (lhelical_test) then
                f(l1:l2,m,n,jf)=forcing_rhs(:,j)
              else
                f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)
              endif
              if (ltestfield_forcing) then
                jf=j+iaatest-1+3*(jtest_aa0-1)
                f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)
              endif
              if (ltestflow_forcing) then
                jf=j+iuutest-1+3*(jtest_uu0-1)
                f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)
              endif
            endif
          enddo
        enddo
      enddo
!
      if (ip<=9) print*,'forcing_hel_both: forcing OK'
!
    endsubroutine forcing_hel_both
!***********************************************************************
    subroutine forcing_chandra_kendall(f)
!
!  Add helical forcing function in spherical polar coordinate system.
!  25-jul-07/dhruba: adapted from forcing_hel
!
      use EquationOfState, only: cs20
      use General, only: random_number_wrapper
      use Sub
!
      real, dimension (mx,my,mz,mfarray) :: f
      !
      real, dimension(3) :: ee
      real, dimension(nx,3) :: capitalT,capitalS,capitalH,psi
      real, dimension(nx,3,3) :: psi_ij,Tij
      integer :: emm,l,j,jf,Legendrel,lmindex,aindex
      real :: a_ell,anum,adenom,jlm_ff,ylm_ff,rphase1,fnorm,alphar,Balpha,&
              psilm,RYlm,IYlm
      real :: rz,rindex,ralpha,rphase2
      real, dimension(mx) :: Z_psi
!
! This is designed for 5 emm values and for each one 5 ell values. Total 25 values.
!
      if (lhelical_test) then
        if (icklist==nlist_ck) &
          call fatal_error("forcing_chandra_kendall","CK testing: no more values in list")
        icklist=icklist+1
        lmindex=icklist
      else
        call random_number_wrapper(rindex,CHANNEL=channel_force)
        lmindex=nint(rindex*(nlist_ck-1))+1
      endif
!
      emm = cklist(lmindex,1)
      Legendrel = cklist(lmindex,2)
!
      call random_number_wrapper(ralpha,CHANNEL=channel_force)
      aindex=nint(ralpha*2)
      Balpha = cklist(lmindex,3+aindex)
!
! Calculate the "potential" for the helical forcing according to Chandrasekhar & Kendall.
! Construct Z_psi(r).
!
      call random_number_wrapper(rphase1,CHANNEL=channel_force)
      rphase1=rphase1*2*pi

      if (lfastCK) then
        do n=1,mz
          do m=1,my
            psilm = RYlm_list(m,n,lmindex)*cos(rphase1) - &
                    IYlm_list(m,n,lmindex)*sin(rphase1)
            psif(:,m,n) = psilm*Zpsi_list(:,lmindex,aindex+1)
            if (ck_equator_gap/=0) psif(:,m,n)=psif(:,m,n)*profy_ampl(m)
          enddo
        enddo

      else
        call sp_bessely_l(anum,Legendrel,Balpha*xyz0(1))
        call sp_besselj_l(adenom,Legendrel,Balpha*xyz0(1))
        a_ell = -anum/adenom
!        write(*,*) 'dhruba:',anum,adenom,Legendrel,Bessel_alpha,xyz0(1)
        do l=1,mx
          alphar=Balpha*x(l)
          call sp_besselj_l(jlm_ff,Legendrel,alphar)
          call sp_bessely_l(ylm_ff,Legendrel,alphar)
          Z_psi(l) = (a_ell*jlm_ff+ylm_ff)
        enddo
!
        do n=1,mz
          do m=1,my
            call sp_harm_real(RYlm,Legendrel,emm,y(m),z(n))
            call sp_harm_imag(IYlm,Legendrel,emm,y(m),z(n))
            psilm = RYlm*cos(rphase1)-IYlm*sin(rphase1)
            psif(:,m,n) = Z_psi*psilm
            if (ck_equator_gap/=0) psif(:,m,n)=psif(:,m,n)*profy_ampl(m)
          enddo
        enddo
      endif
!
! ----- Now calculate the force from the potential and add this to velocity
! get a random unit vector with three components ee_r, ee_theta, ee_phi
! psi at present is just Z_{ell}^m. We next do a sum over random coefficients
! get random psi.
! ----------now generate and add the force ------------
!

      call random_number_wrapper(rz,CHANNEL=channel_force)
      ee(3) = rz
      call random_number_wrapper(rphase2,CHANNEL=channel_force)
      rphase2 = pi*rphase2
      ee(1) = sqrt(1.-rz*rz)*cos(rphase2)
      ee(2) = sqrt(1.-rz*rz)*sin(rphase2)
      fnorm = fpre*cs20*sqrt(1./(cs0*Balpha))*sqrt(dt)
!     write(*,*) 'dhruba:',fnorm*sqrt(dt),dt,ee(1),ee(2),ee(3)

      do n=n1,n2
        do m=m1,m2
          psi(:,1) = psif(l1:l2,m,n)*ee(1)
          psi(:,2) = psif(l1:l2,m,n)*ee(2)
          psi(:,3) = psif(l1:l2,m,n)*ee(3)
          call gij_psi(psif,ee,psi_ij)
          call curl_mn(psi_ij,capitalT,psi)
          call gij_psi_etc(psif,ee,psi,psi_ij,Tij)
          call curl_mn(Tij,capitalS,capitalT)
          capitalS = float(helsign)*(1./Balpha)*capitalS
          if ((y(m)<pi/2.).or.(lsamesign)) then
            capitalH = capitalT + capitalS
          else
            capitalH = capitalT - capitalS
          endif
          do j=1,3
            jf = iuu+j-1
            if (r_ff /= 0.)  capitalH(:,j)=profx_ampl*capitalH(:,j)
            if (lhelical_test) then
              if (lwrite_psi) then
                f(l1:l2,m,n,jf) = psif(l1:l2,m,n)
              else
                f(l1:l2,m,n,jf) = fnorm*capitalH(:,j)
              endif
            else
!
! stochastic euler scheme of integration [sqrt(dt) is already included in fnorm]
!
              f(l1:l2,m,n,jf) = f(l1:l2,m,n,jf) + fnorm*capitalH(:,j)
            endif
          enddo
        enddo
      enddo

    endsubroutine forcing_chandra_kendall
!***********************************************************************
    subroutine forcing_spherical_radial(f)
!
!  Add radial non-helical forcing function in spherical polar coordinate system.
!  16-Aug-2023/MR+Aba: adapted from chandra_kendall forcing.
!
      use EquationOfState, only: cs20
      use General, only: random_number_wrapper
      use Sub
      use Mpicomm
!
      real, dimension (mx,my,mz,mfarray) :: f
      !
      real, dimension(3) :: ee
      real, dimension(nx,3) :: capitalT,capitalS,psi
      real, dimension(nx,3,3) :: psi_ij,Tij
      integer :: emm,l,j,jf,Legendrel,lmindex,aindex,ell_rand,raindex
      real :: a_ell,anum,adenom,jlm_ff,ylm_ff,rphase1,fnorm,alphar,Balpha,psilm,RYlm,IYlm
      real :: rz,rindex,ralpha,rphase2,thphase,sthphase,cthphase,costhprime,phprime,alphap,&
              ralpha2,ellSR,fmaxloc,fmax,gindex,crphase1,srphase1
      real, dimension(mx) :: Z_psi
      real, dimension(nx) :: Smod
!
      if (lhelical_test) then
        if (icklist==nlist_ck) &
          call fatal_error("forcing_spherical_radial","SR testing: no more values in list")
        icklist=icklist+1
        lmindex=icklist
      else
        call random_number_wrapper(rindex,CHANNEL=channel_force)
        lmindex=nint(rindex*(nlist_ck-1))+1
      endif
!
      emm = cklist(lmindex,1)
      Legendrel = cklist(lmindex,2)
!
      call random_number_wrapper(ralpha,CHANNEL=channel_force)
      aindex=nint(ralpha*(ncol-3))
      Balpha = cklist(lmindex,3+aindex)
!
      if (SR_alpha_shift) then
        call random_number_wrapper(ralpha2,CHANNEL=channel_force)
        alphap=alpha_shift*2*ralpha2-alpha_shift
        Balpha=Balpha+alphap
      endif
! 
      if (SR_ell_shift) then
        call random_number_wrapper(ellSR,CHANNEL=channel_force)
        ell_rand=nint(ellSR*2*ell_shift)-ell_shift
        Legendrel=Legendrel+ell_rand 
      endif
!
! Calculate the "potential" for the non-helical forcing according to Chandrasekhar & Kendall.
! Construct Z_psi(r)
!
      call sp_bessely_l(anum,Legendrel,Balpha*xyz0(1))
      call sp_besselj_l(adenom,Legendrel,Balpha*xyz0(1))
      a_ell = -anum/adenom
!
      do l=1,mx
        alphar=Balpha*x(l)
        call sp_besselj_l(jlm_ff,Legendrel,alphar)
        call sp_bessely_l(ylm_ff,Legendrel,alphar)
        Z_psi(l) = (a_ell*jlm_ff+ylm_ff)
      enddo

      if (lisotropize_SR) then

        call random_number_wrapper(rphase1,CHANNEL=channel_force)
        rphase1=rphase1*2*pi

        if (lfastSR) then
          call random_number_wrapper(gindex,CHANNEL=channel_force)
          raindex=nint(gindex*(n_axisrot_angles-1))+1
          crphase1=cos(rphase1)
          srphase1=sin(rphase1)
          do n=1,mz
            do m=1,my
              psilm = RYlm_rot(m,n,lmindex,raindex)*crphase1 - &
                      IYlm_rot(m,n,lmindex,raindex)*srphase1
              psif(:,m,n) = Z_psi*psilm
            enddo
          enddo
        else
!
! Randomize axis of spherical harmonic by rotation about the x direction.
!
          call random_number_wrapper(thphase,CHANNEL=channel_force)
          thphase = pi*thphase
!
          do n=1,mz
            do m=1,my
              if (lisotropize_SR) then
                sthphase=sin(thphase); cthphase=cos(thphase)
                costhprime=sthphase*cosph(n)*sinth(m)+cthphase*costh(m)
                if (abs(costhprime)==1.) then
                  phprime=0.
                else
                  phprime=acos((cthphase*cosph(n)*sinth(m)-sthphase*costh(m))/sqrt(1.-costhprime**2))
                endif
              else
                costhprime=costh(m); phprime=z(n)
              endif
              call sp_harm_real_costh(RYlm,Legendrel,emm,costhprime,phprime)
              call sp_harm_imag_costh(IYlm,Legendrel,emm,costhprime,phprime)
              psilm = RYlm*cos(rphase1)-IYlm*sin(rphase1)
              psif(:,m,n) = Z_psi*psilm
            enddo
          enddo

        endif
!
! Calculate the force from the potential and store it as an auxiliary (except in test branches).
!
        ee = (/1.,0.,0./)
        fnorm = cs20*sqrt(1./(cs0*Balpha))*sqrt(dt)
        fmaxloc=0.
        do n=n1,n2
          do m=m1,m2
            psi(:,1) = psif(l1:l2,m,n)*ee(1)
            call gij_psi(psif,ee,psi_ij)
            call curl_mn(psi_ij,capitalT,psi)
            call gij_psi_etc(psif,ee,psi,psi_ij,Tij)
            call curl_mn(Tij,capitalS,capitalT)
            capitalS(:,1) = (1./Balpha)*capitalS(:,1)
            if (r_ff /= 0.) capitalS(:,1)=profx_ampl*capitalS(:,1)
            if (lhelical_test) then
              if (lwrite_psi) then
                f(l1:l2,m,n,iux) = psif(l1:l2,m,n)
              else
                f(l1:l2,m,n,iux) = fpre*fnorm*capitalS(:,1)
              endif
            else
              fmaxloc=max(fmaxloc,maxval(abs(capitalS(:,1))))
              f(l1:l2,m,n,iff_aux) = fpre*capitalS(:,1)
            endif
          enddo
        enddo
  
        if (.not.lhelical_test) then
          call mpiallreduce_max(fmaxloc,fmax)
!
! Stochastic euler scheme of integration [sqrt(dt) is already included in fnorm].
!
          f(l1:l2,m1:m2,n1:n2,iux) = (1./fmax)*f(l1:l2,m1:m2,n1:n2,iff_aux)
        endif
      endif

    endsubroutine forcing_spherical_radial
!***********************************************************************
    subroutine forcing_GP(f)
!
!  Add Galloway-Proctor forcing function.
!
!  24-jul-06/axel: coded
!
      use Mpicomm, only: mpireduce_sum
      use Sub
      use DensityMethods, only: getrho
      use Diagnostics, only: sum_mn_name
!
      real, dimension (nx) :: ruf,rho
      real, dimension (nx,3) :: variable_rhs,forcing_rhs,force_all
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx) :: cosx,sinx
      real :: cost,sint,cosym,sinym
      integer :: j,jf
      real :: fact
!
!  The default for omega_ff is now (12-jan-2018) changed from 1 to 0.
!
      if (ip<=6) print*,'forcing_GP: t=',t
      cost=cos(omega_ff*t)
      sint=sin(omega_ff*t)
!
!  Normalize ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
      fact=sqrt(1.5)*force*sqrt(dt)
!
!  loop the two cases separately, so we don't check for r_ff during
!  each loop cycle which could inhibit (pseudo-)vectorisation
!  calculate energy input from forcing; must use lout (not ldiagnos)
!
      lfirstpoint = .false.
      do m=m1,m2
        cosx=cos(k1_ff*x+cost)
        sinx=sin(k1_ff*x+cost)
        cosym=cos(k1_ff*y(m)+sint)
        sinym=sin(k1_ff*y(m)+sint)
        forcing_rhs(:,1)=-fact*sinym
        forcing_rhs(:,2)=-fact*cosx(l1:l2)
        forcing_rhs(:,3)=+fact*(sinx(l1:l2)+cosym)
        do n=n1,n2
          variable_rhs=f(l1:l2,m,n,iffx:iffz)
          f(l1:l2,m,n,iffx:iffz)=f(l1:l2,m,n,iffx:iffz)+forcing_rhs
          if (lout) then
            if (idiag_rufm/=0) then
              call getrho(f(:,m,n,ilnrho),rho)
              call multsv_mn(rho/dt,forcing_rhs,force_all)
              call dot_mn(variable_rhs,force_all,ruf)
              call sum_mn_name(ruf,idiag_rufm)
            endif
          endif
          lfirstpoint = .false.
        enddo
      enddo
!
! For printouts, rufm needs to be communicated to other processors
!
      if (lout) then
        if (idiag_rufm/=0) call mpireduce_sum(fname(idiag_rufm)/nwgrid,fname(idiag_rufm))
      endif
!
      if (ip<=9) print*,'forcing_GP: forcing OK'
!
    endsubroutine forcing_GP
!***********************************************************************
    subroutine forcing_GP92(f)
!
!  Add Galloway-Proctor (1992) forcing function.
!
!  23-mar-20/axel: adapted from forcing_GP92
!
      use Mpicomm, only: mpireduce_sum
      use Sub
      use DensityMethods, only: getrho
      use Diagnostics, only: sum_mn_name
!
      real, dimension (nx) :: ruf,rho
      real, dimension (nx,3) :: variable_rhs,forcing_rhs,force_all
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx) :: cosx,sinx
      real :: cost,sint,cosym,sinym
      integer :: j,jf
      real :: fact
!
!  The default for omega_ff is now (12-jan-2018) changed from 1 to 0.
!
      if (ip<=6) print*,'forcing_GP: t=',t
      cost=cos(omega_ff*t)
      sint=sin(omega_ff*t)
!
!  Normalize ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
      fact=sqrt(1.5)*force*sqrt(dt)
!
!  loop the two cases separately, so we don't check for r_ff during
!  each loop cycle which could inhibit (pseudo-)vectorisation
!  calculate energy input from forcing; must use lout (not ldiagnos)
!
      cosx=cos(k1_ff*x+cost)
      sinx=sin(k1_ff*x+cost)
      forcing_rhs(:,2)=-fact*cosx(l1:l2)
      lfirstpoint = .true.
      do m=m1,m2
        cosym=cos(k1_ff*y(m)+sint)
        sinym=sin(k1_ff*y(m)+sint)
        forcing_rhs(:,1)=-fact*sinym
        forcing_rhs(:,3)=+fact*(sinx(l1:l2)+cosym)
        do n=n1,n2
          variable_rhs=f(l1:l2,m,n,iffx:iffz)
          f(l1:l2,m,n,iffx:iffz)=f(l1:l2,m,n,iffx:iffz)+forcing_rhs
          if (lout) then
            if (idiag_rufm/=0) then
              call getrho(f(:,m,n,ilnrho),rho)
              call multsv_mn(rho/dt,forcing_rhs,force_all)
              call dot_mn(variable_rhs,force_all,ruf)
              call sum_mn_name(ruf,idiag_rufm)
            endif
          endif
          lfirstpoint = .false.
        enddo
      enddo
!
! For printouts, rufm needs to be communicated to other processors
!
      if (lout) then
        if (idiag_rufm/=0) call mpireduce_sum(fname(idiag_rufm)/nwgrid,fname(idiag_rufm))
      endif
!
      if (ip<=9) print*,'forcing_GP: forcing OK'
!
    endsubroutine forcing_GP92
!***********************************************************************
    subroutine forcing_TG(f)
!
!  Add Taylor-Green forcing function.
!
!   9-oct-04/axel: coded
!
      use Mpicomm, only: mpireduce_sum
      use Sub
      use DensityMethods, only: getrho
      use Diagnostics, only: sum_mn_name
!
      real, dimension (nx) :: ruf,rho
      real, dimension (nx,3) :: variable_rhs,forcing_rhs,force_all
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx), save :: sinx,cosx
      real, dimension (my), save :: siny,cosy
      real, dimension (mz), save :: cosz
      logical, save :: lfirst_call=.true.
      integer :: j,jf
      real :: fact
!
      if (lfirst_call) then
        if (lroot) print*,'forcing_TG: calculate sinx,cosx,siny,cosy,cosz'
        sinx=sin(k1_ff*x)
        cosx=cos(k1_ff*x)
        siny=sin(k1_ff*y)
        cosy=cos(k1_ff*y)
        cosz=cos(k1_ff*z)
        lfirst_call=.false.
      endif
!
      if (ip<=6) print*,'forcing_TG: dt, lfirst_call=',dt,lfirst_call
!
!  Normalize ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
      fact=2*force*sqrt(dt)
!
!  loop the two cases separately, so we don't check for r_ff during
!  each loop cycle which could inhibit (pseudo-)vectorisation
!  calculate energy input from forcing; must use lout (not ldiagnos)
!
      forcing_rhs(:,3)=0.
      lfirstpoint = .true.
      do n=n1,n2
        do m=m1,m2
          variable_rhs=f(l1:l2,m,n,iffx:iffz)
          forcing_rhs(:,1)=+fact*sinx(l1:l2)*cosy(m)*cosz(n)
          forcing_rhs(:,2)=-fact*cosx(l1:l2)*siny(m)*cosz(n)
          do j=1,3
            if (lactive_dimension(j)) then
              jf=j+ifff-1
              f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)
            endif
          enddo
          if (lout) then
            if (idiag_rufm/=0) then
              call getrho(f(:,m,n,ilnrho),rho)
              call multsv_mn(rho/dt,forcing_rhs,force_all)
              call dot_mn(variable_rhs,force_all,ruf)
              call sum_mn_name(ruf,idiag_rufm)
            endif
          endif
          lfirstpoint = .false.
        enddo
      enddo
!
! For printouts, rufm needs to be communicated to other processors
!
      if (lout) then
        if (idiag_rufm/=0) call mpireduce_sum(fname(idiag_rufm)/nwgrid,fname(idiag_rufm))
      endif
!
      if (ip<=9) print*,'forcing_TG: forcing OK'
!
    endsubroutine forcing_TG
!***********************************************************************
    subroutine forcing_ABC(f)
!
!  Added ABC forcing function
!
!  17-jul-06/axel: coded
!
      use Diagnostics
      use DensityMethods, only: getrho
      use Mpicomm, only: mpireduce_sum
      use Sub
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      real, dimension (nx) :: ruf,rho
      real, dimension (nx,3) :: variable_rhs,forcing_rhs,force_all,bb,fxb
      real, dimension (mx), save :: sinx,cosx
      real, dimension (my), save :: siny,cosy
      real, dimension (mz), save :: sinz,cosz
      logical, save :: lfirst_call=.true.
      integer :: j,jf
      real :: fact
!
!  at the first step, the sin and cos functions are calculated for all
!  x,y,z points and are then saved and used for all subsequent steps
!  and pencils
!
      if (ip<=6) print*,'forcing_ABC: lfirst_call=',lfirst_call
      if (lfirst_call) then
        if (lroot) print*,'forcing_ABC: calculate sinx,cosx,siny,cosy,sinz,cosz'
        sinx=sin(k1_ff*x); cosx=cos(k1_ff*x)
        siny=sin(k1_ff*y); cosy=cos(k1_ff*y)
        sinz=sin(k1_ff*z); cosz=cos(k1_ff*z)
        lfirst_call=.false.
      endif
      if (ip<=6) print*,'forcing_ABC: dt, lfirst_call=',dt,lfirst_call
!
!  Normalize ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
      fact=2*force*sqrt(dt)
!
!  loop the two cases separately, so we don't check for r_ff during
!  each loop cycle which could inhibit (pseudo-)vectorisation
!  calculate energy input from forcing; must use lout (not ldiagnos)
!
      lfirstpoint = .true.
      do n=n1,n2
        forcing_rhs(:,2)=fact*(sinx(l1:l2)+cosz(n))
        do m=m1,m2
          variable_rhs=f(l1:l2,m,n,iffx:iffz)
          forcing_rhs(:,1)=fact*(sinz(n    )+cosy(m)    )
          forcing_rhs(:,3)=fact*(siny(m    )+cosx(l1:l2))
          f(l1:l2,m,n,iffx:iffz)=f(l1:l2,m,n,iffx:iffz)+forcing_rhs
          if (lout) then
            if (idiag_rufm/=0) then
              call getrho(f(:,m,n,ilnrho),rho)
              call multsv_mn(rho/dt,forcing_rhs,force_all)
              call dot_mn(variable_rhs,force_all,ruf)
              call sum_mn_name(ruf,idiag_rufm)
            endif
            if (lmagnetic) then
!             if (idiag_fxbxm/=0.or.idiag_fxbym/=0.or.idiag_fxbzm/=0) then
!               call curl(f,iaa,bb)
!               call cross(forcing_rhs,bb,fxb)
!               call sum_mn_name(fxb(:,1),idiag_fxbxm)
!               call sum_mn_name(fxb(:,2),idiag_fxbym)
!               call sum_mn_name(fxb(:,3),idiag_fxbzm)
!             endif
            endif
          endif
          lfirstpoint = .false.
        enddo
      enddo
!
! For printouts global, summation is needed as not done by finalize_diagnostics.
!
      if (lout) then
        if (idiag_rufm/=0) call mpireduce_sum(fname(idiag_rufm)/nwgrid,fname(idiag_rufm))
!       if (lmagnetic) then
!         if (idiag_fxbxm/=0) call mpireduce_sum(fname(idiag_fxbxm)/nwgrid,fname(idiag_fxbxm))
!         if (idiag_fxbym/=0) call mpireduce_sum(fname(idiag_fxbym)/nwgrid,fname(idiag_fxbym))
!         if (idiag_fxbzm/=0) call mpireduce_sum(fname(idiag_fxbzm)/nwgrid,fname(idiag_fxbzm))
!       endif
      endif
!
      if (ip<=9) print*,'forcing_ABC: forcing OK'
!
    endsubroutine forcing_ABC
!***********************************************************************
    subroutine forcing_mhd_mode(f)
!
!  Added ABC forcing function
!
!  17-jul-06/axel: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      real :: fact, om
      integer :: m,n
!
      fact=force*sqrt(dt)
      do n=n1,n2
        do m=m1,m2
          !f(:,m,n,iuy)=f(:,m,n,iuy)+fact*sin(k1_ff*x)
          !f(:,m,n,iAz)=f(:,m,n,iAz)+fact*sin(k1_ff*x)
!
!  fast mode
!
          om=k1_ff*sqrt(b0_mode**2+1.)
          f(:,m,n,iuy)=f(:,m,n,iuy)+fact*b0_mode*sin(k1_ff*x)/om
          f(:,m,n,iAz)=f(:,m,n,iAz)+fact*b0_mode*sin(k1_ff*x)*(k1_ff/om)**2
        enddo
      enddo
      !
      if (ip<=9) print*,'forcing_mhd: forcing OK'
!
    endsubroutine forcing_mhd_mode
!***********************************************************************
    subroutine forcing_nocos(f)
!
!  Add no-cosine forcing function.
!
!  27-oct-04/axel: coded
!
      use DensityMethods, only: getrho
      use Diagnostics, only: sum_mn_name
      use Mpicomm, only: mpireduce_sum
      use Sub
!
      real, dimension (nx) :: ruf,rho
      real, dimension (nx,3) :: variable_rhs,forcing_rhs,force_all
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx), save :: sinx
      real, dimension (my), save :: siny
      real, dimension (mz), save :: sinz
      logical, save :: lfirst_call=.true.
      integer :: j,jf
      real :: fact
!
      if (lfirst_call) then
        if (lroot) print*,'forcing_nocos: calculate sinx,siny,sinz'
        sinx=sin(k1_ff*x)
        siny=sin(k1_ff*y)
        sinz=sin(k1_ff*z)
        lfirst_call=.false.
      endif
!
      if (ip<=6) print*,'forcing_hel_noshear: dt, lfirst_call=',dt,lfirst_call
!
!  Normalize ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
      fact=force*sqrt(dt)
!
!  loop the two cases separately, so we don't check for r_ff during
!  each loop cycle which could inhibit (pseudo-)vectorisation
!  calculate energy input from forcing; must use lout (not ldiagnos)
!
      forcing_rhs(:,2)=fact*sinx(l1:l2)
      lfirstpoint = .true.
      do n=n1,n2
        forcing_rhs(:,1)=fact*sinz(n)
        do m=m1,m2
          variable_rhs=f(l1:l2,m,n,iffx:iffz)
          forcing_rhs(:,3)=fact*siny(m)
          do j=1,3
            if (lactive_dimension(j)) then
              jf=j+ifff-1
              f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+forcing_rhs(:,j)
            endif
          enddo
          if (lout) then
            if (idiag_rufm/=0) then
              call getrho(f(:,m,n,ilnrho),rho)
              call multsv_mn(rho/dt,forcing_rhs,force_all)
              call dot_mn(variable_rhs,force_all,ruf)
              call sum_mn_name(ruf,idiag_rufm)
            endif
          endif
          lfirstpoint = .false.
        enddo
      enddo
!
! For printouts, rufm needs to be communicated to other processors.
!
      if (lout) then
        if (idiag_rufm/=0) call mpireduce_sum(fname(idiag_rufm)/nwgrid,fname(idiag_rufm))
      endif
!
      if (ip<=9) print*,'forcing_nocos: forcing OK'
!
    endsubroutine forcing_nocos
!***********************************************************************
    subroutine forcing_gaussianpot(f,force_ampl)
!
!  gradient of gaussians as forcing function
!
!  19-dec-05/tony: coded, adapted from forcing_nocos
!  14-jul-10/axel: in less then 3-D, project forcing to computational domain
!
      use DensityMethods, only: getrho
      use Diagnostics, only: sum_mn_name
      use General, only: random_number_wrapper
      use Mpicomm, only: mpireduce_sum
      use Sub, only: multsv_mn,dot_mn,del2v_etc
!
      real, dimension (mx,my,mz,mfarray) :: f
      real :: force_ampl, force_tmp
!
      real, dimension (3) :: fran
      real, dimension (nx) :: radius2,gaussian,gaussian_fact,ruf,rho,rho1,qf
      real, dimension (nx,3) :: variable_rhs,delta,curlo,forcing_rhs,force_all
      integer :: j, jf, ilocation
      real :: fact,width_ff21
!
!  check length of time step
!
      if (ip<=6) print*,'forcing_gaussianpot: dt=',dt
!
!  generate random numbers
!
      if (t>tsforce) then
        if (lrandom_location) then
          call random_number_wrapper(fran,CHANNEL=channel_force)
          location=fran*Lxyz+xyz0
        else
          location=location_fixed(:,1)
        endif
!
!  It turns out that in the presence of shear, and even for weak shear,
!  vorticity is being produced. In order to check whether the shearing
!  periodic boundaries are to blame, we can cut the y extent of forcing
!  locations by half.
!
        if (lforce_cuty) location(2)=location(2)*.5
!
!  reset location(i) to x, y, or z, so forcing does still occur,
!  but it is restricted to the plane or line that is being solved.
!
        if (.not.lactive_dimension(1)) location(1)=x(1)
        if (.not.lactive_dimension(2)) location(2)=y(1)
        if (.not.lactive_dimension(3)) location(3)=z(1)
!
!  write location to file; this is off by default.
!
        if (lroot .and. lwrite_gausspot_to_file) then
          open(1,file=trim(datadir)//'/gaussian_pot_forcing.dat',status='unknown',position='append')
          write(1,'(4f14.7)') t, location
          close (1)
        endif
!
!  set next forcing time
!
        tsforce=t+dtforce
        if (ip<=6) print*,'forcing_gaussianpot: location=',location
      endif   ! if (t>tsforce)
!
!  Set forcing amplitude (same value for each location by default).
!  For iforce_tprofile=='sin^2', we have a temporal sine profile.
!  By default, iforce_tprofile='nothing'.
!
      if (iforce_tprofile=='nothing') then
        force_tmp=force_ampl
      elseif (iforce_tprofile=='sin^2') then
        force_tmp=force_ampl*sin(pi*(tsforce-t)/dtforce)**2
      else
        call fatal_error('forcing_gaussianpot','no such iforce_tprofile: '//trim(iforce_tprofile))
      endif
!
!  Possibility of outputting data at each time step (for testing).
!
      if (lroot .and. lwrite_gausspot_to_file_always) then
        open(1,file=trim(datadir)//'/gaussian_pot_forcing_always.dat', &
            status='unknown',position='append')
        write(1,'(5f14.7)') t, location, force_tmp
        close (1)
      endif
!
!  Let explosion last dtforce_duration or, by default, until next explosion.
!  By default, dtforce_duration=1 (this is a special value, different from zero),
!  but when t>dtforce, forcing is off again, regardless of dtforce_duration.
!
      if ( (dtforce_duration<0.0) .or. (t-(tsforce-dtforce))<=dtforce_duration ) then
!
!  Normalize ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!
!  We multiply the forcing term by dt and add to the right-hand side
!  of the momentum equation for an Euler step, but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference.
!  When dtforce is finite, take dtforce+.5*dt.
!  The 1/2 factor takes care of round-off errors.
!  Also define width_ff21 = 1/width^2.
!  The factor 2 in fact takes care of the 2 in 2*xi/R^2.
!
        width_ff21=1./width_ff**2
        fact=2.*width_ff21*force_tmp*dt*sqrt(cs0*width_ff/max(dtforce+.5*dt,dt))
!
!  Loop the two cases separately, so we don't check for r_ff during
!  each loop cycle which could inhibit (pseudo-)vectorization.
!  Calculate energy input from forcing; must use lout (not ldiagnos).
!
!  loop over all pencils
!
        lfirstpoint = .true.
        do n=n1,n2
          do m=m1,m2
!
!  loop over multiple locations
!
            do ilocation=1,nlocation
!
!  Obtain distance (=delta) to center of blob
!
              if (ilocation==1) then
                delta(:,1)=x(l1:l2)-location(1)
                delta(:,2)=y(m)-location(2)
                delta(:,3)=z(n)-location(3)
              else
                delta(:,1)=x(l1:l2)-location2(1)
                delta(:,2)=y(m)-location2(2)
                delta(:,3)=z(n)-location2(3)
              endif
              do j=1,3
                if (lperi(j)) then
                  if (lforce_peri) then
                    if (j==2) then
                      delta(:,2) = 2*atan(tan(.5*(delta(:,2) + 2.*deltay*atan(1000.*tan(.25* &
                                   (pi+x(l1:l2)-location(1)))))))
                    else
                      delta(:,j) = 2*atan(tan(.5*delta(:,j)))
                    endif
                  else
                    where (delta(:,j) >  Lxyz(j)/2.) delta(:,j)=delta(:,j)-Lxyz(j)
                    where (delta(:,j) < -Lxyz(j)/2.) delta(:,j)=delta(:,j)+Lxyz(j)
                  endif
                endif
                if (.not.lactive_dimension(j)) delta(:,j)=0.
              enddo
!
!  compute gaussian blob and set to forcing function
!
              radius2=delta(:,1)**2+delta(:,2)**2+delta(:,3)**2
              gaussian=exp(-radius2*width_ff21)
              if (lmomentum_ff) then
                if (ldensity_nolog) then
                  rho1=1./f(l1:l2,m,n,irho)
                else
                  rho1=exp(-f(l1:l2,m,n,ilnrho))
                endif
                gaussian_fact=gaussian*fact*rho1
              else
                gaussian_fact=gaussian*fact
              endif
!
!  Forced variable before adding time-integrated (Euler's method) force.
!
              variable_rhs=f(l1:l2,m,n,iffx:iffz)
              if (iphiuu==0) then
                do j=1,3
                  if (lactive_dimension(j)) then
                    jf=j+ifff-1
                    if (iforce_profile=='nothing') then
                      forcing_rhs(:,j)=gaussian_fact*delta(:,j)
                    else
                      forcing_rhs(:,j)=gaussian_fact*delta(:,j)*profx_ampl*profy_ampl(m)*profz_ampl(n)
                    endif
                    f(l1:l2,m,n,jf) = f(l1:l2,m,n,jf)+forcing_rhs(:,j)
                  endif
                enddo
              else
!
!  add possibility of "potential forcing" (but only if iphiuu/=0, i.e., hydro_potential in use)
!
                if (iforce_profile=='nothing') then
                  forcing_rhs(:,1)=gaussian
                else
                  forcing_rhs(:,1)=gaussian*profx_ampl*profy_ampl(m)*profz_ampl(n)
                endif
                forcing_rhs(:,2:3)=0.
                f(l1:l2,m,n,iphiuu)=f(l1:l2,m,n,iphiuu)+forcing_rhs(:,1)
              endif
!
!  test
!
!--           if (icc/=0) f(l1:l2,m,n,icc)=f(l1:l2,m,n,icc)+gaussian
!
!  diagnostics (currently without this possibility of a modulation applied above)
!
              if (lout) then
                if (idiag_rufm/=0) then
                  call getrho(f(:,m,n,ilnrho),rho)
                  call multsv_mn(rho/dt,forcing_rhs,force_all)
                  call dot_mn(variable_rhs,force_all,ruf)
                  call sum_mn_name(ruf,idiag_rufm)
                endif
                if (idiag_qfm/=0) then
                  call del2v_etc(f,iuu,curlcurl=curlo)
                  call dot_mn(curlo,forcing_rhs,qf)
                  call sum_mn_name(qf/dt,idiag_qfm)
                endif
              endif
            enddo
            lfirstpoint = .false.
          enddo
        enddo
      endif    ! if ( (dtforce_duration<0.0) etc.
!
!  For printouts, rufm needs to be communicated to other processors.
!
      if (lout) then
        if (idiag_rufm/=0) call mpireduce_sum(fname(idiag_rufm)/nwgrid,fname(idiag_rufm))
        if (idiag_qfm/=0) call mpireduce_sum(fname(idiag_qfm)/nwgrid,fname(idiag_qfm))
      endif
!
      if (ip<=9) print*,'forcing_gaussianpot: forcing OK'
!
    endsubroutine forcing_gaussianpot
!***********************************************************************
    subroutine forcing_hillrain(f,force_ampl)
!
!  gradient of gaussians as forcing function
!
!  29-sep-15/axel: adapted from forcing_gaussianpot
!
      use DensityMethods, only: getrho
      use Diagnostics, only: sum_mn_name
      use General, only: random_number_wrapper
      use Mpicomm, only: mpireduce_sum
      use Sub
!
      real, dimension (mx,my,mz,mfarray) :: f
      real :: force_ampl, force_tmp
!
      real, dimension (3) :: fran
      real, dimension (nx) :: r, r2, r3, r5, pom2, ruf, rho
      real, dimension (nx,3) :: variable_rhs, force_all, delta
      integer :: j,jf
      real :: fact
      real :: a_hill, a2_hill, a3_hill
!
!  check length of time step
!
      if (ip<=6) print*,'forcing_hillrain: dt=',dt
!
!  generate random numbers
!
      if (t>tsforce) then
        if (lrandom_location) then
          call random_number_wrapper(fran,CHANNEL=channel_force)
          location=fran*Lxyz+xyz0
        else
          location=location_fixed(:,1)
        endif
!
!  It turns out that in the presence of shear, and even for weak shear,
!  vorticitity is being produced. In order to check whether the shearing
!  periodic boundaries are to blame, we can cut the y extent of forcing
!  locations by half.
!
        if (lforce_cuty) location(2)=location(2)*.5
!
!  reset location(i) to x or y, and keep location(3)=0 fixed
!
        if (.not.lactive_dimension(1)) location(1)=x(1)
        if (.not.lactive_dimension(2)) location(2)=y(1)
        location(3)=0.
!
!  write location to file
!
        if (lroot .and. lwrite_gausspot_to_file) then
          open(1,file=trim(datadir)//'/hillrain_forcing.dat',status='unknown',position='append')
          write(1,'(4f14.7)') t, location
          close (1)
        endif
!
!  set next forcing time
!
        tsforce=t+dtforce
        if (ip<=6) print*,'forcing_hillrain: location=',location
      endif
!
!  Set forcing amplitude (same value for each location by default)
!  We multiply the forcing term by dt and add to the right-hand side
!  of the momentum equation for an Euler step, but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference.
!  When dtforce is finite, take dtforce+.5*dt.
!  The 1/2 factor takes care of round-off errors.
!
!
      if (iforce_tprofile=='nothing') then
        force_tmp=force_ampl
        fact=force_tmp*dt*sqrt(cs0*radius_ff/max(dtforce+.5*dt,dt))
      elseif (iforce_tprofile=='sin^2') then
        force_tmp=force_ampl*sin(pi*(tsforce-t)/dtforce)**2
        fact=force_tmp*dt*sqrt(cs0*radius_ff/max(dtforce+.5*dt,dt))
      elseif (iforce_tprofile=='delta') then
        if (tsforce-t <= dt) then
          force_tmp=force_ampl
        else
          force_tmp=0.
        endif
        fact=force_tmp
      else
        call fatal_error('forcing_hillrain','no such iforce_tprofile: '//trim(iforce_tprofile))
      endif
!
!  Let explosion last dtforce_duration or, by default, until next explosion.
!
      if ( force_tmp /= 0. .and. ((dtforce_duration<0.0) .or. &
           (t-(tsforce-dtforce))<=dtforce_duration) ) then
!
!  loop the two cases separately, so we don't check for r_ff during
!  each loop cycle which could inhibit (pseudo-)vectorisation
!  calculate energy input from forcing; must use lout (not ldiagnos)
!
!  loop over all pencils
!
        lfirstpoint = .true.
        do n=n1,n2; do m=m1,m2
!
!  Obtain distance to center of blob
!
          delta(:,1)=x(l1:l2)-location(1)
          delta(:,2)=y(m)-location(2)
          delta(:,3)=z(n)-location(3)
          do j=1,3
            if (lperi(j)) then
              if (lforce_peri) then
                if (j==2) then
                  delta(:,2) = 2*atan(tan(.5*(delta(:,2)+2.*deltay*atan(1000.*tan(.25* &
                               (pi+x(l1:l2)-location(1)))))))
                else
                  delta(:,j) = 2*atan(tan(.5*delta(:,j)))
                endif
              else
                where (delta(:,j) >  Lxyz(j)/2.) delta(:,j)=delta(:,j)-Lxyz(j)
                where (delta(:,j) < -Lxyz(j)/2.) delta(:,j)=delta(:,j)+Lxyz(j)
              endif
            endif
            if (.not.lactive_dimension(j)) delta(:,j)=0.
          enddo
!
!  compute factors for Hill vortex
!
          r2=delta(:,1)**2+delta(:,2)**2+delta(:,3)**2
          pom2=delta(:,1)**2+delta(:,2)**2
          r=sqrt(r2)
          r3=r2*r
          r5=r2*r3
          a_hill=radius_ff
          a2_hill=a_hill**2
          a3_hill=a_hill*a2_hill
!
!  compute Hill vortex
!
          where (r<=a_hill)
            variable_rhs(:,1)=-1.5*delta(:,1)*delta(:,3)/a2_hill
            variable_rhs(:,2)=-1.5*delta(:,2)*delta(:,3)/a2_hill
            variable_rhs(:,3)=-2.5+1.5*(pom2+r2)/a2_hill
          elsewhere
            variable_rhs(:,1)=-1.5*delta(:,1)*delta(:,3)*a3_hill/r5
            variable_rhs(:,2)=-1.5*delta(:,2)*delta(:,3)*a3_hill/r5
            variable_rhs(:,3)=-a3_hill/r3+1.5*pom2*a3_hill/r5
          endwhere
!
!  add to the relevant variable
!
          do j=1,3
            if (lactive_dimension(j)) then
              jf=j+ifff-1
              f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+fact*variable_rhs(:,j)
            endif
          enddo
!
!  passive scalar as a possible test
!
!--         if (icc/=0) f(l1:l2,m,n,icc)=f(l1:l2,m,n,icc)+gaussian
!
!  diagnostics
!
          if (lout) then
            if (idiag_rufm/=0) then
              call getrho(f(:,m,n,ilnrho),rho)
              call multsv_mn(rho/dt,f(l1:l2,m,n,iux:iuz),force_all)
              call dot_mn(variable_rhs,force_all,ruf)
              call sum_mn_name(ruf,idiag_rufm)
            endif
          endif
          lfirstpoint = .false.
        enddo; enddo
!
!  For printouts, rufm needs to be communicated to other processors.
!
        if (lout) then
          if (idiag_rufm/=0) call mpireduce_sum(fname(idiag_rufm)/nwgrid,fname(idiag_rufm))
        endif
      endif
!
      if (ip<=9) print*,'forcing_hillrain: forcing OK'
!
    endsubroutine forcing_hillrain
!***********************************************************************
    subroutine forcing_white_noise(f)
!
!  gradient of gaussians as forcing function
!
!  19-dec-13/axel: added
!
      use DensityMethods, only: getrho
      use Diagnostics, only: sum_mn_name
      use General, only: random_number_wrapper
      use Mpicomm, only: mpireduce_sum
      use Sub
!
      real, dimension (mx,my,mz,mfarray) :: f
      real :: ampl
!
      real, dimension (nx) :: r,p,tmp,rho,ruf
      real, dimension (nx,3) :: force_all,variable_rhs,forcing_rhs
      integer :: j,jf
!
!  check length of time step
!
      if (ip<=6) print*,'forcing_white_noise: dt=',dt
!
!  extent
!
      if (.not.lactive_dimension(1)) location(1)=x(1)
      if (.not.lactive_dimension(2)) location(2)=y(1)
      if (.not.lactive_dimension(3)) location(3)=z(1)
      if (ip<=6) print*,'forcing_white_noise: location=',location
!
!  Normalize ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!
!  We multiply the forcing term by dt and add to the right-hand side
!  of the momentum equation for an Euler step, but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference.
!  When dtforce is finite, take dtforce+.5*dt.
!  The 1/2 factor takes care of round-off errors.
!  Also define width_ff21 = 1/width^2
!
      ampl=force*sqrt(dt*cs0)*cs0
!
!  loop the two cases separately, so we don't check for r_ff during
!  each loop cycle which could inhibit (pseudo-)vectorisation
!  calculate energy input from forcing; must use lout (not ldiagnos)
!
!  loop over all pencils
!
      lfirstpoint = .true.
      do n=n1,n2
        do m=m1,m2
          do j=1,3
            if (lactive_dimension(j)) then
              jf=j+ifff-1
              if (modulo(j-1,2)==0) then
                call random_number_wrapper(r,CHANNEL=channel_force)
                call random_number_wrapper(p,CHANNEL=channel_force)
                tmp=sqrt(-2*log(r))*sin(2*pi*p)
              else
                tmp=sqrt(-2*log(r))*cos(2*pi*p)
              endif
              f(l1:l2,m,n,jf)=f(l1:l2,m,n,jf)+ampl*tmp
              forcing_rhs(:,j)=ampl*tmp
            endif
          enddo
!
!  diagnostics
!
          if (lout) then
            if (idiag_rufm/=0) then
              variable_rhs=f(l1:l2,m,n,iffx:iffz)
              call getrho(f(:,m,n,ilnrho),rho)
              call multsv_mn(rho/dt,forcing_rhs,force_all)
              call dot_mn(variable_rhs,force_all,ruf)
              call sum_mn_name(ruf,idiag_rufm)
            endif
          endif
          lfirstpoint = .false.
        enddo
      enddo
!
!  For printouts, rufm needs to be communicated to other processors
!
      if (lout) then
        if (idiag_rufm/=0) call mpireduce_sum(fname(idiag_rufm)/nwgrid,fname(idiag_rufm))
      endif
!
      if (ip<=9) print*,'forcing_white_noise: forcing OK'
!
    endsubroutine forcing_white_noise
!***********************************************************************
    function calc_force_ampl(f,fx,fy,fz,coef) result(force_ampl)
!
!  calculates the coefficient for a forcing that satisfies
!  <rho*u*f> = constant.
!
!   7-sep-02/axel: coded
!
      use DensityMethods, only: getrho
      use Mpicomm
      use Sub
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,3) :: uu
      real, dimension (nx) :: rho,udotf
      complex, dimension (mx) :: fx
      complex, dimension (my) :: fy
      complex, dimension (mz) :: fz
      complex, dimension (3) :: coef
      real :: rho_uu_ff,force_ampl,fsum_tmp,fsum
      integer :: j,m,n
!
      rho_uu_ff=0.
      do n=n1,n2
        do m=m1,m2
          call getrho(f(:,m,n,ilnrho),rho)
          uu=f(l1:l2,m,n,iffx:iffz)
          udotf=0.
          do j=1,3
            udotf=udotf+uu(:,j)*real(coef(j)*fx(l1:l2)*fy(m)*fz(n))
          enddo
          rho_uu_ff=rho_uu_ff+sum(rho*udotf)
        enddo
      enddo
!
!  on different processors, this result needs to be communicated
!  to other processors
!
      fsum_tmp=rho_uu_ff
      call mpireduce_sum(fsum_tmp,fsum)
      if (lroot) rho_uu_ff=fsum/nwgrid
!      if (lroot) rho_uu_ff=rho_uu_ff/nwgrid
      call mpibcast_real(rho_uu_ff)
!
!  scale forcing function
!  but do this only when rho_uu_ff>0.; never allow it to change sign
!
        if (headt) print*,'calc_force_ampl: divide forcing function by rho_uu_ff=',rho_uu_ff
        !      force_ampl=work_ff/(.1+max(0.,rho_uu_ff))
        force_ampl=work_ff/rho_uu_ff
        if (force_ampl > max_force) force_ampl=max_force
        if (force_ampl < -max_force) force_ampl=-max_force
!
    endfunction calc_force_ampl
!***********************************************************************
    subroutine forcing_hel_noshear(f)
!
!  add helical forcing function, using a set of precomputed wavevectors
!
!  10-apr-00/axel: coded
!  06-dec-13/nishant: made kkx etc allocatable
!
      use General, only: random_number_wrapper
      use Sub
!
      real :: phase,ffnorm
      real, dimension (2) :: fran
      real, dimension (nx) :: radius,tmpx
!      real, dimension (mz) :: tmpz
      real, dimension (mx,my,mz,mfarray) :: f
      complex, dimension (mx) :: fx
      complex, dimension (my) :: fy
      complex, dimension (mz) :: fz
      complex, dimension (3) :: coef
      integer :: ik,j,jf,kx,ky,kz,kex,key,kez,kkex,kkey,kkez
      real :: k2,k,ex,ey,ez,kde,sig,fact
      real, dimension(3) :: e1,e2,ee,kk
      real :: norm,phi
!
!  generate random coefficients -1 < fran < 1
!  ff=force*Re(exp(i(kx+phase)))
!  |k_i| < akmax
!
      call random_number_wrapper(fran,CHANNEL=channel_force)
      phase=pi*(2*fran(1)-1.)
      ik=nk*.9999*fran(2)+1
      if (ip<=6) print*,'force_hel_noshear: ik,phase,kk=',ik,phase,kkx(ik),kky(ik),kkz(ik),dt
!
      kx=kkx(ik)
      ky=kky(ik)
      kz=kkz(ik)
      if (ip<=4) print*, 'force_hel_noshear: kx,ky,kz=',kx,ky,kz
!
      k2=float(kx**2+ky**2+kz**2)
      k=sqrt(k2)
!
! Find e-vector
!
      !
      ! Start with old method (not isotropic) for now.
      ! Pick e1 if kk not parallel to ee1. ee2 else.
      !
      if ((ky==0).and.(kz==0)) then
        ex=0; ey=1; ez=0
      else
        ex=1; ey=0; ez=0
      endif
      if (.not. old_forcing_evector) then
        !
        !  Isotropize ee in the plane perp. to kk by
        !  (1) constructing two basis vectors for the plane perpendicular
        !      to kk, and
        !  (2) choosing a random direction in that plane (angle phi)
        !  Need to do this in order for the forcing to be isotropic.
        !
        kk = (/kx, ky, kz/)
        ee = (/ex, ey, ez/)
        call cross(kk,ee,e1)
        call dot2(e1,norm); e1=e1/sqrt(norm) ! e1: unit vector perp. to kk
        call cross(kk,e1,e2)
        call dot2(e2,norm); e2=e2/sqrt(norm) ! e2: unit vector perp. to kk, e1
        call random_number_wrapper(phi,CHANNEL=channel_force); phi = phi*2*pi
        ee = cos(phi)*e1 + sin(phi)*e2
        ex=ee(1); ey=ee(2); ez=ee(3)
      endif
!
!  k.e
!
      call dot(kk,ee,kde)
!
!  k x e
!
      kex=ky*ez-kz*ey
      key=kz*ex-kx*ez
      kez=kx*ey-ky*ex
!
!  k x (k x e)
!
      kkex=ky*kez-kz*key
      kkey=kz*kex-kx*kez
      kkez=kx*key-ky*kex
!
!  ik x (k x e) + i*phase
!
!  Normalise ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!
!  This does already include the new sqrt(2) factor (missing in B01).
!  So, in order to reproduce the 0.1 factor mentioned in B01
!  we have to set force=0.07.
!
      ffnorm=sqrt(2.)*k*sqrt(k2-kde**2)/sqrt(kav*cs0**3)
      if (ip<=12) then
        print*,'force_hel_noshear: k,kde,ffnorm,kav,dt,cs0=',k,kde,ffnorm,kav,dt,cs0
        print*,'force_hel_noshear: k*sqrt(k2-kde**2)=',k*sqrt(k2-kde**2)
      endif
      !!(debug...) write(21,'(f10.4,3i3,f7.3)') t,kx,ky,kz,phase
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
      fact=force/ffnorm*sqrt(dt)
!
!  The wavevector is for the case where Lx=Ly=Lz=2pi. If that is not the
!  case one needs to scale by 2pi/Lx, etc.
!
      fx=exp(cmplx(0.,2*pi/Lx*kx*x+phase))*fact
      fy=exp(cmplx(0.,2*pi/Ly*ky*y))
      fz=exp(cmplx(0.,2*pi/Lz*kz*z))
!
!  possibly multiply forcing by z-profile
!
!-    if (height_ff/=0.) then
!-      if (lroot) print*,'forcing_hel_noshear: include z-profile'
!-      tmpz=(z/height_ff)**2
!-      fz=fz*exp(-tmpz**5/max(1.-tmpz,1e-5))
!-    endif
!
!  need to discuss with axel
!
!  possibly multiply forcing by sgn(z) and radial profile
!
!      if (r_ff/=0.) then
!        if (lroot) print*,'forcing_hel_noshear: applying sgn(z)*xi(r) profile'
!        !
!        ! only z-dependent part can be done here; radial stuff needs to go
!        ! into the loop
!        !
!        tmpz = tanh(z/width_ff)
!        fz = fz*tmpz
!      endif
!
      if (ip<=5) then
        print*,'force_hel_noshear: fx=',fx
        print*,'force_hel_noshear: fy=',fy
        print*,'force_hel_noshear: fz=',fz
      endif
!
!  prefactor
!
      sig=relhel
      coef(1)=cmplx(k*float(kex),sig*float(kkex))
      coef(2)=cmplx(k*float(key),sig*float(kkey))
      coef(3)=cmplx(k*float(kez),sig*float(kkez))
      if (ip<=5) print*,'force_hel_noshear: coef=',coef
!
! loop the two cases separately, so we don't check for r_ff during
! each loop cycle which could inhibit (pseudo-)vectorisation
!
      if (r_ff == 0) then       ! no radial profile
        do j=1,3
          jf=j+ifff-1
          do n=n1,n2
            do m=m1,m2
              f(l1:l2,m,n,jf) = f(l1:l2,m,n,jf)+real(coef(j)*fx(l1:l2)*fy(m)*fz(n))
            enddo
          enddo
        enddo
      else                      ! with radial profile
        do j=1,3
          jf=j+ifff-1
          do n=n1,n2
!---        sig = relhel*tmpz(n)
!AB: removed tmpz factor
              sig = relhel
            call fatal_error('forcing_hel_noshear','radial profile should be quenched')
            coef(1)=cmplx(k*float(kex),sig*float(kkex))
            coef(2)=cmplx(k*float(key),sig*float(kkey))
            coef(3)=cmplx(k*float(kez),sig*float(kkez))
            do m=m1,m2
              radius = sqrt(x(l1:l2)**2+y(m)**2+z(n)**2)
              tmpx = 0.5*(1.-tanh((radius-r_ff)/width_ff))
              f(l1:l2,m,n,jf) = f(l1:l2,m,n,jf) + real(coef(j)*tmpx*fx(l1:l2)*fy(m)*fz(n))
            enddo
          enddo
        enddo
      endif
!
      if (ip<=12) print*,'force_hel_noshear: forcing OK'
!
    endsubroutine forcing_hel_noshear
!***********************************************************************
    subroutine forcing_roberts(f)
!
!  add some artificial fountain flow
!  (to check for example small scale magnetic helicity loss)
!
!  30-may-02/axel: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: sxx,cxx
      real, dimension (mx) :: sx,cx
      real, dimension (my) :: sy,cy
      real, dimension (mz) :: sz,cz,tmpz,gz,gg,ss,gz1
      real :: kx,ky,kz,ffnorm,fac
!
!  identify ourselves
!
      if (headtt.or.ldebug) print*,'forcing_roberts: ENTER'
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
      kx=kfountain
      ky=kfountain
      kz=1.
!
      sx=sin(kx*x); cx=cos(kx*x)
      sy=sin(ky*y); cy=cos(ky*y)
      sz=sin(kz*z); cz=cos(kz*z)
!
!  abbreviation
!
      sxx=sx(l1:l2)
      cxx=cx(l1:l2)
!
!  g(z) and g'(z)
!  use z-profile to cut off
!
      if (height_ff/=0.) then
        tmpz=(z/height_ff)**2
        gz=sz*exp(-tmpz**5/max(1.-tmpz,1e-5))
!
        fac=1./(60.*dz)
        gg(1:3)=0.; gg(mz-2:mz)=0. !!(border pts to zero)
        gg(4:mz-3)=fac*(45.*(gz(5:mz-2)-gz(3:mz-4)) &
                        -9.*(gz(6:mz-1)-gz(2:mz-5)) &
                           +(gz(7:mz)  -gz(1:mz-6)))
      else
        gz=0
        gg=0
      endif
!
!  make sign antisymmetric
!
      where(z<0)
        ss=-1.
      elsewhere
        ss=1.
      endwhere
      gz1=-ss*gz !!(negative for z>0)
!
!AB: removed nu dependence here. This whole routine is probably not
!AB: needed at the moment, because it is superseded by continuous forcing
!AB: in hydro.f90
!
      !ffnorm=fountain*nu*dt
      ffnorm=fountain*dt
!
!  set forcing function
!
      do n=n1,n2
      do m=m1,m2
        f(l1:l2,m,n,iffx)=f(l1:l2,m,n,iffx)+ffnorm*(+sxx*cy(m)*gz1(n)+cxx*sy(m)*gg(n))
        f(l1:l2,m,n,iffy)=f(l1:l2,m,n,iffy)+ffnorm*(-cxx*sy(m)*gz1(n)+sxx*cy(m)*gg(n))
        f(l1:l2,m,n,iffz)=f(l1:l2,m,n,iffz)+ffnorm*sxx*sy(m)*gz(n)*2.
      enddo
      enddo
!
    endsubroutine forcing_roberts
!***********************************************************************
    subroutine forcing_fountain(f)
!
!  add some artificial fountain flow
!  (to check for example small scale magnetic helicity loss)
!
!  30-may-02/axel: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: sxx,cxx
      real, dimension (mx) :: sx,cx
      real, dimension (my) :: sy,cy
      real, dimension (mz) :: sz,cz,tmpz,gz,gg,ss,gz1
      real :: kx,ky,kz,ffnorm,fac
!
!  identify ourselves
!
      if (headtt.or.ldebug) print*,'forcing_fountain: ENTER'
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
      kx=kfountain
      ky=kfountain
      kz=1.
!
      sx=sin(kx*x); cx=cos(kx*x)
      sy=sin(ky*y); cy=cos(ky*y)
      sz=sin(kz*z); cz=cos(kz*z)
!
!  abbreviation
!
      sxx=sx(l1:l2)
      cxx=cx(l1:l2)
!
!  g(z) and g'(z)
!  use z-profile to cut off
!
      if (height_ff/=0.) then
        tmpz=(z/height_ff)**2
        gz=sz*exp(-tmpz**5/max(1.-tmpz,1e-5))
!
        fac=1./(60.*dz)
        gg(1:3)=0.; gg(mz-2:mz)=0. !!(border pts to zero)
        gg(4:mz-3)=fac*(45.*(gz(5:mz-2)-gz(3:mz-4)) &
                        -9.*(gz(6:mz-1)-gz(2:mz-5)) &
                           +(gz(7:mz)  -gz(1:mz-6)))
      else
        gz=0
        gg=0
      endif
!
!  make sign antisymmetric
!
      where(z<0)
        ss=-1.
      elsewhere
        ss=1.
      endwhere
      gz1=-ss*gz !!(negative for z>0)
!
!AB: removed nu dependence here. This whole routine is probably not
!AB: needed at the moment, because it should be superseded by continuous
!AB: forcing in hydro.f90
!
      !ffnorm=fountain*nu*kfountain**2*dt
      ffnorm=fountain*kfountain**2*dt
!
!  set forcing function
!
      do n=n1,n2
      do m=m1,m2
        f(l1:l2,m,n,iffx)=f(l1:l2,m,n,iffx)+ffnorm*(cxx*sy(m)*gg(n))
        f(l1:l2,m,n,iffy)=f(l1:l2,m,n,iffy)+ffnorm*(sxx*cy(m)*gg(n))
        f(l1:l2,m,n,iffz)=f(l1:l2,m,n,iffz)+ffnorm* sxx*sy(m)*gz(n)*2.
      enddo
      enddo
!
    endsubroutine forcing_fountain
!***********************************************************************
    subroutine forcing_hshear(f)
!
!  add horizontal shear
!
!  19-jun-02/axel+bertil: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx) :: fx
      real, dimension (mz) :: fz
      real :: kx,ffnorm
!
!  need to multiply by dt (for Euler step).
!  Define fz with ghost zones, so fz(n) is the correct position
!  Comment: Brummell et al have a polynomial instead!
!
      kx=2*pi/Lx
      fx=cos(kx*x(l1:l2))
      fz=1./cosh(z/width_ff)**2
      ffnorm=force*dt  !(dt for the timestep)
!
!  add to velocity (here only y-component)
!
      do n=n1,n2
      do m=m1,m2
        f(l1:l2,m,n,iuy)=f(l1:l2,m,n,iuy)+ffnorm*fx*fz(n)
      enddo
      enddo
!
    endsubroutine forcing_hshear
!***********************************************************************
    subroutine forcing_twist(f)
!
!  add circular twisting motion, (ux, 0, uz)
!
!  19-jul-02/axel: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,nz) :: xx,zz,r2,tmp,fx,fz
      real :: ffnorm,ry2,fy,ytwist1,ytwist2
!
!  identifier
!
      if (headt) print*,'forcing_twist: r_ff,width_ff=',r_ff,width_ff
!
!  need to multiply by dt (for Euler step).
!
      ffnorm=force*dt  !(dt for the timestep)
!
!  add to velocity
!  calculate r2=(x^2+z^2)/r^2
!
      xx=spread(x(l1:l2),2,nz)
      zz=spread(z(n1:n2),1,nx)
      r2=(xx**2+zz**2)/r_ff**2
      tmp=exp(-r2/max(1.-r2,1e-5))*ffnorm
      fx=-zz*tmp
      fz=+xx*tmp
!
!  have opposite twists at
!
      ytwist1=y0+0.25*Ly
      ytwist2=y0+0.75*Ly
!
      do m=m1,m2
        !
        ! first twister
        !
        ry2=((y(m)-ytwist1)/width_ff)**2
        fy=exp(-ry2/max(1.-ry2,1e-5))
        f(l1:l2,m,n1:n2,iffx)=f(l1:l2,m,n1:n2,iffx)+fy*fx
        f(l1:l2,m,n1:n2,iffz)=f(l1:l2,m,n1:n2,iffz)+fy*fz
        !
        ! second twister
        !
        ry2=((y(m)-ytwist2)/width_ff)**2
        fy=exp(-ry2/max(1.-ry2,1e-5))
        f(l1:l2,m,n1:n2,iffx)=f(l1:l2,m,n1:n2,iffx)-fy*fx
        f(l1:l2,m,n1:n2,iffz)=f(l1:l2,m,n1:n2,iffz)-fy*fz
      enddo
!
    endsubroutine forcing_twist
!***********************************************************************
    subroutine forcing_diffrot(f,force_ampl)
!
!  Add differential rotation. This routine does not employ continuous
!  forcing, which would be better. It is therefore inferior to the
!  differential rotation procedure implemented directly in hydro.
!
!  26-jul-02/axel: coded
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (nx,nz) :: fx,fz,tmp
      real :: force_ampl,ffnorm,ffnorm2
!
!  identifier
!
      if (headt) print*,'forcing_diffrot: ENTER'
!
!  need to multiply by dt (for Euler step).
!
      ffnorm=force_ampl*dt  !(dt for the timestep)
!
!  prepare velocity, Uy=cosx*cosz
!
      fx=spread(cos(x(l1:l2)),2,nz)
      fz=spread(cos(z(n1:n2)),1,nx)
!
!  this forcing term is balanced by diffusion operator;
!  need to multiply by nu*k^2, but k=sqrt(1+1) for the forcing
!
!AB: removed nu dependence here. This whole routine is probably not
!AB: needed at the moment, because it should be superseded by continuous
!AB: forcing in hydro.f90
!
      !ffnorm2=ffnorm*nu*2
      ffnorm2=ffnorm*2
      tmp=ffnorm2*fx*fz
!
!  add
!
      do m=m1,m2
        f(l1:l2,m,n1:n2,iuy)=f(l1:l2,m,n1:n2,iuy)+tmp
      enddo
!
    endsubroutine forcing_diffrot
!***********************************************************************
    subroutine forcing_blobs(f)
!
!  add blobs in entropy every dtforce time units
!
!  28-jul-02/axel: coded
!
      use Sub
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, save :: tforce=0.
      logical, save :: lfirst_call=.true.
      integer, save :: nforce=0
      logical :: lforce
      character (len=intlen) :: ch
      character (len=fnlen) :: file
!
!  identifier
!
      if (headt) print*,'forcing_blobs: ENTER'
!
!  the last forcing time is recorded in tforce.dat
!
      file=trim(datadir)//'/tforce.dat'
      if (lfirst_call) then
        call read_snaptime(trim(file),tforce,nforce,dtforce,t)
        lfirst_call=.false.
      endif
!
!  Check whether we want to do forcing at this time.
!
      call update_snaptime(file,tforce,nforce,dtforce,t,lforce,ch)
      if (lforce) call blob(force,f,iss,radius_ff,location(1),location(2),location(3))
!
    endsubroutine forcing_blobs
!***********************************************************************
    subroutine forcing_blobHS_random(f)
!
!  add blobs in HS in entropy and density
!  location & time can be random
!
!  08-nov-20/boris: coded
!
      use Sub
      use General, only: random_number_wrapper
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, save :: t_next_blob=1.
      logical, save :: lfirst_call=.true.
      integer, parameter :: t_interval_blobs=50.
      logical :: lforce
      character (len=intlen) :: ch
      character (len=fnlen) :: file
      real, dimension (3) :: fran
      real :: scaled_interval
!
!  identifier
!
      if (headt) print*,'forcing_blobHS_random: ENTER'
!
      if (lfirst_call) then
        t_next_blob=t+1.      !MR: why +1?
        lfirst_call=.false.
      endif
!
!  Check whether we want to do forcing at this time.
!
      if (t>=t_next_blob) then
        if (lrandom_location) then
          call random_number_wrapper(fran,CHANNEL=channel_force)
          location=fran*Lxyz+xyz0
        else
          location=location_fixed(:,1)
        endif
        open (111,file=trim(datadir)//'/tblobs.dat',position="append")
        write (111,'(f12.3,3f8.4)') t_next_blob, location
        close (111)
!
!  Add a blob in HS equilibrium (entropy & density at the same time)
!
        call blob(force,f,iss,radius_ff,location(1),location(2),location(3))
        call blob(-force,f,ilnrho,radius_ff,location(1),location(2),location(3))
!
        if (lrandom_time) then
          call random_number_wrapper(fran)
          scaled_interval=-log(fran(1))*t_interval_blobs
          t_next_blob=t+scaled_interval
        else
          t_next_blob=t+dtforce
        endif
        print*,'t_next_blob=', t_next_blob
      endif
!
    endsubroutine forcing_blobHS_random
!***********************************************************************
    subroutine forcing_hel_smooth(f)
!
      use Diagnostics
      use DensityMethods, only: getrho
      use General, only: random_number_wrapper
      use Mpicomm, only: mpiallreduce_sum
      use Sub
!
!  06-dec-13/nishant: made kkx etc allocatable
!  23-dec-18/axel: forcing_helicity has now similar capabilities
!
      real, dimension (mx,my,mz,mfarray) :: f

      real, dimension (mx,my,mz,3) :: force_vec
      real, dimension (nx) :: ruf,rho
      real, dimension (nx,3) :: variable_rhs,forcing_rhs,force_all
      real :: phase1,phase2,p_weight
      real :: kx01,ky1,kz1,kx02,ky2,kz2
      real :: mulforce_vec
      integer :: ik1,ik2,ik
!
!  Re-calculate forcing wave numbers if necessary
!
      !tsforce is set to -10 in cdata.f90. It should also be saved in a file
      !so that it will be read again on restarts.
      if (t > tsforce) then
        if (tsforce < 0) then
          call random_number_wrapper(fran1,CHANNEL=channel_force)
        else
          fran1=fran2
        endif
        call random_number_wrapper(fran2,CHANNEL=channel_force)
        tsforce=t+dtforce
      endif
      phase1=pi*(2*fran1(1)-1.)
      ik1=nk*.9999*fran1(2)+1
      kx01=kkx(ik1)
      ky1=kky(ik1)
      kz1=kkz(ik1)
      phase2=pi*(2*fran2(1)-1.)
      ik2=nk*.9999*fran2(2)+1
      kx02=kkx(ik2)
      ky2=kky(ik2)
      kz2=kkz(ik2)
!
!  Determine weight parameter
!
      p_weight=(tsforce-t)/dtforce
!
!  Calculate forcing function
!
      force_vec=0.
      call hel_vec(f,kx01,ky1,kz1,phase1,kav,force_vec,   p_weight)
      call hel_vec(f,kx02,ky2,kz2,phase2,kav,force_vec,1.-p_weight)
!
! Find energy input
!
      if (lout.and.idiag_rufm/=0 .or. lwork_ff) then

        lfirstpoint = .true.
        do n=n1,n2
          do m=m1,m2
            forcing_rhs=force_vec(l1:l2,m,n,:)
            variable_rhs=f(l1:l2,m,n,iffx:iffz)!-force_vec(l1:l2,m,n,:)
            call getrho(f(:,m,n,ilnrho),rho)
            call multsv_mn(rho/dt,forcing_rhs,force_all)
            call dot_mn(variable_rhs,force_all,ruf)
            call sum_mn_name(ruf,idiag_rufm)
            lfirstpoint = .false.
          enddo
        enddo
!
! irufm needs to be summed-up globally and communicated to other processors.
!
        if (lout) then
          if (idiag_rufm/=0) call mpiallreduce_sum(fname(idiag_rufm)/nwgrid,fname(idiag_rufm))
        endif
!
! If we want to make energy input constant
!
        if (lwork_ff) then
!
! What should be added to force_vec in order to make the energy
! input equal to work_ff?
!
          mulforce_vec=min(max_force,work_ff/fname(idiag_rufm))
!
!  Add rescaled forcing
!
          f(l1:l2,m1:m2,n1:n2,iffx:iffz) = f(l1:l2,m1:m2,n1:n2,iffx:iffz) &
                                          +force_vec(l1:l2,m1:m2,n1:n2,:)*mulforce_vec
        else
          mulforce_vec=1.
!
!  Add forcing
!
          f(l1:l2,m1:m2,n1:n2,iffx:iffz) = f(l1:l2,m1:m2,n1:n2,iffx:iffz)+force_vec(l1:l2,m1:m2,n1:n2,:)

        endif
!
! Save for printouts.
!
        if (lout.and.idiag_rufm/=0) call save_name(fname(idiag_rufm)*mulforce_vec,idiag_rufm)
      endif
!
    endsubroutine forcing_hel_smooth
!***********************************************************************
    subroutine forcing_tidal(f,force_ampl)
!
!  Added tidal forcing function
!  NB: This is still experimental. Use with care!
!
!  20-Nov-12/simon: coded
!  21-jan-15/MR: changes for use of reference state.
!
      use Diagnostics
      use DensityMethods, only: getrho, getrho1
      use Mpicomm, only: mpireduce_sum
      use Sub
!
      real, dimension (mx,my,mz,mfarray) :: f
!
      real :: force_ampl
      real, dimension (nx) :: ruf,rho,rho1
      real, dimension (nx,3) :: variable_rhs,forcing_rhs,force_all
      real, dimension (nx,3) :: bb,fxb
      integer :: j,jf,l
      real :: fact, dist3
!
!  Normalize ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
      fact=2*force_ampl*sqrt(dt)
!
!  loop the two cases separately, so we don't check for r_ff during
!  each loop cycle which could inhibit (pseudo-)vectorisation
!  calculate energy input from forcing; must use lout (not ldiagnos)
!
      lfirstpoint = .true.
      do n=n1,n2
        do m=m1,m2
          variable_rhs=f(l1:l2,m,n,iffx:iffz)
          if (lmomentum_ff) then
            if (ldensity_nolog) then
              if (lreference_state) then
                rho1=1./(f(l1:l2,m,n,irho)+reference_state(:,iref_rho))
              else
                rho1=1./f(l1:l2,m,n,irho)
              endif
            else
              rho1=exp(-f(l1:l2,m,n,ilnrho))
              call getrho1(f(:,m,n,ilnrho),rho1)
            endif
          else
            rho1=1.
          endif
          do l=l1,l2
            dist3 = sqrt( (R0_tidal*cos(omega_tidal*t)*cos(phi_tidal)-x(l))**2 + &
                          (R0_tidal*sin(omega_tidal*t)               -y(m))**2 + &
                          (R0_tidal*cos(omega_tidal*t)*sin(phi_tidal)-z(n))**2  )**3
            forcing_rhs(l-l1+1,1)=rho1(l)*fact*(R0_tidal*cos(omega_tidal*t)*cos(phi_tidal)-x(l))/dist3
            forcing_rhs(l-l1+1,2)=rho1(l)*fact*(R0_tidal*sin(omega_tidal*t)-y(m))/dist3
            forcing_rhs(l-l1+1,3)=rho1(l)*fact*(R0_tidal*cos(omega_tidal*t)*sin(phi_tidal)-z(n))/dist3
          enddo
          f(l1:l2,m,n,iffx:iffz)=f(l1:l2,m,n,iffx:iffz)+forcing_rhs
          if (lout) then
            if (idiag_rufm/=0) then
              call getrho(f(:,m,n,ilnrho),rho)
              call multsv_mn(rho/dt,forcing_rhs,force_all)
              call dot_mn(variable_rhs,force_all,ruf)
              call sum_mn_name(ruf,idiag_rufm)
            endif
            if (lmagnetic) then
!             if (idiag_fxbxm/=0.or.idiag_fxbym/=0.or.idiag_fxbzm/=0) then
!               call curl(f,iaa,bb)
!               call cross(forcing_rhs,bb,fxb)
!               call sum_mn_name(fxb(:,1),idiag_fxbxm)
!               call sum_mn_name(fxb(:,2),idiag_fxbym)
!               call sum_mn_name(fxb(:,3),idiag_fxbzm)
!             endif
            endif
          endif
          lfirstpoint = .false.
        enddo
      enddo
!
! For printouts, global summation is needed as not done by finalize_diagnostics.
!
      if (lout) then
        if (idiag_rufm/=0) call mpireduce_sum(fname(idiag_rufm)/nwgrid,fname(idiag_rufm))
!       if (lmagnetic) then
!         if (idiag_fxbxm/=0) call mpireduce_sum(fname(idiag_fxbxm)/nwgrid,fname(idiag_fxbxm))
!         if (idiag_fxbym/=0) call mpireduce_sum(fname(idiag_fxbym)/nwgrid,fname(idiag_fxbym))
!         if (idiag_fxbzm/=0) call mpireduce_sum(fname(idiag_fxbzm)/nwgrid,fname(idiag_fxbzm))
!       endif
      endif
!
      if (ip<=9) print*,'forcing_tidal: forcing OK'
!
    endsubroutine forcing_tidal
!***********************************************************************
    subroutine hel_vec(f,kx0,ky,kz,phase,kav,force1,fac)
!
!  Add helical forcing function, using a set of precomputed wavevectors.
!  The relative helicity of the forcing function is determined by the factor
!  sigma, called here also relhel. If it is +1 or -1, the forcing is a fully
!  helical Beltrami wave of positive or negative helicity. For |relhel| < 1
!  the helicity less than maximum. For relhel=0 the forcing is nonhelical.
!  The forcing function is now normalized to unity (also for |relhel| < 1).
!
!  10-apr-00/axel: coded
!   3-sep-02/axel: introduced k1_ff, to rescale forcing function if k1/=1.
!  25-sep-02/axel: preset force_ampl to unity (in case slope is not controlled)
!   9-nov-02/axel: corrected normalization factor for the case |relhel| < 1.
!  17-jan-03/nils: adapted from forcing_hel
!  21-apr-23/MR: changed to cumulative work: NEEDS INITIALIZATION OF force1 BEFORE CALL!
!                added parameter fac for immediate scaling of force1
!
      use General, only: random_number_wrapper
      use Sub
!
      real, dimension (mx,my,mz,mfarray) :: f
      real :: kx0,ky,kz
      real :: phase
      real :: kav
      real :: fac
      real, dimension (mx,my,mz,3) :: force1
!
      real, dimension (nx) :: radius,tmpx
      complex, dimension (mx) :: fx
      complex, dimension (my) :: fy
      complex, dimension (mz) :: fz
      complex, dimension (3) :: coef
      integer :: j,jf
      real :: kx,k2,k,force_ampl,ffnorm
      real :: ex,ey,ez,kde,sig,fact,kex,key,kez,kkex,kkey,kkez
      real, dimension(3) :: e1,e2,ee,kk
      real :: norm,phi
!
!  in the shearing sheet approximation, kx = kx0 - St*k_y.
!  Here, St=-deltay/Lx
!
      if (Sshear==0.) then
        kx=kx0
      else
        kx=kx0+ky*deltay/Lx
      endif
!
      if (headt.or.ip<5) print*, 'hel_vec: kx0,kx,ky,kz=',kx0,kx,ky,kz
      k2=kx**2+ky**2+kz**2
      k=sqrt(k2)
!
! Find e-vector
!
      !
      ! Start with old method (not isotropic) for now.
      ! Pick e1 if kk not parallel to ee1. ee2 else.
      !
      if ((ky==0).and.(kz==0)) then
        ex=0; ey=1; ez=0
      else
        ex=1; ey=0; ez=0
      endif
      if (.not. old_forcing_evector) then
        !
        !  Isotropize ee in the plane perp. to kk by
        !  (1) constructing two basis vectors for the plane perpendicular
        !      to kk, and
        !  (2) choosing a random direction in that plane (angle phi)
        !  Need to do this in order for the forcing to be isotropic.
        !
        kk = (/kx, ky, kz/)
        ee = (/ex, ey, ez/)
        call cross(kk,ee,e1)
        call dot2(e1,norm); e1=e1/sqrt(norm) ! e1: unit vector perp. to kk
        call cross(kk,e1,e2)
        call dot2(e2,norm); e2=e2/sqrt(norm) ! e2: unit vector perp. to kk, e1
        call random_number_wrapper(phi,CHANNEL=channel_force); phi = phi*2*pi
        ee = cos(phi)*e1 + sin(phi)*e2
        ex=ee(1); ey=ee(2); ez=ee(3)
      endif
!
!  k.e
!
      call dot(kk,ee,kde)
!
!  k x e
!
      kex=ky*ez-kz*ey
      key=kz*ex-kx*ez
      kez=kx*ey-ky*ex
!
!  k x (k x e)
!
      kkex=ky*kez-kz*key
      kkey=kz*kex-kx*kez
      kkez=kx*key-ky*kex
!
!  ik x (k x e) + i*phase
!
!  Normalize ff; since we don't know dt yet, we finalize this
!  within timestep where dt is determined and broadcast.
!
!  This does already include the new sqrt(2) factor (missing in B01).
!  So, in order to reproduce the 0.1 factor mentioned in B01
!  we have to set force=0.07.
!
!  Furthermore, for |relhel| < 1, sqrt(2) should be replaced by
!  sqrt(1.+relhel**2). This is done now (9-nov-02).
!  This means that the previous value of force=0.07 (for relhel=0)
!  should now be replaced by 0.05.
!
!  Note: kav is not to be scaled with k1_ff (forcing should remain
!  unaffected when changing k1_ff).
!
      ffnorm=sqrt(1.+relhel**2)*k*sqrt(k2-kde**2)/sqrt(kav*cs0**3)*(k/kav)**slope_ff
      if (ip<=9) then
        print*,'hel_vec: k,kde,ffnorm,kav,dt,cs0=',k,kde,ffnorm,kav,dt,cs0
        print*,'hel_vec: k*sqrt(k2-kde**2)=',k*sqrt(k2-kde**2)
      endif
!
!  need to multiply by dt (for Euler step), but it also needs to be
!  divided by sqrt(dt), because square of forcing is proportional
!  to a delta function of the time difference
!
      fact=force/ffnorm*sqrt(dt)
!
!  The wavevector is for the case where Lx=Ly=Lz=2pi. If that is not the
!  case one needs to scale by 2pi/Lx, etc.
!
      fx=exp(cmplx(0.,kx*k1_ff*x+phase))*fact
      fy=exp(cmplx(0.,ky*k1_ff*y))
      fz=exp(cmplx(0.,kz*k1_ff*z))
!
!  possibly multiply forcing by z-profile
!
!-    if (height_ff/=0.) then
!-      if (lroot) print*,'hel_vec: include z-profile'
!-      tmpz=(z/height_ff)**2
!-      fz=fz*exp(-tmpz**5/max(1.-tmpz,1e-5))
!-    endif
!
!  possibly multiply forcing by sgn(z) and radial profile
!
      if (r_ff/=0.) then
        if (lroot) print*,'hel_vec: applying sgn(z)*xi(r) profile'
        !
        ! only z-dependent part can be done here; radial stuff needs to go
        ! into the loop
        !
        fz = fz*profz_k
      endif
!
      if (ip<=5) then
        print*,'hel_vec: fx=',fx
        print*,'hel_vec: fy=',fy
        print*,'hel_vec: fz=',fz
      endif
!
!  prefactor
!
      coef(1)=cmplx(k*kex,relhel*kkex)
      coef(2)=cmplx(k*key,relhel*kkey)
      coef(3)=cmplx(k*kez,relhel*kkez)
      if (ip<=5) print*,'hel_vec: coef=',coef
!
! loop the two cases separately, so we don't check for r_ff during
! each loop cycle which could inhibit (pseudo-)vectorisation
!
      if (r_ff == 0) then       ! no radial profile
        if (lwork_ff) then
          force_ampl=fac*calc_force_ampl(f,fx,fy,fz,coef)
        else
          force_ampl=fac
        endif
        do j=1,3
          jf=j+ifff-1
          do n=n1,n2
            do m=m1,m2
               force1(l1:l2,m,n,jf) = force1(l1:l2,m,n,jf) + force_ampl*real(coef(j)*fx(l1:l2)*fy(m)*fz(n))
            enddo
          enddo
        enddo
      else                      ! with radial profile
        call fatal_error('hel_vec','radial profile should be quenched')
        do j=1,3
          jf=j+ifff-1
          do n=n1,n2
!--         sig = relhel*tmpz(n)
!AB: removed tmpz factor
            sig = relhel
            coef(1)=cmplx(k*kex,sig*kkex)
            coef(2)=cmplx(k*key,sig*kkey)
            coef(3)=cmplx(k*kez,sig*kkez)
            do m=m1,m2
              radius = sqrt(x(l1:l2)**2+y(m)**2+z(n)**2)
              tmpx = 0.5*(1.-tanh((radius-r_ff)/width_ff))
              force1(l1:l2,m,n,jf) = force1(l1:l2,m,n,jf) + fac*real(coef(j)*tmpx*fx(l1:l2)*fy(m)*fz(n))
            enddo
          enddo
        enddo
      endif
!
      if (ip<=9) print*,'hel_vec: forcing OK'
!
    endsubroutine hel_vec
!***********************************************************************
    subroutine calc_fluxring_cylindrical(force,i)
!
!   4-aug-11/dhruba+axel: adapted from fluxring_cylindrical
!
      real, dimension (nx,3), intent(out):: force
      integer,                intent(in) :: i
!
      real, dimension (nx) :: argum,s
      real ::  b0=1., s0=2., width=.2
!
!  density for the magnetic flux flux ring
!
      s=x(l1:l2)
      argum=(s-s0)/width
      force(:,1)=2.*s*(b0/(width*s0))**2*exp(-2*argum**2)*(width**2-s**2+s*s0)
      force(:,2)=0.
      force(:,3)=0.
      force = fcont_ampl(i)*force
!
    endsubroutine calc_fluxring_cylindrical
!***********************************************************************
    subroutine calc_counter_centrifugal(force,i)
!
!   4-aug-11/dhruba+axel: adapted from fluxring_cylindrical
! Calculates the force required to counter the centrifugal force coming
! from an existing differential rotation.
!
      real, dimension (nx,3), intent(out):: force
      integer,                intent(in) :: i
!
      real, dimension (nx) :: vphi,omega_diffrot
!
      omega_diffrot = ampl_diffrot*x(l1:l2)**(omega_exponent)
      vphi =  x(l1:l2)*omega_diffrot
      force(:,1)=-vphi**2/x(l1:l2)
      force(:,2)=0.
      force(:,3)=0.
      force = fcont_ampl(i)*force
!
    endsubroutine calc_counter_centrifugal
!***********************************************************************
    subroutine calc_GP_TC13(i,sp)
!
!   4-apr-22/hongzhe: The Galloway-Proctor forcing used in Tobias &
! Cattaneo (2013). See also Appendix A of Pongkitiwanichakul+2016.
!
      use General, only: random_number_wrapper
!
      integer,  intent(in) :: i
      character (len=*), intent(in) :: sp
!
      integer :: ii
      real :: k,knorm,omega,tcor,R_GP,xi,eta
!
      sinxt(:,i)=0.
      cosxt(:,i)=0.
      sinyt(:,i)=0.
      cosyt(:,i)=0.
      sinzt(:,i)=0.
      coszt(:,i)=0.
!
      if (nk_GP>100) call not_implemented('calc_GP_TC13','for nk_GP>100')
      do ii=1,nk_GP
        if (nk_GP==1) then
          k = kmin_GP
        else
          k = kmin_GP + (ii-1) * (kmax_GP-kmin_GP)/(nk_GP-1)
        endif
        knorm = k/kmin_GP
        omega = omega_fcont(i)*( knorm**(2-beta_GP) )
        tcor = tcor_GP*( knorm**(beta_GP-2) )
        !
        !  Refresh random parameters
        !
        if (lfirst) then
          call random_number_wrapper(R_GP,CHANNEL=channel_force)
          if ( it==1 .or. R_GP<=(dt/tcor) ) then
            call random_number_wrapper(xi_GP(ii,i), CHANNEL=channel_force)
            call random_number_wrapper(eta_GP(ii,i),CHANNEL=channel_force)
          endif
        endif
        xi = 2*pi*(xi_GP(ii,i)-0.5)
        eta= 2*pi*(eta_GP(ii,i)-0.5)
        !
        !  Compute needed functions
        !
        select case (sp)
        case('xyz')
          sinxt(:,i) = sinxt(:,i) + k*( knorm**(-beta_GP) ) * sin( k*(x-xi +cos(omega*t) ) )
          cosxt(:,i) = cosxt(:,i) + k*( knorm**(-beta_GP) ) * cos( k*(x-xi +cos(omega*t) ) )
          sinyt(:,i) = sinyt(:,i) + k*( knorm**(-beta_GP) ) * sin( k*(y-eta+sin(omega*t) ) )
          cosyt(:,i) = sinyt(:,i) + k*( knorm**(-beta_GP) ) * cos( k*(y-eta+sin(omega*t) ) )
        case('yzx')
          sinzt(:,i) = sinzt(:,i) + k*( knorm**(-beta_GP) ) * sin( k*(z-xi +cos(omega*t) ) )
          coszt(:,i) = coszt(:,i) + k*( knorm**(-beta_GP) ) * cos( k*(z-xi +cos(omega*t) ) )
          sinxt(:,i) = sinxt(:,i) + k*( knorm**(-beta_GP) ) * sin( k*(x-eta+sin(omega*t) ) )
          cosxt(:,i) = sinxt(:,i) + k*( knorm**(-beta_GP) ) * cos( k*(x-eta+sin(omega*t) ) )
        case default
          call fatal_error('calc_GP_TC13','no such sp: '//trim(sp))
        endselect
      enddo
!
    endsubroutine calc_GP_TC13
!***********************************************************************
   subroutine calc_TG_random(i)
!
!  28-sep-22/hongzhe: The Taylor-Green forcing with randomness.
!
      use General, only: random_number_wrapper
!
      integer,  intent(in) :: i
!
      real :: r_TG
!
      call random_number_wrapper(r_TG,CHANNEL=channel_force)
      theta_TG(i)=r_TG*2.*pi
!
    endsubroutine calc_TG_random
!***********************************************************************
    subroutine forcing_after_boundary(f)
!
      real, dimension (mx,my,mz,mfarray),intent(OUT) :: f
!
      if (lforcing_cont) call forcing_cont_after_boundary
!
      if (lff_as_aux) f(l1:l2,m1:m2,n1:n2,ifx:ifz)=0.

    endsubroutine forcing_after_boundary
!***********************************************************************
    subroutine forcing_cont_after_boundary
!
!  precalculate parameters that are new at each timestep,
!  but the same for all pencils
!
      use General, only: random_number_wrapper
!
      integer :: i
      real :: ecost,esint,ecoxt,ecoyt,ecozt
!
!  for the AKA effect, calculate auxiliary functions phi1_ff and phi2_ff
!
      do i=1,n_forcing_cont
        if (iforcing_cont(i)=='AKA') then
          phi1_ff(:,i)=cos(kf_fcont(i)*y+omega_fcont(i)*t)
          phi2_ff(:,i)=cos(kf_fcont(i)*x-omega_fcont(i)*t)
        elseif (iforcing_cont(i)=='GP') then
          ecost=eps_fcont(i)*cos(omega_fcont(i)*t)
          esint=eps_fcont(i)*sin(omega_fcont(i)*t)
          sinxt(:,i)=sin(kf_fcont(i)*x+ecost)
          cosxt(:,i)=cos(kf_fcont(i)*x+ecost)
          sinyt(:,i)=sin(kf_fcont(i)*y+esint)
          cosyt(:,i)=cos(kf_fcont(i)*y+esint)
        elseif (iforcing_cont(i)=='GP_TC13') then
          call calc_GP_TC13(i,'xyz')
        elseif (iforcing_cont(i)=='GP_TC13_yzx') then
          call calc_GP_TC13(i,'yzx')
        elseif (iforcing_cont(i)=='ABCtdep') then
          ecoxt=eps_fcont(i)*cos( omega_fcont(i)*t)
          ecoyt=eps_fcont(i)*cos(omegay_fcont(i)*t)
          ecozt=eps_fcont(i)*cos(omegaz_fcont(i)*t)
          sinxt(:,i)=sin(kf_fcont(i)*x+ecoxt); cosxt(:,i)=cos(kf_fcont(i)*x+ecoxt)
          sinyt(:,i)=sin(kf_fcont(i)*y+ecoyt); cosyt(:,i)=cos(kf_fcont(i)*y+ecoyt)
          sinzt(:,i)=sin(kf_fcont(i)*z+ecozt); coszt(:,i)=cos(kf_fcont(i)*z+ecozt)
        elseif (iforcing_cont(i)=='TG-random-nonhel' .or. &
                iforcing_cont(i)=='TG-random-hel') then
          call calc_TG_random(i)
        endif
      enddo
!
    endsubroutine forcing_cont_after_boundary
!***********************************************************************
    subroutine pencil_criteria_forcing
!
!  All pencils that the Forcing module depends on are specified here.
!
!  24-mar-08/axel: adapted from density.f90
!
      if (lforcing_cont) lpenc_requested(i_fcont)=.true.
      if (iforcing_cont(1)=='(0,cosx*cosz,0)_Lor') lpenc_requested(i_rho1)=.true.
      if (lmomentum_ff) lpenc_requested(i_rho1)=.true.
      if (idiag_qfm/=0) lpenc_diagnos(i_curlo)=.true.
      if (idiag_rufm/=0) lpenc_diagnos(i_rho)=.true.
!
    endsubroutine pencil_criteria_forcing
!***********************************************************************
    subroutine pencil_interdep_forcing(lpencil_in)
!
!  Interdependency among pencils from the Forcing module is specified here.
!
!  24-mar-08/axel: adapted from density.f90
!
      logical, dimension(npencils) :: lpencil_in
!
      call keep_compiler_quiet(lpencil_in)
!
    endsubroutine pencil_interdep_forcing
!***********************************************************************
    subroutine calc_pencils_forcing(f,p)
!
!  Calculate forcing pencils.
!
!  24-mar-08/axel: adapted from density.f90
!   6-feb-09/axel: added epsilon factor in ABC flow (eps_fcont=1. -> nocos)
!
      use Sub, only: multsv_mn
!
      real, dimension (mx,my,mz,mfarray) :: f
      type (pencil_case) :: p
!
      intent(inout) :: f,p
!
      integer :: i
!
!  calculate forcing
!
      if (lpencil(i_fcont)) then
        if (headtt) print*,'forcing: add continuous forcing'

        do i=1,n_forcing_cont
          call forcing_cont(i,p%fcont(:,:,i),rho1=p%rho1)
          ! put force into auxiliary variable, if requested
          if (lff_as_aux) then
            if (i == 1) then
              f(l1:l2,m,n,ifx:ifz) = p%fcont(:,:,i)
            else
              f(l1:l2,m,n,ifx:ifz) = f(l1:l2,m,n,ifx:ifz) + p%fcont(:,:,i)
            endif
          endif
!
!  divide by rho if lmomentum_ff=T
!  MR: better to place it in hydro
!
          if (i==1 .and. lmomentum_ff) call multsv_mn(p%rho1,p%fcont(:,:,1),p%fcont(:,:,1))
        enddo
      endif
!
    endsubroutine calc_pencils_forcing
!***********************************************************************
    subroutine random_isotropic_KS_setup(initpower,kmin,kmax)
!
!  Produces random, isotropic field from energy spectrum following the
!  KS method (Malik and Vassilicos, 1999.)
!
!  More to do; unsatisfactory so far - at least for a steep power-law
!  energy spectrum.
!
!  27-may-05/tony: modified from snod's KS hydro initial
!  03-feb-06/weezy: Attempted rewrite to guarantee periodicity of KS modes.
!
      use Sub, only: cross, dot2
      use General, only: random_number_wrapper
!
      integer :: modeN
!
      real, dimension (3) :: k_unit
      real, dimension (3) :: ee,e1,e2
      real, dimension (6) :: r
      real :: initpower,kmin,kmax
      real, dimension(KS_modes) :: k,dk,energy,ps
      real :: theta,phi,alpha,beta
      real :: ex,ey,ez,norm,a
!
      if (.not.allocated(KS_k)) then
        allocate(KS_k(3,KS_modes))
        allocate(KS_A(3,KS_modes))
        allocate(KS_B(3,KS_modes))
        allocate(KS_omega(KS_modes))
      endif
!
      kmin=2.*pi      !/(1.0*Lxyz(1))
      kmax=128.*pi    !nx*pi
      a=(kmax/kmin)**(1./(KS_modes-1.))
!
!  Loop over all modes.
!
      do modeN=1,KS_modes
!
!  Pick wavenumber.
!
        k=kmin*(a**(modeN-1.))
!
!  Pick 4 random angles for each mode.
!
        call random_number_wrapper(r);
        theta=pi*(r(1) - 0.)
        phi=pi*(2*r(2) - 0.)
        alpha=pi*(2*r(3) - 0.)
        beta=pi*(2*r(4) - 0.)
!
!  Make a random unit vector by rotating fixed vector to random position
!  (alternatively make a random transformation matrix for each k).
!
        k_unit(1)=sin(theta)*cos(phi)
        k_unit(2)=sin(theta)*sin(phi)
        k_unit(3)=cos(theta)
!
        energy=(((k/kmin)**2. +1.)**(-11./6.))*(k**2.)*exp(-0.5*(k/kmax)**2.)
!
!  Make a vector KS_k of length k from the unit vector for each mode.
!
        KS_k(:,modeN)=k*k_unit(:)
        KS_omega(:)=sqrt(energy(:)*(k(:)**3.))
!
!  Construct basis for plane having rr normal to it
!  (bit of code from forcing to construct x', y').
!
      if ((k_unit(2)==0).and.(k_unit(3)==0)) then
          ex=0.; ey=1.; ez=0.
        else
          ex=1.; ey=0.; ez=0.
        endif
        ee = (/ex, ey, ez/)
!
        call cross(k_unit(:),ee,e1)
!  e1: unit vector perp. to KS_k
        call dot2(e1,norm); e1=e1/sqrt(norm)
        call cross(k_unit(:),e1,e2)
!  e2: unit vector perp. to KS_k, e1
        call dot2(e2,norm); e2=e2/sqrt(norm)
!
!  Make two random unit vectors KS_B and KS_A in the constructed plane.
!
        KS_A(:,modeN) = cos(alpha)*e1 + sin(alpha)*e2
        KS_B(:,modeN) = cos(beta)*e1  + sin(beta)*e2
!
!  Make sure dk is set.
!
        call error('random_isotropic_KS_setup', 'Using uninitialized dk')
        dk=0.                     ! to make compiler happy
!
        ps=sqrt(2.*energy*dk)   !/3.0)
!
!  Give KS_A and KS_B length ps.
!
        KS_A(:,modeN)=ps*KS_A(:,modeN)
        KS_B(:,modeN)=ps*KS_B(:,modeN)
!
      enddo
!
!  Form RA = RA x k_unit and RB = RB x k_unit.
!  Note: cannot reuse same vector for input and output.
!
      do modeN=1,KS_modes
        call cross(KS_A(:,modeN),k_unit(:),KS_A(:,modeN))
        call cross(KS_B(:,modeN),k_unit(:),KS_B(:,modeN))
      enddo
!
      call keep_compiler_quiet(initpower)
!
    endsubroutine random_isotropic_KS_setup
!***********************************************************************
    subroutine random_isotropic_KS_setup_test
!
!  Produces random, isotropic field from energy spectrum following the
!  KS method (Malik and Vassilicos, 1999.)
!  This test case only uses 3 very specific modes (useful for comparison
!  with Louise's kinematic dynamo code.
!
!  03-feb-06/weezy: modified from random_isotropic_KS_setup
!
      use Sub, only: cross
      use General, only: random_number_wrapper
!
      integer :: modeN
!
      real, dimension (3,KS_modes) :: k_unit
      real, dimension(KS_modes) :: k,dk,energy,ps
      real :: initpower,kmin,kmax
!
      if (.not.allocated(KS_k)) then
        allocate(KS_k(3,KS_modes))
        allocate(KS_A(3,KS_modes))
        allocate(KS_B(3,KS_modes))
        allocate(KS_omega(KS_modes))
      endif
!
      initpower=-5./3.
      kmin=10.88279619
      kmax=23.50952672
!
      KS_k(1,1)=2.00*pi
      KS_k(2,1)=-2.00*pi
      KS_k(3,1)=2.00*pi
!
      KS_k(1,2)=-4.00*pi
      KS_k(2,2)=0.00*pi
      KS_k(3,2)=2.00*pi
!
      KS_k(1,3)=4.00*pi
      KS_k(2,3)=2.00*pi
      KS_k(3,3)=-6.00*pi
!
      KS_k(1,1)=+1; KS_k(2,1)=-1; KS_k(3,1)=1
      KS_k(1,2)=+0; KS_k(2,2)=-2; KS_k(3,2)=1
      KS_k(1,3)=+0; KS_k(2,3)=-0; KS_k(3,3)=1
!
      k(1)=kmin
      k(2)=14.04962946
      k(3)=kmax
!
      do modeN=1,KS_modes
        k_unit(:,modeN)=KS_k(:,modeN)/k(modeN)
      enddo
!
      kmax=k(KS_modes)
      kmin=k(1)
!
      do modeN=1,KS_modes
        if (modeN==1) dk(modeN)=(k(modeN+1)-k(modeN))/2.
        if (modeN>1.and.modeN<KS_modes) dk(modeN)=(k(modeN+1)-k(modeN-1))/2.
        if (modeN==KS_modes) dk(modeN)=(k(modeN)-k(modeN-1))/2.
      enddo
!
      do modeN=1,KS_modes
         energy(modeN)=((k(modeN)**2 +1.)**(-11./6.))*(k(modeN)**2)*exp(-0.5*(k(modeN)/kmax)**2)
      enddo
!
      ps=sqrt(2.*energy*dk)
!
      KS_A(1,1)=1.00/sqrt(2.00)
      KS_A(2,1)=-1.00/sqrt(2.00)
      KS_A(3,1)=0.00
!
      KS_A(1,2)=1.00/sqrt(3.00)
      KS_A(2,2)=1.00/sqrt(3.00)
      KS_A(3,2)=-1.00/sqrt(3.00)
!
      KS_A(1,3)=-1.00/2.00
      KS_A(2,3)=-1.00/2.00
      KS_A(3,3)=1.00/sqrt(2.00)
!
      KS_B(1,3)=1.00/sqrt(2.00)
      KS_B(2,3)=-1.00/sqrt(2.00)
      KS_B(3,3)=0.00
!
      KS_B(1,1)=1.00/sqrt(3.00)
      KS_B(2,1)=1.00/sqrt(3.00)
      KS_B(3,1)=-1.00/sqrt(3.00)
!
      KS_B(1,2)=-1.00/2.00
      KS_B(2,2)=-1.00/2.00
      KS_B(3,2)=1.00/sqrt(2.00)
!
      do modeN=1,KS_modes
        KS_A(:,modeN)=ps(modeN)*KS_A(:,modeN)
        KS_B(:,modeN)=ps(modeN)*KS_B(:,modeN)
      enddo
!
!  Form RA = RA x k_unit and RB = RB x k_unit.
!
       do modeN=1,KS_modes
         call cross(KS_A(:,modeN),k_unit(:,modeN),KS_A(:,modeN))
         call cross(KS_B(:,modeN),k_unit(:,modeN),KS_B(:,modeN))
       enddo
!
    endsubroutine random_isotropic_KS_setup_test
!***********************************************************************
    subroutine forcing_cont(i,force,rho1)
!
!   9-apr-10/MR: added RobertsFlow_exact forcing, compensates \nu\nabla^2 u
!                and u.grad u for Roberts geometry
!   4-nov-11/MR: now also compensates Coriolis force
!   16-jun-23/Kishore G: Added ability to read forcing function from a file.
!
!  Note: It is not enough to set lforcing_cont = T in input parameters of
!  forcing one must also set  lforcing_cont_uu = T in hydro for the
!  continuous-in-time forcing to be added to velocity.
!  Alternatively, one can set lforcing_cont_bb = T in magnetic.
!
      use Gravity, only: gravz
      use Sub, only: quintic_step, quintic_der_step, step, vortex
      use Viscosity, only: getnu
!
      integer,                intent(in) :: i
      real, dimension (nx,3), intent(out):: force
      real, dimension (nx), optional, intent(in) :: rho1
!
      real, dimension (nx) :: tmp
      real :: fact, fact1, fact2, fpara, dfpara, sqrt21k1
      real :: kf, kx, ky, kz, nu, arg, ecost, esint
      integer :: i2d1=1,i2d2=2,i2d3=3,modeN
      real, dimension(nx) :: kdotxwt, cos_kdotxwt, sin_kdotxwt
!
        select case (iforcing_cont(i))
        case('Fy=const')
          force(:,1)=0.
          force(:,2)=ampl_ff(i)
          force(:,3)=0.
        case('Fz=const')
          force(:,1)=0.
          force(:,2)=0.
          force(:,3)=ampl_ff(i)
        case('ABC')
          fact2=relhel
          fact1=.5*(1.+relhel**2)
          fact=ampl_ff(i)/sqrt(fact1*(ABC_A(i)**2+ABC_B(i)**2+ABC_C(i)**2))
          if (lforcing_osc_double) then
            force(:,1)=fact*(ABC_C(i)*(sin(kf_fcont(i)*z(n)    +omega_ff       *t) &
                       +ampl_double_ff*cos(kf_fcont(i)*z(n)    +omega_double_ff*t)) &
                      +fact2*ABC_B(i)*(cos(kf_fcont(i)*y(m)    +omega_ff       *t) &
                       +ampl_double_ff*sin(kf_fcont(i)*y(m)    +omega_double_ff*t)))
            force(:,2)=fact*(ABC_A(i)*(sin(kf_fcont(i)*x(l1:l2)+omega_ff       *t) &
                       +ampl_double_ff*cos(kf_fcont(i)*x(l1:l2)+omega_double_ff*t)) &
                      +fact2*ABC_C(i)*(cos(kf_fcont(i)*z(n)    +omega_ff       *t) &
                       +ampl_double_ff*sin(kf_fcont(i)*z(n)    +omega_double_ff*t)))
            force(:,3)=fact*(ABC_B(i)*(sin(kf_fcont(i)*y(m)    +omega_ff       *t) &
                       +ampl_double_ff*cos(kf_fcont(i)*y(m)    +omega_double_ff*t)) &
                      +fact2*ABC_A(i)*(cos(kf_fcont(i)*x(l1:l2)+omega_ff       *t) &
                       +ampl_double_ff*sin(kf_fcont(i)*x(l1:l2)+omega_double_ff*t)))
          else
            if (lforcing_osc) then
              fact=fact*(cos(omega_ff*t)+cos(omega_double_ff*t))
            else if (lforcing_osc2) then
              fact=fact*(cos(omega_ff*t)+cos(omega_double_ff*t))**2
            endif
            force(:,1)=fact*(ABC_C(i)*sinz(n    ,i)+fact2*ABC_B(i)*cosy(m    ,i))
            force(:,2)=fact*(ABC_A(i)*sinx(l1:l2,i)+fact2*ABC_C(i)*cosz(n    ,i))
            force(:,3)=fact*(ABC_B(i)*siny(m    ,i)+fact2*ABC_A(i)*cosx(l1:l2,i))
          endif
        case('Schur_nonhelical')
          force(:,1)= ampl_ff(i)*(cosx(l1:l2,i)*siny(m,i))
          force(:,2)=-ampl_ff(i)*(sinx(l1:l2,i)*cosy(m,i))
          force(:,3)= ampl_ff(i)*(sinx(l1:l2,i)*siny(m,i)*cosz(n,i))
        case('Schur_helical')
          force(:,1)=     0.5*ampl_ff(i)*(           cosx(l1:l2,i)*siny(m,i))+ &
                      0.5*0.5*ampl_ff(i)*( sqrt(2.0)*cosx(l1:l2,i)*siny(m,i) +  sinx(l1:l2,i)*cosy(m,i))
          force(:,2)=    -0.5*ampl_ff(i)*(           sinx(l1:l2,i)*cosy(m,i))- &
                      0.5*0.5*ampl_ff(i)*(-sqrt(2.0)*sinx(l1:l2,i)*cosy(m,i) -  cosx(l1:l2,i)*siny(m,i))
          force(:,3)=     0.5*ampl_ff(i)*(           sinx(l1:l2,i)*siny(m,i))+ &
                      0.5*0.5*ampl_ff(i)*( sqrt(2.0)*sinx(l1:l2,i)*siny(m,i) -2*cosx(l1:l2,i)*cosy(m,i))
        case('ABCtdep')
          fact2=relhel
          fact=ampl_ff(i)/sqrt(ABC_A(i)**2+ABC_B(i)**2+ABC_C(i)**2)
          force(:,1)=fact*(ABC_C(i)*sinzt(n    ,i)+fact2*ABC_B(i)*cosyt(m    ,i))
          force(:,2)=fact*(ABC_A(i)*sinxt(l1:l2,i)+fact2*ABC_C(i)*coszt(n    ,i))
          force(:,3)=fact*(ABC_B(i)*sinyt(m    ,i)+fact2*ABC_A(i)*cosxt(l1:l2,i))
        case ('AKA')
          fact=sqrt(2.)*ampl_ff(i)
          force(:,1)=fact*phi1_ff(m    ,i)
          force(:,2)=fact*phi2_ff(l1:l2,i)
          force(:,3)=fact*(phi1_ff(m,i)+phi2_ff(l1:l2,i))
        case ('grav_z')
          force(:,1)=0
          force(:,2)=0
          force(:,3)=gravz*ampl_ff(i)*cos(omega_ff*t)
        case ('uniform_vorticity')
          if (lcylindrical_coords) then
            force(:,1)=-x(l1:l2)*y(m)*ampl_ff(i)*cos(omega_ff*t)
            force(:,2)=x(l1:l2)*ampl_ff(i)*cos(omega_ff*t)
            force(:,3)=0
          else
            force(:,1)=z(n)*ampl_ff(i)*cos(omega_ff*t)
            force(:,2)=0
            force(:,3)=-x(l1:l2)*ampl_ff(i)*cos(omega_ff*t)
          endif
        case('KolmogorovFlow-x')
          fact=ampl_ff(i)
          force(:,1)=0
          force(:,2)=fact*cosx(l1:l2,i)
          force(:,3)=0
        case('KolmogorovFlow-z')
          fact=ampl_ff(i)
          force(:,1)=fact*cosz(n,i)
          force(:,2)=0.
          force(:,3)=0.
!
!  nocos (=Archontis flow), here with temporal modulation
!  if omega_ff or omega_double_ff are non-vanishing.
!  For lforcing_osc2: quadratic time modulation.
!
        case ('nocos')
          fact=ampl_ff(i)
          if (lforcing_osc) then
            fact=fact*(cos(omega_ff*t)+cos(omega_double_ff*t))
          else if (lforcing_osc2) then
            fact=fact*(cos(omega_ff*t)+cos(omega_double_ff*t))**2
          endif
          force(:,1)=fact*sinz(n    ,i)
          force(:,2)=fact*sinx(l1:l2,i)
          force(:,3)=fact*siny(m    ,i)
!
!  Amplitude is given by ampl_ff/nu*k^2, and k^2=2 in a 2-D box
!  of size (2pi)^2, for example.
!
        case('Straining')
          fact=ampl_ff(i)
          fact2=-(dimensionality-1)*ampl_ff(i)
          force(:,1)=fact *sinx(l1:l2,i)*cosy(m,i)*cosz(n,i)
          force(:,2)=fact *cosx(l1:l2,i)*siny(m,i)*cosz(n,i)
          force(:,3)=fact2*cosx(l1:l2,i)*cosy(m,i)*sinz(n,i)
!
!  Amplitude is given by ampl_ff/nu*k^2, and k^2=2 in a 2-D box
!  of size (2pi)^2, for example.
!
        case('StrainingExcact')
          call fatal_error("forcing_cont","'StrainingExcact' not tested yet")
          call getnu(nu_input=nu)
          kx=k1_ff; kz=k1_ff
          kf=sqrt(kx**2+kz**2)
          fact=ampl_ff(i)*kf**2*nu
          fact1=ampl_ff(i)*ampl_ff(i)*kx*kz
          fact2=-(dimensionality-1)*ampl_ff(i)*kf**2*nu
          force(:,1)=fact *sinx(l1:l2,i)*cosy(m,i)*cosz(n,i)-fact1*kz*cosx(l1:l2,i)*cosy(m,i)*sinz(n,i)
          force(:,2)=fact *cosx(l1:l2,i)*siny(m,i)*cosz(n,i)
          force(:,3)=fact2*cosx(l1:l2,i)*cosy(m,i)*sinz(n,i)-fact1*kx*sinx(l1:l2,i)*cosy(m,i)*cosz(n,i)
        case('RobertsFlow')
          fact=ampl_ff(i)
          force(:,1)=-fact*cosx(l1:l2,i)*siny(m,i)
          force(:,2)=+fact*sinx(l1:l2,i)*cosy(m,i)
          force(:,3)=+relhel*fact*cosx(l1:l2,i)*cosy(m,i)*sqrt2
        case('RobertsFlowII')
          fact=ampl_ff(i)
          force(:,1)=-fact*cosx(l1:l2,i)*siny(m,i)
          force(:,2)=+fact*sinx(l1:l2,i)*cosy(m,i)
          force(:,3)=+relhel*fact*sinx(l1:l2,i)*siny(m,i)*sqrt2
        case('RobertsFlowMask')
          tmp=ampl_ff(i)*profx_ampl*profy_ampl(m)
          force(:,1)=-tmp*cosx(l1:l2,i)*siny(m,i)
          force(:,2)=+tmp*sinx(l1:l2,i)*cosy(m,i)
          force(:,3)=+relhel*tmp*cosx(l1:l2,i)*cosy(m,i)*sqrt2
        case('RobertsFlow2d')
          fact=ampl_ff(i)
          i2d1=1;i2d2=2;i2d3=3
          if (l2dxz) then
            i2d1=2;i2d2=1;i2d3=2
          endif
          if (l2dyz) then
            i2d1=3;i2d2=2;i2d3=1
          endif
          force(:,i2d1)=-fact*cos(k2d*x(l1:l2))*sin(k2d*y(m))
          force(:,i2d2)=+fact*sin(k2d*x(l1:l2))*cos(k2d*y(m))
          force(:,i2d3)= 0.
        case ('RobertsFlow_exact')
          kx=kf_fcont(i); ky=kf_fcont(i)
          kf=sqrt(kx*kx+ky*ky)
          call getnu(nu_input=nu)
          fact=ampl_ff(i)*kf*kf*nu
          fact2=ampl_ff(i)*ampl_ff(i)*kx*ky
!
!!print*, 'forcing: kx, ky, kf, fact, fact2=', kx, ky, kf, fact, fact2
!
          force(:,1)=-fact*ky*cosx(l1:l2,i)*siny(m,i) - fact2*ky*sinx(l1:l2,i)*cosx(l1:l2,i)
          force(:,2)=+fact*kx*sinx(l1:l2,i)*cosy(m,i) - fact2*kx*siny(m,i)*cosy(m,i)
          force(:,3)=+fact*relhel*kf*cosx(l1:l2,i)*cosy(m,i)
!
          if ( Omega/=0. .and. theta==0. ) then  ! Obs, only implemented for rotation axis in z direction.
            fact = 2.*ampl_ff(i)*Omega
            force(:,1)= force(:,1)-fact*kx*sinx(l1:l2,i)*cosy(m,i)
            force(:,2)= force(:,2)-fact*ky*cosx(l1:l2,i)*siny(m,i)
          endif
!
        case ('RobertsFlow-zdep')
          if (headtt) print*,'z-dependent Roberts flow; eps_fcont=',eps_fcont
          fpara = quintic_step(z(n),-1.+eps_fcont(i),eps_fcont(i)) &
                 -quintic_step(z(n),+1.-eps_fcont(i),eps_fcont(i))
          dfpara = quintic_der_step(z(n),-1.+eps_fcont(i),eps_fcont(i)) &
                  -quintic_der_step(z(n),+1.-eps_fcont(i),eps_fcont(i))
!
!  abbreviations
!
          sqrt21k1=1./(sqrt2*k1_ff)
!
!  amplitude factor missing in upper lines
!
          force(:,1)=-ampl_ff(i)*cosx(l1:l2,i)*siny(m,i) &
              -dfpara*ampl_ff(i)*sinx(l1:l2,i)*cosy(m,i)*sqrt21k1
          force(:,2)=+ampl_ff(i)*sinx(l1:l2,i)*cosy(m,i) &
              -dfpara*ampl_ff(i)*cosx(l1:l2,i)*siny(m,i)*sqrt21k1
          force(:,3)=+fpara*ampl_ff(i)*cosx(l1:l2,i)*cosy(m,i)*sqrt2
!
        case ('elevator-flow')
          if (headtt) print*,'z-dependent elevator-flow; eps_fcont=',eps_fcont
          fpara=quintic_step(z(n),-1.+eps_fcont(i),eps_fcont(i)) &
               -quintic_step(z(n),+1.-eps_fcont(i),eps_fcont(i))
          dfpara=quintic_der_step(z(n),-1.+eps_fcont(i),eps_fcont(i))&
                -quintic_der_step(z(n),+1.-eps_fcont(i),eps_fcont(i))
!
!  abbreviations
!
          sqrt21k1=1./(sqrt2*k1_ff)
!
!  amplitude factor missing in upper lines
!
          force(:,1)=-dfpara*ampl_ff(i)*sinx(l1:l2,i)*cosy(m,i)*sqrt21k1
          force(:,2)=-dfpara*ampl_ff(i)*cosx(l1:l2,i)*siny(m,i)*sqrt21k1
          force(:,3)=+fpara*ampl_ff(i)*cosx(l1:l2,i)*cosy(m,i)*sqrt2
!
!  f=(sinx,0,0)
!
        case('Roberts-for-SSD')
          fact=ampl_ff(i)
          force(:,1)=-fact*cosx(l1:l2,i)*siny(m,i)
          force(:,2)=+fact*sinx(l1:l2,i)*cosy(m,i)
          force(:,3)= fact*rel_zcomp*(cosx(l1:l2,i)**2-0.5)*(cosy(m,i)**2-0.5)
!
        case ('sinx')
          if (lgentle(i).and.t<tgentle(i)) then
            fact=.5*ampl_ff(i)*(1.-cos(pi*t/tgentle(i)))
          else
            fact=ampl_ff(i)
          endif
          force(:,1)=fact*sinx(l1:l2,i)
          force(:,2)=0.
          force(:,3)=0.
!
!  f=(0,0,cosx)
!
        case ('(0,0,cosx)')
          force(:,1)=0.
          force(:,2)=0.
          force(:,3)=ampl_ff(i)*cosx(l1:l2,i)
!
!  f=(0,0,cosx)
!
        case ('(0,0,cosxcosy)')
          force(:,1)=0.
          force(:,2)=0.
          force(:,3)=ampl_ff(i)*cosx(l1:l2,i)*cosy(m,i)
!
!  f=B=(0,0,cosx*cosy)
!
        case ('B=(0,0,cosxcosy)')
          force(:,1)=-0.5*ampl_ff(i)*cosx(l1:l2,i)*siny(m,i)
          force(:,2)= 0.5*ampl_ff(i)*sinx(l1:l2,i)*cosy(m,i)
          force(:,3)= 0.
!
!  f=(0,x,0)
!
        case ('(0,x,0)')
          force(:,1)=0.
          force(:,2)=ampl_ff(i)*x(l1:l2)
          force(:,3)=0.
!
!  f=(0,sinx*sint,0)
!
        case ('(0,sinxsint,0)')
          force(:,1)=0.
          force(:,2)=ampl_ff(i)*sinx(l1:l2,i)*sin(omega_ff*t)
          force(:,3)=0.
!
!  f=(0,sinx,0)
!
        case ('(0,sinx,0)')
          force(:,1)=0.
          force(:,2)=ampl_ff(i)*sinx(l1:l2,i)
          force(:,3)=0.
!
!  f=(0,cosx*cosz,0), modulated by profy_ampl
!
        case ('(0,cosx*cosz,0)')
          force(:,1)=0.
          force(:,2)=ampl_ff(i)*cosx(l1:l2,i)*cosz(n,i)*profy_ampl(m)
          force(:,3)=0.
!
!  f=(0,sinx*exp(-z^2),0)
!
        case ('(0,sinx*exp(-z^2),0)')
          force(:,1)=0.
          force(:,2)= ampl_ff(i)*sinx(l1:l2,i)*exp(-((z(n)-r_ff)/width_ff)**2) &
                     *(step(x(l1:l2),xminf,2.)-step(x(l1:l2),xmaxf,2.))
          force(:,3)=0.
!
!  f=(0,Aycont_z,0)
!  This ensures vanishing Ay at both boundaries if Bslope=-2Bconst/Lz
!  (making it suitable for perfect conductor boundary condition)
!
        case ('(0,Aycont_z,0)')
          force(:,1)=0.
          force(:,2)=-ampl_ff(i)*((Bconst*(z(n)-xyz0(3)))+(0.5*Bslope*((z(n)-xyz0(3))**2.)))
          force(:,3)=0.
!
!  f=(sinz,cosz,0)
!
        case ('(sinz,cosz,0)')
          force(:,1)=ampl_ff(i)*sinz(n,i)
          force(:,2)=ampl_ff(i)*cosz(n,i)
          force(:,3)=0.
!
!  Taylor-Green forcing
!
        case ('TG')
          fact=2.*ampl_ff(i)
          force(:,1)=+fact*sinx(l1:l2,i)*cosy(m,i)*cosz(n,i)
          force(:,2)=-fact*cosx(l1:l2,i)*siny(m,i)*cosz(n,i)
          force(:,3)=0.
!
!  Taylor-Green forcing with randomness
!
        case ('TG-random-nonhel')
          fact=ampl_ff(i)*sqrt(dt)
          force(:,1)=fact*sin(theta_TG(i)+2.*pi/3)/sqrt(3.)*sinx(l1:l2,i)*cosy(m,i)*cosz(n,i)
          force(:,2)=fact*sin(theta_TG(i)-2.*pi/3)/sqrt(3.)*cosx(l1:l2,i)*siny(m,i)*cosz(n,i)
          force(:,3)=fact*sin(theta_TG(i))/sqrt(3.)*cosx(l1:l2,i)*cosy(m,i)*sinz(n,i)
!
!  Fully helical Taylor-Green forcing with randomness
!
        case ('TG-random-hel')
          fact=ampl_ff(i)*sqrt(dt)
          force(:,1)=fact*(cos(theta_TG(i)+pi/6.)/sqrt(6.)*sinx(l1:l2,i)*cosy(m,i)*cosz(n,i) &
          -(sin(theta_TG(i))+cos(theta_TG(i)-pi/6.))/sqrt(18.)*cosx(l1:l2,i)*siny(m,i)*sinz(n,i))
          force(:,2)=fact*(-cos(theta_TG(i)-pi/6.)/sqrt(6.)*cosx(l1:l2,i)*siny(m,i)*cosz(n,i) &
          +(sin(theta_TG(i))-cos(theta_TG(i)+pi/6.))/sqrt(18.)*sinx(l1:l2,i)*cosy(m,i)*sinz(n,i))
          force(:,3)=fact*(sin(theta_TG(i))/sqrt(6.)*cosx(l1:l2,i)*cosy(m,i)*sinz(n,i) &
          +(cos(theta_TG(i)-pi/6.)+cos(theta_TG(i)+pi/6.))/sqrt(18.)*sinx(l1:l2,i)*siny(m,i)*cosz(n,i))
!
!  Compressive u=-grad(phi) with phi=cos(x+y+z) forcing
!
        case ('cosx*cosy*cosz')
          fact=-ampl_ff(i)
          force(:,1)=fact*sinx(l1:l2,i)*cosy(m,i)*cosz(n,i)
          force(:,2)=fact*cosx(l1:l2,i)*siny(m,i)*cosz(n,i)
          force(:,3)=fact*cosx(l1:l2,i)*cosy(m,i)*sinz(n,i)
!
!  Galloway-Proctor (GP) forcing
!
        case ('GP')
          fact=ampl_ff(i)/sqrt(2.)
          force(:,1)=-fact*                sinyt(m,i)
          force(:,2)=+fact* sinxt(l1:l2,i)
          force(:,3)=-fact*(cosxt(l1:l2,i)+cosyt(m,i))
!
!  Galloway-Proctor (GP) forcing
!
        case ('Galloway-Proctor-92')
          fact=ampl_ff(i)
          ecost=eps_fcont(i)*cos(omega_fcont(i)*t)
          esint=eps_fcont(i)*sin(omega_fcont(i)*t)
          force(:,1)=fact*(sin(kf_fcont(i)*z(n)+esint)+cos(kf_fcont(i)*y(m)+ecost))
          force(:,2)=fact* cos(kf_fcont(i)*z(n)+esint)
          force(:,3)=fact* sin(kf_fcont(i)*y(m)+ecost)
!
!  Galloway-Proctor (GP) forcing in Tobias & Cattaneo (2013)
!
        case ('GP_TC13')
          fact=ampl_ff(i)
          force(:,1) = -fact*sinyt(m,i)
          force(:,2) = -fact*cosxt(l1:l2,i)
          force(:,3) = +fact*( sinxt(l1:l2,i)+cosyt(m,i) )
!
!  Same as GP_TC13 with a rotation of axes y->z->x->y.
!  i.e., the new z axis is the old y axis, etc.
!  This is for the convenience of testfield_z.
!
        case ('GP_TC13_yzx')
          fact=ampl_ff(i)
          force(:,1) = -fact*coszt(n,i)
          force(:,2) = +fact*( sinzt(n,i)+cosxt(l1:l2,i) )
          force(:,3) = -fact*sinxt(l1:l2,i)
!
! Continuous emf required in Induction equation for the Mag Buoy Inst
!
        case('mbi_emf')
          fact=2.*ampl_bb(i)*eta_bb(i)/width_bb(i)**2
          arg=(z(n)-z_bb(i))/width_bb(i)
          force(:,1)=fact*tanh(arg)/(cosh(arg))**2
          force(:,2)=0.0
          force(:,3)=0.0
!
!  Bessel function forcing
!
        case('J0_k1x')
          force(:,1)=0.0
          force(:,2)=profx_ampl1
          force(:,3)=profx_ampl
!
!  Gaussian blob in the z direction
!
        case('gaussian-z')
          force(:,1)=0.0
          force(:,2)=0.0
          force(:,3)=profx_ampl*profy_ampl(m)*profz_ampl(n)
!
!  fluxring_cylindrical
!
        case('fluxring_cylindrical')
          call calc_fluxring_cylindrical(force,i)
        case('counter_centrifugal')
          call calc_counter_centrifugal(force,i)
        case('vortex')
          call vortex(torus,Omega_vortex,force)
!
!  blob-like disturbance (gradient of gaussian)
!
        case('blob')
          fact=ampl_ff(i)/radius_ff**2
          tmp=fact*exp(-.5*((x(l1:l2)-location_fixed(1,1))**2 &
                           +(y(m)    -location_fixed(2,1))**2 &
                           +(z(n)    -location_fixed(3,1))**2)/radius_ff**2)
          force(:,1)=tmp*(x(l1:l2)-location_fixed(1,1))
          force(:,2)=tmp*(y(m)    -location_fixed(2,1))
          force(:,3)=tmp*(z(n)    -location_fixed(3,1))
!
!  blob-like disturbance (gradient of gaussian)
!
        case('zblob')
          tmp=ampl_ff(i)*exp(-.5*(x(l1:l2)**2+y(m)**2+z(n)**2)/radius_ff**2)
          force(:,1)=0.
          force(:,2)=0.
          force(:,3)=tmp
!
!  blob-like vertical field patch
!
        case('vert_field_blob')
          tmp=x(l1:l2)**2+y(m)**2+z(n)**2
          force(:,1)=0.
          where (tmp<=1.)
            force(:,2)=ampl_ff(i)*x(l1:l2)
          elsewhere
            force(:,2)=0.
          endwhere
          force(:,3)=0.
!
!  KS-flow
!
      case ('KS')
        force=0.
        do modeN=1,KS_modes  ! sum over KS_modes modes
          kdotxwt=KS_k(1,modeN)*x(l1:l2)+(KS_k(2,modeN)*y(m)+KS_k(3,modeN)*z(n))+KS_omega(modeN)*t
          cos_kdotxwt=cos(kdotxwt) ;  sin_kdotxwt=sin(kdotxwt)
          force(:,1) = force(:,1) + cos_kdotxwt*KS_A(1,modeN) + sin_kdotxwt*KS_B(1,modeN)
          force(:,2) = force(:,2) + cos_kdotxwt*KS_A(2,modeN) + sin_kdotxwt*KS_B(2,modeN)
          force(:,3) = force(:,3) + cos_kdotxwt*KS_A(3,modeN) + sin_kdotxwt*KS_B(3,modeN)
        enddo
        force=ampl_ff(i)*force
!
!  fz=Exp(-(kf_cont*r)^2), where r=sqrt(x^2+y^2)
!
      case ('exp(-x2-y2)')
        force(:,1) = 0.
        force(:,2) = 0.
        force(:,3) = ampl_ff(i)*expmk2x2(l1:l2,i)*expmk2y2(m,i)
!
!  The following 2 entries were taken verbatim from
!  src/experimental/forcing_cont_sample.f90
!
      case('xz')
        force(:,1)=0.
        force(:,2)=profx_ampl*profz_ampl(n)
        force(:,3)=0.

      case ('1-(4/3)(1-r^2/4)*r^2')
        force(:,1)=0.
        force(:,2)=0.
        force(:,3)=profx_ampl
!
!  Tidal forcing in a convective box; a reference will be added.
!  phi_tidal is used as the colatitude of the box and has to be
!  the same as theta of Omega in the hydro module.
!
      case ('tidal')
        force(:,1) = ampl_ff(i) * ( &
            sinx(l1:l2,i) * cos(2.*phi_tidal) * cos(omega_ff*t) + &
            siny(m,i)     * 2.*cos(phi_tidal) * sin(omega_ff*t) + &
            z(n)          * 0.5*sin(2.*phi_tidal) * cos(omega_ff*t) )
        force(:,2) = ampl_ff(i) * ( &
            sinx(l1:l2,i) * cos(phi_tidal) * sin(omega_ff*t) + &
            siny(m,i)     * (-2.) * cos(omega_ff*t) + &
            z(n)          * sin(phi_tidal) * sin(omega_ff*t) )
        force(:,3) = ampl_ff(i) * ( &
            sinx(l1:l2,i) * sin(2.*phi_tidal) * cos(omega_ff*t) + &
            siny(m,i)     * 2.*sin(phi_tidal) * sin(omega_ff*t) + &
            z(n)          * sin(phi_tidal)**2 * cos(omega_ff*t) )
!
!
!  possibility of putting zero, e.g., for purely magnetic forcings
!
      case('zero')
        force=0.
!
!   Read forcing profile from file. Currently can be used only for one variable
!   (e.g. either uu or aa).
!
      case('from_file')
        force(:,1) = fcont_from_file(1,l1-nghost+ipx*nx:l2-nghost+ipx*nx,m-nghost+ipy*ny,n-nghost+ipz*nz)
        force(:,2) = fcont_from_file(2,l1-nghost+ipx*nx:l2-nghost+ipx*nx,m-nghost+ipy*ny,n-nghost+ipz*nz)
        force(:,3) = fcont_from_file(3,l1-nghost+ipx*nx:l2-nghost+ipx*nx,m-nghost+ipy*ny,n-nghost+ipz*nz)
        force=ampl_ff(i)*force
!
!  nothing 
!
      case ('nothing')
        call warning('forcing_cont',"iforcing_cont='nothing'")
!
      case default
        call fatal_error('forcing_cont','no such iforcing_cont: '//trim(iforcing_cont(i)))
      endselect
!
    endsubroutine forcing_cont
!***********************************************************************
    subroutine calc_diagnostics_forcing(p)
!
!  add a continuous forcing term (used to be in hydro.f90)
!
!  17-sep-06/axel: coded
!
      use Diagnostics
      use Sub
!
      type (pencil_case), intent(IN) :: p

      real, dimension (nx,3) :: fxb
      real, dimension (nx) :: tmp
!
!  diagnostics
!
      if (ldiagnos) then
        if (idiag_rufm/=0) then
          call dot_mn(p%uu,p%fcont(:,:,1),tmp)
          call sum_mn_name(p%rho*tmp,idiag_rufm)
        endif
!
        if (idiag_rufint/=0) then
          call dot_mn(p%uu,p%fcont(:,:,1),tmp)
          call integrate_mn_name(p%rho*tmp,idiag_rufint)
        endif
!
        if (idiag_ufm/=0) then
          call dot_mn(p%uu,p%fcont(:,:,1),tmp)
          call sum_mn_name(tmp,idiag_ufm)
        endif
!
        if (idiag_ofm/=0) then
          call dot_mn(p%oo,p%fcont(:,:,1),tmp)
          call sum_mn_name(tmp,idiag_ofm)
        endif
!
        if (idiag_qfm/=0) then
          call dot_mn(p%curlo,p%fcont(:,:,1),tmp)
          call sum_mn_name(tmp,idiag_qfm)
        endif
!
        if (lmagnetic) then
!         if (idiag_fxbxm/=0.or.idiag_fxbym/=0.or.idiag_fxbzm/=0) then
!           call cross(p%fcont(:,:,2),p%bb,fxb)
!           call sum_mn_name(fxb(:,1),idiag_fxbxm)
!           call sum_mn_name(fxb(:,2),idiag_fxbym)
!           call sum_mn_name(fxb(:,3),idiag_fxbzm)
!         endif
        endif
      endif
!
    endsubroutine calc_diagnostics_forcing
!***********************************************************************
    subroutine read_forcing_run_pars(iostat)
!
      use File_io, only: parallel_unit
!
      integer, intent(out) :: iostat
!
      read(parallel_unit, NML=forcing_run_pars, IOSTAT=iostat)
!
    endsubroutine read_forcing_run_pars
!***********************************************************************
    subroutine write_forcing_run_pars(unit)
!
      integer, intent(in) :: unit
!
      write(unit, NML=forcing_run_pars)
!
    endsubroutine write_forcing_run_pars
!***********************************************************************
    subroutine input_persist_forcing_id(id,done)
!
!  Read in the stored time of the next SNI
!
!  21-dec-05/tony: coded
!  13-Dec-2011/Bourdin.KIS: reworked
!
      use IO, only: read_persist
!
      integer :: id
      logical :: done
!
      select case (id)
        case (id_record_FORCING_LOCATION)
          done=read_persist ('FORCING_LOCATION', location)
          if (lroot.and..not.done) print *, 'input_persist_forcing: ', location
        case (id_record_FORCING_TSFORCE)
          done=read_persist ('FORCING_TSFORCE', tsforce)
          if (lroot.and..not.done) print *, 'input_persist_forcing: ', tsforce
        case (id_record_FORCING_TORUS)
          done=read_persist ('FORCING_TORUS', torus)
      endselect
!
    endsubroutine input_persist_forcing_id
!***********************************************************************
    subroutine input_persist_forcing
!
!  Read in the persistent forcing variables.
!
!  11-Oct-2019/PABourdin: coded
!
      use IO, only: read_persist
!
      logical :: error
!
      error = read_persist ('FORCING_LOCATION', location)
      if (lroot .and. .not. error) print *, 'input_persist_forcing: location: ', location
!print *, 'input_persist_forcing: location: ', location, error
!
      error = read_persist ('FORCING_TSFORCE', tsforce)
      if (lroot .and. .not. error) print *, 'input_persist_forcing: tsforce: ', tsforce
!
      !error = read_persist ('FORCING_TORUS', torus)
      !if (lroot .and. .not. error) print *, 'input_persist_forcing: torus: ', torus
!
    endsubroutine input_persist_forcing
!***********************************************************************
    logical function output_persistent_forcing()
!
!  This is used, for example, for forcing functions with temporal
!  memory, such as in the paper by Mee & Brandenburg (2006, MNRAS)
!
!  21-dec-05/tony: coded
!  13-Dec-2011/Bourdin.KIS: reworked
!
      use IO, only: write_persist
!
      if (ip<=6.and.lroot .and. (tsforce>=0.)) print *,'output_persistent_forcing: ',location,tsforce
!
!  write details
!
      output_persistent_forcing = .true.
!
      if (write_persist ('FORCING_LOCATION', id_record_FORCING_LOCATION, location)) return
      if (write_persist ('FORCING_TSFORCE', id_record_FORCING_TSFORCE, tsforce)) return
      if (torus%thick>0.) then
        if (write_persist ('FORCING_TORUS', id_record_FORCING_TORUS, torus)) return
      endif
!
      output_persistent_forcing = .false.
!
    endfunction output_persistent_forcing
!***********************************************************************
    subroutine rprint_forcing(lreset,lwrite)
!
!  reads and registers print parameters relevant for hydro part
!
!  26-jan-04/axel: coded
!
      use Diagnostics
!
      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
        idiag_rufm=0; idiag_rufint=0; idiag_ufm=0; idiag_ofm=0; idiag_qfm=0; idiag_ffm=0
        idiag_ruxfxm=0; idiag_ruyfym=0; idiag_ruzfzm=0
        idiag_ruxfym=0; idiag_ruyfxm=0
!       idiag_fxbxm=0; idiag_fxbym=0; idiag_fxbzm=0
      endif
!
!  iname runs through all possible names that may be listed in print.in
!
      if (lroot.and.ip<14) print*,'rprint_forcing: run through parse list'
      do iname=1,nname
        call parse_name(iname,cname(iname),cform(iname),'rufm',idiag_rufm)
        call parse_name(iname,cname(iname),cform(iname),'rufint',idiag_rufint)
        call parse_name(iname,cname(iname),cform(iname),'ruxfxm',idiag_ruxfxm)
        call parse_name(iname,cname(iname),cform(iname),'ruxfym',idiag_ruxfym)
        call parse_name(iname,cname(iname),cform(iname),'ruyfxm',idiag_ruyfxm)
        call parse_name(iname,cname(iname),cform(iname),'ruyfym',idiag_ruyfym)
        call parse_name(iname,cname(iname),cform(iname),'ruzfzm',idiag_ruzfzm)
        call parse_name(iname,cname(iname),cform(iname),'ufm',idiag_ufm)
        call parse_name(iname,cname(iname),cform(iname),'ofm',idiag_ofm)
        call parse_name(iname,cname(iname),cform(iname),'qfm',idiag_qfm)
        call parse_name(iname,cname(iname),cform(iname),'ffm',idiag_ffm)
!       call parse_name(iname,cname(iname),cform(iname),'fxbxm',idiag_fxbxm)
!       call parse_name(iname,cname(iname),cform(iname),'fxbym',idiag_fxbym)
!       call parse_name(iname,cname(iname),cform(iname),'fxbzm',idiag_fxbzm)
      enddo
!
!  write column where which forcing variable is stored
!
      if (lwr) then
!
      endif
!
    endsubroutine rprint_forcing
!***********************************************************************
    subroutine forcing_clean_up
!
!  deallocate the variables needed for forcing
!
!   12-aug-09/dhruba: coded
!
      if (allocated(psif)) deallocate(psif)
      if (allocated(cklist)) deallocate(cklist)
      if (allocated(Zpsi_list)) deallocate(Zpsi_list)
      if (allocated(RYlm_list)) deallocate(RYlm_list)
      if (allocated(IYlm_list)) deallocate(IYlm_list)
      if (allocated(RYlm_rot)) deallocate(RYlm_rot)
      if (allocated(IYlm_rot)) deallocate(IYlm_rot)
      if (allocated(random2d_kmodes)) deallocate(random2d_kmodes)
      if (allocated(fcont_from_file)) deallocate(fcont_from_file)
!
    endsubroutine forcing_clean_up
!***********************************************************************
    subroutine pushpars2c(p_par)

    use Syscalls, only: copy_addr

    integer, parameter :: n_pars=9
    integer(KIND=ikind8), dimension(n_pars) :: p_par

    call copy_addr(k1_ff,p_par(1))
    call copy_addr(tforce_stop,p_par(2))
    call copy_addr(iforcing_zsym,p_par(3)) ! int
    call copy_addr(profx_ampl,p_par(4))  ! (nx)
    call copy_addr(profy_ampl,p_par(5))  ! (my)
    call copy_addr(profz_ampl,p_par(6))  ! (mz)
    call copy_addr(profx_hel,p_par(7))   ! (nx)
    call copy_addr(profy_hel,p_par(8))   ! (my)
    call copy_addr(profz_hel,p_par(9))   ! (mz)

    endsubroutine pushpars2c
!*******************************************************************
endmodule Forcing
