Changeset 2292 for palm/trunk

Timestamp:

Jun 20, 2017 9:51:42 AM (7 years ago)

Author:

schwenkel

Message:

implementation of new bulk microphysics scheme

Location:

palm/trunk/SOURCE

Files:

• advec_s_bc.f90 (modified) (3 diffs)
• advec_ws.f90 (modified) (12 diffs)
• average_3d_data.f90 (modified) (2 diffs)
• boundary_conds.f90 (modified) (8 diffs)
• calc_liquid_water_content.f90 (modified) (4 diffs)
• check_parameters.f90 (modified) (5 diffs)
• data_output_2d.f90 (modified) (3 diffs)
• data_output_3d.f90 (modified) (3 diffs)
• data_output_mask.f90 (modified) (3 diffs)
• flow_statistics.f90 (modified) (5 diffs)
• init_3d_model.f90 (modified) (12 diffs)
• init_masks.f90 (modified) (3 diffs)
• microphysics_mod.f90 (modified) (51 diffs)
• modules.f90 (modified) (16 diffs)
• netcdf_interface_mod.f90 (modified) (7 diffs)
• palm.f90 (modified) (3 diffs)
• parin.f90 (modified) (4 diffs)
• pmc_interface_mod.f90 (modified) (24 diffs)
• prognostic_equations.f90 (modified) (5 diffs)
• read_3d_binary.f90 (modified) (3 diffs)
• sum_up_3d_data.f90 (modified) (4 diffs)
• surface_layer_fluxes_mod.f90 (modified) (7 diffs)
• surface_mod.f90 (modified) (40 diffs)
• swap_timelevel.f90 (modified) (6 diffs)
• time_integration.f90 (modified) (5 diffs)
• write_3d_binary.f90 (modified) (4 diffs)

palm/trunk/SOURCE/advec_s_bc.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2101 2017-01-05 16:42:31Z suehring
 !
 ! 2000 2016-08-20 18:09:15Z knoop
@@ -1012 +1017 @@
           ENDIF
 
-       ELSEIF ( sk_char == 'qr' )  THEN
-
-!
-!--       Rain water content boundary condition at the bottom boundary:
+       ELSEIF ( sk_char == 'qc' )  THEN
+
+!
+!--       Cloud water content boundary condition at the bottom boundary:
 !--       Dirichlet (fixed surface rain water content).
           DO  i = nxl, nxr
@@ -1026 +1031 @@
 
 !
+!--       Cloud water content boundary condition at the top boundary: Dirichlet
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
+                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
+             ENDDO
+          ENDDO
+
+       ELSEIF ( sk_char == 'qr' )  THEN
+
+!
+!--       Rain water content boundary condition at the bottom boundary:
+!--       Dirichlet (fixed surface rain water content).
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
+                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
+             ENDDO
+          ENDDO
+
+!
 !--       Rain water content boundary condition at the top boundary: Dirichlet
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
+                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
+             ENDDO
+          ENDDO
+
+       ELSEIF ( sk_char == 'nc' )  THEN
+
+!
+!--       Cloud drop concentration boundary condition at the bottom boundary:
+!--       Dirichlet (fixed surface cloud drop concentration).
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
+                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
+             ENDDO
+          ENDDO
+
+!
+!--       Cloud drop concentration boundary condition at the top boundary: Dirichlet
           DO  i = nxl, nxr
              DO  j = nys, nyn

palm/trunk/SOURCE/advec_ws.f90

@@ -25 +25 @@
 ! -----------------
 ! $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 +244 @@
 
        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 +261 @@
        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 +272 @@
 
        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 +318 @@
              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 +388 @@
                 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 +1040 @@
        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 +1069 @@
           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 +1114 @@
 
        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 +1587 @@
              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 +3146 @@
        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 +3609 @@
                             ) * 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) ) )  &

palm/trunk/SOURCE/average_3d_data.f90

@@ -25 +25 @@
 ! -----------------
 ! $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
@@ -194 +199 @@
              ENDDO
 
+          CASE ( 'nc' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      nc_av(k,j,i) = nc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                   ENDDO
+                ENDDO
+             ENDDO
+
           CASE ( 'nr' )
              DO  i = nxlg, nxrg

palm/trunk/SOURCE/boundary_conds.f90

@@ -25 +25 @@
 ! -----------------
 ! $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
@@ -151 +156 @@
     USE arrays_3d,                                                             &
         ONLY:  c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l,  &
-               dzu, e_p, nr_p, pt, pt_p, q, q_p, qr_p, s, s_p, sa, sa_p,       &
-               u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p,                 &
+               dzu, e_p, nc_p, nr_p, pt, pt_p, q, q_p, qc_p, qr_p, s, s_p, sa, & 
+               sa_p, u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p,           &
                v, vg, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p,                 &
                w, w_p, w_m_l, w_m_n, w_m_r, w_m_s, pt_init
-
+               
     USE control_parameters,                                                    &
         ONLY:  bc_pt_t_val, bc_q_t_val, bc_s_t_val, constant_diffusion,        &
@@ -162 +167 @@
                ibc_sa_t, ibc_uv_b, ibc_uv_t, inflow_l, inflow_n, inflow_r,     &
                inflow_s, intermediate_timestep_count, large_scale_forcing,     &
-               microphysics_seifert, nest_domain, nest_bound_l, nest_bound_s,  &
-               nudging, ocean, outflow_l, outflow_n, outflow_r, outflow_s,     &
-               passive_scalar, tsc, use_cmax
+               microphysics_morrison, microphysics_seifert, nest_domain,       &
+               nest_bound_l, nest_bound_s, nudging, ocean, outflow_l,          &
+               outflow_n, outflow_r, outflow_s, passive_scalar, tsc, use_cmax
 
     USE grid_variables,                                                        &
@@ -393 +398 @@
        ELSEIF ( ibc_q_t == 1 ) THEN
           q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
+       ENDIF
+
+       IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+!             
+!--       Surface conditions cloud water (Dirichlet)
+!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
+!--       the k coordinate belongs to the atmospheric grid point, therefore, set
+!--       qr_p and nr_p at k-1
+          !$OMP PARALLEL DO PRIVATE( i, j, k )
+          DO  m = 1, bc_h(0)%ns
+             i = bc_h(0)%i(m)            
+             j = bc_h(0)%j(m)
+             k = bc_h(0)%k(m)
+             qc_p(k-1,j,i) = 0.0_wp
+             nc_p(k-1,j,i) = 0.0_wp
+          ENDDO
+!
+!--       Top boundary condition for cloud water (Dirichlet)
+          qc_p(nzt+1,:,:) = 0.0_wp
+          nc_p(nzt+1,:,:) = 0.0_wp
+           
        ENDIF
 
@@ -514 +540 @@
        IF ( humidity )  THEN
           q_p(:,nys-1,:) = q_p(:,nys,:)
+          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+             qc_p(:,nys-1,:) = qc_p(:,nys,:)
+             nc_p(:,nys-1,:) = nc_p(:,nys,:)
+          ENDIF
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              qr_p(:,nys-1,:) = qr_p(:,nys,:)
@@ -525 +555 @@
        IF ( humidity )  THEN
           q_p(:,nyn+1,:) = q_p(:,nyn,:)
+          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+             qc_p(:,nyn+1,:) = qc_p(:,nyn,:)
+             nc_p(:,nyn+1,:) = nc_p(:,nyn,:)
+          ENDIF
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
@@ -536 +570 @@
        IF ( humidity )  THEN
           q_p(:,:,nxl-1) = q_p(:,:,nxl)
+          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+             qc_p(:,:,nxl-1) = qc_p(:,:,nxl)
+             nc_p(:,:,nxl-1) = nc_p(:,:,nxl)
+          ENDIF
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
@@ -547 +585 @@
        IF ( humidity )  THEN
           q_p(:,:,nxr+1) = q_p(:,:,nxr)
+          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+             qc_p(:,:,nxr+1) = qc_p(:,:,nxr)
+             nc_p(:,:,nxr+1) = nc_p(:,:,nxr)
+          ENDIF
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              qr_p(:,:,nxr+1) = qr_p(:,:,nxr)

palm/trunk/SOURCE/calc_liquid_water_content.f90

@@ -25 +25 @@
 ! -----------------
 ! $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
@@ -83 +88 @@
 
     USE control_parameters,                                                    &
-        ONLY:  microphysics_seifert
+        ONLY:  microphysics_morrison, microphysics_seifert
 
     USE indices,                                                               &
@@ -132 +137 @@
 !
 !--          Compute the liquid water content
-             IF ( microphysics_seifert )  THEN
-                IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp ) THEN
+             IF ( microphysics_seifert  .AND.  .NOT. microphysics_morrison )   &
+             THEN
+                IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp )  THEN
                    qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) )                  &
                                       * MERGE( 1.0_wp, 0.0_wp,                 &
@@ -147 +153 @@
                                                BTEST( wall_flags_0(k,j,i), 0 ) )
                 ENDIF
+             ELSEIF ( microphysics_morrison )  THEN
+                ql(k,j,i) = qc(k,j,i) + qr(k,j,i)                              &
+                                      * MERGE( 1.0_wp, 0.0_wp,                 &
+                                               BTEST( wall_flags_0(k,j,i), 0 ) )
              ELSE
