Index: palm/trunk/SOURCE/diffusion_e.f90
===================================================================
--- palm/trunk/SOURCE/diffusion_e.f90	(revision 2119)
+++ palm/trunk/SOURCE/diffusion_e.f90	(revision 2232)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adjustments to new topography and surface concept
 ! 
 ! Former revisions:
@@ -129,5 +129,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,        &
+           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, l_wall, tend, zu, zw,&
                   drho_air, rho_air_zw
            
@@ -140,6 +140,6 @@
            
        USE indices,                                                            &
-           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_s_inner,&
-                  nzt
+           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max,    &
+                  nzt, wall_flags_0
            
        USE kinds
@@ -151,10 +151,16 @@
            ONLY:  use_sgs_for_particles, wang_kernel
 
+       USE surface_mod,                                                        &
+          ONLY :  bc_h
+
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i              !< 
-       INTEGER(iwp) ::  j              !< 
-       INTEGER(iwp) ::  k              !< 
+       INTEGER(iwp) ::  i              !< running index x direction
+       INTEGER(iwp) ::  j              !< running index y direction
+       INTEGER(iwp) ::  k              !< running index z direction
+       INTEGER(iwp) ::  m              !< running index surface elements
+ 
        REAL(wp)     ::  dvar_dz        !< 
+       REAL(wp)     ::  flag           !< flag to mask topography
        REAL(wp)     ::  l_stable       !< 
        REAL(wp)     ::  var_reference  !< 
@@ -177,5 +183,5 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt
+                DO  k = nzb+1, nzt
 !
 !--                Calculate the mixing length (for dissipation)
@@ -191,9 +197,7 @@
 !--                Adjustment of the mixing length
                    IF ( wall_adjustment )  THEN
-                      l(k,j)  = MIN( wall_adjustment_factor *          &
-                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+                      l(k,j)  = MIN( wall_adjustment_factor * l_wall(k,j,i),   &
                                      l_grid(k), l_stable )
-                      ll(k,j) = MIN( wall_adjustment_factor *          &
-                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+                      ll(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i),   &
                                      l_grid(k) )
                    ELSE
@@ -208,18 +212,21 @@
 !--          Calculate the tendency terms
              DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt
-
-                    dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
+                DO  k = nzb+1, nzt
+!
+!--                Predetermine flag to mask topography
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
+                   dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
                                        e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
 
-                    tend(k,j,i) = tend(k,j,i)                                  &
+                   tend(k,j,i) = tend(k,j,i)                                   &
                                         + (                                    &
                           ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
                         - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
-                                          ) * ddx2                             &
+                                          ) * ddx2  * flag                     &
                                         + (                                    &
                           ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
                         - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
-                                          ) * ddy2                             &
+                                          ) * ddy2  * flag                     &
                                         + (                                    &
                ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
@@ -227,6 +234,6 @@
              - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
                                                              * rho_air_zw(k-1) &
-                                          ) * ddzw(k) * drho_air(k)            &
-                             - dissipation(k,j)
+                                          ) * ddzw(k) * drho_air(k) * flag     &
+                             - dissipation(k,j) * flag
 
                 ENDDO
@@ -239,6 +246,8 @@
                   collision_turbulence )  THEN
                 DO  j = nys, nyn
