diff --git a/Registry/Registry.EM_COMMON.var b/Registry/Registry.EM_COMMON.var
index 6281c0f8e0..5df5eedcc0 100644
--- a/Registry/Registry.EM_COMMON.var
+++ b/Registry/Registry.EM_COMMON.var
@@ -347,6 +347,20 @@ state    real  QFX              ij      misc        1         -      irh      "Q
 state    real  REGIME           ij      misc        1         -      irh      "REGIME"  "FLAGS: 1=Night/Stable, 2=Mechanical Turbulent, 3=Forced Conv, 4=Free Conv" ""
 state integer  KPBL             ij      misc        1         -      irh      "KPBL"                  "LEVEL OF PBL TOP"                             ""
 
+# Increment Output
+state    real   u_iau        ijk     dyn_em      1         X     ih5       "U_IAU"    "x-wind component inc"           "m s-1"
+state    real   v_iau        ijk     dyn_em      1         Y     ih5       "V_IAU"    "y-wind component inc"           "m s-1"
+state    real   t_iau        ijk     dyn_em      1         -     ih5       "T_IAU"    "potential temp inc"             "K"
+state    real   w_iau        ijk     dyn_em      1         -     ih5       "W_IAU"    "z-wind component inc"           "m s-1"
+state    real   qv_iau       ijk     dyn_em      1         -     ih5       "QV_IAU"   "water water mixing ratio inc"   "kg kg-1"
+state    real   qc_iau       ijk     dyn_em      1         -     ih5       "QC_IAU"   "cloud water mixing ratio inc"   "kg kg-1"
+state    real   qr_iau       ijk     dyn_em      1         -     ih5       "QR_IAU"   "rain water mixing ratio inc"    "kg kg-1"
+state    real   qi_iau       ijk     dyn_em      1         -     ih5       "QI_IAU"   "ice water mixing ratio inc"     "kg kg-1"
+state    real   qs_iau       ijk     dyn_em      1         -     ih5       "QS_IAU"   "snow water mixing ratio inc"    "kg kg-1"
+state    real   qg_iau       ijk     dyn_em      1         -     ih5       "QG_IAU"   "graupel mixing ratio inc"       "kg kg-1"
+state    real   ph_iau       ijk     dyn_em      1         -     ih5       "PH_IAU"   "perturbation geopotential inc"  "m2 s-2"
+state    real   mu_iau        ij     dyn_em      1         -     ih5       "MU_IAU"   "dry air mass inc"     "pa"
+
 #
 #---------------------------------------------------------------------------------------------------------------------------------------
 #                                               
diff --git a/Registry/Registry.wrfvar b/Registry/Registry.wrfvar
index aa3c5bba64..85d694cdef 100644
--- a/Registry/Registry.wrfvar
+++ b/Registry/Registry.wrfvar
@@ -10,7 +10,7 @@ state   real    -              ijkft   g_scalar    1         -     -    -
 
 state   real   landmask            ij    misc          1     -     i012rhd=(interp_fcnm)u=(copy_fcnm)   "LANDMASK"      "LAND MASK (1 FOR LAND, 0 FOR WATER)"  ""
 
-state    real   SST              ij     misc        1         -     i01245rh05d=(interp_mask_water_field:lu_index,iswater)   "SST"              "SEA SURFACE TEMPERATURE" "K"
+state    real   SST              ij     misc        1         -     i01245rd=(interp_mask_water_field:lu_index,iswater)   "SST"              "SEA SURFACE TEMPERATURE" "K"
 
 # Registry entries that are exclusive to Registry.EM
 
diff --git a/Registry/registry.em_shared_collection b/Registry/registry.em_shared_collection
index 484514c482..28a75f175d 100644
--- a/Registry/registry.em_shared_collection
+++ b/Registry/registry.em_shared_collection
@@ -28,4 +28,5 @@ include registry.new3d_wif
 include registry.trad_fields
 include registry.solar_fields
 include registry.diags
+include registry.iau
 include registry.CMAQ
