! Initial condition (density, magnetic field, velocity) ! for a particular configuration of a braided magnetic field ! with magnetic domes with nulls. ! !** 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 ! implicit none ! include '../initial_condition.h' ! ! ampl = amplitude of the magnetic field ! real :: ampl = 1.0 logical :: stretch = .false. ! namelist /initial_condition_pars/ & ampl, stretch ! contains !*********************************************************************** subroutine register_initial_condition() ! ! Configure pre-initialised (i.e. before parameter read) variables ! which should be know to be able to evaluate ! ! 16-december-14/simon: coded ! ! Identify CVS/SVN version information. ! if (lroot) call svn_id( & "$Id$") ! 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. ! ! 16-december-14/simon: coded ! ! Braided magnetic flux tube starting from the lower xy-plane and ! ending at the top plane together with magnetic domes with nulls ! ! Created 2014-12-16 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 ! ! clear the magnetic field to zero f(:,:,:,iax:iaz) = 0. if (stretch .eqv. .false.) then do n = 1, mz, 1 do m = 1, my, 1 do l = 1, mx, 1 f(l,m,n,iax) = 0 f(l,m,n,iay) = ampl * (& -16. * ((x(l) - 0.2) ** 2 + (y(m) - 0.2) ** 2 + (z(n) -& 25.8) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) - 0.2) **& 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l)& - 0.2) ** 2 + (y(m) - 0.2) ** 2 + (z(n) - 27.6) ** 2) **& (-1. / 2.) * (x(l) - 0.2) / ((y(m) - 0.2) ** 2 + (z(n) -& 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l) - 0.2) ** 2 +& (y(m) + 0.2) ** 2 + (z(n) + 25.8) ** 2) ** (-1. / 2.) * (x(l)& - 0.2) / ((y(m) + 0.2) ** 2 + (z(n) + 25.8) ** 2) * (-z(n)& - 25.8) - 18.95 * ((x(l) - 0.2) ** 2 + (y(m) + 0.2) **& 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m)& + 0.2) ** 2 + (z(n) + 27.6) ** 2) * (-z(n) - 27.6) +& 18.95 * ((x(l) + 8.2) ** 2 + (y(m) - 7.8) ** 2 + (z(n) -& 27.6) ** 2) ** (-1. / 2.) * (x(l) + 8.2) / ((y(m) - 7.8) **& 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l)& + 8.2) ** 2 + (y(m) - 8.2) ** 2 + (z(n) + 25.8) ** 2) ** (-1.& / 2.) * (x(l) + 8.2) / ((y(m) - 8.2) ** 2 + (z(n) + 25.8)& ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) + 8.2) ** 2 +& (y(m) - 8.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l)& + 8.2) / ((y(m) - 8.2) ** 2 + (z(n) + 27.6) ** 2) * (-z(n)& - 27.6) - 16. * ((x(l) - 7.8) ** 2 + (y(m) + 8.2) **& 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) - 7.8) / ((y(m)& + 8.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95& * ((x(l) - 7.8) ** 2 + (y(m) + 8.2) ** 2 + (z(n) - 27.6)& ** 2) ** (-1. / 2.) * (x(l) - 7.8) / ((y(m) + 8.2) **& 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l) -& 7.8) ** 2 + (y(m) + 7.8) ** 2 + (z(n) + 25.8) ** 2) ** (-1.& / 2.) * (x(l) - 7.8) / ((y(m) + 7.8) ** 2 + (z(n) + 25.8)& ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 7.8) ** 2 + (y(m)& + 7.8) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l)& - 7.8) / ((y(m) + 7.8) ** 2 + (z(n) + 27.6) ** 2) * (-z(n)& - 27.6) - 16. * ((x(l) - 7.8) ** 2 + (y(m) - 7.8) ** 2& + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) - 7.8) / ((y(m)& - 7.8) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95& * ((x(l) - 7.8) ** 2 + (y(m) - 7.8) ** 2 + (z(n) - 27.6)& ** 2) ** (-1. / 2.) * (x(l) - 7.8) / ((y(m) - 7.8) ** 2& + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + x(l) + 16. * ((x(l)& - 7.8) ** 2 + (y(m) - 8.2) ** 2 + (z(n) + 25.8) ** 2) **& (-1. / 2.) * (x(l) - 7.8) / ((y(m) - 8.2) ** 2 + (z(n) +& 25.8) ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 7.8) ** 2& + (y(m) - 8.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) *& (x(l) - 7.8) / ((y(m) - 8.2) ** 2 + (z(n) + 27.6) ** 2) *& (-z(n) - 27.6) - 16. * ((x(l) - 0.2) ** 2 + (y(m) + 8.2)& ** 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) - 0.2) /& ((y(m) + 8.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8)& + 18.95 * ((x(l) - 0.2) ** 2 + (y(m) + 8.2) ** 2 + (z(n) -& 27.6) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) + 8.2)& ** 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l)& - 0.2) ** 2 + (y(m) + 7.8) ** 2 + (z(n) + 25.8) ** 2) **& (-1. / 2.) * (x(l) - 0.2) / ((y(m) + 7.8) ** 2 + (z(n) + 25.8)& ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 0.2) ** 2 +& (y(m) + 7.8) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l)& - 0.2) / ((y(m) + 7.8) ** 2 + (z(n) + 27.6) ** 2) * (-z(n)& - 27.6) - 16. * ((x(l) + 8.2) ** 2 + (y(m) + 8.2) **& 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) + 8.2) / ((y(m)& + 8.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) +& 18.95 * ((x(l) + 8.2) ** 2 + (y(m) + 8.2) ** 2 + (z(n) - 27.6)& ** 2) ** (-1. / 2.) * (x(l) + 8.2) / ((y(m) + 8.2) **& 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l) +& 8.2) ** 2 + (y(m) + 7.8) ** 2 + (z(n) + 25.8) ** 2) ** (-1.& / 2.) * (x(l) + 8.2) / ((y(m) + 7.8) ** 2 + (z(n) + 25.8)& ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) + 8.2) ** 2 + (y(m)& + 7.8) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l)& + 8.2) / ((y(m) + 7.8) ** 2 + (z(n) + 27.6) ** 2) * (-z(n)& - 27.6) - 16. * ((x(l) + 8.2) ** 2 + (y(m) - 7.8) ** 2& + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) + 8.2) / ((y(m)& - 7.8) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) - 16.& * ((x(l) + 8.2) ** 2 + (y(m) - 0.2) ** 2 + (z(n) - 25.8)& ** 2) ** (-1. / 2.) * (x(l) + 8.2) / ((y(m) - 0.2) ** 2 +& (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l) + 8.2)& ** 2 + (y(m) - 0.2) ** 2 + (z(n) - 27.6) ** 2) ** (-1.& / 2.) * (x(l) + 8.2) / ((y(m) - 0.2) ** 2 + (z(n) - 27.6)& ** 2) * (-z(n) + 27.6) + 16. * ((x(l) + 8.2) ** 2 + (y(m)& + 0.2) ** 2 + (z(n) + 25.8) ** 2) ** (-1. / 2.) * (x(l) +& 8.2) / ((y(m) + 0.2) ** 2 + (z(n) + 25.8) ** 2) * (-z(n)& - 25.8) - 18.95 * ((x(l) + 8.2) ** 2 + (y(m) + 0.2) ** 2 +& (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l) + 8.2) / ((y(m)& + 0.2) ** 2 + (z(n) + 27.6) ** 2) * (-z(n) - 27.6) - 16.& * ((x(l) - 7.8) ** 2 + (y(m) - 0.2) ** 2 + (z(n) - 25.8) **& 2) ** (-1. / 2.) * (x(l) - 7.8) / ((y(m) - 0.2) ** 2 + (z(n)& - 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l) - 7.8)& ** 2 + (y(m) - 0.2) ** 2 + (z(n) - 27.6) ** 2) ** (-1. /& 2.) * (x(l) - 7.8) / ((y(m) - 0.2) ** 2 + (z(n) - 27.6) **& 2) * (-z(n) + 27.6) + 16. * ((x(l) - 7.8) ** 2 + (y(m) +& 0.2) ** 2 + (z(n) + 25.8) ** 2) ** (-1. / 2.) * (x(l) - 7.8)& / ((y(m) + 0.2) ** 2 + (z(n) + 25.8) ** 2) * (-z(n) -& 25.8) - 18.95 * ((x(l) - 7.8) ** 2 + (y(m) + 0.2) ** 2 + (z(n)& + 27.6) ** 2) ** (-1. / 2.) * (x(l) - 7.8) / ((y(m) +& 0.2) ** 2 + (z(n) + 27.6) ** 2) * (-z(n) - 27.6) - 16. *& ((x(l) - 0.2) ** 2 + (y(m) - 7.8) ** 2 + (z(n) - 25.8) **& 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) - 7.8) ** 2 + (z(n)& - 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l) - 0.2)& ** 2 + (y(m) - 7.8) ** 2 + (z(n) - 27.6) ** 2) ** (-1. / 2.)& * (x(l) - 0.2) / ((y(m) - 7.8) ** 2 + (z(n) - 27.6) **& 2) * (-z(n) + 27.6) + 16. * ((x(l) - 0.2) ** 2 + (y(m) - 8.2)& ** 2 + (z(n) + 25.8) ** 2) ** (-1. / 2.) * (x(l) - 0.2)& / ((y(m) - 8.2) ** 2 + (z(n) + 25.8) ** 2) * (-z(n) - 25.8)& - 18.95 * ((x(l) - 0.2) ** 2 + (y(m) - 8.2) ** 2 + (z(n)& + 27.6) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) - 8.2)& ** 2 + (z(n) + 27.6) ** 2) * (-z(n) - 27.6)) f(l,m,n,iaz) = ampl * (& -exp(-((x(l) + 1) ** 2) / 2. - y(m) ** 2 / 2. - ((z(n)& - 20) ** 2) / 4.) * sqrt(2.) + exp(-((x(l) - 1)& ** 2) / 2. - y(m) ** 2 / 2. - ((z(n) - 12) ** 2) / 4.)& * sqrt(2.) - exp(-((x(l) + 1) ** 2) / 2. - y(m) ** 2& / 2. - ((z(n) - 4) ** 2) / 4.) * sqrt(2.) + exp(-((x(l)& - 1) ** 2) / 2. - y(m) ** 2 / 2. - ((z(n) + 4)& ** 2) / 4.) * sqrt(2.) - exp(-((x(l) + 1) ** 2) / 2. -& y(m) ** 2 / 2. - ((z(n) + 12) ** 2) / 4.) * sqrt(2.)& + exp(-((x(l) - 1) ** 2) / 2. - y(m) ** 2 / 2. - ((z(n)& + 20) ** 2) / 4.) * sqrt(2.) - 16. * (((x(l)) -& 0.2) ** 2 + (y(m) + 8.2) ** 2 + ((z(n)) - 25.8) ** 2)& ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 8.2) ** 2 +& ((z(n)) - 25.8) ** 2) * (y(m) + 8.2) + 18.95 * (((x(l))& - 0.2) ** 2 + (y(m) + 8.2) ** 2 + ((z(n)) - 27.6)& ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 8.2)& ** 2 + ((z(n)) - 27.6) ** 2) * (y(m) + 8.2) + 16. * (((x(l))& - 0.2) ** 2 + (y(m) + 7.8) ** 2 + ((z(n))& + 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m)& + 7.8) ** 2 + ((z(n)) + 25.8) ** 2) * (y(m) + 7.8) - 18.95& * (((x(l)) - 0.2) ** 2 + (y(m) + 7.8) ** 2 + ((z(n))& + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) /& ((y(m) + 7.8) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) +& 7.8) - 16. * (((x(l)) + 8.2) ** 2 + (y(m) + 8.2) ** 2& + ((z(n)) - 25.8) ** 2) ** (-1. / 2.) * ((x(l)) +& 8.2) / ((y(m) + 8.2) ** 2 + ((z(n)) - 25.8) ** 2) * (y(m)& + 8.2) + 18.95 * (((x(l)) + 8.2) ** 2 + (y(m) + 8.2)& ** 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l))& + 8.2) / ((y(m) + 8.2) ** 2 + ((z(n)) - 27.6) **& 2) * (y(m) + 8.2) + 16. * (((x(l)) + 8.2) ** 2 + (y(m)& + 7.8) ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.) *& ((x(l)) + 8.2) / ((y(m) + 7.8) ** 2 + ((z(n)) + 25.8)& ** 2) * (y(m) + 7.8) - 18.95 * (((x(l)) + 8.2) **& 2 + (y(m) + 7.8) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1.& / 2.) * ((x(l)) + 8.2) / ((y(m) + 7.8) ** 2 + ((z(n))& + 27.6) ** 2) * (y(m) + 7.8) - 16. * (((x(l)) + 8.2)& ** 2 + (y(m) - 7.8) ** 2 + ((z(n)) - 25.8) ** 2) **& (-1. / 2.) * ((x(l)) + 8.2) / ((y(m) - 7.8) ** 2 + ((z(n))& - 25.8) ** 2) * (y(m) - 7.8) + 18.95 * (((x(l))& + 8.2) ** 2 + (y(m) - 7.8) ** 2 + ((z(n)) - 27.6)& ** 2) ** (-1. / 2.) * ((x(l)) + 8.2) / ((y(m) - 7.8) **& 2 + ((z(n)) - 27.6) ** 2) * (y(m) - 7.8) + 16. * (((x(l))& + 8.2) ** 2 + (y(m) - 8.2) ** 2 + ((z(n)) +& 25.8) ** 2) ** (-1. / 2.) * ((x(l)) + 8.2) / ((y(m) -& 8.2) ** 2 + ((z(n)) + 25.8) ** 2) * (y(m) - 8.2) - 18.95& * (((x(l)) + 8.2) ** 2 + (y(m) - 8.2) ** 2 + ((z(n))& + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) + 8.2) / ((y(m)& - 8.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) - 8.2)& - 16. * (((x(l)) + 8.2) ** 2 + (y(m) - 0.2) ** 2 +& ((z(n)) - 25.8) ** 2) ** (-1. / 2.) * ((x(l)) + 8.2)& / ((y(m) - 0.2) ** 2 + ((z(n)) - 25.8) ** 2) * (y(m)& - 0.2) + 18.95 * (((x(l)) + 8.2) ** 2 + (y(m) - 0.2)& ** 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l))& + 8.2) / ((y(m) - 0.2) ** 2 + ((z(n)) - 27.6) ** 2)& * (y(m) - 0.2) + 16. * (((x(l)) + 8.2) ** 2 + (y(m)& + 0.2) ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l))& + 8.2) / ((y(m) + 0.2) ** 2 + ((z(n)) + 25.8)& ** 2) * (y(m) + 0.2) - 18.95 * (((x(l)) + 8.2) ** 2& + (y(m) + 0.2) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1. /& 2.) * ((x(l)) + 8.2) / ((y(m) + 0.2) ** 2 + ((z(n))& + 27.6) ** 2) * (y(m) + 0.2) - 16. * (((x(l)) - 7.8)& ** 2 + (y(m) - 0.2) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1.& / 2.) * ((x(l)) - 7.8) / ((y(m) - 0.2) ** 2 + ((z(n))& - 25.8) ** 2) * (y(m) - 0.2) + 18.95 * (((x(l))& - 7.8) ** 2 + (y(m) - 0.2) ** 2 + ((z(n)) - 27.6) **& 2) ** (-1. / 2.) * ((x(l)) - 7.8) / ((y(m) - 0.2) **& 2 + ((z(n)) - 27.6) ** 2) * (y(m) - 0.2) + 16. * (((x(l))& - 7.8) ** 2 + (y(m) + 0.2) ** 2 + ((z(n)) + 25.8)& ** 2) ** (-1. / 2.) * ((x(l)) - 7.8) / ((y(m) + 0.2)& ** 2 + ((z(n)) + 25.8) ** 2) * (y(m) + 0.2) - 18.95& * (((x(l)) - 7.8) ** 2 + (y(m) + 0.2) ** 2 + ((z(n))& + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 7.8) / ((y(m)& + 0.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) + 0.2)& - 16. * (((x(l)) - 0.2) ** 2 + (y(m) - 7.8) ** 2 + ((z(n))& - 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2)& / ((y(m) - 7.8) ** 2 + ((z(n)) - 25.8) ** 2) * (y(m)& - 7.8) + 18.95 * (((x(l)) - 0.2) ** 2 + (y(m) - 7.8) **& 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l))& - 0.2) / ((y(m) - 7.8) ** 2 + ((z(n)) - 27.6) ** 2)& * (y(m) - 7.8) + 16. * (((x(l)) - 0.2) ** 2 + (y(m) -& 8.2) ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l))& - 0.2) / ((y(m) - 8.2) ** 2 + ((z(n)) + 25.8)& ** 2) * (y(m) - 8.2) - 18.95 * (((x(l)) - 0.2) ** 2 +& (y(m) - 8.2) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1. / 2.)& * ((x(l)) - 0.2) / ((y(m) - 8.2) ** 2 + ((z(n))& + 27.6) ** 2) * (y(m) - 8.2) - 16. * (((x(l)) - 0.2) **& 2 + (y(m) - 0.2) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1.& / 2.) * ((x(l)) - 0.2) / ((y(m) - 0.2) ** 2 + ((z(n))& - 25.8) ** 2) * (y(m) - 0.2) + 18.95 * (((x(l))& - 0.2) ** 2 + (y(m) - 0.2) ** 2 + ((z(n)) - 27.6) ** 2)& ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) - 0.2) ** 2& + ((z(n)) - 27.6) ** 2) * (y(m) - 0.2) + 16. * (((x(l))& - 0.2) ** 2 + (y(m) + 0.2) ** 2 + ((z(n)) + 25.8)& ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 0.2)& ** 2 + ((z(n)) + 25.8) ** 2) * (y(m) + 0.2) - 18.95 *& (((x(l)) - 0.2) ** 2 + (y(m) + 0.2) ** 2 + ((z(n))& + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m)& + 0.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) + 0.2) -& 16. * (((x(l)) - 7.8) ** 2 + (y(m) + 8.2) ** 2 + ((z(n))& - 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 7.8) /& ((y(m) + 8.2) ** 2 + ((z(n)) - 25.8) ** 2) * (y(m) +& 8.2) + 18.95 * (((x(l)) - 7.8) ** 2 + (y(m) + 8.2) **& 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l))& - 7.8) / ((y(m) + 8.2) ** 2 + ((z(n)) - 27.6) ** 2) *& (y(m) + 8.2) + 16. * (((x(l)) - 7.8) ** 2 + (y(m) + 7.8)& ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l))& - 7.8) / ((y(m) + 7.8) ** 2 + ((z(n)) + 25.8) **& 2) * (y(m) + 7.8) - 18.95 * (((x(l)) - 7.8) ** 2 + (y(m)& + 7.8) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1. / 2.)& * ((x(l)) - 7.8) / ((y(m) + 7.8) ** 2 + ((z(n)) +& 27.6) ** 2) * (y(m) + 7.8) - 16. * (((x(l)) - 7.8) **& 2 + (y(m) - 7.8) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1.& / 2.) * ((x(l)) - 7.8) / ((y(m) - 7.8) ** 2 + ((z(n))& - 25.8) ** 2) * (y(m) - 7.8) + 18.95 * (((x(l)) -& 7.8) ** 2 + (y(m) - 7.8) ** 2 + ((z(n)) - 27.6) ** 2)& ** (-1. / 2.) * ((x(l)) - 7.8) / ((y(m) - 7.8) ** 2 +& ((z(n)) - 27.6) ** 2) * (y(m) - 7.8) + 16. * (((x(l))& - 7.8) ** 2 + (y(m) - 8.2) ** 2 + ((z(n)) + 25.8)& ** 2) ** (-1. / 2.) * ((x(l)) - 7.8) / ((y(m) - 8.2) **& 2 + ((z(n)) + 25.8) ** 2) * (y(m) - 8.2) - 18.95 * (((x(l))& - 7.8) ** 2 + (y(m) - 8.2) ** 2 + ((z(n))& + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 7.8) / ((y(m)& - 8.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) - 8.2)) enddo enddo enddo else do n = 1, mz, 1 do m = 1, my, 1 do l = 1, mx, 1 f(l,m,n,iax) = 0 f(l,m,n,iay) = ampl * (& -16. * ((x(l) - 15.8) ** 2 + (y(m) - 15.8) ** 2 + (z(n)& - 25.8) ** 2) ** (-1. / 2.) * (x(l) - 15.8) / ((y(m) - 15.8)& ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l)& - 15.8) ** 2 + (y(m) - 15.8) ** 2 + (z(n) - 27.6) **& 2) ** (-1. / 2.) * (x(l) - 15.8) / ((y(m) - 15.8) ** 2 +& (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l) - 15.8)& ** 2 + (y(m) - 16.2) ** 2 + (z(n) + 25.8) ** 2) ** (-1.& / 2.) * (x(l) - 15.8) / ((y(m) - 16.2) ** 2 + (z(n) + 25.8)& ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 15.8) ** 2 + (y(m)& - 16.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l)& - 15.8) / ((y(m) - 16.2) ** 2 + (z(n) + 27.6) ** 2) *& (-z(n) - 27.6) + 18.95 * ((x(l) + 16.2) ** 2 + (y(m) + 16.2)& ** 2 + (z(n) - 27.6) ** 2) ** (-1. / 2.) * (x(l) + 16.2)& / ((y(m) + 16.2) ** 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6)& + 16. * ((x(l) + 16.2) ** 2 + (y(m) + 15.8) ** 2 + (z(n)& + 25.8) ** 2) ** (-1. / 2.) * (x(l) + 16.2) / ((y(m) +& 15.8) ** 2 + (z(n) + 25.8) ** 2) * (-z(n) - 25.8) - 18.95& * ((x(l) + 16.2) ** 2 + (y(m) + 15.8) ** 2 + (z(n) + 27.6)& ** 2) ** (-1. / 2.) * (x(l) + 16.2) / ((y(m) + 15.8) **& 2 + (z(n) + 27.6) ** 2) * (-z(n) - 27.6) - 16. * ((x(l) +& 16.2) ** 2 + (y(m) - 15.8) ** 2 + (z(n) - 25.8) ** 2) ** (-1.& / 2.) * (x(l) + 16.2) / ((y(m) - 15.8) ** 2 + (z(n) -& 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l) + 16.2) ** 2& + (y(m) - 15.8) ** 2 + (z(n) - 27.6) ** 2) ** (-1. / 2.)& * (x(l) + 16.2) / ((y(m) - 15.8) ** 2 + (z(n) - 27.6) ** 2)& * (-z(n) + 27.6) + 16. * ((x(l) + 16.2) ** 2 + (y(m) - 16.2)& ** 2 + (z(n) + 25.8) ** 2) ** (-1. / 2.) * (x(l) + 16.2)& / ((y(m) - 16.2) ** 2 + (z(n) + 25.8) ** 2) * (-z(n) -& 25.8) - 18.95 * ((x(l) + 16.2) ** 2 + (y(m) - 16.2) ** 2& + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l) + 16.2) / ((y(m)& - 16.2) ** 2 + (z(n) + 27.6) ** 2) * (-z(n) - 27.6) - 16.& * ((x(l) - 15.8) ** 2 + (y(m) + 16.2) ** 2 + (z(n) - 25.8)& ** 2) ** (-1. / 2.) * (x(l) - 15.8) / ((y(m) + 16.2) **& 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95 * ((x(l)& - 15.8) ** 2 + (y(m) + 16.2) ** 2 + (z(n) - 27.6) ** 2)& ** (-1. / 2.) * (x(l) - 15.8) / ((y(m) + 16.2) ** 2 + (z(n)& - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l) - 15.8) **& 2 + (y(m) + 15.8) ** 2 + (z(n) + 25.8) ** 2) ** (-1. / 2.)& * (x(l) - 15.8) / ((y(m) + 15.8) ** 2 + (z(n) + 25.8) **& 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 15.8) ** 2 + (y(m)& + 15.8) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) * (x(l)& - 15.8) / ((y(m) + 15.8) ** 2 + (z(n) + 27.6) ** 2) * (-z(n)& - 27.6) - 16. * ((x(l) - 15.8) ** 2 + (y(m) - 0.2) ** 2& + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) - 15.8) / ((y(m)& - 0.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8) + 18.95& * ((x(l) - 15.8) ** 2 + (y(m) - 0.2) ** 2 + (z(n) - 27.6)& ** 2) ** (-1. / 2.) * (x(l) - 15.8) / ((y(m) - 0.2) **& 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l)& - 15.8) ** 2 + (y(m) + 0.2) ** 2 + (z(n) + 25.8) ** 2) **& (-1. / 2.) * (x(l) - 15.8) / ((y(m) + 0.2) ** 2 + (z(n) +& 25.8) ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 15.8) ** 2& + (y(m) + 0.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. / 2.) *& (x(l) - 15.8) / ((y(m) + 0.2) ** 2 + (z(n) + 27.6) ** 2)& * (-z(n) - 27.6) - 16. * ((x(l) - 0.2) ** 2 + (y(m) - 15.8)& ** 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) - 0.2)& / ((y(m) - 15.8) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8)& + 18.95 * ((x(l) - 0.2) ** 2 + (y(m) - 15.8) ** 2 + (z(n)& - 27.6) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) - 15.8)& ** 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l)& - 0.2) ** 2 + (y(m) - 16.2) ** 2 + (z(n) + 25.8) **& 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) - 16.2) ** 2 + (z(n)& + 25.8) ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 0.2)& ** 2 + (y(m) - 16.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. /& 2.) * (x(l) - 0.2) / ((y(m) - 16.2) ** 2 + (z(n) + 27.6)& ** 2) * (-z(n) - 27.6) - 16. * ((x(l) - 0.2) ** 2 + (y(m)& + 16.2) ** 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) -& 0.2) / ((y(m) + 16.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n)& + 25.8) + 18.95 * ((x(l) - 0.2) ** 2 + (y(m) + 16.2) ** 2& + (z(n) - 27.6) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m)& + 16.2) ** 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16.& * ((x(l) - 0.2) ** 2 + (y(m) + 15.8) ** 2 + (z(n) + 25.8)& ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) + 15.8) **& 2 + (z(n) + 25.8) ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l)& - 0.2) ** 2 + (y(m) + 15.8) ** 2 + (z(n) + 27.6) ** 2) **& (-1. / 2.) * (x(l) - 0.2) / ((y(m) + 15.8) ** 2 + (z(n) +& 27.6) ** 2) * (-z(n) - 27.6) - 16. * ((x(l) + 16.2) ** 2 +& (y(m) + 16.2) ** 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) *& (x(l) + 16.2) / ((y(m) + 16.2) ** 2 + (z(n) - 25.8) ** 2)& * (-z(n) + 25.8) - 16. * ((x(l) + 16.2) ** 2 + (y(m) - 0.2)& ** 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) + 16.2)& / ((y(m) - 0.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) + 25.8)& + 18.95 * ((x(l) + 16.2) ** 2 + (y(m) - 0.2) ** 2 + (z(n)& - 27.6) ** 2) ** (-1. / 2.) * (x(l) + 16.2) / ((y(m) -& 0.2) ** 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. * ((x(l)& + 16.2) ** 2 + (y(m) + 0.2) ** 2 + (z(n) + 25.8) **& 2) ** (-1. / 2.) * (x(l) + 16.2) / ((y(m) + 0.2) ** 2 + (z(n)& + 25.8) ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) + 16.2)& ** 2 + (y(m) + 0.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. /& 2.) * (x(l) + 16.2) / ((y(m) + 0.2) ** 2 + (z(n) + 27.6)& ** 2) * (-z(n) - 27.6) - 16. * ((x(l) - 0.2) ** 2 + (y(m)& - 0.2) ** 2 + (z(n) - 25.8) ** 2) ** (-1. / 2.) * (x(l) -& 0.2) / ((y(m) - 0.2) ** 2 + (z(n) - 25.8) ** 2) * (-z(n) +& 25.8) + 18.95 * ((x(l) - 0.2) ** 2 + (y(m) - 0.2) ** 2 +& (z(n) - 27.6) ** 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m)& - 0.2) ** 2 + (z(n) - 27.6) ** 2) * (-z(n) + 27.6) + 16. *& ((x(l) - 0.2) ** 2 + (y(m) + 0.2) ** 2 + (z(n) + 25.8) **& 2) ** (-1. / 2.) * (x(l) - 0.2) / ((y(m) + 0.2) ** 2 + (z(n)& + 25.8) ** 2) * (-z(n) - 25.8) - 18.95 * ((x(l) - 0.2)& ** 2 + (y(m) + 0.2) ** 2 + (z(n) + 27.6) ** 2) ** (-1. /& 2.) * (x(l) - 0.2) / ((y(m) + 0.2) ** 2 + (z(n) + 27.6) **& 2) * (-z(n) - 27.6) + x(l)) f(l,m,n,iaz) = ampl * (& -exp(-((x(l) + 1) ** 2) / 2. - y(m) ** 2 / 2. - ((z(n) -& 20) ** 2) / 4.) * sqrt(2.) + exp(-((x(l) - 1) ** 2) / 2.& - y(m) ** 2 / 2. - ((z(n) - 12) ** 2) / 4.) * sqrt(2.) - exp(-((x(l)& + 1) ** 2) / 2. - y(m) ** 2 / 2. - ((z(n) - 4)& ** 2) / 4.) * sqrt(2.) + exp(-((x(l) - 1) ** 2) / 2. - y(m)& ** 2 / 2. - ((z(n) + 4) ** 2) / 4.) * sqrt(2.) - exp(-((x(l)& + 1) ** 2) / 2. - y(m) ** 2 / 2. - ((z(n) + 12) ** 2)& / 4.) * sqrt(2.) + exp(-((x(l) - 1) ** 2) / 2. - y(m) **& 2 / 2. - ((z(n) + 20) ** 2) / 4.) * sqrt(2.) - 16. * (((x(l))& - 15.8) ** 2 + (y(m) - 15.8) ** 2 + ((z(n)) - 25.8) **& 2) ** (-1. / 2.) * ((x(l)) - 15.8) / ((y(m) - 15.8) ** 2& + ((z(n)) - 25.8) ** 2) * (y(m) - 15.8) + 18.95 * (((x(l))& - 15.8) ** 2 + (y(m) - 15.8) ** 2 + ((z(n)) - 27.6) ** 2)& ** (-1. / 2.) * ((x(l)) - 15.8) / ((y(m) - 15.8) ** 2 + ((z(n))& - 27.6) ** 2) * (y(m) - 15.8) + 16. * (((x(l)) - 15.8)& ** 2 + (y(m) - 16.2) ** 2 + ((z(n)) + 25.8) ** 2) ** (-1.& / 2.) * ((x(l)) - 15.8) / ((y(m) - 16.2) ** 2 + ((z(n))& + 25.8) ** 2) * (y(m) - 16.2) - 18.95 * (((x(l)) - 15.8)& ** 2 + (y(m) - 16.2) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1.& / 2.) * ((x(l)) - 15.8) / ((y(m) - 16.2) ** 2 + ((z(n)) +& 27.6) ** 2) * (y(m) - 16.2) - 16. * (((x(l)) + 16.2) ** 2& + (y(m) - 15.8) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1. / 2.)& * ((x(l)) + 16.2) / ((y(m) - 15.8) ** 2 + ((z(n)) - 25.8)& ** 2) * (y(m) - 15.8) + 18.95 * (((x(l)) + 16.2) ** 2 +& (y(m) - 15.8) ** 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) *& ((x(l)) + 16.2) / ((y(m) - 15.8) ** 2 + ((z(n)) - 27.6) **& 2) * (y(m) - 15.8) + 16. * (((x(l)) + 16.2) ** 2 + (y(m)& - 16.2) ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l))& + 16.2) / ((y(m) - 16.2) ** 2 + ((z(n)) + 25.8) ** 2)& * (y(m) - 16.2) - 18.95 * (((x(l)) + 16.2) ** 2 + (y(m) -& 16.2) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1. / 2.) * ((x(l))& + 16.2) / ((y(m) - 16.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m)& - 16.2) - 16. * (((x(l)) - 15.8) ** 2 + (y(m) + 16.2)& ** 2 + ((z(n)) - 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 15.8)& / ((y(m) + 16.2) ** 2 + ((z(n)) - 25.8) ** 2) * (y(m)& + 16.2) + 18.95 * (((x(l)) - 15.8) ** 2 + (y(m) + 16.2) **& 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 15.8)& / ((y(m) + 16.2) ** 2 + ((z(n)) - 27.6) ** 2) * (y(m) + 16.2)& + 16. * (((x(l)) - 15.8) ** 2 + (y(m) + 15.8) ** 2 +& ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 15.8) / ((y(m)& + 15.8) ** 2 + ((z(n)) + 25.8) ** 2) * (y(m) + 15.8)& - 18.95 * (((x(l)) - 15.8) ** 2 + (y(m) + 15.8) ** 2 + ((z(n))& + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 15.8) / ((y(m)& + 15.8) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) + 15.8) - 16.& * (((x(l)) - 15.8) ** 2 + (y(m) - 0.2) ** 2 + ((z(n)) -& 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 15.8) / ((y(m) - 0.2)& ** 2 + ((z(n)) - 25.8) ** 2) * (y(m) - 0.2) + 18.95 * (((x(l))& - 15.8) ** 2 + (y(m) - 0.2) ** 2 + ((z(n)) - 27.6)& ** 2) ** (-1. / 2.) * ((x(l)) - 15.8) / ((y(m) - 0.2) **& 2 + ((z(n)) - 27.6) ** 2) * (y(m) - 0.2) + 16. * (((x(l))& - 15.8) ** 2 + (y(m) + 0.2) ** 2 + ((z(n)) + 25.8) ** 2) **& (-1. / 2.) * ((x(l)) - 15.8) / ((y(m) + 0.2) ** 2 + ((z(n))& + 25.8) ** 2) * (y(m) + 0.2) - 18.95 * (((x(l)) - 15.8)& ** 2 + (y(m) + 0.2) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1.& / 2.) * ((x(l)) - 15.8) / ((y(m) + 0.2) ** 2 + ((z(n)) +& 27.6) ** 2) * (y(m) + 0.2) - 16. * (((x(l)) - 0.2) ** 2 +& (y(m) - 15.8) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1. / 2.) *& ((x(l)) - 0.2) / ((y(m) - 15.8) ** 2 + ((z(n)) - 25.8) **& 2) * (y(m) - 15.8) + 18.95 * (((x(l)) - 0.2) ** 2 + (y(m)& - 15.8) ** 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l))& - 0.2) / ((y(m) - 15.8) ** 2 + ((z(n)) - 27.6) ** 2) *& (y(m) - 15.8) + 16. * (((x(l)) - 0.2) ** 2 + (y(m) - 16.2)& ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2)& / ((y(m) - 16.2) ** 2 + ((z(n)) + 25.8) ** 2) * (y(m)& - 16.2) - 18.95 * (((x(l)) - 0.2) ** 2 + (y(m) - 16.2) **& 2 + ((z(n)) + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) /& ((y(m) - 16.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) - 16.2)& - 16. * (((x(l)) - 0.2) ** 2 + (y(m) + 16.2) ** 2 + ((z(n))& - 25.8) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m)& + 16.2) ** 2 + ((z(n)) - 25.8) ** 2) * (y(m) + 16.2) + 18.95& * (((x(l)) - 0.2) ** 2 + (y(m) + 16.2) ** 2 + ((z(n))& - 27.6) ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 16.2)& ** 2 + ((z(n)) - 27.6) ** 2) * (y(m) + 16.2) + 16. * (((x(l))& - 0.2) ** 2 + (y(m) + 15.8) ** 2 + ((z(n)) + 25.8)& ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 15.8) **& 2 + ((z(n)) + 25.8) ** 2) * (y(m) + 15.8) - 18.95 * (((x(l))& - 0.2) ** 2 + (y(m) + 15.8) ** 2 + ((z(n)) + 27.6) ** 2)& ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 15.8) ** 2 + ((z(n))& + 27.6) ** 2) * (y(m) + 15.8) - 16. * (((x(l)) + 16.2)& ** 2 + (y(m) + 16.2) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1.& / 2.) * ((x(l)) + 16.2) / ((y(m) + 16.2) ** 2 + ((z(n))& - 25.8) ** 2) * (y(m) + 16.2) + 18.95 * (((x(l)) + 16.2)& ** 2 + (y(m) + 16.2) ** 2 + ((z(n)) - 27.6) ** 2) ** (-1.& / 2.) * ((x(l)) + 16.2) / ((y(m) + 16.2) ** 2 + ((z(n)) -& 27.6) ** 2) * (y(m) + 16.2) + 16. * (((x(l)) + 16.2) ** 2& + (y(m) + 15.8) ** 2 + ((z(n)) + 25.8) ** 2) ** (-1. / 2.)& * ((x(l)) + 16.2) / ((y(m) + 15.8) ** 2 + ((z(n)) + 25.8)& ** 2) * (y(m) + 15.8) - 18.95 * (((x(l)) + 16.2) ** 2 +& (y(m) + 15.8) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1. / 2.) *& ((x(l)) + 16.2) / ((y(m) + 15.8) ** 2 + ((z(n)) + 27.6) **& 2) * (y(m) + 15.8) - 16. * (((x(l)) + 16.2) ** 2 + (y(m)& - 0.2) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1. / 2.) * ((x(l))& + 16.2) / ((y(m) - 0.2) ** 2 + ((z(n)) - 25.8) ** 2) *& (y(m) - 0.2) + 18.95 * (((x(l)) + 16.2) ** 2 + (y(m) - 0.2)& ** 2 + ((z(n)) - 27.6) ** 2) ** (-1. / 2.) * ((x(l)) + 16.2)& / ((y(m) - 0.2) ** 2 + ((z(n)) - 27.6) ** 2) * (y(m)& - 0.2) + 16. * (((x(l)) + 16.2) ** 2 + (y(m) + 0.2) ** 2 +& ((z(n)) + 25.8) ** 2) ** (-1. / 2.) * ((x(l)) + 16.2) / ((y(m)& + 0.2) ** 2 + ((z(n)) + 25.8) ** 2) * (y(m) + 0.2) -& 18.95 * (((x(l)) + 16.2) ** 2 + (y(m) + 0.2) ** 2 + ((z(n))& + 27.6) ** 2) ** (-1. / 2.) * ((x(l)) + 16.2) / ((y(m)& + 0.2) ** 2 + ((z(n)) + 27.6) ** 2) * (y(m) + 0.2) + 18.95& * (((x(l)) - 0.2) ** 2 + (y(m) - 0.2) ** 2 + ((z(n)) - 27.6)& ** 2) ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) - 0.2) **& 2 + ((z(n)) - 27.6) ** 2) * (y(m) - 0.2) + 16. * (((x(l))& - 0.2) ** 2 + (y(m) + 0.2) ** 2 + ((z(n)) + 25.8) ** 2)& ** (-1. / 2.) * ((x(l)) - 0.2) / ((y(m) + 0.2) ** 2 + ((z(n))& + 25.8) ** 2) * (y(m) + 0.2) - 18.95 * (((x(l)) - 0.2)& ** 2 + (y(m) + 0.2) ** 2 + ((z(n)) + 27.6) ** 2) ** (-1.& / 2.) * ((x(l)) - 0.2) / ((y(m) + 0.2) ** 2 + ((z(n)) + 27.6)& ** 2) * (y(m) + 0.2) - 16. * (((x(l)) - 0.2) ** 2 + (y(m)& - 0.2) ** 2 + ((z(n)) - 25.8) ** 2) ** (-1. / 2.) * ((x(l))& - 0.2) / ((y(m) - 0.2) ** 2 + ((z(n)) - 25.8) ** 2)& * (y(m) - 0.2)) enddo enddo enddo end if ! 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