Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 2293)
+++ palm/trunk/SOURCE/Makefile	(revision 2296)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Added time_integration_spinup
+# 
+# 2269 2017-06-09 11:57:32Z suehring
 # Add dependency in read_3d_binary
 # 
@@ -377,5 +380,5 @@
 	surface_coupler.f90 surface_layer_fluxes_mod.f90 swap_timelevel.f90 \
 	synthetic_turbulence_generator_mod.f90 temperton_fft_mod.f90 \
-	time_integration.f90 time_to_string.f90 timestep.f90 \
+	time_integration.f90 time_integration_spinup.f90 time_to_string.f90 timestep.f90 \
 	timestep_scheme_steering.f90 transpose.f90 tridia_solver_mod.f90 \
 	urban_surface_mod.f90 \
@@ -619,4 +622,5 @@
         surface_layer_fluxes_mod.o surface_mod.o \
 	urban_surface_mod.o virtual_flight_mod.o wind_turbine_model_mod.o
+time_integration_spinup.o: cpulog_mod.o modules.o land_surface_model_mod.o mod_kinds.o radiation_model_mod.o surface_layer_fluxes_mod.o urban_surface_mod.o 
 time_to_string.o: mod_kinds.o
 timestep.o: modules.o cpulog_mod.o mod_kinds.o microphysics_mod.o \
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 2293)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 2296)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Enabled output of radiation quantities for radiation_scheme = 'constant'
+! 
+! 2292 2017-06-20 09:51:42Z schwenkel
 ! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
 ! includes two more prognostic equations for cloud drop concentration (nc)  
@@ -1087,5 +1090,5 @@
              ENDIF
 
-             IF ( radiation  .AND.  radiation_scheme /= 'constant' )  THEN
+             IF ( radiation )  THEN
                 sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  + rad_net(j,i)
                 sums_l(nzb,100,tn)  = sums_l(nzb,100,tn)  + rad_lw_in(nzb,j,i)
Index: palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2293)
+++ palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2296)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Bugfix in calculation of bare soil heat capacity.
+! Bugfix in calculation of shf
+! Added support for spinups
+! 
+! 2282 2017-06-13 11:38:46Z schwenkel
 ! Bugfix for check of saturation moisture
 ! 
@@ -226,8 +231,10 @@
 
     USE control_parameters,                                                    &
-        ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,    &
+        ONLY:  cloud_physics, coupling_start_time, dt_3d, end_time, humidity,  &
+               intermediate_timestep_count,                                    &
                initializing_actions, intermediate_timestep_count_max,          &
                land_surface, max_masks, precipitation, pt_surface,             &
-               rho_surface, surface_pressure, timestep_scheme, tsc,            &
+               rho_surface, spinup, spinup_pt_mean, spinup_time,               &
+               surface_pressure, timestep_scheme, tsc,                         &
                time_since_reference_point
 
@@ -1230,5 +1237,5 @@
        ENDIF
     ENDDO
-  
+ 
  END SUBROUTINE lsm_check_parameters
  
@@ -1374,5 +1381,5 @@
           c_surface_tmp = (rho_c_soil * ( 1.0_wp - surf%m_sat(nzb_soil,m) )    &
                            + rho_c_water * surf_m_soil%var_2d(nzb_soil,m) )    &
-                           * dz_soil_layer(nzb_soil) * 0.25_wp               
+                           * dz_soil_layer(nzb_soil) * 0.5_wp               
 
        ELSEIF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m)) &
@@ -2069,8 +2076,5 @@
        dz_soil_layer(nzb_soil) = 2.0_wp * zs(nzb_soil)
        zs_layer(nzb_soil)      = 2.0_wp * zs(nzb_soil)
-
-
-       
-       
+  
        DO  k = nzb_soil+1, nzt_soil
           dz_soil_layer(k) = ( zs(k) - zs_layer(k-1) ) * 2.0_wp
@@ -2322,6 +2326,5 @@
                 DO  k = nzb_soil, nzt_soil  
                    t_soil_h%var_2d(k,m)   = soil_temperature(k)
