Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 2233)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 2292)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$ 
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2233 2017-05-30 18:08:54Z suehring
 ! 
 ! 2232 2017-05-30 17:47:52Z suehring
@@ -239,12 +244,14 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  diss_l_e, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr,         &
-                  diss_l_s, diss_l_sa, diss_l_u, diss_l_v, diss_l_w,  flux_l_e,&
-                  flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_s,         &
-                  flux_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_s_e, diss_s_nr,&
-                  diss_s_pt, diss_s_q, diss_s_qr, diss_s_s, diss_s_sa,         &
-                  diss_s_u, diss_s_v, diss_s_w, flux_s_e, flux_s_nr, flux_s_pt,&
-                  flux_s_q, flux_s_qr, flux_s_s, flux_s_sa, flux_s_u, flux_s_v,&
-                  flux_s_w
+           ONLY:  diss_l_e, diss_l_nc, diss_l_nr, diss_l_pt, diss_l_q,         &
+                  diss_l_qc, diss_l_qr, diss_l_s, diss_l_sa, diss_l_u,         &
+                  diss_l_v, diss_l_w, flux_l_e, flux_l_nc, flux_l_nr,          & 
+                  flux_l_pt, flux_l_q, flux_l_qc, flux_l_qr, flux_l_s,         &
+                  flux_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_s_e,           &
+                  diss_s_nc,  diss_s_nr, diss_s_pt, diss_s_q, diss_s_qc,       &
+                  diss_s_qr, diss_s_s, diss_s_sa, diss_s_u, diss_s_v,          &
+                  diss_s_w, flux_s_e, flux_s_nc, flux_s_nr, flux_s_pt,         &
+                  flux_s_q, flux_s_qc, flux_s_qr, flux_s_s, flux_s_sa,         &
+                  flux_s_u, flux_s_v, flux_s_w
 
        USE constants,                                                          &
@@ -254,5 +261,5 @@
        USE control_parameters,                                                 &
            ONLY:  cloud_physics, humidity, loop_optimization,                  &
-                  passive_scalar, microphysics_seifert,                        &
+                  passive_scalar, microphysics_morrison, microphysics_seifert, &
                   ocean, ws_scheme_mom, ws_scheme_sca
 
@@ -265,8 +272,10 @@
 
        USE statistics,                                                         &
-           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,&
-                  sums_wspts_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l,            &
-                  sums_wsss_ws_l, sums_wssas_ws_l,  sums_wsss_ws_l,            &
-                  sums_wsus_ws_l, sums_wsvs_ws_l  
+           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsncs_ws_l,& 
+                  sums_wsnrs_ws_l,sums_wspts_ws_l, sums_wsqcs_ws_l,            &
+                  sums_wsqrs_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l,             &
+                  sums_wssas_ws_l,  sums_wsss_ws_l, sums_wsus_ws_l,            &
+                  sums_wsvs_ws_l
+  
 
 !
@@ -309,4 +318,11 @@
              ALLOCATE( sums_wsss_ws_l(nzb:nzt+1,0:threads_per_task-1) )
              sums_wsss_ws_l = 0.0_wp
+          ENDIF
+
+          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+             ALLOCATE( sums_wsqcs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             ALLOCATE( sums_wsncs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             sums_wsqcs_ws_l = 0.0_wp
+             sums_wsncs_ws_l = 0.0_wp
           ENDIF
 
@@ -372,5 +388,16 @@
                 ALLOCATE( flux_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
                           diss_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
-             ENDIF             
+             ENDIF
+
+             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+                ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1),           &
+                          diss_s_qc(nzb+1:nzt,0:threads_per_task-1),           &
+                          flux_s_nc(nzb+1:nzt,0:threads_per_task-1),           &
+                          diss_s_nc(nzb+1:nzt,0:threads_per_task-1) )
+                ALLOCATE( flux_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                          diss_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                          flux_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
+                          diss_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 
+             ENDIF                  
 
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
@@ -1013,13 +1040,15 @@
        USE control_parameters,                                                 &
            ONLY:  cloud_physics, humidity, passive_scalar, ocean,              &
-                  microphysics_seifert, ws_scheme_mom, ws_scheme_sca
+                  microphysics_morrison, microphysics_seifert, ws_scheme_mom,  & 
+                  ws_scheme_sca
 
        USE kinds
 
        USE statistics,                                                         &
-           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,&
-                  sums_wspts_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l,            &
-                  sums_wsss_ws_l, sums_wssas_ws_l,  sums_wsus_ws_l,            &
-                  sums_wsvs_ws_l  
+           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsncs_ws_l,&
+                  sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqcs_ws_l,           &
+                  sums_wsqrs_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l,             &
+                  sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l      
+                   
 
        IMPLICIT NONE
@@ -1040,4 +1069,8 @@
           IF ( humidity       )  sums_wsqs_ws_l = 0.0_wp
           IF ( passive_scalar )  sums_wsss_ws_l = 0.0_wp
+          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+             sums_wsqcs_ws_l = 0.0_wp
+             sums_wsncs_ws_l = 0.0_wp
+          ENDIF
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              sums_wsqrs_ws_l = 0.0_wp
@@ -1081,7 +1114,8 @@
 
        USE statistics,                                                         &
-           ONLY:  hom, sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqrs_ws_l,      &
-                  sums_wsqs_ws_l, sums_wssas_ws_l, sums_wsss_ws_l,             &
-                  weight_substep
+           ONLY:  hom, sums_wsncs_ws_l, sums_wsnrs_ws_l, sums_wspts_ws_l,      &
+                  sums_wsqcs_ws_l,  sums_wsqrs_ws_l, sums_wsqs_ws_l,           &
+                  sums_wssas_ws_l, sums_wsss_ws_l, weight_substep    
+                  
 
        IMPLICIT NONE
@@ -1553,8 +1587,31 @@
              ENDDO
 
+          CASE ( 'qc' )
+
+             DO  k = nzb, nzt
+                sums_wsqcs_ws_l(k,tn)  = sums_wsqcs_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 ( 'qr' )
 
              DO  k = nzb, nzt
                 sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_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' )
+
+             DO  k = nzb, nzt
+                sums_wsncs_ws_l(k,tn)  = sums_wsncs_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)                 )  &
@@ -3089,6 +3146,8 @@
        USE statistics,                                                        &
            ONLY:  hom, sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,      &
-                  sums_wsqrs_ws_l, sums_wsnrs_ws_l, sums_wsss_ws_l,           &
-                  weight_substep
+                  sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsncs_ws_l,          &
+                  sums_wsnrs_ws_l, sums_wsss_ws_l, weight_substep  
+                  
+
 
        IMPLICIT NONE
@@ -3550,7 +3609,29 @@
                             ) * weight_substep(intermediate_timestep_count)
                     ENDDO
+                 CASE ( 'qc' )
+                    DO  k = nzb, nzt
+                       sums_wsqcs_ws_l(k,tn)  = sums_wsqcs_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 ( 'qr' )
                     DO  k = nzb, nzt
                        sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_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' )
+                    DO  k = nzb, nzt
+                       sums_wsncs_ws_l(k,tn)  = sums_wsncs_ws_l(k,tn)          &
                           + ( flux_t(k)                                        &
                                 / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
