Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 86)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 87)
@@ -1,11 +1,14 @@
 New:
 ---
+Calculation and output of user-defined profiles. New &userpar parameters data_output_pr_user and max_pr_user.
+
+check_parameters, flow_statistics, modules, parin, read_var_list, user_interface, write_var_list
 
 
 Changed:
 -------
-Division through dt_3d replaced by multiplication of the inverse. For performance optimisation, this is done in the loop calculating the divergence instead of using a seperate loop. (pres.f90)
+Division through dt_3d replaced by multiplication of the inverse. For performance optimisation, this is done in the loop calculating the divergence instead of using a seperate loop. (pres.f90) var_hom and var_sum renamed pr_palm.
 
-pres
+data_output_profiles, flow_statistics, init_3d_model, modules, parin, pres, read_var_list, run_control, time_integration
 
 
@@ -16,3 +19,5 @@
 Bugfix: field_chr renamed field_char (user_interface).
 
-poisfft, user_interface
+Bugfix: output of use_upstream_for_tke (header).
+
+header, poisfft, user_interface
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 86)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 87)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Check for user-defined profiles
 !
 ! Former revisions:
@@ -1214,4 +1214,5 @@
           CASE ( 'u', '#u' )
              dopr_index(i) = 1
+             dopr_unit(i)  = 'm/s'
              hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
              IF ( data_output_pr(i)(1:1) == '#' )  THEN
@@ -1223,5 +1224,6 @@
           CASE ( 'v', '#v' )
              dopr_index(i) = 2
-             hom(:,2,2,:) = SPREAD( zu, 2, statistic_regions+1 )
+             dopr_unit(i)  = 'm/s'
+             hom(:,2,2,:)  = SPREAD( zu, 2, statistic_regions+1 )
              IF ( data_output_pr(i)(1:1) == '#' )  THEN
                 dopr_initial_index(i) = 6
@@ -1232,9 +1234,251 @@
           CASE ( 'w' )
              dopr_index(i) = 3
-             hom(:,2,3,:) = SPREAD( zw, 2, statistic_regions+1 )
+             dopr_unit(i)  = 'm/s'
+             hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'pt', '#pt' )
              IF ( .NOT. cloud_physics ) THEN
                 dopr_index(i) = 4
+                dopr_unit(i)  = 'K'
+                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
+                IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                   dopr_initial_index(i) = 7
+                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
+                   hom(nzb,2,7,:)        = 0.0    ! because zu(nzb) is negative
+                   data_output_pr(i)     = data_output_pr(i)(2:)
+                ENDIF
+             ELSE
+                dopr_index(i) = 43
+                dopr_unit(i)  = 'K'
+                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
+                IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                   dopr_initial_index(i) = 28
+                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
+                   hom(nzb,2,28,:)       = 0.0    ! because zu(nzb) is negative
+                   data_output_pr(i)     = data_output_pr(i)(2:)
+                ENDIF
+             ENDIF
+
+          CASE ( 'e' )
+             dopr_index(i)  = 8
+             dopr_unit(i)   = 'm2/s2'
+             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
+             hom(nzb,2,8,:) = 0.0
+
+          CASE ( 'km', '#km' )
+             dopr_index(i)  = 9
+             dopr_unit(i)   = 'm2/s'
+             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
+             hom(nzb,2,9,:) = 0.0
+             IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                dopr_initial_index(i) = 23
+                hom(:,2,23,:)         = hom(:,2,9,:)
+                data_output_pr(i)     = data_output_pr(i)(2:)
+             ENDIF
+
+          CASE ( 'kh', '#kh' )
+             dopr_index(i)   = 10
+             dopr_unit(i)    = 'm2/s'
+             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
+             hom(nzb,2,10,:) = 0.0
+             IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                dopr_initial_index(i) = 24
+                hom(:,2,24,:)         = hom(:,2,10,:)
+                data_output_pr(i)     = data_output_pr(i)(2:)
+             ENDIF
+
+          CASE ( 'l', '#l' )
+             dopr_index(i)   = 11
+             dopr_unit(i)    = 'm'
+             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
+             hom(nzb,2,11,:) = 0.0
+             IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                dopr_initial_index(i) = 25
+                hom(:,2,25,:)         = hom(:,2,11,:)
+                data_output_pr(i)     = data_output_pr(i)(2:)
+             ENDIF
+
+          CASE ( 'w"u"' )
+             dopr_index(i) = 12
+             dopr_unit(i)  = 'm2/s2'
+             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
+             IF ( prandtl_layer )  hom(nzb,2,12,:) = zu(1)
+
+          CASE ( 'w*u*' )
+             dopr_index(i) = 13
+             dopr_unit(i)  = 'm2/s2'
+             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'w"v"' )
+             dopr_index(i) = 14
+             dopr_unit(i)  = 'm2/s2'
+             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
+             IF ( prandtl_layer )  hom(nzb,2,14,:) = zu(1)
+
+          CASE ( 'w*v*' )
+             dopr_index(i) = 15
+             dopr_unit(i)  = 'm2/s2'
+             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'w"pt"' )
+             dopr_index(i) = 16
+             dopr_unit(i)  = 'K m/s'
+             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'w*pt*' )
+             dopr_index(i) = 17
+             dopr_unit(i)  = 'K m/s'
+             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'wpt' )
+             dopr_index(i) = 18
+             dopr_unit(i)  = 'K m/s'
+             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'wu' )
+             dopr_index(i) = 19
+             dopr_unit(i)  = 'm2/s2'
+             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
+             IF ( prandtl_layer )  hom(nzb,2,19,:) = zu(1)
+
+          CASE ( 'wv' )
+             dopr_index(i) = 20
+             dopr_unit(i)  = 'm2/s2'
+             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
+             IF ( prandtl_layer )  hom(nzb,2,20,:) = zu(1)
+
+          CASE ( 'w*pt*BC' )
+             dopr_index(i) = 21
+             dopr_unit(i)  = 'K m/s'
+             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'wptBC' )
+             dopr_index(i) = 22
+             dopr_unit(i)  = 'K m/s'
+             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'u*2' )
+             dopr_index(i) = 30
+             dopr_unit(i)  = 'm2/s2'
+             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
+
+          CASE ( 'v*2' )
+             dopr_index(i) = 31
+             dopr_unit(i)  = 'm2/s2'
+             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
+
+          CASE ( 'w*2' )
+             dopr_index(i) = 32
+             dopr_unit(i)  = 'm2/s2'
+             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'pt*2' )
+             dopr_index(i) = 33
+             dopr_unit(i)  = 'K2'
+             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
+
+          CASE ( 'e*' )
+             dopr_index(i) = 34
+             dopr_unit(i)  = 'm2/s2'
+             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
+
+          CASE ( 'w*2pt*' )
+             dopr_index(i) = 35
+             dopr_unit(i)  = 'K m2/s2'
+             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'w*pt*2' )
+             dopr_index(i) = 36
+             dopr_unit(i)  = 'K2 m/s'
+             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'w*e*' )
+             dopr_index(i) = 37
+             dopr_unit(i)  = 'm3/s3'
+             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'w*3' )
+             dopr_index(i) = 38
+             dopr_unit(i)  = 'm3/s3'
+             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'Sw' )
+             dopr_index(i) = 39
+             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
+
+          CASE ( 'q', '#q' )
+             IF ( .NOT. cloud_physics )  THEN
+                IF ( myid == 0 )  THEN
+                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
+                           data_output_pr(i),                          &
+                           '    is not implemented for cloud_physics = FALSE'
+                ENDIF
+                CALL local_stop
+             ELSE
+                dopr_index(i) = 41
+                dopr_unit(i)  = 'kg/kg'
+                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
+                IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                   dopr_initial_index(i) = 26
+                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
+                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
+                   data_output_pr(i)     = data_output_pr(i)(2:)
+                ENDIF
+             ENDIF
+
+          CASE ( 's', '#s' )
+             IF ( .NOT. passive_scalar )  THEN
+                IF ( myid == 0 )  THEN
+                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
+                           data_output_pr(i),                          &
+                           '    is not implemented for passive_scalar = FALSE'
+                ENDIF
+                CALL local_stop
+             ELSE
+                dopr_index(i) = 41
+                dopr_unit(i)  = 'kg/m3'
+                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
+                IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                   dopr_initial_index(i) = 26
+                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
+                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
+                   data_output_pr(i)     = data_output_pr(i)(2:)
+                ENDIF
+             ENDIF
+
+          CASE ( 'qv', '#qv' )
+             IF ( .NOT. cloud_physics ) THEN
+                dopr_index(i) = 41
+                dopr_unit(i)  = 'kg/kg'
+                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
+                IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                   dopr_initial_index(i) = 26
+                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
+                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
+                   data_output_pr(i)     = data_output_pr(i)(2:)
+                ENDIF
+             ELSE
+                dopr_index(i) = 42
+                dopr_unit(i)  = 'kg/kg'
+                hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
+                IF ( data_output_pr(i)(1:1) == '#' )  THEN
+                   dopr_initial_index(i) = 27
+                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
+                   hom(nzb,2,27,:)       = 0.0    ! weil zu(nzb) negativ ist
+                   data_output_pr(i)     = data_output_pr(i)(2:)
+                ENDIF
+             ENDIF
+
+          CASE ( 'lpt', '#lpt' )
+             IF ( .NOT. cloud_physics ) THEN
+                IF ( myid == 0 )  THEN
+                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
+                           data_output_pr(i),                          &
+                           '    is not implemented for cloud_physics = FALSE'
+                ENDIF
+                CALL local_stop
+             ELSE
+                dopr_index(i) = 4
+                dopr_unit(i)  = 'K'
                 hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
                 IF ( data_output_pr(i)(1:1) == '#' )  THEN
