diff --git a/src/main.f90 b/src/main.f90
index 4f36c1e..000c1d6 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -375,7 +375,6 @@ program cans
       call updatep(pp,p)
       call boundp(cbcpre,n,bcpre,nb,is_bound,dl,dzc,p)
     end if
-    dto = dt
     !
     ! check simulation stopping criteria
     !
@@ -389,6 +388,7 @@ program cans
       tw = (MPI_WTIME()-twi)/3600.
       if(tw    >= tw_max  ) is_done = is_done.or..true.
     end if
+    dto = dt
     if(mod(istep,icheck) == 0) then
       if(myid == 0) print*, 'Calculating maximum velocity to set ACDI gamma parameter...'
       call acdi_set_gamma(n,acdi_gam_factor,u,v,w,gam)
@@ -403,7 +403,6 @@ program cans
         is_done = .true.
         kill = .true.
       end if
-      dto = dt
       dti = 1./dt
       call chkdiv(lo,hi,dli,dzfi,u,v,w,divtot,divmax)
       if(myid == 0) print*, 'Total divergence = ', divtot, '| Maximum divergence = ', divmax
diff --git a/src/mom.f90 b/src/mom.f90
index 2542861..59f3a62 100644
--- a/src/mom.f90
+++ b/src/mom.f90
@@ -826,9 +826,9 @@ subroutine mom_xyz_oth(nx,ny,nz,dxi,dyi,dzci,rho12,beta12,bforce,gacc,sigma,rho0
           factorxp = rhobeta + drhobeta*psixp
           factoryp = rhobeta + drhobeta*psiyp
           factorzp = rhobeta + drhobeta*psizp
-          dudt_aux = dudt_aux - gaccx*factorxp*0.5*(s_pcc+s_ccc)
-          dvdt_aux = dvdt_aux - gaccy*factoryp*0.5*(s_cpc+s_ccc)
-          dwdt_aux = dwdt_aux - gaccz*factorzp*0.5*(s_ccp+s_ccc)
+          dudt_aux = dudt_aux - gaccx*factorxp*0.5*(s_pcc+s_ccc)/rhoxp
+          dvdt_aux = dvdt_aux - gaccy*factoryp*0.5*(s_cpc+s_ccc)/rhoyp
+          dwdt_aux = dwdt_aux - gaccz*factorzp*0.5*(s_ccp+s_ccc)/rhozp
 #endif
           dudt(i,j,k) = dudt_aux
           dvdt(i,j,k) = dvdt_aux
@@ -1101,9 +1101,9 @@ subroutine mom_xyz_all(nx,ny,nz,dxi,dyi,dzci,dzfi,rho12,mu12,beta12,bforce,gacc,
           factorxp = rhobeta + drhobeta*psixp
           factoryp = rhobeta + drhobeta*psiyp
           factorzp = rhobeta + drhobeta*psizp
-          dudt_aux = dudt_aux - gaccx*factorxp*0.5*(s_pcc+s_ccc)
-          dvdt_aux = dvdt_aux - gaccy*factoryp*0.5*(s_cpc+s_ccc)
-          dwdt_aux = dwdt_aux - gaccz*factorzp*0.5*(s_ccp+s_ccc)
+          dudt_aux = dudt_aux - gaccx*factorxp*0.5*(s_pcc+s_ccc)/rhoxp
+          dvdt_aux = dvdt_aux - gaccy*factoryp*0.5*(s_cpc+s_ccc)/rhoyp
+          dwdt_aux = dwdt_aux - gaccz*factorzp*0.5*(s_ccp+s_ccc)/rhozp
 #endif
           dudt(i,j,k) = dudt_aux
           dvdt(i,j,k) = dvdt_aux
diff --git a/src/rk.f90 b/src/rk.f90
index 9a2d1b8..baa2aba 100644
--- a/src/rk.f90
+++ b/src/rk.f90
@@ -51,11 +51,6 @@ subroutine rk(rkpar,n,dli,dzci,dzfi,dt, &
       allocate(dudtrko_t(n(1),n(2),n(3)),dvdtrko_t(n(1),n(2),n(3)),dwdtrko_t(n(1),n(2),n(3)))
       !$acc enter data create(dudtrk_t ,dvdtrk_t ,dwdtrk_t ) async(1)
       !$acc enter data create(dudtrko_t,dvdtrko_t,dwdtrko_t) async(1)
-      !$acc kernels default(present) async(1) ! not really necessary
-      dudtrko_t(:,:,:) = dudtrk_t(:,:,:)
-      dvdtrko_t(:,:,:) = dvdtrk_t(:,:,:)
-      dwdtrko_t(:,:,:) = dwdtrk_t(:,:,:)
-      !$acc end kernels
       dudtrk  => dudtrk_t
       dvdtrk  => dvdtrk_t
       dwdtrk  => dwdtrk_t
@@ -146,16 +141,18 @@ subroutine rk_scal(rkpar,n,dli,dzci,dzfi,dt, &
     factor12 = factor1 + factor2
     rhocp12 = rho12(:)*cp12(:)
     if(is_first) then ! leverage save attribute to allocate these arrays on the device only once
-      is_first = .false.
       allocate(dsdtrk_t(n(1),n(2),n(3)),dsdtrko_t(n(1),n(2),n(3)))
       !$acc enter data create(dsdtrk_t,dsdtrko_t) async(1)
-      !$acc kernels default(present) async(1) ! not really necessary
-      dsdtrko_t(:,:,:) = 0._rp
-      !$acc end kernels
       dsdtrk  => dsdtrk_t
       dsdtrko => dsdtrko_t
     end if
     call scal_ad(n(1),n(2),n(3),dli(1),dli(2),dzci,dzfi,ssource,ka12,rhocp12,psi,u,v,w,s,dsdtrk)
+    if(is_first) then ! use Euler forward
+      !$acc kernels
+      dsdtrko(:,:,:) = dsdtrk(:,:,:)
+      !$acc end kernels
+      is_first = .false.
+    end if
     !$acc parallel loop collapse(3) default(present) async(1)
     do k=1,n(3)
       do j=1,n(2)
@@ -195,16 +192,18 @@ subroutine rk_2fl(rkpar,n,dli,dzci,dzfi,dt,gam,seps,u,v,w,psi)
     factor2 = rkpar(2)*dt
     factor12 = factor1 + factor2
     if(is_first) then ! leverage save attribute to allocate these arrays on the device only once
-      is_first = .false.
       allocate(dpsidtrk_t(n(1),n(2),n(3)),dpsidtrko_t(n(1),n(2),n(3)))
       !$acc enter data create(dpsidtrk_t,dpsidtrko_t) async(1)
-      !$acc kernels default(present) async(1) ! not really necessary
-      dpsidtrko_t(:,:,:) = 0._rp
-      !$acc end kernels
       dpsidtrk  => dpsidtrk_t
       dpsidtrko => dpsidtrko_t
     end if
     call acdi_transport_pf(n(1),n(2),n(3),dli(1),dli(2),dli(3),dzci,dzfi,gam,seps,u,v,w,psi,dpsidtrk)
+    if(is_first) then ! use Euler forward
+      !$acc kernels
+      dpsidtrko(:,:,:) = dpsidtrk(:,:,:)
+      !$acc end kernels
+      is_first = .false.
+    end if
     !$acc parallel loop collapse(3) default(present) async(1)
     !$OMP PARALLEL DO   COLLAPSE(3) DEFAULT(shared)
     do k=1,n(3)