compute_effective_field Subroutine

subroutine compute_effective_field(me)
  type(strand_t), intent(inout) :: me

  real(kind=dp) :: dBdx,dBdy,T,Jop,E0,n,Area,Jc,B0,eps
  real(kind=dp) :: fIntegral,dTheta,dR,Theta,R,xR,yR,f,A,E,JcEq
  real(kind=dp) :: tolerance,Bmin,Bmax,BError,BEff
  real(kind=dp) :: r0,emax,St,Tc0,Bc
  integer :: i,nTheta,nR,iTheta,iR,Iterat
  logical :: converged

  E0   = me%mat_supc%E0
  n    = me%mat_supc%nPow
  Area = me%ct(SUPC)%area
  eps=1.0e-8_dp

  do i=1,me%SC_Prop%NbCells
    T=me%StVar%SCtemp(i) ! Solid strand temperature
    Jop=abs(me%scen%ElCur)/Area !  Current density
    dBdx=me%scen%dBfield(i)
    dBdy=0.0_dp
    B0=me%scen%Bfield(i)

    if(Jop<=0.0_dp) cycle

    St  = me%mat_supc%strain(me%scen%Bfield(i),me%scen%ElCur)
    Tc0 = me%mat_supc%critical_temperature(0.0_dp,St)
    Bc  = me%mat_supc%critical_field(T,St)
    Jc  = me%mat_supc%critical_current_density(T,me%scen%Bfield(i),St,Tc0,Bc)
    if(Jc<=0.0_dp) cycle

    ! Integrate the electric field
    fIntegral = 0.0_dp
    nTheta    = 5
    nR        = 5
    dTheta    = 2.0_dp*Pi_value/dfloat(nTheta)
    dR        = (me%SC_Prop%outer_rad-me%SC_Prop%inner_rad)/dfloat(nR)
    r0=1.0e7_dp
	  emax=r0*Jop
    do iTheta = 1,nTheta
      ! compute angle
      Theta  = 2.0_dp*Pi_value*dfloat(iTheta-1)/dfloat(nTheta-1)
      do iR  = 1,nR
        ! compute radius
        R = dfloat(iR-1)/dfloat(nR-1)
        R = me%SC_Prop%inner_rad*(1.0_dp-R) + me%SC_Prop%outer_rad*R
        ! compute x and y locations in the cross section
        xR = R*cos(Theta)
        yR = R*sin(Theta)
        ! compute B at the location
        me%scen%Bfield(i) = B0 + dBdx*xR + dBdy*yR
        ! compute Jc
        St  = me%mat_supc%strain(me%scen%Bfield(i),me%scen%ElCur)
        Tc0 = me%mat_supc%critical_temperature(0.0_dp,St)
        Bc  = me%mat_supc%critical_field(T,St)
        Jc  = me%mat_supc%critical_current_density(T,me%scen%Bfield(i),St,Tc0,Bc)        
        ! compute the local electric field checking for normal state
        if(Jc>0.0_dp) then
           f = E0*(Jop/Jc)**n
        else
           f=emax
          ! print*, 'critical current calculation failed'
          ! print*, iTheta,iR
          ! print*, Jc,i,T,me%scen%Bfield(i)
          ! print*, 'failed in compute_effective_field'
          ! read(*,*)
        endif
        ! integrate adding contributions of single area differentials
        fIntegral = fIntegral+f*dR*R*dTheta
      enddo
    enddo
    ! average electric field from normalised integral
    A = Pi_value*(me%SC_Prop%outer_rad**2-me%SC_Prop%inner_rad**2)
    E = fIntegral/A

    ! trap limiting case of zero electric field
    if(E<=eps*E0) then
      me%scen%Bfield(i)=B0
      cycle
    endif

    ! find equivalent Jc, corresponding to the computed electric field
    JcEq = Jop/(E/E0)**(1.0/n)

    ! find iteratively equivalent magnetic field
    converged = .false.
    Tolerance = eps
    Bmin      = B0
    Bmax      = B0 + sqrt(dBdx**2+dBdy**2)*me%SC_Prop%outer_rad
    Iterat = 0
    do while(.not.converged)
      Iterat = Iterat + 1
      BEff = 0.5*(BMin+BMax)
      if(Iterat>50) then
        call set_error('convergence failure in compute_effective_field')
        return
      endif
      me%scen%Bfield(i)=BEff
      St  = me%mat_supc%strain(me%scen%Bfield(i),me%scen%ElCur)
      Tc0 = me%mat_supc%critical_temperature(0.0_dp,St)
      Bc  = me%mat_supc%critical_field(T,St)
      Jc  = me%mat_supc%critical_current_density(T,me%scen%Bfield(i),St,Tc0,Bc)
      if(Jc>JcEq) then
         Bmin=BEff
      elseif(Jc<JcEq) then
         Bmax=BEff
      else
         Bmin=BEff
         Bmax=BEff
      endif
      BError    = (Bmax-Bmin)/BEff
      converged = BError<=Tolerance
    enddo
    me%scen%Bfield(i)=BEff
  enddo
end subroutine compute_effective_field