-                IF ( ( q(k,j,i) - q_s ) > 0.0_wp ) THEN
+                IF ( ( q(k,j,i) - q_s ) > 0.0_wp )  THEN
                    qc(k,j,i) = ( q(k,j,i) - q_s )                              &
                                       * MERGE( 1.0_wp, 0.0_wp,                 &

palm/trunk/SOURCE/check_parameters.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2274 2017-06-09 13:27:48Z Giersch
 ! Changed error messages 
 ! 
@@ -983 +988 @@
        microphysics_seifert    = .TRUE.
        microphysics_kessler    = .FALSE.
+       microphysics_morrison  = .FALSE.
        precipitation           = .TRUE.
     ELSEIF ( cloud_scheme == 'kessler' )  THEN
@@ -988 +994 @@
        microphysics_seifert    = .FALSE.
        microphysics_kessler    = .TRUE.
+       microphysics_morrison   = .FALSE.
+       precipitation           = .TRUE.
+    ELSEIF ( cloud_scheme == 'morrison' )  THEN
+       microphysics_sat_adjust = .FALSE.
+       microphysics_seifert    = .TRUE.
+       microphysics_kessler    = .FALSE.
+       microphysics_morrison   = .TRUE.
        precipitation           = .TRUE.
     ELSE
@@ -2766 +2779 @@
              hom(:,2,122,:) = SPREAD( zw, 2, statistic_regions+1 )
 
+          CASE ( 'nc' )
+             IF (  .NOT.  cloud_physics )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_physics = .FALSE.'
+                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
+             ELSEIF ( .NOT.  microphysics_morrison )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_scheme /= morrison'
+                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 89
+                dopr_unit(i)  = '1/m3'
+                hom(:,2,89,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
           CASE ( 'nr' )
              IF (  .NOT.  cloud_physics )  THEN
@@ -3061 +3091 @@
              ENDIF
              unit = 'K'
+
+          CASE ( 'nc' )
+             IF (  .NOT.  cloud_physics )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                         'res cloud_physics = .TRUE.'
+                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
+             ELSEIF ( .NOT.  microphysics_morrison )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                         'res = microphysics morrison '
+                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = '1/m3'
 
           CASE ( 'nr' )

palm/trunk/SOURCE/data_output_2d.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2277 2017-06-12 10:47:51Z kanani
 ! Removed unused variables do2d_xy_n, do2d_xz_n, do2d_yz_n
 ! 
@@ -169 +174 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  dzw, e, nr, p, pt, precipitation_amount, precipitation_rate,    &
+        ONLY:  dzw, e, nc, nr, p, pt, precipitation_amount, precipitation_rate,&
                prr, q, qc, ql, ql_c, ql_v, ql_vp, qr, rho_ocean, s, sa,        &
                tend, u, v, vpt, w, zu, zw
@@ -492 +497 @@
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'nc_xy', 'nc_xz', 'nc_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => nc
+                ELSE
+                   to_be_resorted => nc_av
+                ENDIF
+                IF ( mode == 'xy' )  level_z = zu
 
              CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )

palm/trunk/SOURCE/data_output_3d.f90

@@ -25 +25 @@
 ! -----------------
 ! $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
@@ -155 +160 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa, &
-               tend, u, v, vpt, w
+        ONLY:  e, nc, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, &
+               sa, tend, u, v, vpt, w
         
     USE averaging
@@ -325 +330 @@
              ELSE
                 to_be_resorted => lpt_av
+             ENDIF
+
+          CASE ( 'nc' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => nc
+             ELSE
+                to_be_resorted => nc_av
              ENDIF
 

palm/trunk/SOURCE/data_output_mask.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2101 2017-01-05 16:42:31Z suehring
 !
 ! 2031 2016-10-21 15:11:58Z knoop
@@ -105 +110 @@
 #if defined( __netcdf )
     USE arrays_3d,                                                             &
-        ONLY:  e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa, tend, u,   &
-               v, vpt, w
+        ONLY:  e, nc, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa,  &
+               tend, u, v, vpt, w
     
     USE averaging,                                                             &
-        ONLY:  e_av, lpt_av, nr_av, p_av, pc_av, pr_av, pt_av, q_av, qc_av,    &
-               ql_av, ql_c_av, ql_v_av, ql_vp_av, qv_av, qr_av, rho_ocean_av, s_av,  &
-               sa_av, u_av, v_av, vpt_av, w_av 
+        ONLY:  e_av, lpt_av, nc_av, nr_av, p_av, pc_av, pr_av, pt_av, q_av,    &
+               qc_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qv_av, qr_av,         &
+               rho_ocean_av, s_av, sa_av, u_av, v_av, vpt_av, w_av 
     
     USE cloud_parameters,                                                      &
@@ -243 +248 @@
              ELSE
                 to_be_resorted => lpt_av
+             ENDIF
+
+          CASE ( 'nc' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => nc
+             ELSE
+                to_be_resorted => nc_av
              ENDIF
 

palm/trunk/SOURCE/flow_statistics.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2270 2017-06-09 12:18:47Z maronga
 ! Revised numbering (removed 2 timeseries)
 ! 
@@ -233 +238 @@
     USE arrays_3d,                                                             &
         ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
-               momentumflux_output_conversion, nr, p, prho, prr, pt, q,        &
+               momentumflux_output_conversion, nc, nr, p, prho, prr, pt, q,    &
                qc, ql, qr, rho_air, rho_air_zw, rho_ocean, s,                  &
                sa, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, u,   &
@@ -245 +250 @@
                 dt_3d, g, humidity, kappa, land_surface, large_scale_forcing,  &
                 large_scale_subsidence, max_pr_user, message_string, neutral,  &
-                microphysics_seifert, ocean, passive_scalar, simulated_time,   &
-                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
-                ws_scheme_mom, ws_scheme_sca
+                microphysics_morrison, microphysics_seifert, ocean,            &
+                passive_scalar, simulated_time, use_subsidence_tendencies,     &
+                use_surface_fluxes, use_top_fluxes, ws_scheme_mom,             &
+                ws_scheme_sca
         
     USE cpulog,                                                                &
@@ -1217 +1223 @@
                                                              rmask(j,i,sr) *   &
                                                              flag
+                         IF ( microphysics_morrison )  THEN
+                            sums_l(k,123,tn) = sums_l(k,123,tn) + nc(k,j,i) *  &
+                                                                rmask(j,i,sr) *&
+                                                                flag
+                         ENDIF
                          IF ( microphysics_seifert )  THEN
                             sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
@@ -1759 +1770 @@
        hom(:,1,71,sr) = sums(:,71)     ! prho
        hom(:,1,72,sr) = hyp * 1E-2_wp  ! hyp in hPa
+       hom(:,1,123,sr) = sums(:,123)   ! nc
        hom(:,1,73,sr) = sums(:,73)     ! nr
        hom(:,1,74,sr) = sums(:,74)     ! qr

palm/trunk/SOURCE/init_3d_model.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2277 2017-06-12 10:47:51Z kanani
 ! Removed unused variable sums_up_fraction_l
 ! 
@@ -560 +565 @@
 
           IF ( cloud_physics )  THEN
-
 !
 !--          Liquid water content
@@ -568 +572 @@
              ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 #endif
+
+!
+!--          3D-cloud water content
+             IF ( .NOT. microphysics_morrison )  THEN
+#if defined( __nopointer )
+                ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+#else
+                ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+#endif
+             ENDIF
 !
 !--          Precipitation amount and rate (only needed if output is switched)
@@ -574 +588 @@
 
 !
-!--          3D-cloud water content
-#if defined( __nopointer )
-             ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-#else
-             ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-#endif
-!
 !--          3d-precipitation rate
              ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+             IF ( microphysics_morrison )  THEN
+!
+!--             3D-cloud drop water content, cloud drop concentration arrays
+#if defined( __nopointer )
+                ALLOCATE( nc(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
+                          nc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
+                          qc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          tnc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 
+                          tqc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+#else
+                ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
+                          qc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+#endif
+             ENDIF
 
              IF ( microphysics_seifert )  THEN
@@ -880 +907 @@
           IF ( cloud_physics )  THEN
              ql => ql_1
-             qc => qc_1
+             IF ( .NOT. microphysics_morrison )  THEN
+                qc => qc_1
+             ENDIF
+             IF ( microphysics_morrison )  THEN
+                qc => qc_1;  qc_p  => qc_2;  tqc_m  => qc_3
+                nc => nc_1;  nc_p  => nc_2;  tnc_m  => nc_3
+             ENDIF
              IF ( microphysics_seifert )  THEN
                 qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
@@ -979 +1012 @@
                 ENDDO
              ENDDO
+             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qc(:,j,i) = 0.0_wp
+                      nc(:,j,i) = 0.0_wp
+                   ENDDO
+                ENDDO
+             ENDIF
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
                 DO  i = nxlg, nxrg
@@ -986 +1027 @@
                    ENDDO
                 ENDDO
-
              ENDIF
           ENDIF
+
           IF ( passive_scalar )  THEN
              DO  i = nxlg, nxrg
@@ -1107 +1148 @@
                 ENDDO
              ENDDO
+             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qc(:,j,i) = 0.0_wp
+                      nc(:,j,i) = 0.0_wp
+                   ENDDO
+                ENDDO
+             ENDIF
+
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -1115 +1164 @@
                    ENDDO
                 ENDDO
-
              ENDIF
+
           ENDIF
           
@@ -1319 +1368 @@
           tq_m = 0.0_wp
           q_p = q
+          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+             tqc_m = 0.0_wp
+             qc_p  = qc
+             tnc_m = 0.0_wp
+             nc_p  = nc
+          ENDIF
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              tqr_m = 0.0_wp
@@ -1532 +1587 @@
        IF ( humidity )  THEN
           q_p = q
+          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+             qc_p = qc
+             nc_p = nc
+          ENDIF
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              qr_p = qr
@@ -1547 +1606 @@
        IF ( humidity )  THEN
           tq_m = 0.0_wp
+          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+             tqc_m = 0.0_wp
+             tnc_m = 0.0_wp
+          ENDIF
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
              tqr_m = 0.0_wp

palm/trunk/SOURCE/init_masks.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2271 2017-06-09 12:34:55Z sward
 ! Changed error message
 ! 
@@ -115 +120 @@
                mask_x_loop, mask_xyz_dimension, mask_y, mask_y_loop, mask_z,   &
                mask_z_loop, max_masks,  message_string, mid,                   &
-               microphysics_seifert, passive_scalar, ocean
+               microphysics_morrison, microphysics_seifert, passive_scalar,    &
+               ocean
 
     USE grid_variables,                                                        &
@@ -278 +284 @@
                 unit = 'K'
 
+             CASE ( 'nc' )
+                IF ( .NOT. cloud_physics )  THEN
+                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
+                        '" requires cloud_physics = .TRUE.'
+                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
+                 ELSEIF ( .NOT. microphysics_morrison ) THEN
+                   message_string = 'output of "' // TRIM( var ) // '" requi' //  &
+                         'res  = morrison'
+                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
+                ENDIF
+                unit = '1/m3'
+
              CASE ( 'nr' )
                 IF ( .NOT. cloud_physics )  THEN

palm/trunk/SOURCE/microphysics_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! - The process of activation is parameterized with a simple Twomey 
+!   activion scheme or with considering solution and curvature 
+!   effects (Khvorostyanov and Curry ,2006).
+! - The saturation adjustment scheme is replaced by the parameterization
+!   of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
+! - All other microphysical processes of Seifert and Beheng are used.
+!   Additionally, in those processes the reduction of cloud number concentration 
+!   is considered. 
+! 
+! 2233 2017-05-30 18:08:54Z suehring
 !
 ! 2232 2017-05-30 17:47:52Z suehring
@@ -123 +136 @@
     IMPLICIT NONE
 
-    LOGICAL ::  cloud_water_sedimentation = .FALSE.  !< cloud water sedimentation
-    LOGICAL ::  limiter_sedimentation = .TRUE.       !< sedimentation limiter
-    LOGICAL ::  collision_turbulence = .FALSE.       !< turbulence effects
-    LOGICAL ::  ventilation_effect = .TRUE.          !< ventilation effect
+    LOGICAL ::  cloud_water_sedimentation = .FALSE.       !< cloud water sedimentation
+    LOGICAL ::  curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
+    LOGICAL ::  limiter_sedimentation = .TRUE.            !< sedimentation limiter
+    LOGICAL ::  collision_turbulence = .FALSE.            !< turbulence effects
+    LOGICAL ::  ventilation_effect = .TRUE.               !< ventilation effect
 
     REAL(wp) ::  a_1 = 8.69E-4_wp          !< coef. in turb. parametrization (cm-2 s3)
@@ -147 +161 @@
     REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
     REAL(wp) ::  eps_sb = 1.0E-10_wp       !< threshold in two-moments scheme
+    REAL(wp) ::  eps_mr = 0.0_wp           !< threshold for morrison scheme
     REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
     REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
@@ -161 +176 @@
     REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
     REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
+    REAL(wp) ::  xamin = 5.24E-19_wp       !< average aerosol mass (kg) (~ 0.05µm)
+    REAL(wp) ::  xcmin = 4.18E-15_wp       !< minimum cloud drop size (kg) (~ 1µm)
+    REAL(wp) ::  xcmax = 2.6E-10_wp        !< maximum cloud drop size (kg) (~ 40µm)
     REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
     REAL(wp) ::  xrmax = 5.0E-6_wp         !< maximum rain drop site (kg)
@@ -167 +185 @@
     REAL(wp) ::  dpirho_l                  !< 6.0 / ( pi * rho_l )
     REAL(wp) ::  dt_micro                  !< microphysics time step
+    REAL(wp) ::  na_init = 100.0E6_wp      !< Total particle/aerosol concentration (cm-3)
     REAL(wp) ::  nc_const = 70.0E6_wp      !< cloud droplet concentration
     REAL(wp) ::  dt_precipitation = 100.0_wp !< timestep precipitation (s)
@@ -184 +203 @@
     PUBLIC microphysics_control, microphysics_init
 
-    PUBLIC cloud_water_sedimentation, collision_turbulence, c_sedimentation,   &
-           dt_precipitation, limiter_sedimentation, nc_const, sigma_gc,        &
+    PUBLIC cloud_water_sedimentation, collision_turbulence,                    &
+           curvature_solution_effects_bulk, c_sedimentation, dt_precipitation, &
+           limiter_sedimentation, na_init, nc_const, sigma_gc,                 &
            ventilation_effect
+           
 
     INTERFACE microphysics_control
@@ -198 +219 @@
     END INTERFACE adjust_cloud
 
+    INTERFACE activation
+       MODULE PROCEDURE activation
+       MODULE PROCEDURE activation_ij
+    END INTERFACE activation
+
+    INTERFACE condensation
+       MODULE PROCEDURE condensation
+       MODULE PROCEDURE condensation_ij
+    END INTERFACE condensation
+
     INTERFACE autoconversion
        MODULE PROCEDURE autoconversion
@@ -256 +287 @@
 
        USE control_parameters,                                                 &
-           ONLY:  microphysics_seifert
+           ONLY:  microphysics_morrison, microphysics_seifert
 
        USE indices,                                                            &
@@ -283 +314 @@
 !
 !--    Allocate 1D microphysics arrays
-       ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1),         &
+       ALLOCATE ( pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1),         &
                   qc_1d(nzb:nzt+1) )
+
+       IF ( microphysics_morrison .OR. microphysics_seifert ) THEN
+          ALLOCATE ( nc_1d(nzb:nzt+1) )
+       ENDIF
 
        IF ( microphysics_seifert )  THEN
@@ -290 +325 @@
        ENDIF
 
-!
-!--    Initialze nc_1d with nc_const
-       nc_1d = nc_const
-
     END SUBROUTINE microphysics_init
 
@@ -313 +344 @@
            ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
                   intermediate_timestep_count, large_scale_forcing,            &
-                  lsf_surf, microphysics_kessler, microphysics_seifert,        &
-                  pt_surface, rho_surface,surface_pressure
+                  lsf_surf, microphysics_kessler, microphysics_morrison,       &
+                  microphysics_seifert, pt_surface,                            &
+                  rho_surface,surface_pressure
 
        USE indices,                                                            &
@@ -373 +405 @@
 
           CALL adjust_cloud
+          IF ( microphysics_morrison )  CALL activation
+          IF ( microphysics_morrison )  CALL condensation
           CALL autoconversion
           CALL accretion
@@ -396 +430 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  qr, nr
+           ONLY:  qc, nc, qr, nr
 
        USE cloud_parameters,                                                   &
            ONLY:  hyrho
+
+       USE control_parameters,                                                 &
+           ONLY:  microphysics_morrison
 
        USE cpulog,                                                             &
@@ -434 +471 @@
                    ENDIF
                 ENDIF
+                IF ( microphysics_morrison ) THEN
+                   IF ( qc(k,j,i) <= eps_sb )  THEN
+                      qc(k,j,i) = 0.0_wp
+                      nc(k,j,i) = 0.0_wp
+                   ELSE
+                      IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
+                         nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin *            &
+                                          MERGE( 1.0_wp, 0.0_wp,               &           
+                                              BTEST( wall_flags_0(k,j,i), 0 ) )
+                      ENDIF
+                   ENDIF
+                ENDIF
              ENDDO
           ENDDO
@@ -442 +491 @@
     END SUBROUTINE adjust_cloud
 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate number of activated condensation nucleii after simple activation
+!> scheme of Twomey, 1959.
+!------------------------------------------------------------------------------!
+    SUBROUTINE activation
+
+       USE arrays_3d,                                                          &
+           ONLY:  hyp, nc, nr, pt, q,  qc, qr
+
+       USE cloud_parameters,                                                   &
+           ONLY:  hyrho, l_d_cp, l_d_r, l_v, rho_l, r_v, t_d_pt
+
+       USE constants,                                                          &
+           ONLY:  pi
+
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
+       USE indices,                                                            &
+           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
+
+       USE kinds
+
+       USE control_parameters,                                                 &
+          ONLY: simulated_time
+
+       USE particle_attributes,                                                &
+          ONLY:  molecular_weight_of_solute, molecular_weight_of_water, rho_s, &
+              s1, s2, s3, vanthoff
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+
+       REAL(wp)     ::  activ             !<
+       REAL(wp)     ::  afactor           !<
+       REAL(wp)     ::  alpha             !<
+       REAL(wp)     ::  beta_act          !<
+       REAL(wp)     ::  bfactor           !<
+       REAL(wp)     ::  e_s               !<
+       REAL(wp)     ::  k_act             !<
+       REAL(wp)     ::  n_act             !<
+       REAL(wp)     ::  n_ccn             !<
+       REAL(wp)     ::  q_s               !<
+       REAL(wp)     ::  rd0               !<
+       REAL(wp)     ::  s_0               !<
+       REAL(wp)     ::  sat               !<
+       REAL(wp)     ::  sat_max           !<
+       REAL(wp)     ::  sigma             !<
+       REAL(wp)     ::  t_int             !<
+       REAL(wp)     ::  t_l               !<
+
+       CALL cpu_log( log_point_s(65), 'activation', 'start' )
+
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             DO  k = nzb+1, nzt
+
+!
+!--             Actual liquid water temperature:
+                t_l = t_d_pt(k) * pt(k,j,i)
+
+!
+!--             Calculate actual temperature
+                t_int = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
+!
+!--             Saturation vapor pressure at t_l:
+                e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
+                                       ( t_l - 35.86_wp )                   &
+                                     )
+!
+!--             Computation of saturation humidity:
+                q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
+                alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
+                q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /               &
+                        ( 1.0_wp + alpha * q_s )
+
+!--             Supersaturation:
+                sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
+
+!
+!--             Prescribe parameters for activation
+!--             (see: Bott + Trautmann, 2002, Atm. Res., 64)
+                k_act  = 0.7_wp 
+
+                IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
+!
+!--                Compute the number of activated Aerosols 
+!--                (see: Twomey, 1959, Pure and applied Geophysics, 43)
+                   n_act     = na_init * sat**k_act
+!
+!--                Compute the number of cloud droplets
+!--                (see: Morrison + Grabowski, 2007, JAS, 64)
+!                  activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro
+
+!
+!--                Compute activation rate after Khairoutdinov and Kogan
+!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)   
+                   sat_max = 1.0_wp / 100.0_wp                  
+                   activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN       & 
+                      ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /    & 
+                       dt_micro
+                ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk ) THEN
+!
+!--                Curvature effect (afactor) with surface tension 
+!--                parameterization by Straka (2009)
+                   sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
+                   afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
+!
+!--                Solute effect (bfactor)
+                   bfactor = vanthoff * molecular_weight_of_water *            &
+                       rho_s / ( molecular_weight_of_solute * rho_l )
+
+!
+!--                Prescribe power index that describes the soluble fraction 
+!--                of an aerosol particle (beta) and mean geometric radius of
+!--                dry aerosol spectrum
+!--                (see: Morrison + Grabowski, 2007, JAS, 64)
+                   beta_act = 0.5_wp
+                   rd0      = 0.05E-6_wp
+!
+!--                Calculate mean geometric supersaturation (s_0) with 
+!--                parameterization by Khvorostyanov and Curry (2006)
+                   s_0   = rd0 **(- ( 1.0_wp + beta_act ) ) *                  &
+                       ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
+
+!
+!--                Calculate number of activated CCN as a function of 
+!--                supersaturation and taking Koehler theory into account
+!--                (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)       
+                   n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(              & 
+                      LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(s1) ) ) )      
+                   activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
+                ENDIF
+
+                nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+       CALL cpu_log( log_point_s(65), 'activation', 'stop' )
+
+    END SUBROUTINE activation
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate condensation rate for cloud water content (after Khairoutdinov and  
+!> Kogan, 2000).
+!------------------------------------------------------------------------------!
+    SUBROUTINE condensation
+
+       USE arrays_3d,                                                          &
+           ONLY:  hyp, nr, pt, q,  qc, qr, nc
+
+       USE cloud_parameters,                                                   &
+           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
+
+       USE constants,                                                          &
+           ONLY:  pi
+
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
+       USE indices,                                                            &
+           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
+
+       USE kinds
+
+       USE control_parameters,                                                 &
+          ONLY: simulated_time
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+
+       REAL(wp)     ::  alpha             !<
+       REAL(wp)     ::  cond              !<
+       REAL(wp)     ::  cond_max          !<
+       REAL(wp)     ::  dc                !<
+       REAL(wp)     ::  e_s               !<
+       REAL(wp)     ::  evap              !<
+       REAL(wp)     ::  evap_nc           !<
+       REAL(wp)     ::  g_fac             !<
+       REAL(wp)     ::  nc_0              !<
+       REAL(wp)     ::  q_s               !<
+       REAL(wp)     ::  sat               !<
+       REAL(wp)     ::  t_l               !<
+       REAL(wp)     ::  temp              !<
+       REAL(wp)     ::  xc                !<
+
+       CALL cpu_log( log_point_s(66), 'condensation', 'start' )
+
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             DO  k = nzb+1, nzt
+!
+!--             Actual liquid water temperature:
+                t_l = t_d_pt(k) * pt(k,j,i)
+!
+!--             Saturation vapor pressure at t_l:
+                e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /       &
+                                       ( t_l - 35.86_wp )                      &
+                                     )
+!
+!--             Computation of saturation humidity:
+                q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
+                alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
+                q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /                  &
+                        ( 1.0_wp + alpha * q_s )
+
+!--             Supersaturation:
+                sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
+
+!
+!--             Actual temperature:
+                temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
+
+                g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
+                                    l_v / ( thermal_conductivity_l * temp )    &
+                                    + r_v * temp / ( diff_coeff_l * e_s )      &
+                                    )
+!
+!--             Mean weight of cloud drops
+                IF ( nc(k,j,i) <= 0.0_wp) CYCLE
+                xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)                
+!
+!--             Weight averaged diameter of cloud drops:
+                dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
+!
+!--             Integral diameter of cloud drops
+                nc_0 = nc(k,j,i) * dc                
+!
+!--             Condensation needs only to be calculated in supersaturated regions
+                IF ( sat > 0.0_wp )  THEN
+!
+!--                Condensation rate of cloud water content
+!--                after KK scheme.
+!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
+                   cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
+                   cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
+                   cond      = MIN( cond, cond_max / dt_micro )
+               
+                   qc(k,j,i) = qc(k,j,i) + cond * dt_micro
+                ELSEIF ( sat < 0.0_wp ) THEN
+                   evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
+                   evap      = MAX( evap, -qc(k,j,i) / dt_micro )
+
+                   qc(k,j,i) = qc(k,j,i) + evap * dt_micro
+                ENDIF
+                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
+                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+
+       CALL cpu_log( log_point_s(66), 'condensation', 'stop' )
+
+    END SUBROUTINE condensation
+
 
 !------------------------------------------------------------------------------!
@@ -451 +769 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  diss, dzu, nr, qc, qr
+           ONLY:  diss, dzu, nc, nr, qc, qr
 
        USE cloud_parameters,                                                   &
