Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 49)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 51)
@@ -1,4 +1,6 @@
 New:
 ---
+
+Wall functions for vertical walls now include diabatic conditions. New subroutine wall_fluxes. New 4D-array rif_wall.
 
 new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files
@@ -10,5 +12,7 @@
 samples added to the user interface which show how to add user-define time series quantities.
 
-check_open, check_parameters, header, init_3d_model, modules, netcdf, parin, user_interface
+Makefile, check_open, check_parameters, header, init_3d_model, modules, netcdf, parin, user_interface
+
+New: wall_fluxes
 
 
@@ -20,4 +24,6 @@
 
 d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.
+
+Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.
 
 check_parameters, data_output_ts, flow_statistics, init_3d_model, modules, parin, time_integration
Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 49)
+++ palm/trunk/SOURCE/Makefile	(revision 51)
@@ -6,8 +6,9 @@
 # default suffixes removed from the suffix list to avoid calling of m2c in
 # case of .mod files
+# +wall_fluxes
 #
 # Former revisions:
 # -----------------
-# $Id $
+# $Id$
 # RCS Log replace by Id keyword, revision history cleaned up
 #
@@ -55,5 +56,5 @@
 	swap_timelevel.f90 temperton_fft.f90 time_integration.f90 \
 	time_to_string.f90 timestep.f90 timestep_scheme_steering.f90 \
-	transpose.f90 user_interface.f90 write_3d_binary.f90 \
+	transpose.f90 user_interface.f90 wall_fluxes.f90 write_3d_binary.f90 \
 	write_compressed.f90 write_var_list.f90
 
@@ -84,5 +85,5 @@
 	temperton_fft.o time_integration.o time_to_string.o timestep.o \
 	timestep_scheme_steering.o transpose.o user_interface.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
 
 CC = cc
@@ -216,4 +217,5 @@
 transpose.o: modules.o
 user_interface.o: modules.o user_interface.f90
+wall_fluxes.o: modules.o
 write_3d_binary.o: modules.o random_function.o
 write_compressed.o: modules.o
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 49)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 51)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! "by_user" allowed as initializing action, -data_output_ts
+! "by_user" allowed as initializing action, -data_output_ts,
+! leapfrog with non-flat topography not allowed any more
 !
 ! Former revisions:
@@ -88,4 +89,7 @@
        IF ( momentum_advec /= 'pw-scheme' )  THEN
           WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
+       ENDIF
+       IF ( timestep_scheme(1:8) == 'leapfrog' )  THEN
+          WRITE( action, '(A,A)' )  'timestep_scheme = ', timestep_scheme
        ENDIF
        IF ( psolver == 'multigrid'  .OR.  psolver == 'sor' )  THEN
Index: palm/trunk/SOURCE/diffusion_u.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_u.f90	(revision 49)
+++ palm/trunk/SOURCE/diffusion_u.f90	(revision 51)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! wall functions now include diabatic conditions, call of routine wall_fluxes
 ! 
 ! Former revisions:
@@ -30,4 +30,6 @@
 ! ------------
 ! Diffusion term of the u-component
+! To do: additional damping (needed for non-cyclic bc) causes bad vectorization
+!        and slows down the speed on NEC about 5-10%
 !------------------------------------------------------------------------------!
 
@@ -55,8 +57,9 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp, usvs
+       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
        REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL, DIMENSION(nzb:nzt+1)      ::  usvs
        REAL, DIMENSION(:,:),   POINTER ::  usws
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
@@ -100,7 +103,11 @@
 !--          Wall functions at the north and south walls, respectively
              IF ( wall_u(j,i) /= 0.0 )  THEN
+
+!
+!--             Calculate the horizontal momentum flux u'v'
+                CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),  &
+                                  usvs, 1.0, 0.0, 0.0, 0.0 )
+
                 DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
-                   usvs   = kappa * u(k,j,i) / LOG( 0.5 * dy / z0(j,i) )
-                   usvs   = -usvs * ABS( usvs )
                    kmyp_x = 0.25 * &
                             ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
