GCEoM1_fio_p Subroutine

  subroutine GCEoM1_fio_p(tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
       B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, &
       gradB_R,gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp, &
       ne,ni,nimp,Te,Zeff,flagCon,flagCol,hint)

    USE omp_lib
    IMPLICIT NONE
    
    TYPE(KORC_PARAMS), INTENT(INOUT)                           :: params
    !! Core KORC simulation parameters.
    TYPE(FIELDS), INTENT(IN)      :: F
    TYPE(PROFILES), INTENT(IN)                                 :: P
    REAL(rp),DIMENSION(params%pchunk)  :: Bmag,bhat_R,bhat_PHI,bhat_Z,Bst_R,Bst_PHI
    REAL(rp),DIMENSION(params%pchunk)  :: BstdotE,BstdotgradB,EcrossB_R,EcrossB_PHI,bdotBst
    REAL(rp),DIMENSION(params%pchunk)  :: bcrossgradB_R,bcrossgradB_PHI,bcrossgradB_Z,gamgc
    REAL(rp),DIMENSION(params%pchunk)  :: EcrossB_Z,Bst_Z
    REAL(rp),DIMENSION(params%pchunk)  :: pm,xi,tau_R
    REAL(rp),DIMENSION(params%pchunk)  :: SR_PLL,SR_MU,BREM_PLL,BREM_MU,BREM_P
    REAL(rp),DIMENSION(params%pchunk),INTENT(in) :: gradB_R,gradB_PHI,gradB_Z,curlb_R
    REAL(rp),DIMENSION(params%pchunk),INTENT(in) :: curlb_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z
    REAL(rp),DIMENSION(params%pchunk),INTENT(OUT) :: RHS_R,RHS_PHI,RHS_Z
    REAL(rp),DIMENSION(params%pchunk),INTENT(OUT) :: RHS_PLL,RHS_MU
    REAL(rp),DIMENSION(params%pchunk),INTENT(IN) :: V_PLL,V_MU,Y_R,Y_PHI,Y_Z,curlb_PHI
    REAL(rp),DIMENSION(params%pchunk),INTENT(IN) :: PSIp
    REAL(rp),INTENT(in) :: q_cache,m_cache
    INTEGER(ip)  :: cc,pchunk,ii
    INTEGER(ip),INTENT(IN)  :: tt
    REAL(rp)  :: time,re_cache,alpha_cache
    REAL(rp), DIMENSION(params%pchunk),INTENT(OUT) 		:: ne,Te,ni,nimp,Zeff
    TYPE(C_PTR), DIMENSION(params%pchunk), INTENT(INOUT)  :: hint
    INTEGER(is),DIMENSION(params%pchunk),intent(INOUT) :: flagCon,flagCol
    INTEGER             :: thread_num

    pchunk=params%pchunk

    thread_num = OMP_GET_THREAD_NUM()
    
!    !$OMP SIMD
!    do cc=1_idef,pchunk 
!       if(isnan(B_R(cc))) then
!          write(6,*) thread_num,'0Y',Y_R(cc)*params%cpp%length,Y_PHI(cc),Y_Z(cc)*params%cpp%length
!          write(6,*) thread_num,'0B',B_R(cc),B_PHI(cc),B_Z(cc)
!          write(6,*) thread_num,'0E',E_R(cc),E_PHI(cc),E_Z(cc)
!          write(6,*) thread_num,'0gradB',gradB_R(cc),gradB_PHI(cc),gradB_Z(cc)
!          write(6,*) thread_num,'0curlb',curlb_R(cc),curlb_PHI(cc),curlb_Z(cc)
!          write(6,*) thread_num,'0V',V_PLL(cc),V_MU(cc)
!          write(6,*) 'Exb',Ecrossb_R(cc)
!          write(6,*) 'bxgradB',bcrossgradB_R(cc)
!          write(6,*) 'Bst',Bst_R(cc)
!          write(6,*) 'bdotBst',bdotBst(cc)
!          write(6,*) thread_num,'0gamma',gamgc(cc)

