Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 151)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 152)
@@ -755,4 +755,50 @@
 
 !
+!--    Calculate the initial volume flow at the right and north boundary
+       IF ( conserve_volume_flow )  THEN
+
+          volume_flow_initial_l = 0.0
+          volume_flow_area_l    = 0.0
+  
+          IF ( nxr == nx )  THEN
+             DO  j = nys, nyn
+                DO  k = nzb_2d(j,nx) + 1, nzt
+                   volume_flow_initial_l(1) = volume_flow_initial_l(1) + &
+                                              u(k,j,nx) * dzu(k)
+                   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
+
+          IF ( nyn == ny )  THEN
+             DO  i = nxl, nxr
+                DO  k = nzb_2d(ny,i) + 1, nzt  
+                   volume_flow_initial_l(2) = volume_flow_initial_l(2) + &
+                                              v(k,ny,i) * dzu(k)
+                   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
+
+#if defined( __parallel )
+          CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&
+                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
+          CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1),      &
+                              2, MPI_REAL, MPI_SUM, comm2d, ierr )
+#else
+          volume_flow_initial = volume_flow_initial_l
+          volume_flow_area    = volume_flow_area_l
+#endif  
+       ENDIF
+
+!
 !--    For the moment, perturbation pressure and vertical velocity are zero
        p = 0.0; w = 0.0
@@ -918,58 +964,4 @@
        ENDIF
 
-
-!
-!--    Read binary data from restart file
-          WRITE (9,*) 'before read_3d_binary'
-          CALL local_flush( 9 )
-       CALL read_3d_binary
-          WRITE (9,*) 'after read_3d_binary'
-          CALL local_flush( 9 )
-
-!
-!--    Calculate initial temperature field and other constants used in case
-!--    of a sloping surface
-       IF ( sloping_surface )  CALL init_slope
-
-!
-!--    Initialize new time levels (only done in order to set boundary values
-!--    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 ( ocean )  sa_p = sa
-
-    ELSE
-!
-!--    Actually this part of the programm should not be reached
-       IF ( myid == 0 )  PRINT*,'+++ init_3d_model: unknown initializing ', &
-                                                    'problem'
-       CALL local_stop
-    ENDIF
-
-
-    IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
-!
-!--    Initialize old timelevels needed for radiation boundary conditions
-       IF ( outflow_l )  THEN
-          u_m_l(:,:,:) = u(:,:,1:2)
-          v_m_l(:,:,:) = v(:,:,0:1)
-          w_m_l(:,:,:) = w(:,:,0:1)
-       ENDIF
-       IF ( outflow_r )  THEN
-          u_m_r(:,:,:) = u(:,:,nx-1:nx)
-          v_m_r(:,:,:) = v(:,:,nx-1:nx)
-          w_m_r(:,:,:) = w(:,:,nx-1:nx)
-       ENDIF
-       IF ( outflow_s )  THEN
-          u_m_s(:,:,:) = u(:,0:1,:)
-          v_m_s(:,:,:) = v(:,1:2,:)
-          w_m_s(:,:,:) = w(:,0:1,:)
-       ENDIF
-       IF ( outflow_n )  THEN
-          u_m_n(:,:,:) = u(:,ny-1:ny,:)
-          v_m_n(:,:,:) = v(:,ny-1:ny,:)
-          w_m_n(:,:,:) = w(:,ny-1:ny,:)
-       ENDIF
-
 !
 !--    Calculate the initial volume flow at the right and north boundary
@@ -1016,4 +1008,58 @@
           volume_flow_area    = volume_flow_area_l
 #endif  
+       ENDIF
+
+
+!
+!--    Read binary data from restart file
+          WRITE (9,*) 'before read_3d_binary'
+          CALL local_flush( 9 )
+       CALL read_3d_binary
+          WRITE (9,*) 'after read_3d_binary'
+          CALL local_flush( 9 )
+
+!
+!--    Calculate initial temperature field and other constants used in case
+!--    of a sloping surface
+       IF ( sloping_surface )  CALL init_slope
+
+!
+!--    Initialize new time levels (only done in order to set boundary values
+!--    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 ( ocean )  sa_p = sa
+
+    ELSE
+!
+!--    Actually this part of the programm should not be reached
+       IF ( myid == 0 )  PRINT*,'+++ init_3d_model: unknown initializing ', &
+                                                    'problem'
+       CALL local_stop
+    ENDIF
+
+
+    IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+!
+!--    Initialize old timelevels needed for radiation boundary conditions
+       IF ( outflow_l )  THEN
+          u_m_l(:,:,:) = u(:,:,1:2)
+          v_m_l(:,:,:) = v(:,:,0:1)
+          w_m_l(:,:,:) = w(:,:,0:1)
+       ENDIF
+       IF ( outflow_r )  THEN
+          u_m_r(:,:,:) = u(:,:,nx-1:nx)
+          v_m_r(:,:,:) = v(:,:,nx-1:nx)
+          w_m_r(:,:,:) = w(:,:,nx-1:nx)
+       ENDIF
+       IF ( outflow_s )  THEN
+          u_m_s(:,:,:) = u(:,0:1,:)
+          v_m_s(:,:,:) = v(:,1:2,:)
+          w_m_s(:,:,:) = w(:,0:1,:)
+       ENDIF
+       IF ( outflow_n )  THEN
+          u_m_n(:,:,:) = u(:,ny-1:ny,:)
+          v_m_n(:,:,:) = v(:,ny-1:ny,:)
+          w_m_n(:,:,:) = w(:,ny-1:ny,:)
        ENDIF
 
