Index: /palm/trunk/SCRIPTS/.mrun.config.default
===================================================================
--- /palm/trunk/SCRIPTS/.mrun.config.default	(revision 62)
+++ /palm/trunk/SCRIPTS/.mrun.config.default	(revision 63)
@@ -59,5 +59,5 @@
 %cpumax            1000                                       ibmb parallel
 %numprocs          4                                          ibmb parallel
-#%tmp_data_catalog  /fastfs/work/<replace by your HLRN username>/palm_restart_data    ibmh parallel debug
+#%tmp_data_catalog  /fastfs/work/<replace by your HLRN username>/palm_restart_data    ibmb parallel debug
 #
 #%remote_username   <replace by your DKRZ username>           nech parallel
Index: /palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 62)
+++ /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 63)
@@ -9,4 +9,6 @@
 
 new d3par-parameter dt_max to define the maximum value for the allowed timestep 
+
+new inipar-parameter loop_optimization to control the loop optimization method
 
 new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).
Index: /palm/trunk/SOURCE/advec_particles.f90
===================================================================
--- /palm/trunk/SOURCE/advec_particles.f90	(revision 62)
+++ /palm/trunk/SOURCE/advec_particles.f90	(revision 63)
@@ -2315,5 +2315,5 @@
           ENDIF
 
