anchor.diff:              ! 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
  +            ! if (useDensity) then
  +                ! linear_mask = linear_mask * 2
  +            ! endif
           endif
   
           ! set up nsq_calibration variable
  
positive.code:                                 ( lt_data_m%v_perturb(1+buffer:fftnx-buffer,     buffer:fftny-buffer)         &
                                  + lt_data_m%v_perturb(1+buffer:fftnx-buffer,     1+buffer:fftny-buffer+1)) )) / 2
  
                          call copy_data_to_remote(temporary_u, u_grids, hi_u_LUT, i,j,k, z)
                          call copy_data_to_remote(temporary_v, v_grids, hi_v_LUT, i,j,k, z)
  
  
negative.code:                 end do
              endif
          end do
          ! nsquared(:,:,:) = log(6e-4)
          ! nsquared(:,:,:) = log(1e-7)
          ! nsquared(:,:,:) = log(7.73310239E-06)
          ! !$omp end do
          ! !$omp end parallel
  
  
          ! smooth array has it's own parallelization, so this probably can't go in a critical section
          if (smooth_nsq) then
              call smooth_array(nsquared, winsz, ydim=3)
          endif
          ! if (this_image()==3) print*, ims, ims_u, nxu,"lbound(u3d)",lbound(u3d,1),lbound(u3d,2),lbound(u3d,3)
          ! if (this_image()==3) print*, jms, jms_v, nyv,"ubound(u3d)",ubound(u3d,1),ubound(u3d,2),ubound(u3d,3)
          !$omp parallel firstprivate(ims, ime, jms, jme, ims_u, ime_u, jms_v, jme_v, kms, kme, nx,nxu,ny,nyv,nz), &
          !$omp firstprivate(reverse, vsmooth, winsz, using_blocked_flow), default(none), &
          !$omp private(i,j,k,step, uk, vi, east, west, north, south, top, bottom, u1d, v1d), &
  