Index: palm/trunk/SOURCE/microphysics_mod.f90
===================================================================
--- palm/trunk/SOURCE/microphysics_mod.f90	(revision 2156)
+++ palm/trunk/SOURCE/microphysics_mod.f90	(revision 2232)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Adjustments to new topography and surface concept
 ! 
 ! Former revisions:
@@ -402,5 +402,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
 
        USE kinds
@@ -416,5 +416,5 @@
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
-             DO  k = nzb_s_inner(j,i)+1, nzt
+             DO  k = nzb+1, nzt
                 IF ( qr(k,j,i) <= eps_sb )  THEN
                    qr(k,j,i) = 0.0_wp
@@ -422,7 +422,11 @@
                 ELSE
                    IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
-                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
+                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin *               &
+                                       MERGE( 1.0_wp, 0.0_wp,                  &           
+                                              BTEST( wall_flags_0(k,j,i), 0 ) )
                    ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
-                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
+                      nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax *               &
+                                       MERGE( 1.0_wp, 0.0_wp,                  &           
+                                              BTEST( wall_flags_0(k,j,i), 0 ) )
                    ENDIF
                 ENDIF
@@ -459,5 +463,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
 
        USE kinds
@@ -472,4 +476,5 @@
        REAL(wp)     ::  autocon           !<
        REAL(wp)     ::  dissipation       !<
+       REAL(wp)     ::  flag              !< flag to mask topography grid points
        REAL(wp)     ::  k_au              !<
        REAL(wp)     ::  l_mix             !<
@@ -487,5 +492,8 @@
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
-             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 ) )
 
                 IF ( qc(k,j,i) > eps_sb )  THEN
@@ -554,7 +562,8 @@
                    autocon = MIN( autocon, qc(k,j,i) / dt_micro )
 
-                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro
-                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro
-                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro
+                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag
+                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag
+                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro  &
+                                                              * flag
 
                 ENDIF
@@ -583,5 +592,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nyng, nysg, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
 
        USE kinds
@@ -590,13 +599,25 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i !<
-       INTEGER(iwp) ::  j !<
-       INTEGER(iwp) ::  k !<
+       INTEGER(iwp) ::  i      !<
+       INTEGER(iwp) ::  j      !<
+       INTEGER(iwp) ::  k      !<
+       INTEGER(iwp) ::  k_wall !< topgraphy top index
 
        REAL(wp)    ::  dqdt_precip !<
+       REAL(wp)    ::  flag        !< flag to mask topography grid points
 
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
-             DO  k = nzb_s_inner(j,i)+1, nzt
+!
+!--          Determine vertical index of topography top
+             k_wall = MAXLOC(                                                  &
+                          MERGE( 1, 0,                                         &
+                                 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )    &
+                               ), DIM = 1                                      &
+                            ) - 1
+             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 ) )
 
                 IF ( qc(k,j,i) > ql_crit )  THEN
@@ -606,13 +627,12 @@
                 ENDIF
 
-                qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro
-                q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro
+                qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
+                q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
                 pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp *      &
-                                        pt_d_t(k)
+                                        pt_d_t(k)              * flag
 
 !
 !--             Compute the rain rate (stored on surface grid point)
-                prr(nzb_s_inner(j,i),j,i) = prr(nzb_s_inner(j,i),j,i) +        &
-                                            dqdt_precip * dzw(k)
+                prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
 
              ENDDO
@@ -643,5 +663,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
 
        USE kinds
@@ -654,4 +674,5 @@
 
        REAL(wp)     ::  accr              !<
+       REAL(wp)     ::  flag              !< flag to mask topography grid points
        REAL(wp)     ::  k_cr              !<
        REAL(wp)     ::  phi_ac            !<
@@ -662,5 +683,8 @@
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
-             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 ) )
 
                 IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb ) )  THEN
@@ -690,6 +714,6 @@
                    accr = MIN( accr, qc(k,j,i) / dt_micro )
 
-                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro
-                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro
+                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
+                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
 
                 ENDIF
@@ -724,5 +748,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
 
        USE kinds
