diff --git a/README.md b/README.md
index 26145eb..bd87dc8 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-TEMPO-v3.1.0
+TEMPO-v3.1.1
 
 ## Thompson-Eidhammer Microphysics Parameterization for Operations
 
diff --git a/src/module_mp_tempo_driver.F90 b/src/module_mp_tempo_driver.F90
index 0b5eb1d..b54563b 100644
--- a/src/module_mp_tempo_driver.F90
+++ b/src/module_mp_tempo_driver.F90
@@ -273,7 +273,7 @@ subroutine tempo_run(tempo_cfgs, dt, itimestep, &
     real(wp), dimension(:), allocatable :: thten_swrad1d
 
     integer :: i, j, k, nz
-    logical :: use_temperature 
+    logical :: use_temperature, is_first_step
 
     type(ty_tempo_main_diags) :: tempo_main_diags
     type(ty_tempo_driver_diags), intent(inout) :: tempo_diags
@@ -394,6 +394,8 @@ subroutine tempo_run(tempo_cfgs, dt, itimestep, &
     !   error stop "tempo_run() --- requires either temperature or theta and Exner function"
     endif 
 
+    is_first_step = (itimestep == 1)
+
     ! tempo driver code
     do j = jts, jte
       do i = its, ite
@@ -437,7 +439,7 @@ subroutine tempo_run(tempo_cfgs, dt, itimestep, &
         if (present(land_input)) land1d = land_input(i,j)
 
         ! main call to the 1d tempo microphysics
-        call tempo_main(tempo_cfgs=tempo_cfgs, &
+        call tempo_main(tempo_cfgs=tempo_cfgs, is_first_step=is_first_step, &
           qv1d=qv1d, qc1d=qc1d, qi1d=qi1d, qr1d=qr1d, qs1d=qs1d, qg1d=qg1d, qb1d=qb1d, &
           ni1d=ni1d, nr1d=nr1d, nc1d=nc1d, ng1d=ng1d, nwfa1d=nwfa1d, nifa1d=nifa1d, t1d=t1d, p1d=p1d, &
           w1d=w1d, dz1d=dz1d, land1d=land1d, &
diff --git a/src/module_mp_tempo_main.F90 b/src/module_mp_tempo_main.F90
index b8de91f..8c03185 100644
--- a/src/module_mp_tempo_main.F90
+++ b/src/module_mp_tempo_main.F90
@@ -70,7 +70,7 @@ module module_mp_tempo_main
 
   contains
 
-  subroutine tempo_main(tempo_cfgs, &
+  subroutine tempo_main(tempo_cfgs, is_first_step, &
     qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, ni1d, nr1d, nc1d, ng1d, &
     nwfa1d, nifa1d, t1d, p1d, w1d, dz1d, land1d, &
     qcfrac1d, qifrac1d, qc_bl1d, qcfrac_bl1d, &
@@ -85,6 +85,7 @@ subroutine tempo_main(tempo_cfgs, &
     type(ty_tempo_cfgs), intent(in) :: tempo_cfgs
     integer, intent(in) :: kts, kte, ii, jj
     real(wp), intent(in) :: dt
+    logical, intent(in), optional :: is_first_step
     real(wp), dimension(kts:kte), intent(inout) :: t1d !! 1D temperature \([K]\)
     real(wp), dimension(kts:kte), intent(in) :: p1d !! 1D pressure \([Pa]\)
     real(wp), dimension(kts:kte), intent(inout) :: qv1d !! 1D water vapor mixing ratio \([kg\; kg^{-1}]\)
@@ -148,11 +149,13 @@ subroutine tempo_main(tempo_cfgs, &
     
     ! local variables
     real(wp) :: tempc, tc0, odt
-    logical :: do_micro, supersaturated
-    logical, save :: first_call_main = .true.
+    logical :: do_micro, supersaturated, local_first_step
     integer :: k, nz
 
     ! --------------------------------------------------------------------------------------------
+    local_first_step = .false.
+    if (present(is_first_step)) local_first_step = is_first_step
+
     do_micro = .false.
     supersaturated = .false.
     odt = 1._wp / dt
@@ -323,7 +326,7 @@ subroutine tempo_main(tempo_cfgs, &
     enddo
 
     if (present(nwfa1d)) then
-      if (first_call_main) then
+      if (local_first_step) then
         if (sum(nwfa1d) < eps) call init_water_friendly_aerosols(dz1d, nwfa)
       endif
     else
@@ -331,7 +334,7 @@ subroutine tempo_main(tempo_cfgs, &
     endif
     
     if (present(nifa1d)) then
-      if (first_call_main) then
+      if (local_first_step) then
         if (sum(nifa1d) < eps) call init_ice_friendly_aerosols(dz1d, nifa)
       endif
     else
@@ -381,7 +384,7 @@ subroutine tempo_main(tempo_cfgs, &
       dt=dt, odt=odt)
 
     ! init ng and qb
-    if (first_call_main) then
+    if (local_first_step) then
       if (present(ng1d) .and. present(qb1d)) then
         if (sum(qg1d) > r1 .and. sum(ng1d) < eps .and. sum(qb1d) < eps) then
           call graupel_init(rho, qg1d, ng1d, qb1d)
@@ -410,8 +413,6 @@ subroutine tempo_main(tempo_cfgs, &
       satw, sati, ssatw, ssati, diffu, visco, vsc2, ocp, lvap, tcond, lvt2, &
       supersaturated)
 
-    if (first_call_main) first_call_main = .false.
-
     ! check for hydrometeors or supersaturation --------------------------------------------------
     do_micro = any(l_qc) .or. any(l_qr) .or. any(l_qi) .or. any(l_qs) .or. any(l_qg)
     if (.not. do_micro .and. .not. supersaturated) return