Changeset 1053 for palm

Timestamp:

Nov 13, 2012 5:11:03 PM (12 years ago)

Author:

hoffmann

Message:

two-moment cloud physics implemented

Location:

palm/trunk

Files:

• . (modified) (1 prop)
• EXAMPLES/dvr (modified) (1 prop)
• SOURCE (modified) (1 prop)
• SOURCE/Makefile (modified) (5 diffs)
• SOURCE/advec_ws.f90 (modified) (6 diffs)
• SOURCE/average_3d_data.f90 (modified) (3 diffs)
• SOURCE/boundary_conds.f90 (modified) (3 diffs)
• SOURCE/calc_liquid_water_content.f90 (modified) (2 diffs)
• SOURCE/check_parameters.f90 (modified) (12 diffs)
• SOURCE/data_output_2d.f90 (modified) (5 diffs)
• SOURCE/data_output_3d.f90 (modified) (5 diffs)
• SOURCE/flow_statistics.f90 (modified) (4 diffs)
• SOURCE/init_3d_model.f90 (modified) (14 diffs)
• SOURCE/init_cloud_physics.f90 (modified) (4 diffs)
• SOURCE/microphysics.f90 (copied) (copied from palm/branches/hoffmann/SOURCE/microphysics.f90)
• SOURCE/modules.f90 (modified) (16 diffs)
• SOURCE/netcdf.f90 (modified) (7 diffs)
• SOURCE/parin.f90 (modified) (7 diffs)
• SOURCE/prognostic_equations.f90 (modified) (6 diffs)
• SOURCE/read_3d_binary.f90 (modified) (4 diffs)
• SOURCE/read_var_list.f90 (modified) (12 diffs)
• SOURCE/sum_up_3d_data.f90 (modified) (6 diffs)
• SOURCE/swap_timelevel.f90 (modified) (4 diffs)
• SOURCE/time_integration.f90 (modified) (2 diffs)
• SOURCE/timestep.f90 (modified) (3 diffs)
• SOURCE/user_actions.f90 (modified) (3 diffs)
• SOURCE/write_3d_binary.f90 (modified) (4 diffs)
• SOURCE/write_var_list.f90 (modified) (11 diffs)

palm/trunk/SOURCE/Makefile

@@ -20 +20 @@
 # Current revisions:
 # ------------------
-# 
+# +microphysics for two-moment cloud physics
 #
 # Former revisions:
@@ -159 +159 @@
         user_parin.f90 user_read_restart_data.f90 \
         user_spectra.f90 user_statistics.f90 wall_fluxes.f90 \
-        write_3d_binary.f90 write_compressed.f90 write_var_list.f90
+        write_3d_binary.f90 write_compressed.f90 write_var_list.f90 \
+        microphysics.f90
 
 OBJS =  advec_s_bc.o advec_s_pw.o advec_s_up.o advec_u_pw.o advec_u_up.o \
@@ -205 +206 @@
         user_lpm_init.o user_lpm_set_attributes.o user_module.o user_parin.o \
         user_read_restart_data.o user_spectra.o user_statistics.o \
-        wall_fluxes.o write_3d_binary.o write_compressed.o write_var_list.o
+        wall_fluxes.o write_3d_binary.o write_compressed.o write_var_list.o \
+        microphysics.o
 
 CC = cc
@@ -346 +348 @@
         eqn_state_seawater.o impact_of_latent_heat.o \
         plant_canopy_model.o production_e.o \
-        subsidence.o user_actions.o
+        subsidence.o user_actions.o \
+        microphysics.o
 random_gauss.o: random_function.o
 read_3d_binary.o: modules.o random_function.o
@@ -394 +397 @@
 write_compressed.o: modules.o
 write_var_list.o: modules.o
+microphysics.o: modules.o

palm/trunk/SOURCE/advec_ws.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme:
+! +flux_l_*, flux_s_*, diss_l_*, diss_s_*, sums_ws*s_ws_l 
 !
 ! Former revisions:
@@ -210 +213 @@
           ENDIF
 
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
+             sums_wsqrs_ws_l = 0.0
+             sums_wsnrs_ws_l = 0.0
+          ENDIF
+
           IF ( ocean )  THEN
              ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
@@ -253 +263 @@
 
              IF ( humidity .OR. passive_scalar )  THEN
-                ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1),         &
+                ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1),          &
                           diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
-                ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
+                ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1),  &
                           diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
              ENDIF
+
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1),         &
+                          diss_s_qr(nzb+1:nzt,0:threads_per_task-1),         &
+                          flux_s_nr(nzb+1:nzt,0:threads_per_task-1),         &
+                          diss_s_nr(nzb+1:nzt,0:threads_per_task-1) )
+                ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
+                          diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
+                          flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
+                          diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 
+             END IF
 
              IF ( ocean )  THEN
@@ -297 +318 @@
           sums_wspts_ws_l = 0.0
           IF ( humidity .OR. passive_scalar )  sums_wsqs_ws_l = 0.0
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             sums_wsqrs_ws_l = 0.0
+             sums_wsnrs_ws_l = 0.0
+          ENDIF
           IF ( ocean )  sums_wssas_ws_l = 0.0
 
@@ -307 +332 @@
 ! Scalar advection - Call for grid point i,j
 !------------------------------------------------------------------------------!
-    SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local,  &
-                              swap_diss_y_local, swap_flux_x_local, &
+    SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local,  &
+                              swap_diss_y_local, swap_flux_x_local,  &
                               swap_diss_x_local, i_omp, tn )
 
@@ -727 +752 @@
                                  * 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) + diss_t(k) )                &
+                                 * weight_substep(intermediate_timestep_count)
+             ENDDO
+
+          CASE ( 'nr' )
+
+             DO  k = nzb, nzt
+                sums_wsnrs_ws_l(k,tn)  = sums_wsnrs_ws_l(k,tn) +               &
+                                      ( flux_t(k) + diss_t(k) )                &
+                                 * weight_substep(intermediate_timestep_count)
+             ENDDO
+
 
          END SELECT

palm/trunk/SOURCE/average_3d_data.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! averaging of nr, qr added
 !
 ! Former revisions:
@@ -117 +117 @@
              ENDDO
 
