GCEoM1_p Subroutine

  subroutine GCEoM1_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_Z,q_cache,m_cache,PSIp,ne)
    TYPE(KORC_PARAMS), INTENT(INOUT)                           :: params
    !! Core KORC simulation parameters.
    TYPE(FIELDS), INTENT(IN)      :: F
    TYPE(PROFILES), INTENT(IN)                                 :: P
    REAL(rp),DIMENSION(8)  :: Bmag,bhat_R,bhat_PHI,bhat_Z,Bst_R,Bst_PHI
    REAL(rp),DIMENSION(8)  :: BstdotE,BstdotgradB,EcrossB_R,EcrossB_PHI,bdotBst
    REAL(rp),DIMENSION(8)  :: bcrossgradB_R,bcrossgradB_PHI,bcrossgradB_Z,gamgc
    REAL(rp),DIMENSION(8)  :: EcrossB_Z,Bst_Z
    REAL(rp),DIMENSION(8)  :: pm,xi,tau_R
    REAL(rp),DIMENSION(8)  :: SR_PLL,SR_MU,BREM_PLL,BREM_MU,BREM_P
    REAL(rp),DIMENSION(8),INTENT(in) :: gradB_R,gradB_PHI,gradB_Z,curlb_R
    REAL(rp),DIMENSION(8),INTENT(in) :: curlb_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z
    REAL(rp),DIMENSION(8),INTENT(OUT) :: RHS_R,RHS_PHI,RHS_Z
    REAL(rp),DIMENSION(8),INTENT(OUT) :: RHS_PLL,RHS_MU
    REAL(rp),DIMENSION(8),INTENT(IN) :: V_PLL,V_MU,Y_R,Y_Z,curlb_PHI
    REAL(rp),DIMENSION(8),INTENT(IN) :: PSIp
    REAL(rp),INTENT(in) :: q_cache,m_cache
    INTEGER(ip)  :: cc
    INTEGER(ip),INTENT(IN)  :: tt
    REAL(rp)  :: time,re_cache,alpha_cache
    REAL(rp), DIMENSION(8) 			:: Te,Zeff
    REAL(rp), DIMENSION(8),INTENT(OUT) 		:: ne
    
    !$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,8
       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)

       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)

       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)


       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)

       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

!    !$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(6,*) 'RHS_PLL',RHS_PLL(1)

       re_cache=C_RE/params%cpp%length
       alpha_cache=C_a
       
       time=params%init_time+(params%it-1+tt)*params%dt
       call analytical_profiles_p(time,params,Y_R,Y_Z,P,F,ne,Te,Zeff,PSIp)
       
       !$OMP SIMD
!       !$OMP& aligned(tau_R,Bmag,RHS_PLL,V_PLL,xi,gamgc,RHS_MU,V_MU)
       do cc=1_idef,8
          
          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

!    !$OMP SIMD
!    do cc=1_idef,8
!       if(isnan(ne(cc))) stop 'ne is a NaN'
!       if(isnan(Zeff(cc))) stop 'Zeff is a NaN'
!       if(isnan(gamgc(cc))) stop 'gamgc is a NaN'
!       if(isnan(BREM_P(cc))) stop 'BREM_P is a NaN'
!       if(isnan(BREM_PLL(cc))) stop 'BREM_PLL is a NaN'
!       if(isnan(BREM_MU(cc))) stop 'BREM_MU is a NaN'
!       if(isnan(SR_PLL(cc))) stop 'SR_PLL is a NaN'
!       if(isnan(SR_MU(cc))) stop 'SR_MU is a NaN'
!       if(isnan(RHS_R(cc))) stop 'RHS_R1 is a NaN'
!       if(isnan(RHS_PHI(cc))) stop 'RHS_PHI1 is a NaN'
!       if(isnan(RHS_Z(cc))) stop 'RHS_Z1 is a NaN'
!       if(isnan(RHS_PLL(cc))) stop 'RHS_PLL1 is a NaN'
!       if(isnan(RHS_MU(cc))) stop 'RHS_MU1 is a NaN'

!    end do
!    !$OMP END SIMD

       
!    write(6,*) 'RHS_R: ',RHS_R(1)
!    write(6,*) 'RHS_PHI: ',RHS_PHI(1)
!    write(6,*) 'RHS_Z: ',RHS_Z(1)
!    write(6,*) 'RHS_PLL: ',RHS_PLL(1)
!    write(6,*) 'RHS_MU: ',RHS_MU(1)
    
  end subroutine GCEoM1_p