Index: palm/trunk/SOURCE/advec_s_bc.f90
===================================================================
--- palm/trunk/SOURCE/advec_s_bc.f90	(revision 4)
+++ 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
