Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 4166)
+++ palm/trunk/SOURCE/Makefile	(revision 4167)
@@ -25,4 +25,8 @@
 # -----------------
 # $Id$
+# Remove no longer needed dependencies on surface_mod
+#
+# 
+# 4127 2019-07-30 14:47:10Z suehring
 # Add dependency of data_output_3d on plant_canopy_model_mod 
 # (merge from branch resler)
@@ -930,6 +934,5 @@
 	modules.o \
 	netcdf_interface_mod.o \
-	salsa_mod.o \
-	surface_mod.o
+	salsa_mod.o
 data_output_netcdf4_module.o: \
 	mod_kinds.o
@@ -987,6 +990,5 @@
 diagnostic_output_quantities_mod.o: \
 	mod_kinds.o \
-	modules.o \
-	surface_mod.o
+	modules.o
 diffusion_s.o: \
 	mod_kinds.o \
@@ -1642,5 +1644,4 @@
 	mod_kinds.o \
 	modules.o \
-	surface_mod.o \
 	user_module.o
 user_flight.o: \
Index: palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4166)
+++ palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 4167)
@@ -27,4 +27,9 @@
 ! -----------------
 ! $Id$
+! Changed behaviour of masked output over surface to follow terrain and ignore 
+! buildings (J.Resler, T.Gronemeier)
+! 
+! 
+! 4166 2019-08-16 07:54:21Z resler
 ! Bugfix in decycling
 ! 
@@ -1554,12 +1559,7 @@
     USE control_parameters
 
-    USE surface_mod,                                                                  &
-        ONLY:  get_topography_top_index_ji
-
-
-    CHARACTER(LEN=5)  ::  grid        !< flag to distinquish between staggered grids
     CHARACTER(LEN=16) ::  spec_name
-    CHARACTER(LEN=*)  ::  variable    !< 
-    
+    CHARACTER(LEN=*)  ::  variable    !<
+
     INTEGER(iwp) ::  av              !< flag to control data output of instantaneous or time-averaged data
     INTEGER(iwp) ::  lsp
@@ -1567,6 +1567,9 @@
     INTEGER(iwp) ::  j               !< grid index along y-direction
     INTEGER(iwp) ::  k               !< grid index along z-direction
+    INTEGER(iwp) ::  im              !< loop index for masked variables
+    INTEGER(iwp) ::  jm              !< loop index for masked variables
+    INTEGER(iwp) ::  kk              !< masked output index along z-direction
     INTEGER(iwp) ::  mid             !< masked output running index
-    INTEGER(iwp) ::  topo_top_ind    !< k index of highest horizontal surface
+    INTEGER(iwp) ::  ktt             !< k index of highest terrain surface
     
     LOGICAL ::  found
@@ -1574,11 +1577,12 @@
     REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
               local_pf   !< 
+
+    REAL(wp), PARAMETER ::  fill_value = -9999.0_wp    !< value for the _FillValue attribute
+
 !
 !-- local variables.
 
-    found = .TRUE.
-    grid = 's'
-
     spec_name = TRIM( variable(4:) )
+    found = .FALSE.
 
     DO  lsp=1,nspec
@@ -1607,18 +1611,16 @@
                 DO  i = 1, mask_size_l(mid,1)
                    DO  j = 1, mask_size_l(mid,2)
-!              
-!--                   Get k index of highest horizontal surface
-                      topo_top_ind = get_topography_top_index_ji( &
-                                        mask_j(mid,j),  &
-                                        mask_i(mid,i),  &
-                                        grid                    )
-!              
-!--                   Save output array
+!--                   Get k index of the highest terraing surface
+                      im = mask_i(mid,i)
+                      jm = mask_j(mid,j)
+                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                       DO  k = 1, mask_size_l(mid,3)
-                         local_pf(i,j,k) = chem_species(lsp)%conc( &
-                                              MIN( topo_top_ind+mask_k(mid,k), &
-                                                   nzt+1 ),        &
-                                              mask_j(mid,j),       &
-                                              mask_i(mid,i)      )
+                         kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!--                      Set value if not in building
+                         IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                            local_pf(i,j,k) = fill_value
+                         ELSE
+                            local_pf(i,j,k) = chem_species(lsp)%conc(kk,jm,im)
+                         ENDIF
                       ENDDO
                    ENDDO