-          CALL MPI_SENDRECV( trlp, trlp_count, mpi_particle_type, pleft, 1,    &
+          CALL MPI_SENDRECV( trlp(1), trlp_count, mpi_particle_type, pleft, 1, &
                              particles(number_of_particles+1), trrp_count_recv,&
                              mpi_particle_type, pright, 1,                     &
@@ -2343,5 +2343,6 @@
              ENDIF
 
-             CALL MPI_SENDRECV( trlpt, trlpt_count*tlength, MPI_REAL, pleft, 1,&
+             CALL MPI_SENDRECV( trlpt(1,1,1), trlpt_count*tlength, MPI_REAL,   &
+                                pleft, 1,                                      &
                              particle_tail_coordinates(1,1,number_of_tails+1), &
                                 trrpt_count_recv*tlength, MPI_REAL, pright, 1, &
@@ -2391,5 +2392,5 @@
           ENDIF
 
-          CALL MPI_SENDRECV( trrp, trrp_count, mpi_particle_type, pright, 1,   &
+          CALL MPI_SENDRECV( trrp(1), trrp_count, mpi_particle_type, pright, 1,&
                              particles(number_of_particles+1), trlp_count_recv,&
                              mpi_particle_type, pleft, 1,                      &
@@ -2419,6 +2420,7 @@
              ENDIF
 
-             CALL MPI_SENDRECV( trrpt, trrpt_count*tlength, MPI_REAL, pright,  &
-                          1, particle_tail_coordinates(1,1,number_of_tails+1), &
+             CALL MPI_SENDRECV( trrpt(1,1,1), trrpt_count*tlength, MPI_REAL,   &
+                                pright, 1,                                     &
+                             particle_tail_coordinates(1,1,number_of_tails+1), &
                                 trlpt_count_recv*tlength, MPI_REAL, pleft, 1,  &
                                 comm2d, status, ierr )
@@ -2754,5 +2756,5 @@
           ENDIF
 
-          CALL MPI_SENDRECV( trsp, trsp_count, mpi_particle_type, psouth, 1,   &
+          CALL MPI_SENDRECV( trsp(1), trsp_count, mpi_particle_type, psouth, 1,&
                              particles(number_of_particles+1), trnp_count_recv,&
                              mpi_particle_type, pnorth, 1,                     &
@@ -2782,6 +2784,7 @@
              ENDIF
 
-             CALL MPI_SENDRECV( trspt, trspt_count*tlength, MPI_REAL, psouth,  &
-                          1, particle_tail_coordinates(1,1,number_of_tails+1), &
+             CALL MPI_SENDRECV( trspt(1,1,1), trspt_count*tlength, MPI_REAL,   &
+                                psouth, 1,                                     &
+                             particle_tail_coordinates(1,1,number_of_tails+1), &
                                 trnpt_count_recv*tlength, MPI_REAL, pnorth, 1, &
                                 comm2d, status, ierr )
@@ -2830,5 +2833,5 @@
           ENDIF
 
-          CALL MPI_SENDRECV( trnp, trnp_count, mpi_particle_type, pnorth, 1,   &
+          CALL MPI_SENDRECV( trnp(1), trnp_count, mpi_particle_type, pnorth, 1,&
                              particles(number_of_particles+1), trsp_count_recv,&
                              mpi_particle_type, psouth, 1,                     &
@@ -2858,6 +2861,7 @@
              ENDIF
 
-             CALL MPI_SENDRECV( trnpt, trnpt_count*tlength, MPI_REAL, pnorth,  &
-                          1, particle_tail_coordinates(1,1,number_of_tails+1), &
+             CALL MPI_SENDRECV( trnpt(1,1,1), trnpt_count*tlength, MPI_REAL,   &
+                                pnorth, 1,                                     &
+                             particle_tail_coordinates(1,1,number_of_tails+1), &
                                 trspt_count_recv*tlength, MPI_REAL, psouth, 1, &
                                 comm2d, status, ierr )
Index: /palm/trunk/SOURCE/advec_s_bc.f90
===================================================================
--- /palm/trunk/SOURCE/advec_s_bc.f90	(revision 62)
+++ /palm/trunk/SOURCE/advec_s_bc.f90	(revision 63)
@@ -1,8 +1,8 @@
  SUBROUTINE advec_s_bc( sk, sk_char )
 
-!-------------------------------------------------------------------------------!
+!------------------------------------------------------------------------------!
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation extended for gridpoint nzt
 !
 ! Former revisions:
@@ -29,8 +29,9 @@
 ! stretched grids.
 ! NOTE: This is a provisional, non-optimised version!
-!-------------------------------------------------------------------------------!
+!------------------------------------------------------------------------------!
 
     USE advection
     USE arrays_3d
+    USE control_parameters
     USE cpulog
     USE grid_variables
@@ -39,5 +40,4 @@
     USE pegrid
     USE statistics
-    USE control_parameters
 
     IMPLICIT NONE
@@ -57,5 +57,5 @@
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::  sk_p
 
-#if defined( __t3eb ) || defined( __t3eh ) || defined( __t3ej2 ) || defined( __t3ej5 )
+#if defined( __nec )
     REAL (kind=4) ::  m1n, m1z  !Wichtig: Division
     REAL (kind=4), DIMENSION(:,:), ALLOCATABLE :: m1, sw
@@ -68,5 +68,5 @@
 !
 !-- Array sk_p requires 2 extra elements for each dimension
-    ALLOCATE( sk_p(nzb-2:nzt+2,nys-3:nyn+3,nxl-3:nxr+3) )
+    ALLOCATE( sk_p(nzb-2:nzt+3,nys-3:nyn+3,nxl-3:nxr+3) )
     sk_p = 0.0
 
@@ -89,5 +89,5 @@
     DO  i = nxl-1, nxr+1
        DO  j = nys, nyn
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb, nzt+1
              sk_p(k,j,i) = sk(k,j,i)
           ENDDO
@@ -95,5 +95,5 @@
     ENDDO
 #if defined( __parallel )
-    ngp = 2 * ( nzt - nzb + 5 ) * ( nyn - nys + 7 )
+    ngp = 2 * ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
     CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'start' )
 !
@@ -142,5 +142,5 @@
     DO  i = nxl, nxr
        DO  j = nys, nyn
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              zaehler = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )
              nenner  = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
@@ -160,12 +160,12 @@
 !
 !-- Allocate temporary arrays
-    ALLOCATE( a0(nzb+1:nzt-1,nxl-1:nxr+1),   a1(nzb+1:nzt-1,nxl-1:nxr+1),   &
-              a2(nzb+1:nzt-1,nxl-1:nxr+1),   a12(nzb+1:nzt-1,nxl-1:nxr+1),  &
-              a22(nzb+1:nzt-1,nxl-1:nxr+1),  immb(nzb+1:nzt-1,nxl-1:nxr+1), &
-              imme(nzb+1:nzt-1,nxl-1:nxr+1), impb(nzb+1:nzt-1,nxl-1:nxr+1), &
-              impe(nzb+1:nzt-1,nxl-1:nxr+1), ipmb(nzb+1:nzt-1,nxl-1:nxr+1), &
-              ipme(nzb+1:nzt-1,nxl-1:nxr+1), ippb(nzb+1:nzt-1,nxl-1:nxr+1), &
-              ippe(nzb+1:nzt-1,nxl-1:nxr+1), m1(nzb+1:nzt-1,nxl-2:nxr+2),   &
-              sw(nzb+1:nzt-1,nxl-1:nxr+1)                                   &
+    ALLOCATE( a0(nzb+1:nzt,nxl-1:nxr+1),   a1(nzb+1:nzt,nxl-1:nxr+1),   &
+              a2(nzb+1:nzt,nxl-1:nxr+1),   a12(nzb+1:nzt,nxl-1:nxr+1),  &
+              a22(nzb+1:nzt,nxl-1:nxr+1),  immb(nzb+1:nzt,nxl-1:nxr+1), &
+              imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1), &
+              impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1), &
+              ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1), &
+              ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2),   &
+              sw(nzb+1:nzt,nxl-1:nxr+1)                                 &
             )
     imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
@@ -184,5 +184,5 @@
 !--    Compute polynomial coefficients
        DO  i = nxl-1, nxr+1
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              a12(k,i) = 0.5 * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
              a22(k,i) = 0.5 * ( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) &
@@ -204,5 +204,5 @@
 !--    *VOCL LOOP,UNROLL(2)
        DO  i = nxl, nxr
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              cip  =  MAX( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
              cim  = -MIN( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
@@ -240,5 +240,5 @@
 !--    Compute monitor function m1
        DO  i = nxl-2, nxr+2
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              m1z = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )
              m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
@@ -258,5 +258,5 @@
        sw = 0.0
        DO  i = nxl-1, nxr+1
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              m2 = 2.0 * ABS( a1(k,i) - a12(k,i) ) / &
                   MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35 )
@@ -284,5 +284,5 @@
        CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
        DO  i = nxl, nxr
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
 
 !--          *VOCL STMT,IF(10)
@@ -374,5 +374,5 @@
 !--    Prognostic equation
        DO  i = nxl, nxr
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              fplus  = ( 1.0 - sw(k,i)   ) * ippb(k,i) + sw(k,i)   * ippe(k,i) &
                     - ( 1.0 - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i)
@@ -404,6 +404,6 @@
 !-- Enlarge boundary of local array cyclically in y-direction
 #if defined( __parallel )
-    ngp = ( nzt - nzb + 5 ) * ( nyn - nys + 7 )
-    CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+5), ngp, MPI_REAL, &
+    ngp = ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
+    CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL, &
                           type_xz_2, ierr )
     CALL MPI_TYPE_COMMIT( type_xz_2, ierr )
@@ -423,5 +423,5 @@
 #else
     DO  i = nxl, nxr
-       DO  k = nzb+1, nzt-1
+       DO  k = nzb+1, nzt
           sk_p(k,nys-1,i) = sk_p(k,nyn,i)
           sk_p(k,nys-2,i) = sk_p(k,nyn-1,i)
@@ -439,5 +439,5 @@
     DO  i = nxl, nxr
        DO  j = nys, nyn
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              zaehler = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )
              nenner  = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
@@ -457,12 +457,12 @@
 !
 !-- Allocate temporary arrays
-    ALLOCATE( a0(nzb+1:nzt-1,nys-1:nyn+1),   a1(nzb+1:nzt-1,nys-1:nyn+1),   &
-              a2(nzb+1:nzt-1,nys-1:nyn+1),   a12(nzb+1:nzt-1,nys-1:nyn+1),  &
-              a22(nzb+1:nzt-1,nys-1:nyn+1),  immb(nzb+1:nzt-1,nys-1:nyn+1), &
-              imme(nzb+1:nzt-1,nys-1:nyn+1), impb(nzb+1:nzt-1,nys-1:nyn+1), &
-              impe(nzb+1:nzt-1,nys-1:nyn+1), ipmb(nzb+1:nzt-1,nys-1:nyn+1), &
-              ipme(nzb+1:nzt-1,nys-1:nyn+1), ippb(nzb+1:nzt-1,nys-1:nyn+1), &
-              ippe(nzb+1:nzt-1,nys-1:nyn+1), m1(nzb+1:nzt-1,nys-2:nyn+2),   &
-              sw(nzb+1:nzt-1,nys-1:nyn+1)                                   &
+    ALLOCATE( a0(nzb+1:nzt,nys-1:nyn+1),   a1(nzb+1:nzt,nys-1:nyn+1),   &
+              a2(nzb+1:nzt,nys-1:nyn+1),   a12(nzb+1:nzt,nys-1:nyn+1),  &
+              a22(nzb+1:nzt,nys-1:nyn+1),  immb(nzb+1:nzt,nys-1:nyn+1), &
+              imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1), &
+              impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1), &
+              ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1), &
+              ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2),   &
+              sw(nzb+1:nzt,nys-1:nyn+1)                                 &
             )
     imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
@@ -475,5 +475,5 @@
 !--    Compute polynomial coefficients
        DO  j = nys-1, nyn+1
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              a12(k,j) = 0.5 * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
              a22(k,j) = 0.5 * ( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) &
@@ -495,5 +495,5 @@
 !--    *VOCL LOOP,UNROLL(2)
        DO  j = nys, nyn
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              cip  =  MAX( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
              cim  = -MIN( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
@@ -531,5 +531,5 @@
 !--    Compute monitor function m1
        DO  j = nys-2, nyn+2
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              m1z = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )
              m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
@@ -549,5 +549,5 @@
        sw = 0.0
        DO  j = nys-1, nyn+1
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &
                   MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 )
@@ -575,5 +575,5 @@
        CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
        DO  j = nys, nyn
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
 
 !--          *VOCL STMT,IF(10)
@@ -665,5 +665,5 @@
 !--    Prognostic equation
        DO  j = nys, nyn
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              fplus  = ( 1.0 - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
                     - ( 1.0 - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
@@ -720,7 +720,6 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
-                sk_p(nzb-1,j,i) = sk_p(nzb+1,j,i)
-                sk_p(nzb-2,j,i) = sk_p(nzb+1,j,i)
+                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
              ENDDO
           ENDDO
@@ -730,22 +729,21 @@
 !
 !--    Temperature boundary condition at the top boundary
-       IF ( ibc_pt_t == 0 )  THEN
-!
-!--       Dirichlet
+       IF ( ibc_pt_t == 0  .OR.  ibc_pt_t == 1 )  THEN
+!
+!--       Dirichlet or Neumann (zero gradient)
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                sk_p(nzt+1,j,i)   = sk_p(nzt,j,i)
-                sk_p(nzt+2,j,i)   = sk_p(nzt,j,i)
+                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
+                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
              ENDDO
           ENDDO
 
-       ELSE
-!
-!--       Neumann: dzu(nzt+2) is not defined, thus instead dzu(nzt+1) is used
+       ELSEIF ( ibc_pt_t == 2 )  THEN
+!
+!--       Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                sk_p(nzt,j,i)     = sk_p(nzt-1,j,i) + bc_pt_t_val * dzu(nzt)
-                sk_p(nzt+1,j,i)   = sk_p(nzt,j,i)   + bc_pt_t_val * dzu(nzt+1)
                 sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_pt_t_val * dzu(nzt+1)
+                sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_pt_t_val * dzu(nzt+1)
              ENDDO
           ENDDO
@@ -756,8 +754,8 @@
 
 !
-!--    Temperature boundary condition at the bottom boundary
+!--    Specific humidity boundary condition at the bottom boundary
        IF ( ibc_q_b == 0 )  THEN
 !
-!--       Dirichlet (fixed surface temperature)
+!--       Dirichlet (fixed surface humidity)
           DO  i = nxl, nxr
              DO  j = nys, nyn
@@ -772,7 +770,6 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
-                sk_p(nzb-1,j,i) = sk_p(nzb+1,j,i)
-                sk_p(nzb-2,j,i) = sk_p(nzb+1,j,i)
+                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
              ENDDO
           ENDDO
@@ -781,5 +778,5 @@
 
 !
-!--    Temperature boundary condition at the top boundary
+!--    Specific humidity boundary condition at the top boundary
        IF ( ibc_q_t == 0 )  THEN
 !
@@ -787,6 +784,6 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                sk_p(nzt+1,j,i)   = sk_p(nzt,j,i)
-                sk_p(nzt+2,j,i)   = sk_p(nzt,j,i)
+                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
+                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
              ENDDO
           ENDDO
@@ -794,10 +791,9 @@
        ELSE
 !
-!--       Neumann: dzu(nzt+2) is not defined, thus instead dzu(nzt+1) is used
+!--       Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                sk_p(nzt,j,i)     = sk_p(nzt-1,j,i) + bc_q_t_val * dzu(nzt)
-                sk_p(nzt+1,j,i)   = sk_p(nzt,j,i)   + bc_q_t_val * dzu(nzt+1)
                 sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i) + bc_q_t_val * dzu(nzt+1)
+                sk_p(nzt+3,j,i)   = sk_p(nzt+2,j,i) + bc_q_t_val * dzu(nzt+1)
              ENDDO
           ENDDO
@@ -809,10 +805,8 @@
 !
 !--    TKE boundary condition at bottom and top boundary (generally Neumann)
-       sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
-       sk_p(nzb-1,j,i) = sk_p(nzb+1,j,i)
-       sk_p(nzb-2,j,i) = sk_p(nzb+1,j,i)
-       sk_p(nzt,j,i)   = sk_p(nzt-1,j,i)
-       sk_p(nzt+1,j,i) = sk_p(nzt-1,j,i)
-       sk_p(nzt+2,j,i) = sk_p(nzt-1,j,i)
+       sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
+       sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
+       sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
+       sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
 
     ELSE
@@ -829,5 +823,5 @@
     DO  i = nxl, nxr
        DO  j = nys, nyn
-          DO  k = nzb, nzt
+          DO  k = nzb, nzt+1
              zaehler = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
              nenner  = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) )
@@ -847,12 +841,12 @@
 !
 !-- Allocate temporary arrays
-    ALLOCATE( a0(nzb:nzt,nys:nyn),       a1(nzb:nzt,nys:nyn),       &
-              a2(nzb:nzt,nys:nyn),       a12(nzb:nzt,nys:nyn),      &
-              a22(nzb:nzt,nys:nyn),      immb(nzb+1:nzt-1,nys:nyn), &
-              imme(nzb+1:nzt-1,nys:nyn), impb(nzb+1:nzt-1,nys:nyn), &
-              impe(nzb+1:nzt-1,nys:nyn), ipmb(nzb+1:nzt-1,nys:nyn), &
-              ipme(nzb+1:nzt-1,nys:nyn), ippb(nzb+1:nzt-1,nys:nyn), &
-              ippe(nzb+1:nzt-1,nys:nyn), m1(nzb-1:nzt+1,nys:nyn),   &
-              sw(nzb:nzt,nys:nyn)                                   &
+    ALLOCATE( a0(nzb:nzt+1,nys:nyn),   a1(nzb:nzt+1,nys:nyn),   &
+              a2(nzb:nzt+1,nys:nyn),   a12(nzb:nzt+1,nys:nyn),  &
+              a22(nzb:nzt+1,nys:nyn),  immb(nzb+1:nzt,nys:nyn), &
+              imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn), &
+              impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn), &
+              ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn), &
+              ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn), &
+              sw(nzb:nzt+1,nys:nyn)                             &
             )
     imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
