reinit_SC_E1D_FS Subroutine

public subroutine reinit_SC_E1D_FS(params, F)

Arguments

Type IntentOptional AttributesName
type(KORC_PARAMS), intent(in) :: params
type(FIELDS), intent(inout) :: F

Calls

proc~~reinit_sc_e1d_fs~~CallsGraph proc~reinit_sc_e1d_fs reinit_SC_E1D_FS proc~initialize_sc1d_field_interpolant_fs initialize_SC1D_field_interpolant_FS proc~reinit_sc_e1d_fs->proc~initialize_sc1d_field_interpolant_fs ezspline_init ezspline_init proc~initialize_sc1d_field_interpolant_fs->ezspline_init ezspline_allocated ezspline_allocated proc~initialize_sc1d_field_interpolant_fs->ezspline_allocated ezspline_error ezspline_error proc~initialize_sc1d_field_interpolant_fs->ezspline_error ezspline_setup ezspline_setup proc~initialize_sc1d_field_interpolant_fs->ezspline_setup

Called by

proc~~reinit_sc_e1d_fs~~CalledByGraph proc~reinit_sc_e1d_fs reinit_SC_E1D_FS program~main main program~main->proc~reinit_sc_e1d_fs

Contents

Source Code


Source Code

  subroutine reinit_SC_E1D_FS(params,F)
 
    TYPE(FIELDS), INTENT(INOUT)                 :: F
    TYPE(KORC_PARAMS), INTENT(IN) 		:: params
    real(rp),dimension(F%dim_1D) :: Jsamall,Jexp,dJdt,PSIP_1D
    real(rp),dimension(F%dim_1D) :: a,b,c,u,gam,r,alpha,beta,gamma
    real(rp) :: dPSIP,Isam,bet
    integer :: pp,ii,PSIPind
    INTEGER 				:: mpierr

!    if (params%mpi_params%rank .EQ. 0) then
!       write(output_unit_write,*) 'Calculating SC_E1D'
!    end if

    PSIP_1D=F%PSIP_1D
    dPSIP=PSIP_1D(2)-PSIP_1D(1)
    Jsamall=F%J0_SC_1D%PHI

    Isam=0._rp
    do ii=1_idef,F%dim_1D
       if ((ii.eq.1).or.(ii.eq.F%dim_1D)) then
          Isam=Isam+Jsamall(ii)/2._rp
       else 
          Isam=Isam+Jsamall(ii)
       end if
    end do
    Isam=Isam*dPSIP
!    write(output_unit_write,*) params%mpi_params%rank,'Isam: ',Isam

    F%Ip0=F%Ip_exp/Isam
    
    Jexp=Jsamall*F%Ip0

    F%J3_SC_1D%PHI=Jexp 
    F%J2_SC_1D%PHI=Jexp
    F%J1_SC_1D%PHI=Jexp        
    
    ! Calculating time-derivative of E_phi
    
    dJdt=(3._rp*F%J1_SC_1D%PHI-4._rp*F%J2_SC_1D%PHI+F%J3_SC_1D%PHI)/ &
         (2._rp*F%dt_E_SC)
    
!    write(output_unit_write,*) params%mpi_params%rank,'J(1)',F%J_SC_1D%PHI(1)
    
    ! Solving 1D Poisson equation with tridiagonal matrix solve

    alpha=F%ddMagPsiSqdPsiPSq
    beta=F%dMagPsiSqdPsiP
    gamma=C_MU*dJdt
    
    
    a=-alpha*dPSIP/2._rp+beta
    b=-2._rp*beta
    c=alpha*dPSIP/2._rp+beta
    u=0._rp
    gam=0._rp
!    r=-2*dr**2*C_MU*Jexp
    r=dPSIP**2*gamma

    c(2)=c(2)-a(2)*a(1)/c(1)
    b(2)=b(2)-a(2)*b(1)/c(1)
    r(2)=r(2)-a(2)*r(1)/c(1)
    
    bet=b(2)
    u(2)=r(2)/bet
    do ii=3_idef,F%dim_1D-1
       gam(ii)=c(ii-1)/bet
       bet=b(ii)-a(ii)*gam(ii)
       if (bet.eq.0) then
          stop 'tridiag failed'         
       end if
       u(ii)=(r(ii)-a(ii)*u(ii-1))/bet
    end do

    do ii=F%dim_1D-2,2,-1
       u(ii)=u(ii)-gam(ii+1)*u(ii+1)
    end do

    u(1)=2*u(2)-u(3)


    F%E_SC_1D%PHI=u

    if (params%mpi_params%rank.eq.0) then
       write(output_unit_write,*) 'J1(1)',F%J1_SC_1D%PHI(1)
       write(output_unit_write,*) 'J2(1)',F%J2_SC_1D%PHI(1)
       write(output_unit_write,*) 'J3(1)',F%J3_SC_1D%PHI(1)
       
       write(output_unit_write,*) 'E(1)',F%E_SC_1D%PHI(1)
    end if
       
    ! Normalizing inductive E_phi
    
    F%E_SC_1D%PHI=F%E_SC_1D%PHI/params%cpp%Eo

#ifdef PSPLINE
    call initialize_SC1D_field_interpolant_FS(params,F)
#endif
    
  end subroutine reinit_SC_E1D_FS