Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 939)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 940)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! checks for parameter neutral
 !
 ! Former revisions:
@@ -1039,53 +1039,57 @@
 !
 !--    Compute initial temperature profile using the given temperature gradients
-       i = 1
-       gradient = 0.0
-
-       IF ( .NOT. ocean )  THEN
-
-          pt_vertical_gradient_level_ind(1) = 0
-          DO  k = 1, nzt+1
-             IF ( i < 11 ) THEN
-                IF ( pt_vertical_gradient_level(i) < zu(k)  .AND. &
-                     pt_vertical_gradient_level(i) >= 0.0 )  THEN
-                   gradient = pt_vertical_gradient(i) / 100.0
-                   pt_vertical_gradient_level_ind(i) = k - 1
-                   i = i + 1
+       IF ( .NOT. neutral )  THEN
+
+          i = 1
+          gradient = 0.0
+
+          IF ( .NOT. ocean )  THEN
+
+             pt_vertical_gradient_level_ind(1) = 0
+             DO  k = 1, nzt+1
+                IF ( i < 11 ) THEN
+                   IF ( pt_vertical_gradient_level(i) < zu(k)  .AND. &
+                        pt_vertical_gradient_level(i) >= 0.0 )  THEN
+                      gradient = pt_vertical_gradient(i) / 100.0
+                      pt_vertical_gradient_level_ind(i) = k - 1
+                      i = i + 1
+                   ENDIF
                 ENDIF
-             ENDIF
-             IF ( gradient /= 0.0 )  THEN
-                IF ( k /= 1 )  THEN
-                   pt_init(k) = pt_init(k-1) + dzu(k) * gradient
+                IF ( gradient /= 0.0 )  THEN
+                   IF ( k /= 1 )  THEN
+                      pt_init(k) = pt_init(k-1) + dzu(k) * gradient
+                   ELSE
+                      pt_init(k) = pt_surface   + 0.5 * dzu(k) * gradient
+                   ENDIF
                 ELSE
-                   pt_init(k) = pt_surface   + 0.5 * dzu(k) * gradient
+                   pt_init(k) = pt_init(k-1)
                 ENDIF
-             ELSE
-                pt_init(k) = pt_init(k-1)
-             ENDIF
-          ENDDO
-
-       ELSE
-
-          pt_vertical_gradient_level_ind(1) = nzt+1
-          DO  k = nzt, 0, -1
-             IF ( i < 11 ) THEN
-                IF ( pt_vertical_gradient_level(i) > zu(k)  .AND. &
-                     pt_vertical_gradient_level(i) <= 0.0 )  THEN
-                   gradient = pt_vertical_gradient(i) / 100.0
-                   pt_vertical_gradient_level_ind(i) = k + 1
-                   i = i + 1
+             ENDDO
+
+          ELSE
+
+             pt_vertical_gradient_level_ind(1) = nzt+1
+             DO  k = nzt, 0, -1
+                IF ( i < 11 ) THEN
+                   IF ( pt_vertical_gradient_level(i) > zu(k)  .AND. &
+                        pt_vertical_gradient_level(i) <= 0.0 )  THEN
+                      gradient = pt_vertical_gradient(i) / 100.0
+                      pt_vertical_gradient_level_ind(i) = k + 1
+                      i = i + 1
+                   ENDIF
                 ENDIF
-             ENDIF
-             IF ( gradient /= 0.0 )  THEN
-                IF ( k /= nzt )  THEN
-                   pt_init(k) = pt_init(k+1) - dzu(k+1) * gradient
+                IF ( gradient /= 0.0 )  THEN
+                   IF ( k /= nzt )  THEN
+                      pt_init(k) = pt_init(k+1) - dzu(k+1) * gradient
+                   ELSE
+                      pt_init(k)   = pt_surface - 0.5 * dzu(k+1) * gradient
+                      pt_init(k+1) = pt_surface + 0.5 * dzu(k+1) * gradient
+                   ENDIF
                 ELSE
-                   pt_init(k)   = pt_surface - 0.5 * dzu(k+1) * gradient
-                   pt_init(k+1) = pt_surface + 0.5 * dzu(k+1) * gradient
+                   pt_init(k) = pt_init(k+1)
                 ENDIF
-             ELSE
-                pt_init(k) = pt_init(k+1)
-             ENDIF
-          ENDDO
+             ENDDO
+
+          ENDIF
 
        ENDIF
@@ -1505,4 +1509,21 @@
     IF ( surface_heatflux == 9999999.9 )  constant_heatflux     = .FALSE.
     IF ( top_heatflux     == 9999999.9 )  constant_top_heatflux = .FALSE.
+
+    IF ( neutral )  THEN
+
+       IF ( surface_heatflux /= 0.0  .AND.  surface_heatflux /= 9999999.9 ) &
+       THEN
+          message_string = 'heatflux must not be set for pure neutral flow'
+          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( top_heatflux /= 0.0  .AND.  top_heatflux /= 9999999.9 ) &
+       THEN
+          message_string = 'heatflux must not be set for pure neutral flow'
+          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
+       ENDIF
+
+    ENDIF
+
     IF ( top_momentumflux_u /= 9999999.9  .AND.  &
          top_momentumflux_v /= 9999999.9 )  THEN
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 939)
+++ /palm/trunk/SOURCE/header.f90	(revision 940)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Output in case of simulations for pure neutral stratification (no pt-equation
+! solved)
 ! 
 ! Former revisions:
@@ -349,4 +350,5 @@
        ENDIF
     ENDIF
+    IF ( neutral )  WRITE ( io, 131 )  pt_surface
     IF ( humidity )  THEN
        IF ( .NOT. cloud_physics )  THEN
@@ -1682,4 +1684,6 @@
 129 FORMAT (' --> Additional prognostic equation for the specific humidity')
 130 FORMAT (' --> Additional prognostic equation for the total water content')
+131 FORMAT (' --> No pt-equation solved. Neutral stratification with pt = ', &
+                  F6.2, ' K assumed')
 132 FORMAT ('     Parameterization of long-wave radiation processes via'/ &
             '     effective emissivity scheme')
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 939)
+++ /palm/trunk/SOURCE/modules.f90	(revision 940)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! +neutral
 !
 ! Former revisions:
@@ -591,5 +591,5 @@
                 iso2d_output = .FALSE., large_scale_subsidence = .FALSE., &
                 masking_method = .FALSE., mg_switch_to_pe0 = .FALSE., &