@@ -132,5 +139,5 @@
                                 + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
                                                   )                            &
-                                     + wall_u(j,i) * usvs                      &
+                                     + wall_u(j,i) * usvs(k)                   &
                                    ) * ddy
                 ENDDO
@@ -204,8 +211,9 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp, usvs
+       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
        REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL, DIMENSION(nzb:nzt+1)      ::  usvs
        REAL, DIMENSION(:,:),   POINTER ::  usws
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
@@ -246,7 +254,11 @@
 !--    Wall functions at the north and south walls, respectively
        IF ( wall_u(j,i) .NE. 0.0 )  THEN
+
+!
+!--       Calculate the horizontal momentum flux u'v'
+          CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i),  &
+                            usvs, 1.0, 0.0, 0.0, 0.0 )
+
           DO  k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i)
-             usvs   = kappa * u(k,j,i) / LOG( 0.5 * dy / z0(j,i) )
-             usvs   = -usvs * ABS( usvs )
              kmyp_x = 0.25 * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
              kmym_x = 0.25 * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
@@ -276,5 +288,5 @@
                                 + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx     &
                                                   )                            &
-                                     + wall_u(j,i) * usvs                      &
+                                     + wall_u(j,i) * usvs(k)                   &
                                    ) * ddy
           ENDDO
Index: palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_v.f90	(revision 49)
+++ palm/trunk/SOURCE/diffusion_v.f90	(revision 51)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! wall functions now include diabatic conditions, call of routine wall_fluxes
 !
 ! Former revisions:
@@ -55,8 +55,9 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp, vsus
+       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
        REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL, DIMENSION(nzb:nzt+1)      ::  vsus
        REAL, DIMENSION(:,:),   POINTER ::  vsws
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
@@ -100,7 +101,11 @@
 !--          Wall functions at the left and right walls, respectively
              IF ( wall_v(j,i) /= 0.0 )  THEN
+
+!
+!--             Calculate the horizontal momentum flux v'u'
+                CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), &
+                                  vsus, 0.0, 1.0, 0.0, 0.0 )
+
                 DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
-                   vsus   = kappa * v(k,j,i) / LOG( 0.5 * dx / z0(j,i))
-                   vsus   = -vsus * ABS( vsus )
                    kmxp_x = 0.25 * &
                             ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
@@ -132,5 +137,5 @@
                                 + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
                                                   )                            &
-                                     + wall_v(j,i) * vsus                      &
+                                     + wall_v(j,i) * vsus(k)                   &
                                    ) * ddx
                 ENDDO
@@ -204,8 +209,9 @@
 
        INTEGER ::  i, j, k
-       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp, vsus
+       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
        REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL, DIMENSION(nzb:nzt+1)      ::  vsus
        REAL, DIMENSION(:,:),   POINTER ::  vsws
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
@@ -245,7 +251,11 @@
 !--    Wall functions at the left and right walls, respectively
        IF ( wall_v(j,i) /= 0.0 )  THEN
+
+!
+!--       Calculate the horizontal momentum flux v'u'
+          CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), &
+                            vsus, 0.0, 1.0, 0.0, 0.0 )
+
           DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
-             vsus   = kappa * v(k,j,i) / LOG( 0.5 * dx / z0(j,i))
-             vsus   = -vsus * ABS( vsus )
              kmxp_x = 0.25 * &
                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
@@ -277,5 +287,5 @@
                                 + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy     &
                                                   )                            &
-                                     + wall_v(j,i) * vsus                      &
+                                     + wall_v(j,i) * vsus(k)                   &
                                    ) * ddx
           ENDDO
Index: palm/trunk/SOURCE/diffusion_w.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_w.f90	(revision 49)
+++ palm/trunk/SOURCE/diffusion_w.f90	(revision 51)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! wall functions now include diabatic conditions, call of routine wall_fluxes
 !
 ! Former revisions:
@@ -54,9 +54,10 @@
        INTEGER ::  i, j, k
        REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