!          stop 'B_R is a NaN'
!       endif
!    end do
!    !$OMP END SIMD

    
    !$OMP SIMD
    !    !$OMP& aligned(gradB_R,gradB_PHI,gradB_Z,curlb_R,curlb_Z, &
    !    !$OMP& B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, &
    !    !$OMP& V_PLL,V_MU,Y_R,curlb_PHI,tau_R)
    do cc=1_idef,pchunk

       ne(cc)=1._rp
       Te(cc)=1._rp
       Zeff(cc)=1._rp
       
       Bmag(cc) = SQRT(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+B_Z(cc)*B_Z(cc))
       
       bhat_R(cc) = B_R(cc)/Bmag(cc)
       bhat_PHI(cc) = B_PHI(cc)/Bmag(cc)
       bhat_Z(cc) = B_Z(cc)/Bmag(cc)

       Bst_R(cc)=q_cache*B_R(cc)+V_PLL(cc)*curlb_R(cc)
       Bst_PHI(cc)=q_cache*B_PHI(cc)+V_PLL(cc)*curlb_PHI(cc)
       Bst_Z(cc)=q_cache*B_Z(cc)+V_PLL(cc)*curlb_Z(cc)

      ! write(output_unit_write,*) 'bmag',Bmag(cc),'bhat',bhat_R(cc),bhat_PHI(cc),bhat_Z(cc),'Bst',Bst_R(cc),Bst_PHI(cc),Bst_Z(cc)
       
       bdotBst(cc)=bhat_R(cc)*Bst_R(cc)+bhat_PHI(cc)*Bst_PHI(cc)+ &
            bhat_Z(cc)*Bst_Z(cc)
       BstdotE(cc)=Bst_R(cc)*E_R(cc)+Bst_PHI(cc)*E_PHI(cc)+Bst_Z(cc)*E_Z(cc)   
       BstdotgradB(cc)=Bst_R(cc)*gradB_R(cc)+Bst_PHI(cc)*gradB_PHI(cc)+ &
            Bst_Z(cc)*gradB_Z(cc)

       !write(output_unit_write,*) 'bdotBst',bdotBst(cc),BstdotE(cc),BstdotgradB(cc)
       
       Ecrossb_R(cc)=E_PHI(cc)*bhat_Z(cc)-E_Z(cc)*bhat_PHI(cc)
       Ecrossb_PHI(cc)=E_Z(cc)*bhat_R(cc)-E_R(cc)*bhat_Z(cc)
       Ecrossb_Z(cc)=E_R(cc)*bhat_PHI(cc)-E_PHI(cc)*bhat_R(cc)

       !write(output_unit_write,*) 'Ecrossb',Ecrossb_R(cc),Ecrossb_PHI(cc),Ecrossb_Z(cc)
       
       bcrossgradB_R(cc)=bhat_PHI(cc)*gradB_Z(cc)-bhat_Z(cc)*gradB_PHI(cc)
       bcrossgradB_PHI(cc)=bhat_Z(cc)*gradB_R(cc)-bhat_R(cc)*gradB_Z(cc)
       bcrossgradB_Z(cc)=bhat_R(cc)*gradB_PHI(cc)-bhat_PHI(cc)*gradB_R(cc)

      ! write(output_unit_write,*) 'bcrossgradB',bcrossgradB_R(cc),bcrossgradB_PHI(cc),bcrossgradB_Z(cc)
       
       gamgc(cc)=sqrt(1+V_PLL(cc)*V_PLL(cc)+2*V_MU(cc)*Bmag(cc))
       
       pm(cc)=sqrt(gamgc(cc)**2-1)
       xi(cc)=V_PLL(cc)/pm(cc)

       RHS_R(cc)=(q_cache*Ecrossb_R(cc)+(m_cache*V_MU(cc)* &
            bcrossgradB_R(cc)+V_PLL(cc)*Bst_R(cc))/(m_cache*gamgc(cc)))/ &
            bdotBst(cc)
       RHS_PHI(cc)=(q_cache*Ecrossb_PHI(cc)+(m_cache*V_MU(cc)* &
            bcrossgradB_PHI(cc)+V_PLL(cc)*Bst_PHI(cc))/(m_cache*gamgc(cc)))/ &
            (Y_R(cc)*bdotBst(cc))
       RHS_Z(cc)=(q_cache*Ecrossb_Z(cc)+(m_cache*V_MU(cc)* &
            bcrossgradB_Z(cc)+V_PLL(cc)*Bst_Z(cc))/(m_cache*gamgc(cc)))/ &
            bdotBst(cc)
       RHS_PLL(cc)=(q_cache*BstdotE(cc)-V_MU(cc)*BstdotgradB(cc)/gamgc(cc))/ &
            bdotBst(cc)
       RHS_MU(cc)=0._rp

    end do
    !$OMP END SIMD

    !write(output_unit_write,*) 'bmag',Bmag(1),'bhat',bhat_R(1),bhat_PHI(1),bhat_Z(1),'Bst',Bst_R(1),Bst_PHI(1),Bst_Z(1)
    !write(output_unit_write,*) 'bdotBst',bdotBst(1),BstdotE(1),BstdotgradB(1)
    !write(output_unit_write,*) 'Ecrossb',Ecrossb_R(1),Ecrossb_PHI(1),Ecrossb_Z(1)
    !write(output_unit_write,*) 'bcrossgradB',bcrossgradB_R(1),bcrossgradB_PHI(1),bcrossgradB_Z(1)

    !    !$OMP SIMD
    !    do cc=1_idef,8
    !       if(isnan(RHS_R(cc))) stop 'RHS_R0 is a NaN'
    !       if(isnan(RHS_PHI(cc))) stop 'RHS_PHI0 is a NaN'
    !       if(isnan(RHS_Z(cc))) stop 'RHS_Z0 is a NaN'
    !       if(isnan(RHS_PLL(cc))) stop 'RHS_PLL0 is a NaN'
    !       if(isnan(RHS_MU(cc))) stop 'RHS_MU0 is a NaN'
    !    end do
    !    !$OMP END SIMD

    if (params%radiation.and.(params%GC_rad_model.eq.'SDE')) then

       !       write(output_unit_write,*) 'RHS_PLL',RHS_PLL(1)

       re_cache=C_RE/params%cpp%length
       alpha_cache=C_a

       call get_fio_profile_p(params,P,Y_R,Y_PHI,Y_Z, &
            ne,Te,flagCon,hint)

       call get_fio_ion_p(params,P,Y_R,Y_PHI,Y_Z, &
            ne,ni,nimp,Zeff,flagCon,hint)
          
       !$OMP SIMD
       !       !$OMP& aligned(tau_R,Bmag,RHS_PLL,V_PLL,xi,gamgc,RHS_MU,V_MU)
       do cc=1_idef,pchunk

          tau_R(cc)=6*C_PI*E0/(Bmag(cc)*Bmag(cc))

          SR_PLL(cc)=V_PLL(cc)*(1._rp-xi(cc)*xi(cc))/tau_R(cc)* &
               (1._rp/gamgc(cc)-gamgc(cc))
          SR_MU(cc)=-2._rp*V_MU(cc)/tau_R(cc)* &
               (gamgc(cc)*(1-xi(cc)*xi(cc))+xi(cc)*xi(cc)/gamgc(cc))
          
          !Normalizations done here
          BREM_P(cc)=-4._rp*re_cache**2*ne(cc)* &
               Zeff(cc)*(Zeff(cc)+1._rp)*alpha_cache* &
               (gamgc(cc)-1._rp)*(log(2._rp*gamgc(cc))-1._rp/3._rp)
          BREM_PLL(cc)=xi(cc)*BREM_P(cc)
          BREM_MU(cc)=(1._rp-xi(cc)*xi(cc))*V_PLL(cc)/ &
               (Bmag(cc)*xi(cc))*BREM_P(cc)

          RHS_PLL(cc)=RHS_PLL(cc)+SR_PLL(cc)+BREM_PLL(cc)
          RHS_MU(cc)=SR_MU(cc)+BREM_MU(cc)

       end do
       !$OMP END SIMD


    end if