-                netcdf_output = .FALSE., ocean = .FALSE., &
+                netcdf_output = .FALSE., neutral = .FALSE., ocean = .FALSE., &
                 outflow_l = .FALSE., outflow_n = .FALSE., outflow_r = .FALSE., &
                 outflow_s = .FALSE., passive_scalar = .FALSE., &
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 939)
+++ /palm/trunk/SOURCE/parin.f90	(revision 940)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-!
+! +neutral in inipar
 !
 ! Former revisions:
@@ -165,6 +165,6 @@
              leaf_surface_concentration, long_filter_factor, &
              loop_optimization, masking_method, mg_cycles, &
-             mg_switch_to_pe0_level, &
-             mixing_length_1d, momentum_advec, netcdf_precision, ngsrb, nsor, &
+             mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, &
+             netcdf_precision, neutral, ngsrb, nsor, &
              nsor_ini, nx, ny, nz, ocean, omega, omega_sor, &
              outflow_damping_width, overshoot_limit_e, overshoot_limit_pt, &
Index: /palm/trunk/SOURCE/poisfft.f90
===================================================================
--- /palm/trunk/SOURCE/poisfft.f90	(revision 939)
+++ /palm/trunk/SOURCE/poisfft.f90	(revision 940)
@@ -4,4 +4,6 @@
 ! Current revisions:
 ! -----------------
+! special handling of tri-array as an argument in tridia_1dd routines switched
+! off because it caused segmentation faults with intel 12.1 compiler
 !
 ! Former revisions:
@@ -1381,7 +1383,8 @@
 ! tridia)
 !
-! Attention:  when using the intel compiler, array tri must be passed as an
-!             argument to the contained subroutines. Otherwise addres faults
-!             will occur.
+! Attention:  when using the intel compilers older than 12.0, array tri must
+!             be passed as an argument to the contained subroutines. Otherwise
+!             addres faults will occur. This feature can be activated with
+!             cpp-switch __intel11
 !             On NEC, tri should not be passed (except for routine substi_1dd)
 !             because this causes very bad performance.
@@ -1425,5 +1428,5 @@
 
        IF ( j <= nnyh )  THEN
-#if defined( __lc )
+#if defined( __intel11 )
           CALL maketri_1dd( j, tri )
 #else
@@ -1431,5 +1434,5 @@
 #endif
        ELSE
-#if defined( __lc )
+#if defined( __intel11 )
           CALL maketri_1dd( ny+1-j, tri )
 #else
@@ -1437,5 +1440,5 @@
 #endif
        ENDIF
-#if defined( __lc )
+#if defined( __intel11 )
        CALL split_1dd( tri )
 #else
@@ -1446,5 +1449,5 @@
     CONTAINS
 
-#if defined( __lc )
+#if defined( __intel11 )
        SUBROUTINE maketri_1dd( j, tri )
 #else
@@ -1465,5 +1468,5 @@
           REAL, DIMENSION(0:nx) ::  l
 
-#if defined( __lc )
+#if defined( __intel11 )
           REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
 #endif
@@ -1512,5 +1515,5 @@
 
 
-#if defined( __lc )
+#if defined( __intel11 )
        SUBROUTINE split_1dd( tri )
 #else
@@ -1526,5 +1529,5 @@
           INTEGER ::  i, k
 
-#if defined( __lc )
+#if defined( __intel11 )
           REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
 #endif
Index: /palm/trunk/SOURCE/production_e.f90
===================================================================
--- /palm/trunk/SOURCE/production_e.f90	(revision 939)
+++ /palm/trunk/SOURCE/production_e.f90	(revision 940)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! TKE production by buoyancy can be switched off in case of runs with pure
+! neutral stratification
 !
 ! Former revisions:
@@ -127,8 +128,9 @@
 !       ENDIF
 
-!
-!--    Calculate TKE production by shear
+
        DO  i = nxl, nxr
 
+!
+!--       Calculate TKE production by shear
           DO  j = nys, nyn
              DO  k = nzb_diff_s_outer(j,i), nzt
@@ -453,41 +455,88 @@
 
 !
-!--       Calculate TKE production by buoyancy
-          IF ( .NOT. humidity )  THEN
-
-             IF ( use_reference )  THEN
-
-                IF ( ocean )  THEN
-!
-!--                So far in the ocean no special treatment of density flux in
-!--                the bottom and top surface layer
-                   DO  j = nys, nyn
-                      DO  k = nzb_s_inner(j,i)+1, nzt
-                         tend(k,j,i) = tend(k,j,i) +                   &
-                                       kh(k,j,i) * g / rho_reference * &
-                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
+!--       If required, calculate TKE production by buoyancy
+          IF ( .NOT. neutral )  THEN
+
+             IF ( .NOT. humidity )  THEN
+
+                IF ( use_reference )  THEN
+
+                   IF ( ocean )  THEN
+!
+!--                   So far in the ocean no special treatment of density flux
+!--                   in the bottom and top surface layer
+                      DO  j = nys, nyn
+                         DO  k = nzb_s_inner(j,i)+1, nzt
+                            tend(k,j,i) = tend(k,j,i) +                     &
+                                          kh(k,j,i) * g / rho_reference *   &
+                                          ( rho(k+1,j,i) - rho(k-1,j,i) ) * &
+                                          dd2zu(k)
+                         ENDDO
                       ENDDO
-                   ENDDO
+
+                   ELSE
+
+                      DO  j = nys, nyn
+                         DO  k = nzb_diff_s_inner(j,i), nzt_diff
+                            tend(k,j,i) = tend(k,j,i) -                   &
+                                          kh(k,j,i) * g / pt_reference *  &
+                                          ( pt(k+1,j,i) - pt(k-1,j,i) ) * &
+                                          dd2zu(k)
+                         ENDDO
+
+                         IF ( use_surface_fluxes )  THEN
+                            k = nzb_diff_s_inner(j,i)-1
+                            tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
+                                                        shf(j,i)
+                         ENDIF
+
+                         IF ( use_top_fluxes )  THEN
+                            k = nzt
+                            tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
+                                                        tswst(j,i)
+                         ENDIF
+                      ENDDO
+
+                   ENDIF
 
                 ELSE
 
