anchor.diff:          class(domain_t),  intent(inout) :: this
           type(options_t), intent(in)     :: options
   
  -        real, allocatable :: forcing_terrain_u(:,:), forcing_terrain_v(:,:), delta_terrain(:,:), delta_dzdx_lc(:,:,:)
  -        real, allocatable :: zf_interface(:,:,:), dzf_interface(:,:,:), zf(:,:,:), dzf_mass(:,:,:), dzfdx(:,:,:), dzfdy(:,:,:), delta_dzdx(:,:,:)
  +        real, allocatable :: forcing_terrain_u(:,:), forcing_terrain_v(:,:), delta_terrain(:,:), delta_dzdx_sc(:,:,:), delta_dzdy_sc(:,:,:)
  +        real, allocatable :: zf_interface(:,:,:), dzf_interface(:,:,:), zf(:,:,:), dzf_mass(:,:,:), dzfdx(:,:,:), dzfdy(:,:,:)!, delta_dzdx(:,:,:)
           real :: s
           integer :: i
   
  
positive.code:         class(domain_t), intent(inout)  :: this
          type(options_t), intent(in)     :: options
  
          real, allocatable :: temp(:,:,:), terrain_u(:,:), terrain_v(:,:)
          integer :: i
          real :: smooth_height, H, s, n
          logical :: SLEVE  
  
negative.code: 
  
                  ! write and compare 
          ! if ((this_image()==1).and.(options%parameters%debug)) then  ! Print some diagnostics. USeful for development.         
            if ((this_image()==1)) then  ! Print some diagnostics. USeful for development.         
            call io_write("forcing_terrain.nc", "forcing_terrain", forcing_terrain(:,:) ) ! check in plot
            call io_write("terrain.nc", "terrain", terrain(:,:) ) ! check in plot        
            call io_write("delta_dzdx_sc.nc", "delta_dzdx_sc", delta_dzdx_sc(:,:,:) )
            ! call io_write("delta_dzdx_lc.nc", "delta_dzdx_lc", delta_dzdx(:,:,:) )
            call io_write("dzdx.nc", "dzdx", dzdx(:,:,:) )
            call io_write("zfr_u.nc", "zfr_u", zfr_u(:,:,:) ) ! check in plot
          endif
  