@@ -736,4 +760,5 @@
        REAL(wp)     ::  breakup           !<
        REAL(wp)     ::  dr                !<
+       REAL(wp)     ::  flag              !< flag to mask topography grid points
        REAL(wp)     ::  phi_br            !<
        REAL(wp)     ::  selfcoll          !<
@@ -743,5 +768,9 @@
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
-             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 ) )
+
                 IF ( qr(k,j,i) > eps_sb )  THEN
 !
@@ -763,5 +792,5 @@
 
                    selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
-                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro
+                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
 
                 ENDIF
@@ -796,5 +825,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
 
        USE kinds
@@ -812,4 +841,5 @@
        REAL(wp)     ::  evap_nr           !<
        REAL(wp)     ::  f_vent            !<
+       REAL(wp)     ::  flag              !< flag to mask topography grid points
        REAL(wp)     ::  g_evap            !<
        REAL(wp)     ::  lambda_r          !<
@@ -828,5 +858,9 @@
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
-             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 ) )
+
                 IF ( qr(k,j,i) > eps_sb )  THEN
 !
@@ -916,6 +950,6 @@
                                      -nr(k,j,i) / dt_micro )
 
-                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro
-                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro
+                      qr(k,j,i) = qr(k,j,i) + evap    * dt_micro * flag
+                      nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro * flag
 
                    ENDIF
@@ -951,5 +985,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
 
        USE kinds
@@ -965,5 +999,6 @@
        INTEGER(iwp) ::  k             !<
 
-       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<
+       REAL(wp)                       ::  flag   !< flag to mask topography grid points
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
 
        CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
@@ -973,9 +1008,13 @@
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
-             DO  k = nzt, nzb_s_inner(j,i)+1, -1
+             DO  k = nzt, nzb+1, -1
+!
+!--             Predetermine flag to mask topography
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
 
                 IF ( qc(k,j,i) > eps_sb )  THEN
                    sed_qc(k) = sed_qc_const * nc_const**( -2.0_wp / 3.0_wp ) * &
-                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp )
+                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
+                                                                           flag
                 ELSE
                    sed_qc(k) = 0.0_wp
@@ -984,13 +1023,13 @@
                 sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
                                             dt_micro + sed_qc(k+1)             &
-                               )
+                               ) * flag
 
                 q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * dt_micro
+                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
                 qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * dt_micro
+                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
                 pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
                                         ddzu(k+1) / hyrho(k) * l_d_cp *        &
-                                        pt_d_t(k) * dt_micro
+                                        pt_d_t(k) * dt_micro            * flag
 
 !
@@ -998,7 +1037,8 @@
                 IF ( call_microphysics_at_all_substeps )  THEN
                    prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
-                                * weight_substep(intermediate_timestep_count)
+                                * weight_substep(intermediate_timestep_count)  &
+                                * flag
                 ELSE
-                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
+                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
                 ENDIF
 
@@ -1032,5 +1072,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzb_s_inner, nzt
+           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
 
        USE kinds
@@ -1039,10 +1079,17 @@
            ONLY:  weight_substep
 
+       USE surface_mod,                                                        &
+           ONLY :  bc_h
+       
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i                          !<
-       INTEGER(iwp) ::  j                          !<
-       INTEGER(iwp) ::  k                          !<
-       INTEGER(iwp) ::  k_run                      !<
+       INTEGER(iwp) ::  i             !< running index x direction
+       INTEGER(iwp) ::  j             !< running index y direction
+       INTEGER(iwp) ::  k             !< running index z direction
+       INTEGER(iwp) ::  k_run         !< 
+       INTEGER(iwp) ::  l             !< running index of surface type
+       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)     ::  c_run                      !<
@@ -1052,4 +1099,5 @@
        REAL(wp)     ::  dr                         !<
        REAL(wp)     ::  flux                       !<
+       REAL(wp)     ::  flag                       !< flag to mask topography grid points
        REAL(wp)     ::  lambda_r                   !<
        REAL(wp)     ::  mu_r                       !<
@@ -1071,5 +1119,9 @@
        DO  i = nxlg, nxrg
           DO  j = nysg, nyng
