anchor.diff:          write(*,*) "Max U perturb:",maxval(u_perturbation), "Max V perturb:",maxval(v_perturbation)
           if (use_spatial_linear_fields) then
               if ((.not.allocated(hi_u_LUT) .and. (.not.reverse)) .or. ((.not.allocated(rev_u_LUT)) .and. reverse)) then
  -
  +                ! Need linear_contribution to be 1 for initialize_spatial_winds, it will be applied at "real" time. 
  +                ! should probably set linear_mask to 1 too... for not linear_mask is set at LUT creation. 
  +                linear_mask = linear_mask / linear_contribution
  +                saved_linear_contribution = linear_contribution
  +                linear_contribution = 1.0
  +                
                   write(*,*) "Generating a spatially variable linear perturbation look up table"
                   call initialize_spatial_winds(domain,options,reverse,useDensity)
  +                
  +                linear_contribution = saved_linear_contribution
  +                linear_mask = linear_mask * linear_contribution
               ! else
               !     write(*,*) "Skipping spatial wind field for presumed domain repeat"
  +            endif
  +            if (minval(linear_mask) > linear_contribution) then
  +                write(*,*) " "
  +                write(*,*) "  ---------------  !!!WARNING!!!  --------------- "
  +                write(*,*) " Spatially variable linear mask is fixed at LUT creation"
  +                write(*,*) " If the LUT was read from disk it could be 'wrong' "
  +                write(*,*) "  ---------------  !!!WARNING!!!  --------------- "
  +                write(*,*) " "
               endif
           endif
   
  
positive.code: 
      end subroutine initialize_spatial_winds
  
  ! Code to test that the direction calculation works.
  !     if ((k==n_spd_values).and.(j==1)) then
  !         print*, "Direction:", &
  !             dir_values(i), "Calc=", calc_direction(domain%u(1,1,1), domain%v(1,1,1)), &
  !             "  U=",domain%u(1,1,1), "  V=",domain%v(1,1,1)
  !         if (abs(dir_values(i) - calc_direction(domain%u(1,1,1), domain%v(1,1,1)))>1e-5) then
  !             print*, "ERROR ERROR ERROR ERROR ERROR ERROR ERROR ERROR ERROR ERROR ERROR "
  !             print*, "Error with direction:", &
  !                 dir_values(i), "Calc=", calc_direction(domain%u(1,1,1), domain%v(1,1,1)), &
  !                 "  U=",domain%u(1,1,1), "  V=",domain%v(1,1,1)
  !         endif
  !     endif
  
  
      function calc_weight(indata, bestpos, nextpos, match) result(weight)
          ! simply calculate the weights between the positions bestpos and nextpos
  
negative.code:                         u_perturbation(i,j,k)=u_perturbation(i,j,k) * (1-linear_update_fraction) &
                                      + linear_update_fraction * (sweight*wind_first + (1-sweight)*wind_second)
  
                          domain%u(i,j,k) = domain%u(i,j,k) + u_perturbation(i,j,k) * linear_contribution
                      endif
                      if (i<=nx) then
                          wind_first =      nweight  * (dweight * v_LUT(spos,dpos,npos,i,j,k)  + (1-dweight) * v_LUT(spos,nextd,npos,i,j,k))    &
  