-                   DO  j = nys, nyn
-                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                         tend(k,j,i) = tend(k,j,i) -                  &
-                                       kh(k,j,i) * g / pt_reference * &
-                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
+                   IF ( ocean )  THEN
+!
+!--                   So far in the ocean no special treatment of density flux
+!--                   in the bottom and top surface layer
+                      DO  j = nys, nyn
+                         DO  k = nzb_s_inner(j,i)+1, nzt
+                            tend(k,j,i) = tend(k,j,i) +                     &
+                                          kh(k,j,i) * g / rho(k,j,i) *      &
+                                          ( rho(k+1,j,i) - rho(k-1,j,i) ) * &
+                                          dd2zu(k)
+                         ENDDO
                       ENDDO
 
-                      IF ( use_surface_fluxes )  THEN
-                         k = nzb_diff_s_inner(j,i)-1
-                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
-                      ENDIF
-
-                      IF ( use_top_fluxes )  THEN
-                         k = nzt
-                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
-                                                     tswst(j,i)
-                      ENDIF
-                   ENDDO
+                   ELSE
+
+                      DO  j = nys, nyn
+                         DO  k = nzb_diff_s_inner(j,i), nzt_diff
+                            tend(k,j,i) = tend(k,j,i) -                   &
+                                          kh(k,j,i) * g / pt(k,j,i) *     &
+                                          ( pt(k+1,j,i) - pt(k-1,j,i) ) * &
+                                          dd2zu(k)
+                         ENDDO
+
+                         IF ( use_surface_fluxes )  THEN
+                            k = nzb_diff_s_inner(j,i)-1
+                            tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
+                                                        shf(j,i)
+                         ENDIF
+
+                         IF ( use_top_fluxes )  THEN
+                            k = nzt
+                            tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
+                                                        tswst(j,i)
+                         ENDIF
+                      ENDDO
+
+                   ENDIF
 
                 ENDIF
@@ -495,132 +544,96 @@
              ELSE
 
-                IF ( ocean )  THEN
-!
-!--                So far in the ocean no special treatment of density flux in
-!--                the bottom and top surface layer
-                   DO  j = nys, nyn
-                      DO  k = nzb_s_inner(j,i)+1, nzt
-                         tend(k,j,i) = tend(k,j,i) +                &
-                                       kh(k,j,i) * g / rho(k,j,i) * &
-                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
-                      ENDDO
-                   ENDDO
-
-                ELSE
-
-                   DO  j = nys, nyn
-                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                         tend(k,j,i) = tend(k,j,i) -               &
-                                       kh(k,j,i) * g / pt(k,j,i) * &
-                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
-                      ENDDO
-
-                      IF ( use_surface_fluxes )  THEN
-                         k = nzb_diff_s_inner(j,i)-1
-                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
-                      ENDIF
-
-                      IF ( use_top_fluxes )  THEN
-                         k = nzt
-                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
-                      ENDIF
-                   ENDDO
-
-                ENDIF
-
-             ENDIF
-
-          ELSE
-
-             DO  j = nys, nyn
-
-                DO  k = nzb_diff_s_inner(j,i), nzt_diff
-
-                   IF ( .NOT. cloud_physics )  THEN
-                      k1 = 1.0 + 0.61 * q(k,j,i)
-                      k2 = 0.61 * pt(k,j,i)
-                   ELSE
-                      IF ( ql(k,j,i) == 0.0 )  THEN
+                DO  j = nys, nyn
+
+                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
+
+                      IF ( .NOT. cloud_physics )  THEN
                          k1 = 1.0 + 0.61 * q(k,j,i)
                          k2 = 0.61 * pt(k,j,i)
                       ELSE
-                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                         temp  = theta * t_d_pt(k)
-                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                    ( q(k,j,i) - ql(k,j,i) ) *          &
-                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
-                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                         IF ( ql(k,j,i) == 0.0 )  THEN
+                            k1 = 1.0 + 0.61 * q(k,j,i)
+                            k2 = 0.61 * pt(k,j,i)
+                         ELSE
+                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                            temp  = theta * t_d_pt(k)
+                            k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+                                       ( q(k,j,i) - ql(k,j,i) ) *          &
+                                 ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+                                 ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                            k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                         ENDIF
                       ENDIF
-                   ENDIF
-
-                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
-                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
-                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
-                                      ) * dd2zu(k)
+
+                      tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
+                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
+                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
+                                         ) * dd2zu(k)
+                   ENDDO
+
                 ENDDO
 
-             ENDDO
-
-             IF ( use_surface_fluxes )  THEN
-
-                DO  j = nys, nyn
-
-                   k = nzb_diff_s_inner(j,i)-1
-
-                   IF ( .NOT. cloud_physics )  THEN
-                      k1 = 1.0 + 0.61 * q(k,j,i)
-                      k2 = 0.61 * pt(k,j,i)
-                   ELSE
-                      IF ( ql(k,j,i) == 0.0 )  THEN
+                IF ( use_surface_fluxes )  THEN
+
+                   DO  j = nys, nyn
+
+                      k = nzb_diff_s_inner(j,i)-1
+
+                      IF ( .NOT. cloud_physics )  THEN
                          k1 = 1.0 + 0.61 * q(k,j,i)
                          k2 = 0.61 * pt(k,j,i)
                       ELSE
-                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                         temp  = theta * t_d_pt(k)
-                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                    ( q(k,j,i) - ql(k,j,i) ) *          &
-                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
-                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                         IF ( ql(k,j,i) == 0.0 )  THEN
+                            k1 = 1.0 + 0.61 * q(k,j,i)
+                            k2 = 0.61 * pt(k,j,i)
+                         ELSE
+                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                            temp  = theta * t_d_pt(k)
+                            k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+                                       ( q(k,j,i) - ql(k,j,i) ) *          &
+                                 ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+                                 ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                            k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                         ENDIF
                       ENDIF
-                   ENDIF
-
-                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
-                                         ( k1* shf(j,i) + k2 * qsws(j,i) )
-                ENDDO
-
-             ENDIF
-
-             IF ( use_top_fluxes )  THEN
-
-                DO  j = nys, nyn
-
-                   k = nzt
-
-                   IF ( .NOT. cloud_physics )  THEN
-                      k1 = 1.0 + 0.61 * q(k,j,i)
-                      k2 = 0.61 * pt(k,j,i)
-                   ELSE
-                      IF ( ql(k,j,i) == 0.0 )  THEN
+
+                      tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
+                                            ( k1* shf(j,i) + k2 * qsws(j,i) )
+                   ENDDO
+
+                ENDIF
+
+                IF ( use_top_fluxes )  THEN
+
+                   DO  j = nys, nyn
+
+                      k = nzt
+
+                      IF ( .NOT. cloud_physics )  THEN
                          k1 = 1.0 + 0.61 * q(k,j,i)
                          k2 = 0.61 * pt(k,j,i)
                       ELSE
