Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2214)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2232)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Adjustments to new topography concept
 ! 
 ! Former revisions:
@@ -118,9 +118,9 @@
  
     USE arrays_3d,                                                             &
-        ONLY:  dzu, dzw, e, q, s, shf, tend, u, v, w, zu, zw 
+        ONLY:  dzu, dzw, e, q, s, tend, u, v, w, zu, zw 
 
     USE indices,                                                               &
         ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
-               nz, nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
+               nz, nzb, nzb_max, nzt, wall_flags_0
 
     USE kinds
@@ -549,8 +549,8 @@
        USE control_parameters,                                                 &
            ONLY:  coupling_char, dz, humidity, io_blocks, io_group,            &
-                  message_string, ocean, passive_scalar 
-
-       USE control_parameters,                                                 &
-           ONLY:  urban_surface
+                  message_string, ocean, passive_scalar, urban_surface
+
+       USE surface_mod,                                                        &
+           ONLY: surf_def_h, surf_lsm_h, surf_usm_h
 
        IMPLICIT NONE
@@ -562,4 +562,5 @@
        INTEGER(iwp) ::  j   !< running index
        INTEGER(iwp) ::  k   !< running index
+       INTEGER(iwp) ::  m   !< running index
 
        REAL(wp) ::  int_bpdf        !< vertical integral for lad-profile construction
@@ -639,5 +640,5 @@
           canopy_height = pch_index * dz
 
-          DO k = nzb, pch_index
+          DO k = 0, pch_index
              int_bpdf = int_bpdf +                                             &
                       ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *  &
@@ -649,5 +650,5 @@
 !
 !--       Preliminary lad profile (defined on w-grid)
-          DO k = nzb, pch_index
+          DO k = 0, pch_index
              pre_lad(k) =  lai_beta *                                          &
                         ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) )  &
@@ -662,5 +663,5 @@
 !--       when calculating the canopy tendencies)
           lad(0) = pre_lad(0)
-          DO k = nzb+1, pch_index
+          DO k = 1, pch_index
              lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
           ENDDO          
@@ -761,27 +762,45 @@
           ENDDO
 
+!            
+!--       In areas with canopy the surface value of the canopy heat 
+!--       flux distribution overrides the surface heat flux (shf) 
+!--       Start with default surface type
+          DO  m = 1, surf_def_h(0)%ns
+             k = surf_def_h(0)%k(m)
+             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
+                surf_def_h(0)%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
+          ENDDO
+!
+!--       Natural surfaces
+          DO  m = 1, surf_lsm_h%ns
+             k = surf_lsm_h%k(m)
+             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
+                surf_lsm_h%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
+          ENDDO
+!
+!--       Urban surfaces
+          DO  m = 1, surf_usm_h%ns
+             k = surf_usm_h%k(m)
+             IF ( cum_lai_hf(0,j,i) /= 0.0_wp )                                &
+                surf_usm_h%shf(m) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
+          ENDDO
+!
 !
 !--       Calculation of the heating rate (K/s) within the different layers of 
-!--       the plant canopy
+!--       the plant canopy. Calculation is only necessary in areas covered with
+!--       canopy. 
+!--       Within the different canopy layers the plant-canopy heating
+!--       rate (pc_heating_rate) is calculated as the vertical 
+!--       divergence of the canopy heat fluxes at the top and bottom
+!--       of the respective layer
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-!
-!--             Calculation only necessary in areas covered with canopy
-                IF ( cum_lai_hf(0,j,i) /= 0.0_wp )  THEN
-!--             
-!--                In areas with canopy the surface value of the canopy heat 
-!--                flux distribution overrides the surface heat flux (shf) 
-                   shf(j,i) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
-!
-!--                Within the different canopy layers the plant-canopy heating
-!--                rate (pc_heating_rate) is calculated as the vertical 
-!--                divergence of the canopy heat fluxes at the top and bottom
-!--                of the respective layer
-                   DO  k = 1, pch_index
-                      pc_heating_rate(k,j,i) = cthf *                         &
-                                                ( exp(-ext_coef*cum_lai_hf(k,j,i)) -  &
-                                                  exp(-ext_coef*cum_lai_hf(k-1,j,i)) ) / dzw(k)
-                   ENDDO
-                ENDIF
+                DO  k = 1, pch_index
+                   IF ( cum_lai_hf(0,j,i) /= 0.0_wp )  THEN
+                      pc_heating_rate(k,j,i) = cthf *                             &
+                                ( exp(-ext_coef*cum_lai_hf(k,j,i)) -              &
+                                  exp(-ext_coef*cum_lai_hf(k-1,j,i) ) ) / dzw(k)
+                   ENDIF
+                ENDDO
              ENDDO
           ENDDO
@@ -966,4 +985,5 @@
        INTEGER(iwp) ::  j         !< running index
        INTEGER(iwp) ::  k         !< running index
+       INTEGER(iwp) ::  k_wall    !< vertical index of topography top
        INTEGER(iwp) ::  kk        !< running index for flat lad arrays
 
@@ -987,7 +1007,14 @@
              DO  i = nxlu, nxr
                 DO  j = nys, nyn
-                   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
+!
+!--                Determine topography-top index on u-grid
+                   k_wall = MAXLOC(                                            &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 14 )    &
+                               ), DIM = 1                                      &
+                                  ) - 1
+                   DO  k = k_wall+1, k_wall+pch_index
+
+                      kk = k - k_wall   !- lad arrays are defined flat
 !
 !--                   In order to create sharp boundaries of the plant canopy, 
@@ -1048,7 +1075,14 @@
              DO  i = nxl, nxr
                 DO  j = nysv, nyn