-             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 ) )
+
                 IF ( qr(k,j,i) > eps_sb )  THEN
 !
@@ -1093,5 +1145,5 @@
                                                   ( mu_r + 1.0_wp ) )          &
                                               )                                &
-                                )
+                                ) * flag
 
                    w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
@@ -1101,5 +1153,5 @@
                                                   ( mu_r + 4.0_wp ) )          &
                                              )                                 &
-                                )
+                                ) * flag
                 ELSE
                    w_nr(k) = 0.0_wp
@@ -1108,18 +1160,39 @@
              ENDDO
 !
-!--          Adjust boundary values
-             w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
-             w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
+!--          Adjust boundary values using surface data type. 
+!--          Upward-facing
+             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)
+                w_nr(k-1) = w_nr(k)
+                w_qr(k-1) = w_qr(k)
+             ENDDO
+!
+!--          Downward-facing
+             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)
+                w_nr(k+1) = w_nr(k)
+                w_qr(k+1) = w_qr(k)
+             ENDDO
+!
+!--          Model top boundary value
              w_nr(nzt+1) = 0.0_wp
              w_qr(nzt+1) = 0.0_wp
 !
 !--          Compute Courant number
-             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 ) )
+
                 c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
                                       2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
-                          dt_micro * ddzu(k)
+                          dt_micro * ddzu(k) * flag
                 c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
                                       2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
-                          dt_micro * ddzu(k)
+                          dt_micro * ddzu(k) * flag
              ENDDO
 !
@@ -1127,5 +1200,9 @@
              IF ( limiter_sedimentation )  THEN
 
-                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 ) )
+
                    d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
                    d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
@@ -1134,5 +1211,6 @@
                    qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
                                                               2.0_wp * d_max,  &
-                                                              ABS( d_mean ) )
+                                                              ABS( d_mean ) )  &
+                                                      * flag
 
                    d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
@@ -1156,7 +1234,10 @@
 !
 !--          Compute sedimentation flux
-             DO  k = nzt, nzb_s_inner(j,i)+1, -1
-!
-!--             Sum up all rain drop number densities which contribute to the flux
+             DO  k = nzt, nzb+1, -1
+!
+!--             Predetermine flag to mask topography
+                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+!
+!--             Sum up all rain drop number densities which contribute to the flux 
 !--             through k-1/2
                 flux  = 0.0_wp
@@ -1167,8 +1248,10 @@
                    flux  = flux + hyrho(k_run) *                               &
                            ( nr(k_run,j,i) + nr_slope(k_run) *                 &
-                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)
-                   z_run = z_run + dzu(k_run)
-                   k_run = k_run + 1
-                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
+                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)  &
+                                              * flag
+                   z_run = z_run + dzu(k_run) * flag
+                   k_run = k_run + 1          * flag
+                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )    &
+                                              * flag
                 ENDDO
 !
@@ -1180,7 +1263,7 @@
                           )
 
-                sed_nr(k) = flux / dt_micro
+                sed_nr(k) = flux / dt_micro * flag
                 nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * dt_micro
+                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
 !
 !--             Sum up all rain water content which contributes to the flux
@@ -1195,8 +1278,9 @@
                    flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
                                   qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
-                                  0.5_wp ) * c_run * dzu(k_run)
-                   z_run = z_run + dzu(k_run)
-                   k_run = k_run + 1
-                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
+                                  0.5_wp ) * c_run * dzu(k_run) * flag
+                   z_run = z_run + dzu(k_run)                   * flag
+                   k_run = k_run + 1                            * flag
+                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )    &
+                                                                * flag
 
                 ENDDO
@@ -1209,20 +1293,21 @@
                           )
 
-                sed_qr(k) = flux / dt_micro
+                sed_qr(k) = flux / dt_micro * flag
 
                 qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * dt_micro
+                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
                 q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
-                                        ddzu(k+1) / hyrho(k) * dt_micro
+                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
                 pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
                                         ddzu(k+1) / hyrho(k) * l_d_cp *        &
