Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 2298)
+++ /palm/trunk/SOURCE/header.f90	(revision 2299)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Modified output for spinups
+! 
+! 2298 2017-06-29 09:28:18Z raasch
 ! MPI2 related parts removed
 ! 
@@ -507,5 +510,5 @@
     ENDIF
 #if defined( __parallel )
-    IF ( coupling_start_time /= 0.0_wp )  THEN
+    IF ( coupling_start_time /= 0.0_wp  .AND. .NOT. spinup )  THEN
        IF ( coupling_start_time > simulated_time_at_begin )  THEN
           WRITE ( io, 109 )
Index: /palm/trunk/SOURCE/init_1d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_1d_model.f90	(revision 2298)
+++ /palm/trunk/SOURCE/init_1d_model.f90	(revision 2299)
@@ -21,8 +21,11 @@
 ! -----------------
 ! 
-!
+! 
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Removed german text
+!
+! 2101 2017-01-05 16:42:31Z suehring
 !
 ! 2059 2016-11-10 14:20:40Z maronga
@@ -851,5 +854,5 @@
 !
 !-- formats
-100 FORMAT (///'1D-Zeitschrittkontrollausgaben:'/ &
+100 FORMAT (///'1D run control output:'/ &
               &'------------------------------'// &
            &'ITER.  HH:MM:SS    DT      UMAX   VMAX    U*   ALPHA   ENERG.'/ &
Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2298)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2299)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Removed pt_p from USE statement. Adjusted call to lsm_soil_model to allow
+! spinups without soil moisture prediction
+! 
+! 2298 2017-06-29 09:28:18Z raasch
 ! type of write_binary changed from CHARACTER to LOGICAL
 ! 
@@ -228,5 +232,5 @@
  
     USE arrays_3d,                                                             &
-        ONLY:  hyp, pt, pt_p, prr, q, q_p, ql, vpt, u, v, w
+        ONLY:  hyp, pt, prr, q, q_p, ql, vpt, u, v, w
 
     USE cloud_parameters,                                                      &
@@ -2884,5 +2888,5 @@
 !> temperature and water content.
 !------------------------------------------------------------------------------!
-    SUBROUTINE lsm_soil_model( horizontal, l )
+    SUBROUTINE lsm_soil_model( horizontal, l, calc_soil_moisture )
 
 
@@ -2892,4 +2896,6 @@
        INTEGER(iwp) ::  l       !< surface-data type index indication facing
        INTEGER(iwp) ::  m       !< running index
+
+       LOGICAL, INTENT(IN) ::  calc_soil_moisture !< flag indicating whether soil moisture shall be calculated or not.
 
        LOGICAL      ::  horizontal !< flag indication horizontal wall, required to set pointer accordingly
@@ -3059,4 +3065,14 @@
 
              ENDDO
+
+          ENDIF
+
+       ENDDO
+
+
+       DO  m = 1, surf%ns
+
+          IF (  .NOT.  surf%water_surface(m)  .AND.  calc_soil_moisture )  THEN
+
 
 !
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2298)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2299)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Improved syntax layout
+! 
+! 2298 2017-06-29 09:28:18Z raasch
 ! type of write_binary changed from CHARACTER to LOGICAL
 ! 
@@ -505,6 +508,6 @@
 !
 !-- Public variables and constants / NEEDS SORTING
-    PUBLIC dots_rad, dt_radiation, force_radiation_call,                       &
-           rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
+    PUBLIC decl_1, decl_2, decl_3, dots_rad, dt_radiation, force_radiation_call,&
+           lat, lon, rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
            rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0,       &
            rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
@@ -632,5 +635,6 @@
 !> Check data output of profiles for radiation model
 !------------------------------------------------------------------------------!  
-    SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, dopr_unit )
+    SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit,      &
+               dopr_unit )
  
        USE arrays_3d,                                                          &
@@ -1237,5 +1241,6 @@
           DO j = nysg, nyng
 !
