anchor.diff:      !! compute FFT(terrain), and dzdx,dzdy components
       !!
       !!----------------------------------------------------------
  -    subroutine setup_linwinds(domain,options,reverse,useDensity)
  +    subroutine setup_linwinds(domain, options, reverse, useDensity)
           implicit none
           class(linearizable_type),intent(inout)  :: domain
           type(options_type),      intent(in)     :: options
  
positive.code:             endif
              u_LUT=>hi_u_LUT
              v_LUT=>hi_v_LUT
          endif
  
          if (options%debug) then
  
negative.code:                                                                 +domain%qrain(i,j,k)+domain%qsnow(i,j,k))
                          domain%nsquared(i,j,k) = max(min(domain%nsquared(i,j,k) * nsq_calibration(i,k), max_stability), min_stability)
                      else
                          ! Low-res boundary condition variables will be in a different array format.  It should be 
                          ! easy enough to call calc_stability after e.g. transposing z and y dimension, but some
                          ! e.g. pii will not be set in the forcing data, so this may need a little thought. 
                          domain%nsquared(i,j,k) = 3e-6
                      endif
                  end do
  