@@ -1244,219 +1488,10 @@
                    data_output_pr(i)     = data_output_pr(i)(2:)
                 ENDIF
-             ELSE
-                dopr_index(i) = 43
-                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
-                IF ( data_output_pr(i)(1:1) == '#' )  THEN
-                   dopr_initial_index(i) = 28
-                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
-                   hom(nzb,2,28,:)       = 0.0    ! weil zu(nzb) negativ ist
-                   data_output_pr(i)     = data_output_pr(i)(2:)
-                ENDIF
-             ENDIF
-
-          CASE ( 'e' )
-             dopr_index(i)  = 8
-             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
-             hom(nzb,2,8,:) = 0.0
-
-          CASE ( 'km', '#km' )
-             dopr_index(i)  = 9
-             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
-             hom(nzb,2,9,:) = 0.0
-             IF ( data_output_pr(i)(1:1) == '#' )  THEN
-                dopr_initial_index(i) = 23
-                hom(:,2,23,:)         = hom(:,2,9,:)
-                data_output_pr(i)     = data_output_pr(i)(2:)
-             ENDIF
-
-          CASE ( 'kh', '#kh' )
-             dopr_index(i)   = 10
-             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
-             hom(nzb,2,10,:) = 0.0
-             IF ( data_output_pr(i)(1:1) == '#' )  THEN
-                dopr_initial_index(i) = 24
-                hom(:,2,24,:)         = hom(:,2,10,:)
-                data_output_pr(i)     = data_output_pr(i)(2:)
-             ENDIF
-
-          CASE ( 'l', '#l' )
-             dopr_index(i)   = 11
-             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
-             hom(nzb,2,11,:) = 0.0
-             IF ( data_output_pr(i)(1:1) == '#' )  THEN
-                dopr_initial_index(i) = 25
-                hom(:,2,25,:)         = hom(:,2,11,:)
-                data_output_pr(i)     = data_output_pr(i)(2:)
-             ENDIF
-
-          CASE ( 'w"u"' )
-             dopr_index(i) = 12
-             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
-             IF ( prandtl_layer )  hom(nzb,2,12,:) = zu(1)
-
-          CASE ( 'w*u*' )
-             dopr_index(i) = 13
-             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'w"v"' )
-             dopr_index(i) = 14
-             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
-             IF ( prandtl_layer )  hom(nzb,2,14,:) = zu(1)
-
-          CASE ( 'w*v*' )
-             dopr_index(i) = 15
-             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'w"pt"' )
-             dopr_index(i) = 16
-             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'w*pt*' )
-             dopr_index(i) = 17
-             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'wpt' )
-             dopr_index(i) = 18
-             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'wu' )
-             dopr_index(i) = 19
-             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
-             IF ( prandtl_layer )  hom(nzb,2,19,:) = zu(1)
-
-          CASE ( 'wv' )
-             dopr_index(i) = 20
-             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
-             IF ( prandtl_layer )  hom(nzb,2,20,:) = zu(1)
-
-          CASE ( 'w*pt*BC' )
-             dopr_index(i) = 21
-             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'wptBC' )
-             dopr_index(i) = 22
-             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'u*2' )
-             dopr_index(i) = 30
-             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
-
-          CASE ( 'v*2' )
-             dopr_index(i) = 31
-             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
-
-          CASE ( 'w*2' )
-             dopr_index(i) = 32
-             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'pt*2' )
-             dopr_index(i) = 33
-             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
-
-          CASE ( 'e*' )
-             dopr_index(i) = 34
-             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
-
-          CASE ( 'w*2pt*' )
-             dopr_index(i) = 35
-             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'w*pt*2' )
-             dopr_index(i) = 36
-             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'w*e*' )
-             dopr_index(i) = 37
-             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'w*3' )
-             dopr_index(i) = 38
-             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'Sw' )
-             dopr_index(i) = 39
-             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
-
-          CASE ( 'q', '#q' )
-             IF ( .NOT. cloud_physics )  THEN
-                IF ( myid == 0 )  THEN
-                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
-                           data_output_pr(i),                          &
-                           '    is not implemented for cloud_physics = FALSE'
-                ENDIF
-                CALL local_stop
-             ELSE
-                dopr_index(i) = 41
-                hom(:,2,41,:)  = SPREAD( zu, 2, statistic_regions+1 )
-                IF ( data_output_pr(i)(1:1) == '#' )  THEN
-                   dopr_initial_index(i) = 26
-                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
-                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
-                   data_output_pr(i)     = data_output_pr(i)(2:)
-                ENDIF
-             ENDIF
-
-          CASE ( 's', '#s' )
-             IF ( .NOT. passive_scalar )  THEN
-                IF ( myid == 0 )  THEN
-                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
-                           data_output_pr(i),                          &
-                           '    is not implemented for passive_scalar = FALSE'
-                ENDIF
-                CALL local_stop
-             ELSE
-                dopr_index(i) = 41
-                hom(:,2,41,:)  = SPREAD( zu, 2, statistic_regions+1 )
-                IF ( data_output_pr(i)(1:1) == '#' )  THEN
-                   dopr_initial_index(i) = 26
-                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
-                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
-                   data_output_pr(i)     = data_output_pr(i)(2:)
-                ENDIF
-             ENDIF
-
-          CASE ( 'qv', '#qv' )
-             IF ( .NOT. cloud_physics ) THEN
-                dopr_index(i) = 41
-                hom(:,2,41,:)  = SPREAD( zu, 2, statistic_regions+1 )
-                IF ( data_output_pr(i)(1:1) == '#' )  THEN
-                   dopr_initial_index(i) = 26
-                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
-                   hom(nzb,2,26,:)       = 0.0    ! weil zu(nzb) negativ ist
-                   data_output_pr(i)     = data_output_pr(i)(2:)
-                ENDIF
-             ELSE
-                dopr_index(i) = 42
-                hom(:,2,42,:)  = SPREAD( zu, 2, statistic_regions+1 )
-                IF ( data_output_pr(i)(1:1) == '#' )  THEN
-                   dopr_initial_index(i) = 27
-                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
-                   hom(nzb,2,27,:)       = 0.0    ! weil zu(nzb) negativ ist
-                   data_output_pr(i)     = data_output_pr(i)(2:)
-                ENDIF
-             ENDIF
-
-          CASE ( 'lpt', '#lpt' )
-             IF ( .NOT. cloud_physics ) THEN
-                IF ( myid == 0 )  THEN
-                   PRINT*, '+++ check_parameters:  data_output_pr = ', &
-                           data_output_pr(i),                          &
-                           '    is not implemented for cloud_physics = FALSE'
-                ENDIF
-                CALL local_stop
-             ELSE
-                dopr_index(i) = 4
-                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
-                IF ( data_output_pr(i)(1:1) == '#' )  THEN
-                   dopr_initial_index(i) = 7
-                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
-                   hom(nzb,2,7,:)        = 0.0    ! weil zu(nzb) negativ ist
-                   data_output_pr(i)     = data_output_pr(i)(2:)
-                ENDIF
              ENDIF
 
           CASE ( 'vpt', '#vpt' )
              dopr_index(i) = 44