@@ -457 +775 @@
 
        USE control_parameters,                                                 &
-           ONLY:  rho_surface
+           ONLY:  microphysics_morrison, rho_surface
 
        USE cpulog,                                                             &
@@ -482 +800 @@
        REAL(wp)     ::  k_au              !<
        REAL(wp)     ::  l_mix             !<
+       REAL(wp)     ::  nc_auto           !<
        REAL(wp)     ::  nu_c              !<
        REAL(wp)     ::  phi_au            !<
@@ -500 +819 @@
                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 
+                nc_auto = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
+
                 IF ( qc(k,j,i) > eps_sb )  THEN
 
@@ -518 +839 @@
 !
 !--                Mean weight of cloud droplets:
-                   xc = hyrho(k) * qc(k,j,i) / nc_const
+                   xc = hyrho(k) * qc(k,j,i) / nc_auto
 !
 !--                Parameterized turbulence effects on autoconversion (Seifert,
@@ -569 +890 @@
                    nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro  &
                                                               * flag
+                   IF ( microphysics_morrison ) THEN
+                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *         & 
+                                    autocon / x0 * hyrho(k) * dt_micro * flag )
+                   ENDIF
 
                 ENDIF
@@ -654 +979 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  diss, qc, qr
+           ONLY:  diss, qc, qr, nc
 
        USE cloud_parameters,                                                   &
@@ -660 +985 @@
 
        USE control_parameters,                                                 &
-           ONLY:  rho_surface
+           ONLY:  microphysics_morrison, rho_surface
 
        USE cpulog,                                                             &
@@ -679 +1004 @@
        REAL(wp)     ::  flag              !< flag to mask topography grid points
        REAL(wp)     ::  k_cr              !<
+       REAL(wp)     ::  nc_accr           !<
        REAL(wp)     ::  phi_ac            !<
        REAL(wp)     ::  tau_cloud         !<
+       REAL(wp)     ::  xc                !<
+
 
        CALL cpu_log( log_point_s(56), 'accretion', 'start' )
@@ -691 +1019 @@
                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 
-                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb ) )  THEN
+                nc_accr = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
+
+                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )    &
+                                             .AND.  ( nc_accr > eps_mr ) ) THEN
 !
 !--                Intern time scale of coagulation (Seifert and Beheng, 2006):
@@ -699 +1030 @@
 !--                Beheng, 2001):
                    phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
+
+!
+!--                Mean weight of cloud drops
+                   xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)      
 !
 !--                Parameterized turbulence effects on autoconversion (Seifert,
@@ -719 +1054 @@
                    qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
                    qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
+                   IF ( microphysics_morrison )  THEN
+                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i),                  &
+                                         accr / xc * hyrho(k) * dt_micro * flag)
+                   ENDIF
 
                 ENDIF
@@ -976 +1315 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  ddzu, dzu, pt, prr, q, qc
+           ONLY:  ddzu, dzu, nc, pt, prr, q, qc
 
        USE cloud_parameters,                                                   &
@@ -982 +1321 @@
 
        USE control_parameters,                                                 &
-           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
+           ONLY:  call_microphysics_at_all_substeps,                           & 
+                  intermediate_timestep_count, microphysics_morrison
 
        USE cpulog,                                                             &
@@ -1002 +1342 @@
        INTEGER(iwp) ::  k             !<
 
-       REAL(wp)                       ::  flag   !< flag to mask topography grid points
+       REAL(wp) ::  flag              !< flag to mask topography grid points
+       REAL(wp) ::  nc_sedi           !<
+
        REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc !<
+
 
        CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
@@ -1015 +1359 @@
 !--             Predetermine flag to mask topography
                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+                nc_sedi = MERGE ( nc(k,j,i), nc_const, microphysics_morrison )
+
+!
+!--             Sedimentation fluxes for number concentration are only calculated
+!--             for cloud_scheme = 'morrison'
+                IF ( microphysics_morrison ) THEN
+                   IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
+                      sed_nc(k) = sed_qc_const *                               & 
+                              ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
+                              ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
+                   ELSE
+                      sed_nc(k) = 0.0_wp
+                   ENDIF
+
+                   sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *           &
+                                    nc(k,j,i) / dt_micro + sed_nc(k+1)         &
+                                  ) * flag
+
+                   nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *       &
+                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
+                ENDIF
 
                 IF ( qc(k,j,i) > eps_sb )  THEN
-                   sed_qc(k) = sed_qc_const * nc_const**( -2.0_wp / 3.0_wp ) * &
+                   sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *  &
                                ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
                                                                            flag
@@ -1389 +1754 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  hyp, nr, pt, pt_init, prr, q, qc, qr, zu
+           ONLY:  hyp, nc, nr, pt, pt_init, prr, q, qc, qr, zu
 
        USE cloud_parameters,                                                   &
@@ -1397 +1762 @@
            ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
                   intermediate_timestep_count, large_scale_forcing,            &
-                  lsf_surf, microphysics_seifert, microphysics_kessler,        &
-                  pt_surface, rho_surface, surface_pressure
+                  lsf_surf, microphysics_morrison, microphysics_seifert,       &
+                  microphysics_kessler, pt_surface, rho_surface,               &
+                  surface_pressure
 
        USE indices,                                                            &
@@ -1448 +1814 @@
        pt_1d(:) = pt(:,j,i)
        qc_1d(:) = qc(:,j,i)
-       nc_1d(:) = nc_const
        IF ( microphysics_seifert )  THEN
           qr_1d(:) = qr(:,j,i)
           nr_1d(:) = nr(:,j,i)
        ENDIF
+       IF ( microphysics_morrison )  THEN
+          nc_1d(:) = nc(:,j,i)
+       ENDIF
+
 
 !
@@ -1468 +1837 @@
 
           CALL adjust_cloud( i,j )
+          IF ( microphysics_morrison ) CALL activation( i,j )
+          IF ( microphysics_morrison ) CALL condensation( i,j )
           CALL autoconversion( i,j )
           CALL accretion( i,j )
@@ -1483 +1854 @@
        q(:,j,i)  = q_1d(:)
        pt(:,j,i) = pt_1d(:)
+       IF ( microphysics_morrison ) THEN
+          qc(:,j,i) = qc_1d(:)
+          nc(:,j,i) = nc_1d(:)
+       ENDIF
        IF ( microphysics_seifert )  THEN
           qr(:,j,i) = qr_1d(:)
@@ -1503 +1878 @@
        USE cloud_parameters,                                                   &
            ONLY:  hyrho
+
+       USE control_parameters,                                                 &
+           ONLY: microphysics_morrison
 
        USE indices,                                                            &
@@ -1538 +1916 @@
           ENDIF
 
+          IF ( microphysics_morrison ) THEN
+             IF ( qc_1d(k) <= eps_sb )  THEN
+                qc_1d(k) = 0.0_wp
+                nc_1d(k) = 0.0_wp
+             ELSE
+                IF ( nc_1d(k) * xcmin > qc_1d(k) * hyrho(k) )  THEN
+                   nc_1d(k) = qc_1d(k) * hyrho(k) / xamin * flag
+                ENDIF
+             ENDIF
+          ENDIF
+
        ENDDO
 
     END SUBROUTINE adjust_cloud_ij
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate number of activated condensation nucleii after simple activation
+!> scheme of Twomey, 1959.
+!------------------------------------------------------------------------------!
+    SUBROUTINE activation_ij( i, j )
+
+       USE arrays_3d,                                                          &
+           ONLY:  hyp, nr, pt, q,  qc, qr, nc
+
+       USE cloud_parameters,                                                   &
+           ONLY:  hyrho, l_d_cp, l_d_r, l_v, rho_l, r_v, t_d_pt
+
+       USE constants,                                                          &
+           ONLY:  pi
+
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
+       USE indices,                                                            &
+           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
+
+       USE kinds
+
+       USE control_parameters,                                                 &
+          ONLY: simulated_time
+
+       USE particle_attributes,                                                &
+          ONLY:  molecular_weight_of_solute, molecular_weight_of_water, rho_s, &
+              s1, s2, s3, vanthoff
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+
+       REAL(wp)     ::  activ             !<
+       REAL(wp)     ::  afactor           !<
+       REAL(wp)     ::  alpha             !<
+       REAL(wp)     ::  beta_act          !<
+       REAL(wp)     ::  bfactor           !<
+       REAL(wp)     ::  e_s               !<
+       REAL(wp)     ::  k_act             !<
+       REAL(wp)     ::  n_act             !<
+       REAL(wp)     ::  n_ccn             !<
+       REAL(wp)     ::  q_s               !<
+       REAL(wp)     ::  rd0               !<
+       REAL(wp)     ::  s_0               !<
+       REAL(wp)     ::  sat               !<
+       REAL(wp)     ::  sat_max           !<
+       REAL(wp)     ::  sigma             !<
+       REAL(wp)     ::  t_int             !<
+       REAL(wp)     ::  t_l               !<
+
+       DO  k = nzb+1, nzt
+!
+!--       Actual liquid water temperature:
+          t_l = t_d_pt(k) * pt_1d(k)
+
+!
+!--       Calculate actual temperature
+          t_int = pt_1d(k) * ( hyp(k) / 100000.0_wp )**0.286_wp
+!
+!--             Saturation vapor pressure at t_l:
+          e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /             &
+                                 ( t_l - 35.86_wp )                            &
+                                 )
+!
+!--       Computation of saturation humidity:
+          q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
+          alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
+          q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) /                         &
+                        ( 1.0_wp + alpha * q_s )
+
+!--       Supersaturation:
+          sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
+
+!
+!--       Prescribe parameters for activation
+!--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
+          k_act  = 0.7_wp 
+
+          IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
+!
+!--          Compute the number of activated Aerosols 
+!--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
+             n_act     = na_init * sat**k_act
+!
+!--          Compute the number of cloud droplets
+!--          (see: Morrison + Grabowski, 2007, JAS, 64)
+!            activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
+
+!
+!--          Compute activation rate after Khairoutdinov and Kogan
+!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)   
+             sat_max = 0.8_wp / 100.0_wp                  
+             activ   = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN              & 
+                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) /          & 
+                  dt_micro
+
+             nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
+          ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk )  THEN
+!
+!--          Curvature effect (afactor) with surface tension 
+!--          parameterization by Straka (2009)
+             sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
+             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
+!
+!--          Solute effect (bfactor)
+             bfactor = vanthoff * molecular_weight_of_water *                  &
+                  rho_s / ( molecular_weight_of_solute * rho_l )
+
+!
+!--          Prescribe power index that describes the soluble fraction 
+!--          of an aerosol particle (beta) and mean geometric radius of
+!--          dry aerosol spectrum
+!--          (see: Morrison + Grabowski, 2007, JAS, 64)
+             beta_act = 0.5_wp
+             rd0      = 0.05E-6_wp
+!
+!--          Calculate mean geometric supersaturation (s_0) with 
+!--          parameterization by Khvorostyanov and Curry (2006)
+             s_0   = rd0 **(- ( 1.0_wp + beta_act ) ) *                        &
+               ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
+
+!
+!--          Calculate number of activated CCN as a function of 
+!--          supersaturation and taking Koehler theory into account
+!--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
+             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                     & 
+                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(s1) ) ) )      
+             activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
+                                      
+             nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)                   
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE activation_ij
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculate condensation rate for cloud water content (after Khairoutdinov and  
+!> Kogan, 2000).
+!------------------------------------------------------------------------------!
+    SUBROUTINE condensation_ij( i, j )
+
+       USE arrays_3d,                                                          &
+           ONLY:  hyp, nr, pt, q,  qc, qr, nc
+
+       USE cloud_parameters,                                                   &
+           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
+
+       USE constants,                                                          &
+           ONLY:  pi
+
+       USE cpulog,                                                             &
+           ONLY:  cpu_log, log_point_s
+
+       USE indices,                                                            &
+           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
+
+       USE kinds
+
+       USE control_parameters,                                                 &
+          ONLY: simulated_time
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i                 !<
+       INTEGER(iwp) ::  j                 !<
+       INTEGER(iwp) ::  k                 !<
+
+       REAL(wp)     ::  alpha             !<
+       REAL(wp)     ::  cond              !<
+       REAL(wp)     ::  cond_max          !<
+       REAL(wp)     ::  dc                !<
+       REAL(wp)     ::  e_s               !<
+       REAL(wp)     ::  evap              !<
+       REAL(wp)     ::  evap_nc           !<
+       REAL(wp)     ::  g_fac             !<
+       REAL(wp)     ::  nc_0              !<
+       REAL(wp)     ::  q_s               !<
+       REAL(wp)     ::  sat               !<
+       REAL(wp)     ::  t_l               !<
+       REAL(wp)     ::  temp              !<
+       REAL(wp)     ::  xc                !<
+
+
+       DO  k = nzb+1, nzt
+!
+!--       Actual liquid water temperature:
+          t_l = t_d_pt(k) * pt_1d(k)
+!
+!--             Saturation vapor pressure at t_l:
+          e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
+                                 ( t_l - 35.86_wp )                   &
+                                 )
+!
+!--       Computation of saturation humidity:
+          q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
+          alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
+          q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) /               &
+                        ( 1.0_wp + alpha * q_s )
+
+!--       Supersaturation:
+          sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
+
+
+!
+!--       Actual temperature:
+          temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
+
+          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
+                              l_v / ( thermal_conductivity_l * temp )    &
+                              + r_v * temp / ( diff_coeff_l * e_s )      &
+                            )
+!
+!--       Mean weight of cloud drops
+          IF ( nc_1d(k) <= 0.0_wp) CYCLE
+          xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin)                
+!
+!--       Weight averaged diameter of cloud drops:
+          dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
+!
+!--       Integral diameter of cloud drops
+          nc_0 = nc_1d(k) * dc                
+!
+!--       Condensation needs only to be calculated in supersaturated regions
+          IF ( sat > 0.0_wp )  THEN
+!
+!--          Condensation rate of cloud water content
+!--          after KK scheme.
+!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
+             cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
+             cond_max  = q_1d(k) - q_s - qc_1d(k) - qr_1d(k)
+             cond      = MIN( cond, cond_max / dt_micro )
+               
+             qc_1d(k) = qc_1d(k) + cond * dt_micro
+          ELSEIF ( sat < 0.0_wp ) THEN
+             evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
+             evap      = MAX( evap, -qc_1d(k) / dt_micro )
+
+             qc_1d(k) = qc_1d(k) + evap * dt_micro
+          ENDIF
+       ENDDO
+
+    END SUBROUTINE condensation_ij
 
 
@@ -1557 +2198 @@
 
        USE control_parameters,                                                 &
-           ONLY:  rho_surface
+           ONLY:  microphysics_morrison, rho_surface
 
        USE grid_variables,                                                     &
@@ -1579 +2220 @@
        REAL(wp)     ::  k_au              !<
        REAL(wp)     ::  l_mix             !<
+       REAL(wp)     ::  nc_auto           !<
        REAL(wp)     ::  nu_c              !<
        REAL(wp)     ::  phi_au            !<
@@ -1592 +2234 @@
 !--       Predetermine flag to mask topography
           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+          nc_auto = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
 
           IF ( qc_1d(k) > eps_sb )  THEN
@@ -1610 +2253 @@
 !
 !--          Mean weight of cloud droplets:
-             xc = hyrho(k) * qc_1d(k) / nc_1d(k)
+             xc = hyrho(k) * qc_1d(k) / nc_auto
 !
 !--          Parameterized turbulence effects on autoconversion (Seifert,
@@ -1658 +2301 @@
              qc_1d(k) = qc_1d(k) - autocon * dt_micro                 * flag
              nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag
+             IF ( microphysics_morrison ) THEN
+                nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp *                & 
+                                    autocon / x0 * hyrho(k) * dt_micro * flag )
+             ENDIF
 
           ENDIF
@@ -1738 +2385 @@
 
        USE control_parameters,                                                 &
-           ONLY:  rho_surface
+           ONLY:  microphysics_morrison, rho_surface
 
        USE indices,                                                            &
@@ -1754 +2401 @@
        REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
        REAL(wp)     ::  k_cr              !<
+       REAL(wp)     ::  nc_accr           !<
        REAL(wp)     ::  phi_ac            !<
        REAL(wp)     ::  tau_cloud         !<
+       REAL(wp)     ::  xc                !<
+
 
        DO  k = nzb+1, nzt
@@ -1761 +2411 @@
 !--       Predetermine flag to mask topography
           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
-
-          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
+          nc_accr = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
+
+          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb )  .AND.  &
+               ( nc_accr > eps_mr ) )  THEN
 !
 !--          Intern time scale of coagulation (Seifert and Beheng, 2006):
@@ -1770 +2422 @@
 !--          (Seifert and Beheng, 2001):
              phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