-                   kmyp_z, wsus, wsvs
+                   kmyp_z
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
                    km_damp_y(nys-1:nyn+1)
        REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL, DIMENSION(nzb:nzt+1)      ::  wsus, wsvs
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
 
@@ -113,17 +114,24 @@
 !--          Wall functions at all vertical walls, where necessary
              IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
+
+!
+!--             Calculate the horizontal momentum fluxes w'u' and/or w'v'
+                IF ( wall_w_x(j,i) /= 0.0 )  THEN
+                   CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1,              &
+                                     nzb_w_outer(j,i), wsus, 0.0, 0.0, 0.0, &
+                                     1.0 )
+                ELSE
+                   wsus = 0.0
+                ENDIF
+
+                IF ( wall_w_y(j,i) /= 0.0 )  THEN
+                   CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1,              &
+                                     nzb_w_outer(j,i), wsvs, 0.0, 0.0, 1.0, &
+                                     0.0 )
+                ELSE
+                   wsvs = 0.0
+                ENDIF
+
                 DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
-                   IF ( wall_w_x(j,i) /= 0.0 )  THEN
-                      wsus = kappa * w(k,j,i) / LOG( 0.5 * dx / z0(j,i))
-                      wsus = -wsus * ABS( wsus )
-                   ELSE
-                      wsus = 0.0
-                   ENDIF
-                   IF ( wall_w_y(j,i) /= 0.0 )  THEN
-                      wsvs = kappa * w(k,j,i) / LOG( 0.5 * dy / z0(j,i))
-                      wsvs = -wsvs * ABS( wsvs )
-                   ELSE
-                      wsvs = 0.0
-                   ENDIF
 !
 !--                Interpolate eddy diffusivities on staggered gridpoints
@@ -163,5 +171,5 @@
                           + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
                                                    )                           &
-                                     + wall_w_x(j,i) * wsus                    &
+                                     + wall_w_x(j,i) * wsus(k)                 &
                                    ) * ddx                                     &
                                  + (   fwyp(j,i) * (                           &
@@ -173,5 +181,5 @@
                           + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
                                                    )                           &
-                                     + wall_w_y(j,i) * wsvs                    &
+                                     + wall_w_y(j,i) * wsvs(k)                 &
                                    ) * ddy                                     &
                                  + 2.0 * (                                     &
@@ -202,9 +210,10 @@
        INTEGER ::  i, j, k
        REAL    ::  kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, &
-                   kmyp_z, wsus, wsvs
+                   kmyp_z
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
                    km_damp_y(nys-1:nyn+1)
        REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+       REAL, DIMENSION(nzb:nzt+1)      ::  wsus, wsvs
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
 
@@ -255,17 +264,22 @@
 !--    Wall functions at all vertical walls, where necessary
        IF ( wall_w_x(j,i) /= 0.0  .OR.  wall_w_y(j,i) /= 0.0 )  THEN
+
+!
+!--       Calculate the horizontal momentum fluxes w'u' and/or w'v'
+          IF ( wall_w_x(j,i) /= 0.0 )  THEN
+             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i), &
+                               wsus, 0.0, 0.0, 0.0, 1.0 )
+          ELSE
+             wsus = 0.0
+          ENDIF
+
+          IF ( wall_w_y(j,i) /= 0.0 )  THEN
+             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),  &
+                               wsvs, 0.0, 0.0, 1.0, 0.0 )
+          ELSE
+             wsvs = 0.0
+          ENDIF
+
           DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
-             IF ( wall_w_x(j,i) /= 0.0 )  THEN
-                wsus = kappa * w(k,j,i) / LOG( 0.5 * dx / z0(j,i))
-                wsus = -wsus * ABS( wsus )
-             ELSE
-                wsus = 0.0
-             ENDIF
-             IF ( wall_w_y(j,i) /= 0.0 )  THEN
-                wsvs = kappa * w(k,j,i) / LOG( 0.5 * dy / z0(j,i))
-                wsvs = -wsvs * ABS( wsvs )
-             ELSE
-                wsvs = 0.0
-             ENDIF
 !
 !--          Interpolate eddy diffusivities on staggered gridpoints
