anchor.diff:          if (smooth_nsq) then
               call smooth_array(nsquared, winsz, ydim=3)
           endif
  -        ! if (this_image()==3) print*, ims, ims_u, nxu,"lbound(u3d)",lbound(u3d,1),lbound(u3d,2),lbound(u3d,3)
  -        ! if (this_image()==3) print*, jms, jms_v, nyv,"ubound(u3d)",ubound(u3d,1),ubound(u3d,2),ubound(u3d,3)
  +
           !$omp parallel firstprivate(ims, ime, jms, jme, ims_u, ime_u, jms_v, jme_v, kms, kme, nx,nxu,ny,nyv,nz), &
           !$omp firstprivate(reverse, vsmooth, winsz, using_blocked_flow), default(none), &
           !$omp private(i,j,k,step, uk, vi, east, west, north, south, top, bottom, u1d, v1d), &
  
positive.code:                             ! if (reverse) then
                              !     u3d(i,j,k) = u3d(i,j,k) - u_perturbation(i,j,k) * linear_contribution
                              ! else
                                  ! if (this_image()==3) then
                                  !     if ((j==4).and.(i>13).and.(k>16)) then
                                  !         print*, "i,j",i,k, u_perturbation(i,j,k)
                                  !     endif
                                  ! endif
                                  u3d(i+ims_u-1,j,k+jms-1) = u3d(i+ims_u-1,j,k+jms-1) + u_perturbation(i,j,k) *linear_contribution ! * linear_mask(min(nx,i),min(ny,k)) * (1-blocked)
                              ! endif
                          endif
  
negative.code:             this_n = stop_pos-start_pos+1
              ! $omp do
  
              do ijk = start_pos, stop_pos
                  ! loop over the combined ijk space to improve parallelization (more granular parallelization)
                  ! because it is one combined loop, we have to calculate the i,k indicies from the combined ik variable
  