diff --git a/Makefile.in b/Makefile.in
index df8c35c..9485658 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -40,7 +40,7 @@ HASH_OBJ = @HASH_OBJ@
 
 OBJS = array.o cmdparse.o conf.o local.o logging.o main.o memory.o mkdirpp.o \
        reference.o regress.o rtc.o sched.o sources.o sourcestats.o stubs.o \
-       sys.o tempcomp.o util.o $(HASH_OBJ)
+       sys.o smooth.o tempcomp.o util.o $(HASH_OBJ)
 
 EXTRA_OBJS=@EXTRA_OBJECTS@
 
diff --git a/conf.c b/conf.c
index a07a80e..4af99a5 100644
--- a/conf.c
+++ b/conf.c
@@ -73,6 +73,7 @@ static void parse_peer(char *);
 static void parse_pool(char *);
 static void parse_refclock(char *);
 static void parse_server(char *);
+static void parse_smoothtime(char *);
 static void parse_tempcomp(char *);
 
 /* ================================================== */
@@ -185,6 +186,10 @@ static IPAddr bind_cmd_address4, bind_cmd_address6;
  * chronyds being started. */
 static char *pidfile;
 
+/* Smoothing constants */
+static double smooth_max_freq = 0.0; /* in ppm */
+static double smooth_max_wander = 0.0; /* in ppm/s */
+
 /* Temperature sensor, update interval and compensation coefficients */
 static char *tempcomp_sensor_file = NULL;
 static char *tempcomp_point_file = NULL;