-                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                         temp  = theta * t_d_pt(k)
-                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                                    ( q(k,j,i) - ql(k,j,i) ) *          &
-                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
-                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                         IF ( ql(k,j,i) == 0.0 )  THEN
+                            k1 = 1.0 + 0.61 * q(k,j,i)
+                            k2 = 0.61 * pt(k,j,i)
+                         ELSE
+                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                            temp  = theta * t_d_pt(k)
+                            k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+                                       ( q(k,j,i) - ql(k,j,i) ) *          &
+                                 ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+                                 ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                            k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                         ENDIF
                       ENDIF
-                   ENDIF
-
-                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
-                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
-                ENDDO
+
+                      tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
+                                            ( k1* tswst(j,i) + k2 * qswst(j,i) )
+                   ENDDO
+
+                ENDIF
 
              ENDIF
@@ -950,33 +963,71 @@
 
 !
-!--    Calculate TKE production by buoyancy
-       IF ( .NOT. humidity )  THEN
-
-          IF ( use_reference )  THEN
-
-             IF ( ocean )  THEN
-!
-!--             So far in the ocean no special treatment of density flux in the
-!--             bottom and top surface layer
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho_reference * &
-                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
-                ENDDO
+!--    If required, calculate TKE production by buoyancy
+       IF ( .NOT. neutral )  THEN
+
+          IF ( .NOT. humidity )  THEN
+
+             IF ( use_reference )  THEN
+
+                IF ( ocean )  THEN
+!
+!--                So far in the ocean no special treatment of density flux in
+!--                the bottom and top surface layer
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tend(k,j,i) = tend(k,j,i) +                   &
+                                    kh(k,j,i) * g / rho_reference * &
+                                    ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
+                   ENDDO
+
+                ELSE
+
+                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
+                      tend(k,j,i) = tend(k,j,i) -                  &
+                                    kh(k,j,i) * g / pt_reference * &
+                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
+                   ENDDO
+
+                   IF ( use_surface_fluxes )  THEN
+                      k = nzb_diff_s_inner(j,i)-1
+                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
+                   ENDIF
+
+                   IF ( use_top_fluxes )  THEN
+                      k = nzt
+                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
+                   ENDIF
+
+                ENDIF
 
              ELSE
 
-                DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
-                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
-                ENDDO
-
-                IF ( use_surface_fluxes )  THEN
-                   k = nzb_diff_s_inner(j,i)-1
-                   tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
-                ENDIF
-
-                IF ( use_top_fluxes )  THEN
-                   k = nzt
-                   tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
+                IF ( ocean )  THEN
+!
+!--                So far in the ocean no special treatment of density flux in
+!--                the bottom and top surface layer
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tend(k,j,i) = tend(k,j,i) +                &
+                                    kh(k,j,i) * g / rho(k,j,i) * &
+                                    ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
+                   ENDDO
+
+                ELSE
+
+                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
+                      tend(k,j,i) = tend(k,j,i) -               &
+                                    kh(k,j,i) * g / pt(k,j,i) * &
+                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
+                   ENDDO
+
+                   IF ( use_surface_fluxes )  THEN
+                      k = nzb_diff_s_inner(j,i)-1
+                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
+                   ENDIF
+
+                   IF ( use_top_fluxes )  THEN
+                      k = nzt
+                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
+                   ENDIF
+
                 ENDIF
 
@@ -985,113 +1036,83 @@
           ELSE
 
-             IF ( ocean )  THEN
-!
-!--             So far in the ocean no special treatment of density flux in the
-!--             bottom and top surface layer
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho(k,j,i) * &
-                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
-                ENDDO
-
-             ELSE
-
-                DO  k = nzb_diff_s_inner(j,i), nzt_diff
-                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
-                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
-                ENDDO
-
-                IF ( use_surface_fluxes )  THEN
-                   k = nzb_diff_s_inner(j,i)-1
-                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
-                ENDIF
-
-                IF ( use_top_fluxes )  THEN
-                   k = nzt
-                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
-                ENDIF
-
-             ENDIF
-
-          ENDIF
-
-       ELSE
-
-          DO  k = nzb_diff_s_inner(j,i), nzt_diff
-
-             IF ( .NOT. cloud_physics )  THEN
-                k1 = 1.0 + 0.61 * q(k,j,i)
-                k2 = 0.61 * pt(k,j,i)
-             ELSE
-                IF ( ql(k,j,i) == 0.0 )  THEN
+             DO  k = nzb_diff_s_inner(j,i), nzt_diff
+
+                IF ( .NOT. cloud_physics )  THEN
                    k1 = 1.0 + 0.61 * q(k,j,i)
                    k2 = 0.61 * pt(k,j,i)
                 ELSE
-                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                   temp  = theta * t_d_pt(k)
-                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                              ( q(k,j,i) - ql(k,j,i) ) *          &
-                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
-                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
-                ENDIF
-             ENDIF
-
-             tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *      &
-                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
-                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
-                                      ) * dd2zu(k)
-          ENDDO
-
-          IF ( use_surface_fluxes )  THEN
-             k = nzb_diff_s_inner(j,i)-1
-
-             IF ( .NOT. cloud_physics ) THEN
-                k1 = 1.0 + 0.61 * q(k,j,i)
-                k2 = 0.61 * pt(k,j,i)
-             ELSE
-                IF ( ql(k,j,i) == 0.0 )  THEN
+                   IF ( ql(k,j,i) == 0.0 )  THEN
+                      k1 = 1.0 + 0.61 * q(k,j,i)
+                      k2 = 0.61 * pt(k,j,i)
+                   ELSE
+                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                      temp  = theta * t_d_pt(k)
+                      k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+                                 ( q(k,j,i) - ql(k,j,i) ) *          &
+                           ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+                           ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                      k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                   ENDIF
+                ENDIF
+
+                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *      &
+                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
+                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
+                                         ) * dd2zu(k)
+             ENDDO
+
+             IF ( use_surface_fluxes )  THEN
+                k = nzb_diff_s_inner(j,i)-1
+
+                IF ( .NOT. cloud_physics ) THEN
                    k1 = 1.0 + 0.61 * q(k,j,i)
                    k2 = 0.61 * pt(k,j,i)
                 ELSE