-                   m_soil_h%var_2d(k,m)   = MAX(soil_moisture(k),surf_lsm_h%m_wilt(k,m))
-                   soil_moisture(k)       = MAX(soil_moisture(k),wilting_point)           
+                   m_soil_h%var_2d(k,m)   = soil_moisture(k)
                 ENDDO
                 t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
@@ -2343,6 +2346,5 @@
                 DO  k = nzb_soil, nzt_soil                        
                    t_soil_h%var_2d(k,m)   = soil_temperature(k)
-                   m_soil_h%var_2d(k,m)   = MAX(soil_moisture(k),surf_lsm_h%m_wilt(k,m))
-                   soil_moisture(k)       = MAX(soil_moisture(k),wilting_point)         
+                   m_soil_h%var_2d(k,m)   = soil_moisture(k)       
                 ENDDO 
                 t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
@@ -2411,7 +2413,6 @@
                    surf_lsm_v(l)%z0q(m)         = z0h_pavement  
                    DO  k = nzb_soil, nzt_soil  
-                      t_soil_v(l)%var_2d(k,m)      = soil_temperature(k)
-                      m_soil_v(l)%var_2d(k,m)   = MAX(soil_moisture(k),surf_lsm_v(l)%m_wilt(k,m))
-                      soil_moisture(k)          = MAX(soil_moisture(k),wilting_point)   
+                      t_soil_v(l)%var_2d(k,m)   = soil_temperature(k)
+                      m_soil_v(l)%var_2d(k,m)   = soil_moisture(k) 
                    ENDDO
                    t_soil_v(l)%var_2d(nzt_soil+1,m)  = soil_temperature(nzt_soil+1)
@@ -2433,6 +2434,5 @@
                    DO  k = nzb_soil, nzt_soil                        
                       t_soil_v(l)%var_2d(k,m)   = soil_temperature(k)
-                      m_soil_v(l)%var_2d(k,m)   = MAX(soil_moisture(k),surf_lsm_v(l)%m_wilt(k,m))
-                      soil_moisture(k)          = MAX(soil_moisture(k),wilting_point)         
+                      m_soil_v(l)%var_2d(k,m)   = soil_moisture(k)       
                    ENDDO 
                    t_soil_v(l)%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
@@ -2857,4 +2857,16 @@
        land_surface = .TRUE.
 
+!
+!--    Activate spinup
+       IF ( spinup_time > 0.0_wp )  THEN
+          coupling_start_time = spinup_time
+          end_time = end_time + spinup_time
+          IF ( spinup_pt_mean == 9999999.9_wp )  THEN
+             spinup_pt_mean = pt_surface
+          ENDIF
+          spinup = .TRUE.
+       ENDIF
+
+
  10    CONTINUE
        
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 2293)
+++ palm/trunk/SOURCE/modules.f90	(revision 2296)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Added parameters for model spinup
+! 
+! 2292 2017-06-20 09:51:42Z schwenkel
 ! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
 ! includes two more prognostic equations for cloud drop concentration (nc)  
@@ -1144,4 +1147,5 @@
     LOGICAL ::  constant_waterflux = .TRUE.                  !< water flux at all surfaces constant?
     LOGICAL ::  create_disturbances = .TRUE.                 !< namelist parameter
+    LOGICAL ::  data_output_during_spinup = .FALSE.          !< namelist parameter
     LOGICAL ::  data_output_2d_on_each_pe = .TRUE.           !< namelist parameter
     LOGICAL ::  disturbance_created = .FALSE.                !< flow disturbance imposed?
@@ -1197,4 +1201,5 @@
     LOGICAL ::  sloping_surface = .FALSE.                    !< use sloped surface?
                                                              !< (namelist parameter alpha_surface)
+    LOGICAL ::  spinup = .FALSE.                             !< perform model spinup without atmosphere code?
     LOGICAL ::  stop_dt = .FALSE.                            !< (Siggi add short description)
     LOGICAL ::  synchronous_exchange = .FALSE.               !< namelist parameter
