Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 1917)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 1918)
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! setting of a fixed reference state ('single_value') for ocean runs removed
 ! 
 ! Former revisions:
@@ -1526,11 +1526,4 @@
 
 !
-!-- Ocean runs always use reference values in the buoyancy term
-    IF ( ocean )  THEN
-       reference_state = 'single_value'
-       use_single_reference_value = .TRUE.
-    ENDIF
-
-!
 !-- Sign of buoyancy/stability terms
     IF ( ocean )  atmos_ocean_sign = -1.0_wp
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 1917)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 1918)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! in case of Wicker-Skamarock scheme, calculate disturbance kinetic energy here,
+! if flow_statistics is called before the first initial time step
 ! 
 ! Former revisions:
@@ -630,5 +631,6 @@
 !--    Computation of statistics when ws-scheme is not used. Else these 
 !--    quantities are evaluated in the advection routines.
-       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 )  THEN
+       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
+       THEN
           !$OMP DO
           DO  i = nxl, nxr
@@ -1623,6 +1625,7 @@
 
 !
-!-- If required, sum up horizontal averages for subsequent time averaging
-    IF ( do_sum )  THEN
+!-- If required, sum up horizontal averages for subsequent time averaging.
+!-- Do not sum, if flow statistics is called before the first initial time step.
+    IF ( do_sum  .AND.  simulated_time /= 0.0_wp )  THEN
        IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
        hom_sum = hom_sum + hom(:,1,:,:)
@@ -2230,5 +2233,6 @@
 !--    Computation of statistics when ws-scheme is not used. Else these
 !--    quantities are evaluated in the advection routines.
-       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0 )  THEN
+       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
+       THEN
 
           !$OMP DO
@@ -3566,5 +3570,6 @@
 !
 !-- If required, sum up horizontal averages for subsequent time averaging
-    IF ( do_sum )  THEN
+!-- Do not sum, if flow statistics is called before the first initial time step.
+    IF ( do_sum  .AND.  simulated_time /= 0.0_wp )  THEN
        IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
        hom_sum = hom_sum + hom(:,1,:,:)
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 1917)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 1918)
@@ -19,5 +19,7 @@
 ! Current revisions:
 ! ------------------
-! 
+! intermediate_timestep_count is set 0 instead 1 for first call of pres,
+! bugfix: initialization of local sum arrays are moved to the beginning of the
+!         routine because otherwise results from pres are overwritten
 ! 
 ! Former revisions:
@@ -709,7 +711,20 @@
     weight_substep = 1.0_wp
     weight_pres    = 1.0_wp
-    intermediate_timestep_count = 1  ! needed when simulated_time = 0.0
+    intermediate_timestep_count = 0  ! needed when simulated_time = 0.0
        
     CALL location_message( 'finished', .TRUE. )
+
+!
+!-- Initialize local summation arrays for routine flow_statistics.
+!-- This is necessary because they may not yet have been initialized when they
+!-- are called from flow_statistics (or - depending on the chosen model run -
+!-- are never initialized)
+    sums_divnew_l      = 0.0_wp
+    sums_divold_l      = 0.0_wp
+    sums_l_l           = 0.0_wp
+    sums_up_fraction_l = 0.0_wp
+    sums_wsts_bc_l     = 0.0_wp
+
+
 !
 !-- Initialize model variables
@@ -1821,15 +1836,4 @@
        ENDDO
     ENDIF
-
-!
-!-- Initialize local summation arrays for routine flow_statistics.
-!-- This is necessary because they may not yet have been initialized when they
-!-- are called from flow_statistics (or - depending on the chosen model run -
-!-- are never initialized)
-    sums_divnew_l      = 0.0_wp
-    sums_divold_l      = 0.0_wp
-    sums_l_l           = 0.0_wp
-    sums_up_fraction_l = 0.0_wp
-    sums_wsts_bc_l     = 0.0_wp
 
 !
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 1917)
+++ palm/trunk/SOURCE/modules.f90	(revision 1918)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! ------------------
-! 
+! default timestep switched from -1.0 to +1.0 in order to avoid wrong sign of
+! initially calculated divergence
 ! 
 ! Former revisions:
@@ -751,5 +752,5 @@
                  dt_max = 20.0_wp, &
                  dt_restart = 9999999.9_wp, &
-                 dt_run_control = 60.0_wp, dt_3d = -1.0_wp, dz = -1.0_wp, &
+                 dt_run_control = 60.0_wp, dt_3d = 1.0_wp, dz = -1.0_wp, &
                  dz_max = 9999999.9_wp, dz_stretch_factor = 1.08_wp, &
                  dz_stretch_level = 100000.0_wp, e_init = 0.0_wp, e_min = 0.0_wp, &
