anchor.diff:                              dpos=step
                           endif
                       end do
  +                    
  +                    spos=1
  +                    curspd=calc_speed(domain%u(i,j,uk), domain%v(vi,j,k))
  +                    do step=1,n_spd_values
  +                        if (curspd>spd_values(step)) then
  +                            spos=step
  +                        endif
  +                    end do
  +                    
                       npos=1
  -                    curnsq = calc_stability(domain%th(i,j,k), domain%th(i,top,k),  &
  -                                            domain%pii(i,j,k),domain%pii(i,top,k), &
  -                                            domain%z(i,j,k),  domain%z(i,top,k),   &
  -                                            domain%qv(i,j,k), domain%qv(i,top,k),  &
  -                                            domain%cloud(i,j,k)+domain%ice(i,j,k)+domain%qrain(i,j,k)+domain%qsnow(i,j,k))
  +                    curnsq = calc_stability(domain%th(vi,bottom,uk), domain%th(vi,top,uk),  &
  +                                            domain%pii(vi,bottom,uk),domain%pii(vi,top,uk), &
  +                                            domain%z(vi,bottom,uk),  domain%z(vi,top,uk),   &
  +                                            domain%qv(vi,bottom,uk), domain%qv(vi,top,uk),  &
  +                                            domain%cloud(vi,j,uk)+domain%ice(vi,j,uk)+domain%qrain(vi,j,uk)+domain%qsnow(vi,j,uk))
                       curnsq = log(curnsq)
                       do step=1,n_nsq_values
                           if (curnsq>nsq_values(step)) then
  
positive.code:             write(*,*) i/real(n_dir_values)*100," %"
              ! 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...
              domain%u=calc_u(dir_values(i),1.0)
              domain%v=calc_v(dir_values(i),1.0)
              do j=1,n_nsq_values
                  
                  ! 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.
                  do ii=1,nz
                      u_LUT(i,j,:,ii,:)=(domain%u(:,ii,:)-calc_u(dir_values(i),1.0))
                      v_LUT(i,j,:,ii,:)=(domain%v(:,ii,:)-calc_v(dir_values(i),1.0))
                  end do
              end do
          end do
  
negative.code:                     ! calculate the weights and the "next" u/v position
                      ! "next" usually = pos+1 but for edge cases next = 1 or n
                      dweight=calc_weight(dir_values, dpos,nextd,curdir)
                      sweight=calc_weight(spd_values, spos,nexts,curspd)
                      nweight=calc_weight(nsq_values, npos,nextn,curnsq)
  !                     if ((i==50).and.(k==50)) then
  !                         print*, i,j,k, exp(curnsq)
  !                         print*, domain%th(vi,bottom,uk),domain%p(vi,bottom,uk),domain%z(vi,bottom,uk),domain%qv(vi,bottom,uk)
  !                         print*, u_LUT(spos,dpos,npos,i,j,k), v_LUT(spos,dpos,npos,i,j,k)
  !                     endif
                      if (k<=ny) then
                          ! perform linear interpolation between LUT values
                          wind_first =      nweight  * (dweight * u_LUT(spos,dpos,npos,i,j,k)  + (1-dweight) * u_LUT(spos,nextd,npos,i,j,k))    &
                                      +  (1-nweight) * (dweight * u_LUT(spos,dpos,nextn,i,j,k) + (1-dweight) * u_LUT(spos,nextd,nextn,i,j,k))
                          wind_second=      nweight  * (dweight * u_LUT(nexts,dpos,npos,i,j,k) + (1-dweight) * u_LUT(nexts,nextd,npos,i,j,k))    &
                                      +  (1-nweight) * (dweight * u_LUT(nexts,dpos,nextn,i,j,k)+ (1-dweight) * u_LUT(nexts,nextd,nextn,i,j,k))
                          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)
                      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))    &
                                      +  (1-nweight) * (dweight * v_LUT(spos,dpos,nextn,i,j,k) + (1-dweight) * v_LUT(spos,nextd,nextn,i,j,k))
                          wind_second=      nweight  * (dweight * v_LUT(nexts,dpos,npos,i,j,k) + (1-dweight) * v_LUT(nexts,nextd,npos,i,j,k))    &
                                      +  (1-nweight) * (dweight * v_LUT(nexts,dpos,nextn,i,j,k)+ (1-dweight) * v_LUT(nexts,nextd,nextn,i,j,k))
                          v_perturbation(i,j,k)=v_perturbation(i,j,k) * (1-linear_update_fraction) &
                                      + linear_update_fraction * (sweight*wind_first + (1-sweight)*wind_second)
                                      
                          domain%v(i,j,k) = domain%v(i,j,k) + v_perturbation(i,j,k)
                      endif
                      
  