@@ -1273,4 +1278,5 @@
     REAL(wp) ::  dt_restart = 9999999.9_wp                     !< namelist parameter
     REAL(wp) ::  dt_run_control = 60.0_wp                      !< namelist parameter
+    REAL(wp) ::  dt_spinup = 60.0_wp                           !< namelist parameter
     REAL(wp) ::  dt_3d = 1.0_wp                                !< time step
     REAL(wp) ::  dz = -1.0_wp                                  !< namelist parameter
@@ -1330,4 +1336,7 @@
     REAL(wp) ::  skip_time_do2d_yz = 9999999.9_wp              !< namelist parameter
     REAL(wp) ::  skip_time_do3d = 9999999.9_wp                 !< namelist parameter
+    REAL(wp) ::  spinup_pt_amplitude = 9999999.9_wp            !< namelist parameter
+    REAL(wp) ::  spinup_pt_mean = 9999999.9_wp                 !< namelist parameter
+    REAL(wp) ::  spinup_time = 0.0_wp                          !< namelist parameter
     REAL(wp) ::  surface_heatflux = 9999999.9_wp               !< namelist parameter
     REAL(wp) ::  surface_pressure = 1013.25_wp                 !< namelist parameter
Index: palm/trunk/SOURCE/palm.f90
===================================================================
--- palm/trunk/SOURCE/palm.f90	(revision 2293)
+++ palm/trunk/SOURCE/palm.f90	(revision 2296)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Added call to new spinup routine
+! 
+! 2292 2017-06-20 09:51:42Z schwenkel
 ! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
 ! includes two more prognostic equations for cloud drop concentration (nc)  
@@ -184,5 +187,5 @@
                message_string, microphysics_morrison, microphysics_seifert,    &
                nest_domain, neutral, nudging, passive_scalar, runnr,           &
-               simulated_time, simulated_time_chr, urban_surface,              &
+               simulated_time, simulated_time_chr, spinup, urban_surface,      &
                user_interface_current_revision,                                &
                user_interface_required_revision, version, wall_heatflux,       &
@@ -392,4 +395,11 @@
 
     CALL cpu_log( log_point(2), 'initialisation', 'stop' )
+
+!
+!-- Integration of the non-atmospheric equations (land surface model, urban
+!-- surface model)
+    IF ( spinup )  THEN
+       CALL time_integration_spinup
+    ENDIF
 
 !
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 2293)
+++ palm/trunk/SOURCE/parin.f90	(revision 2296)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Added parameters for model spinup
+! 
+! 2292 2017-06-20 09:51:42Z schwenkel
 ! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
 ! includes two more prognostic equations for cloud drop concentration (nc)  
@@ -357,7 +360,8 @@
              coupling_start_time, curvature_solution_effects_bulk,             &
              cycle_mg, damp_level_1d,                                          &
+             data_output_during_spinup,                                        &
              dissipation_1d,                                                   &
              dp_external, dp_level_b, dp_smooth, dpdxy,                        &
-             dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, dz_max,              &
+             dt, dt_pr_1d, dt_run_control_1d, dt_spinup, dx, dy, dz, dz_max,   &
              dz_stretch_factor, dz_stretch_level, end_time_1d,                 &
              ensemble_member_nr, e_init, e_min, fft_method,                    &
@@ -386,4 +390,5 @@
              sa_vertical_gradient, sa_vertical_gradient_level, scalar_advec,   &
              scalar_rayleigh_damping,                                          &
+             spinup_time, spinup_pt_amplitude, spinup_pt_mean,                 &
              statistic_regions, subs_vertical_gradient,                        &
              subs_vertical_gradient_level, surface_heatflux, surface_pressure, &
@@ -547,5 +552,4 @@
              initializing_actions = 'set_constant_profiles'
           ENDIF
-
             
 !
Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2293)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2296)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Added output of rad_sw_out for radiation_scheme = 'constant'
+! 
+! 2270 2017-06-09 12:18:47Z maronga
 ! Numbering changed (2 timeseries removed)
 ! 