Index: palm/trunk/SOURCE/pres.f90
===================================================================
--- palm/trunk/SOURCE/pres.f90	(revision 1917)
+++ palm/trunk/SOURCE/pres.f90	(revision 1918)
@@ -19,5 +19,11 @@
 ! Current revisions:
 ! ------------------
-! 
+! sum of divergence is also calculated when pres is called before the initial
+! first time step,
+! stearing is modified by using intermediate_timestep_count = 0 in order to
+! determine if pres is called before the first initial timestep,
+! bugfix: in case of Neumann conditions for pressure both at bottom and top,
+!         mean vertical velocity is also removed for the first time step
+! bugfix for calculating divergences
 ! 
 ! Former revisions:
@@ -114,6 +120,6 @@
                nest_domain, nest_bound_l, nest_bound_n, nest_bound_r,          &
                nest_bound_s, on_device, outflow_l, outflow_n, outflow_r,       &
-               outflow_s, psolver, simulated_time, subdomain_size, topography, &
-               volume_flow, volume_flow_area, volume_flow_initial
+               outflow_s, psolver, subdomain_size, topography, volume_flow,    &
+               volume_flow_area, volume_flow_initial
 
     USE cpulog,                                                                &
@@ -149,7 +155,8 @@
 
     REAL(wp)     ::  ddt_3d         !<
+    REAL(wp)     ::  d_weight_pres  !<
     REAL(wp)     ::  localsum       !<
     REAL(wp)     ::  threadsum      !<
-    REAL(wp)     ::  d_weight_pres  !<
+    REAL(wp)     ::  weight_pres_l  !<
 
     REAL(wp), DIMENSION(1:3)   ::  volume_flow_l       !<
@@ -162,6 +169,16 @@
 
 
+!
+!-- Calculate quantities to be used locally
     ddt_3d = 1.0_wp / dt_3d
-    d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
+    IF ( intermediate_timestep_count == 0 )  THEN
+!
+!--    If pres is called before initial time step
+       weight_pres_l = 1.0_wp
+       d_weight_pres = 1.0_wp
+    ELSE
+       weight_pres_l = weight_pres(intermediate_timestep_count)
+       d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
+    ENDIF
 
 !
@@ -187,5 +204,5 @@
        ENDIF
        
-    ELSEIF ( psolver == 'sor'  .AND.  intermediate_timestep_count == 1 )  THEN
+    ELSEIF ( psolver == 'sor'  .AND.  intermediate_timestep_count <= 1 )  THEN
 
 !
@@ -283,33 +300,35 @@
 !
 !-- Remove mean vertical velocity in case that Neumann conditions are
-!-- used both at bottom and top boundary, and if not a nested domain
-    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.  .NOT. nest_domain)  THEN
-       IF ( simulated_time > 0.0_wp )  THEN ! otherwise nzb_w_inner not yet known
-          w_l = 0.0_wp;  w_l_l = 0.0_wp
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb_w_inner(j,i)+1, nzt
-                   w_l_l(k) = w_l_l(k) + w(k,j,i)
-                ENDDO
-             ENDDO
-          ENDDO
+!-- used both at bottom and top boundary, and if not a nested domain.
+!-- This cannot be done before the first initial time step because ngp_2dh_outer
+!-- is not yet known then.
+    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.  .NOT. nest_domain  .AND.    &
+         intermediate_timestep_count /= 0 )                                    &
+    THEN
+       w_l = 0.0_wp;  w_l_l = 0.0_wp
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_w_inner(j,i)+1, nzt
+                w_l_l(k) = w_l_l(k) + w(k,j,i)
+             ENDDO
+          ENDDO
+       ENDDO
 #if defined( __parallel )   
-          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
-                              comm2d, ierr )
+       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
+       CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
+                           comm2d, ierr )
 #else
-          w_l = w_l_l  
+       w_l = w_l_l
 #endif
-          DO  k = 1, nzt
-             w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
-          ENDDO
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-                DO  k = nzb_w_inner(j,i)+1, nzt
-                   w(k,j,i) = w(k,j,i) - w_l(k)
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
+       DO  k = 1, nzt
+          w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
+       ENDDO
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
+             DO  k = nzb_w_inner(j,i)+1, nzt
+                w(k,j,i) = w(k,j,i) - w_l(k)
+             ENDDO
+          ENDDO
+       ENDDO
     ENDIF
 
@@ -363,7 +382,7 @@
     ENDDO
 