@@ -1645,18 +1647,16 @@
                 DO  i = 1, mask_size_l(mid,1)
                    DO  j = 1, mask_size_l(mid,2)
-!              
-!--                   Get k index of highest horizontal surface
-                      topo_top_ind = get_topography_top_index_ji( &
-                                        mask_j(mid,j),  &
-                                        mask_i(mid,i),  &
-                                        grid                    )
-!              
-!--                   Save output array
+!--                   Get k index of the highest terraing surface
+                      im = mask_i(mid,i)
+                      jm = mask_j(mid,j)
+                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                       DO  k = 1, mask_size_l(mid,3)
-                         local_pf(i,j,k) = chem_species(lsp)%conc_av(  &
-                                              MIN( topo_top_ind+mask_k(mid,k), &
-                                                   nzt+1 ),            &
-                                              mask_j(mid,j),           &
-                                              mask_i(mid,i)         )
+                         kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!--                      Set value if not in building
+                         IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                            local_pf(i,j,k) = fill_value
+                         ELSE
+                            local_pf(i,j,k) = chem_species(lsp)%conc_av(kk,jm,im)
+                         ENDIF
                       ENDDO
                    ENDDO
@@ -1666,5 +1666,6 @@
 
           ENDIF
-          found = .FALSE.
+          found = .TRUE.
+          EXIT
        ENDIF
     ENDDO
Index: palm/trunk/SOURCE/data_output_mask.f90
===================================================================
--- palm/trunk/SOURCE/data_output_mask.f90	(revision 4166)
+++ palm/trunk/SOURCE/data_output_mask.f90	(revision 4167)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Changed behaviour of masked output over surface to follow terrain and ignore 
+! buildings (J.Resler, T.Gronemeier)
+! 
+! 4069 2019-07-01 14:05:51Z Giersch
 ! Masked output running index mid has been introduced as a local variable to 
 ! avoid runtime error (Loop variable has been modified) in time_integration 
@@ -194,5 +198,5 @@
 
     USE indices,                                                               &
-        ONLY:  nbgp, nxl, nxr, nyn, nys, nzb, nzt
+        ONLY:  nbgp, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
 
     USE kinds
@@ -204,6 +208,6 @@
     
     USE netcdf_interface,                                                      &
-        ONLY:  id_set_mask, id_var_domask, id_var_time_mask, nc_stat,          &
-               netcdf_data_format, netcdf_handle_error
+        ONLY:  fill_value, id_set_mask, id_var_domask, id_var_time_mask,       &
+               nc_stat, netcdf_data_format, netcdf_handle_error
     
     USE particle_attributes,                                                   &
@@ -218,7 +222,5 @@
     USE salsa_mod,                                                             &
         ONLY:  salsa_data_output_mask     
-        
-    USE surface_mod,                                                           &
-        ONLY :  get_topography_top_index_ji
+
 
     IMPLICIT NONE
@@ -232,4 +234,6 @@
     INTEGER(iwp) ::  j                       !< loop index
     INTEGER(iwp) ::  k                       !< loop index
+    INTEGER(iwp) ::  im                      !< loop index for masked variables
+    INTEGER(iwp) ::  jm                      !< loop index for masked variables
     INTEGER(iwp) ::  kk                      !< vertical index
     INTEGER(iwp) ::  mid                     !< masked output running index
@@ -237,5 +241,5 @@
     INTEGER(iwp) ::  netcdf_data_format_save !< value of netcdf_data_format
     INTEGER(iwp) ::  sender                  !< PE id of sending PE
-    INTEGER(iwp) ::  topo_top_ind            !< k index of highest horizontal surface
+    INTEGER(iwp) ::  ktt                     !< k index of highest terrain surface
     INTEGER(iwp) ::  ind(6)                  !< index limits (lower/upper bounds) of array 'local_2d'
 
@@ -309,8 +313,4 @@
        grid = 's'
 !
-!--    Set flag to steer output of radiation, land-surface, or user-defined
-!--    quantities
-       found = .FALSE.
-!
 !--    Store the variable chosen.
        resorted = .FALSE.
@@ -370,14 +370,16 @@
                    DO  i = 1, mask_size_l(mid,1)
                       DO  j = 1, mask_size_l(mid,2)