-                   DO  k = nzb_s_inner(j,i)+1, nzt
-                      diss(k,j,i) = dissipation(k,j)
+                   DO  k = nzb+1, nzt
+                      diss(k,j,i) = dissipation(k,j) *                         &
+                                        MERGE( 1.0_wp, 0.0_wp,                 &
+                                               BTEST( wall_flags_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -251,5 +260,5 @@
           DO  i = nxl, nxr
              DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt
+                DO  k = nzb+1, nzt
 !
 !--                Calculate the mixing length (for dissipation)
@@ -265,9 +274,7 @@
 !--                Adjustment of the mixing length
                    IF ( wall_adjustment )  THEN
-                      l(k,j)  = MIN( wall_adjustment_factor *          &
-                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+                      l(k,j)  = MIN( wall_adjustment_factor * l_wall(k,j,i),   &
                                      l_grid(k), l_stable )
-                      ll(k,j) = MIN( wall_adjustment_factor *          &
-                                     ( zu(k) - zw(nzb_s_inner(j,i)) ), &
+                      ll(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i),   &
                                      l_grid(k) )
                    ELSE
@@ -282,18 +289,21 @@
 !--          Calculate the tendency terms
              DO  j = nys, nyn
-                DO  k = nzb_s_inner(j,i)+1, nzt
-
-                    dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
+                DO  k = nzb+1, nzt
+!
+!--                Predetermine flag to mask topography
+                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
+                   dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) )   *&
                                              e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
 
-                    tend(k,j,i) = tend(k,j,i)                                  &
+                   tend(k,j,i) = tend(k,j,i)                                   &
                                         + (                                    &
                           ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
                         - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
-                                          ) * ddx2                             &
+                                          ) * ddx2  * flag                     &
                                         + (                                    &
                           ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
                         - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
-                                          ) * ddy2                             &
+                                          ) * ddy2  * flag                     &
                                         + (                                    &
                ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
@@ -301,6 +311,6 @@
              - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
                                                              * rho_air_zw(k-1) &
-                                          ) * ddzw(k) * drho_air(k)            &
-                             - dissipation(k,j)
+                                          ) * ddzw(k) * drho_air(k) * flag     &
+                             - dissipation(k,j) * flag
 
                 ENDDO
@@ -313,6 +323,8 @@
                   collision_turbulence )  THEN
                 DO  j = nys, nyn
-                   DO  k = nzb_s_inner(j,i)+1, nzt
-                      diss(k,j,i) = dissipation(k,j)
+                   DO  k = nzb+1, nzt
+                      diss(k,j,i) = dissipation(k,j) *                         &
+                                        MERGE( 1.0_wp, 0.0_wp,                 &
+                                               BTEST( wall_flags_0(k,j,i), 0 ) )
                    ENDDO
                 ENDDO
@@ -324,12 +336,24 @@
 
 !
-!--    Boundary condition for dissipation
+!--    Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
        IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  &
             collision_turbulence )  THEN
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
-             ENDDO
-          ENDDO
+!
+!--       Upward facing surfaces
+          DO  m = 1, bc_h(0)%ns
+             i = bc_h(0)%i(m)            
+             j = bc_h(0)%j(m)
+             k = bc_h(0)%k(m)
+             diss(k-1,j,i) = diss(k,j,i)
+          ENDDO
+!
+!--       Downward facing surfaces
+          DO  m = 1, bc_h(1)%ns
+             i = bc_h(1)%i(m)            
+             j = bc_h(1)%j(m)
+             k = bc_h(1)%k(m)
+             diss(k+1,j,i) = diss(k,j,i)
+          ENDDO
+
        ENDIF
 
@@ -345,5 +369,5 @@
 
        USE arrays_3d,                                                          &
-           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw,        &
+           ONLY:  dd2zu, ddzu, ddzw, diss, e, km, l_grid, l_wall, tend, zu, zw,&
                   drho_air, rho_air_zw
           
@@ -356,5 +380,5 @@
            
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
            
        USE kinds
@@ -366,10 +390,18 @@
            ONLY:  use_sgs_for_particles, wang_kernel
 
+       USE surface_mod,                                                        &
+          ONLY :  bc_h
+
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i              !< 
-       INTEGER(iwp) ::  j              !< 
-       INTEGER(iwp) ::  k              !< 
+       INTEGER(iwp) ::  i              !< running index x direction
+       INTEGER(iwp) ::  j              !< running index y direction
+       INTEGER(iwp) ::  k              !< running index z direction
+       INTEGER(iwp) ::  m              !< running index surface elements
+       INTEGER(iwp) ::  surf_e         !< End index of surface elements at (j,i)-gridpoint
+       INTEGER(iwp) ::  surf_s         !< Start index of surface elements at (j,i)-gridpoint
+
        REAL(wp)     ::  dvar_dz        !< 
