anchor.diff:          debug=.False.
           
       end subroutine linear_perturb
  -end module linear_theory_winds+end module linear_theory_winds
  
positive.code:         endif
          if (reverse) then
              write(*,*) "ERROR: reversing linear winds not set up for parallel fftw computation yet"
              stop
          endif
  
  
negative.code:             ! set the domain wide U and V values to the current u and v values
              ! this could use u/v_perturbation, but those would need to be put in a linearizable structure...
              do k=1,n_spd_values
                  do j=1,n_nsq_values
                      domain%u=calc_u(dir_values(i),spd_values(k))
                      domain%v=calc_v(dir_values(i),spd_values(k))
                      ! calculate the linear wind field for the current u and v values
                      call linear_winds(domain,exp(nsq_values(j)), 0, reverse,useDensity,debug=debug)
                      debug=.False.
                      u_LUT(k,i,j,:,:,:)=(domain%u-calc_u(dir_values(i),spd_values(k)))
                      v_LUT(k,i,j,:,:,:)=(domain%v-calc_v(dir_values(i),spd_values(k)))
                      print*, i,j,k
                      print*, "dir=",dir_values(i),"  nsq=",exp(nsq_values(j)),"  spd=",spd_values(k)
                      print*, "u=",calc_u(dir_values(i),spd_values(k)), "  v=",calc_v(dir_values(i),spd_values(k))
                      print*, "ulut=",u_LUT(k,i,j,100,1,100), "  vlut=",v_LUT(k,i,j,100,1,100)
                      call io_write3d("u_LUT_dir"//trim(str(i))//"_spd"//trim(str(k))//"_nsq"//trim(str(j))//".nc","data", u_LUT(k,i,j,:,:,:) )
                      call io_write3d("v_LUT_dir"//trim(str(i))//"_spd"//trim(str(k))//"_nsq"//trim(str(j))//".nc","data", v_LUT(k,i,j,:,:,:) )
                      print*, "---------------------------------------------"
                  end do
              end do
          end do
  