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 Bϕ<0 and Bθ<0. Sign convention in analytical fields corresponds to DIII-D fields with Bϕ<0 and Bθ<0.
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 redge as measured from the magnetic axis. Radial position of the magnetic axis R0 Safety factor at the plasma edge. Safety factor at the magnetic axis q0. 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