Index: palm/trunk/SOURCE/time_integration_spinup.f90
===================================================================
--- palm/trunk/SOURCE/time_integration_spinup.f90	(revision 2297)
+++ 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