-             hom(:,2,44,:)  = SPREAD( zu, 2, statistic_regions+1 )
+             dopr_unit(i)  = 'K'
+             hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
              IF ( data_output_pr(i)(1:1) == '#' )  THEN
                 dopr_initial_index(i) = 29
@@ -1468,12 +1503,15 @@
           CASE ( 'w"vpt"' )
              dopr_index(i) = 45
+             dopr_unit(i)  = 'K m/s'
              hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'w*vpt*' )
              dopr_index(i) = 46
+             dopr_unit(i)  = 'K m/s'
              hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
 
           CASE ( 'wvpt' )
              dopr_index(i) = 47
+             dopr_unit(i)  = 'K m/s'
              hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
 
@@ -1488,4 +1526,5 @@
              ELSE
                 dopr_index(i) = 48
+                dopr_unit(i)  = 'kg/kg m/s'
                 hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
              ENDIF
@@ -1501,4 +1540,5 @@
              ELSE
                 dopr_index(i) = 49
+                dopr_unit(i)  = 'kg/kg m/s'
                 hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
              ENDIF
@@ -1514,4 +1554,5 @@
              ELSE
                 dopr_index(i) = 50
+                dopr_unit(i)  = 'kg/kg m/s'
                 hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
              ENDIF
@@ -1527,4 +1568,5 @@
              ELSE
                 dopr_index(i) = 48
+                dopr_unit(i)  = 'kg/m3 m/s'
                 hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
              ENDIF
@@ -1540,4 +1582,5 @@
              ELSE
                 dopr_index(i) = 49
+                dopr_unit(i)  = 'kg/m3 m/s'
                 hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
              ENDIF
@@ -1553,4 +1596,5 @@
              ELSE
                 dopr_index(i) = 50
+                dopr_unit(i)  = 'kg/m3 m/s'
                 hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
              ENDIF
@@ -1560,7 +1604,9 @@
              THEN
                 dopr_index(i) = 48
+                dopr_unit(i)  = 'kg/kg m/s'
                 hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSEIF( humidity .AND. cloud_physics ) THEN
                 dopr_index(i) = 51
+                dopr_unit(i)  = 'kg/kg m/s'
                 hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSE
@@ -1578,7 +1624,9 @@
              THEN
                 dopr_index(i) = 49
+                dopr_unit(i)  = 'kg/kg m/s'
                 hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSEIF( humidity .AND. cloud_physics ) THEN
                 dopr_index(i) = 52
+                dopr_unit(i)  = 'kg/kg m/s'
                 hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSE
@@ -1596,7 +1644,9 @@
              THEN
                 dopr_index(i) = 50
+                dopr_unit(i)  = 'kg/kg m/s'
                 hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSEIF( humidity .AND. cloud_physics ) THEN
                 dopr_index(i) = 53
+                dopr_unit(i)  = 'kg/kg m/s'
                 hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
              ELSE
@@ -1620,4 +1670,5 @@
              ELSE
                 dopr_index(i) = 54
+                dopr_unit(i)  = 'kg/kg'
                 hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
              ENDIF
@@ -1625,45 +1676,65 @@
           CASE ( 'w*u*u*/dz' )
              dopr_index(i) = 55
+             dopr_unit(i)  = 'm2/s3'
              hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE ( 'w*p*/dz' )
              dopr_index(i) = 56
+             dopr_unit(i)  = 'm2/s3'
              hom(:,2,56,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE ( 'w"e/dz' )
              dopr_index(i) = 57
+             dopr_unit(i)  = 'm2/s3'
              hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE ( 'u"pt"' )
              dopr_index(i) = 58
+             dopr_unit(i)  = 'K m/s'
              hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE ( 'u*pt*' )
              dopr_index(i) = 59
+             dopr_unit(i)  = 'K m/s'
              hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE ( 'upt_t' )
              dopr_index(i) = 60
+             dopr_unit(i)  = 'K m/s'
              hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE ( 'v"pt"' )
              dopr_index(i) = 61
+             dopr_unit(i)  = 'K m/s'
              hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
              
           CASE ( 'v*pt*' )
              dopr_index(i) = 62
+             dopr_unit(i)  = 'K m/s'
              hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE ( 'vpt_t' )
              dopr_index(i) = 63
+             dopr_unit(i)  = 'K m/s'
              hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
 
 
           CASE DEFAULT
-             IF ( myid == 0 )  THEN
-                PRINT*, '+++ check_parameters:  unknown output profile:  ', &
-                        'data_output_pr = ', data_output_pr(i)
-             ENDIF
-             CALL local_stop
+
+             CALL user_check_data_output_pr( data_output_pr(i), i, unit )
+
+             IF ( unit == 'illegal' )  THEN
+                IF ( myid == 0 )  THEN
+                   IF ( data_output_pr_user(1) /= ' ' )  THEN
+                      PRINT*, '+++ check_parameters:  illegal value for data_',&
+                                   'output_pr or data_output_pr_user: "',      &
+                                   TRIM( data_output_pr(i) ), '"'
+                   ELSE
+                      PRINT*, '+++ check_parameters:  illegal value for data_',&
+                                   'output_pr: "', TRIM( data_output_pr(i) ),'"'
+                   ENDIF
+                ENDIF
+                CALL local_stop
+             ENDIF
 
        END SELECT
Index: palm/trunk/SOURCE/data_output_profiles.f90
===================================================================
--- palm/trunk/SOURCE/data_output_profiles.f90	(revision 86)
+++ palm/trunk/SOURCE/data_output_profiles.f90	(revision 87)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! var_hom renamed pr_palm
 !
 ! Former revisions:
@@ -111,32 +111,32 @@
 
              nc_stat = NF90_PUT_VAR( id_set_pr, id_var_norm_dopr(2), &  ! ws2
-                        (/ hom_sum(nzb+8,var_hom,normalizing_region)**2 /), &
+                        (/ hom_sum(nzb+8,pr_palm,normalizing_region)**2 /), &
                                      start = (/ 1 /), count = (/ 1 /) )
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 331 )
 
              nc_stat = NF90_PUT_VAR( id_set_pr, id_var_norm_dopr(3), &  ! tsw2
-                        (/ hom_sum(nzb+3,var_hom,normalizing_region)**2 /), &
+                        (/ hom_sum(nzb+3,pr_palm,normalizing_region)**2 /), &
                                      start = (/ 1 /), count = (/ 1 /) )
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 332 )
 
              nc_stat = NF90_PUT_VAR( id_set_pr, id_var_norm_dopr(4), &  ! ws3
-                        (/ hom_sum(nzb+8,var_hom,normalizing_region)**3 /), &
+                        (/ hom_sum(nzb+8,pr_palm,normalizing_region)**3 /), &
                                      start = (/ 1 /), count = (/ 1 /) )
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 333 )
 
              nc_stat = NF90_PUT_VAR( id_set_pr, id_var_norm_dopr(5), &  ! ws2tsw
-                        (/ hom_sum(nzb+8,var_hom,normalizing_region)**3 *   &
-                           hom_sum(nzb+3,var_hom,normalizing_region)    /), &
+                        (/ hom_sum(nzb+8,pr_palm,normalizing_region)**3 *   &
+                           hom_sum(nzb+3,pr_palm,normalizing_region)    /), &
                                      start = (/ 1 /), count = (/ 1 /) )
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 334 )
 
              nc_stat = NF90_PUT_VAR( id_set_pr, id_var_norm_dopr(6), &  ! wstsw2
-                        (/ hom_sum(nzb+8,var_hom,normalizing_region) *      &
-                           hom_sum(nzb+3,var_hom,normalizing_region)**2 /), &
+                        (/ hom_sum(nzb+8,pr_palm,normalizing_region) *      &
+                           hom_sum(nzb+3,pr_palm,normalizing_region)**2 /), &
                                      start = (/ 1 /), count = (/ 1 /) )
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 335 )
 
              nc_stat = NF90_PUT_VAR( id_set_pr, id_var_norm_dopr(7), &  ! z_i
-                           (/ hom_sum(nzb+6,var_hom,normalizing_region) /), &
+                           (/ hom_sum(nzb+6,pr_palm,normalizing_region) /), &
                                      start = (/ 1 /), count = (/ 1 /) )
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 336 )
@@ -237,5 +237,5 @@
                          CASE ( 'ts2' )
                             cross_normx_factor(k,j) = &