+          CASE ( 'nr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      nr_av(k,j,i) = nr_av(k,j,i) / REAL( average_count_3d )
+                   ENDDO
+                ENDDO
+             ENDDO
+
           CASE ( 'p' )
              DO  i = nxlg, nxrg
@@ -207 +216 @@
              ENDDO
 
+          CASE ( 'qr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      qr_av(k,j,i) = qr_av(k,j,i) / REAL( average_count_3d )
+                   ENDDO
+                ENDDO
+             ENDDO
+
           CASE ( 'qv' )
              DO  i = nxlg, nxrg

palm/trunk/SOURCE/boundary_conds.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! boundary conditions for the two new prognostic equations (nr, qr) of the 
+! two-moment cloud scheme
 !  
 ! Former revisions:
@@ -216 +217 @@
 !--       Top boundary
           q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
+
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+!             
+!--          Surface conditions for constant_humidity_flux
+             IF ( ibc_qr_b == 0 ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qr_p(nzb_s_inner(j,i),j,i) = qr(nzb_s_inner(j,i),j,i)
+                   ENDDO
+                ENDDO
+             ELSE
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qr_p(nzb_s_inner(j,i),j,i) = qr_p(nzb_s_inner(j,i)+1,j,i)
+                   ENDDO
+                ENDDO
+             ENDIF
+!
+!--          Top boundary
+             qr_p(nzt+1,:,:) = qr_p(nzt,:,:) + bc_qr_t_val * dzu(nzt+1)
+!             
+!--          Surface conditions for constant_humidity_flux
+             IF ( ibc_nr_b == 0 ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      nr_p(nzb_s_inner(j,i),j,i) = nr(nzb_s_inner(j,i),j,i)
+                   ENDDO
+                ENDDO
+             ELSE
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      nr_p(nzb_s_inner(j,i),j,i) = nr_p(nzb_s_inner(j,i)+1,j,i)
+                   ENDDO
+                ENDDO
+             ENDIF
+!
+!--          Top boundary
+             nr_p(nzt+1,:,:) = nr_p(nzt,:,:) + bc_nr_t_val * dzu(nzt+1)
+          ENDIF
 
        ENDIF
@@ -242 +282 @@
           pt_p(:,nys-1,:)     = pt_p(:,nys,:)
           IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
-          IF ( humidity .OR. passive_scalar )  q_p(:,nys-1,:) = q_p(:,nys,:)
+          IF ( humidity .OR. passive_scalar )  THEN
+             q_p(:,nys-1,:) = q_p(:,nys,:)
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                qr_p(:,nys-1,:) = qr_p(:,nys,:)
+                nr_p(:,nys-1,:) = nr_p(:,nys,:)
+             ENDIF
+          ENDIF
        ELSEIF ( outflow_n )  THEN
           pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
           IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
-          IF ( humidity .OR. passive_scalar )  q_p(:,nyn+1,:) = q_p(:,nyn,:)
+          IF ( humidity .OR. passive_scalar )  THEN
+             q_p(:,nyn+1,:) = q_p(:,nyn,:)
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
+                nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
+             ENDIF
+          ENDIF
        ELSEIF ( outflow_l )  THEN
           pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
           IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
-          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxl-1) = q_p(:,:,nxl)
+          IF ( humidity .OR. passive_scalar )  THEN
+             q_p(:,:,nxl-1) = q_p(:,:,nxl)
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
+                nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
+             ENDIF
+          ENDIF
        ELSEIF ( outflow_r )  THEN
           pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
           IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
-          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxr+1) = q_p(:,:,nxr)
+          IF ( humidity .OR. passive_scalar )  THEN
+             q_p(:,:,nxr+1) = q_p(:,:,nxr)
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
+                nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
+             ENDIF
+          ENDIF
        ENDIF
 

palm/trunk/SOURCE/calc_liquid_water_content.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
+! description expanded to the two-moment cloud scheme
 !
 ! Former revisions:
@@ -51 +52 @@
 ! Description:
 ! ------------
-! Calculation of the liquid water content (0%-or-100%-scheme)
+! Calculation of the liquid water content (0%-or-100%-scheme). This scheme is 
+! used by the one and the two moment cloud physics scheme. Using the two moment 
+! scheme, this calculation results in the cloud water content.
 !------------------------------------------------------------------------------!
 

palm/trunk/SOURCE/check_parameters.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary changes for the new two-moment cloud physics scheme added:
+! - check cloud physics scheme (Kessler or Seifert and Beheng)
+! - plant_canopy is not allowed
+! - currently, only cache loop_optimization is allowed 
+! - initial profiles of nr, qr
+! - boundary condition of nr, qr
+! - check output quantities (qr, nr, prr)
 !
 ! Former revisions:
@@ -617 +623 @@
 
     ENDIF
-
+!
+!-- Check cloud scheme
+    IF ( cloud_scheme == 'seifert_beheng' )  THEN
+       icloud_scheme = 0
+    ELSEIF ( cloud_scheme == 'kessler' )  THEN
+       icloud_scheme = 1
+    ELSE
+       message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
+                        TRIM( cloud_scheme ) // '"'
+       CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
+    ENDIF
 !
 !-- Check whether there are any illegal values
@@ -853 +869 @@
     ENDIF 
 
+    IF ( plant_canopy  .AND.  cloud_physics  .AND.  icloud_scheme == 0 ) THEN
+       message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' //  &
+                        ' seifert_beheng'
+       CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
+    ENDIF 
+
+    IF ( loop_optimization /= 'cache' .AND.  cloud_physics  .AND.            &
+         icloud_scheme == 0 ) THEN
+       message_string = 'cloud_scheme = seifert_beheng requires ' // &
+                        'loop_optimization = cache'
+       CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
+    ENDIF 
+
 !
 !-- In case of no model continuation run, check initialising parameters and
@@ -861 +890 @@
 !--    Initial profiles for 1D and 3D model, respectively (u,v further below)
        pt_init = pt_surface
-       IF ( humidity )        q_init  = q_surface
+       IF ( humidity )  THEN
+          q_init  = q_surface
+!
+!--       It is not allowed to choose initial profiles of rain water content 
+!--       and rain drop concentration. They are set to 0.0.
+          IF ( cloud_physics .AND. icloud_scheme == 0 ) THEN
+             qr_init = 0.0
+             nr_init = 0.0
+          ENDIF
+       ENDIF
        IF ( ocean )           sa_init = sa_surface
        IF ( passive_scalar )  q_init  = s_surface
@@ -1167 +1205 @@
           ENDIF
 
-!
-!--       Store humidity gradient at the top boundary for possile Neumann
-!--       boundary condition
-          bc_q_t_val = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+
+             i = 1
+             gradient = 0.0
+             qr_vertical_gradient_level_ind(1) = 0
+             DO  k = 1, nzt+1
+                IF ( i < 11 ) THEN
+                   IF ( qr_vertical_gradient_level(i) < zu(k)  .AND. &
+                        qr_vertical_gradient_level(i) >= 0.0 )  THEN
+                      gradient = qr_vertical_gradient(i) / 100.0
+                      qr_vertical_gradient_level_ind(i) = k - 1
+                      i = i + 1
+                   ENDIF
+                ENDIF
+                IF ( gradient /= 0.0 )  THEN
+                   IF ( k /= 1 )  THEN
+                      qr_init(k) = qr_init(k-1) + dzu(k) * gradient
+                   ELSE
+                      qr_init(k) = qr_init(k-1) + 0.5 * dzu(k) * gradient
+                   ENDIF
+                ELSE
+                   qr_init(k) = qr_init(k-1)
+                ENDIF
+!
+!--             Avoid negative rain water content
+                IF ( qr_init(k) < 0.0 )  THEN
+                   qr_init(k) = 0.0
+                ENDIF
+             ENDDO
+!
+!--          In case of no given rain water content gradients, choose zero gradient
+!--          conditions
+             IF ( qr_vertical_gradient_level(1) == -1.0 )  THEN
+                qr_vertical_gradient_level(1) = 0.0
+             ENDIF
+
+             i = 1
+             gradient = 0.0
+             nr_vertical_gradient_level_ind(1) = 0
+             DO  k = 1, nzt+1
+                IF ( i < 11 ) THEN
+                   IF ( nr_vertical_gradient_level(i) < zu(k)  .AND. &
+                        nr_vertical_gradient_level(i) >= 0.0 )  THEN
+                      gradient = nr_vertical_gradient(i) / 100.0
+                      nr_vertical_gradient_level_ind(i) = k - 1
+                      i = i + 1
+                   ENDIF
+                ENDIF
+                IF ( gradient /= 0.0 )  THEN
+                   IF ( k /= 1 )  THEN
+                      nr_init(k) = nr_init(k-1) + dzu(k) * gradient
+                   ELSE
+                      nr_init(k) = nr_init(k-1) + 0.5 * dzu(k) * gradient
+                   ENDIF
+                ELSE
+                   nr_init(k) = nr_init(k-1)
+                ENDIF
+!
+!--             Avoid negative rain water content
+                IF ( nr_init(k) < 0.0 )  THEN
+                   nr_init(k) = 0.0
+                ENDIF
+             ENDDO
+!
+!--          In case of no given rain water content gradients, choose zero gradient
+!--          conditions
+             IF ( nr_vertical_gradient_level(1) == -1.0 )  THEN
+                nr_vertical_gradient_level(1) = 0.0
+             ENDIF
+
+          ENDIF
+!
+!--       Store humidity, rain water content and rain drop concentration 
+!--       gradient at the top boundary for possile Neumann boundary condition
+          bc_q_t_val  = ( q_init(nzt+1) - q_init(nzt) ) / dzu(nzt+1)
+            
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             bc_qr_t_val = ( qr_init(nzt+1) - qr_init(nzt) ) / dzu(nzt+1)
+             bc_nr_t_val = ( nr_init(nzt+1) - nr_init(nzt) ) / dzu(nzt+1)
+          ENDIF
 
        ENDIF
@@ -1344 +1458 @@
 !-- Set wind speed in the Galilei-transformed system
     IF ( galilei_transformation )  THEN
-       IF ( use_ug_for_galilei_tr .AND.                &
-            ug_vertical_gradient_level(1) == 0.0 .AND. & 
-            vg_vertical_gradient_level(1) == 0.0 )  THEN
+       IF ( use_ug_for_galilei_tr .AND.                  &
+            ug_vertical_gradient_level(1) == 0.0  .AND.  &
+            ug_vertical_gradient(1) == 0.0  .AND.        & 
+            vg_vertical_gradient_level(1) == 0.0  .AND.  &
+            vg_vertical_gradient(1) == 0.0 )  THEN
           u_gtrans = ug_surface * 0.6
           v_gtrans = vg_surface * 0.6
-       ELSEIF ( use_ug_for_galilei_tr .AND.                &
-                ug_vertical_gradient_level(1) /= 0.0 )  THEN
+       ELSEIF ( use_ug_for_galilei_tr  .AND.                  &
+                ( ug_vertical_gradient_level(1) /= 0.0  .OR.  &
+                ug_vertical_gradient(1) /= 0.0 ) )  THEN
           message_string = 'baroclinicity (ug) not allowed simultaneously' // &
                            ' with galilei transformation'
           CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
-       ELSEIF ( use_ug_for_galilei_tr .AND.                &
-                vg_vertical_gradient_level(1) /= 0.0 )  THEN
+       ELSEIF ( use_ug_for_galilei_tr  .AND.                  &
+                ( vg_vertical_gradient_level(1) /= 0.0  .OR.  &
+                vg_vertical_gradient(1) /= 0.0 ) )  THEN
           message_string = 'baroclinicity (vg) not allowed simultaneously' // &
                            ' with galilei transformation'
@@ -1650 +1768 @@
           CALL message( 'check_parameters', 'PA0074', 1, 2, 0, 6, 0 )
        ENDIF
-       
+
+       IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+          IF ( bc_qr_b == 'dirichlet' )  THEN
+             ibc_qr_b = 0
+          ELSEIF ( bc_qr_b == 'neumann' )  THEN
+             ibc_qr_b = 1
+          ELSE
+             message_string = 'unknown boundary condition: bc_qr_b ="' // &
+                              TRIM( bc_qr_b ) // '"'
+             CALL message( 'check_parameters', 'PA0352', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( bc_qr_t == 'dirichlet' )  THEN
+             ibc_qr_t = 0
+          ELSEIF ( bc_qr_t == 'neumann' )  THEN
+             ibc_qr_t = 1
+          ELSE
+             message_string = 'unknown boundary condition: bc_qr_t ="' // &
+                              TRIM( bc_qr_t ) // '"'
+             CALL message( 'check_parameters', 'PA0353', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( bc_nr_b == 'dirichlet' )  THEN
+             ibc_nr_b = 0
+          ELSEIF ( bc_nr_b == 'neumann' )  THEN
+             ibc_nr_b = 1
+          ELSE
+             message_string = 'unknown boundary condition: bc_nr_b ="' // &
+                              TRIM( bc_nr_b ) // '"'
+             CALL message( 'check_parameters', 'PA0355', 1, 2, 0, 6, 0 )
+          ENDIF
+          IF ( bc_nr_t == 'dirichlet' )  THEN
+             ibc_nr_t = 0
+          ELSEIF ( bc_nr_t == 'neumann' )  THEN
+             ibc_nr_t = 1
+          ELSE
+             message_string = 'unknown boundary condition: bc_nr_t ="' // &
+                              TRIM( bc_nr_t ) // '"'
+             CALL message( 'check_parameters', 'PA0356', 1, 2, 0, 6, 0 )
+          ENDIF
+       ENDIF       
+
     ENDIF
 
@@ -2483 +2640 @@
              hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
 
+          CASE ( 'nr' )
+             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 ( icloud_scheme /= 0 )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_scheme /= seifert_beheng'
+                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 73
+                dopr_unit(i)  = '1/m3'
+                hom(:,2,73,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'qr' )
+             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 ( icloud_scheme /= 0 )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_scheme /= seifert_beheng'
+                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 74
+                dopr_unit(i)  = 'kg/kg'
+                hom(:,2,74,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'qc' )
+             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 ( icloud_scheme /= 0 )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_scheme /= seifert_beheng'
+                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 75
+                dopr_unit(i)  = 'kg/kg'
+                hom(:,2,75,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'prr' )
+             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 ( icloud_scheme /= 0 )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for cloud_scheme /= seifert_beheng'
+                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
+             ELSEIF ( .NOT. precipitation )  THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for precipitation = .FALSE.'
+                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
+
+             ELSE
+                dopr_index(i) = 76
+                dopr_unit(i)  = 'kg/kg m/s'
+                hom(:,2,76,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
           CASE DEFAULT
 
@@ -2572 +2803 @@
              unit = 'K'
 
+          CASE ( 'nr' )
+             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 ( icloud_scheme /= 0 )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                         'res cloud_scheme = seifert_beheng'
+                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = '1/m3'
+
           CASE ( 'pc', 'pr' )
              IF ( .NOT. particle_advection )  THEN
@@ -2581 +2824 @@
              IF ( TRIM( var ) == 'pr' )  unit = 'm'
 
+          CASE ( 'prr' )
+             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 ( icloud_scheme /= 0 )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                         'res cloud_scheme = seifert_beheng'
+                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
+             ELSEIF ( .NOT. precipitation )  THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                                 'res precipitation = .TRUE.'
+                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'kg/kg m/s'
+
           CASE ( 'q', 'vpt' )
              IF ( .NOT. humidity )  THEN
@@ -2590 +2849 @@
              IF ( TRIM( var ) == 'vpt' )  unit = 'K'
 
+          CASE ( 'qc' )
+             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 ( icloud_scheme /= 0 ) THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                         'res cloud_scheme = seifert_beheng'
+                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'kg/kg'
+
           CASE ( 'ql' )
              IF ( .NOT. ( cloud_physics  .OR.  cloud_droplets ) )  THEN
@@ -2607 +2878 @@
              IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
              IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
+
+          CASE ( 'qr' )
+             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 ( icloud_scheme /= 0 ) THEN
+                message_string = 'output of "' // TRIM( var ) // '" requi' // &
+                         'res cloud_scheme = seifert_beheng'
+                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
+             ENDIF
+             unit = 'kg/kg'
 
           CASE ( 'qv' )

palm/trunk/SOURCE/data_output_2d.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! +qr, nr, qc and cross sections
 !
 ! Former revisions:
@@ -290 +290 @@
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => nr
+                ELSE
+                   to_be_resorted => nr_av
+                ENDIF
+                IF ( mode == 'xy' )  level_z = zu
 
              CASE ( 'p_xy', 'p_xz', 'p_yz' )
@@ -376 +384 @@
 
              CASE ( 'prr*_xy' )        ! 2d-array
-                IF ( av == 0 )  THEN
-                   CALL exchange_horiz_2d( precipitation_rate )
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) =  precipitation_rate(j,i)
-                      ENDDO
-                   ENDDO
-                ELSE
-                   CALL exchange_horiz_2d( precipitation_rate_av )
-                   DO  i = nxlg, nxrg
-                      DO  j = nysg, nyng
-                         local_pf(i,j,nzb+1) =  precipitation_rate_av(j,i)
-                      ENDDO
-                   ENDDO
+                IF ( icloud_scheme == 1 )  THEN
+                   IF ( av == 0 )  THEN
+                      CALL exchange_horiz_2d( precipitation_rate )
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            local_pf(i,j,nzb+1) =  precipitation_rate(j,i)
+                         ENDDO
+                      ENDDO
+                   ELSE
+                      CALL exchange_horiz_2d( precipitation_rate_av )
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            local_pf(i,j,nzb+1) =  precipitation_rate_av(j,i)
+                         ENDDO
+                      ENDDO
+                   ENDIF
+                ELSE
+                   IF ( av == 0 )  THEN
+                      CALL exchange_horiz_2d( prr(nzb+1,:,:) )
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            local_pf(i,j,nzb+1) = prr(nzb+1,j,i) * hyrho(nzb+1)
+                         ENDDO
+                      ENDDO
+                   ELSE
+                      CALL exchange_horiz_2d( prr_av(nzb+1,:,:) )
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            local_pf(i,j,nzb+1) = prr_av(nzb+1,j,i) * hyrho(nzb+1)
+                         ENDDO
+                      ENDDO
+                   ENDIF
                 ENDIF
                 resorted = .TRUE.
                 two_d = .TRUE.
                 level_z(nzb+1) = zu(nzb+1)
+
+             CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
+                IF ( av == 0 )  THEN
+                   CALL exchange_horiz( prr, nbgp )
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            local_pf(i,j,k) = prr(k,j,i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ELSE
+                   CALL exchange_horiz( prr_av, nbgp )
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nzt+1
+                            local_pf(i,j,k) = prr_av(k,j,i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                ENDIF
+                resorted = .TRUE.
+                IF ( mode == 'xy' )  level_z = zu
 
              CASE ( 'pt_xy', 'pt_xz', 'pt_yz' )
@@ -424 +473 @@
                 IF ( mode == 'xy' )  level_z = zu
 
+             CASE ( 'qc_xy', 'qc_xz', 'qc_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => ql
+                ELSE
+                   to_be_resorted => ql_av
+                ENDIF
+                IF ( mode == 'xy' )  level_z = zu
+
              CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
                 IF ( av == 0 )  THEN
-                   to_be_resorted => ql
-                ELSE
-                   to_be_resorted => ql_av
+                   IF ( icloud_scheme == 0 )  THEN
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            DO  k = nzb, nz_do3d
+                               local_pf(i,j,k) = ql(k,j,i) + qr(k,j,i)
+                            ENDDO
+                         ENDDO
+                      ENDDO
+                      resorted = .TRUE.
+                   ELSE
+                      to_be_resorted => ql
+                   ENDIF
+                ELSE
+                   IF ( icloud_scheme == 0 )  THEN
+                      DO  i = nxlg, nxrg
+                         DO  j = nysg, nyng
+                            DO  k = nzb, nz_do3d
+                               local_pf(i,j,k) = ql_av(k,j,i) + qr_av(k,j,i)
+                            ENDDO
+                         ENDDO
+                      ENDDO
+                      resorted = .TRUE.
+                   ELSE
+                      to_be_resorted => ql_av
+                   ENDIF
                 ENDIF
                 IF ( mode == 'xy' )  level_z = zu
@@ -478 +557 @@
                    CALL exchange_horiz( ql_vp_av, nbgp )
                    to_be_resorted => ql_vp
+                ENDIF
+                IF ( mode == 'xy' )  level_z = zu
+
+             CASE ( 'qr_xy', 'qr_xz', 'qr_yz' )
+                IF ( av == 0 )  THEN
+                   to_be_resorted => qr
+                ELSE
+                   to_be_resorted => qr_av
                 ENDIF
                 IF ( mode == 'xy' )  level_z = zu

palm/trunk/SOURCE/data_output_3d.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! +nr, qr, prr, qc and averaged quantities
 !
 ! Former revisions:
@@ -207 +207 @@
              ENDIF
 
+          CASE ( 'nr' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => nr
+             ELSE
+                to_be_resorted => nr_av
+             ENDIF
+
           CASE ( 'p' )
              IF ( av == 0 )  THEN
@@ -270 +277 @@
              ENDIF
 
+          CASE ( 'prr' )
+             IF ( av == 0 )  THEN
+                CALL exchange_horiz( prr, nbgp )
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nzt+1
+                         local_pf(i,j,k) = prr(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ELSE
+                CALL exchange_horiz( prr_av, nbgp )
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      DO  k = nzb, nzt+1
+                         local_pf(i,j,k) = prr_av(k,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
+             resorted = .TRUE.
+
           CASE ( 'pt' )
              IF ( av == 0 )  THEN
@@ -297 +326 @@
              ENDIF
 
+          CASE ( 'qc' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => ql
+             ELSE
+                to_be_resorted => ql_av
+             ENDIF
+
           CASE ( 'ql' )
              IF ( av == 0 )  THEN
-                to_be_resorted => ql
-             ELSE
-                to_be_resorted => ql_av
+                IF ( icloud_scheme == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nz_do3d
+                            local_pf(i,j,k) = ql(k,j,i) + qr(k,j,i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                   resorted = .TRUE.
+                ELSE
+                   to_be_resorted => ql
+                ENDIF
+             ELSE
+                IF ( icloud_scheme == 0 )  THEN
+                   DO  i = nxlg, nxrg
+                      DO  j = nysg, nyng
+                         DO  k = nzb, nz_do3d
+                            local_pf(i,j,k) = ql_av(k,j,i) + qr_av(k,j,i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+                   resorted = .TRUE.
+                ELSE
+                   to_be_resorted => ql_av
+                ENDIF
              ENDIF
 
@@ -344 +402 @@
                 CALL exchange_horiz( ql_vp_av, nbgp )
                 to_be_resorted => ql_vp_av
+             ENDIF
+
+          CASE ( 'qr' )
+             IF ( av == 0 )  THEN
+                to_be_resorted => qr
+             ELSE
+                to_be_resorted => qr_av
              ENDIF
 

palm/trunk/SOURCE/flow_statistics.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary additions for two-moment cloud physics scheme:
+! +nr, qr, qc, prr
 !
 ! Former revisions:
@@ -748 +749 @@
 
 !
-!--             Buoyancy flux, water flux, humidity flux and liquid water
-!--             content
+!--             Buoyancy flux, water flux, humidity flux, liquid water
+!--             content, rain drop concentration and rain water content
                 IF ( humidity )  THEN
                    IF ( cloud_physics .OR. cloud_droplets )  THEN
-                      pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
+                      pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +        &
                                     vpt(k+1,j,i) - hom(k+1,1,44,sr) )
                       sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
                                                           rmask(j,i,sr)
-                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
-                                                          rmask(j,i,sr)
+                      IF ( .NOT. cloud_droplets )  THEN
+                         pts = 0.5 *                                          &
+                              ( ( q(k,j,i) - ql(k,j,i) ) -                    &
+                              hom(k,1,42,sr) +                                &
+                              ( q(k+1,j,i) - ql(k+1,j,i) ) -                  &
+                              hom(k+1,1,42,sr) )
+                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
+                                                             rmask(j,i,sr)
+                         IF ( icloud_scheme == 0  )  THEN
+                            sums_l(k,54,tn) = sums_l(k,54,tn) + ( ql(k,j,i) + &
+                                                                qr(k,j,i) ) * &
+                                                                rmask(j,i,sr)
+                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *   &
+                                                                rmask(j,i,sr)
+                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *   &
+                                                                rmask(j,i,sr)
+                            sums_l(k,75,tn) = sums_l(k,75,tn) + ql(k,j,i) *   &
+                                                                rmask(j,i,sr)
+                            sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *  &
+                                                                rmask(j,i,sr)
+                         ELSE
+                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *   &
+                                                                rmask(j,i,sr)
+                         ENDIF
+                      ELSE
+                         sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *      &
+                                                             rmask(j,i,sr)
+                      ENDIF
                    ELSE
                       IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
@@ -765 +792 @@
                                                              rmask(j,i,sr)
                       ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
-                         sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * &
+                         sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) *  &
                                              sums_l(k,17,tn) +      &
                                           0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
@@ -1105 +1132 @@
        hom(:,1,71,sr) = sums(:,71)     ! prho
        hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
+       hom(:,1,73,sr) = sums(:,73)     ! nr
+       hom(:,1,74,sr) = sums(:,74)     ! qr
+       hom(:,1,75,sr) = sums(:,75)     ! qc
+       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
 
        hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)

palm/trunk/SOURCE/init_3d_model.f90

@@ -23 +23 @@
 ! Current revisions:
 ! ------------------
-!
+! allocation and initialisation of necessary data arrays for the two-moment 
+! cloud physics scheme the two new prognostic equations (nr, qr):
+! +dr, lambda_r, mu_r, sed_*, xr, *s, *sws, *swst, *, *_p, t*_m, *_1, *_2, *_3, 
+! +tend_*, prr
 !
 ! Former revisions:
@@ -353 +356 @@
 
           IF ( cloud_physics ) THEN
+
 !
 !--          Liquid water content
@@ -364 +368 @@
              ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg), &
                        precipitation_rate(nysg:nyng,nxlg:nxrg) )
+
+             IF ( icloud_scheme == 0 )  THEN
+!
+!--             1D-rain sedimentation fluxes and rain drop size distribution
+!--             properties 
+                ALLOCATE ( dr(nzb:nzt+1), lambda_r(nzb:nzt+1),  &
+                           mu_r(nzb:nzt+1), sed_q(nzb:nzt+1),   &
+                           sed_qr(nzb:nzt+1), sed_nr(nzb:nzt+1),& 
+                           xr(nzb:nzt+1) ) 
+!
+!--             2D-rain water content and rain drop concentration arrays
+                ALLOCATE ( qrs(nysg:nyng,nxlg:nxrg),   &
+                   qrsws(nysg:nyng,nxlg:nxrg),         &
+                   qrswst(nysg:nyng,nxlg:nxrg),        &
+                   nrs(nysg:nyng,nxlg:nxrg),           &
+                   nrsws(nysg:nyng,nxlg:nxrg),         &
+                   nrswst(nysg:nyng,nxlg:nxrg) )
+!
+!--             3D-rain water content, rain drop concentration arrays
+#if defined( __nopointer )
+                ALLOCATE( nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),         &
+                          nr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),         &
+                          qr_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          tnr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg),      &
+                          tqr_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+#else
+                ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
+                          qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+#endif
+!
+!--             3D-tendency arrays
+                ALLOCATE( tend_nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
+                          tend_pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg),    &
+                          tend_q(nzb:nzt+1,nysg:nyng,nxlg:nxrg),     &
+                          tend_qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+!
+!--             3d-precipitation rate
+                IF ( precipitation )  THEN
+                   ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+
+             ENDIF
           ENDIF
 
@@ -508 +559 @@
     IF ( humidity  .OR.  passive_scalar )  THEN
        q => q_1;  q_p => q_2;  tq_m => q_3
-       IF ( humidity )        vpt  => vpt_1
-       IF ( cloud_physics )   ql   => ql_1
+       IF ( humidity )  THEN
+          vpt  => vpt_1    
+          IF ( cloud_physics )  THEN
+             ql => ql_1
+             IF ( icloud_scheme == 0 )  THEN
+                qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
+                nr => nr_1;  nr_p  => nr_2;  tnr_m  => nr_3
+             ENDIF
+          ENDIF
+       ENDIF
        IF ( cloud_droplets )  THEN
           ql   => ql_1
@@ -564 +623 @@
                 ENDDO
              ENDDO
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qr(:,j,i) = qr_init
+                      nr(:,j,i) = nr_init
+                   ENDDO
+                ENDDO
+             ENDIF
           ENDIF
 
@@ -606 +673 @@
 !--       This could actually be computed more accurately in the 1D model.
 !--       Update when opportunity arises!
-          IF ( humidity  .OR.  passive_scalar )  qs = 0.0
+          IF ( humidity  .OR.  passive_scalar )  THEN
+             qs = 0.0
+             IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+                qrs = 0.0
+                nrs = 0.0
+             ENDIF
+          ENDIF
 
 !
@@ -671 +744 @@
                 ENDDO
              ENDDO
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      qr(:,j,i) = qr_init
+                      nr(:,j,i) = nr_init
+                   ENDDO
+                ENDDO
+             ENDIF
           ENDIF
 
@@ -811 +892 @@
 !--       Determine the near-surface water flux
           IF ( humidity  .OR.  passive_scalar )  THEN
+             IF ( cloud_physics .AND. icloud_scheme )  THEN
+                IF ( constant_waterflux_qr )  THEN
+                   qrsws = surface_waterflux_qr
+                ENDIF
+                IF (constant_waterflux_nr )  THEN
+                   nrsws = surface_waterflux_nr
+                ENDIF
+             ENDIF
              IF ( constant_waterflux )  THEN
                 qsws   = surface_waterflux
@@ -842 +931 @@
              tswst = top_heatflux
 
-             IF ( humidity  .OR.  passive_scalar )  qswst = 0.0
+             IF ( humidity  .OR.  passive_scalar )  THEN
+                qswst = 0.0
+                IF ( cloud_physics  .AND.  icloud_scheme == 0 ) THEN
+                   nrswst = 0.0
+                   qrswst = 0.0
+                ENDIF
+             ENDIF
 
              IF ( ocean )  THEN
@@ -877 +972 @@
           IF ( humidity  .OR.  passive_scalar )  THEN
              IF ( .NOT. constant_waterflux )  qsws   = 0.0
+             IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+                IF ( .NOT. constant_waterflux_qr )  THEN
+                   qrsws = 0.0
+                ENDIF
+                IF ( .NOT. constant_waterflux_nr )  THEN
+                   nrsws = 0.0
+                ENDIF
+             ENDIF
           ENDIF
 
@@ -925 +1028 @@
             q_surface_initial_change /= 0.0 )  THEN
           q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             IF ( qr_surface_initial_change /= 0.0 )  THEN
+                qr(nzb,:,:) = qr(nzb,:,:) + qr_surface_initial_change
+             ELSEIF ( nr_surface_initial_change /= 0.0 ) THEN
+                nr(nzb,:,:) = nr(nzb,:,:) + nr_surface_initial_change
+             ENDIF
+          ENDIF
        ENDIF
 
@@ -939 +1049 @@
           tq_m = 0.0
           q_p = q
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             tqr_m = 0.0
+             qr_p = qr
+             tnr_m = 0.0
+             nr_p = nr
+          ENDIF
        ENDIF
 
@@ -1109 +1225 @@
 !--    including ghost points)
        e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
-       IF ( humidity  .OR.  passive_scalar )  q_p = q
+       IF ( humidity  .OR.  passive_scalar )  THEN
+          q_p = q
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             qr_p = qr
+             nr_p = nr
+          ENDIF
+       ENDIF
        IF ( ocean )  sa_p = sa
 
@@ -1117 +1239 @@
 !--    there before they are set. 
        te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
-       IF ( humidity  .OR.  passive_scalar )  tq_m = 0.0
+       IF ( humidity  .OR.  passive_scalar )  THEN
+          tq_m = 0.0
+          IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+             tqr_m = 0.0
+             tnr_m = 0.0
+          ENDIF
+       ENDIF
        IF ( ocean )  tsa_m = 0.0
 

palm/trunk/SOURCE/init_cloud_physics.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! calculation of the maximum timestep according to the terminal velocity of rain 
+! drops in the two moment cloud scheme
+!
+! calculation of frequently used constants (pirho_l, dpirho_l, schmidt_p_1d3, 
+! hyrho) 
 !
 ! Former revisions:
@@ -66 +70 @@
     REAL    ::  t_surface
 
-    ALLOCATE( hyp(nzb:nzt+1), pt_d_t(nzb:nzt+1), t_d_pt(nzb:nzt+1) )
+    ALLOCATE( hyp(nzb:nzt+1), pt_d_t(nzb:nzt+1), t_d_pt(nzb:nzt+1),  &
+              hyrho(nzb:nzt+1) )
 
 !
@@ -73 +78 @@
     l_d_r  = l_v / r_d
     l_d_rv = l_v / r_v
+
+    schmidt_p_1d3 = schmidt**( 1.0 / 3.0 )
+
+    pirho_l  = pi * rho_l / 6.0
+    dpirho_l = 1.0 / pirho_l 
+!
+!-- Calculate timestep according to precipitation
+    IF ( icloud_scheme == 0  .AND.  precipitation )  THEN
+       dt_precipitation = MINVAL( dzu(nzb+2:nzt) ) / w_precipitation
+    ENDIF
 !
 !-- Calculate factor used in equation for droplet growth by condensation
@@ -94 +109 @@
                    ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0/0.286)
        pt_d_t(k) = ( 100000.0 / hyp(k) )**0.286
-       t_d_pt(k) = 1.0 / pt_d_t(k)       
+       t_d_pt(k) = 1.0 / pt_d_t(k)
+       hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
     ENDDO
 

palm/trunk/SOURCE/modules.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme:
+! +*_init, flux_l_*, diss_l_*, flux_s_*, diss_s_*, *sws, *swst, tend_*, *, *_p
+! +t*_m, *_1, *_2, *_3, *_av, bc_*_b, bc_*_t, ibc_*_b, ibc_*_t, bc_*_t_val, 
+! +*_vertical_gradient, *_surface_initial_change, *_vertical_gradient_level, 
+! +*_vertical_gradient_level_ind, *_surface, constant_waterflux_*,  
+! +cloud_scheme, icloud_scheme, surface_waterflux_*, sums_ws*s_ws_l, wall_*flux
+!
+! constants for the two-moment scheme:
+! +a_vent, a_term, b_vent, b_term, c_evap, c_term, cof, eps_sb, k_cc, k_cr, k_rr, 
+! +k_br, kappa_rr, kin_vis_air, mu_constant_value, nc, pirho_l, dpirho_l, rho_1,
+! +schmidt, schmidt_p_1d3, stp, x0, xmin, xmax, dt_precipitation, w_precipitation
+!
+! steering parameters for the two_moment scheme:
+! +mu_constant, ventilation_effect
 !
 ! Former revisions:
@@ -369 +383 @@
           c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l, ddzu, ddzu_pres,     &
           dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, lad, l_grid,      &
-          ptdf_x, ptdf_y, pt_init, q_init, rdf, rdf_sc, sa_init, ug, u_init,   &
-          u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs, zu, zw 
+          nr_init, ptdf_x, ptdf_y, pt_init, q_init, qr_init, rdf, rdf_sc,      &
+          sa_init, ug, u_init, u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, &
+          w_subs, zu, zw
 
     REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
-          c_u, c_v, c_w, diss_s_e, diss_s_pt, diss_s_q, diss_s_sa, diss_s_u,   &
-          diss_s_v, diss_s_w, dzu_mg, dzw_mg, flux_s_e, flux_s_pt, flux_s_q,   &
-          flux_s_sa, flux_s_u, flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg,        &
-          mean_inflow_profiles, pt_slope_ref, qs, qsws, qswst, qswst_remote,   &
-          rif, saswsb, saswst, shf, total_2d_a, total_2d_o, ts, tswst, us,     &
-          usws, uswst, vsws, vswst, z0, z0h
+          c_u, c_v, c_w, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q, diss_s_qr,  &
+          diss_s_sa, diss_s_u, diss_s_v, diss_s_w, dzu_mg, dzw_mg, flux_s_e,   &
+          flux_s_nr, flux_s_pt, flux_s_q, flux_s_qr, flux_s_sa, flux_s_u,      &
+          flux_s_v, flux_s_w, f1_mg, f2_mg, f3_mg, mean_inflow_profiles, nrs,  &
+          nrsws, nrswst, pt_slope_ref, qs, qsws, qswst, qswst_remote, qrs,     &
+          qrsws, qrswst, rif, saswsb, saswst, shf, total_2d_a, total_2d_o, ts, &
+          tswst, us, usws, uswst, vsws, vswst, z0, z0h
+          
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
           canopy_heat_flux, cdc, d, de_dx, de_dy, de_dz, diss, diss_l_e,       &
-          diss_l_pt, diss_l_q, diss_l_sa, diss_l_u, diss_l_v, diss_l_w,        &
-          flux_l_e, flux_l_pt, flux_l_q, flux_l_sa, flux_l_u, flux_l_v,        &
-          flux_l_w, kh, km, lad_s, lad_u, lad_v, lad_w, lai, l_wall, p_loc,    &
-          sec, sls, tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r,     &
-          v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
+          diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr, 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_sa, flux_l_u, flux_l_v, flux_l_w, kh, km, lad_s,   &
+          lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, tend_pt,    &
+          tend_nr, tend_q, tend_qr, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n,  &
+          v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
+            
 
 #if defined( __nopointer )
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
-          e, e_p, p, prho, pt, pt_p, q, q_p, ql, ql_c, ql_v, ql_vp, rho, sa,   &
-          sa_p, te_m, tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p,    &
-          vpt, w, w_p
+          e, e_p, nr, nr_p, p, prho, pt, pt_p, q, q_p, ql, ql_c, ql_v, ql_vp,  &
+          qr, qr_p, rho, sa, sa_p, te_m, tnr_m, tpt_m, tq_m, tqr_m, tsa_m,     &
+          tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p
 #else
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
-          e_1, e_2, e_3, p, prho_1, pt_1, pt_2, pt_3, q_1, q_2, q_3, ql_v,     &
-          ql_vp, ql_1, ql_2, rho_1, sa_1, sa_2, sa_3, u_1, u_2, u_3,           &
-          v_1, v_2, v_3, vpt_1, w_1, w_2, w_3
+          e_1, e_2, e_3, p, prho_1, nr_1, nr_2, nr_3, pt_1, pt_2, pt_3, q_1,   &
+          q_2, q_3, ql_v, ql_vp, ql_1, ql_2, qr_1, qr_2, qr_3, rho_1, sa_1,    &
+          sa_2, sa_3, u_1, u_2, u_3, v_1, v_2, v_3, vpt_1, w_1, w_2, w_3
 
     REAL, DIMENSION(:,:,:), POINTER ::                                         &
-          e, e_p, prho, pt, pt_p, q, q_p, ql, ql_c, rho, sa, sa_p, te_m,       &
-          tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, u, u_p, v, v_p, vpt, w, w_p
+          e, e_p, nr, nr_p, prho, pt, pt_p, q, q_p, ql, ql_c, qr, qr_p, rho,   &
+          sa, sa_p, te_m, tnr_m, tpt_m, tq_m, tqr_m, tsa_m, tu_m, tv_m, tw_m,  &
+          u, u_p, v, v_p, vpt, w, w_p
 #endif
 
@@ -429 +449 @@
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
-          e_av, lpt_av, p_av, pc_av, pr_av, pt_av, q_av, ql_av, ql_c_av,     &
-          ql_v_av, ql_vp_av, qv_av, rho_av, s_av, sa_av, u_av, v_av, vpt_av, &
-          w_av
+          e_av, lpt_av, nr_av, p_av, pc_av, pr_av, prr_av, pt_av, q_av, ql_av, &
+          ql_c_av, ql_v_av, ql_vp_av, qr_av, qv_av, rho_av, s_av, sa_av, u_av, &
+          v_av, vpt_av, w_av
+  
  END MODULE averaging
 
@@ -444 +465 @@
 !------------------------------------------------------------------------------!
 
-    LOGICAL ::  curvature_solution_effects = .FALSE.
-
-    REAL  ::  bfactor, cp = 1005.0, diff_coeff_l = 0.23E-4,                    &
-              effective_coll_efficiency,                                       &
-              eps_ros = 1.0E-4,             &  ! accuracy of Rosenbrock method
-              l_d_cp, l_d_r, l_d_rv, l_v = 2.5E+06,                            &
-              mass_of_solute = 1.0E-17,              & ! soluted NaCl in kg
-              molecular_weight_of_solute = 0.05844,  & ! mol. mass NaCl (kg/mol)
-              molecular_weight_of_water = 0.01801528,& ! mol. mass H2O (kg/mol)
-              prec_time_const = 0.001, ql_crit = 0.0005, rho_l = 1.0E3,        &
-              r_d = 287.0, r_v = 461.51, thermal_conductivity_l = 2.43E-2,     &
-              vanthoff = 2.0                        ! van't Hoff factor (NaCl)
-
-    REAL, DIMENSION(:), ALLOCATABLE   ::  pt_d_t, t_d_pt
-
-    REAL, DIMENSION(:,:), ALLOCATABLE ::  precipitation_amount, &
-                                          precipitation_rate
+    LOGICAL ::  curvature_solution_effects = .FALSE., &
+                ventilation_effect = .FALSE., &
+                mu_constant = .FALSE.
+
+    REAL ::  a_vent = 0.78,      & ! coef. for ventilation effect 
+             a_term = 9.65,      & ! coef. for terminal velocity (m s-1)
+             b_vent = 0.308,     & ! coef. for ventilation effect
+             b_term = 9.8,       & ! coef. for terminal velocity (m s-1)
+             bfactor,            &
+             c_evap = 0.7,       & ! constant in evaporation
+             c_term = 600.0,     & ! coef. for terminal velocity (m-1)
+             cof(6) = (/ 76.18009172947146,      & ! coefficients in the 
+                         -86.50532032941677,     & ! numerical
+                         24.01409824083091,      & ! calculation of the
+                         -1.231739572450155,     & ! gamma function
+                         0.1208650973866179E-2,  &
+                         -0.5395239384953E-5 /), &
+             cp = 1005.0,        & ! heat capacity of dry air (J kg-1 K-1) 
+             diff_coeff_l = 0.23E-4, & ! diffusivity of water vapor (m2 s-1)
+             effective_coll_efficiency, &
+             eps_ros = 1.0E-4,   & ! accuracy of Rosenbrock method
+             eps_sb = 1.0E-20,   & ! threshold in two-moments scheme
+             k_cc = 4.44E09,     & ! const. rain-rain kernel (m3 kg-2 s-1) 
+             k_cr = 5.25,        & ! const. cloud-rain kernel (m3 kg-1 s-1) 
+             k_rr = 7.12,        & ! const. rain-rain kernel (m3 kg-1 s-1) 
+             k_br = 1000.,       & ! const. in breakup parametrization (m-1)
+             kappa_rr = 60.7,    & ! const. in collision kernel (kg-1/3)
+             kin_vis_air = 1.4086E-5, & ! kin. viscosity of air (m2 s-1)
+             l_v = 2.5E+06,      & ! latent heat of vaporization (J kg-1)
+             l_d_cp, l_d_r, l_d_rv, & ! l_v / cp, l_v / r_d, l_v / r_v
+             mass_of_solute = 1.0E-17, & ! soluted NaCl (kg)
+             molecular_weight_of_solute = 0.05844, & ! mol. m. NaCl (kg mol-1)
+             molecular_weight_of_water = 0.01801528, & ! mol. m. H2O (kg mol-1)
+             mu_constant_value = 0.0, & ! shape param. of gamma distribution
+             nc = 70.0E6,        & ! cloud droplet concentration 
+             prec_time_const = 0.001, & !coef. in Kessler scheme
+             pirho_l, dpirho_l,  & ! pi * rho_l / 6.0; 6.0 / ( pi * rho_l )
+             rho_l = 1.0E3,      & ! density of water (kg m-3)
+             ql_crit = 0.0005,   & ! coef. in Kessler scheme
+             r_d = 287.0,        & ! sp. gas const. dry air (J kg-1 K-1)
+             r_v = 461.51,       & ! sp. gas const. water vapor (J kg-1 K-1)
+             schmidt = 0.71,     & ! Schmidt number
+             schmidt_p_1d3,      & ! schmidt**( 1.0 / 3.0 )
+             stp = 2.5066282746310005, & ! parameter in gamma function
+             thermal_conductivity_l = 2.43E-2, & ! therm. cond. air (J m-1 s-1 K-1)
+             vanthoff = 2.0,     & ! van't Hoff factor for NaCl
+             x0 = 2.6E-10,       & ! separating drop mass (kg)
+             xrmin = 2.6E-10,    & ! minimum rain drop size (kg)
+             xrmax = 5.0E-6,     & ! maximum rain drop site (kg)
+             dt_precipitation = 100.0, & ! timestep precipitation (s)
+             w_precipitation = 9.65      ! maximum terminal velocity (m s-1)
+
+    REAL, DIMENSION(:), ALLOCATABLE     ::  dr, hyrho, lambda_r, mu_r, pt_d_t, &
+                                            sed_nr, sed_q, sed_qr, t_d_pt, xr  
+
+    REAL, DIMENSION(:,:), ALLOCATABLE   ::  precipitation_amount, &
+                                            precipitation_rate
+!
+!-- 3D array of precipitation rate
+    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  prr
 
     SAVE
@@ -529 +593 @@
                              scalar_advec = 'ws-scheme'
     CHARACTER (LEN=20)   ::  bc_e_b = 'neumann', bc_lr = 'cyclic', &
+                             bc_nr_b = 'dirichlet', bc_nr_t = 'neumann', &
                              bc_ns = 'cyclic', bc_p_b = 'neumann', &
                              bc_p_t = 'dirichlet', bc_pt_b = 'dirichlet', &
                              bc_pt_t = 'initial_gradient', &
                              bc_q_b = 'dirichlet', bc_q_t = 'neumann', &
+                             bc_qr_b = 'dirichlet', bc_qr_t = 'neumann',&
                              bc_s_b = 'dirichlet', bc_s_t = 'neumann', &
                              bc_sa_t = 'neumann', &
                              bc_uv_b = 'dirichlet', bc_uv_t = 'dirichlet', &
                              canopy_mode = 'block', &
+                             cloud_scheme = 'seifert_beheng', &
                              coupling_mode = 'uncoupled', &
                              coupling_mode_remote = 'uncoupled', &
@@ -544 +611 @@
                              random_generator = 'numerical-recipes', &
                              return_addres, return_username, &
-                             timestep_scheme = 'runge-kutta-3'
+                             timestep_scheme = 'runge-kutta-3'                             
     CHARACTER (LEN=40)   ::  avs_data_file, output_format_netcdf, &
                              topography = 'flat'
@@ -577 +644 @@
                 dp_level_ind_b = 0, dvrp_filecount = 0, &
                 dz_stretch_level_index, gamma_mg, gathered_size, &
-                grid_level, ibc_e_b, ibc_p_b, ibc_p_t, ibc_pt_b, ibc_pt_t, &
-                ibc_q_b, ibc_q_t, ibc_sa_t, ibc_uv_b, ibc_uv_t, &
+                grid_level, ibc_e_b, ibc_nr_b, ibc_nr_t, ibc_p_b, ibc_p_t, &
+                ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_qr_b, ibc_qr_t, &
+                ibc_sa_t, ibc_uv_b, ibc_uv_t, icloud_scheme, &
                 inflow_disturbance_begin = -1, inflow_disturbance_end = -1, &
                 intermediate_timestep_count, intermediate_timestep_count_max, &
@@ -600 +668 @@
                 mask_size(max_masks,3) = -1, mask_size_l(max_masks,3) = -1, &
                 mask_start_l(max_masks,3) = -1, &
+                nr_vertical_gradient_level_ind(10) = -9999, &
                 pt_vertical_gradient_level_ind(10) = -9999, &
                 q_vertical_gradient_level_ind(10) = -9999, &
+                qr_vertical_gradient_level_ind(10) = -9999, &
                 sa_vertical_gradient_level_ind(10) = -9999, &
                 section(100,3), section_xy(100) = -9999, &
@@ -631 +701 @@
                 constant_top_momentumflux = .FALSE., &
                 constant_top_salinityflux = .TRUE., &
-                constant_waterflux = .TRUE., create_disturbances = .TRUE., &
+                constant_waterflux = .TRUE., constant_waterflux_nr = .TRUE., &
+                constant_waterflux_qr = .TRUE., create_disturbances = .TRUE., &
                 data_output_2d_on_each_pe = .TRUE., & 
                 dissipation_control = .FALSE., disturbance_created = .FALSE., &
                 do2d_at_begin = .FALSE., do3d_at_begin = .FALSE., &
                 do3d_compress = .FALSE., do_sum = .FALSE., &
-                dp_external = .FALSE., dp_smooth = .FALSE., dt_fixed = .FALSE., &
+                dp_external = .FALSE., dp_smooth = .FALSE., &
+                drizzle = .TRUE., dt_fixed = .FALSE., &
                 dt_3d_reached, dt_3d_reached_l, exchange_mg = .FALSE., &
                 first_call_lpm = .TRUE., &
@@ -668 +740 @@
              averaging_interval = 0.0, averaging_interval_pr = 9999999.9, &
              averaging_interval_sp = 9999999.9, bc_pt_t_val, bc_q_t_val, &
-             bottom_salinityflux = 0.0, &
+             bc_qr_t_val, bc_nr_t_val, bottom_salinityflux = 0.0, &
              building_height = 50.0, building_length_x = 50.0, &
              building_length_y = 50.0, building_wall_left = 9999999.9, &
@@ -699 +771 @@
              maximum_cpu_time_allowed = 0.0,  &
              molecular_viscosity = 1.461E-5, &
+             nr_surface = 0.0, nr_surface_initial_change = 0.0, &
              old_dt = 1.0E-10, omega = 7.29212E-5, omega_sor = 1.8, &
              particle_maximum_age = 9999999.9, &
@@ -707 +780 @@
              pt_surface = 300.0, pt_surface_initial_change = 0.0, &
              q_surface = 0.0, q_surface_initial_change = 0.0, &
+             qr_surface = 0.0, qr_surface_initial_change = 0.0, &
              rayleigh_damping_factor = -1.0, rayleigh_damping_height = -1.0, &
              recycling_width = 9999999.9, residual_limit = 1.0E-4, &
@@ -719 +793 @@
              surface_heatflux = 9999999.9, surface_pressure = 1013.25, &
              surface_scalarflux = 9999999.9, surface_waterflux = 9999999.9, &
+             surface_waterflux_nr = 0.0, surface_waterflux_qr = 0.0, &
              s_surface = 0.0, s_surface_initial_change = 0.0, &
              termination_time_needed = -1.0, time_coupling = 0.0, &
@@ -740 +815 @@
              lad_vertical_gradient_level(10) = -9999999.9, &
              mask_scale(3), &
+             nr_vertical_gradient(10) = 0.0, &
+             nr_vertical_gradient_level(10) = -1.0, &
              pt_vertical_gradient(10) = 0.0, &
              pt_vertical_gradient_level(10) = -9999999.9, &
              q_vertical_gradient(10) = 0.0, &
              q_vertical_gradient_level(10) = -1.0, &
+             qr_vertical_gradient(10) = 0.0, &
+             qr_vertical_gradient_level(10) = -1.0, &
              s_vertical_gradient(10) = 0.0, &
              s_vertical_gradient_level(10) = -1.0, &
@@ -759 +838 @@
              volume_flow(1:2) = 0.0, volume_flow_area(1:2) = 0.0, &
              volume_flow_initial(1:2) = 0.0, wall_heatflux(0:4) = 0.0, &
-             wall_humidityflux(0:4) = 0.0, wall_qflux(0:4) = 0.0, &
+             wall_humidityflux(0:4) = 0.0, wall_nrflux(0:4) = 0.0, & 
+             wall_qflux(0:4) = 0.0, wall_qrflux(0:4) = 0.0, &
              wall_salinityflux(0:4) = 0.0, wall_scalarflux(0:4) = 0.0, &
              subs_vertical_gradient(10) = 0.0, &
@@ -1532 +1612 @@
                                               sums_us2_ws_l, sums_vs2_ws_l,   &
                                               sums_ws2_ws_l,                  &
+                                              sums_wsnrs_ws_l,                &
                                               sums_wspts_ws_l,                &
-                                              sums_wssas_ws_l,sums_wsqs_ws_l
+                                              sums_wssas_ws_l,                &
+                                              sums_wsqs_ws_l,                 &
+                                              sums_wsqrs_ws_l
+                                              
     REAL, DIMENSION(:,:,:), ALLOCATABLE   ::  hom_sum, rmask, spectrum_x, &
                                               spectrum_y, sums_l, sums_l_l, &

palm/trunk/SOURCE/netcdf.f90

@@ -23 +23 @@
 ! Current revisions:
 ! ------------------
-!
+! +qr, nr, prr
 !
 ! Former revisions:
@@ -466 +466 @@
 !
 !--             Most variables are defined on the scalar grid
-                CASE ( 'e', 'lpt', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c', &
-                       'ql_v', 'ql_vp', 'qv', 'rho', 's', 'sa', 'vpt' )
+                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',&
+                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',    &
+                       'rho', 's', 'sa', 'vpt' )
 
                    grid_x = 'x'
@@ -994 +995 @@
 !
 !--             Most variables are defined on the scalar grid
-                CASE ( 'e', 'lpt', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c', &
-                       'ql_v', 'ql_vp', 'qv', 'rho', 's', 'sa', 'vpt' )
+                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',   &
+                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho',&
+                       's', 'sa', 'vpt' )
 
                    grid_x = 'x'
@@ -1556 +1558 @@
 !
 !--                   Most variables are defined on the zu grid
-                      CASE ( 'e_xy', 'lpt_xy', 'p_xy', 'pc_xy', 'pr_xy', 'pt_xy',&
-                             'q_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy', 'ql_vp_xy',  &
-                             'qv_xy', 'rho_xy', 's_xy', 'sa_xy', 'vpt_xy' )
+                      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_xy', 's_xy', 'sa_xy', 'vpt_xy' )
 
                          grid_x = 'x'
@@ -2142 +2145 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_xz', 'lpt_xz', 'p_xz', 'pc_xz', 'pr_xz', 'pt_xz', &
-                          'q_xz', 'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz',   &
-                          'qv_xz', 'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' )
+                   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_xz', 's_xz', 'sa_xz', 'vpt_xz' )
 
                       grid_x = 'x'
@@ -2728 +2732 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_yz', 'lpt_yz', 'p_yz', 'pc_yz', 'pr_yz', 'pt_yz',&
-                          'q_yz', 'ql_yz', 'ql_c_yz', 'ql_v_yz', 'ql_vp_yz',  &
-                          'qv_yz', 'rho_yz', 's_yz', 'sa_yz', 'vpt_yz' )
+                   CASE ( 'e_yz', 'lpt_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_yz', 's_yz', 'sa_yz', 'vpt_yz' )
 
                       grid_x = 'x'
