check_if_in_fields_domain_p Subroutine

  subroutine check_if_in_fields_domain_p(pchunk,F,Y_R,Y_PHI,Y_Z,flag)
    !! @note Subrotuine that checks if particles in the simulation are within
    !! the spatial domain where interpolants and fields are known. @endnote
    !! External fields and interpolants can have different spatial domains where
    !! they are defined. Therefore, it is necessary to
    !! check if a given particle has left these spatial domains to
    !! stop following it, otherwise this will cause an error in the simulation.
    INTEGER, INTENT(IN)  :: pchunk
    TYPE(FIELDS), INTENT(IN)                                   :: F
    REAL(rp), DIMENSION(pchunk),  INTENT(IN)      :: Y_R,Y_PHI,Y_Z    
    INTEGER(is), DIMENSION(pchunk), INTENT(INOUT)  :: flag
    !! Flag that determines whether particles are followed in the
    !! simulation (flag=1), or not (flag=0).
    INTEGER                                                :: IR
    !! Variable used to localize the grid cell in the \((R,\phi,Z)\)
    !! or \((R,Z)\) grid containing the fields data that corresponds
    !! to the radial position of the particles.
    INTEGER                                                :: IPHI
    !! Variable used to localize the grid cell in the \((R,\phi,Z)\)
    !! or \((R,Z)\) grid containing the fields data that corresponds
    !! to the azimuthal position of the particles.
    INTEGER                                                :: IZ
    !! Variable used to localize the grid cell in the \((R,\phi,Z)\)
    !! or \((R,Z)\) grid containing the fields data that corresponds
    !! to the vertical position of the particles.
    INTEGER(ip)                                            :: pp
    !! Particle iterator.
    INTEGER(ip)                                            :: ss
    !! Species iterator.
    REAL(rp) :: Rwall

    
!    write(output_unit_write,'("YR:",E17.10)') Y_R
!    write(output_unit_write,'("YPHI:",E17.10)') Y_PHI
!    write(output_unit_write,'("YZ:",E17.10)') Y_Z

!    write(output_unit_write,'("Ro:",E17.10)') fields_domain%Ro
!    write(output_unit_write,'("Zo:",E17.10)') fields_domain%Zo
!    write(output_unit_write,'("DR:",E17.10)') fields_domain%DR
    !    write(output_unit_write,'("DZ:",E17.10)') fields_domain%DZ
!    write(output_unit_write,'("DT:",E17.10)') fields_domain%DT


    
    if (ALLOCATED(fields_domain%FLAG3D)) then
       if (F%Dim2x1t) then
          !$OMP SIMD
          !       !$OMP&  aligned(IR,IPHI,IZ)
          do pp=1_idef,pchunk

             IR = INT(FLOOR((Y_R(pp)  - fields_domain%Ro + &
                  0.5_rp*fields_domain%DR)/fields_domain%DR) + 1.0_rp,idef)

             IZ = INT(FLOOR((Y_Z(pp)  + ABS(fields_domain%Zo) + &
                  0.5_rp*fields_domain%DZ)/fields_domain%DZ) + 1.0_rp,idef)

!             write(output_unit_write,'("IR: ",I16)') IR
!             write(output_unit_write,'("IPHI: ",I16)') IPHI
!             write(output_unit_write,'("IZ: ",I16)') IZ
             
             if ((fields_domain%FLAG3D(IR,1,IZ).NE.1_is).OR. &
                  ((IR.GT.bfield_2X1T%NR).OR.(IZ.GT.bfield_2X1T%NZ))) then
                if (F%Analytic_IWL.eq.'NONE') then
                   flag(pp) = 0_is
                else if (F%Analytic_IWL.eq.'D3D') then
                   if ((IR.lt.floor(bfield_2d%NR/6._rp)).and. &
                        (IZ.gt.floor(bfield_2d%NZ/5._rp)).and. &
                        (IZ.lt.floor(4._rp*bfield_2d%NZ/5._rp))) then

                      Rwall=F%circumradius*cos(C_PI/F%ntiles)/ &
                           (cos((mod(Y_PHI(pp),2*C_PI/F%ntiles))-C_PI/F%ntiles))
                      if (Y_R(pp).lt.Rwall) flag(pp) = 0_is

                   else
                      flag(pp) = 0_is
                   endif
                endif

                !write(output_unit_write,'("YR:",E17.10)') Y_R(pp)
                !write(output_unit_write,'("YPHI:",E17.10)') Y_PHI(pp)
                !write(output_unit_write,'("YZ:",E17.10)') Y_Z(pp)

                !write(output_unit_write,'("IR: ",I16)') IR
                !write(output_unit_write,'("IPHI: ",I16)') IPHI
                !write(output_unit_write,'("IZ: ",I16)') IZ

                !call KORC_ABORT()

             end if

             !write(output_unit_write,'("IPHI: ",I16)') IPHI
             !write(output_unit_write,'("flag: ",I16)') flag(pp)


          end do
          !$OMP END SIMD
       else
          !$OMP SIMD
          !       !$OMP&  aligned(IR,IPHI,IZ)
          do pp=1_idef,pchunk

             IR = INT(FLOOR((Y_R(pp)  - fields_domain%Ro + &
                  0.5_rp*fields_domain%DR)/fields_domain%DR) + 1.0_rp,idef)
             IPHI = INT(FLOOR((Y_PHI(pp)  + 0.5_rp*fields_domain%DPHI)/ &
                  fields_domain%DPHI) + 1.0_rp,idef)
             IZ = INT(FLOOR((Y_Z(pp)  + ABS(fields_domain%Zo) + &
                  0.5_rp*fields_domain%DZ)/fields_domain%DZ) + 1.0_rp,idef)

             if ((fields_domain%FLAG3D(IR,IPHI,IZ).NE.1_is).OR. &
                  ((IR.GT.bfield_3d%NR).OR.(IZ.GT.bfield_3d%NZ))) then
                flag(pp) = 0_is

                !write(output_unit_write,'("YR:",E17.10)') Y_R
                !write(output_unit_write,'("YPHI:",E17.10)') Y_PHI
                !write(output_unit_write,'("YZ:",E17.10)') Y_Z

                !write(output_unit_write,'("IR: ",I16)') IR
                !write(output_unit_write,'("IPHI: ",I16)') IPHI
                !write(output_unit_write,'("IZ: ",I16)') IZ

                !call KORC_ABORT()

             end if

             !write(output_unit_write,'("IPHI: ",I16)') IPHI
             !write(output_unit_write,'("flag: ",I16)') flag(pp)


          end do
          !$OMP END SIMD
       end if
    else if (ALLOCATED(fields_domain%FLAG2D)) then
       !$OMP SIMD