-!--          Obtain vertical index of topography top. So far it is identical to nzb.
+!--          Obtain vertical index of topography top. So far it is identical to
+!--          nzb.
              k = MAXLOC(                                                       &
                         MERGE( 1, 0,                                           &
@@ -1314,11 +1319,11 @@
 
  1 FORMAT ('    Geograph. longitude            :   lambda = ',F4.1,' degr')
- 2 FORMAT ('    Day of the year at model start :   day_init = ',I3   &
+ 2 FORMAT ('    Day of the year at model start :   day_init = ',I3             &
             /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
  3 FORMAT (//' Radiation model information:'/                                  &
               ' ----------------------------'/)
- 4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,        &
+ 4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,     &
            // 'W/m**2')
- 5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',   &
+ 5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',&
                    ' default)')
  6 FORMAT ('    --> RRTMG scheme is used')
@@ -1493,8 +1498,8 @@
 !--                   Calculate cloud droplet effective radius
                       IF ( cloud_physics )  THEN
-                         rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)        &
-                                           * rho_surface                          &
-                                           / ( 4.0_wp * pi * nc_const * rho_l )   &
-                                           )**0.33333333333333_wp                 &
+                         rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)     &
+                                           * rho_surface                       &
+                                           / ( 4.0_wp * pi * nc_const * rho_l )&
+                                           )**0.33333333333333_wp              &
                                            * EXP( LOG( sigma_gc )**2 )
 
@@ -1636,5 +1641,5 @@
 !--    Calculate current day and time based on the initial values and simulation
 !--    time
-       day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point)    &
+       day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point)&
                                / 86400.0_wp ), KIND=iwp)
        time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
@@ -1648,5 +1653,5 @@
 !
 !--    Calculate cosine of solar zenith angle
-       zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)      &
+       zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)   &
                                             * COS(hour_angle)
        zenith(0) = MAX(0.0_wp,zenith(0))
@@ -1655,6 +1660,10 @@
 !--    Calculate solar directional vector
        IF ( sun_direction )  THEN
+
+!
 !--       Direction in longitudes equals to sin(solar_azimuth) * sin(zenith)
           sun_dir_lon(0) = -SIN(hour_angle) * COS(declination)
+
+!
 !--       Direction in latitues equals to cos(solar_azimuth) * sin(zenith)
           sun_dir_lat(0) = SIN(declination) * COS(lat) - COS(hour_angle) &
@@ -1942,5 +1951,5 @@
           rrtm_h2ovmr(0,k) = q_snd(k)
        ENDDO
-       rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                 &
+       rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                  &
                                 - rrtm_tlay(0,nzt_rad-1)
        DO k = nzt+9, nzt_rad+1
@@ -2058,8 +2067,8 @@
 
 
-       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  p_mls,         & !< pressure levels for the absorbers
-                                                 rrtm_play_tmp, & !< temporary array for pressure zu-levels
-                                                 rrtm_plev_tmp, & !< temporary array for pressure zw-levels
-                                                 trace_path_tmp   !< temporary array for storing trace gas path data
+       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  p_mls,          & !< pressure levels for the absorbers
+                                                 rrtm_play_tmp,  & !< temporary array for pressure zu-levels
+                                                 rrtm_plev_tmp,  & !< temporary array for pressure zw-levels
+                                                 trace_path_tmp    !< temporary array for storing trace gas path data
 
        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  trace_mls,      & !< array for storing the absorber amounts
@@ -2497,5 +2506,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) + rad_lw_out(k,j,i)
+                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i)              &
+                                             + rad_lw_out(k,j,i)
                    ENDDO
                 ENDDO
@@ -2506,5 +2516,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) + rad_lw_cs_hr(k,j,i)
+                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i)          &
+                                               + rad_lw_cs_hr(k,j,i)
                    ENDDO
                 ENDDO
@@ -2515,5 +2526,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) + rad_lw_hr(k,j,i)
+                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i)                &
+                                            + rad_lw_hr(k,j,i)
                    ENDDO
                 ENDDO
@@ -2524,5 +2536,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) + rad_sw_in(k,j,i)
+                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i)                &
+                                            + rad_sw_in(k,j,i)
                    ENDDO
                 ENDDO
@@ -2533,5 +2546,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) + rad_sw_out(k,j,i)
+                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i)              &
+                                             + rad_sw_out(k,j,i)
                    ENDDO
                 ENDDO
@@ -2542,5 +2556,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) + rad_sw_cs_hr(k,j,i)
+                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i)          &
+                                               + rad_sw_cs_hr(k,j,i)
                    ENDDO
                 ENDDO
