Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 1499)
+++ /palm/trunk/SOURCE/Makefile	(revision 1500)
@@ -20,5 +20,5 @@
 # Current revisions:
 # ------------------
-# 
+# Bugfix: missing adjustments for land surface model and radiation model
 # 
 # Former revisions:
@@ -298,5 +298,5 @@
 check_open.o: modules.o mod_kinds.o mod_particle_attributes.o
 check_parameters.o: modules.o mod_kinds.o subsidence.o land_surface_model.o\
-	plant_canopy_model.o
+	plant_canopy_model.o radiation_model.o
 close_file.o: modules.o mod_kinds.o
 compute_vpt.o: modules.o mod_kinds.o
@@ -398,5 +398,5 @@
 poisfft.o: modules.o cpulog.o fft_xy.o mod_kinds.o tridia_solver.o
 poismg.o: modules.o cpulog.o mod_kinds.o
-prandtl_fluxes.o: modules.o mod_kinds.o
+prandtl_fluxes.o: modules.o mod_kinds.o land_surface_model.o
 pres.o: modules.o cpulog.o mod_kinds.o poisfft.o
 print_1d.o: modules.o cpulog.o mod_kinds.o
@@ -422,5 +422,5 @@
 sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o
 surface_coupler.o: modules.o cpulog.o mod_kinds.o
-swap_timelevel.o: modules.o cpulog.o mod_kinds.o
+swap_timelevel.o: modules.o cpulog.o mod_kinds.o land_surface_model.o
 temperton_fft.o: modules.o mod_kinds.o
 time_integration.o: modules.o advec_ws.o buoyancy.o calc_mean_profile.o \
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 1499)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 1500)
@@ -958,12 +958,14 @@
        ENDIF   
 
-!      Dirichlet boundary conditions are allowed at the moment for testing 
-!      purposes
-!        IF ( bc_pt_b == 'dirichlet' .OR. bc_q_b == 'dirichlet' )  THEN
-!           message_string = 'lsm requires setting of'//                         &
-!                            'bc_pt_b = "neumann" and '//                        &
-!                            'bc_q_b  = "neumann"'
-!           CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
-!        ENDIF
+!
+!--    Dirichlet boundary conditions are required as the surface fluxes are
+!--    calculated from the temperature/humidity gradients in the land surface
+!--    model
+       IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' )  THEN
+          message_string = 'lsm requires setting of'//                         &
+                           'bc_pt_b = "dirichlet" and '//                        &
+                           'bc_q_b  = "dirichlet"'
+          CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
+       ENDIF
 
        IF ( .NOT. prandtl_layer )  THEN
Index: /palm/trunk/SOURCE/land_surface_model.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model.f90	(revision 1499)
+++ /palm/trunk/SOURCE/land_surface_model.f90	(revision 1500)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! -----------------
-! 
+! Corrected calculation of aerodynamic resistance (r_a).
+! Precipitation is now added to liquid water reservoir using LE_liq.
+! Added support for dry runs.
 ! 
 ! Former revisions:
@@ -37,4 +39,15 @@
 ! H-TESSEL. The implementation is based on the formulation implemented in the 
 ! DALES model.
+!
+! To do list:
+! -----------
+! - Add support for binary I/O support
+! - Add support for lsm data output
+! - Check for time step criterion
+! - Check use with RK-2 and Euler time-stepping
+! - Adaption for use with cloud physics (liquid water potential temperature)
+! - Check reaction of plants at wilting point and at atmospheric saturation
+! - Consider partial absorption of the net shortwave radiation by the skin layer
+! - Allow for water surfaces, check performance for bare soils 
 !------------------------------------------------------------------------------!
      USE arrays_3d,                                                            &
@@ -42,10 +55,10 @@
 
      USE cloud_parameters,                                                     &
-         ONLY: cp, l_d_r, l_v, rho_l, r_d, r_v
+         ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
 
      USE control_parameters,                                                   &
          ONLY: dt_3d, humidity, intermediate_timestep_count,                   &
-               intermediate_timestep_count_max, pt_surface, rho_surface,       &
-               surface_pressure, timestep_scheme, tsc
+               intermediate_timestep_count_max, precipitation, pt_surface,     &
+               rho_surface, surface_pressure, timestep_scheme, tsc
 
      USE indices,                                                              &
