Index: palm/trunk/SOURCE/salsa_mod.f90
===================================================================
--- palm/trunk/SOURCE/salsa_mod.f90	(revision 4131)
+++ 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