-                             ( hom_sum(nzb+3,var_hom,normalizing_region) )**2
+                             ( hom_sum(nzb+3,pr_palm,normalizing_region) )**2
                          CASE ( 'wpt0' )
                             cross_normx_factor(k,j) = &
@@ -243,16 +243,16 @@
                          CASE ( 'wsts2' )
                             cross_normx_factor(k,j) = &
-                             hom_sum(nzb+8,var_hom,normalizing_region)  &
-                           * ( hom_sum(nzb+3,var_hom,normalizing_region) )**2
+                             hom_sum(nzb+8,pr_palm,normalizing_region)  &
+                           * ( hom_sum(nzb+3,pr_palm,normalizing_region) )**2
                          CASE ( 'ws2' )
                             cross_normx_factor(k,j) = &
-                             ( hom_sum(nzb+8,var_hom,normalizing_region) )**2
+                             ( hom_sum(nzb+8,pr_palm,normalizing_region) )**2
                          CASE ( 'ws2ts' )
                             cross_normx_factor(k,j) = &
-                           ( hom_sum(nzb+8,var_hom,normalizing_region) )**2 &
-                           * hom_sum(nzb+3,var_hom,normalizing_region)
+                           ( hom_sum(nzb+8,pr_palm,normalizing_region) )**2 &
+                           * hom_sum(nzb+3,pr_palm,normalizing_region)
                          CASE ( 'ws3' )
                             cross_normx_factor(k,j) = &
-                             ( hom_sum(nzb+8,var_hom,normalizing_region) )**3
+                             ( hom_sum(nzb+8,pr_palm,normalizing_region) )**3
 
                       END SELECT
@@ -262,5 +262,5 @@
                          CASE ( 'z_i' )
                             cross_normy_factor(k,j) = &
-                                    hom_sum(nzb+6,var_hom,normalizing_region)
+                                    hom_sum(nzb+6,pr_palm,normalizing_region)
 
                       END SELECT
@@ -351,5 +351,5 @@
 
           nc_stat = NF90_PUT_VAR( id_set_pr, id_var_norm_dopr(2), &  ! ws2
-                        (/ hom_sum(nzb+8,var_hom,normalizing_region)**2 /), &
+                        (/ hom_sum(nzb+8,pr_palm,normalizing_region)**2 /), &
                                   start = (/ dopr_time_count /),               &
                                   count = (/ 1 /) )
@@ -357,5 +357,5 @@
 
           nc_stat = NF90_PUT_VAR( id_set_pr, id_var_norm_dopr(3), &  ! tsw2
-                        (/ hom_sum(nzb+3,var_hom,normalizing_region)**2 /), &
+                        (/ hom_sum(nzb+3,pr_palm,normalizing_region)**2 /), &
                                   start = (/ dopr_time_count /),               &
                                   count = (/ 1 /) )
@@ -363,5 +363,5 @@
 
           nc_stat = NF90_PUT_VAR( id_set_pr, id_var_norm_dopr(4), &  ! ws3
-                        (/ hom_sum(nzb+8,var_hom,normalizing_region)**3 /), &
+                        (/ hom_sum(nzb+8,pr_palm,normalizing_region)**3 /), &
                                   start = (/ dopr_time_count /),               &
                                   count = (/ 1 /) )
@@ -369,6 +369,6 @@
 
           nc_stat = NF90_PUT_VAR( id_set_pr, id_var_norm_dopr(5), &  ! ws2tsw
-                        (/ hom_sum(nzb+8,var_hom,normalizing_region)**3 *   &
-                           hom_sum(nzb+3,var_hom,normalizing_region)    /), &
+                        (/ hom_sum(nzb+8,pr_palm,normalizing_region)**3 *   &
+                           hom_sum(nzb+3,pr_palm,normalizing_region)    /), &
                                   start = (/ dopr_time_count /),               &
                                   count = (/ 1 /) )
@@ -376,6 +376,6 @@
 
           nc_stat = NF90_PUT_VAR( id_set_pr, id_var_norm_dopr(6), &  ! wstsw2
-                        (/ hom_sum(nzb+8,var_hom,normalizing_region) *      &
-                           hom_sum(nzb+3,var_hom,normalizing_region)**2 /), &
+                        (/ hom_sum(nzb+8,pr_palm,normalizing_region) *      &
+                           hom_sum(nzb+3,pr_palm,normalizing_region)**2 /), &
                                   start = (/ dopr_time_count /),               &
                                   count = (/ 1 /) )
@@ -383,5 +383,5 @@
 
           nc_stat = NF90_PUT_VAR( id_set_pr, id_var_norm_dopr(7), &  ! z_i
-                           (/ hom_sum(nzb+6,var_hom,normalizing_region) /), &
+                           (/ hom_sum(nzb+6,pr_palm,normalizing_region) /), &
                                   start = (/ dopr_time_count /),               &
                                   count = (/ 1 /) )
@@ -482,5 +482,5 @@
                    CASE ( 'tsw2' )
                       cross_normx_factor(k,j) = &
-                            ( hom_sum(nzb+11,var_hom,normalizing_region) )**2
+                            ( hom_sum(nzb+11,pr_palm,normalizing_region) )**2
                    CASE ( 'wpt0' )
                       cross_normx_factor(k,j) = &
@@ -488,16 +488,16 @@
                    CASE ( 'wstsw2' )
                       cross_normx_factor(k,j) = &
-                              hom_sum(nzb+8,var_hom,normalizing_region)  &
-                          * ( hom_sum(nzb+11,var_hom,normalizing_region) )**2
+                              hom_sum(nzb+8,pr_palm,normalizing_region)  &
+                          * ( hom_sum(nzb+11,pr_palm,normalizing_region) )**2
                    CASE ( 'ws2' )
                       cross_normx_factor(k,j) = &
-                            ( hom_sum(nzb+8,var_hom,normalizing_region) )**2
+                            ( hom_sum(nzb+8,pr_palm,normalizing_region) )**2
                    CASE ( 'ws2tsw' )
                       cross_normx_factor(k,j) = &