-                                        pt_d_t(k) * dt_micro
+                                        pt_d_t(k) * dt_micro            * flag
 !
 !--             Compute the rain rate
                 IF ( call_microphysics_at_all_substeps )  THEN
                    prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
-                                * weight_substep(intermediate_timestep_count)
+                                * weight_substep(intermediate_timestep_count) &
+                                * flag
                 ELSE
-                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
+                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
                 ENDIF
 
@@ -1256,13 +1341,19 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxr, nys, nyn, nzb_s_inner
+           ONLY:  nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0
 
        USE kinds
 
+       USE surface_mod,                                                        &
+           ONLY :  bc_h
+
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i                          !:
-       INTEGER(iwp) ::  j                          !:
-
+       INTEGER(iwp) ::  i             !< running index x direction
+       INTEGER(iwp) ::  j             !< running index y direction
+       INTEGER(iwp) ::  k             !< running index y 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
 
        IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
@@ -1270,14 +1361,15 @@
             intermediate_timestep_count == intermediate_timestep_count_max ) ) &
        THEN
-
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-
-                precipitation_amount(j,i) = precipitation_amount(j,i) +        &
-                                            prr(nzb_s_inner(j,i)+1,j,i) *      &
-                                            hyrho(nzb_s_inner(j,i)+1) * dt_3d
-
-             ENDDO
+!
+!--       Run over all upward-facing surface elements, i.e. non-natural, 
+!--       natural and urban
+          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)
+             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
+                                               prr(k,j,i) * hyrho(k) * dt_3d
           ENDDO
+
        ENDIF
 
@@ -1410,5 +1502,5 @@
 
        USE indices,                                                            &
-           ONLY:  nzb_s_inner, nzt
+           ONLY:  nzb, nzt, wall_flags_0
 
        USE kinds
@@ -1420,5 +1512,10 @@
        INTEGER(iwp) ::  k                 !<
 
-       DO  k = nzb_s_inner(j,i)+1, nzt
+       REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
+
+       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 ) )
 
           IF ( qr_1d(k) <= eps_sb )  THEN
@@ -1431,7 +1528,7 @@
 !--          too big weights of rain drops (Stevens and Seifert, 2008).
              IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) )  THEN
-                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin
+                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin * flag
              ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) )  THEN
-                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax
+                nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax * flag
              ENDIF
 
@@ -1463,5 +1560,5 @@
 
        USE indices,                                                            &
-           ONLY:  nzb_s_inner, nzt
+           ONLY:  nzb, nzt, wall_flags_0
 
        USE kinds
@@ -1476,4 +1573,5 @@
        REAL(wp)     ::  autocon           !<
        REAL(wp)     ::  dissipation       !<
+       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
        REAL(wp)     ::  k_au              !<
        REAL(wp)     ::  l_mix             !<
@@ -1487,5 +1585,8 @@
        REAL(wp)     ::  xc                !<
 
-       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 ) )
 
           IF ( qc_1d(k) > eps_sb )  THEN
@@ -1551,7 +1652,7 @@
              autocon = MIN( autocon, qc_1d(k) / dt_micro )
 
-             qr_1d(k) = qr_1d(k) + autocon * dt_micro
-             qc_1d(k) = qc_1d(k) - autocon * dt_micro
-             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
+             qr_1d(k) = qr_1d(k) + autocon * dt_micro                 * flag
+             qc_1d(k) = qc_1d(k) - autocon * dt_micro                 * flag
+             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag
 
           ENDIF
@@ -1575,5 +1676,5 @@
 
        USE indices,                                                            &
-           ONLY:  nzb_s_inner, nzt
+           ONLY:  nzb, nzb_max, nzt, wall_flags_0
 
        USE kinds
@@ -1582,11 +1683,23 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i !<
-       INTEGER(iwp) ::  j !<
-       INTEGER(iwp) ::  k !<
+       INTEGER(iwp) ::  i      !<
+       INTEGER(iwp) ::  j      !<
+       INTEGER(iwp) ::  k      !<
+       INTEGER(iwp) ::  k_wall !< topography top index
 
        REAL(wp)    ::  dqdt_precip !<