@@ -1254,4 +1257,6 @@
                                   + rad_lw_out(0,j,i) )                        &
                                   / ( 1.0_wp - alpha(j,i) )
+
+             rad_sw_out(0,j,i) =  alpha(j,i) * rad_sw_in(0,j,i)
 
           ENDDO
Index: palm/trunk/SOURCE/time_integration_spinup.f90
===================================================================
--- palm/trunk/SOURCE/time_integration_spinup.f90	(revision 2296)
+++ palm/trunk/SOURCE/time_integration_spinup.f90	(revision 2296)
@@ -0,0 +1,366 @@
+!> @file time_integration_spinup.f90
+!------------------------------------------------------------------------------!
+! This file is part of PALM.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2017 Leibniz Universitaet Hannover
+!------------------------------------------------------------------------------!
+!
+! Current revisions:
+! ------------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
+! Initial revision
+! 
+! 
+!
+!
+! Description:
+! ------------
+!> Integration in time of the non-atmospheric model components such as land
+!> surface model and urban surface model
+!------------------------------------------------------------------------------!
+ SUBROUTINE time_integration_spinup
+ 
+    USE arrays_3d,                                                             &
+        ONLY:  pt, pt_p
+
+    USE control_parameters,                                                    &
+        ONLY:  averaging_interval_pr, constant_diffusion, constant_flux_layer, coupling_start_time, current_timestep_number, data_output_during_spinup, disturbance_created, dopr_n, do_sum, dt_averaging_input_pr, dt_dopr, dt_dots, dt_run_control, dt_spinup, humidity, intermediate_timestep_count, intermediate_timestep_count_max, land_surface, nr_timesteps_this_run, simulated_time, simulated_time_chr, skip_time_dopr, spinup, spinup_pt_amplitude, spinup_pt_mean, spinup_time, timestep_count, timestep_scheme, time_dopr, time_dopr_av, time_dots, time_run_control, time_since_reference_point, urban_surface
+
+    USE constants,                                                             &
+        ONLY:  pi
+
+    USE cpulog,                                                                &
+        ONLY:  cpu_log, log_point, log_point_s
+
+    USE indices,                                                               &
+        ONLY:  nbgp, nzb, nzt, nysg, nyng, nxlg, nxrg
+
+
+    USE land_surface_model_mod,                                                &
+        ONLY:  lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel,         &
+               skip_time_do_lsm
+
+    USE pegrid
+
+    USE kinds
+
+    USE radiation_model_mod,                                                   &
+        ONLY:  dt_radiation, force_radiation_call, radiation,                  &
+               radiation_control, skip_time_do_radiation, time_radiation
+
+    USE statistics,                                                            &
+        ONLY:  flow_statistics_called
+
+    USE surface_layer_fluxes_mod,                                              &
+        ONLY:  surface_layer_fluxes
+
+    USE surface_mod,                                                        &
+        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
+                surf_usm_v
+
+    USE urban_surface_mod,                                                     &
+        ONLY:  usm_material_heat_model, usm_material_model,                    &
+               usm_radiation, usm_surface_energy_balance, usm_swap_timelevel
+
+
+
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=9) ::  time_to_string          !< 
+  
+    INTEGER(iwp) :: i
+    INTEGER(iwp) :: j
+    INTEGER(iwp) :: k
+    INTEGER(iwp) :: l
+    INTEGER(iwp) :: m
+  
+    REAL(wp) ::  pt_spinup   !< temporary storage of temperature
+                  
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pt_save   !< temporary storage of temperature
+
+    REAL(wp) ::  day_length = 43200.0_wp !! must be calculated from time_utc_init, day_init, and latitude/longitude
+
+    ALLOCATE( pt_save(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+    CALL exchange_horiz( pt_p, nbgp )   
+    pt_save = pt_p
+
+    CALL location_message( 'starting spinup-sequence', .TRUE. )
+!
+!-- Start of the time loop
+    DO  WHILE ( simulated_time < spinup_time )
+
+       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'start' )
+   
+!
+!--    Start of intermediate step loop
+       intermediate_timestep_count = 0
+       DO  WHILE ( intermediate_timestep_count < &
+                   intermediate_timestep_count_max )
+
+          intermediate_timestep_count = intermediate_timestep_count + 1
+
+!
+!--       Set the steering factors for the prognostic equations which depend
+!--       on the timestep scheme
+          CALL timestep_scheme_steering
+
+
+!!!!      Set new values of pt_p here
+          pt_spinup = spinup_pt_mean + spinup_pt_amplitude * SIN( pi* (time_since_reference_point/day_length) - pi)
+
+          IF ( land_surface )  THEN
+             DO  m = 1, surf_lsm_h%ns
+                i   = surf_lsm_h%i(m)            
+                j   = surf_lsm_h%j(m)
+                k   = surf_lsm_h%k(m)
+                pt_p(k,j,i) = pt_spinup
+             ENDDO
+
+             DO  l = 0, 3
+                DO  m = 1, surf_lsm_v(l)%ns
+                   i   = surf_lsm_v(l)%i(m)            
+                   j   = surf_lsm_v(l)%j(m)
+                   k   = surf_lsm_v(l)%k(m)
+                   pt_p(k,j,i) = pt_spinup
+                ENDDO
+             ENDDO
+          ENDIF
+
+          IF ( urban_surface )  THEN
+             DO  m = 1, surf_usm_h%ns
+                i   = surf_usm_h%i(m)            
+                j   = surf_usm_h%j(m)
+                k   = surf_usm_h%k(m)
+                pt_p(k,j,i) = pt_spinup
+             ENDDO
+
+             DO  l = 0, 3
+                DO  m = 1, surf_usm_v(l)%ns
+                   i   = surf_usm_v(l)%i(m)            
+                   j   = surf_usm_v(l)%j(m)
+                   k   = surf_usm_v(l)%k(m)
+                   pt_p(k,j,i) = pt_spinup
+                ENDDO
+             ENDDO
+          ENDIF
+
+!
+!--       Swap the time levels in preparation for the next time step.
+          timestep_count = timestep_count + 1
+     
+          IF ( land_surface )  THEN
+              CALL lsm_swap_timelevel ( 0 )
+          ENDIF
+
+          IF ( urban_surface )  THEN
+             CALL usm_swap_timelevel ( 0 )
+          ENDIF
+
+          IF ( land_surface )  THEN
+             CALL lsm_swap_timelevel ( MOD( timestep_count, 2) )
+          ENDIF
+
+          IF ( urban_surface )  THEN
+             CALL usm_swap_timelevel ( MOD( timestep_count, 2) )
+          ENDIF
+         
+!
+!--       If required, compute virtual potential temperature 
+          IF ( humidity )  THEN 
+             CALL compute_vpt 
+          ENDIF 
+
+!
+!--       Compute the diffusion quantities
+          IF ( .NOT. constant_diffusion )  THEN
+
+!
+!--          First the vertical (and horizontal) fluxes in the surface 
+!--          (constant flux) layer are computed
+             IF ( constant_flux_layer )  THEN
+                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'start' )
+                CALL surface_layer_fluxes
+                CALL cpu_log( log_point(19), 'surface_layer_fluxes', 'stop' )
+             ENDIF
+
+!
+!--          If required, solve the energy balance for the surface and run soil 
+!--          model. Call for horizontal as well as vertical surfaces
+             IF ( land_surface .AND. simulated_time > skip_time_do_lsm)  THEN
+
+                CALL cpu_log( log_point(54), 'land_surface', 'start' )
+!
+!--             Call for horizontal upward-facing surfaces
+                CALL lsm_energy_balance( .TRUE., -1 )
+                CALL lsm_soil_model( .TRUE., -1 )
+!
+!--             Call for northward-facing surfaces
+                CALL lsm_energy_balance( .FALSE., 0 )
+                CALL lsm_soil_model( .FALSE., 0 )
+!
+!--             Call for southward-facing surfaces
+                CALL lsm_energy_balance( .FALSE., 1 )
+                CALL lsm_soil_model( .FALSE., 1 )
+!
+!--             Call for eastward-facing surfaces
+                CALL lsm_energy_balance( .FALSE., 2 )
+                CALL lsm_soil_model( .FALSE., 2 )
+!
+!--             Call for westward-facing surfaces
+                CALL lsm_energy_balance( .FALSE., 3 )
+                CALL lsm_soil_model( .FALSE., 3 )
+
+                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
+             ENDIF
+
+!
+!--          If required, solve the energy balance for urban surfaces and run 
+!--          the material heat model
+             IF (urban_surface) THEN
+                CALL cpu_log( log_point(74), 'urban_surface', 'start' )
+                CALL usm_surface_energy_balance
+                IF ( usm_material_model )  THEN
+                   CALL usm_material_heat_model
+                ENDIF
+                CALL cpu_log( log_point(74), 'urban_surface', 'stop' )
+             ENDIF
+
+          ENDIF
+
+!
+!--       If required, calculate radiative fluxes and heating rates
+          IF ( radiation .AND. intermediate_timestep_count                     &
+               == intermediate_timestep_count_max .AND. simulated_time >    &
+               skip_time_do_radiation )  THEN
+
+               time_radiation = time_radiation + dt_spinup
+
+             IF ( time_radiation >= dt_radiation .OR. force_radiation_call )   &
+             THEN
+
+                CALL cpu_log( log_point(50), 'radiation', 'start' )
+
+                IF ( .NOT. force_radiation_call )  THEN
+                   time_radiation = time_radiation - dt_radiation
+                ENDIF
+
+                CALL radiation_control
+
+                CALL cpu_log( log_point(50), 'radiation', 'stop' )
+
+                IF (urban_surface)  THEN
+                   CALL cpu_log( log_point(75), 'usm_radiation', 'start' )
+                   CALL usm_radiation
+                   CALL cpu_log( log_point(75), 'usm_radiation', 'stop' )
+                ENDIF
+             ENDIF
+          ENDIF
+
+       ENDDO   ! Intermediate step loop
+
+!
+!--    Increase simulation time and output times
+       nr_timesteps_this_run      = nr_timesteps_this_run + 1
+       current_timestep_number    = current_timestep_number + 1
+       simulated_time             = simulated_time   + dt_spinup
+       simulated_time_chr         = time_to_string( simulated_time )
+       time_since_reference_point = simulated_time - coupling_start_time
+
+       IF ( data_output_during_spinup )  THEN
+          time_dots          = time_dots        + dt_spinup
+          IF ( simulated_time >= skip_time_dopr )  THEN
+             time_dopr       = time_dopr        + dt_spinup
+          ENDIF
+          time_run_control   = time_run_control + dt_spinup
+
+!
+!--       Carry out statistical analysis and output at the requested output times.
+!--       The MOD function is used for calculating the output time counters (like
+!--       time_dopr) in order to regard a possible decrease of the output time
+!--       interval in case of restart runs
+
+!
+!--       Set a flag indicating that so far no statistics have been created
+!--       for this time step
+          flow_statistics_called = .FALSE.
+
+!
+!--       If required, call flow_statistics for averaging in time
+          IF ( averaging_interval_pr /= 0.0_wp  .AND.                          &
+             ( dt_dopr - time_dopr ) <= averaging_interval_pr  .AND.           &
+             simulated_time >= skip_time_dopr )  THEN
+             time_dopr_av = time_dopr_av + dt_spinup
+             IF ( time_dopr_av >= dt_averaging_input_pr )  THEN
+                do_sum = .TRUE.
+                time_dopr_av = MOD( time_dopr_av,                              &
+                               MAX( dt_averaging_input_pr, dt_spinup ) )
+             ENDIF
+          ENDIF
+          IF ( do_sum )  CALL flow_statistics
+
+!
+!--       Output of profiles
+          IF ( time_dopr >= dt_dopr )  THEN
+             IF ( dopr_n /= 0 )  CALL data_output_profiles
+             time_dopr = MOD( time_dopr, MAX( dt_dopr, dt_spinup ) )
+             time_dopr_av = 0.0_wp    ! due to averaging (see above)
+          ENDIF
+
+!
+!--       Output of time series
+          IF ( time_dots >= dt_dots )  THEN
+             CALL data_output_tseries
+             time_dots = MOD( time_dots, MAX( dt_dots, dt_spinup ) )
+          ENDIF
+
+       ENDIF
+
+!
+!--    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_spinup ) )
+          ENDIF
+       ENDIF
+
+       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
+
+       IF ( myid == 0 )  THEN
+          PRINT*, time_since_reference_point
+       ENDIF
+
+    ENDDO   ! time loop
+
+!
+!-- Write back saved temperature to the 3D arrays
+    pt(:,:,:)   = pt_save
+    pt_p(:,:,:) = pt_save
+
+    DEALLOCATE(pt_save)
+
+    CALL location_message( 'finished time-stepping spinup', .TRUE. )
+
+ END SUBROUTINE time_integration_spinup
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 2293)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 2296)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! Bugfix: adjust output of shf. Added support for spinups
+! 
+! 2287 2017-06-15 16:46:30Z suehring
 ! Bugfix in determination topography-top index
 ! 