@@ -2551,5 +2566,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) + rad_sw_hr(k,j,i)
+                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i)                &
+                                            + rad_sw_hr(k,j,i)
                    ENDDO
                 ENDDO
@@ -2568,5 +2584,6 @@
              DO  i = nxlg, nxrg
                 DO  j = nysg, nyng
-                   rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
+                   rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, &
+                                     KIND=wp )
                 ENDDO
              ENDDO
@@ -2576,5 +2593,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i)                &
+                                            / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -2585,5 +2603,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i)              &
+                                             / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -2594,5 +2613,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i)          &
+                                             / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -2603,5 +2623,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i)                &
+                                            / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -2612,5 +2633,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i)                &
+                                            / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -2621,5 +2643,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i)              &
+                                             / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -2630,5 +2653,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i)          &
+                                             / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -2639,5 +2663,6 @@
                 DO  j = nysg, nyng
                    DO  k = nzb, nzt+1
-                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
+                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i)                &
+                                            / REAL( average_count_3d, KIND=wp )
                    ENDDO
                 ENDDO
@@ -3412,8 +3437,8 @@
 
 
-SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file,   &
-                                     nynfa, nyn_on_file, nysfa, nys_on_file,   &
-                                     offset_xa, offset_ya, overlap_count,      &
-                                     tmp_2d, tmp_3d )
+SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa,          &
+                                        nxr_on_file, nynfa, nyn_on_file, nysfa,&
+                                        nys_on_file, offset_xa, offset_ya,     &
+                                        overlap_count, tmp_2d, tmp_3d )
  
 
@@ -3491,6 +3516,6 @@
                    ENDIF  
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
-                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                   rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =         &
+                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
                 CASE ( 'rad_net_av' )
@@ -3499,6 +3524,6 @@
                    ENDIF  
                    IF ( k == 1 )  READ ( 13 )  tmp_2d
-                   rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
-                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                   rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =      &
+                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
                 CASE ( 'rad_lw_in' )
                    IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
@@ -3515,5 +3540,5 @@
                          READ ( 13 )  tmp_3d2
                          rad_lw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
-                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
+                            tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
                       ELSE
                          READ ( 13 )  tmp_3d
@@ -3602,5 +3627,5 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                   rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -3618,5 +3643,5 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                   rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -3626,5 +3651,5 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -3722,5 +3747,5 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                   rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -3738,5 +3763,5 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                   rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -3746,5 +3771,5 @@
                    ENDIF
                    IF ( k == 1 )  READ ( 13 )  tmp_3d
-                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
+                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
                            tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
 
@@ -3753,6 +3778,6 @@
                                         TRIM( field_char ), '" found in',      &
                                         '&data from prior run on PE ', myid
-                  CALL message( 'radiation_read_restart_data', 'PA0302', 1, 2, 0, 6, &
-                                 0 )
+                  CALL message( 'radiation_read_restart_data', 'PA0302', 1, 2, &
+                                0, 6, 0 )
 
             END SELECT
Index: /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 2298)
+++ /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 2299)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Adjusted for allow separate spinups of LSM and atmosphere code
+! 
+! 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)  
@@ -203,5 +206,5 @@
                most_method, neutral, passive_scalar, pt_surface, q_surface,    &
                run_coupled, surface_pressure, simulated_time, terminate_run,   &
-               urban_surface, zeta_max, zeta_min
+               time_since_reference_point, urban_surface, zeta_max, zeta_min
 
     USE grid_variables,                                                        &
@@ -1921,7 +1924,8 @@
 !
 !--       Compute the vertical kinematic heat flux
-          IF (  .NOT.  constant_heatflux  .AND.  ( simulated_time <=           &
-               skip_time_do_lsm  .OR.  .NOT.  land_surface )  .AND.            &
-               .NOT.  urban_surface  .AND.  .NOT. downward )  THEN
+          IF (  .NOT.  constant_heatflux  .AND.  ( ( time_since_reference_point&
+               <=  skip_time_do_lsm  .AND. simulated_time > 0.0_wp ) .OR.      &
+               .NOT.  land_surface )  .AND.  .NOT.  urban_surface  .AND.       &
+               .NOT. downward )  THEN
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, surf%ns 
@@ -1935,6 +1939,7 @@
 !--       Compute the vertical water flux
           IF (  .NOT.  constant_waterflux  .AND.  humidity  .AND.              &
-                ( simulated_time <= skip_time_do_lsm                           &
-               .OR.  .NOT.  land_surface )  .AND.  .NOT. downward )  THEN
+               ( ( time_since_reference_point <= skip_time_do_lsm  .AND.       &
+               simulated_time > 0.0_wp ) .OR.  .NOT.  land_surface )  .AND.    &
+               .NOT. downward )  THEN
              !$OMP PARALLEL DO PRIVATE( i, j, k )
              DO  m = 1, surf%ns 
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 2298)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 2299)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Call of soil model adjusted
+! 
+! 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)  
@@ -770,5 +773,5 @@
 !--          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
