Subroutine that initializes the analytical or externally calculated electric and magnetic fields.
In this subroutine we load the parameters of the electric and magnetic fields from the namelists 'analytical_fields_params' and 'externalPlasmaModel' in the input file. Sign convention in analytical fields corresponds to DIII-D fields with and . Sign convention in analytical fields corresponds to DIII-D fields with and .
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(KORC_PARAMS), | intent(inout) | :: | params | Core KORC simulation parameters. |
||
| type(FIELDS), | intent(out) | :: | F | An instance of the KORC derived type FIELDS. Magnetic field at magnetic axis for an 'ANALYTICAL' magnetic field, or the magnitude of the magnetic field for a 'UNFIROM' plasma. Plasma edge as measured from the magnetic axis. Radial position of the magnetic axis Safety factor at the plasma edge. Safety factor at the magnetic axis . String with information about the direction of the plasma current, 'PARALLEL' or 'ANTI-PARALLEL' to the toroidal magnetic field. Electric field at the magnetic axis. Logical variable that specifies if the electric field is going to be used on in a given simulation. Logical variable that specifies if the magnetic field is going to be used on in a given simulation. Logical variable that specifies if the poloidal magnetic flux is going to be used on in a given simulation. Logical variable that specifies if the plasma is axisymmetric. |
subroutine initialize_fields(params,F)
!! @note Subroutine that initializes the analytical or externally
!! calculated electric and magnetic fields. @endnote
!! In this subroutine we load the parameters of the electric and
!! magnetic fields from the namelists 'analytical_fields_params' and
!! 'externalPlasmaModel' in the input file.
TYPE(KORC_PARAMS), INTENT(INOUT) :: params
!! Core KORC simulation parameters.
TYPE(FIELDS), INTENT(OUT) :: F
!! An instance of the KORC derived type FIELDS.
!REAL(rp) :: Bo
!! Magnetic field at magnetic axis for an 'ANALYTICAL' magnetic field,
!! or the magnitude of the magnetic field for a 'UNFIROM' plasma.
!REAL(rp) :: minor_radius
!! Plasma edge \(r_{edge}\) as measured from the magnetic axis.
!REAL(rp) :: major_radius
!! Radial position of the magnetic axis \(R_0\)
!REAL(rp) :: qa
!! Safety factor at the plasma edge.
!REAL(rp) :: qo
!! Safety factor at the magnetic axis \(q_0\).
!CHARACTER(MAX_STRING_LENGTH) :: current_direction
!! String with information about the direction of the plasma current,
!! 'PARALLEL' or 'ANTI-PARALLEL' to the toroidal magnetic field.
!CHARACTER(MAX_STRING_LENGTH) :: E_model
!REAL(rp) :: Eo,E_dyn,E_pulse,E_width
!! Electric field at the magnetic axis.
!LOGICAL :: Efield
!! Logical variable that specifies if the electric field is
!! going to be used on in a given simulation.
!LOGICAL :: dBfield
!LOGICAL :: Bfield
!! Logical variable that specifies if the magnetic field is
!! going to be used on in a given simulation.
!LOGICAL :: Bflux
!LOGICAL :: Bflux3D
!LOGICAL :: Dim2x1t
!LOGICAL :: E_2x1t,ReInterp_2x1t
!! Logical variable that specifies if the poloidal magnetic
!! flux is going to be used on in a given simulation.
!LOGICAL :: axisymmetric_fields
!! Logical variable that specifies if the plasma is axisymmetric.
INTEGER :: ii
!! Iterators for creating mesh for GC model with analytic fields
INTEGER :: kk
!! Iterators for creating mesh for GC model with analytic fields
!INTEGER :: nR
!! Number of mesh points in R for grid in GC model of analytical field
!INTEGER :: nZ,nPHI
!! Number of mesh points in Z for grid in GC model of analytical field
real(rp) :: rm
!! Minor radius at each position in the grid for
!! GC model of analytical field
real(rp) :: qr
!! Safety factor at each position in the grid for
!! GC model of analytical field
real(rp) :: theta
!! Poloidal angle at each position in the grid for
!! GC model of analytical field
logical :: test
!integer :: res_double
real(rp) :: RMAX,RMIN,ZMAX,ZMIN
!integer :: dim_1D,ind0_2x1t
!real(rp) :: dt_E_SC,Ip_exp,PSIp_lim,PSIp_0
!real(rp) :: t0_2x1t
!NAMELIST /analytical_fields_params/ Bo,minor_radius,major_radius,&
! qa,qo,Eo,current_direction,nR,nZ,nPHI,dim_1D,dt_E_SC,Ip_exp, &
! E_dyn,E_pulse,E_width
!NAMELIST /externalPlasmaModel/ Efield, Bfield, Bflux,Bflux3D,dBfield, &
! axisymmetric_fields, Eo,E_dyn,E_pulse,E_width,res_double, &
! dim_1D,dt_E_SC,Ip_exp,PSIp_lim,Dim2x1t,t0_2x1t,E_2x1t,ReInterp_2x1t, &
! ind0_2x1t,PSIp_0
#ifdef FIO
F%FIO_B = -1
F%FIO_E = -1
#endif
if (params%mpi_params%rank .EQ. 0) then
write(output_unit_write,'(/,"* * * * * * * * INITIALIZING FIELDS * * * * * * * *")')
end if
! SELECT CASE (TRIM(params%field_model))
if (params%field_model(1:10).eq.'ANALYTICAL') then
! CASE('ANALYTICAL')
! Load the parameters of the analytical magnetic field
!open(unit=default_unit_open,file=TRIM(params%path_to_inputs), &
! status='OLD',form='formatted')
!read(default_unit_open,nml=analytical_fields_params)
!close(default_unit_open)
F%AB%Bo = Bo
F%AB%a = minor_radius
F%AB%Ro = major_radius
F%Ro = major_radius
F%Zo = 0.0_rp
F%AB%qa = qa
F%AB%qo = qo
F%AB%lambda = F%AB%a/SQRT(qa/qo - 1.0_rp)
F%AB%Bpo = F%AB%lambda*F%AB%Bo/(F%AB%qo*F%AB%Ro)
F%AB%current_direction = TRIM(current_direction)
SELECT CASE (TRIM(F%AB%current_direction))
CASE('PARALLEL')
F%AB%Bp_sign = 1.0_rp
CASE('ANTI-PARALLEL')
F%AB%Bp_sign = -1.0_rp
CASE DEFAULT
END SELECT
F%Eo = Eo
F%Bo = F%AB%Bo
F%E_dyn = E_dyn
F%E_pulse = E_pulse
F%E_width = E_width
F%PSIp_lim=PSIp_lim
F%AB%Ero=Ero
F%AB%rmn=rmn
F%AB%sigmamn=sigmamn
!write(output_unit_write,*) E_dyn,E_pulse,E_width
F%res_double=res_double
if (params%mpi_params%rank .EQ. 0) then
write(output_unit_write,'("ANALYTIC")')
write(output_unit_write,'("Magnetic field: ",E17.10)') F%Bo
write(output_unit_write,'("Electric field: ",E17.10)') F%Eo
end if
if (params%field_eval.eq.'interp') then
F%dims(1) = nR
F%dims(2) = nPHI
F%dims(3) = nZ
if (params%field_model(12:14).eq.'PSI') then
F%axisymmetric_fields = .TRUE.
F%Bfield=.TRUE.
F%Efield=.TRUE.
F%Bflux=.TRUE.
call ALLOCATE_2D_FIELDS_ARRAYS(params,F,F%Bfield,F%Bflux, &
.false.,F%Efield,.FALSE.,.FALSE.)
do ii=1_idef,F%dims(1)
F%X%R(ii)=(F%Ro-F%AB%a)+(ii-1)*2*F%AB%a/(F%dims(1)-1)
end do
do ii=1_idef,F%dims(3)
F%X%Z(ii)=(F%Zo-F%AB%a)+(ii-1)*2*F%AB%a/(F%dims(3)-1)
end do
!write(6,*) F%X%R
!write(6,*) F%X%Z
do ii=1_idef,F%dims(1)
do kk=1_idef,F%dims(3)
rm=sqrt((F%X%R(ii)-F%Ro)**2+(F%X%Z(kk)-F%Zo)**2)
qr=F%AB%qo*(1+(rm/F%AB%lambda)**2)
theta=atan2(F%X%Z(kk)-F%Zo,F%X%R(ii)-F%Ro)
F%B_2D%R(ii,kk)=(rm/F%X%R(ii))* &
(F%AB%Bo/qr)*sin(theta)
F%B_2D%PHI(ii,kk)=-(F%Ro/F%X%R(ii))*F%AB%Bo
F%B_2D%Z(ii,kk)=-(rm/F%X%R(ii))* &
(F%AB%Bo/qr)*cos(theta)
F%E_2D%R(ii,kk)=0.0_rp
F%E_2D%PHI(ii,kk)=-(F%Ro/F%X%R(ii))*F%Eo
F%E_2D%Z(ii,kk)=0.0_rp
F%PSIp(ii,kk)=F%X%R(ii)*F%AB%lambda**2*F%Bo/ &
(2*F%AB%qo*(F%Ro+rm*cos(theta)))* &
log(1+(rm/F%AB%lambda)**2)
!! Sign convention in analytical fields corresponds to
!! DIII-D fields with \(B_\phi<0\) and \(B_\theta<0\).
F%FLAG2D=1.
end do
end do
F%FLAG2D(1:2,:)=0.
F%FLAG2D(F%dims(1)-1:F%dims(1),:)=0.
F%FLAG2D(:,1:2)=0.
F%FLAG2D(:,F%dims(3)-1:F%dims(3))=0.
if (F%Bflux) F%PSIP_min=minval(F%PSIp)
F%Bfield=.FALSE.
else if (params%field_model(12:13).eq.'3D') then
F%axisymmetric_fields = .FALSE.
F%Bfield=.TRUE.
F%Efield=.TRUE.
call ALLOCATE_3D_FIELDS_ARRAYS(params,F,F%Bfield,F%Efield,.false.)
do ii=1_idef,F%dims(1)
F%X%R(ii)=(F%Ro-F%AB%a)+(ii-1)*2*F%AB%a/(F%dims(1)-1)
end do
do ii=1_idef,F%dims(2)
F%X%PHI(ii)=0._rp+(ii-1)*2*C_PI/(F%dims(1)-1)
end do
do ii=1_idef,F%dims(3)
F%X%Z(ii)=(F%Zo-F%AB%a)+(ii-1)*2*F%AB%a/(F%dims(3)-1)
end do
!write(output_unit_write,*) size(F%B_3D%R)
do ii=1_idef,F%dims(1)
do kk=1_idef,F%dims(3)
!write(output_unit_write,*) ii,kk
rm=sqrt((F%X%R(ii)-F%Ro)**2+(F%X%Z(kk)-F%Zo)**2)
qr=F%AB%qo*(1+(rm/F%AB%lambda)**2)
theta=atan2(F%X%Z(kk)-F%Zo,F%X%R(ii)-F%Ro)
F%B_3D%R(ii,:,kk)=(rm/F%X%R(ii))* &
(F%AB%Bo/qr)*sin(theta)
F%B_3D%PHI(ii,:,kk)=-(F%Ro/F%X%R(ii))*F%AB%Bo
!write(output_unit_write,*) F%B_3D%PHI(ii,1,kk)
F%B_3D%Z(ii,:,kk)=-(rm/F%X%R(ii))* &
(F%AB%Bo/qr)*cos(theta)
F%E_3D%R(ii,:,kk)=0.0_rp
F%E_3D%PHI(ii,:,kk)=-(F%Ro/F%X%R(ii))*F%Eo
F%E_3D%Z(ii,:,kk)=0.0_rp
!! Sign convention in analytical fields corresponds to
!! DIII-D fields with \(B_\phi<0\) and \(B_\theta<0\).
F%FLAG3D=1.
end do
end do
F%FLAG3D(1:2,:,:)=0.
F%FLAG3D(F%dims(1)-1:F%dims(1),:,:)=0.
F%FLAG3D(:,:,1:2)=0.
F%FLAG3D(:,:,F%dims(3)-1:F%dims(3))=0.
end if
if (params%orbit_model(3:5).eq.'pre') then
if (params%mpi_params%rank .EQ. 0) then
write(output_unit_write,'("Initializing GC fields from analytic EM fields")')
end if
if (params%field_model(12:13).eq.'2D') then
call initialize_GC_fields(F)
else if (params%field_model(12:13).eq.'3D') then
call initialize_GC_fields_3D(F)
end if
end if
!F%Bfield= .FALSE.
!F%axisymmetric_fields = .TRUE.
!F%Bflux=.TRUE.
!F%Efield=.FALSE.
end if
if (params%SC_E) then
F%dim_1D=dim_1D
F%dt_E_SC=dt_E_SC
F%Ip_exp=Ip_exp
ALLOCATE(F%E_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%A1_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%A2_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%A3_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%J1_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%J2_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%J3_SC_1D%PHI(F%dim_1D))
ALLOCATE(F%r_1D(F%dim_1D))
F%E_SC_1D%PHI=0._rp
F%A1_SC_1D%PHI=0._rp
F%A2_SC_1D%PHI=0._rp
F%A3_SC_1D%PHI=0._rp
F%J1_SC_1D%PHI=0._rp
F%J2_SC_1D%PHI=0._rp
F%J3_SC_1D%PHI=0._rp
F%r_1D=0._rp
do ii=1_idef,F%dim_1D
F%r_1D(ii)=(ii-1)*F%AB%a/(F%dim_1D-1)
end do
end if
! CASE('EXTERNAL')
else if (params%field_model(1:8).eq.'EXTERNAL') then
! Load the magnetic field from an external HDF5 file
!open(unit=default_unit_open,file=TRIM(params%path_to_inputs), &
! status='OLD',form='formatted')
!read(default_unit_open,nml=externalPlasmaModel)
!close(default_unit_open)
F%Bfield = Bfield
F%B1field = B1field
F%dBfield = dBfield
F%Bflux = Bflux
F%Bflux3D = Bflux3D
F%Efield = Efield
F%E1field = E1field
F%axisymmetric_fields = axisymmetric_fields
F%Dim2x1t = Dim2x1t
F%ReInterp_2x1t = ReInterp_2x1t
F%t0_2x1t = t0_2x1t
F%ind0_2x1t = ind0_2x1t
F%psip_conv = psip_conv
F%MARS_AMP_Scale = MARS_AMP_Scale
F%Analytic_IWL=Analytic_IWL
F%ntiles=ntiles
F%circumradius=circumradius
F%AORSA_AMP_Scale=AORSA_AMP_Scale
F%AORSA_freq=AORSA_freq
F%AORSA_nmode=AORSA_nmode
F%useLCFS = useLCFS
if (params%proceed.and.F%ReInterp_2x1t) then
call load_prev_iter(params)
F%ind_2x1t=params%prev_iter_2x1t
else
F%ind_2x1t=F%ind0_2x1t
end if
F%E_2x1t = E_2x1t
F%E_profile = E_profile
F%E_dyn = E_dyn
F%E_pulse = E_pulse
F%E_width = E_width
F%AB%a = minor_radius
F%AB%Ro = major_radius
F%PSIp_lim=PSIp_lim
F%res_double=res_double
!write(output_unit_write,'("E_dyn: ",E17.10)') E_dyn
! write(output_unit_write,'("E_pulse: ",E17.10)') E_pulse
! write(output_unit_write,'("E_width: ",E17.10)') E_width
call load_dim_data_from_hdf5(params,F)
!sets F%dims for 2D or 3D data
!write(6,*) F%dims
call which_fields_in_file(params,F%Bfield_in_file,F%Efield_in_file, &
F%Bflux_in_file,F%dBfield_in_file,F%B1field_in_file, &
F%E1field_in_file)
if (F%Bflux.AND..NOT.F%Bflux_in_file) then
write(output_unit_write,'("ERROR: Magnetic flux to be used but no data in file!")')
call KORC_ABORT(18)
end if
if (F%Bfield.AND..NOT.F%Bfield_in_file) then
write(output_unit_write,'("ERROR: Magnetic field to be used but no data in file!")')
call KORC_ABORT(18)
end if
if (F%B1field.AND..NOT.F%B1field_in_file) then
write(output_unit_write,'("ERROR: Magnetic perturbation field to be used but no data in file!")')
call KORC_ABORT(18)
end if
if (F%E1field.AND..NOT.F%E1field_in_file) then
write(output_unit_write,'("ERROR: Electric perturbation field to be used but no data in file!")')
call KORC_ABORT(18)
end if
if (F%dBfield.AND..NOT.F%dBfield_in_file) then
write(output_unit_write,'("ERROR: differential Magnetic field to be used &
but no data in file!")')
call KORC_ABORT(18)
end if
if (F%Efield.AND..NOT.F%Efield_in_file) then
if (params%mpi_params%rank.EQ.0_idef) then
write(output_unit_write,'(/,"* * * * * * * * * * FIELDS * * * * * * * * * *")')
write(output_unit_write,'("MESSAGE: Analytical electric field will be used.")')
write(output_unit_write,'("* * * * * * * * * * * * ** * * * * * * * * * * *",/)')
flush(output_unit_write)
end if
end if
if (F%axisymmetric_fields) then
if (F%Dim2x1t) then
call ALLOCATE_2D_FIELDS_ARRAYS(params,F,F%Bfield, &
F%Bflux,F%dBfield,F%Efield.AND.F%Efield_in_file,F%B1field, &
F%E1field)
call ALLOCATE_3D_FIELDS_ARRAYS(params,F,F%Bfield, &
F%Efield,F%dBfield)
else if (params%orbit_model(1:2).eq.'FO') then
call ALLOCATE_2D_FIELDS_ARRAYS(params,F,F%Bfield, &
.TRUE.,F%dBfield,F%Efield,F%B1field,F%E1field)
else
call ALLOCATE_2D_FIELDS_ARRAYS(params,F,F%Bfield, &
F%Bflux,F%dBfield,F%Efield,F%B1field,F%E1field)
end if
else if ((params%field_model(10:13).eq.'MARS').OR. &
(params%field_model(10:14).eq.'AORSA')) then
call ALLOCATE_2D_FIELDS_ARRAYS(params,F,F%Bfield, &
F%Bflux,F%dBfield,F%Efield,F%B1field,F%E1field)
else
call ALLOCATE_3D_FIELDS_ARRAYS(params,F,F%Bfield,F%Efield,F%dBfield)
end if
!allocates 2D or 3D data arrays (fields and spatial)
call load_field_data_from_hdf5(params,F)
!write(output_unit_write,*) F%PSIp
!write(output_unit_write,*) F%E_3D%PHI(:,F%ind0_2x1t,:)
!write(6,*) F%B1Re_2D%R(:,200)
!flush(6)
! end if
if (F%Bflux.and..not.(params%field_model(10:13).eq.'MARS')) then
F%PSIP_min=minval(F%PSIp)
else if(F%Bflux3D) then
F%PSIP_min=minval(F%PSIp3D(:,1,:))
end if
if ((.not.F%Efield_in_file).and.(.not.F%Dim2x1t).and.F%Efield) then
if (F%Eo.eq.0) F%Eo = Eo
if (F%axisymmetric_fields) then
F%E_2D%R=0._rp
do ii=1_idef,F%dims(1)
F%E_2D%PHI(ii,:)=F%Eo*F%Ro/F%X%R(ii)
end do
F%E_2D%Z=0._rp
else
F%E_3D%R=0._rp
do ii=1_idef,F%dims(1)
F%E_3D%PHI(ii,:,:)=F%Eo*F%Ro/F%X%R(ii)
end do
F%E_3D%Z=0._rp
end if
end if
if(F%dBfield.and..not.F%dBfield_in_file) then
if (F%axisymmetric_fields) then
F%dBdR_2D%R=0._rp
F%dBdR_2D%PHI=0._rp
F%dBdR_2D%Z=0._rp
F%dBdPHI_2D%R=0._rp
F%dBdPHI_2D%PHI=0._rp
F%dBdPHI_2D%Z=0._rp
F%dBdZ_2D%R=0._rp
F%dBdZ_2D%PHI=0._rp
F%dBdZ_2D%Z=0._rp
else
F%dBdR_3D%R=0._rp
F%dBdR_3D%PHI=0._rp
F%dBdR_3D%Z=0._rp
F%dBdPHI_3D%R=0._rp
F%dBdPHI_3D%PHI=0._rp
F%dBdPHI_3D%Z=0._rp
F%dBdZ_3D%R=0._rp
F%dBdZ_3D%PHI=0._rp
F%dBdZ_3D%Z=0._rp
end if
end if
if (params%mpi_params%rank .EQ. 0) then
write(output_unit_write,'("EXTERNAL")')
write(output_unit_write,'("Magnetic field: ",E17.10)') F%Bo
write(output_unit_write,'("Electric field: ",E17.10)') F%Eo
end if
if (params%SC_E) then
F%dim_1D=dim_1D
F%dt_E_SC=dt_E_SC
F%Ip_exp=Ip_exp
write(output_unit_write,*) 'dt_E_SC',F%dt_E_SC,'Ip_exp',Ip_exp
call allocate_1D_FS_arrays(params,F)
call load_1D_FS_from_hdf5(params,F)
! write(output_unit_write,*) F%PSIP_1D
end if
! test=.true.
! if (F%Bflux.and.(.not.test)) then
! call initialize_fields_interpolant(params,F)
! F%Bfield=.TRUE.
! F%Efield=.TRUE.
! F%Efield_in_file=.TRUE.
! RMIN=F%X%R(1)
! RMAX=F%X%R(F%dims(1))
! ZMIN=F%X%Z(1)
! ZMAX=F%X%Z(F%dims(3))
! do ii=1_idef,res_double
! F%dims(1)=2*F%dims(1)-1
! F%dims(3)=2*F%dims(3)-1
! end do
! if (res_double>0) then
! DEALLOCATE(F%X%R)
! DEALLOCATE(F%X%Z)
! DEALLOCATE(F%PSIp)
! end if
! call ALLOCATE_2D_FIELDS_ARRAYS(params,F,F%Bfield, &
! F%Bflux,F%Efield.AND.F%Efield_in_file)
! do ii=1_idef,F%dims(1)
! F%X%R(ii)=RMIN+REAL(ii-1)/REAL(F%dims(1)-1)*(RMAX-RMIN)
! end do
! do ii=1_idef,F%dims(3)
! F%X%Z(ii)=ZMIN+REAL(ii-1)/REAL(F%dims(3)-1)*(ZMAX-ZMIN)
! end do
! call calculate_initial_magnetic_field(F)
! F%E_2D%R=0._rp
! do ii=1_idef,F%dims(1)
! F%E_2D%PHI(ii,:)=F%Eo*F%Ro/F%X%R(ii)
! end do
! F%E_2D%Z=0._rp
! end if
! if (F%Bflux.and.test) then
! F%Bfield=.TRUE.
! F%Efield=.TRUE.
! F%Efield_in_file=.TRUE.
! call ALLOCATE_2D_FIELDS_ARRAYS(params,F,F%Bfield, &
! F%Bflux,F%Efield.AND.F%Efield_in_file)
! B
! edge nodes at minimum R,Z
! F%B_2D%Z(1,:)=-(F%PSIp(2,:)-F%PSIp(1,:))/(F%X%R(2)-F%X%R(1))/F%X%R(1)
! F%B_2D%R(:,1)=(F%PSIp(:,2)-F%PSIp(:,1))/(F%X%Z(2)-F%X%Z(1))/F%X%R(:)
! edge nodes at maximum R,Z
! F%B_2D%Z(F%dims(1),:)=-(F%PSIp(F%dims(1),:)-F%PSIp(F%dims(1)-1,:))/ &
! (F%X%R(F%dims(1))-F%X%R(F%dims(1)-1))/F%X%R(F%dims(1))
! F%B_2D%R(:,F%dims(3))=(F%PSIp(:,F%dims(3))-F%PSIp(:,F%dims(3)-1))/ &
! (F%X%Z(F%dims(3))-F%X%Z(F%dims(3)-1))/F%X%R(:)
! do ii=2_idef,F%dims(1)-1
! central difference over R for interior nodes for BZ
! F%B_2D%Z(ii,:)=-(F%PSIp(ii+1,:)-F%PSIp(ii-1,:))/ &
! (F%X%R(ii+1)-F%X%R(ii-1))/F%X%R(ii)
! end do
! do ii=2_idef,F%dims(3)-1
! central difference over Z for interior nodes for BR
! F%B_2D%R(:,ii)=(F%PSIp(:,ii+1)-F%PSIp(:,ii-1))/ &
! (F%X%Z(ii+1)-F%X%Z(ii-1))/F%X%R(:)
! end do
! do ii=1_idef,F%dims(1)
! F%B_2D%PHI(ii,:)=-F%Bo*F%Ro/F%X%R(ii)
! end do
! F%E_2D%R=0._rp
! do ii=1_idef,F%dims(1)
! F%E_2D%PHI(ii,:)=F%Eo*F%Ro/F%X%R(ii)
! end do
! F%E_2D%Z=0._rp
! F%Bfield=.FALSE.
if (params%mpi_params%rank.EQ.0) then
if (F%axisymmetric_fields) then
if (F%Bflux) then
write(output_unit_write,'("PSIp(r=0)",E17.10)') F%PSIp(F%dims(1)/2,F%dims(3)/2)
write(output_unit_write,'("BPHI(r=0)",E17.10)') F%Bo
write(output_unit_write,'("EPHI(r=0)",E17.10)') F%Eo
else if (F%Bflux3D) then
write(output_unit_write,'("PSIp(r=0)",E17.10)') F%PSIp3D(F%dims(1)/2,1,F%dims(3)/2)
else
write(output_unit_write,'("BR(r=0)",E17.10)') F%B_2D%R(F%dims(1)/2,F%dims(3)/2)
write(output_unit_write,'("BPHI(r=0)",E17.10)') &
F%B_2D%PHI(F%dims(1)/2,F%dims(3)/2)
write(output_unit_write,'("BZ(r=0)",E17.10)') F%B_2D%Z(F%dims(1)/2,F%dims(3)/2)
write(output_unit_write,'("EPHI(r=0)",E17.10)') &
F%E_2D%PHI(F%dims(1)/2,F%dims(3)/2)
end if
end if
end if
if (params%orbit_model(3:5).EQ.'pre') then
if (params%mpi_params%rank.eq.0) then
write(output_unit_write,'("Initializing GC fields from external EM fields")')
end if
if (params%field_model(10:12).eq.'2DB') then
if (F%axisymmetric_fields) then
call initialize_GC_fields(F)
else
call initialize_GC_fields_3D(F)
end if
end if
end if
! write(output_unit_write,'("gradBR",E17.10)') F%gradB_2D%R(F%dims(1)/2,F%dims(3)/2)
! write(output_unit_write,'("gradBPHI",E17.10)') F%gradB_2D%PHI(F%dims(1)/2,F%dims(3)/2)
! write(output_unit_write,'("gradBZ",E17.10)') F%gradB_2D%Z(F%dims(1)/2,F%dims(3)/2)
else
F%Bo=Bo
F%Eo=Eo
if (params%mpi_params%rank.EQ.0) then
write(output_unit_write,'("UNIFORM")')
write(output_unit_write,'("BPHI(r=0)",E17.10)') F%Bo
write(output_unit_write,'("EPHI(r=0)",E17.10)') F%Eo
end if
end if
if (params%mpi_params%rank.eq.0) then
write(output_unit_write,'("* * * * * * * * * * * * * * * * * * * * * * * * *",/)')
end if
end subroutine initialize_fields