Index: palm/trunk/DOC/app/appendix_a.html
===================================================================
--- palm/trunk/DOC/app/appendix_a.html	(revision 410)
+++ palm/trunk/DOC/app/appendix_a.html	(revision 411)
@@ -3535,4 +3535,11 @@
       </td>
     </tr>
+    <tr>
+      <td style="vertical-align: top;"><a name="PA0324"></a>PA0324<br>
+      </td>
+      <td style="vertical-align: top;">Applying large scale vertical motion is not 
+                          allowed for ocean runs<br>
+      </td>
+    </tr>
 <tr>
 
Index: palm/trunk/DOC/app/chapter_4.1.html
===================================================================
--- palm/trunk/DOC/app/chapter_4.1.html	(revision 410)
+++ palm/trunk/DOC/app/chapter_4.1.html	(revision 411)
@@ -17249,6 +17249,104 @@
 
 
-
- 
+ </tr>
+
+ <tr>
+
+ <td style="vertical-align: top;"> 
+     
+  <p><a name="ws_vertical_gradient"></a><b>ws_vertical_gradient</b></p>
+
+ </td>
+
+  <td style="vertical-align: top;">R (10)</td>
+
+  <td style="vertical-align: top;"><i>10 * 0.0</i></td>
+
+  <td style="vertical-align: top;"> 
+      
+    <p>Gradient(s) of the profile for the large scale subsidence/ascent velocity
+       (in (m/s) / 100 m).&nbsp; </p>
+      
+      <p>This gradient holds starting from the height&nbsp;
+        level defined by <a href="#ws_vertical_gradient_level">ws_vertical_gradient_level</a>
+       (precisely: for all uv levels k where zu(k) &gt;
+       ws_vertical_gradient_level,
+       w_subs(k) is set: w_subs(k) = w_subs(k-1) + dzu(k) * <b>ws_vertical_gradient</b>)
+       up to the top boundary or up to the next height level defined
+       by <a href="#ws_vertical_gradient_level">ws_vertical_gradient_level</a>.
+       A total of 10 different gradients for 11 height intervals (10 intervals
+       if <a href="#ws_vertical_gradient_level">ws_vertical_gradient_level</a>(1)
+       = <i>0.0</i>) can be assigned. &nbsp;
+      </p>
+      
+      <p>Example:&nbsp; </p>
+      
+      <ul>
+      
+        <p><b>ws_vertical_gradient</b>
+= <i>-0.002</i>, <i>0.0</i>,&nbsp; <br>
+
+        <b>ws_vertical_gradient_level</b> = <i>0.0</i>,
+        <i>1000.0</i>,</p>
+      
+      </ul>   
+      
+      <p>That defines the subsidence/ascent profile to be linear up
+         to z = 1000.0 m with a surface value of 0 m/s.
+         For z &gt; 1000.0 m up to the top boundary the gradient is <i>0.0</i> (m/s) / 100 m
+         (it is assumed that the assigned height levels correspond
+         with uv levels).</p>
+            
+       <p>With an appropriate construction of w_subs the height of the boundary layer z_i 
+        can be kept approximately constant.</p>
+      
+   <p><span style="font-weight: bold;">Attention:</span><br>
+
+      The large scale vertical motion is only applied to the prognostic equation for the potential 
+      temperature because it cannot be applied to the momentum equations due to incompressibility. 
+      Thus, the model is not mass consistent.</p>
+
+
+ </td>
+
+ </tr>
+
+
+ </tr>
+
+ <tr>
+
+ <td style="vertical-align: top;">   
+      
+      <p><a name="ws_vertical_gradient_level"></a><b>ws_vertical_gradient</b>
+      <br><b>_level</b></p>
+
+ </td>
+
+ <td style="vertical-align: top;">R (10)</td>
+
+ <td style="vertical-align: top;"> 
+         
+      <p><i>10 *</i>&nbsp;
+      <span style="font-style: italic;">0.0</span><br>
+      </p>
+
+ </td>
+
+
+ <td style="vertical-align: top;">
+           
+      <p>Height level from which on the gradient for the subsidence/ascent velocity  defined by
+      <a href="#ws_vertical_gradient">ws_vertical_gradient</a>
+      is effective (in m).&nbsp; </p>
+   
+      <p>The height levels have to be assigned in ascending order. The
+default values result in a profile which is zero everywhere regardless of the
+values of <a href="#ws_vertical_gradient">ws_vertical_gradient</a>.
+For the piecewise construction of temperature profiles see <a href="#ws_vertical_gradient">ws_vertical_gradient</a>.</p>
+      </td>
+
+ </tr>
+
   
   
Index: palm/trunk/DOC/app/chapter_4.6.html
===================================================================
--- palm/trunk/DOC/app/chapter_4.6.html	(revision 410)
+++ palm/trunk/DOC/app/chapter_4.6.html	(revision 411)
@@ -15149,4 +15149,33 @@
 
  </tr>