-!
-!--                      Get k index of highest horizontal surface
-                         topo_top_ind =  &
-                            get_topography_top_index_ji( mask_j(mid,j),  &
-                                                         mask_i(mid,i),  &
-                                                         grid )
+!--                      Get k index of the highest terraing surface
+                         im = mask_i(mid,i)
+                         jm = mask_j(mid,j)
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
-                            kk = MIN( topo_top_ind+mask_k(mid,k), nzt+1 )
-                            local_pf(i,j,k) =  &
-                               tend(kk,mask_j(mid,j),mask_i(mid,i))
+                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!--                         Set value if not in building
+                            IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                               local_pf(i,j,k) = fill_value
+                            ELSE
+                               local_pf(i,j,k) =  tend(kk,jm,im)
+                            ENDIF
                          ENDDO
                       ENDDO
@@ -436,14 +438,16 @@
                    DO  i = 1, mask_size_l(mid,1)
                       DO  j = 1, mask_size_l(mid,2)
-!
-!--                      Get k index of highest horizontal surface
-                         topo_top_ind =  &
-                            get_topography_top_index_ji( mask_j(mid,j),  &
-                                                         mask_i(mid,i),  &
-                                                         grid )
+!--                      Get k index of the highest terraing surface
+                         im = mask_i(mid,i)
+                         jm = mask_j(mid,j)
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
-                            kk = MIN( topo_top_ind+mask_k(mid,k), nzt+1 )
-                            local_pf(i,j,k) =  &
-                               tend(kk,mask_j(mid,j),mask_i(mid,i))
+                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!--                         Set value if not in building
+                            IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                               local_pf(i,j,k) = fill_value
+                            ELSE
+                               local_pf(i,j,k) =  tend(kk,jm,im)
+                            ENDIF
                          ENDDO
                       ENDDO
@@ -477,16 +481,16 @@
                       DO  i = 1, mask_size_l(mid,1)
                          DO  j = 1, mask_size_l(mid,2)
-!
-!--                         Get k index of highest horizontal surface
-                            topo_top_ind =  &
-                               get_topography_top_index_ji( mask_j(mid,j),  &
-                                                            mask_i(mid,i),  &
-                                                            grid )
+!--                         Get k index of the highest terraing surface
+                            im = mask_i(mid,i)
+                            jm = mask_j(mid,j)
+                            ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                             DO  k = 1, mask_size_l(mid,3)
-                               kk = MIN( topo_top_ind+mask_k(mid,k), nzt+1 )
-                               local_pf(i,j,k) =  &
-                                    pt(kk,mask_j(mid,j),mask_i(mid,i) ) &
-                                    + lv_d_cp * d_exner(kk) *           &
-                                      ql(kk,mask_j(mid,j),mask_i(mid,i))
+                               kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!--                            Set value if not in building
+                               IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                                  local_pf(i,j,k) = fill_value
+                               ELSE
+                                  local_pf(i,j,k) = pt(kk,jm,im) + lv_d_cp * d_exner(kk) * ql(kk,jm,im)
+                               ENDIF
                             ENDDO
                          ENDDO
@@ -571,14 +575,16 @@
                    DO  i = 1, mask_size_l(mid,1)
                       DO  j = 1, mask_size_l(mid,2)
-!
-!--                      Get k index of highest horizontal surface
-                         topo_top_ind =  &
-                            get_topography_top_index_ji( mask_j(mid,j),  &
-                                                         mask_i(mid,i),  &
-                                                         grid )
+!--                      Get k index of the highest terraing surface
+                         im = mask_i(mid,i)
+                         jm = mask_j(mid,j)
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
-                            kk = MIN( topo_top_ind+mask_k(mid,k), nzt+1 )
-                            local_pf(i,j,k) =  &
-                               tend(kk,mask_j(mid,j),mask_i(mid,i))
+                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!--                         Set value if not in building
+                            IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                               local_pf(i,j,k) = fill_value
+                            ELSE
+                               local_pf(i,j,k) = tend(kk,jm,im)
+                            ENDIF
                          ENDDO
                       ENDDO
@@ -608,15 +614,16 @@
                    DO  i = 1, mask_size_l(mid,1)
                       DO  j = 1, mask_size_l(mid,2)