-
-       DO  k = nzb_s_inner(j,i)+1, nzt
+       REAL(wp)    ::  flag              !< flag to indicate first grid level above surface
+
+!
+!--    Determine vertical index of topography top
+       k_wall = MAXLOC(                                                        &
+                        MERGE( 1, 0,                                           &
+                               BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 )      &
+                             ), DIM = 1                                        &
+                      ) - 1
+       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 ) )
 
           IF ( qc_1d(k) > ql_crit )  THEN
@@ -1596,12 +1709,11 @@
           ENDIF
 
-          qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro
-          q_1d(k)  = q_1d(k)  - dqdt_precip * dt_micro
-          pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k)
+          qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro * flag
+          q_1d(k)  = q_1d(k)  - dqdt_precip * dt_micro * flag
+          pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k) * flag
 
 !
 !--       Compute the rain rate (stored on surface grid point)
-          prr(nzb_s_inner(j,i),j,i) = prr(nzb_s_inner(j,i),j,i) +              &
-                                      dqdt_precip * dzw(k)
+          prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
 
        ENDDO
@@ -1626,5 +1738,5 @@
 
        USE indices,                                                            &
-           ONLY:  nzb_s_inner, nzt
+           ONLY:  nzb, nzt, wall_flags_0
 
        USE kinds
@@ -1637,9 +1749,14 @@
 
        REAL(wp)     ::  accr              !<
+       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
        REAL(wp)     ::  k_cr              !<
        REAL(wp)     ::  phi_ac            !<
        REAL(wp)     ::  tau_cloud         !<
 
-       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 ) )
+
           IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
 !
@@ -1667,6 +1784,6 @@
              accr = MIN( accr, qc_1d(k) / dt_micro )
 
-             qr_1d(k) = qr_1d(k) + accr * dt_micro
-             qc_1d(k) = qc_1d(k) - accr * dt_micro
+             qr_1d(k) = qr_1d(k) + accr * dt_micro * flag
+             qc_1d(k) = qc_1d(k) - accr * dt_micro * flag
 
           ENDIF
@@ -1691,5 +1808,5 @@
 
        USE indices,                                                            &
-           ONLY:  nzb_s_inner, nzt
+           ONLY:  nzb, nzt, wall_flags_0
 
        USE kinds
@@ -1703,8 +1820,13 @@
        REAL(wp)     ::  breakup           !<
        REAL(wp)     ::  dr                !<
+       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
        REAL(wp)     ::  phi_br            !<
        REAL(wp)     ::  selfcoll          !<
 
-       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 ) )
+
           IF ( qr_1d(k) > eps_sb )  THEN
 !
@@ -1724,5 +1846,5 @@
 
              selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro )
-             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
+             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro * flag
 
           ENDIF
@@ -1750,5 +1872,5 @@
 
        USE indices,                                                            &
-           ONLY:  nzb_s_inner, nzt
+           ONLY:  nzb, nzt, wall_flags_0
 
        USE kinds
@@ -1766,4 +1888,5 @@
        REAL(wp)     ::  evap_nr           !<
        REAL(wp)     ::  f_vent            !<
+       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
        REAL(wp)     ::  g_evap            !<
        REAL(wp)     ::  lambda_r          !<
@@ -1778,5 +1901,9 @@
        REAL(wp)     ::  xr                !<
 
-       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 ) )
+
           IF ( qr_1d(k) > eps_sb )  THEN
 !
@@ -1862,6 +1989,6 @@
                                -nr_1d(k) / dt_micro )
 
-                qr_1d(k) = qr_1d(k) + evap    * dt_micro
-                nr_1d(k) = nr_1d(k) + evap_nr * dt_micro
+                qr_1d(k) = qr_1d(k) + evap    * dt_micro * flag
+                nr_1d(k) = nr_1d(k) + evap_nr * dt_micro * flag
 
              ENDIF
@@ -1891,5 +2018,5 @@
 
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nzb, nzb, nzt, wall_flags_0
 
        USE kinds
@@ -1904,12 +2031,17 @@
        INTEGER(iwp) ::  k             !<
 