+
+  <tr>
+
+ <td><a href="chapter_4.1.html#ws_vertical_gradient"><span style="font-weight: bold;">ws_vertical_gradient</span></a></td>
+
+ <td>I</td>
+
+ <td>R(10)</td>
+
+ <td><span style="font-style: italic;">10 * 0.0</span></td>
+
+ <td>Gradient(s) of the profile for the large scale subsidence/ascent velocity in ((m/s) / 100 m)</td>
+
+ </tr>
+
+  <tr>
+
+ <td><a href="chapter_4.1.html#ws_vertical_gradient_level"><span style="font-weight: bold;">ws_vertical_gradient_level</span></a></td>
+
+ <td>I</td>
+
+ <td>R(10)</td>
+
+ <td><span style="font-style: italic;">10 * 0.0</span></td>
+
+ <td>Height level from which on the gradient for the subsidence/ascent velocity defined by 
+     ws_vertical_gradient is effective (in m).</td>
+
+ </tr>
  
  <tr>
Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 410)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 411)
@@ -5,5 +5,10 @@
 (check_parameters, init_3d_model, parin)
 
-check parameters, init_3d_model, parin
+Large scale vertical motion (subsidence/ascent) can be applied to the 
+prognostic equation for the potential temperature. (check_parameters, header, 
+Makefile, modules, parin, prognostic_equations, read_var_list, subsidence,
+write_var_list)
+
+check parameters, header, init_3d_model, Makefile, modules, parin prognostic equations, read_var_list, subsidence, write_var_list
 
 Changed:
Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 410)
+++ palm/trunk/SOURCE/Makefile	(revision 411)
@@ -5,4 +5,5 @@
 # -----------------
 # +set_particle_attributes, set_slicer_attributes_dvrp
+# +subsidence
 #
 # Branch revisions:
@@ -81,5 +82,6 @@
 	set_particle_attributes.f90 set_slicer_attributes_dvrp.f90 \
         singleton.f90 sor.f90 spline_x.f90 \
-	spline_y.f90 spline_z.f90 sum_up_3d_data.f90 surface_coupler.f90 \
+	spline_y.f90 spline_z.f90 subsidence.f90 \
+        sum_up_3d_data.f90 surface_coupler.f90 \
 	swap_timelevel.f90 temperton_fft.f90 time_integration.f90 \
         time_to_string.f90 timestep.f90 timestep_scheme_steering.f90 \
@@ -124,5 +126,6 @@
 	set_particle_attributes.o set_slicer_attributes_dvrp.o \
         singleton.o sor.o spline_x.o spline_y.o \
-        spline_z.o sum_up_3d_data.o surface_coupler.o swap_timelevel.o \
+        spline_z.o subsidence.o sum_up_3d_data.o surface_coupler.o \
+        swap_timelevel.o \
         temperton_fft.o time_integration.o time_to_string.o \
         timestep.o timestep_scheme_steering.o transpose.o \
@@ -188,5 +191,5 @@
 check_for_restart.o: modules.o
 check_open.o: modules.o
