Your modification on the code is correct.

PS: A more elegant way is adding your function as a new real space function in functions.f90, and link it into user-defined function ("function userfunc(x,y,z)"); in your function code simply set the returned value to zero if it is found to be >=0.5. In this case you can integrate your function in basin analysis module in terms of integrating user-defined function, any modification of basin.f90 is not needed.

I'm faced with the task of integrating existing Bader basins, but in a specific way. The integration needs to be done in a mixed-type grid. For my task, I need to take only those voxels in which the value of my function (below, MYFUNC) is >= 0.5. Am I correct in understanding that I need to add code in two places in the basin.f90 file, namely:
1)
...
!Use DFT integration algorithm to integrate the region inside the trust radius
...
!$OMP do schedule(DYNAMIC)
do ipt=1,nintgrid
ptx=gridatt(ipt)%x
pty=gridatt(ipt)%y
ptz=gridatt(ipt)%z
MYFUNC = f(ptx, pty, ptz, "f") #my function
if (MYFUNC < 0.5D0) cycle
rx=ptx-CPpos(1,numcp) !The relative distance between the current point and the corresponding attractor
ry=pty-CPpos(2,numcp)
rz=ptz-CPpos(3,numcp)
!Calculate switching function
dist=dsqrt(rx*rx+ry*ry+rz*rz)
....
2)
...
!--------- Integrating uniform grids, basin boundary grids will be calculated in the later stage
!$OMP do schedule(DYNAMIC)
do iz=izlow,izup
do iy=iylow,iyup
do ix=ixlow,ixup
if ((itype==2.or.itype==3).and.interbasgrid(ix,iy,iz)) cycle !If refine boundary grid at next stage, we don't calculate them at present stage
call getgridxyz(ix,iy,iz,rnowx,rnowy,rnowz)
MYFUNC = f(rnowx, rnowy, rnowz, "f")
if (MYFUNC < 0.5D0) cycle
iatt=gridbas(ix,iy,iz)
! do icp=1,numcp
! dist=dsqrt( (rnowx-CPpos(1,icp))**2+(rnowy-CPpos(2,icp))**2+(rnowz-CPpos(3,icp))**2 )
! if (disttest<trustrad(icp)) cycle cycix !The function inside trust radius is integrated by DFT integration
! end do
rx=rnowx-CPpos(1,iatt) !The relative distance between the current point and the exact position of the corresponding attractor
ry=rnowy-CPpos(2,iatt)
...

I would be very grateful if you could tell me if I'm on the right track and if I need to add an "if MYFUNC..." condition somewhere else in the code (When integrating the basin boundary mb?)