subroutine hpf_fill_in_data(target, w, h, ccr, cci, cstep, nmin, nmax)
       integer, intent(in)     :: w, h
       byte, intent(out)       :: target(w,h)
       real*8, intent(in)      :: ccr, cci, cstep
       integer, intent(in)     :: nmin, nmax
 !hpf$ distribute target(*, cyclic)

       integer :: cx, cy
       cx = w/2
       cy = h/2

       forall(ix = 1:w, iy = 1:h)                                  &
                  target(ix,iy) = mandel_val(CMPLX(ccr + ((ix-cx)*cstep),   &
                                          cci + ((iy-cx)*cstep),   &
                                          KIND=KIND(0.0D0)),       &
                                    nmin, nmax)                                                                             
 contains

     pure byte function mandel_val(x, nmin, nmax)
       complex(KIND=KIND(0.0D0)), intent(in) :: x
       integer, intent(in)                   :: nmin, nmax

       integer                         :: n

       real(kind=KIND(0.0D0)) :: xorgr, xorgi, xr, xi, xr2, xi2, rad2
       logical                :: keepgoing

       n = -1
       xorgr = REAL(x)
       xorgi = AIMAG(x)
       xr = xorgr
       xi = xorgi

       do
          n = n + 1
          xr2 = xr*xr
          xi2 = xi*xi
          xi = 2*(xr*xi) + xorgi
          keepgoing = n < nmax
          rad2 = xr2 + xi2
          xr = xr2 - xi2 + xorgr
          if (keepgoing .AND. (rad2 <= 4.0)) cycle
          exit
       end do

       if (n >= nmax) then
          mandel_val = nmax-nmin
       else
          mandel_val = MOD(n, nmax-nmin)
       end if

     end function mandel_val

 end subroutine hpf_fill_in_data