anchor.diff:                              ! 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
  
positive.code:         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), &
  
negative.code:                         sweight = calc_weight(spd_values, spos, nexts, curspd)
                          nweight = calc_weight(nsq_values, npos, nextn, curnsq)
  
                          ! perform linear interpolation between LUT values
                          if (k<=ny) then
                              wind_first =      nweight  * (dweight * hi_u_LUT(spos, dpos,npos, i,j,k) + (1-dweight) * hi_u_LUT(spos, nextd,npos, i,j,k))   &
  