-                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                   temp  = theta * t_d_pt(k)
-                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                              ( q(k,j,i) - ql(k,j,i) ) *          &
-                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
-                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
-                ENDIF
+                   IF ( ql(k,j,i) == 0.0 )  THEN
+                      k1 = 1.0 + 0.61 * q(k,j,i)
+                      k2 = 0.61 * pt(k,j,i)
+                   ELSE
+                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                      temp  = theta * t_d_pt(k)
+                      k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+                                 ( q(k,j,i) - ql(k,j,i) ) *          &
+                           ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+                           ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                      k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                   ENDIF
+                ENDIF
+
+                tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
+                                            ( k1* shf(j,i) + k2 * qsws(j,i) )
              ENDIF
 
-             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
-                                         ( k1* shf(j,i) + k2 * qsws(j,i) )
-          ENDIF
-
-          IF ( use_top_fluxes )  THEN
-             k = nzt
-
-             IF ( .NOT. cloud_physics ) THEN
-                k1 = 1.0 + 0.61 * q(k,j,i)
-                k2 = 0.61 * pt(k,j,i)
-             ELSE
-                IF ( ql(k,j,i) == 0.0 )  THEN
+             IF ( use_top_fluxes )  THEN
+                k = nzt
+
+                IF ( .NOT. cloud_physics ) THEN
                    k1 = 1.0 + 0.61 * q(k,j,i)
                    k2 = 0.61 * pt(k,j,i)
                 ELSE
-                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
-                   temp  = theta * t_d_pt(k)
-                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
-                              ( q(k,j,i) - ql(k,j,i) ) *          &
-                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
-                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
-                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
-                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
-                ENDIF
+                   IF ( ql(k,j,i) == 0.0 )  THEN
+                      k1 = 1.0 + 0.61 * q(k,j,i)
+                      k2 = 0.61 * pt(k,j,i)
+                   ELSE
+                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                      temp  = theta * t_d_pt(k)
+                      k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+                                 ( q(k,j,i) - ql(k,j,i) ) *          &
+                           ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+                           ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                      k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                   ENDIF
+                ENDIF
+
+                tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
+                                            ( k1* tswst(j,i) + k2 * qswst(j,i) )
              ENDIF
 
-             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
-                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
           ENDIF
 
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 939)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 940)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! temperature equation can be switched off
 !
 ! Former revisions:
@@ -159,7 +159,7 @@
 !-- Calculate those variables needed in the tendency terms which need
 !-- global communication
-    CALL calc_mean_profile( pt, 4 )
-    IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
-    IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
+    IF ( .NOT. neutral )  CALL calc_mean_profile( pt, 4 )
+    IF ( ocean         )  CALL calc_mean_profile( rho, 64 )
+    IF ( humidity      )  CALL calc_mean_profile( vpt, 44 )
     IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
          intermediate_timestep_count == 1 )  CALL ws_statistics
@@ -205,5 +205,7 @@
           ENDIF
           CALL coriolis( i, j, 1 )
-          IF ( sloping_surface )  CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
+          IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
+             CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
+          ENDIF
 
 !
@@ -375,11 +377,14 @@
           ENDIF
           CALL coriolis( i, j, 3 )
-          IF ( ocean )  THEN
-             CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
-          ELSE
-             IF ( .NOT. humidity )  THEN
-                CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
-             ELSE
-                CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
+
+          IF ( .NOT. neutral )  THEN
+             IF ( ocean )  THEN
+                CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
+             ELSE
+                IF ( .NOT. humidity )  THEN
+                   CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
+                ELSE
+                   CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
+                ENDIF
              ENDIF
           ENDIF
@@ -422,119 +427,126 @@
 
 !
-!-- potential temperature
-    CALL cpu_log( log_point(13), 'pt-equation', 'start' )
-
-!
-!-- pt-tendency terms with communication
-    sat = tsc(1)
-    sbt = tsc(2)
-    IF ( scalar_advec == 'bc-scheme' )  THEN
-
-       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
-!
-!--       Bott-Chlond scheme always uses Euler time step when leapfrog is
-!--       switched on. Thus:
-          sat = 1.0
-          sbt = 1.0
-       ENDIF
-       tend = 0.0
-       CALL advec_s_bc( pt, 'pt' )
-    ELSE
-       IF ( tsc(2) /= 2.0  .AND.  scalar_advec == 'ups-scheme' )  THEN
+!-- If required, compute prognostic equation for potential temperature
+    IF ( .NOT. neutral )  THEN
+
+       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
+
+   !
+   !-- pt-tendency terms with communication
+       sat = tsc(1)
+       sbt = tsc(2)
+       IF ( scalar_advec == 'bc-scheme' )  THEN
+
+          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+   !
+   !--       Bott-Chlond scheme always uses Euler time step when leapfrog is
+   !--       switched on. Thus:
+             sat = 1.0
+             sbt = 1.0
+          ENDIF
           tend = 0.0
-          CALL advec_s_ups( pt, 'pt' )
-       ENDIF
-    ENDIF
-
-!
-!-- pt-tendency terms with no communication
-    DO  i = nxl, nxr
-       DO  j = nys, nyn
-!
-!--       Tendency terms
-          IF ( scalar_advec == 'bc-scheme' )  THEN
-             CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
-                               wall_heatflux, tend )
-          ELSE
-             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
-                tend(:,j,i) = 0.0
-                IF ( ws_scheme_sca )  THEN
-                   CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, &
-                                 diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn )
-                ELSE
-                    CALL advec_s_pw( i, j, pt )
-                ENDIF
-             ELSE
-                IF ( scalar_advec /= 'ups-scheme' )  THEN
-                   tend(:,j,i) = 0.0
-                   CALL advec_s_up( i, j, pt )
-                ENDIF
-             ENDIF
-             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
-             THEN
-                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
-                                  tswst_m, wall_heatflux, tend )
-             ELSE
+          CALL advec_s_bc( pt, 'pt' )
+       ELSE
+          IF ( tsc(2) /= 2.0  .AND.  scalar_advec == 'ups-scheme' )  THEN
+             tend = 0.0
+             CALL advec_s_ups( pt, 'pt' )
+          ENDIF
+       ENDIF
+
+   !
+   !-- pt-tendency terms with no communication
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+   !
+   !--       Tendency terms
+             IF ( scalar_advec == 'bc-scheme' )  THEN
                 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
                                   wall_heatflux, tend )