-!
-!--                      Get k index of highest horizontal surface
-                         topo_top_ind =  &
-                            get_topography_top_index_ji( mask_j(mid,j),  &
-                                                         mask_i(mid,i),  &
-                                                         grid )
+!--                      Get k index of the highest terraing surface
+                         im = mask_i(mid,i)
+                         jm = mask_j(mid,j)
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
-                            kk = MIN( topo_top_ind+mask_k(mid,k), nzt+1 )
-                            local_pf(i,j,k) =  &
-                                 q(kk,mask_j(mid,j),mask_i(mid,i)) -  &
-                                 ql(kk,mask_j(mid,j),mask_i(mid,i))
+                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!--                         Set value if not in building
+                            IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                               local_pf(i,j,k) = fill_value
+                            ELSE
+                               local_pf(i,j,k) = q(kk,jm,im) - ql(kk,jm,im)
+                            ENDIF
                          ENDDO
                       ENDDO
@@ -686,13 +693,16 @@
 
           CASE DEFAULT
-
+!
+!--          Set flag to steer output of radiation, land-surface, or user-defined
+!--          quantities
+             found = .FALSE.
 !
 !--          Radiation quantity
-             IF ( radiation )  THEN
+             IF ( .NOT. found  .AND. radiation )  THEN
                 CALL radiation_data_output_mask(av, domask(mid,av,ivar), found,&
                                                 local_pf, mid )
              ENDIF
 
-             IF ( air_chemistry )  THEN
+             IF ( .NOT. found  .AND. air_chemistry )  THEN
                 CALL chem_data_output_mask(av, domask(mid,av,ivar), found,     &
                                            local_pf, mid )
@@ -700,9 +710,11 @@
 !
 !--          Check for diagnostic quantities
-             CALL doq_output_mask( av, domask(mid,av,ivar), found, local_pf,   &
+             IF ( .NOT. found )  THEN
+                CALL doq_output_mask( av, domask(mid,av,ivar), found, local_pf,   &
                                    mid)
+             ENDIF
 !
 !--          SALSA quantities
-             IF (  salsa )  THEN
+             IF ( .NOT. found .AND. salsa )  THEN
                 CALL salsa_data_output_mask( av, domask(mid,av,ivar), found,   &
                                              local_pf, mid )
@@ -745,17 +757,16 @@
              DO  i = 1, mask_size_l(mid,1)
                 DO  j = 1, mask_size_l(mid,2)
-!
-!--                Get k index of highest horizontal surface
-                   topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), &
-                                                               mask_i(mid,i), &
-                                                               grid )
-!
-!--                Save output array
+!--                Get k index of the highest terraing surface
+                   im = mask_i(mid,i)
+                   jm = mask_j(mid,j)
+                   ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                    DO  k = 1, mask_size_l(mid,3)
-                      local_pf(i,j,k) = to_be_resorted(                       &
-                                             MIN( topo_top_ind+mask_k(mid,k), &
-                                                  nzt+1 ),                    &
-                                             mask_j(mid,j),                   &
-                                             mask_i(mid,i)                     )
+                      kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!--                   Set value if not in building
+                      IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                         local_pf(i,j,k) = fill_value
+                      ELSE
+                         local_pf(i,j,k) = to_be_resorted(kk,jm,im)
+                      ENDIF
                    ENDDO
                 ENDDO
Index: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90
===================================================================
--- palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4166)
+++ palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90	(revision 4167)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Changed behaviour of masked output over surface to follow terrain and ignore 
+! buildings (J.Resler, T.Gronemeier)
+! 
+! 4157 2019-08-14 09:19:12Z suehring
 ! Initialization restructured, in order to work also when data output during
 ! spin-up is enabled. 
@@ -86,5 +90,5 @@
 ! 
     USE kinds
-! 
+!
 !     USE land_surface_model_mod,                                                &
 !         ONLY:  zs
@@ -694,23 +698,23 @@
         
     USE indices