+             IF ( land_surface .AND. time_since_reference_point > skip_time_do_lsm)  THEN
 
                 CALL cpu_log( log_point(54), 'land_surface', 'start' )
@@ -776,21 +779,21 @@
 !--             Call for horizontal upward-facing surfaces
                 CALL lsm_energy_balance( .TRUE., -1 )
-                CALL lsm_soil_model( .TRUE., -1 )
+                CALL lsm_soil_model( .TRUE., -1, .TRUE. )
 !
 !--             Call for northward-facing surfaces
                 CALL lsm_energy_balance( .FALSE., 0 )
-                CALL lsm_soil_model( .FALSE., 0 )
+                CALL lsm_soil_model( .FALSE., 0, .TRUE. )
 !
 !--             Call for southward-facing surfaces
                 CALL lsm_energy_balance( .FALSE., 1 )
-                CALL lsm_soil_model( .FALSE., 1 )
+                CALL lsm_soil_model( .FALSE., 1, .TRUE. )
 !
 !--             Call for eastward-facing surfaces
                 CALL lsm_energy_balance( .FALSE., 2 )
-                CALL lsm_soil_model( .FALSE., 2 )
+                CALL lsm_soil_model( .FALSE., 2, .TRUE. )
 !
 !--             Call for westward-facing surfaces
                 CALL lsm_energy_balance( .FALSE., 3 )
-                CALL lsm_soil_model( .FALSE., 3 )
+                CALL lsm_soil_model( .FALSE., 3, .TRUE. )
 
                 CALL cpu_log( log_point(54), 'land_surface', 'stop' )
@@ -829,5 +832,5 @@
 !--       If required, calculate radiative fluxes and heating rates
           IF ( radiation .AND. intermediate_timestep_count                     &
-               == intermediate_timestep_count_max .AND. simulated_time >    &
+               == intermediate_timestep_count_max .AND. time_since_reference_point >    &
                skip_time_do_radiation )  THEN
 
Index: /palm/trunk/SOURCE/time_integration_spinup.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration_spinup.f90	(revision 2298)
+++ /palm/trunk/SOURCE/time_integration_spinup.f90	(revision 2299)
@@ -25,10 +25,15 @@
 ! -----------------
 ! $Id$
+! Call of soil model adjusted to avoid prognostic equation for soil moisture
+! during spinup.
+! Better representation of diurnal cycle of near-surface temperature.
+! Excluded prognostic equation for soil moisture during spinup.
+! Added output of run control data for spinup.
+! 
+! 2297 2017-06-28 14:35:57Z scharf
 ! bugfixes
 ! 
 ! 2296 2017-06-28 07:53:56Z maronga
 ! Initial revision
-! 
-! 
 !
 !
@@ -50,5 +55,5 @@
                dt_spinup, humidity, intermediate_timestep_count,               &
                intermediate_timestep_count_max, land_surface,                  &
-               nr_timesteps_this_run, simulated_time, simulated_time_chr,      &
+               simulated_time, simulated_time_chr,      &
                skip_time_dopr, spinup, spinup_pt_amplitude, spinup_pt_mean,    &
                spinup_time, timestep_count, timestep_scheme, time_dopr,        &
@@ -67,8 +72,8 @@
 
     USE land_surface_model_mod,                                                &
-        ONLY:  lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel,         &
-               skip_time_do_lsm
-
-    USE pegrid
+        ONLY:  lsm_energy_balance, lsm_soil_model, lsm_swap_timelevel
+
+    USE pegrid,                                                                &
+        ONLY:  myid
 
     USE kinds