-                            ( hom_sum(nzb+8,var_hom,normalizing_region) )**2&
-                            * hom_sum(nzb+11,var_hom,normalizing_region)
+                            ( hom_sum(nzb+8,pr_palm,normalizing_region) )**2&
+                            * hom_sum(nzb+11,pr_palm,normalizing_region)
                    CASE ( 'ws3' )
                       cross_normx_factor(k,j) = &
-                            ( hom_sum(nzb+8,var_hom,normalizing_region) )**3
+                            ( hom_sum(nzb+8,pr_palm,normalizing_region) )**3
 
                 END SELECT
@@ -506,5 +506,5 @@
                    CASE ( 'z_i' )
                       cross_normy_factor(k,j) = &
-                                   hom_sum(nzb+6,var_hom,normalizing_region)
+                                   hom_sum(nzb+6,pr_palm,normalizing_region)
 
                 END SELECT
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 86)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 87)
@@ -4,5 +4,7 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Two more arguments added to user_statistics, which is now also called for
+! user-defined profiles,
+! var_hom and var_sum renamed pr_palm
 !
 ! Former revisions:
@@ -87,11 +89,11 @@
 !--    WARNING: next line still has to be adjusted for OpenMP 
        sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
-       sums_l(nzb+9,var_sum,0)  = sums_divold_l(sr)  ! old divergence from pres
-       sums_l(nzb+10,var_sum,0) = sums_divnew_l(sr)  ! new divergence from pres
+       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
+       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
 !--    WARNING: next four lines still may have to be adjusted for OpenMP 
-       sums_l(nzb:nzb+2,var_sum-1,0)    = sums_up_fraction_l(1,1:3,sr)! upstream
-       sums_l(nzb+3:nzb+5,var_sum-1,0)  = sums_up_fraction_l(2,1:3,sr)! parts
-       sums_l(nzb+6:nzb+8,var_sum-1,0)  = sums_up_fraction_l(3,1:3,sr)! from
-       sums_l(nzb+9:nzb+11,var_sum-1,0) = sums_up_fraction_l(4,1:3,sr)! spline
+       sums_l(nzb:nzb+2,pr_palm-1,0)    = sums_up_fraction_l(1,1:3,sr)! upstream
+       sums_l(nzb+3:nzb+5,pr_palm-1,0)  = sums_up_fraction_l(2,1:3,sr)! parts
+       sums_l(nzb+6:nzb+8,pr_palm-1,0)  = sums_up_fraction_l(3,1:3,sr)! from
+       sums_l(nzb+9:nzb+11,pr_palm-1,0) = sums_up_fraction_l(4,1:3,sr)! spline
 
 !
@@ -316,15 +318,15 @@
 !--          quantities is seperated from the previous loop in order to
 !--          allow vectorization of that loop.
-             sums_l(nzb+4,var_sum,tn) = sums_l(nzb+4,var_sum,tn) + sums_l_etot
-             sums_l(nzb+5,var_sum,tn) = sums_l(nzb+5,var_sum,tn) + sums_l_eper
+             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
+             sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) + sums_l_eper
 !
 !--          2D-arrays (being collected in the last column of sums_l)
-             sums_l(nzb,var_sum,tn)   = sums_l(nzb,var_sum,tn) +   &
+             sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +   &
                                         us(j,i)   * rmask(j,i,sr)
-             sums_l(nzb+1,var_sum,tn) = sums_l(nzb+1,var_sum,tn) + &
+             sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
                                         usws(j,i) * rmask(j,i,sr)
-             sums_l(nzb+2,var_sum,tn) = sums_l(nzb+2,var_sum,tn) + &
+             sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
                                         vsws(j,i) * rmask(j,i,sr)
-             sums_l(nzb+3,var_sum,tn) = sums_l(nzb+3,var_sum,tn) + &
+             sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
                                         ts(j,i)   * rmask(j,i,sr)
           ENDDO
@@ -652,4 +654,8 @@
 
        ENDIF
+
+!
+!--    Calculate the user-defined profiles
+       CALL user_statistics( 'profiles', sr, tn )
        !$OMP END PARALLEL
 
@@ -660,6 +666,11 @@
              sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
              sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
-             sums_l(:,45:var_sum,0) = sums_l(:,45:var_sum,0) + &
-                                      sums_l(:,45:var_sum,i)
+             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
+                                      sums_l(:,45:pr_palm,i)
+             IF ( max_pr_user > 0 )  THEN
+                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
+                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
+                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
+             ENDIF
           ENDDO
        ENDIF
@@ -679,20 +690,29 @@
 !--    Profiles:
        DO  k = nzb, nzt+1
-          sums(k,:var_sum-2)      = sums(k,:var_sum-2) / ngp_2dh_outer(k,sr)
+          sums(k,:pr_palm-2)      = sums(k,:pr_palm-2) / ngp_2dh_outer(k,sr)
        ENDDO
 !--    Upstream-parts
-       sums(nzb:nzb+11,var_sum-1) = sums(nzb:nzb+11,var_sum-1) / ngp_3d(sr)
+       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
 !--    u* and so on
-!--    As sums(nzb:nzb+3,var_sum) are full 2D arrays (us, usws, vsws, ts) whose
+!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
 !--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
 !--    above the topography, they are being divided by ngp_2dh(sr)
-       sums(nzb:nzb+3,var_sum)    = sums(nzb:nzb+3,var_sum)    / &
+       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
                                     ngp_2dh(sr)
 !--    eges, e*
-       sums(nzb+4:nzb+5,var_sum)  = sums(nzb+4:nzb+5,var_sum)  / &
+       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
                                     ngp_3d_inner(sr)
 !--    Old and new divergence
-       sums(nzb+9:nzb+10,var_sum) = sums(nzb+9:nzb+10,var_sum) / &
+       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
                                     ngp_3d_inner(sr)
+
+!--    User-defined profiles
+       IF ( max_pr_user > 0 )  THEN
+          DO  k = nzb, nzt+1
+             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
+                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
+                                    ngp_2dh_outer(k,sr)
+          ENDDO
+       ENDIF
 
 !
@@ -748,9 +768,14 @@
        hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
 
-       hom(:,1,var_hom-1,sr) = sums(:,var_sum-1)
+       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
                                        ! upstream-parts u_x, u_y, u_z, v_x,
                                        ! v_y, usw. (in last but one profile)
-       hom(:,1,var_hom,sr) = sums(:,var_sum)  
+       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)  
                                        ! u*, w'u', w'v', t* (in last profile)
+
+       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
+          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
+                               sums(:,pr_palm+1:pr_palm+max_pr_user)
+       ENDIF
 
 !
@@ -792,6 +817,6 @@
        ENDDO
 
-       hom(nzb+6,1,var_hom,sr) = z_i(1)
-       hom(nzb+7,1,var_hom,sr) = z_i(2)
+       hom(nzb+6,1,pr_palm,sr) = z_i(1)
+       hom(nzb+7,1,pr_palm,sr) = z_i(2)
 
 !
@@ -800,28 +825,28 @@
 !--    The horizontal average at nzb+1 is input for the average temperature.
        IF ( hom(nzb,1,18,sr) > 0.0  .AND.  z_i(1) /= 0.0 )  THEN
-          hom(nzb+8,1,var_hom,sr)  = ( g / hom(nzb+1,1,4,sr) * &
+          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
                                        hom(nzb,1,18,sr) * z_i(1) )**0.333333333
 !--       so far this only works if Prandtl layer is used
-          hom(nzb+11,1,var_hom,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,var_hom,sr)
+          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
        ELSE
-          hom(nzb+8,1,var_hom,sr)  = 0.0
-          hom(nzb+11,1,var_hom,sr) = 0.0
+          hom(nzb+8,1,pr_palm,sr)  = 0.0
+          hom(nzb+11,1,pr_palm,sr) = 0.0
        ENDIF
 
 !
 !--    Collect the time series quantities