+       REAL(wp)     ::  flag           !< flag to mask topography
        REAL(wp)     ::  l_stable       !<
        REAL(wp)     ::  var_reference  !<
@@ -384,8 +416,11 @@
        REAL(wp), DIMENSION(nzb+1:nzt) ::  ll           !<
 
-
 !
 !--    Calculate the mixing length (for dissipation)
-       DO  k = nzb_s_inner(j,i)+1, nzt
+       DO  k = nzb+1, nzt
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
           dvar_dz = atmos_ocean_sign * &
                     ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
@@ -404,9 +439,7 @@
 !--       Adjustment of the mixing length
           IF ( wall_adjustment )  THEN
-             l(k)  = MIN( wall_adjustment_factor *                     &
-                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
-                          l_stable )
-             ll(k) = MIN( wall_adjustment_factor *                     &
-                          ( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
+             l(k)  = MIN( wall_adjustment_factor * l_wall(k,j,i),              &
+                          l_grid(k), l_stable )
+             ll(k) = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
           ELSE
              l(k)  = MIN( l_grid(k), l_stable )
@@ -415,23 +448,23 @@
 !
 !--       Calculate the tendency term
-          dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * &
+          dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) *   &
                                  SQRT( e(k,j,i) ) / l(k)
 
-          tend(k,j,i) = tend(k,j,i)                                           &
-                                       + (                                    &
-                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
-                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
-                                         ) * ddx2                             &
-                                       + (                                    &
-                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
-                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
-                                         ) * ddy2                             &
-                                       + (                                    &
-              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
-                                                            * rho_air_zw(k)   &
-            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
-                                                            * rho_air_zw(k-1) &
-                                         ) * ddzw(k) * drho_air(k)            &
-                                       - dissipation(k)
+          tend(k,j,i) = tend(k,j,i)                                            &
+                                       + (                                     &
+                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )   &
+                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )   &
+                                         ) * ddx2  * flag                      &
+                                       + (                                     &
+                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )   &
+                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )   &
+                                         ) * ddy2  * flag                      &
+                                       + (                                     &
+              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1)  &
+                                                            * rho_air_zw(k)    &
+            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)    &
+                                                            * rho_air_zw(k-1)  &
+                                         ) * ddzw(k) * drho_air(k) * flag      &
+                                       - dissipation(k) * flag
 
        ENDDO
@@ -441,10 +474,27 @@
        IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.                     &
             collision_turbulence )  THEN
-          DO  k = nzb_s_inner(j,i)+1, nzt
-             diss(k,j,i) = dissipation(k)
-          ENDDO
-!
-!--       Boundary condition for dissipation
-          diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
+          DO  k = nzb+1, nzt
+             diss(k,j,i) = dissipation(k) *                                    &
+                                        MERGE( 1.0_wp, 0.0_wp,                 &
+                                               BTEST( wall_flags_0(k,j,i), 0 ) )
+          ENDDO
+!
+!--       Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
+!--       For each surface type determine start and end index (in case of elevated
+!--       toopography several up/downward facing surfaces may exist.
+          surf_s = bc_h(0)%start_index(j,i)   
+          surf_e = bc_h(0)%end_index(j,i)   
+          DO  m = surf_s, surf_e
+             k             = bc_h(0)%k(m)
+             diss(k-1,j,i) = diss(k,j,i)
+          ENDDO
+!
+!--       Downward facing surfaces
+          surf_s = bc_h(1)%start_index(j,i)   
+          surf_e = bc_h(1)%end_index(j,i)   
+          DO  m = surf_s, surf_e
+             k             = bc_h(1)%k(m)
+             diss(k+1,j,i) = diss(k,j,i)
+          ENDDO
        ENDIF
 