@@ -76,5 +81,5 @@
     USE radiation_model_mod,                                                   &
         ONLY:  dt_radiation, force_radiation_call, radiation,                  &
-               radiation_control, skip_time_do_radiation, time_radiation
+               radiation_control, rad_sw_in, time_radiation, time_utc_init
 
     USE statistics,                                                            &
@@ -99,20 +104,22 @@
     CHARACTER (LEN=9) ::  time_to_string          !< 
   
-    INTEGER(iwp) :: i
-    INTEGER(iwp) :: j
-    INTEGER(iwp) :: k
-    INTEGER(iwp) :: l
-    INTEGER(iwp) :: m
+    INTEGER(iwp) ::  i !< running index
+    INTEGER(iwp) ::  j !< running index
+    INTEGER(iwp) ::  k !< running index
+    INTEGER(iwp) ::  l !< running index
+    INTEGER(iwp) ::  m !< running index
+
+    INTEGER(iwp) :: current_timestep_number_spinup = 0  !< number if timestep during spinup
   
+    LOGICAL :: run_control_header_spinup = .FALSE.  !< flag parameter for steering whether the header information must be output
+
     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 exchange_horiz( pt, nbgp )   
+    pt_save = pt
 
     CALL location_message( 'starting spinup-sequence', .TRUE. )
@@ -137,7 +144,15 @@
 
 
-!!!!      Set new values of pt_p here
-          pt_spinup = spinup_pt_mean + spinup_pt_amplitude * SIN( pi* (time_since_reference_point/day_length) - pi)
-
+!
+!--       Estimate a near-surface air temperature based on the position of the 
+!--       sun and user input about mean temperature and amplitude. The time is
+!--       shifted by one hour to simulate a lag between air temperature and
+!--       incoming radiation
+          pt_spinup = spinup_pt_mean + spinup_pt_amplitude                     &
+             * solar_angle (time_utc_init + time_since_reference_point - 3600.0)
+
+!
+!--       Map air temperature to all grid points in the vicinity of a surface
+!--       element
           IF ( land_surface )  THEN
              DO  m = 1, surf_lsm_h%ns
@@ -145,5 +160,5 @@
                 j   = surf_lsm_h%j(m)
                 k   = surf_lsm_h%k(m)
-                pt_p(k,j,i) = pt_spinup
+                pt(k,j,i) = pt_spinup
              ENDDO
 
@@ -153,5 +168,5 @@
                    j   = surf_lsm_v(l)%j(m)
                    k   = surf_lsm_v(l)%k(m)
-                   pt_p(k,j,i) = pt_spinup
+                   pt(k,j,i) = pt_spinup
                 ENDDO
              ENDDO
@@ -163,5 +178,5 @@
                 j   = surf_usm_h%j(m)
                 k   = surf_usm_h%k(m)
-                pt_p(k,j,i) = pt_spinup
+                pt(k,j,i) = pt_spinup
              ENDDO
 
@@ -171,5 +186,5 @@
                    j   = surf_usm_v(l)%j(m)
                    k   = surf_usm_v(l)%k(m)
-                   pt_p(k,j,i) = pt_spinup
+                   pt(k,j,i) = pt_spinup
                 ENDDO
              ENDDO
@@ -217,6 +232,7 @@
 !
 !--          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
+!--          model. Call for horizontal as well as vertical surfaces.
+!--          The prognostic equation for soil moisure is switched off
+             IF ( land_surface )  THEN
 
                 CALL cpu_log( log_point(54), 'land_surface', 'start' )
@@ -224,21 +240,21 @@
 !--             Call for horizontal upward-facing surfaces
                 CALL lsm_energy_balance( .TRUE., -1 )
-                CALL lsm_soil_model( .TRUE., -1 )
+                CALL lsm_soil_model( .TRUE., -1, .FALSE. )
 !
 !--             Call for northward-facing surfaces
                 CALL lsm_energy_balance( .FALSE., 0 )
-                CALL lsm_soil_model( .FALSE., 0 )
+                CALL lsm_soil_model( .FALSE., 0, .FALSE. )
 !
 !--             Call for southward-facing surfaces
                 CALL lsm_energy_balance( .FALSE., 1 )