-    
-    USE surface_mod,                                                           &
-        ONLY:  get_topography_top_index_ji
- 
+
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  variable  !< 
-    CHARACTER (LEN=5) ::  grid      !< flag to distinquish between staggered grids
-
-    INTEGER(iwp) ::  av           !< index indicating averaged or instantaneous output
-    INTEGER(iwp) ::  flag_nr      !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
-    INTEGER(iwp) ::  i            !< index variable along x-direction
-    INTEGER(iwp) ::  j            !< index variable along y-direction
-    INTEGER(iwp) ::  k            !< index variable along z-direction
-    INTEGER(iwp) ::  mid          !< masked output running index 
-    INTEGER(iwp) ::  topo_top_ind !< k index of highest horizontal surface
-
-    LOGICAL ::  found             !< true if variable is in list
-    LOGICAL ::  resorted          !< true if array is resorted
+    CHARACTER (LEN=*) ::  variable   !<
+    CHARACTER (LEN=5) ::  grid       !< flag to distinquish between staggered grids
+
+    INTEGER(iwp) ::  av              !< index indicating averaged or instantaneous output
+    INTEGER(iwp) ::  flag_nr         !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
+    INTEGER(iwp) ::  i               !< index variable along x-direction
+    INTEGER(iwp) ::  j               !< index variable along y-direction
+    INTEGER(iwp) ::  k               !< index variable along z-direction
+    INTEGER(iwp) ::  im              !< loop index for masked variables
+    INTEGER(iwp) ::  jm              !< loop index for masked variables
+    INTEGER(iwp) ::  kk              !< masked output index variable along z-direction
+    INTEGER(iwp) ::  mid             !< masked output running index
+    INTEGER(iwp) ::  ktt             !< k index of highest horizontal surface
+
+    LOGICAL      ::  found           !< true if variable is in list
+    LOGICAL      ::  resorted        !< true if array is resorted
 
     REAL(wp),                                                                  &
@@ -718,4 +722,6 @@
           local_pf   !< 
     REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
+
+    REAL(wp), PARAMETER   ::  fill_value = -9999.0_wp       !< value for the _FillValue attribute
 
     flag_nr  = 0
@@ -788,16 +794,18 @@
              DO  j = 1, mask_size_l(mid,2)
 !
-!--             Get k index of highest horizontal surface
-                topo_top_ind = get_topography_top_index_ji( mask_j(mid,j),     &
-                                                            mask_i(mid,i),     &
-                                                            grid )
-!
-!--             Save output array
+!--             Get k index of the highest terraing surface
+                im = mask_i(mid,i)
+                jm = mask_j(mid,j)
+                ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), &
+                              DIM = 1 ) - 1
                 DO  k = 1, mask_size_l(mid,3)
-                   local_pf(i,j,k) = to_be_resorted(                           &
-                                          MIN( topo_top_ind+mask_k(mid,k),     &
-                                               nzt+1 ),                        &
-                                          mask_j(mid,j),                       &
-                                          mask_i(mid,i)                     )
+                   kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
+!--                Set value if not in building
+                   IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                      local_pf(i,j,k) = fill_value
+                   ELSE
+                      local_pf(i,j,k) = to_be_resorted(kk,jm,im)
+                   ENDIF
                 ENDDO
              ENDDO
Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4166)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4167)
@@ -28,4 +28,8 @@
 ! -----------------
 ! $Id$
+! Changed behaviour of masked output over surface to follow terrain and ignore 
+! buildings (J.Resler, T.Gronemeier)
+! 
+! 4157 2019-08-14 09:19:12Z suehring
 ! Give informative message on raytracing distance only by core zero
 ! 
@@ -741,4 +745,6 @@
 
     CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
+
+    REAL(wp), PARAMETER ::  fill_value = -9999.0_wp       !< value for the _FillValue attribute
 
 !
@@ -10137,6 +10143,4 @@
     LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
 
-    REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
-
     REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
 
@@ -10536,6 +10540,4 @@
     LOGICAL      ::  found       !<
 
-    REAL(wp)     ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
-
     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
 
@@ -11105,12 +11107,13 @@
     CHARACTER (LEN=*) ::  variable   !< 
 
-    CHARACTER(LEN=5) ::  grid        !< flag to distinquish between staggered grids
-
     INTEGER(iwp) ::  av              !< 
     INTEGER(iwp) ::  i               !< 
     INTEGER(iwp) ::  j               !< 
-    INTEGER(iwp) ::  k               !< 
+    INTEGER(iwp) ::  k               !<
+    INTEGER(iwp) ::  im              !< loop index for masked variables
+    INTEGER(iwp) ::  jm              !< loop index for masked variables
+    INTEGER(iwp) ::  kk              !< 
     INTEGER(iwp) ::  mid             !< masked output running index
