diff --git a/arcos_f90/src/cdr_advect_diffu_vec.f90 b/arcos_f90/src/cdr_advect_diffu_vec.f90
new file mode 100644
index 0000000000000000000000000000000000000000..691d4b3470a69a2b0d22070d385a99bbe732e3ec
--- /dev/null
+++ b/arcos_f90/src/cdr_advect_diffu_vec.f90
@@ -0,0 +1,414 @@
+!> @file cdr_advect_diffu_vec.f90
+!! @brief This code replaces the original c code of cdr_advect_diffu in file
+!! cdr.c
+!!
+!! @author Margreet Nool
+!!
+!! The arrays er_array and ez_array are in C stored as
+!! 1D arrays. In the wrapper routine cdr_advect_diffu_wrap.c it is
+!! checked of their sizes, er_len, ez_len are equal to rmax-by-zmax.
+!!
+!! Note that these are NOT the r0, z0 we received. The idea here is that
+!! RZ(array, r0, z0) will produce the correct result, where the allocated
+!! array looks like (in the picture, r0, z0 ARE the received ones).
+!! Note: in the picture there is only one buffer cell, while we actually
+!! allocate two.
+!!
+!! +--------------+--------------+...+--------------+--------------+
+!! |r0 - 1, z0 - 1| | | | |
+!! +--------------+--------------+...+--------------+--------------+
+!! | | r0, z0 | | | |
+!! +--------------+--------------+...+--------------+--------------+
+!! ... |...| ...
+!! +--------------+--------------+...+--------------+--------------+
+!! | | | |r1 - 1, z1 - 1| |
+!! +--------------+--------------+...+--------------+--------------+
+!! | | | | | r1, z1 |
+!! +--------------+--------------+...+--------------+--------------+
+!!
+!! but the rows z0 +/- 1 and the columns r0 +/- 1 do not belong to the
+!! physical space and are used only to specify boundary conditions.
+!!
+!! In the C routines r0, and z0 are zero, and the calculations start
+!! at the third element (r0-2, r0-1, r0, ...., r1, r1+1), and
+!! (z0-2, z0-1, r0, ...., z1, z1+1), respectively.
+!! rmax = r1 - r0 + 4 and zmax = z1 - z0 + 4.
+!!-----------------------------------------------
+
+! subroutine cdr_advect_vec
+! subroutine mat_efield_r
+! subroutine mat_efield_z
+! subroutine mat_cdr_f_ad_r
+! subroutine mat_cdr_f_ad_z
+! subroutine make_vecs_er
+
+ subroutine cdr_advect_vec(mass,charge,dens_array,d_dens_array, &
+ er_array,ez_array,diffusion_coeff,dr,dz, &
+ sprite_module,r0,r1,z0,z1)
+ implicit none
+! ..
+! .. Parameters ..
+ double precision,parameter :: ONE = 1.0D0
+! ..
+! .. Dummy Variables ..
+ logical,intent(in) :: sprite_module
+ integer,intent(in) :: r0,r1,z0,z1
+ double precision,intent(in) :: charge,dr,dz,mass,diffusion_coeff
+! ..
+! .. Dummy Arrays ..
+ double precision,dimension(r0-2:r1+1,z0-2:z1+1),intent(in) :: &
+ dens_array,er_array,ez_array
+ double precision,dimension(r0-2:r1+1,z0-2:z1+1),intent(out) :: &
+ d_dens_array
+! ..
+! .. Local Variables ..
+ integer :: ir,ishift
+ double precision :: cm0,cm,gdz_inv,d_iso0,d_iso
+! ..
+! .. Local arrays ..
+ double precision,dimension(r0:r1-1) :: er_min
+ double precision,dimension(r0+1:r1) :: er
+ double precision,dimension(z1-z0) :: d_iso_r
+ double precision,dimension(z1-z0+1) :: d_iso_z
+ double precision,dimension(r1-r0+1,z1-z0 ) :: mat_r
+ double precision,dimension(r1-r0 ,z1-z0+1) :: mat_z
+ double precision,dimension(r1-r0+1,z1-z0 ) :: efield_r
+ double precision,dimension(r0:r1-1,z0-1:z1-1) :: efield_z
+! ..
+! .. External subroutines ..
+ external :: mat_efield_r,mat_efield_z,mat_cdr_f_ad_r,mat_cdr_f_ad_z
+! ..
+
+ if (mass <= 0) return
+
+ !
+ d_iso0 = diffusion_coeff / mass
+ !
+
+ gdz_inv = ONE / dz
+
+ cm0 = charge / mass
+
+ if (.not.sprite_module) then
+ cm = cm0
+ d_iso = d_iso0
+ end if
+
+!-----------------------------------------------
+! Increase the values r0, r1, z0 and z1 with 3:
+! because of 2 buffer cells at the left side and to compensate the C 0
+!-----------------------------------------------
+ ishift = 3
+
+ ! Compute the electric field in radial direction
+ call mat_efield_r(er_array(r0-1:r1-1,z0:z1-1), &
+ efield_r,d_iso_r,cm0,d_iso0, &
+ r0,r1,dz,z0,z1,ishift,dr,sprite_module)
+ ! Calculate the radial advection from left and right
+ call mat_cdr_f_ad_r(dens_array(r0-2:r1+1,z0:z1-1), &
+ efield_r(1:r1-r0+1,1:z1-z0),d_iso_r, &
+ mat_r,r0,r1,z0,z1)
+ ! Add the radial advection
+ call make_vecs_er (er_min,er,r0,r1,dr)
+ do ir = r0,r1-1
+ d_dens_array(ir,z0:z1-1) = d_dens_array(ir,z0:z1-1) + &
+ er_min(ir) * mat_r(ir-r0+1,1:z1-z0) - &
+ er(ir+1) * mat_r(ir-r0+2,1:z1-z0)
+ end do
+
+ ! Compute the electric field in axial direction
+ call mat_efield_z(ez_array(r0 :r1-1,z0-1:z1-1), &
+ efield_z(r0:r1-1,z0-1:z1-1),d_iso_z,cm0,d_iso0, &
+ r0,r1,dz,z0,z1,ishift,dz,sprite_module)
+ ! Calculate the axial advection from the north and south
+ call mat_cdr_f_ad_z(dens_array(r0 :r1-1,z0-2:z1+1), &
+ efield_z(r0:r1-1,z0-1:z1-1),d_iso_z, &
+ mat_z,r0,r1,z0,z1)
+ ! Add the axial advection
+ d_dens_array(r0:r1-1,z0:z1-1) = d_dens_array(r0:r1-1,z0:z1-1) + &
+ gdz_inv * (mat_z(1:,1:z1-z0) - &
+ mat_z(1:,2:z1-z0+1))
+
+ end subroutine cdr_advect_vec
+
+ subroutine mat_efield_r(er_array,efield,d_iso,cm0,d_iso0, &
+ r0,r1,dz,z0,z1,z_shift,dx,sprite_module)
+! ..
+ implicit none
+! ..
+! .. Parameters ..
+ double precision,parameter :: HALF = 0.5D0,ONE = 1.0D0
+! ..
+! .. Dummy Variables ..
+ logical,intent(in) :: sprite_module
+ integer,intent(in) :: r0,r1,z0,z1,z_shift
+ double precision,intent(in) :: cm0,d_iso0,dz,dx
+! ..
+! .. Dummy Arrays ..
+ double precision,dimension(r1-r0+1,z1-z0),intent(in) :: er_array
+ double precision,dimension(r1-r0+1,z1-z0),intent(out) :: efield
+ double precision,dimension(z1-z0),intent(out) :: d_iso
+! ..
+! .. Local Variables ..
+ integer :: iz
+! ..
+! .. Local arrays ..
+ double precision,dimension(z1-z0) :: back_dens_inv,cm
+!-----------------------------------------------
+
+ ! If the densities are varying...
+ if (sprite_module) then
+ do iz = z0,z1
+ ! spr_density_at(iz) = (iz-z_shift + half) * dz
+ ! back_dens_inv = ONE / spr_density_at
+ back_dens_inv(iz-z0+1) = ONE / ((iz-z_shift + half) * dz)
+ end do
+ cm = cm0 * back_dens_inv
+ d_iso = d_iso0 * back_dens_inv
+ else
+ back_dens_inv = ONE
+ cm = cm0
+ d_iso = d_iso0
+ end if
+
+ d_iso = d_iso / dx
+
+ ! efield = cm * er_array
+ do iz = 1,z1-z0
+ efield(:,iz) = cm(iz) * er_array(:,iz)
+ end do
+
+ end subroutine mat_efield_r
+
+ subroutine mat_efield_z (ez_array,efield,d_iso,cm0,d_iso0, &
+ r0,r1,dz,z0,z1,z_shift,dx,sprite_module)
+! ..
+ implicit none
+! ..
+! .. Parameters ..
+ double precision,parameter :: HALF = 0.5D0,ONE = 1.0D0
+! ..
+! .. Dummy Variables ..
+ logical,intent(in) :: sprite_module
+ integer,intent(in) :: r0,r1,z0,z1,z_shift
+ double precision,intent(in) :: cm0,d_iso0,dz,dx
+! ..
+! .. Dummy Arrays ..
+ double precision,dimension(r0:r1-1,z0-1:z1-1),intent(in) :: ez_array
+ double precision,dimension(r0:r1-1,z0-1:z1-1),intent(out) :: efield
+ double precision,dimension(z0-1:z1-1),intent(out) :: d_iso
+! ..
+! .. Local Variables ..
+ integer :: iz
+! ..
+! .. Local arrays ..
+ double precision,dimension(z0-1:z1-1) :: back_dens_inv,cm
+! ..
+!-----------------------------------------------
+
+ ! If the densities are varying...
+ if (sprite_module) then
+ do iz = z0-1,z1+1
+ ! spr_density_at(iz) = (iz-z_shift + half) * dz
+ ! back_dens_inv = ONE / spr_density_at
+ back_dens_inv(iz) = ONE / ((iz-z_shift - half) * dz)
+ end do
+ cm = cm0 * back_dens_inv
+ d_iso = d_iso0 * back_dens_inv
+ else
+ back_dens_inv = ONE
+ cm = cm0
+ d_iso = d_iso0
+ end if
+
+ d_iso = d_iso / dx
+
+ ! efield = cm * ez_array
+ do iz = z0-1,z1-1
+ efield(:,iz) = cm(iz) * ez_array(:,iz)
+ end do
+
+ end subroutine mat_efield_z
+
+ subroutine mat_cdr_f_ad_r(data_array,efield,d_iso,mat_r,r0,r1,z0,z1)
+! ..
+ implicit none
+! ..
+! .. Parameters ..
+ double precision,parameter :: ZERO = 0.0D0,HALF = 0.5D0, &
+ ONE = 1.0D0,THREE = 3.0D0, &
+ verysmall = 1.0e-20
+! ..
+! .. Dummy Variables ..
+ integer,intent(in) :: r0,r1,z0,z1
+! ..
+! .. Dummy Arrays ..
+ double precision,dimension(z1-z0),intent(in) :: d_iso
+ double precision,dimension(r1-r0+1,z1-z0),intent(in) :: efield
+ double precision,dimension(r1-r0+4,z1-z0),intent(in) :: data_array
+ double precision,dimension(r1-r0+1,z1-z0),intent(out) :: mat_r
+! ..
+! .. Local Variables ..
+ integer :: iz,length
+! ..
+! .. Local arrays ..
+ double precision,dimension(r1-r0+1,z1-z0) :: aux,psi_p
+ double precision,dimension(r1-r0+3,z1-z0) :: mat_diff
+! ..
+!-----------------------------------------------
+ length = r1 - r0 + 1
+ mat_diff(1:length+2,:) = data_array(2:length+3,:) - &
+ data_array(1:length+2,:)
+
+ where (abs(mat_diff(2:length+1,:)) > verysmall)
+ where (efield>=ZERO)
+ ! p_min =data_array(ir-k1,iz-k2)
+ ! p =data_array(ir,iz)
+ ! p_plus =data_array(ir+k1,iz+k2)
+
+ psi_p = mat_diff(1:length,:) / mat_diff(2:length+1,:)
+ elsewhere (efield= 4)
+ psi_p = ONE
+ elsewhere (psi_p >= 0.4)
+ psi_p = (ONE + HALF * psi_p) / THREE
+ end where
+ elsewhere ! abs(mat_diff(2:length+1,:)) < verysmall
+ psi_p = -ONE
+ end where
+
+ ! sigmasigma is (sigma_{i+1,j} - sigma_{i,j})*/
+ ! sigmasigma = p - p_plus
+ aux = psi_p * mat_diff(2:length+1,:)
+ ! aux = p + psi_p * sigmasigma
+ where (efield>=ZERO)
+ aux = data_array(2:length+1,:) + aux
+ elsewhere
+ aux = data_array(3:length+2,:) - aux
+ end where
+
+ ! aux = efield * (p + psi_p * sigmasigma)
+ aux = aux * efield
+ do iz = 1,z1-z0
+ ! mat_r = d_iso * sigmasigma
+ mat_r(1:length,iz) = d_iso(iz) * mat_diff(2:length+1,iz)
+ end do
+ ! mat_r = efield * (p + psi_p * sigmasigma) + d_iso * sigmasigma
+ mat_r = aux - mat_r
+
+ end subroutine mat_cdr_f_ad_r
+
+ subroutine mat_cdr_f_ad_z(data_array,efield,d_iso,mat_z,r0,r1,z0,z1)
+! ..
+ implicit none
+! ..
+! .. Parameters ..
+ double precision,parameter :: ZERO = 0.0D0,HALF = 0.5D0, &
+ ONE = 1.0D0,THREE = 3.0D0, &
+ verysmall = 1.0e-20
+! ..
+! .. Dummy Variables ..
+ integer,intent(in) :: r0,r1,z0,z1
+! ..
+! .. Dummy Arrays ..
+ double precision,dimension(z1-z0+1),intent(in) :: d_iso
+ double precision,dimension(r0:r1-1,z0-1:z1-1),intent(in) :: efield
+ double precision,dimension(r1-r0,z1-z0+4),intent(in) :: data_array
+ double precision,dimension(r1-r0,z1-z0+1),intent(out) :: mat_z
+! ..
+! .. Local Variables ..
+ integer :: iz,length
+! ..
+! .. Local Arrays ..
+ double precision,dimension(r1-r0,z1-z0+1) :: aux,psi_p,pij
+ double precision,dimension(r1-r0,z1-z0+3) :: mat_diff
+!-----------------------------------------------
+ length = z1 - z0 + 1
+ mat_diff(:,1:length+2) = data_array(:,2:length+3) - &
+ data_array(:,1:length+2)
+
+ where (abs(mat_diff(:,2:length+1)) > verysmall)
+ ! psi_p = (p - p_min) / (p_plus- p)
+ where (efield>ZERO)
+ ! p_min =data_array(ir-k1,iz-k2)
+ ! p =data_array(ir,iz)
+ ! p_plus =data_array(ir+k1,iz+k2)
+
+ psi_p = mat_diff(:,1:length) / mat_diff(:,2:length+1)
+ elsewhere ! efield= 0.4) .and. (psi_p < 4))
+ pij = (ONE + HALF * psi_p) / THREE
+ elsewhere (psi_p >= 4)
+ pij = ONE
+ else where (psi_p <= ZERO)
+ pij = ZERO
+ end where
+ elsewhere ! abs(mat_diff(:,2:length+1)) < verysmall
+ pij = ONE
+ end where
+
+ ! sigmasigma is (sigma_{i+1,j} - sigma_{i,j})*/
+ ! sigmasigma = p - p_plus
+ aux = pij * mat_diff(:,2:length+1)
+ ! aux = p + pij * sigmasigma
+ where (efield>=0)
+ aux = data_array(:,2:length+1) + aux
+ elsewhere
+ aux = data_array(:,3:length+2) - aux
+ end where
+
+ ! aux = efield * (p + pij * sigmasigma)
+ aux = aux * efield
+ do iz = 1,length
+ ! mat_z = d_iso * sigmasigma
+ mat_z(:,iz) = d_iso(iz) * mat_diff(:,iz+1)
+ end do
+ ! mat_z = efield * (p + pij * sigmasigma) + d_iso * sigmasigma
+ mat_z = aux - mat_z
+
+ end subroutine mat_cdr_f_ad_z
+
+ subroutine make_vecs_er (er_min,er,r0,r1,dr)
+! ..
+ implicit none
+!
+! .. Parameters ..
+ double precision,parameter :: HALF = 0.5D0,ONE = 1.0D0
+! ..
+! .. Dummy Variables ..
+ integer,intent(in) :: r0,r1
+ double precision,intent(in) :: dr
+! ..
+! .. Dummy Arrays ..
+ double precision,dimension(r0:r1-1),intent(out) :: er_min
+ double precision,dimension(r0+1:r1),intent(out) :: er
+! ..
+! .. Local Variables ..
+ integer :: ir
+ double precision :: r_inv
+! ..
+! .............................................................
+ do ir = r0,r1-1
+ r_inv = ONE / ((ir+ HALF) * dr)
+ er_min(ir) = ir * r_inv
+ er(ir+1) = (ir+1) * r_inv
+ end do
+
+ end subroutine make_vecs_er