File diff 000000000000 → d6faa5ffcedf
arcos_f90/src/cdr_advect_diffu_vec.f90
Show inline comments
 
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