Index: palm/trunk/SOURCE/plant_canopy_model.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model.f90	(revision 1683)
+++ palm/trunk/SOURCE/plant_canopy_model.f90	(revision 1721)
@@ -19,5 +19,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! bugfixes: shf is reduced in areas covered with canopy only,
+!           canopy is set on top of topography
 ! 
 ! Former revisions:
@@ -357,7 +358,14 @@
 
 !
-!--       The surface heat flux is set to the surface value of the calculated
-!--       in-canopy heat flux distribution  
-          shf(:,:) = canopy_heat_flux(0,:,:)
+!--       In areas covered with canopy, the surface heat flux is set to
+!--       the surface value of the above calculated in-canopy heat flux
+!--       distribution
+          DO  i = nxlg,nxrg
+             DO  j = nysg, nyng
+                IF ( canopy_heat_flux(0,j,i) /= cthf )  THEN
+                   shf(j,i) = canopy_heat_flux(0,j,i)
+                ENDIF
+             ENDDO
+          ENDDO
 
        ENDIF
@@ -407,4 +415,5 @@
        INTEGER(iwp) ::  j         !< running index
        INTEGER(iwp) ::  k         !< running index
+       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
 
        REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
@@ -427,6 +436,7 @@
              DO  i = nxlu, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb_u_inner(j,i)+1, pch_index
-
+                   DO  k = nzb_u_inner(j,i)+1, nzb_u_inner(j,i)+pch_index
+
+                      kk = k - nzb_u_inner(j,i)  !- lad arrays are defined flat
 !
 !--                   In order to create sharp boundaries of the plant canopy, 
@@ -437,8 +447,8 @@
 !--                   For the same reason, the lad at the rightmost(i+1)canopy 
 !--                   boundary on the u-grid equals lad_s(k,j,i).
-                      lad_local = lad_s(k,j,i)
-                      IF ( lad_local == 0.0_wp  .AND.                          &
-                           lad_s(k,j,i-1) > 0.0_wp )  THEN
-                         lad_local = lad_s(k,j,i-1)
+                      lad_local = lad_s(kk,j,i)
+                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )&
+                      THEN
+                         lad_local = lad_s(kk,j,i-1)
                       ENDIF
 
@@ -487,6 +497,7 @@
              DO  i = nxl, nxr
                 DO  j = nysv, nyn
-                   DO  k = nzb_v_inner(j,i)+1, pch_index
-
+                   DO  k = nzb_v_inner(j,i)+1, nzb_v_inner(j,i)+pch_index
+
+                      kk = k - nzb_v_inner(j,i)  !- lad arrays are defined flat
 !
 !--                   In order to create sharp boundaries of the plant canopy, 
@@ -497,8 +508,8 @@
 !--                   For the same reason, the lad at the northmost(j+1) canopy 
 !--                   boundary on the v-grid equals lad_s(k,j,i).
-                      lad_local = lad_s(k,j,i)
-                      IF ( lad_local == 0.0_wp  .AND.                          &
-                           lad_s(k,j-1,i) > 0.0_wp )  THEN
-                         lad_local = lad_s(k,j-1,i)
+                      lad_local = lad_s(kk,j,i)
+                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )&
+                      THEN
+                         lad_local = lad_s(kk,j-1,i)
                       ENDIF
 
@@ -547,5 +558,7 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb_w_inner(j,i)+1, pch_index-1
+                   DO  k = nzb_w_inner(j,i)+1, nzb_w_inner(j,i)+pch_index-1
+
+                      kk = k - nzb_w_inner(j,i)  !- lad arrays are defined flat
 
                       pre_tend = 0.0_wp
@@ -555,5 +568,5 @@
                       pre_tend = - cdc *                                       &
                                    (0.5_wp *                                   &
-                                      ( lad_s(k+1,j,i) + lad_s(k,j,i) )) *     &
+                                      ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *   &
                                    SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
                                                        u(k,j,i+1) +            &
@@ -593,8 +606,9 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb_s_inner(j,i)+1, pch_index
+                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
+                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
                       tend(k,j,i) = tend(k,j,i) +                              &
-                                       ( canopy_heat_flux(k,j,i) -             &
-                                         canopy_heat_flux(k-1,j,i) ) / dzw(k)
+                                       ( canopy_heat_flux(kk,j,i) -            &
+                                         canopy_heat_flux(kk-1,j,i) ) / dzw(k)
                    ENDDO
                 ENDDO
@@ -606,8 +620,9 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb_s_inner(j,i)+1, pch_index
+                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
+                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
                       tend(k,j,i) = tend(k,j,i) -                              &
                                        lsec *                                  &