@@ -865,5 +859,5 @@
 !--    Compute polynomial coefficients
        DO  j = nys, nyn
-          DO  k = nzb, nzt
+          DO  k = nzb, nzt+1
              a12(k,j) = 0.5 * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
              a22(k,j) = 0.5 * ( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) &
@@ -885,5 +879,5 @@
 !--    *VOCL LOOP,UNROLL(2)
        DO  j = nys, nyn
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              cip  =  MAX( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
              cim  = -MIN( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
@@ -921,5 +915,5 @@
 !--    Compute monitor function m1
        DO  j = nys, nyn
-          DO  k = nzb-1, nzt+1
+          DO  k = nzb-1, nzt+2
              m1z = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
              m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
@@ -939,5 +933,5 @@
        sw = 0.0
        DO  j = nys, nyn
-          DO  k = nzb, nzt
+          DO  k = nzb, nzt+1
              m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &
                   MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 )
@@ -965,5 +959,5 @@
        CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
        DO  j = nys, nyn
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
 
 !--          *VOCL STMT,IF(10)
@@ -1055,5 +1049,5 @@
 !--    Prognostic equation
        DO  j = nys, nyn
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              fplus  = ( 1.0 - sw(k,j)   ) * ippb(k,j) + sw(k,j)   * ippe(k,j) &
                     - ( 1.0 - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
@@ -1082,5 +1076,5 @@
           DO  sr = 0, statistic_regions
              DO  j = nys, nyn
-                DO  k = nzb+1, nzt-1
+                DO  k = nzb+1, nzt
                    sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) + &
                                           m1(k,j) * rmask(j,i,sr)
@@ -1103,5 +1097,5 @@
     DO  i = nxl, nxr
        DO  j = nys, nyn
-          DO  k = nzb+1, nzt-1
+          DO  k = nzb+1, nzt
              tend(k,j,i) = tend(k,j,i) + ( sk_p(k,j,i) - sk(k,j,i) ) / dt_3d
           ENDDO
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 62)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 63)
@@ -5,6 +5,6 @@
 ! -----------------
 ! "by_user" allowed as initializing action, -data_output_ts,
