Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 3298)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 3302)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Stokes drift velocity added in scalar advection
+! 
+! 3298 2018-10-02 12:21:11Z kanani
 ! Add formatted note from ketelsen about statistics for chemistry
 !
@@ -1147,5 +1150,6 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air_zw
+           ONLY:  ddzw, drho_air, rho_air_zw, tend, u, u_stokes_zu, v,         &
+                  v_stokes_zu, w
 
        USE control_parameters,                                                 &
@@ -1227,5 +1231,5 @@
              ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
 
-             v_comp                  = v(k,j,i) - v_gtrans
+             v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
              swap_flux_y_local(k,tn) = v_comp *         (                     &
                                                ( 37.0_wp * ibit5 * adv_sca_5  &
@@ -1263,5 +1267,5 @@
           DO  k = nzb_max+1, nzt
 
-             v_comp                  = v(k,j,i) - v_gtrans
+             v_comp                  = v(k,j,i) - v_gtrans + v_stokes_zu(k)
              swap_flux_y_local(k,tn) = v_comp * (                             &
                                     37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )   &
@@ -1288,5 +1292,5 @@
              ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
 
-             u_comp                     = u(k,j,i) - u_gtrans
+             u_comp                     = u(k,j,i) - u_gtrans + u_stokes_zu(k)
              swap_flux_x_local(k,j,tn) = u_comp * (                           &
                                                ( 37.0_wp * ibit2 * adv_sca_5  &
@@ -1323,5 +1327,5 @@
           DO  k = nzb_max+1, nzt
 
-             u_comp                 = u(k,j,i) - u_gtrans
+             u_comp                 = u(k,j,i) - u_gtrans + u_stokes_zu(k)
              swap_flux_x_local(k,j,tn) = u_comp * (                           &
                                       37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) ) &
@@ -1355,5 +1359,5 @@
           ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
 
-          u_comp    = u(k,j,i+1) - u_gtrans
+          u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
           flux_r(k) = u_comp * (                                              &
                      ( 37.0_wp * ibit2 * adv_sca_5                            &
@@ -1390,5 +1394,5 @@
           ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
 
-          v_comp    = v(k,j+1,i) - v_gtrans
+          v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
           flux_n(k) = v_comp * (                                              &
                      ( 37.0_wp * ibit5 * adv_sca_5                            &
@@ -1514,5 +1518,5 @@
        DO  k = nzb_max+1, nzt
 
-          u_comp    = u(k,j,i+1) - u_gtrans
+          u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
           flux_r(k) = u_comp * (                                              &
                       37.0_wp * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
@@ -1524,5 +1528,5 @@
                     +           ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
 
-          v_comp    = v(k,j+1,i) - v_gtrans
+          v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
           flux_n(k) = v_comp * (                                              &
                       37.0_wp * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
@@ -3192,22 +3196,23 @@
     SUBROUTINE advec_s_ws( sk, sk_char )
 
-       USE arrays_3d,                                                         &
-           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air_zw
-
-       USE control_parameters,                                                &
+       USE arrays_3d,                                                          &
+           ONLY:  ddzw, drho_air, rho_air_zw, tend, u, u_stokes_zu, v,         &
+                  v_stokes_zu, w
+
+       USE control_parameters,                                                 &
            ONLY:  intermediate_timestep_count, u_gtrans, v_gtrans 
 
-       USE grid_variables,                                                    &
+       USE grid_variables,                                                     &
            ONLY:  ddx, ddy
 
-       USE indices,                                                           &
-           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max,   &
+       USE indices,                                                            &
+           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max,    &
                   nzt, advc_flags_1
            
        USE kinds
        
-       USE statistics,                                                        &
-           ONLY:  hom, sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,      &
-                  sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsncs_ws_l,          &
+       USE statistics,                                                         &
+           ONLY:  hom, sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,       &
+                  sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsncs_ws_l,           &
                   sums_wsnrs_ws_l, sums_wsss_ws_l, weight_substep  
                   
@@ -3272,5 +3277,5 @@
              ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)
 
-             u_comp                 = u(k,j,i) - u_gtrans
+             u_comp                 = u(k,j,i) - u_gtrans + u_stokes_zu(k)
              swap_flux_x_local(k,j) = u_comp * (                              &
                                              ( 37.0_wp * ibit2 * adv_sca_5    &
@@ -3307,5 +3312,5 @@
           DO  k = nzb_max+1, nzt
 
-             u_comp                 = u(k,j,i) - u_gtrans
+             u_comp                 = u(k,j,i) - u_gtrans + u_stokes_zu(k)
              swap_flux_x_local(k,j) = u_comp * (                              &
                                       37.0_wp * ( sk(k,j,i)   + sk(k,j,i-1) ) &
@@ -3333,5 +3338,5 @@
              ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)
 
-             v_comp               = v(k,j,i) - v_gtrans
+             v_comp               = v(k,j,i) - v_gtrans + v_stokes_zu(k)
              swap_flux_y_local(k) = v_comp * (                                &
                                              ( 37.0_wp * ibit5 * adv_sca_5    &
@@ -3369,5 +3374,5 @@
           DO  k = nzb_max+1, nzt
 
-             v_comp               = v(k,j,i) - v_gtrans
+             v_comp               = v(k,j,i) - v_gtrans + v_stokes_zu(k)
              swap_flux_y_local(k) = v_comp * (                               &
                                     37.0_wp * ( sk(k,j,i)   + sk(k,j-1,i) )  &
@@ -3375,5 +3380,5 @@
                                   +           ( sk(k,j+2,i) + sk(k,j-3,i) )  &
                                              ) * adv_sca_5
-              swap_diss_y_local(k) = -ABS( v_comp ) * (                      &
+             swap_diss_y_local(k) = -ABS( v_comp ) * (                       &
                                     10.0_wp * ( sk(k,j,i)   - sk(k,j-1,i) )  &
                                   -  5.0_wp * ( sk(k,j+1,i) - sk(k,j-2,i) )  &
@@ -3396,5 +3401,5 @@
                 ibit0 = IBITS(advc_flags_1(k,j,i),0,1)
 
-                u_comp    = u(k,j,i+1) - u_gtrans
+                u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
                 flux_r(k) = u_comp * (                                        &
                           ( 37.0_wp * ibit2 * adv_sca_5                       &
@@ -3431,5 +3436,5 @@
                 ibit3 = IBITS(advc_flags_1(k,j,i),3,1)
 
-                v_comp    = v(k,j+1,i) - v_gtrans
+                v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
                 flux_n(k) = v_comp * (                                        &
                           ( 37.0_wp * ibit5 * adv_sca_5                       &
@@ -3548,5 +3553,5 @@
              DO  k = nzb_max+1, nzt
 
-                u_comp    = u(k,j,i+1) - u_gtrans
+                u_comp    = u(k,j,i+1) - u_gtrans + u_stokes_zu(k)
                 flux_r(k) = u_comp * (                                        &
                       37.0_wp * ( sk(k,j,i+1) + sk(k,j,i)   )                 &
@@ -3558,5 +3563,5 @@
                     +           ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
 
-                v_comp    = v(k,j+1,i) - v_gtrans
+                v_comp    = v(k,j+1,i) - v_gtrans + v_stokes_zu(k)
                 flux_n(k) = v_comp * (                                        &
                       37.0_wp * ( sk(k,j+1,i) + sk(k,j,i)   )                 &
