anchor.diff:          class(domain_t), intent(inout)  :: this
           type(options_t), intent(in)     :: options
   
  -        real, allocatable :: temp(:,:,:)
  +        real, allocatable :: temp(:,:,:), terrain_u(:,:), terrain_v(:,:)
           integer :: i, max_level
  -        real :: smooth_height, z_sleve(:,:,:), z_s_u(:,:,:), z_s_v(:,:,:) ! here???
  +        real :: smooth_height, H, s
           logical :: SLEVE
   
   
  
positive.code:                   zr_u                  => this%zr_u,                           &
                    zr_v                  => this%zr_v)
  
            SLEVE = .True.
            
            if (SLEVE==.True.) then  
              ! ______ SLEVE simple Implementation  _______________________
  
              smooth_height = sum(dz(1:nz)) ! this overrides any flat_z_height, will parameterize later
  
              i=kms
  
              ! First calculate 2 adjacent SLEVE z levels:
              H = smooth_height  ! lexicon
              s = smooth_height / 3  ! should become namlist option
              h = terrain(:,:) 
  
              ! This = z_interface, i think..... 
              ! z_sleve(:,i,:) = dz(i) + h * (SINH( (H - dz(i))/s ) / np.sinh(H/s) )
              z_sleve(:,i,:) = terrain
              z_sleve(:,i+1,:) = dz(i+1) + h * (np.sinh( (H - dz(i+1))/s ) / np.sinh(H/s) )
              
              dz_interface(:,i,:)  =  z_sleve(:,i+1,:) - z_sleve(:,i,:)  ! = 'dz_sleve'
              dz_interface(:,i,:)  =  z_sleve(:,i+1,:) - terrain(:,:)  ! = same as above
              z_level_ratio(:,i,:) = dz_interface(:,i,:)  / dz(i)
              dz_mass(:,i,:) = dz_interface(:,i,:) / 2
  
  
  
              zr_u(:,i,:) = 
              
              ! similar for zr_u and zr_v ? But now h is zr_u ?
              ! z_s_u(:,i,:) = 
  
              ! for the u and v grids, z(1) was already initialized with terrain.
              ! but the first level needs to be offset, and the rest of the levels need to be created
              z_u(:,i,:)          = z_u(:,i,:) + dz(i) / 2 * zr_u(:,i,:)
              z_v(:,i,:)          = z_v(:,i,:) + dz(i) / 2 * zr_v(:,i,:)
  
  
  
              ! ____ end SLEVE simple Implementation  _______
              ! ############################################
  
              ! dz_mass(:,i,:)      = dz(i) / 2 * z_level_ratio(:,i,:)
              ! dz_interface(:,i,:) = dz(i) * z_level_ratio(:,i,:)
              ! z(:,i,:)            = terrain + dz_mass(:,i,:)
              ! z_interface(:,i,:)  = terrain
  
  
  
  
  
  
  
            else  
              i = this%grid%kms
  
negative.code:             enddo
  
  
              ! ! # # # # # #   BK debug wind field:    # # # # # # # # #
              ! if (options%parameters%debug ) then
              !     print*, ""
              !     print*, "z_level_ratio(150,i,150) : ", z_level_ratio(150,:,150)
              !     print*, "dzdx at  i,j,k = 150,150,k: " !, k, " : "
              !     print*, dzdx(150,:,150)
              !     print*, "dzdy at i,j,k = 150,150,k: "!, k, " : "
              !     print*, dzdy(150,:,150)
  
              ! endif  
          endif
  
  
          end associate
  