Loading [MathJax]/jax/output/HTML-CSS/jax.js

initialize_fields Subroutine

public subroutine initialize_fields(params, F)

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.

Arguments

Type IntentOptional AttributesName
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.


Calls

proc~~initialize_fields~~CallsGraph proc~initialize_fields initialize_fields proc~allocate_2d_fields_arrays ALLOCATE_2D_FIELDS_ARRAYS proc~initialize_fields->proc~allocate_2d_fields_arrays proc~allocate_3d_fields_arrays ALLOCATE_3D_FIELDS_ARRAYS proc~initialize_fields->proc~allocate_3d_fields_arrays proc~load_prev_iter load_prev_iter proc~initialize_fields->proc~load_prev_iter proc~allocate_1d_fs_arrays allocate_1D_FS_arrays proc~initialize_fields->proc~allocate_1d_fs_arrays proc~load_1d_fs_from_hdf5 load_1D_FS_from_hdf5 proc~initialize_fields->proc~load_1d_fs_from_hdf5 proc~load_dim_data_from_hdf5 load_dim_data_from_hdf5 proc~initialize_fields->proc~load_dim_data_from_hdf5 proc~which_fields_in_file which_fields_in_file proc~initialize_fields->proc~which_fields_in_file proc~korc_abort korc_abort proc~initialize_fields->proc~korc_abort proc~load_field_data_from_hdf5 load_field_data_from_hdf5 proc~initialize_fields->proc~load_field_data_from_hdf5 proc~initialize_gc_fields initialize_GC_fields proc~initialize_fields->proc~initialize_gc_fields proc~initialize_gc_fields_3d initialize_GC_fields_3D proc~initialize_fields->proc~initialize_gc_fields_3d proc~allocate_v_field_2d ALLOCATE_V_FIELD_2D proc~allocate_2d_fields_arrays->proc~allocate_v_field_2d proc~allocate_v_field_2dx ALLOCATE_V_FIELD_2DX proc~allocate_2d_fields_arrays->proc~allocate_v_field_2dx proc~allocate_v_field_3d ALLOCATE_V_FIELD_3D proc~allocate_3d_fields_arrays->proc~allocate_v_field_3d proc~load_prev_iter->proc~korc_abort mpi_bcast mpi_bcast proc~load_prev_iter->mpi_bcast h5fclose_f h5fclose_f proc~load_prev_iter->h5fclose_f h5fopen_f h5fopen_f proc~load_prev_iter->h5fopen_f interface~load_from_hdf5 load_from_hdf5 proc~load_prev_iter->interface~load_from_hdf5 proc~allocate_1d_fs_arrays->h5fclose_f proc~allocate_1d_fs_arrays->h5fopen_f proc~allocate_1d_fs_arrays->interface~load_from_hdf5 proc~load_1d_fs_from_hdf5->h5fclose_f interface~load_array_from_hdf5 load_array_from_hdf5 proc~load_1d_fs_from_hdf5->interface~load_array_from_hdf5 proc~load_1d_fs_from_hdf5->h5fopen_f proc~load_dim_data_from_hdf5->h5fclose_f proc~load_dim_data_from_hdf5->h5fopen_f proc~load_dim_data_from_hdf5->interface~load_from_hdf5 proc~which_fields_in_file->h5fclose_f proc~which_fields_in_file->h5fopen_f mpi_abort mpi_abort proc~korc_abort->mpi_abort proc~load_field_data_from_hdf5->h5fclose_f proc~load_field_data_from_hdf5->interface~load_array_from_hdf5 proc~load_field_data_from_hdf5->h5fopen_f proc~load_field_data_from_hdf5->interface~load_from_hdf5 proc~rload_1d_array_from_hdf5 rload_1d_array_from_hdf5 interface~load_array_from_hdf5->proc~rload_1d_array_from_hdf5 proc~rload_2d_array_from_hdf5 rload_2d_array_from_hdf5 interface~load_array_from_hdf5->proc~rload_2d_array_from_hdf5 proc~rload_3d_array_from_hdf5 rload_3d_array_from_hdf5 interface~load_array_from_hdf5->proc~rload_3d_array_from_hdf5 proc~iload_from_hdf5 iload_from_hdf5 interface~load_from_hdf5->proc~iload_from_hdf5 proc~rload_from_hdf5 rload_from_hdf5 interface~load_from_hdf5->proc~rload_from_hdf5 proc~iload_from_hdf5->proc~korc_abort h5dclose_f h5dclose_f proc~iload_from_hdf5->h5dclose_f h5dread_f h5dread_f proc~iload_from_hdf5->h5dread_f proc~rload_from_hdf5->proc~korc_abort proc~rload_from_hdf5->h5dclose_f proc~rload_from_hdf5->h5dread_f proc~rload_1d_array_from_hdf5->proc~korc_abort proc~rload_1d_array_from_hdf5->h5dclose_f proc~rload_1d_array_from_hdf5->h5dread_f proc~rload_2d_array_from_hdf5->proc~korc_abort proc~rload_2d_array_from_hdf5->h5dclose_f proc~rload_2d_array_from_hdf5->h5dread_f proc~rload_3d_array_from_hdf5->proc~korc_abort proc~rload_3d_array_from_hdf5->h5dclose_f proc~rload_3d_array_from_hdf5->h5dread_f

Called by

proc~~initialize_fields~~CalledByGraph proc~initialize_fields initialize_fields program~main main program~main->proc~initialize_fields

Contents

Source Code


Source Code

  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