-                CALL lsm_soil_model( .FALSE., 1 )
+                CALL lsm_soil_model( .FALSE., 1, .FALSE. )
 !
 !--             Call for eastward-facing surfaces
                 CALL lsm_energy_balance( .FALSE., 2 )
-                CALL lsm_soil_model( .FALSE., 2 )
+                CALL lsm_soil_model( .FALSE., 2, .FALSE. )
 !
 !--             Call for westward-facing surfaces
                 CALL lsm_energy_balance( .FALSE., 3 )
-                CALL lsm_soil_model( .FALSE., 3 )
+                CALL lsm_soil_model( .FALSE., 3, .FALSE. )
 
                 CALL cpu_log( log_point(54), 'land_surface', 'stop' )
@@ -262,6 +278,5 @@
 !--       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
+               == intermediate_timestep_count_max )  THEN
 
                time_radiation = time_radiation + dt_spinup
@@ -292,6 +307,5 @@
 !
 !--    Increase simulation time and output times
-       nr_timesteps_this_run      = nr_timesteps_this_run + 1
-       current_timestep_number    = current_timestep_number + 1
+       current_timestep_number_spinup = current_timestep_number_spinup + 1
        simulated_time             = simulated_time   + dt_spinup
        simulated_time_chr         = time_to_string( simulated_time )
@@ -350,19 +364,36 @@
 !--    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
+!        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' )
+
+
+!
+!--    Run control output
+       IF ( myid == 0 )  THEN
+!
+!--       If necessary, write header
+          IF ( .NOT. run_control_header_spinup )  THEN
+             CALL check_open( 15 )
+             WRITE ( 15, 100 )
+             run_control_header_spinup = .TRUE.
+          ENDIF
+!
+!--       Write some general information about the spinup in run control file
+          WRITE ( 15, 101 )  current_timestep_number_spinup, simulated_time_chr, dt_spinup, pt_spinup, rad_sw_in(0,nysg,nxlg)
+!
+!--       Write buffer contents to disc immediately
+          FLUSH( 15 )
        ENDIF
 
-       CALL cpu_log( log_point_s(15), 'timesteps spinup', 'stop' )
-
-       IF ( myid == 0 )  THEN
-          PRINT*, time_since_reference_point
-       ENDIF
+
 
     ENDDO   ! time loop
@@ -375,5 +406,61 @@
     DEALLOCATE(pt_save)
 
-    CALL location_message( 'finished time-stepping spinup', .TRUE. )
+    CALL location_message( 'finished spinup-sequence', .TRUE. )
+
+
+!
+!-- Formats
+100 FORMAT (///'Spinup control output:'/  &
+            '----------------------------------------'// &
+            'ITER.  HH:MM:SS    DT   PT(z_MO)     SWD'/   &
+            '----------------------------------------')
+101 FORMAT (I5,2X,A9,1X,F6.2,3X,F6.2,2X,F6.2)
+
+ CONTAINS
+
+!
+!-- Returns the cosine of the solar zenith angle at a given time. This routine
+!-- is similar to that for calculation zenith (see radiation_model_mod.f90)
+    FUNCTION solar_angle( local_time ) 
+
+       USE constants,                                                          &
+           ONLY:  pi
+       
+       USE kinds
+
+       USE radiation_model_mod,                                                &
+           ONLY:  day_init, decl_1, decl_2, decl_3, lat, lon
+
+       IMPLICIT NONE 
+
+
+       REAL(wp) ::  solar_angle     !< cosine of the solar zenith angle
+
+       REAL(wp) ::  day             !< day of the year
+       REAL(wp) ::  declination     !< solar declination angle
+       REAL(wp) ::  hour_angle      !< solar hour angle
+       REAL(wp) ::  time_utc        !< current time in UTC
+       REAL(wp), INTENT(IN) :: local_time
+!
+!--    Calculate current day and time based on the initial values and simulation
+!--    time
+
+       day = day_init + INT(FLOOR( local_time / 86400.0_wp ), KIND=iwp)
+       time_utc = MOD(local_time, 86400.0_wp)
+
+
+!
+!--    Calculate solar declination and hour angle   
+       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
+       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
+
+!
+!--    Calculate cosine of solar zenith angle
+       solar_angle = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &
+                     * COS(hour_angle)
+
+
+    END FUNCTION solar_angle
+
 
  END SUBROUTINE time_integration_spinup