-    IF ( intermediate_timestep_count == intermediate_timestep_count_max )  THEN
-       localsum = localsum + threadsum * dt_3d * &
-                             weight_pres(intermediate_timestep_count)
+    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
+         intermediate_timestep_count == 0 )  THEN
+       localsum = localsum + threadsum * dt_3d * weight_pres_l
     ENDIF
     !$OMP END PARALLEL
@@ -390,5 +409,6 @@
 !-- Compute possible PE-sum of divergences for flow_statistics. Carry out
 !-- computation only at last Runge-Kutta substep.
-    IF ( intermediate_timestep_count == intermediate_timestep_count_max )  THEN
+    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
+         intermediate_timestep_count == 0 )  THEN
        !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
        !$OMP DO SCHEDULE( STATIC )
@@ -402,6 +422,5 @@
        ENDDO
        !$acc end parallel loop
-       localsum = localsum + threadsum * dt_3d * &
-                             weight_pres(intermediate_timestep_count)
+       localsum = localsum + threadsum * dt_3d * weight_pres_l
        !$OMP END PARALLEL
     ENDIF
@@ -411,6 +430,8 @@
 !-- For completeness, set the divergence sum of all statistic regions to those
 !-- of the total domain
-    IF ( intermediate_timestep_count == intermediate_timestep_count_max )      &
+    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
+         intermediate_timestep_count == 0 )  THEN
        sums_divold_l(0:statistic_regions) = localsum
+    ENDIF
 
     CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
@@ -564,5 +585,5 @@
 !-- Store perturbation pressure on array p, used for pressure data output.
 !-- Ghost layers are added in the output routines (except sor-method: see below)
-    IF ( intermediate_timestep_count == 1 )  THEN
+    IF ( intermediate_timestep_count <= 1 )  THEN
        !$OMP PARALLEL PRIVATE (i,j,k)
        !$OMP DO
@@ -616,5 +637,5 @@
     !$OMP PARALLEL PRIVATE (i,j,k)
     !$OMP DO
-    !$acc kernels present( ddzu, nzb_u_inner, nzb_v_inner, nzb_w_inner, tend, u, v, w, weight_pres )
+    !$acc kernels present( ddzu, nzb_u_inner, nzb_v_inner, nzb_w_inner, tend, u, v, w )
     !$acc loop independent
     DO  i = nxl, nxr   
@@ -626,5 +647,5 @@
                 w(k,j,i) = w(k,j,i) - dt_3d *                                 &
                            ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) *      &
-                           weight_pres(intermediate_timestep_count)
+                           weight_pres_l
              ENDIF
           ENDDO
@@ -634,5 +655,5 @@
                 u(k,j,i) = u(k,j,i) - dt_3d *                                 &
                            ( tend(k,j,i) - tend(k,j,i-1) ) * ddx *            &
-                           weight_pres(intermediate_timestep_count)
+                           weight_pres_l
              ENDIF
           ENDDO
@@ -642,5 +663,5 @@
                 v(k,j,i) = v(k,j,i) - dt_3d *                                 &
                            ( tend(k,j,i) - tend(k,j-1,i) ) * ddy *            &
-                           weight_pres(intermediate_timestep_count)
+                           weight_pres_l
              ENDIF
           ENDDO                                                         
@@ -737,5 +758,6 @@
 !-- A possible PE-sum is computed in flow_statistics. Carry out computation
 !-- only at last Runge-Kutta step.
-    IF ( intermediate_timestep_count == intermediate_timestep_count_max )  THEN 
+    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
+         intermediate_timestep_count == 0 )  THEN
        CALL cpu_log( log_point_s(1), 'divergence', 'start' )
        sums_divnew_l = 0.0_wp
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 1917)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 1918)
@@ -19,5 +19,7 @@
 ! Current revisions:
 ! ------------------
-! 
+! determination of time step moved to the end of the time step loop,
+! the first time step is now always calculated before the time step loop (i.e.
+! also in case of restart runs)
 ! 
 ! Former revisions:
@@ -260,8 +262,10 @@
     USE indices,                                                               &
         ONLY:  i_left, i_right, j_north, j_south, nbgp, nx, nxl, nxlg, nxr,    &
-               nxrg, nyn, nys, nzb, nzt, nzb_u_inner, nzb_v_inner 
+               nxrg, nyn, nyng, nys, nysg, nzb, nzt, nzb_u_inner, nzb_v_inner
 
     USE interaction_droplets_ptq_mod,                                          &
         ONLY:  interaction_droplets_ptq
+
+    USE interfaces
 
     USE kinds
@@ -311,5 +315,6 @@
 
     USE statistics,                                                            &
-        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l
+        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l, u_max,         &
+               u_max_ijk, v_max, v_max_ijk, w_max, w_max_ijk
 
     USE surface_layer_fluxes_mod,                                              &