@@ -139,12 +142,12 @@
     
     USE control_parameters,                                                    &
-        ONLY:  dz, topography, dt_3d, intermediate_timestep_count,             &
-               initializing_actions, intermediate_timestep_count_max,          &
-               simulated_time, end_time, timestep_scheme, tsc,                 &
-               coupling_char, io_blocks, io_group, message_string,             &
-               time_since_reference_point, surface_pressure,                   &
-               g, pt_surface, large_scale_forcing, lsf_surf,                   &
-               time_do3d, dt_do3d, average_count_3d, varnamelength,            &
-               urban_surface
+        ONLY:  coupling_start_time, dz, topography, dt_3d,                     &
+               intermediate_timestep_count, initializing_actions,              &
+               intermediate_timestep_count_max, simulated_time, end_time,      &
+               timestep_scheme, tsc, coupling_char, io_blocks, io_group,       &
+               message_string, time_since_reference_point, surface_pressure,   &
+               g, pt_surface, large_scale_forcing, lsf_surf, spinup,           &
+               spinup_pt_mean, spinup_time, time_do3d, dt_do3d,                &
+               average_count_3d, varnamelength, urban_surface
 
     USE cpulog,                                                                &
@@ -3489,4 +3492,15 @@
 !--    Set flag that indicates that the land surface model is switched on
        urban_surface = .TRUE.
