tmpar <- par( mar = c(3.1,3.1,0,0), mfrow = c(1,2) )
on.exit( {par( tmpar ); rm( data, ypos, ylim, xlab, ylab, lsm, tmpar )} )

xlab <- c("(b) control","(a) tree")
ylab <- c("","score means")
ypos <- c(2,2)

lsm <- lsmean( trees.all, factors = c("soil","date") )
lsm$type <- factor( !is.na( match( lsm$soil,
   c("dtac2","nc11004","ne-332","nm-6") ) ) )
ylim <- range( lsm$pred )

lsm <- split( lsm, lsm$type )
# following values are SE for means taken for SAS output
lsm[[1]]$se <- rep( .84086871, nrow( lsm[[1]] ) )
lsm[[2]]$se <- rep( .63433993, nrow( lsm[[2]] ) )

data <- split( trees, trees$type )

for ( i in 2:1 ) { # "FALSE" comes before "TRUE"
   int.plot( trees.all, data[[i]], c("date","soil"), lsm[[i]],
      bar.plot = "none",
      xlim = c(1,7.7), ylim = ylim, xaxt = "n", xlab = "", ylab = "" )
   se.bar( 4.5, ypos[i], lsm[[i]]$se[1] )
   axis(1,1:6)
   mtext( paste( xlab[i], "soils by date" ), 1, 2 )
   mtext( ylab[i], 2, 2 )
   tmp <- lsm[[i]][lsm[[i]]$date==6,]
   text( rep( 6.2, nrow( tmp ) ), uncollide( tmp$pred, .1 ),
      as.character( tmp$soil ), adj = 0 )
}