-    INTEGER(iwp) ::  topo_top_ind    !< k index of highest horizontal surface
+    INTEGER(iwp) ::  ktt             !< k index of highest terrain surface
 
     LOGICAL ::  found                !< true if output array was found
@@ -11124,8 +11127,6 @@
     REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which needs to be resorted for output
 
-
-    found    = .TRUE.
-    grid     = 's'
     resorted = .FALSE.
+    found = .TRUE.
 
     SELECT CASE ( TRIM( variable ) )
@@ -11214,16 +11215,18 @@
              DO  j = 1, mask_size_l(mid,2)
 !
-!--             Get k index of highest horizontal surface
-                topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), &
-                                                            mask_i(mid,i), &
-                                                            grid )
-!
-!--             Save output array
+!--             Get k index of the highest terraing surface
+                im = mask_i(mid,i)
+                jm = mask_j(mid,j)
+                ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), &
+                              DIM = 1 ) - 1
                 DO  k = 1, mask_size_l(mid,3)
-                   local_pf(i,j,k) = to_be_resorted(                       &
-                                          MIN( topo_top_ind+mask_k(mid,k), &
-                                               nzt+1 ),                    &
-                                          mask_j(mid,j),                   &
-                                          mask_i(mid,i)                     )
+                   kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
+!--                Set value if not in building
+                   IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                      local_pf(i,j,k) = fill_value
+                   ELSE
+                      local_pf(i,j,k) = to_be_resorted(kk,jm,im)
+                   ENDIF
                 ENDDO
              ENDDO
Index: palm/trunk/SOURCE/salsa_mod.f90
===================================================================
--- palm/trunk/SOURCE/salsa_mod.f90	(revision 4166)
+++ palm/trunk/SOURCE/salsa_mod.f90	(revision 4167)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Changed behaviour of masked output over surface to follow terrain and ignore 
+! buildings (J.Resler, T.Gronemeier)
+! 
+! 4131 2019-08-02 11:06:18Z monakurppa
 ! - Add "salsa_" before each salsa output variable
 ! - Add a possibility to output the number (salsa_N_UFP) and mass concentration
@@ -256,4 +260,6 @@
                                              !< 2 = autumn (no harvest yet), 3 = late autumn
                                              !< (already frost), 4 = winter, 5 = transitional spring
+
+    REAL(wp), PARAMETER ::  fill_value = -9999.0_wp    !< value for the _FillValue attribute
 !
 !-- Universal constants
@@ -10393,5 +10399,4 @@
                                           !< depositing in the alveolar (or tracheobronchial)
                                           !< region of the lung. Depends on the particle size
-    REAL(wp) ::  fill_value = -9999.0_wp  !< value for the _FillValue attribute
     REAL(wp) ::  mean_d                   !< Particle diameter in micrometres
     REAL(wp) ::  temp_bin                 !< temporary array for calculating output variables
@@ -10808,5 +10813,4 @@
                                           !< depositing in the alveolar (or tracheobronchial)
                                           !< region of the lung. Depends on the particle size
-    REAL(wp) ::  fill_value = -9999.0_wp  !< value for the _FillValue attribute
     REAL(wp) ::  mean_d                   !< Particle diameter in micrometres
     REAL(wp) ::  temp_bin                 !< temporary array for calculating output variables
@@ -11174,7 +11178,4 @@
         ONLY:  mask_i, mask_j, mask_k, mask_size_l, mask_surface, nz_do3d
 
-    USE surface_mod,                                                                               &
-        ONLY:  get_topography_top_index_ji
-
     IMPLICIT NONE
 
@@ -11191,6 +11192,9 @@
     INTEGER(iwp) ::  j              !< loop index in y-direction
     INTEGER(iwp) ::  k              !< loop index in z-direction
+    INTEGER(iwp) ::  im             !< loop index for masked variables
+    INTEGER(iwp) ::  jm             !< loop index for masked variables
+    INTEGER(iwp) ::  kk             !< loop index for masked output in z-direction
     INTEGER(iwp) ::  mid            !< masked output running index
-    INTEGER(iwp) ::  topo_top_ind   !< k index of highest horizontal surface
+    INTEGER(iwp) ::  ktt            !< k index of highest terrain surface
 
     LOGICAL ::  found      !<