-! leapfrog with non-flat topography not allowed any more, pt_reference is
-! checked
+! leapfrog with non-flat topography not allowed any more, loop_optimization
+! and pt_reference are checked
 !
 ! Former revisions:
@@ -82,4 +82,23 @@
 
 !
+!-- Check the general loop optimization method
+    IF ( loop_optimization == 'default' )  THEN
+       IF ( host(1:3) == 'nec' )  THEN
+          loop_optimization = 'vector'
+       ELSE
+          loop_optimization = 'cache'
+       ENDIF
+    ENDIF
+    IF ( loop_optimization /= 'noopt'  .AND.  loop_optimization /= 'cache' &
+         .AND.  loop_optimization /= 'vector' )  THEN
+       IF ( myid == 0 )  THEN
+          PRINT*, '+++ check_parameters:'
+          PRINT*, '    illegal value given for loop_optimization: ', &
+                  TRIM( loop_optimization )
+       ENDIF
+       CALL local_stop
+    ENDIF
+
+!
 !-- Check topography setting (check for illegal parameter combinations)
     IF ( topography /= 'flat' )  THEN
@@ -257,5 +276,5 @@
     END SELECT
 
-    IF ( scalar_advec /= 'pw-scheme'  .AND.  timestep_scheme(1:5) == 'runge' ) &
+    IF ( scalar_advec == 'ups-scheme'  .AND.  timestep_scheme(1:5) == 'runge' )&
     THEN
        IF ( myid == 0 )  THEN
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 62)
+++ /palm/trunk/SOURCE/header.f90	(revision 63)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! Output of netcdf_64bit_3d, particles-package is now part of the default code.
+! Output of netcdf_64bit_3d, particles-package is now part of the default code,
+! output of the loop optimization method.
 !
 ! Former revisions:
@@ -189,8 +190,7 @@
        WRITE ( io, 118 )
     ENDIF
-    IF ( momentum_advec /= 'ups-scheme' .AND. scalar_advec /= 'ups-scheme' &
-         .AND. scalar_advec /= 'bc-scheme'  .AND.  host(1:3) /= 'nec' )  THEN
-       WRITE ( io, 139 )
-    ENDIF
+
+    WRITE ( io, 139 )  TRIM( loop_optimization )
+
     IF ( galilei_transformation )  THEN
        IF ( use_ug_for_galilei_tr )  THEN
@@ -1182,6 +1182,5 @@
                   I3,')')
 138 FORMAT ('     Using hybrid version for 1d-domain-decomposition')
-139 FORMAT (' --> Prognostic equations are solved in one single loop (fast', &
-                  ' method)')
+139 FORMAT (' --> Loop optimization method: ',A)
 140 FORMAT ('     maximum residual allowed:                ',E10.3)
 141 FORMAT ('     fixed number of multigrid cycles:        ',I4)
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 62)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 63)
@@ -756,5 +756,5 @@
 !
 !-- If required, initialize particles
-    CALL init_particles
+    IF ( particle_advection )  CALL init_particles
 
 !
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 62)
+++ /palm/trunk/SOURCE/modules.f90	(revision 63)
@@ -6,6 +6,7 @@
 ! -----------------
 ! +array rif_wall
