!> @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