anchor.diff:  
       end subroutine spatial_winds
   
  +    !>----------------------------------------------------------
  +    !! Copy options from the options data structure into module level variables
  +    !!
  +    !!----------------------------------------------------------
       subroutine set_module_options(options)
           implicit none
           type(options_type), intent(in) :: options
   
  -        original_buffer = options%lt_options%buffer
  -        variable_N = options%lt_options%variable_N
  -        smooth_nsq = options%lt_options%smooth_nsq
  +        original_buffer       = options%lt_options%buffer
  +        variable_N            = options%lt_options%variable_N
  +        smooth_nsq            = options%lt_options%smooth_nsq
   
           stability_window_size = options%lt_options%stability_window_size        ! window to average nsq over
  -        max_stability = options%lt_options%max_stability                        ! limits on the calculated Brunt Vaisala Frequency
  -        min_stability = options%lt_options%min_stability                        ! these may need to be a little narrower.
  -        linear_contribution = options%lt_options%linear_contribution            ! multiplier on uhat,vhat before adding to u,v
  +        max_stability         = options%lt_options%max_stability                ! limits on the calculated Brunt Vaisala Frequency
  +        min_stability         = options%lt_options%min_stability                ! these may need to be a little narrower.
  +        linear_contribution   = options%lt_options%linear_contribution          ! multiplier on uhat,vhat before adding to u,v
                                                                                   ! =fractional contribution of linear perturbation to wind field
           ! these are defined per call to permit different values for low res and high res domains
  -        ! N_squared = options%lt_options%N_squared                              ! static Brunt Vaisala Frequency (N^2) to use
  -        ! rm_N_squared = options%lt_options%rm_N_squared                        ! static BV Frequency (N^2) to use in removing linear wind field
  +        ! N_squared              = options%lt_options%N_squared                 ! static Brunt Vaisala Frequency (N^2) to use
  +        ! rm_N_squared           = options%lt_options%rm_N_squared              ! static BV Frequency (N^2) to use in removing linear wind field
           ! rm_linear_contribution = options%lt_options%rm_linear_contribution    ! fractional contribution of linear perturbation to remove wind field
  -        ! remove_lowres_linear = options%lt_options%remove_lowres_linear        ! remove the linear mountain wave from low res forcing model
  +        ! remove_lowres_linear   = options%lt_options%remove_lowres_linear      ! remove the linear mountain wave from low res forcing model
   
           linear_update_fraction    = options%lt_options%linear_update_fraction   ! controls the rate at which the linearfield updates
                                                                                   ! =fraction of linear perturbation to add each time step
  
positive.code:         !$omp end do
          !$omp barrier
          !$omp do
          do k=1,ny+1
              do j=1,nz
                  do i=1,nx+1
                      uk=min(k,ny)
                      vi=min(i,nx)
  
                      dpos=1
                      curdir=calc_direction(domain%u(i,j,uk), domain%v(vi,j,k))
                      do step=1,n_dir_values
                          if (curdir>dir_values(step)) then
                              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
  
                      west  = max(i-winsz,1)
                      east  = min(i+winsz,nx)
                      bottom= max(j-winsz,1)
                      top   = min(j+winsz,nz)
                      south = max(k-winsz,1)
                      north = min(k+winsz,ny)
                      n = (((east-west)+1) * ((top-bottom)+1) * ((north-south)+1))
  
                      curnsq = sum(log(domain%nsquared(west:east, bottom:top,south:north)))
                      curnsq = curnsq / n
  
                      npos=1
                      do step=1,n_nsq_values
                          if (curnsq > nsq_values(step)) then
                              npos = step
                          endif
                      end do
  
                      top = min(j+1, nz)
                      if (top == j) then
                          bottom = j-1
                      else
                          bottom = j
                      endif
  
                      ! 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)
                      ! perform linear interpolation between LUT values
                      if (k<=ny) then
                          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)
  
negative.code: 
      end subroutine set_module_options
  
      !>----------------------------------------------------------
      !! Called from linear_perturb the first time perturb is called
      !! compute FFT(terrain), and dzdx,dzdy components
      !!
      !!----------------------------------------------------------
      subroutine setup_linwinds(domain,options,reverse,useDensity)
          implicit none
          class(linearizable_type),intent(inout)::domain
  