-             ENDIF
-          ENDIF
-
-!
-!--       If required compute heating/cooling due to long wave radiation
-!--       processes
-          IF ( radiation )  THEN
-             CALL calc_radiation( i, j )
-          ENDIF
-
-!
-!--       If required compute impact of latent heat due to precipitation
-          IF ( precipitation )  THEN
-             CALL impact_of_latent_heat( i, j )
-          ENDIF
-
-!
-!--       Consideration of heat sources within the plant canopy
-          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
-
-          CALL user_actions( i, j, 'pt-tendency' )
-
-!
-!--       Prognostic equation for potential temperature
-          DO  k = nzb_s_inner(j,i)+1, nzt
-             pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
-                           dt_3d * (                                           &
-                                     sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
-                                   ) -                                         &
-                           tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
+             ELSE
+                IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
+                THEN
+                   tend(:,j,i) = 0.0
+                   IF ( ws_scheme_sca )  THEN
+                      CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt,       &
+                                       diss_s_pt, flux_l_pt, diss_l_pt, &
+                                       i_omp_start, tn )
+                   ELSE
+                       CALL advec_s_pw( i, j, pt )
+                   ENDIF
+                ELSE
+                   IF ( scalar_advec /= 'ups-scheme' )  THEN
+                      tend(:,j,i) = 0.0
+                      CALL advec_s_up( i, j, pt )
+                   ENDIF
+                ENDIF
+                IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
+                THEN
+                   CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
+                                     tswst_m, wall_heatflux, tend )
+                ELSE
+                   CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
+                                     wall_heatflux, tend )
+                ENDIF
+             ENDIF
+
+   !
+   !--       If required compute heating/cooling due to long wave radiation
+   !--       processes
+             IF ( radiation )  THEN
+                CALL calc_radiation( i, j )
+             ENDIF
+
+   !
+   !--       If required compute impact of latent heat due to precipitation
+             IF ( precipitation )  THEN
+                CALL impact_of_latent_heat( i, j )
+             ENDIF
+
+   !
+   !--       Consideration of heat sources within the plant canopy
+             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
+
+             CALL user_actions( i, j, 'pt-tendency' )
+
+   !
+   !--       Prognostic equation for potential temperature
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
+                              dt_3d * (                                           &
+                                       sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
+                                      ) -                                         &
+                              tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
+             ENDDO
+
+   !
+   !--       Calculate tendencies for the next Runge-Kutta step
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( intermediate_timestep_count == 1 )  THEN
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tpt_m(k,j,i) = tend(k,j,i)
+                   ENDDO
+                ELSEIF ( intermediate_timestep_count < &
+                         intermediate_timestep_count_max )  THEN
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                     5.3125 * tpt_m(k,j,i)
+                   ENDDO
+                ENDIF
+             ENDIF
+
           ENDDO
-
-!
-!--       Calculate tendencies for the next Runge-Kutta step
-          IF ( timestep_scheme(1:5) == 'runge' )  THEN
-             IF ( intermediate_timestep_count == 1 )  THEN
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   tpt_m(k,j,i) = tend(k,j,i)
-                ENDDO
-             ELSEIF ( intermediate_timestep_count < &
-                      intermediate_timestep_count_max )  THEN
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
-                ENDDO
-             ENDIF
-          ENDIF
-
        ENDDO
-    ENDDO
-
-    CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
+
+       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
+
+    ENDIF
 
 !
@@ -714,6 +726,6 @@
 !
 !--          If required compute influence of large-scale subsidence/ascent
-             IF ( large_scale_subsidence ) THEN
-                CALL subsidence ( i, j, tend, q, q_init )
+             IF ( large_scale_subsidence )  THEN
+                CALL subsidence( i, j, tend, q, q_init )
              ENDIF
 
@@ -936,7 +948,7 @@
 !-- Calculate those variables needed in the tendency terms which need
 !-- global communication
-    CALL calc_mean_profile( pt, 4 )
-    IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
-    IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
+    IF ( .NOT. neutral )  CALL calc_mean_profile( pt, 4 )
+    IF ( ocean         )  CALL calc_mean_profile( rho, 64 )
+    IF ( humidity      )  CALL calc_mean_profile( vpt, 44 )
     IF ( .NOT. constant_diffusion )  CALL production_e_init
     IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
@@ -989,6 +1001,7 @@
              ENDIF
              CALL coriolis( i, j, 1 )
-             IF ( sloping_surface )  CALL buoyancy( i, j, pt, pt_reference, 1, &
-                                                    4 )
+             IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
+                CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
+             ENDIF
 
 !
@@ -1121,11 +1134,14 @@
           ENDIF
           CALL coriolis( i, j, 3 )
-          IF ( ocean )  THEN
-             CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
-          ELSE
-             IF ( .NOT. humidity )  THEN
-                CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
-             ELSE
-                CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
+
+          IF ( .NOT. neutral )  THEN
+             IF ( ocean )  THEN
+                CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
+             ELSE
+                IF ( .NOT. humidity )  THEN
+                   CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
+                ELSE
+                   CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
+                ENDIF
              ENDIF
           ENDIF
@@ -1163,78 +1179,82 @@
 
 !