-! +netcdf_64bit_3d, zu_s_inner, zw_w_inner, id_var_zusi_*, id_var_zwwi_*,
-! ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc, pt_reference, use_pt_reference,
+! +loop_optimization, netcdf_64bit_3d, zu_s_inner, zw_w_inner, id_var_zusi_*,
+! id_var_zwwi_*, ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc, pt_reference,
+! use_pt_reference,
 ! +age_m in particle_type
 ! -data_output_ts, dots_n
@@ -211,5 +212,6 @@
     CHARACTER (LEN=9)   ::  simulated_time_chr
     CHARACTER (LEN=12)  ::  version = 'PALM   3.1c'
-    CHARACTER (LEN=16)  ::  momentum_advec = 'pw-scheme', &
+    CHARACTER (LEN=16)  ::  loop_optimization = 'default', &
+                            momentum_advec = 'pw-scheme', &
                             psolver = 'poisfft', &
                             scalar_advec = 'pw-scheme'
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 62)
+++ /palm/trunk/SOURCE/parin.f90	(revision 63)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! +dt_max, netcdf_64bit_3d in d3par, +pt_reference in inipar, -data_output_ts
+! +dt_max, netcdf_64bit_3d in d3par, +loop_optimization, pt_reference in
+! inipar, -data_output_ts
 !
 ! Former revisions:
