|
new file 100644
|
|
|
!> @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
|
|
|
|
|
|
! <ISOTROPIC DIFFUSSION>
|
|
|
d_iso0 = diffusion_coeff / mass
|
|
|
! </ISOTROPIC DIFFUSION>
|
|
|
|
|
|
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<ZERO)
|
|
|
! p =data_array(ir,iz)
|
|
|
! p_plus =data_array(ir+ k1,iz+ k2)
|
|
|
! p_2plus =data_array(ir+2*k1,iz+2*k2)
|
|
|
|
|
|
psi_p = mat_diff(3:length+2,:) / mat_diff(2:length+1,:)
|
|
|
end where
|
|
|
where (psi_p <= ZERO)
|
|
|
! psi_p = psi_MN (psi_p)
|
|
|
psi_p = ZERO
|
|
|
elsewhere (psi_p >= 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=<ZERO
|
|
|
! p =data_array(ir,iz)
|
|
|
! p_plus =data_array(ir+ k1,iz+ k2)
|
|
|
! p_2plus =data_array(ir+2*k1,iz+2*k2)
|
|
|
|
|
|
psi_p = mat_diff(:,3:length+2) / mat_diff(:,2:length+1)
|
|
|
end where
|
|
|
! psi_p = psi_MN (psi_p)
|
|
|
pij = psi_p
|
|
|
where ((psi_p > 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
|