-       REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc  !<
+       REAL(wp)                       ::  flag    !< flag to indicate first grid level above surface
+       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc  !<
 
        sed_qc(nzt+1) = 0.0_wp
 
-       DO  k = nzt, nzb_s_inner(j,i)+1, -1
+       DO  k = nzt, nzb+1, -1
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+
           IF ( qc_1d(k) > eps_sb )  THEN
              sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) *       &
-                         ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
+                         ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )  * flag
           ELSE
              sed_qc(k) = 0.0_wp
@@ -1918,20 +2050,20 @@
           sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) /          &
                                       dt_micro + sed_qc(k+1)                   &
-                         )
+                         ) * flag
 
           q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
-                                hyrho(k) * dt_micro
+                                hyrho(k) * dt_micro * flag
           qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
-                                hyrho(k) * dt_micro
+                                hyrho(k) * dt_micro * flag
           pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
-                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
+                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
 
 !
 !--       Compute the precipitation rate of cloud (fog) droplets
           IF ( call_microphysics_at_all_substeps )  THEN
-             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *               &
-                             weight_substep(intermediate_timestep_count)
+             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) *                  &
+                              weight_substep(intermediate_timestep_count) * flag
           ELSE
-             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k)
+             prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
           ENDIF
 
@@ -1959,5 +2091,5 @@
 
        USE indices,                                                            &
-           ONLY:  nzb, nzb_s_inner, nzt
+           ONLY:  nzb, nzb, nzt, wall_flags_0
 
        USE kinds
@@ -1966,10 +2098,16 @@
            ONLY:  weight_substep
 
+       USE surface_mod,                                                        &
+           ONLY :  bc_h
+       
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i                          !<
-       INTEGER(iwp) ::  j                          !<
-       INTEGER(iwp) ::  k                          !<
-       INTEGER(iwp) ::  k_run                      !<
+       INTEGER(iwp) ::  i             !< running index x direction
+       INTEGER(iwp) ::  j             !< running index y direction
+       INTEGER(iwp) ::  k             !< running index z direction
+       INTEGER(iwp) ::  k_run         !<
+       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)     ::  c_run                      !<
@@ -1979,4 +2117,5 @@
        REAL(wp)     ::  dr                         !<
        REAL(wp)     ::  flux                       !<
+       REAL(wp)     ::  flag                       !< flag to indicate first grid level above surface
        REAL(wp)     ::  lambda_r                   !<
        REAL(wp)     ::  mu_r                       !<
@@ -1994,5 +2133,9 @@
 !
 !--    Compute velocities
-       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 ) )
+
           IF ( qr_1d(k) > eps_sb )  THEN
 !
@@ -2013,5 +2156,5 @@
                                             ( mu_r + 1.0_wp ) )                &
                                         )                                      &
-                          )
+                          ) * flag
              w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                              &
                                          a_term - b_term * ( 1.0_wp +          &
@@ -2019,5 +2162,5 @@
                                             ( mu_r + 4.0_wp ) )                &
                                        )                                       &
-                          )
+                          ) * flag
           ELSE
              w_nr(k) = 0.0_wp
@@ -2026,16 +2169,37 @@
        ENDDO
 !
-!--    Adjust boundary values
-       w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1)
-       w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1)
+!--    Adjust boundary values using surface data type. 
+!--    Upward facing non-natural
+       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)
+          w_nr(k-1) = w_nr(k)
+          w_qr(k-1) = w_qr(k)
+       ENDDO
+!
+!--    Downward facing non-natural
+       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)
+          w_nr(k+1) = w_nr(k)
+          w_qr(k+1) = w_qr(k)
+       ENDDO
+!
+!--    Neumann boundary condition at model top
        w_nr(nzt+1) = 0.0_wp
        w_qr(nzt+1) = 0.0_wp
 !
 !--    Compute Courant number
-       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 ) )
+
           c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) *   &
-                    dt_micro * ddzu(k)
+                    dt_micro * ddzu(k) * flag
           c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
-                    dt_micro * ddzu(k)
+                    dt_micro * ddzu(k) * flag
        ENDDO
 !
@@ -2043,5 +2207,9 @@
        IF ( limiter_sedimentation )  THEN
 