@@ -119,5 +132,5 @@
               root_fraction = (/0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp/), & !: distribution of root surface area to the individual soil layers
               soil_level = (/0.07_wp, 0.28_wp, 1.00_wp,  2.89_wp/),   & !: soil layer depths (m)
-              soil_moisture    = 0.0_wp       !: soil moisture content (m3/m3)
+              soil_moisture = 0.0_wp          !: soil moisture content (m3/m3)
 
     REAL(wp), DIMENSION(0:soil_layers) ::   &
@@ -532,5 +545,11 @@
        G       = 0.0_wp
        H       = rho_cp * shf
-       LE      = rho_l * l_v * qsws
+
+       IF ( humidity )  THEN
+          LE = rho_l * l_v * qsws
+       ELSE
+          LE = 0.0_wp
+       ENDIF
+
        LE_veg  = 0.0_wp
        LE_soil = LE
@@ -658,4 +677,8 @@
           DO  j = nysg, nyng
              k = nzb_s_inner(j,i)
+!
+!--          Assure that r_a cannot be zero at model start
+             IF ( pt(k+1,j,i) == pt(k,j,i) )  pt(k+1,j,i) = pt(k+1,j,i) + 1.0E-10_wp
+
              us(j,i) = 0.1_wp
              ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i)
@@ -717,5 +740,5 @@
        DO i = nxlg, nxrg
           DO j = nysg, nyng
-
+             k = nzb_s_inner(j,i)
 
 !
@@ -726,12 +749,11 @@
                 lambda_skin = lambda_skin_u(j,i)
              ENDIF
-!
-!--          First step: calculate aerodyamic resistance. As pt(0), us, ts
-!--          are not available for the current time step, data from the last
-!--          time step is used here.
-             k = nzb_s_inner(j,i)
-
-!            r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20)
-             r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / - (shf(j,i) + 1.0E-20)
+
+!
+!--          First step: calculate aerodyamic resistance. As pt, us, ts
+!--          are not available for the prognostic time step, data from the last
+!--          time step is used here. Note that this formulation is the
+!--          equivalent to the ECMWF formulation using drag coefficients
+             r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20)
 
 !
@@ -793,8 +815,13 @@
 
              q_s = 0.622_wp * e_s / surface_pressure
-             IF ( q_s .LE. q_p(k+1,j,i))  THEN
-!               PRINT*, "dew fall at (before)", time_since_reference_point
-                r_canopy(j,i) = 0.0_wp
-                r_soil(j,i) = 0.0_wp
+
+!
+!--          In case of dew fall, set resistances to zero.
+!--          To do: what does that physically reasoning behind this?
+             IF ( humidity )  THEN
+                IF ( q_s .LE. q_p(k+1,j,i) )  THEN
+                   r_canopy(j,i) = 0.0_wp
+                   r_soil(j,i) = 0.0_wp
+                ENDIF
              ENDIF 
 
@@ -808,9 +835,10 @@
 
 
-!            Plant cannot transpirate below wilting point. here, r_canopy
-!            should go to infinity
-!              IF ( m_soil(k,j,i) .LT. m_wilt(j,i) )  THEN
-!                 f_LE_veg(j,i) = 0.0
-!              ENDIF
+!
+!--          If soil moisture is below wilting point, plants do no longer
+!--          transpirate.
+             IF ( m_soil(k,j,i) .LT. m_wilt(j,i) )  THEN
+                f_LE_veg = 0.0
+             ENDIF
 
              f_H  = rho_cp / r_a(j,i)
@@ -831,14 +859,30 @@
              Rn(j,i) = Rn(j,i) + sigma_SB * T_0(j,i) ** 4
 