@@ -301,5 +315,5 @@
                           + kmxm_z * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1) &
                                                    )                           &
-                                     + wall_w_x(j,i) * wsus                    &
+                                     + wall_w_x(j,i) * wsus(k)                 &
                                    ) * ddx                                     &
                                  + (   fwyp(j,i) * (                           &
@@ -311,5 +325,5 @@
                           + kmym_z * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1) &
                                                    )                           &
-                                     + wall_w_y(j,i) * wsvs                    &
+                                     + wall_w_y(j,i) * wsvs(k)                 &
                                    ) * ddy                                     &
                                  + 2.0 * (                                     &
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 49)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 51)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Collection of time series quantities moved from routine flow_statistics to
+! here, routine user_statistics is called for each statistic region
 !
 ! Former revisions:
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 49)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 51)
@@ -8,5 +8,7 @@
 ! -----------------
 ! New initializing action "by_user" calls user_init_3d_model,
-! ts_value is allocated
+! ts_value is allocated, +module netcdf_control,
+! initial velocities at nzb+1 are regarded for volume
+! flow control in case they have been set zero before (to avoid small timesteps)
 !
 ! Former revisions:
@@ -44,4 +46,5 @@
     USE interfaces
     USE model_1d
+    USE netcdf_control
     USE particle_attributes
     USE pegrid
@@ -183,4 +186,19 @@
 
 !
+!-- 4D-array for storing the Rif-values at vertical walls
+    IF ( topography /= 'flat' )  THEN
+       ALLOCATE( rif_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1,1:4) )
+       rif_wall = 0.0
+    ENDIF
+
+!
+!-- Velocities at nzb+1 needed for volume flow control
+    IF ( conserve_volume_flow )  THEN
+       ALLOCATE( u_nzb_p1_for_vfc(nys:nyn), v_nzb_p1_for_vfc(nxl:nxr) )
+       u_nzb_p1_for_vfc = 0.0
+       v_nzb_p1_for_vfc = 0.0
+    ENDIF
+
+!
 !-- Initial assignment of the pointers
     IF ( timestep_scheme(1:5) /= 'runge' )  THEN
@@ -382,8 +400,9 @@
           ENDDO
 !
-!--       Set initial horizontal velocities at the lowest grid levels to zero
-!--       in order to avoid too small time steps caused by the diffusion
+!--       Set initial horizontal velocities at the lowest computational grid levels
+!--       to zero in order to avoid too small time steps caused by the diffusion
 !--       limit in the initial phase of a run (at k=1, dz/2 occurs in the
-!--       limiting formula!)
+!--       limiting formula!). The original values are stored to be later used for
+!--       volume flow control.
           DO  i = nxl-1, nxr+1
              DO  j = nys-1, nyn+1
@@ -392,4 +411,18 @@
              ENDDO
           ENDDO
+          IF ( conserve_volume_flow )  THEN
+             IF ( nxr == nx )  THEN
+                DO  j = nys, nyn
+                   k = nzb_u_inner(j,nx) + 1
+                   u_nzb_p1_for_vfc(j) = u_init(k) * dzu(k)
+                ENDDO
+             ENDIF
+             IF ( nyn == ny )  THEN
+                DO  i = nxl, nxr
+                   k = nzb_v_inner(ny,i) + 1
+                   v_nzb_p1_for_vfc(i) = v_init(k) * dzu(k)
+                ENDDO
+             ENDIF
+          ENDIF
 
           IF ( moisture  .OR.  passive_scalar )  THEN
@@ -568,4 +601,8 @@
                    volume_flow_area_l(1)    = volume_flow_area_l(1) + dzu(k)
                 ENDDO
+!
+!--             Correction if velocity at nzb+1 has been set zero further above
+                volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
+                                           u_nzb_p1_for_vfc(j)
              ENDDO
           ENDIF
@@ -578,4 +615,8 @@
                    volume_flow_area_l(2)    = volume_flow_area_l(2) + dzu(k)
                 ENDDO