@@ -11231,10 +11235,18 @@
                 DO  i = 1, mask_size_l(mid,1)
                    DO  j = 1, mask_size_l(mid,2)
-                      topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i),    &
-                                                                  grid )
+!
+!--                   Get k index of the highest terraing surface
+                      im = mask_i(mid,i)
+                      jm = mask_j(mid,j)
+                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                       DO  k = 1, mask_size_l(mid,3)
-                         local_pf(i,j,k) = aerosol_number(ib)%conc(MIN( topo_top_ind+mask_k(mid,k),&
-                                                                        nzt+1 ),                   &
-                                                                   mask_j(mid,j), mask_i(mid,i) )
+                         kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
+!--                      Set value if not in building
+                         IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                            local_pf(i,j,k) = fill_value
+                         ELSE
+                            local_pf(i,j,k) = aerosol_number(ib)%conc(kk,jm,im)
+                         ENDIF
                       ENDDO
                    ENDDO
@@ -11277,9 +11289,18 @@
                 DO  i = 1, mask_size_l(mid,1)
                    DO  j = 1, mask_size_l(mid,2)
-                      topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i),    &
-                                                                  grid )
+!
+!--                   Get k index of the highest terraing surface
+                      im = mask_i(mid,i)
+                      jm = mask_j(mid,j)
+                      ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                       DO  k = 1, mask_size_l(mid,3)
-                         local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ),         &
-                                                 mask_j(mid,j), mask_i(mid,i) )
+                         kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
+!--                      Set value if not in building
+                         IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                            local_pf(i,j,k) = fill_value
+                         ELSE
+                            local_pf(i,j,k) = tend(kk,jm,im)
+                         ENDIF
                       ENDDO
                    ENDDO
@@ -11347,9 +11368,18 @@
                    DO  i = 1, mask_size_l(mid,1)
                       DO  j = 1, mask_size_l(mid,2)
-                         topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), &
-                                                                     grid )
+!
+!--                      Get k index of the highest terraing surface
+                         im = mask_i(mid,i)
+                         jm = mask_j(mid,j)
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
-                            local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ),      &
-                                                    mask_j(mid,j), mask_i(mid,i) )
+                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
+!--                         Set value if not in building
+                            IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                               local_pf(i,j,k) = fill_value
+                            ELSE
+                               local_pf(i,j,k) = tend(kk,jm,im)
+                            ENDIF
                          ENDDO
                       ENDDO
@@ -11387,9 +11417,18 @@
                    DO  i = 1, mask_size_l(mid,1)
                       DO  j = 1, mask_size_l(mid,2)
-                         topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), &
-                                                                     grid )
+!
+!--                      Get k index of the highest terraing surface
+                         im = mask_i(mid,i)
+                         jm = mask_j(mid,j)
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
-                            local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ),      &
-                                                    mask_j(mid,j), mask_i(mid,i) )
+                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
+!--                         Set value if not in building
+                            IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                               local_pf(i,j,k) = fill_value
+                            ELSE
+                               local_pf(i,j,k) = tend(kk,jm,im)
+                            ENDIF
                          ENDDO
                       ENDDO
@@ -11425,9 +11464,18 @@
                    DO  i = 1, mask_size_l(mid,1)
                       DO  j = 1, mask_size_l(mid,2)
-                         topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), &
-                                                                     grid )
+!
+!--                      Get k index of the highest terraing surface
+                         im = mask_i(mid,i)
+                         jm = mask_j(mid,j)
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
-                            local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ),      &
-                                                    mask_j(mid,j), mask_i(mid,i) )
+                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
+!--                         Set value if not in building
+                            IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                               local_pf(i,j,k) = fill_value
+                            ELSE
+                               local_pf(i,j,k) = tend(kk,jm,im)
+                            ENDIF
                          ENDDO
                       ENDDO
@@ -11467,9 +11515,18 @@
                    DO  i = 1, mask_size_l(mid,1)
                       DO  j = 1, mask_size_l(mid,2)
-                         topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), &
-                                                                     grid )
+!
+!--                      Get k index of the highest terraing surface
+                         im = mask_i(mid,i)
+                         jm = mask_j(mid,j)
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
-                            local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ),      &
-                                                    mask_j(mid,j), mask_i(mid,i) )
+                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
+!--                         Set value if not in building
+                            IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                               local_pf(i,j,k) = fill_value
+                            ELSE
+                               local_pf(i,j,k) = tend(kk,jm,im)
+                            ENDIF
                          ENDDO
                       ENDDO