@@ -3915 +3920 @@
 !
 !--             Most variables are defined on the zu levels
-                CASE ( 'e', 'lpt', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c',  &
-                       'ql_v', 'ql_vp', 'qv', 'rho', 's', 'sa', 'u', 'v', 'vpt' )
+                CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',   &
+                       'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho',&
+                       's', 'sa', 'u', 'v', 'vpt' )
 
                    grid_z = 'zu'

palm/trunk/SOURCE/parin.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme and steering parameters:
+! +*_init, *_surface, *_surface_initial_change, *_vertical_gradient, 
+! +*_vertical_gradient_level, surface_waterflux_*, 
+! +cloud_scheme, drizzle, mu_constant, mu_constant_value, ventilation_effect
 !
 ! Former revisions:
@@ -194 +198 @@
              canyon_width_x, canyon_width_y, canyon_wall_left, &
              canyon_wall_south, cfl_factor, cloud_droplets, cloud_physics, &
-             collective_wait, conserve_volume_flow, conserve_volume_flow_mode, &
-             coupling_start_time, cthf, curvature_solution_effects, &
-             cycle_mg, damp_level_1d, dissipation_1d, & !dissipation_control, &
+             cloud_scheme, collective_wait, conserve_volume_flow, &
+             conserve_volume_flow_mode, coupling_start_time, cthf, &
+             curvature_solution_effects, cycle_mg, damp_level_1d, &
+             dissipation_1d, & !dissipation_control, &
              dp_external, dp_level_b, dp_smooth, dpdxy, drag_coefficient, &
