! Initial condition (density, magnetic field, velocity)
! for a particular configuration of magnetic carpets.
!
!** 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 :: linitial_condition = .true.
!
!***************************************************************
module InitialCondition
!
use Cparam
use Cdata
use General, only: keep_compiler_quiet
use Mpicomm
use Messages
use Boundcond ! for the core boundary communication
!
implicit none
!
include '../initial_condition.h'
!
! ampl = amplitude of the magnetic field
!
real :: ampl = 1.0
character (len=labellen) :: config='parasitic_polarities'
!
namelist /initial_condition_pars/ &
ampl, config
!
contains
!***********************************************************************
subroutine register_initial_condition()
!
! Configure pre-initialised (i.e. before parameter read) variables
! which should be know to be able to evaluate
!
! 30-april-15/simon: coded
!
endsubroutine register_initial_condition
!***********************************************************************
subroutine initial_condition_uu(f)
!
! Initialize the velocity field.
!
! 07-may-09/wlad: coded
!
real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
call keep_compiler_quiet(f)
!
endsubroutine initial_condition_uu
!***********************************************************************
subroutine initial_condition_lnrho(f)
!
! Initialize logarithmic density. init_lnrho
! will take care of converting it to linear
! density if you use ldensity_nolog
!
! 07-may-09/wlad: coded
!
real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
call keep_compiler_quiet(f)
!
endsubroutine initial_condition_lnrho
!***********************************************************************
subroutine initial_condition_aa(f)
!
! Initialize the magnetic vector potential.
!
! 30-april-15/simon: coded
!
! Magnetic carpet consisting of three spots.
!
! Created 2015-04-30 by Simon Candelaresi (Iomsn)
!
use Mpicomm, only: stop_it
use Poisson
use Sub
!
real, dimension (mx,my,mz,mfarray) :: f
!
! The next variables are used for the uncurling.
integer :: l, j, ju, k
real, dimension (nx,ny,nz,3) :: jj, tmpJ ! This is phi for poisson.f90
!
! clear the magnetic field to zero
f(:,:,:,iax:iaz) = 0.
!
if (config == 'homogeneous') then
do m = 1, my, 1
do l = 1, mx, 1
f(l,m,:,iax) = -ampl*y(m)/2
f(l,m,:,iay) = ampl*x(l)/2
enddo
enddo
endif
if (config == 'parasitic_polarities') then
do n = 1, mz, 1
do m = 1, my, 1
do l = 1, mx, 1
f(l,m,n,iax) = ampl * (&
2. * ((x(l) - 10.) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2)&
** (-3.0 / 2.) * y(m) + 2. * ((x(l) + 8.) ** 2 + y(m) **&
2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * y(m) + 2. * ((x(l)&
+ 10.) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
/ 2.) * y(m) + 2. * ((x(l) + 6.) ** 2 + y(m) ** 2 + (z(n)&
+ 0.85) ** 2) ** (-3.0 / 2.) * y(m) + 2. * (x(l) ** 2 + y(m)&
** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * y(m) + 2. *&
((x(l) + 2.) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
/ 2.) * y(m) + 2. * ((x(l) - 2.) ** 2 + y(m) ** 2 + (z(n)&
+ 0.85) ** 2) ** (-3.0 / 2.) * y(m) + 2. * ((x(l) - 8.)&
** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 /&
2.) * (y(m) - 8.) + 2. * ((x(l) - 6.) ** 2 + (y(m) - 8.)&
** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) - 8.) +&
2. * ((x(l) - 10.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85)&
** 2) ** (-3.0 / 2.) * (y(m) - 8.) + 2. * ((x(l) - 8.) **&
2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.)&
* (y(m) + 8.) + 2. * ((x(l) - 6.) ** 2 + (y(m) + 8.) ** 2&
+ (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) + 8.) + 2. *&
((x(l) - 10.) ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) **&
2) ** (-3.0 / 2.) * (y(m) + 8.) + 2. * ((x(l) + 8.) ** 2 +&
(y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m)&
- 8.) + 2. * ((x(l) + 10.) ** 2 + (y(m) - 8.) ** 2 +&
(z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) - 8.) + 2. * ((x(l)&
+ 6.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) **&
(-3.0 / 2.) * (y(m) - 8.) + 2. * ((x(l) + 8.) ** 2 + (y(m)&
+ 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m)&
+ 8.) + 2. * ((x(l) + 10.) ** 2 + (y(m) + 8.) ** 2 + (z(n)&
+ 0.85) ** 2) ** (-3.0 / 2.) * (y(m) + 8.) + 2. * ((x(l)&
+ 6.) ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
/ 2.) * (y(m) + 8.) + 2. * (x(l) ** 2 + (y(m) + 8.) **&
2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) + 8.) + 2.&
* ((x(l) + 2.) ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) **&
2) ** (-3.0 / 2.) * (y(m) + 8.) + 2. * ((x(l) - 2.) ** 2&
+ (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) *&
(y(m) + 8.) + 2. * (x(l) ** 2 + (y(m) - 8.) ** 2 + (z(n) +&
0.85) ** 2) ** (-3.0 / 2.) * (y(m) - 8.) + 2. * ((x(l) +&
2.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
/ 2.) * (y(m) - 8.) + 2. * ((x(l) - 2.) ** 2 + (y(m) - 8.)&
** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) - 8.)&
+ 2. * ((x(l) - 8.) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2)&
** (-3.0 / 2.) * y(m) + 2. * ((x(l) - 6.) ** 2 + y(m) **&
2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * y(m))
f(l,m,n,iay) = ampl * (&
-2. * ((x(l) - 10.) ** 2 + y(m) ** 2 + (z(n) + 0.85) **&
2) ** (-3.0 / 2.) * (x(l) - 10.) - 2. * ((x(l) + 8.) ** 2&
+ y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) +&
8.) - 2. * ((x(l) + 10.) ** 2 + y(m) ** 2 + (z(n) + 0.85)&
** 2) ** (-3.0 / 2.) * (x(l) + 10.) - 2. * ((x(l) + 6.) **&
2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l)&
+ 6.) - 2. * (x(l) ** 2 + y(m) ** 2 + (z(n) + 0.85) **&
2) ** (-3.0 / 2.) * x(l) - 2. * ((x(l) + 2.) ** 2 + y(m) **&
2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) + 2.) - 2.&
* ((x(l) - 2.) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2) **&
(-3.0 / 2.) * (x(l) - 2.) - 2. * ((x(l) - 8.) ** 2 + (y(m)&
- 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l)&
- 8.) - 2. * ((x(l) - 6.) ** 2 + (y(m) - 8.) ** 2 + (z(n)&
+ 0.85) ** 2) ** (-3.0 / 2.) * (x(l) - 6.) - 2. * ((x(l) -&
10.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
/ 2.) * (x(l) - 10.) - 2. * ((x(l) - 8.) ** 2 + (y(m) +&
8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) - 8.)&
- 2. * ((x(l) - 6.) ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85)&
** 2) ** (-3.0 / 2.) * (x(l) - 6.) - 2. * ((x(l) - 10.)&
** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
/ 2.) * (x(l) - 10.) - 2. * ((x(l) + 8.) ** 2 + (y(m) - 8.)&
** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) + 8.)&
- 2. * ((x(l) + 10.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85)&
** 2) ** (-3.0 / 2.) * (x(l) + 10.) - 2. * ((x(l) + 6.)&
** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 /&
2.) * (x(l) + 6.) - 2. * ((x(l) + 8.) ** 2 + (y(m) + 8.) **&
2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) - 2.&
* ((x(l) + 10.) ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85)&
** 2) ** (-3.0 / 2.) * (x(l) + 10.) - 2. * ((x(l) + 6.) **&
2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.)&
* (x(l) + 6.) - 2. * (x(l) ** 2 + (y(m) + 8.) ** 2 + (z(n)&
+ 0.85) ** 2) ** (-3.0 / 2.) * x(l) - 2. * ((x(l) + 2.)&
** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.)&
* (x(l) + 2.) - 2. * ((x(l) - 2.) ** 2 + (y(m) + 8.) **&
2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) - 2.) - 2.&
* (x(l) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) **&
(-3.0 / 2.) * x(l) - 2. * ((x(l) + 2.) ** 2 + (y(m) - 8.)&
** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) + 2.) -&
2. * ((x(l) - 2.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85)&
** 2) ** (-3.0 / 2.) * (x(l) - 2.) - 2. * ((x(l) - 8.) **&
2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l)&
- 8.) - 2. * ((x(l) - 6.) ** 2 + y(m) ** 2 + (z(n) + 0.85)&
** 2) ** (-3.0 / 2.) * (x(l) - 6.) + x(l))
f(l,m,n,iaz) = 0
enddo
enddo
enddo
elseif (config == 'single_parasitic_polarity') then
do n = 1, mz, 1
do m = 1, my, 1
do l = 1, mx, 1
f(l,m,n,iax) = ampl * (&
2. * (x(l) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
/ 2.) * y(m) + 2. * ((x(l) + 8.) ** 2 + y(m) ** 2 + (z(n)&
+ 0.85) ** 2) ** (-3.0 / 2.) * y(m) + 2. * ((x(l) - 8.)&
** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * y(m)&
+ 2. * (x(l) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) **&
2) ** (-3.0 / 2.) * (y(m) - 8.) + 2. * (x(l) ** 2 + (y(m)&
+ 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m)&
+ 8.) + 2. * ((x(l) + 8.) ** 2 + (y(m) + 8.) ** 2 + (z(n)&
+ 0.85) ** 2) ** (-3.0 / 2.) * (y(m) + 8.) + 2. * ((x(l) +&
8.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
/ 2.) * (y(m) - 8.) + 2. * ((x(l) - 8.) ** 2 + (y(m) + 8.)&
** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (y(m) + 8.)&
+ 2. * ((x(l) - 8.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85)&
** 2) ** (-3.0 / 2.) * (y(m) - 8.))
f(l,m,n,iay) = ampl * (&
x(l) - 2. * (x(l) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2)&
** (-3.0 / 2.) * x(l) - 2. * ((x(l) + 8.) ** 2 + y(m) **&
2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) - 2.&
* ((x(l) - 8.) ** 2 + y(m) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
/ 2.) * (x(l) - 8.) - 2. * (x(l) ** 2 + (y(m) - 8.) **&
2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * x(l) - 2. * (x(l)&
** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0&
/ 2.) * x(l) - 2. * ((x(l) + 8.) ** 2 + (y(m) + 8.) ** 2 +&
(z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) - 2. * ((x(l)&
+ 8.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.85) ** 2)&
** (-3.0 / 2.) * (x(l) + 8.) - 2. * ((x(l) - 8.) ** 2 + (y(m)&
+ 8.) ** 2 + (z(n) + 0.85) ** 2) ** (-3.0 / 2.) * (x(l)&
- 8.) - 2. * ((x(l) - 8.) ** 2 + (y(m) - 8.) ** 2 + (z(n)&
+ 0.85) ** 2) ** (-3.0 / 2.) * (x(l) - 8.))
f(l,m,n,iaz) = 0
enddo
enddo
enddo
elseif (config == 'dominant_polarities') then
do n = 1, mz, 1
do m = 1, my, 1
do l = 1, mx, 1
f(l,m,n,iax) = ampl * (&
-0.3 * ((x(l) - 105) ** 2 + ((y(m) + 8) ** 2) + (z(n) +&
0.3) ** 2) ** (-3.0 / 2.) * (y(m) + 8) - 0.3 * ((x(l) - 8.)&
** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 /&
2.) * (y(m) + 8) - 0.3 * ((x(l) + 8.) ** 2 + ((y(m) - 8)&
** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m) - 8) - 0.3&
* ((x(l) + 105) ** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3)&
** 2) ** (-3.0 / 2.) * (y(m) - 8) - 0.3 * ((x(l) + 5.5) **&
2 + ((y(m) - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.)&
* (y(m) - 8) - 0.3 * ((x(l) + 8.) ** 2 + ((y(m) + 8) **&
2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m) + 8) - 0.3&
* ((x(l) + 105) ** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3) **&
2) ** (-3.0 / 2.) * (y(m) + 8) - 0.3 * ((x(l) + 5.5) ** 2&
+ ((y(m) + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) *&
(y(m) + 8) - 0.3 * (x(l) ** 2 + ((y(m) + 8) ** 2) + (z(n)&
+ 0.3) ** 2) ** (-3.0 / 2.) * (y(m) + 8) - 0.3 * ((x(l) +&
2.5) ** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
/ 2.) * (y(m) + 8) - 0.3 * ((x(l) - 2.5) ** 2 + ((y(m)&
+ 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m) + 8)&
- 0.3 * (x(l) ** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3) **&
2) ** (-3.0 / 2.) * (y(m) - 8) - 0.3 * ((x(l) + 2.5) ** 2&
+ ((y(m) - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) *&
(y(m) - 8) - 0.3 * ((x(l) - 2.5) ** 2 + ((y(m) - 8) ** 2)&
+ (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m) - 8) - 0.3 *&
((x(l) - 8.) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
/ 2.) * (y(m)) - 0.3 * ((x(l) - 5.5) ** 2 + (y(m) ** 2)&
+ (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m)) - 0.3 * ((x(l)&
- 105) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
/ 2.) * (y(m)) - 0.3 * ((x(l) + 8.) ** 2 + (y(m) ** 2) +&
(z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m)) - 0.3 * ((x(l)&
+ 105) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 /&
2.) * (y(m)) - 0.3 * ((x(l) + 5.5) ** 2 + (y(m) ** 2) + (z(n)&
+ 0.3) ** 2) ** (-3.0 / 2.) * (y(m)) - 0.3 * ((x(l) -&
2.5) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.)&
* (y(m)) - 0.3 * (x(l) ** 2 + (y(m) ** 2) + (z(n) + 0.3)&
** 2) ** (-3.0 / 2.) * (y(m)) - 0.3 * ((x(l) + 2.5) ** 2&
+ (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m))&
- 0.3 * ((x(l) - 8.) ** 2 + ((y(m) - 8) ** 2) + (z(n) +&
0.3) ** 2) ** (-3.0 / 2.) * (y(m) - 8) - 0.3 * ((x(l) - 5.5)&
** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
/ 2.) * (y(m) - 8) - 0.3 * ((x(l) - 105) ** 2 + ((y(m) - 8)&
** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (y(m) - 8) -&
0.3 * ((x(l) - 5.5) ** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3)&
** 2) ** (-3.0 / 2.) * (y(m) + 8))
f(l,m,n,iay) = ampl * (&
0.3 * ((x(l) - 8.) ** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3)&
** 2) ** (-3.0 / 2.) * (x(l) - 8.) + 0.3 * ((x(l) - 5.5)&
** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 /&
2.) * (x(l) - 5.5) + 0.3 * ((x(l) - 105) ** 2 + ((y(m) -&
8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 105)&
+ 0.3 * ((x(l) - 8.) ** 2 + ((y(m) + 8) ** 2) + (z(n) +&
0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 8.) + 0.3 * ((x(l) - 5.5)&
** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
/ 2.) * (x(l) - 5.5) + 0.3 * ((x(l) - 105) ** 2 + ((y(m)&
+ 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 105)&
+ 0.3 * ((x(l) + 8.) ** 2 + ((y(m) - 8) ** 2) + (z(n)&
+ 0.3) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) + 0.3 * ((x(l) +&
105) ** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
/ 2.) * (x(l) + 105) + 0.3 * ((x(l) + 5.5) ** 2 + ((y(m)&
- 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) +&
5.5) + 0.3 * ((x(l) + 8.) ** 2 + ((y(m) + 8) ** 2) + (z(n)&
+ 0.3) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) + 0.3 * ((x(l)&
+ 105) ** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
/ 2.) * (x(l) + 105) + 0.3 * ((x(l) + 5.5) ** 2 + ((y(m)&
+ 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l)&
+ 5.5) + 0.3 * (x(l) ** 2 + ((y(m) + 8) ** 2) + (z(n) + 0.3)&
** 2) ** (-3.0 / 2.) * x(l) + 0.3 * ((x(l) + 2.5) ** 2&
+ ((y(m) + 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) *&
(x(l) + 2.5) + 0.3 * ((x(l) - 2.5) ** 2 + ((y(m) + 8) **&
2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 2.5) + 0.5D-1&
* x(l) + 0.3 * (x(l) ** 2 + ((y(m) - 8) ** 2) + (z(n)&
+ 0.3) ** 2) ** (-3.0 / 2.) * x(l) + 0.3 * ((x(l) + 2.5)&
** 2 + ((y(m) - 8) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.)&
* (x(l) + 2.5) + 0.3 * ((x(l) - 2.5) ** 2 + ((y(m) - 8)&
** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 2.5)&
+ 0.3 * ((x(l) - 5.5) ** 2 + (y(m) ** 2) + (z(n) + 0.3) **&
2) ** (-3.0 / 2.) * (x(l) - 5.5) + 0.3 * ((x(l) - 105) **&
2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l)&
- 105) + 0.3 * ((x(l) - 8.) ** 2 + (y(m) ** 2) + (z(n)&
+ 0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 8.) + 0.3 * ((x(l) +&
8.) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.)&
* (x(l) + 8.) + 0.3 * ((x(l) + 105) ** 2 + (y(m) ** 2) +&
(z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) + 105) + 0.3 *&
((x(l) + 5.5) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0&
/ 2.) * (x(l) + 5.5) + 0.3 * (x(l) ** 2 + (y(m) ** 2)&
+ (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * x(l) + 0.3 * ((x(l)&
+ 2.5) ** 2 + (y(m) ** 2) + (z(n) + 0.3) ** 2) ** (-3.0 /&
2.) * (x(l) + 2.5) + 0.3 * ((x(l) - 2.5) ** 2 + (y(m) **&
2) + (z(n) + 0.3) ** 2) ** (-3.0 / 2.) * (x(l) - 2.5))
f(l,m,n,iaz) = 0
enddo
enddo
enddo
elseif (config == 'single_dominant_polarity') then
do n = 1, mz, 1
do m = 1, my, 1
do l = 1, mx, 1
f(l,m,n,iax) = ampl * (&
-1. * (x(l) ** 2 + y(m) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0&
/ 2.) * y(m) - 1. * ((x(l) + 8.) ** 2 + y(m) ** 2 + (z(n)&
+ 0.5) ** 2) ** (-3.0 / 2.) * y(m) - 1. * ((x(l) - 8.)&
** 2 + y(m) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0 / 2.) * y(m)&
- 1. * (x(l) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.5) ** 2)&
** (-3.0 / 2.) * (y(m) - 8.) - 1. * (x(l) ** 2 + (y(m) +&
8.) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0 / 2.) * (y(m) + 8.)&
- 1. * ((x(l) + 8.) ** 2 + (y(m) + 8.) ** 2 + (z(n) + 0.5)&
** 2) ** (-3.0 / 2.) * (y(m) + 8.) - 1. * ((x(l) + 8.)&
** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0 / 2.)&
* (y(m) - 8.) - 1. * ((x(l) - 8.) ** 2 + (y(m) + 8.) **&
2 + (z(n) + 0.5) ** 2) ** (-3.0 / 2.) * (y(m) + 8.) - 1. *&
((x(l) - 8.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.5) ** 2)&
** (-3.0 / 2.) * (y(m) - 8.))
f(l,m,n,iay) = ampl * (&
x(l) + 1. * (x(l) ** 2 + y(m) ** 2 + (z(n) + 0.5) ** 2)&
** (-3.0 / 2.) * x(l) + 1. * ((x(l) + 8.) ** 2 + y(m) ** 2&
+ (z(n) + 0.5) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) + 1. *&
((x(l) - 8.) ** 2 + y(m) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0&
/ 2.) * (x(l) - 8.) + 1. * (x(l) ** 2 + (y(m) - 8.) ** 2&
+ (z(n) + 0.5) ** 2) ** (-3.0 / 2.) * x(l) + 1. * (x(l) **&
2 + (y(m) + 8.) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0 / 2.)&
* x(l) + 1. * ((x(l) + 8.) ** 2 + (y(m) + 8.) ** 2 + (z(n)&
+ 0.5) ** 2) ** (-3.0 / 2.) * (x(l) + 8.) + 1. * ((x(l)&
+ 8.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.5) ** 2) ** (-3.0&
/ 2.) * (x(l) + 8.) + 1. * ((x(l) - 8.) ** 2 + (y(m) + 8.)&
** 2 + (z(n) + 0.5) ** 2) ** (-3.0 / 2.) * (x(l) - 8.)&
+ 1. * ((x(l) - 8.) ** 2 + (y(m) - 8.) ** 2 + (z(n) + 0.5)&
** 2) ** (-3.0 / 2.) * (x(l) - 8.))
f(l,m,n,iaz) = 0
enddo
enddo
enddo
endif
!
endsubroutine initial_condition_aa
!***********************************************************************
subroutine read_initial_condition_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=initial_condition_pars, IOSTAT=iostat)
!
endsubroutine read_initial_condition_pars
!***********************************************************************
subroutine write_initial_condition_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=initial_condition_pars)
!
endsubroutine write_initial_condition_pars
!***********************************************************************
!***********************************************************************
!************ DO NOT DELETE THE FOLLOWING **************
!********************************************************************
!** This is an automatically generated include file that creates **
!** copies dummy routines from noinitial_condition.f90 for any **
!** InitialCondition routines not implemented in this file **
!** **
include '../initial_condition_dummies.inc'
!********************************************************************
endmodule InitialCondition