@@ -11509,9 +11566,18 @@
                    DO  i = 1, mask_size_l(mid,1)
                       DO  j = 1, mask_size_l(mid,2)
-                         topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), &
-                                                                     grid )
+!
+!--                      Get k index of the highest terraing surface
+                         im = mask_i(mid,i)
+                         jm = mask_j(mid,j)
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
-                            local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ),      &
-                                                    mask_j(mid,j), mask_i(mid,i) )
+                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
+!--                         Set value if not in building
+                            IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                               local_pf(i,j,k) = fill_value
+                            ELSE
+                               local_pf(i,j,k) = tend(kk,jm,im)
+                            ENDIF
                          ENDDO
                       ENDDO
@@ -11551,9 +11617,18 @@
                    DO  i = 1, mask_size_l(mid,1)
                       DO  j = 1, mask_size_l(mid,2)
-                         topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), &
-                                                                     grid )
+!
+!--                      Get k index of the highest terraing surface
+                         im = mask_i(mid,i)
+                         jm = mask_j(mid,j)
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
-                            local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ),      &
-                                                    mask_j(mid,j), mask_i(mid,i) )
+                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
+!--                         Set value if not in building
+                            IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                               local_pf(i,j,k) = fill_value
+                            ELSE
+                               local_pf(i,j,k) = tend(kk,jm,im)
+                            ENDIF
                          ENDDO
                       ENDDO
@@ -11594,9 +11669,18 @@
                    DO  i = 1, mask_size_l(mid,1)
                       DO  j = 1, mask_size_l(mid,2)
-                         topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), &
-                                                                     grid )
+!
+!--                      Get k index of the highest terraing surface
+                         im = mask_i(mid,i)
+                         jm = mask_j(mid,j)
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
-                            local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ),      &
-                                                    mask_j(mid,j), mask_i(mid,i) )
+                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
+!--                         Set value if not in building
+                            IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                               local_pf(i,j,k) = fill_value
+                            ELSE
+                               local_pf(i,j,k) = tend(kk,jm,im)
+                            ENDIF
                          ENDDO
                       ENDDO
@@ -11641,9 +11725,18 @@
                    DO  i = 1, mask_size_l(mid,1)
                       DO  j = 1, mask_size_l(mid,2)
-                         topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), &
-                                                                     grid )
+!
+!--                      Get k index of the highest terraing surface
+                         im = mask_i(mid,i)
+                         jm = mask_j(mid,j)
+                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                          DO  k = 1, mask_size_l(mid,3)
-                            local_pf(i,j,k) = tend( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ),      &
-                                                    mask_j(mid,j), mask_i(mid,i) )
+                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!
+!--                         Set value if not in building
+                            IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                               local_pf(i,j,k) = fill_value
+                            ELSE
+                               local_pf(i,j,k) =  tend(kk,jm,im)
+                            ENDIF
                          ENDDO
                       ENDDO
@@ -11661,5 +11754,5 @@
     ENDIF
 
-    IF ( .NOT. resorted )  THEN
+    IF ( found  .AND.  .NOT. resorted )  THEN
        IF ( .NOT. mask_surface(mid) )  THEN
 !
@@ -11677,12 +11770,16 @@
           DO  i = 1, mask_size_l(mid,1)
              DO  j = 1, mask_size_l(mid,2)
-!
-!--             Get k index of highest horizontal surface
-                topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), mask_i(mid,i), grid )
-!
-!--             Save output array
+!--             Get k index of the highest terraing surface
+                im = mask_i(mid,i)
+                jm = mask_j(mid,j)
+                ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
                 DO  k = 1, mask_size_l(mid,3)
-                   local_pf(i,j,k) = to_be_resorted( MIN( topo_top_ind+mask_k(mid,k), nzt+1 ),     &
-                                                     mask_j(mid,j), mask_i(mid,i) )
+                   kk = MIN( ktt+mask_k(mid,k), nzt+1 )
+!--                Set value if not in building
+                   IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) )  THEN
+                      local_pf(i,j,k) = fill_value
+                   ELSE
+                      local_pf(i,j,k) = to_be_resorted(kk,jm,im)
+                   ENDIF
                 ENDDO
              ENDDO