-             dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, dz_max, & 
+             drizzle, dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, dz_max, & 
              dz_stretch_factor, dz_stretch_level, e_init, e_min, end_time_1d, &
              fft_method, galilei_transformation, humidity, &
@@ -208 +213 @@
              loop_optimization, masking_method, mg_cycles, &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, &
-             netcdf_precision, neutral, ngsrb, nsor, &
+             mu_constant, mu_constant_value, &
+             netcdf_precision, neutral, ngsrb, nr_surface, &
+             nr_surface_initial_change, nr_vertical_gradient, &
+             nr_vertical_gradient_level, nsor, &
              nsor_ini, nx, ny, nz, ocean, omega, omega_sor, &
              passive_scalar, pch_index, phi, plant_canopy, prandtl_layer, &
@@ -215 +223 @@
              pt_surface_initial_change, pt_vertical_gradient, &
              pt_vertical_gradient_level, q_surface, q_surface_initial_change, &
-             q_vertical_gradient, q_vertical_gradient_level, radiation, &
+             q_vertical_gradient, q_vertical_gradient_level, &
+             qr_surface, qr_surface_initial_change, &
+             qr_vertical_gradient, qr_vertical_gradient_level, radiation, &
              random_generator, random_heatflux, rayleigh_damping_factor, &
              rayleigh_damping_height, recycling_width, residual_limit, &