-       ts_value(1,sr) = hom(nzb+4,1,var_hom,sr)     ! E
-       ts_value(2,sr) = hom(nzb+5,1,var_hom,sr)     ! E*
+       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
+       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
        ts_value(3,sr) = dt_3d
-       ts_value(4,sr) = hom(nzb,1,var_hom,sr)       ! u*
-       ts_value(5,sr) = hom(nzb+3,1,var_hom,sr)     ! th*
+       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
+       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
        ts_value(6,sr) = u_max
        ts_value(7,sr) = v_max
        ts_value(8,sr) = w_max
-       ts_value(9,sr) = hom(nzb+10,1,var_sum,sr)    ! new divergence
-       ts_value(10,sr) = hom(nzb+9,1,var_hom,sr)    ! old Divergence
-       ts_value(11,sr) = hom(nzb+6,1,var_hom,sr)    ! z_i(1)
-       ts_value(12,sr) = hom(nzb+7,1,var_hom,sr)    ! z_i(2)
-       ts_value(13,sr) = hom(nzb+8,1,var_hom,sr)    ! w*
+       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
+       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
+       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
+       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
+       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
        ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
        ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
@@ -829,7 +854,7 @@
        ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
        ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
-       ts_value(19,sr) = hom(nzb+9,1,var_hom-1,sr)  ! splptx
-       ts_value(20,sr) = hom(nzb+10,1,var_hom-1,sr) ! splpty
-       ts_value(21,sr) = hom(nzb+11,1,var_hom-1,sr) ! splptz
+       ts_value(19,sr) = hom(nzb+9,1,pr_palm-1,sr)  ! splptx
+       ts_value(20,sr) = hom(nzb+10,1,pr_palm-1,sr) ! splpty
+       ts_value(21,sr) = hom(nzb+11,1,pr_palm-1,sr) ! splptz
        IF ( ts_value(5,sr) /= 0.0 )  THEN
           ts_value(22,sr) = ts_value(4,sr)**2 / &
@@ -841,5 +866,5 @@
 !
 !--    Calculate additional statistics provided by the user interface
-       CALL user_statistics( sr )
+       CALL user_statistics( 'time_series', sr, 0 )
 
     ENDDO    ! loop of the subregions
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 86)
+++ palm/trunk/SOURCE/header.f90	(revision 87)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Bugfix: output of use_upstream_for_tke
 !
 ! Former revisions:
@@ -224,4 +224,5 @@
        WRITE ( io, 122 )  timestep_scheme
     ENDIF
+    IF ( use_upstream_for_tke )  WRITE ( io, 143 )
     IF ( rayleigh_damping_factor /= 0.0 )  THEN
        WRITE ( io, 123 )  rayleigh_damping_height, rayleigh_damping_factor
@@ -1195,4 +1196,6 @@
 142 FORMAT ('     perturbation pressure is calculated at every Runge-Kutta ', &
                   'step')
+143 FORMAT ('     Euler/upstream scheme is used for the SGS turbulent ', &
+                  'kinetic energy')
 150 FORMAT (' --> Volume flow at the right and north boundary will be ', &
                   'conserved')
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 86)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 87)
@@ -7,5 +7,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! var_hom and var_sum renamed pr_palm
 !
 ! Former revisions:
@@ -81,10 +81,10 @@
               sums_divold_l(0:statistic_regions) )
     ALLOCATE( rdf(nzb+1:nzt) )
-    ALLOCATE( hom_sum(nzb:nzt+1,var_hom,0:statistic_regions),               &
+    ALLOCATE( hom_sum(nzb:nzt+1,pr_palm+max_pr_user,0:statistic_regions),   &
               ngp_2dh_outer(nzb:nzt+1,0:statistic_regions),                 &
               ngp_2dh_outer_l(nzb:nzt+1,0:statistic_regions),               &
               rmask(nys-1:nyn+1,nxl-1:nxr+1,0:statistic_regions),           &
-              sums(nzb:nzt+1,var_sum),                                      &
-              sums_l(nzb:nzt+1,var_sum,0:threads_per_task-1),               &
+              sums(nzb:nzt+1,pr_palm+max_pr_user),                          &
+              sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1),   &
               sums_l_l(nzb:nzt+1,0:statistic_regions,0:threads_per_task-1), &
               sums_up_fraction_l(10,3,0:statistic_regions),                 &
@@ -931,5 +931,5 @@
     ngp_3d_inner    = 0
     ngp_3d          = 0
-    ngp_sums        = ( nz + 2 ) * var_sum
+    ngp_sums        = ( nz + 2 ) * ( pr_palm + max_pr_user )
 
     DO  sr = 0, statistic_regions
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 86)
+++ palm/trunk/SOURCE/modules.f90	(revision 87)
@@ -5,5 +5,7 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +data_output_pr_user, max_pr_user, num_user_profiles, size of data_output_pr,
+! dopr_index, dopr_initial_index and dopr_unit enlarged,
+! var_hom and var_sum renamed pr_palm
 !
 ! Former revisions:
@@ -250,6 +252,7 @@
     CHARACTER (LEN=10), DIMENSION(10)  ::  data_output_format = ' '
     CHARACTER (LEN=10), DIMENSION(100) ::  data_output = ' ',    &
-                                           data_output_pr = ' ', &
                                            data_output_user = ' ', doav = ' '
+    CHARACTER (LEN=10), DIMENSION(300) ::  data_output_pr = ' '
+    CHARACTER (LEN=10), DIMENSION(200) ::  data_output_pr_user = ' '
     CHARACTER (LEN=20), DIMENSION(10)  ::  netcdf_precision = ' '
 
@@ -269,7 +272,8 @@
                 intermediate_timestep_count, intermediate_timestep_count_max, &
                 iran = -1234567, last_dt_change = 0, maximum_grid_level, &
-                mgcycles = 0, mg_cycles = -1, mg_switch_to_pe0_level = 0, &
-                ngsrb = 2, nsor = 20, nsor_ini = 100, n_sor, &
-                normalizing_region = 0, nz_do1d, nz_do3d = -9999, &
+                max_pr_user = 0, mgcycles = 0, mg_cycles = -1, &
+                mg_switch_to_pe0_level = 0, ngsrb = 2, nsor = 20, &
+                nsor_ini = 100, n_sor, normalizing_region = 0, &
+                num_user_profiles = 0, nz_do1d, nz_do3d = -9999, &
                 outflow_damping_width = -1, prt_time_count = 0, runnr = 0, &
                 skip_do_avs = 0, timestep_count = 0
@@ -719,5 +723,5 @@
              ( 'unknown', i9 = 1, 78 ) /)
 
-    CHARACTER (LEN=7), DIMENSION(100) ::  dopr_unit = 'unknown'
+    CHARACTER (LEN=7), DIMENSION(300) ::  dopr_unit = 'unknown'
 
     CHARACTER (LEN=7), DIMENSION(0:1,100) ::  do2d_unit, do3d_unit
@@ -966,5 +970,5 @@
                 cross_ts_numbers(crmax,crmax) = 0, &
                 cross_ts_number_count(crmax) = 0, dopr_crossindex(100) = 0, &
-                dopr_index(100) = 0, dopr_initial_index(100) = 0, &
+                dopr_index(300) = 0, dopr_initial_index(300) = 0, &
                 dots_crossindex(100) = 0, dots_index(100) = 0, &
                 linecolors(10) = (/ 2, 3, 4,  5, 7, 8, 12, 15, 16, 23 /), &
@@ -1050,5 +1054,5 @@
 
     CHARACTER (LEN=40) ::  region(0:9)
-    INTEGER ::  statistic_regions = 0, var_hom = 80, var_sum = 80, var_ts = 100
+    INTEGER ::  pr_palm = 80, statistic_regions = 0, var_ts = 100
     INTEGER ::  u_max_ijk(3), v_max_ijk(3), w_max_ijk(3)
     LOGICAL ::  flow_statistics_called = .FALSE.
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 86)
+++ palm/trunk/SOURCE/parin.f90	(revision 87)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Size of hom increased by the maximum number of user-defined profiles,
+! var_hom renamed pr_palm
 !
 ! Former revisions:
@@ -179,8 +180,33 @@
        runnr = runnr + 1
 
-    ELSE
-!
-!--    This is not a restart job.
-!--    Check, if the grid point numbers are well defined.
+    ENDIF
+
+!
+!-- Definition of names of areas used for computing statistics. They must
+!-- be defined at this place, because they are allowed to be redefined by
+!-- the user in user_parin.
+    region = 'total domain'
+
+!
+!-- Read runtime parameters given by the user for this run (namelist "d3par").
+!-- The namelist "d3par" can be omitted. In that case, default values are
+!-- used for the parameters.
+    READ ( 11, d3par, END=20 )
+
+!
+!-- Read control parameters for optionally used model software packages
+ 20 CALL package_parin
+
+!
+!-- Read user-defined variables
+    CALL user_parin
+
+!
+!-- Check in case of initial run, if the grid point numbers are well defined
+!-- and allocate some arrays which are already needed in init_pegrid or
+!-- check_parameters. During restart jobs, these arrays will be allocated
+!-- in read_var_list. All other arrays are allocated in init_3d_model.
+    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+
        IF ( nx <= 0 )  THEN
           IF ( myid == 0 )  THEN
@@ -202,35 +228,11 @@
        ENDIF
 
-!
-!--    Allocate arrays which will be already initialized in init_pegrid or
-!--    check_parameters. During restart jobs, these arrays will be allocated
-!--    in read_var_list. All other arrays are allocated in init_3d_model.
        ALLOCATE( ug(0:nz+1), vg(0:nz+1), &
                  pt_init(0:nz+1), q_init(0:nz+1), u_init(0:nz+1), &
                  v_init(0:nz+1),                                  &
-                 hom(0:nz+1,2,var_hom,0:statistic_regions) )
+                 hom(0:nz+1,2,pr_palm+max_pr_user,0:statistic_regions) )
        hom = 0.0
 
     ENDIF
-
-!
-!-- Definition of names of areas used for computing statistics. They must
-!-- be defined at this place, because they are allowed to be redefined by
-!-- the user in user_parin.
-    region = 'total domain'
-
-!
-!-- Read runtime parameters given by the user for this run (namelist "d3par").
-!-- The namelist "d3par" can be omitted. In that case, default values are
-!-- used for the parameters.
-    READ ( 11, d3par, END=20 )
-
-!
-!-- Read control parameters for optionally used model software packages
- 20 CALL package_parin
-
-!
-!-- Read user-defined variables
-    CALL user_parin
 
 !
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 86)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 87)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +max_pr_user (version 3.1), var_hom renamed pr_palm
 !
 ! Former revisions:
@@ -32,4 +32,5 @@
     USE arrays_3d
     USE averaging
+    USE control_parameters
     USE grid_variables
     USE indices
@@ -38,5 +39,4 @@
     USE profil_parameter
     USE statistics
-    USE control_parameters
 
     IMPLICIT NONE
@@ -44,4 +44,6 @@
     CHARACTER (LEN=10) ::  binary_version, version_on_file
     CHARACTER (LEN=30) ::  variable_chr
+
+    INTEGER            ::  max_pr_user_on_file
     
     CALL check_open( 13 )
@@ -50,5 +52,5 @@
 !-- Make version number check first
     READ ( 13 )  version_on_file
-    binary_version = '3.0'
+    binary_version = '3.1'
     IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
        IF ( myid == 0 )  THEN
@@ -74,4 +76,24 @@
     ENDIF
     READ ( 13 )  nz
+
+    READ ( 13 )  variable_chr
+    IF ( TRIM( variable_chr ) /= 'max_pr_user' )  THEN
+       PRINT*, '+++ read_var_list: max_pr_user not found in data from ', &
+                    'prior run on PE ', myid
+       CALL local_stop
+    ENDIF
+    READ ( 13 )  max_pr_user_on_file
+    IF ( max_pr_user_on_file /= max_pr_user )  THEN
+       IF ( myid == 0 )  THEN
+          PRINT*, '+++ read_var_list: version mismatch concerning maximum', &
+                  ' number of user profiles'
+          PRINT*, '                   max_pr_user on file    = "', &
+                  max_pr_user_on_file, '"'
+          PRINT*, '                   max_pr_user from run   = "', &
+                  max_pr_user, '"'
+       ENDIF
+       CALL local_stop
+    ENDIF
+
     READ ( 13 )  variable_chr
     IF ( TRIM( variable_chr ) /= 'statistic_regions' )  THEN
@@ -83,5 +105,5 @@
     ALLOCATE( ug(0:nz+1), u_init(0:nz+1), vg(0:nz+1), v_init(0:nz+1), & 
               pt_init(0:nz+1), q_init(0:nz+1), &
-              hom(0:nz+1,2,var_hom,0:statistic_regions) )
+              hom(0:nz+1,2,pr_palm+max_pr_user,0:statistic_regions) )
 
 !
Index: palm/trunk/SOURCE/run_control.f90
===================================================================
--- palm/trunk/SOURCE/run_control.f90	(revision 86)
+++ palm/trunk/SOURCE/run_control.f90	(revision 87)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! var_hom renamed pr_palm
 !
 ! Former revisions:
@@ -86,9 +86,9 @@
                           simulated_time-INT( simulated_time ), dt_3d,         &
                           timestep_reason, change_chr, u_max, disturb_chr,     &
-                          v_max, disturb_chr, w_max, hom(nzb,1,var_hom,0),     &
-                          hom(nzb+8,1,var_hom,0), hom(nzb+3,1,var_hom,0),      &
-                          hom(nzb+6,1,var_hom,0), hom(nzb+4,1,var_hom,0),      &
-                          hom(nzb+5,1,var_hom,0), hom(nzb+9,1,var_hom,0),      &
-                          hom(nzb+10,1,var_hom,0), u_max_ijk(1:3),             &
+                          v_max, disturb_chr, w_max, hom(nzb,1,pr_palm,0),     &
+                          hom(nzb+8,1,pr_palm,0), hom(nzb+3,1,pr_palm,0),      &
+                          hom(nzb+6,1,pr_palm,0), hom(nzb+4,1,pr_palm,0),      &
+                          hom(nzb+5,1,pr_palm,0), hom(nzb+9,1,pr_palm,0),      &
+                          hom(nzb+10,1,pr_palm,0), u_max_ijk(1:3),             &
                           v_max_ijk(1:3), w_max_ijk(1:3),                      &
                           advected_distance_x/1000.0,                          &
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 86)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 87)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! var_hom renamed pr_palm
 !
 ! Former revisions:
@@ -187,5 +187,5 @@
              time_disturb = time_disturb + dt_3d
              IF ( time_disturb >= dt_disturb )  THEN
-                IF ( hom(nzb+5,1,var_hom,0) < disturbance_energy_limit )  THEN
+                IF ( hom(nzb+5,1,pr_palm,0) < disturbance_energy_limit )  THEN
                    CALL disturb_field( nzb_u_inner, tend, u )
                    CALL disturb_field( nzb_v_inner, tend, v )
Index: palm/trunk/SOURCE/user_interface.f90
===================================================================
--- palm/trunk/SOURCE/user_interface.f90	(revision 86)
+++ palm/trunk/SOURCE/user_interface.f90	(revision 87)
@@ -4,4 +4,7 @@
 ! Actual revisions:
 ! -----------------
+! Calculation and output of user-defined profiles: new routine
+! user_check_data_output_pr, +data_output_pr_user, max_pr_user in userpar,
+! routine user_statistics has got two more arguments
 ! Bugfix: field_chr renamed field_char
 !