+
+!
+!--          Mean weight of cloud drops
+             xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin)      
 !
 !--          Parameterized turbulence effects on autoconversion (Seifert,
@@ -1784 +2440 @@
 !
 !--          Accretion rate (Seifert and Beheng, 2006):
-             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
+             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac *                      &
+                    SQRT( rho_surface * hyrho(k) )
              accr = MIN( accr, qc_1d(k) / dt_micro )
 
              qr_1d(k) = qr_1d(k) + accr * dt_micro * flag
              qc_1d(k) = qc_1d(k) - accr * dt_micro * flag
+             IF ( microphysics_morrison )  THEN
+                nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), accr / xc *               &
+                                             hyrho(k) * dt_micro * flag        &
+                                         )
+             ENDIF
+
 
           ENDIF
@@ -2018 +2681 @@
 
        USE control_parameters,                                                 &
-           ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
+           ONLY:  call_microphysics_at_all_substeps,                           &
+                  intermediate_timestep_count, microphysics_morrison
 
        USE indices,                                                            &
@@ -2034 +2698 @@
        INTEGER(iwp) ::  k             !<
 
-       REAL(wp)                       ::  flag    !< flag to indicate first grid level above surface
+       REAL(wp)     ::  flag    !< flag to indicate first grid level above surface
+       REAL(wp)     ::  nc_sedi !<
+
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc  !<
        REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc  !<
 
@@ -2043 +2710 @@
 !--       Predetermine flag to mask topography
           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+          nc_sedi = MERGE( nc_1d(k), nc_const, microphysics_morrison )
+!
+!--       Sedimentation fluxes for number concentration are only calculated
+!--       for cloud_scheme = 'morrison'
+          IF ( microphysics_morrison ) THEN
+             IF ( qc_1d(k) > eps_sb  .AND.  nc_1d(k) > eps_mr )  THEN
+                sed_nc(k) = sed_qc_const *                                     & 
+                            ( qc_1d(k) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *     &
+                            ( nc_1d(k) )**( 1.0_wp / 3.0_wp )
+             ELSE
+                sed_nc(k) = 0.0_wp
+             ENDIF
+
+             sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *                 &
+                              nc_1d(k) / dt_micro + sed_nc(k+1)                &
+                            ) * flag
+
+             nc_1d(k) = nc_1d(k) + ( sed_nc(k+1) - sed_nc(k) ) *               &
+                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
+          ENDIF
 
           IF ( qc_1d(k) > eps_sb )  THEN
-             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) *       &
+             sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *        &
                          ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )  * flag
           ELSE

palm/trunk/SOURCE/modules.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2277 2017-06-12 10:47:51Z kanani
 ! Added doxygen comments for variables/parameters,
 ! removed unused variables dissipation_control, do2d_xy_n, do2d_xz_n, do2d_yz_n, 
@@ -543 +548 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  c_w                   !< phase speed of w-velocity component
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_e              !< artificial numerical dissipation flux at south face of grid box - subgrid-scale TKE
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_nc             !< artificial numerical dissipation flux at south face of grid box - clouddrop-number concentration    
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_nr             !< artificial numerical dissipation flux at south face of grid box - raindrop-number concentration    
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_pt             !< artificial numerical dissipation flux at south face of grid box - potential temperature 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_q              !< artificial numerical dissipation flux at south face of grid box - mixing ratio 
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qc             !< artificial numerical dissipation flux at south face of grid box - cloudwater 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qr             !< artificial numerical dissipation flux at south face of grid box - rainwater 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_s              !< artificial numerical dissipation flux at south face of grid box - passive scalar
@@ -555 +562 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dzw_mg                !< vertical grid size (w-grid) for multigrid pressure solver
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_e              !< 6th-order advective flux at south face of grid box - subgrid-scale TKE
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_nc             !< 6th-order advective flux at south face of grid box - clouddrop-number concentration 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_nr             !< 6th-order advective flux at south face of grid box - raindrop-number concentration 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_pt             !< 6th-order advective flux at south face of grid box - potential temperature 
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_q              !< 6th-order advective flux at south face of grid box - mixing ratio 
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_qc             !< 6th-order advective flux at south face of grid box - cloudwater
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_qr             !< 6th-order advective flux at south face of grid box - rainwater
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_s              !< 6th-order advective flux at south face of grid box - passive scalar
@@ -594 +603 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss       !< TKE dissipation rate
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_e   !< artificial numerical dissipation flux at left face of grid box - subgrid-scale TKE
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_nc  !< artificial numerical dissipation flux at left face of grid box - clouddrop-number concentration
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_nr  !< artificial numerical dissipation flux at left face of grid box - raindrop-number concentration
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_pt  !< artificial numerical dissipation flux at left face of grid box - potential temperature
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_q   !< artificial numerical dissipation flux at left face of grid box - mixing ratio
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_qc  !< artificial numerical dissipation flux at left face of grid box - cloudwater
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_qr  !< artificial numerical dissipation flux at left face of grid box - rainwater
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_s   !< artificial numerical dissipation flux at left face of grid box - passive scalar
@@ -604 +615 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_w   !< artificial numerical dissipation flux at left face of grid box - w-component
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_e   !< 6th-order advective flux at south face of grid box - subgrid-scale TKE
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_nc  !< 6th-order advective flux at south face of grid box - clouddrop-number concentration
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_nr  !< 6th-order advective flux at south face of grid box - raindrop-number concentration
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_pt  !< 6th-order advective flux at south face of grid box - potential temperature
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_q   !< 6th-order advective flux at south face of grid box - mixing ratio
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_qc  !< 6th-order advective flux at south face of grid box - cloudwater
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_qr  !< 6th-order advective flux at south face of grid box - rainwater
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_s   !< 6th-order advective flux at south face of grid box - passive scalar
@@ -636 +649 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  e          !< subgrid-scale turbulence kinetic energy (sgs tke)
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  e_p        !< prognostic value of sgs tke
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc         !< cloud drop number density
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_p       !< prognostic value of cloud drop number density
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr         !< rain drop number density
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_p       !< prognostic value of rain drop number density
@@ -646 +661 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  q_p        !< prognostic value of specific humidity 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc         !< cloud water content
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_p       !< cloud water content
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql         !< liquid water content
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_c       !< change in liquid water content due to 
@@ -659 +675 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  sa_p       !< prognostic value of ocean salinity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  te_m       !< weighted tendency of e for previous sub-timestep (Runge-Kutta)
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  tnc_m      !< weighted tendency of nc for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  tnr_m      !< weighted tendency of nr for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  tpt_m      !< weighted tendency of pt for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  tq_m       !< weighted tendency of q for previous sub-timestep (Runge-Kutta)
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  tqc_m      !< weighted tendency of qc for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  tqr_m      !< weighted tendency of qr for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ts_m       !< weighted tendency of s for previous sub-timestep (Runge-Kutta)
@@ -681 +699 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  p       !< pointer: perturbation pressure
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  prho_1  !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_1    !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_2    !< pointer for swapping of timelevels for respective quantity 
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_3    !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_1    !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_2    !< pointer for swapping of timelevels for respective quantity 
@@ -691 +712 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  q_3     !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_1    !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_2    !< pointer for swapping of timelevels for respective quantity
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_3    !< pointer for swapping of timelevels for respective quantity
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_v    !< pointer: volume of liquid water
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_vp   !< pointer: liquid water weighting factor
@@ -718 +741 @@
     REAL(wp), DIMENSION(:,:,:), POINTER ::  e          !< pointer: subgrid-scale turbulence kinetic energy (sgs tke)
     REAL(wp), DIMENSION(:,:,:), POINTER ::  e_p        !< pointer: prognostic value of sgs tke
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  nc         !< pointer: cloud drop number density
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  nc_p       !< pointer: prognostic value of cloud drop number density
     REAL(wp), DIMENSION(:,:,:), POINTER ::  nr         !< pointer: rain drop number density
     REAL(wp), DIMENSION(:,:,:), POINTER ::  nr_p       !< pointer: prognostic value of rain drop number density
@@ -726 +751 @@
     REAL(wp), DIMENSION(:,:,:), POINTER ::  q_p        !< pointer: prognostic value of specific humidity
     REAL(wp), DIMENSION(:,:,:), POINTER ::  qc         !< pointer: cloud water content
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  qc_p       !< pointer: cloud water content
     REAL(wp), DIMENSION(:,:,:), POINTER ::  ql         !< pointer: liquid water content
     REAL(wp), DIMENSION(:,:,:), POINTER ::  ql_c       !< pointer: change in liquid water content due to 
@@ -737 +763 @@
     REAL(wp), DIMENSION(:,:,:), POINTER ::  sa_p       !< pointer: prognostic value of ocean salinity
     REAL(wp), DIMENSION(:,:,:), POINTER ::  te_m       !< pointer: weighted tendency of e for previous sub-timestep (Runge-Kutta)
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  tnc_m      !< pointer: weighted tendency of nc for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), POINTER ::  tnr_m      !< pointer: weighted tendency of nr for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), POINTER ::  tpt_m      !< pointer: weighted tendency of pt for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), POINTER ::  tq_m       !< pointer: weighted tendency of q for previous sub-timestep (Runge-Kutta)
+    REAL(wp), DIMENSION(:,:,:), POINTER ::  tqc_m      !< pointer: weighted tendency of qc for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), POINTER ::  tqr_m      !< pointer: weighted tendency of qr for previous sub-timestep (Runge-Kutta)
     REAL(wp), DIMENSION(:,:,:), POINTER ::  ts_m       !< pointer: weighted tendency of s for previous sub-timestep (Runge-Kutta)
@@ -800 +828 @@
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  e_av          !< avg. subgrid-scale tke
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  lpt_av        !< avg. liquid water potential temperature
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_av         !< avg. cloud drop number density
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_av         !< avg. rain drop number density
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  p_av          !< avg. perturbation pressure
@@ -1144 +1173 @@
     LOGICAL ::  microphysics_sat_adjust = .FALSE.            !< use saturation adjust bulk scheme?
     LOGICAL ::  microphysics_kessler = .FALSE.               !< use kessler bulk scheme?
-    LOGICAL ::  microphysics_seifert = .FALSE.               !< use 2-moment Seifert and Beheng scheme?
+    LOGICAL ::  microphysics_morrison = .FALSE.              !< use 2-moment Morrison (add. prog. eq. for nc and qc)
+    LOGICAL ::  microphysics_seifert = .FALSE.               !< use 2-moment Seifert and Beheng scheme
     LOGICAL ::  mg_switch_to_pe0 = .FALSE.                   !< (Siggi add short description)
     LOGICAL ::  nest_bound_l = .FALSE.                       !< nested boundary on left side?
@@ -1980 +2010 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_vs2_ws_l    !< subdomain sum of horizontal momentum flux v'v' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_ws2_ws_l    !< subdomain sum of vertical momentum flux w'w' (5th-order advection scheme only)
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsncs_ws_l  !< subdomain sum of vertical clouddrop-number concentration flux w'nc' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsnrs_ws_l  !< subdomain sum of vertical raindrop-number concentration flux w'nr' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wspts_ws_l  !< subdomain sum of vertical sensible heat flux w'pt' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqs_ws_l   !< subdomain sum of vertical latent heat flux w'q' (5th-order advection scheme only)
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqcs_ws_l  !< subdomain sum of vertical cloudwater flux w'qc' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqrs_ws_l  !< subdomain sum of vertical rainwater flux w'qr' (5th-order advection scheme only)
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wssas_ws_l  !< subdomain sum of vertical salinity flux w'sa' (5th-order advection scheme only)

palm/trunk/SOURCE/netcdf_interface_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2270 2017-06-09 12:18:47Z maronga
 ! Removed 2 timeseries (shf_eb + qsws_eb). Removed _eb suffixes
 ! 
@@ -844 +849 @@
 !
 !--             Most variables are defined on the scalar grid
-                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
-                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',        &
+                CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
+                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
                        'rho_ocean', 's', 'sa', 'vpt'  )
 
@@ -1376 +1381 @@
 !
 !--             Most variables are defined on the scalar grid
-                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
-                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho_ocean', &
-                       's', 'sa', 'vpt' )
+                CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
+                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
+                       'rho_ocean', 's', 'sa', 'vpt' )
 
                    grid_x = 'x'
@@ -2027 +2032 @@
 !
 !--                   Most variables are defined on the zu grid
-                      CASE ( 'e_xy', 'lpt_xy', 'nr_xy', 'p_xy', 'pc_xy',       &
-                             'pr_xy', 'prr_xy', 'pt_xy', 'q_xy', 'qc_xy',      &
-                             'ql_xy', 'ql_c_xy', 'ql_v_xy', 'ql_vp_xy',        &
-                             'qr_xy', 'qv_xy', 'rho_ocean_xy', 's_xy',         &
-                             'sa_xy', 'vpt_xy' )
+                      CASE ( 'e_xy', 'lpt_xy', 'nc_xy', 'nr_xy', 'p_xy',       &
+                             'pc_xy', 'pr_xy', 'prr_xy', 'pt_xy', 'q_xy',      &
+                             'qc_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy',           &
+                             'ql_vp_xy', 'qr_xy', 'qv_xy', 'rho_ocean_xy',     &
+                             's_xy', 'sa_xy', 'vpt_xy' )
 
                          grid_x = 'x'
@@ -2726 +2731 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_xz', 'lpt_xz', 'nr_xz', 'p_xz', 'pc_xz', 'pr_xz', &
-                          'prr_xz', 'pt_xz', 'q_xz', 'qc_xz', 'ql_xz',         &
-                          'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz', 'qv_xz',  &
-                          'rho_ocean_xz', 's_xz', 'sa_xz', 'vpt_xz' )
+                   CASE ( 'e_xz', 'lpt_xz', 'nc_xz', 'nr_xz', 'p_xz', 'pc_xz', &
+                          'pr_xz', 'prr_xz', 'pt_xz', 'q_xz', 'qc_xz',         &
+                          'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz',  &
+                          'qv_xz', 'rho_ocean_xz', 's_xz', 'sa_xz', 'vpt_xz' )
 
                       grid_x = 'x'
@@ -3380 +3385 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_yz', 'lpt_yz', 'nr_yz', 'p_yz', 'pc_yz', 'pr_yz', &
-                          'prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz',         &
+                   CASE ( 'e_yz', 'lpt_yz', 'nc_yz', 'nr_yz', 'p_yz', 'pc_yz', &
+                          'pr_yz','prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz', &
                           'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz',  &
                           'rho_ocean_yz', 's_yz', 'sa_yz', 'vpt_yz' )
@@ -4557 +4562 @@
 !
 !--             Most variables are defined on the zu levels
-                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
-                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',        &
+                CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
+                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
                        'rho_ocean', 's', 'sa', 'u', 'v', 'vpt' )
 

palm/trunk/SOURCE/palm.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2261 2017-06-08 14:25:57Z raasch
 ! output of run number for mrun to create unified cycle numbers
 ! 
@@ -177 +182 @@
                do2d_at_begin, do3d_at_begin, humidity, initializing_actions,   &
                io_blocks, io_group, land_surface, large_scale_forcing,         &
-               message_string, microphysics_seifert, nest_domain, neutral,     &
-               nudging, passive_scalar, runnr, simulated_time,                 &
-               simulated_time_chr, urban_surface,                              &
+               message_string, microphysics_morrison, microphysics_seifert,    &
+               nest_domain, neutral, nudging, passive_scalar, runnr,           &
+               simulated_time, simulated_time_chr, urban_surface,              &
                user_interface_current_revision,                                &
                user_interface_required_revision, version, wall_heatflux,       &
@@ -364 +369 @@
              IF ( humidity )  THEN
                 CALL exchange_horiz( q, nbgp )
+                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+                  CALL exchange_horiz( qc, nbgp )
+                  CALL exchange_horiz( nc, nbgp )
+                ENDIF
                 IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-!                   CALL exchange_horiz( qc, nbgp )
                    CALL exchange_horiz( qr, nbgp ) 
-!                   CALL exchange_horiz( nc, nbgp )
                    CALL exchange_horiz( nr, nbgp )
                 ENDIF

palm/trunk/SOURCE/parin.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2267 2017-06-09 09:33:25Z gronemeier
 ! Bugfix: removed skipping of reading namelists in case of omitted d3par
 ! 
@@ -292 +297 @@
     USE microphysics_mod,                                                      &
         ONLY:  c_sedimentation, cloud_water_sedimentation,                     &
-               collision_turbulence, limiter_sedimentation, nc_const,          &
-               ventilation_effect
+               collision_turbulence, curvature_solution_effects_bulk,          &
+               limiter_sedimentation ,na_init, nc_const, ventilation_effect
 
     USE model_1d,                                                              &
@@ -350 +355 @@
              collective_wait, collision_turbulence, conserve_volume_flow,      &
              conserve_volume_flow_mode, constant_flux_layer,                   &
-             coupling_start_time,                                              &
+             coupling_start_time, curvature_solution_effects_bulk,             &
              cycle_mg, damp_level_1d,                                          &
              dissipation_1d,                                                   &
@@ -366 +371 @@
              loop_optimization, lsf_exception, masking_method, mg_cycles,      &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
-             most_method, nc_const, netcdf_precision, neutral, ngsrb,          &
+             most_method, na_init, nc_const, netcdf_precision, neutral, ngsrb, &
              nsor, nsor_ini, nudging, nx, ny, nz, ocean, omega, omega_sor,     &
              outflow_source_plane, passive_scalar, phi,                        &

palm/trunk/SOURCE/pmc_interface_mod.f90