!       !$OMP& aligned(IR,IZ)
       do pp=1_idef,pchunk

          !write(output_unit_write,*) Y_R(pp),Y_Z(pp),fields_domain%Ro,fields_domain%DR,fields_domain%Zo,fields_domain%DZ
          
          IR = INT(FLOOR((Y_R(pp)  - fields_domain%Ro + &
               0.5_rp*fields_domain%DR)/fields_domain%DR) + 1.0_rp,idef)
          IZ = INT(FLOOR((Y_Z(pp)  + ABS(fields_domain%Zo) + &
               0.5_rp*fields_domain%DZ)/fields_domain%DZ) + 1.0_rp,idef)

!          write(output_unit_write,*) pp

!          write(output_unit_write,'("Size of fields_domain R: ",I16)') &
!               size(fields_domain%FLAG2D,1)
!          write(output_unit_write,'("Size of fields_domain Z: ",I16)') &
!               size(fields_domain%FLAG2D,2)   
          
!          if ((IR.lt.0).or.(IZ.lt.0)) then
!             write(output_unit_write,'("YR:",E17.10)') Y_R(pp)
!             write(output_unit_write,'("YZ:",E17.10)') Y_Z(pp)
!             write(output_unit_write,'("IR: ",I16)') IR
!             write(output_unit_write,'("IZ: ",I16)') IZ
!          end if

#if DBG_CHECK    
          if ((IR.lt.1).or.(IZ.lt.1).or.(IR.GT. &
               bfield_2d%NR).OR.(IZ.GT.bfield_2d%NZ)) then
          
             !write(6,'("YR:",E17.10)') Y_R(pp)
             !write(6,'("YZ:",E17.10)') Y_Z(pp)
             !write(6,'("IR: ",I16)') IR
             !write(6,'("IZ: ",I16)') IZ
             !write(6,*) 'NR',bfield_2d%NR
             !write(6,*) 'NZ',bfield_2d%NZ
             !write(6,*) 'FLAG',fields_domain%FLAG2D(IR,IZ)
             !call KORC_ABORT(23)
          end if
#endif
          
          if ((IR.lt.1).or.(IZ.lt.1).or. &
               (IR.GT.bfield_2d%NR).OR.(IZ.GT.bfield_2d%NZ).or. &
               (fields_domain%FLAG2D(IR,IZ).NE.1_is)) then
             if (F%Analytic_IWL.eq.'NONE') then
                flag(pp) = 0_is
             else if (F%Analytic_IWL.eq.'D3D') then
                if ((IR.lt.floor(bfield_2d%NR/6._rp)).and. &
                     (IZ.gt.floor(bfield_2d%NZ/5._rp)).and. &
                     (IZ.lt.floor(4._rp*bfield_2d%NZ/5._rp))) then

                   Rwall=F%circumradius*cos(C_PI/F%ntiles)/ &
                        (cos((mod(Y_PHI(pp),2*C_PI/F%ntiles))-C_PI/F%ntiles))
                   if (Y_R(pp).lt.Rwall) flag(pp) = 0_is

                else
                   flag(pp) = 0_is
                endif
             endif

!             write(output_unit_write,'("Shit''s fucked.")')
          end if
       end do      
       !$OMP END SIMD
!       write(output_unit_write,'("Shit''s not fucked.")')
    end if
  end subroutine check_if_in_fields_domain_p