-check_parameters.o: modules.o
+check_parameters.o: modules.o subsidence.o
 close_file.o: modules.o
 compute_vpt.o: modules.o
@@ -260,5 +263,5 @@
 	eqn_state_seawater.o impact_of_latent_heat.o \
         plant_canopy_model.o production_e.o \
-	user_actions.o
+	subsidence.o user_actions.o
 random_gauss.o: random_function.o
 read_3d_binary.o: modules.o random_function.o
@@ -272,4 +275,5 @@
 spline_y.o: modules.o
 spline_z.o: modules.o
+subsidence.o: modules.o
 sum_up_3d_data.o: modules.o
 surface_coupler.o: modules.o
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 410)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 411)
@@ -5,4 +5,5 @@
 ! -----------------
 ! Enabled passive scalar/humidity wall fluxes for non-flat topography
+! Initialization of large scale vertical motion (subsidence/ascent)
 !
 ! Branch revisions:
@@ -114,4 +115,5 @@
     USE pegrid
     USE profil_parameter
+    USE subsidence_mod
     USE statistics
     USE transpose_indices
@@ -983,4 +985,12 @@
          
     ENDIF
+
+!
+!-- Initialize large scale subsidence if required
+    IF ( ws_vertical_gradient_level(1) /= -9999999.9 )  THEN
+       large_scale_subsidence = .TRUE.
+       CALL init_w_subsidence
+    END IF
+ 
              
 !
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 410)
+++ palm/trunk/SOURCE/header.f90	(revision 411)
@@ -4,4 +4,5 @@
 ! Current revisions:
 ! -----------------
+! +large scale vertical motion (subsidence/ascent)
 ! 
 !
@@ -84,5 +85,5 @@
 ! end of the run, then writing additional information about CPU-usage on file
 ! header.
-!------------------------------------------------------------------------------!
+!-----------------------------------------------------------------------------!
 
     USE arrays_3d
@@ -96,4 +97,5 @@
     USE particle_attributes
     USE pegrid
+    USE subsidence_mod
     USE spectrum
 
@@ -323,4 +325,8 @@
           WRITE ( io, 152 )  dpdxy, dp_level_b, '.'
        ENDIF
+    ENDIF
+    IF ( large_scale_subsidence )  THEN
+        WRITE ( io, 153 )
+        WRITE ( io, 154 )
     ENDIF
     WRITE ( io, 99 )
@@ -1357,4 +1363,35 @@
 
        WRITE ( io, 425 )  TRIM( coordinates ), TRIM( temperatures ), &
+                          TRIM( gradients ), TRIM( slices )
+    ENDIF
+
+!
+!-- Profile for the large scale vertial velocity
+!-- Building output strings, starting with surface value
+    IF ( large_scale_subsidence )  THEN
+       temperatures = '   0.0'
+       gradients = '------'
+       slices = '     0'
+       coordinates = '   0.0'
+       i = 1
+       DO  WHILE ( ws_vertical_gradient_level_ind(i) /= -9999 )
+
+          WRITE (coor_chr,'(E10.2,7X)')  &
+                                w_subs(ws_vertical_gradient_level_ind(i))
+          temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(E10.2,7X)')  ws_vertical_gradient(i)
+          gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(I10,7X)')  ws_vertical_gradient_level_ind(i)
+          slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F10.2,7X)')  ws_vertical_gradient_level(i)
+          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
+
+          i = i + 1
+       ENDDO
+
+       WRITE ( io, 426 )  TRIM( coordinates ), TRIM( temperatures ), &
                           TRIM( gradients ), TRIM( slices )
     ENDIF
@@ -1568,4 +1605,7 @@
            /'     ',2(1X,E12.5),'Pa/m in x/y direction', &
            /'     starting from dp_level_b =', F8.3, 'm', A /)
+153 FORMAT (' --> Large-scale vertical motion is used in the ', &
+                  'prognostic equation for')
+154 FORMAT ('     the potential temperature')
 200 FORMAT (//' Run time and time step information:'/ &
              ' ----------------------------------'/)