@@ -26 +26 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2285 2017-06-15 13:15:41Z suehring
 ! Consider multiple pure-vertical nest domains in overlap check
 ! 
@@ -154 +159 @@
 #if defined( __nopointer )
     USE arrays_3d,                                                             &
-        ONLY:  dzu, dzw, e, e_p, nr, pt, pt_p, q, q_p, qr, u, u_p, v, v_p,     &
-               w, w_p, zu, zw
+        ONLY:  dzu, dzw, e, e_p, nc, nr, pt, pt_p, q, q_p, qc, qr, u, u_p, v,  &
+               v_p, w, w_p, zu, zw
 #else
    USE arrays_3d,                                                              &
-        ONLY:  dzu, dzw, e, e_p, e_1, e_2, nr, nr_2, nr_p, pt, pt_p, pt_1,     &
-               pt_2, q, q_p, q_1, q_2, qr, qr_2, s, s_2, u, u_p, u_1, u_2, v,  &
-               v_p, v_1, v_2, w, w_p, w_1, w_2, zu, zw
+        ONLY:  dzu, dzw, e, e_p, e_1, e_2, nc, nc_2, nc_p, nr, nr_2, nr_p, pt, &
+               pt_p, pt_1, pt_2, q, q_p, q_1, q_2, qc, qc_2, qr, qr_2, s, s_2, &
+               u, u_p, u_1, u_2, v, v_p, v_1, v_2, w, w_p, w_1, w_2, zu, zw
 #endif
 
     USE control_parameters,                                                     &
         ONLY:  cloud_physics, coupling_char, dt_3d, dz, humidity,               &
-               message_string, microphysics_seifert, nest_bound_l, nest_bound_r,&
-               nest_bound_s, nest_bound_n, nest_domain, neutral, passive_scalar,& 
-               roughness_length, simulated_time, topography, volume_flow
+               message_string, microphysics_morrison, microphysics_seifert,     &
+               nest_bound_l, nest_bound_r, nest_bound_s, nest_bound_n,          &
+               nest_domain, neutral, passive_scalar, roughness_length,          &
+               simulated_time, topography, volume_flow
 
     USE cpulog,                                                                 &
@@ -281 +287 @@
     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wc   !:
     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  q_c  !:
-!     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qcc  !:
+    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qcc  !:
     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qrc  !:
     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nrc  !:
-!     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ncc  !:
+    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ncc  !:
     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  sc   !:
 
@@ -966 +972 @@
           CALL pmc_set_dataarray_name( 'coarse', 'q'  ,'fine', 'q',  ierr )
 
+          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+            CALL pmc_set_dataarray_name( 'coarse', 'qc'  ,'fine', 'qc',  ierr )  
+            CALL pmc_set_dataarray_name( 'coarse', 'nc'  ,'fine', 'nc',  ierr ) 
+          ENDIF
+
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-!              CALL pmc_set_dataarray_name( 'coarse', 'qc'  ,'fine', 'qc',  ierr )  
              CALL pmc_set_dataarray_name( 'coarse', 'qr'  ,'fine', 'qr',  ierr )
-!             CALL pmc_set_dataarray_name( 'coarse', 'nc'  ,'fine', 'nc',  ierr ) 
              CALL pmc_set_dataarray_name( 'coarse', 'nr'  ,'fine', 'nr',  ierr ) 
-
           ENDIF
      
@@ -2990 +2998 @@
     IF ( TRIM(name) == "pt" )  p_3d => pt
     IF ( TRIM(name) == "q"  )  p_3d => q
-!     IF ( TRIM(name) == "qc" )  p_3d => qc
+    IF ( TRIM(name) == "qc" )  p_3d => qc
     IF ( TRIM(name) == "qr" )  p_3d => qr
     IF ( TRIM(name) == "nr" )  p_3d => nr
-!     IF ( TRIM(name) == "nc" )  p_3d => nc
+    IF ( TRIM(name) == "nc" )  p_3d => nc
     IF ( TRIM(name) == "s"  )  p_3d => s
 !
@@ -3026 +3034 @@
     IF ( TRIM(name) == "pt" )  p_3d_sec => pt_2
     IF ( TRIM(name) == "q"  )  p_3d_sec => q_2
-!     IF ( TRIM(name) == "qc" )  p_3d_sec => qc_2
+    IF ( TRIM(name) == "qc" )  p_3d_sec => qc_2
     IF ( TRIM(name) == "qr" )  p_3d_sec => qr_2
     IF ( TRIM(name) == "nr" )  p_3d_sec => nr_2
-!     IF ( TRIM(name) == "nc" )  p_3d_sec => nc_2
+    IF ( TRIM(name) == "nc" )  p_3d_sec => nc_2
     IF ( TRIM(name) == "s"  )  p_3d_sec => s_2
 
@@ -3102 +3110 @@
        IF ( .NOT. ALLOCATED( q_c ) ) ALLOCATE( q_c(0:nzc+1, js:je, is:ie) )
        p_3d => q_c
-!     ELSEIF ( TRIM( name ) == "qc")  THEN
-!        IF ( .NOT. ALLOCATED( qcc ) ) ALLOCATE( qcc(0:nzc+1, js:je, is:ie) )
-!        p_3d => qcc
+    ELSEIF ( TRIM( name ) == "qc")  THEN
+       IF ( .NOT. ALLOCATED( qcc ) ) ALLOCATE( qcc(0:nzc+1, js:je, is:ie) )
+       p_3d => qcc
     ELSEIF ( TRIM( name ) == "qr")  THEN
        IF ( .NOT. ALLOCATED( qrc ) ) ALLOCATE( qrc(0:nzc+1, js:je, is:ie) )
@@ -3111 +3119 @@
        IF ( .NOT. ALLOCATED( nrc ) ) ALLOCATE( nrc(0:nzc+1, js:je, is:ie) )
        p_3d => nrc
-!     ELSEIF ( TRIM( name ) == "nc")  THEN
-!        IF ( .NOT. ALLOCATED( ncc ) ) ALLOCATE( ncc(0:nzc+1, js:je, is:ie) )
-!        p_3d => ncc
+    ELSEIF ( TRIM( name ) == "nc")  THEN
+       IF ( .NOT. ALLOCATED( ncc ) ) ALLOCATE( ncc(0:nzc+1, js:je, is:ie) )
+       p_3d => ncc
     ELSEIF ( TRIM( name ) == "s")  THEN
        IF ( .NOT. ALLOCATED( sc ) ) ALLOCATE( sc(0:nzc+1, js:je, is:ie) )
@@ -3224 +3232 @@
                                       r2yo, r1zo, r2zo, 's' )
 
+          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+             CALL pmci_interp_tril_all ( qc, qcc, ico, jco, kco, r1xo, r2xo,    &
+                                          r1yo, r2yo, r1zo, r2zo, 's' ) 
+             CALL pmci_interp_tril_all ( nc, ncc, ico, jco, kco, r1xo, r2xo,    &
+                                         r1yo, r2yo, r1zo, r2zo, 's' )   
+          ENDIF
+
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-!              CALL pmci_interp_tril_all ( qc, qcc, ico, jco, kco, r1xo, r2xo,    &
-!                                          r1yo, r2yo, r1zo, r2zo, 's' ) 
              CALL pmci_interp_tril_all ( qr, qrc, ico, jco, kco, r1xo, r2xo,    &
                                          r1yo, r2yo, r1zo, r2zo, 's' )
-!             CALL pmci_interp_tril_all ( nc, ncc, ico, jco, kco, r1xo, r2xo,    &
-!                                         r1yo, r2yo, r1zo, r2zo, 's' )   
              CALL pmci_interp_tril_all ( nr, nrc, ico, jco, kco, r1xo, r2xo,    &
                                          r1yo, r2yo, r1zo, r2zo, 's' )
@@ -3955 +3966 @@
                                           nzt_topo_nestbc_l, 'l', 's' )
 
+                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+                   CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo, r2xo,&
+                                             r1yo, r2yo, r1zo, r2zo,            &
+                                             logc_u_l,                          &
+                                             logc_ratio_u_l, nzt_topo_nestbc_l, &
+                                             'l', 's' )  
+
+                   CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo, r2xo,&
+                                             r1yo, r2yo, r1zo, r2zo,            &
+                                             logc_u_l,                          &
+                                             logc_ratio_u_l, nzt_topo_nestbc_l, &
+                                             'l', 's' )          
+                ENDIF
 
                 IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-!                    CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo, r2xo,&
-!                                              r1yo, r2yo, r1zo, r2zo,            &
-!                                              logc_u_l,                          &
-!                                              logc_ratio_u_l, nzt_topo_nestbc_l, &
-!                                              'l', 's' )  
-
                    CALL pmci_interp_tril_lr( qr, qrc, ico, jco, kco, r1xo, r2xo,&
                                              r1yo, r2yo, r1zo, r2zo,            &
@@ -3968 +3986 @@
                                              logc_ratio_u_l, nzt_topo_nestbc_l, &
                                              'l', 's' ) 
-
-!                    CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo, r2xo,&
-!                                              r1yo, r2yo, r1zo, r2zo,            &
-!                                              logc_u_l,                          &
-!                                              logc_ratio_u_l, nzt_topo_nestbc_l, &
-!                                              'l', 's' ) 
 
                    CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo, r2xo,&
@@ -4004 +4016 @@
                    CALL pmci_extrap_ifoutflow_lr( q, 'l', 's' )
 
+                   IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+
+                      CALL pmci_extrap_ifoutflow_lr( qc, 'l', 's' )
+                      CALL pmci_extrap_ifoutflow_lr( nc, 'l', 's' )
+
+                   ENDIF
+
                    IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
 
-!                       CALL pmci_extrap_ifoutflow_lr( qc, 'l', 's' )
                       CALL pmci_extrap_ifoutflow_lr( qr, 'l', 's' )
-!                      CALL pmci_extrap_ifoutflow_lr( nc, 'l', 's' )
                       CALL pmci_extrap_ifoutflow_lr( nr, 'l', 's' )
 
@@ -4061 +4078 @@
                                           nzt_topo_nestbc_r, 'r', 's' )
 
+                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+
+                   CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo,     &
+                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
+                                             logc_u_r,                         &
+                                             logc_ratio_u_r, nzt_topo_nestbc_r,&
+                                             'r', 's' ) 
+     
+                   CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo,     &
+                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
+                                             logc_u_r,                         &
+                                             logc_ratio_u_r, nzt_topo_nestbc_r,&
+                                             'r', 's' ) 
+
+
+                ENDIF
 
                 IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
 
-!                    CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo,     &
-!                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
-!                                              logc_u_r,                         &
-!                                              logc_ratio_u_r, nzt_topo_nestbc_r,&
-!                                              'r', 's' ) 
      
                    CALL pmci_interp_tril_lr( qr, qrc, ico, jco, kco, r1xo,     &
@@ -4076 +4104 @@
                                              'r', 's' ) 
 
-!                    CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo,     &
-!                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
-!                                              logc_u_r,                         &
-!                                              logc_ratio_u_r, nzt_topo_nestbc_r,&
-!                                              'r', 's' ) 
-
                    CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo,     &
                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
@@ -4113 +4135 @@
                    CALL pmci_extrap_ifoutflow_lr( q, 'r', 's' )
 
+                   IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+                      CALL pmci_extrap_ifoutflow_lr( qc, 'r', 's' )
+                      CALL pmci_extrap_ifoutflow_lr( nc, 'r', 's' )  
+                   ENDIF
+
                    IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-!                       CALL pmci_extrap_ifoutflow_lr( qc, 'r', 's' )
                       CALL pmci_extrap_ifoutflow_lr( qr, 'r', 's' )
-!                      CALL pmci_extrap_ifoutflow_lr( nc, 'r', 's' )  
                       CALL pmci_extrap_ifoutflow_lr( nr, 'r', 's' )
                    ENDIF
@@ -4162 +4187 @@
                                           nzt_topo_nestbc_s, 's', 's' )
 
+                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+
+                   CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
+                                             r2xo, r1yo,r2yo, r1zo, r2zo,      &
+                                             logc_u_s,                         &
+                                             logc_ratio_u_s, nzt_topo_nestbc_s,&
+                                             's', 's' )
+
+                   CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
+                                             r2xo, r1yo,r2yo, r1zo, r2zo,      &
+                                             logc_u_s,                         &
+                                             logc_ratio_u_s, nzt_topo_nestbc_s,&
+                                             's', 's' )
+
+                ENDIF
+
                 IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-
-!                    CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
-!                                              r2xo, r1yo,r2yo, r1zo, r2zo,      &
-!                                              logc_u_s,                         &
-!                                              logc_ratio_u_s, nzt_topo_nestbc_s,&
-!                                              's', 's' )
 
                    CALL pmci_interp_tril_sn( qr, qrc, ico, jco, kco, r1xo,     &
@@ -4175 +4210 @@
                                              logc_ratio_u_s, nzt_topo_nestbc_s,&
                                              's', 's' )
-
-!                    CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
-!                                              r2xo, r1yo,r2yo, r1zo, r2zo,      &
-!                                              logc_u_s,                         &
-!                                              logc_ratio_u_s, nzt_topo_nestbc_s,&
-!                                              's', 's' )
 
                    CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo,     &
@@ -4212 +4241 @@
                    CALL pmci_extrap_ifoutflow_sn( q,  's', 's' )
 
+                   IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+                      CALL pmci_extrap_ifoutflow_sn( qc, 's', 's' ) 
+                      CALL pmci_extrap_ifoutflow_sn( nc, 's', 's' ) 
+                   ENDIF
+
                    IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-!                       CALL pmci_extrap_ifoutflow_sn( qc, 's', 's' ) 
                       CALL pmci_extrap_ifoutflow_sn( qr, 's', 's' )      
-!                       CALL pmci_extrap_ifoutflow_sn( nc, 's', 's' ) 
                       CALL pmci_extrap_ifoutflow_sn( nr, 's', 's' ) 
 
@@ -4266 +4298 @@
                                           nzt_topo_nestbc_n, 'n', 's' )
 
+                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+
+                   CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
+                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
+                                             logc_u_n,                         &
+                                             logc_ratio_u_n, nzt_topo_nestbc_n,& 
+                                             'n', 's' )
+
+                   CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
+                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
+                                             logc_u_n,                         &
+                                             logc_ratio_u_n, nzt_topo_nestbc_n,& 
+                                             'n', 's' )
+
+                ENDIF
+
                 IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-
-!                    CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
-!                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
-!                                              logc_u_n,                         &
-!                                              logc_ratio_u_n, nzt_topo_nestbc_n,& 
-!                                              'n', 's' )
 
                    CALL pmci_interp_tril_sn( qr, qrc, ico, jco, kco, r1xo,     &
@@ -4279 +4321 @@
                                              logc_ratio_u_n, nzt_topo_nestbc_n,& 
                                              'n', 's' )
-
-!                    CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
-!                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
-!                                              logc_u_n,                         &
-!                                              logc_ratio_u_n, nzt_topo_nestbc_n,& 
-!                                              'n', 's' )
 
                    CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo,     &
@@ -4317 +4353 @@
 
                    IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-!                       CALL pmci_extrap_ifoutflow_sn( qc, 'n', 's' )
+                      CALL pmci_extrap_ifoutflow_sn( qc, 'n', 's' )
+                      CALL pmci_extrap_ifoutflow_sn( nc, 'n', 's' )
+                   ENDIF
+
+                   IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
                       CALL pmci_extrap_ifoutflow_sn( qr, 'n', 's' )
-!                       CALL pmci_extrap_ifoutflow_sn( nc, 'n', 's' )
                       CALL pmci_extrap_ifoutflow_sn( nr, 'n', 's' )
                    ENDIF
@@ -4356 +4395 @@
                                    r2yo, r1zo, r2zo, 's' )
 
+          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+
+             CALL pmci_interp_tril_t( qc, qcc, ico, jco, kco, r1xo, r2xo, r1yo,&
+                                      r2yo, r1zo, r2zo, 's' )
+
+             CALL pmci_interp_tril_t( nc, ncc, ico, jco, kco, r1xo, r2xo, r1yo,&
+                                      r2yo, r1zo, r2zo, 's' )
+
+          ENDIF
+
           IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
 
-!              CALL pmci_interp_tril_t( qc, qcc, ico, jco, kco, r1xo, r2xo, r1yo,&
-!                                       r2yo, r1zo, r2zo, 's' )
 
              CALL pmci_interp_tril_t( qr, qrc, ico, jco, kco, r1xo, r2xo, r1yo,&
                                       r2yo, r1zo, r2zo, 's' )
-
-!              CALL pmci_interp_tril_t( nc, ncc, ico, jco, kco, r1xo, r2xo, r1yo,&
-!                                       r2yo, r1zo, r2zo, 's' )
 
              CALL pmci_interp_tril_t( nr, nrc, ico, jco, kco, r1xo, r2xo, r1yo,&
@@ -4394 +4438 @@
              CALL pmci_extrap_ifoutflow_t( q, 's' )
 
+             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+                CALL pmci_extrap_ifoutflow_t( qc, 's' )
+                CALL pmci_extrap_ifoutflow_t( nc, 's' )
+             ENDIF
+
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-!                 CALL pmci_extrap_ifoutflow_t( qc, 's' )
                 CALL pmci_extrap_ifoutflow_t( qr, 's' )
-!                 CALL pmci_extrap_ifoutflow_t( nc, 's' )
                 CALL pmci_extrap_ifoutflow_t( nr, 's' )
 
@@ -4438 +4485 @@
                                   kfuo, ijfc_s, kfc_s, 's' )
 