-!
-!--          Numerator of the prognostic equation
-             coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H / exn  &
-                      * T_1 + f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT * T_0(j,i) &
-                                     ) + lambda_skin * T_soil(0,j,i)
-
-!
-!--          Denominator of the prognostic equation
-             coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 + f_LE * dq_s_dT +     &
-                      lambda_skin + f_H / exn
+             IF ( humidity )  THEN
+
+!
+!--             Numerator of the prognostic equation
+                coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H     &
+                         / exn * T_1 + f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT   &
+                         * T_0(j,i) ) + lambda_skin * T_soil(0,j,i)
+
+!
+!--             Denominator of the prognostic equation
+                coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 + f_LE * dq_s_dT    &
+                         + lambda_skin + f_H / exn
+
+             ELSE
+
+!
+!--             Numerator of the prognostic equation
+                coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H /   &
+                         exn * T_1 + lambda_skin * T_soil(0,j,i)
+
+!
+!--             Denominator of the prognostic equation
+                coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3                     &
+                         + lambda_skin + f_H / exn
+
+             ENDIF
 
              tend = 0.0_wp
@@ -877,14 +921,16 @@
              G(j,i)         = lambda_skin * (T_0_p(j,i) - T_soil(0,j,i))
              H(j,i)         = - f_H  * ( pt_p(k+1,j,i) - pt_p(k,j,i) )
-             LE(j,i)        = - f_LE      * ( q_p(k+1,j,i) - q_s + dq_s_dT *   &
-                                T_0(j,i) - dq_s_dT * T_0_p(j,i) )
-
-             LE_veg(j,i)    = - f_LE_veg  * ( q_p(k+1,j,i) - q_s + dq_s_dT *   &
-                                T_0(j,i) - dq_s_dT * T_0_p(j,i) )
-             LE_soil(j,i)   = - f_LE_soil * ( q_p(k+1,j,i) - q_s + dq_s_dT *   &
-                                T_0(j,i) - dq_s_dT * T_0_p(j,i) )
-             LE_liq(j,i)    = - f_LE_liq  * ( q_p(k+1,j,i) - q_s + dq_s_dT *   &
-                                T_0(j,i) - dq_s_dT * T_0_p(j,i) )
-
+
+             IF ( humidity )  THEN
+                LE(j,i)        = - f_LE      * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
+                                   * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
+
+                LE_veg(j,i)    = - f_LE_veg  * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
+                                   * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
+                LE_soil(j,i)   = - f_LE_soil * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
+                                   * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
+                LE_liq(j,i)    = - f_LE_liq  * ( q_p(k+1,j,i) - q_s + dq_s_dT  &
+                                   * T_0(j,i) - dq_s_dT * T_0_p(j,i) )
+             ENDIF
 
 !              IF ( i == 1 .AND. j == 1 )  THEN
@@ -898,7 +944,7 @@
 !              ENDIF
 
-
+!
+!--          Calculate the true surface resistance
              IF ( LE(j,i) .EQ. 0.0 )  THEN
-!               PRINT*, "+++ Evapotranspiration -> 0"
                 r_s(j,i) = 1.0E10
              ELSE
@@ -908,62 +954,77 @@
 
 !
-!--          Calculate change in liquid water reservoir due to dew fall or 
-!--          evaporation of liquid water (to do: add interception from rainfall)
-             IF ( q_s .LE. q_p(k+1,j,i))  THEN
-!
-!--             Check if reservoir is full (avoid values > m_liq_max)
-!--             In that case, LE_liq goes to LE_soil. In this case
-!--             LE_veg is zero anyway (because c_liq = 1), so that tend is
-!--             zero and no further check is needed
-                IF ( m_liq(j,i) .EQ. m_liq_max )  THEN
-                   LE_soil(j,i) = LE_soil(j,i) + LE_liq(j,i)
-                   LE_liq(j,i) = 0.0_wp
-                ENDIF
-
-!
-!--             In case LE_veg becomes negative (unphysical behavior), let
-!--             the water enter the liquid water reservoir as dew on the
-!--             plant
-                IF ( LE_veg(j,i) .LT. 0.0_wp )  THEN
-                   LE_liq(j,i) = LE_liq(j,i) + LE_veg(j,i)
-                   LE_veg(j,i) = 0.0_wp
-                ENDIF
-             ENDIF                    
-  
-             tend = - LE_liq(j,i) * drho_l_lv
-  
-
-             m_liq_p(j,i) = m_liq(j,i) + dt_3d * ( tsc(2) * tend               &
-                                                   + tsc(3) * tm_liq_m(j,i) )
-
-!
-!--          Check if reservoir is overfull -> reduce to maximum
-!--          (conservation of water is violated here)
-             m_liq_p(j,i) = MIN(m_liq_p(j,i),m_liq_max)
-
-!
-!--          Check if reservoir is empty (avoid values < 0.0)
-!--          (conservation of water is violated here)
-             m_liq_p(j,i) = MAX(m_liq_p(j,i),0.0_wp)
-
-
-!
-!--          Calculate m_liq tendencies for the next Runge-Kutta step
-             IF ( timestep_scheme(1:5) == 'runge' )  THEN
-                IF ( intermediate_timestep_count == 1 )  THEN
-                   tm_liq_m(j,i) = tend
-                ELSEIF ( intermediate_timestep_count <                         &
-                         intermediate_timestep_count_max )  THEN
-                   tm_liq_m(j,i) = -9.5625_wp * tend + 5.3125_wp * tm_liq_m(j,i)
-                ENDIF
-             ENDIF
-
-!
 !--          Calculate fluxes in the atmosphere
              shf(j,i) = H(j,i) / rho_cp