@@ -512,6 +517,8 @@ CNF_ParseLine(const char *filename, int number, char *line)
     parse_int(p, &sched_priority);
   } else if (!strcasecmp(command, "server")) {
     parse_server(p);
+  } else if (!strcasecmp(command, "smoothtime")) {
+    parse_smoothtime(p);
   } else if (!strcasecmp(command, "stratumweight")) {
     parse_double(p, &stratum_weight);
   } else if (!strcasecmp(command, "tempcomp")) {
@@ -1164,6 +1171,17 @@ parse_broadcast(char *line)
 
 /* ================================================== */
 
+static void
+parse_smoothtime(char *line)
+{
+  check_number_of_args(line, 2);
+  if (sscanf(line, "%lf %lf", &smooth_max_freq, &smooth_max_wander) != 2) {
+    smooth_max_freq = 0.0;
+    command_parse_error();
+  }
+}
+
+/* ================================================== */
 static void
 parse_tempcomp(char *line)
 {
@@ -1719,6 +1737,15 @@ CNF_GetLockMemory(void)
 
 /* ================================================== */
 
+void
+CNF_GetSmooth(double *max_freq, double *max_wander)
+{
+  *max_freq = smooth_max_freq;
+  *max_wander = smooth_max_wander;
+}
+
+/* ================================================== */
+
 void
 CNF_GetTempComp(char **file, double *interval, char **point_file, double *T0, double *k0, double *k1, double *k2)
 {
diff --git a/conf.h b/conf.h
index ab415fa..4d51d24 100644
--- a/conf.h
+++ b/conf.h
@@ -96,6 +96,7 @@ extern void CNF_SetupAccessRestrictions(void);
 extern int CNF_GetSchedPriority(void);
 extern int CNF_GetLockMemory(void);
 
+extern void CNF_GetSmooth(double *max_freq, double *max_wander);
 extern void CNF_GetTempComp(char **file, double *interval, char **point_file, double *T0, double *k0, double *k1, double *k2);
 
 extern char *CNF_GetUser(void);
diff --git a/local.c b/local.c
index 6c337b2..f3dd15f 100644
--- a/local.c
+++ b/local.c
@@ -36,6 +36,7 @@
 #include "local.h"
 #include "localp.h"
 #include "memory.h"
+#include "smooth.h"
 #include "util.h"
 #include "logging.h"
 
@@ -488,6 +489,9 @@ LCL_ApplyStepOffset(double offset)
 
   (*drv_apply_step_offset)(offset);
 
+  /* Reset smoothing on all clock steps */
+  SMT_Reset(&cooked);
+
   /* Dispatch to all handlers */
   invoke_parameter_change_handlers(&raw, &cooked, 0.0, offset, LCL_ChangeStep);
 }
diff --git a/logging.h b/logging.h
index aa221fa..6dfb774 100644
--- a/logging.h
+++ b/logging.h
@@ -100,7 +100,8 @@ typedef enum {
   LOGF_SysWinnt,
   LOGF_TempComp,
   LOGF_RtcLinux,
-  LOGF_Refclock
+  LOGF_Refclock,
+  LOGF_Smooth,
 } LOG_Facility;
 
 /* Init function */
diff --git a/main.c b/main.c
index 3d78a72..f63ea1c 100644
--- a/main.c
+++ b/main.c
@@ -49,6 +49,7 @@
 #include "refclock.h"
 #include "clientlog.h"
 #include "nameserv.h"
+#include "smooth.h"
 #include "tempcomp.h"
 
 /* ================================================== */
@@ -88,6 +89,7 @@ MAI_CleanupAndExit(void)
   /* Don't update clock when removing sources */
   REF_SetMode(REF_ModeIgnore);
 
+  SMT_Finalise();
   TMC_Finalise();
   MNL_Finalise();
   CLG_Finalise();
@@ -495,6 +497,7 @@ int main
   CLG_Initialise();
   MNL_Initialise();
   TMC_Initialise();
+  SMT_Initialise();
 
   /* From now on, it is safe to do finalisation on exit */
   initialised = 1;
diff --git a/ntp_core.c b/ntp_core.c
index c62aafc..a24b516 100644
--- a/ntp_core.c
+++ b/ntp_core.c
@@ -36,6 +36,7 @@
 #include "sched.h"
 #include "reference.h"
 #include "local.h"
+#include "smooth.h"
 #include "sources.h"
 #include "util.h"
 #include "conf.h"
@@ -757,14 +758,14 @@ transmit_packet(NTP_Mode my_mode, /* The mode this machine wants to be */
 {
   NTP_Packet message;
   int leap, auth_len, length, ret;
-  struct timeval local_transmit;
+  struct timeval local_receive, local_transmit;
 
   /* Parameters read from reference module */
-  int are_we_synchronised, our_stratum;
+  int are_we_synchronised, our_stratum, smooth_time;
   NTP_Leap leap_status;
   uint32_t our_ref_id, ts_fuzz;
   struct timeval our_ref_time;
-  double our_root_delay, our_root_dispersion;
+  double our_root_delay, our_root_dispersion, smooth_offset;
 
   /* Don't reply with version higher than ours */
   if (version > NTP_VERSION) {
@@ -781,6 +782,15 @@ transmit_packet(NTP_Mode my_mode, /* The mode this machine wants to be */
                          &our_ref_id, &our_ref_time,
                          &our_root_delay, &our_root_dispersion);
 
+  /* Get current smoothing offset when sending packet to a client */
+  if (SMT_IsEnabled() && (my_mode == MODE_SERVER || my_mode == MODE_BROADCAST)) {
+    smooth_time = 1;
+    smooth_offset = SMT_GetOffset(&local_transmit);
+  } else {
+    smooth_time = 0;
+    smooth_offset = 0.0;
+  }
+
   if (are_we_synchronised) {
     leap = (int) leap_status;
   } else {
@@ -807,6 +817,14 @@ transmit_packet(NTP_Mode my_mode, /* The mode this machine wants to be */
   message.reference_id = htonl((NTP_int32) our_ref_id);
 
   /* Now fill in timestamps */
+
+  if (smooth_time) {
+    UTI_AddDoubleToTimeval(&our_ref_time, smooth_offset, &our_ref_time);
+    UTI_AddDoubleToTimeval(local_rx, smooth_offset, &local_receive);
+  } else {
+    local_receive = *local_rx;
+  }
+
   UTI_TimevalToInt64(&our_ref_time, &message.reference_ts, 0);
 
   /* Originate - this comes from the last packet the source sent us */
@@ -816,7 +834,7 @@ transmit_packet(NTP_Mode my_mode, /* The mode this machine wants to be */
      This timestamp will have been adjusted so that it will now look to
      the source like we have been running on our latest estimate of
      frequency all along */
-  UTI_TimevalToInt64(local_rx, &message.receive_ts, 0);
+  UTI_TimevalToInt64(&local_receive, &message.receive_ts, 0);
 
   /* Prepare random bits which will be added to the transmit timestamp. */
   ts_fuzz = UTI_GetNTPTsFuzz(message.precision);
@@ -826,6 +844,9 @@ transmit_packet(NTP_Mode my_mode, /* The mode this machine wants to be */
      from the source we're sending to now. */
   LCL_ReadCookedTime(&local_transmit, NULL);
 
+  if (smooth_time)
+    UTI_AddDoubleToTimeval(&local_transmit, smooth_offset, &local_transmit);
+
   length = NTP_NORMAL_PACKET_LENGTH;
 
   /* Authenticate */
diff --git a/reference.c b/reference.c
index df2f74f..2b1bfea 100644
--- a/reference.c
+++ b/reference.c
@@ -1241,6 +1241,14 @@ REF_GetOurStratum(void)
 
 /* ================================================== */
 
+double
+REF_GetSkew(void)
+{
+  return our_skew;
+}
+
+/* ================================================== */
+
 void
 REF_ModifyMaxupdateskew(double new_max_update_skew)
 {
diff --git a/reference.h b/reference.h
index d7f4874..2634c77 100644
--- a/reference.h
+++ b/reference.h
@@ -162,6 +162,9 @@ REF_SetUnsynchronised(void);
    synchronised */
 extern int REF_GetOurStratum(void);
 
+/* Return the current skew */
+extern double REF_GetSkew(void);
+
 /* Modify the setting for the maximum skew we are prepared to allow updates on (in ppm). */
 extern void REF_ModifyMaxupdateskew(double new_max_update_skew);
 
diff --git a/smooth.c b/smooth.c
new file mode 100644
index 0000000..dc145da
--- /dev/null
+++ b/smooth.c
@@ -0,0 +1,267 @@
+/*
+  chronyd/chronyc - Programs for keeping computer clocks accurate.
+
+ **********************************************************************
+ * Copyright (C) Miroslav Lichvar  2015
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of version 2 of the GNU General Public License as
+ * published by the Free Software Foundation.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License along
+ * with this program; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ * 
+ **********************************************************************
+
+  =======================================================================
+
+  Routines implementing time smoothing.
+
+  */
+
+#include "config.h"
+
+#include "sysincl.h"
+
+#include "conf.h"
+#include "local.h"
+#include "logging.h"
+#include "reference.h"
+#include "smooth.h"
+#include "util.h"
+
+/*
+  Time smoothing determines an offset that needs to be applied to the cooked
+  time to make it smooth for external observers.  Observed offset and frequency
+  change slowly and there are no discontinuities.  This can be used on an NTP
+  server to make it easier for the clients to track the time and keep their
+  clocks close together even when large offset or frequency corrections are
+  applied to the server's clock (e.g. after being offline for longer time).
+
+  Accumulated offset and frequency are smoothed out in three stages.  In the
+  first stage, the frequency is changed at a constant rate (wander) up to a
+  maximum, in the second stage the frequency stays at the maximum for as long
+  as needed and in the third stage the frequency is brought back to zero.
+
+              |
+    max_freq  +-------/--------\-------------
+              |      /|        |\
+        freq  |     / |        | \
+              |    /  |        |  \
+              |   /   |        |   \
+           0  +--/----+--------+----\--------
+              | /     |        |    |    time
+              |/      |        |    |
+
+        stage     1       2      3
+
+  Integral of this function is the smoothed out offset.  It's a continuous
+  piecewise polynomial with two quadratic parts and one linear.
+*/
+
+struct stage {
+  double wander;
+  double length;
+};
+
+#define NUM_STAGES 3
+
+static struct stage stages[NUM_STAGES];
+
+/* Enabled/disabled smoothing */
+static int enabled;
+
+/* Maximum skew/max_wander ratio to start updating offset and frequency */
+#define UNLOCK_SKEW_WANDER_RATIO 10000
+
+static int locked;
+
+/* Maximum wander and frequency offset */
+static double max_wander;
+static double max_freq;
+
+/* Frequency offset, time offset and the time of the last smoothing update */
+static double smooth_freq;
+static double smooth_offset;
+static struct timeval last_update;
+
+
+static void
+get_offset_freq(struct timeval *now, double *offset, double *freq)
+{
+  double elapsed, length;
+  int i;
+
+  UTI_DiffTimevalsToDouble(&elapsed, now, &last_update);
+
+  *offset = smooth_offset;
+  *freq = smooth_freq;
+
+  for (i = 0; i < NUM_STAGES; i++) {
+    if (elapsed <= 0.0)
+      break;
+
+    length = stages[i].length;
+    if (length >= elapsed)
+      length = elapsed;
+
+    *offset -= length * (2.0 * *freq + stages[i].wander * length) / 2.0;
+    *freq += stages[i].wander * length;
+    elapsed -= length;
+  }
+
+  if (elapsed > 0.0)
+    *offset -= elapsed * *freq;
+}
+
+static void
+update_stages(void)
+{
+  double s1, s2, s, l1, l2, l3, lc, f, f2;
+  int i, dir;
+
+  /* Prepare the three stages so that the integral of the frequency offset
+     is equal to the offset that should be smoothed out */
+
+  s1 = smooth_offset / max_wander;
+  s2 = smooth_freq * smooth_freq / (2.0 * max_wander * max_wander);
+  
+  l1 = l2 = l3 = 0.0;
+
+  /* Calculate the lengths of the 1st and 3rd stage assuming there is no
+     frequency limit.  If length of the 1st stage comes out negative, switch
+     its direction. */
+  for (dir = -1; dir <= 1; dir += 2) {
+    s = dir * s1 + s2;
+    if (s >= 0.0) {
+      l3 = sqrt(s);
+      l1 = l3 - dir * smooth_freq / max_wander;
+      if (l1 >= 0.0)
+        break;
+    }
+  }
+
+  assert(dir <= 1 && l1 >= 0.0 && l3 >= 0.0);
+
+  /* If the limit was reached, shorten 1st+3rd stages and set a 2nd stage */
+  f = dir * smooth_freq + l1 * max_wander - max_freq;
+  if (f > 0.0) {
+    lc = f / max_wander;
+
+    /* No 1st stage if the frequency is already above the maximum */
+    if (lc > l1) {
+      lc = l1;
+      f2 = dir * smooth_freq;
+    } else {
+      f2 = max_freq;
+    }
+
+    l2 = lc * (2.0 + f / f2);
+    l1 -= lc;
+    l3 -= lc;
+  }
+
+  stages[0].wander = dir * max_wander;
+  stages[0].length = l1;
+  stages[1].wander = 0.0;
+  stages[1].length = l2;
+  stages[2].wander = -dir * max_wander;
+  stages[2].length = l3;
+
+  for (i = 0; i < NUM_STAGES; i++) {
+    DEBUG_LOG(LOGF_Smooth, "Smooth stage %d wander %e length %f",
+              i + 1, stages[i].wander, stages[i].length);
+  }
+}
+
+static void
+update_smoothing(struct timeval *now, double offset, double freq)
+{
+  /* Don't accept offset/frequency until the clock has stabilized */
+  if (locked) {
+    if (REF_GetSkew() / max_wander < UNLOCK_SKEW_WANDER_RATIO) {
+      LOG(LOGS_INFO, LOGF_Smooth, "Time smoothing activated");
+      locked = 0;
+    }
+    return;
+  }
+
+  get_offset_freq(now, &smooth_offset, &smooth_freq);
+  smooth_offset += offset;
+  smooth_freq = (smooth_freq - freq) / (1.0 - freq);
+  last_update = *now;
+
+  update_stages();
+
+  DEBUG_LOG(LOGF_Smooth, "Smooth offset %e freq %e", smooth_offset, smooth_freq);
+}
+
+static void
+handle_slew(struct timeval *raw, struct timeval *cooked, double dfreq,
+            double doffset, LCL_ChangeType change_type, void *anything)
+{
+  double delta;
+
+  if (change_type == LCL_ChangeAdjust)
+    update_smoothing(cooked, doffset, dfreq);
+
+  UTI_AdjustTimeval(&last_update, cooked, &last_update, &delta, dfreq, doffset);
+}
+
+void SMT_Initialise(void)
+{
+  CNF_GetSmooth(&max_freq, &max_wander);
+  if (max_freq <= 0.0 || max_wander <= 0.0) {
+      enabled = 0;
+      return;
+  }
+
+  enabled = 1;
+  locked = 1;
+
+  /* Convert from ppm */
+  max_freq *= 1e-6;
+  max_wander *= 1e-6;
+
+  LCL_AddParameterChangeHandler(handle_slew, NULL);
+}
+
+void SMT_Finalise(void)
+{
+}
+
+int SMT_IsEnabled(void)
+{
+  return enabled;
+}
+
+double
+SMT_GetOffset(struct timeval *now)
+{
+  double offset, freq;
+
+  if (!enabled)
+    return 0.0;
+  
+  get_offset_freq(now, &offset, &freq);
+
+  return offset;
+}
+
+void
+SMT_Reset(struct timeval *now)
+{
+  if (!enabled)
+    return;
+
+  locked = 1;
+  smooth_offset = 0.0;
+  smooth_freq = 0.0;
+  last_update = *now;
+}
diff --git a/smooth.h b/smooth.h
new file mode 100644
index 0000000..c2877d6
--- /dev/null
+++ b/smooth.h
@@ -0,0 +1,40 @@
+/*
+  chronyd/chronyc - Programs for keeping computer clocks accurate.
+
+ **********************************************************************
+ * Copyright (C) Miroslav Lichvar  2015
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of version 2 of the GNU General Public License as
+ * published by the Free Software Foundation.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License along
+ * with this program; if not, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ * 
+ **********************************************************************
+
+  =======================================================================
+
+  This module implements time smoothing.
+  */
+
+#ifndef GOT_SMOOTH_H
+#define GOT_SMOOTH_H
+
+extern void SMT_Initialise(void);
+
+extern void SMT_Finalise(void);
+
+extern int SMT_IsEnabled(void);
+
+extern double SMT_GetOffset(struct timeval *now);
+
+extern void SMT_Reset(struct timeval *now);
+
+#endif