+         IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+
+            CALL pmci_anterp_tophat( qc, qcc, kctu, iflo, ifuo, jflo, jfuo,    &
+                                     kflo, kfuo, ijfc_s, kfc_s, 's' )
+
+            CALL pmci_anterp_tophat( nc, ncc, kctu, iflo, ifuo, jflo, jfuo,    &
+                                     kflo, kfuo, ijfc_s, kfc_s,  's' )
+
+         ENDIF
+
          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-
-!             CALL pmci_anterp_tophat( qc, qcc, kctu, iflo, ifuo, jflo, jfuo,    &
-!                                      kflo, kfuo, ijfc_s, kfc_s, 's' )
 
             CALL pmci_anterp_tophat( qr, qrc, kctu, iflo, ifuo, jflo, jfuo,    &
                                      kflo, kfuo, ijfc_s, kfc_s, 's' )
-
-!             CALL pmci_anterp_tophat( nc, ncc, kctu, iflo, ifuo, jflo, jfuo,    &
-!                                      kflo, kfuo, ijfc_s, kfc_s,  's' )
 
             CALL pmci_anterp_tophat( nr, nrc, kctu, iflo, ifuo, jflo, jfuo,    &

palm/trunk/SOURCE/prognostic_equations.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2261 2017-06-08 14:25:57Z raasch
 ! bugfix for r2232: openmp directives removed
 ! 
@@ -227 +232 @@
 
     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_s_e, diss_s_nr, diss_s_pt, diss_s_q,  &
-               diss_s_qr, diss_s_s, diss_s_sa, e, e_p, flux_s_e, flux_s_nr,    &
-               flux_s_pt, flux_s_q, flux_s_qr, flux_s_s, flux_s_sa, flux_l_e,  &
-               flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_s, flux_l_sa, &
-               nr, nr_p, pt, ptdf_x, ptdf_y, pt_init, pt_p,                    &
-               prho, q, q_init, q_p, qr, qr_p, rdf,                            &
-               rdf_sc, ref_state, rho_ocean, s, s_init, s_p, sa, sa_init, sa_p,&
-               tend, te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m, tsa_m, tu_m, tv_m, &
-               tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
+        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_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, &
+               e, e_p, 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_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, nc, nc_p, nr, nr_p, pt, ptdf_x, ptdf_y, pt_init,     &
+               pt_p, prho, q, q_init, q_p, qc, qc_p, qr, qr_p, rdf, rdf_sc,    &
+               ref_state, rho_ocean, s,  s_init, s_p, sa, sa_init, sa_p, tend, &
+               te_m, tnc_m,  tnr_m, tpt_m, tq_m, tqc_m, tqr_m, ts_m, tsa_m,    &
+               tu_m, tv_m, tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p,  &
+               w, w_p
 
     USE control_parameters,                                                    &
@@ -244 +251 @@
                inflow_l, intermediate_timestep_count,                          &
                intermediate_timestep_count_max, large_scale_forcing,           &
-               large_scale_subsidence, microphysics_seifert,                   &
+               large_scale_subsidence, microphysics_morrison, microphysics_seifert, &
                microphysics_sat_adjust, neutral, nudging, ocean, outflow_l,    &
                outflow_s, passive_scalar, prho_reference, prho_reference,      &
@@ -892 +899 @@
 
 !
+!--          If required, calculate prognostic equations for cloud water content
+!--          and cloud drop concentration
+             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+!
+!--             Calculate prognostic equation for cloud water content
+                tend(:,j,i) = 0.0_wp
+                IF ( timestep_scheme(1:5) == 'runge' ) &
+                THEN
+                   IF ( ws_scheme_sca )  THEN
+                      CALL advec_s_ws( i, j, qc, 'qc', flux_s_qc,       &
+                                       diss_s_qc, flux_l_qc, diss_l_qc, &
+                                       i_omp_start, tn )
+                   ELSE
+                      CALL advec_s_pw( i, j, qc )
+                   ENDIF
+                ELSE
+                   CALL advec_s_up( i, j, qc )
+                ENDIF
+                CALL diffusion_s( i, j, qc,                                   &
+                                  surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,   &
+                                  surf_def_h(2)%qcsws,                        &
+                                  surf_lsm_h%qcsws,    surf_usm_h%qcsws,      &  
+                                  surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,   &
+                                  surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,   &
+                                  surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,   &
+                                  surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,   &
+                                  surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,   &
+                                  surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
+
+!
+!--             Prognostic equation for cloud water content
+                DO  k = nzb+1, nzt
+                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d *                         &
+                                                      ( tsc(2) * tend(k,j,i) + &
+                                                        tsc(3) * tqc_m(k,j,i) )&
+                                                      - tsc(5) * rdf_sc(k)     &
+                                                               * qc(k,j,i)     &
+                                             )                                 &
+                                       * MERGE( 1.0_wp, 0.0_wp,                &
+                                                BTEST( wall_flags_0(k,j,i), 0 )&
+                                              ) 
+                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
+                ENDDO
+!
+!--             Calculate tendencies for the next Runge-Kutta step
+                IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( intermediate_timestep_count == 1 )  THEN
+                      DO  k = nzb+1, nzt
+                         tqc_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ELSEIF ( intermediate_timestep_count < &
+                            intermediate_timestep_count_max )  THEN
+                      DO  k = nzb+1, nzt
+                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
+                                           5.3125_wp * tqc_m(k,j,i)
+                      ENDDO
+                   ENDIF
+                ENDIF
+
+!
+!--             Calculate prognostic equation for cloud drop concentration.
+                tend(:,j,i) = 0.0_wp
+                IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( ws_scheme_sca )  THEN
+                      CALL advec_s_ws( i, j, nc, 'nc', flux_s_nc,    &
+                                    diss_s_nc, flux_l_nc, diss_l_nc, &
+                                    i_omp_start, tn )
+                   ELSE
+                      CALL advec_s_pw( i, j, nc )
+                   ENDIF
+                ELSE
+                   CALL advec_s_up( i, j, nc )
+                ENDIF
+                CALL diffusion_s( i, j, nc,                                    &
+                                  surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,    &
+                                  surf_def_h(2)%ncsws,                         &
+                                  surf_lsm_h%ncsws,    surf_usm_h%ncsws,       &
+                                  surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,    &
+                                  surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,    &
+                                  surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,    &
+                                  surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,    &
+                                  surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,    &
+                                  surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
+
+!
+!--             Prognostic equation for cloud drop concentration
+                DO  k = nzb+1, nzt
+                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d *                         &
+                                                      ( tsc(2) * tend(k,j,i) + &
+                                                        tsc(3) * tnc_m(k,j,i) )&
+                                                      - tsc(5) * rdf_sc(k)     &
+                                                               * nc(k,j,i)     &
+                                             )                                 &
+                                       * MERGE( 1.0_wp, 0.0_wp,                &
+                                                BTEST( wall_flags_0(k,j,i), 0 )&
+                                              )
+                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
+                ENDDO
+!
+!--             Calculate tendencies for the next Runge-Kutta step
+                IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( intermediate_timestep_count == 1 )  THEN
+                      DO  k = nzb+1, nzt
+                         tnc_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ELSEIF ( intermediate_timestep_count < &
+                            intermediate_timestep_count_max )  THEN
+                      DO  k = nzb+1, nzt
+                         tnc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
+                                           5.3125_wp * tnc_m(k,j,i)
+                      ENDDO
+                   ENDIF
+                ENDIF
+
+             ENDIF
+!
 !--          If required, calculate prognostic equations for rain water content
 !--          and rain drop concentration
@@ -1808 +1931 @@
 
 !
+!--    If required, calculate prognostic equations for cloud water content
+!--    and cloud drop concentration
+       IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+
+          CALL cpu_log( log_point(67), 'qc-equation', 'start' )
+
+!
+!--       Calculate prognostic equation for cloud water content
+          sbt = tsc(2)
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0_wp
+             ENDIF
+             tend = 0.0_wp
+             CALL advec_s_bc( qc, 'qc' )
+
+          ENDIF
+
+!
+!--       qc-tendency terms with no communication
+          IF ( scalar_advec /= 'bc-scheme' )  THEN
+             tend = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( qc, 'qc' )
+                ELSE
+                   CALL advec_s_pw( qc )
+                ENDIF
+             ELSE
+                CALL advec_s_up( qc )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( qc,                                                &
+                            surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,          &
+                            surf_def_h(2)%qcsws,                               &
+                            surf_lsm_h%qcsws,    surf_usm_h%qcsws,             &
+                            surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,          &
+                            surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,          &
+                            surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,          &
+                            surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,          &
+                            surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,          &
+                            surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
+
+!
+!--       Prognostic equation for cloud water content
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
+                                                      tsc(3) * tqc_m(k,j,i) )  &
+                                                    - tsc(5) * rdf_sc(k) *     &
+                                                               qc(k,j,i)       &
+                                             )                                 &
+                                    * MERGE( 1.0_wp, 0.0_wp,                   &
+                                             BTEST( wall_flags_0(k,j,i), 0 )   &
+                                          )
+                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb+1, nzt
+                         tqc_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb+1, nzt
+                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
+                                         + 5.3125_wp * tqc_m(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+
+          CALL cpu_log( log_point(67), 'qc-equation', 'stop' )
+          CALL cpu_log( log_point(68), 'nc-equation', 'start' )
+
+!
+!--       Calculate prognostic equation for cloud drop concentration
+          sbt = tsc(2)
+          IF ( scalar_advec == 'bc-scheme' )  THEN
+
+             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--             Bott-Chlond scheme always uses Euler time step. Thus:
+                sbt = 1.0_wp
+             ENDIF
+             tend = 0.0_wp
+             CALL advec_s_bc( nc, 'nc' )
+
+          ENDIF
+
+!
+!--       nc-tendency terms with no communication
+          IF ( scalar_advec /= 'bc-scheme' )  THEN
+             tend = 0.0_wp
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( ws_scheme_sca )  THEN
+                   CALL advec_s_ws( nc, 'nc' )
+                ELSE
+                   CALL advec_s_pw( nc )
+                ENDIF
+             ELSE
+                CALL advec_s_up( nc )
+             ENDIF
+          ENDIF
+
+          CALL diffusion_s( nc,                                                &
+                            surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,          &
+                            surf_def_h(2)%ncsws,                               &
+                            surf_lsm_h%ncsws,    surf_usm_h%ncsws,             & 
+                            surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,          &
+                            surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,          &
+                            surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,          &
+                            surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,          &
+                            surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,          &
+                            surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
+
+!
+!--       Prognostic equation for cloud drop concentration
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb+1, nzt
+                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
+                                                      tsc(3) * tnc_m(k,j,i) )  &
+                                                    - tsc(5) * rdf_sc(k) *     &
+                                                               nc(k,j,i)       &
+                                             )                                 &
+                                   * MERGE( 1.0_wp, 0.0_wp,                    &
+                                             BTEST( wall_flags_0(k,j,i), 0 )   &
+                                          )
+                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       Calculate tendencies for the next Runge-Kutta step
+          IF ( timestep_scheme(1:5) == 'runge' )  THEN
+             IF ( intermediate_timestep_count == 1 )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb+1, nzt
+                         tnc_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSEIF ( intermediate_timestep_count < &
+                      intermediate_timestep_count_max )  THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      DO  k = nzb+1, nzt
+                         tnc_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
+                                         + 5.3125_wp * tnc_m(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+          ENDIF
+
+          CALL cpu_log( log_point(68), 'nc-equation', 'stop' )
+
+       ENDIF
+!
 !--    If required, calculate prognostic equations for rain water content
 !--    and rain drop concentration

palm/trunk/SOURCE/read_3d_binary.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2269 2017-06-09 11:57:32Z suehring
 ! Enable restart runs for urban_surface_mod
 ! 
@@ -127 +132 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, kh, km, p, pt, q, ql, qc, nr, prr, precipitation_amount, qr, &
-               s, sa, u, u_m_l, u_m_n, u_m_r, u_m_s, v, v_m_l, v_m_n, v_m_r,   &
-               v_m_s, vpt, w, w_m_l, w_m_n, w_m_r, w_m_s
+        ONLY:  e, kh, km, p, pt, q, ql, qc, nc, nr, prr, precipitation_amount, &
+               qr, s, sa, u, u_m_l, u_m_n, u_m_r, u_m_s, v, v_m_l, v_m_n,      &
+               v_m_r, v_m_s, vpt, w, w_m_l, w_m_n, w_m_r, w_m_s
 
     USE averaging
@@ -496 +501 @@
                                   tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+                CASE ( 'nc' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   nc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'nc_av' )
+                   IF ( .NOT. ALLOCATED( nc_av ) )  THEN
+                      ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   nc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+
                 CASE ( 'nr' )
                    IF ( k == 1 )  READ ( 13 )  tmp_3d

palm/trunk/SOURCE/sum_up_3d_data.f90

@@ -25 +25 @@
 ! -----------------
 ! $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
@@ -135 +140 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  dzw, e, nr, p, pt, precipitation_rate, q, qc, ql, ql_c, ql_v,   &
-               qr, rho_ocean, s, sa, u, v, vpt, w
+        ONLY:  dzw, e, nc, nr, p, pt, precipitation_rate, q, qc, ql, ql_c,     &
+               ql_v, qr, rho_ocean, s, sa, u, v, vpt, w
 
     USE averaging,                                                             &
-        ONLY:  e_av, lpt_av, lwp_av, nr_av, ol_av, p_av, pc_av, pr_av, prr_av, &
-               precipitation_rate_av, pt_av, q_av, qc_av, ql_av, ql_c_av,      &
-               ql_v_av, ql_vp_av, qr_av, qsws_av, qv_av, rho_ocean_av, s_av, sa_av,  &
-               shf_av, ssws_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, &
-               z0h_av, z0q_av
+        ONLY:  e_av, lpt_av, lwp_av, nc_av, nr_av, ol_av, p_av, pc_av, pr_av,  &
+               prr_av, precipitation_rate_av, pt_av, q_av, qc_av, ql_av,       &
+               ql_c_av, ql_v_av, ql_vp_av, qr_av, qsws_av, qv_av, rho_ocean_av,&
+               s_av, sa_av, shf_av, ssws_av, ts_av, u_av, us_av, v_av, vpt_av, &
+               w_av, z0_av, z0h_av, z0q_av
 
     USE cloud_parameters,                                                      &
@@ -230 +235 @@
                 lwp_av = 0.0_wp
 
+             CASE ( 'nc' )
+                IF ( .NOT. ALLOCATED( nc_av ) )  THEN
+                   ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                nc_av = 0.0_wp
+
              CASE ( 'nr' )
                 IF ( .NOT. ALLOCATED( nr_av ) )  THEN
@@ -483 +494 @@
                    lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i)            &
                                                * dzw(1:nzt+1) ) * rho_surface