+!
+!--             Correction if velocity at nzb+1 has been set zero further above
+                volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
+                                           v_nzb_p1_for_vfc(i)
              ENDDO
           ENDIF
@@ -704,9 +745,5 @@
 !
 !-- If required, initialize dvrp-software
-!    WRITE ( 9, * ) '*** myid=', myid, ' vor init_dvrp'
-!    CALL FLUSH_( 9 )
     IF ( dt_dvrp /= 9999999.9 )  CALL init_dvrp
-!    WRITE ( 9, * ) '*** myid=', myid, ' nach init_dvrp'
-!    CALL FLUSH_( 9 )
 
 !
@@ -719,9 +756,5 @@
 !
 !-- If required, initialize particles
-!    WRITE ( 9, * ) '*** myid=', myid, ' vor init_particles'
-!    CALL FLUSH_( 9 )
     CALL init_particles
-!    WRITE ( 9, * ) '*** myid=', myid, ' nach init_particles'
-!    CALL FLUSH_( 9 )
 
 !
@@ -825,6 +858,16 @@
 
 !
-!-- User-defined initializing actions
+!-- User-defined initializing actions. Check afterwards, if maximum number
+!-- of allowed timeseries is not exceeded
     CALL user_init
+
+    IF ( dots_num > dots_max )  THEN
+       IF ( myid == 0 )  THEN
+          PRINT*, '+++ user_init: number of time series quantities exceeds', &
+                  ' its maximum of dots_max = ', dots_max
+          PRINT*, '    Please increase dots_max in modules.f90.'
+       ENDIF
+       CALL local_stop
+    ENDIF
 
 !
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 49)
+++ palm/trunk/SOURCE/modules.f90	(revision 51)
@@ -5,6 +5,7 @@
 ! Actual revisions:
 ! -----------------
+! +array rif_wall
 ! +netcdf_64bit_3d, zu_s_inner, zw_w_inner, id_var_zusi_*, id_var_zwwi_*,
-! ts_value
+! ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc
 ! -data_output_ts, dots_n
 ! arrays dots_label and dots_unit now dimensioned with dots_max
@@ -77,6 +78,6 @@
     REAL, DIMENSION(:), ALLOCATABLE ::                                         &
           ddzu, dd2zu, dzu, ddzw, dzw, km_damp_x, km_damp_y, l_grid, pt_init,  &
-          q_init, rdf, ug, u_init, uvmean_outflow, uvmean_outflow_l, vg,       &
-          v_init, zu, zw
+          q_init, rdf, ug, uvmean_outflow, uvmean_outflow_l, u_init,           &
+          u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, zu, zw
 
     REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
@@ -106,4 +107,6 @@
           ql_c, te_m, tpt_m, tq_m, tu_m, tv_m, tw_m, u, u_m, u_p, v, v_m, v_p, &
           vpt, vpt_m, w, w_m, w_p
+
+    REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall
                                                     
 
Index: palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/read_3d_binary.f90	(revision 49)
+++ palm/trunk/SOURCE/read_3d_binary.f90	(revision 51)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +rif_wall
 !
 ! Former revisions:
@@ -267,4 +267,6 @@
           CASE ( 'rif_m' )
              READ ( 13 )  rif_m
+          CASE ( 'rif_wall' )
+             READ ( 13 )  rif_wall
           CASE ( 's_av' )
              ALLOCATE( s_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
Index: palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- palm/trunk/SOURCE/write_3d_binary.f90	(revision 49)
+++ palm/trunk/SOURCE/write_3d_binary.f90	(revision 51)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +rif_wall
 !
 ! Former revisions:
@@ -137,4 +137,7 @@
     IF ( timestep_scheme(1:5) /= 'runge' )  THEN
        WRITE ( 14 )  'rif_m               ';  WRITE ( 14 )  rif_m
+    ENDIF
+    IF ( topography /= 'flat' )  THEN
+       WRITE ( 14 )  'rif_wall            ';  WRITE ( 14 )  rif_wall
     ENDIF
     IF ( ALLOCATED( s_av ) )  THEN