#if DBG_CHECK    
    !$OMP SIMD
    do cc=1_idef,pchunk 
       if(isnan(gamgc(cc)).and. &
            ((flagCon(cc)==1).and.(flagCol(cc)==1))) then           
          write(6,*) 'gamgc is NaN'
          write(6,*) 'V_PLL',V_PLL(cc)
          write(6,*) 'V_MU',V_MU(cc)
          write(6,*) 'Bmag',Bmag(cc)
          
          stop 'gamgc is a NaN'
       endif
       if(isnan(RHS_R(cc)).and. &
            ((flagCon(cc)==1).and.(flagCol(cc)==1))) then  

          write(6,*) thread_num,'flags',flagCon(cc),flagCol(cc)
          write(6,*) thread_num,'Y',Y_R(cc)*params%cpp%length,Y_PHI(cc),Y_Z(cc)*params%cpp%length
          write(6,*) thread_num,'B',B_R(cc),B_PHI(cc),B_Z(cc)
          write(6,*) thread_num,'E',E_R(cc),E_PHI(cc),E_Z(cc)
          write(6,*) thread_num,'gradB',gradB_R(cc),gradB_PHI(cc),gradB_Z(cc)
          write(6,*) thread_num,'curlb',curlb_R(cc),curlb_PHI(cc),curlb_Z(cc)
          write(6,*) thread_num,'V',V_PLL(cc),V_MU(cc)
          write(6,*) thread_num,'Exb',Ecrossb_R(cc)
          write(6,*) thread_num,'bxgradB',bcrossgradB_R(cc)
          write(6,*) thread_num,'Bst',Bst_R(cc)
          write(6,*) thread_num,'bdotBst',bdotBst(cc)
          write(6,*) thread_num,'gamma',gamgc(cc)
          
          stop 'RHS_R1 is a NaN'
       endif
       if(isnan(RHS_PHI(cc)).and. &
            ((flagCon(cc)==1).and.(flagCol(cc)==1))) then
          stop 'RHS_PHI1 is a NaN'
       end if
       if(isnan(RHS_Z(cc)).and. &
            ((flagCon(cc)==1).and.(flagCol(cc)==1))) then
          stop 'RHS_Z1 is a NaN'
       endif
       if(isnan(RHS_PLL(cc)).and. &
            ((flagCon(cc)==1).and.(flagCol(cc)==1))) then
          stop 'RHS_PLL1 is a NaN'
       endif
       if(isnan(RHS_MU(cc)).and. &
            ((flagCon(cc)==1).and.(flagCol(cc)==1))) then
          stop 'RHS_MU1 is a NaN'
       endif
    end do
    !$OMP END SIMD

    if (params%radiation.and.(params%GC_rad_model.eq.'SDE')) then
       !$OMP SIMD
       do cc=1_idef,pchunk
          if(isnan(ne(cc)).and. &
               ((flagCon(cc)==1).and.(flagCol(cc)==1))) then
             stop 'ne is a NaN'
          endif
          if(isnan(Zeff(cc)).and. &
               ((flagCon(cc)==1).and.(flagCol(cc)==1))) then
             stop 'Zeff is a NaN'
          endif
          if(isnan(BREM_P(cc)).and. &
               ((flagCon(cc)==1).and.(flagCol(cc)==1))) then
             stop 'BREM_P is a NaN'
          endif
          if(isnan(BREM_PLL(cc)).and. &
               ((flagCon(cc)==1).and.(flagCol(cc)==1))) then
             stop 'BREM_PLL is a NaN'
          endif
          if(isnan(BREM_MU(cc)).and. &
               ((flagCon(cc)==1).and.(flagCol(cc)==1))) then
             stop 'BREM_MU is a NaN'
          endif
          if(isnan(SR_PLL(cc)).and. &
               ((flagCon(cc)==1).and.(flagCol(cc)==1))) then
             stop 'SR_PLL is a NaN'
          endif
          if(isnan(SR_MU(cc)).and. &
               ((flagCon(cc)==1).and.(flagCol(cc)==1))) then
             stop 'SR_MU is a NaN'
          endif
       end do
       !$OMP END SIMD
    end if
#endif
    !    write(output_unit_write,*) 'RHS_R: ',RHS_R(1)
    !    write(output_unit_write,*) 'RHS_PHI: ',RHS_PHI(1)
    !    write(output_unit_write,*) 'RHS_Z: ',RHS_Z(1)
    !    write(output_unit_write,*) 'RHS_PLL: ',RHS_PLL(1)
    !    write(output_unit_write,*) 'RHS_MU: ',RHS_MU(1)

  end subroutine GCEoM1_fio_p