@@ -223 +233 @@
              statistic_regions, subs_vertical_gradient, &
              subs_vertical_gradient_level, surface_heatflux, surface_pressure, &
-             surface_scalarflux, surface_waterflux, s_surface, &
+             surface_scalarflux, surface_waterflux, surface_waterflux_nr, &
+             surface_waterflux_qr, s_surface, &
              s_surface_initial_change, s_vertical_gradient, &
              s_vertical_gradient_level, timestep_scheme, &
@@ -232 +243 @@
              use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke, &
              uv_heights, u_bulk, u_profile, vg_surface, vg_vertical_gradient, &
-             vg_vertical_gradient_level, v_bulk, v_profile, wall_adjustment, &
-             wall_heatflux, wall_humidityflux, wall_scalarflux, z0h_factor
-
-
+             vg_vertical_gradient_level, v_bulk, v_profile, ventilation_effect, &
+             wall_adjustment, wall_heatflux, wall_humidityflux, wall_scalarflux, &
+             z0h_factor
+      
     NAMELIST /d3par/  averaging_interval, averaging_interval_pr, &
              create_disturbances, &
@@ -391 +402 @@
 !--                  also check the allocate statement in routine read_var_list
              ALLOCATE( lad(0:nz+1),pt_init(0:nz+1), q_init(0:nz+1),           &
+                       qr_init(0:nz+1), nr_init(0:nz+1),                      &
                        sa_init(0:nz+1), ug(0:nz+1), u_init(0:nz+1),           &
                        v_init(0:nz+1), vg(0:nz+1),                            &

palm/trunk/SOURCE/prognostic_equations.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! implementation of two new prognostic equations for rain drop concentration (nr)
+! and rain water content (qr)
+!
+! currently, only available for cache loop optimization
 !
 ! Former revisions:
@@ -154 +157 @@
     USE diffusion_w_mod
     USE impact_of_latent_heat_mod
+    USE microphysics_mod
     USE plant_canopy_model_mod
     USE production_e_mod
@@ -409 +413 @@
 
 !
+!--      If required, calculate tendencies for total water content, rain water 
+!--      content, rain drop concentration and liquid temperature
+
+         IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+
+            tend_q(:,j,i)  = 0.0
+            tend_qr(:,j,i) = 0.0
+            tend_nr(:,j,i) = 0.0
+            tend_pt(:,j,i) = 0.0
+!
+!--         Droplet size distribution (dsd) properties are needed for the 
+!--         computation of selfcollection, breakup, evaporation and 
+!--         sedimentation of rain.  
+            IF ( precipitation )  THEN
+               CALL dsd_properties( i,j ) 
+               CALL autoconversion( i,j )
+               CALL accretion( i,j )
+               CALL selfcollection_breakup( i,j )
+               CALL evaporation_rain( i,j )
+               CALL sedimentation_rain( i,j )
+            ENDIF
+
+            IF ( drizzle )  CALL sedimentation_cloud( i,j )
+
+         ENDIF
+
+!
 !--       If required, compute prognostic equation for potential temperature
           IF ( .NOT. neutral )  THEN
@@ -433 +464 @@
              ENDIF
 
-!
-!--          If required compute impact of latent heat due to precipitation
-             IF ( precipitation )  THEN
-                CALL impact_of_latent_heat( i, j )
+!--          Using microphysical tendencies (latent heat)
+             IF ( cloud_physics )  THEN
+                IF ( icloud_scheme == 0 )  THEN
+                   tend(:,j,i) = tend(:,j,i) + tend_pt(:,j,i)
+                ELSEIF ( icloud_scheme == 1 .AND. precipitation)  THEN
+                   CALL impact_of_latent_heat( i, j ) 
+                ENDIF
              ENDIF
 
@@ -555 +589 @@
              ENDIF
              CALL diffusion_s( i, j, q, qsws, qswst, wall_qflux )
-       
-!
-!--          If required compute decrease of total water content due to
-!--          precipitation
-             IF ( precipitation )  THEN
-                CALL calc_precipitation( i, j )
-             ENDIF
-
+     
+!
+!--          Using microphysical tendencies
+             IF ( cloud_physics )  THEN
+                IF ( icloud_scheme == 0 )  THEN
+                   tend(:,j,i) = tend(:,j,i) + tend_q(:,j,i)
+                ELSEIF ( icloud_scheme == 1 .AND. precipitation )  THEN
+                   CALL calc_precipitation( i, j )
+                ENDIF
+             ENDIF
 !
 !--          Sink or source of scalar concentration due to canopy elements
              IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
 
+!
 !--          If required compute influence of large-scale subsidence/ascent
              IF ( large_scale_subsidence )  THEN
@@ -598 +635 @@
                    ENDDO
                 ENDIF
+             ENDIF
+
+!
+!--          If required, calculate prognostic equations for rain water content 
+!--          and rain drop concentration
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+!
+!--             Calculate prognostic equation for rain water content
+                tend(:,j,i) = 0.0
+                IF ( timestep_scheme(1:5) == 'runge' ) &
+                THEN
+                   IF ( ws_scheme_sca )  THEN
+                      CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       & 
+                                       diss_s_qr, flux_l_qr, diss_l_qr, &
+                                       i_omp_start, tn )
+                   ELSE 
+                      CALL advec_s_pw( i, j, qr )
+                   ENDIF
+                ELSE
+                   CALL advec_s_up( i, j, qr )
+                ENDIF
+                CALL diffusion_s( i, j, qr, qrsws, qrswst, wall_qrflux )
+
+!
+!--             Using microphysical tendencies (autoconversion, accretion, 
+!--             evaporation; if required: sedimentation)
+                tend(:,j,i) = tend(:,j,i) + tend_qr(:,j,i)
+
+!
+!--             If required, compute influence of large-scale subsidence/ascent
+                IF ( large_scale_subsidence )  THEN
+                   CALL subsidence( i, j, tend, qr, qr_init )
+                ENDIF
+
+!              CALL user_actions( i, j, 'qr-tendency' )
+
+!
+!--             Prognostic equation for rain water content
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +    &
+                                                     tsc(3) * tqr_m(k,j,i) )     &
+                                         - tsc(5) * rdf_sc(k) *                  &
+                                        ( qr(k,j,i) - qr_init(k) )
+                   IF ( qr_p(k,j,i) < 0.0 )  qr_p(k,j,i) = 0.1 * qr(k,j,i)
+                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_s_inner(j,i)+1, nzt
+                         tqr_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ELSEIF ( intermediate_timestep_count < &
+                            intermediate_timestep_count_max )  THEN
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tqr_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                        5.3125 * tqr_m(k,j,i)
+                      ENDDO
+                   ENDIF
+                ENDIF
+
+!
+!--             Calculate prognostic equation for rain drop concentration.
+                tend(:,j,i) = 0.0
+                IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( ws_scheme_sca )  THEN
+                      CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    & 
+                                    diss_s_nr, flux_l_nr, diss_l_nr, &
+                                    i_omp_start, tn )
+                   ELSE 
+                      CALL advec_s_pw( i, j, nr )
+                   ENDIF
+                ELSE
+                   CALL advec_s_up( i, j, nr )
+                ENDIF
+                CALL diffusion_s( i, j, nr, nrsws, nrswst, wall_nrflux )
+
+!--             Using microphysical tendencies (autoconversion, accretion, 
+!--             selfcollection, breakup, evaporation; 
+!--             if required: sedimentation)
+                tend(:,j,i) = tend(:,j,i) + tend_nr(:,j,i)
+
+!
+!--             If required, compute influence of large-scale subsidence/ascent
+                IF ( large_scale_subsidence )  THEN
+                   CALL subsidence( i, j, tend, nr, nr_init )
+                ENDIF
+
+!                CALL user_actions( i, j, 'nr-tendency' )
+
+!
+!--             Prognostic equation for rain drop concentration
+                DO  k = nzb_s_inner(j,i)+1, nzt
+                   nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
+                                                  tsc(3) * tnr_m(k,j,i) )     &
+                                      - tsc(5) * rdf_sc(k) *                  &
+                                        ( nr(k,j,i) - nr_init(k) )
+                   IF ( nr_p(k,j,i) < 0.0 )  nr_p(k,j,i) = 0.1 * nr(k,j,i)
+                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_s_inner(j,i)+1, nzt
+                         tnr_m(k,j,i) = tend(k,j,i)
+                      ENDDO
+                   ELSEIF ( intermediate_timestep_count < &
+                            intermediate_timestep_count_max )  THEN
+                      DO  k = nzb_s_inner(j,i)+1, nzt
+                         tnr_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                        5.3125 * tnr_m(k,j,i)
+                      ENDDO
+                   ENDIF
+                ENDIF
+
              ENDIF
 

