anchor.diff:              do j=1,nz
                   do i=1,nx
   
  -                    ! look up vsmooth gridcells up to nz at the maximum
  -                    top = min(j+vsmooth, nz)
  -                    ! if (top-j)/=vsmooth, then look down enough layers to make the window vsmooth in size
  -                    bottom = max(1, j - (vsmooth - (top-j)))
  -
  -                    if (.not.reverse) then
  -                        nsquared(i+ims-1,j+kms-1,k+jms-1) = calc_stability(domain%potential_temperature%data_3d(i+ims-1,bottom+kms-1,k+jms-1),                              &
  -                                                         domain%potential_temperature%data_3d(i+ims-1,top+kms-1,k+jms-1),                                                   &
  -                                                         domain%exner%data_3d(i+ims-1,bottom+kms-1,k+jms-1),domain%exner%data_3d(i+ims-1,top+kms-1,k+jms-1),                &
  -                                                         domain%z%data_3d(i+ims-1,bottom+kms-1,k+jms-1),  domain%z%data_3d(i+ims-1,top+kms-1,k+jms-1),                      &
  -                                                         domain%water_vapor%data_3d(i+ims-1,bottom+kms-1,k+jms-1), domain%water_vapor%data_3d(i+ims-1,top+kms-1,k+jms-1),   &
  -                                                         domain%cloud_water_mass%data_3d(i+ims-1,j+kms-1,k+jms-1)     &
  -                                                         +domain%cloud_ice_mass%data_3d(i+ims-1,j+kms-1,k+jms-1)      &
  -                                                         +domain%rain_mass%data_3d(i+ims-1,j+kms-1,k+jms-1)           &
  -                                                         +domain%snow_mass%data_3d(i+ims-1,j+kms-1,k+jms-1))
  -
  -                        nsquared(i+ims-1,j+kms-1,k+jms-1) = max(min_stability, min(max_stability, &
  -                                                nsquared(i+ims-1,j+kms-1,k+jms-1) * nsq_calibration(i,k)))
  +                    if (variable_N) then
  +                        ! look up vsmooth gridcells up to nz at the maximum
  +                        top = min(j+vsmooth, nz)
  +                        ! if (top-j)/=vsmooth, then look down enough layers to make the window vsmooth in size
  +                        bottom = max(1, j - (vsmooth - (top-j)))
  +
  +                        hydrometeors = 0
  +                        if (associated(domain%cloud_water_mass%data_3d)) &
  +                            hydrometeors = hydrometeors + domain%cloud_water_mass%data_3d(i+ims-1,j+kms-1,k+jms-1)
  +                        if (associated(domain%cloud_ice_mass%data_3d)) &
  +                            hydrometeors = hydrometeors + domain%cloud_ice_mass%data_3d(i+ims-1,j+kms-1,k+jms-1)
  +                        if (associated(domain%rain_mass%data_3d)) &
  +                            hydrometeors = hydrometeors + domain%rain_mass%data_3d(i+ims-1,j+kms-1,k+jms-1)
  +                        if (associated(domain%snow_mass%data_3d)) &
  +                            hydrometeors = hydrometeors + domain%snow_mass%data_3d(i+ims-1,j+kms-1,k+jms-1)
  +
  +
  +                        if (.not.reverse) then
  +                            nsquared(i+ims-1,j+kms-1,k+jms-1) = calc_stability(domain%potential_temperature%data_3d(i+ims-1,bottom+kms-1,k+jms-1),                              &
  +                                                             domain%potential_temperature%data_3d(i+ims-1,top+kms-1,k+jms-1),                                                   &
  +                                                             domain%exner%data_3d(i+ims-1,bottom+kms-1,k+jms-1),domain%exner%data_3d(i+ims-1,top+kms-1,k+jms-1),                &
  +                                                             domain%z%data_3d(i+ims-1,bottom+kms-1,k+jms-1),  domain%z%data_3d(i+ims-1,top+kms-1,k+jms-1),                      &
  +                                                             domain%water_vapor%data_3d(i+ims-1,bottom+kms-1,k+jms-1), domain%water_vapor%data_3d(i+ims-1,top+kms-1,k+jms-1),   &
  +                                                             hydrometeors)
  +
  +                            nsquared(i+ims-1,j+kms-1,k+jms-1) = max(min_stability, min(max_stability, &
  +                                                    nsquared(i+ims-1,j+kms-1,k+jms-1)))! * nsq_calibration(i,k)))
  +                        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.
  +                            nsquared(i+ims-1,j+kms-1,k+jms-1) = 3e-6
  +                        endif
                       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.
  -                        nsquared(i+ims-1,j+kms-1,k+jms-1) = 3e-6
  +                        nsquared(i+ims-1,j+kms-1,k+jms-1) = N_squared
                       endif
                   end do
                   ! look up table is computed in log space
  
positive.code:             return
          endif
  
          start_z = minval(z_bottom)
          end_z = maxval(z_top)
  
          allocate(internal_z_top(size(lt_data%u_perturb,1),size(lt_data%u_perturb,2)))
          internal_z_top = maxval(z_top)
          internal_z_top(buffer:buffer+size(z_top,1)-1, buffer:buffer+size(z_top,2)-1) = z_top(:,:)
  
          allocate(internal_z_bottom(size(lt_data%u_perturb,1),size(lt_data%u_perturb,2)), source=maxval(z_bottom))
          internal_z_bottom(buffer:buffer+size(z_bottom,1)-1, buffer:buffer+size(z_bottom,2)-1) = z_bottom
  
          allocate(layer_count(size(lt_data%u_perturb,1), size(lt_data%u_perturb,2)))
  
negative.code:             ! If we are using density in the advection calculations, modify the linear perturbation
              ! to get the vertical velocities closer to what they would be without density (boussinesq)
              ! need to check if this makes the most sense when close to the surface
              ! if (useDensity) then
                  ! linear_mask = linear_mask * 2
              ! endif
          endif
  
          ! set up nsq_calibration variable
  