Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 4697)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 4742)
@@ -24,4 +24,7 @@
 ! -----------------
 ! $Id$
+! Implement snow and graupel (bulk microphysics)
+!
+! 4697 2020-09-25 08:24:29Z suehring
 ! To avoid numerical oscillations which may lead to a built-up of passive scalar near non-cyclic
 ! boundaries, employ a first-order scheme at the 3 grid points next to the boundary for passive
@@ -252,4 +255,8 @@
                sums_wsqis_ws_l,                                                                    &
                sums_wsnis_ws_l,                                                                    &
+               sums_wsqgs_ws_l,                                                                    &
+               sums_wsngs_ws_l,                                                                    &
+               sums_wsqss_ws_l,                                                                    &
+               sums_wsnss_ws_l,                                                                    &
                sums_wssas_ws_l,                                                                    &
                sums_wsus_ws_l,                                                                     &
@@ -1828,4 +1835,14 @@
           ENDDO
 
+       CASE ( 'qg' )
+
+          DO  k = nzb, nzt
+             sums_wsqgs_ws_l(k,tn)  = sums_wsqgs_ws_l(k,tn) +                                      &
+                                      ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
+                                                  * ( w(k,j,i) - hom(k,1,3,0)                 )    &
+                                      + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )    &
+                                                  *   ABS( w(k,j,i) - hom(k,1,3,0)            )    &
+                                      ) * weight_substep(intermediate_timestep_count)
+          ENDDO
 
        CASE ( 'qi' )
@@ -1851,4 +1868,15 @@
           ENDDO
 
+       CASE ( 'qs' )
+
+          DO  k = nzb, nzt
+             sums_wsqss_ws_l(k,tn)  = sums_wsqss_ws_l(k,tn) +                                      &
+                                      ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
+                                                  * ( w(k,j,i) - hom(k,1,3,0)                 )    &
+                                      + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )    &
+                                                  *   ABS( w(k,j,i) - hom(k,1,3,0)            )    &
+                                      ) * weight_substep(intermediate_timestep_count)
+          ENDDO
+
        CASE ( 'nc' )
 
@@ -1862,4 +1890,16 @@
           ENDDO
 
+       CASE ( 'ng' )
+
+          DO  k = nzb, nzt
+             sums_wsngs_ws_l(k,tn)  = sums_wsngs_ws_l(k,tn) +                                      &
+                                      ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
+                                                  * ( w(k,j,i) - hom(k,1,3,0)                 )    &
+                                      + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )    &
+                                                  *   ABS( w(k,j,i) - hom(k,1,3,0)            )    &
+                                      ) * weight_substep(intermediate_timestep_count)
+          ENDDO
+
+
        CASE ( 'ni' )
 
@@ -1877,4 +1917,15 @@
           DO  k = nzb, nzt
              sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +                                      &
+                                      ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
+                                                  * ( w(k,j,i) - hom(k,1,3,0)                 )    &
+                                      + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )    &
+                                                  *   ABS( w(k,j,i) - hom(k,1,3,0)            )    &
+                                      ) * weight_substep(intermediate_timestep_count)
+          ENDDO
+
+       CASE ( 'ns' )
+
+          DO  k = nzb, nzt
+             sums_wsnss_ws_l(k,tn)  = sums_wsnss_ws_l(k,tn) +                                      &
                                       ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )    &
                                                   * ( w(k,j,i) - hom(k,1,3,0)                 )    &
@@ -3450,4 +3501,12 @@
         CASE ( 'qi' )
            sk_num = 11
+        CASE ( 'ng' )
+           sk_num = 12
+        CASE ( 'qg' )
+           sk_num = 13
+        CASE ( 'ns' )
+           sk_num = 14
+        CASE ( 'qs' )
+           sk_num = 15
         CASE DEFAULT
            sk_num = 0
@@ -3477,4 +3536,6 @@
     !$ACC PRESENT(sums_wsnrs_ws_l, sums_wsss_ws_l) &
     !$ACC PRESENT(sums_wsnis_ws_l, sums_wsqis_ws_l) &
+    !$ACC PRESENT(sums_wsngs_ws_l, sums_wsqgs_ws_l) &
+    !$ACC PRESENT(sums_wsnss_ws_l, sums_wsqss_ws_l) &
     !$ACC PRESENT(sums_salsa_ws_l)
     DO  i = nxl, nxr
@@ -4123,4 +4184,46 @@
                                                 *   ABS(w(k,j,i) - hom(k,1,3,0)             )      &
                                             ) * weight_substep(intermediate_timestep_count)
+                CASE ( 12 )
+                   !$ACC ATOMIC
+                   sums_wsngs_ws_l(k,tn)  = sums_wsngs_ws_l(k,tn) +                                &
+                                            ( flux_t(k)                                            &
+                                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )      &
+                                                * ( w(k,j,i) - hom(k,1,3,0)                 )      &
+                                            + diss_t(k)                                            &
+                                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )      &
+                                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )      &
+                                            ) * weight_substep(intermediate_timestep_count)
+                CASE ( 13 )
+                   !$ACC ATOMIC
+                   sums_wsqgs_ws_l(k,tn)  = sums_wsqgs_ws_l(k,tn) +                                &
+                                            ( flux_t(k)                                            &
+                                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )      &
+                                                * ( w(k,j,i) - hom(k,1,3,0)                 )      &
+                                            + diss_t(k)                                            &
+                                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )      &
+                                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )      &
+                                            ) * weight_substep(intermediate_timestep_count)
+
+                CASE ( 14 )
+                   !$ACC ATOMIC
+                   sums_wsnss_ws_l(k,tn)  = sums_wsnss_ws_l(k,tn) +                                &
+                                            ( flux_t(k)                                            &
+                                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )      &
+                                                * ( w(k,j,i) - hom(k,1,3,0)                 )      &
+                                            + diss_t(k)                                            &
+                                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )      &
+                                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )      &
+                                            ) * weight_substep(intermediate_timestep_count)
+                CASE ( 15 )
+                   !$ACC ATOMIC
+                   sums_wsqss_ws_l(k,tn)  = sums_wsqss_ws_l(k,tn) +                                &
+                                            ( flux_t(k)                                            &
+                                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )      &
+                                                * ( w(k,j,i) - hom(k,1,3,0)                 )      &
+                                            + diss_t(k)                                            &
+                                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )      &
+                                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )      &
+                                            ) * weight_substep(intermediate_timestep_count)
+
 
              END SELECT