-                   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
+!
+!--                Determine topography-top index on v-grid
+                   k_wall = MAXLOC(                                            &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 16 )    &
+                               ), DIM = 1                                      &
+                                  ) - 1
+                   DO  k = k_wall+1, k_wall+pch_index
+
+                      kk = k - k_wall   !- lad arrays are defined flat
 !
 !--                   In order to create sharp boundaries of the plant canopy, 
@@ -1109,7 +1143,14 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   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
+!
+!--                Determine topography-top index on w-grid
+                   k_wall = MAXLOC(                                            &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 18 )    &
+                               ), DIM = 1                                      &
+                                  ) - 1
+                   DO  k = k_wall+1, k_wall+pch_index-1
+
+                      kk = k - k_wall   !- lad arrays are defined flat
 
                       pre_tend = 0.0_wp
@@ -1157,6 +1198,14 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   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
+!
+!--                Determine topography-top index on scalar-grid
+                   k_wall = MAXLOC(                                            &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                                  ) - 1
+                   DO  k = k_wall+1, k_wall+pch_index
+
+                      kk = k - k_wall   !- lad arrays are defined flat
                       tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
                    ENDDO
@@ -1169,6 +1218,14 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   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
+!
+!--                Determine topography-top index on scalar-grid
+                   k_wall = MAXLOC(                                            &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                                  ) - 1
+                   DO  k = k_wall+1, k_wall+pch_index
+
+                      kk = k - k_wall   !- lad arrays are defined flat
                       tend(k,j,i) = tend(k,j,i) -                              &
                                        lsec *                                  &
@@ -1194,6 +1251,14 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   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
+!
+!--                Determine topography-top index on scalar-grid
+                   k_wall = MAXLOC(                                            &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                                  ) - 1
+                   DO  k = k_wall+1, k_wall+pch_index
+
+                      kk = k - k_wall   !- lad arrays are defined flat
                       tend(k,j,i) = tend(k,j,i) -                              &
                                        2.0_wp * cdc *                          &
@@ -1218,6 +1283,14 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   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
+!
+!--                Determine topography-top index on scalar-grid
+                   k_wall = MAXLOC(                                            &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                                  ) - 1
+                   DO  k = k_wall+1, k_wall+pch_index
+
+                      kk = k - k_wall   !- lad arrays are defined flat
                       tend(k,j,i) = tend(k,j,i) -                              &
                                        lsec *                                  &
@@ -1288,4 +1361,5 @@
        INTEGER(iwp) ::  j         !< running index
        INTEGER(iwp) ::  k         !< running index
+       INTEGER(iwp) ::  k_wall    !< vertical index of topography top
        INTEGER(iwp) ::  kk        !< running index for flat lad arrays
 
@@ -1307,7 +1381,14 @@
 !--       u-component
           CASE ( 1 )
-             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
+!
+!--          Determine topography-top index on u-grid
+             k_wall = MAXLOC(                                                  &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 14 )    &
+                               ), DIM = 1                                      &
+                            ) - 1
+             DO  k = k_wall+1, k_wall+pch_index
+
+                kk = k - k_wall  !- lad arrays are defined flat
 !
 !--             In order to create sharp boundaries of the plant canopy, 
@@ -1363,7 +1444,14 @@
 !--       v-component
           CASE ( 2 )
-             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
+!
+!--          Determine topography-top index on v-grid
+             k_wall = MAXLOC(                                                  &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 16 )    &
+                               ), DIM = 1                                      &
+                            ) - 1
+             DO  k = k_wall+1, k_wall+pch_index
+
+                kk = k - k_wall  !- lad arrays are defined flat
 !
 !--             In order to create sharp boundaries of the plant canopy, 
@@ -1419,7 +1507,14 @@
 !--       w-component
           CASE ( 3 )
-             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
+!
+!--          Determine topography-top index on w-grid
+             k_wall = MAXLOC(                                                  &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 18 )    &
+                               ), DIM = 1                                      &
+                            ) - 1
+             DO  k = k_wall+1, k_wall+pch_index-1
+
+                kk = k - k_wall  !- lad arrays are defined flat
 
                 pre_tend = 0.0_wp
@@ -1462,6 +1557,13 @@
 !--       potential temperature
           CASE ( 4 )
-             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
+!
+!--          Determine topography-top index on scalar grid
+             k_wall = MAXLOC(                                                  &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                            ) - 1
+             DO  k = k_wall+1, k_wall+pch_index
+                kk = k - k_wall  !- lad arrays are defined flat
                 tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
              ENDDO
@@ -1471,6 +1573,14 @@
 !--       humidity
           CASE ( 5 )
-             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
+!
+!--          Determine topography-top index on scalar grid
+             k_wall = MAXLOC(                                                  &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                            ) - 1
+             DO  k = k_wall+1, k_wall+pch_index
+
+                kk = k - k_wall
                 tend(k,j,i) = tend(k,j,i) -                                    &
                                  lsec *                                        &
@@ -1492,6 +1602,14 @@
 !--       sgs-tke
           CASE ( 6 )
-             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
+!
+!--          Determine topography-top index on scalar grid
+             k_wall = MAXLOC(                                                  &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                            ) - 1
+             DO  k = k_wall+1, k_wall+pch_index
+
+                kk = k - k_wall
                 tend(k,j,i) = tend(k,j,i) -                                    &
                                  2.0_wp * cdc *                                &
@@ -1513,6 +1631,14 @@
 !--       scalar concentration 
           CASE ( 7 )
-             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
+!
+!--          Determine topography-top index on scalar grid
+             k_wall = MAXLOC(                                                  &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                            ) - 1
+             DO  k = k_wall+1, k_wall+pch_index
+
+                kk = k - k_wall
                 tend(k,j,i) = tend(k,j,i) -                                    &
                                  lsec *                                        &
