Index: palm/trunk/SOURCE/prandtl_fluxes.f90
===================================================================
--- palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 484)
+++ palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 667)
@@ -4,5 +4,9 @@
 ! Current revisions:
 ! -----------------
-! 
+!
+! Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
+! therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
+!
+! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
 !
 ! Former revisions:
@@ -53,5 +57,5 @@
     REAL    ::  a, b, e_q, rifm, uv_total, z_p
 
-!
+! 
 !-- Compute theta*
     IF ( constant_heatflux )  THEN
@@ -60,6 +64,6 @@
 !--    for u* use the value from the previous time step
        !$OMP PARALLEL DO
-       DO  i = nxl-1, nxr+1
-          DO  j = nys-1, nyn+1
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              ts(j,i) = -shf(j,i) / ( us(j,i) + 1E-30 )
 !
@@ -76,6 +80,6 @@
 !--    (the Richardson number is still the one from the previous time step)
        !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
-       DO  i = nxl-1, nxr+1
-          DO  j = nys-1, nyn+1
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
 
              k   = nzb_s_inner(j,i)
@@ -108,6 +112,6 @@
     IF ( .NOT. humidity )  THEN
        !$OMP PARALLEL DO PRIVATE( k, z_p )
-       DO  i = nxl-1, nxr+1
-          DO  j = nys-1, nyn+1
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              k   = nzb_s_inner(j,i)
              z_p = zu(k+1) - zw(k)
@@ -126,6 +130,6 @@
     ELSE
        !$OMP PARALLEL DO PRIVATE( k, z_p )
-       DO  i = nxl-1, nxr+1
-          DO  j = nys-1, nyn+1
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              k   = nzb_s_inner(j,i)
              z_p = zu(k+1) - zw(k)
@@ -155,8 +159,11 @@
 
 !
-!--       Compute the absolute value of the horizontal velocity
-          uv_total = SQRT( ( 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) ) )**2 + &
-                           ( 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) ) )**2   &
-                         )
+!--       Compute the absolute value of the horizontal velocity 
+!--       (relative to the surface)
+          uv_total = SQRT( ( 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1)        &
+                                   - u(k,j,i)   - u(k,j,i+1) ) )**2 + &
+                           ( 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i)        &
+                                   - v(k,j,i)   - v(k,j+1,i) ) )**2 )    
+
 
           IF ( rif(j,i) >= 0.0 )  THEN
@@ -203,5 +210,5 @@
 !
 !--          Stable stratification
-             usws(j,i) = kappa * u(k+1,j,i) / (                           &
+             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )/ (              &
                                      LOG( z_p / z0(j,i) ) +               &
                                      5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
@@ -213,5 +220,5 @@
              b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
 
-             usws(j,i) = kappa * u(k+1,j,i) / (                           &
+             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) ) / (            &
                          LOG( z_p / z0(j,i) ) -                           &
                          LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
@@ -240,5 +247,5 @@
 !
 !--          Stable stratification
-             vsws(j,i) = kappa * v(k+1,j,i) / (                           &
+             vsws(j,i) = kappa * ( v(k+1,j,i) -  v(k,j,i) ) / (           &
                                      LOG( z_p / z0(j,i) ) +               &
                                      5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
@@ -250,5 +257,5 @@
              b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
 
-             vsws(j,i) = kappa * v(k+1,j,i) / (                           &
+             vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / (            &
                          LOG( z_p / z0(j,i) ) -                           &
                          LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
@@ -268,6 +275,6 @@
 !--       For a given water flux in the Prandtl layer:
           !$OMP PARALLEL DO
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
                 qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30 )
              ENDDO
@@ -276,6 +283,6 @@
        ELSE          
           !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
-          DO  i = nxl-1, nxr+1
-             DO  j = nys-1, nyn+1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
 
                 k   = nzb_s_inner(j,i)
@@ -325,6 +332,6 @@
     IF ( .NOT. constant_heatflux )  THEN
        !$OMP PARALLEL DO
-       DO  i = nxl-1, nxr+1
-          DO  j = nys-1, nyn+1
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              shf(j,i) = -ts(j,i) * us(j,i)
           ENDDO
@@ -336,6 +343,6 @@
     IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar ) ) THEN
        !$OMP PARALLEL DO
-       DO  i = nxl-1, nxr+1
-          DO  j = nys-1, nyn+1
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              qsws(j,i) = -qs(j,i) * us(j,i)
           ENDDO
@@ -347,6 +354,6 @@
     IF ( ibc_e_b == 2 )  THEN
        !$OMP PARALLEL DO
-       DO  i = nxl-1, nxr+1
-          DO  j = nys-1, nyn+1
+       DO  i = nxlg, nxrg
+          DO  j = nysg, nyng
              e(nzb_s_inner(j,i)+1,j,i) = ( us(j,i) / 0.1 )**2
 !
