kernel smoothing in openmp

This code takes a list of positions x,y,z and creates a kernal smoothed mesh using smoothing lengths hsml.
Basically, I want to just to loop through i as fast a possible as every calculation is independent of one another. It just has to populate mesh_out which is shared. kern_result is a function representing the smoothing kernel.

!$OMP PARALLEL DO DEFAULT(private) SHARED(mesh_out)
do i = 1,Ntot
    xtmp = x0(i)
    ytmp = y0(i)
    ztmp = z0(i)
    hsmltmp = hsml(i)

    if (h.lt.hmin) then
        hsmltmp = hmin
    elseif (h.gt.hmax) then
        hsmltmp = hmax
    end if

    Xrel = xtmp-Xc+Lx/2;
    Yrel = ytmp-Yc+Ly/2;
    Zrel = ztmp-Zc+Lz/2;

    if (Xrel.gt.-hsmltmp.and.Yrel.gt.-hsmltmp.and.Xrel-Lx.lt.hsmltmp &
      .and.Yrel-Ly.lt.hsmltmp.and.Zrel.gt.-hsmltmp.and.Zrel-Lz.lt.hsmltmp) then
        binned_particles = binned_particles + 1

        pixelX = (FLOOR(Xrel/pixelsizeX)+.5)*pixelsizeX
        pixelY = (FLOOR(Yrel/pixelsizeY)+.5)*pixelsizeY
        pixelZ = (FLOOR(Zrel/pixelsizeZ)+.5)*pixelsizeZ

        !write(*,*)"pixelXYZ",pixelX,pixelY,pixelZ

        nx = INT(CEILING(hsmltmp/pixelsizeX))
        ny = INT(CEILING(hsmltmp/pixelsizeY))
        nz = INT(CEILING(hsmltmp/pixelsizeZ))
        !write(*,*)"Nxyz",nx,ny,nz
        normalization = 0

        do dx = -nx,nx
            do dy = -ny,ny
                do dz = -nz,nz
                    distX = Xrel-(pixelX-REAL(dx)*pixelsizeX)
                    distY = Yrel-(pixelY-REAL(dy)*pixelsizeY)
                    distZ = Zrel-(pixelZ-REAL(dz)*pixelsizeZ)
                    r2 = distX**2+distY**2+distZ**2;
                    if (r2.lt.hsmltmp**2) then
                        kern_result = kernel(hsmltmp,r2)
                        normalization = normalization + kern_result
                        xloc = CEILING(pixelX/pixelsizeX-REAL(dx));
                        yloc = CEILING(pixelY/pixelsizeY-REAL(dy));
                        zloc = CEILING(pixelZ/pixelsizeZ-REAL(dz));
                        !write(*,*)"XYZloc",xloc,yloc,zloc
                        if (xloc.gt.0.and.xloc.le.num_pixelsX.and.yloc.gt.0.and.yloc.le.num_pixelsY.and. &
                             zloc.gt.0.and.zloc.le.num_pixelsZ) then
                            kern_result = kernel(hsmltmp,r2)
                            !$OMP ATOMIC UPDATE
                            mesh_out(INT(zloc),INT(yloc),INT(xloc)) = &
                                    mesh_out(INT(zloc),INT(yloc),INT(xloc)) &
                                                                     + kern_result/normalization

                            !write(*,*)"kern",hsmltmp,r2,kern_result
                        end if
                    end if
                end do
            end do
        end do
    end if
end do
!$OMP END PARALLEL DO


RECURSIVE FUNCTION kernel(hdum,rdum) RESULT(kern)
real*8 rdum
real*8 hdum
real*8 kern

! Kernal Coefficients
real*8, parameter :: K1 = 2.546479089470
real*8, parameter :: K2 = 15.278874536822
real*8, parameter :: K3 = 5.092958178941
real*8 u

u = sqrt(rdum)/hdum;

if (u.lt.0.5) then
    kern = K1+K2*(u-1)*u**2;
else
    kern = K3*(1-u)**3;
end if
END FUNCTION kernel