-                                       lad_s(k,j,i) *                          &
+                                       lad_s(kk,j,i) *                         &
                                        SQRT( ( 0.5_wp * ( u(k,j,i) +           &
                                                           u(k,j,i+1) )         &
@@ -630,8 +645,9 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb_s_inner(j,i)+1, pch_index
+                   DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
+                      kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
                       tend(k,j,i) = tend(k,j,i) -                              &
                                        2.0_wp * cdc *                          &
-                                       lad_s(k,j,i) *                          &
+                                       lad_s(kk,j,i) *                         &
                                        SQRT( ( 0.5_wp * ( u(k,j,i) +           &
                                                           u(k,j,i+1) )         &
@@ -698,4 +714,5 @@
        INTEGER(iwp) ::  j         !< running index
        INTEGER(iwp) ::  k         !< running index
+       INTEGER(iwp) ::  kk        !< running index for flat lad arrays
 
        REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
@@ -716,6 +733,7 @@
 !--       u-component
           CASE ( 1 )
-             DO  k = nzb_u_inner(j,i)+1, pch_index
-
+             DO  k = nzb_u_inner(j,i)+1, nzb_u_inner(j,i)+pch_index
+
+                kk = k - nzb_u_inner(j,i)  !- lad arrays are defined flat
 !
 !--             In order to create sharp boundaries of the plant canopy, 
@@ -726,8 +744,7 @@
 !--             For the same reason, the lad at the rightmost(i+1)canopy 
 !--             boundary on the u-grid equals lad_s(k,j,i).
-                lad_local = lad_s(k,j,i)
-                IF ( lad_local == 0.0_wp  .AND.                                &
-                     lad_s(k,j,i-1) > 0.0_wp )  THEN
-                   lad_local = lad_s(k,j,i-1)
+                lad_local = lad_s(kk,j,i)
+                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )  THEN
+                   lad_local = lad_s(kk,j,i-1)
                 ENDIF
 
@@ -772,6 +789,7 @@
 !--       v-component
           CASE ( 2 )
-             DO  k = nzb_v_inner(j,i)+1, pch_index
-
+             DO  k = nzb_v_inner(j,i)+1, nzb_v_inner(j,i)+pch_index
+
+                kk = k - nzb_v_inner(j,i)  !- lad arrays are defined flat
 !
 !--             In order to create sharp boundaries of the plant canopy, 
@@ -782,8 +800,7 @@
 !--             For the same reason, the lad at the northmost(j+1)canopy 
 !--             boundary on the v-grid equals lad_s(k,j,i).
-                lad_local = lad_s(k,j,i)
-                IF ( lad_local == 0.0_wp  .AND.                                &
-                     lad_s(k,j-1,i) > 0.0_wp )  THEN
-                   lad_local = lad_s(k,j-1,i)
+                lad_local = lad_s(kk,j,i)
+                IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )  THEN
+                   lad_local = lad_s(kk,j-1,i)
                 ENDIF
 
@@ -828,5 +845,7 @@
 !--       w-component
           CASE ( 3 )
-             DO  k = nzb_w_inner(j,i)+1, pch_index-1
+             DO  k = nzb_w_inner(j,i)+1, nzb_w_inner(j,i)+pch_index-1
+
+                kk = k - nzb_w_inner(j,i)  !- lad arrays are defined flat
 
                 pre_tend = 0.0_wp
@@ -836,5 +855,5 @@
                 pre_tend = - cdc *                                             &
                              (0.5_wp *                                         &
-                                ( lad_s(k+1,j,i) + lad_s(k,j,i) )) *           &
+                                ( lad_s(kk+1,j,i) + lad_s(kk,j,i) )) *         &
                              SQRT( ( 0.25_wp * ( u(k,j,i)    +                 &  
                                                  u(k,j,i+1)  +                 &
@@ -869,8 +888,9 @@
 !--       potential temperature
           CASE ( 4 )
-             DO  k = nzb_s_inner(j,i)+1, pch_index
+             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
+                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
                 tend(k,j,i) = tend(k,j,i) +                                    &
-                                 ( canopy_heat_flux(k,j,i) -                   &
-                                   canopy_heat_flux(k-1,j,i) ) / dzw(k)
+                                 ( canopy_heat_flux(kk,j,i) -                  &
+                                   canopy_heat_flux(kk-1,j,i) ) / dzw(k)
              ENDDO
 
@@ -879,8 +899,9 @@
 !--       scalar concentration 
           CASE ( 5 )
-             DO  k = nzb_s_inner(j,i)+1, pch_index
+             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
+                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
                 tend(k,j,i) = tend(k,j,i) -                                    &
                                  lsec *                                        &
-                                 lad_s(k,j,i) *                                &
+                                 lad_s(kk,j,i) *                               &
                                  SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
                                                     u(k,j,i+1) )               &
@@ -899,8 +920,9 @@
 !--       sgs-tke
           CASE ( 6 )
-             DO  k = nzb_s_inner(j,i)+1, pch_index   
+             DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
+                kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
                 tend(k,j,i) = tend(k,j,i) -                                    &
                                  2.0_wp * cdc *                                &
-                                 lad_s(k,j,i) *                                &
+                                 lad_s(kk,j,i) *                               &
                                  SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
                                                     u(k,j,i+1) )               &