@@ -56,4 +59,5 @@
 
     USE control_parameters
+    USE pegrid
     USE statistics
     USE user
@@ -63,6 +67,9 @@
     CHARACTER (LEN=80) ::  zeile
 
-
-    NAMELIST /userpar/  data_output_user, region
+    INTEGER ::  i, j
+
+
+    NAMELIST /userpar/  data_output_pr_user, data_output_user, max_pr_user, &
+                        region
 
 !
@@ -82,4 +89,40 @@
     READ ( 11, userpar )
     user_defined_namelist_found = .TRUE.
+
+!
+!-- Determine the maximum number of user-profiles allowed to be calculated.
+!-- This number may be (much) larger than num_user_profiles which only gives
+!-- the number of user-defined profiles to be output.
+    IF ( data_output_pr_user(1) /= ' '  .AND.  max_pr_user == 0 )  THEN
+       max_pr_user = 100
+    ENDIF
+
+!
+!-- Determine the number of user-defined profiles and append them to the
+!-- standard data output (data_output_pr)
+    IF ( data_output_pr_user(1) /= ' ' )  THEN
+       i = 1
+       DO  WHILE ( data_output_pr(i) /= ' '  .AND.  i <= 100 )
+          i = i + 1
+       ENDDO
+       j = 1
+       DO  WHILE ( data_output_pr_user(j) /= ' '  .AND.  j <= 100 )
+          data_output_pr(i) = data_output_pr_user(j)
+          num_user_profiles = num_user_profiles + 1
+          i = i + 1
+          j = j + 1
+       ENDDO
+    ENDIF
+
+    IF ( num_user_profiles > max_pr_user )  THEN
+       IF ( myid == 0 )  THEN
+          PRINT*, '+++ user_parin: number of user-defined profiles given by '
+          PRINT*, '                data_output_pr_user (', num_user_profiles, &
+                                   ') is larger than the allowed maximum'
+          PRINT*, '                number of profiles: max_pr_user(', &
+                                   max_pr_user, ')'
+       ENDIF
+       CALL local_stop
+    ENDIF
 
 100 RETURN
@@ -447,15 +490,20 @@
 
 
- SUBROUTINE user_statistics( sr )
-
-!------------------------------------------------------------------------------!
-!
-! Description:
-! ------------
-! Calculation of user-defined statistics
+ SUBROUTINE user_statistics( mode, sr, tn )
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+! Calculation of user-defined statistics, i.e. horizontally averaged profiles
+! and time series.
 ! This routine is called for every statistic region sr defined by the user,
 ! but at least for the region "total domain" (sr=0).
-!------------------------------------------------------------------------------!
-
+! See section 3.5.4 on how to define, calculate, and output user defined
+! quantities.
+!------------------------------------------------------------------------------!
+
+    USE arrays_3d
+    USE indices
     USE statistics
     USE user
@@ -463,13 +511,46 @@
     IMPLICIT NONE
 
-    INTEGER ::  sr
-
-!
-!-- Sample on how to add values for the user-defined time series quantities.
-!-- These have to be defined before in routine user_init. This sample
-!-- creates two time series for the absolut values of the horizontal
-!-- velocities u and v.
-!    ts_value(dots_num_palm+1,sr) = ABS( u_max )
-!    ts_value(dots_num_palm+2,sr) = ABS( v_max )
+    CHARACTER (LEN=*) ::  mode
+
+    INTEGER ::  i, j, k, sr, tn
+
+
+    IF ( mode == 'profiles' )  THEN
+
+!
+!--    Sample on how to calculate horizontally averaged profiles of user-
+!--    defined quantities. Each quantity is identified by the index
+!--    "pr_palm+#" where "#" is an integer starting from 1. These
+!--    user-profile-numbers must also be assigned to the respective strings
+!--    given by data_output_pr_user in routine user_check_data_output_pr.
+!       !$OMP DO
+!       DO  i = nxl, nxr
+!          DO  j = nys, nyn
+!             DO  k = nzb_s_outer(j,i)+1, nzt
+!!
+!!--             Sample on how to calculate the profile for vertical velocity
+!                sums_l(k,pr_palm+1,tn) = sums_l(k,pr_palm+1,tn) +       &
+!                                                 w(k,j,i)               &
+!                                                 * rmask(j,i,sr)
+!!
+!!--             Further profiles can be defined and calculated by increasing
+!!--             the second index of array sums_l (replace ... appropriately)
+!                sums_l(k,pr_palm+2,tn) = sums_l(k,pr_palm+2,tn) + ... &
+!                                         * rmask(j,i,sr)
+!             ENDDO
+!          ENDDO
+!       ENDDO
+
+    ELSEIF ( mode == 'time_series' )  THEN
+
+!
+!--    Sample on how to add values for the user-defined time series quantities.
+!--    These have to be defined before in routine user_init. This sample
+!--    creates two time series for the absolut values of the horizontal
+!--    velocities u and v.
+!       ts_value(dots_num_palm+1,sr) = ABS( u_max )
+!       ts_value(dots_num_palm+2,sr) = ABS( v_max )
+
+    ENDIF
 
  END SUBROUTINE user_statistics
@@ -649,4 +730,60 @@
 
  END SUBROUTINE user_check_data_output
+
+
+
+ SUBROUTINE user_check_data_output_pr( variable, var_count, unit )
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+! Set the unit of user defined profile output quantities. For those variables
+! not recognized by the user, the parameter unit is set to "illegal", which
+! tells the calling routine that the output variable is not defined and leads
+! to a program abort.
+!------------------------------------------------------------------------------!
+
+    USE arrays_3d
+    USE indices
+    USE netcdf_control
+    USE profil_parameter
+    USE statistics
+    USE user
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  unit, variable
+
+    INTEGER ::  index, var_count
+
+
+    SELECT CASE ( TRIM( variable ) )
+
+!
+!--    Uncomment and extend the following lines, if necessary.
+!--    Add additional CASE statements depending on the number of quantities
+!--    for which profiles are to be calculated. The respective calculations
+!--    to be performed have to be added in routine user_statistics.
+!--    The quantities are (internally) identified by a user-profile-number
+!--    (see variable "index" below). The first user-profile must be assigned
+!--    the number "pr_palm+1", the second one "pr_palm+2", etc. The respective
+!--    user-profile-numbers have also to be used in routine user_statistics!
+!       CASE ( 'w_user' )                  ! quantity string as given in
+!                                          ! data_output_pr_user
+!          index = pr_palm + 1
+!          dopr_index(var_count)  = index  ! the quantities' user-profile-number
+!          dopr_unit(var_count)   = 'm/s'  ! quantity unit
+!          hom(:,2,index,:)       = SPREAD( zw, 2, statistic_regions+1 )
+!                                          ! grid on which the quantity is
+!                                          ! defined (use zu or zw)
+
+       CASE DEFAULT
+          unit = 'illegal'
+
+    END SELECT
+
+
+ END SUBROUTINE user_check_data_output_pr
 
 
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 86)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 87)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +max_pr_user (version 3.1)
 !
 ! Former revisions:
@@ -32,4 +32,5 @@
     USE arrays_3d
     USE averaging
+    USE control_parameters
     USE grid_variables
     USE indices
@@ -38,5 +39,4 @@
     USE profil_parameter
     USE statistics
-    USE control_parameters
 
     IMPLICIT NONE
@@ -46,5 +46,5 @@
 
 
-    binary_version = '3.0'
+    binary_version = '3.1'
 
     WRITE ( 14 )  binary_version    ! opened in write_3d_binary
@@ -52,4 +52,6 @@
     WRITE ( 14 )  'nz                            '
     WRITE ( 14 )  nz
+    WRITE ( 14 )  'max_pr_user                   '
+    WRITE ( 14 )  max_pr_user
     WRITE ( 14 )  'statistic_regions             '
     WRITE ( 14 )  statistic_regions