@@ -60,7 +61,8 @@
                        inflow_disturbance_end, initializing_actions, &
                        km_constant, km_damp_max, long_filter_factor, &
-                       mixing_length_1d, moisture, momentum_advec, &
-                       netcdf_precision, npex, npey, nsor_ini, nx, ny, nz, &
-                       omega, outflow_damping_width, overshoot_limit_e, &
+                       loop_optimization, mixing_length_1d, moisture, &
+                       momentum_advec, netcdf_precision, npex, npey, nsor_ini, &
+                       nx, ny, nz, omega, outflow_damping_width, &
+                       overshoot_limit_e, &
                        overshoot_limit_pt, overshoot_limit_u, &
                        overshoot_limit_v, overshoot_limit_w, passive_scalar, &
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 62)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 63)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! z0 removed from arguments in calls of diffusion_u/v/w
+! z0 removed from arguments in calls of diffusion_u/v/w,
+! subroutine names changed to .._noopt, .._cache, and .._vector
 !
 ! Former revisions:
@@ -63,18 +64,18 @@
 
     PRIVATE
-    PUBLIC prognostic_equations, prognostic_equations_fast, &
-           prognostic_equations_vec
-
-    INTERFACE prognostic_equations
-       MODULE PROCEDURE prognostic_equations
-    END INTERFACE prognostic_equations
-
-    INTERFACE prognostic_equations_fast
-       MODULE PROCEDURE prognostic_equations_fast
-    END INTERFACE prognostic_equations_fast
-
-    INTERFACE prognostic_equations_vec
-       MODULE PROCEDURE prognostic_equations_vec
-    END INTERFACE prognostic_equations_vec
+    PUBLIC prognostic_equations_noopt, prognostic_equations_cache, &
+           prognostic_equations_vector
+
+    INTERFACE prognostic_equations_noopt
+       MODULE PROCEDURE prognostic_equations_noopt
+    END INTERFACE prognostic_equations_noopt
+
+    INTERFACE prognostic_equations_cache
+       MODULE PROCEDURE prognostic_equations_cache
+    END INTERFACE prognostic_equations_cache
+
+    INTERFACE prognostic_equations_vector
+       MODULE PROCEDURE prognostic_equations_vector
+    END INTERFACE prognostic_equations_vector
 
 
@@ -82,5 +83,5 @@
 
 
- SUBROUTINE prognostic_equations
+ SUBROUTINE prognostic_equations_noopt
 
 !------------------------------------------------------------------------------!
@@ -615,8 +616,8 @@
 
 
- END SUBROUTINE prognostic_equations
-
-
- SUBROUTINE prognostic_equations_fast
+ END SUBROUTINE prognostic_equations_noopt
+
+
+ SUBROUTINE prognostic_equations_cache
 
 !------------------------------------------------------------------------------!
@@ -996,8 +997,8 @@
 
 
- END SUBROUTINE prognostic_equations_fast
-
-
- SUBROUTINE prognostic_equations_vec
+ END SUBROUTINE prognostic_equations_cache
+
+
+ SUBROUTINE prognostic_equations_vector
 
 !------------------------------------------------------------------------------!
@@ -1242,14 +1243,18 @@
 !
 !-- pt-tendency terms with communication
+    sat = tsc(1)
+    sbt = tsc(2)
     IF ( scalar_advec == 'bc-scheme' )  THEN