@@ -1814,4 +1854,9 @@
             '       Gradient:   ',A,'  psu/100m'/ &
             '       Gridpoint:  ',A)
+426 FORMAT (/'    Characteristic levels of the subsidence/ascent profile:'// &
+            '       Height:      ',A,'  m'/ &
+            '       w_subs:      ',A,'  m/s'/ &
+            '       Gradient:    ',A,'  (m/s)/100m'/ &
+            '       Gridpoint:   ',A)
 450 FORMAT (//' LES / Turbulence quantities:'/ &
               ' ---------------------------'/)
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 410)
+++ palm/trunk/SOURCE/modules.f90	(revision 411)
@@ -8,4 +8,6 @@
 ! initial data assignments to some dvrp arrays changed due to error messages
 ! from gfortran compiler
+! +large_scale_subsidence, ws_vertical_gradient, ws_vertical_gradient_level,
+! ws_vertical_gradient_level_ind, w_subs
 !
 ! Branch revisions:
@@ -155,5 +157,5 @@
           ddzu, dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, km_damp_x,  &
           km_damp_y, lad, l_grid, pt_init, q_init, rdf, sa_init, ug, u_init,   &
-          u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, zu, zw
+          u_nzb_p1_for_vfc, vg, v_init, v_nzb_p1_for_vfc, w_subs, zu, zw
 
     REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
@@ -398,5 +400,6 @@
                 section_xz(100) = -9999, section_yz(100) = -9999, &
                 ug_vertical_gradient_level_ind(10) = -9999, &
-                vg_vertical_gradient_level_ind(10) = -9999
+                vg_vertical_gradient_level_ind(10) = -9999, &
+                ws_vertical_gradient_level_ind(10) = -9999
 
     INTEGER, DIMENSION(:), ALLOCATABLE ::  grid_level_count
@@ -426,5 +429,6 @@
                 humidity = .FALSE., humidity_remote = .FALSE., &
                 inflow_l = .FALSE., inflow_n = .FALSE., &
-                inflow_r = .FALSE., inflow_s = .FALSE., iso2d_output = .FALSE.,&
+                inflow_r = .FALSE., inflow_s = .FALSE., &
+                iso2d_output = .FALSE., large_scale_subsidence = .FALSE., &
                 mg_switch_to_pe0 = .FALSE., &
                 netcdf_output = .FALSE., netcdf_64bit = .FALSE., &
@@ -545,5 +549,7 @@
              volume_flow_initial(1:2) = 0.0, wall_heatflux(0:4) = 0.0, &
              wall_humidityflux(0:4) = 0.0, wall_qflux(0:4) = 0.0, &
-             wall_salinityflux(0:4) = 0.0, wall_scalarflux(0:4) = 0.0
+             wall_salinityflux(0:4) = 0.0, wall_scalarflux(0:4) = 0.0, &
+             ws_vertical_gradient(10) = 0.0, &
+             ws_vertical_gradient_level(10) = -9999999.9
 
     REAL, DIMENSION(:), ALLOCATABLE ::  dp_smooth_factor
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 410)
+++ palm/trunk/SOURCE/parin.f90	(revision 411)
@@ -5,4 +5,5 @@
 ! -----------------
 ! +wall_humidityflux, wall_scalarflux
+! +ws_vertical_gradient, ws_vertical_gradient_level
 !
 ! Branch revisions:
@@ -137,5 +138,6 @@
              use_upstream_for_tke, vg_surface, vg_vertical_gradient, &
              vg_vertical_gradient_level, v_bulk, wall_adjustment, &
-             wall_heatflux, wall_humidityflux, wall_scalarflux
+             wall_heatflux, wall_humidityflux, wall_scalarflux, &
+             ws_vertical_gradient, ws_vertical_gradient_level
 
 
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 410)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 411)
@@ -4,4 +4,6 @@
 ! Actual revisions:
 ! -----------------
+! add call of subsidence in the prognostic equation for
+! the potential temperature
 ! 
 !
@@ -94,4 +96,5 @@
     USE plant_canopy_model_mod
     USE production_e_mod