+                ENDDO
+             ENDDO
+
+          CASE ( 'nc' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      nc_av(k,j,i) = nc_av(k,j,i) + nc(k,j,i)
+                   ENDDO
                 ENDDO
              ENDDO

palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2281 2017-06-13 11:34:50Z suehring
 ! Clean-up unnecessary index access to surface type
 ! 
@@ -178 +183 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, kh, nr, pt, q, ql, qr, s, u, v, vpt, w, zu, zw, drho_air_zw, &
-               rho_air_zw
+        ONLY:  e, kh, nc, nr, pt, q, ql, qc, qr, s, u, v, vpt, w, zu, zw,      &
+               drho_air_zw, rho_air_zw
 
     USE cloud_parameters,                                                      &
@@ -195 +200 @@
                intermediate_timestep_count, intermediate_timestep_count_max,   &
                land_surface, large_scale_forcing, lsf_surf,                    &
-               message_string, microphysics_seifert, most_method, neutral,     &
-               passive_scalar, pt_surface, q_surface, run_coupled,             &
-               surface_pressure, simulated_time, terminate_run,                &
+               message_string, microphysics_morrison, microphysics_seifert,    &
+               most_method, neutral, passive_scalar, pt_surface, q_surface,    &
+               run_coupled, surface_pressure, simulated_time, terminate_run,   &
                urban_surface, zeta_max, zeta_min
 
@@ -1526 +1531 @@
 
 !
-!-- Calculate the other MOST scaling parameters theta*, q*, (qr*, nr*)
+!-- Calculate the other MOST scaling parameters theta*, q*, (qc*, qr*, nc*, nr*)
     SUBROUTINE calc_scaling_parameters
 
@@ -1749 +1754 @@
        ENDIF
 
+!
+!--    If required compute qc* and nc*
+       IF ( cloud_physics  .AND.  microphysics_morrison  .AND.                 &
+            .NOT. surf_vertical )  THEN
+          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
+          DO  m = 1, surf%ns    
+
+             i    = surf%i(m)            
+             j    = surf%j(m)
+             k    = surf%k(m)
+
+             z_mo = surf%z_mo(m)
+
+             surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) )                 &
+                                 / ( LOG( z_mo / surf%z0q(m) )                 &
+                                 - psi_h( z_mo / surf%ol(m) )                  &
+                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
+
+             surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) )                 &
+                                 / ( LOG( z_mo / surf%z0q(m) )                 &
+                                 - psi_h( z_mo / surf%ol(m) )                  &
+                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
+          ENDDO
+
+       ENDIF
 
 !
@@ -1781 +1811 @@
 
 !
-!-- Calculate surface fluxes usws, vsws, shf, qsws, (qrsws, nrsws)
+!-- Calculate surface fluxes usws, vsws, shf, qsws, (qcsws, qrsws, ncsws, nrsws)
     SUBROUTINE calc_surface_fluxes
 
@@ -1927 +1957 @@
 
              ENDDO
-          ENDIF       
+          ENDIF   
+!
+!--       Compute (turbulent) fluxes of cloud water content and cloud drop conc.
+          IF ( cloud_physics  .AND.  microphysics_morrison  .AND.              &
+               .NOT. downward)  THEN
+             !$OMP PARALLEL DO PRIVATE( i, j )
+             DO  m = 1, surf%ns   
+
+                i    = surf%i(m)            
+                j    = surf%j(m)
+
+                surf%qcsws(m) = -surf%qcs(m) * surf%us(m)
+                surf%ncsws(m) = -surf%ncs(m) * surf%us(m)
+             ENDDO
+          ENDIF    
 !
 !--       Compute (turbulent) fluxes of rain water content and rain drop conc.
-          IF ( cloud_physics  .AND.  microphysics_seifert  .AND.              &
+          IF ( cloud_physics  .AND.  microphysics_seifert  .AND.               &
                .NOT. downward)  THEN
              !$OMP PARALLEL DO PRIVATE( i, j )

palm/trunk/SOURCE/surface_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2270 2017-06-09 12:18:47Z maronga
 ! Parameters removed/added due to changes in the LSM
 ! 
@@ -104 +109 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  qs        !< scaling parameter humidity
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  ss        !< scaling parameter passive scalar
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qcs       !< scaling parameter qc
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ncs       !< scaling parameter nc
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  qrs       !< scaling parameter qr
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrs       !< scaling parameter nr
@@ -124 +131 @@
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws      !< surface flux latent heat 
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  ssws      !< surface flux passive scalar
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qcsws     !< surface flux qc
+       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ncsws     !< surface flux nc
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  qrsws     !< surface flux qr
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrsws     !< surface flux nr
@@ -719 +728 @@
 !
 !--       
+       IF ( cloud_physics .AND. microphysics_morrison)  THEN
+          ALLOCATE ( surfaces%qcs(1:surfaces%ns)   )
+          ALLOCATE ( surfaces%ncs(1:surfaces%ns)   )
+          ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
+          ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
+       ENDIF
+!
+!--       
        IF ( cloud_physics .AND. microphysics_seifert)  THEN
           ALLOCATE ( surfaces%qrs(1:surfaces%ns)   )
@@ -775 +792 @@
           ALLOCATE ( surfaces%ssws(1:surfaces%ns) ) 
        ENDIF 
+!
+!--       
+       IF ( cloud_physics .AND. microphysics_morrison)  THEN
+          ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
+          ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
+       ENDIF
 !
 !--       
@@ -877 +900 @@
           ALLOCATE ( surfaces%ssws(1:surfaces%ns) ) 
        ENDIF 
+
+       IF ( cloud_physics .AND. microphysics_seifert)  THEN
+          ALLOCATE ( surfaces%qcs(1:surfaces%ns)   )
+          ALLOCATE ( surfaces%ncs(1:surfaces%ns)   )
+          ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
+          ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
+       ENDIF
 
        IF ( cloud_physics .AND. microphysics_seifert)  THEN
@@ -1308 +1338 @@
              IF ( humidity )  THEN
                 surf%qs(num_h)   = 0.0_wp
+                IF ( cloud_physics .AND. microphysics_morrison)  THEN
+                   surf%qcs(num_h) = 0.0_wp
+                   surf%ncs(num_h) = 0.0_wp
+   
+                   surf%qcsws(num_h) = 0.0_wp
+                   surf%ncsws(num_h) = 0.0_wp
+
+                ENDIF
                 IF ( cloud_physics .AND. microphysics_seifert)  THEN
                    surf%qrs(num_h) = 0.0_wp
@@ -1447 +1485 @@
              IF ( humidity )  THEN
              surf%qsws = 0.0_wp
+                IF ( cloud_physics  .AND.  microphysics_morrison ) THEN
+                   surf%ncsws = 0.0_wp
+                   surf%qcsws = 0.0_wp
+                ENDIF
                 IF ( cloud_physics  .AND.  microphysics_seifert ) THEN
                    surf%nrsws = 0.0_wp
@@ -1568 +1610 @@
 !
 !--             Following wall fluxes are assumed to be zero 
+                IF ( cloud_physics .AND. microphysics_morrison)  THEN
+                   surf%qcs(num_v) = 0.0_wp
+                   surf%ncs(num_v) = 0.0_wp
+   
+                   surf%qcsws(num_v) = 0.0_wp
+                   surf%ncsws(num_v) = 0.0_wp
+                ENDIF
                 IF ( cloud_physics .AND. microphysics_seifert)  THEN
                    surf%qrs(num_v) = 0.0_wp
@@ -1656 +1705 @@
                    IF ( ALLOCATED( surf_def_h(l)%ss ) )                        &
                       surf_h(l)%ss(mm(l))      = surf_def_h(l)%ss(m)
+                   IF ( ALLOCATED( surf_def_h(l)%qcs ) )                       &
+                      surf_h(l)%qcs(mm(l))     = surf_def_h(l)%qcs(m)
+                   IF ( ALLOCATED( surf_def_h(l)%ncs ) )                       &
+                      surf_h(l)%ncs(mm(l))     = surf_def_h(l)%ncs(m)
                    IF ( ALLOCATED( surf_def_h(l)%qrs ) )                       &
                       surf_h(l)%qrs(mm(l))     = surf_def_h(l)%qrs(m)
@@ -1674 +1727 @@
                    IF ( ALLOCATED( surf_def_h(l)%ssws ) )                      &
                       surf_h(l)%qsws(mm(l))    = surf_def_h(l)%ssws(m)
+                   IF ( ALLOCATED( surf_def_h(l)%ncsws ) )                     &
+                      surf_h(l)%ncsws(mm(l))   = surf_def_h(l)%ncsws(m)
                    IF ( ALLOCATED( surf_def_h(l)%nrsws ) )                     &
                       surf_h(l)%nrsws(mm(l))   = surf_def_h(l)%nrsws(m)
@@ -1693 +1748 @@
                       IF ( ALLOCATED( surf_lsm_h%ss ) )                        &
                          surf_h(0)%ss(mm(0))      = surf_lsm_h%ss(m)
+                      IF ( ALLOCATED( surf_lsm_h%qcs ) )                       &
+                         surf_h(0)%qcs(mm(0))     = surf_lsm_h%qcs(m)
+                      IF ( ALLOCATED( surf_lsm_h%ncs ) )                       &
+                         surf_h(0)%ncs(mm(0))     = surf_lsm_h%ncs(m)
                       IF ( ALLOCATED( surf_lsm_h%qrs ) )                       &
                          surf_h(0)%qrs(mm(0))     = surf_lsm_h%qrs(m)
@@ -1711 +1770 @@
                       IF ( ALLOCATED( surf_lsm_h%ssws ) )                      &
                          surf_h(0)%qsws(mm(0))    = surf_lsm_h%ssws(m)
+                      IF ( ALLOCATED( surf_lsm_h%ncsws ) )                     &
+                         surf_h(0)%ncsws(mm(0))   = surf_lsm_h%ncsws(m)
                       IF ( ALLOCATED( surf_lsm_h%nrsws ) )                     &
                          surf_h(0)%nrsws(mm(0))   = surf_lsm_h%nrsws(m)
@@ -1730 +1791 @@
                       IF ( ALLOCATED( surf_usm_h%ss ) )                        &
                          surf_h(0)%ss(mm(0))      = surf_usm_h%ss(m)
+                      IF ( ALLOCATED( surf_usm_h%qcs ) )                       &
+                         surf_h(0)%qcs(mm(0))     = surf_usm_h%qcs(m)
+                      IF ( ALLOCATED( surf_usm_h%ncs ) )                       &
+                         surf_h(0)%ncs(mm(0))     = surf_usm_h%ncs(m)
                       IF ( ALLOCATED( surf_usm_h%qrs ) )                       &
                          surf_h(0)%qrs(mm(0))     = surf_usm_h%qrs(m)
@@ -1748 +1813 @@
                       IF ( ALLOCATED( surf_usm_h%ssws ) )                      &
                          surf_h(0)%qsws(mm(0))    = surf_usm_h%ssws(m)
+                      IF ( ALLOCATED( surf_usm_h%ncsws ) )                     &
+                         surf_h(0)%ncsws(mm(0))   = surf_usm_h%ncsws(m)
                       IF ( ALLOCATED( surf_usm_h%nrsws ) )                     &
                          surf_h(0)%nrsws(mm(0))   = surf_usm_h%nrsws(m)
@@ -1791 +1858 @@
                    IF ( ALLOCATED( surf_def_v(l)%ss ) )                        &
                       surf_v(l)%ss(mm(l))      = surf_def_v(l)%ss(m)
+                   IF ( ALLOCATED( surf_def_v(l)%qcs ) )                       &
+                      surf_v(l)%qcs(mm(l))     = surf_def_v(l)%qcs(m)
+                   IF ( ALLOCATED( surf_def_v(l)%ncs ) )                       &
+                      surf_v(l)%ncs(mm(l))     = surf_def_v(l)%ncs(m)
                    IF ( ALLOCATED( surf_def_v(l)%qrs ) )                       &
                       surf_v(l)%qrs(mm(l))     = surf_def_v(l)%qrs(m)
@@ -1805 +1876 @@
                    IF ( ALLOCATED( surf_def_v(l)%ssws ) )                      &
                       surf_v(l)%qsws(mm(l))    = surf_def_v(l)%ssws(m)
+                   IF ( ALLOCATED( surf_def_v(l)%ncsws ) )                     &
+                      surf_v(l)%ncsws(mm(l))   = surf_def_v(l)%ncsws(m)
                    IF ( ALLOCATED( surf_def_v(l)%nrsws ) )                     &
                       surf_v(l)%nrsws(mm(l))   = surf_def_v(l)%nrsws(m)
@@ -1829 +1902 @@
                    IF ( ALLOCATED( surf_lsm_v(l)%ss ) )                        &
                       surf_v(l)%ss(mm(l))      = surf_lsm_v(l)%ss(m)
+                   IF ( ALLOCATED( surf_lsm_v(l)%qcs ) )                       &
+                      surf_v(l)%qcs(mm(l))     = surf_lsm_v(l)%qcs(m)
+                   IF ( ALLOCATED( surf_lsm_v(l)%ncs ) )                       &
+                      surf_v(l)%ncs(mm(l))     = surf_lsm_v(l)%ncs(m)
                    IF ( ALLOCATED( surf_lsm_v(l)%qrs ) )                       &
                       surf_v(l)%qrs(mm(l))     = surf_lsm_v(l)%qrs(m)
@@ -1847 +1924 @@
                    IF ( ALLOCATED( surf_lsm_v(l)%ssws ) )                      &
                       surf_v(l)%qsws(mm(l))    = surf_lsm_v(l)%ssws(m)
+                   IF ( ALLOCATED( surf_lsm_v(l)%ncsws ) )                     &
+                      surf_v(l)%ncsws(mm(l))   = surf_lsm_v(l)%ncsws(m)
                    IF ( ALLOCATED( surf_lsm_v(l)%nrsws ) )                     &
                       surf_v(l)%nrsws(mm(l))   = surf_lsm_v(l)%nrsws(m)
@@ -1871 +1950 @@
                    IF ( ALLOCATED( surf_usm_v(l)%ss ) )                        &
                       surf_v(l)%ss(mm(l))      = surf_usm_v(l)%ss(m)
+                   IF ( ALLOCATED( surf_usm_v(l)%qcs ) )                       &
+                      surf_v(l)%qcs(mm(l))     = surf_usm_v(l)%qcs(m)
+                   IF ( ALLOCATED( surf_usm_v(l)%ncs ) )                       &
+                      surf_v(l)%ncs(mm(l))     = surf_usm_v(l)%ncs(m)
                    IF ( ALLOCATED( surf_usm_v(l)%qrs ) )                       &
                       surf_v(l)%qrs(mm(l))     = surf_usm_v(l)%qrs(m)
@@ -1889 +1972 @@
                    IF ( ALLOCATED( surf_usm_v(l)%ssws ) )                      &
                       surf_v(l)%qsws(mm(l))    = surf_usm_v(l)%ssws(m)
+                   IF ( ALLOCATED( surf_usm_v(l)%ncsws ) )                     &
+                      surf_v(l)%ncsws(mm(l))   = surf_usm_v(l)%ncsws(m)
                    IF ( ALLOCATED( surf_usm_v(l)%nrsws ) )                     &
                       surf_v(l)%nrsws(mm(l))   = surf_usm_v(l)%nrsws(m)
@@ -1946 +2031 @@
              WRITE ( 14 )  surf_h(l)%ss
           ENDIF
+          WRITE ( 14 )  'surf_h(' // dum // ')%qcs                 '
+          IF ( ALLOCATED ( surf_h(l)%qcs ) )  THEN  
+             WRITE ( 14 )  surf_h(l)%qcs
+          ENDIF
+          WRITE ( 14 )  'surf_h(' // dum // ')%ncs                 '  
+          IF ( ALLOCATED ( surf_h(l)%ncs ) )  THEN
+             WRITE ( 14 )  surf_h(l)%ncs
+          ENDIF
           WRITE ( 14 )  'surf_h(' // dum // ')%qrs                 '
           IF ( ALLOCATED ( surf_h(l)%qrs ) )  THEN  
@@ -1981 +2074 @@
           IF ( ALLOCATED ( surf_h(l)%ssws ) )  THEN 
              WRITE ( 14 )  surf_h(l)%ssws
+          ENDIF
+          WRITE ( 14 )  'surf_h(' // dum // ')%qcsws               '  
+          IF ( ALLOCATED ( surf_h(l)%qcsws ) )  THEN
+             WRITE ( 14 )  surf_h(l)%qcsws
+          ENDIF
+          WRITE ( 14 )  'surf_h(' // dum // ')%ncsws               '  
+          IF ( ALLOCATED ( surf_h(l)%ncsws ) )  THEN
+             WRITE ( 14 )  surf_h(l)%ncsws
           ENDIF
           WRITE ( 14 )  'surf_h(' // dum // ')%qrsws               '  
@@ -2021 +2122 @@
              WRITE ( 14 )  surf_v(l)%ss
           ENDIF
+          WRITE ( 14 )  'surf_v(' // dum // ')%qcs                 '  
+          IF ( ALLOCATED ( surf_v(l)%qcs ) )  THEN
+             WRITE ( 14 )  surf_v(l)%qcs
+          ENDIF
+          WRITE ( 14 )  'surf_v(' // dum // ')%ncs                 ' 
+          IF ( ALLOCATED ( surf_v(l)%ncs ) )  THEN 
+             WRITE ( 14 )  surf_v(l)%ncs
+          ENDIF
           WRITE ( 14 )  'surf_v(' // dum // ')%qrs                 '  
           IF ( ALLOCATED ( surf_v(l)%qrs ) )  THEN
@@ -2048 +2157 @@
           IF ( ALLOCATED ( surf_v(l)%ssws ) )  THEN
              WRITE ( 14 )  surf_v(l)%ssws
+          ENDIF
+          WRITE ( 14 )  'surf_v(' // dum // ')%qcsws               '  
+          IF ( ALLOCATED ( surf_v(l)%qcsws ) )  THEN
+             WRITE ( 14 )  surf_v(l)%qcsws
+          ENDIF
+          WRITE ( 14 )  'surf_v(' // dum // ')%ncsws               '  
+          IF ( ALLOCATED ( surf_v(l)%ncsws ) )  THEN
+             WRITE ( 14 )  surf_v(l)%ncsws
           ENDIF
           WRITE ( 14 )  'surf_v(' // dum // ')%qrsws               '  
@@ -2240 +2357 @@
                       IF ( ALLOCATED( surf_h(0)%ss )  .AND.  kk == 1 )         &
                          READ ( 13 )  surf_h(0)%ss
+                   CASE ( 'surf_h(0)%qcs' )         
+                      IF ( ALLOCATED( surf_h(0)%qcs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_h(0)%qcs
+                   CASE ( 'surf_h(0)%ncs' )         
+                      IF ( ALLOCATED( surf_h(0)%ncs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_h(0)%ncs
                    CASE ( 'surf_h(0)%qrs' )         
                       IF ( ALLOCATED( surf_h(0)%qrs )  .AND.  kk == 1 )        &
@@ -2267 +2390 @@
                       IF ( ALLOCATED( surf_h(0)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_h(0)%ssws
+                   CASE ( 'surf_h(0)%qcsws' )         
+                      IF ( ALLOCATED( surf_h(0)%qcsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_h(0)%qcsws
+                   CASE ( 'surf_h(0)%ncsws' )         
+                      IF ( ALLOCATED( surf_h(0)%ncsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_h(0)%ncsws
                    CASE ( 'surf_h(0)%qrsws' )         
                       IF ( ALLOCATED( surf_h(0)%qrsws )  .AND.  kk == 1 )      &
@@ -2296 +2425 @@
                       IF ( ALLOCATED( surf_h(1)%ss )  .AND.  kk == 1 )         &
                          READ ( 13 )  surf_h(1)%ss
+                   CASE ( 'surf_h(1)%qcs' )         
+                      IF ( ALLOCATED( surf_h(1)%qcs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_h(1)%qcs
+                   CASE ( 'surf_h(1)%ncs' )         
+                      IF ( ALLOCATED( surf_h(1)%ncs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_h(1)%ncs
                    CASE ( 'surf_h(1)%qrs' )         
                       IF ( ALLOCATED( surf_h(1)%qrs )  .AND.  kk == 1 )        &
@@ -2323 +2458 @@
                       IF ( ALLOCATED( surf_h(1)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_h(1)%ssws
+                   CASE ( 'surf_h(1)%qcsws' )         
+                      IF ( ALLOCATED( surf_h(1)%qcsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_h(1)%qcsws
+                   CASE ( 'surf_h(1)%ncsws' )         
+                      IF ( ALLOCATED( surf_h(1)%ncsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_h(1)%ncsws
                    CASE ( 'surf_h(1)%qrsws' )         
                       IF ( ALLOCATED( surf_h(1)%qrsws )  .AND.  kk == 1 )      &
@@ -2352 +2493 @@
                       IF ( ALLOCATED( surf_h(2)%ss )  .AND.  kk == 1 )         &
                          READ ( 13 )  surf_h(2)%ss
+                   CASE ( 'surf_h(2)%qcs' )         
+                      IF ( ALLOCATED( surf_h(2)%qcs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_h(2)%qcs
+                   CASE ( 'surf_h(2)%ncs' )         
+                      IF ( ALLOCATED( surf_h(2)%ncs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_h(2)%ncs
                    CASE ( 'surf_h(2)%qrs' )         
                       IF ( ALLOCATED( surf_h(2)%qrs )  .AND.  kk == 1 )        &
@@ -2379 +2526 @@
                       IF ( ALLOCATED( surf_h(2)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_h(2)%ssws
+                   CASE ( 'surf_h(2)%qcsws' )         
+                      IF ( ALLOCATED( surf_h(2)%qcsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_h(2)%qcsws
+                   CASE ( 'surf_h(2)%ncsws' )         
+                      IF ( ALLOCATED( surf_h(2)%ncsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_h(2)%ncsws
                    CASE ( 'surf_h(2)%qrsws' )         
                       IF ( ALLOCATED( surf_h(2)%qrsws )  .AND.  kk == 1 )      &
@@ -2410 +2563 @@
                       IF ( ALLOCATED( surf_v(0)%ss )  .AND.  kk == 1 )         &
                          READ ( 13 )  surf_v(0)%ss
+                   CASE ( 'surf_v(0)%qcs' )         
+                      IF ( ALLOCATED( surf_v(0)%qcs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_v(0)%qcs
+                   CASE ( 'surf_v(0)%ncs' )         
+                      IF ( ALLOCATED( surf_v(0)%ncs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_v(0)%ncs
                    CASE ( 'surf_v(0)%qrs' )         
                       IF ( ALLOCATED( surf_v(0)%qrs )  .AND.  kk == 1 )        &
@@ -2431 +2590 @@
                       IF ( ALLOCATED( surf_v(0)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_v(0)%ssws
+                   CASE ( 'surf_v(0)%qcsws' )         
+                      IF ( ALLOCATED( surf_v(0)%qcsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_v(0)%qcsws
+                   CASE ( 'surf_v(0)%ncsws' )         
+                      IF ( ALLOCATED( surf_v(0)%ncsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_v(0)%ncsws
                    CASE ( 'surf_v(0)%qrsws' )         
                       IF ( ALLOCATED( surf_v(0)%qrsws )  .AND.  kk == 1 )      &
@@ -2469 +2634 @@
                       IF ( ALLOCATED( surf_v(1)%ss )  .AND.  kk == 1 )         &
                          READ ( 13 )  surf_v(1)%ss
+                   CASE ( 'surf_v(1)%qcs' )         
+                      IF ( ALLOCATED( surf_v(1)%qcs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_v(1)%qcs
+                   CASE ( 'surf_v(1)%ncs' )         
+                      IF ( ALLOCATED( surf_v(1)%ncs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_v(1)%ncs
                    CASE ( 'surf_v(1)%qrs' )         
                       IF ( ALLOCATED( surf_v(1)%qrs )  .AND.  kk == 1 )        &
@@ -2490 +2661 @@
                       IF ( ALLOCATED( surf_v(1)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_v(1)%ssws
+                   CASE ( 'surf_v(1)%qcsws' )         
+                      IF ( ALLOCATED( surf_v(1)%qcsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_v(1)%qcsws
+                   CASE ( 'surf_v(1)%ncsws' )         
+                      IF ( ALLOCATED( surf_v(1)%ncsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_v(1)%ncsws
                    CASE ( 'surf_v(1)%qrsws' )         
                       IF ( ALLOCATED( surf_v(1)%qrsws )  .AND.  kk == 1 )      &
@@ -2528 +2705 @@
                       IF ( ALLOCATED( surf_v(2)%ss )  .AND.  kk == 1 )         &
                          READ ( 13 )  surf_v(2)%ss
+                   CASE ( 'surf_v(2)%qcs' )         
+                      IF ( ALLOCATED( surf_v(2)%qcs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_v(2)%qcs
+                   CASE ( 'surf_v(2)%ncs' )         
+                      IF ( ALLOCATED( surf_v(2)%ncs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_v(2)%ncs
                    CASE ( 'surf_v(2)%qrs' )         
                       IF ( ALLOCATED( surf_v(2)%qrs )  .AND.  kk == 1 )        &
@@ -2549 +2732 @@
                       IF ( ALLOCATED( surf_v(2)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_v(2)%ssws
+                   CASE ( 'surf_v(2)%qcsws' )         
+                      IF ( ALLOCATED( surf_v(2)%qcsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_v(2)%qcsws
+                   CASE ( 'surf_v(2)%ncsws' )         
+                      IF ( ALLOCATED( surf_v(2)%ncsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_v(2)%ncsws
                    CASE ( 'surf_v(2)%qrsws' )         
                       IF ( ALLOCATED( surf_v(2)%qrsws )  .AND.  kk == 1 )      &
@@ -2587 +2776 @@
                       IF ( ALLOCATED( surf_v(3)%ss )  .AND.  kk == 1 )         &
                          READ ( 13 )  surf_v(3)%ss
+                   CASE ( 'surf_v(3)%qcs' )         
+                      IF ( ALLOCATED( surf_v(3)%qcs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_v(3)%qcs
+                   CASE ( 'surf_v(3)%ncs' )         
+                      IF ( ALLOCATED( surf_v(3)%ncs )  .AND.  kk == 1 )        &
+                         READ ( 13 )  surf_v(3)%ncs
                    CASE ( 'surf_v(3)%qrs' )         
                       IF ( ALLOCATED( surf_v(3)%qrs )  .AND.  kk == 1 )        &
@@ -2608 +2803 @@
                       IF ( ALLOCATED( surf_v(3)%ssws )  .AND.  kk == 1 )       &
                          READ ( 13 )  surf_v(3)%ssws
+                   CASE ( 'surf_v(3)%qcsws' )         
+                      IF ( ALLOCATED( surf_v(3)%qcsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_v(3)%qcsws
+                   CASE ( 'surf_v(3)%ncsws' )         
+                      IF ( ALLOCATED( surf_v(3)%ncsws )  .AND.  kk == 1 )      &
+                         READ ( 13 )  surf_v(3)%ncsws
                    CASE ( 'surf_v(3)%qrsws' )         
                       IF ( ALLOCATED( surf_v(3)%qrsws )  .AND.  kk == 1 )      &
@@ -2814 +3015 @@
              ENDIF
 
+             IF ( SCAN( TRIM( field_chr ), '%qcs' ) /= 0 )  THEN 
+                IF ( ALLOCATED( surf_target%qcs )  .AND.                       &
+                     ALLOCATED( surf_file%qcs   ) )                            & 
+                  surf_target%qcs(m_target) = surf_file%qcs(m_file)
+             ENDIF
+
+             IF ( SCAN( TRIM( field_chr ), '%qcsws' ) /= 0 )  THEN 
+                IF ( ALLOCATED( surf_target%qcsws )  .AND.                     &
+                     ALLOCATED( surf_file%qcsws   ) )                          & 
+                   surf_target%qcsws(m_target) = surf_file%qcsws(m_file)
+             ENDIF
+
+             IF ( SCAN( TRIM( field_chr ), '%ncs' ) /= 0 )  THEN 
+                IF ( ALLOCATED( surf_target%ncs )  .AND.                       &
+                     ALLOCATED( surf_file%ncs   ) )                            & 
+                   surf_target%ncs(m_target) = surf_file%ncs(m_file)
+             ENDIF
+
+             IF ( SCAN( TRIM( field_chr ), '%ncsws' ) /= 0 )  THEN 
+                IF ( ALLOCATED( surf_target%ncsws )  .AND.                     &
+                     ALLOCATED( surf_file%ncsws   ) )                          & 
+                   surf_target%ncsws(m_target) = surf_file%ncsws(m_file)
+             ENDIF
+
              IF ( SCAN( TRIM( field_chr ), '%qrs' ) /= 0 )  THEN 
                 IF ( ALLOCATED( surf_target%qrs )  .AND.                       &

palm/trunk/SOURCE/swap_timelevel.f90

@@ -25 +25 @@
 ! -----------------
 ! $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
@@ -103 +108 @@
 #if defined( __nopointer )
     USE arrays_3d,                                                             &
-        ONLY:  e, e_p, nr, nr_p, pt, pt_p, q, q_p, qr, qr_p, s, s_p, sa, sa_p, &
-               u, u_p, v, v_p, w, w_p
+        ONLY:  e, e_p, nc, nc_p, nr, nr_p, pt, pt_p, q, q_p, qc, qc_p qr, qr_p,&
+               s, s_p, sa, sa_p, u, u_p, v, v_p, w, w_p
 #else
     USE arrays_3d,                                                             &
-        ONLY:  e, e_1, e_2, e_p, nr, nr_1, nr_2, nr_p, pt, pt_1, pt_2, pt_p, q,&
-               q_1, q_2, q_p, qr, qr_1, qr_2, qr_p, s, s_1, s_2, s_p, sa, sa_1,&
-               sa_2, sa_p, u, u_1, u_2, u_p, v, v_1, v_2, v_p, w, w_1, w_2, w_p
+        ONLY:  e, e_1, e_2, e_p, nc, nc_1, nc_2, nc_p, nr, nr_1, nr_2, nr_p,   &
+               pt, pt_1, pt_2, pt_p, q, q_1, q_2, q_p, qc, qc_1, qc_2, qc_p,   &
+               qr, qr_1, qr_2, qr_p, s, s_1, s_2, s_p, sa, sa_1, sa_2, sa_p,   &
+               u, u_1, u_2, u_p, v, v_1, v_2, v_p, w, w_1, w_2, w_p
 
 #endif
@@ -121 +127 @@
     USE control_parameters,                                                    &
         ONLY:  cloud_physics, constant_diffusion, humidity, land_surface,      &
-               microphysics_seifert, neutral, ocean, passive_scalar,           &
-               timestep_count, urban_surface
+               microphysics_morrison, microphysics_seifert, neutral, ocean,    &
+               passive_scalar, timestep_count, urban_surface
 
     USE indices,                                                               &
@@ -174 +180 @@
 
     IF ( humidity )  THEN
-       q = q_p             
+       q = q_p
+       IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+          qc = qc_p
+          nc = nc_p
+       ENDIF             
        IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
           qr = qr_p
@@ -214 +224 @@
           IF ( humidity )  THEN
              q => q_1;    q_p => q_2
+             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+                qc => qc_1;    qc_p => qc_2
+                nc => nc_1;    nc_p => nc_2
+             ENDIF
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
                 qr => qr_1;    qr_p => qr_2
@@ -241 +255 @@
           IF ( humidity )  THEN
              q => q_2;    q_p => q_1
+             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+                qc => qc_2;    qc_p => qc_1
+                nc => nc_2;    nc_p => nc_1
+             ENDIF
              IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
                 qr => qr_2;    qr_p => qr_1

palm/trunk/SOURCE/time_integration.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
+! includes two more prognostic equations for cloud drop concentration (nc)  
+! and cloud water content (qc). 
+! 
+! 2271 2017-06-09 12:34:55Z sward
 ! Start timestep message changed
 ! 
@@ -267 +272 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  diss, dzu, e, e_p, nr, nr_p, prho, pt, pt_p, pt_init, q_init, q,&
-               ql, ql_c, ql_v, ql_vp, qr, qr_p, q_p, ref_state, rho_ocean,     &
-               s, s_p, sa_p, tend, u, u_p, v, vpt, v_p, w, w_p
+        ONLY:  diss, dzu, e, e_p, nc, nc_p, nr, nr_p, prho, pt, pt_p, pt_init, &
+               q_init, q, qc, qc_p, ql, ql_c, ql_v, ql_vp, qr, qr_p, q_p,      &
+               ref_state, rho_ocean, s, s_p, sa_p, tend, u, u_p, v, vpt,       &
+               v_p, w, w_p
 
     USE calc_mean_profile_mod,                                                 &
@@ -291 +297 @@
                land_surface, large_scale_forcing,                              &
                loop_optimization, lsf_surf, lsf_vert, masks,                   &
-               microphysics_seifert, mid, nest_domain,                         &
+               microphysics_morrison, microphysics_seifert, mid, nest_domain,  &
                neutral, nr_timesteps_this_run, nudging,                        &
                ocean, passive_scalar,                                          &
@@ -585 +591 @@
           IF ( humidity )  THEN
              CALL exchange_horiz( q_p, nbgp )
+             IF ( cloud_physics .AND. microphysics_morrison )  THEN
+                CALL exchange_horiz( qc_p, nbgp )
+                CALL exchange_horiz( nc_p, nbgp )
+             ENDIF
              IF ( cloud_physics .AND. microphysics_seifert )  THEN
                 CALL exchange_horiz( qr_p, nbgp )
@@ -637 +647 @@
                    CALL exchange_horiz( q, nbgp )
 
+                   IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
+                       CALL exchange_horiz( qc, nbgp )
+                       CALL exchange_horiz( nc, nbgp )
+                   ENDIF
                    IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
-!                        CALL exchange_horiz( qc, nbgp )
                        CALL exchange_horiz( qr, nbgp )
-!                        CALL exchange_horiz( nc, nbgp )
                        CALL exchange_horiz( nr, nbgp )
                    ENDIF

palm/trunk/SOURCE/write_3d_binary.f90

@@ -25 +25 @@
 ! -----------------
 ! $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
@@ -116 +121 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, kh, km, p, pt, q, ql, qc, nr, prr, precipitation_amount, qr, &
-               s, sa, u, u_m_l, u_m_n, u_m_r, u_m_s, v, v_m_l, v_m_n, v_m_r,   &
-               v_m_s, vpt, w, w_m_l, w_m_n, w_m_r, w_m_s
+        ONLY:  e, kh, km, p, pt, q, ql, qc, nc, nr, prr, precipitation_amount, &
+               qr, s, sa, u, u_m_l, u_m_n, u_m_r, u_m_s, v, v_m_l, v_m_n,      &
+               v_m_r, v_m_s, vpt, w, w_m_l, w_m_n, w_m_r, w_m_s
         
     USE averaging
@@ -124 +129 @@
     USE control_parameters,                                                    &
         ONLY:  iran, humidity, passive_scalar, cloud_physics, cloud_droplets,  &
-               microphysics_seifert, ocean, topography
+               microphysics_morrison, microphysics_seifert, ocean, topography
                 
     USE indices,                                                               &
@@ -222 +227 @@
              IF ( ALLOCATED( qc_av ) )  THEN
                 WRITE ( 14 )  'qc_av               ';  WRITE ( 14 )  qc_av
+             ENDIF
+             IF ( microphysics_morrison )  THEN
+                WRITE ( 14 )  'nc                  ';  WRITE ( 14 )  nc
+                IF ( ALLOCATED( nc_av ) )  THEN
+                   WRITE ( 14 )  'nc_av               ';  WRITE ( 14 )  nc_av
+                ENDIF
              ENDIF
              IF ( microphysics_seifert )  THEN