-!
-!--    Bott-Chlond scheme always uses Euler time step. Thus:
-       sat = 1.0
-       sbt = 1.0
+
+       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
-       sat = tsc(1)
-       sbt = tsc(2)
        IF ( tsc(2) /= 2.0  .AND.  scalar_advec == 'ups-scheme' )  THEN
           tend = 0.0
@@ -1340,14 +1345,18 @@
 !
 !--    Scalar/q-tendency terms with communication
+       sat = tsc(1)
+       sbt = tsc(2)
        IF ( scalar_advec == 'bc-scheme' )  THEN
-!
-!--       Bott-Chlond scheme always uses Euler time step. Thus:
-          sat = 1.0
-          sbt = 1.0
+
+          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( q, 'q' )
        ELSE
-          sat = tsc(1)
-          sbt = tsc(2)
           IF ( tsc(2) /= 2.0 )  THEN
              IF ( scalar_advec == 'ups-scheme' )  THEN
@@ -1438,15 +1447,20 @@
 !--    TKE-tendency terms with communication
        CALL production_e_init
+
+       sat = tsc(1)
+       sbt = tsc(2)
        IF ( .NOT. use_upstream_for_tke )  THEN
           IF ( scalar_advec == 'bc-scheme' )  THEN
-!
-!--          Bott-Chlond scheme always uses Euler time step. Thus:
-             sat = 1.0
-             sbt = 1.0
+
+             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( e, 'e' )
           ELSE
-             sat = tsc(1)
-             sbt = tsc(2)
              IF ( tsc(2) /= 2.0 )  THEN
                 IF ( scalar_advec == 'ups-scheme' )  THEN
@@ -1550,5 +1564,5 @@
 
 
- END SUBROUTINE prognostic_equations_vec
+ END SUBROUTINE prognostic_equations_vector
 
 
Index: /palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/read_var_list.f90	(revision 62)
+++ /palm/trunk/SOURCE/read_var_list.f90	(revision 63)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! +pt_reference
+! +loop_optimization, pt_reference
 !
 ! Former revisions:
@@ -208,4 +208,6 @@
           CASE ( 'long_filter_factor' )
              READ ( 13 )  long_filter_factor
+          CASE ( 'loop_optimization' )
+             READ ( 13 )  loop_optimization
           CASE ( 'mixing_length_1d' )
              READ ( 13 )  mixing_length_1d
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 62)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 63)
@@ -5,5 +5,7 @@
 ! -----------------
 ! Move call of user_actions( 'after_integration' ) below increment of times
-! and counters
+! and counters,
+! calls of prognostic_equations_.. changed to .._noopt, .._cache, and
+! .._vector, these calls are now controlled by switch loop_optimization
 !
 ! Former revisions:
@@ -98,6 +100,6 @@
 !--       in the other versions a good vectorization is prohibited due to
 !--       inlining problems.
-          IF ( host(1:3) == 'nec' )  THEN
-             CALL prognostic_equations_vec
+          IF ( loop_optimization == 'vector' )  THEN
+             CALL prognostic_equations_vector
           ELSE
              IF ( momentum_advec == 'ups-scheme'  .OR.  &
@@ -105,7 +107,7 @@
                   scalar_advec == 'bc-scheme' )        &
              THEN
-                CALL prognostic_equations
+                CALL prognostic_equations_noopt
              ELSE
-                CALL prognostic_equations_fast
+                CALL prognostic_equations_cache
              ENDIF
           ENDIF
@@ -114,5 +116,6 @@
 !--       Particle advection (only once during intermediate steps, because
 !--       it uses an Euler-step)
-          IF ( simulated_time >= particle_advection_start  .AND. &
+          IF ( particle_advection  .AND.                         &
+               simulated_time >= particle_advection_start  .AND. &
                intermediate_timestep_count == 1 )  THEN
              CALL advec_particles
Index: /palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/write_var_list.f90	(revision 62)
+++ /palm/trunk/SOURCE/write_var_list.f90	(revision 63)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! +pt_refrence
+! +loop_optimization, pt_refrence
 !
 ! Former revisions:
@@ -174,4 +174,6 @@
     WRITE ( 14 )  'long_filter_factor            '
     WRITE ( 14 )  long_filter_factor
+    WRITE ( 14 )  'loop_optimization             '
+    WRITE ( 14 )  loop_optimization
     WRITE ( 14 )  'mixing_length_1d              '
     WRITE ( 14 )  mixing_length_1d