-             qsws(j,i) = LE(j,i) / rho_lv
-
-             ENDDO
+
+!
+!--          Calculate change in liquid water reservoir due to dew fall or 
+!--          evaporation of liquid water
+             IF ( humidity )  THEN
+!
+!--             If precipitation is activated, add rain water to LE_liq.
+!--             precipitation_rate is given in mm.
+                IF ( precipitation )  THEN
+                   LE_liq(j,i) = LE_liq(j,i) + precipitation_rate(j,i)         &
+                                               * 0.001_wp * rho_l * l_v
+                ENDIF
+!
+!--             If the air is saturated, check the reservoir water level
+                IF ( q_s .LE. q_p(k+1,j,i))  THEN
+!
+!--                Check if reservoir is full (avoid values > m_liq_max)
+!--                In that case, LE_liq goes to LE_soil. In this case
+!--                LE_veg is zero anyway (because c_liq = 1), so that tend is
+!--                zero and no further check is needed
+                   IF ( m_liq(j,i) .EQ. m_liq_max )  THEN
+                      LE_soil(j,i) = LE_soil(j,i) + LE_liq(j,i)
+                      LE_liq(j,i) = 0.0_wp
+                   ENDIF
+
+!
+!--                In case LE_veg becomes negative (unphysical behavior), let
+!--                the water enter the liquid water reservoir as dew on the
+!--                plant
+                   IF ( LE_veg(j,i) .LT. 0.0_wp )  THEN
+                      LE_liq(j,i) = LE_liq(j,i) + LE_veg(j,i)
+                      LE_veg(j,i) = 0.0_wp
+                   ENDIF
+                ENDIF                    
+  
+                tend = - LE_liq(j,i) * drho_l_lv
+
+                m_liq_p(j,i) = m_liq(j,i) + dt_3d * ( tsc(2) * tend            &
+                                                   + tsc(3) * tm_liq_m(j,i) )
+
+!
+!--             Check if reservoir is overfull -> reduce to maximum
+!--             (conservation of water is violated here)
+                m_liq_p(j,i) = MIN(m_liq_p(j,i),m_liq_max)
+
+!
+!--             Check if reservoir is empty (avoid values < 0.0)
+!--             (conservation of water is violated here)
+                m_liq_p(j,i) = MAX(m_liq_p(j,i),0.0_wp)
+
+
+!
+!--             Calculate m_liq tendencies for the next Runge-Kutta step
+                IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( intermediate_timestep_count == 1 )  THEN
+                      tm_liq_m(j,i) = tend
+                   ELSEIF ( intermediate_timestep_count <                      &
+                            intermediate_timestep_count_max )  THEN
+                      tm_liq_m(j,i) = -9.5625_wp * tend + 5.3125_wp            &
+                                      * tm_liq_m(j,i)
+                   ENDIF
+                ENDIF
+
+!
+!--             Calculate moisture flux in the atmosphere
+                qsws(j,i) = LE(j,i) / rho_lv
+
+             ENDIF
+
           ENDDO
+       ENDDO
 
 
@@ -1221,5 +1282,4 @@
        DO i = nxlg, nxrg    
           DO j = nysg, nyng
-
              k = nzb_s_inner(j,i)
 !