+    USE subsidence_mod
     USE user_actions_mod
 
@@ -452,4 +455,10 @@
           IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
              CALL plant_canopy_model( i, j, 4 )
+          ENDIF
+
+!
+!--       If required compute influence of large-scale subsidence/ascent
+          IF ( large_scale_subsidence ) THEN
+             CALL subsidence ( i, j, tend, pt, pt_init )
           ENDIF
 
@@ -1094,4 +1103,11 @@
           ENDIF
 
+
+!--       If required compute influence of large-scale subsidence/ascent
+          IF ( large_scale_subsidence ) THEN
+             CALL subsidence ( i, j, tend, pt, pt_init )
+          ENDIF
+
+
           CALL user_actions( i, j, 'pt-tendency' )
 
@@ -1680,4 +1696,8 @@
     ENDIF
 
+!--If required compute influence of large-scale subsidence/ascent
+   IF ( large_scale_subsidence ) THEN
+      CALL subsidence ( tend, pt, pt_init )
+   ENDIF
 
     CALL user_actions( 'pt-tendency' )
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 410)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 411)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! +large_scale_subsidence, ws_vertical_gradient, ws_vertical_gradient_level,
+! ws_vertical_gradient_level_ind
 !
 ! Former revisions:
@@ -346,4 +347,6 @@
           CASE ( 'last_dt_change' )
              READ ( 13 )  last_dt_change
+           CASE ( 'large_scale_subsidence' )
+             READ ( 13 )  large_scale_subsidence
           CASE ( 'leaf_surface_concentration' )
              READ ( 13 )  leaf_surface_concentration
@@ -593,4 +596,10 @@
           CASE ( 'wall_adjustment' )
              READ ( 13 )  wall_adjustment
+          CASE ( 'ws_vertical_gradient' )
+             READ ( 13 )  ws_vertical_gradient
+          CASE ( 'ws_vertical_gradient_level' )
+             READ ( 13 )  ws_vertical_gradient_level
+          CASE ( 'ws_vertical_gradient_level_ind' )
+             READ ( 13 )  ws_vertical_gradient_level_ind
           CASE ( 'w_max' )
              READ ( 13 )  w_max