@@ -326,8 +331,10 @@
     CHARACTER (LEN=9) ::  time_to_string          !< 
 
-!
-!-- At the beginning of a simulation determine the time step as well as
-!-- determine and print out the run control parameters
-    IF ( simulated_time == 0.0_wp )  CALL timestep
+    REAL(wp) ::  dt_3d_old  !< temporary storage of timestep to be used for
+                            !< steering of run control output interval
+
+!
+!-- At beginning determine the first time step
+    CALL timestep
 
 !
@@ -340,4 +347,9 @@
     ENDIF
 
+!
+!-- Determine and print out the run control quantities before the first time
+!-- step of this run. For the initial run, some statistics (e.g. divergence)
+!-- need to be determined first.
+    IF ( simulated_time == 0.0_wp )  CALL flow_statistics
     CALL run_control
 
@@ -349,9 +361,10 @@
 !--    In case of model termination initiated by the local model the coupler 
 !--    must not be called because this would again cause an MPI hang. 
-       DO WHILE ( time_coupling >= dt_coupling .AND. terminate_coupled == 0 )
+       DO WHILE ( time_coupling >= dt_coupling  .AND.  terminate_coupled == 0 )
           CALL surface_coupler
           time_coupling = time_coupling - dt_coupling
        ENDDO
-       IF (time_coupling == 0.0_wp .AND. time_since_reference_point < dt_coupling)&
+       IF (time_coupling == 0.0_wp  .AND.                                      &
+           time_since_reference_point < dt_coupling )                          &
        THEN
           time_coupling = time_since_reference_point
@@ -372,18 +385,4 @@
 
        CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
-
-!
-!--    Determine size of next time step
-       IF ( simulated_time /= 0.0_wp )  THEN
-          CALL timestep
-!
-!--       Synchronize the timestep in case of nested run.
-          IF ( nested_run )  THEN
-!
-!--          Synchronize by unifying the time step.
-!--          Global minimum of all time-steps is used for all.      
-             CALL pmci_synchronize
-          ENDIF
-       ENDIF
 
 !
@@ -1026,17 +1025,4 @@
 
 !
-!--    Computation and output of run control parameters.
-!--    This is also done whenever perturbations have been imposed
-       IF ( time_run_control >= dt_run_control  .OR.                     &
-            timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created ) &
-       THEN
-          CALL run_control
-          IF ( time_run_control >= dt_run_control )  THEN
-             time_run_control = MOD( time_run_control, &
-                                     MAX( dt_run_control, dt_3d ) )
-          ENDIF
-       ENDIF
-
-!
 !--    Profile output (ASCII) on file
        IF ( time_dopr_listing >= dt_dopr_listing )  THEN
@@ -1149,4 +1135,46 @@
 
 !
+!--    Determine size of next time step. Save timestep dt_3d because it is
+!--    newly calculated in routine timestep, but required further below for
+!--    steering the run control output interval
+       dt_3d_old = dt_3d
+       CALL timestep
+
+!
+!--    Computation and output of run control parameters.
+!--    This is also done whenever perturbations have been imposed
+       IF ( time_run_control >= dt_run_control  .OR.                     &
+            timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created ) &
+       THEN
+!          IF ( current_timestep_number == 1 )  THEN
+!             IF ( nxl < 7 .AND.  nxr > 7 .AND. nys < 7 .AND. nyn > 7 )  THEN
+!                u(10,7,7) = 0.55
+!             ENDIF
+!             PRINT*, 'calculating minmax'
+!             CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u,       &
+!                                  'abs', 0.0_wp, u_max, u_max_ijk )
+!             CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v,       &
+!                                  'abs', 0.0_wp, v_max, v_max_ijk )
+!             CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w,       &
+!                                  'abs', 0.0_wp, w_max, w_max_ijk )
+!             PRINT*, 'calculated u_max = ', u_max, '  myid = ', myid
+!          ENDIF
+          CALL run_control
+          IF ( time_run_control >= dt_run_control )  THEN
+             time_run_control = MOD( time_run_control, &
+                                     MAX( dt_run_control, dt_3d_old ) )
+          ENDIF
+       ENDIF
+
+!
+!--    Synchronize the timestep in case of nested run.
+       IF ( nested_run )  THEN
+!
+!--       Synchronize by unifying the time step.
+!--       Global minimum of all time-steps is used for all.
+          CALL pmci_synchronize
+       ENDIF
+
+!
 !--    Output elapsed simulated time in form of a progress bar on stdout
        IF ( myid == 0 )  CALL output_progress_bar
