anchor.diff:          h2_v =  h_v  - h1_v
   
           ! In case one wants to see how the terrain is split by smoothing, activate the block below and run in debug:
  -        if ((this_image()==1).and.(options%parameters%debug)) then
  -          call io_write("terrain_smooth_h1.nc", "h1", h1(:,:) )
  -          call io_write("terrain_smooth_h2.nc", "h2", h2(:,:) )
  +        ! if ((this_image()==1).and.(options%parameters%debug)) then
  +        !   call io_write("terrain_smooth_h1.nc", "h1", h1(:,:) )
  +        !   call io_write("terrain_smooth_h2.nc", "h2", h2(:,:) )
           !   call io_write("h1_u.nc", "h1_u", h1_u(:,:) )
           !   call io_write("h2_u.nc", "h2_u", h2_u(:,:) )
  -        endif
  +        ! endif
   
           if (this_image()==1) then
              print*, "       Max of full topography", MAXVAL(global_terrain )
  
positive.code:                     zr_u(:,i,:)  = (z_u(:,i,:) - z_u(:,i-1,:)) / (dz_scl(i)/2 + dz_scl(i-1)/2 )  ! if dz_scl(i-1) = 0 (and no error)  k=1 can be included
                      zr_v(:,i,:)  = (z_v(:,i,:) - z_v(:,i-1,:)) / (dz_scl(i)/2 + dz_scl(i-1)/2 )
  
                      ! TODO does this need to go the parallelized grid (i.e ims:ime like line 1060??)
  
                  else ! above the flat_z_height
  
                      zr_u(:,i,:) = 1
  
negative.code:                           - MAXVAL(h2)/s2 * COSH(smooth_height/s2)/SINH(smooth_height/s2)
  
              ! with the new (leuenberger et al 2010) Sleve formulation, the inveribiltiy criterion is as follows:
              ! ( Although an argument could be made to calculate this on the offset (u/v) grid b/c that is most
              !   relevant for advection? In reality this is probably a sufficient approximation, as long as we
              !   aren't pushing the gamma factor too close to zero )
              allocate(gamma_n(this%kds : this%kde+1))
              i=kms                          
              gamma_n(i) =  1                                                     &    
  