palm/trunk/SOURCE/read_3d_binary.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme:
+! +prr, prr_av, *, *_av, *s, *sws, *swst
 !
 ! Former revisions:
@@ -279 +281 @@
 !--    First compare the version numbers
        READ ( 13 )  version_on_file
-       binary_version = '3.2'
+       binary_version = '3.3'
        IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
           WRITE( message_string, * ) 'version mismatch concerning data ',      &
@@ -457 +459 @@
                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
+                CASE ( 'prr' )
+                   IF ( .NOT. ALLOCATED( prr ) )  THEN
+                      ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   prr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'prr_av' )
+                   IF ( .NOT. ALLOCATED( prr_av ) )  THEN
+                      ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   prr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
                 CASE ( 'precipitation_amount' )
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
@@ -508 +526 @@
                    ql_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
                                        tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qr' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   qr(:,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
+                   nr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qr_av' )
+                   IF ( .NOT. ALLOCATED( qr_av ) )  THEN
+                      ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   qr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'nr_av' )
+                   IF ( .NOT. ALLOCATED( nr_av ) )  THEN
+                      ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                   ENDIF
+                   IF ( k == 1 )  READ ( 13 )  tmp_3d
+                   nr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qrs' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qrs(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qrsws' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qrsws(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'qrswst' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   qrswst(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'nrs' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   nrs(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'nrsws' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   nrsws(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+
+                CASE ( 'nrswst' )
+                   IF ( k == 1 )  READ ( 13 )  tmp_2d
+                   nrswst(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
+                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'ql_c_av' )

palm/trunk/SOURCE/read_var_list.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-!
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme:
+! +bc_*_b, +bc_*_t, +bc_*_t_val, *_init, *_surface, *_surface_initial_change,
+! +*_vertical_gradient, +*_vertical_gradient_level, *_vertical_gradient_level_ind,
+! +surface_waterflux_*
+!
+! in addition, steering parameters parameters of the two-moment cloud physics 
+! scheme:   
+! +cloud_scheme, +drizzle, +mu_constant, +mu_constant_value, +ventilation_effect
 !
 ! Former revisions:
@@ -169 +177 @@
 !-- Make version number check first
     READ ( 13 )  version_on_file
-    binary_version = '3.5'
+    binary_version = '3.6'
     IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
        WRITE( message_string, * ) 'version mismatch concerning control ', &
@@ -235 +243 @@
        ALLOCATE( lad(0:nz+1), ug(0:nz+1), u_init(0:nz+1), vg(0:nz+1),    &
                  v_init(0:nz+1), pt_init(0:nz+1), q_init(0:nz+1),        &
-                 sa_init(0:nz+1),                                        &
+                 qr_init(0:nz+1), nr_init(0:nz+1), sa_init(0:nz+1),      &
                  hom(0:nz+1,2,pr_palm+max_pr_user,0:statistic_regions),  &
                  hom_sum(0:nz+1,pr_palm+max_pr_user,0:statistic_regions) )
@@ -268 +276 @@
           CASE ( 'bc_ns' )
              READ ( 13 )  bc_ns
+          CASE ( 'bc_nr_b' )
+             READ ( 13 )  bc_nr_b
+          CASE ( 'bc_nr_t' )
+             READ ( 13 )  bc_nr_t
+          CASE ( 'bc_nr_t_val' )
+             READ ( 13 )  bc_nr_t_val
           CASE ( 'bc_p_b' )
              READ ( 13 )  bc_p_b
@@ -284 +298 @@
           CASE ( 'bc_q_t_val' )
              READ ( 13 )  bc_q_t_val
+          CASE ( 'bc_qr_b' )
+             READ ( 13 )  bc_qr_b
+          CASE ( 'bc_qr_t' )
+             READ ( 13 )  bc_qr_t
+          CASE ( 'bc_qr_t_val' )
+             READ ( 13 )  bc_qr_t_val
           CASE ( 'bc_s_b' )
              READ ( 13 )  bc_s_b
@@ -326 +346 @@
           CASE ( 'cloud_physics' )
              READ ( 13 )  cloud_physics
+          CASE ( 'cloud_scheme' )
+             READ ( 13 )  cloud_scheme
           CASE ( 'collective_wait' )
              READ ( 13 )  collective_wait
@@ -356 +378 @@
           CASE ( 'drag_coefficient' )
              READ ( 13 )  drag_coefficient
+          CASE ( 'drizzle' )
+             READ ( 13 )  drizzle
           CASE ( 'dt_pr_1d' )
              READ ( 13 )  dt_pr_1d
@@ -438 +462 @@
           CASE ( 'momentum_advec' )
              READ ( 13 )  momentum_advec
+          CASE ( 'mu_constant' )
+             READ ( 13 )  mu_constant
+          CASE ( 'mu_constant_value' )
+             READ ( 13 )  mu_constant_value
           CASE ( 'netcdf_precision' )
              READ ( 13 )  netcdf_precision
@@ -444 +472 @@
           CASE ( 'ngsrb' )
              READ ( 13 )  ngsrb
+          CASE ( 'nr_init' )
+             READ ( 13 )  nr_init
+          CASE ( 'nr_surface' )
+             READ ( 13 )  nr_surface
+          CASE ( 'nr_surface_initial_change' )
+             READ ( 13 )  nr_surface_initial_change
+          CASE ( 'nr_vertical_gradient' )
+             READ ( 13 )  nr_vertical_gradient
+          CASE ( 'nr_vertical_gradient_level' )
+             READ ( 13 )  nr_vertical_gradient_level
+          CASE ( 'nr_vertical_gradient_level_ind' )
+             READ ( 13 )  nr_vertical_gradient_level_ind
           CASE ( 'nsor' )
              READ ( 13 )  nsor
@@ -510 +550 @@
           CASE ( 'q_vertical_gradient_level_ind' )
              READ ( 13 )  q_vertical_gradient_level_ind
+          CASE ( 'qr_init' )
+             READ ( 13 )  qr_init
+          CASE ( 'qr_surface' )
+             READ ( 13 )  qr_surface
+          CASE ( 'qr_surface_initial_change' )
+             READ ( 13 )  qr_surface_initial_change
+          CASE ( 'qr_vertical_gradient' )
+             READ ( 13 )  qr_vertical_gradient
+          CASE ( 'qr_vertical_gradient_level' )
+             READ ( 13 )  qr_vertical_gradient_level
+          CASE ( 'qr_vertical_gradient_level_ind' )
+             READ ( 13 )  qr_vertical_gradient_level_ind
           CASE ( 'radiation' )
              READ ( 13 )  radiation
@@ -555 +607 @@
              READ ( 13 )  surface_scalarflux             
           CASE ( 'surface_waterflux' )
-             READ ( 13 )  surface_waterflux             
+             READ ( 13 )  surface_waterflux      
+          CASE ( 'surface_waterflux_nr' )
+             READ ( 13 )  surface_waterflux_nr    
+          CASE ( 'surface_waterflux_qr' )
+             READ ( 13 )  surface_waterflux_qr    
           CASE ( 's_surface' )
              READ ( 13 )  s_surface
@@ -656 +712 @@
           CASE ( 'v_max_ijk' )
              READ ( 13 )  v_max_ijk
+          CASE ( 'ventilation_effect' )
+             READ ( 13 )  ventilation_effect
           CASE ( 'vg' )
              READ ( 13 )  vg

palm/trunk/SOURCE/sum_up_3d_data.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! +nr, prr, qr
 !
 ! Former revisions:
@@ -115 +115 @@
                 lwp_av = 0.0
 
+             CASE ( 'nr' )
+                IF ( .NOT. ALLOCATED( nr_av ) )  THEN
+                   ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                nr_av = 0.0
+
              CASE ( 'p' )
                 IF ( .NOT. ALLOCATED( p_av ) )  THEN
@@ -133 +139 @@
                 pr_av = 0.0
 
+             CASE ( 'prr' )
+                IF ( .NOT. ALLOCATED( prr_av ) )  THEN
+                   ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                prr_av = 0.0
+
              CASE ( 'prr*' )
                 IF ( .NOT. ALLOCATED( precipitation_rate_av ) )  THEN
@@ -175 +187 @@
                 ql_vp_av = 0.0
 
+             CASE ( 'qr' )
+                IF ( .NOT. ALLOCATED( qr_av ) )  THEN
+                   ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                ENDIF
+                qr_av = 0.0
+
              CASE ( 'qsws*' )
                 IF ( .NOT. ALLOCATED( qsws_av ) )  THEN
@@ -301 +319 @@
                    lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i) * &
                                                     dzw(1:nzt+1) )
+                ENDDO
+             ENDDO
+
+          CASE ( 'nr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      nr_av(k,j,i) = nr_av(k,j,i) + nr(k,j,i)
+                   ENDDO
                 ENDDO
              ENDDO
@@ -423 +450 @@
              ENDDO
 
+          CASE ( 'qr' )
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   DO  k = nzb, nzt+1
+                      qr_av(k,j,i) = qr_av(k,j,i) + qr(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
           CASE ( 'qsws*' )
              DO  i = nxlg, nxrg

palm/trunk/SOURCE/swap_timelevel.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
+! swap of timelevels for nr, qr added
 ! 
-!
 ! Former revisions:
 ! -----------------
@@ -91 +91 @@
     ENDIF
     IF ( humidity  .OR.  passive_scalar )  THEN
-       q = q_p
+       q = q_p             
+       IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+          qr = qr_p
+          nr = nr_p
+       ENDIF
     ENDIF
 
@@ -116 +120 @@
           IF ( humidity  .OR.  passive_scalar )  THEN
              q => q_1;    q_p => q_2
+             IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+                qr => qr_1;    qr_p => qr_2
+                nr => nr_1;    nr_p => nr_2
+             ENDIF
           ENDIF
 
@@ -135 +143 @@
           IF ( humidity  .OR.  passive_scalar )  THEN
              q => q_2;    q_p => q_1
+             IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+                qr => qr_2;    qr_p => qr_1
+                nr => nr_2;    nr_p => nr_1
+             ENDIF
           ENDIF
 

palm/trunk/SOURCE/time_integration.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! exchange of ghost points for nr, qr added
 !
 ! Former revisions:
@@ -240 +240 @@
              CALL exchange_horiz( prho, nbgp )
           ENDIF
-          IF (humidity  .OR.  passive_scalar)  CALL exchange_horiz( q_p, nbgp )
+          IF (humidity  .OR.  passive_scalar)  THEN
+             CALL exchange_horiz( q_p, nbgp )
+             IF ( cloud_physics .AND. icloud_scheme == 0 )  THEN
+                CALL exchange_horiz( qr_p, nbgp )
+                CALL exchange_horiz( nr_p, nbgp )
+             ENDIF
+          ENDIF
           IF ( cloud_droplets )  THEN
              CALL exchange_horiz( ql, nbgp )

palm/trunk/SOURCE/timestep.f90

@@ -20 +20 @@
 ! Current revisions:
 ! ------------------
-!
+! timestep is reduced in two-moment cloud scheme according to the maximum 
+! terminal velocity of rain drops
 !
 ! Former revisions:
@@ -75 +76 @@
 
     USE arrays_3d
+    USE cloud_parameters
     USE control_parameters
     USE cpulog
@@ -301 +303 @@
 !--    step minus a reduction (cfl_factor) to be on the safe side.
 !--    The time step must not exceed the maximum allowed value.
-       dt_3d = cfl_factor * MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w )
+       dt_3d = cfl_factor * MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w,   &
+                                 dt_precipitation )
        dt_3d = MIN( dt_3d, dt_max )
 

palm/trunk/SOURCE/user_actions.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
+! +qr-tendency, nr-tendency
 !
 ! Former revisions:
@@ -136 +137 @@
 
 
+          CASE ( 'qr-tendency' )
+
+
+          CASE ( 'nr-tendency' )
+
+
           CASE DEFAULT
              message_string = 'unknown location "' // location // '"'
@@ -188 +195 @@
 
           CASE ( 'q-tendency' )
+
+
+          CASE ( 'qr-tendency' )
+
+
+          CASE ( 'nr-tendency' )
 
 

palm/trunk/SOURCE/write_3d_binary.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme:
+! +prr, prr_av, *, *_av, *s, *sws, *swst
 !
 ! Former revisions:
@@ -110 +112 @@
 !
 !-- Write arrays.
-    binary_version = '3.2'
+    binary_version = '3.3'
 
     WRITE ( 14 )  binary_version
@@ -145 +147 @@
        WRITE ( 14 )  'pr_av               ';  WRITE ( 14 )  pr_av
     ENDIF
+    IF ( ALLOCATED( prr ) )  THEN
+       WRITE ( 14 )  'prr                 ';  WRITE ( 14 )  prr
+    ENDIF
+    IF ( ALLOCATED( prr_av ) )  THEN
+       WRITE ( 14 )  'prr_av              ';  WRITE ( 14 )  prr_av
+    ENDIF
     IF ( ALLOCATED( precipitation_amount ) )  THEN
        WRITE ( 14 )  'precipitation_amount';  WRITE ( 14 )  precipitation_amount
@@ -166 +174 @@
              WRITE ( 14 )  'ql_av               ';  WRITE ( 14 )  ql_av
           ENDIF
+          IF ( icloud_scheme == 0 )  THEN
+             WRITE ( 14 )  'qr                  ';  WRITE ( 14 ) qr
+             WRITE ( 14 )  'nr                  ';  WRITE ( 14 ) nr
+             IF ( ALLOCATED( qr_av ) )  THEN
+                WRITE ( 14 )  'qr_av               ';  WRITE ( 14 )  qr_av
+             ENDIF
+             IF ( ALLOCATED( nr_av ) )  THEN
+                WRITE ( 14 )  'nr_av               ';  WRITE ( 14 )  nr_av
+             ENDIF
+             WRITE ( 14 )  'qrs                 ';  WRITE ( 14 ) qrs
+             WRITE ( 14 )  'qrsws               ';  WRITE ( 14 ) qrsws
+             WRITE ( 14 )  'qrswst              ';  WRITE ( 14 ) qrswst
+             WRITE ( 14 )  'nrs                 ';  WRITE ( 14 ) nrs
+             WRITE ( 14 )  'nrsws               ';  WRITE ( 14 ) nrsws
+             WRITE ( 14 )  'nrswst              ';  WRITE ( 14 ) nrswst
+          ENDIF
        ENDIF
        WRITE ( 14 )  'qs                  ';  WRITE ( 14 ) qs

palm/trunk/SOURCE/write_var_list.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-!
+! necessary expansions according to the two new prognostic equations (nr, qr) 
+! of the two-moment cloud physics scheme:
+! +bc_*_b, bc_*_t, bc_*_t_val, *_init, *_surface, *_surface_initial_change,
+! +*_vertical_gradient, *_vertical_gradient_level, *_vertical_gradient_level_ind,
+! +surface_waterflux_*
+!
+! in addition, steering parameters parameters of the two-moment cloud physics 
+! scheme:   
+! +cloud_scheme, +drizzle, +mu_constant, +mu_constant_value, +ventilation_effect
 !
 ! Former revisions:
@@ -151 +159 @@
 
 
-    binary_version = '3.5'
+    binary_version = '3.6'
 
     WRITE ( 14 )  binary_version
@@ -190 +198 @@
     WRITE ( 14 )  'bc_ns                         '
     WRITE ( 14 )  bc_ns
+    WRITE ( 14 )  'bc_nr_b                       '
+    WRITE ( 14 )  bc_nr_b
+    WRITE ( 14 )  'bc_nr_t                       '
+    WRITE ( 14 )  bc_nr_t
+    WRITE ( 14 )  'bc_nr_t_val                   '
+    WRITE ( 14 )  bc_nr_t_val
     WRITE ( 14 )  'bc_p_b                        '
     WRITE ( 14 )  bc_p_b
@@ -206 +220 @@
     WRITE ( 14 )  'bc_q_t_val                    '
     WRITE ( 14 )  bc_q_t_val
+    WRITE ( 14 )  'bc_qr_b                       '
+    WRITE ( 14 )  bc_qr_b
+    WRITE ( 14 )  'bc_qr_t                       '
+    WRITE ( 14 )  bc_qr_t
+    WRITE ( 14 )  'bc_qr_t_val                   '
+    WRITE ( 14 )  bc_qr_t_val
     WRITE ( 14 )  'bc_s_b                        '
     WRITE ( 14 )  bc_s_b
@@ -248 +268 @@
     WRITE ( 14 )  'cloud_physics                 '
     WRITE ( 14 )  cloud_physics
+    WRITE ( 14 )  'cloud_scheme                  '
+    WRITE ( 14 )  cloud_scheme
     WRITE ( 14 )  'collective_wait               '
     WRITE ( 14 )  collective_wait
@@ -278 +300 @@
     WRITE ( 14 )  'drag_coefficient              '
     WRITE ( 14 )  drag_coefficient
+    WRITE ( 14 )  'drizzle                       '
+    WRITE ( 14 )  drizzle
     WRITE ( 14 )  'dt_pr_1d                      '
     WRITE ( 14 )  dt_pr_1d
@@ -358 +382 @@
     WRITE ( 14 )  'momentum_advec                '
     WRITE ( 14 )  momentum_advec
+    WRITE ( 14 )  'mu_constant                   '
+    WRITE ( 14 )  mu_constant
+    WRITE ( 14 )  'mu_constant_value             '
+    WRITE ( 14 )  mu_constant_value
     WRITE ( 14 )  'netcdf_precision              '
     WRITE ( 14 )  netcdf_precision
@@ -364 +392 @@
     WRITE ( 14 )  'ngsrb                         '
     WRITE ( 14 )  ngsrb
+    WRITE ( 14 )  'nr_init                       '
+    WRITE ( 14 )  nr_init
+    WRITE ( 14 )  'nr_surface                    '
+    WRITE ( 14 )  nr_surface
+    WRITE ( 14 )  'nr_surface_initial_change     '
+    WRITE ( 14 )  nr_surface_initial_change
+    WRITE ( 14 )  'nr_vertical_gradient          '
+    WRITE ( 14 )  nr_vertical_gradient
+    WRITE ( 14 )  'nr_vertical_gradient_level    '
+    WRITE ( 14 )  nr_vertical_gradient_level
+    WRITE ( 14 )  'nr_vertical_gradient_level_ind'
+    WRITE ( 14 )  nr_vertical_gradient_level_ind
     WRITE ( 14 )  'nsor                          '
     WRITE ( 14 )  nsor
@@ -428 +468 @@
     WRITE ( 14 )  'q_vertical_gradient_level_ind '
     WRITE ( 14 )  q_vertical_gradient_level_ind
+    WRITE ( 14 )  'qr_init                       '
+    WRITE ( 14 )  qr_init
+    WRITE ( 14 )  'qr_surface                    '
+    WRITE ( 14 )  qr_surface
+    WRITE ( 14 )  'qr_surface_initial_change     '
+    WRITE ( 14 )  qr_surface_initial_change
+    WRITE ( 14 )  'qr_vertical_gradient          '
+    WRITE ( 14 )  qr_vertical_gradient
+    WRITE ( 14 )  'qr_vertical_gradient_level    '
+    WRITE ( 14 )  qr_vertical_gradient_level
+    WRITE ( 14 )  'qr_vertical_gradient_level_ind'
+    WRITE ( 14 )  qr_vertical_gradient_level_ind
     WRITE ( 14 )  'radiation                     '
     WRITE ( 14 )  radiation
@@ -474 +526 @@
     WRITE ( 14 )  'surface_waterflux             '
     WRITE ( 14 )  surface_waterflux    
+    WRITE ( 14 )  'surface_waterflux_nr          '
+    WRITE ( 14 )  surface_waterflux_nr    
+    WRITE ( 14 )  'surface_waterflux_qr          '
+    WRITE ( 14 )  surface_waterflux_qr    
     WRITE ( 14 )  's_surface                     '
     WRITE ( 14 )  s_surface
@@ -574 +630 @@
     WRITE ( 14 )  'v_max_ijk                     '
     WRITE ( 14 )  v_max_ijk
+    WRITE ( 14 )  'ventilation_effect            '
+    WRITE ( 14 )  ventilation_effect
     WRITE ( 14 )  'vg                            '
     WRITE ( 14 )  vg