Index: palm/trunk/SOURCE/subsidence.f90
===================================================================
--- palm/trunk/SOURCE/subsidence.f90	(revision 411)
+++ palm/trunk/SOURCE/subsidence.f90	(revision 411)
@@ -0,0 +1,255 @@
+ MODULE subsidence_mod
+
+!-----------------------------------------------------------------------------!
+! Current revisions:
+! -----------------
+! 
+!
+! Former revisions:
+! -----------------
+! $Id$
+!
+!
+! Initial revision
+!
+! Description:
+! ------------
+! Impact of large-scale subsidence or ascent as tendency term for use 
+! in the prognostic equation of potential temperature. This enables the 
+! construction of a constant boundary layer height z_i with time.
+!-----------------------------------------------------------------------------!
+
+
+    IMPLICIT NONE
+
+    PRIVATE
+    PUBLIC  init_w_subsidence, subsidence
+
+    INTERFACE init_w_subsidence
+       MODULE PROCEDURE init_w_subsidence
+    END INTERFACE init_w_subsidence
+
+    INTERFACE subsidence
+       MODULE PROCEDURE subsidence
+       MODULE PROCEDURE subsidence_ij
+    END INTERFACE subsidence
+
+ CONTAINS
+
+    SUBROUTINE init_w_subsidence 
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE pegrid
+       USE statistics 
+
+       IMPLICIT NONE
+
+       INTEGER :: i, k
+       REAL    :: gradient, ws_surface
+
+       IF ( .NOT. ALLOCATED( w_subs )) THEN
+          ALLOCATE( w_subs(nzb:nzt+1) )
+          w_subs = 0.0
+       ENDIF 
+
+      IF ( ocean )  THEN
+          message_string = 'Applying large scale vertical motion is not ' // &
+                           'allowed for ocean runs'
+          CALL message( 'init_w_subsidence', 'PA0324', 2, 2, 0, 6, 0 )
+       ENDIF
+
+!
+!--   Compute the profile of the subsidence/ascent velocity 
+!--   using the given gradients
+      i = 1
+      gradient = 0.0
+      ws_surface = 0.0
+      
+
+      ws_vertical_gradient_level_ind(1) = 0
+      DO  k = 1, nzt+1
+         IF ( i < 11 ) THEN
+            IF ( ws_vertical_gradient_level(i) < zu(k)  .AND. &
+                 ws_vertical_gradient_level(i) >= 0.0 )  THEN
+               gradient = ws_vertical_gradient(i) / 100.0
+               ws_vertical_gradient_level_ind(i) = k - 1
+               i = i + 1
+            ENDIF
+         ENDIF
+         IF ( gradient /= 0.0 )  THEN
+            IF ( k /= 1 )  THEN
+               w_subs(k) = w_subs(k-1) + dzu(k) * gradient
+            ELSE
+               w_subs(k) = ws_surface   + 0.5 * dzu(k) * gradient
+            ENDIF
+         ELSE
+            w_subs(k) = w_subs(k-1)
+         ENDIF
+      ENDDO
+
+!
+!--   In case of no given gradients for the subsidence/ascent velocity,
+!--   choose zero gradient
+      IF ( ws_vertical_gradient_level(1) == -9999999.9 )  THEN
+         ws_vertical_gradient_level(1) = 0.0
+      ENDIF
+
+    END SUBROUTINE init_w_subsidence
+
+
+    SUBROUTINE subsidence( tendency, var, var_init ) 
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE pegrid
+       USE statistics 
+
+       IMPLICIT NONE
+ 
+       INTEGER :: i, j, k
+
+       REAL :: tmp_grad
+    
+       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: var, tendency
+       REAL, DIMENSION(nzb:nzt+1) :: var_init, var_mod
+
+       var_mod = pt_surface
+
+!
+!--    Influence of w_subsidence on the current tendency term
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_s_inner(j,i)+1, nzt 
+                IF ( w_subs(k) < 0.0 ) THEN    ! large-scale subsidence
+                   tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
+                                     ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1)
+                ELSE                           ! large-scale ascent
+                   tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
+                                     ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)
+                ENDIF
+             ENDDO
+          ENDDO
+       ENDDO
+
+!
+!--    Shifting of the initial profile is especially necessary with Rayleigh 
+!--    damping switched on
+
+       DO  k = nzb, nzt
+          IF ( w_subs(k) < 0.0 ) THEN      ! large-scale subsidence
+             var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
+                               ( var_init(k+1) - var_init(k) ) * ddzu(k+1)
+          ENDIF
+       ENDDO
+!
+!--   At the upper boundary, the initial profile is shifted with aid of 
+!--   the gradient tmp_grad. (This is ok if the gradients are linear.)
+      IF ( w_subs(nzt) < 0.0 ) THEN
+         tmp_grad = ( var_init(nzt+1) - var_init(nzt) ) * ddzu(nzt+1)
+         var_mod(nzt+1) = var_init(nzt+1) -  &
+                              dt_3d * w_subs(nzt+1) * tmp_grad
+      ENDIF
+       
+
+      DO  k = nzt+1, nzb+1, -1
+         IF ( w_subs(k) >= 0.0 ) THEN  ! large-scale ascent
+            var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
+                               ( var_init(k) - var_init(k-1) ) * ddzu(k+1) 
+         ENDIF
+      ENDDO
+!
+!--   At the lower boundary shifting is not necessary because the 
+!--   subsidence velocity w_subs(nzb) vanishes.
+
+
+      IF ( w_subs(nzb+1) >= 0.0 ) THEN
+         var_mod(nzb) = var_init(nzb)
+      ENDIF
+
+      var_init = var_mod
+
+
+ END SUBROUTINE subsidence
+
+ SUBROUTINE subsidence_ij( i, j, tendency, var, var_init ) 
+
+       USE arrays_3d
+       USE control_parameters
+       USE grid_variables
+       USE indices
+       USE pegrid
+       USE statistics 
+
+       IMPLICIT NONE
+ 
+       INTEGER :: i, j, k
+
+       REAL :: tmp_grad
+    
+       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: var, tendency
+       REAL, DIMENSION(nzb:nzt+1) :: var_init, var_mod
+
+       var_mod = pt_surface
+
+!
+!--    Influence of w_subsidence on the current tendency term
+       DO  k = nzb_s_inner(j,i)+1, nzt 
+          IF ( w_subs(k) < 0.0 ) THEN      ! large-scale subsidence
+             tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
+                               ( var(k+1,j,i) - var(k,j,i) ) * ddzu(k+1)
+          ELSE                             ! large-scale ascent
+             tendency(k,j,i) = tendency(k,j,i) - w_subs(k) *  &
+                               ( var(k,j,i) - var(k-1,j,i) ) * ddzu(k)
+          ENDIF
+       ENDDO
+
+
+!
+!--    Shifting of the initial profile is especially necessary with Rayleigh 
+!--    damping switched on
+       IF ( i == nxl .AND. j == nys ) THEN ! shifting only once per PE 
+
+          DO  k = nzb, nzt
+             IF ( w_subs(k) < 0.0 ) THEN      ! large-scale subsidence
+                var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
+                                  ( var_init(k+1) - var_init(k) ) * ddzu(k+1)
+             ENDIF
+          ENDDO
+!
+!--       At the upper boundary, the initial profile is shifted with aid of 
+!--       the gradient tmp_grad. (This is ok if the gradients are linear.)
+          IF ( w_subs(nzt) < 0.0 ) THEN
+             tmp_grad = ( var_init(nzt+1) - var_init(nzt) ) * ddzu(nzt+1)
+             var_mod(nzt+1) = var_init(nzt+1) -  &
+                                  dt_3d * w_subs(nzt+1) * tmp_grad
+          ENDIF
+       
+
+          DO  k = nzt+1, nzb+1, -1
+             IF ( w_subs(k) >= 0.0 ) THEN  ! large-scale ascent
+                var_mod(k) = var_init(k) - dt_3d * w_subs(k) *  &
+                                   ( var_init(k) - var_init(k-1) ) * ddzu(k+1)
+             ENDIF
+          ENDDO
+!
+!--       At the lower boundary shifting is not necessary because the 
+!--       subsidence velocity w_subs(nzb) vanishes.
+
+
+          IF ( w_subs(nzb+1) >= 0.0 ) THEN
+             var_mod(nzb) = var_init(nzb)
+          ENDIF
+
+          var_init = var_mod(:) 
+
+       ENDIF
+
+ END SUBROUTINE subsidence_ij
+
+
+ END MODULE subsidence_mod
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 410)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 411)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +large_scale_subsidence, ws_vertical_gradient, ws_vertical_gradient_level,
+! ws_vertical_gradient_level_ind
 !
 ! Former revisions:
@@ -267,4 +268,6 @@
     WRITE ( 14 )  'last_dt_change                '
     WRITE ( 14 )  last_dt_change
+    WRITE ( 14 )  'large_scale_subsidence        '
+    WRITE ( 14 )  large_scale_subsidence
     WRITE ( 14 )  'leaf_surface_concentration    '
     WRITE ( 14 )  leaf_surface_concentration
@@ -511,4 +514,10 @@
     WRITE ( 14 )  'wall_adjustment               '
     WRITE ( 14 )  wall_adjustment
+    WRITE ( 14 )  'ws_vertical_gradient          '
+    WRITE ( 14 )  ws_vertical_gradient
+    WRITE ( 14 )  'ws_vertical_gradient_level    '
+    WRITE ( 14 )  ws_vertical_gradient_level
+    WRITE ( 14 )  'ws_vertical_gradient_level_ind'
+    WRITE ( 14 )  ws_vertical_gradient_level_ind
     WRITE ( 14 )  'w_max                         '
     WRITE ( 14 )  w_max