-!--       Tendency terms for potential temperature
-          tend(:,j,i) = 0.0
-          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
-                IF ( ws_scheme_sca )  THEN
-       !            CALL local_diss( i, j, pt )
-                   CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, &
-                             diss_s_pt, flux_l_pt, diss_l_pt, i_omp_start, tn )
-                ELSE
-                   CALL advec_s_pw( i, j, pt )
-                ENDIF
-          ELSE
-             CALL advec_s_up( i, j, pt )
-          ENDIF
-          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
-          THEN
-             CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
-                               tswst_m, wall_heatflux, tend )
-          ELSE
-             CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
-                               wall_heatflux, tend )
-          ENDIF
-
-!
-!--       If required compute heating/cooling due to long wave radiation
-!--       processes
-          IF ( radiation )  THEN
-             CALL calc_radiation( i, j )
-          ENDIF
-
-!
-!--       If required compute impact of latent heat due to precipitation
-          IF ( precipitation )  THEN
-             CALL impact_of_latent_heat( i, j )
-          ENDIF
-
-!
-!--       Consideration of heat sources within the plant canopy
-          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
-
-
-          CALL user_actions( i, j, 'pt-tendency' )
-
-!
-!--       Prognostic equation for potential temperature
-          DO  k = nzb_s_inner(j,i)+1, nzt
-             pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
-                           dt_3d * (                                        &
-                               tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
-                                   ) -                                      &
-                           tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
-          ENDDO
-
-!
-!--       Calculate tendencies for the next Runge-Kutta step
-          IF ( timestep_scheme(1:5) == 'runge' )  THEN
-             IF ( intermediate_timestep_count == 1 )  THEN
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   tpt_m(k,j,i) = tend(k,j,i)
-                ENDDO
-             ELSEIF ( intermediate_timestep_count < &
-                      intermediate_timestep_count_max )  THEN
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
-                                   5.3125 * tpt_m(k,j,i)
-                ENDDO
-             ENDIF
+!--       If required, compute prognostic equation for potential temperature
+          IF ( .NOT. neutral )  THEN
+!
+!--          Tendency terms for potential temperature
+             tend(:,j,i) = 0.0
+             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
+                   IF ( ws_scheme_sca )  THEN
+                      CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt, &
+                                       flux_l_pt, diss_l_pt, i_omp_start, tn )
+                   ELSE
+                      CALL advec_s_pw( i, j, pt )
+                   ENDIF
+             ELSE
+                CALL advec_s_up( i, j, pt )
+             ENDIF
+             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
+             THEN
+                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
+                                  tswst_m, wall_heatflux, tend )
+             ELSE
+                CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
+                                  wall_heatflux, tend )
+             ENDIF
+
+!
+!--          If required compute heating/cooling due to long wave radiation
+!--          processes
+             IF ( radiation )  THEN
+                CALL calc_radiation( i, j )
+             ENDIF
+
+!
+!--          If required compute impact of latent heat due to precipitation
+             IF ( precipitation )  THEN
+                CALL impact_of_latent_heat( i, j )
+             ENDIF
+
+!
+!--          Consideration of heat sources within the plant canopy
+             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
+
+
+             CALL user_actions( i, j, 'pt-tendency' )
+
+!
+!--          Prognostic equation for potential temperature
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
+                              dt_3d * (                                        &
+                                  tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
+                                      ) -                                      &
+                              tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
+             ENDDO
+
+!
+!--          Calculate tendencies for the next Runge-Kutta step
+             IF ( timestep_scheme(1:5) == 'runge' )  THEN
+                IF ( intermediate_timestep_count == 1 )  THEN
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tpt_m(k,j,i) = tend(k,j,i)
+                   ENDDO
+                ELSEIF ( intermediate_timestep_count < &
+                         intermediate_timestep_count_max )  THEN
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                      5.3125 * tpt_m(k,j,i)
+                   ENDDO
+                ENDIF
+             ENDIF
+
           ENDIF
 
@@ -1337,6 +1357,6 @@
 
 !--          If required compute influence of large-scale subsidence/ascent
-             IF ( large_scale_subsidence ) THEN
-                CALL subsidence ( i, j, tend, q, q_init )
+             IF ( large_scale_subsidence )  THEN
+                CALL subsidence( i, j, tend, q, q_init )
              ENDIF
 
@@ -1484,7 +1504,7 @@
 !-- Calculate those variables needed in the tendency terms which need
 !-- global communication
-    CALL calc_mean_profile( pt, 4 )
-    IF ( ocean    )  CALL calc_mean_profile( rho, 64 )
-    IF ( humidity )  CALL calc_mean_profile( vpt, 44 )
+    IF ( .NOT. neutral )  CALL calc_mean_profile( pt, 4 )
+    IF ( ocean         )  CALL calc_mean_profile( rho, 64 )
+    IF ( humidity      )  CALL calc_mean_profile( vpt, 44 )
     IF ( ( ws_scheme_mom .OR. ws_scheme_sca )  .AND.  &
          intermediate_timestep_count == 1 )  CALL ws_statistics
@@ -1523,5 +1543,7 @@
     ENDIF
     CALL coriolis( 1 )
-    IF ( sloping_surface )  CALL buoyancy( pt, pt_reference, 1, 4 )
+    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
+       CALL buoyancy( pt, pt_reference, 1, 4 )
+    ENDIF
 
 !
@@ -1706,11 +1728,14 @@
     ENDIF
     CALL coriolis( 3 )
-    IF ( ocean )  THEN
-       CALL buoyancy( rho, rho_reference, 3, 64 )
-    ELSE
-       IF ( .NOT. humidity )  THEN
-          CALL buoyancy( pt, pt_reference, 3, 4 )
+
+    IF ( .NOT. neutral )  THEN
+       IF ( ocean )  THEN
+          CALL buoyancy( rho, rho_reference, 3, 64 )
        ELSE
-          CALL buoyancy( vpt, pt_reference, 3, 44 )
+          IF ( .NOT. humidity )  THEN
+             CALL buoyancy( pt, pt_reference, 3, 4 )
+          ELSE
+             CALL buoyancy( vpt, pt_reference, 3, 44 )
+          ENDIF
        ENDIF
     ENDIF
@@ -1761,122 +1786,128 @@
     CALL cpu_log( log_point(7), 'w-equation', 'stop' )
 
-!
-!-- potential temperature
-    CALL cpu_log( log_point(13), 'pt-equation', 'start' )
-
-!
-!-- pt-tendency terms with communication
-    sat = tsc(1)
-    sbt = tsc(2)
-    IF ( scalar_advec == 'bc-scheme' )  THEN
-
-       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
-!
-!--       Bott-Chlond scheme always uses Euler time step when leapfrog is
-!--       switched on. Thus:
-          sat = 1.0
-          sbt = 1.0
-       ENDIF
-       tend = 0.0
-       CALL advec_s_bc( pt, 'pt' )
-    ELSE
-       IF ( tsc(2) /= 2.0  .AND.  scalar_advec == 'ups-scheme' )  THEN
+
+!
+!-- If required, compute prognostic equation for potential temperature
+    IF ( .NOT. neutral )  THEN
+
+       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
+
+!
+!--    pt-tendency terms with communication
+       sat = tsc(1)
+       sbt = tsc(2)
+       IF ( scalar_advec == 'bc-scheme' )  THEN
+
+          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+!
+!--          Bott-Chlond scheme always uses Euler time step when leapfrog is
+!--          switched on. Thus:
+             sat = 1.0
+             sbt = 1.0
+          ENDIF
           tend = 0.0