-          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 ) )
+
              d_mean = 0.5_wp * ( qr_1d(k+1) - qr_1d(k-1) )
              d_min  = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) )
@@ -2050,5 +2218,5 @@
              qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
                                                         2.0_wp * d_max,        &
-                                                        ABS( d_mean ) )
+                                                        ABS( d_mean ) ) * flag
 
              d_mean = 0.5_wp * ( nr_1d(k+1) - nr_1d(k-1) )
@@ -2058,5 +2226,5 @@
              nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,        &
                                                         2.0_wp * d_max,        &
-                                                        ABS( d_mean ) )
+                                                        ABS( d_mean ) ) * flag
           ENDDO
 
@@ -2072,7 +2240,10 @@
 !
 !--    Compute sedimentation flux
-       DO  k = nzt, nzb_s_inner(j,i)+1, -1
-!
-!--       Sum up all rain drop number densities which contribute to the flux
+       DO  k = nzt, nzb+1, -1
+!
+!--       Predetermine flag to mask topography
+          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
+!
+!--       Sum up all rain drop number densities which contribute to the flux 
 !--       through k-1/2
           flux  = 0.0_wp
@@ -2083,8 +2254,8 @@
              flux  = flux + hyrho(k_run) *                                     &
                      ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) *   &
-                     0.5_wp ) * c_run * dzu(k_run)
-             z_run = z_run + dzu(k_run)
-             k_run = k_run + 1
-             c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )
+                     0.5_wp ) * c_run * dzu(k_run) * flag
+             z_run = z_run + dzu(k_run)            * flag
+             k_run = k_run + 1                     * flag
+             c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) * flag
           ENDDO
 !
@@ -2094,7 +2265,7 @@
                       hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro )
 
-          sed_nr(k) = flux / dt_micro
+          sed_nr(k) = flux / dt_micro * flag
           nr_1d(k)  = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) /     &
-                                    hyrho(k) * dt_micro
+                                    hyrho(k) * dt_micro * flag
 !
 !--       Sum up all rain water content which contributes to the flux
@@ -2109,8 +2280,8 @@
              flux  = flux + hyrho(k_run) *                                     &
                      ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) *   &
-                     0.5_wp ) * c_run * dzu(k_run)
-             z_run = z_run + dzu(k_run)
-             k_run = k_run + 1
-             c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )
+                     0.5_wp ) * c_run * dzu(k_run) * flag
+             z_run = z_run + dzu(k_run)            * flag
+             k_run = k_run + 1                     * flag
+             c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) * flag
 
           ENDDO
@@ -2120,19 +2291,19 @@
                       hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro )
 
-          sed_qr(k) = flux / dt_micro
+          sed_qr(k) = flux / dt_micro * flag
 
           qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
-                                hyrho(k) * dt_micro
+                                hyrho(k) * dt_micro * flag
           q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
-                                hyrho(k) * dt_micro
+                                hyrho(k) * dt_micro * flag
           pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
-                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
+                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro * flag
 !
 !--       Compute the rain rate
           IF ( call_microphysics_at_all_substeps )  THEN
              prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)                    &
-                          * weight_substep(intermediate_timestep_count)
+                          * weight_substep(intermediate_timestep_count) * flag
           ELSE
-             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)
+             prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
           ENDIF
 
@@ -2162,13 +2333,19 @@
 
        USE indices,                                                            &
-           ONLY:  nzb_s_inner
+           ONLY:  nzb, nzt, wall_flags_0
 
        USE kinds
 
+       USE surface_mod,                                                        &
+           ONLY :  bc_h
+
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  i                          !:
-       INTEGER(iwp) ::  j                          !:
-
+       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
 
        IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
@@ -2177,7 +2354,12 @@
        THEN
 
-          precipitation_amount(j,i) = precipitation_amount(j,i) +              &
-                                      prr(nzb_s_inner(j,i)+1,j,i) *            &
-                                      hyrho(nzb_s_inner(j,i)+1) * dt_3d
+          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)
+             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
+                                               prr(k,j,i) * hyrho(k) * dt_3d
+          ENDDO
+
        ENDIF
 