diff --git a/Registry/registry.iau b/Registry/registry.iau
new file mode 100644
index 0000000000..882dfba9ce
--- /dev/null
+++ b/Registry/registry.iau
@@ -0,0 +1,38 @@
+# IAU variables
+
+state   character iau_time     -         -        -      -     i{15}r  "TIME_IAU"     "  "                                   "   "
+state   real   mu_iau         ij      misc        1      -     i{15}r  "MU_IAU"       "mu analysis increments array"         "   "
+state   real   u_iau          ikj     misc        1      -     i{15}r  "U_IAU"        "u analysis increments array"          "   "
+state   real   v_iau          ikj     misc        1      -     i{15}r  "V_IAU"        "v analysis increments array"          "   "
+state   real   w_iau          ikj     misc        1      -     i{15}r  "W_IAU"        "w analysis increments array"          "   "
+state   real   p_iau          ikj     misc        1      -     i{15}r  "P_IAU"        "p analysis increments array"          "   "
+state   real   t_iau          ikj     misc        1      -     i{15}r  "T_IAU"        "t analysis increments array"          "   "
+state   real   ph_iau         ikj     misc        1      -     i{15}r  "PH_IAU"       "ph analysis increments array"         "   "
+state   real   qv_iau         ikj     misc        1      -     i{15}r  "QV_IAU"       "qv analysis increments array"         "   "
+state   real   qr_iau         ikj     misc        1      -     i{15}r  "QR_IAU"       "qr analysis increments array"         "   "
+state   real   qc_iau         ikj     misc        1      -     i{15}r  "QC_IAU"       "qc analysis increments array"         "   "
+state   real   qs_iau         ikj     misc        1      -     i{15}r  "QS_IAU"       "qs analysis increments array"         "   "
+state   real   qi_iau         ikj     misc        1      -     i{15}r  "QI_IAU"       "qice analysis increments array"       "   "
+state   real   qg_iau         ikj     misc        1      -     i{15}r  "QG_IAU"       "qgraupel analysis increments array"   "   "
+
+state   real   RUIAUTEN       ikj     misc        1      X      r   "RUIAUTEN"     "X WIND TENDENCY DUE TO IAU"  "m s-2" 
+state   real   RVIAUTEN       ikj     misc        1      Y      r   "RVIAUTEN"     "Y WIND TENDENCY DUE TO IAU"  "m s-2"
+state   real   RTHIAUTEN      ikj     misc        1      -      r   "RTHIAUTEN"    "THETA TENDENCY DUE TO IAU"   "K s-1"
+state   real   RPHIAUTEN      ikj     misc        1      -      r   "RPHIAUTEN"    "GEOPOTENTIAL TENDENCY DUE TO IAU"   "m2 s-3" 
+state   real   RQVIAUTEN      ikj     misc        1      -      r   "RQVIAUTEN"    "Q_V TENDENCY DUE TO IAU"     "kg kg-1 s-1" 
+state   real   RQCIAUTEN      ikj     misc        1      -      r   "RQCIAUTEN"    "Q_C TENDENCY DUE TO IAU"     "kg kg-1 s-1" 
+state   real   RQRIAUTEN      ikj     misc        1      -      r   "RQRIAUTEN"    "Q_R TENDENCY DUE TO IAU"     "kg kg-1 s-1" 
+state   real   RQIIAUTEN      ikj     misc        1      -      r   "RQIIAUTEN"    "Q_I TENDENCY DUE TO IAU"     "kg kg-1 s-1" 
+state   real   RQSIAUTEN      ikj     misc        1      -      r   "RQSIAUTEN"    "Q_S TENDENCY DUE TO IAU"     "kg kg-1 s-1" 
+state   real   RQGIAUTEN      ikj     misc        1      -      r   "RQGIAUTEN"    "Q_G TENDENCY DUE TO IAU"     "kg kg-1 s-1" 
+state   real   RMUIAUTEN      ij      misc        1      -      r   "RMUIAUTEN"    "MU TENDENCY DUE TO IAU"      "Pa s-1" 
+
+# IAU namelist options
+
+rconfig   integer iau                   namelist,time_control     max_domains    0      irh   "analysis increments read" "0/1 ACTIVATE FOR ANALYSIS INCREMENTS UPDATES" ""
+rconfig   real    iau_time_window_sec   namelist,time_control     max_domains    3600.  irh   "iau_time_window_sec"      "TIME WINDOW OF INCREMENTS ANALYSIS UPDATES" "SECONDS"
+
+# IAU packages
+
+package   noiau        iau==0   -   -
+package   iau          iau==1   -   state:u_iau,v_iau,w_iau,p_iau,t_iau,ph_iau,qqv_iau,qqr_iau,qqc_iau,qqs_iau,qqi_iau,qqg_iau,ruiauten,rviauten,rthiauten,rqviauten,rqciauten,rqriauten,rqiiauten,rqsiauten,rqgiauten
diff --git a/dyn_em/module_em.F b/dyn_em/module_em.F
index 56df890f90..62abb8b7c3 100644
--- a/dyn_em/module_em.F
+++ b/dyn_em/module_em.F
@@ -2089,6 +2089,14 @@ SUBROUTINE calculate_phy_tend (config_flags,c1,c2,                     &
                      RQISHTEN,RQSSHTEN,RQGSHTEN,                       &
                      RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,        &
                      RMUNDGDTEN,                                       &
+                     ruiauten,rviauten,rthiauten,                      &
+                     rqviauten,rqciauten,rqriauten,                    &
+                     rqiiauten,rqsiauten,rqgiauten,                    &
+                     rmuiauten,rphiauten,                              &
+                     u_iau,v_iau,t_iau,ph_iau,                         &
+                     mu_iau,qv_iau,qc_iau,qr_iau,                      &
+                     qi_iau,qs_iau,qg_iau,                             &
+                     itimestep,                                        &
                      scalar, scalar_tend, num_scalar,                  &
                      tracer, tracer_tend, num_tracer,                  &
                      ids,ide, jds,jde, kds,kde,                        &
@@ -2163,6 +2171,38 @@ SUBROUTINE calculate_phy_tend (config_flags,c1,c2,                     &
       REAL,     DIMENSION( ims:ime, jms:jme )               , &
                 INTENT(INOUT)   ::                           RMUNDGDTEN
 
+! iau
+      REAL,     DIMENSION( ims:ime , kms:kme , jms:jme ),              &
+                INTENT(INOUT)   ::                           RUIAUTEN, &
+                                                             RVIAUTEN, &
+                                                            RTHIAUTEN, &
+                                                            RPHIAUTEN, &
+                                                            RQVIAUTEN, &
+                                                            RQCIAUTEN, &
+                                                            RQRIAUTEN, &
+                                                            RQIIAUTEN, &
+                                                            RQSIAUTEN, &
+                                                            RQGIAUTEN   
+
+      REAL,     DIMENSION( ims:ime, jms:jme ),                         &
+                INTENT(INOUT)   ::                          RMUIAUTEN
+
+      REAL,     DIMENSION( ims:ime , kms:kme , jms:jme ),              &
+                INTENT(IN   )   ::                                     &
+                                                                u_iau, &
+                                                                v_iau, &
+                                                                t_iau, &
+                                                               ph_iau, &
+                                                               qv_iau, &
+                                                               qc_iau, &
+                                                               qr_iau, &
+                                                               qs_iau, &
+                                                               qi_iau, &
+                                                               qg_iau   
+
+      REAL,     DIMENSION( ims:ime, jms:jme ),                         &
+                INTENT(IN   )   ::                             mu_iau
+
 ! 4d arrays
 
     REAL    ,DIMENSION(ims:ime,kms:kme,jms:jme,num_tracer),INTENT(INOUT)   :: tracer
@@ -2170,9 +2210,13 @@ SUBROUTINE calculate_phy_tend (config_flags,c1,c2,                     &
     REAL    ,DIMENSION(ims:ime,kms:kme,jms:jme,num_scalar),INTENT(INOUT)   :: scalar
     REAL    ,DIMENSION(ims:ime,kms:kme,jms:jme,num_scalar),INTENT(INOUT)   :: scalar_tend
 
+      INTEGER :: itimestep
       INTEGER :: i,k,j, im
       INTEGER :: itf,ktf,jtf,itsu,jtsv
 
+! local
+    real   :: wgt_iau
+
 !-----------------------------------------------------------------------
 
 !<DESCRIPTION>
@@ -2448,6 +2492,88 @@ SUBROUTINE calculate_phy_tend (config_flags,c1,c2,                     &
          ENDDO
    ENDDO
 
+! IAU (Incremental Analysis Updates)
+
+   IF (config_flags%iau .gt. 0) THEN
+
+      wgt_iau = iau_coef(config_flags%iau_time_window_sec,config_flags%dt,itimestep)
+
+      DO J=jts,jtf
+      DO K=kts,ktf
+      DO I=its,itf
+         RUIAUTEN(I,K,J) =(c1(k)*mut(I,J)+c2(k))*U_IAU(I,K,J)*WGT_IAU
+         RVIAUTEN(I,K,J) =(c1(k)*mut(I,J)+c2(k))*V_IAU(I,K,J)*WGT_IAU
+         RTHIAUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*T_IAU(I,K,J)*WGT_IAU
+         RQVIAUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*QV_IAU(I,K,J)*WGT_IAU
+      ENDDO
+      ENDDO
+      ENDDO
+
+      DO J=jts,jtf
+      DO K=kts,ktf+1
+      DO I=its,itf
+         RPHIAUTEN(I,K,J)=PH_IAU(I,K,J)*WGT_IAU
+      ENDDO 
+      ENDDO
+      ENDDO
+
+      DO J=jts,jtf
+      DO I=its,itf
+         RMUIAUTEN(I,J) =MU_IAU(I,J)*WGT_IAU
+      ENDDO
+      ENDDO
+
+      IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
+         DO J=jts,jtf
+         DO K=kts,ktf
+         DO I=its,itf
+            RQCIAUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*QC_IAU(I,K,J)*WGT_IAU
+         ENDDO
+         ENDDO
+         ENDDO
+      ENDIF
+
+      IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
+         DO J=jts,jtf
+         DO K=kts,ktf
+         DO I=its,itf
+            RQRIAUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*QR_IAU(I,K,J)*WGT_IAU
+         ENDDO
+         ENDDO
+         ENDDO
+      ENDIF
+
+      IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
+         DO J=jts,jtf
+         DO K=kts,ktf
+         DO I=its,itf
+            RQIIAUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*QI_IAU(I,K,J)*WGT_IAU
+         ENDDO
+         ENDDO
+         ENDDO
+      ENDIF
+
+      IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
+         DO J=jts,jtf
+         DO K=kts,ktf
+         DO I=its,itf
+            RQSIAUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*QS_IAU(I,K,J)*WGT_IAU
+         ENDDO
+         ENDDO
+         ENDDO
+      ENDIF
+
+      IF(P_QG .ge. PARAM_FIRST_SCALAR)THEN
+         DO J=jts,jtf
+         DO K=kts,ktf
+         DO I=its,itf
+            RQGIAUTEN(I,K,J)=(c1(k)*mut(I,J)+c2(k))*QG_IAU(I,K,J)*WGT_IAU
+         ENDDO
+         ENDDO
+         ENDDO
+      ENDIF
+
+   ENDIF
 
 END SUBROUTINE calculate_phy_tend
 
@@ -3018,4 +3144,22 @@ subroutine trajmapproj (grid,config_flags,ts_proj)
 
 end subroutine trajmapproj
 
+FUNCTION iau_coef (iau_time_window_sec, dt, itimestep)  result(wgt_iau)
+
+  implicit none
+
+! This function returns the coeficient for IAU
+
+      INTEGER :: itimestep
+      INTEGER :: nsteps_iau
+      REAL    :: wgt_iau, iau_time_window_sec, dt
+
+      wgt_iau = 0.
+      nsteps_iau = nint(iau_time_window_sec / dt)
+      if (itimestep <= nsteps_iau) then
+         wgt_iau = 1.0/iau_time_window_sec
+      endif
+
+END FUNCTION iau_coef
+
 END MODULE module_em
diff --git a/dyn_em/module_first_rk_step_part2.F b/dyn_em/module_first_rk_step_part2.F
index 8f092ce794..656f063f7c 100644
--- a/dyn_em/module_first_rk_step_part2.F
+++ b/dyn_em/module_first_rk_step_part2.F
@@ -402,6 +402,14 @@ SUBROUTINE first_rk_step_part2 (   grid , config_flags        &
                      grid%rqishten,grid%rqsshten,grid%rqgshten,        &
                      grid%RUNDGDTEN,grid%RVNDGDTEN,grid%RTHNDGDTEN,grid%RQVNDGDTEN, &
                      grid%RMUNDGDTEN,                                  &
+                     grid%ruiauten,grid%rviauten,grid%rthiauten,       &
+                     grid%rqviauten,grid%rqciauten,grid%rqriauten,     &
+                     grid%rqiiauten,grid%rqsiauten,grid%rqgiauten,     &
+                     grid%rmuiauten,grid%rphiauten,                    &
+                     grid%u_iau,grid%v_iau,grid%t_iau,grid%ph_iau,     &
+                     grid%mu_iau,grid%qv_iau,grid%qc_iau,grid%qr_iau,  &
+                     grid%qs_iau,grid%qi_iau,grid%qg_iau,              &
+                     grid%itimestep,                                   &
                      scalar, scalar_tend, num_scalar,                  &
                      tracer, tracer_tend, num_tracer,                  &
                      ids,ide, jds,jde, kds,kde,                        &
@@ -807,6 +815,10 @@ SUBROUTINE first_rk_step_part2 (   grid , config_flags        &
                            grid%RVNDGDTEN,grid%RTHNDGDTEN,grid%RPHNDGDTEN,  &
                            grid%RQVNDGDTEN,grid%RMUNDGDTEN,                 &
                            grid%rthfrten,grid%rqvfrten,                     &  ! fire
+                           grid%ruiauten,grid%rviauten,grid%rthiauten,      &
+                           grid%rqviauten,grid%rqciauten,grid%rqriauten,    &
+                           grid%rqiiauten,grid%rqsiauten,grid%rqgiauten,    &
+                           grid%rphiauten,grid%rmuiauten,                   &
                            num_moist,num_scalar,config_flags,rk_step,       &
                            grid%adv_moist_cond,                             &
                            ids, ide, jds, jde, kds, kde,                    &
diff --git a/phys/module_physics_addtendc.F b/phys/module_physics_addtendc.F
index 5e15a32b83..a0a360deda 100644
--- a/phys/module_physics_addtendc.F
+++ b/phys/module_physics_addtendc.F
@@ -38,6 +38,10 @@ SUBROUTINE update_phy_ten(rph_tendf,rt_tendf,ru_tendf,rv_tendf,moist_tendf, &
                       RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RPHNDGDTEN,            &
                       RQVNDGDTEN,RMUNDGDTEN,                                &
                       rthfrten,rqvfrten,                                    & !fire
+                      ruiauten,rviauten,rthiauten,                          &
+                      rqviauten,rqciauten,rqriauten,                        &
+                      rqiiauten,rqsiauten,rqgiauten,                        &
+                      rphiauten,rmuiauten,                                  &
                       n_moist,n_scalar,config_flags,rk_step,adv_moist_cond, &
                       ids, ide, jds, jde, kds, kde,                         &
                       ims, ime, jms, jme, kms, kme,                         &
@@ -103,9 +107,20 @@ SUBROUTINE update_phy_ten(rph_tendf,rt_tendf,ru_tendf,rv_tendf,moist_tendf, &
                                                      RPHNDGDTEN, &
                                                      RQVNDGDTEN, &
                                                       RUNDGDTEN, &
-                                                      RVNDGDTEN
-
-   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN  ) :: RMUNDGDTEN
+                                                      RVNDGDTEN, &
+                                                       RUIAUTEN, &
+                                                       RVIAUTEN, &
+                                                      RPHIAUTEN, &
+                                                      RTHIAUTEN, &
+                                                      RQVIAUTEN, &
+                                                      RQCIAUTEN, &
+                                                      RQRIAUTEN, &
+                                                      RQIIAUTEN, &
+                                                      RQSIAUTEN, &
+                                                      RQGIAUTEN
+
+   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN  ) :: RMUNDGDTEN, &
+                                                       RMUIAUTEN
 
    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN  ) ::   & ! fire
                                                        rthfrten, &
@@ -164,6 +179,18 @@ SUBROUTINE update_phy_ten(rph_tendf,rt_tendf,ru_tendf,rv_tendf,moist_tendf, &
                       ims, ime, jms, jme, kms, kme,            &
                       its, ite, jts, jte, kts, kte             )
 
+   if (config_flags%iau .gt. 0)                                &
+      CALL phy_iau_ten(config_flags,rk_step,n_moist,           &
+                      rph_tendf,rt_tendf,ru_tendf,rv_tendf,    &
+                      mu_tendf, moist_tendf,                   &
+                      RUIAUTEN,RVIAUTEN,RTHIAUTEN,             &
+                      RPHIAUTEN,RMUIAUTEN,                     &
+                      RQVIAUTEN,RQCIAUTEN,RQRIAUTEN,           &
+                      RQIIAUTEN,RQSIAUTEN,RQGIAUTEN,           &
+                      ids, ide, jds, jde, kds, kde,            &
+                      ims, ime, jms, jme, kms, kme,            &
+                      its, ite, jts, jte, kts, kte             )
+
    if (config_flags%ifire .gt. 0)                              & ! fire
       CALL phy_fr_ten(config_flags,rk_step,n_moist,            &
                       rt_tendf,ru_tendf,rv_tendf,              &
@@ -1992,6 +2019,129 @@ SUBROUTINE phy_fg_ten(config_flags,rk_step,n_moist,            &
 
 END SUBROUTINE phy_fg_ten
 
+!=================================================================
+SUBROUTINE phy_iau_ten(config_flags,rk_step,n_moist,             &
+                      rph_tendf,rt_tendf,ru_tendf,rv_tendf,      &
+                      mu_tendf, moist_tendf,                     &
+                      RUIAUTEN,RVIAUTEN,RTHIAUTEN,               &
+                      RPHIAUTEN,RMUIAUTEN,                       &
+                      RQVIAUTEN,RQCIAUTEN,RQRIAUTEN,             &
+                      RQIIAUTEN,RQSIAUTEN,RQGIAUTEN,             &
+                      ids, ide, jds, jde, kds, kde,              &
+                      ims, ime, jms, jme, kms, kme,              &
+                      its, ite, jts, jte, kts, kte               )
+!-----------------------------------------------------------------
+   IMPLICIT NONE
+!-----------------------------------------------------------------
+   TYPE(grid_config_rec_type) ,     INTENT(IN   ) :: config_flags
+
+   INTEGER , INTENT(IN)        ::  ids, ide, jds, jde, kds, kde, &
+                                   ims, ime, jms, jme, kms, kme, &
+                                   its, ite, jts, jte, kts, kte, &
+                                   n_moist, rk_step
+
+   REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_moist),         &
+          INTENT(INOUT)     ::                      moist_tendf
+
+   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN  ) ::   &
+                                                       RPHIAUTEN, &
+                                                       RTHIAUTEN, &
+                                                        RUIAUTEN, &
+                                                        RVIAUTEN, &
+                                                       RQVIAUTEN
+   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN  ) ::    &
+                                                       RQCIAUTEN, &
+                                                       RQRIAUTEN, &
+                                                       RQIIAUTEN, &
+                                                       RQSIAUTEN, &
+                                                       RQGIAUTEN
+
+   REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN  ) ::  RMUIAUTEN
+
+   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT)::   &
+                                                       rph_tendf,&
+                                                       rt_tendf, &
+                                                       ru_tendf, &
+                                                       rv_tendf
+
+   REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT)::  mu_tendf
+
+! LOCAL VARS
+
+   INTEGER :: i,j,k,IBGN,IEND,JBGN,JEND
+
+!-----------------------------------------------------------------
+
+           CALL add_a2a(rt_tendf,RTHIAUTEN,config_flags,         &
+                ids,ide, jds, jde, kds, kde,                     &
+                ims, ime, jms, jme, kms, kme,                    &
+                its, ite, jts, jte, kts, kte                     )
+
+           CALL add_a2a_ph(rph_tendf,RPHIAUTEN,config_flags,     &
+                ids,ide, jds, jde, kds, kde,                     &
+                ims, ime, jms, jme, kms, kme,                    &
+                its, ite, jts, jte, kts, kte                     )
+
+! note fdda u and v tendencies are staggered
+           CALL add_c2c_u(ru_tendf,RUIAUTEN,config_flags,        &
+                ids,ide, jds, jde, kds, kde,                     &
+                ims, ime, jms, jme, kms, kme,                    &
+                its, ite, jts, jte, kts, kte                     )
+
+           CALL add_c2c_v(rv_tendf,RVIAUTEN,config_flags,        &
+                ids,ide, jds, jde, kds, kde,                     &
+                ims, ime, jms, jme, kms, kme,                    &
+                its, ite, jts, jte, kts, kte                     )
+
+           CALL add_a2a(mu_tendf,RMUIAUTEN,config_flags,         &
+                ids,ide, jds, jde, kds, kds,                     &
+                ims, ime, jms, jme, kms, kms,                    &
+                its, ite, jts, jte, kts, kts                     )
+
+        if (P_QV .ge. PARAM_FIRST_SCALAR)                        &
+           CALL add_a2a(moist_tendf(ims,kms,jms,P_QV),RQVIAUTEN, &
+                config_flags,                                    &
+                ids,ide, jds, jde, kds, kde,                     &
+                ims, ime, jms, jme, kms, kme,                    &
+                its, ite, jts, jte, kts, kte                     )
+
+        if (P_QC .ge. PARAM_FIRST_SCALAR)                        &
+           CALL add_a2a(moist_tendf(ims,kms,jms,P_QC),RQCIAUTEN, &
+                config_flags,                                    &
+                ids,ide, jds, jde, kds, kde,                     &
+                ims, ime, jms, jme, kms, kme,                    &
+                its, ite, jts, jte, kts, kte                     )
+
+        if (P_QR .ge. PARAM_FIRST_SCALAR)                        &
+           CALL add_a2a(moist_tendf(ims,kms,jms,P_QR),RQRIAUTEN, &
+                config_flags,                                    &
+                ids,ide, jds, jde, kds, kde,                     &
+                ims, ime, jms, jme, kms, kme,                    &
+                its, ite, jts, jte, kts, kte                     )
+
+        if (P_QI .ge. PARAM_FIRST_SCALAR)                        &
+           CALL add_a2a(moist_tendf(ims,kms,jms,P_QI),RQIIAUTEN, &
+                config_flags,                                    &
+                ids,ide, jds, jde, kds, kde,                     &
+                ims, ime, jms, jme, kms, kme,                    &
+                its, ite, jts, jte, kts, kte                     )
+
+        if (P_QS .ge. PARAM_FIRST_SCALAR)                        &
+           CALL add_a2a(moist_tendf(ims,kms,jms,P_QS),RQSIAUTEN, &
+                config_flags,                                    &
+                ids,ide, jds, jde, kds, kde,                     &
+                ims, ime, jms, jme, kms, kme,                    &
+                its, ite, jts, jte, kts, kte                     )
+
+        if (P_QG .ge. PARAM_FIRST_SCALAR)                        &
+           CALL add_a2a(moist_tendf(ims,kms,jms,P_QG),RQGIAUTEN, &
+                config_flags,                                    &
+                ids,ide, jds, jde, kds, kde,                     &
+                ims, ime, jms, jme, kms, kme,                    &
+                its, ite, jts, jte, kts, kte                     )
+
+END SUBROUTINE phy_iau_ten
+
 !=================================================================
 SUBROUTINE phy_fr_ten(config_flags,rk_step,n_moist,            &
                       rt_tendf,ru_tendf,rv_tendf,              &
diff --git a/var/build/da.make b/var/build/da.make
index 898d57464c..522a013f78 100644
--- a/var/build/da.make
+++ b/var/build/da.make
@@ -445,10 +445,10 @@ da_setup_structures.o :
          fi
 	if $(FGREP) '!$$OMP' $*.f ; then \
           if [ -n "$(OMP)" ] ; then echo COMPILING $*.f90 WITH OMP ; fi ; \
-	  $(FC) -c $(FCFLAGS) $(OMP) $(PROMOTION) $(CRTM_SRC) $(RTTOV_SRC) $(HDF5_INC) $*.f ; \
+	  $(FC) -c $(FCFLAGS) $(OMP) $(PROMOTION) $(CRTM_SRC) $(RTTOV_SRC) $(HDF5_INC) -I$(NETCDF)/include  $*.f ; \
         else \
 	if [ -n "$(OMP)" ] ; then echo COMPILING $*.f90 WITHOUT OMP ; fi ; \
-	  $(FC) -c $(FCFLAGS) $(PROMOTION) $(CRTM_SRC) $(RTTOV_SRC) $(HDF5_INC) $*.f ; \
+	  $(FC) -c $(FCFLAGS) $(PROMOTION) $(CRTM_SRC) $(RTTOV_SRC) $(HDF5_INC) -I$(NETCDF)/include $*.f ; \
         fi
 
 da_obs_io.o \
diff --git a/var/build/depend.txt b/var/build/depend.txt
index 3632d9e20e..8ad937e19c 100644
--- a/var/build/depend.txt
+++ b/var/build/depend.txt
@@ -380,7 +380,7 @@ da_wavelet.o : da_wavelet.f90 da_transform_through_wavelet_adj.inc da_transform_
 da_wrf_interfaces.o : da_wrf_interfaces.f90 module_configure.o module_domain.o 
 da_wrfvar_esmf.o : da_wrfvar_esmf.f90 
 da_wrfvar_esmf_super.o : da_wrfvar_esmf_super.f90 da_wrfvar_interface.inc da_esmf_finalize.inc da_esmf_run.inc da_esmf_init.inc 
-da_wrfvar_io.o : copyfile.c da_wrfvar_io.f90 da_med_initialdata_output_lbc.inc da_med_initialdata_output.inc da_med_initialdata_input.inc da_update_firstguess.inc da_4dvar.o da_tracing.o da_reporting.o da_control.o module_io_domain.o module_domain.o module_date_time.o module_configure.o module_domain_type.o 
+da_wrfvar_io.o : copyfile.c da_wrfvar_io.f90 da_med_initialdata_output_lbc.inc da_med_initialdata_output.inc da_med_initialdata_input.inc da_update_firstguess.inc da_write_anaincrements.inc da_4dvar.o da_tracing.o da_reporting.o da_control.o module_io_domain.o module_domain.o module_date_time.o module_configure.o module_domain_type.o 
 da_wrfvar_main.o : da_wrfvar_main.f90 da_4dvar.o da_wrfvar_top.o da_wrf_interfaces.o da_tracing.o da_control.o module_symbols_util.o 
 da_wrfvar_top.o : da_wrfvar_top.f90 da_solve_init.inc da_solve_dual_res_init.inc da_solve.inc da_wrfvar_finalize.inc da_wrfvar_interface.inc da_wrfvar_run.inc da_wrfvar_init2.inc da_wrfvar_init1.inc da_wrf_interfaces.o da_rain.o da_synop.o da_ssmi.o da_sound.o da_ships.o da_satem.o da_radar.o da_lightning.o da_mtgirs.o da_qscat.o da_profiler.o da_polaramv.o da_pilot.o da_metar.o da_gpsref.o da_gpspw.o da_geoamv.o da_buoy.o da_bogus.o da_airsr.o da_airep.o da_crtm.o da_tools.o da_vtox_transforms.o da_transfer_model.o da_tracing.o da_tools_serial.o da_test.o da_setup_structures.o da_reporting.o da_varbc.o da_radiance1.o da_physics.o da_par_util.o da_obs_io.o da_obs.o da_minimisation.o da_define_structures.o da_control.o module_comm_dm.o module_dm.o module_tiles.o module_state_description.o module_radiance.o da_wrfvar_io.o da_4dvar.o module_symbols_util.o module_driver_constants.o module_domain.o module_configure.o module_io_domain.o da_netcdf_interface.o da_gpseph.o da_varbc_tamdar.o module_io_wrf.o da_chem_sfc.o 
 decode_airs.o : decode_airs.f90 module_read_airs.o 
diff --git a/var/da/da_main/da_solve.inc b/var/da/da_main/da_solve.inc
index a2ce91ccba..9e2caa2548 100644
--- a/var/da/da_main/da_solve.inc
+++ b/var/da/da_main/da_solve.inc
@@ -1129,6 +1129,9 @@
        (config_flags%real_data_init_type == 3)) then
       if ( .not. offline_varbc ) &
       call da_update_firstguess(input_grid)
+   ! output analysis increments
+      call da_write_anaincrements (grid, config_flags)
+
 #ifdef VAR4D
       !if (var4d) call da_med_initialdata_output_lbc (head_grid , config_flags)
       if ( var4d_lbc ) then
diff --git a/var/da/da_main/da_wrfvar_io.f90 b/var/da/da_main/da_wrfvar_io.f90
index eca072a404..c95952656b 100644
--- a/var/da/da_main/da_wrfvar_io.f90
+++ b/var/da/da_main/da_wrfvar_io.f90
@@ -6,6 +6,7 @@ module da_wrfvar_io
    use module_domain, only : domain, get_ijk_from_grid
    use module_io_domain, only : open_r_dataset,close_dataset, &
       input_input, open_w_dataset,output_input, &
+      output_auxhist5, &
       input_boundary, output_boundary, output_auxhist4, &
       input_auxhist6, input_auxhist4
    use module_io, only: wrf_get_dom_ti_integer
@@ -30,5 +31,6 @@ module da_wrfvar_io
 #include "da_med_initialdata_output.inc"
 #include "da_med_initialdata_output_lbc.inc"
 #include "da_update_firstguess.inc"
+#include "da_write_anaincrements.inc"
 
 end module da_wrfvar_io
diff --git a/var/da/da_main/da_wrfvar_top.f90 b/var/da/da_main/da_wrfvar_top.f90
index 240ceb5be0..3d0ccd78a9 100644
--- a/var/da/da_main/da_wrfvar_top.f90
+++ b/var/da/da_main/da_wrfvar_top.f90
@@ -95,7 +95,7 @@ module da_wrfvar_top
    use da_vtox_transforms, only : da_transform_vtox, da_transform_xtoxa, &
       da_transform_xtoxa_adj, da_copy_xa, da_add_xa, da_transform_vpatox, &
       da_transform_vtox_inv
-   use da_wrfvar_io, only : da_med_initialdata_input, da_update_firstguess
+   use da_wrfvar_io, only : da_med_initialdata_input, da_update_firstguess,da_write_anaincrements
    use da_tools, only : da_set_randomcv, da_get_julian_time
 
    use da_tools, only : map_info,map_info_ens,proj_merc, proj_ps,proj_lc,proj_latlon, &
diff --git a/var/da/da_main/da_write_anaincrements.inc b/var/da/da_main/da_write_anaincrements.inc
new file mode 100644
index 0000000000..7a351ccc52
--- /dev/null
+++ b/var/da/da_main/da_write_anaincrements.inc
@@ -0,0 +1,76 @@
+subroutine da_write_anaincrements(grid, config_flags)
+
+  use module_domain, only : get_ijk_from_grid, program_name
+  use da_control,    only : use_radarobs, use_rad, crtm_cloud, &
+                            use_radar_rhv, use_radar_rqv
+  use module_state_description, only : p_qv, p_qc, p_qr, p_qi, &
+                                       p_qs, p_qg, &
+                                       f_qc, f_qr, f_qi, f_qs, f_qg
+  use module_model_constants, only: R_d, R_v, T0
+#if (WRF_CHEM == 1)
+  use module_domain_type, only : fieldlist
+  use da_control,    only : use_chemic_surfobs, stdout
+  use module_state_description, only : PARAM_FIRST_SCALAR, num_chem
+#endif
+
+  implicit none
+
+  INTERFACE
+    integer(c_int32_t) function copyfile(ifile, ofile) bind(c)
+      import :: c_int32_t, C_CHAR
+      CHARACTER(KIND=C_CHAR), DIMENSION(*), intent(in) :: ifile, ofile
+    END function copyfile
+  END INTERFACE
+
+  include 'netcdf.inc'
+
+  type(domain), intent(inout)             :: grid
+  type(grid_config_rec_type),intent(inout) :: config_flags
+
+! Declare local parameters
+  character(len=120) :: file_name,filnam
+  character(len=19)  :: DateStr1
+  character(len=4)   :: staggering=' N/A'
+  character(len=3)   :: ordering
+  character(len=80), dimension(3)  ::  dimnames
+  character(len=80)  :: rmse_var
+  integer            :: dh1,dh0
+  integer            :: i,j,k
+  integer            :: ndim1
+  integer            :: WrfType
+  integer            :: it, ierr, Status, Status_next_time
+  integer            :: wrf_real
+  integer            :: nlon_regional,nlat_regional,nsig_regional
+  integer            :: ids, ide, jds, jde, kds, kde,    &
+                        ims, ime, jms, jme, kms, kme,    &
+                        ips, ipe, jps, jpe, kps, kpe    
+  integer, dimension(4)           :: start_index, end_index1
+  real, dimension(:), allocatable :: globbuf
+  real*4,allocatable :: field3(:,:,:),field2(:,:)
+  real*4,allocatable :: field3u(:,:,:),field3v(:,:,:),field3ph(:,:,:)
+  character(len=4) ::  fgname
+  integer :: julyr, julday
+  real    :: gmt
+  real*4  :: gmt4
+  real    :: qvf
+
+      if ( grid%auxhist5_oid .NE. 0 ) then
+         call close_dataset ( grid%auxhist5_oid , config_flags , "DATASET=AUXHIST5" )
+      endif
+
+      call open_w_dataset (grid%auxhist5_oid, trim(config_flags%auxhist5_outname), grid, config_flags, &
+                            output_auxhist5, "DATASET=AUXHIST5", ierr)
+      if ( ierr .NE. 0 ) CALL wrf_error_fatal('Error opening '//trim(filnam))
+
+      start_date=current_date
+
+      call geth_julgmt(julyr, julday, gmt)
+      config_flags%gmt = gmt
+      config_flags%julyr = julyr
+      config_flags%julday = julday
+
+      call output_auxhist5 (grid%auxhist5_oid, grid , config_flags , ierr)
+      if ( ierr .NE. 0 ) CALL wrf_error_fatal('Error writing Gradient in auxhist5')
+      call close_dataset (grid%auxhist5_oid, config_flags, "DATASET=AUXHIST5")
+
+end subroutine da_write_anaincrements
diff --git a/var/da/da_transfer_model/da_transfer_xatowrf.inc b/var/da/da_transfer_model/da_transfer_xatowrf.inc
index 1c8b64856e..85328935dd 100644
--- a/var/da/da_transfer_model/da_transfer_xatowrf.inc
+++ b/var/da/da_transfer_model/da_transfer_xatowrf.inc
@@ -265,6 +265,7 @@ subroutine da_transfer_xatowrf(grid, config_flags)
             ! The analysis perturbation = Hydro_ph - base_ph + nonhydro_xb_ph:
             grid%ph_2(i,j,k+1) = ph_full - grid%phb(i,j,k+1) &
                             + (grid%xb%hf(i,j,k+1)*gravity - ph_xb_hd)
+            grid%t_iau(i,j,k) = grid%xa%t(i,j,k)
          end do
       end do
    end do
@@ -303,6 +304,7 @@ subroutine da_transfer_xatowrf(grid, config_flags)
           phm = mu_full*c3h(k)   + c4h(k)   + grid%p_top
           ph(k+1) = ph(k) + ald(k)*phm*LOG(pfd/pfu) 
           grid%ph_2(i,j,k+1) = ph(k+1) - grid%phb(i,j,k+1)
+          grid%t_iau(i,j,k) = grid%xa%t(i,j,k)
        END DO
 
       ! Update geopotential perturbation
@@ -323,6 +325,7 @@ subroutine da_transfer_xatowrf(grid, config_flags)
      do j=jts,jte
         do i=its,ite
            ph_cgrid(i,j,k) = grid%ph_2(i,j,k) - ph_cgrid(i,j,k)
+           grid%ph_iau(i,j,k) = ph_cgrid(i,j,k)
         end do
      end do
    end do
@@ -429,6 +432,7 @@ subroutine da_transfer_xatowrf(grid, config_flags)
          grid%mu_2(i,j) = grid%mu_2(i,j) + mu_cgrid(i,j)
          grid%w_2(i,j,kte+1)=  grid%w_2(i,j,kte+1) + grid%xa%w(i,j,kte+1)
          grid%psfc(i,j) = grid%psfc(i,j) + grid%xa%psfc(i,j)
+         grid%mu_iau(i,j) = mu_cgrid(i,j)
       end do
 
       do k=kts,kte
@@ -437,6 +441,9 @@ subroutine da_transfer_xatowrf(grid, config_flags)
             grid%u_2(i,j,k) = grid%u_2(i,j,k) + u_cgrid(i,j,k)
             grid%v_2(i,j,k) = grid%v_2(i,j,k) + v_cgrid(i,j,k)
 #endif
+            grid%u_iau(i,j,k) = u_cgrid(i,j,k)
+            grid%v_iau(i,j,k) = v_cgrid(i,j,k)
+            grid%w_iau(i,j,k) = grid%xa%w(i,j,k)
             grid%w_2(i,j,k) = grid%w_2(i,j,k) + grid%xa%w(i,j,k)
 
             ! (xb%q+xa%q in specific humidity) >= 0.0
@@ -448,19 +455,25 @@ subroutine da_transfer_xatowrf(grid, config_flags)
             else
                grid%moist(i,j,k,P_QV) = grid%moist(i,j,k,P_QV)+q_cgrid(i,j,k)
             end if
+            grid%qv_iau(i,j,k) = q_cgrid(i,j,k)
 
             if (size(grid%moist,dim=4) >= 4 .and. cloud_cv_options >= 1) then
                grid%moist(i,j,k,p_qc) = max(grid%moist(i,j,k,p_qc) + grid%xa%qcw(i,j,k), 0.0)
                grid%moist(i,j,k,p_qr) = max(grid%moist(i,j,k,p_qr) + grid%xa%qrn(i,j,k), 0.0)
+               grid%qc_iau(i,j,k) = grid%xa%qcw(i,j,k)
+               grid%qr_iau(i,j,k) = grid%xa%qrn(i,j,k)
             end if
 
             if (size(grid%moist,dim=4) >= 6 .and. cloud_cv_options >= 2) then
                grid%moist(i,j,k,p_qi) = max(grid%moist(i,j,k,p_qi) + grid%xa%qci(i,j,k), 0.0)
                grid%moist(i,j,k,p_qs) = max(grid%moist(i,j,k,p_qs) + grid%xa%qsn(i,j,k), 0.0)
+               grid%qi_iau(i,j,k) = grid%xa%qci(i,j,k)
+               grid%qs_iau(i,j,k) = grid%xa%qsn(i,j,k)
             end if
 
             if (size(grid%moist,dim=4) >= 7 .and. cloud_cv_options >= 2) then
                grid%moist(i,j,k,p_qg) = max(grid%moist(i,j,k,p_qg) + grid%xa%qgr(i,j,k), 0.0)
+               grid%qg_iau(i,j,k) = grid%xa%qgr(i,j,k)
             end if
          end do
       end do