-          CALL advec_s_ups( pt, 'pt' )
-       ENDIF
-    ENDIF
-          
-!
-!-- pt-tendency terms with no communication
-    IF ( scalar_advec == 'bc-scheme' )  THEN
-       CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
-                         tend )
-    ELSE
-       IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
-          tend = 0.0
-          IF ( ws_scheme_sca )  THEN
-             CALL advec_s_ws( pt, 'pt' )
-          ELSE
-             CALL advec_s_pw( pt )
-          ENDIF
+          CALL advec_s_bc( pt, 'pt' )
        ELSE
-          IF ( scalar_advec /= 'ups-scheme' )  THEN
+          IF ( tsc(2) /= 2.0  .AND.  scalar_advec == 'ups-scheme' )  THEN
              tend = 0.0
-             CALL advec_s_up( pt )
-          ENDIF
-       ENDIF
-       IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
-          CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, &
-                            wall_heatflux, tend )
-       ELSE
+             CALL advec_s_ups( pt, 'pt' )
+          ENDIF
+       ENDIF
+
+!
+!--    pt-tendency terms with no communication
+       IF ( scalar_advec == 'bc-scheme' )  THEN
           CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
                             tend )
-       ENDIF
-    ENDIF
-
-!
-!-- If required compute heating/cooling due to long wave radiation
-!-- processes
-    IF ( radiation )  THEN
-       CALL calc_radiation
-    ENDIF
-
-!
-!-- If required compute impact of latent heat due to precipitation
-    IF ( precipitation )  THEN
-       CALL impact_of_latent_heat
-    ENDIF
-
-!
-!-- Consideration of heat sources within the plant canopy
-    IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
-       CALL plant_canopy_model( 4 )
-    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' )
-
-!
-!-- Prognostic equation for potential temperature
-    DO  i = nxl, nxr
-       DO  j = nys, nyn
-          DO  k = nzb_s_inner(j,i)+1, nzt
-             pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) +       &
-                           dt_3d * (                                           &
-                                     sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
-                                   ) -                                         &
-                           tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
+       ELSE
+          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
+             tend = 0.0
+             IF ( ws_scheme_sca )  THEN
+                CALL advec_s_ws( pt, 'pt' )
+             ELSE
+                CALL advec_s_pw( pt )
+             ENDIF
+          ELSE
+             IF ( scalar_advec /= 'ups-scheme' )  THEN
+                tend = 0.0
+                CALL advec_s_up( pt )
+             ENDIF
+          ENDIF
+          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
+             CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, &
+                               wall_heatflux, tend )
+          ELSE
+             CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
+                               tend )
+          ENDIF
+       ENDIF
+
+!
+!--    If required compute heating/cooling due to long wave radiation processes
+       IF ( radiation )  THEN
+          CALL calc_radiation
+       ENDIF
+
+!
+!--    If required compute impact of latent heat due to precipitation
+       IF ( precipitation )  THEN
+          CALL impact_of_latent_heat
+       ENDIF
+
+!
+!--    Consideration of heat sources within the plant canopy
+       IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
+          CALL plant_canopy_model( 4 )
+       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' )
+
+!
+!--    Prognostic equation for potential temperature
+       DO  i = nxl, nxr
+          DO  j = nys, nyn
+             DO  k = nzb_s_inner(j,i)+1, nzt
+                pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) +       &
+                              dt_3d * (                                           &
+                                       sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
+                                      ) -                                         &
+                              tsc(5) * rdf_sc(k) * ( pt(k,j,i) - pt_init(k) )
+             ENDDO
           ENDDO
        ENDDO
-    ENDDO
-
-!
-!-- Calculate tendencies for the next Runge-Kutta step
-    IF ( timestep_scheme(1:5) == 'runge' )  THEN
-       IF ( intermediate_timestep_count == 1 )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   tpt_m(k,j,i) = tend(k,j,i)
-                ENDDO
-             ENDDO
-          ENDDO
-       ELSEIF ( intermediate_timestep_count < &
-                intermediate_timestep_count_max )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt
-                   tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-    ENDIF
-
-    CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
+
+!
+!--    Calculate tendencies for the next Runge-Kutta step
+       IF ( timestep_scheme(1:5) == 'runge' )  THEN
+          IF ( intermediate_timestep_count == 1 )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tpt_m(k,j,i) = tend(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSEIF ( intermediate_timestep_count < &
+                   intermediate_timestep_count_max )  THEN
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, nzt
+                      tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
+                                      5.3125 * tpt_m(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+       ENDIF
+
+       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
+
+    ENDIF
 
 !
@@ -2054,6 +2085,6 @@
 !
 !--    If required compute influence of large-scale subsidence/ascent
-       IF ( large_scale_subsidence ) THEN
-         CALL subsidence ( tend, q, q_init )
+       IF ( large_scale_subsidence )  THEN
+         CALL subsidence( tend, q, q_init )
        ENDIF
 
Index: /palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/read_var_list.f90	(revision 939)
+++ /palm/trunk/SOURCE/read_var_list.f90	(revision 940)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! +neutral
 !
 ! Former revisions:
@@ -421,4 +421,6 @@
           CASE ( 'netcdf_precision' )
              READ ( 13 )  netcdf_precision
+          CASE ( 'neutral' )
+             READ ( 13 )  neutral
           CASE ( 'ngsrb' )
              READ ( 13 )  ngsrb
Index: /palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/write_var_list.f90	(revision 939)
+++ /palm/trunk/SOURCE/write_var_list.f90	(revision 940)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! +neutral
 !
 ! Former revisions:
@@ -335,4 +335,6 @@
     WRITE ( 14 )  'netcdf_precision              '
     WRITE ( 14 )  netcdf_precision
+    WRITE ( 14 )  'neutral                       '
+    WRITE ( 14 )  neutral
     WRITE ( 14 )  'ngsrb                         '
     WRITE ( 14 )  ngsrb