+
+!
+!--    Activate spinup
+       IF ( spinup_time > 0.0_wp )  THEN
+          coupling_start_time = spinup_time
+          end_time = end_time + spinup_time
+          IF ( spinup_pt_mean == 9999999.9_wp )  THEN
+             spinup_pt_mean = pt_surface
+          ENDIF
+          spinup = .TRUE.
+       ENDIF
 
  10    CONTINUE
@@ -5161,5 +5175,5 @@
 !--        store kinematic surface heat fluxes for utilization in other processes
 !--        diffusion_s, surface_layer_fluxes,...
-           surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / rho_cp
+           surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / cp
 
        ENDDO
@@ -5282,5 +5296,5 @@
 !--           store kinematic surface heat fluxes for utilization in other processes
 !--           diffusion_s, surface_layer_fluxes,...
-              surf_usm_v(l)%shf(m) = surf_usm_v(l)%wshf_eb(m) / rho_cp
+              surf_usm_v(l)%shf(m) = surf_usm_v(l)%wshf_eb(m) / cp
 
            ENDDO
@@ -5316,5 +5330,5 @@
 !--              Instead of this, we can adjust shf in case AH only at surface
                  surf_usm_h%shf(m) = surf_usm_h%shf(m) +                       &
-                                   aheat(j,i) * acoef * ddx * ddy / rho_cp
+                                   aheat(j,i) * acoef * ddx * ddy / cp
               ENDIF
            ENDDO
